基于2023面状人口估算更新2020栅格人口数据的技术求助

问题说明 #

目标是基于2020年栅格人口估算数据和2023年县级面状人口结果，生成2023年更新后的栅格人口数据，要求保留原有人口分布密度。已完成前4步计算得到包含县ID、单元编号、新旧值的表格，但执行values(r2) <- new_cell_tbl$new_cell_value时出现错误：

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

即使验证过有效单元数与新值数量一致。

错误原因 #

r2是经过mask处理后的栅格，其总单元数包含mask范围外的NA单元，而new_cell_tbl仅包含mask范围内的有效单元。values()函数要求赋值的向量长度必须等于栅格总单元数，两者长度不匹配导致报错。

解决方案 #

通过单元编号精准定位赋值，仅更新mask范围内的有效单元，保留原栅格的NA结构。修改后的关键代码如下：

# 复制原栅格结构，保留原NA值
updated_raster <- r2

# 利用单元编号（new_cell_tbl$cell）将新值填充到对应位置
updated_raster[new_cell_tbl$cell] <- new_cell_tbl$new_cell_value

# 验证结果：检查新栅格总和是否等于县级新人口值
cellStats(updated_raster, 'sum', na.rm = TRUE) == new_tbl$new_sum

完整可复现代码 #

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

# 创建示例栅格
ras <- rast(nrows=100, ncols=80, xmn=0, xmx=1000, ymn=0, ymx=800)
values(ras) <- runif(ncell(ras))

# 创建示例多边形
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))

# 裁剪栅格到多边形范围
r2 <- mask(ras, spdf)

# 县级新人口数据
new_tbl <- tibble(X1 = 1, new_sum = 10000)

# 1. 汇总县级栅格人口总和
poly_pop <- exact_extract(
  r2, 
  spdf, 
  fun = 'sum',
  weights = 'area',
  append_cols = T
) 

# 2. 提取栅格单元值及编号
r2_cell_values <- terra::extract(
  r2, 
  spdf, 
  df = TRUE, 
  cells = TRUE, 
  ID = TRUE,  
  exact = TRUE
)

# 3. 计算单元占县级人口比例
cell_props <- r2_cell_values %>%
  filter(!is.na(layer)) %>%
  left_join(poly_pop, by = c("ID" = "X1")) %>%
  mutate(cell_prop = layer / sum)

# 4. 计算单元新人口值
new_cell_tbl <- cell_props %>%
  left_join(new_tbl, by = c("ID" = "X1")) %>%
  mutate(new_cell_value = cell_prop * new_sum)

# 5. 正确更新栅格
updated_raster <- r2
updated_raster[new_cell_tbl$cell] <- new_cell_tbl$new_cell_value

# 可视化验证
plot(updated_raster)
lines(spdf)

补充提示 #

• 建议统一使用terra包的rast()创建栅格，避免混用raster包函数，减少兼容性问题
• 多县级单元场景下，只需确保new_cell_tbl包含所有待更新单元的编号和新值即可，赋值逻辑不变

内容的提问来源于stack exchange，提问作者RaBe