Exécution dans le cloud : ce notebook peut aussi tourner dans le cloud (Colab, Kaggle, Renku). Avant de lancer les cellules, télécharge tp3.gpkg et meuse.csv depuis le dossier OneDrive du cours et dépose-les dans le répertoire de travail du notebook (Colab : panneau 📁 → Upload ; Kaggle : File → Upload ; Renku : déjà dans le projet). Procédure complète : docs/cloud-badges.md.

TP Python 3 : Géotraitement vectoriel & modèles non-linéaires¶

Ce TP Python est le pendant numérique du TP QGIS 3. Tu vas reproduire les mêmes opérations de géotraitement vectoriel en Python, puis appliquer des modèles de régression non-linéaires sur des données géographiques historiques.

🔗 Lien avec le TP QGIS 3 : les parties 1–3 reproduisent en Python le même résultat que dans QGIS — la zone d’interface habitat-forêt (WUI) du district bernois de Frutigen-Niedersimmental. La partie 4 modélise la pollution des sols (jeu Meuse, Pays-Bas) : régression polynomiale, GLM et GAM avec lisseur spatial, pour prédire la concentration en zinc à partir de variables de terrain et de la position géographique.

Remplacer des valeurs avec .map()

La couche districts ne couvre pas seulement la Suisse : la colonne ICC contient le code du pays de chaque district (CH, AT, DE, FR, IT). Pour rendre ces codes lisibles, on peut les remplacer par le nom du pays à l’aide d’un dictionnaire et de la méthode serie.map(dictionnaire) :

• le dictionnaire associe à chaque clé (ici, un code de pays) une valeur (ici, le nom du pays) ;

• pour chaque ligne de la colonne, .map() cherche la clé correspondante dans le dictionnaire et la remplace par la valeur associée ;

• le résultat est une nouvelle colonne, de même longueur que la colonne de départ. Une valeur absente du dictionnaire devient NaN (valeur manquante).

Au TP2 (section 6.7), .map() recevait une fonction plutôt qu’un dictionnaire : le principe est le même, chaque valeur de la colonne est transformée une à une.

Partie 1 — Sélection et extraction¶

Dans le TP QGIS 3, nous avons commencé par sélectionner les districts bernois de Frutigen-Niedersimmental et Obersimmental-Saanen, les avons fusionnés pour définir le périmètre d’étude, puis avons découpé les autres couches dans ce périmètre.

En Python avec geopandas, les équivalents sont :

1.1 Sélection des districts et fusion (étapes 2a/2b)¶

On sélectionne les districts par expression (équivalent du filtre SQL "NAME" ILIKE '...' dans QGIS), puis on les fusionne en une seule entité (dissolve) pour définir le périmètre global de la zone d’étude.

Exercice 1 : sélectionner et fusionner les districts

Complète la cellule à trous ci-dessous pour sélectionner les deux districts de la zone d’étude, Frutigen-Niedersimmental et Obersimmental-Saanen, puis pour les fusionner en une seule entité study_area.

Pour cela, tu peux utiliser :

• serie.str.contains(motif) teste, pour chaque ligne d’une colonne de texte, si elle contient motif. Elle renvoie une série de True/False, qui sert de filtre entre crochets : gdf[gdf['col'].str.contains('Lausanne')] ne garde que les lignes dont la colonne col contient « Lausanne ». Tu as déjà rencontré cette méthode au TP2, dans une requête .query() (section 2.4).

• Le motif de str.contains est une expression régulière, dans laquelle | signifie « ou » : 'Lausanne|Morges' garde les lignes qui contiennent l’un ou l’autre nom. Enfin, case=False ignore les majuscules et na=False écarte les valeurs manquantes.

• gdf.dissolve() fusionne les entités d’un GeoDataFrame, comme l’outil Regrouper de QGIS. Avec l’argument by='colonne', elle fusionne les entités qui partagent la même valeur dans cette colonne ; sans argument, elle fusionne toutes les entités en une seule.

1.2 Découpage des couches (étape 2c)¶

Une fois le périmètre défini, on découpe toutes les couches avec gpd.clip(), l’équivalent de l’outil Couper (Clip) dans QGIS. Cela garantit que les géotraitements suivants ne portent que sur la zone d’étude.

Partie 2 — Zones tampons (Buffer)¶

Les zones tampons permettent de définir des zones de proximité autour d’entités. Dans QGIS, nous avons utilisé l’outil Zone tampon pour :

1. Zone Densément Bâtie (ZDB) — buffer +75 m sur les bâtiments (regroupé), puis buffer −75 m pour conserver uniquement les zones densément bâties

2. Emprise des routes — buffer de 6 m (= moitié de la largeur de 12 m)

2.1 ZDB — Zone Densément Bâtie (étape 3a)¶

Exercice 2 : effet de la distance de regroupement

La ZDB dépend de la distance choisie : avec 75 m, deux bâtiments distants de moins de 150 m sont regroupés. Recalcule la ZDB avec une autre distance d, en enchaînant sur une seule ligne les deux zones tampons de l’étape 3a (+d regroupée, puis −d). Essaie plusieurs valeurs, par exemple 10, 25, 150 et 300 m : comment évoluent la surface de la ZDB et le nombre de groupes de bâtiments ? Pourquoi ?

Pour cela, tu peux utiliser :

• geom.buffer(distance) renvoie la zone tampon d’une géométrie, ou de chaque géométrie d’une GeoSeries. La distance s’exprime dans l’unité du système de coordonnées, ici le mètre. Une distance négative rétrécit la géométrie au lieu de l’agrandir.

• serie.union_all() regroupe toutes les géométries d’une GeoSeries en une seule géométrie, comme l’option Dissoudre le résultat de l’outil Zone tampon dans QGIS.

• Le chaînage de méthodes : chaque méthode s’applique au résultat de la précédente, de gauche à droite. Ainsi, a.methode1(x).methode2(y) applique methode1 à a, puis methode2 au résultat obtenu.

Partie 3 — Superposition (Union, Différence, Intersection)¶

Ces outils permettent de combiner plusieurs couches pour créer de nouvelles entités. C’est l’étape centrale du TP pour produire la WUI (Wildland-Urban Interface).

3.1 Zone Urbaine (ZU) = ZDB ∪ routes (étape 4a)¶

3.2 Zone à risque = Buffer(ZU, 80 m) − ZU (étape 4b)¶

3.3 WUI = zone à risque ∩ Forêt (étape 4c)¶

Exercice 3 : construire la WUI

Traduis en code les trois opérations des titres 3.1 à 3.3 ci-dessus. Les géométries de départ sont zdb_geom (ZDB), roads_buf_geom (routes bufferisées) et foret_union (forêt regroupée). Une fois la carte finale affichée, compare-la avec ta couche WUI du TP QGIS.

Pour cela, tu peux utiliser les opérations de superposition de shapely. Elles s’appliquent à une géométrie a, prennent une autre géométrie b en argument, et renvoient une nouvelle géométrie, que tu peux stocker dans une variable pour l’utiliser à l’étape suivante :

• a.union(b) renvoie la surface couverte par a ou par b (a ∪ b).

• a.difference(b) renvoie la partie de a qui n’est pas dans b (a − b). L’ordre compte : b.difference(a) donne un tout autre résultat.

• a.intersection(b) renvoie la partie commune à a et à b (a ∩ b).

• a.buffer(distance) renvoie la zone tampon de a, comme en partie 2.

Partie 4 — Régression non-linéaire : polynômes, GLM et GAM¶

Pourquoi aller au-delà de la régression linéaire ?¶

La régression linéaire (OLS) suppose que l’effet de chaque variable explicative est constant et proportionnel, et que la variable cible peut prendre n’importe quelle valeur réelle. Ces deux hypothèses sont souvent trop restrictives pour des données réelles, en particulier géographiques : l’effet d’une variable peut se courber (s’atténuer, s’inverser), et une variable physiquement bornée (une concentration, un pourcentage, un compte) ne se prête pas à une droite qui peut, en théorie, prendre n’importe quelle valeur.

Cette partie relâche ces hypothèses progressivement, en gardant le même jeu de données (Meuse) et la même variable cible (zinc) comme fil conducteur :

Ces fonctions lisses (GAM) sont approximées par des B-splines (polynômes par morceaux raccordés) avec une pénalité de lissage pour éviter le surapprentissage.

Les données : la Meuse (pollution des sols, Pays-Bas)¶

Le jeu Meuse est un classique de la géostatistique : 155 sites de sol échantillonnés dans la plaine inondable de la rivière Meuse, près de Stein (Pays-Bas). Deux sites n’ont pas de mesure de matière organique et sont écartés : 153 sites sont utilisés pour la modélisation. Les crues déposent des sédiments chargés en métaux lourds ; on mesure notamment la concentration en zinc.

Question : peut-on expliquer la concentration en zinc à partir de ces variables de terrain ? La position géographique (proximité de la rivière, effet des crues) joue-t-elle un rôle propre, au-delà des variables mesurées ?

4.1 Exploration : carte et corrélations¶

Avant de modéliser, on visualise la distribution spatiale du zinc et les corrélations entre variables. La carte révèle d’emblée un fort gradient géographique, les concentrations les plus élevées longent la rivière, qui motivera plus loin l’ajout d’un lisseur spatial.

4.2 Référence : régression linéaire (OLS)¶

On commence par un modèle de référence — la régression linéaire multiple OLS — qui suppose une relation linéaire entre chaque variable explicative et le zinc. Son R² sera la baseline à améliorer.

4.3 Régression polynomiale¶

La régression polynomiale remplace la droite $\hat{y} = \beta_0 + \beta_1 x$ par un polynôme

Malgré son nom, c’est toujours un modèle linéaire au sens statistique : linéaire par rapport aux paramètres $\beta_j$ (on l’ajuste avec np.polyfit, ou une régression OLS sur les colonnes $x, x^2, \ldots, x^d$). Seule la relation entre $x$ et $y$ devient non-linéaire.

Le degré $d$ contrôle la flexibilité. Ici, la relation zinc ~ dist est réellement courbée : le zinc décroît vite près de la rivière puis se stabilise (décroissance de type distance-décroissance). C’est le cas idéal pour voir les deux faces du polynôme :

• une droite (degré 1) sous-ajuste — elle rate la courbure et laisse une structure nette dans les résidus ;

• un polynôme de degré modéré (2–3) épouse la courbure et améliore la prédiction hors-échantillon — c’est le vrai intérêt du polynôme, celui qu’une droite ne peut pas offrir ;

• un degré élevé (8) sur-ajuste : il oscille pour coller au bruit et généralise moins bien qu’une droite.

On l’illustre en ajustant sur un petit échantillon de 40 sites (pour rendre le surapprentissage visible) et en évaluant sur les 113 restants, jamais vus lors de l’ajustement.

Exercice 4 : ajuster un polynôme

La cellule précédente a tiré 40 sites d’entraînement (x_tr, y_tr) ; les autres sites (x_te, y_te) ne servent qu’à l’évaluation. Ajuste un polynôme du degré de ton choix sur les sites d’entraînement, puis calcule son R² sur ces mêmes sites et sur les sites jamais vus. Essaie plusieurs degrés : lequel prédit le mieux les sites jamais vus ? La cellule suivante compare ensuite les degrés 1, 2, 3 et 8.

Pour cela, tu peux utiliser trois fonctions déjà rencontrées au TP2, dans la section 6.4 sur le surapprentissage :

• np.polyfit(x, y, deg) ajuste par moindres carrés un polynôme de degré deg aux points (x, y). Elle renvoie ses coefficients, du plus haut degré au terme constant.

• np.polyval(coeffs, x) évalue le polynôme défini par coeffs en chaque valeur de x. Le résultat, ce sont les prédictions du modèle en ces points.

• r2_score(y_vrai, y_predit) compare les valeurs observées aux prédictions et renvoie le R². Pour évaluer le polynôme sur un groupe de sites, il faut donc d’abord calculer ses prédictions pour ces sites.

4.4 GLM — un modèle adapté à une grandeur positive¶

zinc est une concentration : elle est strictement positive par nature. Rien ne garantit pourtant qu’une droite OLS y reste positive. Pire, ici elle ne l’est même pas sur les données observées : la droite zinc ~ dist croise zéro dès dist ≈ 0.63 et prédit donc des concentrations négatives pour les sites les plus éloignés de la rivière — ce qui n’a aucun sens physique.

Un GLM (Generalized Linear Model) règle ce problème en ne modélisant plus directement le zinc, mais une transformation de celui-ci. Avec le lien log et une loi Gamma (adaptée à une grandeur positive continue), on modélise :

Le membre de gauche — le log de la moyenne prédite — peut prendre n’importe quelle valeur réelle, donc la combinaison linéaire à droite s’ajuste comme dans un OLS. Mais dès qu’on inverse le lien pour revenir à l’échelle du zinc ($\exp(\cdot)$), le résultat est mathématiquement garanti positif, quels que soient les prédicteurs.

Exercice 5 : ajuster le GLM

Ajuste le GLM décrit ci-dessus : mêmes variables que l’OLS de la section 4.2, loi Gamma et lien logarithmique. La cellule calcule ensuite son R² pour le comparer à celui de l’OLS.

Pour cela, tu peux utiliser :

• smf.glm(formule, data=..., family=...) s’utilise comme smf.ols (section 4.2). La formule 'cible ~ variable1 + variable2' décrit le modèle, data donne le tableau de données, et .fit() lance l’ajustement. L’argument supplémentaire family précise la loi de la variable cible et la fonction de lien.

• sm.families.Loi(sm.families.links.Lien()) construit cette famille. Loi est le nom de la loi, par exemple Gaussian, Poisson ou Gamma (voir la liste des lois). Lien est le nom de la fonction de lien, par exemple Identity, Log ou Logit (voir la liste des liens). Les deux noms commencent par une majuscule.

4.5 GAM thématique — splines sur les variables de terrain¶

Un GAM remplace la droite de régression de chaque prédicteur par une courbe lisse (B-spline). Cela permet de capturer des relations non-linéaires — par exemple un effet de la distance qui s’atténue au-delà d’un certain seuil.

La pénalité de lissage (alpha) contrôle la souplesse de la courbe : plus alpha est grand, plus la courbe se rapproche d’une droite. En statsmodels, le GAM utilise des B-splines pénalisées (P-splines).

Le graphe d’effets partiels (plot_partial) montre la contribution de chaque variable à la prédiction, toutes choses égales par ailleurs.

4.6 GAM spatial, ajouter un lisseur géographique¶

La carte du zinc (§4.1) montrait un gradient spatial marqué (concentrations élevées le long de la rivière). Une partie de ce gradient n’est pas capturée par les seules variables de terrain.

On l’intègre en ajoutant un lisseur spatial 2D : $f(x, y)$ — une surface lisse sur les coordonnées des sites. Cela revient à modéliser l’autocorrélation spatiale résiduelle.

Le lisseur spatial est ajusté avec une pénalité plus faible (alpha plus petit) pour lui permettre de capturer des gradients larges.

4.7 Validation croisée K-fold¶

Les R² vus jusqu’ici (OLS, polynomiale, GLM, GAM) sont pour la plupart calculés sur les mêmes données que celles utilisées pour l’ajustement : ils sont optimistes. Un modèle plus complexe peut avoir un R² d’entraînement élevé tout en généralisant mal à de nouvelles observations (surapprentissage) — on vient de le voir avec le polynôme de degré 8 (§4.3).

La validation croisée à K plis (K-fold CV) fournit une estimation honnête, applicable uniformément aux cinq modèles :

1. Diviser les 153 sites en K=5 groupes (plis)

2. Pour chaque pli : entraîner sur les 4 autres, évaluer sur celui-ci

3. Moyenner les R² des 5 plis

Pour que la CV fonctionne avec les GAM, on utilise un pipeline sklearn équivalent : SplineTransformer (base B-spline, knots sur quantiles) + Ridge (pénalisation L₂ ≈ pénalité de lissage). Le GLM n’a pas de wrapper scikit-learn standard : sa CV est effectuée manuellement, pli par pli, avec statsmodels.

Exercice 6 : validation croisée de l’OLS

Complète la cellule à trous ci-dessous pour découper les 153 sites en K = 5 plis, puis calculer le R² de l’OLS en validation croisée, à partir des variables de terrain X_th et du zinc y. Les cellules suivantes appliquent la même démarche aux autres modèles.

Pour cela, tu peux utiliser :

• KFold(n_splits=K, shuffle=True, random_state=...) définit un découpage des données en K plis. shuffle=True mélange les sites avant de les répartir, et random_state fixe ce mélange, pour obtenir le même découpage à chaque exécution.

• cross_val_score(modele, X, y, cv=decoupage, scoring='r2') entraîne modele sur tous les plis sauf un, puis l’évalue sur le pli restant, et recommence pour chaque pli. Elle renvoie un tableau avec un R² par pli.

• LinearRegression() est le modèle de régression linéaire (OLS) de scikit-learn, déjà utilisé au TP2. Inutile de l’entraîner avec .fit() : cross_val_score s’en charge pour chaque pli.

4.8 Incertitude liée au découpage — validation croisée répétée¶

Le tableau précédent donne une seule valeur de R² CV par modèle — celle obtenue avec random_state=42. Mais avec seulement 153 sites, répartis aléatoirement en 5 plis de ~31 observations, le découpage lui-même est une source de variabilité : un autre tirage donnerait un R² CV légèrement différent, simplement parce que les plis contiendraient des sites différents.

On répète la validation croisée 20 fois, avec un random_state différent à chaque fois, pour visualiser la distribution du R² CV plutôt qu’un seul chiffre. Plus le jeu de données est petit, plus cette variabilité est grande — un avertissement à garder en tête chaque fois qu’on compare deux R² CV calculés une seule fois.

⚠️ Ce n’est pas du bootstrap, même si les deux techniques « répètent une procédure aléatoire plusieurs fois » : ici, on ré-utilise toujours les 153 mêmes sites, seul le découpage en plis change. Le bootstrap (§4.9) tire un nouvel échantillon de sites à chaque répétition — une question différente.

4.9 Bootstrap, intervalle de confiance par ré-échantillonnage¶

La CV répétée (§4.8) fait varier le découpage en plis, mais utilise toujours les mêmes 153 sites. Le bootstrap répond à une question différente : si on avait observé un autre échantillon de sites dans la plaine de la Meuse (même processus, tirage différent), quelle serait la variabilité du R² ?

Le principe : tirer, $B$ fois, un nouvel échantillon de 153 sites avec remise depuis les données observées (certains sites apparaissent plusieurs fois, d’autres pas du tout), ré-ajuster le modèle sur chaque tirage, et regarder la distribution des R² obtenus. Pour que chaque réplication reste une évaluation hors échantillon (et non optimiste comme un R² d’entraînement), on évalue chaque modèle sur les sites jamais tirés dans cette réplication: les points out-of-bag (OOB), qui représentent en moyenne environ 37 % de l’échantillon original.

Bootstrap et CV répétée sont deux outils complémentaires, pas concurrents, pour quantifier l’incertitude d’un R² : la CV répétée mesure l’incertitude due au découpage train/test ; le bootstrap mesure l’incertitude due à l’échantillonnage de la population source. Dans un projet réel, les deux questions sont légitimes, et rien n’empêche de les utiliser ensemble.

Récapitulatif¶

Parties 1–3 : Géotraitement vectoriel (parallèle du TP QGIS 3)¶

Partie 4 : Régression non-linéaire : polynômes, GLM et GAM (jeu Meuse)¶

Concepts clés :

🎯 Bonnes pratiques ML : comparer et quantifier¶

Règle clé : le R² d’entraînement ne dit pas si le modèle généralisera, et un R² CV calculé une seule fois n’est lui-même qu’une estimation bruitée (§4.8) — le bootstrap (§4.9) en donne une seconde estimation, généralement plus large. Le GAM spatial capture un effet géographique réel (le gradient de la rivière) qui reste robuste en validation croisée.

📝 Quiz Moodle¶

Une fois ce TP terminé, teste tes connaissances avec le quiz Moodle du TP3 :

👉 Ouvrir le quiz du TP3 sur Moodle (accès réservé aux étudiant·e·s UNIL)