Réaliser une carte avec R (5)

Suite aux parties 1, 2, 3, et 4, du tutoriel de cartographie, nous allons maintenant aborder la cartographie en 3D.

Vous pouvez télécharger le code source complet du notebook R de ce tutoriel ici et le fichier contenant uniquement le code R ici.

Avant toute chose

Autant le dire tout desuite, la cartographie en 3D avec R est plus complexe et moins souple que la cartographie en 2D. Un inconvénient est l’absence de prise en charge de la 3D par ggplot2, qui oblige à faire appel à d’autres modules de visualisation. Toutefois, ponctuellement, la 3D peut avoir son utilité et valoir le coût de son investissement.

Commençons par récupérer les données d’altimétrie pour le Japon. Ensuite nous pouvons découper ce raster au cadre de Kyūshū avec la fonction crop().

library(raster)
jpAlt <- getData("alt", country = "JPN")
xlims <- c(128.5,132)
ylims <- c(31,35)
kyushuAlt <- crop(jpAlt,extent(c(xlims,ylims)))

Vérifions le résultat en visualisant les deux rasters avec la fonction de base plot().

plot(jpAlt)

plot(kyushuAlt)

Visualisation en perspective

On peut alors se servir de la fonction persp pour visualiser nos données de terrain en 3D.

persp(kyushuAlt)

Piètre résultat! Mais pas de panique, il suffit de régler les options de persp() (voir la documentation pour un liste complète) en commençant par expand, qui règle l’échelle de la dimension z (altitude).

persp(kyushuAlt, expand = .025)

Ensuite, il vaut mieux ne pas afficher les bordures de polygones (border) et ajouter à la place une couleur de remplissage (col).

persp(kyushuAlt, expand = .025, border = NA, col = "grey90")

Le relief qui semble avoir disparu réapparait en réglant l’ombrage (shade).

persp(kyushuAlt, expand = .025, border = NA, col = "grey90", shade = .5)

L’orientation et l’inclinaison se règlent avec les options theta et phi.

persp(kyushuAlt, expand = .025, border = NA, col = "grey90", shade = .5, theta = 0, phi = 30)

Le rendu peut être amélioré en augmentant le nombre maximum de pixels (maxpixels).

persp(kyushuAlt, expand = .025, border = NA, col = "grey90", shade = .5, theta = 0, phi = 30, maxpixels = ncell(kyushuAlt))

Enfin, on peut supprimer le cadre autour de la carte.

persp(kyushuAlt, expand = .025, border = NA, col = "grey90", shade = .5, theta = 0, phi = 30, maxpixels = ncell(kyushuAlt), box = FALSE)

Améliorer la visualisation et ajouter des couleurs

Évidemment, on souhaiterait colorier la carte en perspective en fonction de l’altitude. Étrangement, ce n’est pas si simple que ça puisqu’il faut effectuer une transformation préalable.

Il est en effet nécessaire de transformer d’abord notre raster en matrice.

kyushuAltM <- as.matrix(kyushuAlt)

Ensuite, il faut faire appel à la fonction drape.plot du module fields en spécifiant quelles données sont utilisées pour les coordonnées x, y et z.

require(fields)
drape.plot(x = 1:nrow(kyushuAltM), y = 1:ncol(kyushuAltM), z = kyushuAltM, expand = .025, border = NA, col = "grey90", shade = .5, theta = 90, phi = 30, box = FALSE)

Pour obtenir une belle coloration en fonction de l’altitude, il suffit d’utiliser dans l’option col la fonction terrain.colors() et de lui donner un nombre de couleurs.

require(fields)
drape.plot(x = 1:nrow(kyushuAltM), y = 1:ncol(kyushuAltM), z = kyushuAltM, expand = .025, border = NA, col = terrain.colors(128), shade = .5, theta = 90, phi = 30, box = FALSE)

Si l’on désire lisser le relief, on peut modifier le raster au préalable.

kyushuAlt2 <- disaggregate(kyushuAlt, fact = 2, dissolve = TRUE)
kyushuAlt2 <- focal(kyushuAlt2, w = matrix(1, 5, 5), mean)
kyushuAltM2 <- as.matrix(kyushuAlt2)
drape.plot(x = 1:nrow(kyushuAltM2), y = 1:ncol(kyushuAltM2), z = kyushuAltM2, expand = .025, border = NA, col = terrain.colors(128), shade = .5, theta = 90, phi = 30, box = FALSE)

Perspective et polygones

Malheureusement je n’ai pas trouvé de méthode simple et efficace pour projeter des polygones sur une surface en perspective. La méthode ci-dessous est donc loin d’être utilisable pour des besoins réels.

Il faut en premier lieu récupérer nos polygones et on extrait la préfecture de Kagoshima. Ici, comme nous n’utilisons pas ggpolt2, pas besoin d’utiliser fortify().

require(rgdal)
jpn1 <- readOGR(dsn = path.expand("~/repos/cartes/jpn/adm/jpn"), layer = "JPN_adm1")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/tpellard/repos/cartes/jpn/adm/jpn", layer: "JPN_adm1"
## with 47 features
## It has 9 fields
kag <- subset(jpn1, NAME_1 == "Kagoshima")

Ensuite on crée un masque du raster, c’est à dire que nous allons effacer les données hors des polygones de Kagoshima. On peut visualiser rapidement le résultat avec plot().

kagAlt <- mask(kyushuAlt, kag)
plot(kagAlt)

On peut maintenant superposer nos deux raster et colorier Kagoshima en rouge, en prennant soin d’utiliser les mêmes valeurs pour les options.

persp(kyushuAlt, expand = .025, border = NA, col = "grey90", shade = .5, theta = 0, phi = 30, maxpixels = ncell(kyushuAlt), box = FALSE)
par(new=TRUE)
persp(kagAlt, expand = .025, border = NA, col = "red", shade = .5, theta = 0, phi = 30, maxpixels = ncell(kagAlt), box = FALSE)

On peut superposer également Kagoshima à une carte coloriée selon l’altitude, même si cela présente peut d’intérêt.

drape.plot(x = 1:nrow(kyushuAltM), y = 1:ncol(kyushuAltM), z = kyushuAltM, expand = .025, border = NA, col = terrain.colors(128), shade = .5, theta = 90, phi = 30, box = FALSE)
par(new=TRUE)
persp(kagAlt, expand = .025, border = NA, col = "red", shade = .5, theta = 0, phi = 30, maxpixels = ncell(kagAlt), box = FALSE)

À l’inverse, on peut afficher uniquement l’altitude en couleur de Kagoshima.

kagAltM <- as.matrix(kagAlt)
persp(kyushuAlt, expand = .025, border = NA, col = "grey90", shade = .5, theta = 0, phi = 30, maxpixels = ncell(kyushuAlt), box = FALSE)
par(new=TRUE)
drape.plot(x = 1:nrow(kagAltM), y = 1:ncol(kagAltM), z = kagAltM, expand = .025, border = NA, col = terrain.colors(128), shade = .5, theta = 90, phi = 30, box = FALSE)

Malheureusement, je ne connais pas de méthode pour projeter automatiquement un grand nombre de polygones.

Aller encore plus loin

Le module rasterVis permet de manipuler des raster en 3D de manière interactive et créer des animations.

require(rgl)
require(rasterVis)
plot3D(kyushuAlt, zfac = .25)
snapshot3d("3d.png")

À suivre dans la sixième partie.

Ce billet peut être cité de la façon suivante: Thomas Pellard, « Réaliser une carte avec R (5) », Cipanglossia, 26 juin 2017, https://cipanglo.hypotheses.org/466.

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *