#1. Introduction et préparation des données

Ce document présente une analyse géostatistique des sols du Jura suisse. L’objectif est d’étudier la variabilité spatiale de la concentration en zinc (Zn) sur 259 points de mesures pour réaliser une interpolation par krigeage.

1.1. Acquisition et description des données

Le jeu de données utilisé pour cette étude provient directement du package R gstat. Il est issu d’une campagne d’échantillonnage réalisée dans le massif du Jura suisse, couvrant une surface d’environ 14,5 km². Les mesures portent sur la concentration en divers métaux lourds dans la couche superficielle du sol (0 à 20 cm de profondeur).

Pour notre analyse, nous nous focalisons exclusivement sur le zinc (Zn), mesuré sur un semis de 259 stations d’observation (prediction.dat). L’intérêt majeur de ce jeu de données est qu’il fournit non seulement les points d’échantillonnage avec leurs coordonnées spatiales (locales, exprimées en kilomètres), mais également une grille de prédiction régulière (jura.grid) qui servira de support géométrique pour notre interpolation spatiale finale.

# Chargement des packages
library(gstat)
library(sp)
library(sf)
library(mapsf)
library(bestNormalize)
library(tmap)
library(tmap.mapgl)
library(terra)

# Chargement du jeu de données jura
data(jura)
donnees_zn <- prediction.dat[, c("Xloc", "Yloc", "Zn")]

# Export du jeu de données en CSV
write.csv(donnees_zn, "donnees_jura_zinc.csv", row.names = FALSE)

# Conversion en objet spatial sf pour la cartographie
pts_zn_sf <- st_as_sf(donnees_zn, coords = c("Xloc", "Yloc"), remove = FALSE)
# On lui assigne un CRS générique local car les coordonnées sont relatives
st_crs(pts_zn_sf) <- NA

#2. Cartographie exploratoire

La carte ci-dessous illustre la répartition spatiale des teneurs en zinc.

# Création d'un polygone de fond pour la zone d'étude
grille_pts <- st_as_sf(jura.grid[, c("Xloc", "Yloc")], coords = c("Xloc", "Yloc"))
st_crs(grille_pts) <- NA
zone_etude_sf <- st_sf(geometry = st_buffer(st_convex_hull(st_union(grille_pts)), dist = 0.4))

# Initialisation de la carte
mf_init(zone_etude_sf)
mf_map(zone_etude_sf, col = "ivory", border = "gray40", lwd = 1.2, add = TRUE)

# Cercles proportionnels
mf_map(
  x = pts_zn_sf,
  var = "Zn",
  type = "prop",
  inches = 0.15,
  col = "cadetblue",
  border = "white",
  lwd = 1,
  leg_title = "Concentration en Zinc\n(mg/kg)",
  leg_pos = "topright",
  add = TRUE
)

mf_title("Distribution spatiale des teneurs en zinc (Jura suisse)")
mf_scale(size = 2, pos = "bottomleft")
mf_arrow(pos = "topleft")

La carte exploratoire met en évidence une forte hétérogénéité spatiale des concentrations en zinc. On observe un gradient géographique très net : les teneurs les plus élevées (cercles les plus larges) se concentrent massivement dans la partie sud-ouest de la zone d’étude, tandis que le secteur nord-est présente des valeurs significativement plus faibles. Cette forte structuration spatiale justifie pleinement le recours à une modélisation par krigeage.

3. Normalisation et Variographie

Avant de réaliser le variogramme expérimental et le krigeage, les données asymétriques doivent être normalisées pour stabiliser la variance.

# Normalisation avec bestNormalize
bn_result <- bestNormalize(donnees_zn$Zn)
donnees_zn$Zn_trans <- bn_result$x.t

# Conversion en SpatialPointsDataFrame (sp) pour gstat
pts_zn_sp <- donnees_zn
coordinates(pts_zn_sp) <- ~Xloc+Yloc

# Variogramme expérimental avec tendance linéaire (Xloc, Yloc)
vario_exp <- variogram(Zn_trans ~ Xloc + Yloc, pts_zn_sp)

# Ajustement du modèle sphérique (nugget estimé à 0.2, palier à 1, portée à 1.5)
modele_ajuste <- fit.variogram(vario_exp, vgm(psill = 1, model = "Sph", range = 1.5, nugget = 0.2))

# Affichage
plot(vario_exp, modele_ajuste, main = "Modèle sphérique ajusté (avec tendance linéaire)")

Les données environnementales présentant une forte asymétrie, la transformation via bestNormalize a permis de stabiliser la variance et de s’approcher d’une distribution normale, une condition optimale pour le krigeage. Le semi-variogramme expérimental, calculé après l’extraction de la tendance linéaire (pour neutraliser le gradient Sud-Ouest / Nord-Est), révèle une autocorrélation spatiale claire. L’ajustement du modèle sphérique permet de quantifier un effet pépite (traduisant la micro-variabilité locale ou l’erreur de mesure) et une portée, définissant la distance au-delà de laquelle les teneurs en zinc ne sont plus spatialement corrélées.

4. Interpolation par Krigeage et Validation

Nous utilisons la grille de prédiction fournie pour estimer les valeurs de zinc via un Krigeage Universel. Une validation croisée est effectuée pour évaluer la qualité du modèle.

# Préparation de la grille
grille_sp <- jura.grid[, c("Xloc", "Yloc")]
coordinates(grille_sp) <- ~Xloc+Yloc
gridded(grille_sp) <- TRUE

# Krigeage
krigeage_res <- krige(Zn_trans ~ Xloc + Yloc, locations = pts_zn_sp, newdata = grille_sp, model = modele_ajuste)
## [using universal kriging]
# Transformation inverse (retour aux mg/kg)
krigeage_res$Pred_Zn <- predict(bn_result, krigeage_res$var1.pred, inverse = TRUE)

# Validation croisée
validation <- krige.cv(Zn_trans ~ Xloc + Yloc, locations = pts_zn_sp, model = modele_ajuste)
validation$residus_vrais <- donnees_zn$Zn - predict(bn_result, validation$var1.pred, inverse = TRUE)

# Affichage des résidus
hist(validation$residus_vrais, breaks = 30, col = "lightblue", border = "white", main = "Histogramme des résidus de la validation croisée", xlab = "Erreur (mg/kg)")

Le krigeage universel exploite conjointement le modèle de dépendance spatiale (le variogramme) et la dérive régionale (la tendance linéaire) pour estimer les concentrations sur l’ensemble de la grille de prédiction. L’histogramme issu de la validation croisée montre des résidus (différences entre les valeurs observées et les valeurs prédites) globalement centrés sur zéro et d’allure symétrique. Ce résultat démontre que le modèle d’interpolation est robuste, non biaisé, et qu’il ne sous-estime ni ne surestime systématiquement les teneurs en zinc.

5. Cartographie finale interactive

Les coordonnées du jeu de données étant locales, nous les convertissons en mètres et les décalons vers les véritables coordonnées du massif du Jura (en projection Lambert 93) pour permettre l’affichage des fonds de carte interactifs.

# 1. Extraction des prédictions et décalage géographique
df_map <- as.data.frame(krigeage_res)

# Conversion des km en mètres (* 1000) et positionnement dans le Jura (EPSG:2154)
df_map$X_reel <- (df_map$Xloc * 1000) + 950000
df_map$Y_reel <- (df_map$Yloc * 1000) + 6650000

# 2. Création du raster géoréférencé avec terra
krigeage_raster <- rast(df_map[, c("X_reel", "Y_reel", "Pred_Zn")], type = "xyz", crs = "EPSG:2154")

# 3. Création du masque de détourage (enveloppe convexe autour des points réels)
pts_sf <- st_as_sf(df_map, coords = c("X_reel", "Y_reel"), crs = 2154)
zone_masque <- st_buffer(st_convex_hull(st_union(pts_sf)), dist = 150)

# Application du détourage sur le raster
krigeage_raster_masque <- mask(krigeage_raster, vect(zone_masque))

# 4. Affichage interactif des deux cartes
tmap_mode("view")

# Carte 1 : Fond de carte classique
carte_classique <- tm_basemap("OpenStreetMap") +
  tm_shape(krigeage_raster_masque) +
  tm_raster(
    col = "Pred_Zn",
    style = "jenks",
    n = 6,
    palette = "-YlOrRd",
    title = "Prédiction Zn (mg/kg)",
    alpha = 0.7
  ) +
  tm_shape(zone_masque) +
  tm_borders(col = "black", lwd = 2)

# Carte 2 : Fond de carte satellite
carte_satellite <- tm_basemap("Esri.WorldImagery") +
  tm_shape(krigeage_raster_masque) +
  tm_raster(
    col = "Pred_Zn",
    style = "jenks",
    n = 6,
    palette = "-YlOrRd",
    title = "Prédiction Zn (mg/kg)",
    alpha = 0.7
  ) +
  tm_shape(zone_masque) +
  tm_borders(col = "white", lwd = 2.5)

# Impression des deux cartes dans le document
carte_classique
carte_satellite

Ces deux cartes interactives, repositionnées géographiquement et détourées sur l’emprise réelle, offrent une modélisation continue et lissée de la zone étudiée. La classification met parfaitement en valeur le pic de concentration au sud-ouest. Le basculement vers l’imagerie satellite apporte un contexte d’occupation des sols essentiel : il permet de confronter la modélisation mathématique à la réalité topographique et environnementale du terrain pour mieux interpréter l’origine de ces métaux lourds.

6. Discussion

L’approche géostatistique adoptée dans cette étude a permis d’obtenir des résultats cohérents, mais soulève quelques points de discussion.

Premièrement, l’utilisation d’un krigeage universel (intégrant une dérive linéaire sur les coordonnées X et Y) s’est avérée indispensable pour corriger le fort gradient spatial du zinc et modéliser correctement la structure spatiale résiduelle. Toutefois, le jeu de données jura contient d’autres variables environnementales très corrélées au zinc (comme le cadmium, le plomb ou le cuivre) ainsi que des informations sur l’usage des sols. L’intégration de ces informations secondaires via un co-krigeage (ou un krigeage avec dérive externe) aurait potentiellement permis de réduire la variance d’estimation et d’affiner la précision des prédictions.

Deuxièmement, d’un point de vue cartographique, les coordonnées d’origine du jeu de données sont exprimées dans un référentiel local arbitraire. Leur repositionnement pour la cartographie interactive satellite (Partie 5) a nécessité une translation mathématique (offset) vers le système EPSG:2154. Bien que visuellement très pertinente pour comprendre le contexte paysager, cette superposition conserve une marge d’incertitude quant à son alignement géographique absolu.

7. Conclusion

Cette analyse a permis de modéliser avec succès la distribution spatiale des teneurs en zinc dans les sols de cette zone du Jura suisse. La chaîne de traitement mise en place (normalisation de la variable, modélisation variographique avec prise en compte d’une tendance, puis krigeage universel) s’est révélée robuste, ce qui a été confirmé par l’analyse des résidus de la validation croisée.

La cartographie finale met en évidence une concentration accrue de zinc dans le secteur sud-ouest de la zone d’étude, s’estompant progressivement vers le nord-est. Le passage d’un semis de points irrégulier à une surface continue lissée offre une vision globale du phénomène. Ce type de modèle prédictif interactif constitue un excellent outil d’aide à la décision pour le diagnostic environnemental, permettant par exemple d’identifier les zones dépassant certains seuils de pollution et de cibler efficacement de futures campagnes de prélèvements.