Évolution du littoral de la grotte Cosquer

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 des données libres.

rayshader 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 animations1.

En complément de rayshader, quelques packages supplémentaires sont nécessaires : notamment raster pour le traitement des données raster et magick pour le traitement d’images.

## Installer les packages nécessaires
# pkg <- c("ggplot2", "khroma", "magick","magrittr", "raster", "rayshader")
# install.packages(pkg)

## Charger les packages
library(ggplot2)
library(khroma)
library(magick)
library(magrittr)
library(raster)
library(rayshader)

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 Litto3D produites et diffusées par l’IGN (2017) et le SHOM (2015) sous licence ouverte etalab.

Les données BD ALTI sont produites avec un pas de 75 m. Ces dernières utilisent le système de coordonnées projetées RGF93-Lambert 93. Trois dalles sont nécessaires :

  • 0825-6300
  • 0900-6225
  • 0900-6300
## Lire les données raster BD ALTI
BDALTI_1 <- raster::raster("BDALTIV2_75M_FXX_0825_6300_MNT_LAMB93_IGN69.asc")
BDALTI_2 <- raster::raster("BDALTIV2_75M_FXX_0900_6225_MNT_LAMB93_IGN69.asc")
BDALTI_3 <- raster::raster("BDALTIV2_75M_FXX_0900_6300_MNT_LAMB93_IGN69.asc")

## Fusionner les trois dalles
BDALTI <- raster::merge(BDALTI_1, BDALTI_2, BDALTI_3)

## Préciser le système de coordonnées
## RGF93 - Lambert 93 (EPSG:2154)
raster::crs(BDALTI) <- "+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

BDALTI
## class      : RasterLayer 
## dimensions : 2000, 2000, 4e+06  (nrow, ncol, ncell)
## resolution : 75, 75  (x, y)
## extent     : 824962.5, 974962.5, 6150038, 6300038  (xmin, xmax, ymin, ymax)
## crs        : +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : -47.17, 1466  (min, max)

Les données Litto3D 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 Litto3D
Litto3d <- raster::raster("MNT_MED100m_GDL-CA_HOMONIM_WGS84_NM_ZNEG.asc")

## Préciser le système de coordonnées
## WGS 84 (EPSG:4326)
raster::crs(Litto3d) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

Litto3d
## class      : RasterLayer 
## dimensions : 2701, 5001, 13507701  (nrow, ncol, ncell)
## resolution : 0.001, 0.001  (x, y)
## extent     : 2.8995, 7.9005, 41.6995, 44.4005  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : MNT_MED100m_GDL.CA_HOMONIM_WGS84_NM_ZNEG 
## values     : -2844.2, 9.99  (min, max)

Les données BD ALTI et Litto3D 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.

La fonction projectRaster() permet de spécifier un objet cible (to) dont les paramètres seront utilisés pour la projection de l’object de départ (from). 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).

## Projection des données Litto3D
## WGS 84 vers RGF93-Lambert 93
Litto3d <- raster::projectRaster(from = Litto3d, to = BDALTI, 
                                 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 découper selon une emprise carrée de 70 km de côté, approximativement centrée sur la grotte Cosquer, et à les fusionner en un seul objet.

## Définir l'emprise du découpage (xmin, xmax, ymin, ymax)
extent_2154 <- raster::extent(860000, 930000, 6193000, 6263000)

## Découper les données
BDALTI_crop <- raster::crop(x = BDALTI, y = extent_2154)
Litto3d_crop <- raster::crop(x = Litto3d, y = extent_2154)

## Fusion des données BD ALTI et Litto3D
cosquer_merge <- raster::merge(BDALTI_crop, Litto3d_crop)
cosquer_merge
## class      : RasterLayer 
## dimensions : 933, 934, 871422  (nrow, ncol, ncell)
## resolution : 75, 75  (x, y)
## extent     : 859987.5, 930037.5, 6193012, 6262988  (xmin, xmax, ymin, ymax)
## crs        : +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : -2360.077, 1144.46  (min, max)

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
raster::cellStats(cosquer_merge, stat = "countNA")
## [1] 41
## Estimer les valeurs manquantes
cosquer_focal <- raster::focal(
  cosquer_merge,
  w = matrix(1, nrow = 3, ncol = 3),
  fun = mean,
  NAonly = TRUE,
  na.rm = TRUE
)

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

## Convertir les données raster en data.frame
cosquer_point <- cosquer_focal %>%
  raster::rasterToPoints() %>% 
  as.data.frame()

## Visualiser
ggplot2::ggplot(cosquer_point) +
  ggplot2::aes(x = x, y = y, fill = layer) + 
  ggplot2::geom_raster() +
  ggplot2::coord_fixed(expand = FALSE) +
  ggplot2::theme_bw() +
  ggplot2::labs(
    x = "", y = "",
    fill = "Élévation (m)",
    caption = "IGN BD ALTI 2.0 - SHOM Litto3D"
  ) +
  khroma::scale_fill_sunset()

Terrain et hillshade

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
par(mfrow = c(1, 2))

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

Tyller Morgan-Wall expose en détail ces différentes méthodes.

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
par(mfrow = c(1, 2))

## 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)
## Visualiser avec les ombrages
par(mfrow = c(1, 2))

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) %>%
  rayshader::plot_map()

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) %>%
  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
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) %>%
  rayshader::plot_3d(
    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(1000, 700),
  )

## Ajouter des étiquettes
rayshader::render_label(
  cosquer,
  text = "Cosquer Cave",
  lat = 6236818,
  long = 899155,
  altitude = 12000,
  zscale = 37,
  extent = extent_2154,
  clear_previous = TRUE
)
rayshader::render_label(
  cosquer,
  text = "Marseille",
  lat = 6246798,
  long = 892378,
  altitude = 7000,
  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)
sea_stack <- read.table("spratt2016.txt", sep = "\t", dec = ".",
                        header = TRUE, skip = 95)
## Replacer le niveau zéro à 0 ka BP
sea_stack$SeaLev_shortPC1 <- sea_stack$SeaLev_shortPC1 - 8.49

## Conserver uniquement les 130 premières valeurs
## (Pléistocène supérieur et Holocène)
sea <- sea_stack[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 = "Sea level (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
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) %>%
  rayshader::plot_3d(
    cosquer,
    zscale = 75,
    water = TRUE,
    waterdepth = -130,
    watercolor = "lightblue",
    wateralpha = 0.5,
    windowsize = c(1000, 700),
  )

## Ajouter une étiquette
rayshader::render_label(
  cosquer,
  text = "Cosquer Cave",
  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()

Pour produire l’animation, on va générer une vue par pas de 1000 ans en faisant varier le niveau marin. Sur chaque vue2, 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 peinture3.

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

## Itérer sur les données de Spratt et Lisiecki 2016
## (limitées à 129 - 0 ka cal BP)
for (i in 1: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
  cosquer %>%
    rayshader::sphere_shade(zscale = 1) %>%
    rayshader::add_shadow(ray, max_darken = 0.5) %>%
    rayshader::add_shadow(lam, max_darken = 0.7) %>%
    rayshader::add_shadow(amb, max_darken = 0.1) %>%
    rayshader::plot_3d(
      cosquer,
      zscale = 37,
      water = TRUE,
      waterdepth = sea_level,
      watercolor = "lightblue",
      wateralpha = 0.5,
      theta = -40,
      phi = 30,
      fov = 0,
      zoom = 0.7,
      windowsize = c(1000, 700),
    )
  
  ## Ajouter les étiquettes
  rayshader::render_label(
    cosquer,
    text = "Cosquer Cave",
    lat = 6236818,
    long = 899155,
    altitude = 12000,
    zscale = 37,
    extent = extent_2154,
    clear_previous = TRUE
  )
  rayshader::render_label(
    cosquer,
    text = "Marseille",
    lat = 6246798,
    long = 892378,
    altitude = 7000,
    zscale = 37,
    extent = extent_2154
  )
  
  ## Enregistrer la vue
  img_file <- sprintf("temp/render%i.png", nrow(sea) - i)
  rayshader::render_snapshot(filename = img_file)
  ## Fermer la fenêtre 3D
  rgl::rgl.close()
  
  ## 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 = 660, labels = "Cosquer Cave, France", 
       cex = 3, adj = c(0, 0.5))
  text(x = 20, y = 620, labels = sprintf("%d cal ka BP", sea_age), 
       cex = 2, adj = c(0, 0.5))
  text(x = 20, y = 590, 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 Litto3D",
       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 vue générées sont combinées en une vidéo à l’aide de la librairie ffmpeg.

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

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.

Session

## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.1 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] rayshader_0.19.2 raster_3.4-5     sp_1.4-4         magrittr_1.5    
## [5] magick_2.4.0     khroma_1.4.0     ggplot2_3.3.2   
## 
## loaded via a namespace (and not attached):
##  [1] progress_1.2.2          tidyselect_1.1.0        xfun_0.19              
##  [4] purrr_0.3.4             lattice_0.20-41         colorspace_1.4-1       
##  [7] vctrs_0.3.3             generics_0.0.2          miniUI_0.1.1.1         
## [10] htmltools_0.5.0         yaml_2.2.1              rayimage_0.3.1         
## [13] rlang_0.4.7             manipulateWidget_0.10.1 pillar_1.4.6           
## [16] later_1.1.0.1           glue_1.4.2              withr_2.2.0            
## [19] foreach_1.5.0           lifecycle_0.2.0         stringr_1.4.0          
## [22] munsell_0.5.0           blogdown_0.21           gtable_0.3.0           
## [25] htmlwidgets_1.5.1       codetools_0.2-17        evaluate_0.14          
## [28] labeling_0.3            knitr_1.29              fastmap_1.0.1          
## [31] doParallel_1.0.15       httpuv_1.5.4            crosstalk_1.1.0.1      
## [34] parallel_4.0.3          Rcpp_1.0.5              xtable_1.8-4           
## [37] scales_1.1.1            promises_1.1.1          webshot_0.5.2          
## [40] jsonlite_1.7.0          farver_2.0.3            mime_0.9               
## [43] png_0.1-7               hms_0.5.3               digest_0.6.25          
## [46] stringi_1.4.6           bookdown_0.20           dplyr_1.0.2            
## [49] shiny_1.5.0             grid_4.0.3              rgdal_1.5-18           
## [52] tools_4.0.3             rgl_0.100.54            tibble_3.0.3           
## [55] crayon_1.3.4            pkgconfig_2.0.3         ellipsis_0.3.1         
## [58] prettyunits_1.1.1       rmarkdown_2.5           iterators_1.0.12       
## [61] R6_2.4.1                compiler_4.0.3

Bibliographie

Bache, Stefan Milton, et Hadley Wickham. 2014. magrittr: A Forward-Pipe Operator for R. https://CRAN.R-project.org/package=magrittr.

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. 2020. raster: Geographic Data Analysis and Modeling. https://CRAN.R-project.org/package=raster.

IGN. 2017. BD ALTI 75 m.

Morgan-Wall, Tyler. 2020. rayshader: Create Maps and Visualize Data in 2D and 3D. https://CRAN.R-project.org/package=rayshader.

Ooms, Jeroen. 2020. magick: Advanced Graphics and Image-Processing in R. https://CRAN.R-project.org/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.

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


  1. Un résultat similaire pourrait être obtenu avec QGIS et blender, mais R et rayshader offrent la possibilité de produire des figures reproductibles.↩︎

  2. Inspiré de A 3D tour over Lake Geneva with rayshader.↩︎

  3. Très grossièrement approximée à partir des données de Clottes et al. (1992).↩︎

Nicolas Frerebeau
Nicolas Frerebeau
Archaeologist

Technological choices in ancient societies.

Related