R 栅格包中的提取和重采样函数：面积加权值
Extract and Resample Functions in R raster package: Area-Weighted Values
我也在 other forum 中做了这个 post，但是因为我真的需要回复，所以我再次 post 在这里。
我在 R 中工作，想计算从栅格的相交像元派生的多边形的值。该值应考虑每个相交单元格的权重。当我尝试使用样本栅格和多边形 运行 "extract" 函数时，我得到的权重与我手动计算的权重不同，导致最终值不同。
这是我的示例代码：
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。
难道是因为lat/lon中的投影？但即使我以米为单位分配投影，我也会得到相同的结果。如何获得我感兴趣的 2 的值？我也尝试使用 "resample" 函数，但我也得到了不同的结果。
我的最终目标是创建一个与原始栅格具有不同分辨率和范围的新栅格，并根据与新栅格像元相交的原始栅格像元的权重分配值。但似乎 resample 和 extract 函数都没有给出预期的结果。
这是我根据 this post.
的回复设法做到的
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)
如果有人可以建议更好的 and/or 更快的代码，请告诉我，因为我不太喜欢我的方法。
让我们假设我们有一个栅格 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（右侧）与值为 c(1:8, 10:17) 的 4x4 栅格的交点。 B 有 12 个交点，C 有 8 个交点。请再次注意，由于对称性，B 的均值恰好为 9。
这应该适用于任何 SpatialPolygons 对象。请务必对您放入 intersect.
的对象使用相同的 CRS
我也在 other forum 中做了这个 post，但是因为我真的需要回复，所以我再次 post 在这里。
我在 R 中工作，想计算从栅格的相交像元派生的多边形的值。该值应考虑每个相交单元格的权重。当我尝试使用样本栅格和多边形 运行 "extract" 函数时，我得到的权重与我手动计算的权重不同，导致最终值不同。
这是我的示例代码：
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。
难道是因为lat/lon中的投影？但即使我以米为单位分配投影，我也会得到相同的结果。如何获得我感兴趣的 2 的值？我也尝试使用 "resample" 函数，但我也得到了不同的结果。
我的最终目标是创建一个与原始栅格具有不同分辨率和范围的新栅格，并根据与新栅格像元相交的原始栅格像元的权重分配值。但似乎 resample 和 extract 函数都没有给出预期的结果。
这是我根据 this post.
的回复设法做到的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)
如果有人可以建议更好的 and/or 更快的代码，请告诉我，因为我不太喜欢我的方法。
让我们假设我们有一个栅格 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（右侧）与值为 c(1:8, 10:17) 的 4x4 栅格的交点。 B 有 12 个交点，C 有 8 个交点。请再次注意，由于对称性，B 的均值恰好为 9。
这应该适用于任何 SpatialPolygons 对象。请务必对您放入 intersect.