分区统计R（光栅/多边形）
在R中，与包“光栅”中的“提取”相比，在包“空间生态”的函数“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偏差的原因是什么？
共有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\u extract）。
每个算法都使用不同数量的像素来计算分区统计信息；因此，差异可能是由（z1高于z2）引起的。根据这个小例子，我可以推断出地带性。统计数据的限制性小于提取数据。因此，可能是带状的。统计信息将考虑多边形内光栅的每个值；但是，“提取”（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)
空间：地带性。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
-
我想将光栅数据聚合到自定义形状文件中的每个多边形。 在这种情况下，我想获得撒哈拉以南非洲次国家区域城市化的平均程度。 我的sf如下所示： 或绘制： 另一方面，光栅数据采用以下形式： 这些比整个星球所需的要细得多。为了加速计算，我首先聚合光栅，然后将其转换为shapefile，剩余的每个光栅像素都转换为shapefile中的点几何形状。然后，这个shapefile可以聚合到我的区域边界。诚然，这不是
-
我想计算区域内的分区统计数据，以符合等宽。 得到区段列表后，想法是分配块号，因为栅格::区段函数需要一个带有表示区段代码的栅格层。 当我尝试用分区编号填充范围时，填充的区域与范围不对应（请参见图）。为什么会这样？
-
下面的代码在我的图像上生成两个框。我正计划进一步分析这些框内的像素。 在下面的例子中，在红色方块的情况下，我不想继续下去，因为它的右上角有黑色像素。而我想继续在绿色方块的情况下，因为它没有一个黑色像素沿着它的边缘。
-
我正在尝试在R中设置一个randomForest，以便根据其他光栅图像对光栅图像进行分类。我的训练数据是一个完全填充的光栅图像，我想训练许多其他光栅，以尝试基于初始光栅创建光栅输出。代码示例如下： <代码>rf1 ...其中，是我的光栅格式的实际已知值，而到是我想用来预测trainingRaster1是什么的其他光栅图像。我知道您将使用向量或点的训练类来训练一系列光栅，但在我的情况下，我希望使用光
-
我已经为此挣扎了几个小时。我有一个包含177个多边形（即177个县）的shapefile（称为“shp”）。这个shapefile覆盖在光栅上。我的光栅（称为“ras”）由具有不同污染值的像素组成。 现在我想提取每个多边形的所有像素值及其出现次数。 这正是QGIS功能“分区直方图”所做的。但我想在R中做同样的事情。 我尝试了提取（）函数，并设法获得了每个县的平均值，这已经是第一步，但我想制作像素分
-
我想知道是否有办法为光栅图层对象创建分区统计数据，特别是R中给定单元值（例如土地使用类别）的计数，而无需重新分类整个光栅。解决方案应具有内存效率，以便处理大型光栅文件，即不需要将值提取到R中的矩阵中。 下面是一个到目前为止我是如何处理它的示例。在这种情况下，我将原始光栅重新分类为仅保留1的利息值和所有其他值的缺失值。 我提出的解决方案创建了冗余数据和额外的处理步骤，以实现我的初始目标。我认为类似于