Interpolation spatiale

• L’interpolation spatiale est une méthode géostatistique permettant d’estimer les valeurs d’une variable à des emplacements non mesurés en utilisant des données ponctuelles existantes. Cette technique repose sur l’hypothèse fondamentale selon laquelle les valeurs proches dans l’espace sont plus similaires que les valeurs éloignées (principe de Tobler).

• L’interpolation spatiale a de nombreuses applications, notamment dans la modélisation de phénomènes environnementaux tels que les précipitations, les températures, ou encore la pollution de l’air. Elle permet de générer une surface continue à partir de points discrets. Le choix de la méthode d’interpolation dépend des caractéristiques de la donnée et des objectifs de l’analyse.

Différentes méthodes

• Méthodes disponibles :

  • Pondération par distance inverse (IDW)

  • Polygones de Thiessen

  • Kriging

  • Analyse de régression

    • Modèles linéaires

    • Modèles bayésiens

Chargement des données

L’objectif de cette étape est de préparer l’environnement pour l’analyse. Les bibliothèques utilisées incluent des outils pour :

• tmap : Création de cartes interactives ou statiques.

• terra : Gestion et visualisation des rasters et vecteurs. Ce package possède également des méthodes d’interpolation faciles à utiliser.

• gstat : pour les outils d’interpolations spatiales.

Vous êtes habitués !

# Charger les bibliothèques nécessaires
library(tmap)      # Pour la cartographie
library(tidyverse) # Pour le traitement des données 
library(terra)     # Pour les rasters et les vecteurs
library(gstat)     # Pour les méthodes d'interpolation
library(sf)        # Pour faire communiquer terra et gstat
library(yardstick) # Pour les métriques de cross-validation

Et pour les données

# Charger les données de précipitations
dta_precip <- vect(x = "Data/TP/precipitation_texas.gpkg")
projected_precip <- project(x = dta_precip, y = "EPSG:3081")
# Ajouter deux colonnes avec les coordonées des points 
projected_precip_coords <- cbind(projected_precip, crds(projected_precip))
# Charger la carte des limites du Texas
texas_admin <- vect(x = "Data/TP/texas_admin.gpkg")
projected_texas <- project(x = texas_admin, y = "EPSG:3081")

Késako ?

Les données importées comprennent :

• Un jeu de données de précipitations au Texas (valeurs ponctuelles).

• Les limites administratives du Texas (polygones).

Visualisation des données

Cette étape permet de visualiser les données de base pour vérifier leur contenu et leur distribution spatiale.

• Les points représentent les précipitations mesurées à des emplacements spécifiques.

• Le fond polygonal délimite le Texas, qui constitue notre zone d’étude.

• La carte nous aide à comprendre les relations spatiales et la couverture des données échantillonnées.

Une bonne visualisation initiale est essentielle pour identifier les éventuels biais spatiaux ou lacunes dans les données.

Méthode 1 : Polygones de Thiessen / Voronoi

Principe : Les polygones de Thiessen (ou diagrammes de Voronoï) attribuent à chaque point un espace où il est considéré comme le plus proche. Cette méthode crée des zones homogènes autour des points d’échantillonnage.

• Théorie sous-jacente :
  Chaque polygone est une zone d’influence unique où tout point intérieur est plus proche de son point d’échantillonnage que de tout autre. Cette méthode repose sur l’idée que chaque point représente une zone homogène, sans variation spatiale interne.

• Exemple concret :
  Considérez un réseau de stations météorologiques. Chaque station “représente” la zone qui l’entoure. Si vous avez une station mesurant les précipitations, les zones environnantes qui lui sont les plus proches auront des valeurs similaires.

Implémentation

• Cette méthode, aussi appelée interpolation par proximité, peut être réalisée avec la fonction voronoi du package terra.

# Create a tessellated surface
tesselated_surf <- voronoi(x = projected_precip_coords)
# Clip to the extent of Texas 
tesselated_texas <- intersect(x = tesselated_surf, y = projected_texas)

La cartographie les résultats !

Méthode 2 : Pondération par distance inverse (IDW)

Principe : L’IDW (Inverse Distance Weighting) applique un modèle où les valeurs estimées dépendent directement de la distance aux points mesurés : plus un point est proche, plus il contribue à l’estimation.

• Théorie sous-jacente :
  L’effet de distance est modulé par un paramètre de puissance (power). Une valeur power = 2 indique que l’influence d’un point sur son voisinage décroît avec le carré de la distance. Plus un point est éloigné, moins il influence l’estimation.

• Exemple concret :
  Dans l’estimation des précipitations, les stations proches auront une plus grande influence sur l’estimation des valeurs à des points non échantillonnés que celles situées à des centaines de kilomètres. Ce modèle suppose que la variabilité dans les données est fortement influencée par la proximité géographique.

Création d’une grille

# Create an empty grid where resolution is the size of the cells
grd <- rast(
    crs = "EPSG:3081",
    resolution = c(1e04, 1e04),
    vals = NA,
    extent = projected_texas
)

Faire communiquer gstat et terra

• Pour faire communiquer ces deux packages, nous devons d’abord créer un objet gstat qui contient les différentes informations pour notre modèle.

gstat_idw <- gstat(
    id = "Precip_in",
    formula = Precip_in ~ 1,
    data = st_as_sf(x = projected_precip_coords),
    set = list(idp = 2)
)

• Nous devons également créer une fonction qui va nous permettre de faire communiquer les différents packages entre eux !

interpolate_gstat <- function(model, x, crs, ...) {
    v <- st_as_sf(x = x, coords = c("x", "y"), crs = crs)
    p <- predict(model, v, ...)
    return(as.data.frame(p)[, 1:2])
}

Interpolation avec l’IDW

• Nous pouvons ensuite utiliser la commande interpolate pour notre interpolation

P.idw <- interpolate(
    object = grd,
    model = gstat_idw, 
    fun = interpolate_gstat, 
    crs = crs(grd)
)

[inverse distance weighted interpolation]
[inverse distance weighted interpolation]

# Mask to Texas
P.idw <- crop(x = P.idw$Precip_in.pred, y = projected_texas, mask = T, touches = T)

Cartographie

Le lien avec la régression linéaire

Principe : Cette méthode utilise un modèle statistique pour estimer les valeurs inconnues en fonction des coordonnées spatiales et, éventuellement, d’autres variables explicatives.

• Théorie sous-jacente :
  La régression linéaire suppose que la variable étudiée suit une relation fonctionnelle simple avec les prédicteurs. Elle utilise une équation de la forme : \(P = \alpha + \beta_{1} \cdot X + \beta_{2} \cdot Y\) où X et Y sont les coordonnées spatiales.

• Exemple concret :
  Modéliser les précipitations en fonction de l’altitude ou de la latitude, qui sont connues pour influencer ce phénomène.

Importation des données nécessaires

elevation <- rast(x = "Data/TP/elevation_texas.tif")

projected_elev <- project(x = elevation, y = grd, method = "bilinear") %>%
    crop(y = projected_texas, mask = T, touches = T)
names(projected_elev) <- "Elev"

Comment faire ?

# Add X and Y to P
dta_reg_lin <- data.frame(
    Precip_in = projected_precip_coords$Precip_in,
    Elev = projected_precip_coords$elevation_texas, 
    x = projected_precip_coords$x, 
    y = projected_precip_coords$y
)

# Run the regression model
reg_lin_interpolation <- lm(formula = Precip_in ~ x + y + Elev, 
           data = dta_reg_lin)

# Use the regression model output to interpolate the surface
interpolation_lm <- interpolate(object = c(grd, projected_elev),
                       model = reg_lin_interpolation,
                       xyNames = c("x", "y"), 
                       fun = predict)

Et les cartes !

Juste avec l’élévation

# Run the regression model
reg_lin_interpolation <- lm(formula = Precip_in ~ Elev, 
           data = dta_reg_lin)

# Use the regression model output to interpolate the surface
interpolation_lm <- interpolate(object = c(grd, projected_elev),
                       model = reg_lin_interpolation,
                       # xyNames = c("x", "y"), 
                       fun = predict)

En carte !

Et le code

# Leave-One-Out Cross-Validation Loop
IDW.out <- map(
    .x = 1:nrow(projected_precip_coords),
    .f = function(i) {
        # Creation of a gstat object with i - 1 observation
        tmp_gstat_idw <- gstat(
            id = "Precip_in",
            formula = Precip_in ~ 1,
            data = st_as_sf(projected_precip_coords)[-i, ],
            set = list(idp = 2)
        )
        # IDW interpolation excluding the i-th observation
        tmp_idw_result <- interpolate(
            object = grd,
            model = tmp_gstat_idw,
            fun = interpolate_gstat,
            crs = crs(grd),
            debug.level = 0
        )
        
        # Extract the interpolated value at the i-th location
        tmp_extracted <- extract(x = tmp_idw_result, 
                                 y = projected_precip_coords[i, ], 
                                 raw = FALSE) %>%
            mutate(observed = projected_precip_coords$Precip_in[i])  # Store observed value
        return(tmp_extracted)
    },
    .progress = T
)

Calcul des métriques

# Convert list to a data frame
validation.df <- bind_rows(IDW.out)

# Compute RMSE
rmse <- rmse(data = validation.df, truth = observed, estimate = Precip_in.pred)
# Compute R²
r_squared <- rsq_trad(data = validation.df, truth = observed, estimate = Precip_in.pred)

cat("LOOCV RMSE:", round(rmse$.estimate, 2), "\n")

LOOCV RMSE: 6.98 

cat("LOOCV R²:", round(r_squared$.estimate, 2), "\n")

LOOCV R²: 0.59 

Et le graphe !

Pour aller plus loin…

• Le kriging ! Pas matière d’examen : je ne vous demanderai pas de faire cette méthode d’interpolation !

• Tout d’abord, nous devons trouver la bonne formulation du variogramme théorique

• Le variogramme empirique s’estime avec la fonction variogram.

var_empi <- variogram(
    object = Precip_in ~ 1,
    data = st_as_sf(projected_precip_coords)
)

Le variogramme empirique

plot(var_empi)

Le variogramme théorique

• Grâce au variogramme théorique, nous pourrons estimer les valeurs pour toute la zone, peu importe sa distance aux points connus !
  • Plusieurs possibilités dans la forme du variogramme (vgm) : "Exp", "Sph", "Gau" et "Mat"
  • Pour choisir : optimisation :)

var_model <- fit.variogram(
    object = var_empi,
    model = vgm(model = c("Gau", "Exp", "Sph", "Mat")),
    fit.sills = T,
    fit.ranges = T,
    fit.method = T,
    fit.kappa = T, 
    warn.if.neg = T
)

Empirique et théorique

plot(var_empi, var_model)

Pour interpoler…

krige_model <- gstat(
    id = "Precip_in",
    formula = Precip_in ~ 1,
    data = st_as_sf(x = projected_precip_coords),
    model = var_model
)

Ce qui donne…

• Nous pouvons ensuite utiliser la fonction interpolate pour interpoler nos valeurs avec le modèle de kriging !

krige_surface <- interpolate(
    object = grd,
    model = krige_model,
    fun = interpolate_gstat,
    crs = crs(grd)
) %>%
    crop(y = projected_texas, touches = T, mask = T)

[using ordinary kriging]
[using ordinary kriging]

• Nous avons également confirmation que nous utilisons bien du kriging ordinaire.
• Vous pourriez également cross-valider votre résultat sur la même méthode que décrite plus haut. Faites bien attention à changer les données pour le variogramme empirique ! Ma RMSE = 5.53

La carte

Et sa déviation standard !

Question 1

L’ensemble de données sur lequel nous travaillons est un ensemble de données sur la qualité de l’air obtenu auprès de l’Agence Européenne pour l’Environnement (AEE) en Allemagne.

Votre travail est de comparer les méthodes d’interpolation vues lors de ce TP grâce aux données de concentration \(NO_2\). Pour l’interpolation par IDW, vous comparerez également deux valeurs différentes du paramètre idp. Pour la régression linéaire, j’ai ajouté la pente et l’élévation comme prédicteurs.

En guise d’exercice de dépassement, vous réaliserez une validation croisée sur le modèle de régression linéaire et l’interpréterez.

Si vous souhaitez, vous pouvez utiliser du kriging universel, c’est à dire utiliser également les valeurs de pente et d’élévation dans votre modèle ! Le kriging ne sera pas matière d’examen ;).