问如何从光栅中提取面积加权和为R中的多边形？
EN

Stack Overflow用户

提问于 2017-02-20 23:31:15

回答 1查看 1.9K关注 0票数 1

我有一个栅格(属于“RasterLayer”类)，它可以计数人，N.我有代表管理边界的多边形(类"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中运行此命令时，我会收到以下命令：

“警告信息:在.local(x，y，.)中：”weights=TRUE“改为”weights=TRUE“；”weights=TRUE“时不能使用其他函数

当使用权重时，r强制更改为fun=mean。我需要使用权重，以说明栅格在多边形中的比例，但是一个平均值不会让我计算我想要的数量。我需要找到一种方法，把栅格中的计数N加起来，用位于多边形内的栅格面积来加权(加权和而不是加权平均值)。有人知道我怎么能做到这一点吗？

我很高兴尝试使用其他空间数据类来实现这一目标。我尝试将光栅表示为SpatialPixels、SpatialPoints和SpatialGrid，并使用sp::over()函数，但无法计算出所需的加权和。帮助是非常感谢的！

r-raster

TDSQL-CMySQL免费体验

云原生数据库TDSQL-C MySQL 2核4G 免费体验15天

回答 1

Stack Overflow用户

回答已采纳

发布于 2017-02-21 00:30:30

首先，请举出一个可复制的例子，这将有助于我们帮助你；)

*编辑--这对于版本raster v2.5-8有效，但不适用于raster v2.3-12__！

weights的解释如下：

如果为TRUE和normalizeWeights=FALSE，则该函数为每个多边形返回一个矩阵，该矩阵包含单元格值和每个单元格被多边形覆盖的近似分数(四舍五入为1/100)。如果为TRUE和normalizeWeights=TRUE，则对权重进行规范化，使它们相加为1。权重可用于平均；参见示例。如果多边形相对于Raster*对象的单元格大小很小，则此选项可能有用(但速度慢)。

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

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