Les données spatiales

class: right, middle, title-slide

# Les données spatiales
## <code>{sf}</code>  <code>{ggplot2}</code>
### Nina SCHIETTEKATTE (Nicolas CASAJUS)
### .inst[Mardi 3 novembre 2020]

---

class: inverse, center, middle

## Introduction

---

## Qu'est-ce qu'un objet spatial ?

&nbsp;Deux catégories :

**1)** Les données vectorielles

.center[![:scale 100%](img/vectors.png)]

---

## Qu'est-ce qu'un objet spatial ?

&nbsp;Deux catégories :

**2)** Les données matricielles

.center[![:scale 100%](img/rasters.png)]

---

## Le système de coordonnées

&nbsp;Un objet spatial se représente
dans l'espace selon un référentiel spatial : le **système de coordonnées** (**CRS** pour _Coordinates Reference System_)

&nbsp;Deux types de systèmes de coordonnées existent :

**1)** Les systèmes géographiques (ou non projetés) [en degrés]

.center[![:scale 100%](img/projections01.png)]

---

## Le système de coordonnées

&nbsp;Deux types de systèmes de coordonnées existent :

**2)** Les systèmes projetés [en mètres]

.center[![:scale 100%](img/projections02.png)]

&nbsp;Le choix du CRS peut être crucial et dépend souvent de l'information que l'on souhaite représenter

---

## Le système de coordonnées

Sous le CRS s'exprime selon le standard `proj4string` défini par le projet [**PROJ**](https://proj.org/)

.small[`+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0`]

Chaque CRS possède un **SRID** (_Spatial Reference Identifier_) ou **EPSG**
(_European Petroleum Survey Group_) : un identifiant unique qui nous évite de
devoir écrire à la main le CRS au format `proj4string`

`+init=epsg:4326` est l'identifiant du CRS ci-dessus

.small[`+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0`]

&nbsp;Où trouver le bon CRS ? Réponse : [**Spatial Reference**](https://spatialreference.org/)

---
class: inverse, center, middle

## Objets vectoriels

---

## Les packages - Vecteurs

Package           | Fonctionnalités
------------------|------------------
`{sp}`            | Classes et méthodes pour couches vectorielles
`{rgdal}`         | Lecture et écriture de couches spatiales
`{rgeos}`         | Opérations spatiales sur couches vectorielles
`{maptools}`      | Opérations spatiales sur couches vectorielles
`{sf}`            | Le petit dernier - Une vraie bête !

&nbsp;`{sf}` est devenu la "norme" : `{sf} = {sp} + {rgdal} + {rgeos}`

---

## Pourquoi `{sf}` ?

- Implémentation sous des _Simple Features_
 - Standard reposant sur la norme ISO 19125 (stockage et accès)
 - Remplace le format d'ESRI.small[TM], vous savez le shapefile !?

- 3 packages en 1

- Syntaxe (_snake_case_) des fonctions simples (`{st_*()}`)

- Objets stockés dans des `data.frame` ou des `tibble`
  - avec la géométrie stockée dans des `list-column`

- Très grandes performances comparées à `{sp}`

- Pipe-friendly et intégration au tidyverse (dont `{ggplot2}`)

&nbsp;Edzer Pebesma le recommende 😄

---

## Structure d'un objet `sf`

.center[![:scale 100%](img/str_sf.png)]

---

## Les _Simple feature geometry_ `sfg`

.center[![:scale 60%](img/sf_obj.png)]

---

## Installation

```r
## install.packages("sf")
```

---

## Création de points spatiaux

Soit les villes avec les coordonnées suivantes :

```r
foix <- c(1.6053807, 42.9638998)
carcassonne <- c(2.3491069, 43.2130358)
rodez <- c(2.5728419, 44.3516207)
```

Convertissons ces vecteurs en points spatiaux :

```r
foix <- sf::st_point(foix)
carcassonne <- sf::st_point(carcassonne)
rodez <- sf::st_point(rodez)
```

Et regroupons-les en un seul _spatial column_ en définissant le CRS :

```r
villes <- sf::st_sfc(
 list(foix, carcassonne, rodez),
 crs = 4326
)
class(villes)
```

```
## [1] "sfc_POINT" "sfc"
```

---

## Création de points spatiaux

Ajoutons une table d'attributs et convertissons l'objet en _simple  feature_ :

```r
datas <- data.frame(
 ville = c("Foix", "Carcassonne", "Rodez"),
 population = c(9613, 45895, 23739)
)

villes <- sf::st_sf(datas, geom = villes)
class(villes)
```

```
## [1] "sf"         "data.frame"
```

&nbsp;Quid des lignes et polygones ?

&nbsp;Voir les fonctions `st_linestring()` et `st_polygon()`

---

## Création de points spatiaux

Bien souvent, vous importerez directement un `data.frame` que vous souhaiterez convertir en _simple feature_

```r
## URL et nom du fichier
url <- "https://raw.githubusercontent.com/FRBCesab/datatoolbox/master/data/"
filename <- "occitanie_prefectures.csv"

## Téléchargement du fichier
download.file(
  url      = paste0(url, filename),
  destfile = filename
)

## Extraction du ZIP
unzip(zipfile = filename)
```

---

## Création de points spatiaux

```r
## Ouverture du csv
(prefectures <- readr::read_csv("occitanie_prefectures.csv"))
```

```
## # A tibble: 13 × 5
## departement prefecture population longitude latitude
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 Ariège Foix 9613 1.61 43.0
## 2 Aude Carcassonne 45895 2.35 43.2
## 3 Aveyron Rodez 23739 2.57 44.4
## 4 Gard Nîmes 151001 4.36 43.8
## 5 Haute-Garonne Toulouse 475438 1.44 43.6
## 6 Gers Auch 21618 0.585 43.6
## 7 Hérault Montpellier 281613 3.88 43.6
## 8 Lot Cahors 19405 1.44 44.4
## 9 Lozère Mende 11860 3.50 44.5
## 10 Hautes-Pyrénées Tarbes 40318 0.0781 43.2
## 11 Pyrénées-Orientales Perpignan 121875 2.90 42.7
## 12 Tarn Albi 49024 2.15 43.9
## 13 Tarn-et-Garonne Montauban 60444 1.35 44.0
```

---

## Création de points spatiaux

```r
library("magrittr")
prefectures <- sf::st_as_sf(
 x = prefectures,
 coords = c("longitude", "latitude"),
 crs = 4326
)
prefectures %>% head(6)
```

```
## Simple feature collection with 6 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 0.5850507 ymin: 42.9639 xmax: 4.360069 ymax: 44.35162
## Geodetic CRS: WGS 84
## # A tibble: 6 × 4
## departement prefecture population geometry
## <chr> <chr> <dbl> <POINT [°]>
## 1 Ariège Foix 9613 (1.605381 42.9639)
## 2 Aude Carcassonne 45895 (2.349107 43.21304)
## 3 Aveyron Rodez 23739 (2.572842 44.35162)
## 4 Gard Nîmes 151001 (4.360069 43.83742)
## 5 Haute-Garonne Toulouse 475438 (1.444247 43.60446)
## 6 Gers Auch 21618 (0.5850507 43.64636)
```

---

## Une première carte

Les packages `{rnaturalearth}` et `{rnaturalearthdata}` proposent différentes
couches vectorielles, dont les limites administratives des pays du monde.

Installons ces packages :

```r
## install.packages("rnaturalearth")
## install.packages("rnaturalearthdata")
```

Et chargeons le fond de carte du monde entier :

```r
world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
class(world)
```

```
## [1] "sf"         "data.frame"
```

---

## Une première carte

```r
ggplot(data = world) +
  geom_sf() +
  xlab("Longitude") + ylab("Latitude") + ggtitle("Carte du monde")
```

---

## Une première carte

```r
ggplot(data = world) +
  geom_sf(fill = "#49847b", colour = "#e1ddc0") +
  xlab("Longitude") + ylab("Latitude") + ggtitle("Carte du monde")
```

---

## Une première carte

```r
ggplot(data = world) +
  geom_sf(fill = "#49847b", colour = "#e1ddc0") +
  xlab("Longitude") + ylab("Latitude") + ggtitle("Carte du monde") +
  geom_sf(data = prefectures, colour = "black")
```

---

## Importation d'une couche vectorielle

Dans la plupart des cas, vous devrez importer une couche spatiale déjà prête.

Téléchargeons les départements de la région Occitanie

```r
## URL et nom du fichier
url <- "https://raw.githubusercontent.com/FRBCesab/datatoolbox/master/data/"
filename <- "occitanie.zip"

## Téléchargement du fichier
download.file(
  url      = paste0(url, filename),
  destfile = filename
)

## Extraction du ZIP
unzip(zipfile = filename)
```

---

## Importation d'une couche vectorielle

Que contient la couche ?

```r
sf::st_layers("occitanie.shp")
```

```
## Driver: ESRI Shapefile 
## Available layers:
##   layer_name geometry_type features fields
## 1  occitanie       Polygon       13     13
```

Importons-la :

```r
occitanie <- sf::st_read(dsn = ".", layer = "occitanie")
```

```
## Reading layer 'occitanie' from data source '/Users/nicolascasajus/Desktop/...'
## Simple feature collection with 13 features and 13 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -0.326875 ymin: 42.33344 xmax: 4.845344 ymax: 45.04467
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
```

---

## Importation d'une couche vectorielle

```r
occitanie %>% head(3)
```

```
## Simple feature collection with 3 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 0.8266816 ymin: 42.57232 xmax: 3.450981 ymax: 44.94122
## Geodetic CRS: WGS 84
## GID_0 NAME_0 GID_1 NAME_1 NL_NAME_1 GID_2 NAME_2 VARNAME_2 NL_NAME_2 TYPE_2 ENGTYPE_2
## 1 FRA France FRA.11_1 Occitanie <NA> FRA.11.1_1 Ariège <NA> <NA> Département Department
## 2 FRA France FRA.11_1 Occitanie <NA> FRA.11.2_1 Aude <NA> <NA> Département Department
## 3 FRA France FRA.11_1 Occitanie <NA> FRA.11.3_1 Aveyron <NA> <NA> Département Department
## CC_2 HASC_2 geometry
## 1 09 FR.AG MULTIPOLYGON (((1.494104 42...
## 2 11 FR.AD MULTIPOLYGON (((3.027638 42...
## 3 12 FR.AV MULTIPOLYGON (((3.236248 43...
```

---

```r
ggplot(data = occitanie) +
  geom_sf(fill = "#49847b", colour = "#e1ddc0") +
  xlab("Longitude") + ylab("Latitude") + ggtitle("OCCITANIE")
```

---

```r
ggplot(data = occitanie) +
  geom_sf(fill = "#49847b", colour = "#e1ddc0") +
  xlab("Longitude") + ylab("Latitude") + ggtitle("OCCITANIE") +
  geom_sf(data = prefectures, colour = "#3f3f3f")
```

---

```r
ggplot(data = occitanie) +
  geom_sf(fill = "#49847b", colour = "#e1ddc0") +
  xlab("Longitude") + ylab("Latitude") + ggtitle("OCCITANIE") +
  geom_sf(data = prefectures, colour = "#3f3f3f") +
  geom_sf_label(aes(label = NAME_2))
```

---

---

## Exportation d'une couche vectorielle

```r
sf::st_drivers()
```

```
          name                       long_name write  copy is_raster is_vector   vsi
        PCIDSK            PCIDSK Database File  TRUE FALSE      TRUE      TRUE  TRUE
        netCDF      Network Common Data Format  TRUE  TRUE      TRUE      TRUE FALSE
           PDF                  Geospatial PDF  TRUE  TRUE      TRUE      TRUE FALSE
       MBTiles                         MBTiles  TRUE  TRUE      TRUE      TRUE  TRUE
          EEDA           Earth Engine Data API FALSE FALSE     FALSE      TRUE FALSE
ESRI Shapefile                  ESRI Shapefile  TRUE FALSE     FALSE      TRUE  TRUE
  MapInfo File                    MapInfo File  TRUE FALSE     FALSE      TRUE  TRUE
       UK .NTF                         UK .NTF FALSE FALSE     FALSE      TRUE  TRUE
      OGR_SDTS                            SDTS FALSE FALSE     FALSE      TRUE  TRUE
           S57                  IHO S-57 (ENC)  TRUE FALSE     FALSE      TRUE  TRUE
           DGN                Microstation DGN  TRUE FALSE     FALSE      TRUE  TRUE
       OGR_VRT        VRT - Virtual Datasource FALSE FALSE     FALSE      TRUE  TRUE
           REC                    EPIInfo .REC FALSE FALSE     FALSE      TRUE FALSE
        Memory                          Memory  TRUE FALSE     FALSE      TRUE FALSE
           BNA                       Atlas BNA  TRUE FALSE     FALSE      TRUE  TRUE
           CSV    Comma Separated Value (.csv)  TRUE FALSE     FALSE      TRUE  TRUE
           GML Geography Markup Language (GML)  TRUE FALSE     FALSE      TRUE  TRUE
           GPX                             GPX  TRUE FALSE     FALSE      TRUE  TRUE
           KML   Keyhole Markup Language (KML)  TRUE FALSE     FALSE      TRUE  TRUE
       GeoJSON                         GeoJSON  TRUE FALSE     FALSE      TRUE  TRUE
           ...                             ...   ...   ...       ...       ...   ...
```

---

## Exportation d'une couche vectorielle

Au format ESRI.small[TM] :

```r
sf::st_write(
  obj    = prefectures,
  dsn    = ".",
  layer  = "prefectures",
  driver = "ESRI Shapefile"
)
```

---

## Le système de coordonnées

Le CRS de notre couche vectorielle (monde) est-il défini ?

```r
world
```
```
Simple feature collection with 241 features and 63 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
CRS:            +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
...
```

Supprimons-le :

```r
world <- world %>% sf::st_set_crs(NA)
```

---

## Le système de coordonnées

Et redéfinissons-le :

```r
world <- world %>% sf::st_set_crs(4326)
```
```
Simple feature collection with 241 features and 63 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
geographic CRS: WGS 84
...
```

Projetons notre couche dans un autre système (par .ex, la projection Lambert
azimuthal equal area centrée sur la France)

```r
prj <- paste(
 "+proj=laea +lat_0=48 +lon_0=5",
 "+x_0=4321000 +y_0=3210000",
 "+ellps=GRS80 +units=m +no_defs"
)

world_ortho <- world %>% sf::st_transform(crs = prj)
```

---

```r
ggplot(data = world_ortho) +
  geom_sf(fill = "#49847b", colour = "#e1ddc0")
```

---

```r
ggplot(data = world) +
  geom_sf(fill = "#49847b", colour = "#e1ddc0") +
  coord_sf(crs = prj)
```

---

## Manipulations spatiales : sf cheat sheet

.center[![:scale 90%](img/geom_ops.png)]

---

## Manipulations spatiales : sf cheat sheet

Mesures géométriques

Fonction R              | Détails
------------------------|------------------
`sf::st_length(x)`      | Longueur de lignes spatiales
`sf::st_area(x)`        | Superficie de polygones spatiaux
`sf::st_distance(x, y)` | Distance entre deux couches
`sf::st_bbox(x)`        | Etentue spatiale (emprise) d'une couche

Opérations géométriques unitaires

Fonction R              | Détails
------------------------|------------------
`sf::st_buffer(x)`      | Buffers autour d'une géométrie
`sf::st_boundary(x)`    | Limites d'une géométrie
`sf::st_convex_hull(x)` | Convex hull d'un ensemble de points
`sf::st_centroid(x)`    | Centroïde d'une géométrie

---

## Manipulations spatiales : sf cheat sheet

Opérations géométriques binaires

- `sf::st_intersects(x, y)`
- `sf::st_disjoint(x, y)`
- `sf::st_contains(x, y)`
- `sf::st_overlaps(x, y)`
- `...`

Autres opérations

- `sf::st_crop()`
- `sf::st_intersection()`
- `sf::st_difference()`
- `sf::st_contains()`
- `sf::st_union()`

```r
ls("package:sf")
```

---

## Interactivité

<div id="htmlwidget-6a6b458327077b27e559" style="width:100%;height:432px;" class="leaflet html-widget"></div>
<script type="application/json" data-for="htmlwidget-6a6b458327077b27e559">{"x":{"options":{"minZoom":1,"maxZoom":52,"crs":{"crsClass":"L.CRS.EPSG3857","code":null,"proj4def":null,"projectedBounds":null,"options":{}},"preferCanvas":false,"bounceAtZoomLimits":false,"maxBounds":[[[-90,-370]],[[90,370]]]},"calls":[{"method":"addProviderTiles","args":["CartoDB.Positron","CartoDB.Positron","CartoDB.Positron",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["CartoDB.DarkMatter","CartoDB.DarkMatter","CartoDB.DarkMatter",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["OpenStreetMap","OpenStreetMap","OpenStreetMap",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["Esri.WorldImagery","Esri.WorldImagery","Esri.WorldImagery",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["OpenTopoMap","OpenTopoMap","OpenTopoMap",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"createMapPane","args":["point",440]},{"method":"addFlatGeoBuf","args":["prefectures","prefectures",null,true,"mvFeatureId",{"radius":6,"stroke":true,"color":"#333333","weight":1,"opacity":0.9,"fill":true,"fillColor":"#6666FF","fillOpacity":0.6},{"className":"","pane":"point"},"mapview-popup",{"radius":{"to":[3,15],"from":[3,15]},"weight":{"to":[1,10],"from":[1,10]},"opacity":{"to":[0,1],"from":[0,1]},"fillOpacity":{"to":[0,1],"from":[0,1]}}]},{"method":"addScaleBar","args":[{"maxWidth":100,"metric":true,"imperial":true,"updateWhenIdle":true,"position":"bottomleft"}]},{"method":"addHomeButton","args":[0.0781021,42.6985304,4.3600687,44.5180226,true,"prefectures","Zoom to prefectures"," prefectures <\/strong>","bottomright"]},{"method":"addLayersControl","args":[["CartoDB.Positron","CartoDB.DarkMatter","OpenStreetMap","Esri.WorldImagery","OpenTopoMap"],"prefectures",{"collapsed":true,"autoZIndex":true,"position":"topleft"}]},{"method":"addLegend","args":[{"colors":["#6666FF"],"labels":["prefectures"],"na_color":null,"na_label":"NA","opacity":1,"position":"topright","type":"factor","title":"","extra":null,"layerId":null,"className":"info legend","group":"prefectures"}]}],"fitBounds":[42.6985304,0.0781021,44.5180226,4.3600687,[]]},"evals":[],"jsHooks":{"render":[{"code":"function(el, x, data) {\n return (\n function(el, x, data) {\n // get the leaflet map\n var map = this; //HTMLWidgets.find('#' + el.id);\n // we need a new div element because we have to handle\n // the mouseover output separately\n // debugger;\n function addElement () {\n // generate new div Element\n var newDiv = $(document.createElement('div'));\n // append at end of leaflet htmlwidget container\n $(el).append(newDiv);\n //provide ID and style\n newDiv.addClass('lnlt');\n newDiv.css({\n 'position': 'relative',\n 'bottomleft': '0px',\n 'background-color': 'rgba(255, 255, 255, 0.7)',\n 'box-shadow': '0 0 2px #bbb',\n 'background-clip': 'padding-box',\n 'margin': '0',\n 'padding-left': '5px',\n 'color': '#333',\n 'font': '9px/1.5 \"Helvetica Neue\", Arial, Helvetica, sans-serif',\n 'z-index': '700',\n });\n return newDiv;\n }\n\n\n // check for already existing lnlt class to not duplicate\n var lnlt = $(el).find('.lnlt');\n\n if(!lnlt.length) {\n lnlt = addElement();\n\n // grab the special div we generated in the beginning\n // and put the mousmove output there\n\n map.on('mousemove', function (e) {\n if (e.originalEvent.ctrlKey) {\n if (document.querySelector('.lnlt') === null) lnlt = addElement();\n lnlt.text(\n ' lon: ' + (e.latlng.lng).toFixed(5) +\n ' | lat: ' + (e.latlng.lat).toFixed(5) +\n ' | zoom: ' + map.getZoom() +\n ' | x: ' + L.CRS.EPSG3857.project(e.latlng).x.toFixed(0) +\n ' | y: ' + L.CRS.EPSG3857.project(e.latlng).y.toFixed(0) +\n ' | epsg: 3857 ' +\n ' | proj4: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs ');\n } else {\n if (document.querySelector('.lnlt') === null) lnlt = addElement();\n lnlt.text(\n ' lon: ' + (e.latlng.lng).toFixed(5) +\n ' | lat: ' + (e.latlng.lat).toFixed(5) +\n ' | zoom: ' + map.getZoom() + ' ');\n }\n });\n\n // remove the lnlt div when mouse leaves map\n map.on('mouseout', function (e) {\n var strip = document.querySelector('.lnlt');\n if( strip !==null) strip.remove();\n });\n\n };\n\n //$(el).keypress(67, function(e) {\n map.on('preclick', function(e) {\n if (e.originalEvent.ctrlKey) {\n if (document.querySelector('.lnlt') === null) lnlt = addElement();\n lnlt.text(\n ' lon: ' + (e.latlng.lng).toFixed(5) +\n ' | lat: ' + (e.latlng.lat).toFixed(5) +\n ' | zoom: ' + map.getZoom() + ' ');\n var txt = document.querySelector('.lnlt').textContent;\n console.log(txt);\n //txt.innerText.focus();\n //txt.select();\n setClipboardText('\"' + txt + '\"');\n }\n });\n\n }\n ).call(this.getMap(), el, x, data);\n}","data":null},{"code":"function(el, x, data) {\n return (function(el,x,data){\n var map = this;\n\n map.on('keypress', function(e) {\n console.log(e.originalEvent.code);\n var key = e.originalEvent.code;\n if (key === 'KeyE') {\n var bb = this.getBounds();\n var txt = JSON.stringify(bb);\n console.log(txt);\n\n setClipboardText('\\'' + txt + '\\'');\n }\n })\n }).call(this.getMap(), el, x, data);\n}","data":null}]}}</script>

```r
mapview::mapview(prefectures)
```

---

## Interactivité

<div id="htmlwidget-f4ed25483e186b4b2f5c" style="width:100%;height:432px;" class="leaflet html-widget"></div>
<script type="application/json" data-for="htmlwidget-f4ed25483e186b4b2f5c">{"x":{"options":{"minZoom":1,"maxZoom":52,"crs":{"crsClass":"L.CRS.EPSG3857","code":null,"proj4def":null,"projectedBounds":null,"options":{}},"preferCanvas":false,"bounceAtZoomLimits":false,"maxBounds":[[[-90,-370]],[[90,370]]]},"calls":[{"method":"addProviderTiles","args":["CartoDB.Positron","CartoDB.Positron","CartoDB.Positron",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["CartoDB.DarkMatter","CartoDB.DarkMatter","CartoDB.DarkMatter",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["OpenStreetMap","OpenStreetMap","OpenStreetMap",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["Esri.WorldImagery","Esri.WorldImagery","Esri.WorldImagery",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"addProviderTiles","args":["OpenTopoMap","OpenTopoMap","OpenTopoMap",{"errorTileUrl":"","noWrap":false,"detectRetina":false,"pane":"tilePane"}]},{"method":"createMapPane","args":["polygon",420]},{"method":"addFlatGeoBuf","args":["Département","Département",null,true,"NAME_2",{"radius":6,"stroke":true,"color":"#333333","weight":0.5,"opacity":0.9,"fill":true,"fillColor":null,"fillOpacity":0.6},{"className":"","pane":"polygon"},"mapview-popup",{"radius":{"to":[3,15],"from":[3,15]},"weight":{"to":[1,10],"from":[1,10]},"opacity":{"to":[0,1],"from":[0,1]},"fillOpacity":{"to":[0,1],"from":[0,1]}}]},{"method":"addScaleBar","args":[{"maxWidth":100,"metric":true,"imperial":true,"updateWhenIdle":true,"position":"bottomleft"}]},{"method":"addHomeButton","args":[-0.32687503,42.33344269,4.84534407,45.04466629,true,"Département","Zoom to Département"," Département <\/strong>","bottomright"]},{"method":"addLegend","args":[{"colors":["#4B0055","#46226A","#353E7C","#00588B","#007094","#008798","#009B95","#00AE8C","#00BE7D","#53CC67","#96D84B","#CDE030","#FDE333"],"labels":["Ariège","Aude","Aveyron","Gard","Gers","Haute-Garonne","Hautes-Pyrénées","Hérault","Lot","Lozère","Pyrénées-Orientales","Tarn","Tarn-et-Garonne"],"na_color":null,"na_label":"NA","opacity":1,"position":"topright","type":"factor","title":"Département","extra":null,"layerId":null,"className":"info legend","group":"Département"}]},{"method":"createMapPane","args":["point",440]},{"method":"addFlatGeoBuf","args":["prefectures","prefectures",null,true,"mvFeatureId",{"radius":"population","stroke":true,"color":"#333333","weight":1,"opacity":0.9,"fill":true,"fillColor":"#3F3F3F","fillOpacity":0.6},{"className":"","pane":"point"},"mapview-popup",{"radius":{"to":[3,15],"from":[9613,475438]},"weight":{"to":[1,10],"from":[1,10]},"opacity":{"to":[0,1],"from":[0,1]},"fillOpacity":{"to":[0,1],"from":[0,1]}}]},{"method":"addHomeButton","args":[0.0781021,42.6985304,4.3600687,44.5180226,true,"prefectures","Zoom to prefectures"," prefectures <\/strong>","bottomright"]},{"method":"addLayersControl","args":[["CartoDB.Positron","CartoDB.DarkMatter","OpenStreetMap","Esri.WorldImagery","OpenTopoMap"],["Département","prefectures"],{"collapsed":true,"autoZIndex":true,"position":"topleft"}]},{"method":"addLegend","args":[{"colors":["#3F3F3F"],"labels":["prefectures"],"na_color":null,"na_label":"NA","opacity":1,"position":"topright","type":"factor","title":"","extra":null,"layerId":null,"className":"info legend","group":"prefectures"}]},{"method":"addHomeButton","args":[-0.32687503,42.33344269,4.84534407,45.04466629,true,null,"Zoom to full extent","Zoom full<\/strong>","bottomleft"]}],"fitBounds":[42.33344269,-0.32687503,45.04466629,4.84534407,[]]},"evals":[],"jsHooks":{"render":[{"code":"function(el, x, data) {\n return (\n function(el, x, data) {\n // get the leaflet map\n var map = this; //HTMLWidgets.find('#' + el.id);\n // we need a new div element because we have to handle\n // the mouseover output separately\n // debugger;\n function addElement () {\n // generate new div Element\n var newDiv = $(document.createElement('div'));\n // append at end of leaflet htmlwidget container\n $(el).append(newDiv);\n //provide ID and style\n newDiv.addClass('lnlt');\n newDiv.css({\n 'position': 'relative',\n 'bottomleft': '0px',\n 'background-color': 'rgba(255, 255, 255, 0.7)',\n 'box-shadow': '0 0 2px #bbb',\n 'background-clip': 'padding-box',\n 'margin': '0',\n 'padding-left': '5px',\n 'color': '#333',\n 'font': '9px/1.5 \"Helvetica Neue\", Arial, Helvetica, sans-serif',\n 'z-index': '700',\n });\n return newDiv;\n }\n\n\n // check for already existing lnlt class to not duplicate\n var lnlt = $(el).find('.lnlt');\n\n if(!lnlt.length) {\n lnlt = addElement();\n\n // grab the special div we generated in the beginning\n // and put the mousmove output there\n\n map.on('mousemove', function (e) {\n if (e.originalEvent.ctrlKey) {\n if (document.querySelector('.lnlt') === null) lnlt = addElement();\n lnlt.text(\n ' lon: ' + (e.latlng.lng).toFixed(5) +\n ' | lat: ' + (e.latlng.lat).toFixed(5) +\n ' | zoom: ' + map.getZoom() +\n ' | x: ' + L.CRS.EPSG3857.project(e.latlng).x.toFixed(0) +\n ' | y: ' + L.CRS.EPSG3857.project(e.latlng).y.toFixed(0) +\n ' | epsg: 3857 ' +\n ' | proj4: +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs ');\n } else {\n if (document.querySelector('.lnlt') === null) lnlt = addElement();\n lnlt.text(\n ' lon: ' + (e.latlng.lng).toFixed(5) +\n ' | lat: ' + (e.latlng.lat).toFixed(5) +\n ' | zoom: ' + map.getZoom() + ' ');\n }\n });\n\n // remove the lnlt div when mouse leaves map\n map.on('mouseout', function (e) {\n var strip = document.querySelector('.lnlt');\n if( strip !==null) strip.remove();\n });\n\n };\n\n //$(el).keypress(67, function(e) {\n map.on('preclick', function(e) {\n if (e.originalEvent.ctrlKey) {\n if (document.querySelector('.lnlt') === null) lnlt = addElement();\n lnlt.text(\n ' lon: ' + (e.latlng.lng).toFixed(5) +\n ' | lat: ' + (e.latlng.lat).toFixed(5) +\n ' | zoom: ' + map.getZoom() + ' ');\n var txt = document.querySelector('.lnlt').textContent;\n console.log(txt);\n //txt.innerText.focus();\n //txt.select();\n setClipboardText('\"' + txt + '\"');\n }\n });\n\n }\n ).call(this.getMap(), el, x, data);\n}","data":null},{"code":"function(el, x, data) {\n return (function(el,x,data){\n var map = this;\n\n map.on('keypress', function(e) {\n console.log(e.originalEvent.code);\n var key = e.originalEvent.code;\n if (key === 'KeyE') {\n var bb = this.getBounds();\n var txt = JSON.stringify(bb);\n console.log(txt);\n\n setClipboardText('\\'' + txt + '\\'');\n }\n })\n }).call(this.getMap(), el, x, data);\n}","data":null}]}}</script>

```r
mapview::mapview(occitanie, zcol = "NAME_2", layer.name = "Département") +
mapview::mapview(prefectures, col.regions = "#3f3f3f", cex = "population")
```

---

## Quelques bonnes ressources

- [**Spatial Data Science with R**](https://rspatial.org/raster/index.html)

- [**sf vignettes**](https://r-spatial.github.io/sf/articles/sf1.html)

- [**Introduction to Geospatial Raster and Vector Data with R**](https://datacarpentry.org/r-raster-vector-geospatial/)

- [**StatnMap**](https://statnmap.com/fr/categories/spatial/)

- [**Raster analysis in R**](https://mgimond.github.io/megug2017/)

- [**Geocomputation with R**](https://geocompr.robinlovelace.net/)

- [**gglot2 et sf**](https://www.r-spatial.org/r/2018/10/25/ggplot2-sf.html)

- [**Très bonne introduction d'Olivier**](https://oliviergimenez.github.io/introspatialR/#1)