区域统计 R (raster/polygon)
Zonal Statistics R (raster/polygon)
在 R 中，与包 'raster' 中的 'extract' 相比，包 'spatialEco' 中的函数 'zonal.stats' 计算平均值存在偏差。对于两者，我都使用多边形作为区域字段，并使用栅格作为值。
这是一个例子：
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出现偏差的原因是什么？
每种算法都使用不同数量的像素来计算区域统计数据；因此，差异可能是由此引起的（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)
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
感谢@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 是更好的选择。
在 R 中，与包 'raster' 中的 'extract' 相比，包 'spatialEco' 中的函数 'zonal.stats' 计算平均值存在偏差。对于两者，我都使用多边形作为区域字段，并使用栅格作为值。
这是一个例子：
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出现偏差的原因是什么？
每种算法都使用不同数量的像素来计算区域统计数据；因此，差异可能是由此引起的（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)
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
感谢@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 是更好的选择。