如何将栅格中的面积加权和提取到 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 来实现它。我尝试将栅格表示为 SpatialPixels
、SpatialPoints
和 SpatialGrid
并使用 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
[警告:尝试将 normalizeWeights
与 raster v2.3-12
一起使用,它既没有崩溃也没有抛出错误,但没有完成这项工作,因此请注意并更新您的 raster
版本!]
我有一个栅格(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 来实现它。我尝试将栅格表示为 SpatialPixels
、SpatialPoints
和 SpatialGrid
并使用 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
[警告:尝试将 normalizeWeights
与 raster v2.3-12
一起使用,它既没有崩溃也没有抛出错误,但没有完成这项工作,因此请注意并更新您的 raster
版本!]