1
votes

Statistiques zonales R (raster / polygone)

Dans R, il y a des écarts calculant la moyenne de la fonction 'zonal.stats' du package 'spatialEco' par rapport à 'extract' dans le package 'raster'. Pour les deux, j'ai utilisé un polygone comme champ de zone et un raster pour les valeurs.

Voici un exemple:

library(raster)
library(spatialEco)
library(sp)

#Create raster
ras <- raster(nrows=100, ncols=80, xmn=0, xmx=1000, ymn=0, ymx=800)
val <- runif(ncell(ras))
values(ras) <- val

#Create polygon within raster extent
xym <- cbind(runif(3,0,1000), runif(3,0,800))
p <- Polygons(list(Polygon(xym)),1)
sp <- SpatialPolygons(list(p))
spdf <- SpatialPolygonsDataFrame(sp, data=data.frame(1))

#z1 zonal statistics using "spatialECO"
z1 <- zonal.stats(spdf, ras, stats="mean")

#z2 zonal statistics using "raster"
z2 <- extract(ras, spdf, fun=mean)

Qu'est-ce qui cause l'écart de z2 et z1?

r statistics polygon raster spatial

0 commentaires

3 Réponses :

3
votes

Chaque algorithme utilise un nombre différent de pixels pour calculer les statistiques zonales; ainsi, la différence est probablement causée par cela (z1 est supérieur à z2). D'après ce petit exemple, j'inférerais que zonal.stats est moins restrictif que extract . Ainsi, zonal.stats prend probablement en compte chaque valeur du raster qui se trouve à l'intérieur du polygone; cependant, extract ne prend en compte que les pixels dont le centre se trouve à l'intérieur du polygone (consultez la documentation de la fonction).

# Create a function to count the number of pixels used to calculate the zonal stats
counter <- function(x, na.rm = T) { 
  length(x)
} 

#z1 zonal statistics using "spatialECO"
z1 <- zonal.stats(spdf, ras, stats="counter")

#z2 zonal statistics using "raster"
z2 <- extract(ras, spdf, fun=counter)

1 commentaires

Merci de clarifier. En effet, la fonction «extraire» fournit une description dans la documentation du nombre de pixels choisis et offre des possibilités de modification. «zonal.stats» reste un peu une boîte noire. Je m'en tiendrai alors à «extraire» (même si c'est beaucoup plus lent ..)

1
votes

spatialEco :: zonal.stats utilise exactextractr (je n'ai pas vérifié le code, mais il m'a dit de l'installer pour pouvoir utiliser zonal. stats ) qui devrait être plus précis si vous envisagez des polygones (le package raster les transforme d'abord en rasters, voir zonal ci-dessous). Cependant, l'exemple ci-dessous (ce n'est qu'un cas) suggère que spatialEco est moins précis.

Exemple (évitez les nombres aléatoires, mais si vous les utilisez, utilisez set.seed ). Je commence par de très grandes cellules de grille.

library(terra)    
r <- rast(ras)
v <- vect(sp)
extract(r, v, "mean")    
#     ID    layer
#[1,]  1 5.271039

Avec raster, vous pouvez également obtenir une estimation plus précise comme celle-ci

ras <- disaggregate(ras, 10)
zonal.stats(sp, ras, stats="mean")
#mean.layer
#1   5.299915
extract(ras, sp, small=TRUE, fun=mean)
#         [,1]
#[1,] 5.271039


zonal.stats(sp, ras, stats="counter")
# counter.layer
#1         17695

extract(ras, sp, fun=function(x,...)length(x))
#      [,1]
#[1,] 17289

Pour voir combien de cellules sont utilisées

ras <- disaggregate(ras, 10)
zonal.stats(sp, ras, stats="mean")
#  mean.layer
#1        5.5
extract(ras, sp, fun=mean)
#         [,1]
#[1,] 5.245614

zonal.stats(sp, ras, stats="counter")
#  counter.layer
#1           218
extract(ras, sp, fun=function(x,...)length(x))
#     [,1]
#[1,]  171

Une façon de voir cela est de créer des données raster de résolution plus élevée.

Résolution 10x plus élevée, mêmes valeurs

zonal.stats(sp, ras, stats="counter")
#  counter.layer
#1             6
extract(ras, sp, fun=function(x,...)length(x))
#     [,1]
#[1,]    3

Résolution 100x plus élevée, mêmes valeurs

e <- extract(ras, sp, weights=T)[[1]]
weighted.mean(e[,1], e[,2])
#[1] 5.269565

À la résolution la plus élevée, les valeurs moyennes sont similaires (et la différence relative le nombre de cellules est petit); mais le raster était plus proche de la valeur correcte (quoi que ce soit, exactement) à une résolution inférieure (et aussi avec la moyenne pondérée). C'est inattendu.

Pour une meilleure vitesse, il existe maintenant aussi le package terra

library(raster)
library(spatialEco)

ras <- raster(nrows=4, ncols=4, xmn=0, xmx=1000, ymn=0, ymx=800)
values(ras) <- 1:ncell(ras)
set.seed(1)
xy <- cbind(runif(3,0,1000), runif(3,0,800))
xy <- rbind(xy, xy[1,])
sp <- spPolygons(xy, attr=data.frame(x=1))


### zonal statistics using "spatialECO"
zonal.stats(sp, ras, stats="mean")
#  mean.layer
#1          7
### zonal statistics using "raster"
extract(ras, sp, fun=mean)
#     [,1]
#[1,]    6
### same as 
# x <- rasterize(sp, ras)
# zonal(ras, x, "mean")

0 commentaires

1
votes

Merci à @Robert Hijmans pour vos idées. Sur la base de votre exemple, j'ai fait quelques comparaisons supplémentaires de 1) précision et 2) temps de calcul en calculant les moyennes zonales pour différentes fonctions et résolutions:

library(raster)
library(spatialEco)
library(terra)
library(exactextractr)
library(sf)

ras <- raster(nrows=4, ncols=4, xmn=0, xmx=1000, ymn=0, ymx=800)
values(ras) <- 1:ncell(ras)
set.seed(1)
xy <- cbind(runif(3,0,1000), runif(3,0,800))
xy <- rbind(xy, xy[1,])
sp <- spPolygons(xy, attr=data.frame(x=1))

mn <- data.frame(matrix(ncol=6, nrow=4))
colnames(mn) <- c("disagr", "raster", "raster_weight", "spatialEco", "exactextractr", "terra")
mn[,1] <- c(2,10,50,250)
tm <- mn

for (i in 1:5){
  d <- mn[i,1]
  rasd <- disaggregate(ras, d)
  
  on <- Sys.time()
  mn[i,2] <- raster::extract(rasd, sp, fun=mean)
  off <- Sys.time()
  tm[i,2] <- off - on
  
  on <- Sys.time()
  mn[i,3] <- raster::extract(rasd, sp, fun=mean, weights=T)[[1]]
  off <- Sys.time()
  tm[i,3] <- off - on
  
  on <- Sys.time()
  mn[i,4] <- spatialEco::zonal.stats(sp, rasd, stats="mean")
  off <- Sys.time()
  tm[i,4] <- off - on
  
  on <- Sys.time()
  mn[i,5] <- exactextractr::exact_extract(rasd, st_as_sf(sp), fun="mean")
  off <- Sys.time()
  tm[i,5] <- off - on

  on <- Sys.time()
  val <- terra::extract(rast(rasd), vect(sp))
  mn[i,6] <- mean(val[,2])
  off <- Sys.time()
  tm[i,6] <- off - on
  print(i)
}
mn  # arithmetic mean
  disagr   raster raster_weight spatialEco exactextractr    terra
1      2 5.333333      5.269565   6.647059      5.271303 5.333333
2     10 5.245614      5.271039   5.500000      5.271303 5.245614
3     50 5.272370      5.271328   5.325525      5.271303 5.272370
4    250 5.271314      5.271303   5.282827      5.271303 5.271314
tm   # computing time in seconds
  disagr     raster raster_weight  spatialEco exactextractr       terra
1      2 0.03998685    0.03598809 0.003000021   0.002999067 0.008996964
2     10 0.09783196    0.04598618 0.003997803   0.003000021 0.008984089
3     50 0.37189507    0.40886688 0.004998922   0.003998041 0.021993160
4    250 4.10671687    8.29134583 0.035988092   0.019991875 0.336881876

Sur la base de cet exemple et lors de l'utilisation de polygones comme zones exact_extract est le choix préférable.

0 commentaires