r - R 栅格包中的提取和重采样函数 : Area-Weighted Values

Question

我在 other forum 发表了这篇文章同样，但由于我确实需要回复，所以我再次将其发布在这里。

我正在使用 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和两个 SpatialPolygon 对象 [B, C]不是矩形(在本例中为六边形)。为了演示目的，六边形的中心 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 (右)带有 4x4栅格值 c(1:8, 10:17) 。对于 B有 12 个交叉路口，C 8. 再次注意 B 的平均值由于对称性，正好是 9。

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