该课程已经讲解完，感兴趣的小伙伴可以点击文末的阅读原文跳转到 RStata 平台观看～
今天给大家分享使用 R 语言测算中国各城市多中心度与集聚程度的方法。
该指标在附件中的文献「多中心空间结构对中国城市创新绩效的影响_乔艺波」及「中国城市空间结构演化特征和和路径_袁超君」有提及：
接下来我们就按照论文介绍的过程进行计算。
首先加载所需的 R 包：
library(tidyverse)
library(sf)
library(spdep)
library(raster)
其中 spdep 用于计算邻接矩阵、空间滞后和局部莫兰指数。
读取市级行政区划数据：
read_sf("2021行政区划/市.shp") %>%
dplyr::select(-contains("类型")) %>%
st_make_valid() -> citymap
以广州市为例：
citymap %>%
filter(市 == "广州市") -> citytemp
使用 raster 包读取栅格数据：
raster::raster("pop-tif/2000.tif") -> rst
虽然现在比较推崇使用 terra 包进行栅格数据计算，不过 spdep 包不支持 terra 包的对象，所以这里还是使用 raster 包。
裁剪出北京市的：
rst %>%
raster::crop(citytemp) %>%
raster::mask(citytemp) -> rsttemp
更改成一个统一的名字：
names(rsttemp) <- "popden"
转换成 sf 包的 polygon 对象：
rsttemp %>%
raster::rasterToPolygons() %>%
st_as_sf() -> rstpoly
rstpoly
#> Simple feature collection with 9011 features and 1 field
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 112.95 ymin: 22.6 xmax: 114.0583 ymax: 23.93333
#> Geodetic CRS: +proj=longlat +datum=WGS84 +no_defs
#> First 10 features:
#> popden geometry
#> 1 25 POLYGON ((113.8583 23.93333...
#> 2 19 POLYGON ((113.8667 23.93333...
#> 3 24 POLYGON ((113.8833 23.93333...
#> 4 22 POLYGON ((113.8917 23.93333...
#> 5 22 POLYGON ((113.9167 23.93333...
#> 6 32 POLYGON ((113.925 23.93333,...
#> 7 40 POLYGON ((113.9417 23.93333...
#> 8 36 POLYGON ((113.95 23.93333, ...
#> 9 21 POLYGON ((113.9583 23.93333...
#> 10 21 POLYGON ((113.9667 23.93333...
计算各个多边形的相邻关系：
poly2nb(rstpoly, row.names = row.names(rstpoly), queen = TRUE) -> rstnb
rstnb
#> Neighbour list object:
#> Number of regions: 9011
#> Number of nonzero links: 66560
#> Percentage nonzero weights: 0.08197234
#> Average number of links: 7.386528
#> 3 regions with no links:
#> 657, 1753, 5027
#> 16 disjoint connected subgraphs
这里的 queen = TRUE 表示使用 queen 规则。
QUEEN 规则（共点或共边即视为邻接）；如果设置为 FALSE，则使用 ROOK 规则（仅共边视为邻接）。
将 nb 对象转换为空间权重矩阵：
w_mat <- nb2listw(rstnb, style = "W", zero.policy = TRUE)
w_mat
#> Characteristics of weights list object:
#> Neighbour list object:
#> Number of regions: 9011
#> Number of nonzero links: 66560
#> Percentage nonzero weights: 0.08197234
#> Average number of links: 7.386528
#> 3 regions with no links:
#> 657, 1753, 5027
#> 16 disjoint connected subgraphs
#>
#> Weights style: W
#> Weights constants summary:
#> n nn S0 S1 S2
#> W 9008 81144064 9008 2519.087 36210.48
计算局部莫兰指数：
localmoran(rstpoly$popden, w_mat) %>%
as_tibble() -> local_moran
local_moran
#> # A tibble: 9,011 × 5
#> Ii E.Ii Var.Ii Z.Ii `Pr(z != E(Ii))`
#> <localmrn> <localmrn> <localmrn> <localmrn> <localmrn>
#> 1 0.04770805 -5.298516e-06 0.023869687 0.3088279 0.7574524
#> 2 0.04875499 -5.382952e-06 0.024250069 0.3131198 0.7541896
#> 3 0.04802408 -5.312542e-06 0.011963781 0.4391096 0.6605822
#> 4 0.04759623 -5.340651e-06 0.012027080 0.4340514 0.6642511
#> 5 0.04748529 -5.340651e-06 0.012027080 0.4330398 0.6649859
#> 6 0.04710991 -5.200850e-06 0.011712253 0.4353515 0.6633073
#> 7 0.03292813 -5.090345e-06 0.011463398 0.3075937 0.7583915
#> 8 0.03598082 -5.145449e-06 0.009268964 0.3737813 0.7085670
#> 9 0.04775099 -5.354733e-06 0.009645963 0.4862485 0.6267910
#> 10 0.04757324 -5.354733e-06 0.009645963 0.4844387 0.6280746
#> # ℹ 9,001 more rows
莫兰指数的取值范围一般在 -1 到 1 之间。
当莫兰指数为正值时，表示数据呈正相关性，即相似的数据点（高值或低值）在空间上趋于聚集。 当莫兰指数为负值时，表示数据呈负相关性，即相似的数据点在空间上趋于分散。 莫兰指数接近 0 时，表示数据之间不存在明显的空间相关性，即数据在空间上随机分布。
计算空间滞后：
rstpoly$lpopden <- lag.listw(w_mat, rstpoly$popden)
合并 rstpoly 和 局部莫兰指数计算结果：
bind_cols(rstpoly, local_moran) -> rstpoly
根据局部莫兰指数及其显著性程度，挑选出市域范围内所有人口密度和人口密度的空间滞后值均高于该城市的相应平均值，且在 0.1 的显著性水平上显著的栅格。
rstpoly %>%
mutate(mean = mean(popden, na.rm = T),
lmean = mean(lpopden, na.rm = T)) %>%
filter(popden >= mean & lpopden >= lmean &
`Pr(z != E(Ii))` <= 0.1) -> rstpoly_main
把邻近的融合：
rstpoly_main %>%
st_union() -> rstpoly_main
拆分 polygons：
rstpoly_main %>%
st_cast("POLYGON") -> rstpoly_main
# plot(rstpoly_main)
rstpoly_main %>%
st_sf() %>%
mutate(id = row_number()) %>%
dplyr::select(id, everything()) %>%
st_make_valid() -> rmdf
rmdf
#> Simple feature collection with 19 features and 1 field
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 113.1917 ymin: 22.925 xmax: 113.8417 ymax: 23.55833
#> Geodetic CRS: +proj=longlat +datum=WGS84 +no_defs
#> First 10 features:
#> id geometry
#> 1 1 POLYGON ((113.4917 23.09167...
#> 2 2 POLYGON ((113.4417 22.98333...
#> 3 3 POLYGON ((113.3833 22.975, ...
#> 4 4 POLYGON ((113.4 22.95833, 1...
#> 5 5 POLYGON ((113.375 22.975, 1...
#> 6 6 POLYGON ((113.4 22.94167, 1...
#> 7 7 POLYGON ((113.3833 22.94167...
#> 8 8 POLYGON ((113.3667 22.95, 1...
#> 9 9 POLYGON ((113.5 23.08333, 1...
#> 10 10 POLYGON ((113.4 23.14167, 1...
下面我们需要找到距离在 2km 范围内的区域再合并起来。为此首先计算各个 polygon 间的距离、过滤出距离小于 2km 的：
st_distance(rmdf, rmdf) %>%
as_tibble() %>%
mutate_all(as.numeric) %>%
mutate(id1 = row_number()) %>%
dplyr::select(id1, everything()) %>%
gather(contains("V"), key = "id2", value = "dist") %>%
mutate(id2 = str_remove(id2, "V"),
id2 = as.numeric(id2)) %>%
filter(dist <= 2000 & id1 != id2 & id2 > id1) %>%
dplyr::select(-dist) -> idpair
idpair
#> # A tibble: 22 × 2
#> id1 id2
#> <int> <dbl>
#> 1 3 4
#> 2 3 5
#> 3 4 5
#> 4 4 6
#> 5 5 6
#> 6 4 7
#> 7 5 7
#> 8 6 7
#> 9 4 8
#> 10 5 8
#> # ℹ 12 more rows
这里我们就得到了距离在 2km 内的区域对，下面我们借助网络分析的方法把这些合并起来：
library(igraph)
g <- graph_from_data_frame(idpair, directed = FALSE)
# plot(g)
我们需要把相连接的这些区域合并起来。
找到图的连通分量：
components <- components(g)
components$membership -> membership
tibble(
id = names(membership),
group = membership
) %>%
mutate(id = as.numeric(id)) -> groupdf
groupdf
#> # A tibble: 15 × 2
#> id group
#> <dbl> <dbl>
#> 1 3 1
#> 2 4 1
#> 3 5 1
#> 4 6 1
#> 5 7 1
#> 6 1 2
#> 7 10 2
#> 8 8 1
#> 9 11 2
#> 10 12 2
#> 11 16 3
#> 12 9 2
#> 13 15 1
#> 14 17 2
#> 15 18 3
然后就可以分组聚合了：
rmdf %>%
left_join(groupdf) %>%
mutate(groupmax = max(group, na.rm = T)) %>%
mutate(group = if_else(is.na(group), id + groupmax, group)) %>%
group_by(group) %>%
summarise() %>%
mutate(group = row_number()) -> rmdf2
rmdf2
#> Simple feature collection with 7 features and 1 field
#> Geometry type: GEOMETRY
#> Dimension: XY
#> Bounding box: xmin: 113.1917 ymin: 22.925 xmax: 113.8417 ymax: 23.55833
#> Geodetic CRS: +proj=longlat +datum=WGS84 +no_defs
#> # A tibble: 7 × 2
#> group geometry
#> * <int> <GEOMETRY [°]>
#> 1 1 MULTIPOLYGON (((113.375 22.975, 113.3667 22.975, 113.3667 22.96667, 113…
#> 2 2 MULTIPOLYGON (((113.2417 23.05833, 113.25 23.05833, 113.2583 23.05833, …
#> 3 3 MULTIPOLYGON (((113.8333 23.3, 113.8417 23.3, 113.8417 23.30833, 113.83…
#> 4 4 POLYGON ((113.4417 22.98333, 113.45 22.98333, 113.45 22.99167, 113.4417…
#> 5 5 POLYGON ((113.6083 23.125, 113.6 23.125, 113.6 23.11667, 113.5917 23.11…
#> 6 6 POLYGON ((113.5917 23.54167, 113.5917 23.55, 113.5917 23.55833, 113.583…
#> 7 7 POLYGON ((113.2167 23.375, 113.2167 23.38333, 113.2083 23.38333, 113.20…
再检查有没有 > 2km 的：
as.numeric(st_distance(rmdf2, rmdf2)) -> distmat
sum(distmat < 2000 & distmat != 0)
#> [1] 0
看来是没有了。
依据面积大于 1 km2 和总人口大于 5 万人的标准，挑选出最终的城市区域，即市域的主中心及次中心。
计算每个区域的总人口：
raster::extract(raster::area(rsttemp) * rsttemp, rmdf2, fun = sum, na.rm = T) %>%
as_tibble() -> popres
popres
#> # A tibble: 7 × 1
#> V1
#> <dbl>
#> 1 125111.
#> 2 3025257.
#> 3 59911.
#> 4 763.
#> 5 6413.
#> 6 39301.
#> 7 86020.
计算面积：
rmdf2 %>%
st_area() %>%
units::set_units(km2) %>%
as.numeric() -> area
rmdf2 %>%
mutate(pop = popres$V1 / 10000, area_km2 = area) %>%
filter(pop >= 5 & area_km2 >= 1) -> rmdf3
# plot(rmdf3[1])
其中主中心为：
rmdf3 %>%
slice_max(n = 1, order_by = pop) -> maindf
n 为每个城市市域范围内的中心个数 (包括主中心及次中心)：
n <- nrow(rmdf3)
IMP：为每个中心的重要性 = POP x AREA x DIS：
rmdf3 %>%
anti_join(maindf %>% st_drop_geometry() %>% dplyr::select(group)) -> subdf
subdf %>%
mutate(DIS_km = st_distance(., maindf)[,1] %>%
units::set_units("km") %>%
as.numeric()) %>%
mutate(IMP = pop * area_km2 * DIS_km) -> rmdfa
对于主中心而言，其重要性为主中心人口、所有中心中的最大面积（通常即主中心本身的面积）和距主中心最远的次中心与主中心之间的距离三者之乘积：
maindf %>%
mutate(DIS_km = st_distance(maindf, subdf) %>%
max(na.rm = T) %>%
units::set_units("km") %>%
as.numeric()) %>%
mutate(IMP = pop * area_km2 * DIS_km) -> rmdfb
# 合并两部分结果
bind_rows(rmdfa, rmdfb) -> rmdf4
σ_Cobs 为市域各中心重要性的标准差：
rmdf4 %>%
mutate(sigma_Cobs = sd(IMP, na.rm = T)) %>%
mutate(sigma_Cmax = sd(c(max(pop, na.rm = T) * max(area_km2, na.rm = T) * max(DIS_km, na.rm = T), 0), na.rm = T)) %>%
mutate(polycentricity = 1 - (sigma_Cobs / sigma_Cmax)) %>%
st_drop_geometry() -> rmdf5
rmdf5
#> # A tibble: 4 × 8
#> group pop area_km2 DIS_km IMP sigma_Cobs sigma_Cmax polycentricity
#> * <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 12.5 19.8 11.9 2945. 1073098. 1518783. 0.293
#> 2 3 5.99 3.94 38.8 916. 1073098. 1518783. 0.293
#> 3 7 8.60 7.88 17.8 1206. 1073098. 1518783. 0.293
#> 4 2 303. 183. 38.8 2147883. 1073098. 1518783. 0.293
主城区人口在中心区域的集聚程度为主城区中心 区域人口与主城区总人口之比：
maindf$pop / sum(rmdf5$pop, na.rm = T)
#> [1] 0.9177738
最后我们就可以把上面的代码封装成一个函数了：get_polycentricity.R。使用：
source("get_polycentricity.R", local = T)
get_polycentricity("北京市", "pop-tif/2000.tif")
#> [1] "pop-tif/2000.tif"
#> # A tibble: 1 × 5
#> city rst polycentricity N_center pop_cluster
#> <chr> <chr> <dbl> <int> <dbl>
#> 1 北京市 pop-tif/2000.tif 0.294 4 0.905
由于该函数里面没有包的加载和矢量数据读取的代码，所以需要确保加载了相关的 R 包和准备好相关的矢量数据才可以使用。
再例如计算广州市所有年份的：
"广州市" %>%
tidyr::crossing(as.character(fs::dir_ls("pop-tif"))) %>%
set_names("cityname", "rstpath") %>%
mutate(res = map2(cityname, rstpath, get_polycentricity)) -> dfres
dfres %>%
unnest(res)
#> # A tibble: 23 × 7
#> cityname rstpath city rst polycentricity N_center pop_cluster
#> <chr> <chr> <chr> <chr> <dbl> <int> <dbl>
#> 1 广州市 pop-tif/2000.tif 广州市 pop-tif… 0.293 4 0.918
#> 2 广州市 pop-tif/2001.tif 广州市 pop-tif… 0.293 4 0.933
#> 3 广州市 pop-tif/2002.tif 广州市 pop-tif… 0.293 4 0.938
#> 4 广州市 pop-tif/2003.tif 广州市 pop-tif… 0.293 4 0.952
#> 5 广州市 pop-tif/2004.tif 广州市 pop-tif… 0.293 4 0.951
#> 6 广州市 pop-tif/2005.tif 广州市 pop-tif… 0.293 4 0.951
#> 7 广州市 pop-tif/2006.tif 广州市 pop-tif… 0.293 4 0.937
#> 8 广州市 pop-tif/2007.tif 广州市 pop-tif… 0.368 5 0.921
#> 9 广州市 pop-tif/2008.tif 广州市 pop-tif… 0.293 4 0.930
#> 10 广州市 pop-tif/2009.tif 广州市 pop-tif… 0.293 4 0.930
#> # ℹ 13 more rows
绘图展示：
dfres %>%
unnest(res) %>%
mutate(year = str_extract(rstpath, "\\d{4}")) %>%
type_convert() %>%
ggplot(aes(year, polycentricity)) +
geom_line() +
geom_point() -> p
# ggsave("unnamed-chunk-33-1.png", device = png, width = 10, height = 6)
这样我们就得到了广州历年的计算结果。
直播信息
为了让大家更好的理解本文内容，欢迎各位名师讲堂会员参加明晚 8 点的直播课：「使用 R 语言测算中国各城市多中心度与集聚程度」。
直播地址：腾讯会议(需要报名 RStata 培训班参加) 讲义材料：需要购买 RStata 名师讲堂会员，详情可阅读：一起来学习 R 语言和 Stata 啦！学习过程中遇到的问题也可以随时提问！
更多关于 RStata 会员的更多信息可添加微信号 r_stata 咨询：
附件下载（点击文末的阅读原文即可跳转）：
https://rstata.duanshu.com/#/brief/course/4db63085cd024c3eb7a1746af2068498
