R：光栅化具有复杂权重的多边形

如何解决R：光栅化具有复杂权重的多边形

想象一下地球表面有一个规则的 0.5° 网格。该网格的 3x3 子集如下所示。作为我正在处理的一个程式化的例子，假设我有三个多边形——黄色、橙色和蓝色——为了简单起见，它们的面积都是 1 个单位。这些多边形具有属性 Population 和 Value，您可以在图例中看到：

我想将这些多边形转换为 0.5° 栅格（具有全局范围），其值基于多边形的加权平均值。棘手的部分是我想根据多边形的值而不是根据它们的人口来加权多边形的值，而是根据它们包含的人口。

我知道——理论上——我想做什么，下面已经为中心网格单元做了。

1. 将人口乘以包含（网格单元中包含的多边形面积）以获得人口。包括。 （假设人口均匀分布在整个多边形中，这是可以接受的。）
2. 将每个多边形的 Included_pop 除以所有多边形的 Included_pop 总和 (32) 以获得权​​重。
3. 将每个多边形的值乘以权重以获得结果。
4. 对所有多边形的结果求和以获得中心网格单元的值 (0.31)。

	人口	价值	碎裂。包括	流行音乐。包括	重量	结果
黄色	24	0.8	0.25	6	0.1875	0.15
橙色	16	0.4	0.5	8	0.25	0.10
蓝色	18	0.1	1	18	0.5625	0.06
				32		0.31

我有一个想法如何在 R 中实现这一点，如下所述。在可能的情况下，我已经填写了我认为会做我想做的代码。我的问题：我该如何进行第 2 步和第 3 步？或者有没有更简单的方法来做到这一点？如果您想尝试一下，我已经将 old_polygons 作为 .rds 文件 here 上传。

library("sf")
library("raster")

1. 计算每个多边形的面积：old_polygons$area <- as.numeric(st_area(old_polygons))
2. 将全局 0.5° 网格生成为某种空间对象。
3. 按网格分割多边形，生成 new_polygons。
4. 计算新多边形的面积：new_polygons$new_area <- as.numeric(st_area(new_polygons))
5. 计算每个新多边形的分数：new_polygons$frac_included <- new_polygons$new_area / new_polygons$old_area
6. 计算新多边形中的“包含人口”：new_polygons$pop_included <- new_polygons$pop * new_polygons$frac_included
7. 为每个多边形计算一个新属性，即其值乘以包含的人口。 new_polygons$tmp <- new_polygons$Value * new_polygons$frac_included
8. 为后续步骤设置一个空栅格：empty_raster <- raster(nrows=360,ncols=720,xmn=-180,xmx=180,ymn=-90,ymx=90)
9. 通过将每个网格单元内的新属性相加来栅格化多边形。 tmp_raster <- rasterize(new_polygons,empty_raster,"tmp",fun = "sum")
10. 创建另一个栅格，它只是每个网格单元中的总人口：pop_raster <- rasterize(new_polygons,"pop_included",fun = "sum")
11. 将第一个栅格除以第二个得到我想要的：

output_raster <- empty_raster
values(output_raster) <- getValues(tmp_raster) / getValues(pop_raster)

任何帮助将不胜感激！

解决方法

示例数据：

library(terra)
f <- system.file("ex/lux.shp",package="terra")
v <- vect(f)
values(v) <- data.frame(population=1:12,value=round(c(2:13)/14,2))
r <- rast(ext(v)+.05,ncols=4,nrows=6,names="cell")

说明数据

p <- as.polygons(r)
plot(p,lwd=2,col="gray",border="light gray")
lines(v,col=rainbow(12),lwd=2)
txt <- paste0(v$value," (",v$population,")") 
text(v,txt,cex=.8,halo=TRUE)

解决方案：

# area of the polygons
v$area1 <- expanse(v)
# intersect with raster cell boundaries
values(r) <- 1:ncell(r) 
p <- as.polygons(r)
pv <- intersect(p,v)

# area of the polygon parts
pv$area2 <- expanse(pv)
pv$frac <- pv$area2 / pv$area1

现在我们只使用带有多边形属性的 data.frame 来计算多边形覆盖加权人口加权值。

z <- values(pv)
a <- aggregate(z[,"frac",drop=FALSE],z[,"cell",sum)
names(a)[2] <- 'fsum'
z <- merge(z,a)
z$weight <- z$population * z$frac / z$fsum
z$wvalue <- z$value * z$weight
b <- aggregate(z[,c("wvalue","weight")],sum)
b$bingo <- b$wvalue / b$weight

将值分配回栅格单元

x <- rast(r)
x[b$cell] <- b$bingo

检查结果

plot(x)
lines(v)
text(x,digits=2,halo=TRUE,cex=.9)
text(v,"value",col="red",halo=TRUE)

这可能无法很好地扩展到大型数据集，但您或许可以分块进行。