GeoViz2018_R_espace.Rmd

---
title: "Espace et Temps avec R"
subtitle: '1 - Manipulation de données spatiales'
author: "Robin Cura, Hélène Mathian & Lise Vaudor"
date: ["16/10/2018<br/><br /> École Thématique GeoViz 2018 "]
output:
  xaringan::moon_reader:
    chakra: "www/libs/remark.min.js"
    lib_dir: "www/libs"
    css: ["www/xaringan-themer.css", "www/style.css"]
    nature:
        beforeInit: "www/macros.js"
        highlightLines: true
editor_options: 
  chunk_output_type: console
---
```{r xaringan-themer, include = FALSE}
library(tidyverse)
library(sf)
library(xaringan)
library(xaringanthemer) # devtools::install_github("gadenbuie/xaringanthemer")
duo(
  primary_color = "#1F4257",
  secondary_color = "#F97B64",
  header_font_google = google_font("Josefin Sans"),
  text_font_google   = google_font("Montserrat", "300", "300i"),
  code_font_google   = google_font("Droid Mono"),
  outfile = "www/xaringan-themer.css"
)
library(knitr)
opts_chunk$set(cache = TRUE,
               cache.lazy = FALSE)

hook_source <- knitr::knit_hooks$get('source')
knitr::knit_hooks$set(source = function(x, options) {
  x <- stringr::str_replace(x, "^[[:blank:]]?([^*].+?)[[:blank:]]*#<<[[:blank:]]*$", "*\\1")
  hook_source(x, options)
})
```

## Sommaire
.pull-left[
Charger des données spatiales
- [Les Simples Features et R](#objets-sf)
- [Lecture de fichiers spatiaux](#lecture-sf)
- [Conversion de tableaux en sf](#conversion-spatiale)

Manipuler des données spatiales
- [Éditer des données](#modifier-sf)
- [Systèmes de projections](#projections-sf)
- [Agrégations spatiales](#agregations-sf)
- [Jointures spatiales](#jointures-spatiales)
- [Opérations géométriques](#operations-spatiales)


- [Quelques exemples d'analyse spatiale avec R](#analyse-spatiale)
]

.pull-right[
Exploration visuelle de données spatiales
- [Afficher un objet avec plot](#plot-sf)
- [Exploration rapide avec mapview](#plot-mapview)
- [Exploration thématique avec mapview](#carto-mapview)

Cartographier avec R
- [Cartographier avec ggplot2](#carto-ggplot2)
- [Cartographie statique avec cartography](#carto-cartography)
- [Cartographie dynamique avec leaflet](#carto-leaflet)

]


---
## Les données spatiales dans R



Deux formats de données spatiales coexistent :
  - `sp`, le format le plus ancien et répandu
  - `sf`, un format plus récent, puissant, mais pas encore universel
  
.pull-left[
.smaller[
```{r, message = FALSE}
library(rgdal)
library(sp)
iris_Paris_sp <- readOGR("data/ParisIris2017.shp",
                         stringsAsFactors = FALSE)
```
]

![](www/img/spatial_str_sp.png)
]

.pull-right[
.smaller[
```{r}
library(sf)
iris_Paris_sf <- st_read("data/ParisIris2017.shp",
                         stringsAsFactors = FALSE)
str(iris_Paris_sf)
```
]
]

---
## Les données spatiales dans R

.pull-left[
.smaller[
```{r}
plot(iris_Paris_sp)
```
]
]


.pull-right[
.smaller[
```{r}
plot(iris_Paris_sf)
```
]
]

---
## Les données spatiales dans R

On va se contenter de mobiliser le *package* `sf`, qui permet une manipulation d'objets simples, semblables à des `data.frame`/`tibbles`.

**N.B.** : Pour une très grande majorité d'opérations, `sf` suffit. Il peut toutefois être nécessaire de passer de `sf` à `sp` dès lors qu'on souhaite utiliser des fonctions spatiales avancées, par exemple pour mobiliser des méthodes d'analyse spatiale.


---
name: objets-sf
## `sf` : Les *Simples Features*


.small[
.center[
![:scale 60%](www/img/sf-classes.png)

R. Lovelace, J. Nowosad & J. Muenchow (2018), Geocomputation with R - https://geocompr.robinlovelace.net/
]
]

---
## `sf` : Structure d'un objet `sf`

.small[
```{r}
iris_Paris_sf
```
]

- $\rightarrow$ Un `data.frame`, doté d'une colonne de `geometry` et d'un **système de coordonnées/projection**.

```{r, echo = FALSE}
rm(iris_Paris_sf, iris_Paris_sp)
```

---
name: lecture-sf
## Lecture de fichiers géographiques

#### Toutes les fonctions de `sf` s'inspirent de la syntaxe des fonctions `PostGIS` : `st_OPERATION` :

- Lecture d'un shapefile
.small[
```{r}
irisParis <- st_read(dsn = "data/ParisIris2017.shp", 
                       stringsAsFactors = FALSE) # Comme pour la lecture des data.frame
```
]
---
name: projections-sf
## Projections / transformations

#### Comme dans tout SIG, le système de projection (`st_crs`) est récuperé à la lecture d'un objet :

.small[
```{r}
st_crs(irisParis)
```
]

Ici, le "code de projection" est bien reconnu, mais le code EPSG n'est pas renseigné. On va le spécifier pour clarifier (même si ce n'est pas utile) :

.small[
```{r}
irisParis <- irisParis %>%
  st_set_crs(2154) # SRID/EPSG du Lambert 93
head(irisParis)
```
]

---
## Projections / transformations

#### On peut re-projeter un objet `sf` avec la fonction `st_transform(SRID/CRS)` :

.pull-left[
.small[
```{r}
irisParis %>%
  select(geometry) %>% # Lambert 93
  plot()
```
]
]

.pull-right[
.small[
```{r}
GallPeters <- "+proj=cea +lon_0=0 +lat_ts=45 \
  +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"
irisParis %>%
  select(geometry) %>%
  st_transform(crs = GallPeters) %>% #<<
  plot()
```
]
]

---
name: plot-sf
## Visualisation exploratoire : `plot()`

On peut utiliser la fonction de base `plot()`, qui affiche alors l'ensemble des attributs du jeu de données :

.small[
```{r, fig.width = 12, fig.height = 5}
plot(irisParis)
```
]

---
## Visualisation exploratoire : `plot()`

Pour "cartographier" le contenu d'une seule variable :

.small[
```{r, fig.width = 12, fig.height = 5}
plot(irisParis["NOM_C"])
```
]

---
name: plot-mapview
## Visualisation exploratoire : `mapview`

Le *package* `mapview` permet de mener rapidement une exploration des données spatiales et attributaires :

.small[
```{r, message = FALSE, fig.width=10, fig.height = 5}
library(mapview) # install.packages("mapview")
mapview(irisParis)
```
]

---
## Visualisation exploratoire : `mapview`

Le *package* `mapview` permet de mener rapidement une exploration des données spatiales et attributaires.

On peut aussi spécifier un attribut à observer avec l'instruction `zcol` :

.small[
```{r, message = FALSE, fig.width=10, fig.height = 5}
mapview(irisParis, zcol = "NOM_C")
```
]

---
name: modifier-sf
## Opérations attributaires

Comme l'objet `sf` est fondamentalement un `data.frame`, on peut lui appliquer les manipulations vues hier : filtrage (`filter`), réductions de variables (`select`), renommage de variables (`rename`) et création/modification de variables (`mutate`) :

.small[
```{r, fig.width=11, fig.height = 3.5}
irisParis %>%
  select(INSEE, CODE_) %>% # La sélection de la colonne `geometry` est implicite #<<
  mutate(ARRONDISSEMENT = as.numeric(INSEE) - 100) %>% #<<
  filter(ARRONDISSEMENT >= 75001, ARRONDISSEMENT <= 75007) %>%#<<
  mapview(zcol = "ARRONDISSEMENT")
```
]

---
name: agregations-sf
## Agrégations

Comme avec un `tibble`, on peut réaliser des opérations d'agrégation.
La géométrie est alors, elle-aussi, agrégée en conséquence :

.small[
```{r}
arrondissements_Paris <- irisParis %>%
  group_by(INSEE) %>%
  summarise()

plot(arrondissements_Paris)
```
]

---
## Jointures attributaires

Comme pour tous les `data.frame`, on peut réaliser des jointures attributaires :

.smaller[
```{r}
df_dmr  <- readRDS("dans_ma_rue_clean.RDS")
joinData <-  df_dmr %>%
  group_by(CODE_POSTAL, ANNEE_DECLARATION) %>%
  summarise(NbIncidents = n()) %>%
  ungroup()

donnees_incidents <- arrondissements_Paris %>%
  mutate(CODGEO = as.character(as.numeric(INSEE) - 100)) %>%
  left_join(joinData, by = c("CODGEO" = "CODE_POSTAL")) #<<

incidents2015 <- filter(donnees_incidents, ANNEE_DECLARATION == 2015)
incidents2017 <- filter(donnees_incidents, ANNEE_DECLARATION == 2017)

commonBreaks <- joinData %>%
  filter(ANNEE_DECLARATION %in% c(2015, 2017)) %>%
  pull(NbIncidents) %>%
  quantile(probs = seq(from = 0, to = 1, by = 0.1))
```

.pull-left[
```{r, fig.height = 3.5}
plot(incidents2015[,"NbIncidents"],
     main = "Incidents 2015",
     breaks = commonBreaks)
```
]

.pull-right[
```{r, fig.height = 3.5}
plot(incidents2017[,"NbIncidents"],
     main = "Incidents 2017",
     breaks = commonBreaks)
```
]
]
---
name: conversion-spatiale
## Conversion en `sf` et jointures spatiales

Les objets `sf` sont des objets spatiaux, on peut donc aussi effectuer des jointures spatiales (entre deux objets `sf`) :

- Conversion du jeu de données "Dans ma Rue" en objet spatial avec la fonction `st_as_sf` :

.small[
```{r}
dmr_spatial <- df_dmr %>%
  filter(ANNEE_DECLARATION %in% c(2015, 2017)) %>%
  st_as_sf(coords = c("Long", "Lat"), crs = 4326) %>% #<<
  st_transform(2154) # Pour une opération géométrique, les objets doivent avoir le même SRID
```
]


---
name: jointures-spatiales
## Jointures spatiales


- Jointure spatiale pour récuperer les IRIS contenant les points :

.small[
```{r}
dmr_augmente <- dmr_spatial %>%
  st_join(irisParis %>% select(INSEE), join = st_within) #<<
```

.pull-left[
```{r blob, eval = FALSE}
dmr_augmente %>%
  group_by(INSEE, ANNEE_DECLARATION) %>%
  summarise(NbIncidents = n()) %>%
  ungroup() %>%
  head()
```
]
]
.smaller[
.pull-right[
```{r, eval = TRUE, echo = FALSE, ref.label="blob"}
```
]
]
---
## Opérations spatiales
name: operations-spatiales

`sf` permet de réaliser la quasi-totalité des opérations attendues d'un SIG, avec des fonctions dédiées reprenant le vocabulaire classique (inspiré par `PostGIS`) :

`st_agr`, `st_area`, `st_bbox`, `st_bind_cols`, `st_boundary`, `st_buffer`, `st_cast`, `st_centroid`, `st_combine`, `st_contains`, `st_contains_properly`, `st_convex_hull`, `st_covered_by`, `st_covers`, `st_crop`, `st_crosses`, `st_difference`, `st_disjoint`, `st_distance`, `st_equals`, `st_equals_exact`, `st_graticule`, `st_interpolate_aw`, `st_intersection`, `st_intersects`, `st_is_simple`, `st_is_valid`, `st_is_within_distance`, `st_jitter`, `st_join`, `st_layers`, `st_length`, `st_line_merge`, `st_line_sample`, `st_make_grid`, `st_overlaps`, `st_point_on_surface`, `st_polygonize`, `st_relate`, `st_sample`, `st_segmentize`, `st_simplify`, `st_snap`, `st_sym_difference`, `st_touches`, `st_transform`, `st_triangulate`, `st_union`, `st_viewport`, `st_voronoi`, `st_within`, `st_wrap_dateline`, `st_write`, `st_write_db`, `st_zm`, `write_sf`

Voir les exemples illustrées dans les *vignettes* du *package* `sf` : [Manipulating Simple Feature Geometries](https://cran.r-project.org/web/packages/sf/vignettes/sf3.html) et
[Manipulating Simple Features](https://cran.r-project.org/web/packages/sf/vignettes/sf4.html)


---
## Opérations spatiales : un exemple de traitement

-> On peut chercher à créer une carte des densités d'incidents déclarés par IRIS, potentiellement selon les années :

.small[
```{r}
nbIncidentsIris <- dmr_spatial %>%
  st_join(irisParis %>% select(INSEE), join = st_within) %>%
  st_set_geometry(NULL) %>% # On n'a plus besoin que ce soit un objet spatial
  group_by(INSEE, ANNEE_DECLARATION) %>%
  summarise(NbIncidents = n()) %>% #<<
  ungroup()

irisParis_Incidents <- irisParis %>%
  mutate(surface = st_area(.)) %>% #<<
  left_join(nbIncidentsIris, by = "INSEE")

densite_incidents <- irisParis_Incidents %>%
  mutate(densiteIncidents = NbIncidents / (surface / 1E6)) %>% # Pour avoir une densité / km² #<<
  select(INSEE, CODE_, ANNEE_DECLARATION, NbIncidents, densiteIncidents)

# Les densités sont exprimées avec des "unités", qui complexifient le traitement.
# On les enlève donc :

densite_incidents <- densite_incidents %>%
  mutate(densiteIncidents = units::drop_units(densiteIncidents))
```
]

---
name: carto-mapview
## Cartographie exploratoire : `mapview`

.small[
```{r, eval = FALSE}
densite2015 <- filter(densite_incidents, ANNEE_DECLARATION == 2015)
densite2017 <- filter(densite_incidents, ANNEE_DECLARATION == 2017)

map2015 <- mapview(densite2015, zcol = "densiteIncidents", legend = TRUE)
map2017 <- mapview(densite2017, zcol = "densiteIncidents", legend = TRUE)

latticeView(map2013, map2017, ncol = 2)

```
]

---
## Cartographie exploratoire : `mapview`

.small[
```{r, eval = TRUE, echo = FALSE}
densite2015 <- filter(densite_incidents, ANNEE_DECLARATION == 2015)
densite2017 <- filter(densite_incidents, ANNEE_DECLARATION == 2017)

map2015 <- mapview(densite2015, zcol = "densiteIncidents", legend = TRUE)
map2017 <- mapview(densite2017, zcol = "densiteIncidents", legend = TRUE)

latticeView(map2015, map2017, ncol = 2)

```
]


---
name: carto-ggplot2
## Cartographie statique : `ggplot2`

On peut bien sûr afficher ces éléments avec le *package* `ggplot2`, en faisant appel à la géométrie dédiée aux objets de type `sf` : `geom_sf`

.small[
```{r, fig.width = 12, fig.height = 4.5}
ggplot(densite_incidents) +
  geom_sf(data = densite_incidents, aes(fill = densiteIncidents), lwd = 0) +
  # On est obligé de remettre la couche sf en data #<<
  facet_wrap(~ANNEE_DECLARATION) +
  scale_fill_viridis_c(name = "Incidents / Km²") +
  guides(fill = guide_colourbar(title.position = "top")) +
  coord_sf() + # On s'assure que les coordonnées du graphique respectent le CRS #<< 
  theme_minimal() +
  theme(legend.position = "bottom", legend.key.width = unit(2,"cm")) 
```
]

---
## Cartographie statique : `ggplot2`

- Certains *packages* (`ggsn`, `ggmap`) proposent des foncitonnalités permettant d'améliorer le rendu cartographique (échelles graphiques, flèches d'orientations, fonds de cartes...) ;

- mais `ggplot` est avant tout un outil de visualisation **générique** :

#### - pour créer des cartographies de qualité, mieux vaut se tourner vers des *packages* dédiés : `tmap` (qu'on ne verra pas ici) et `cartography` 

---
name: carto-cartography
## Cartographie statique : `cartography`

#### Exemple d'une carte en symboles proportionnels

.small[
```{r, fig.show="hold", eval = FALSE}
library(cartography)
# On commence par afficher un fond de carte
par(mar=c(0,0,0,0))
plot(st_geometry(arrondissements_Paris),
     col = "#1F4257",
     border = "#B6C1C8",
     bg = "#708694",
     lwd = 2)

# On ajoute des symboles proportionnels pour les incidents
propSymbolsLayer(x = densite_incidents %>% filter(ANNEE_DECLARATION == 2017),
                 var = "NbIncidents",
                 col = "#F97B64",
                 inches = 0.1, 
                 legend.title.txt = "Nombre de déclarations par IRIS en 2017")

# On ajoute les labels des arrondissements
labelLayer(x = arrondissements_Paris %>% mutate(CP = as.numeric(INSEE) - 100),
           txt = "CP",
           col = "#B6C1C8",
           bg = "#1F4257",
           halo = TRUE,
           overlap = FALSE, show.lines = FALSE)

```
]

---
## Cartographie statique : `cartography`

#### Exemple d'une carte en symboles proportionnels
.center[
.small[
```{r, fig.show="hold", eval = TRUE, echo = FALSE, fig.height = 5.5}
library(cartography)
# On commence par afficher un fond de carte
par(mar=c(0,0,0,0))
plot(st_geometry(arrondissements_Paris),
     col = "#1F4257",
     border = "#B6C1C8",
     bg = "#708694",
     lwd = 2)

# On ajoute des symboles proportionnels pour les incidents
propSymbolsLayer(x = densite_incidents %>% filter(ANNEE_DECLARATION == 2017),
                 var = "NbIncidents",
                 col = "#F97B64",
                 inches = 0.1, 
                 legend.title.txt = "Nombre de déclarations par IRIS en 2017")

# On ajoute les labels des arrondissements
labelLayer(x = arrondissements_Paris %>% mutate(CP = as.numeric(INSEE) - 100),
           txt = "CP",
           col = "#B6C1C8",
           bg = "#1F4257",
           halo = TRUE,
           overlap = FALSE, show.lines = FALSE)

```
]
]
---
name: carto-leaflet
## Cartographie dynamique : `leaflet`

#### - Pour créer des cartes "dynamiques", en plus de `mapview` (qui sert surtout à l'exploration de données), on peut utiliser le *package* `leaflet`, qui génère une carte en HTML+JavaScript que l'on pourra placer sur une page web quelconque :

.small[
```{r, eval = FALSE, echo = TRUE}
library(leaflet)
densite_incidents2017 <- densite_incidents %>%
  filter(ANNEE_DECLARATION == 2017) %>%
  st_transform(4326) # Leaflet requiert une couche en WGS84

seuils <- quantile(densite_incidents2017$densiteIncidents,
                   probs = seq(from = 0, to = 1, by = 0.2))

colorPalette <- colorBin("YlOrRd",
                         domain = densite_incidents2017$densiteIncidents,
                         bins = seuils)

infosPopup <- sprintf("<strong>%s</strong><br/>%s incidents  (%.1f / km²)",
                      densite_incidents2017$CODE_,
                      densite_incidents2017$NbIncidents,
                      densite_incidents2017$densiteIncidents) %>% lapply(htmltools::HTML)

leaflet(data = densite_incidents2017) %>%
  addProviderTiles(providers$CartoDB.DarkMatterNoLabels) %>%
  addPolygons(
    fillColor = ~colorPalette(densiteIncidents), fillOpacity = 0.7,
    color = "white", weight = .5, opacity = 1, dashArray = "3",
    label = infosPopup) %>%
  addLegend(pal = colorPalette, values = ~densiteIncidents, opacity = 0.7,
            title = "Densités d'incidents<br/>[incident/km²]", position = "topright")
```
]

---
## Cartographie dynamique : `leaflet`

#### - Pour créer des cartes "dynamiques", en plus de `mapview` (qui sert surtout à l'exploration de données), on peut utiliser le *package* `leaflet`, qui génère une carte en HTML+JavaScript que l'on pourra placer sur une page web quelconque :

.small[
```{r, eval = TRUE, echo = FALSE, message=FALSE, warning=FALSE, fig.width = 10, fig.height = 5.5}
library(leaflet)
densite_incidents2017 <- densite_incidents %>%
  filter(ANNEE_DECLARATION == 2017) %>%
  st_transform(4326) # Leaflet requiert une couche en WGS84

seuils <- quantile(densite_incidents2017$densiteIncidents, probs = seq(from = 0, to = 1, by = 0.2))
colorPalette <- colorBin("YlOrRd", domain = densite_incidents2017$densiteIncidents, bins = seuils)
infosPopup <- sprintf("<strong>%s</strong><br/>%s incidents  (%.1f / km²)",
                      densite_incidents2017$CODE_,
                      densite_incidents2017$NbIncidents,
                      densite_incidents2017$densiteIncidents) %>% lapply(htmltools::HTML)

leaflet(data = densite_incidents2017) %>%
  addProviderTiles(providers$CartoDB.DarkMatterNoLabels) %>%
  addPolygons(
    fillColor = ~colorPalette(densiteIncidents), fillOpacity = 0.7,
    color = "white", weight = .5, opacity = 1, dashArray = "3",
    label = infosPopup) %>%
  addLegend(pal = colorPalette, values = ~densiteIncidents, opacity = 0.7,
            title = "Densités d'incidents<br/>[incident/km²]",
            position = "topright")
```
]