名师讲堂｜使用 R 语言基于LandScan数据测算城市多中心指标（商玉萍版本）

今天给大家介绍如何使用 R 语言基于 LandScan 全球人口栅格数据，按照商玉萍（2022）论文中的方法，同时测算城市多中心结构的 5 大指标，包括中心数量、帕累托指数、多中心指数（含距离）以及去中心化指标。

一、指标来源与计算原理

1.1 数据来源

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

1.2 文献来源

5 大指标均来自以下论文：

商玉萍. 中国城市多中心空间战略的创新绩效研究——基于集聚经济与舒适度的视角 [J]. 经济学（季刊）, 2022.

该论文从集聚经济与舒适度双重视角，考察城市多中心空间战略对创新绩效的影响，所使用的多中心测量指标体系是目前文献中最为系统的之一。

1.3 五大指标定义

变量名	中文名	含义
center	城市中心数量	识别出的有效城市中心（高密度聚集区）数量
pareto	帕累托指数	各中心人口规模的秩-规模幂律回归系数（绝对值），越小表示多中心越均衡
poly	多中心指数	纳入距离的人口规模标准差指数，越小表示分布越均衡
sub3	去中心化指标（3 km）	CBD 3 千米以外的人口占城市总人口的比例
sub5	去中心化指标（5 km）	CBD 5 千米以外的人口占城市总人口的比例

1.4 各指标计算方法

1.5 计算流程总览

整个计算分为以下步骤：

1. 坐标系转换：将 LandScan 栅格数据投影到等面积坐标系（Albers 投影），确保距离计算准确；
2. 裁剪与掩膜：按城市行政边界裁剪栅格，获取该城市范围内的人口格点；
3. 栅格转点：将栅格像元转为空间点数据；
4. 构建点邻接：使用 1000 m 距离阈值构建邻接关系，生成空间权重矩阵；
5. 局部 Moran’s I 检验：计算每个格点的局部莫兰指数，识别显著的 HH 聚集区格点；
6. 聚类成中心：将 HH 格点按空间邻接分组，形成连续的”高密度区块”；
7. 筛选有效中心：要求每个区块格点数 ≥ 3、总人口 ≥ 100,000；
8. 计算 5 大指标：center、pareto、poly、sub3、sub5。

1.6 参数设定说明

ANALYSIS_PARAMS <- list(
  sig_level = 0.05,   # 局部莫兰指数显著性阈值
  min_cells = 3,      # 中心最小格点数（过滤噪点）
  min_pop   = 100000, # 中心最小人口（单位：人），论文基准为 10 万
  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
  units      # 单位转换（用于 st_distance 去单位）
)

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  = "pop-tif2",           # 已重投影到 mycrs 的 tif 文件夹
  city_shp = "2021行政区划/市.shp"  # 2021 年市级行政区划
)

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

数据预处理提示：如果你的 LandScan tif 文件还是 WGS84 坐标系，需要先转换到 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_all <- read_sf(PATH$city_shp) %>%
  st_make_valid() %>%
  filter(!st_is_empty(.)) %>%
  st_transform(mycrs) %>%
  filter(市代码 %in% c(130100, 110000, 120000))

# 筛选北京市
city_demo <- city_all %>% filter(市代码 == 110000)
city_demo

# 读取 2020 年人口栅格
pop_file <- "pop-tif2/2020.tif"
r_demo <- rast(pop_file)
r_demo

3.2 裁剪与栅格转点

# 按北京边界裁剪并掩膜
r_clip <- r_demo %>%
  crop(vect(city_demo)) %>%
  mask(vect(city_demo))

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

pts

3.3 构建空间权重矩阵

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

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

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

参数说明

• style = “W”：行标准化权重。若一个格点有 4 个邻居，则每个邻居权重 = 1/4；有 2 个邻居则权重 = 1/2。
• zero.policy = TRUE：允许无邻居的孤立点（如城市边缘格点）存在。

3.4 局部 Moran’s I 检验

局部 Moran’s I 的核心作用：找出人口高密度且周围也是高密度的区域（HH 聚集区），这正是城市中心的候选位置。

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

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

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

3.5 LISA 空间类型分类

pts <- pts %>%
  mutate(
    p_val   = lmo[, 5],
    sig     = p_val < ANALYSIS_PARAMS$sig_level,
    med     = median(pop),
    lag_val = lag.listw(lw, pop),            # 空间滞后：邻居的加权平均人口
    self    = if_else(pop > med, "H", "L"),
    nbr     = if_else(lag_val > med, "H", "L"),
    type    = case_when(
      sig & self == "H" & nbr == "H" ~ "HH",
      TRUE ~ "other"
    )
  )
pts

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

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

# 筛选 HH 格点
hh <- filter(pts, type == "HH")
hh

# 对 HH 格点再次做空间邻接，将连续区域标为同一 cluster
coords_hh  <- st_coordinates(hh)
nb_hh      <- dnearneigh(coords_hh, 0, ANALYSIS_PARAMS$dist_nb)
hh$cluster <- n.comp.nb(nb_hh)$comp.id

# 按 cluster 汇总，筛选有效中心（人口≥10万，格点数≥3）
centers <- hh %>%
  group_by(cluster) %>%
  summarise(
    pop = sum(pop, na.rm = TRUE),
    n   = n(),
    x   = mean(st_coordinates(geometry)[, 1]),
    y   = mean(st_coordinates(geometry)[, 2]),
    .groups = "drop"
  ) %>%
  filter(n >= ANALYSIS_PARAMS$min_cells, pop >= ANALYSIS_PARAMS$min_pop)

centers

3.7 计算 5 大指标

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

---

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

4.1 封装核心计算函数

将以上步骤封装为函数，接收城市代码和人口 tif 文件路径，返回该城市该年份的 5 大指标：

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

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 = 市)

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

pop_files

4.3 启动多线程并行计算

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

plan(multisession, workers = 4)
message("并行核心数：", nbrOfWorkers())

# 创建结果目录
dir_create("res2")

# 所有城市并行计算 → 每个城市单独保存到 res2/
future_walk(city_all$city_code, function(code) {
  tryCatch({
    # 串行计算该城市所有年份
    res <- map_dfr(pop_files, ~ calc_paper_indicators(code, .x))
    if (nrow(res) > 0) {
      write_csv(res, str_glue("res2/{code}.csv"))
    }
    message("完成城市：", code)
  }, error = function(e) {
    message("失败城市：", code)
  })
})

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

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

---

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

5.1 为什么可以改进？

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

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

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

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

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

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

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

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

改进版函数从外部接收 lw 参数，不在函数内部重新计算权重矩阵：

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

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

dir.create("resb")

# 对每个城市：读取预计算好的 lw，循环所有年份
calc_paper_indicators3 <- function(citycode) {
  city_demo <- city_all %>% filter(市代码 == citycode)
  lw        <- read_rds(paste0("lwres/", citycode, ".rds"))

  for (f in fs::dir_ls("pop-tif2")) {
    res <- calc_paper_indicators2(city_demo, f, lw)
    if (nrow(res) > 0) {
      write_csv(
        res,
        paste0("resb/", citycode, "_", str_extract(f, "\\d{4}"), ".csv")
      )
    }
  }
}

# 并行计算所有城市（城市间并行，城市内串行遍历年份）
plan(multisession, workers = 15)
future_imap(city_all$市代码, ~calc_paper_indicators3(.x))

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

5.5 合并所有结果并保存

library(parallel)

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

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

stopCluster(cl)

# 预览结果
df1

# 保存为 dta 格式
df1 %>%
  haven::write_dta(
    "2000～2024年各城市多中心指标（商玉萍版本）.dta",
    label = "数据处理：微信公众号 RStata"
  )

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

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

变量名	说明
year	年份（2000～2024）
city_code	行政区划代码（2021 年版）
city_name	城市名称
center	有效城市中心数量
pareto	帕累托指数，越小越均衡
poly	含距离的多中心指数，越小越均衡
sub3	CBD 3 千米以外的人口占比
sub5	CBD 5 千米以外的人口占比
total_pop	城市总人口

通常，pareto、poly 值越小，代表各中心的人口分布越均衡（均等）；sub3、sub5 越大，代表城市人口越去中心化。

---

六、稳健性检验：放宽中心人口门槛至 1 万人

商玉萍（2022）论文中提供了一项稳健性检验：将”总人口在 10 万人以上”的条件改为”总人口在 1 万人以上“，重新确定每个城市的中心数量（2center）作为稳健性指标。

这一操作的实现方式非常简单，只需将 ANALYSIS_PARAMS$min_pop 改为 10000，其余代码完全不变：

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

---

关联课程与数据

• 名师讲堂｜使用 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数据测算城市多中心指标（商玉萍版本）