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

发布于05月06日

我从WorldPop下载了100 m分辨率的人口网格数据，通过跨投影和裁剪，得到了一个城市的人口空间分布网格数据，此外还有该市县级的行政边界形状文件数据，目的是根据当年的人口普查(县级)数据进行人口网格修正.

##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,])

我计划分别裁剪每个县的人口网格，然后通过计算其修正因子来进行修正，最后合并这些网格.不过，这种方式很耗时，不知道有没有更简单的方法？

推荐答案

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

创建新版本的WorldPop格栅，其中每个像元代表其所在的WorldPop县人口，例如pop1;
使用原始WorldPop格栅作为参考，使用rasterize()创建具有相应县$census_data人口值的县形状文件的格栅版本，例如county_r;
使用county_r和pop1定义比例以调整pop中的单元格值

请注意，此示例使用tigris个数据的默认CRS，并且为了速度和简单性，其分辨率低于WorldPop数据.

Step 0: Load packages and create example data

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

Step 1: Create population by county raster from WorldPop data

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

Step 2: 100 county shapefile with county$census population values to corresponding cells

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")

Step 3: Adjust original WorldPop values using county_r and pop1

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