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

# Les opérations spatiales sur les données


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}
library(sf)
library(tidyverse)
library(mapview)
library(ggplot2)
```

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}
load("data/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 = c("NOM_DEP", "NOM_EPCI"), legend = F)
```

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 (`, ,`), on peut rajouter, comme dans une opération `[` classique de R les colonnes que l'on souhaite garder.

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.

```{block2, type='rmdnote'}
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), layer.name = c("Loire-Atlantique avec un buffer de 5 km", "Loire-Atlantique"), zcol = c("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 = c("NOM_DEP", "NOM_EPCI"), legend = F)
```

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'Ile 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 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 = c("NOM_DEP", "NOM_EPCI"), legend = F)
```

## 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 test.
Nous allons utiliser un polygone (a), des lignes (l) et des points (p).

```{r}
# polygone (a)
a_poly <- st_polygon(list(rbind(c(-1, -1), c(1, -1), c(1, 1), c(-1, -1))))
a <- st_sfc(a_poly)
# ligne (l)
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)

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

```{r,echo = F}
ggplot() + geom_sf(data = a, color = "red") + geom_sf(data = p) + geom_sf(data = l, color = "green") + theme_minimal()
```

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 vont nous permettre de répondre à ces questions. `sf` contient une liste de fonctions qui permettent chacune de répondre à 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 travail 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=F`.

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

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

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

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 = F)
```

`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 = F)
```

Equivalent à :

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

`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 = F)
```

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()`

### Exercices

- Créer un objet des points de p qui intersectent avec le polygone a


## Les jointures spatiales

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

Les jointures spatiales se basent sur un appariement sur un espace geographique 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 polygone 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 sur ce département et regardons dans quelle partie du territoire elles se trouvent.

```{r}
load("data/sirene.RData")
sirene44_sel <- sirene44 %>%
  filter(APET700 == "0893Z")

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

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 = c("NOM_DEP", "NOM_EPCI", "NOM_EPCI"), legend = F)
```

```{block2, type='rmdnote'}
Une jointure entre deux couches de données géographique 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 dans lesquels il se retrouve, 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))

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=T` 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 = T)

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

### Exercice

Le but de cet exercice va être d'exploiter les données *DVF* sur les transactions immobilières dans l'ancien et la carte des quartiers de Nantes pour obtenir un prix moyen des transactions par quartier.
On va utiliser pour DVF l'API mise en place par Christian Quest.

- Données DVF : http://api.cquest.org/dvf 

- Contour des quartiers de Nantes : https://data.nantesmetropole.fr/explore/dataset/244400404_quartiers-nantes/information/?disjunctive.nom

On veut produire les infos suivantes par quartier et année :

- Volume de ventes
- Pourcentage de maisons dans les ventes
- Prix moyen au m2 par type de bien

## Les calculs de distance

### Matrice de distances

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

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

```{r,warning = F}
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 = T`. Dans ce cas le résultat est un vecteur.

### Identification du plus proche voisin

Un besoin fréquent en traitement 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 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
)
```

`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
  )
```