Accueil Écorégions Données et reproductibilité Scripts de données et de cartes

Scripts de données et de cartes

À quoi sert cette page ? Elle publie les scripts qui enrichissent les unités après le découpage : population par zonage de référence, zones climatiques de Köppen, une classification des climats du monde, fragmentation des nappes, puis auto-documentation du géopackage et cartographie. Le lecteur voit comment chaque diagnostic et chaque carte se produisent.

Après le découpage viennent les traitements qui enrichissent les unités et produisent les diagnostics et les cartes. Tous partent du même fichier pivot ecoregions-2027.gpkg et de la couche commune d'ADMIN EXPRESS.

Population par zonage de référence

La même méthode d'affectation, un point représentatif de commune joint à la zone qui le contient, se réutilise pour toute couche de zonage. Une fonction générique traite ainsi les 7 circonscriptions de bassin, puis les 22 hydro-écorégions de niveau 1 et les 114 de niveau 2. Le point représentatif, garanti à l'intérieur du polygone contrairement au centroïde, évite les erreurs d'affectation aux frontières.

#!/usr/bin/env python3
# Population INSEE par zonage : circonscriptions de bassin, HER-1, HER-2.
import geopandas as gpd, pandas as pd
from pathlib import Path
WORK = Path(__file__).resolve().parents[1] / "ecoregions-2027.gpkg"
ADMIN = ".../ADMIN-EXPRESS.gpkg"

com = gpd.read_file(ADMIN, layer="commune")
com["population"] = pd.to_numeric(com["population"], errors="coerce").fillna(0).astype(int)
pts = com.copy(); pts["geometry"] = com.representative_point()
pts = pts"code_insee", "population", "geometry"

def affecter(niveau, cd, nom):
    zone = gpd.read_file(WORK, layer=niveau).to_crs(pts.crs)
    j = gpd.sjoin(pts, zonecd, nom, "geometry", predicate="within", how="left")
    ok = j.dropna(subset=[cd])
    agg = (ok.groupby([cd, nom])
             .agg(population=("population", "sum"),
                  n_communes=("code_insee", "count")).reset_index())
    res = zone.merge(agg, left_on=[cd, nom], right_on=[cd, nom], how="left")
    res["population"] = res["population"].fillna(0).astype(int)
    res["surface_km2"] = (res.geometry.area / 1e6).round(1)
    res["dens_hab_km2"] = (res["population"] / res["surface_km2"]).round(1)
    res.to_file(WORK, layer=f"{niveau}_pop", driver="GPKG")
    return res

affecter("bassins_agences", "NumCircAdminBassin", "NomCircAdminBassin")
affecter("her1", "CdHER1", "NomHER1")
affecter("her2", "CdHER2", "NomHER2")

Zones climatiques Köppen-Geiger

Le climat de chaque écorégion s'appuie sur la classification Köppen-Geiger (Beck, 2023, période 1991-2020, résolution 1 km). Le raster est polygonisé sous GDAL, dissous par classe, découpé sur le masque France métropolitaine et Corse, simplifié à 500 m, puis étiqueté par sa légende climatique.

#!/usr/bin/env python3
# Vectorise Köppen-Geiger (Beck 2023) 1991-2020 sur France métropole + Corse, L93.
from osgeo import gdal, ogr, osr
gdal.UseExceptions()
import geopandas as gpd
from shapely.geometry import shape
from shapely.ops import unary_union

RASTER = "/tmp/kg_fr_l93.tif"; ADE = ".../ADMIN-EXPRESS.gpkg"
OUT = "/tmp/climat_koppen_l93.gpkg"

# 1. polygonisation du raster (valeur de bande -> champ DN)
src = gdal.Open(RASTER); band = src.GetRasterBand(1)
srs = osr.SpatialReference(); srs.ImportFromWkt(src.GetProjection())
mem = ogr.GetDriverByName("MEMORY").CreateDataSource("mem")
lyr = mem.CreateLayer("poly", srs=srs, geom_type=ogr.wkbPolygon)
lyr.CreateField(ogr.FieldDefn("DN", ogr.OFTInteger))
gdal.Polygonize(band, band.GetMaskBand(), lyr, 0, [], callback=None)
import json
feats = [{"DN": f.GetField("DN"),
          "geometry": shape(json.loads(f.GetGeometryRef().ExportToJson()))}
         for f in lyr if f.GetGeometryRef() is not None]
gdf = gpd.GeoDataFrame(feats, crs="EPSG:2154")
gdf = gdf[gdf["DN"] > 0].copy(); gdf["geometry"] = gdf.buffer(0)

# 2. dissolution par classe, découpe France, simplification 500 m
diss = gdf.dissolve(by="DN", as_index=False)
reg = gpd.read_file(ADE, layer="region")
metro = {"11","24","27","28","32","44","52","53","75","76","84","93","94"}
france = unary_union(reg[reg["code_insee"].isin(metro)].geometry.buffer(0).values)
gfr = gpd.GeoDataFrame(geometry=[france], crs="EPSG:2154")
clip = gpd.overlay(diss, gfr, how="intersection", keep_geom_type=True)
clip["geometry"] = clip.simplify(500, preserve_topology=True).buffer(0)
clip = clip.dissolve(by="DN", as_index=False)
clip.to_file(OUT, layer="climat_koppen", driver="GPKG")

Climat des écorégions

Nappes souterraines et fragmentation

Les masses d'eau souterraines DCE sont importées, filtrées sur la métropole et reprojetées, puis croisées avec les circonscriptions de bassin. Chaque intersection nappe-bassin de plus de 1 km² est comptée ; une nappe qui recoupe au moins deux bassins impose une cogestion.

#!/usr/bin/env python3
# MESO DCE cycle 2022-2027 : import + fragmentation par bassin d'agence.
import geopandas as gpd, pandas as pd
SRC_SHP = ".../MasseDEauSouterraine_FRA.shp"; GPKG = ".../ecoregions-2027.gpkg"
METRO_EU = {"FRA","FRB1","FRB2","FRC","FRD","FRE","FRF","FRG","FRH"}

meso = gpd.read_file(SRC_SHP)
meso = meso[meso["CdEuBassin"].isin(METRO_EU)].copy().to_crs(2154)
meso["geometry"] = meso.geometry.buffer(0)
meso = meso.rename(columns={"CdMasseDEa": "code_meso", "NomMasseDE": "nom_meso"})
meso.to_file(GPKG, layer="meso", driver="GPKG")

bassins = gpd.read_file(GPKG, layer="bassins_agences").to_crs(2154)
bassins["geometry"] = bassins.geometry.buffer(0)
bassins = bassins"NomCircAdminBassin", "geometry".rename(
    columns={"NomCircAdminBassin": "nom_bassin"})
inter = gpd.overlay(meso"code_meso", "nom_meso", "geometry", bassins,
                    how="intersection", keep_geom_type=True)
inter["surf_km2"] = inter.geometry.area / 1e6
inter = inter[inter["surf_km2"] > 1.0]

rows = []
for code, sub in inter.groupby("code_meso"):
    bs = sorted(sub["nom_bassin"].unique().tolist())
    rows.append({"code_meso": code, "nom_meso": sub["nom_meso"].iloc[0],
                 "n_bassins_recoupes": len(bs), "liste_bassins": " | ".join(bs),
                 "surf_totale_km2": round(sub["surf_km2"].sum(), 1)})
frag = pd.DataFrame(rows)
n_cheval = int((frag["n_bassins_recoupes"] >= 2).sum())

Auto-documentation et cartographie

Un script rend le GeoPackage communicable : il inscrit une description lisible pour chaque couche dans gpkg_contents, construit la table catalogue_couches avec producteur, millésime et licence, et embarque des styles QGIS par défaut. Les cartes de synthèse et les cartes thématiques (population, forêts, nappes à cheval) partagent une charte cartographique commune sous matplotlib, rendue en PNG à 200 dpi et en SVG, en projection Lambert-93. Le code complet de ces deux étapes est repris tel quel dans la note de scripts du chantier ; la logique en est directe et sans surprise une fois les couches précédentes produites.

Les contrôles de cohérence montrent comment vérifier que tout ce qui précède se tient.

Ce qu'il faut en retenir

Une même méthode d'affectation, un point représentatif de commune joint à la zone qui le contient, se réutilise pour tous les zonages, ce qui rend les comparaisons homogènes. Le climat vient d'un raster de Köppen à haute résolution vectorisé sous GDAL, les nappes d'un croisement avec les bassins. Le géopackage embarque enfin ses propres métadonnées et ses styles, si bien qu'il se relit sans documentation externe.