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

Question

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

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

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

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

	Population	Value	Frac. included	Pop. included	Weight	Result
Yellow	24	0.8	0.25	6	0.1875	0.15
Orange	16	0.4	0.5	8	0.25	0.10
Blue	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, empty_raster, "pop_included", fun = "sum")
11. 将第一个栅格除以第二个栅格得到我想要的：

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

任何帮助将不胜感激！

Answer 1

示例数据：

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",drop=FALSE], 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")], z[, "cell", drop=FALSE], 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", cex=.8, col="red", halo=TRUE)

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

Answer 2

这是快速且可扩展的：

library(data.table)
library(terra)

# make the 3 polygons with radius = 5km
center_points <- data.frame(lon = c(0.5, 0.65, 1),
                            lat = c(0.75, 0.65, 1),
                            Population = c(16, 18, 24),
                            Value = c(0.4, 0.1, 0.8))
polygon <- vect(center_points, crs = "EPSG:4326")
polygon <- buffer(polygon, 5000)

# make the raster
my_raster <- rast(nrow = 3, ncol = 3, xmin = 0, xmax = 1.5, ymin = 0, ymax = 1.5, crs = "EPSG:4326")
my_raster[] <- 0 # set the value to 0 for now

# find the fractions of cells in each polygon 
# "cells" gives you the cell ID and "weights" (or "exact") gives you the cell fraction in the polygon
# using "exact" instead of "weights" is more accurate
my_Table <- extract(my_raster, polygon, cells = TRUE, weights = TRUE)
setDT(my_Table) # convert to datatable

# merge the polygon attributes to "my_Table"
poly_Table <- setDT(as.data.frame(polygon))
poly_Table[, ID := 1:nrow(poly_Table)] # add the IDs which are the row numbers
merged_Table <- merge(my_Table, poly_Table, by = "ID")

# find Frac_included
merged_Table[, Frac_included := weight / sum(weight), by = ID]

# find Pop_included
merged_Table[, Pop_included := Frac_included * Population]

# find Weight, to avoid confusion with "weight" produced above, I call this "my_Weight"
merged_Table[, my_Weight := Pop_included / sum(Pop_included), by = cell]

# final results
Result <- merged_Table[, .(Result = sum(Value * my_Weight)), by = cell]

# add the values to the raster
my_raster[Result$cell] <- Result$Result
plot(my_raster)