如何将栅格中的面积加权和提取到 R 中的多边形中?

How can I extract an area weighted sum from a raster into a polygon in R?

我有一个栅格(class "RasterLayer"),它包含人数 N。我有代表行政边界的多边形(class "SpatialPolygonsDataFrame")。

我想将栅格图层中的值提取到多边形中,方法是根据覆盖多边形的栅格图块的面积比例进行加权。假设有四个栅格图块,其值(4、8、12、16)与单个多边形重叠。如果每个栅格的 25% 与多边形重叠,我希望它们对多边形总数的贡献为 (1, 2, 3, 4),并且提取到该多边形中的总加权和为 10。

基本上,我想完成以下任务:

 some_polygons$n_people = extract(count_raster$n_people, 
                                   some_polygons, 
                                   fun=sum, 
                                   weights=T, 
                                   na.rm=T)

然而,当我在 R 中 运行 这个命令时,我收到以下内容:

"Warning message: In .local(x, y, ...) : "fun" was changed to "mean"; other functions cannot be used when "weights=TRUE""

R 在使用权重时强制变为 fun=mean。我需要使用权重来说明落在多边形中的栅格的比例,但是平均值不会让我计算我想要的数量。我需要找到一种方法来将栅格中的计数 N 相加,加权栅格中落在多边形内的面积(加权和不是加权平均值)。有谁知道我该怎么做?

我很高兴尝试使用其他空间数据 classes 来实现它。我尝试将栅格表示为 SpatialPixelsSpatialPointsSpatialGrid 并使用 sp::over() 函数,但无法弄清楚如何计算我的加权和需要。非常感谢您的帮助!

首先,请做一个可重现的例子,这将有助于我们帮助你;-)

* 编辑 这对版本 raster v2.5-8 有效,但对 raster v2.3-12!

无效

其中weights解释如下:

If TRUE and normalizeWeights=FALSE, the function returns, for each polygon, a matrix with the cell values and the approximate fraction of each cell that is covered by the polygon(rounded to 1/100). If TRUE and normalizeWeights=TRUE the weights are normalized such that they add up to one. The weights can be used for averaging; see examples. This option can be useful (but slow) if the polygons are small relative to the cells size of the Raster* object

所以你可以用它来计算你的加权和:

library(raster)
library(sp)

# Reproducible example
set.seed(13)
N = raster(nrows=4, ncol=4, xmn=0, xmx=4, ymn=0, ymx=4)
N[] = runif(16, 50, 100)

Ps1 = Polygons(list(Polygon(cbind(c(0.5,2.8,3.7,1.2,0.5), c(2,3.4,1.2,1.1,2)))), 
               ID = "a")
SPs = SpatialPolygons(list(Ps1))
poly = SpatialPolygonsDataFrame(SPs, data.frame(onecol = c("one"), 
                                                row.names = c("a")))
# See plot below
plot(N)
plot(poly, add=T, border='red')

# Extract with the arguments to get the appropriate weights 
# Check the version of your raster package for normalizeWeights!
myextr = as.data.frame(extract(N, poly, weights=T, normalizeWeights=F))
#      value weight
# 1 69.48172   0.16
# 2 98.10323   0.08
# 3 50.54667   0.61
# 4 78.71476   0.99
# 5 88.21990   0.17
# 6 93.66912   0.17
# 7 52.05317   0.87
# 8 83.05608   0.85
# 9 93.91854   0.43

# compute your weighted sum
mywsum = sum(myextr$value * myextr$weight)
# [1] 314.9164

* 编辑前

对于 raster v2.3-12,参数 normalizeWeights 显然不存在,因此我不得不手动完成这项工作,并警告栅格的分辨率与多边形的大小相比(即,需要一个完全封闭的单元,否则就需要进行调整)。因此,此答案下方的前两条评论。

Ps2 = Polygons(list(Polygon(cbind(c(0.5,2.8,3.7,1.2,0.5), c(2,3.4,1.2,0.2,2)))), 
               ID = "a")
SPs2 = SpatialPolygons(list(Ps2))
poly2 = SpatialPolygonsDataFrame(SPs2, data.frame(onecol = c("one"),
                                                  row.names = c("a")))
plot(poly2, add=T, border='blue', lty=3)

myextr2 = as.data.frame(extract(N, poly2, weights=T))

# compute the weight (cells fully enclosed = 1)
myextr2$weight2 = myextr2$weight/max(myextr2$weight)
#       value      weight weight2
# 1  69.48172 0.027777778    0.16
# 2  98.10323 0.013888889    0.08
# 3  50.54667 0.105902778    0.61
# 4  78.71476 0.171875000    0.99
# 5  88.21990 0.029513889    0.17
# 6  93.66912 0.052083333    0.30
# 7  52.05317 0.173611111    1.00
# 8  83.05608 0.173611111    1.00
# 9  93.91854 0.090277778    0.52
# 10 94.52795 0.003472222    0.02
# 11 78.31402 0.107638889    0.62
# 12 79.67737 0.048611111    0.28
# 13 68.22573 0.001736111    0.01

mywsum2 = sum(myextr2$value * myextr2$weight2)
# [1] 428.2086

这表明 raster 包已经很棒了,但仍在改进:-D

[警告:尝试将 normalizeWeightsraster v2.3-12 一起使用,它既没有崩溃也没有抛出错误,但没有完成这项工作,因此请注意并更新您的 raster 版本!]