我从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,])
我计划分别裁剪出每个县的人口栅格,然后通过计算其校正因子进行校正,最后合并栅格。不过这样比较费时间,不知道有没有更简单的方法呢?
如果您的目的是使用 WorldPop 中的县人口与 County$census_data 之间的比率来调整您的 WorldPop 像元值,您可以使用栅格算术来完成此操作:
rasterize()
创建具有相应 County$census_data 人口值的县形状文件的栅格版本,例如county_r;请注意,此示例使用
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 an example 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