On cherche ici à produire carte géologique 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 du BRGM pour la géologie et de l’IGN pour l’élévation.

Pour y parvenir, cinq packages R sont nécessaires en plus de ggplot2 (Wickham 2016) :

sf (Pebesma 2018; Pebesma et Bivand 2023) pour le traitement des données vectorielles,
terra (Hijmans 2025) pour le traitement des données raster,
rayshader (Morgan-Wall 2024) et rayimage (Morgan-Wall 2025) pour la visualisation en 2D et 3D des données d’élévation,
rgl (Murdoch et Adler 2025) pour le rendu des scènes 3D.

## Installer les packages nécessaires
pkg <- c("sf", "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
box <- c(xmin = 861000, xmax = 911000, ymin = 6310000, ymax = 6360000)
extent_ventoux <- sf::st_bbox(box)

Préparation des données

Les données d’élévation nécessaires sont issues de la BD ALTI, un modèle numérique de terrain (MNT) à la résolution de 25 m produit par l’IGN. Les données ont été préparées lors de la première étape.

## Lire le MNT
BDALTI_crop <- terra::rast(here::here("posts/r-ventoux-elevation/ventoux.tif"))

Les cartes géologiques vectorisées harmonisées à l’échelle 1:50000 sont distribuées par le BRGM sous licence ouverte etalab. Les cartes au format shapefile peuvent être téléchargées par département sur le site du BRGM ou sur data.gouv.fr.

La zone d’intérêt est à cheval sur trois départements : les Alpes-de-Haute-Provence, la Drôme et le Vaucluse :

## Téléchargement depuis le site du BRGM
BRGM_url <- "http://infoterre.brgm.fr/telechargements/BDCharm50/"

## Alpes-de-Haute-Provence (04), Drôme (26) et Vaucluse (84)
dpt <- c("GEO050K_HARM_004.zip", "GEO050K_HARM_026.zip", "GEO050K_HARM_084.zip")
for (i in dpt) {
  ## Téléchargement des données
  utils::download.file(paste0(BRGM_url, i), destfile = i, 
                       method = "wget", quiet = TRUE)
  ## Décompression des fichiers
  utils::unzip(i, exdir = tools::file_path_sans_ext(i))
}

Le BRGM met également à disposition des données supplémentaires relatives à la symbologie, notamment les polices de caractères utilisées pour les symboles ponctuels (voir Janjou (2004) pour la documentation).

## Téléchargement des données
utils::download.file(paste0(BRGM_url, "Outils.zip"), destfile = "Outils.zip")

## Décompression des fichiers
utils::unzip("Outils.zip", exdir = "./")

## Installation des polices de caractères (sous Ubuntu)
system("cp -a OUTILS/Point_Symbols/. ~/.fonts")

Pour chaque département, on importe trois couches et on les fusionne par département :

S_FGEOL : formations géologiques ;
L_STRUCT : objets linéaires structuraux ;
P_STRUCT : éléments ponctuels structuraux.

## Lecture des couches
S_FGEOL_04 <- sf::st_read("GEO050K_HARM_004/GEO050K_HARM_004_S_FGEOL_2154.shp")
S_FGEOL_26 <- sf::st_read("GEO050K_HARM_026/GEO050K_HARM_026_S_FGEOL_2154.shp")
S_FGEOL_84 <- sf::st_read("GEO050K_HARM_084/GEO050K_HARM_084_S_FGEOL_2154.shp")

L_STRUCT_04 <- sf::st_read("GEO050K_HARM_004/GEO050K_HARM_004_L_STRUCT_2154.shp")
L_STRUCT_26 <- sf::st_read("GEO050K_HARM_026/GEO050K_HARM_026_L_STRUCT_2154.shp")
L_STRUCT_84 <- sf::st_read("GEO050K_HARM_084/GEO050K_HARM_084_L_STRUCT_2154.shp")

P_STRUCT_04 <- sf::st_read("GEO050K_HARM_004/GEO050K_HARM_004_P_STRUCT_2154.shp")
P_STRUCT_26 <- sf::st_read("GEO050K_HARM_026/GEO050K_HARM_026_P_STRUCT_2154.shp")
P_STRUCT_84 <- sf::st_read("GEO050K_HARM_084/GEO050K_HARM_084_P_STRUCT_2154.shp")

## Union des couches
L_STRUCT <- rbind(L_STRUCT_04, L_STRUCT_26, L_STRUCT_84)
P_STRUCT <- rbind(P_STRUCT_04, P_STRUCT_26, P_STRUCT_84)
S_FGEOL <- rbind(S_FGEOL_04, S_FGEOL_26, S_FGEOL_84)

## Unions des polygones par attributs
S_FGEOL <- sf::st_union(S_FGEOL, by_feature = TRUE)

## Nettoyage des polygones
S_FGEOL <- sf::st_buffer(S_FGEOL, dist = 0)

Toutes les couches sont projetées avec le même système de coordonnées (RGF93-Lambert 93), il reste à les découper selon la même emprise.

## Découper les données
S_FGEOL_crop <- sf::st_crop(x = S_FGEOL, y = extent_ventoux)
L_STRUCT_crop <- sf::st_crop(x = L_STRUCT, y = extent_ventoux)
P_STRUCT_crop <- sf::st_crop(x = P_STRUCT, y = extent_ventoux)

La table attibutaire de la couche S_FGEOL contient les couleurs associées à chaque formation géologique. Ces dernières sont données en coordonnées CMYB dans les colonnes C_FOND, M_FOND, Y_FOND et B_FOND. On ajoute ainsi à la table attributaire une colonne HEX_FOND contenant le code couleur hexadécimal correspondant aux coordonnées CMYB.

## Fonction de conversion CMYB vers RGB
CMYK2RGB <- function(C, M, Y, K, max = 100) {
  CMY <- matrix(data = c(C, M, Y), ncol = 3, byrow = FALSE) / max
  K <- K / max
  RGB <- 255 * (1 - CMY) * (1 - K)
  colnames(RGB) <- c("R", "G", "B")
  return(RGB)
}

## Conversion des coordonnées CMYB en RGB
RGB <- CMYK2RGB(
  C = S_FGEOL_crop$C_FOND,
  M = S_FGEOL_crop$M_FOND,
  Y = S_FGEOL_crop$J_FOND,
  K = S_FGEOL_crop$N_FOND
)

## Conversion en hex
geol_col <- grDevices::rgb(RGB, maxColorValue = 255)
names(geol_col) <- geol_col

## Ajouter comme nouvelle colone
S_FGEOL_crop$HEX_FOND <- geol_col

Pour afficher correctement les éléments ponctuels structuraux, on ajoute une colonne SYMB à la table attributaire de la couche P_STRUCT. Cette colonne contient le caractère correspondant à la valeur de la colonne CODE (cf. lexique national des éléments ponctuels de nature structurale) et permet d’afficher le symbole de l’élément ponctuel en utilisant la police de caractère BRGM_PStruct.ttf¹.

P_STRUCT_crop$SYMB <- factor(
  x = P_STRUCT_crop$CODE, 
  labels = c("!", "'", "#", "$", ">", "?")
)

Visualisation

Géologie

La nouvelle colonne HEX_FOND peut être utilisée avec ggplot2 comme paramètre esthétique de remplissage et utilisée comme palette de couleur grâce à la fonction scale_fill_identity().

GEO050K <- ggplot2::ggplot() +
  ggplot2::aes(fill = HEX_FOND) +
  ggplot2::geom_sf(data = S_FGEOL_crop, color = "black", size = 0.1) +
  ggplot2::scale_fill_identity() +
  ggplot2::scale_x_continuous(expand = c(0, 0)) +
  ggplot2::scale_y_continuous(expand = c(0, 0))

GEO050K +
  ggplot2::coord_sf(datum = sf::st_crs(2154)) +
  ggplot2::labs(
    x = "", y = "",
    caption = "Carte géologique harmonisée (BRGM)"
  ) +
  ggplot2::theme_bw()

On peut ensuite compléter la carte avec les couches L_STRUCT et P_STRUCT. La couche P_STRUCT est ajoutée à la carte comme texte avec geom_sf_text(), en passant les colonnes SYMB et AZIMUT en paramètres esthétiques.

GEO050K + 
  ggplot2::geom_sf(
    data = L_STRUCT_crop, 
    mapping = ggplot2::aes(size = WT_SYMB),
    inherit.aes = FALSE
  ) +
  ggplot2::geom_sf_text(
    data = P_STRUCT_crop, 
    mapping = ggplot2::aes(label = SYMB, angle = AZIMUT),
    family = "BRGM_PStruct", # Police de caractère contenant les symboles
    size = grid::unit(2.5, "mm"),
    inherit.aes = FALSE
  ) +
  ggplot2::coord_sf(datum = sf::st_crs(2154)) +
  ggplot2::scale_size_identity() +
  ggplot2::labs(
    x = "", y = "",
    caption = "Carte géologique harmonisée (BRGM)"
  ) +
  ggplot2::theme_bw()

Terrain

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

## Calcul des ombrages
## 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)

Enfin, on peut créer les couches pour habiller² le modèle avec les fonctions generate_line_overlay() et generate_polygon_overlay() :

## Créer un habillage de lignes
geo_line <- rayshader::generate_line_overlay(
  geometry = L_STRUCT_crop,
  extent = extent_ventoux,
  linewidth = 2,
  color = "black", 
  heightmap = ventoux
)

## Créer un habillage de polygones
geo_poly <- rayshader::generate_polygon_overlay(
  geometry = S_FGEOL_crop,
  extent = extent_ventoux,
  heightmap = ventoux,
  data_column_fill = "HEX_FOND",
  linecolor = "grey80",
  palette = geol_col,
  linewidth = 1
)

Si les ombrages réalisés sont satisfaisant, on peut les combiner avec les données d’élévation et visualiser le terrain et sa géologie :

## Visualiser le modèle 2D
ventoux_geol <- ventoux |> 
  rayshader::sphere_shade(colorintensity = 5, zscale = 1) |> 
  rayimage::render_image_overlay(geo_poly) |>
  rayimage::render_image_overlay(geo_line) |> 
  rayshader::add_shadow(lam, max_darken = 0.5) |>
  rayshader::add_shadow(amb, max_darken = 0.5) |>
  rayshader::add_shadow(tex, max_darken = 0.2)

ventoux_geol |>
  rayshader::plot_map()

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 = ventoux_geol,
  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()

Animation

Pour produire une vue animée, on réalise une rotation complète du modèle en générant une vue par position de la caméra :

rayshader::plot_3d(
  hillshade = ventoux_geol,
  heightmap = ventoux,
  zscale = 12, # Exagération verticale
  shadowdepth = -45,
  theta = -45,
  phi = 30,
  fov = 45,
  zoom = 0.45,
  windowsize = c(1024, 768),
)

## Créer un dossier où générer les vues
## (le nettoyer s'il existe déjà)
if (!dir.exists("tmp")) dir.create("tmp")
file.remove(list.files("tmp", pattern = "render", full.names = TRUE))

angles <- seq(0, 360, length.out = 1441)[-1]
for(i in seq_len(1440)) {
  rayshader::render_camera(theta = -45 + angles[i])
  rayshader::render_snapshot(
    filename = sprintf("tmp/render%i.png", i),
    title_text = "Mont-Ventoux, France | Source : BRGM, IGN",
    title_size = 25,
    title_bar_color = "#999999",
    title_color = "white",
    title_bar_alpha = 1
  )
}
rgl::close3d()

Les vues générées sont combinées en une vidéo à l’aide de la librairie ffmpeg :

system("ffmpeg -framerate 60 -i tmp/render%d.png -pix_fmt yuv420p ventoux.mp4")

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] generics_0.1.4     class_7.3-23       KernSmooth_2.23-26 lattice_0.22-7    
 [5] extrafontdb_1.1    rgl_1.3.24         hms_1.1.3          digest_0.6.37     
 [9] magrittr_2.0.4     evaluate_1.0.5     grid_4.5.1         RColorBrewer_1.1-3
[13] iterators_1.0.14   fastmap_1.2.0      foreach_1.5.2      doParallel_1.0.17 
[17] rprojroot_2.1.1    jsonlite_2.0.0     progress_1.2.3     e1071_1.7-16      
[21] DBI_1.2.3          scales_1.4.0       codetools_0.2-20   cli_3.6.5         
[25] crayon_1.5.3       rlang_1.1.6        units_0.8-7        base64enc_0.1-3   
[29] withr_3.0.2        yaml_2.3.10        raster_3.6-32      tools_4.5.1       
[33] parallel_4.5.1     rayimage_0.15.1    rayshader_0.37.3   dplyr_1.1.4       
[37] ggplot2_4.0.0      here_1.0.2         png_0.1-8          vctrs_0.6.5       
[41] R6_2.6.1           magick_2.9.0       proxy_0.4-27       lifecycle_1.0.4   
[45] classInt_0.4-11    htmlwidgets_1.6.4  pkgconfig_2.0.3    terra_1.8-70      
[49] pillar_1.11.1      gtable_0.3.6       glue_1.8.0         Rcpp_1.1.0        
[53] sf_1.0-21          xfun_0.53          tibble_3.3.0       tidyselect_1.2.1  
[57] rstudioapi_0.17.1  knitr_1.50         dichromat_2.0-0.1  extrafont_0.20    
[61] farver_2.1.2       htmltools_0.5.8.1  rmarkdown_2.30     Rttf2pt1_1.3.14   
[65] compiler_4.5.1     prettyunits_1.2.0  S7_0.2.0           sp_2.2-0

Les références

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

Janjou, D. 2004. « Descriptif des cartes géologiques à 1/50 000 format "vecteurs" ». Vol. RP–53473–FR. BRGM.

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

———. 2025. rayimage: Image Processing for Simulated Cameras. https://doi.org/10.32614/CRAN.package.rayimage.

Murdoch, Duncan, et Daniel Adler. 2025. rgl: 3D Visualization Using OpenGL. https://doi.org/10.32614/CRAN.package.rgl.

Pebesma, Edzer. 2018. « Simple Features for R: Standardized Support for Spatial Vector Data ». The R Journal 10 (1): 439‑46. https://doi.org/10.32614/RJ-2018-009.

Pebesma, Edzer, et Roger Bivand. 2023. Spatial Data Science: With applications in R. Chapman and Hall/CRC. https://doi.org/10.1201/9780429459016.

Wickham, Hadley. 2016. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.

Notes de bas de page

La même stratégie peut être utilisée avec la colonne NOTATION de la couche S_FGEOL et la police de caractère BRGM_NOT.ttf.↩︎
Voir https://www.tylermw.com/adding-open-street-map-data-to-rayshader-maps-in-r 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-10-01},
  url = {https://carnets.archeosciences-bordeaux.fr/outils/r-ventoux-geology/},
  langid = {fr-FR}
}

Veuillez citer ce travail comme suit :

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