今天给大家介绍如何使用 R 语言基于 LandScan 全球人口栅格数据，按照王峤（2024）论文中的方法测算各城市多中心指标。

一、指标来源与计算原理

1.1 数据来源

本方法使用的人口数据为美国能源部橡树岭国家实验室提供的 LandScan 全球人口密度栅格数据，空间分辨率约为 1 km × 1 km（坐标参考系转换后约 820 m）。

1.2 指标定义

多中心指数（Polycentric Index）源自论文：

王峤. 空间结构、城市规模与中国城市的创新绩效 [J]. 中国工业经济, 2024.

其核心定义为：

指数 = 0：城市完全单中心，仅识别出一个中心；
指数越大（趋近于 1）：城市多中心特征越明显，次中心与主中心旗鼓相当。

1.3 计算流程

整个计算分为以下步骤：

坐标系转换：将 LandScan 栅格数据投影到等面积坐标系（Albers 投影），确保距离计算准确；
裁剪与掩膜：按城市行政边界裁剪栅格，获取该城市范围内的人口格点；
栅格转点：将栅格像元转为空间点数据；
构建点邻接：使用距离阈值（1000 m）构建邻接关系（边邻接规则，仅选取上下左右相邻格点），生成空间权重矩阵；
局部 Moran’s I 检验：计算每个格点的局部莫兰指数，识别统计显著的高-高聚类区域（HH 格点）；
聚类成中心：将 HH 格点按空间邻接分组，形成连续的”高密度区块”；
筛选有效中心：要求每个区块格点数 ≥ 3、总人口 ≥ 100,000；
计算多中心指数：按主中心（最大人口区块）与次中心之比计算。

1.4 参数设定说明

ANALYSIS_PARAMS <- list(
  sig_level = 0.05,   # 局部莫兰指数显著性阈值
  min_cells = 3,      # 中心最小格点数（过滤噪点）
  min_pop   = 100000, # 中心最小人口（单位：人）
  dist_nb   = 1000    # 邻接距离阈值（米）
)

关于 dist_nb = 1000 的设定：LandScan 数据经 Albers 投影后分辨率约为 820 m。相邻格点（上下左右）的距离约为 820 m，对角线方向约为 820 × √2 ≈ 1159 m。因此将阈值设为 1000 m，可精确选取边邻接格点，不会错误地把对角线邻居算进来。

二、环境配置与数据准备

2.1 加载 R 包

pacman::p_load(
  terra,     # 栅格数据读取、裁剪、掩膜
  sf,        # 矢量数据处理
  spdep,     # 邻接矩阵、空间权重矩阵、局部莫兰指数
  tidyverse, # 数据处理
  fs,        # 文件系统操作
  furrr,     # 基于 future 的并行 map
  tictoc     # 计时
)

2.2 全局路径与参数配置

# 等面积投影（Albers，适合中国范围计算）
mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +ellps=krass +units=m +no_defs"

PATH <- list(
  pop_tif    = path("pop-tif2"),           # 已重投影到 mycrs 的 tif 文件夹
  city_shp   = path("2021行政区划/市.shp"), # 2021 年市级行政区划
  output_csv = path("多中心指数_2000_2024.csv")
)

ANALYSIS_PARAMS <- list(
  sig_level = 0.05,
  min_cells = 3,
  min_pop   = 100000,
  dist_nb   = 1000
)

数据预处理提示：如果你的 LandScan tif 文件还是 WGS84 坐标系，需要先转换到 r mycrs：
dir.create("pop-tif2")
fs::dir_ls("pop-tif") %>%
  lapply(function(x){
    rast(x) %>%
      project(mycrs) %>%
      writeRaster(str_replace_all(x, "pop-tif", "pop-tif2"))
  })

三、单城市演示计算（北京市，2020 年）

3.1 读取城市边界与人口栅格

# 读取北京市行政边界
city_demo <- read_sf(PATH$city_shp) %>%
  st_make_valid() %>%
  st_transform(mycrs) %>%
  dplyr::filter(市 == "北京市")
city_demo

# 读取 2020 年人口栅格
pop_demo <- rast(path(PATH$pop_tif, "2020.tif"))
pop_demo

3.2 裁剪与栅格转点

# 裁剪并掩膜
r_clip <- pop_demo %>%
  crop(vect(city_demo)) %>%
  mask(vect(city_demo))

# 栅格转点，去掉 NA
pop_points <- as.points(r_clip) %>%
  st_as_sf() %>%
  rename(pop = 1) %>%
  drop_na(pop)

pop_points

3.3 构建空间权重矩阵

# 提取坐标
coords <- st_coordinates(pop_points)

# 构建距离邻接关系（边邻接，1000 m）
nb <- dnearneigh(coords, d1 = 0, d2 = ANALYSIS_PARAMS$dist_nb)
nb

# 转换为行标准化空间权重矩阵
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

参数说明

style = “W”：行标准化权重。若一个格点有 4 个邻居，则每个邻居权重 = 1/4；有 2 个邻居则权重 = 1/2。

zero.policy = TRUE：允许无邻居的孤立点（如城市边缘格点）存在；设为 FALSE 时遇到孤立点会报错，但在城市内部使用 FALSE 也是合理的选择。

3.4 局部 Moran’s I 检验

局部 Moran’s I 是用来找 “高值围着高值、低值围着低值” 的空间集聚热点的工具。在这里，它的唯一作用就是：找出人口高密度、且周围也是高密度的区域 → 也就是城市中心（HH 集聚区）。

# 计算局部莫兰指数
local_moran <- localmoran(pop_points$pop, lw, zero.policy = TRUE)
local_moran %>% as_tibble()

局部莫兰指数结果包含 5 列：

列名	含义
`Ii`	局部莫兰指数值
`E.Ii`	期望值
`Var.Ii`	方差
`Z.Ii`	标准化 Z 值
`Pr(z != E(Ii))`	双侧 p 值

3.5 空间类型分类（LISA 象限）

pop_points <- pop_points %>%
  mutate(
    p_value   = local_moran[, 5],
    sig       = p_value < ANALYSIS_PARAMS$sig_level,
    pop_med   = median(pop, na.rm = TRUE),
    pop_lag   = lag.listw(lw, pop),          # 空间滞后（邻居的加权平均人口）
    self_type = if_else(pop > pop_med, "H", "L"),
    nb_type   = if_else(pop_lag > pop_med, "H", "L"),
    cluster_type = case_when(
      sig & self_type == "H" & nb_type == "H" ~ "HH",
      TRUE ~ "其他"
    )
  )
pop_points

这里只关心 HH 类型（自身人口高 + 邻居人口高 + 统计显著），这类格点代表真正的人口密集聚集区，是城市中心的候选位置。

3.6 提取 HH 格点并聚类成中心

# 筛选 HH 格点
hh_points <- pop_points %>% filter(cluster_type == "HH")

# 对 HH 格点再次做空间邻接，将连续区域标为同一 cluster
if (nrow(hh_points) >= ANALYSIS_PARAMS$min_cells) {
  coords_hh <- st_coordinates(hh_points)
  nb_hh     <- dnearneigh(coords_hh, 0, ANALYSIS_PARAMS$dist_nb)
  hh_points$cluster_id <- n.comp.nb(nb_hh)$comp.id
}
pop_points

3.7 筛选有效中心并计算多中心指数

# 按 cluster 汇总，筛选有效中心
city_centers <- hh_points %>%
  group_by(cluster_id) %>%
  summarise(
    cell_count = n(),
    total_pop  = sum(pop, na.rm = TRUE),
    .groups    = "drop"
  ) %>%
  filter(
    cell_count >= ANALYSIS_PARAMS$min_cells,
    total_pop  >= ANALYSIS_PARAMS$min_pop
  )

# 计算多中心指数
if (nrow(city_centers) >= 1) {
  city_centers <- city_centers %>% arrange(desc(total_pop))
  main_pop     <- first(city_centers$total_pop)
  sub_pop      <- sum(city_centers$total_pop[-1])
  total_pop    <- main_pop + sub_pop
  poly_index   <- sub_pop / total_pop

  message("城市：北京市 | 年份：2020")
  message("有效中心数量：", nrow(city_centers))
  message("主中心人口：",   format(main_pop, big.mark = ","))
  message("次中心总人口：", format(sub_pop,  big.mark = ","))
  message("多中心指数：",   round(poly_index, 4))
}

四、批量并行计算（全国所有城市 × 所有年份）

4.1 封装核心计算函数

此处代码需下载讲义材料查看～

4.2 加载全部城市与年份文件

# 加载全部城市
# 这里仅选择三个城市为例进行计算
city_all <- read_sf(PATH$city_shp) %>%
  st_make_valid() %>%
  filter(!st_is_empty(.)) %>%
  st_transform(mycrs) %>%
  filter(市代码 %in% c(130100, 110000, 120000)) %>%
  select(city_code = 市代码, city_name = 市)

city_all

# 索引所有年份 tif 文件，以年份为名称
pop_files <- dir_ls(PATH$pop_tif, glob = "*.tif") %>%
  as.character() %>%
  set_names(str_extract(basename(.), "\\d{4}"))

pop_files

4.3 启动多线程并行计算

这部分就不要运行了，需要耗费大量的时间。

# 启动多线程（根据电脑核心数调整 workers）
plan(multisession, workers = 4)
message("使用核心数：", nbrOfWorkers())

# 创建结果保存目录
dir.create("res")

# 筛选尚未计算的城市（支持断点续算）
`%w/o%` <- function (x, y) {
  x[!x %in% y]
}
city_all$city_code %w/o% str_extract(fs::dir_ls("res"), "\\d{6}") -> newlist

# 对每个城市，并行遍历所有年份
walk(newlist, function(code) {
  message("正在计算城市代码：", code)

  result <- future_map_dfr(pop_files, ~calculate_polycentric(
    city_code = code,
    pop_file  = .x
  ))

  # 按城市代码保存结果 CSV
  if (nrow(result) > 0) {
    write_csv(result, path(str_c("res/", code, ".csv")))
  }
})

message("======== 全部计算完成！结果保存在 res/ 文件夹 ========")

说明：这里的外层循环（城市）用串行 walk，内层（年份）用 future_map_dfr 并行，兼顾了断点续算和资源利用效率。结果按城市代码逐个保存为 CSV，避免内存积累。

五、改进算法：预计算空间权重矩阵

5.1 为什么可以改进？

在上面的批量计算中，对同一个城市，每个年份都重新计算了一次空间权重矩阵。但实际上，空间权重矩阵只取决于城市边界的形状——它与年份无关，只要城市行政边界不变（本项目使用 2021 年固定边界），同一城市所有年份的空间权重矩阵完全相同。

因此，可以先把所有城市的权重矩阵计算并保存好，然后在计算各年数据时直接读取，显著减少重复计算。

对于一个有 N 个城市、T 个年份的数据集：

方法	权重矩阵计算次数
原始方法	N×T
改进方法	N（预计算一次）

当 T=25（2000～2024 年）时，改进方法可将权重矩阵的计算量缩减为原来的 1/25。

5.2 预计算并保存所有城市的权重矩阵

# 创建权重矩阵保存目录
dir.create("lwres")

city_all <- read_sf(PATH$city_shp) %>%
  st_make_valid() %>%
  filter(!st_is_empty(.)) %>%
  st_transform(mycrs) %>%
  filter(市代码 %in% c(130100, 110000, 120000))

# 计算单个城市权重矩阵的函数（以任意一个年份的人口文件为模板）
compute_lw <- function(city_code_x) {
  city_temp <- read_sf(PATH$city_shp) %>%
    st_make_valid() %>%
    st_transform(mycrs) %>%
    filter(市代码 == city_code_x)

  # 使用任意一年（如 2020 年）的栅格确定格点位置
  pop_demo <- rast(path(PATH$pop_tif, "2020.tif"))

  r_clip <- pop_demo %>%
    crop(vect(city_temp)) %>%
    mask(vect(city_temp))

  pop_points <- as.points(r_clip) %>%
    st_as_sf() %>%
    rename(pop = 1) %>%
    drop_na(pop)

  coords <- st_coordinates(pop_points)
  nb     <- dnearneigh(coords, d1 = 0, d2 = ANALYSIS_PARAMS$dist_nb)
  lw     <- nb2listw(nb, style = "W", zero.policy = TRUE)

  # 保存权重矩阵为 rds 文件
  lw %>% write_rds(paste0("lwres/", city_code_x, ".rds"))
}

# 并行计算所有城市的权重矩阵
plan(multisession, workers = 15)
future_imap(city_all$市代码, ~compute_lw(.x))

message("======== 所有城市权重矩阵计算完成，保存在 lwres/ ========")

5.3 使用预计算权重矩阵的改进版计算函数

此处代码需下载讲义材料查看～

5.4 使用改进算法批量计算全部数据

此处代码需下载讲义材料查看～

5.5 合并所有结果并保存

library(parallel)

# 使用多线程并行读取所有 CSV
makeCluster(16) -> cl

fs::dir_ls("resa") %>%
  parLapply(cl, ., read_csv) %>%
  bind_rows() -> df1

stopCluster(cl)

# 预览结果
df1

# 保存为 dta 格式（附加变量标签）
df1 %>%
  haven::write_dta(
    "sample_2000～2024年各城市多中心指标（王峤版本）.dta",
    label = "数据处理：微信公众号 RStata"
  )

message("最终数据已保存：2000～2024年各城市多中心指标（王峤版本）.dta")

最终结果数据集包含如下变量：

变量名	说明
`year`	年份（2000～2024）
`city_name`	城市名称
`city_code`	行政区划代码（2021 年版）
`center_count`	有效城市中心数量
`main_pop`	主中心人口数量（人）
`sub_pop`	次中心总人口数量（人）
`polycentric_index`	多中心指数 = 次中心人口 / (主中心人口 + 次中心人口)

关联课程与数据

名师讲堂｜使用 R 语言测算中国各城市多中心度与集聚程度：https://rstata.duanshu.com/#/course/4db63085cd024c3eb7a1746af2068498
2000～2024 年各城市多中心程度和人口集聚程度面板数据（Landscan 来源）：https://tidyfriday.cn/rsdb2/
2000～2024 年各省市区县、乡镇人口密度与人口数量面板数据（Landscan 来源）：https://tidyfriday.cn/rsdb2/
栅格数据的裁剪与面积汇总——基于 R 语言的方法：https://rstata.duanshu.com
如何将栅格数据处理成面板数据或时间序列数据（含并行计算内容）：https://rstata.duanshu.com

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 R 语言基于 LandScan 数据测算城市多中心指标

一、指标来源与计算原理

1.1 数据来源

1.2 指标定义

1.3 计算流程

1.4 参数设定说明

二、环境配置与数据准备

2.1 加载 R 包

2.2 全局路径与参数配置

三、单城市演示计算（北京市，2020 年）

3.1 读取城市边界与人口栅格

3.2 裁剪与栅格转点

3.3 构建空间权重矩阵

3.4 局部 Moran’s I 检验

3.5 空间类型分类（LISA 象限）

3.6 提取 HH 格点并聚类成中心

3.7 筛选有效中心并计算多中心指数

四、批量并行计算（全国所有城市 × 所有年份）

4.1 封装核心计算函数

此处代码需下载讲义材料查看～

4.2 加载全部城市与年份文件

4.3 启动多线程并行计算

五、改进算法：预计算空间权重矩阵

5.1 为什么可以改进？

5.2 预计算并保存所有城市的权重矩阵

5.3 使用预计算权重矩阵的改进版计算函数

此处代码需下载讲义材料查看～

5.4 使用改进算法批量计算全部数据

此处代码需下载讲义材料查看～

5.5 合并所有结果并保存

关联课程与数据

相关文章

评论