如何從R中的柵格提取面積加權總和爲多邊形？

Question

我有一個柵格（類「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此命令，收到以下：

「警告消息： 在。本地（X，Y，。 ..）： 「fun」已更改爲「mean」;其他函數在「weights = TRUE」時不能使用「

R當使用權重時強制更改爲fun = mean。我需要使用權重來說明落在多邊形中的柵格比例，但平均值不會讓我計算出我想要的數量。我需要找到一種方法來將光柵中的計數N加起來，這是由多邊形內的光柵面積（加權總和不是加權平均值）加權得出的。有誰知道我能做到這一點？

我很高興嘗試使用其他空間數據類來獲得此功能。我試圖用SpatialPixels，SpatialPoints和SpatialGrid來表示柵格，並使用sp::over()函數，但不知道如何計算我需要的加權和。非常感謝幫助！

Answer 1

首先，請重複的例子，這將幫助我們來幫助你;-)

*編輯 這是有效的版本raster v2.5-8但不是爲raster v2.3-12！

的weights解釋如下：

如果TRUE和normalizeWeights = FALSE，則該函數返回時，對於每個 多邊形，與小區值的矩陣和的 每個單元的近似分數即由多邊形覆蓋（四捨五入爲1/100）。如果TRUE爲 且normalizeWeights = TRUE，則將權重標準化，使得它們的 加起來爲1。權重可以用於平均;看例子。 此選項很有用（但慢），如果多邊形小 相對於柵格的像元尺寸*對象

所以，你可以用它來計算你加權和：

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！]