区域统计 R（栅格/多边形）

Question

在 R 中，计算包 'spatialEco' 的函数 'zonal.stats' 的平均值与包 'raster' 中的 'extract' 相比存在偏差。对于这两者，我都使用多边形作为区域字段，并使用栅格作为值。

这是一个例子：

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)

z2和z1的偏差是什么原因造成的？

Answer 1

spatialEco::zonal.stats 使用exactextractr（我没有检查代码，但它告诉我安装它以便能够使用zonal.stats）如果您正在考虑多边形应该更准确（光栅包将它们变成首先进入光栅，见下面的zonal）。然而，下面的例子（这只是一种情况）表明 spatialEco 不太精确。

示例（避免使用随机数，但如果您确实使用它们，请使用set.seed）。我从非常大的网格单元开始。

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")

使用栅格，您还可以像这样获得更精确的估计

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

查看使用了多少个单元格

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

查看此问题的一种方法是创建更高分辨率的栅格数据。

分辨率提高 10 倍，数值相同

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

分辨率提高 100 倍，值相同

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

在最高分辨率下，平均值相似（细胞数量的相对差异很小）；但是栅格在较低的分辨率（以及加权平均值）下更接近正确的值（无论是什么，确切地说）。这是出乎意料的。

为了更快的速度，现在还有terra 包

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

Answer 2

每种算法都使用不同数量的像素来计算区域统计信息；因此，差异可能是由（z1高于z2）引起的。根据这个小例子，我会推断zonal.stats 的限制比extract 的限制要少。因此，zonal.stats 可能会考虑落在多边形内的栅格的每个值；但是，extract 仅考虑其中心位于多边形内的像素（查看函数文档）。

# 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)

Answer 3

感谢@Robert Hijmans 的见解。根据您的示例，我通过计算不同功能和分辨率的区域均值对 1）精度和 2）计算时间进行了进一步比较：

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

基于此示例，当使用多边形作为区域时，exact_extract 是首选。