Réaliser une carte avec R (4)

Dans les parties 1, 2 et 3, nous avons vu comment réaliser des cartes descriptive simples et nous avons abordé la cartographie de données quantitatives. Nous allons maintenant améliorer nos cartes en leur donnant un peu de relief.

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.

Récupérer et afficher des données de relief

Pour éviter d’avoir à manipuler des fichiers de trop grande taille, nous allons nous concentrer sur la région de l’ile de Kyūshū. Commençons par récupérer les données de bathymétrie et d’altimétrie de cette région avec la fonction getNOAA.bathy() du module marmap. L’option res indique le niveau de résolution des données (attention le fichier devient vite très lourd). L’option keep permet d’enregistrer le fichier automatiquement sur le disque dur pour le réutiliser.

xlims <- c(128.5,132)
ylims <- c(31,35)
require(marmap)
require(ggplot2)
kyushuBath <- getNOAA.bathy(min(xlims), max(xlims), min(ylims), max(ylims), res = 1, keep = TRUE)
## Querying NOAA database ...
## This may take seconds to minutes, depending on grid size
## Building bathy matrix ...
# kyushuBath <- read.bathy("marmap_coord_128.5;31;132;35_res_1.csv", header = T) # Pour lire le fichier enregistré

On peut alors cartographier ces données. Commençons par tracer les données sous forme de ligne de contour, après, comme d’habitude, avoir utilisé fortify().

kyushuBathF <- fortify(kyushuBath)
ggplot() + 
  geom_contour(data = kyushuBathF, aes(x = x, y = y, z = z)) + 
  coord_quickmap() + 
  theme_classic()

On n’y voit pas grand chose. Rajoutons les polygones de la terre ferme en gris.

require(maptools)
require(rgdal)
jpn0 <- readOGR(dsn = path.expand("~/repos/cartes/jpn/adm/jpn"), layer = "JPN_adm0")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/tpellard/repos/cartes/jpn/adm/jpn", layer: "JPN_adm0"
## with 1 features
## It has 24 fields
jpn0F <- fortify(jpn0)
ggplot() + 
  geom_polygon(data = jpn0F, aes(x = long, y = lat, group = group), fill = "grey80", color = "NA") +
  geom_contour(data = kyushuBathF, aes(x = x, y = y, z = z)) + 
  coord_quickmap(xlim = xlims, ylim = ylims) +
  theme_classic()

Afficher les données de relief sous forme de raster

Reconnaissons que les lignes de contour ne sont pas le mode de visualisation de l’altitude le plus efficace. Il vaut mieux cartographier nos données comme une couche de type image raster.

ggplot() + 
  geom_raster(data = kyushuBathF, aes(x = x, y = y, fill = z)) + 
  coord_quickmap(xlim = xlims, ylim = ylims) +
  theme_classic()

La visualisation devient beaucoup plus claire si l’on change la palette de couleurs.

ggplot() + 
  geom_raster(data = kyushuBathF, aes(x = x, y = y, fill = z)) + 
  scale_fill_etopo() +
  coord_quickmap(xlim = xlims, ylim = ylims) +
  theme_classic()

L’option interpolate permet de lisser l’image afin de limiter la pixelisation.

ggplot() + 
  geom_raster(data = kyushuBathF, aes(x = x, y = y, fill = z), interpolate = TRUE) + 
  scale_fill_etopo() +
  coord_quickmap(xlim = xlims, ylim = ylims) +
  theme_classic()

Créer un relief ombré

Si l’on peut visualiser le terrain avec les couleurs, notre carte manque néanmoins encore de relief. Il suffit de créer ce relief avec la suite de manipulations suivantes.

require(raster)
kyushuRast <- marmap::as.raster(kyushuBath)
kyushuSlope <- terrain(kyushuRast, opt="slope")
kyushuAspect <- terrain(kyushuRast, opt="aspect")
kyushuShade <- hillShade(kyushuSlope, kyushuAspect)
kyushuShadeSpdf <- as(kyushuShade, "SpatialPixelsDataFrame")
kyushuShadeDf <- as.data.frame(kyushuShadeSpdf)
ggplot() + 
  geom_raster(data = kyushuShadeDf, aes(x = x, y = y, fill = layer), interpolate = TRUE) +
  coord_quickmap(xlim = xlims, ylim = ylims) +
  theme_classic()

Je trouve la visualisation meilleure en niveaux de gris.

ggplot() + 
  geom_raster(data = kyushuShadeDf, aes(x = x, y = y, fill = layer), interpolate = TRUE) + 
  scale_fill_gradient(low = "black", high = "white") +
  coord_quickmap(xlim = xlims, ylim = ylims) +
  theme_classic()

On peut jouer sur le paramètres de la fonction hillshadepour obtenir une illumination du relief différente.

kyushuShade <- hillShade(kyushuSlope, kyushuAspect, angle = 80, direction = 0)
kyushuShadeSpdf <- as(kyushuShade, "SpatialPixelsDataFrame")
kyushuShadeDf <- as.data.frame(kyushuShadeSpdf)
ggplot() + 
  geom_raster(data = kyushuShadeDf, aes(x = x, y = y, fill = layer), interpolate = TRUE) + 
  scale_fill_gradient(low = "black", high = "white") +
  coord_quickmap(xlim = xlims, ylim = ylims) +
  theme_classic()

On peut alors superposer par transparence (alpha) la couche d’ombres de relief et la couche de couleurs.

ggplot() +
  geom_raster(data = kyushuShadeDf, aes(x = x, y = y, alpha = 1 - layer), fill = "gray10", interpolate = TRUE) +
  geom_raster(data = kyushuBathF, aes(x = x, y = y, fill = z), alpha = .75, interpolate = TRUE) +
  scale_fill_etopo() +
  scale_alpha_continuous(guide = "none") +
  coord_quickmap(xlim = xlims, ylim = ylims) +
  theme_classic()

On peut jouer sur la transparence pour accentuer le relief.

ggplot() +
  geom_raster(data = kyushuShadeDf, aes(x = x, y = y, alpha = 1 - layer), fill = "gray10", interpolate = TRUE) +
  geom_raster(data = kyushuBathF, aes(x = x, y = y, fill = z), alpha = .45, interpolate = TRUE) +
  scale_fill_etopo() +
  scale_alpha_continuous(guide = "none") +
  coord_quickmap(xlim = xlims, ylim = ylims) +
  theme_classic()

On peut bien sûr rajouter des décorations dans des couches supplémentaires.

jpn1 <- readOGR(dsn = path.expand("~/repos/cartes/jpn/adm/jpn"), layer = "JPN_adm1")
jpn1F <- fortify(jpn1, region = "NAME_1")
ggplot() +
  geom_raster(data = kyushuShadeDf, aes(x = x, y = y, alpha = 1 - layer), fill = "gray10", interpolate = TRUE) +
  geom_raster(data = kyushuBathF, aes(x = x, y = y, fill = z), alpha = .45, interpolate = TRUE) +
  geom_polygon(data = subset(jpn1F, id == "Kagoshima"), aes(x = long, y = lat, group = group), fill = "red", alpha = .5) +
  scale_fill_etopo() +
  scale_alpha_continuous(guide = "none") +
  coord_quickmap(xlim = xlims, ylim = ylims) +
  theme_classic()

Aller plus loin

On peut aussi préférer séparer la bathymétrie de l’altimétrie, et cartographier seulement cette dernière en relief ombré, en affichant la bathymétrie sous forme de contour en arrière plan uniquement.

jpAlt <- getData("alt", country = "JPN")
jpAltSpdf <- as(jpAlt, "SpatialPixelsDataFrame")
jpAltDf <- as.data.frame(jpAltSpdf)
jpSlope <- terrain(jpAlt, opt="slope")
jpAspect <- terrain(jpAlt, opt="aspect")
jpShade <- hillShade(jpSlope, jpAspect, angle = 80, direction = 0)
jpShadeSpdf <- as(jpShade, "SpatialPixelsDataFrame")
jpShadeDf <- as.data.frame(jpShadeSpdf)
ggplot() + 
  geom_raster(data = kyushuBathF, aes(x = x, y = y, fill = z), interpolate = TRUE) +
  geom_raster(data=jpShadeDf, aes(x = x, y = y, alpha = 1 - layer), fill = "gray10", interpolate = TRUE) +
  geom_raster(data = jpAltDf, aes(x = x, y = y, fill = JPN_msk_alt), alpha = .45, interpolate = TRUE) +
  scale_fill_etopo() +
  scale_alpha_continuous(guide = "none") +
  coord_quickmap(xlim = xlims, ylim = ylims) + 
  theme_classic()

On peut préférer les fichiers raster colorisés et mis en relief de Natural Earth.

# http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/raster/HYP_HR_SR_OB_DR.zip (385.56 MB)
neHy <- stack("~/Téléchargements/HYP_HR_SR_OB_DR/HYP_HR_SR_OB_DR.tif") # 700 MB une fois décompressée
kyushuHyp <- crop(neHy,extent(c(xlims,ylims)))
kyushuHypSpdf <- as(kyushuHyp, "SpatialPixelsDataFrame")
kyushuHypDf <- as.data.frame(kyushuHypSpdf)
kyushuHypTb <- data.frame(xyFromCell(kyushuHyp, 1:ncell(kyushuHyp)), getValues(kyushuHyp/255))
kyushuHypTb$rgb <- rgb(red = kyushuHypTb[,3], green = kyushuHypTb[,4], blue = kyushuHypTb[,5], alpha = 1)
ggplot() + 
  geom_raster(data = kyushuHypTb, aes(x = x, y = y), fill = kyushuHypTb$rgb, interpolate = TRUE) +
  coord_equal() + 
  theme_classic()

L’inconvénient de ce genre de carte en relief est la taille importante des fichiers et la difficulté à obtenir des données de qualité pour des régions de petite échelle.

À suivre dans la cinquième partie.

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

Laisser un commentaire

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