如何使用行政边界形状文件中的人口普查数据调整栅格数据集中的人口数据

Question

我从WorldPop下载了100m分辨率的人口栅格数据，通过转投影和裁剪，得到了一个城市的人口空间分布栅格数据，另外，我还有这个城市县级的行政边界shapefile数据，我的目的是根据当年的人口普查（县级）数据修正人口栅格。

##Mask and crop one county from the population data
popc <- trim(mask(pop, county[1, ]))
##Calculate the correction factor and make correction
POPc <- popc*(county[1, ]$census_data/global(popc, sum, na.rm = T)[1,])

我计划分别裁剪出每个县的人口栅格，然后通过计算其校正因子进行校正，最后合并栅格。不过这样比较费时间，不知道有没有更简单的方法呢？

Answer 1

如果您的目的是使用 WorldPop 中的县人口与 County$census_data 之间的比率来调整您的 WorldPop 像元值，您可以使用栅格算术来完成此操作：

1. 创建新版本的 WorldPop 栅格，其中每个像元代表其所属的 WorldPop 县人口，例如pop1;
2. 使用原始 WorldPop 栅格作为参考，使用

   rasterize()

   创建具有相应 County$census_data 人口值的县形状文件的栅格版本，例如county_r；
3. 使用 county_r 和 pop1 定义调整 pop
中单元格值的比率

请注意，此示例使用

tigris

数据中的默认 CRS，并且出于速度和简单性的考虑，其分辨率低于 WorldPop 数据。

第0步：加载包并创建示例数据

library(sf) # 1.0-16 
library(terra) # 1.7.71
library(tigris) # 2.1
library(dplyr) # 1.1.4

# Create example city by county sf data using tigris package
# Add ID to your actual data for left_join() in Step 1
county <- counties("Maine", cb = TRUE) %>%
  mutate(ID = 1:n())

# Get extent of county and use values to create a proxy WorldPop SpatRaster
st_bbox(county)

#      xmin      ymin      xmax      ymax 
# -71.08392  42.97776 -66.94989  47.45969

set.seed(1)
pop <- rast(resolution = 100/11113.9, nlyrs = 1,
            xmin = -71.08392, xmax = -66.94989, 
            ymin = 42.97776, ymax = 47.45969, 
            names = "population",
            crs = crs(county)) %>%
  init("cell") %>%
  crop(county, mask = TRUE)

pop[] <- sample(1:100, nrow(pop) * ncol(pop), replace = TRUE)

pop
# class       : SpatRaster
# dimensions  : 498, 459, 1  (nrow, ncol, nlyr)
# resolution  : 0.008997742, 0.008997742  (x, y)
# extent      : -71.08392, -66.95396, 42.97776, 47.45864  (xmin, xmax, ymin, ymax)
# coord. ref. : lon/lat NAD83 (EPSG:4269)
# source(s)   : memory
# name        : population
# min value   :          1
# max value   :        100

第 1 步：根据 WorldPop 数据按县栅格创建人口

pop1 <- extract(pop, county, fun = sum, na.rm = TRUE) %>%
  left_join(county, ., by = "ID") %>%
  rasterize(., pop, field = "population")

第 2 步：

rasterize()

县形状文件，其中包含对应单元格的县$人口普查人口值

county_r <- extract(pop, county, fun = sum, na.rm = TRUE) %>%
  left_join(county, ., by = "ID") %>%
  mutate(census_data = ceiling(population * 1.1)) %>% # For illustrative purposes, delete this line for your actual data
  rasterize(., pop, field = "census_data")

第3步：使用county_r和pop1调整原始WorldPop值

POPc <- pop * (county_r / pop1)

POPc
# class       : SpatRaster
# dimensions  : 498, 459, 1  (nrow, ncol, nlyr)
# resolution  : 0.008997742, 0.008997742  (x, y)
# extent      : -71.08392, -66.95396, 42.97776, 47.45864  (xmin, xmax, ymin, ymax)
# coord. ref. : lon/lat NAD83 (EPSG:4269)
# source(s)   : memory
# name        : population
# min value   :     1.1000
# max value   :   110.0018