Géologie du Mont-Ventoux

Réaliser une carte géologique avec R.

R
Geomorphology
Geology
Auteur·rice
Affiliation
N. Frerebeau

IRAMAT-CRP2A

Date de publication

19 décembre 2020

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 du programme STRM de la NASA pour l’élévation.

Pour y parvenir, quatre packages R sont nécessaires en plus de ggplot2 (Wickham 2016) : raster (Hijmans 2020) et sf (Pebesma 2018) pour le traitement des données raster et vectorielles, ainsi que rayshader (Morgan-Wall 2020), qui permet de visualiser en 2D et 3D des données d’élévation.

## Installer les packages nécessaires
## La version 0.21.0 ou supérieure de rayshader est nécessaire
pkg <- c("ggplot2", "magrittr", "raster", "rayshader", "sf")
install.packages(pkg)
## Charger les packages
library(magrittr)

La zone d’intérêt est un carré de 40 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_RGF93 <- raster::extent(861000, 901000, 6320000, 6360000)

Géologie

Préparation des données

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, la Drôme, le Gard et le Vaucluse :

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

## Départements de la Drôme (26), du Gard (30) et du Vaucluse (84)
dpt <- c("GEO050K_HARM_026.zip", "GEO050K_HARM_030.zip", "GEO050K_HARM_084.zip")
for (i in dpt) {
  ## Téléchargement des données
  utils::download.file(paste0(BRGM_url, i), destfile = i)
  ## 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_26 <- sf::st_read("GEO050K_HARM_026/GEO050K_HARM_026_S_FGEOL_2154.shp")
S_FGEOL_30 <- sf::st_read("GEO050K_HARM_030/GEO050K_HARM_030_S_FGEOL_2154.shp")
S_FGEOL_84 <- sf::st_read("GEO050K_HARM_084/GEO050K_HARM_084_S_FGEOL_2154.shp")

L_STRUCT_26 <- sf::st_read("GEO050K_HARM_026/GEO050K_HARM_026_L_STRUCT_2154.shp")
L_STRUCT_30 <- sf::st_read("GEO050K_HARM_030/GEO050K_HARM_030_L_STRUCT_2154.shp")
L_STRUCT_84 <- sf::st_read("GEO050K_HARM_084/GEO050K_HARM_084_L_STRUCT_2154.shp")

P_STRUCT_26 <- sf::st_read("GEO050K_HARM_026/GEO050K_HARM_026_P_STRUCT_2154.shp")
P_STRUCT_30 <- sf::st_read("GEO050K_HARM_030/GEO050K_HARM_030_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_26, L_STRUCT_30, L_STRUCT_84)
P_STRUCT <- rbind(P_STRUCT_26, P_STRUCT_30, P_STRUCT_84)
S_FGEOL <- rbind(S_FGEOL_26, S_FGEOL_30, S_FGEOL_84)

## Unions des polygones par attributs
S_FGEOL <- aggregate(
  x = S_FGEOL_84, 
  by = list(S_FGEOL_84$NOTATION), 
  FUN = function(x) x[1]
)
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_RGF93)
L_STRUCT_crop <- sf::st_crop(x = L_STRUCT, y = extent_RGF93)
P_STRUCT_crop <- sf::st_crop(x = P_STRUCT, y = extent_RGF93)

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

## Attribution des couleurs
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.ttf1.

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

Cartographie

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()

Élévation

Préparation des données

Les données d’élévation nécessaires sont issues du STRM v3 (NASA JPL 2013) et peuvent être téléchargées librement via le portail Earthdata de la NASA. Les données utilisées ici sont produites avec une résolution spatiale de 1 arcseconde (env. 30 m) et utilisent le système de coordonnées WGS 84. Deux dalles sont nécessaires :

  • N43E005
  • N44E005

La première étape consiste à fusionner les deux dalles, puis à projeter les données en utilisant le système de coordonnées RGF93-Lambert 93. La projection est réalisée avec la fonction projectRaster() en précisant la résolution de sortie (ici 30 m) et en laissant la méthode utilisée à sa valeur par défaut (l’interpolation bilinéaire est adaptée aux données continues, ce qui est le cas des données d’élévation).

## Lire les données raster SRTM
N44E005 <- raster::raster("N44E005.hgt")
N43E005 <- raster::raster("N43E005.hgt")

## Fusion des raster
SRTM <- raster::merge(N44E005, N43E005)

## Découper grossièrement le raster fusionné
## (évite d'avoir à projeter toutes les données)
extent_WGS84 <- raster::extent(5.0, 5.6, 43.5, 44.5)
SRTM_WGS84 <- raster::crop(x = SRTM, y = extent_WGS84)

## Projection des données SRTM
## WGS 84 vers RGF93-Lambert 93 (EPSG:2154)
RGF93 <- "+init=epsg:2154"
SRTM_RGF93 <- raster::projectRaster(from = SRTM_WGS84, res = 30, crs = RGF93)

## Découper les données
SRTM_crop <- raster::crop(x = SRTM_RGF93, y = extent_RGF93)

Ombrages du terrain

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ère2.

## Convertir les données raster en matrice
ventoux <- rayshader::raster_to_matrix(SRTM_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 = 15)

## 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).

## Visualiser les ombrages
par(mfrow = c(2, 2))

## Hillshading
rayshader::plot_map(lam)

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

## Texture shading
rayshader::plot_map(tex)

## Sphere shading
rayshader::plot_map(sph)

## Elévation et ombrages
ventoux_map <- 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)

ventoux_map %>% 
  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
ventoux_map %>% 
  rayshader::plot_3d(
    ventoux,
    zscale = 15, # Exagération verticale
    theta = -45,
    phi = 30,
    fov = 0,
    zoom = 0.6,
    windowsize = c(1000, 700),
  )

## Ajouter des étiquettes
# rayshader::render_label(
#   ventoux,
#   text = "Mont-Ventoux",
#   lat = 6344344,
#   long = 882196,
#   altitude = 3000,
#   zscale = 15,
#   extent = extent_RGF93,
#   clear_previous = TRUE
# )
# rayshader::render_label(
#   ventoux,
#   text = "Mont-Serein",
#   lat = 6345506,
#   long = 880718,
#   altitude = 3000,
#   zscale = 15,
#   extent = extent_RGF93,
#   clear_previous = FALSE
# )

## Enregistrer la vue
rayshader::render_snapshot(clear = TRUE)

Terrain et géologie

Enfin, on peut créer les couches pour habiller3 le modèle avec les fonctions generate_*_overlay() et add_overlay() :

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

## Créer un habillage de polygones
geo_poly <- rayshader::generate_polygon_overlay(
  geometry = S_FGEOL_crop,
  extent = extent_RGF93,
  heightmap = ventoux,
  data_column_fill = "HEX_FOND",
  linecolor = "grey80",
  palette = geol_col,
  linewidth = 1
)
## Visualiser le modèle 2D
ventoux_geol <- ventoux %>%
  rayshader::sphere_shade() %>%
  rayshader::add_overlay(geo_poly) %>%
  rayshader::add_overlay(geo_line) %>% 
  rayshader::add_shadow(lam, max_darken = 0.3) %>%
  rayshader::add_shadow(amb, max_darken = 0.5) %>%
  rayshader::add_shadow(tex, max_darken = 0.3)

ventoux_geol %>%
  rayshader::plot_map()

## Visualiser le modèle 3D
ventoux_geol %>%
  rayshader::plot_3d(
    ventoux,
    zscale = 15, # Exagération verticale
    theta = -45,
    phi = 30,
    fov = 0,
    zoom = 0.6,
    windowsize = c(1000, 700),
  )

## Enregistrer la vue
rayshader::render_snapshot(clear = TRUE)

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.

ventoux_geol %>%
  rayshader::plot_3d(
    ventoux,
    zscale = 15, # Exagération verticale
    shadowdepth = -50,
    theta = -45,
    phi = 30,
    fov = 45,
    zoom = 0.45,
    windowsize = c(1000, 700),
  )

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

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

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

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

Session

R version 4.2.1 (2022-06-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.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       

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

other attached packages:
[1] magrittr_2.0.3

loaded via a namespace (and not attached):
 [1] rgl_0.109.6        Rcpp_1.0.9         here_1.0.1         lattice_0.20-45   
 [5] png_0.1-7          prettyunits_1.1.1  class_7.3-20       assertthat_0.2.1  
 [9] rprojroot_2.0.3    digest_0.6.29      foreach_1.5.2      utf8_1.2.2        
[13] R6_2.5.1           evaluate_0.16      e1071_1.7-11       ggplot2_3.3.6     
[17] pillar_1.8.1       rlang_1.0.5        progress_1.2.2     rstudioapi_0.14   
[21] extrafontdb_1.0    raster_3.5-29      magick_2.7.3       rmarkdown_2.16    
[25] extrafont_0.18     stringr_1.4.1      htmlwidgets_1.5.4  munsell_0.5.0     
[29] proxy_0.4-27       compiler_4.2.1     xfun_0.32          base64enc_0.1-3   
[33] pkgconfig_2.0.3    htmltools_0.5.3    tidyselect_1.1.2   tibble_3.1.8      
[37] codetools_0.2-18   fansi_1.0.3        crayon_1.5.1       dplyr_1.0.10      
[41] rayshader_0.28.2   sf_1.0-8           wk_0.6.0           grid_4.2.1        
[45] Rttf2pt1_1.3.10    jsonlite_1.8.0     gtable_0.3.1       lifecycle_1.0.1   
[49] DBI_1.1.3          units_0.8-0        scales_1.2.1       KernSmooth_2.23-20
[53] cli_3.4.0          stringi_1.7.8      farver_2.1.1       doParallel_1.0.17 
[57] sp_1.5-0           ellipsis_0.3.2     generics_0.1.3     vctrs_0.4.1       
[61] s2_1.1.0           iterators_1.0.14   tools_4.2.1        rayimage_0.6.5    
[65] glue_1.6.2         purrr_0.3.4        hms_1.1.2          parallel_4.2.1    
[69] fastmap_1.1.0      yaml_2.3.5         colorspace_2.0-3   terra_1.6-7       
[73] classInt_0.4-7     knitr_1.40

Les références

Brown, Leland H. 2014. « TEXTURE SHADING : A NEW TECHNIQUE FOR DEPICTING TERRAIN RELIEF ». In.
Hijmans, Robert J. 2020. raster: Geographic Data Analysis and Modeling. https://CRAN.R-project.org/package=raster.
Janjou, D. 2004. « Descriptif des cartes géologiques à 1/50 000 format "vecteurs" ». Vol. RP–53473–FR. BRGM.
Morgan-Wall, Tyler. 2020. rayshader: Create Maps and Visualize Data in 2D and 3D. https://CRAN.R-project.org/package=rayshader.
NASA JPL. 2013. NASA Shuttle Radar Topography Mission Global 1 arc second number [Data set]. NASA EOSDIS Land Processes DAAC. https://doi.org/10.5067/MEaSUREs/SRTM/SRTMGL1N.003.
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.
Wickham, Hadley. 2016. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.

Notes de bas de page

  1. 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.↩︎

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

  3. Voir https://www.tylermw.com/adding-open-street-map-data-to-rayshader-maps-in-r pour plus de détails.↩︎

Réutilisation

CC BY-NC-SA

Citation

BibTeX
@online{frerebeau2021,
  author = {N. Frerebeau},
  title = {Géologie du Mont-Ventoux},
  date = {2021-07-12},
  url = {https://github.com/nfrerebeau/website/posts/2020-12-19-r-ventoux},
  langid = {fr-FR}
}
Veuillez citer ce travail comme suit :
N. Frerebeau. 2021. “Géologie du Mont-Ventoux.” July 12, 2021. https://github.com/nfrerebeau/website/posts/2020-12-19-r-ventoux.