06-les-operations-sur-donnees-spatiales.Rmd

# Les opérations spatiales sur les données {#les-operations-sur-donnees-spatiales}


Les opérations spatiales sont des opérations prenant nos données en entrée pour en sortir un résultat dépendant de leur composante spatiale (forme, localisation).

Dans ce chapitre, nous allons utiliser les packages suivants.

```{r pkg_chap6_op_spat}
# CRAN
library(mapview)
library(sf)
library(tidyverse)
```

Nous utiliserons les données de la table des régions de la France métropolitaine et des établissements publics de coopération intercommunale (EPCI)^[https://fr.wikipedia.org/wiki/%C3%89tablissement_public_de_coop%C3%A9ration_intercommunale] de la France Métropolitaine


```{r, eval = FALSE}
load("extdata/admin_express.RData")
```

```{r}
glimpse(epci_geo)
```

```{r}
glimpse(departements_geo)
```

```{r}
glimpse(regions_geo)
```

## Filtrer

Nous souhaitons par exemple filtrer nos EPCI sur les EPCI du département de Loire-Atlantique.

```{r}
departement_44 <- departements_geo %>%
  filter(INSEE_DEP == "44")

epci_d44 <- epci_geo[departement_44, , op = st_within]

mapview(list(departement_44, epci_d44), zcol = list("NOM_DEP", "NOM_EPCI"), legend = FALSE)
```

L'opération de filtre sur les données spatiales fonctionne en prenant la table en entrée (`epci_geo`), la table avec laquelle on souhaite définir les lignes à garder (`departement_44`), et l'opérateur qui va définir le test entre les deux géométries. Ici cet opérateur est `st_within(x, y)`, qui renvoie `TRUE` si la géométrie de `x` est contenue à l'intérieur de celle de `y`.

On peut spécifier différents prédicats spatiaux pour réaliser ce filtre.

En deuxième argument des `[ , , ]`, on peut ajouter, comme dans une opération `[` classique de R, les colonnes que l'on souhaite garder.

On peut aussi écrire en *tidyverse* :

```{r}
# 1er méthode
epci_d44 <- epci_geo %>% 
  st_filter(departement_44, .predicate = st_within)

# 2e méthode
epci_d44 <- epci_geo %>%
  filter(st_within(., departement_44, sparse = FALSE))
```
Il faut noter la présence du paramètre `sparse = FALSE`. 
L'option `sparse` est abordée dans la partie [Prédicats spatiaux]. Il faut retenir que pour utiliser les prédicats à l'intérieur de la fonction `filter`, il faut utiliser l'option `sparse = FALSE`.

On voit ici que le résultat n'est pas très concluant :  il manque 3 EPCI du département, ceux qui sortent des frontières de celui-ci.
Prenons un buffer autour du département.

> Qu'est ce qu'un buffer ? C'est un tampon qui va nous permettre d'appliquer une transformation sur un objet vectoriel.
>
> A partir d'une couche de départ de type ponctuel, linéaire ou polygonal, le buffer va créer une nouvelle couche vectorielle. La géométrie de cette couche représente des objets surfaciques dont les frontières sont positionnées à une distance euclidienne, définie par l'utilisateur, des limites des objets vectoriels de la couche de départ.

La fonction qui permet de faire cela avec `sf` s'appelle `st_buffer()`.

`st_buffer()` prend en paramètre : 

- un objet de classe *sf* 
- une distance dont l'unité est définie par celle de l'objet `sf`, que l'on peut obtenir comme ceci `st_crs(x)$units`.

```{r}
departement_44_buffer <- departement_44 %>%
  st_buffer(dist = 5000)

mapview(list(departement_44_buffer, departement_44), legend = FALSE, 
        layer.name = c("Loire-Atlantique avec un buffer de 5 km", "Loire-Atlantique"), 
        zcol = list("NOM_DEP", "NOM_DEP"), col.regions = list("#440154FF", "#FDE725FF"))
```

```{r}
epci_d44_buffer <- epci_geo[departement_44_buffer, , op = st_within]

mapview(list(departement_44_buffer, epci_d44_buffer), zcol = list("NOM_DEP", "NOM_EPCI"), legend = FALSE)
```

On récupère 2 des 3 EPCI manquant ainsi. Celui qui manque est l'EPCI de Redon qui est à cheval sur la Loire-Atlantique, le Morbihan et l'Île et Vilaine.
Une méthode pour le récupérer est de prendre l'opérateur de filtre `st_intersect` au lieu de `st_within` en utilisant un buffer légèrement négatif de notre département pour ne pas récupérer tous les EPCI limitrophes.

```{r}
departement_44_buffer_negatif <- departement_44 %>%
  st_buffer(dist = -2000)

epci_d44 <- epci_geo[departement_44_buffer_negatif, , op = st_intersects]

mapview(list(departement_44, epci_d44), zcol = list("NOM_DEP", "NOM_EPCI"), legend = FALSE)
```

On aurait pu obtenir le même résultat en prenant les EPCI à l'intérieur du département 44 (avec `st_within`) ainsi que les EPCI qui chevauchent (avec `st_overlaps`)  le département 44.

```{r}
epci_d44 <- epci_geo %>%
  filter(st_within(., departement_44, sparse = FALSE) | st_overlaps(., departement_44, sparse = FALSE))

mapview(list(departement_44, epci_d44), zcol = c("NOM_DEP", "NOM_EPCI"), legend = FALSE)
```

## Prédicats spatiaux

Les prédicats spatiaux décrivent les relations spatiales entre objets. Pour bien les illustrer, on va utiliser quelques données de notre création.

Nous allons utiliser un polygone (a), des lignes (l) et des points (p), que nous créons à partir de leur coordonnées.

```{r}
# polygone (a - orange)
a_poly <- st_polygon(list(rbind(c(-1, -1), c(1, -1), c(1, 1), c(-1, -1))))
a <- st_sfc(a_poly)

# lignes (l - bleues)
l1 <- st_multilinestring(list(rbind(c(0.5, -1), c(-0.5, 1))))
l2 <- st_multilinestring(list(rbind(c(.9, -.9), c(.5, 0))))
l <- st_sfc(l1, l2)

# multi-points (p - noirs)
p_matrix <- matrix(c(0.5, 1, -1, 0, 0, 1, 0.5, 1), ncol = 2)
p_multi <- st_multipoint(x = p_matrix)
p_multi
p <- st_cast(st_sfc(p_multi), "POINT")
```

La fonction `st_cast` sert à passer d'un type de géométrie à l'autre, ici on transforme un multi-points en points.

```{r, echo=FALSE}
geom_pred <- ggplot() + 
  geom_sf(data = a, color = "salmon3", fill = "salmon2") + 
  geom_sf(data = p) + 
  geom_sf(data = l, color = "dodgerblue4", size = 1) + 
  geom_sf_text(data = rowid_to_column(st_sf(p)), mapping = aes(label = rowid), nudge_x = 0.07) +
  theme_minimal() + labs(x = "", y = "")
```

A partir de ces objets, on peut se poser les questions suivantes :

- Quels sont les points de `p` contenus dans le triangle `a` ?

- Quels sont les points de `p` qui ne sont pas contenus dans le triangle `a` ?

- Quels sont les points de `p` qui touchent le triangle `a` ?

- Quelles sont les lignes de `l` contenues dans `a` ?


Les prédicats spatiaux servent à répondre à ces questions et `sf` contient une liste de fonctions dédiée à l'une ou l'autre de ces questions.

`st_intersects()` permet de répondre à la première question, à savoir quels points de `p` sont dans `a`.

```{r}
st_intersects(p, a)
```

L'opposé de `st_intersects()` est `st_disjoint()` :  `st_disjoint(x,y)` renvoie `TRUE` pour les objets de `x` non reliés à `y`.

```{r}
st_disjoint(p, a)
```

Le résultat de cette opération est une liste. Par défaut, la fonction `st_intersect()` renvoie une *matrice creuse*^[https://fr.wikipedia.org/wiki/Matrice_creuse]. Cette structure permet d'économiser de la mémoire en n'enregistrant que les relations qui existent. Sur une opération de ce type, le gain est peu évident, mais quand on travaille sur des objets plus complexes, le gain est appréciable.

Si on souhaite mieux utiliser cette information, on peut vouloir privilégier la *matrice dense*, qui renvoie une matrice de booléen pour chaque relation possible.

Pour cela on peut utiliser l'option `sparse = FALSE`, 'sparse' signifiant en anglais 'clairsemé'.

```{r}
st_intersects(p, a, sparse = FALSE)
```

`st_within()` est une variante de `st_intersect()` qui ne renvoie `TRUE` que pour les points strictement *à l'intérieur* du polygone.

```{r}
st_within(p, a, sparse = FALSE)
```

Une variante de `st_within()` permet d'ajouter un critère de distance pour intégrer des points *presque* dans le polygone, `st_is_within_distance()`.

```{r}
st_is_within_distance(p, a, dist = 0.8)
```

`st_touches()` permet de récupérer les points qui *touchent* le polygone sans sans être à l'intérieur du polygone.

```{r}
st_touches(p, a, sparse = FALSE)
```

`st_contains(x,y)` est équivalent à `st_within(y,x)`. Par exemple si nous voulons savoir lesquelles de nos lignes `l` sont contenues dans `a`.

```{r}
st_contains(a, l, sparse = FALSE)
```

Equivalent à :

```{r}
st_within(l, a, sparse = FALSE)
```

`st_crosses()`  renvoie TRUE si l'intersection des deux géométries est une géométrie de dimension n-1 ou n est le maximum des dimensions des deux objets et si l'intersection est à l'intérieur des deux objets.

```{r}
st_crosses(l, a, sparse = FALSE)
```

Il existent encore d'autres prédicats qu'on ne détaillera pas ici : 

- `st_covers()`

- `st_covered_by()`

- `st_equals()` et `st_equals_exact()`

- `st_contains_properly()`

- `st_overlaps()`

### Exercice 1 - manipuler des objets sf
```{r mod7_exo1, child=charge_exo("m7", "exo1.rmd"), echo=FALSE}
```

Visualisation graphique :
```{r, message=FALSE, warning=FALSE, echo = FALSE}
res <- st_sf(p)[a, , op = st_intersects]
library(ggplot2)

# graphique
ggplot() + 
  geom_sf(data = a, color = "orange") + 
  geom_sf(data = p) + 
  geom_sf(data = l, color = "blue") +
  geom_sf(data = res, color = "cyan") +
  theme_minimal()
```


## Les jointures spatiales

Les jointures *attributaires*  se basent sur un appariement entre une liste des variables présentes dans les deux tables.

Les jointures spatiales se basent sur un appariement géographique, c'est à dire sur un espace géo commun.

### Jointure de points avec des polygones

Ce cas est relativement simple, une jointure spatiale entre une liste de points et une liste de polygones va attribuer pour chaque point le ou les polygones auquel il appartient.

On va utiliser ici le fichier sirene du département de Loire Atlantique géocodé par Christian Quest^[http://data.cquest.org/geo_sirene/].
Prenons les entreprises de production de sel (code APE '0893Z') sur ce département et regardons dans quelle partie du territoire elles se trouvent.

```{r, eval = FALSE}
load("extdata/sirene.RData")
```


```{r}
sirene44_sel <- sirene44 %>%
  filter(APET700 == "0893Z")

mapview(list(departement_44, epci_d44, sirene44_sel), zcol = list("NOM_DEP", "NOM_EPCI", "NOMEN_LONG"), legend = FALSE)
```

Nous allons réaliser une jointure spatiale pour récupérer le code sirene de l'EPCI où se trouve chaque entreprise.

```{r}
sirene44_sel_avec_code_epci <- sirene44_sel %>%
  st_join(epci_geo)
```


```{r}
mapview(list(departement_44, epci_d44, sirene44_sel_avec_code_epci), zcol = list("NOM_DEP", "NOM_EPCI", "NOM_EPCI"), legend = FALSE)
```


> Une jointure entre deux couches de données géographiques demande à ce que celles-ci partagent la même projection.


### Jointure de polygones avec des polygones

A la différence des appariements entre points et polygones, la jointure spatiales entre deux couches de polygones nécessite quelques critères complémentaires : souhaite-t-on joindre deux polygones dès qu'ils s'intersectent ? Souhaite-t-on joindre à un polygone de la première couche à celui de la deuxième avec lequel il partage le plus de surface en commun ?

Par exemple, imaginons que nous voulions joindre notre couche des EPCI avec celle des départements, souhaite-t-on que l'EPCI de Redon se retrouve apparié avec tous les départements auxquels il appartient, ou seulement le département dans lequel il est principalement situé ?

```{r}
epci_d44_avec_departement <- epci_d44 %>%
  st_join(departements_geo %>% st_buffer(dist = -1000))

names(epci_d44_avec_departement)
```


```{r}
epci_d44_avec_departement %>%
  select(NOM_EPCI, NOM_DEP) %>%
  group_by(NOM_EPCI) %>%
  tally() %>%
  arrange(-n)
```

Une jointure classique va donc rattacher 3 EPCI à plus de 1 département.

Avec l'option `largest = TRUE` la jointure va attribuer aux EPCI le département avec lequel il partage le plus de surface.
On voit ici que tout les EPCI adhérents à la Loire-Atlantique se retrouvent alors rattachés à la Loire-Atlantique.

```{r}
epci_d44_avec_departement <- epci_d44 %>%
  st_join(departements_geo %>% st_buffer(dist = -1000), largest = TRUE)

mapview(list(departement_44, epci_d44, sirene44_sel_avec_code_epci), zcol = list("NOM_DEP", "NOM_EPCI", "NOM_EPCI"), 
        legend = FALSE)
```

> `st_join()` utilise par défaut le prédicat st_intersects, dans le paramètre `join`.
> 
> Après une jointure spatiale, le dataframe résultat contient la géométrie du jeu de données de gauche (càd le 1er jeu de données passé dans les arguments de `st_join()`).


### Exercice 2 : exploitation des données DVF en API

```{r mod7_exo2, child=charge_exo("m7", "exo2.rmd"), echo=FALSE}
``` 

## Les mesures

Contrairement aux opérations précédentes qui sont binaires, les opérations de mesure sont continues.

### Calcul de surface  

Le calcul de surface se fait à l'aide de la fonction `st_area()`. 
Le résultat sera exprimé dans l'unité du crs. 

```{r st_area}
st_area(epci_d44)
```


### Matrice de distances

Les distances se calculent avec la fonction `st_distance()`.

```{r, warning=FALSE}
centres_departements_pdl <- st_centroid(departements_geo) %>%
  filter(INSEE_REG == "52")

st_distance(centres_departements_pdl)
```

Trois choses à noter sur le résultat :

- `st_distance()` retourne une *matrice*...

- ... contenant toute les distances calculables 2 à 2...

- ...et qui a un paramètre `Units` nous donnant l'unité de mesure des distances calculées.

Ici on calcule notre matrice sur un seul objet. Vous pouvez calculer des distances entre deux objets `x` et `y` de classe `sf`. 
Dans ce cas il fera le calcul des distances pour toutes les combinaisons possibles d'objets de `x` et de `y`. 
Une option de `st_distance()` vous permet de limiter le résultat aux calculs 2 à 2 : `by_element = TRUE`. 
Dans ce cas le résultat est un vecteur.

### Identification du plus proche voisin

Un besoin fréquent en géomatique est d'identifier l'objet le plus proche d'un autre.
La fonction qui permet cela est `st_nearest_feature()`.

Prenons l'ensemble des départements français, et trouvons celui de la région qui leur est le plus proche. On va utiliser les centroïdes pour alléger le calcul.

```{r}
index_dep_pdl <- st_nearest_feature(
  departements_geo,
  centres_departements_pdl
)

index_dep_pdl 
```

`st_nearest_feature()` renvoie un vecteur d'index en résultat. 

Pour visualiser cet index, vous pouvez utiliser ensuite la fonction `st_nearest_point()` qui va permettre de faire un lien entre les départements et le département ligérien le plus proche.

`st_nearest_point()` permet en effet de renvoyer pour deux géométries la ligne rejoignant les 2 points les plus proches.

```{r}
liens <- st_nearest_points(departements_geo,
  centres_departements_pdl[index_dep_pdl, ],
  pairwise = TRUE
)

ggplot() +
  geom_sf(data = departements_geo) +
  geom_sf(data = liens)
```

On peut utiliser aussi `st_nearest_feature()` comme un mode de jointure des données.

```{r}
departements_join <- st_join(departements_geo,
  centres_departements_pdl,
  join = st_nearest_feature
)

ggplot() +
  geom_sf(data = departements_join, aes(fill = NOM_DEP.y)) +
  labs(
    title = "Département ligérien le plus proche de chaque département français",
    fill = NULL
  )
```


## Carroyage

### Créer une grille régulière

La fonction `st_make_grid()` permet de créer une grille  régulière à partir d'un objet géo. 
Elle produit un objet rectangulaire de type `sfc`, il faut ensuite utiliser la fonction `st_sf()` pour transformer
cet objet `sfc` en objet `sf`. 

Lors de cette transformation nous rajoutons ici une colonne d'identifiants uniques.   

```{r grid}
grille <- st_make_grid(epci_d44, cellsize = 10000) %>% 
  st_as_sf() %>% 
  rowid_to_column("id")

mapview(grille) + mapview(epci_d44, col.regions = "white", alpha.regions = 0.5)
```

On constate que plusieurs carreaux sont inutiles. On peut réduire leur nombre avec `st_intersects()` comme vu précédemment.
```{r}
grille_ajustee <- grille[epci_d44, , op = st_intersects]
mapview(grille_ajustee) + mapview(epci_d44, col.regions = "white", alpha.regions = 0.5)
```


### Compter des points dans un polygone

On cherche ici à compter par carreau, le nombre d'entreprises de type "Boulangerie et boulangerie-pâtisserie" (code APE 1071C) du fichier sirene 44 .
```{r}
boul_44 <- filter(sirene44, APET700 == "1071C") %>% select(SIREN, NOMEN_LONG)
mapview(grille_ajustee) + mapview(epci_d44, col.regions = "white", alpha.regions = 0.5) + mapview(boul_44, col.regions = "grey")
```

La grille ajustée comprend `r nrow(grille_ajustee)` carreaux. 
Le 44 comprend `r nrow(boul_44)` boulangeries.

```{r}
intersection <- st_intersects(grille_ajustee, boul_44, sparse = TRUE)

intersection[1:10]

```

L'objet **intersection** est une liste de la longueur de l'objet **grille_ajustee** et chaque élément de la liste contient le ou les index des éléments de l'objet **boul_44** qu'il intersecte.
Exemple avec le 46e carreau :

```{r}
intersection[[46]]

```
Le 46e carreau comprend les boulangeries des lignes `r intersection[[46]] %>% paste(collapse = ", ")` de la table `boul_44`.

Pour compter le nombre de boulangeries par carreau, on peut parcourir la liste et reporter la longueur de chacun de ces éléments. 

```{r intersect3.4}
grille_boul <- mutate(grille_ajustee, nb_boulang = map_dbl(intersection, length))
select(grille_boul, nb_boulang) %>% 
  plot
```


Une autre solution est de passer par les jointures spatiales :

```{r}
nb_boul_par_carreau <- st_join(x = grille_ajustee, y = boul_44) %>% 
  group_by(id) %>% 
  summarise(nb_boulang = length(SIREN))

select(nb_boul_par_carreau, nb_boulang) %>% plot
```