library(sf)
library(tidyverse)

Référence

Ce petit tutoriel décrit l’usage basique de sf, les exemples sont extraits de ce tutoriel complet (qui n’utilise cependant pas la logique tidyverse).

Données

Les données utilisées sont issues de ce tutoriel sur la géomatique avec R et accessible sur leur répertoire Github.

communes <- st_read("https://github.com/rCarto/geomatique_avec_r/blob/main/data/lot46.gpkg?raw=true",
               layer="commune")

## Reading layer `commune' from data source 
##   `https://github.com/rCarto/geomatique_avec_r/blob/main/data/lot46.gpkg?raw=true' 
##   using driver `GPKG'
## Simple feature collection with 313 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 539668.5 ymin: 6346290 xmax: 637380.9 ymax: 6439668
## Projected CRS: RGF93 / Lambert-93

route <- st_read("https://github.com/rCarto/geomatique_avec_r/blob/main/data/lot46.gpkg?raw=true", 
                 layer = "route", quiet = TRUE)

Objets `sf`

La première étape est de caractériser ces objets spatiaux. Ces objets sont de la classe sf.

is(communes)

## [1] "sf"       "oldClass"

is(route)

## [1] "sf"       "oldClass"

Un objet sf est un tableau qui se manipulera comme un data.frame. Cependant, l’objet sf a la caractéristique de contenir une colonne geometry qui donne les coordonnées de l’objet géographique. Cet objet géographique peut être

un point;
une ligne;
un polygone;
un ensemble de polygones (MULTIPOLYGON);
etc… (on trouvera le descriptif des formats ici)

communes

## Simple feature collection with 313 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 539668.5 ymin: 6346290 xmax: 637380.9 ymax: 6439668
## Projected CRS: RGF93 / Lambert-93
## First 10 features:
##    INSEE_COM         NOM_COM         STATUT POPULATION     AGR_H    AGR_F
## 1      46001           Albas Commune simple        522  4.978581 0.000000
## 2      46002          Albiac Commune simple         67  0.000000 9.589041
## 3      46003        Alvignac Commune simple        706 10.419682 0.000000
## 4      46004         Anglars Commune simple        219  0.000000 0.000000
## 5      46005 Anglars-Juillac Commune simple        329  4.894895 4.894895
## 6      46006   Anglars-Nozac Commune simple        377  4.840849 0.000000
## 7      46007        Arcambal Commune simple        988 10.126317 0.000000
## 8      46008      Les Arques Commune simple        203  0.000000 0.000000
## 9      46009          Assier Commune simple        642 10.194444 0.000000
## 10     46010          Aujols Commune simple        367  0.000000 0.000000
##        IND_H     IND_F     BTP_H    BTP_F     TER_H      TER_F
## 1   4.936153  0.000000  9.957527 0.000000 44.917145  34.681799
## 2   0.000000  0.000000  4.794521 0.000000  4.794521   9.589041
## 3  10.419682  5.209841 10.419682 0.000000 57.308249  78.147612
## 4  20.000000 15.000000 10.000000 0.000000 20.000000  20.000000
## 5   4.894895  0.000000  0.000000 0.000000 29.369369  29.369369
## 6   0.000000  0.000000  9.681698 4.840849 43.567639  38.726790
## 7  60.757902 10.126317 40.505268 5.063158 55.694743 121.515804
## 8   4.924242  0.000000  0.000000 0.000000  9.848485   9.848485
## 9  49.006029 14.996466 14.649874 0.000000 19.692951  68.227761
## 10  9.837838  0.000000 24.594595 0.000000 24.594595  63.945946
##                              geom
## 1  MULTIPOLYGON (((559262 6371...
## 2  MULTIPOLYGON (((605540.7 64...
## 3  MULTIPOLYGON (((593707.7 64...
## 4  MULTIPOLYGON (((613211.3 64...
## 5  MULTIPOLYGON (((556744.9 63...
## 6  MULTIPOLYGON (((576667.2 64...
## 7  MULTIPOLYGON (((581404 6370...
## 8  MULTIPOLYGON (((558216 6389...
## 9  MULTIPOLYGON (((612729.6 63...
## 10 MULTIPOLYGON (((581404 6370...

route

## Simple feature collection with 16096 features and 17 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 539670.7 ymin: 6346010 xmax: 637504.4 ymax: 6440193
## Projected CRS: RGF93 / Lambert-93
## First 10 features:
##                          ID           VOCATION NB_CHAUSSE NB_VOIES   ETAT ACCES
## 1  BDCTRORO0000000021289990     Liaison locale 1 chaussée   1 voie Revêtu Libre
## 2  BDCTRORO0000000021290017     Liaison locale 1 chaussée   1 voie Revêtu Libre
## 3  BDCTRORO0000000021326471     Liaison locale 1 chaussée   1 voie Revêtu Libre
## 4  BDCTRORO0000000021328689     Liaison locale 1 chaussée   1 voie Revêtu Libre
## 5  BDCTRORO0000000202866717  Liaison régionale 1 chaussée   1 voie Revêtu Libre
## 6  BDCTRORO0000000021272608     Liaison locale 1 chaussée   1 voie Revêtu Libre
## 7  BDCTRORO0000000021272607 Liaison principale 1 chaussée  2 voies Revêtu Libre
## 8  BDCTRORO0000000021337059  Liaison régionale 1 chaussée   1 voie Revêtu Libre
## 9  BDCTRORO0000000021346024     Liaison locale 1 chaussée   1 voie Revêtu Libre
## 10 BDCTRORO0000000021328392     Liaison locale 1 chaussée   1 voie Revêtu Libre
##    POS_SOL RES_VERT        SENS NB_VOIES_M NB_VOIES_D TOPONYME
## 1   Au sol      Non Double sens Sans objet Sans objet     <NA>
## 2   Au sol      Non Double sens Sans objet Sans objet     <NA>
## 3   Au sol      Non Double sens Sans objet Sans objet     <NA>
## 4   Au sol      Non Double sens Sans objet Sans objet     <NA>
## 5   Au sol      Non Double sens Sans objet Sans objet     <NA>
## 6   Au sol      Non Sens direct Sans objet Sans objet     <NA>
## 7   Au sol      Oui Double sens Sans objet Sans objet     <NA>
## 8   Au sol      Non Double sens Sans objet Sans objet     <NA>
## 9   Au sol      Non Double sens Sans objet Sans objet     <NA>
## 10  Au sol      Non Double sens Sans objet Sans objet     <NA>
##                        USAGE       DATE NUM_ROUTE      CLASS_ADM GEST_ROUTE
## 1  Logique et cartographique 2021-05-07      D258 Départementale         46
## 2  Logique et cartographique 2021-05-07      D205 Départementale         46
## 3  Logique et cartographique 2021-05-07      D158 Départementale         46
## 4  Logique et cartographique 2021-05-07      D158 Départementale         46
## 5  Logique et cartographique 2021-05-07       D13 Départementale         46
## 6  Logique et cartographique 2021-05-07      D105 Départementale         46
## 7  Logique et cartographique 2021-05-07      D820 Départementale         46
## 8  Logique et cartographique 2021-05-07       D19 Départementale         46
## 9  Logique et cartographique 2021-05-07       D79 Départementale         46
## 10 Logique et cartographique 2021-05-07       D18 Départementale         46
##                              geom
## 1  LINESTRING (622828.5 640887...
## 2  LINESTRING (614776.4 640763...
## 3  LINESTRING (542094.5 638394...
## 4  LINESTRING (541643.5 638172...
## 5  LINESTRING (578592.6 638806...
## 6  LINESTRING (577875.7 641786...
## 7  LINESTRING (577923 6417959,...
## 8  LINESTRING (606816.7 637497...
## 9  LINESTRING (608857.1 636793...
## 10 LINESTRING (616840.6 638151...

Chaque champ de cette colonne geometry est un ensemble de valeurs numériques exprimées dans un système de référence, le coordinate reference system (CRS). Pour trouver le CRS adapté à sa zone, on pourra consulter cette page.

Pour la France, le CRS 2154 est une bonne (au sens du respect des distances après projection) approximation à l’échelle du territoire métropolitaine.

Mise au format `sf`

Les objets précédents étant nativement des “objets géographiques”, on les chargeait directement avec st_read (formats .shp, .gpkg, possibles). Cependant, on peut transformer un tableau ayant deux colonnes de coordonnées grâce à la fonction st_as_sf.

# Tableau classique
population_brute <- read.csv("https://github.com/rCarto/geomatique_avec_r/raw/main/data/pop.csv")
head(population_brute) # Les colonnes x et y sont déjà en CRS(2154)

##         x       y  Ind
## 1 3630500 2383500 23.0
## 2 3629500 2384500 25.0
## 3 3630500 2384500 18.0
## 4 3631500 2384500 11.5
## 5 3644500 2384500  5.0
## 6 3645500 2384500  9.0

population <- st_as_sf(population_brute, 
                       coords = c("x", "y"), # Spécification des coordonnées
                       crs = st_crs(2154)) # On utilise st_crs pour specifier un CRS
population

## Simple feature collection with 4595 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 3604500 ymin: 2383500 xmax: 3702500 ymax: 2475500
## Projected CRS: RGF93 / Lambert-93
## First 10 features:
##     Ind                geometry
## 1  23.0 POINT (3630500 2383500)
## 2  25.0 POINT (3629500 2384500)
## 3  18.0 POINT (3630500 2384500)
## 4  11.5 POINT (3631500 2384500)
## 5   5.0 POINT (3644500 2384500)
## 6   9.0 POINT (3645500 2384500)
## 7  16.0 POINT (3646500 2384500)
## 8  67.0 POINT (3648500 2384500)
## 9   8.0 POINT (3629500 2385500)
## 10 11.0 POINT (3630500 2385500)

Manipulation et visualisation

La manipulation et la visualisation sont calquées sur celles des data.frame. On pourra utiliser tous les verbes de dplyr.

Par exemple, si on veut obtenir la commune ayant la plus grande population

filter(communes, POPULATION == max(POPULATION)) %>% 
  select(NOM_COM, POPULATION) # On note que la colonne geometry est conservé!

## Simple feature collection with 1 feature and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 570470.4 ymin: 6368056 xmax: 581068.4 ymax: 6380468
## Projected CRS: RGF93 / Lambert-93
##   NOM_COM POPULATION                           geom
## 1  Cahors      19907 MULTIPOLYGON (((580994.5 63...

Pour la représentation graphique, on pourra utiliser le plot natif ou, si l’on préfère, se servir de ggplot2 et de la fonction dédiée geom_sf.

plot(communes)

# Graphique basique des polygones
ggplot(communes) + 
  geom_sf()

# Ajout d'un attribut du tableau
ggplot(communes) +
  geom_sf(aes(fill = POPULATION)) +
  scale_fill_viridis_c()

# Ajout d'un 2e objet sf
ggplot(communes) +
  geom_sf(fill = "lightblue") +
  geom_sf(data = route, color = "red")

Manipulations spatiales

sf propose une série de fonctions (ou verbes commençant par le préfixe st_). Par exemple, st_bbox extrait le rectangle entourant un objet géométrique.

`st_intersects`

Ici, si on veut isoler les routes qui traversent la commune de Gramat et la représenter. On utilise ici st_intersects qui établit si un un ensemble de geometry (dans l’argument x) intersecte un ou plusieurs geometry.

gramat <- filter(communes, NOM_COM == "Gramat") # On extrait la commune Gramat
limites_gramat <- st_bbox(gramat) # On en extrait ses limites
# On s'intéresse aux routes qui traversent Gramat, où ici, st_intersects renvoie 
# un vecteur de booleens (c'est le role de sparse = FALSE) indiquant pour chaque
# route si elle traverse ou non Gramat.
routes_gramat <- filter(route, st_intersects(route, gramat, sparse = FALSE))

# On peut faire la représentation graphique.
ggplot(communes) +
  geom_sf(color = "black", size = 2)  +
  geom_sf(data = route, aes(color = "Hors Gramat")) +
  geom_sf(data = routes_gramat, aes(color = "Traversant Gramat")) +
  labs(color = "", title = "Commune de Gramat") +
  coord_sf(xlim = limites_gramat[c("xmin", "xmax")], 
           ylim = limites_gramat[c("ymin", "ymax")])

`st_within`

Si on ne souhaite que les routes à l’intérieur de la commune, on utilise st_within

routes_dans_gramat <- filter(route, st_within(route, gramat, sparse = FALSE))

# On peut faire la représentation graphique.
ggplot(communes) +
  geom_sf(color = "black", size = 2)  +
  geom_sf(data = route, aes(color = "Hors Gramat")) +
  geom_sf(data = routes_dans_gramat, aes(color = "Traversant Gramat")) +
  labs(color = "", title = "Commune de Gramat") +
  coord_sf(xlim = limites_gramat[c("xmin", "xmax")], 
           ylim = limites_gramat[c("ymin", "ymax")])

Unions

restaurants_complet <- st_read("https://github.com/rCarto/geomatique_avec_r/blob/main/data/lot46.gpkg?raw=true", 
                 layer = "restaurant", quiet = TRUE)

On a importé la liste des restaurants de la région. On voit ci-dessous que nombre d’entre eux sont en dehors du Lot.

ggplot(communes) +
  geom_sf() +
  geom_sf(data = restaurants_complet, color = "red") +
  labs(title = "Positions des restaurants")

Ainsi, on souhaite ne conserver que ceux sont qui sont dans notre ensemble de communes. On utilisera à nouveau st_intersects, mais ici on ne veut pas que le test se fasse pour chaque commune! Ainsi, on regroupera toutes les communes en un seul polygone avec st_union.

restaurant <- filter(restaurants_complet,
                     st_intersects(restaurants_complet,
                                   st_union(communes),
                                   sparse = FALSE))

ggplot(communes) +
  geom_sf() +
  geom_sf(data = restaurant, color = "red") +
  labs(title = "Positions des restaurants du Lot")

Barycentres et buffers

On peut obtenir le barycentre d’un polygon et/ou faire un buffer autour d’un objet géométrique.

buffer_gramat <- st_centroid(gramat) %>% # On extrait le centre de gravité 
  st_buffer(dist = 15000) # On fait le disque centré en ce point, de rayon 15km 
# Représentation
ggplot(buffer_gramat) +
  geom_sf(data = communes) +
  geom_sf(alpha = .3, fill = "red")

Exercices

Afin de s’entraîner, nous avons répondu aux exercices proposés dans le tutoriel Géomatique avec R cité plus haut.

Quelles communes ont plus de 10 restaurants et moins de 1000 habitants ?

# Villages de France
VdF <- communes %>% # Dans communes
  filter(POPULATION < 1000) %>% # On extrait celles ayant < de 1000 habitants
  mutate(NbResto = st_intersects(., restaurant) %>%
           sapply(length)) %>% 
  filter(NbResto >= 10)
pull(VdF, NOM_COM)

## [1] "Rocamadour"         "Saint-Cirq-Lapopie"

Les représenter sur une carte

ggplot(communes) +
  geom_sf(fill = "gray") +
  geom_sf(data = VdF, fill = "red") +
  geom_sf_label(data = VdF, aes(label = NOM_COM), nudge_x = 20000)

Intersecter les routes par les communes Découper les routes par les communes. Chaque segment de route récupère l’identifiant de la commune dans laquelle il se trouve. Rajouter une colonne longueur qui contient la longueur de chacun des segments.

routes_communes <- st_intersection(route, communes) %>% 
  mutate(ID == paste(NUM_ROUTE, NOM_COM)) %>% 
  select(ID, NUM_ROUTE, NOM_COM) %>% 
  mutate(longueur = st_length(.))

Calculer la longueur totale de route par commune;
Joindre ses longueurs à la table des communes;
Calculer les superficies des communes;
Convertir
Calculer les 2 indicateurs suivants:
- le rapport entre la longueur des routes et la population des communes.
- le rapport entre la longueur des routes et la superficie des communes.

indicateur_communes <- routes_communes %>% 
  group_by(NOM_COM) %>% 
  summarise(longueur_totale_route = sum(longueur) * 10^(-3)) %>% # Mise en km 
  st_drop_geometry() %>%  # Ici, on veut faire une jointure, donc on se débarasse de la geometry
  inner_join(x = communes) %>% # On met le sf en 1er argument pour que la sortie soit un sf
  mutate(surface = st_area(.) * 10^(-6),# Création des nouvelles
         longueur_par_habitant = longueur_totale_route / POPULATION, 
         longueur_par_km2 = longueur_totale_route / surface) %>% 
  select(NOM_COM, longueur_totale_route, longueur_par_habitant, longueur_par_km2)

Représenter les deux indicateurs

p1 <- ggplot(indicateur_communes) + 
  geom_sf(aes(fill = as.numeric(longueur_par_habitant))) +
  scale_fill_viridis_c() +
  geom_sf(data = route, color = "red", size = .5, alpha = .2) +
  labs(title = "Longueur de route par habitant") +
  theme(legend.position = "none")
p2 <- ggplot(indicateur_communes) + 
  geom_sf(aes(fill = as.numeric(longueur_par_km2))) +
  scale_fill_viridis_c() +
  geom_sf(data = route, color = "red", size = .5, alpha = .2) +
  labs(title = "Longueur de route par km2") +
  theme(legend.position = "none")
cowplot::plot_grid(p1, p2)

SIG avec SF

Pierre Barbillon, Anna Bonnet, Pierre Gloaguen

2022-08-23

Référence

Données

Objets `sf`

Mise au format `sf`

Manipulation et visualisation

Manipulations spatiales

`st_intersects`

`st_within`

Unions

Barycentres et buffers

Exercices

SIG avec SF

Pierre Barbillon, Anna Bonnet, Pierre Gloaguen

2022-08-23

Référence

Données

Objets sf

Mise au format sf

Manipulation et visualisation

Manipulations spatiales

st_intersects

st_within

Unions

Barycentres et buffers

Exercices

Objets `sf`

Mise au format `sf`

`st_intersects`

`st_within`