Géocomputation d’un indice de géodiversité - Application au Calvados

1 La géodiversité : définition et intérêt scientifique

Ce script vise à automatiser, de manière (presque) complètement reproductible, la géocomputation d’un indice de géodiversité dans le Calvados. La géodiversité, concept inventé dans les années 1990 par les géologues, résume la variété géologique, géomorphologique, pédologique et hydrologique d’un territoire. On suppose que la géodiversité est positivement corrélée à la biodiversité, ce qui en fait un proxy à moindre frais très utile pour établir des recommandations de conservation.

Illustration d’une faible et forte géodiversité, d’après Maliniemi et al. (2024)

Il y a plusieurs manière d’évaluer la géodiversité. Nous allons ici choisir la méthode quantitative, consistant à calculer des indices de Shannon pour chaque sous-dimension (hydrologie, pédologie, etc.), puis à sommer ces indices, comme illustré ci-dessous.

Calcul de l’indice de géodiversité : une somme d’indices de Shannon, d’après Gakou et al. (in prep.)

Pour rappel, l’indice de Shannon est un indice de diversité qui prend en compte à la fois la richesse (le nombre d’espèces ou de catégories d’une variable qualitative présent dans une unité spatiale) et l’abondance relative (le nombre d’individus dans chaque catégorie/espèce de cette unité). L’indice varie entre 0 (une seule catégorie présente) et ln(S), S étant la richesse. Plus l’indice est grand, plus la diversité est grande.

L’indice est défini comme :

\[ H' = - \sum_{i=1}^{S} p_i \log(p_i) \]

où :

• \(S\) est le nombre total de classes ou de catégories,

• \(p_i\) est la proportion d’observations dans la classe \(i\).

2 Préparation de l’analyse

L’indice sera calculé au sein d’une grille à mailles régulières, que nous allons paramétrer dès maintenant, puis importer les données. A ce stade, nous allons également définir l’emprise de la zone d’étude (ici le département du Calvados).

2.1 Activation des packages

Mais avant cela, activons les librairies nécessaires aux analyses à venir :

library(dplyr)
library(tidyr)
library(sf)
library(terra)
library(tmap)
library(vegan)

2.2 Création de la grille et défintion de l’emprise

Deux paramètres généraux peuvent être définis ici, en amont des analyses subséquentes : la résolution de la grille et l’emprise :

cellsize <- 5000 # En mètres si EPSG:2154
emprise <- st_read("data/DEP14.gpkg", quiet = TRUE) # Fichier disponible ici : https://gitlab.huma-num.fr/tfeuillet/cours/-/raw/main/data/DEP14.gpkg?ref_type=heads

Créons maintenant la grille selon notre emprise :

grid <- st_make_grid(x = emprise, cellsize = cellsize) %>%
    st_as_sf() %>%
    st_intersection(emprise) %>%
    mutate(ID = row_number())  # Ajout d'un identifiant unique à chaque maille
plot(st_geometry(grid))

3 Calcul de l’indice de diversité lithologique

La lithologie est disponible par dalle départementale via la BD Charm-50. Pour le Calvados et l’exécution de ce script, vous pourrez télécharger la donnée ici : https://gitlab.huma-num.fr/tfeuillet/cours/-/raw/main/data/geol_14.gpkg?ref_type=heads.

# Import de la donnée
litho <- st_read("data/geol_14.gpkg", quiet = TRUE)

# Intersection entre les mailles et la lithologie
  litho_grid <- st_intersection(litho, grid) %>%
    mutate(surface_intersection = as.numeric(st_area(.)))  # Surface après intersection
  
  # Agrégation par maille et lithologie
  part_surface_litho <- litho_grid %>%
    st_drop_geometry() %>%
    group_by(ID, CODE) %>% # CODE = ID de chaque couche lithologique
    summarise(surface_litho = sum(surface_intersection), .groups = "drop")
  
  # Calcul des pourcentages de chaque lithologie par maille
  part_litho <- part_surface_litho %>%
    group_by(ID) %>%
    mutate(surface_totale = sum(surface_litho),
           part_litho = (surface_litho / surface_totale) * 100) %>%
    ungroup() %>%
    select(ID, CODE, part_litho) %>%
    pivot_wider(names_from = CODE, values_from = part_litho, values_fill = 0)
  
  # Calcul de l'indice de Shannon
  part_litho$shannon_litho <- diversity(part_litho %>% select(-ID), index = "shannon")
  
  # Jointure des résultats avec la grille
  grid <- grid %>% left_join(part_litho, by = "ID") %>% 
    mutate(across(where(is.numeric), ~ replace_na(.x, 0)))
  
  # Visualisation cartographique
  tmap_mode("view")
  tm_shape(grid) +
    tm_basemap("OpenTopoMap") + 
    tm_polygons(col = "shannon_litho", style = "quantile", n = 6, palette = "viridis", title = paste("Diversité lithologique"))