6.3 Les couches “raster”

Jusqu’ici, nous avons travaillé exclusivement avec des données vectorielles, c’est-à-dire, des données représentées par des points, des traits, des polygones dont les sommets sont définis par des positions géographiques. Tous ces objets sont matérialisés en R sous forme sf et nous avons vu que geom_sf() permet de manipuler tout cela pratiquement comme d’habitude à l’aide de chart() et le moteur graphique {ggplot2}. Toutefois, il y a une autre catégorie de données que sont les rasters. Il s’agit ici de données (numériques la plupart du temps) présentées sur une grille rectangulaire. Le principe est le même que pour les images numériques. La scène est découpée en petits carrés appelés pixels et une valeur est associée à chaque pixel.

Une image numérique non vectorielle, tout comme une couche d’information géographique, est un “raster” constitué d’une grille de petits carrés (les pixels) auxquels on attribue une valeur (une couleur pour l’image). Ainsi l’image d’un chat à gauche est encodée à droite sous forme d’une grille de pixels bien visible à la basse résolution choisie ici.

Les couches “raster” sont utilisées pour représenter, par exemple, des modèles terrain reprenant l’élévation du sol, des images satellitaires, ou toute autre caractéristique mesurée ou calculée sur une grille rectangulaire dense. Nous allons illustrer cela en utilisant un raster de type modèle terrain et en calculant un autre (pluviométrie sur un pays) dans la suite de ce module.

6.3.1 Modèle terrain du Maroc

Le Maroc présente un relief contrasté avec, au nord-ouest, la chaîne montagneuse de l’Atlas, et au sud, le désert du Sahara. Le modèle terrain avec une grille carrée d’un kilomètre environ de côté est disponible dans le package {aurelhy}. Les données dans {aurelhy} sont dans un format propriétaire et aucune fonction de conversion en raster ou sf n’y est proposée. Nous allons donc utiliser nos propres fonctions de conversion ici (elles seront intégrées dans la prochaine version du package {aurelhy}). Il n’est pas indispensable de comprendre le contenu de ce code en première lecture.

# Transformation de geomat, geotm & geomask en raster
st_as_raster <- function(x, ...)
  UseMethod("st_as_raster")
st_as_raster.geomat <- function(x, crs = "+proj=longlat +datum=WGS84", ...) {
  x_coords <- aurelhy::coords(x)
  raster::raster(t(unclass(x))[ncol(x):1, ],
  xmn = x_coords["x1"], xmx = x_coords["x2"],
  ymn = x_coords["y1"], ymx = x_coords["y2"],
  crs = crs)
}

# Transformation de geoshapes en sf, sfc_MULTIPOLYGONS
library(sf)

st_as_sf.geoshapes <- function(x, name = "P1", crs = "+proj=longlat +datum=WGS84", ...) {
  x2 <- dplyr::bind_rows(as.list(unclass(x)), .id = "id")
  points <- st_as_sf(x2, coords = c("x", "y"))
  polygons <- st_sf(aggregate(points$geometry, list(points$id),
    function(g) st_cast(st_combine(g), "POLYGON")))
  st_crs(polygons) <- crs
  st_as_sf(tibble::tibble(name = name,
    geometry = st_combine(polygons)))
}

# Transformation de geopoints en sf, sfc_POINTS
st_as_sf.geopoints <- function(x, crs = "+proj=longlat +datum=WGS84", ...) {
  class(x) <- "data.frame"
  x_sf <- st_as_sf(x, coords = c("x", "y"))
  st_crs(x_sf) <- crs
  x_sf
}

À présent, nous allons récupérer les données raster du modèle terrain à partir de morocco et les données vectorielles des frontières du Maroc à partir de mbord, et nous les convertissons en raster et sf, respectivement.

library(raster)

# Le chargement a nécessité le package : sp

# 
# Attachement du package : 'raster'

# L'objet suivant est masqué depuis 'package:ggsn':
# 
#     scalebar

# L'objet suivant est masqué depuis 'package:dplyr':
# 
#     select

# L'objet suivant est masqué depuis 'package:MASS':
# 
#     select

morocco <- read("morocco", package = "aurelhy")
m_tm <- st_as_raster(morocco)
mbord <- read("mbord", package = "aurelhy")
mbord[[1]] <- mbord[[1]][-(1:28), ] # Données inutiles au début à éliminer
m_bord <- st_as_sf(mbord, name = "Morocco")

Vérifions tout d’abord que les frontières (vectorielles) du Maroc ont été converties correctement :

chart(data = m_bord) +
  geom_sf()

Maintenant, notre objet m_tm.

class(m_tm)

# [1] "RasterLayer"
# attr(,"package")
# [1] "raster"

length(m_tm)

# [1] 3752244

res(m_tm)

# [1] 0.008333334 0.008333334

m_tm

# class      : RasterLayer 
# dimensions : 1878, 1998, 3752244  (nrow, ncol, ncell)
# resolution : 0.008333334, 0.008333334  (x, y)
# extent     : -17.34825, -0.6982474, 20.54881, 36.19881  (xmin, xmax, ymin, ymax)
# crs        : +proj=longlat +datum=WGS84 
# source     : memory
# names      : layer 
# values     : 1, 4057  (min, max)

Notre objet est précisément un RasterLayer (une couche raster). Il contient plus de 3,7 millions de pixels. La résolution d’un pixel (carré) est d’environ 0.008° en latitude comme en longitude. Cela représente environ 924km à cet endroit de la terre. Les valeurs (élévations) vont de 1m à 4057m.

La plupart des rasters que vous importerez proviendront de logiciels SIG et seront enregistrés soit en une forme particulière d’image TIFF appelée GeoTIFF avec extension .tif, soit au format NetCDF avec extenstion .nc. Il suffit d’utiliser raster("filename") pour lire ces données et les inclure dans un objet RasterLayer en une seule étape. L’objet supporte également des rasters tellement volumineux qu’ils ne tiennent pas entièrement en mémoire !

Pour exporter vos rasters, vous utiliserez la fonction writeRaster(). Elle supporte de nombreux formats différents dont les deux cités ci-dessus, voir ?writeFormats.

Un raster de plusieurs millions de points n’est pas inhabituel. Cependant, pour représenter l’élévation du sol sur tout le territoire marocain sur une carte à basse ou moyenne résolution, nous n’avons pas forcément besoin d’autant d’information. Nous pouvons utiliser aggregate() pour réduire cette résolution. Nous allons la diviser par cinq. Mais comme cela se fait dans les deux dimensions, l’objet sera 5 * 5 = 25 fois plus petit !

m_tm2 <- aggregate(m_tm, fact = 5)
length(m_tm2)

# [1] 150400

res(m_tm2)

# [1] 0.04166667 0.04166667

Notre objet m_tm2ne contient plus qu’environ 150.000 pixels. Il existe de nombreuses façons de représenter visuellement cet objet sur une carte dans R, mais si nous suivons l’approche choisie jusqu’ici d’utiliser chart() et geom_sf() avec des objets vectoriels sf, nous devons le convertir en un objet compatible pour pouvoir le rajouter comme couche supplémentaire. C’est possible en passant par un autre objet appelé stars avec st_as_stars() et en utilisant geom_stars(, sf = TRUE)… mais déjà ici, ce n’est plus très simple à manipuler. Voici comment cela serait fait (mais le graphique n’est pas généré… explications juste en dessous)…

library(stars)

chart(data = m_bord) +
  # Conversion en stars et ajoute comme couche 'sf-compatible'
  geom_stars(data = st_as_stars(m_tm2), sf = TRUE) +
  # Tracé des frontières du pays
  geom_sf()

Si vous exécutez ce code, vous vous rendrez compte qu’il faut un temps énorme (plusieurs dizaines de minutes sur un PC puissant et doté de beaucoup de mémoire vive) pour réaliser ce graphique. En fait, notre raster est converti en un objet vectoriel constitué d’autant de petits polygones colorés qu’il y a de pixels dans le raster. Cela prend énormément de temps à calculer, et ensuite à tracer sur le graphique.

Clairement, la combinaison de rasters avec sf sur des graphiques {ggplot2}, si elle est techniquement possible, doit être strictement limitée à des rasters de très basse résolution. Dans la section suivante, nous découvrirons un moteur de cartographie dans R qui ne présente pas ces limitations et qui reste assez semblable à {ggplot2} du point de vue de son interface : {tmap}.