25  Manipuler des données spatiales

25.1 Tâches concernées et recommandations

L’utilisateur souhaite traiter avec R des données spatiales (données géolocalisées, polygones…). Si vous ne savez pas si cette fiche répond à votre besoin, les données spatiales sont définies dans le paragraphe Définition des données spatiales.

Tâche concernée et recommandation
  • Il est recommandé d’utiliser le package sf qui couvre les principaux besoins (lecture des formats de données géographiques, traitements des données spatiales, représentations graphiques et cartographiques) ;
  • Sauf cas particuliers, il est recommandé de ne pas utiliser le package sp, qui est rendu obsolète par sf ;
  • Selon le besoin exact, plusieurs autres packages proposent des fonctionnalités complémentaires (voir section Pour aller plus loin).

25.2 Définition des données spatiales

Le terme “données spatiales” désigne les données qui portent sur les caractéristiques géographiques des objets (localisation, contours, liens). Les caractéristiques géographiques des objets sont décrites à l’aide d’un système de coordonnées qui permettent une représentation dans un espace euclidien (\((x,y)\)). Le passage de l’espace réel (la Terre, qui est une sphère) à l’espace plan se fait grâce à un système de projection. Voici quelques exemples de données spatiales :

  • Une table décrivant des bâtiments, avec les coordonnées géographiques de chaque bâtiment ;
  • Le découpage communal du territoire, avec le contour du territoire de chaque commune ;
  • Les routes terrestres, avec les coordonnées décrivant leur parcours.

Les données spatiales rassemblent classiquement deux types de données :

  1. des données géographiques (ou géométries) : objets géométriques tels que des points, des lignes, des polygones, ou des maillages (raster). Exemple : la forme de chaque commune, les coordonnées d’un bâtiment ;
  2. des données attributaires (ou attributs) : des mesures et des caractéristiques associés aux objets géométriques. Exemple : la population de chaque commune, le nombre de fenêtres et le nombre d’étages d’un bâtiment.

Les données spatiales sont fréquemment traitées à l’aide d’un système d’information géographique (SIG), c’est-à-dire un système d’information capable de stocker, d’organiser et de présenter des données alphanumériques spatialement référencées par des coordonnées dans un système de référence (CRS). R dispose de fonctionnalités lui permettant de réaliser les mêmes tâches qu’un SIG (traitement de données spatiales, représentations cartographiques).

Les systèmes de projection font l’objet de standards internationaux et sont souvent désignés par des codes dits codes EPSG . Ce site est un bon aide-mémoire. Les plus fréquents, pour les utilisateurs français, sont les suivants (plus d’infos ici) :

  • 2154 : système de projection Lambert 93. Il s’agit du système de projection officiel : la plupart des données diffusées par l’administration pour la métropole sont disponibles dans ce système de projection.
  • 4326 : WGS 84 ou système de pseudo-Mercator. C’est le système de projection des données GPS.
  • 27572 : Lambert II étendu. Il s’agit de l’ancien système de projection officiel. Les données spatiales anciennes peuvent être dans ce format.

25.3 Comment utiliser le package sf

Le package sf est une boîte à outils conçue pour faciliter la manipulation de données spatiales. La grande force de sf est qu’il permet de manipuler des données spatiales comme s’il s’agissait de données traditionnelles, car il repose sur le standard ISO 19125 simple feature access défini conjointement par l’Open Geospatial Consortium (OGC) et l’International Organization for Standardization (ISO). Cette fiche illustre l’usage de sf avec le jeu de données martinique du package cartography.

Les fonctions de sf sont pour la plupart préfixées par st_ (Spatial Type) et permettent notamment de :

  • Calculer des distances et des surfaces ;
  • Agréger rapidement des zonages (regrouper les communes en département par exemple) ;
  • Trouver dans quelle commune se trouve un bâtiment à partir de ses coordonnées géographiques ;
  • Recalculer des coordonnées dans un autre système de projection.

Le package sf introduit un objet géographique particulier : la table de données spatiale appelée sf. Un sf est simplement une table de données R traditionnelle (un data.frame), enrichie d’une colonne supplémentaire geometry qui contient l’information géographique. Par conséquent, toutes les fonctions qui s’appliquent à un data.frame s’appliquent exactement de la même façon aux attributs des objets sf.

Note

Le package sf est une extension du package dplyr pour les objets géographiques. On peut donc utiliser le pipe (%>%) pour chaîner des opérations avec sf, ce qui est souvent pratique (voir la fiche Manipuler des données avec le tidyverse).

25.3.1 Créer une table de données spatiales

Créer une table de données spatiales avec sf est très simple. Il y a deux méthodes, selon que les données spatiales sont directement disponibles dans un format de données géographiques (par un exemple un fichier shapefile) ou dans un fichier de données tabulées (par exemple un fichier .csv).

25.3.1.1 Importer des fichiers de données géographiques

La fonction sf::st_read() permet de lire différents formats de données géographiques. Les paramètres attendus par la fonction sont :

  • dsn, le nom de la source de donnée (nom du fichier ou chemin vers le répertoire contenant les données géographiques). Le format le plus commun est le .shp mais sf peut lire d’autres formats (.map, .geojson…) ;
  • layer, le nom de la couche. Ce paramètre est facultatif si les données ne contiennent qu’une seule couche ;
  • crs, pour définir le système de projection. Ce paramètre est facultatif. Si le système de projection a été défini à l’écriture du fichier, sf peut le repérer automatiquement à la lecture. Sinon, il faut le définir à la lecture avec l’argument crs.
# Installation du package cartography
# install.packages("cartography")

# Charger les données sur la Martinique
martinique <- sf::st_read(
  dsn = system.file("gpkg/mtq.gpkg", package="cartography"),
  quiet = TRUE)

25.3.1.2 Reconstituer des données géographiques à partir de données brutes

Il peut arriver que les données géographiques ne soient pas stockées dans un format adapté. Par exemple, vous pouvez disposer d’un fichier .csv qui contient une table géolocalisée de logements avec leurs caractéristiques. En ce cas, vous pouvez créer une table de données géographiques en deux temps :

  • Importer les données brutes (ou tabulées, cf + haut) dans R voir la fiche [Importer des fichiers plats (.csv, .tsv, .txt)] ;

  • Transformer le data.frame en objet sf grâce à la commande sf::st_set_geometry(). Par exemple, si les variables lon et lat contiennent des données GPS (code epsg 4326) alors on peut faire :

    df_spatial <- sf::st_as_sf(x = df,                         
             coords = c("lon", "lat"),
             crs = 4326)

25.3.2 Structure d’un objet sf

La commande head() fait apparaître la structure de l’objet sf martinique. L’en-tête précise le type de geometry (polygon, en l’occurrence le contour des communes), et le système de projection (WGS 84 / UTM zone 20N). La table de données contient la liste des communes avec leurs caractéristiques (nom, population…). Enfin, la dernière colonne (geometry) contient l’information géographique (le contour des communes).

head(martinique, 3)
Simple feature collection with 3 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 697602 ymin: 1598818 xmax: 710462 ymax: 1645182
Projected CRS: WGS 84 / UTM zone 20N
  INSEE_COM              STATUS            LIBGEO  POP   MED CHOM  ACT
1     97201 Simple municipality L'Ajoupa-Bouillon 1902 13633  254  801
2     97202 Simple municipality Les Anses-d'Arlet 3737 14558  425 1659
3     97203 Simple municipality      Basse-Pointe 3357 14456  400 1395
                            geom
1 MULTIPOLYGON (((699261 1637...
2 MULTIPOLYGON (((709840 1599...
3 MULTIPOLYGON (((697602 1638...

Pour revenir à un dataframe traditionnel, c’est-à-dire sans la dimension géographique, la manière la plus simple est de rendre nulle la géométrie grâce à la commande sf::drop_geometry. Exemple : martinique <- martinique %>% sf::st_drop_geometry().

On peut produire très rapidement une représentation cartographique d’un attribut avec la fonction plot() :

# Produire une carte de la population par commune
plot(martinique['POP'])

25.3.3 Effectuer des opérations sur un objet sf

25.3.3.1 Effectuer des opérations sur les attributs

Il est possible d’appliquer aux attributs des objets sf les fonctions de la grammaire dplyr comme si l’objet était un data.frame standard. On peut ainsi sélectionner un sous-champ de la base avec dplyr::filter(), faire du tri avec dplyr::arrange() et créer de nouvelles variables avec dplyr::mutate() :

grandes_villes_martinique <- martinique %>% 
  dplyr::filter(POP > 10000) %>% 
  dplyr::arrange(POP) %>% 
  dplyr::mutate(log_pop = log(POP))

On peut également appliquer des fonctions de statistiques descriptives avec dplyr::summarise() et construire des groupes avec dplyr::group_by(). Le group_by() sur un objet spatial va également agréger les géométries : il va créer une nouvelle géométrie recoupant l’ensemble des points, des lignes ou des polygones appartenant au groupe. Par exemple, on peut construire très facilement une géométrie plus agrégée par regroupement des géométries plus fines. Imaginons que la table de données sur la Martinique comprend deux groupes : les villes de plus de 10 000 habitants et celles de moins de 10 000 habitants. On désire faire la somme de la variable POP pour chaque groupe. Dans ce cas, après avoir créé les groupes, on va effectuer un group_by() puis un summarise() :

martinique_agg <- martinique %>%
  dplyr::mutate(group = (POP < 10000)) %>%
  dplyr::group_by(group) %>%
  dplyr::summarise(POP = sum(POP))

head(martinique_agg)
Simple feature collection with 2 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 690574 ymin: 1592536 xmax: 735940.2 ymax: 1645660
Projected CRS: WGS 84 / UTM zone 20N
# A tibble: 2 × 3
  group    POP                                                              geom
  <lgl>  <dbl>                                                <MULTIPOLYGON [m]>
1 FALSE 281471 (((722317 1613892, 722287 1613877, 722225 1613781, 722215 161375…
2 TRUE   99406 (((703149 1624633, 703169 1624533, 703588 1623940, 703848 162345…

Le résultat est analogue à celui qu’on obtiendrait avec undata.frame, à l’exception de la colonne geometry. L’agrégation d’une géométrie de nature POLYGON aboutit à une colonne de format MULTIPOLYGON qui agrège des polygones. La représentation graphique rend plus évidente cette agrégation :

plot(martinique['POP'])

plot(martinique_agg['POP'])

Tip

Les manipulations de données sur un sf sont nettement plus lentes que sur un data.frame traditionnel (car R doit gérer les informations géographiques pendant la manipulation des données). Lorsque vous manipulez des données de grandes dimensions, il peut être préférable d’effectuer les opérations sur les données avant de joindre une géométrie à celles-ci.

25.3.4 Effectuer des opérations sur les géométries d’un objet géométrique

sf propose de nombreuses fonctions pour manipuler la dimension spatiale. On en trouve un certain nombre ici et ici.

Par exemple, on peut calculer la superficie d’un polygone avec sf::st_area(). Si l’unité ne convient pas, il est possible de convertir cette superficie en \(km^2\) :

martinique %>%
  dplyr::mutate(area = sf::st_area(.)) %>%
  dplyr::mutate(x = units::set_units(area, km^2)) %>%
  head()
Simple feature collection with 6 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 695444 ymin: 1598818 xmax: 717731 ymax: 1645182
Projected CRS: WGS 84 / UTM zone 20N
  INSEE_COM              STATUS            LIBGEO  POP   MED CHOM  ACT
1     97201 Simple municipality L'Ajoupa-Bouillon 1902 13633  254  801
2     97202 Simple municipality Les Anses-d'Arlet 3737 14558  425 1659
3     97203 Simple municipality      Basse-Pointe 3357 14456  400 1395
4     97204 Simple municipality         Le Carbet 3683 18847  285 1568
5     97205 Simple municipality       Case-Pilote 4458 21005  347 2096
6     97206 Simple municipality        Le Diamant 5976 19704  430 2654
                            geom           area               x
1 MULTIPOLYGON (((699261 1637... 12176271 [m^2] 12.17627 [km^2]
2 MULTIPOLYGON (((709840 1599... 25783065 [m^2] 25.78307 [km^2]
3 MULTIPOLYGON (((697602 1638... 27753371 [m^2] 27.75337 [km^2]
4 MULTIPOLYGON (((702229 1628... 17806777 [m^2] 17.80678 [km^2]
5 MULTIPOLYGON (((698805 1621... 18609342 [m^2] 18.60934 [km^2]
6 MULTIPOLYGON (((709840 1599... 28042394 [m^2] 28.04239 [km^2]

25.3.5 Joindre des données géographiques et attributaires

25.3.5.1 Jointure sur des attributs

On peut effectuer une jointure sur un objet (ou deux) objet(s) sf en fonction d’une ou plusieurs variables communes. Dans ce cas, on parle de jointure attributaire. Ces opérations peuvent être effectuées comme si les objets sf étaient des data.frames traditionnels avec les fonctions dplyr::_*join.

25.3.5.2 Jointure à partir de la géométrie

Deux objets géométriques peuvent également être associés en fonction de leur dimension spatiale. Par exemple, les fonctions ci-dessous renvoient des TRUE/FALSE pour chaque couple de géométrie des bases x et y (il y en a beaucoup d’autres, voir ?sf::geos_binary_pred) :

Fonction Opération
st_intersects() Quelles géométries de x intersectent celles de y ?
st_contains() Quelles géométries de x contiennent celles de y ?
st_disjoint() Quelles géométries de x sont disjointes à celles de y ?
st_is_within_distance() Quelles géométries de x est à moins de \(X\) m/km de celles de y ?

La fonction sf::st_join() permet d’appliquer ces fonctions pour réaliser la jointure à partir de la dimension spatiale. Deux arguments sont importants :

  • l’argument join = indique la méthode à utiliser pour la jointure. Par défaut, la méthode est st_intersects, ce qui signifie que la table de sortie contient une observation pour chaque couple de géométrie qui a une intersection non vide ;
  • On utilise l’argument left = TRUE si on désire effectuer une jointure sur la gauche (par défaut, un inner join est effectué).

La vignette du package permet d’aller plus loin sur le sujet.

Tip

Les jointures spatiales peuvent être très gourmandes en ressources (car il peut être nécessaire de croiser toutes les géométries de x avec toutes les géométries de y). Voici deux conseils qui peuvent vous aider :

  • Il est préférable de tester les jointures géographiques sur un petit échantillon de données, pour estimer le temps et les ressources nécessaires à la réalisation de la jointure.
  • Il est parfois possible d’écrire une fonction qui réduit la taille du problème. Exemple : vous voulez déterminer dans quelle commune se situe un logement dont vous connaissez les coordonnées et le département ; vous pouvez écrire une fonction qui réalise pour chaque département une jointure spatiale entre les logements situés dans ce département et les communes de ce département, puis empiler les 101 tables de sorties.
Tip

Une application possible de la jointure spatiale est l’obtention de données sur un zonage à façon.

Il est, par exemple, possible d’utiliser les données carroyées mises en ligne sur le site de l’Insee pour obtenir des données sur chaque zone du zonage à façon en transformant des agrégats au niveau carré dans la granularité spatiale désirée.

Une fonction R réalisant une jointure spatiale entre carreaux et zones est déjà disponible sur le site de l’Insee. Le principe de la fonction consiste à déterminer, pour une zone donnée, l’ensemble des carreaux qui la recouvrent puis à calculer les agrégats sur cet ensemble de carreaux. En plus des agrégats, elle permet de calculer, pour chaque zone, la part de la population vivant dans des carreaux imputés, ce qui donne une idée de la fiabilité des résultats obtenus.

25.3.6 Gérer le système de projection

Le système de projection est fondamental pour que la dimension spatiale soit bien interprétée par R. Si un système de projection est défini, il s’affiche dans la ligne Projected CRS lorsqu’on applique la fonction head() à un objet sf. Pour vérifier le système de projection d’une base de données, on peut utiliser sf::st_crs :

sf::st_crs(martinique)
Coordinate Reference System:
  User input: WGS 84 / UTM zone 20N 
  wkt:
PROJCRS["WGS 84 / UTM zone 20N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 20N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-63,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Between 66°W and 60°W, northern hemisphere between equator and 84°N, onshore and offshore. Anguilla. Antigua and Barbuda. Bermuda. Brazil. British Virgin Islands. Canada - New Brunswick; Labrador; Nova Scotia; Nunavut; Prince Edward Island; Quebec. Dominica. Greenland. Grenada. Guadeloupe. Guyana. Martinique. Montserrat. Puerto Rico. St Kitts and Nevis. St Barthelemy. St Lucia. St Maarten, St Martin. St Vincent and the Grenadines. Trinidad and Tobago. Venezuela. US Virgin Islands."],
        BBOX[0,-66,84,-60]],
    ID["EPSG",32620]]

Les deux principales fonctions pour définir le système de projection utilisé sont :

  • sf::st_set_crs : cette commande sert à préciser quel est le système de projection utilisé, c’est-à-dire comment les coordonnées (x,y) sont reliées à la surface terrestre. Cette commande ne doit pas être utilisée pour transformer le système de coordonnées, seulement pour le définir.
  • sf::st_transform : cette commande sert à projeter les points d’une géométrie dans une autre, c’est-à-dire à recalculer les coordonnées selon un autre système de projection. Par exemple, si on désire produire une carte avec un fond openstreetmaps ou une carte dynamique leaflet à partir de données projetées en Lambert 93, il est nécessaire de re-projeter les données dans le système WGS 84 (code EPSG 4326).

La définition du système de projection se fait de la manière suivante :

martinique <- martinique %>% sf::st_set_crs(32620)

Alors que la reprojection s’obtient de la manière suivante :

martinique <- martinique %>% sf::st_transform(4326)

25.4 Pour aller plus loin

Selon le besoin exact, plusieurs autres packages proposent des fonctionnalités complémentaires :

  • le package raster pour gérer le format de données de type maillage ;
  • les packages pour l’analyse statistique spatiale :
    • spdep pour l’économétrie spatiale : relations de voisinage entre objets spatiaux, indices d’autocorrélation spatiale, … ;
    • spatstat pour l’économétrie spatiale ;
    • gstat et geoR pour la géostatistique ;
  • btb (produit par l’Insee) pour du lissage ;
  • le package rpostgis pour interfacer R à une base de données PostGIS.

25.5 Où trouver des données géographiques et des shapefiles ?

Il s’agit du domaine où la diffusion de données par l’open data a permis une grande diversité des sources disponibles. Certaines données sont disponibles sous forme de shapefiles, d’autres avec des identifiants géographiques comme le code commune à apparier à un shapefile proposant la même variable.

Les limites administratives de référence sont disponibles sur le site Admin Express de l’IGN (anciennement GeoFla).

Le code officiel géographique (COG), qui tient à jour les entités administratives (codes communes, régions, etc.) est disponible sur le site de l’Insee et via l’API Métadonnées. La fiche API indique comment accéder à des données via une API.

Tip

Contrairement à ce qui pourrait être pensé, la géographie et le COG sont régulièrement modifiés, pour prendre notamment en compte des fusions et scissions de communes. Deux ensembles apparemment identiques de codes communes au 1er janvier 2022 et 1er janvier 2021 peuvent ainsi différer, et il convient de retraiter ses données pour s’assurer qu’elles sont toutes définies dans une même géographie.

Le package COGugaison fournit un ensemble d’outils répondant à ce besoin. Il permet de réaliser de nombreuses modifications utiles au chargé d’études (remplacement des codes arrondissements dans Paris, Lyon, Marseille, identification du millésime géographique des données, visualisation des changements de géographie, transformation d’un millésime à un autre pour les communes, les zonages standards d’études de l’Insee, et les zonages à façon, etc.) sans avoir à mobiliser le COG.

Note

Ce package n’étant pas disponible sur le CRAN, dans un environnement connecté à internet, il est nécessaire de l’installer depuis Github:

remotes::install_github("antuki/COGugaison")
Spécificité Insee

Ce package n’est pas disponible sur le CRAN. Sur AUS, on peut l’installer avec la commande install.packages("COGugaison", repos = "https://nexus.insee.fr/repository/r-public").

Il est également possible de trouver des données géographiques sur data.gouv, insee.fr ou sur d’autres ressources d’open data, qu’ils soient génériques comme opendatasoft ou plus institutionnels comme opendata.paris.fr/

Spécificité Insee

L’Insee propose un outil pour sélectionner et télécharger des données géographiques via l’application creacartes.

25.6 Pour en savoir plus

  • la documentation du package sf :
  • un tutoriel détaillé sur les données spatiales et sur la cartographie avec R (en français).