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 2020) 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 (Hijmans 2020) pour le traitement des données raster et magick (Ooms 2020) 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(magrittr)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 -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) <- sf::st_crs(2154)
BDALTILe chargement a nécessité le package : rasterLe chargement a nécessité le package : spclass : 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) <- sf::st_crs(4326)
Litto3dclass : 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_mergeclass : 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 (Wickham 2016) 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_vik()
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ère2.
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 = "Grotte Cosquer",
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(clear = TRUE)
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 = "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 vue3, 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 peinture4.
## 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 = "Grotte Cosquer, 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")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] raster_3.5-29 sp_1.5-0 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] prettyunits_1.1.1 png_0.1-7 assertthat_0.2.1 rprojroot_2.0.3
[9] digest_0.6.29 foreach_1.5.2 utf8_1.2.2 R6_2.5.1
[13] evaluate_0.16 ggplot2_3.3.6 pillar_1.8.1 rlang_1.0.5
[17] progress_1.2.2 rstudioapi_0.14 extrafontdb_1.0 magick_2.7.3
[21] rmarkdown_2.16 labeling_0.4.2 rgdal_1.5-32 extrafont_0.18
[25] stringr_1.4.1 htmlwidgets_1.5.4 munsell_0.5.0 compiler_4.2.1
[29] xfun_0.32 pkgconfig_2.0.3 base64enc_0.1-3 htmltools_0.5.3
[33] tidyselect_1.1.2 tibble_3.1.8 khroma_1.9.0 codetools_0.2-18
[37] fansi_1.0.3 crayon_1.5.1 dplyr_1.0.10 rayshader_0.28.2
[41] grid_4.2.1 jsonlite_1.8.0 Rttf2pt1_1.3.10 gtable_0.3.1
[45] lifecycle_1.0.1 DBI_1.1.3 scales_1.2.1 cli_3.4.0
[49] stringi_1.7.8 farver_2.1.1 doParallel_1.0.17 ellipsis_0.3.2
[53] generics_0.1.3 vctrs_0.4.1 iterators_1.0.14 tools_4.2.1
[57] glue_1.6.2 purrr_0.3.4 hms_1.1.2 parallel_4.2.1
[61] fastmap_1.1.0 yaml_2.3.5 colorspace_2.0-3 terra_1.6-7
[65] knitr_1.40 Les références
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.↩︎
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
Citation
@online{frerebeau2020,
author = {N. Frerebeau},
title = {Évolution du littoral de la grotte Cosquer},
date = {2020-05-08},
url = {https://github.com/nfrerebeau/website/posts/2020-08-05-r-cosquer},
langid = {fr-FR}
}