Chapitre 6 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.

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)⁶ de la France Métropolitaine

load("data/admin_express.RData")

glimpse(epci_geo)

Rows: 1,232
Columns: 5
$ ID        <chr> "EPCI00000000000000000001", "EPCI00000000000000000002", "EPCI00000000000000000003", "EPC...
$ CODE_EPCI <chr> "200000172", "200000438", "200000545", "200000628", "200000800", "200000925", "200000933...
$ NOM_EPCI  <chr> "CC Faucigny-Glières", "CC du Pays de Pontchâteau Saint-Gildas-des-Bois", "CC des Portes...
$ TYPE_EPCI <chr> "CC", "CC", "CC", "CC", "CC", "CC", "CC", "CC", "CC", "CC", "CA", "CC", "CC", "CC", "CC"...
$ geometry  <MULTIPOLYGON [m]> MULTIPOLYGON (((964684.4 65..., MULTIPOLYGON (((316650.4 67..., MULTIPOLYGO...

glimpse(departements_geo)

Rows: 96
Columns: 8
$ ID        <chr> "DEP000000000000000000001", "DEP000000000000000000002", "DEP000000000000000000003", "DEP...
$ NOM_DEP   <chr> "Ain", "Aisne", "Allier", "Alpes-de-Haute-Provence", "Hautes-Alpes", "Alpes-Maritimes", ...
$ NOM_DEP_M <chr> "AIN", "AISNE", "ALLIER", "ALPES DE HAUTE PROVENCE", "HAUTES ALPES", "ALPES MARITIMES", ...
$ INSEE_DEP <chr> "01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12", "13", "14", "15"...
$ INSEE_REG <chr> "84", "32", "84", "93", "93", "93", "84", "44", "76", "44", "76", "76", "93", "28", "84"...
$ CHF_DEP   <chr> "01053", "02408", "03190", "04070", "05061", "06088", "07186", "08105", "09122", "10387"...
$ geometry  <MULTIPOLYGON [m]> MULTIPOLYGON (((838383.2 65..., MULTIPOLYGON (((708719 6956..., MULTIPOLYGO...
$ AREA      [m^2] 5774273633 [m^2], 7418609601 [m^2], 7365633524 [m^2], 6994210238 [m^2], 5685029372 [m^2]...

glimpse(regions_geo)

Rows: 13
Columns: 6
$ ID        <chr> "REG000000000000000000001", "REG000000000000000000002", "REG000000000000000000003", "REG...
$ NOM_REG   <chr> "Auvergne-Rhône-Alpes", "Bourgogne-Franche-Comté", "Bretagne", "Centre-Val de Loire", "C...
$ NOM_REG_M <chr> "AUVERGNE RHONE ALPES", "BOURGOGNE FRANCHE COMTE", "BRETAGNE", "CENTRE VAL DE LOIRE", "C...
$ INSEE_REG <chr> "84", "27", "53", "24", "94", "44", "32", "28", "75", "76", "52", "93", "11"
$ CHF_REG   <chr> "69123", "21231", "35238", "45234", "2A004", "67482", "59350", "76540", "33063", "31555"...
$ geometry  <MULTIPOLYGON [m]> MULTIPOLYGON (((985059.5 65..., MULTIPOLYGON (((880572.7 67..., MULTIPOLYGO...

6.1 Filtrer

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

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.

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.

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

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.

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)

6.2 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).

# 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")

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.

st_intersects(p, a)

Sparse geometry binary predicate list of length 4, where the predicate was `intersects'
 1: 1
 2: 1
 3: (empty)
 4: (empty)

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

st_disjoint(p, a)

Sparse geometry binary predicate list of length 4, where the predicate was `disjoint'
 1: (empty)
 2: (empty)
 3: 1
 4: 1

Le résultat de cette opération est une liste. Par défaut, la fonction st_intersect() renvoie une 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.

st_intersects(p, a, sparse = F)

      [,1]
[1,]  TRUE
[2,]  TRUE
[3,] FALSE
[4,] FALSE

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

st_within(p, a, sparse = F)

      [,1]
[1,]  TRUE
[2,] FALSE
[3,] FALSE
[4,] 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().

st_is_within_distance(p, a, dist = 0.8)

Sparse geometry binary predicate list of length 4, where the predicate was `is_within_distance'
 1: 1
 2: 1
 3: (empty)
 4: 1

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

st_touches(p, a, sparse = F)

      [,1]
[1,] FALSE
[2,]  TRUE
[3,] FALSE
[4,] 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.

st_contains(a, l, sparse = F)

      [,1] [,2]
[1,] FALSE TRUE

Equivalent à :

st_within(l, a, sparse = F)

      [,1]
[1,] FALSE
[2,]  TRUE

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.

st_crosses(l, a, sparse = F)

      [,1]
[1,]  TRUE
[2,] 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()

6.2.1 Exercices

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

6.3 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.

6.3.1 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⁸. Prenons les entreprises de production de sel sur ce département et regardons dans quelle partie du territoire elles se trouvent.

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.

sirene44_sel_avec_code_epci <- sirene44_sel %>%
  st_join(epci_geo)

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

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

6.3.2 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é ?

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)

# A tibble: 17 x 3
   NOM_EPCI                                    n                                                       geometry
   <chr>                                   <int>                                             <MULTIPOLYGON [m]>
 1 CA Redon Agglomération                      3 (((306623.5 6741727, 306242.5 6741883, 305602.7 6742118, 3055~
 2 CA de la Presqu'île de Guérande Atlant~     2 (((276740.4 6716295, 276731.6 6716306, 276724.3 6716291, 2767~
 3 CC du Pays d'Ancenis                        2 (((368752 6717819, 368723.8 6717821, 368302 6717946, 368143.4~
 4 CA Clisson Sèvre et Maine Agglo             1 (((377846 6675063, 377875.1 6675064, 377905.6 6675050, 377960~
 5 CA de la Région Nazairienne et de l'Es~     1 (((296085.7 6697957, 296023.1 6697783, 296007.6 6697761, 2957~
 6 CA Pornic Agglo Pays de Retz                1 (((304025.3 6684737, 304057.2 6684904, 304108.9 6684991, 3041~
 7 CC Châteaubriant-Derval                     1 (((352473.9 6745979, 352463.6 6746000, 352442.9 6746125, 3524~
 8 CC d'Erdre et Gesvres                       1 (((344439 6712389, 344500.3 6712440, 344533.6 6712504, 344640~
 9 CC de Grand Lieu                            1 (((350701.8 6678585, 350710.7 6678630, 350792.7 6678700, 3507~
10 CC de la Région de Blain                    1 (((347277.1 6718325, 347241.9 6718337, 347188.6 6718338, 3471~
11 CC de Nozay                                 1 (((348977.2 6729779, 349279.8 6730100, 349309.8 6730099, 3493~
12 CC du Pays de Pontchâteau Saint-Gildas~     1 (((316650.4 6727168, 316716 6727700, 316721.5 6727722, 316814~
13 CC du Sud-Estuaire                          1 (((312929 6700671, 313038.5 6700721, 313134.2 6700753, 313313~
14 CC Estuaire et Sillon                       1 (((321954.7 6705420, 321910 6705825, 321141.1 6706100, 320580~
15 CC Sèvre et Loire                           1 (((369365.5 6686430, 369203.9 6686525, 368865.3 6686806, 3687~
16 CC Sud Retz Atlantique                      1 (((338804.8 6658661, 338783.6 6658652, 338773.2 6658622, 3387~
17 Nantes Métropole                            1 (((343448 6696675, 343551 6697017, 343478.2 6697111, 343478 6~

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.

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)

6.3.3 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

6.4 Les calculs de distance

6.4.1 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().

centres_departements_pdl <- st_centroid(departements_geo) %>%
  filter(INSEE_REG == "52")

st_distance(centres_departements_pdl)

Units: [m]
          [,1]     [,2]      [,3]      [,4]      [,5]
[1,]      0.00 84534.14 115940.12 159007.12  81693.02
[2,]  84534.14     0.00  84436.05  89197.31  97167.70
[3,] 115940.12 84436.05      0.00  67721.05 170454.77
[4,] 159007.12 89197.31  67721.05      0.00 186162.73
[5,]  81693.02 97167.70 170454.77 186162.73      0.00

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.

6.4.2 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.

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.

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.

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
  )