Aller au contenu

5- Indicateurs de diversité écologique

Danger

Ce TP est évalué. Votre rendu, en binôme, est consistué :

  • du code déposé sur un GIT
  • d'un rapport illustrant l'ensemble des étapes du projet.

Vous devez renseigner d'une manière ou d'une autre les packages nécessaires pour construire l'environnement virtuel nécessaire pour exécuter vos scripts. Vous pouvez par exemple le documenter dans le README.md de votre dépôt, ou bien fournir un fichier requirements.txt qui spécifie ces dépendances.

La directive cadre sur l'eau (DCE) est une directive Européenne qui impose aux états Européens des objectifs pour améliorer la qualité des masses d'eau de surface et souterraine. Pour suivre l'atteinte de ces objectifs, elle impose de réaliser un suivi d'un certain nombre de facteurs, qu'ils soient de nature physico-chimiques (température, pH, salinité, azote, nitrates, etc...) ou liés au recensement de la faune et de la flore. Concernant la faune, plusieurs compartiments sont considérés :

  • les diatomées benthiques qui sont des petites algues unicellulaires,
  • les poissons,
  • les macro-inverterbrés (visibles à l'oeil nu)

Dans la suite de ce sujet, nous allons nous intéresser aux diatomées benthiques en exploitant les mesures mises à disposition par le réseau de surveillance et accessible via l'API Hubeau. Ces données sont collectées dans plusieurs stations réparties sur tout le territoire Français. Pour disposer d'une vision d'ensemble des mesures, nous aggrégerons dans un premier temps ces observations à l'échelle des départements et des régions.

Les codes de région/département associés aux stations de mesure sont facilement accessibles dans les données, mais, d'un point de vue écologique ainsi que dans le cadre de la DCE, l'aggrégation au niveau des départements/régions n'est pas la plus pertinente. L'aggrégation de la plus pertinente se fait au niveau des hydro-éco-régions qui sont des territoires géographiquement définis qui regroupent des cours d'eau aux caractéristiques physiques et biologiques similaires. Je propopose de reporter l'aggrégation au niveau des hydro-écrorégions en fin de sujet.

Je tiens à remercier Philippe Le Noac'h pour son aide pour la compréhension des données Naïades et ses conseils.

Code préléminaire : tracé d'un fond de carte avec geopandas

A plusieurs reprises dans ce sujet, nous allons avoir besoin de tracer des fonds de carte, avec le tracé des limites de département, de région, etc... Pour ce faire, nous allons utiliser la librairie geopandas. Avec cette librairie, il est possible de charger des fichiers de description de polygones sous la forme de shapefile.

L'administration Française met à disposition des fichiers shapefile, issus du projet OpenStreetMap pour les départements et les régions :

Par ailleurs, l'INSEE identifie chaque région Française par un code. Ce même code est celui utilisé par Naïades lorsqu'il faut identifier la région dans laquelle se trouve une station de mesure.

Code INSEE Nom de la région
1 Guadeloupe
2 Martinique
3 Guyane
4 La Réunion
6 Mayotte
11 Île-de-France
24 Centre-Val de Loire
27 Bourgogne-Franche-Comté
28 Normandie
32 Hauts-de-France
44 Grand Est
52 Pays de la Loire
53 Bretagne
75 Nouvelle-Aquitaine
76 Occitanie
84 Auvergne-Rhône-Alpes
93 Provence-Alpes-Côte d'Azur
94 Corse

Je vous invite à récupérer et extraire les fichiers pour les départements et les régions. Le code ci-dessous permet alors de générer une carte à partir des frontières de département. Le code a été testé avec pandas==2.3.3, geopandas==1.1.1 et matplotlib==3.10.7.

# coding: utf-8

# External imports
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt

gdf = gpd.read_file("departements/departements-20180101.shp")
print(gdf.columns)
print(gdf.head())

# Un exemple de dictionnaire associant une valeur à plusieurs codes de départements
values = {
    "75": 10, "13": 7, "33": 4, "69": 9, "59": 6,
}

# Construction d'une dataframe pandas à partir de ces valeurs
df_values = pd.DataFrame(list(values.items()), columns=["code_insee", "value"])

# On fusionne les deux dataframes sur la colonne commune "code_insee"
gdf_merged = gdf.merge(df_values, on="code_insee", how="left")

# On utilise matplotlib pour tracer le fond de carte
fig, ax = plt.subplots(figsize=(8, 10))

gdf_merged.plot(
    column="value",
    cmap="OrRd",
    linewidth=0.5,
    edgecolor="black",
    legend=True,
    ax=ax
)

# On calcule le centroïde de chaque département (centre géométrique)
# pour y placer les noms de département
gdf_merged["centroid"] = gdf_merged.geometry.centroid

for _, row in gdf_merged.iterrows():
    if not pd.isna(row["value"]):
        x, y = row["centroid"].x, row["centroid"].y
        label = row["nom"]  
        ax.text(
            x, y,
            label,
            fontsize=12,
            ha='center', va='center',
            color="black"
        )

ax.set_title("Quelques non départements de la région Grand-Est", fontsize=14)
ax.axis("off")

plt.tight_layout()
plt.savefig('illustration_geopandas.png')
print("Image saved to illustration_geopandas.png")

Ce qui produit la carte :

Illustration d'un fond de carte de quelques départements de la région Grand-Est

Le code ci-dessus fonctionne également avec le fichier shapefile des régions :

Illustration d'un fond de carte de quelques régions de France Métropolitaine

Exploration des données

On souhaite calculer et analyser des indicateurs de biodiversité au niveau de la France métropolitaine. Pour ce faire, nous allons utiliser l'API Hubeau. Dans le précédent TP, vous avez commencer à la prendre en main, en particulier au travers du module python hubeau.py qui permet d'interfacer du code python avec quelques endpoints de l'API Hub'Eau.

Question

Mettez en place votre projet, ce qui consiste à :

  • créer un dépôt GIT et inviter vos profs encadrants
  • télécharger le script hubeau.py permettant d'accéder à l'API depuis python
  • adder/comiter/pusher ce code sur votre dépôt GIT

L'analyse de la biodiversité va exploiter l'API Hydrobiologie https://hubeau.eaufrance.fr/page/api-hydrobiologie, en particulier le endpoint taxons qui permet d'accéder aux recensements des poissons, macro-invertébrés, diatomées et macrophytes de 1971 à 2018 avec plus de 3M de taxons répertoriés.

Sur le site de l'API Hydrobiologie https://hubeau.eaufrance.fr/page/api-hydrobiologie, il est possible :

  1. de tester l'API directement dans le navigateur
  2. de voir le format de la réponse, tout en bas de la page, dans la section "Models / Résultat d'une requête ..."

Expérimentation de l'API en ligne

Vous pouvez explorer quelques requêtages de l'API via le portail Web https://hubeau.eaufrance.fr/page/api-hydrobiologie.

Question

Via le portail Web d'accès à l'API Hydrobiologie, pour le endpoint taxons, faites une requête permettant d'obtenir :

  • les effectifs (code_type_resultat=1)
  • des diatomées benthiques (code_support=10)
  • pour la région Grand-Est (code_region et se réferrer à la table en haut de page)
  • pour le mois de Juillet 2015 (date_debut_prelevement et date_fin_prelevement)
  • pour les mesures de bonne qualité (code_qualification=1)

Dans le résultat de la requête :

  • libelle_appel_taxon et code_appel_taxon identifient le taxon dénombré
  • resultat_taxon corresponds à l'abondance (nombre d'individus observés)
  • code_prelevement corresponds à un identifiant unique du prélèvement, associé à une station de mesure pour une date donnée
  • code_station_hydrobio, libelle_station_hydrobio corresponds à la station de mesure
  • date_prelevement corresponds à la date de prélèvement

Combien de mesures totales sont disponibles ?

Récupération de données

Le code du script hubeau.py vous permet d'accéder aux données d'Hydrobiologie, en spécifiant les différents paramètres tels qu'indiqués sur la page de l'API hydrobiologie. Attention, l'année est gérée de manière un peu singulière dans le script hubeau.py et c'est le seul paramètre qui ne correspond pas directement à un paramètre de l'API. Voyons quelques exemples :

import hubeau

hydrobio_client = hubeau.Hydrobiologie()

# Récupération des données de l'année 2010 pour les départements
# de Moselle et Meurthe-et-Moselle
df = hydrobio_client.taxons(year=2010, 
                            code_departement=[54, 57]) 
print(df)

# On peut limiter les champs retournés par la requête 
# grâce au paramètre fields
df = hydrobio_client.taxons(year=2010, 
                            code_departement=[54, 57],
                            fields=["date_prelevement", "resultat_taxon"])

Danger

Il peut arriver que les appels à l'API Hubeau restent bloqués. Dans ce cas, il suffit d'interrompre l'exécution du script et de le relancer.

Question

A l'aide du script hubeau.py (mais dans un script à part !), écrivez une fonction permettant de requêter l'API Hubeau hydrobiologie pour les données. Je vous suggère d'écrire une fonction test_api() dans laquelle vous implémenterez la requête suivante :

  • sur les diatomées benthiques (voir en particulier le champ code_support documenté dans le document sur les codes support SANDRE,
  • pour l'année 2010,
  • pour les départements des Vosges (88), Meurthe et Moselle (54), Moselle (57), Meuse (55) (voir code_departement),
  • pour les comptages d'effectifs (voir code_type_resultat),
  • uniquement pour les relevés qui ont une qualification "correcte" (code_qualification=1),
  • en ne conservant que les colonnes du :
    • code de département (code_departement),
    • code de la station de prélèvement (code_station_hydrobio),
    • les coordonnées de la station (longitude, latitude),
    • identifiant unique du prélèvement (code_prelevement)
    • date de prélèvement (date_prelevement).
    • résultat (resultat_taxon),
    • libellé du taxon (libelle_appel_taxon),

Maintenant que vous avez les données, il faut également ajuster les types des différentes colonnes.

Question

Commencez par observer le type des différentes colonnes. Rappelez vous l'attribute dtypes de la DataFrame.

Ecrivez et appelez une fonction def fix_datatypes(df) pour ajuster les types des différentes colonnes (pd.as_numeric, astype(), astype('category'), pd.to_datetime, ...)

Je vous suggère de partir du modèle : 

def fix_datatypes(df):
    """
    Fix the datatypes of our df columns
    It operates in place
    """

    # Fix numeric columns
    num_columns = ["longitude", "latitude"]
    for col in num_columns:
        if not col in df:
            continue
        df[col] = pd.to_numeric(df[col])

La garde conditionnelle if col in df permet de s'assurer que la conversion ne sera appliquée que si la colonne est effectivement présente dans la dataframe. On sera amené, par la suite, à appliquer cette fonction en l'absence de certaines colonnes.

Question

Combien y'a t'il de stations qui ont réalisé des prélèvements de diatomées benthiques par département du Grand-Est sur l'année 2010 ? Le groupement (voir TP 4 - Grouper et agréger des données) peut être très utile.

En auriez vous trouvé \(17\) en Meuse ?

Visualisation des stations hydrobiologiques

On se propose de tracer la localisation des stations hydrobiologiques. La localisation des stations est accessible via l'API Hydrobiologie, le endpoint stations_hydrobio.

Je vous conseille d'écrire une fonction spécifique, par exemple 

def stations_hydrobio():
    ....

et de l'invoquer en fin de script dans le bloc 

if __name__ == '__main__':
    stations_hydrobio()

Question

Vous devez tracer la localisation des stations hydrobiologique pour l'année 2015 pour la région Grand-Est.

L'accès au endpoint stations_hydrobio se fait de manière similaire au endpoint taxons, à ceci près qu'il n'y a pas d'année à spécifier.

hydrobio_client = hubeau.Hydrobiologie()
df = hydrobio_client.stations_hydrobio(code_region=...,
                            fields=[...])

Vous pouvez récupérer les longitudes et latitudes des stations et les exploiter comme nous l'avons fait dans le TP4 : Localiser des données

L'image ci-dessous représente la position des stations où des prélèvements hydrobiologiques ont eu lieu dans la région Normandie.

Localisation des stations de mesure hydrobiologique dans la région Grand-Est

Effort de prospection

Pour évaluer dans quelle mesure la DCE est mise en oeuvre, nous allons compter le nombre d'analyse sur les diatomées réalisées par département, ce que nous appellerons l'effort de propsection.

L'objectif est de produire une carte de la France Métropolitaine, sur laquelle figure le nombre de campagnes de mesures réalisées sur chaque département.

Conseil de mise en oeuvre

Je vous conseille de répondre à cette partie par une fonction python dédiée qui va collecter les données, les filtrer et réaliser les tracés.

L'objectif étant de mesurer un effort de prospection par département, je vous conseille également de déléguer la mesure de l'effort de prospection à une fonction dédiée. 

def get_effort_prospection(departement):
    ....

def effort_prospection():

    # Collecte des efforts de prospection
    for departement in ... :
        effort = get_effort_prospection(departement)
        ...

    # Tracé de l'effort de propsection par département
    ...

if __name__ == "__main__":
    effort_prospection()

Voyons tout ça par étape.

Aggrégats des mesures par département

Pour commencer, nous travaillons sur la récolte de l'effort de prospection par la fonction get_effort_prospection. Les prélèvements sont identifiés de manière unique par un code de prélèvement.

Cela se déroule de la manière suivante :

  • vous devez collecter les prélèvements d'abondance (code_type_resultat=1) des diatomées benthiques (code_support=10), de bonne qualité (code_qualification=1)
  • vous devez fixer les types des colonnes (numérique, catégorielle, etc...)
  • regrouper et aggréger les prélèvements par département (voir TP 4 - Grouper et agréger des données).

Danger

Attention, les abondances des différentes diatomées sont renseignées par des entrées différentes dans votre dataframe : une ligne par diatomée. En d'autres termes, il y a des doublons de code de prélèvements.

Question

Complétez votre fonction get_effort_prospection pour retourner une dataframe dénombrant le nombre de prélèvements réalisés par département sur une année donnée.

Tracer de l'effort de prospection

On va maintenant passer au tracé en utilisant le code introduit en début de sujet.

Question

Tracer l'effort de propsection pour l'année 2020 pour les régions Grand-Est, Ile-de-France et Hauts-de-France.

A titre d'illustration, vous trouverez ci-dessous l'effort de prospection pour l'année 2015 pour un certain nombre de départements du nord de la France.

Effort de prospection sur les diatomées benthiques pour l'année 2015

Calculs d'indicateurs de biodiversité

On se propose maintenant de passer au calcul d'indicateurs de biodiversité. Plusieurs indicateurs de qualité d'eau sont construit à partir du dénombrement de taxons dans différents compartiments. Par exemple, l'abondance des diatomées, identifiées au niveau de l'éspèce permet de calculer l'indice biologique diatomée (IBD). L'Indice Poisson Rivière (IPR) est calculé à partir du rencesement des abondances des poissons capturés par pêche électrique (avant réintroduction dans le milieu). L'Indice Macrophyte en rivière repose sur le recensement des taxons de végétaux aquatiques.

Ce sont des indicateurs qui ne sont pas forcément très simples à calculer, d'autant qu'ils nécessitent quelques variables environnementales en plus des dénombrements de taxon. Les variables environnementales sont accessibles par l'API Hubeau mais on va se limiter dans cette partie au calcul d'un indice beaucoup plus simple : la richesse spécifique.

La richesse spécifique est définit comme le nombre de taxons observés dans un compartiment donné. Par exemple, pour les diatomées, on peut calculer le nombre d'éspèces de diatomées différentes dont l'abondance est strictement positive.

Question

Je vous propose de tracer la richesse spécifique, moyennée par département, en se basant sur les mesures d'abondances des diatomées.

Je vous conseille de procéder de la manière suivante :

  1. Collectez les abondances via l'endpoint taxons, en filtrant le bon code support, la bonne qualification et s'assurant de bien disposer des abondances (resultat_taxon), le nom du taxon identifié (libelle_appel_taxon), des départements et des codes de prélèvement,
  2. Ajustez les types des colonnes
  3. Nettoyez les noms des taxons (voir ci-dessous)
  4. Transformez la table grâce à un pivot pour obtenir, pour chaque prélèvement (index) les abondances des taxons dans des colonnes spécifiques à ces taxons
  5. Comptez, pour chaque prélèvement, le nombre de taxon dont l'abondance est strictement positive ce qui permet de calculer la richesse spécifique par prélèvement
  6. Moyennez la richesse au niveau des départements

Les noms des taxons sont parfois suffixés de " abnormal form" qu'il faudra supprimer pour éviter de considérer que "Achnanthes" et "Achnanthes abnormal form" soient deux taxons différents ! Il faudra également supprimer les lignes pour lesquelles le "libelle_appel_taxon" est "Abnormal diatom valve" ou "Appellation de Taxon inconnue".

A titre d'exemple, vous trouverez ci-dessous une représentation de la richesse spécifique pour l'année 2012 :

Illustration de la richesse spécifique moyennée par département pour l'année 2012

Pour aller plus loin : aggrégation au niveau des hydro-éco-regions

Les aggrégations considérées jusqu'à maintenant sont réalisées au niveau des départements ou des régions ce qui n'est pas le plus pertinent en terme d'écologie ni même dans le cadre de la directive DCE. Un niveau d'aggrégation plus pertinent est celui des hydro-éco-régions qui sont des zones géographiques qui regroupent des masses d'eau partageant des propriétés physiques et biologiques.

Vous pouvez récupérer un fichier shapefile définissant ces hydro-éco-régions sur le site https://www.data.gouv.fr/datasets/hydroecoregions-de-niveau-1-her-1/. Tout comme pour les régions et les départements, le shapefile des hydro-éco-régions est chargeable avec geopandas.

Le script ci-dessous permet de charger et tracer les polygones des hydro-éco-régions. Attention, il faut bien s'assurer de disposer du fichier .shp et .shx même si le script ne semble faire référence qu'au fichier .shp.

# coding: utf-8

# External imports
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import numpy as np

def plot_her():
    gdf = gpd.read_file("her/Hydroecoregion1.shp")
    print(gdf.head())

    # Créer une palette de couleurs basée sur le gid
    unique_gids = gdf["gid"].unique()
    cmap = plt.colormaps['tab20']
    gid_to_color = {gid: cmap(i / len(unique_gids)) for i, gid in enumerate(unique_gids)}
    gdf['color'] = gdf['gid'].map(gid_to_color)

    # Tracé des polygones colorés par type de l'éco région
    gdf.plot(linewidth=0.5, 
             edgecolor="black",
             color=gdf["color"],
             )

    # Créer une légende avec les couleurs et les noms des hydro-éco-régions
    legend_patches = []
    for _, row in gdf.iterrows():
        color = gdf[gdf['gid'] == row['gid']]['color'].values[0]
        label = row['NomHER1']
        patch = mpatches.Patch(facecolor=color, edgecolor='black', linewidth=0.5, label=label)
        legend_patches.append(patch)

    plt.legend(handles=legend_patches, loc='center left', bbox_to_anchor=(1, 0.5), fontsize=9)

    plt.title("Hydro-éco-régions de niveau 1", fontsize=14)
    plt.gca().axis("off")
    plt.tight_layout()
    plt.savefig('illustration_her.png')
    print("Image saved to illustration_her.png")


if __name__ == "__main__":
    plot_her()

Ce script produit l'illustration ci-dessous

Représentation des hydro-éco-régions de niveau 1

Vous pouvez maintenant revenir sur vos codes précédents pour ajouter une colonne à vos dataframes indiquant à quelle hydro-éco-région votre station appartient et réaliser vos aggrégations par hydro-éco-régions plutôt que par départements ou régions.

Références