Cartographie et géomatique avec R (2023)

Deux modules R principaux servent aujourd’hui pour le traitement de données géomatiques:

• le module sf (simple features) dédié aux géodonnées de type vecteur, qui rassemble les développements antérieurs réalisés notamment dans les modules sp, rgdal et rgeos.
• Le module terra, dédié aux géodonnées de type raster.

Pour comprendre ces modules R en détail, je vous recommande l’excellente documentation de Timothée Giraud & Hugues Pecout. Le présent exercice vise à vous familiariser avec les éléments de base de la cartographie vectorielle avec R. Le prérequis pour comprendre les instructions ci-dessous est de connaître les bases de l’utilisation de R et de RStudio. Vous devriez aussi savoir ce que sont les projections cartographiques. Comme il s’agit d’un exercice avancé, je ne donne pas l’ensemble des solutions, afin que vous appreniez à utiliser la documentation de modules de R. Suivez les hyperliens pour mieux comprendre chaque fonction mentionnée. Le chapitre 3 de la documentation de Giroud & Pecout peut aussi vous servir de repère et de référence. Je vous recommande également de garder ouvert, dans un onglet du navigateur l’index de la page de documentation du module sf pour trouvez facilement des informations sur les fonctions utilisées.

Installation et chargement des modules nécessaires

• Installez le module sf à l’aide de RStudio.
• Activez-le, ainsi que d’autres modules qui nous seront utiles par la suite. Installez ceux qui vous manqueraient:

library(sf)
library(readxl) # permet de lire les anciens fichiers Excel
library(magrittr) # permet d'utiliser l'opérateur "tuyau" (<em>pipe</em>) %>% 
library(ggplot2) # permet de visualiser des donnéesCode language: PHP (php)

Chargement et inspection des données

• Récupérez les données d’exercice unine_exercice1.zip depuis le site du cours, stockez-les dans un répertoire du disque dur de votre ordinateur et décompressez-les. Trouvez le fichier “Poly.shp” dans le dossier décompressé et notez bien son adresse.
• Importez le fichier Poly.shp dans une variable nommée compoly à l’aide de la fonction st_read. [p. ex. compoly <- st_read("Poly.shp")]
• Vérifiez le contenu de vos géodonnées en écrivant compoly dans le terminal, puis en tapant Enter. Vous devriez voir ceci:

Comme vous le voyez, les données que vous avez chargées en mémoire contiennent 2715 objets spatiaux (features), de type “multipolygones” (des polygones complexes composés chacun d’un ou plusieurs polygones). La projection géométrique (Coordinate Reference System, CRS) est inconnue (NA). Nous y remédierons de suite. Notons, avant ça, que nous avons trois colonnes de données: ID0 (l’identifiant numérique de la commune), ID1 (le nom de la commune) et geometry. Cette dernière colonne contient la géométrie du territoire de la commune, stockée sous forme d’une expression du langage WKT (Well Known Text). La géométrie de la commune de Neuchâtel, par exemple, est décrite comme suit:

MULTIPOLYGON (((559008 206103, 560028 206623, 560659 207107, 563408 209215, 565496 212356, 566048 211292, 563943 208665, 563363 208160, 564112 206640, 564283 206294, 563630 205740, 563586 205431, 563268 205128, 562770 205123, 562638 204836, 561734 204520, 561293 204263, 560962 204351, 560169 203910, 560036 203690, 559595 203381, 558613 203283, 558528 203753, 559008 206103)))

L’image ci-dessous permet de comprendre que cette expression WKT génère un seul polygone en commençant au point (559008, 206103), ajoutant les points de la liste dans les sens des aiguille d’une montre, et en finissant au même point (559008, 206103):

Le grand avantage du langage de balisage WKT et de pouvoir faire tenir une forme géométrique complexe dans une seule cellule d’un tableau de données. Cela veut dire que vous pourriez stocker une carte vectorielle complète même dans un fichier Excel. Cette approche vous permet de traiter vos géodonnées comme une data frame R quelconque, tout en conservant la possibilité de générer des cartes.

Nul besoin, en outre, de vous occuper de l’interprétation de la géométrie WKT lorsque vous utilisez le module sf; il suffit d’utiliser, par exemple, les commandes ggplot() et geom_sf() fournies par le module ggplot2 pour mapper les toutes les communes listées dans notre tableau compoly:

ggplot(compoly) + geom_sf() Code language: R (r)

Revenons au problème de la projection CRS. J’ai déjà évoqué que sa définition manquait. Vérifiez-le avec la commande st_crs(compoly) qui retourne, en effet, “Coordinate Reference System: NA“. Ceci nous posera un problème plus tard si nous essayons de combiner notre carte avec des géodonnées issues d’autres sources. Mais nous pouvons déclarer le bon CRS. Le graticule de l’image ci-dessus suggère que nous avons affaire à un système de coordonnées suisse. Connaissant la source de ces données, je sais qu’il s’agit de la projection CH1903 / LV03 dont le numéro EPSG international est 21781. Nous pouvons l’assigner ainsi:

compoly <- st_set_crs(compoly,21781)

Dès lors, lorsque nous redemandons à R de dessiner la carte avec la commande ggplot(compoly) + geom_sf(), nous obtenons une carte des communes suisses correctement placée sur la graticule globale des longitudes et des latitudes:

Ajout de données à partir d’une autre source tabulaire

Le moment est venu d’associer d’autres données à notre carte. Commençons par quelques données tabulaires issues du recensement fédéral de la population 2000.

• Importez les données Excel donnees_communes.xls dans la variable data à l’aide du module readxl et sa fonction read_excel().
• Joignez les données désormais stockées dans la variable data à votre géométrie stockée dans la variable compoly à l’aide de la fonction merge():

compoly <- merge(x=compoly,y=data, by.x="ID0", by.y="GMDE")Code language: R (r)

Si vous inspectez compoly, vous devriez voir ceci:

Vous pouvez à présent utiliser les variables que vous venez d’ajouter à votre géométrie. Si vous consultez la documentation du dataset dans le fichier “donnees_communes.pdf”, vous verrez que la variable P00B23 contient le nombre d’italophones et que la variable P00BTOT contient la population totale d’une commune. Calculons la proportion d’italophones et montrons-la sur une carte.

compoly$p_P00B23 <- compoly$P00BTOT / compoly$P00B23
ggplot(compoly) + geom_sf(aes(fill=p_P00B23)) Code language: R (r)

La carte révèle des proportions élevées d’italophones en dehors du Tessin et des Grisons. Découvrez de quelles communes il s’agit avec la commande suivante:

compoly[which(compoly$p_P00B23 > 0.9 & ! compoly$KT %in% c(18,21)),c("KT","ID1","p_P00B23")]Code language: PHP (php)

La force du module sf tient dans le fait qu’il permet d’accomplir également des opérations spatiales.

• Calculez la superficie des communes à l’aide de la fonction st_area().
• Servez-vous du résultat et de la variable P00BTOT pour calculer leur densité de la population par km^2. (compoly$densite)
• Essayez d’afficher le résultat sous forme de carte

Il est plus que probable que cette tentative vous renvoie une erreur. Comme d’habitude en programmation et, quel que soit le langage utilisé, lisez le message d’erreur! Ces messages vous donnent souvent la solution aux problèmes que vous rencontrez.

library(units)
ggplot(compoly) + geom_sf(aes(fill=densite))

Hélas, dans ce cas, la solution proposée ne résout pas le problème:

Il ne faut pas désespérer, mais tenter de comprendre ce qui arrive. Examinons la variable compoly$densite:

Vous voyez que cette variable comporte une précision sur les unités utilisées: [1/m^2]. En effet, nous avons affaire au nombre de résidents par mètre carré. Mais la commande geom_sf(aes(fill=densite)) ne parvient pas à traiter cette indication en l’état actuel des modules “ggplot2” (version 3.4.0) et “units” (version 0.8-1). De telles choses arrivent, en attendant que les modules soient mis à jour pour une meilleure compatibilité. Pour contourner le problème, il suffit de convertir la valeur de la densité en simple valeur numérique dépourvue d’unités:

  ggplot(compoly) + geom_sf(aes(fill=densite%>%as.numeric)) Code language: JavaScript (javascript)

Nous avons pu résoudre le problème des unités, mais la carte résultante est peu différenciée. De toute évidence, peu de communes suisses sont aussi densément peuplées que Genève, ce qui les met toutes au même niveau sur une échelle de couleurs linéaire. Passons à une échelle logarithmique avec l’instruction ggplot scale_fill_gradient(trans="log"). Transformons au passage les “habitants par mètre carré” en “habitants par km carré” en multiplant la densité par 1000000:

ggplot(compoly) + 
	geom_sf(aes(fill=(densite%>%as.numeric) * 1000000)) + 
	scale_fill_gradient(trans="log", name="résidents / km^2") Code language: R (r)

Notez que le fond de cartes à notre disposition limite l’étendue des communes à l’espace situé en deçà des 1000 mètres d’altitude. Lorsque vous réfléchissez à la densité de peuplement d’un territoire, il est souvent pertinent de focaliser cette réflexion sur la partie de territoire où la population a tendance à se concentrer. Visitez le village de Zermatt pour vous rendre compte que sa densité peut être plutôt élevée, malgré les sommets alpins qui l’entourent…

Ajout de données à partir d’OpenStreetMap

Données vectorielles

Ajoutons, maintenant, d’autres géodonnées en puisant dans des ressources en ligne. Focalisons-nous, ce faisant, sur le canton du Tessin, dont le numéro de canton (KT) est 21 (vous pouvez vous en rendre compte en examinant vos données).

compoly.tessin <- compoly[which(compoly$KT==21),]
ggplot(compoly.tessin) +
	geom_sf(aes(fill=(densite%>%as.numeric) * 1000000)) + 
	scale_fill_gradient(trans="log", name="résidents / km^2") Code language: PHP (php)

Nous allons ajouter des données directement depuis OpenStreetMap (OSM).

• Pour cela, installez le module “osmdata” à l’aide de RStudio et activez-le.

La base de données OSM couvre le monde entier. Au mois de mars 2023, la taille du fichier qui les contient toutes atteint 67.7 GB ! Vous n’avez évidemment pas envie de le télécharger dans sa totalité. Avant d’entamer le téléchargement de quoi que ce soit depuis le serveur d’OSM, il faut définir la zone d’intérêt à l’aide d’une “bounding box“. Nous utiliserons les fonctions st_bbox et opq. Profitons-en également pour étendre le temps que nous voulons laisser à nos requêtes au serveur; par défaut, ce temps est limité à 25 secondes, ce qui ne suffira pas toujours pour récupérer des géodonnées d’un canton entier. Pensez par exemple aux géométries de tous les bâtiments du Tessin!

compoly.tessin.mercator <- st_transform(compoly.tessin,4326) 
bbox <- st_bbox(compoly.tessin.mercator)
q <- opq(bbox, timeout=60) 

Notez que nous avons, au préalable créée une nouvelle couche de données compoly.tessin.mercator en reprojetant la couche précédemment créée, compoly.tessin, dans le système de projection (CRS) 4326. Ceci est nécessaire car la fonction st_bbox() a besoin que les données soient dans le système de coordonnées EPSG:4326 pour rendre des valeurs correctement interprétables par la fonction opq() issue du module osmdata. Ce module va en effet chercher en ligne des données OpenStreetMap stockées avec ce système de coordonnées précis. Si vous oubliez de convertir le CRS, opq() essaiera probablement de télécharger le monde entier…

À présent, récupérez les données des stations de train:

req <- add_osm_feature(opq = q, key = 'railway', value = "station")
res <- osmdata_sf(req) %>% unique_osmdataCode language: JavaScript (javascript)

Examinez le résultat:

Vous voyez que les données contiennent 74 points, 8 polygones et 8 polylignes. Stockons les points dans une nouvelle variable:

stations.train <- res$osm_pointsCode language: PHP (php)

Récupérons aussi les lignes ferroviaires:

req <- add_osm_feature(opq = q, key = 'route', value = "train")
res <- osmdata_sf(req) %>% unique_osmdata
routes.train <- res$osm_multilinesCode language: R (r)

Et cartographions le tout:

ggplot() + 
	geom_sf(data=compoly.tessin, aes(fill=(densite%>%as.numeric) * 1000000)) + 
	scale_fill_gradient(trans="log", name="résidents / km^2\n(communes du Tessin)") + 
	geom_sf(data=stations.train, color="red") +
	geom_sf(data=routes.train, color="red") Code language: R (r)

Notre problème est que, même si la zone de requête est limitée à la bounding box du Tessin, certaines lignes s’étendent jusqu’à Genève, Francfort, ou Milan. Il faut donc couper les données à l’aide de la fonction st_crop.

routes.train.cropped <- st_crop(routes.train,bbox)
ggplot() + 
	geom_sf(data=compoly.tessin, aes(fill=(densite%>%as.numeric) * 1000000)) + 
	scale_fill_gradient(trans="log", name="résidents / km^2\n(communes du Tessin)") + 
	geom_sf(data=stations.train, color="red") +
	geom_sf(data=routes.train.cropped, color="red") Code language: JavaScript (javascript)

• Ajoutez à présent par vous-mêmes d’autres données depuis le serveur OSM. Par exemple les parcs nationaux, ou les écoles.

Fonds de carte

La base de données d’OpenStreetMap permet aussi de générer des fonds de carte. Des services en ligne comme Stamen Maps proposent des fonds de carte déjà générés et prêts à l’usage. Vous pouvez les intégrer dans une carte préparée avec R à l’aide du module ggmap.

library(ggmap)
fondsCarte <- st_buffer(compoly.tessin.mercator,1000) %>% 
	st_bbox() %>% 
	setNames(c("left", "bottom", "right", "top")) %>% 
	get_stamenmap(maptype="terrain",zoom=11)
ggmap(fondsCarte) + 
	geom_sf(
		data=compoly.tessin.webmercator, 
		inherit.aes = FALSE, 
		fill="white",
		color="white",
		alpha=0.5
	) +
	geom_sf(data=stations.train, inherit.aes = FALSE)  Code language: PHP (php)

Ce bout de code demande quelques explications. Notez d’abord que nous utilisons les données reprojetées compoly.tessin.mercator. À l’instar de la fonction opq() du module osmdata, la fonction get_stammenmap() du module ggmap s’attend en effet à recevoir une bonding box exprimées en coordonnées EPSG:4326.

La fonction st_buffer() sert à créer une marge autour de notre géométrie. De cette manière, les polygones représentant les communes du Tessin ne colleront pas au bord de la carte. Nous avons déjà vu la fonction st_bbox qui fournit une bounding box sous forme d’un “vecteur nommé”:

     xmin      ymin      xmax      ymax 
 8.459629 45.808788  9.135888 46.563882 Code language: CSS (css)

Hélas, get_stamenmap() a besoin que les coordonnées de la bounding box portent d’autres noms, ce que nous pouvons heureusement vite régler en appliquant

st_bbox() %>%  setNames(c("left", "bottom", "right", "top"))Code language: JavaScript (javascript)

pour obtenir

     left    bottom     right       top 
 8.459629 45.808788  9.135888 46.563882 Code language: CSS (css)

Notons enfin que nous introduisons le paramètre inherit.aes=FALSE dans les deux occurrences de l’appel à la fonction geom_sf(). Cela remédie à un problème d’incompatibilité partielle entre les modules sf et ggmap. C’est ainsi que nous obtenons une carte dotée d’un fond.

Exercices complémentaires

Le module sf vous permet d’accomplir des opérations spatiales, comme de calculer la moyenne des valeurs d’une variable attachée à tous les points contenus dans un polygone. Les possibilités sont nombreuses et je vous recommande de consulter la documentation ou le tutoriel de Giraud & Pecout. Dans cette section, je réponds aux questions spécifiques posées par les étudiantes de l’Université de Neuchâtel.

Un polygone à partir de points (enveloppes convexes et concaves)

Il peut être utile de dessiner un polygone qui enveloppe un ensemble de points. Définissons d’abord quelques points. Pour cet exemple, j’ai utilisé le site geojson.io pour aisément créer une chaine de caractères GeoJSON contenant les coordonnées de 4 sommets au canton du Tessin. Je vous suggère de créer votre propre collection de points pour rendre cet exercice plus intéressant.

{ "type": "FeatureCollection", "features": [
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "coordinates": [9.007992484256647,46.0299246522398],
        "type": "Point"
      }
    },
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "coordinates": [8.9478934338552,45.977147786636834],
        "type": "Point"
      }
    },
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "coordinates": [8.8329890317913,46.04030106008443],
        "type": "Point"
      }
    },
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "coordinates": [8.950073186262301,45.913272474209464],
        "type": "Point"
      }
    }
]}Code language: JavaScript (javascript)

Pour intégrer cette série de points dans un script R, vous pourriez la stocker dans un fichier de texte brut avec l’extension .geojson et l’importer à l’aide de la fonction st_read(). Comme il y a peu de points, cependant, il vous suffit de copier la chaîne de caractères GeoJSON dans votre code pour la lire. Vous pouvez l’abréger en enlevant toutes les properties, vu que nos points en sont dépourvus:

mespoints <- st_read('{ "type": "FeatureCollection", "features": [
    { "type": "Feature", "geometry": {
        "coordinates": [9.007992484256647,46.0299246522398],
        "type": "Point"
    }},
    { "type": "Feature", "geometry": {
        "coordinates": [8.9478934338552,45.977147786636834],
        "type": "Point"
    }},
    {"type": "Feature", "geometry": {
        "coordinates": [8.8329890317913,46.04030106008443],
        "type": "Point"
    }},
    {"type": "Feature", "geometry": {
        "coordinates": [8.950073186262301,45.913272474209464],
        "type": "Point"
    }}
]}')Code language: R (r)

Voyons ces points sur la carte:

fondsCarte <- st_buffer(mespoints,1000) %>% st_bbox() %>% setNames(c("left", "bottom", "right", "top")) %>% get_stamenmap(maptype = "terrain",zoom=13)
ggmap(fondsCarte) + geom_sf(data=mespoints, inherit.aes = FALSE,color="red",size=3)Code language: PHP (php)

Deux techniques permettent de dessiner un polygone contenant un ensemble de points: l’enveloppe convexe et l’enveloppe concave. Le module sf fournit toutes deux, avec les fonctions st_convex_hull() et st_concave_hull(). Notez bien que leur bon fonctionnement exige que les points pris en cosidération soit réunis au prélable dans un seul objet de type “multipoint” à l’aide de la fonction st_union().

enveloppeConvexe <- st_union(mespoints) %>% st_convex_hull()
ggmap(fondsCarte) + 
	geom_sf(data=mespoints, inherit.aes = FALSE,color="red",size=3) +
	geom_sf(data=enveloppeConvexe, inherit.aes = FALSE,fill="blue",alpha=0.5)Code language: PHP (php)

Pour les enveloppes concaves, il faut préciser le paramètre “ratio” un ratio=1 donnera une envoloppe convexe. Pour toute valeur inférieure à 1, plus le ratio est bas, plus le polyone résultant présentera de concavités.

enveloppeConcave <- st_union(mespoints) %>% st_concave_hull(ratio=0.5)
ggmap(fondsCarte) + 
	geom_sf(data=mespoints, inherit.aes = FALSE,color="red",size=3) +
	geom_sf(data=enveloppeConcave, inherit.aes = FALSE,fill="yellow",alpha=0.5)Code language: PHP (php)

Une carte thermique (heatmap) à partir de points

À partir des points créés dans la section précédente, vous pouvez aussi créer une carte thermique (heatmap) à l’aide de la fonction stat_density2d_filled().

ggmap(fondsCarte) + 
	geom_sf(
		data=mespoints, 
		inherit.aes = FALSE,
		color="red",
		size=3
	) +
	stat_density2d_filled(
		data=st_coordinates(mespoints) %>% as.data.frame, 
		aes(x=X,y=Y),
		inherit.aes = FALSE,
		alpha=0.5
	)Code language: R (r)

Notez que stat_density2d_filled() est une fonction a priori pas conçue pour traiter des géodonnées, mais rappelons-nous que nous travaillons avec de coordonées projetées. Si vous testez

st_coordinates(mespoints) %>% as.data.frame Code language: R (r)

Vous obtenez:

         X        Y
1 9.007992 46.02992
2 8.947893 45.97715
3 8.832989 46.04030
4 8.950073 45.91327Code language: CSS (css)

stat_density2d_filled() peut très bien opérer avec ces coordonnées vous pouvez donc vous en servir pour superposer une heatmap sur une carte créé avec ggmap()

Cette fonction possède en outre un paramètre h vous permettant d’ajuster la “bandwidth“, que vous pouvez comprendre comme rayon de diffusion à partir d’un point dans les deux directions (horizontale et verticale). En utilisant des valeurs basses de h, par exemple, la zone de chaleur sera davantage concentrée autour des points d’intérêt:

ggmap(fondsCarte) + 
	geom_sf(data=mespoints, inherit.aes = FALSE,color="red",size=3) +
	stat_density2d_filled(
		aes(x=X,y=Y),
		data=st_coordinates(mespoints) %>% as.data.frame, 
		inherit.aes = FALSE, 
		alpha=0.5,
		h=c(0.05,0.05) # ajuste la largeur de la zone de chaleur
) Code language: PHP (php)

Autres ressources

• Gimond, M. Intro to GIS and Spatial Analysis. Consulté 31 mai 2023, à l’adresse https://mgimond.github.io/Spatial/index.html