R：栅格化具有复杂权重的多边形
R: Rasterizing polygons with complicated weighting
想象一下地球表面有一个规则的 0.5° 网格。该网格的 3x3 子集如下所示。作为我正在使用的程式化示例，假设我有三个多边形——黄色、橙色和蓝色——为了简单起见，它们的面积都是 1 个单位。这些多边形具有属性 Population 和 Value，您可以在图例中看到它们：
我想将这些多边形转换为 0.5° 栅格（具有全局范围），其值基于多边形的加权平均值。棘手的部分是我想根据它们的 included 人口来加权多边形的值。
我知道——理论上——我想做什么，下面已经为中心网格完成了。
- 将 Population 乘以 Included（包含在网格单元中的多边形区域）得到 Pop。包括。 （假设人口在整个面内均匀分布，这是可以接受的。）
- 将每个多边形的 Included_pop 除以所有多边形的总和 Included_pop (32) 得到权重。
- 将每个多边形的值乘以权重得到结果。
- 对所有多边形的结果求和以获得中心网格单元 (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")
- 计算每个多边形的面积：
old_polygons$area <- as.numeric(st_area(old_polygons))
- 将全局 0.5° 网格生成为某种空间对象。
- 通过网格分割多边形，生成
new_polygons。
- 计算新多边形的面积：
new_polygons$new_area <- as.numeric(st_area(new_polygons))
- 计算每个新多边形包含的分数：
new_polygons$frac_included <- new_polygons$new_area / new_polygons$old_area
- 计算新多边形中的“包含人口”：
new_polygons$pop_included <- new_polygons$pop * new_polygons$frac_included
- 为每个多边形计算一个新属性，该属性只是它们的值乘以它们包含的人口。
new_polygons$tmp <- new_polygons$Value * new_polygons$frac_included
- 为后续步骤设置空栅格：
empty_raster <- raster(nrows=360, ncols=720, xmn=-180, xmx=180, ymn=-90, ymx=90)
- 通过在每个网格单元内将这个新属性加在一起来栅格化多边形。
tmp_raster <- rasterize(new_polygons, empty_raster, "tmp", fun = "sum")
- 创建另一个栅格，它只是每个网格单元中的总人口：
pop_raster <- rasterize(new_polygons, empty_raster, "pop_included", fun = "sum")
- 将第一个栅格除以第二个得到我想要的：
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",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)
这可能无法很好地扩展到大型数据集，但您或许可以分块进行。
这是快速且可扩展的：
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)
想象一下地球表面有一个规则的 0.5° 网格。该网格的 3x3 子集如下所示。作为我正在使用的程式化示例，假设我有三个多边形——黄色、橙色和蓝色——为了简单起见，它们的面积都是 1 个单位。这些多边形具有属性 Population 和 Value，您可以在图例中看到它们：
我想将这些多边形转换为 0.5° 栅格（具有全局范围），其值基于多边形的加权平均值。棘手的部分是我想根据它们的 included 人口来加权多边形的值。
我知道——理论上——我想做什么，下面已经为中心网格完成了。
- 将 Population 乘以 Included（包含在网格单元中的多边形区域）得到 Pop。包括。 （假设人口在整个面内均匀分布，这是可以接受的。）
- 将每个多边形的 Included_pop 除以所有多边形的总和 Included_pop (32) 得到权重。
- 将每个多边形的值乘以权重得到结果。
- 对所有多边形的结果求和以获得中心网格单元 (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")
- 计算每个多边形的面积：
old_polygons$area <- as.numeric(st_area(old_polygons)) - 将全局 0.5° 网格生成为某种空间对象。
- 通过网格分割多边形，生成
new_polygons。 - 计算新多边形的面积：
new_polygons$new_area <- as.numeric(st_area(new_polygons)) - 计算每个新多边形包含的分数：
new_polygons$frac_included <- new_polygons$new_area / new_polygons$old_area - 计算新多边形中的“包含人口”：
new_polygons$pop_included <- new_polygons$pop * new_polygons$frac_included - 为每个多边形计算一个新属性，该属性只是它们的值乘以它们包含的人口。
new_polygons$tmp <- new_polygons$Value * new_polygons$frac_included - 为后续步骤设置空栅格：
empty_raster <- raster(nrows=360, ncols=720, xmn=-180, xmx=180, ymn=-90, ymx=90) - 通过在每个网格单元内将这个新属性加在一起来栅格化多边形。
tmp_raster <- rasterize(new_polygons, empty_raster, "tmp", fun = "sum") - 创建另一个栅格，它只是每个网格单元中的总人口：
pop_raster <- rasterize(new_polygons, empty_raster, "pop_included", fun = "sum") - 将第一个栅格除以第二个得到我想要的：
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",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)
这可能无法很好地扩展到大型数据集，但您或许可以分块进行。
这是快速且可扩展的：
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)