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 de données de type vecteur 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.
• 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 que j’ai affaire à 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 momement et 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

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’entammer 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 reqê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!

bbox <- st_bbox(compoly.tessin %>% st_transform(4326))
q <- opq(bbox, timeout=60) Code language: HTML, XML (xml)

À 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ésulat:

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 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 exmple les parcs nationaux, ou les écoles.