On cherche ici à produire carte du Mont-Ventoux (Vaucluse, France) et de ses alentours, ainsi qu’une vue en 3D des structures géologiques à l’aide d’un modèle d’élévation. Le tout avec R et des données libres.

Les données utilisées sont celles de la BD ALTI, il s’agit d’un modèle numérique de terrain (MNT) à la résolution de 25 m produit par l’IGN.

Pour y parvenir, deux packages R sont nécessaires : terra (Hijmans 2025) pour le traitement des données raster et rayshader (Morgan-Wall 2024) qui permet de visualiser les données d’élévation.

## Installer les packages nécessaires
pkg <- c("khroma", "terra", "rayshader")
pkg_missing <- !(pkg %in% installed.packages())
install.packages(pkg[pkg_missing])

La zone d’intérêt est un carré de 50 km de côté, grossièrement centrée sur le Mont-Ventoux :

## Définir l'emprise du découpage (xmin, xmax, ymin, ymax)
## Système de coordonnées projetées RGF93-Lambert 93
extent_ventoux <- terra::ext(861000, 911000, 6310000, 6360000)

Préparation des données

Les données d’élévation nécessaires sont issues de la BDALTI et peuvent être téléchargées librement depuis le catalogue de l’IGN. Les données sont produites avec une résolution spatiale de 25 m et utilisent le système de coordonnées projeté RGF93 (EPSG:2154). Neuf dalles sont nécessaires :

BDALTIV2_25M_FXX_0850_6325_MNT_LAMB93_IGN69
BDALTIV2_25M_FXX_0850_6350_MNT_LAMB93_IGN69
BDALTIV2_25M_FXX_0850_6375_MNT_LAMB93_IGN69
BDALTIV2_25M_FXX_0875_6325_MNT_LAMB93_IGN69
BDALTIV2_25M_FXX_0875_6350_MNT_LAMB93_IGN69
BDALTIV2_25M_FXX_0875_6375_MNT_LAMB93_IGN69
BDALTIV2_25M_FXX_0900_6325_MNT_LAMB93_IGN69
BDALTIV2_25M_FXX_0900_6350_MNT_LAMB93_IGN69
BDALTIV2_25M_FXX_0900_6375_MNT_LAMB93_IGN69

La première étape consiste à regrouper les neuf dalles, puis à découper les données selon l’emprise souhaitée.

## Lister les fichiers
## (les fichiers .asc sont placés dans le répertoire de travail)
BDALTI_files <- list.files(path = ".", pattern = ".asc", full.names = TRUE)

## Créer un raster virtuel temporaire
vrt_tmp <- paste0(tempfile(), ".vrt")
BDALTI_vrt <- terra::vrt(BDALTI_files, filename = vrt_tmp)

## Préciser le système de coordonnées
terra::crs(BDALTI_vrt) <- "epsg:2154"

## Découper les données
BDALTI_crop <- raster::crop(x = BDALTI_vrt, y = extent_ventoux)

## Enregistrer le raster découpé pour un usage futur
terra::writeRaster(BDALTI_crop, filename = "ventoux.tif")

Visualisation

Élévation

Pour visualiser les données d’élévation, il peut être préférable de les discrétiser : les regrouper selon différents intervalles. Pour réaliser ce découpage, il est utile de commencer par explorer la distribution des données à l’aide d’un histogramme (l’argument breaks de la fonction hist() permet de spécifier le nombre de classes souhaitées) :

## Distribution des valeurs d'élévation
elevation <- terra::hist(
  x = BDALTI_crop,
  breaks = 10, # 10 classes de même amplitude
  main = NULL,
  xlab = "Élévation (m)"
)

## Récupérer les bornes des classes de l'histogramme
breaks <- elevation$breaks

Figure 1: Histogramme des valeurs d’élévation.

Les classes obtenues à partir de l’histogramme sont de même amplitude. Une autre approche peut être d’utiliser les déciles (ou les quartiles, etc.) :

## Calcul des déciles
deciles <- terra::global(BDALTI_crop, fun = quantile, probs = seq(0, 1, 0.1))

## Arrondir à la dizaine pour faciliter la lecture
deciles <- round(deciles, digits = -1)

On peut comparer les deux approches :

terra::plot(
  x = BDALTI_crop, 
  breaks = breaks, # Classes de l'histogramme
  col = khroma::color("bamako")(10),
  mar = c(3, 3, 1, 6) + 0.1
)

Figure 2: Modèle numérique de terrain : classes de l’histogramme.

terra::plot(
  x = BDALTI_crop, 
  breaks = deciles, # Déciles
  col = khroma::color("bamako")(10),
  mar = c(3, 3, 1, 6) + 0.1
)

Figure 3: Modèle numérique de terrain : déciles.

Ombrages

Au lieu d’utiliser l’élévation absolue du modèle de terrain, les fonctions sphere_shade() et lamb_shade() du package rayshader permettent de colorer la surface en fonction de la direction et de la déclivité des pentes et de la position de la source de lumière (hillshading).

Ces méthodes peuvent être combinées avec les fonctions texture_shade() qui permet de calculer les ombrage selon la méthode de Brown (2014) et ambient_shade() qui permet de calculer les ombres liées à la diffusion de la lumière par l’atmosphère¹.

## Convertir les données raster en matrice
ventoux <- rayshader::raster_to_matrix(BDALTI_crop)

## Calcul des ombrages
## Sphere shading
sph <- rayshader::sphere_shade(ventoux, texture = "desert", 
                               colorintensity = 5, zscale = 3)

## Hillshading
lam <- rayshader::lamb_shade(ventoux, sunaltitude = 30, zscale = 12)

## Occlusion atmosphérique
amb <- rayshader::ambient_shade(ventoux)

## Texture Shading
tex <- rayshader::texture_shade(ventoux, detail = 2/3, 
                                contrast = 5, brightness = 4)

On peut visualiser les différentes couches d’ombrage et ajuster au besoin les paramètres zscale (facteur d’échelle entre coordonnées horizontales et verticale), sunaltitude et sunangle (position du soleil).

## Hillshading
rayshader::plot_map(lam)

## Occlusion atmosphérique
rayshader::plot_map(amb)

## Texture shading
rayshader::plot_map(tex)

## Sphere shading
rayshader::plot_map(sph)

Si les ombrages réalisés sont satisfaisant, on peut les combiner avec les données d’élévation pour visualiser le relief :

## Elévation et ombrages
relief <- ventoux |>
  rayshader::height_shade() |>
  rayshader::add_overlay(sph, alphalayer = 0.5) |>
  rayshader::add_shadow(lam, max_darken = 0.5) |>
  rayshader::add_shadow(amb, max_darken = 0.5) |> 
  rayshader::add_shadow(tex, max_darken = 0.2)

Le résultat de la fonction height_shade() est un array de valeurs RGB. Ce tableau à trois dimensions a les mêmes dimensions et la même résolution que le MNT d’origine, il peut aisément être géoréférencé pour être complété (barre d’échelle, flèche vers le nord…) :

## Convertir en raster spatial
ventoux_relief <- terra::rast(relief)

## Préciser que les couches correspondent aux valeurs de rouge, vert et bleu
terra::RGB(ventoux_relief) <- c(1, 2, 3)

## Préciser l'étendue géographique
terra::ext(ventoux_relief) <- extent_ventoux

## Préciser le système de coordonnées
terra::crs(ventoux_relief) <- "epsg:2154"

## Visualiser
terra::plotRGB(ventoux_relief, scale = 1)

## Ajouter une barre d'échelle
terra::sbar(d = 10000, label = c(0, "km", 10), type = "bar", col = "white")

## Ajouter une flèche vers le nord
terra::north(type = 2, col = "white")

Figure 8: Relief de la région du Mont-Ventoux.

Le résultat obtenu avec rayshader peut également être exportée au format GeoTiff pour d’autres usages (dans un logiciel de SIG comme QGIS par exemple) :

## Exporter au format GeoTiff
raster::writeRaster(ventoux_relief, filename = "ventoux_relief.tif")

Finalement, on peut générer un modèle 3D (l’échelle verticale est exagérée d’un facteur 2) :

## Visualiser le modèle 3D
rayshader::plot_3d(
  hillshade = relief,
  heightmap = ventoux,
  zscale = 12, # Exagération verticale
  theta = -45,
  phi = 30,
  fov = 0,
  zoom = 0.6,
  windowsize = c(2048, 1152)
)

## Enregistrer la vue
rayshader::render_snapshot()

Retour au sommet

Session

R version 4.5.1 (2025-06-13)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    
 [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
 [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/Paris
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] rgl_1.3.24         jsonlite_2.0.0     compiler_4.5.1     rayshader_0.37.3  
 [5] crayon_1.5.3       Rcpp_1.1.0         magick_2.9.0       dichromat_2.0-0.1 
 [9] parallel_4.5.1     png_0.1-8          progress_1.2.3     scales_1.4.0      
[13] yaml_2.3.10        fastmap_1.2.0      here_1.0.2         R6_2.6.1          
[17] khroma_1.16.0      knitr_1.50         iterators_1.0.14   htmlwidgets_1.6.4 
[21] rprojroot_2.1.1    RColorBrewer_1.1-3 rlang_1.1.6        Rttf2pt1_1.3.14   
[25] terra_1.8-70       xfun_0.53          doParallel_1.0.17  cli_3.6.5         
[29] magrittr_2.0.4     digest_0.6.37      foreach_1.5.2      grid_4.5.1        
[33] rstudioapi_0.17.1  base64enc_0.1-3    hms_1.1.3          lifecycle_1.0.4   
[37] prettyunits_1.2.0  vctrs_0.6.5        extrafont_0.20     rayimage_0.15.1   
[41] evaluate_1.0.5     glue_1.8.0         farver_2.1.2       extrafontdb_1.1   
[45] codetools_0.2-20   rmarkdown_2.30     tools_4.5.1        pkgconfig_2.0.3   
[49] htmltools_0.5.8.1

Les références

Brown, Leland H. 2014. « TEXTURE SHADING : A NEW TECHNIQUE FOR DEPICTING TERRAIN RELIEF ». In.

Hijmans, Robert J. 2025. terra: Spatial Data Analysis. https://doi.org/10.32614/CRAN.package.terra.

Morgan-Wall, Tyler. 2024. rayshader: Create Maps and Visualize Data in 2D and 3D. https://doi.org/10.32614/CRAN.package.rayshader.

Notes de bas de page

Voir https://github.com/tylermorganwall/MusaMasterclass pour plus de détails.↩︎

Réutilisation

CC BY 4.0

Citation

BibTeX

@online{frerebeau2024,
  author = {Frerebeau, N.},
  title = {Géographie physique du Mont-Ventoux},
  date = {2024-08-06},
  url = {https://carnets.archeosciences-bordeaux.fr/outils/r-ventoux-elevation/},
  langid = {fr-FR}
}

Veuillez citer ce travail comme suit :

Frerebeau, N. 2024. “Géographie physique du Mont-Ventoux.” August 6, 2024. https://carnets.archeosciences-bordeaux.fr/outils/r-ventoux-elevation/.