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

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 的值？我也尝试了“重采样”功能，但也得到了不同的结果。

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

Answer 1

这是我根据this帖子的回复设法做到的。

require(raster)
require(rgeos)
r <- raster(nrow=2, ncol=2, xmn=-180, xmx=60, ymn=-30, ymx=90)    
r[] <- c(1,2,4,5)    
r <- stack(r, r*2, r^2)
s <- raster(xmn=-120, xmx=-40, ymn=20, ymx=60, nrow=1, ncol=1)    
s.pl <- as(s, 'SpatialPolygons')    
r.s <- as(r, 'SpatialPolygonsDataFrame')
pi1 <- gIntersection(r.s, s.pl, byid = T)
areas1 <- data.frame(area=sapply(pi1@polygons, FUN=function(x) {slot(x, 'area')}))
row.names(areas1) <- sapply(pi1@polygons, FUN=function(x) {slot(x, 'ID')})
areas1$Pol.old <- as.numeric(vapply(strsplit(rownames(areas1), " "), `[`, 1, FUN.VALUE=character(1)))
areas1$pol.new <- as.numeric(vapply(strsplit(rownames(areas1), " "), `[`, 2, FUN.VALUE=character(1)))
f <- r.s@data
seqs <- match(areas1$Pol.old, rownames(f))
ar <- cbind(areas1, f[seqs,])
ar[,-(1:3)] <- ar[,-(1:3)]*ar$area
f <- aggregate.data.frame(ar, by=list(ar$pol.new), FUN=sum)
f[,-(1:4)] <- f[,-(1:4)]/f$area  
ar.v <- as.matrix(f[, -c(1:4)])
s2 <- stack(s)
s1 <- setValues(s2, ar.v)

如果有人可以建议更好和/或更快的代码，请告诉我，因为我不太喜欢我的方法。

Answer 2

假设我们有一个栅格 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。