Analyse spatiale avec Python et GeoPandas

GeoPandas étend la célèbre bibliothèque pandas avec une colonne de géométries. Résultat : tout ce que vous faites avec un tableau (filtrer, grouper, joindre) fonctionne aussi avec des données géographiques, plus les opérations spatiales. Idéal pour automatiser des traitements que l'on ferait à la main dans QGIS.

Installation

python -m venv .venv
# Windows : .venv\Scripts\activate    —   macOS/Linux : source .venv/bin/activate
pip install geopandas matplotlib mapclassify folium

Lire et explorer

import geopandas as gpd

communes = gpd.read_file("communes.gpkg")   # GeoPackage, GeoJSON, Shapefile, zip…
print(communes.shape)       # (nb_lignes, nb_colonnes)
print(communes.crs)         # système de coordonnées, ex. EPSG:2154
print(communes.head())
communes.plot(figsize=(8, 8), edgecolor="white", linewidth=0.2)

Un GeoDataFrame est un DataFrame classique avec une colonne spéciale geometry contenant des objets Shapely.

Reprojeter

# Toujours travailler en projeté pour mesurer !
communes = communes.to_crs(epsg=2154)

communes["surface_km2"] = communes.area / 1e6
communes["densite"] = communes["population"] / communes["surface_km2"]

# Définir un CRS manquant (ne transforme PAS les coordonnées)
# gdf = gdf.set_crs(epsg=4326)

Créer des géométries depuis un CSV

import pandas as pd

df = pd.read_csv("stations.csv")   # colonnes lon, lat, nom…
stations = gpd.GeoDataFrame(
    df,
    geometry=gpd.points_from_xy(df["lon"], df["lat"]),
    crs="EPSG:4326",
).to_crs(2154)

Jointure spatiale

Associer à chaque station la commune qui la contient, puis compter :

jointure = gpd.sjoin(stations, communes[["nom", "geometry"]], predicate="within", how="left")

nb_par_commune = (
    jointure.groupby("nom").size()
    .rename("nb_stations")
    .reset_index()
)
communes = communes.merge(nb_par_commune, on="nom", how="left").fillna({"nb_stations": 0})

Zones tampons et superposition

# Zone de 500 m autour de chaque station
tampons = stations.copy()
tampons["geometry"] = stations.buffer(500)

# Fusionner les tampons qui se chevauchent
zone_desservie = tampons.dissolve()

# Part de chaque commune couverte
inter = gpd.overlay(communes, zone_desservie[["geometry"]], how="intersection")
inter["surf_couverte"] = inter.area
couverture = inter.groupby("nom")["surf_couverte"].sum() / communes.set_index("nom").area

Cartographier

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(9, 9))
communes.plot(
    column="densite", scheme="NaturalBreaks", k=6,
    cmap="viridis", legend=True, ax=ax,
    legend_kwds={"title": "Habitants / km²", "loc": "lower left"},
)
ax.set_axis_off()
ax.set_title("Densité de population", fontsize=16)
plt.savefig("densite.png", dpi=200, bbox_inches="tight")

# Carte web interactive en une ligne (Leaflet via folium)
communes.explore(column="densite", scheme="NaturalBreaks", k=6, cmap="viridis").save("carte.html")

Exporter

communes.to_file("resultat.gpkg", layer="communes_densite", driver="GPKG")
communes.to_crs(4326).to_file("resultat.geojson", driver="GeoJSON")
communes.to_parquet("resultat.parquet")   # GeoParquet, très rapide

Démonstration

Les opérations de ce tutoriel rejouées pas à pas dans le navigateur avec Turf.js, l'équivalent JavaScript de Shapely / GeoPandas : vraies communes, stations fictives.

À essayer : suivez les étapes 1 à 5. La jointure spatiale compte les stations par commune, puis viennent les tampons de 500 m, leur fusion en une zone desservie, et enfin la part de chaque commune couverte, calculée par intersection. Le code Python correspondant s'affiche à chaque étape ; cliquez une commune pour lire son résultat.