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

Auteur

Thierry Feuillet (Université de Caen, UMR IDEES)

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

4 Calcul de l’indice de diversité hydrologique

Nous allons capture la diverité hydrologique à travers deux métriques : la diversité hiérarchique des affluents (selon les ordres de Strahler) et la diversité de la nature des plans d’eau.

4.1 Diversité des ordres de Strahler

La classification de Strahler consiste à hiérarchiser l’ensemble des drains d’un réseau hydrographique de manière à ce que l’ordre soit proportionnel au nombre de confluences. Ainsi, un drain d’ordre 1 signifie signifie qu’il n’a pas d’affluent, un drain d’ordre 2 qu’il en a 1, etc.

Classification des drains d’un réseau hydrographique selon Strahler (source wikipédia)

La manière la plus simple de produire cette information est d’utiliser un algorithme de SAGA GIS disponible via QGIS. Cet algorithme calcule les ordres de Strahler d’un réseau uniquement sur la base d’un MNT en entrée.

Il convient donc dans un premier temps de télécharger les dalles de MNT sur votre emprise, depuis le site de l’IGN (ici la BD ALTI® 25 m : https://geoservices.ign.fr/bdalti). Une fois les fichiers .asc décompressés dans un dossier (par exemple “DATA/MNT”), vous pouvez les fusionner et les croper selon votre emprise de cette manière (utilisation du package terra) :

# Code non reprodictible (car les dalles MNT sont lourdes), simplement pour illustrer la procédure
f <-list.files(path = "data/MNT", pattern = '.*\\.asc$', full.names = TRUE) # Ajustez votre path
r <- lapply(f, rast) 
mnt <- do.call("merge",r)
mnt <- crop(mnt, emprise, mask = TRUE)

# Puis il faut exporter ce mnt
writeRaster(mnt, "data/MNT/MNT14.tif")

Le fichier en résultant peut être téléchargé ici :

mnt <- rast("data/MNT/MNT14.tif")
plot(mnt)

Puis importer le MNT dans QGIS et exécuter l’algo Strahler de SAGA GIS (disponible dans la boite à traitements après avoir installé SAGA Next Gen dans les extensions).

Capture d’écran de l’algo de SAGA (dans QGIS) permettant de calculer les ordres de Strahler)

Exprtez le raster résultant du calcul (par exemple : “strahler_14.tif”), puis importez-le dans R :

strahler <- rast("data/strahler_14.tif")

On peut maintenant passer au calcul de l’indice de Shannon des ordres de Strahler. Le raisonnement est le même que pour la lithologie, sauf qu’il s’agit ici de coupler un raster avec un vecteur, et non deux vecteurs. Ce couplage est basé sur la fonction terra::extract().

 # Conversion de la grille en SpatVector pour extraction des valeurs raster
  grid_vect <- vect(grid)

  # Extraction des valeurs du raster de lithologie par maille
  strahler_values <- extract(strahler, grid_vect)
  
  # Calcul des proportions des classes lithologiques par maille
  strahler_counts <- strahler_values %>%
    group_by(ID, strahler_14) %>%  # `litho` est la valeur du raster (ID de la lithologie)
    summarise(n_pixels = n(), .groups = "drop") %>%  # Nombre de pixels de chaque litho
    mutate(total_pixels = sum(n_pixels)) %>%
    mutate(part_strahler = (n_pixels / total_pixels) * 100) %>%
    select(ID, strahler_14, part_strahler)
  
  # Restructuration en format large (pivot_wider)
  part_strahler <- strahler_counts %>%
    pivot_wider(names_from = strahler_14, values_from = part_strahler, values_fill = 0)
  
  # Calcul de l'indice de Shannon
  part_strahler$shannon_strahler <- diversity(part_strahler %>% select(-ID), index = "shannon")
  
  # Jointure avec la grille spatiale
  grid <- grid %>%
    left_join(part_strahler, by = "ID") %>%
    mutate(across(where(is.numeric), ~ replace_na(.x, 0)))
  
tm_shape(grid) +
    tm_basemap("OpenTopoMap") + 
    tm_polygons(col = "shannon_strahler", style = "quantile", n = 6, palette = "viridis", title = paste("Diversité de Strahler"))
── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
`tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style', 'n', 'palette' (rename to 'values') to
  'tm_scale_intervals(<HERE>)'
[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'

4.2 Diversité des plans d’eau

Dans la BD TOPAGE® ou dans la couche HYDROGRAPHIE de la BD TOPO®, nous disposons d’une couche PLAN_D_EAU.shp dont un attribut renseigne la nature des plans d’eau (mare, estuaire, marais, etc.). Nous allons utiliser ces catégories pour calculer un nouvel indice de Shannon.

eau <- st_read("data/PLAN_D_EAU.shp", quiet = TRUE) %>% rename(ID_EAU = ID) 
# Renommer la variable ID pour éviter la confusion avec l'ID de grid

# Intersection entre les mailles et les plans d'eau
  eau_grid <- st_intersection(eau, grid) %>%
    mutate(surface_intersection = as.numeric(st_area(.)))  # Surface après intersection
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
  # Agrégation par maille et plans d'eau
  part_surface_eau <- eau_grid %>%
    st_drop_geometry() %>%
    group_by(ID, NATURE) %>% # NATURE = ID de chaque type de plan d'eau
    summarise(surface_eau = sum(surface_intersection), .groups = "drop")
  
  # Calcul des pourcentages de chaque plan d'eau par maille
  part_eau <- part_surface_eau %>%
    group_by(ID) %>%
    mutate(surface_totale = sum(surface_eau),
           part_eau = (surface_eau / surface_totale) * 100) %>%
    ungroup() %>%
    select(ID, NATURE, part_eau) %>%
    pivot_wider(names_from = NATURE, values_from = part_eau, values_fill = 0)
  
  # Calcul de l'indice de Shannon
  part_eau$shannon_eau <- diversity(part_eau %>% select(-ID), index = "shannon")
  
  # Jointure des résultats avec la grille
  grid <- grid %>% left_join(part_eau, by = "ID") %>% 
    mutate(across(where(is.numeric), ~ replace_na(.x, 0)))

4.3 Diversité hydrologique totale

Il ne reste plus qu’à sommer les deux indices pour produire l’indice de Shannon pour l’hydrologie :

grid$shannon_hydro <- grid$shannon_eau + grid$shannon_strahler

tm_shape(grid) +
    tm_basemap("OpenTopoMap") + 
    tm_polygons(col = "shannon_hydro", style = "quantile", n = 6, palette = "viridis", title = paste("Diversité hydrologique"))
── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
`tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style', 'n', 'palette' (rename to 'values') to
  'tm_scale_intervals(<HERE>)'
[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'

5 Géodiversité totale

grid$geodiversite <- grid$shannon_litho + grid$shannon_hydro

tm_shape(grid) +
    tm_basemap("OpenTopoMap") + 
    tm_polygons(col = "geodiversite", style = "quantile", n = 6, palette = "viridis", title = paste("Géodiversité"))
── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
`tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style', 'n', 'palette' (rename to 'values') to
  'tm_scale_intervals(<HERE>)'
[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'