R 栅格包中的提取和重采样函数：面积加权值

Question

我在其他论坛也发了这个帖子，但因为我真的需要回复，所以我再发一次。

我在 R 中工作并想计算从栅格的相交单元格派生的多边形的值。该值应考虑每个相交单元格的权重。当我尝试使用样本栅格和多边形运行“提取”函数时，我得到的权重与我手动计算的权重不同，从而导致最终值不同。

这是我的示例代码：

require(raster)
r <- raster(nrow=2, ncol=2, xmn=-180, xmx=60, ymn=-30, ymx=90)   
r[] <- c(1,2,4,5)    
s <- raster(xmn=-120, xmx=-40, ymn=20, ymx=60, nrow=1, ncol=1)    
s.pl <- as(s, 'SpatialPolygons')    
w <- raster::extract(r, s.pl, method="simple",weights=T, normalizeWeights=F)    
mean.value <- raster::extract(r, s.pl, method="simple",weights=T, fun=mean) 

我得到的值是 2.14，但根据单元格的实际权重，它应该是 2。更具体地说，对于与不同单元格相交的多边形的每个部分，数据是：

 Area Value
 1800 1
  600 2
  600 4
  200 5

所以基于上面的多边形的最终值应该是2。

可能是因为经纬度的投影吗？但即使我以米为单位分配投影，我也会得到相同的结果。如何获得我感兴趣的 2 的值？我也尝试过“重新采样”功能，但我也得到了不同的结果。

我的最终目标是创建一个与原始栅格具有不同分辨率和范围的新栅格，并根据与新栅格像元相交的原始栅格像元的权重分配值。但似乎重新采样和提取功能都没有给出预期的结果。

Answer 1

假设我们有一个栅格A和两个[B, C]非矩形的 SpatialPolygon 对象（在本例中为六边形）。为了演示目的，六边形的中心B被定义为栅格的中心A（见下面的左图）。六边形C沿水平轴向右移动。

require(raster)
require(scales)

A    <- raster(nrow=2, ncol=2, xmn=-180, xmx=180, ymn=-180, ymx=180)
A[]  <- c(1,2,4,5)    
A.pl <- as(A, 'SpatialPolygons')
B    <- SpatialPolygons(list(Polygons(list(Polygon(cbind(c(0, 100, 100, 0, -100, -100, 0), 
                                                         c(100, 50, -50, -100, -50, 50, 100)))), 'B')))
C    <- SpatialPolygons(list(Polygons(list(Polygon(cbind(c(40, 140, 140, 40, -60, -60, 40), 
                                                         c(100, 50, -50, -100, -50, 50, 100)))), 'C')))

对象B

由于六边形B位于中心，因此权重应全部等于 0.25。我们可以很容易地从图中得出六边形的面积为 30000（想象一个六边形适合的正方形（40000）并减去 2 个矩形（-10000），每个矩形由您必须切除的 4 个角中的 2 个组成） 。因此，每个交叉区域的大小为 7500 并且7500/30000 = 0.25

# get intersections
intsct.B <- raster::intersect(B, A.pl)
intsct.C <- raster::intersect(C, A.pl)

### B
area.B <- B@polygons[[1]]@area

weights <- unlist(lapply(intsct.B@polygons, function(x) {
  slot(x, 'area')/area.B
}))
weights
> [1] 0.25 0.25 0.25 0.25

现在我们获取每个相交多边形所在单元格的值并计算平均值。

vals    <- unlist(lapply(intsct.B@polygons, function(x) { 
  extract(A, data.frame(t(slot(x, 'labpt'))))
}))

sum(weights * vals)
> [1] 3

正如我们所期望的， 的平均值c(1, 2, 4, 5)是3。

对象C

现在让我们对对象做同样的事情C

### C
area.C <- C@polygons[[1]]@area

weights <- unlist(lapply(intsct.C@polygons, function(x) {
  slot(x, 'area')/area.C
}))
weights
> [1] 0.13 0.37 0.13 0.37

vals    <- unlist(lapply(intsct.C@polygons, function(x) { 
  extract(A, data.frame(t(slot(x, 'labpt'))))
}))

sum(weights * vals)
> [1] 3.24

同样，正如我们所期望的，平均值更大（因为值为 2 和 5 的单元格的权重更高）。另外，由于我们仅沿一个轴移动六边形，因此 2 个权重出现两次是有意义的。

具有较多像元的栅格

B下图显示了(左侧) 和C(rhs) 与4x4值为 的栅格的交集c(1:8, 10:17)。因为B有 12 个交集，而 为8 个。再次注意，由于对称性，C的平均值恰好为 9。B

这应该适用于任何SpatialPolygons对象。请务必对您扔入的对象使用相同的 CRS intersect。