On cherche ici à produire une vue du littoral de la grotte Cosquer (dont l’entrée est actuellement située à 37 m sous la surface de la mer Méditerranée) et à animer les variations du niveau marin depuis le Pléistocène supérieur. Le tout avec R et des données libres.

rayshader (Morgan-Wall 2024) est un package R qui permet de visualiser en 2D et 3D des données d’élévation. Ce package permet d’explorer les modèles 3D de façon interactive et de créer des animations¹.

En complément de rayshader, quelques packages supplémentaires sont nécessaires : notamment terra (Hijmans 2025) pour le traitement des données raster, rgl² pour le rendu des scènes 3D (Murdoch et Adler 2025) et magick (Ooms 2025) pour le traitement d’images.

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

Préparation des données

Pour commencer, il faut récupérer les données d’élévation terrestres et sous-marines. Ces dernières sont librement accessibles grâce aux ressources BD ALTI 2.0 et MNT bathymétrique produites et diffusées par l’IGN (2023) et le SHOM (2015) sous licence ouverte etalab.

Les données BD ALTI sont produites avec un pas de 25 m. Ces dernières utilisent le système de coordonnées projetées RGF93-Lambert 93. Six dalles du département des Bouches-du-Rhône (13) sont nécessaires :

0850-6250
0850-6275
0875-6250
0875-6275
0900-6250
0900-6275

## Lire les données raster BD ALTI
## (les fichiers .asc sont placés dans le répertoire de travail)
BDALTI_files <- list.files(
  path = ".",
  pattern = "MNT_LAMB93_IGN69.asc",
  full.names = TRUE
)

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

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

## Visualisation
# terra::plot(BDALTI_vrt)

Les données du MNT bathymétrique Golfe du Lion – Côte d’Azur sont produites avec une résolution de 0.001° (environ 111 m) et utilisent le système de coordonnées géographiques WGS 84.

## Lire les données raster MNT bathymétrique
MNTbathy <- terra::rast("MNT_MED100m_GDL-CA_HOMONIM_WGS84_NM_ZNEG.asc")

## Préciser le système de coordonnées
## WGS 84 (EPSG:4326)
terra::crs(MNTbathy) <- "epsg:4326"

## Visualisation
# terra::plot(MNTbathy)

Les deux jeux de données utilisent des systèmes de coordonnées différents. Il faut donc projeter l’ensemble des données dans le même référentiel et les découper selon une emprise carrée de 75 km de côté, approximativement centrée sur la grotte Cosquer.

La fonction project() permet de spécifier un objet cible (y) dont les paramètres seront utilisés pour la projection de l’object de départ (x). Le résultat est un raster projeté dans le nouveau système de coordonées et rééchantillonné (il a la même résolution que l’objet cible). La méthode utilisée est laissé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).

## Définir l'emprise du découpage (xmin, xmax, ymin, ymax)
extent_2154 <- terra::ext(850000, 925000, 6190000, 6265000)

## Projection du MNT bathymétrique
## WGS 84 vers RGF93-Lambert 93
MNTbathy_crop <- terra::project(
  x = MNTbathy,
  y = terra::rast(extent_2154, crs = "epsg:2154", res = 75),
  method = "bilinear"
)

Les deux jeux de données sont désormais projetés dans le même référentiel et ont la même résolution (75 m). Il reste à les fusionner en un seul objet.

## Découper les données BD ALTI
BDALTI_crop <- terra::crop(x = BDALTI_vrt, y = extent_2154)

## Fusion des données BD ALTI et MNT bathymétrique
cosquer_merge <- terra::merge(BDALTI_crop, MNTbathy_crop, first = FALSE)

Dernière étape : s’il existe des valeurs manquantes, elles sont remplacées par une estimation. Cette estimation correspond à la moyenne des huit valeurs immédiatement voisines.

## Compter le nombre de valeurs manquantes
terra::global(cosquer_merge, fun = "isNA")

                 isNA
file39fd18adeb12 1001

## Estimer les valeurs manquantes
cosquer_focal <- terra::focal(
  cosquer_merge,
  w = matrix(1, nrow = 3, ncol = 3),
  fun = mean,
  na.policy = "only",
  na.rm = TRUE
)

Enfin, on peut visualiser le résultat en colorant la surface en fonction de l’élévation absolue.

## Distribution des valeurs d'élévation
## Calcul de l'histogramme
h <- terra::hist(
  x = cosquer_focal, 
  xlab = "Élévation (m)",
  ylab = "Fréquence", 
  main = ""
)

## Échelle de couleur en fonction des classes de l'histogramme
pal <- khroma::color("bukavu")
cols <- khroma::palette_color_continuous(pal, midpoint = 0)(h$breaks)

## Visualiser
terra::plot(cosquer_focal, breaks = h$breaks, col = cols)

Visualisation

Terrain et ombrage

Au lieu d’utiliser l’élévation absolue du modèle de terrain, la fonction sphere_shade() du package rayshader permet 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).

Il est également possible d’ajouter de l’eau grâce à la fonction add_water(). Pour cela il suffit de créer une matrice d’incidence indiquant la présence d’eau (c’est-à-dire partout où l’élévation est inférieure ou égale à zéro).

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

## Visualiser
## Données d'élévation
cosquer |>
  rayshader::sphere_shade() |>
  rayshader::plot_map()
## Ajouter de l'eau lorsque l'altitude est inférieure à 0
## (il faut inverser l'ordre des lignes dans la matrice d'incidence)
cosquer_water <- cosquer[nrow(cosquer):1, ] <= 0
cosquer |>
  rayshader::sphere_shade() |>
  rayshader::add_water(cosquer_water) |> 
  rayshader::plot_map()

Cette méthode simple de hillshading (dite lambertienne) peut être combinée avec les fonctions ray_shade() qui permet de calculer les ombrage par ray tracing et ambient_shade() qui permet de calculer les ombres liées à la diffusion de la lumière par l’atmophère³.

Le paramètre zscale correspond au facteur d’échelle entre coordonnées horizontales et verticale. La position du soleil peut également être ajustée au sein des fonction ray_shade() et lamb_shade() à l’aide des arguments sunaltitude et sunangle.

## Visualiser les ombrages
## Hillshading
cosquer |> 
  rayshader::lamb_shade(zscale = 75) |>
  rayshader::plot_map()
## Occlusion atmosphérique
cosquer |> 
  rayshader::ambient_shade() |> 
  rayshader::plot_map()

On peut visualiser des données d’élévation avec l’ensemble des couches d’ombrage.

## Calcul des ombrages
## Ray tracing
ray <- rayshader::ray_shade(cosquer, zscale = 75, lambert = FALSE)

## Hillshading
lam <- rayshader::lamb_shade(cosquer, zscale = 75)

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

## Combiner les ombrages
cosquer_shade <- cosquer |>
  rayshader::sphere_shade() |>
  rayshader::add_shadow(ray, max_darken = 0.5) |>
  rayshader::add_shadow(lam, max_darken = 0.7) |>
  rayshader::add_shadow(amb, max_darken = 0.1)

## Visualiser avec les ombrages
cosquer_shade |>
  rayshader::plot_map()
cosquer_shade |>
  rayshader::add_water(cosquer_water) |> 
  rayshader::plot_map()

Scène 3D

Désormais on peut facilement 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(
  cosquer_shade,
  cosquer,
  zscale = 37, # Exagération verticale
  water = TRUE,
  waterdepth = 0,
  watercolor = "lightblue",
  wateralpha = 0.5,
  waterlinecolor = "white",
  waterlinealpha = 1,
  theta = -40,
  phi = 30,
  fov = 0,
  zoom = 0.6,
  windowsize = c(2048, 1152),
)

## Ajouter des étiquettes
rayshader::render_label(
  cosquer,
  text = "Grotte Cosquer",
  textsize = 2,
  lat = 6236818,
  long = 899155,
  altitude = 13000,
  zscale = 37,
  extent = extent_2154,
  clear_previous = TRUE
)
rayshader::render_label(
  cosquer,
  text = "Marseille",
  textsize = 1.5,
  lat = 6246798,
  long = 892378,
  altitude = 9000,
  zscale = 37,
  extent = extent_2154
)

## Enregistrer la vue
rayshader::render_snapshot()

Animation du niveau marin

Dernière étape, on veut animer les variations du niveau marin au cours du temps. On dispose pour cela de la reconstruction du niveau marin global par Spratt et Lisiecki (2016).

## Importer les données de Spratt et Lisiecki 2016
## (les 95 premières lignes contiennent les métadonnées)
data("spratt2016", package = "folio")

## Replacer le niveau zéro à 0 ka BP
spratt2016$SeaLev_shortPC1 <- spratt2016$SeaLev_shortPC1 - 8.49

## Conserver uniquement les 130 premières valeurs
## (Pléistocène supérieur et Holocène)
sea <- spratt2016[1:130, 1:2]

## Tracer les données
par(mar = c(4, 4, 0, 0) + 0.1, las = 1)
plot(
  x = sea$age_calkaBP,
  y = sea$SeaLev_shortPC1,
  type = "l",
  xlab = "Age (cal ka BP)",
  ylab = "Niveau marin (m)",
  xlim = rev(range(sea$age_calkaBP)),
  bty = "n"
)

Le niveau le plus bas (-130 m) est atteint à 24 ka cal BP (approximativement synchrone de la réalisation des peintures). On peut illustrer le littoral à ce moment (sans exagération verticale cette fois) :

## Visualiser le modèle 3D
rayshader::plot_3d(
  cosquer_shade,
  cosquer,
  zscale = 75,
  water = TRUE,
  waterdepth = -130,
  watercolor = "lightblue",
  wateralpha = 0.5,
  zoom = 0.7,
  windowsize = c(2048, 1536),
)

## Ajouter une étiquette
rayshader::render_label(
  cosquer,
  text = "Grotte Cosquer",
  linecolor = "white",
  textcolor = "white",
  textsize = 2,
  linewidth = 3,
  lat = 6236818,
  long = 899155,
  altitude = 3000,
  zscale = 75,
  extent = extent_2154,
  clear_previous = TRUE
)

## Repositionner la caméra
rayshader::render_camera(theta = -15, phi = 20, zoom = 0.2, fov = 0)

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

Pour produire l’animation, on va générer une vue par pas de 1000 ans en faisant varier le niveau marin. Sur chaque vue⁴, on va également afficher l’âge et la valeur du niveau marin, ainsi que la reconstruction du niveau marin sur 129 ka en surlignant la période de réalisation des peinture⁵.

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

## Itérer sur les données de Spratt et Lisiecki 2016
## (limitées à 129 - 0 ka cal BP)
for (i in seq_len(nrow(sea))) {
  ## Extraire l'âge et le niveau marin
  sea_age <- sea$age_calkaBP[[i]]
  sea_level <- sea$SeaLev_shortPC1[[i]]
  
  ## Générer la vue 3D
  rayshader::plot_3d(
    cosquer_shade,
    cosquer,
    zscale = 37,
    water = TRUE,
    waterdepth = sea_level,
    watercolor = "lightblue",
    wateralpha = 0.5,
    theta = -40,
    phi = 30,
    fov = 0,
    zoom = 0.7,
    windowsize = c(1024, 768),
  )
  
  ## Ajouter les étiquettes
  rayshader::render_label(
    cosquer,
    text = "Cosquer Cave",
    lat = 6236818,
    long = 899155,
    altitude = 15000,
    zscale = 37,
    extent = extent_2154,
    clear_previous = TRUE
  )
  rayshader::render_label(
    cosquer,
    text = "Marseille",
    lat = 6246798,
    long = 892378,
    altitude = 8000,
    zscale = 37,
    extent = extent_2154
  )
  
  ## Enregistrer la vue
  img_file <- sprintf("tmp/render%i.png", nrow(sea) - i)
  rayshader::render_snapshot(filename = img_file)
  ## Fermer la fenêtre 3D
  rgl::close3d()
  
  ## Générer le graphique niveau marin vs âge
  img_plot <- magick::image_graph(width = 350, height = 220)
  par(mar = c(4, 4, 0, 0) + 0.1, bg = "transparent", las = 1)
  plot(
    x = sea$age_calkaBP,
    y = sea$SeaLev_shortPC1,
    type = "l",
    xlab = "Age (cal ka BP)",
    ylab = "Sea level (m)",
    xlim = rev(range(sea$age_calkaBP)),
    bty = "n"
  )
  ## Surligner la période de réalisation des peintures
  polygon(x = c(32, 32, 22, 22), y = c(-140, 0, 0, -140), 
          col = adjustcolor("red", alpha.f = 0.5), border = NA)
  lines(x = sea_age, y = sea_level, type = "p", pch = 16, cex = 1.5)
  arrows(x0 = 45, y0 = -10, x1 = 33, y1 = -10, length = 0.1)
  text(x = 45, y = -10, labels = "paintings", pos = 2)
  dev.off()
  
  ## Lire la vue enregistrée précédemment
  img <- magick::image_read(img_file)
  ## Ajouter le graphique
  img <- magick::image_composite(img, img_plot, offset = "+630+10")
  img <- magick::image_draw(
    img, 
    xlim = c(0, magick::image_info(img)[["width"]]), 
    ylim = c(0, magick::image_info(img)[["height"]])
  )
  ## Ajouter le texte
  text(x = 20, y = 740, labels = "Grotte Cosquer, France", 
       cex = 3, adj = c(0, 0.5))
  text(x = 20, y = 700, labels = sprintf("%d cal ka BP", sea_age), 
       cex = 2, adj = c(0, 0.5))
  text(x = 20, y = 670, labels = sprintf("Sea level: %g m", sea_level), 
       cex = 2, adj = c(0, 0.5))
  text(x = 20, y = 40, labels = "Elevation data: IGN BD ALTI 2.0, Shom MNT bathymétrique",
       cex = 1, adj = c(0, 0.5))
  text(x = 20, y = 20, labels = "Sea level data: Spratt & Lisiecki, 2016", 
       cex = 1, adj = c(0, 0.5))
  dev.off()
  
  ## Enregistrer la vue finale
  magick::image_write(img, path = img_file)
}

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

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

Note

Cette animation est uniquement illustrative. Elle suppose que la morphologie du terrain est restée inchangée depuis le Pléistocène supérieur et utilise une estimation globale du niveau marin pour représenter des variations locales.

Retour au sommet

Session

R version 4.5.2 (2025-10-31)
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.31         jsonlite_2.0.0     compiler_4.5.2     rayshader_0.37.3  
 [5] crayon_1.5.3       Rcpp_1.1.0         magick_2.9.0       parallel_4.5.2    
 [9] dichromat_2.0-0.1  progress_1.2.3     scales_1.4.0       png_0.1-8         
[13] yaml_2.3.10        fastmap_1.2.0      R6_2.6.1           khroma_1.17.0     
[17] knitr_1.50         iterators_1.0.14   htmlwidgets_1.6.4  RColorBrewer_1.1-3
[21] rlang_1.1.6        Rttf2pt1_1.3.14    terra_1.8-80       xfun_0.54         
[25] doParallel_1.0.17  cli_3.6.5          magrittr_2.0.4     digest_0.6.39     
[29] foreach_1.5.2      grid_4.5.2         rstudioapi_0.17.1  base64enc_0.1-3   
[33] hms_1.1.4          lifecycle_1.0.4    prettyunits_1.2.0  vctrs_0.6.5       
[37] extrafont_0.20     evaluate_1.0.5     glue_1.8.0         rayimage_0.15.1   
[41] farver_2.1.2       extrafontdb_1.1    codetools_0.2-20   rmarkdown_2.30    
[45] tools_4.5.2        pkgconfig_2.0.3    htmltools_0.5.8.1

Les références

Clottes, J., J. Courtin, H. Valladas, H. Cachier, N. Mercier, et M. Arnold. 1992. « La grotte Cosquer datée ». Bulletin de la Société préhistorique française 89 (8): 230‑34. https://doi.org/10.3406/bspf.1992.9527.

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

IGN. 2023. BD ALTI 75 m.

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

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

Ooms, Jeroen. 2025. magick: Advanced Graphics and Image-Processing in R. https://doi.org/10.32614/CRAN.package.magick.

SHOM. 2015. MNT Bathymétrique de façade Golfe du Lion - Côte d’Azur (Projet Homonim). https://doi.org/10.17183/MNT_MED100m_GDL_CA_HOMONIM_WGS84.

Spratt, Rachel M., et Lorraine E. Lisiecki. 2016. « A Late Pleistocene Sea Level Stack ». Climate of the Past 12 (4): 1079‑92. https://doi.org/10.5194/cp-12-1079-2016.

Notes de bas de page

Un résultat similaire pourrait être obtenu avec QGIS et blender, mais R et rayshader offrent la possibilité de produire des figures reproductibles.↩︎
Sous Ubuntu, OpenGL doit être installé au préalable (sudo apt install libgl1-mesa-dev libglu1-mesa-dev freeglut3-dev).↩︎
Voir https://github.com/tylermorganwall/MusaMasterclass pour le détail de ces méthodes.↩︎
Inspiré de A 3D tour over Lake Geneva with rayshader.↩︎
Très grossièrement approximée à partir des données de Clottes et al. (1992).↩︎

Réutilisation

CC BY 4.0

Citation

BibTeX

@online{frerebeau2020,
  author = {Frerebeau, N.},
  title = {Évolution du littoral de la grotte Cosquer},
  date = {2020-08-05},
  url = {https://carnets.archeosciences-bordeaux.fr/outils/r-cosquer/},
  langid = {fr-FR}
}

Veuillez citer ce travail comme suit :

Frerebeau, N. 2020. “Évolution du littoral de la grotte Cosquer.” August 5, 2020. https://carnets.archeosciences-bordeaux.fr/outils/r-cosquer/.