Chapitre 4 Lire et écrire des données spatiales avec R

Dans ce chapitre, nous allons utiliser les packages suivants :

# CRAN
library(DT)
library(remotes)
library(RPostgres)
library(sf)
library(spData)
library(tidyverse)

# Gitlab-forge
library(datalibaba)

# Github
library(COGiter)
library(CARTElette)

Le package {sf} s’appuie sur la librairie externe GDAL (Geospatial Data Abstraction Library) pour lire et écrire les données spatiales. GDAL est une librairie permettant de lire et d’écrire des données vecteur et raster de n’importe quel format de fichier ou base de données.

Les fonctions de lecture et d’écriture de {sf} s’appellent st_read() et st_write().

Pour lire des données spatiales, st_read() a besoin :

  • d’une source de donnée : un fichier, un répertoire, une base de données ;
  • d’un layer: une table de données spatiale spécifique du fichier, répertoire, ou de la base de données ;

On peut avoir besoin de rajouter des options spécifiques au format de la source. Par exemple, sur des données en csv, il faut pouvoir spécifier la composante spatiale des données (les champs X, Y ou longitude, latitude).

4.1 Lire des fichiers plats

st_read() permet de lire des données spatiales disponibles en fichier plat.

Le format de fichier plat le plus populaire en géomatique est ESRI Shapefile. Ce format, en plus de ne pas être un format ouvert, a des limites bien documentées3.

Avec l’avènement du web, le format GeoJSON se développe beaucoup, bien qu’il soit aussi limité.

Le format considéré comme le plus prometteur est le format OGC GeoPackage promu par l’Open Geospatial Consortium4.

Mais la liste des formats lisibles par {sf} est bien plus vaste.

Pour l’obtenir, on peut utiliser st_drivers(), qui liste les drivers utilisables par {sf} :

DT::datatable(st_drivers() %>% arrange(name))

Exemple de lecture d’une donnée au format GeoPackage: les données sur la Leucémie à New York disponibles dans le package spData5.

geo_fichier <- system.file("shapes/NY8_bna_utm18.gpkg", package = "spData")
geo_fichier # contient l'adresse vers le jeu de données
NY_leukemia <- st_read(dsn = geo_fichier)

4.2 Ecrire des fichiers plats

Ecrire des fichiers plats va se faire avec la fonction st_write().

st_write(obj = NY_leukemia, dsn = "extdata/NY_leukemia.gpkg")

Notez que si vous cherchez à exporter une nouvelle fois ce fichier, vous aurez un message d’erreur.

Pour pouvoir écraser ce fichier, vous avez deux options :

  • Utiliser le paramètre append qui sert à d’ajouter des données à la couche existante quand il est à TRUE, et la remplace sinon.
st_write(
  obj = NY_leukemia,
  dsn = "extdata/NY_leukemia.gpkg",
  append = FALSE
)
  • Utiliser le paramètre delete_layer = TRUE de st_write() qui vous permet d’écraser les layers avant sauvegarde (paramètre qui ne dépend pas du driver utilisé).
st_write(
  obj = NY_leukemia,
  dsn = "extdata/NY_leukemia.gpkg",
  delete_layer = TRUE
)

A noter que st_write() peut exporter des tables géographiques au format csv :

st_write(
  obj = NY_leukemia %>% st_centroid(),
  dsn = "extdata/NY_leukemia.csv",
  layer_options = "GEOMETRY=AS_XY",
  delete_layer = TRUE
)

Il faut juste lui préciser comment on souhaite écrire la géométrie dans la table exportée, avec layer_options = "GEOMETRY=AS_WKT" ou layer_options = "GEOMETRY=AS_XY".

Le paramètre layer_option permet d’inclure des options propres à chaque driver utilisé.

4.3 Lire des données spatiales du web

4.3.1 Lire des données geojson

sf peut lire un geojson qu’il soit hébergé sur un serveur, compressé ou non, ou encore disponible par api.

4.3.1.1 Exemple 1 (geojson classique):

st_read("https://france-geojson.gregoiredavid.fr/repo/regions.geojson") %>% 
  filter(code > "10") %>% 
  select(code) %>% 
  plot()
Reading layer `regions' from data source 
  `https://france-geojson.gregoiredavid.fr/repo/regions.geojson' 
  using driver `GeoJSON'
Simple feature collection with 18 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -61.80596 ymin: -21.3869 xmax: 55.83663 ymax: 51.08854
Geodetic CRS:  WGS 84

4.3.1.2 Exemple 2 (api)

L’api de découpage administratif permet de charger des contours administratifs. Dans l’exemple suivant, nous chargeons les contours des epci du département 44 en complétant l’URL de base avec des paramètre de requête.

url <- 'https://geo.api.gouv.fr/epcis?codeDepartement=44&format=geojson&geometry=contour'
st_read(url) %>%
  select(code) %>% 
  plot()
Reading layer `OGRGeoJSON' from data source 
  `https://geo.api.gouv.fr/epcis?codeDepartement=44&format=geojson&geometry=contour' 
  using driver `GeoJSON'
Simple feature collection with 17 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -2.625436 ymin: 46.86008 xmax: -0.896293 ymax: 47.87169
Geodetic CRS:  WGS 84

4.3.1.3 Exemple 3 (geojson compressé)

sf repose sur GDAL et permet également de lire des geojsons compressés. On pourrait faire de même avec un shapefile compressé en .zip. Cela évite de télécharger le fichier et de le décompresser. Voici un exemple pour récupérer les bâtiments de la commune de Nantes provenant du PCI redistribué par Etalab.

dsn <- '/vsigzip//vsicurl/https://cadastre.data.gouv.fr/data/etalab-cadastre/latest/geojson/communes/44/44109/cadastre-44109-batiments.json.gz'
batiment_nantes <- st_read(dsn = dsn, layer = "cadastre-44109-batiments.json") 
Reading layer `cadastre-44109-batiments.json' from data source 
  `/vsigzip//vsicurl/https://cadastre.data.gouv.fr/data/etalab-cadastre/latest/geojson/communes/44/44109/cadastre-44109-batiments.json.gz' 
  using driver `GeoJSON'
Simple feature collection with 81539 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -1.627776 ymin: 47.18209 xmax: -1.485699 ymax: 47.29455
Geodetic CRS:  WGS 84
mapview(batiment_nantes %>% head(1000))

Pour aller plus loin sur les Virtual File Systems : https://gdal.org/user/virtual_file_systems.html

4.3.2 Lire des données WFS

sf permet aussi de lire des données en flux WFS.

crte <- sf::st_read(dsn = 'https://datacarto.sigloire.fr/wfs?REQUEST=getCapabilities&service=WFS&VERSION=2.0.0', layer = "ms:r_portrait_crte_r52")
Reading layer `ms:r_portrait_crte_r52' from data source 
  `https://datacarto.sigloire.fr/wfs?REQUEST=getCapabilities&service=WFS&VERSION=2.0.0' 
  using driver `WFS'
Simple feature collection with 72 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 276123.7 ymin: 6582709 xmax: 545050 ymax: 6834664
Projected CRS: RGF93 v1 / Lambert-93
crte %>% select(nom_territ) %>% plot()

Autres exemple avec la géoplateforme de l’IGN :

pref_dep <- st_read("https://data.geopf.fr/wfs/ows?SERVICE=WFS&VERSION=2.0.0&REQUEST=GetCapabilities", 
                    layer = "LIMITES_ADMINISTRATIVES_EXPRESS.LATEST:chef_lieu_de_departement") %>% 
  select(nom_du_chef_lieu)
Reading layer `LIMITES_ADMINISTRATIVES_EXPRESS.LATEST:chef_lieu_de_departement' from data source `https://data.geopf.fr/wfs/ows?SERVICE=WFS&VERSION=2.0.0&REQUEST=GetCapabilities' 
  using driver `WFS'
Simple feature collection with 101 features and 10 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -61.72203 ymin: -20.87582 xmax: 55.44641 ymax: 50.63436
Geodetic CRS:  WGS 84
mapview(pref_dep, legend = FALSE)

On indique au niveau du paramètre dsn l’url du webservice WFS de la plateforme et au niveau du parametre layer le nom de la couche. Pour aller plus loin dans l’exploitation des flux WFS dans R : https://inbo.github.io/tutorials/tutorials/spatial_wfs_services/.

4.4 Lire/écrire des données spatiales sur PostGIS

PostGIS, c’est l’extension géomatique de PostgreSQL permettant de stocker des données géo sur un serveur de données et de les exploiter.

R vous permet de vous connecter simplement à une base PostGIS. Pour cela vous aurez besoin du package {DBI} (package permettant de s’interfacer à des bases de données) et du package {RPostgres} (interface à PostgreSQL en particulier).

Rstudio a développé tout un site sur les interactions entre les bases de données et R, vous pouvez le trouver à l’adresse suivante : db.rstudio.com.

Rstudio possède même depuis la version 1.1 un onglet de connexion vous permettant de visualiser vos connexions aux bases.

Pour se connecter à une base, il faut définir un driver (à quelle type de base on veut se connecter et selon quel protocole) et un connecteur (les informations de connexion).

drv <- dbDriver("Postgres")
con <- dbConnect(drv,
  dbname = "nom_de_ma_db",
  host = "adresse_ip_du_serveur",
  port = numero_du_port,
  user = "nom_utilisateur",
  password = "mot_de_passe_super_secret"
)

Pour lire des données géographique sur le serveur, il faut utiliser encore {sf} et sa fonction st_read() en lui définissant 2 paramètres : le connecteur et la requête que l’on veut réaliser.

ma_table <- st_read(con, query = "SELECT * FROM le_schema.ma_table")

L’avantage ici est que vous pouvez faire travailler le serveur directement en SQL sans avoir à faire travailler votre poste de travail.

Vous pouvez très bien, dans cette requête sql, réaliser quelques filtres, sélections et agrégations. Celles-ci seront alors réalisées par le serveur posgreSQL et non par votre session R qui ne récupérera que le résultat.

Vous pouvez écrire vos données ensuite de la même façon avec st_write()

st_write(ma_table,
  dsn = con,
  layer = Id(schema = "schema", table = "ma_table")
)

Un package, encore expérimental, de la DREAL Pays de la Loire, {datalibaba} simplifie l’écriture des instructions de lecture/écriture de données vers ou depuis un SGBD Postgresql/postgis.

Il propose de stocker vos identifiants de connexion dans vos variables d’environnement, c’est à dire en dehors de votre script, afin d’en préserver la confidentialité, mais aussi de vous éviter de les réécrire dans chaque script.

L’utilisateur n’a plus à se préoccuper du driver de connexion ni de la fonction de lecture ou de comment placer le nom du schéma par rapport au nom de la table.

Les types propres à R, comme les facteurs, sont préservés au ré-import.

Les instruction précédentes deviennent :

remotes::install_gitlab('dreal-pdl/csd/datalibaba', host = "gitlab-forge.din.developpement-durable.gouv.fr")
library(datalibaba)

poster_data(data = regions_geo, schema = "admin_express", table = "n_regions", 
            pk = 'REG', db = 'referentiels')

station <- importer_data(table = "station", schema = "pesticides", base = "production")

Pour tester ces fonctions, le SSP Cloud permet de lancer des instances de postgresql/postgis et de les interoger depuis une session RStudio.

4.5 Convertir un dataframe

Une autre façon d’obtenir un spatial dataframe est de convertir un objet existant en objet sf. Par exemple, partir d’un dataframe sur lequel on va définir à la fois la ou les colonnes contenant la dimension géographique et le crs de référence. Ou encore convertir un objet de type sp en objet de type sf.

Pour cela, la fonction st_as_sf() permet de convertir un objet en objet sf en définissant la composante spatiale.

La fonction st_set_crs() permet de définir le crs de notre objet.

Exemple, nous allons lire les coordonnées géographiques des préfectures de région.

prefectures <- read_csv2("extdata/prefecture.csv")
prefectures_geo <- st_as_sf(prefectures, coords = c("x", "y"), crs = 2154)

4.6 Données administratives contenues dans des packages

Depuis le développement de {sf}, de nombreux packages proposent des données géographiques.

Pour la France, les packages {CARTElette} et {COGiter} propose des fonds de cartes aux contours administratifs.

Ils s’installent depuis github :

remotes::install_github("MaelTheuliere/COGiter")
remotes::install_github("antuki/CARTElette/CARTElette@RPackage")

{CARTElette} propose de télécharger les différents fonds de plan administratifs d’admin-Express COG de l’IGN, de tous les millésimes jusque 2021.

charger_carte(COG = 2015, nivsupra = "REG", geometrie_simplifiee = FALSE) %>% 
  st_geometry() %>% 
  plot()
Reading layer `REG_2015_CARTElette' from data source 
  `/tmp/RtmpWJvFEO/REG_2015_CARTElette.shp' using driver `ESRI Shapefile'
Simple feature collection with 28 features and 2 fields (with 1 geometry empty)
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 99217.1 ymin: 5988964 xmax: 1242417 ymax: 7110480
Projected CRS: Lambert_Conformal_Conic

{COGiter} contient le dernier millésime des fonds de plan administratifs. Ces couches sont obtenues par simplification des couches admin-express COG de l’IGN et sont adaptées à la visualisation d’indicateurs statistiques.

Il comprend pour la métropole et chaque DROM un fond de carte pour les niveaux communes, EPCI, départements, régions :

ls("package:COGiter", pattern = "_geo$")
 [1] "communes_971_geo"       "communes_972_geo"       "communes_973_geo"      
 [4] "communes_974_geo"       "communes_976_geo"       "communes_geo"          
 [7] "communes_metro_geo"     "departements_971_geo"   "departements_972_geo"  
[10] "departements_973_geo"   "departements_974_geo"   "departements_976_geo"  
[13] "departements_geo"       "departements_metro_geo" "epci_971_geo"          
[16] "epci_972_geo"           "epci_973_geo"           "epci_974_geo"          
[19] "epci_976_geo"           "epci_geo"               "epci_metro_geo"        
[22] "filtrer_cog_geo"        "regions_971_geo"        "regions_972_geo"       
[25] "regions_973_geo"        "regions_974_geo"        "regions_976_geo"       
[28] "regions_geo"            "regions_metro_geo"

La fonction filtrer_cog_geo() vous permet d’obtenir une liste de spatial dataframes des différentes échelles, centrée sur une partie du territoire (que ce soit une commune, un EPCI, un département ou une région) dès lors que vous connaissez son code officiel.

Exemple sur la Normandie :

normandie <- filtrer_cog_geo(reg = "28")
glimpse(normandie)
List of 4
 $ communes    : sf [2,644 × 3] (S3: sf/tbl_df/tbl/data.frame)
  ..$ DEPCOM  : chr [1:2644] "14380" "27111" "14381" "14184" ...
  ..$ AREA    : num [1:2644] 3100000 3800000 5500000 3800000 2400000 5600000 4200000 4300000 4700000 7200000 ...
  ..$ geometry:sfc_GEOMETRY of length 2644; first list element: List of 1
  .. ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
  ..- attr(*, "sf_column")= chr "geometry"
  ..- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA
  .. ..- attr(*, "names")= chr [1:2] "DEPCOM" "AREA"
 $ epci        : sf [71 × 3] (S3: sf/tbl_df/tbl/data.frame)
  ..$ EPCI    : Factor w/ 1254 levels "200000172","200000438",..: 18 38 95 96 113 118 157 226 231 244 ...
  ..$ AREA    : Units: [m^2] num [1:71] 5.74e+08 6.64e+08 2.79e+08 3.65e+08 5.68e+08 ...
  ..$ geometry:sfc_POLYGON of length 71; first list element: List of 1
  .. ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
  ..- attr(*, "sf_column")= chr "geometry"
  ..- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA
  .. ..- attr(*, "names")= chr [1:2] "EPCI" "AREA"
 $ departements: sf [5 × 3] (S3: sf/tbl_df/tbl/data.frame)
  ..$ DEP     : Factor w/ 101 levels "01","02","03",..: 14 26 51 62 77
  ..$ AREA    : Units: [m^2] num [1:5] 5.53e+09 6.01e+09 5.95e+09 6.10e+09 6.28e+09
  ..$ geometry:sfc_POLYGON of length 5; first list element: List of 1
  .. ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
  ..- attr(*, "sf_column")= chr "geometry"
  ..- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA
  .. ..- attr(*, "names")= chr [1:2] "DEP" "AREA"
 $ regions     : sf [1 × 3] (S3: sf/tbl_df/tbl/data.frame)
  ..$ REG     : Factor w/ 18 levels "01","02","03",..: 9
  ..$ AREA    : Units: [m^2] num 2.99e+10
  ..$ geometry:sfc_POLYGON of length 1; first list element: List of 1
  .. ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
  ..- attr(*, "sf_column")= chr "geometry"
  ..- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA
  .. ..- attr(*, "names")= chr [1:2] "REG" "AREA"
plot(normandie$epci)

Le package COGiter propose également :

  • les tables du dernier millesime du COG de l’Insee,
  • une table de passage des COG historiques vers le COG à jour,
  • des fonctions d’aide à au passage de jeux de données vers le millésime actuel du COG,
  • des fonctions de calculs d’agrégats aux différentes échelles territoriales,
  • des fonctions d’aide à la sélection des différents fond de cartes nécessaire à la mise en page de cartes statistiques.

Son utilisation est désormais vue lors du module 2.

Il est traditionnellement mis à jour courant avril chaque année.

  1. voir http://switchfromshapefile.org/↩︎

  2. https://www.geopackage.org/↩︎

  3. Issues de Waller and Gotway (2004) Applied Spatial Statistics for Public Health Data.↩︎