Packages nécessaires pour l’analyse

# Libraries
library(terra)
library(dismo)
library(tmap)
library(mapview)
library(mapedit)

Comparaison entre les jeux de données de biomes

Pour l’exemple, nous utilisons des couches raster des jeux de données de biomes inclus dans le package {dismo}. Bien que ces jeux de données partagent la même étendue, origine et résolution, nous allons transformer grossièrement l’un d’eux pour simuler une étendue, une origine et une résolution différentes. C’est une situation courante lorsque l’on travaille avec des données spatiales, et probablement celle à laquelle vous serez confronté. Nous allons ensuite découper, agréger et rééchantillonner les deux rasters pour avoir la même étendue, origine et résolution. Cela nous permettra de les comparer.

Nous utilisons bio1 = Température annuelle moyenne (°C * 10) et bio12 = Précipitations annuelles totales (mm). Nous allons transformer bio12 pour simuler une étendue, une origine et une résolution différentes.

# get list of predictor variables in dismo
fnames <- list.files(
  path = paste(system.file(package = "dismo"), "/ex", sep = ""),
  pattern = "grd", full.names = TRUE
)

# Get Bio1 and Bio12
bio1 <- rast(fnames[1])
bio12 <- rast(fnames[2])

# Transform Bio12 to simulate a different extent, origin and resolution
bio12_trans <- aggregate(bio12, 3)

Visualisation des motifs spatiaux

Des visualisations claires et engageantes sont cruciales dans l’analyse spatiale. Le package {tmap} facilite cela tout en offrant une qualité prête pour la publication. Cartographions nos variables côte à côte.

# Temperature and precipiation maps
tm1 <- tm_shape(bio1) +
  tm_raster(
    title = "Temperature (°C)", palette = "-RdBu",
    midpoint = mean(values(bio1), na.rm = TRUE)
  ) +
  tm_layout(title = "Mean annual temperature")

tm2 <- tm_shape(bio12_trans) +
  tm_raster(title = "Precipitation", palette = "Greens") +
  tm_layout(title = "Total annual precipitation (mm)")

# Arrange maps
tmap_arrange(tm1, tm2)

Ces cartes mettent en évidence la variation spatiale des deux variables et suggèrent des corrélations potentielles.

Rééchantillonnage et agrégation des rasters pour permettre la comparaison

Nous devons nous assurer que les deux rasters ont des résolutions compatibles. Nous utilisons la plus petite résolution des deux rasters pour rééchantillonner l’autre. Nous allons également empiler les deux rasters pour une analyse plus approfondie.

# Aggregate and resample chlorophyll raster
bio1_r <- resample(bio1, bio12_trans)

# Stack rasters for analysis
bio_1_12 <- c(bio1_r, bio12_trans)
names(bio_1_12) <- c("temperature", "precipitation")

# Temperature and precipiation maps
tm1 <- tm_shape(bio1_r) +
  tm_raster(
    title = "Temperature (°C)", palette = "-RdBu",
    midpoint = mean(values(bio1_r), na.rm = TRUE)
  ) +
  tm_layout(title = "Mean annual temperature")

tm2 <- tm_shape(bio12_trans) +
  tm_raster(title = "Precipitation", palette = "Greens") +
  tm_layout(title = "Total annual precipitation (mm)")

# Arrange maps
tmap_arrange(tm1, tm2)

Exploration des relations entre les variables

L’un de nos objectifs principaux est d’examiner la relation entre la température et les précipitations. Nous pouvons commencer par calculer la corrélation globale entre les deux rasters.

La corrélation globale est une opération habituelle utilisant la fonction cor() entre les valeurs des deux rasters. Cela renvoie une seule valeur.

# Correlation between layers
cor(
  values(bio_1_12)[, 1],
  values(bio_1_12)[, 2],
  use = "na.or.complete"
)

## [1] 0.6220131

Estimer un modèle linéaire

Une autre façon de voir s’il y a une corrélation entre les deux couches est de définir un modèle linéaire (ici lm(precipitation ~ temperature)). Pour examiner l’hétérogénéité spatiale, nous pouvons tracer la distribution des résidus du modèle. (N’oubliez pas de récupérer les résidus manquants des valeurs NA.). Cette carte montre où les observations s’écartent du modèle.

lm1 <- lm(values(bio_1_12)[, 2] ~ values(bio_1_12)[, 1])
summary(lm1)

## 
## Call:
## lm(formula = values(bio_1_12)[, 2] ~ values(bio_1_12)[, 1])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1509.9  -421.8   -31.5   368.6  3904.2 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           -271.2116    68.5785  -3.955 8.23e-05 ***
## values(bio_1_12)[, 1]    8.0648     0.3304  24.407  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 638.2 on 944 degrees of freedom
##   (3022 observations effacées parce que manquantes)
## Multiple R-squared:  0.3869,	Adjusted R-squared:  0.3863 
## F-statistic: 595.7 on 1 and 944 DF,  p-value: < 2.2e-16

# Create and plot residual raster
resid_r <- rast(bio_1_12, 1)
values(resid_r) <- NA
values(resid_r)[
  !is.na(values(bio_1_12)[, 1]) &
    !is.na(values(bio_1_12)[, 2])
] <- lm1$residuals

tm_shape(resid_r) +
  tm_raster(title = "Residuals", palette = "-RdBu", midpoint = 0) +
  tm_layout(title = "Model Residuals")

Calculer la corrélation locale en utilisant focal() sur deux rasters en même temps

La corrélation locale est une manière différente de voir la relation entre deux covariables dans l’espace. L’idée est simple. Autour de chaque cellule des rasters, nous définissons une zone focale (carrée, ronde, …). Pour l’exemple, j’ai choisi un carré de 5*5. Nous calculons la corrélation entre les 25 valeurs de chaque raster dans ce carré en utilisant la fonction cor(). Cette corrélation locale est enregistrée pour la cellule centrale. Cette opération est répétée sur l’ensemble du raster grâce à la fonction focal() de la bibliothèque {raster}.
Cependant, nous devons utiliser une astuce pour utiliser la fonction sur deux rasters en même temps. Nous définissons un troisième raster (bio_1_12_nb) qui enregistre les positions des valeurs de 1 à ncell(). focal() extraira ainsi les positions pour lesquelles nous pouvons obtenir les valeurs des deux rasters (comme ils ont maintenant la même étendue, origine et résolution).

Attention si vous manipulez de grands rasters, cela peut utiliser plus de RAM.

bio_1_12_nb <- rast(bio_1_12, 1)
values(bio_1_12_nb) <- 1:ncell(bio_1_12)

matrix_bio_1_12 <- values(bio_1_12) # stack as raster [MW]

focal_cor <- focal(
  x = bio_1_12_nb,
  w = matrix(1, 5, 5),
  fun = function(x, y = matrix_bio_1_12) {
    cor(y[x, 1], y[x, 2],
      use = "na.or.complete"
    )
  },
  filename = file.path(extraWD, "focal_cor.tif"),
  overwrite = TRUE
)

La carte résultante montre les zones où la corrélation entre les précipitations et la température est négative, nulle ou positive. C’est une manière différente de voir la relation entre les covariables.

focal_cor <- raster(file.path(extraWD, "focal_cor.tif"))

# Plot local correlation
tm_shape(focal_cor) +
  tm_raster(title = "Local Correlation", palette = "-RdYlGn", midpoint = 0) +
  tm_layout(title = "Local Correlation Map")

Ce billet de blog est une réécriture d’un billet de blog précédent que j’ai écrit en 2018 : Corrélation spatiale entre rasters. J’ai utilisé ChatGPT avec cette commande : Je voudrais utiliser le package 'terra' au lieu de 'raster'. Et utiliser 'tmap' au lieu de 'ggplot2' et le code résultant n’a pas eu besoin de modifications supplémentaires. Ainsi, si vous aussi, avez besoin de mettre à jour votre code utilisant le package {raster}, vous pouvez utiliser ChatGPT pour vous aider…

Le script R de cet article est disponible sur mon dépôt GitHub Blog tips