用更新的聚合多边形值替换栅格中的所有像元值（即更新网格人口估计）

Question

目标：根据 2020 年网格人口栅格和 2023 年面人口估计生成 2023 年更新的网格人口栅格。

为此，我采取了尝试使用按像元编号存储在表中的新像元值来更新栅格中的所有像元值的方法。

我有一个大型栅格，其中包含 2020 年的网格人口估计值、县多边形以及按县划分的 2023 年人口表。我想以相同的分布/密度更新 2020 年栅格中的像元值，但反映 2023 年的人口。

我已从以下方面开始：

1. 使用

   exactextractr::exact_extract

   对每个县多边形内的 2020 年像元值（人口）求和。
2. 使用

   terra:extract

   提取 2020 年栅格中的所有像元值
3. 连接1和2来计算每个单元格的县人口比例（即单元格值/县人口）。
4. 将 2023 年县人口估计值加入到 3，并计算反映 2023 年人口的新单元格值（2023 年县人口 * 每个单元格中县人口的比例）
5. 我似乎不知道如何将新的像元值放入 2020 年栅格中。

我有一个包含县 ID、像元编号、旧像元值和新像元值的表格，但我不确定如何使用新值更新 2020 年栅格。如果这是一个数据框，我将按单元格编号加入，或者直接分配新值，假设行顺序相同。我看过以下 SO 文章（以及更多文章），但它们似乎只解决了条件值替换。

• 使用基于 R terra 中的多边形字段的向量或值矩阵更新栅格值
• 用另一个栅格中的值替换栅格中的值
• https://rdrr.io/cran/raster/man/replacement.html
• https://rdrr.io/cran/terra/man/subst.html

下面是一个可重现的示例（感谢 Zonal Statistics R（栅格/多边形） 提供通用栅格和多边形）。细胞数量在

extract

步骤中的显示方式与真实数据的显示方式不同，但在其他方面非常相似。理想情况下，我想通过单元格号码加入，以便安心，而不是依赖于行顺序。

library(terra)
library(exactextractr)
library(tidyverse)


# Create raster
ras <- raster(nrows=100, ncols=80, xmn=0, xmx=1000, ymn=0, ymx=800)
val <- runif(ncell(ras))
values(ras) <- val

# Create polygon within raster extent
xym <- cbind(runif(3,0,1000), runif(3,0,800))
p <- Polygons(list(Polygon(xym)),1)
sp <- SpatialPolygons(list(p))
spdf <- SpatialPolygonsDataFrame(sp, data=data.frame(1))

# Mask raster to within polygon
r2 <- mask(ras, spdf)

plot(r2)
lines(spdf)

#Create table of for new polygon value
new_tbl <- 
  tibble(X1 = 1,
         new_sum = 10000)

#1. Use ```extact_extractr::exact_extract``` to sum the 2020 cell values within each polygon.

poly_pop <- exact_extract(
  r2, 
  spdf, 
  fun = 'sum',
  weights = 'area',
  append_cols = T
) 

#2. Use ```terra:extract``` to extract all cell values in the old raster

r2_cell_values <- 
  terra::extract(
  r2, # raster file
  spdf, # polygons
  df = TRUE, 
  cells = TRUE, # Include cell numbers
  ID = TRUE,  # Include polygon IDs
  exact = TRUE # Include the proportion of the cell covered by the polygon
)

#3. Joined 1 and 2 to calculate the proportion of the value for each cell (i.e. cell value / polygon sum).

cell_props <-
r2_cell_values %>%
  filter(!is.na(layer)) %>% # Remove NA values
  left_join(poly_pop, by = c("ID" = "X1")) %>%
  mutate(
    cell_prop = layer / sum
  )

#4. Joined thenew values to 3 and calculated new cell values reflecting the mew polygon value (new polygon value * proportion of old polygon value in each cell)

new_cell_tbl <- 
cell_props %>%
  left_join(new_tbl, by = c("ID" = "X1")) %>%
  mutate(new_cell_value = cell_prop * new_sum)

#5. I can't seem to figure out how to get the new cell values into the 2020 raster. 

x <- rasterize(spdf, r2, 1)
cellStats(x, 'sum') == nrow(new_cell_tbl) # Check whether number of values in r2 and new_cell_tbl are the same

values(r2) <- new_cell_tbl$new_cell_value

我最终得到了这个错误消息：

Error in setValues(x, value) :  length(values) is not equal to ncell(x), or to 1

，即使我尝试测试之前步骤中的栅格单元数量和新向量的长度是否相同。

我对此相当陌生，因此任何关于如何根据 2023 年管理区域值更新 2020 年单元格的指导将不胜感激。谢谢！

Answer 1

示例数据

library(terra)
pop <- rast(system.file("ex/elev.tif", package="terra"))
names(pop) <- "pop2020"
adm <- vect(system.file("ex/lux.shp", package="terra"))

计算“2020”每个行政单位的相对人口

adm <- extract(pop, adm, "sum", na.rm=TRUE, bind=TRUE)
pop_agg <- rasterize(adm, pop, "pop2020")
rel_pop <- pop / pop_agg

使用相对分布计算“2030”的栅格像元计数

#example data
adm$pop2030 <- adm$pop2020 * (1:12) / 2
# solution
new_pop <- rasterize(adm, pop, "pop2030") * rel_pop