该课程已经讲解完，感兴趣的小伙伴可以文末的阅读原文跳转到 RStata 平台观看录播～
方颖、白秀叶（2022）一文中使用了各城市空间形态的紧凑度的指标，该指标的计算原理是这样的：
这一操作可以使用 R 语言完全实现。步骤如下：
使用市辖区的矢量数据裁剪土地利用栅格数据得到每个城市的城镇建成区数据； 提取面积最大的； 在该区域内随机抽点，计算任意两点间的平均距离，记为 dist1； 构造与该区域面积相当的圆，同样计算圆内任意两点间的平均距离，记为 dist2； dist1 / dist2 即为城市空间形态的紧凑度指标。
裁剪各城市的城镇建成区
附件中给大家准备了 2021 年各级行政区划矢量数据以及中国土地利用数据（1980-2015），论文中使用的是 2010 年的，下面我们会循环计算所有年份的。
首先加载相关的 R 包：
library(tidyverse)
library(sf)
library(geosphere)
library(raster)
由于我们需要使用 raster::rasterToPolygons() 函数，所以这次就使用 raster 包，而不是 terra 包。
读取 2015 年的土地利用数据：
raster("中国土地利用数据1980-2015/2015年中国土地利用现状遥感监测数据/lucc2015/w001001.adf") -> rst
rst
#> class : RasterLayer
#> dimensions : 4083, 5116, 20888628 (nrow, ncol, ncell)
#> resolution : 1000, 1000 (x, y)
#> extent : -2745867, 2370133, 1871283, 5954283 (xmin, xmax, ymin, ymax)
#> crs : +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
#> source : w001001.adf
#> names : w001001
#> values : 11, 99 (min, max)
#> attributes :
#> ID COUNT
#> from: 11 464950
#> to : 99 1037
plot(rst)
该数据的坐标参考系是：
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"
因此裁剪使用的矢量数据也需要是该坐标系的。城市市辖区的数据可以通过下面的方法构造：
read_sf("2021行政区划/县.shp") %>%
filter(县类型 == "市辖区") %>%
dplyr::select(-contains("类型")) %>%
group_by(省, 省代码, 市, 市代码) %>%
summarise() %>%
ungroup() %>%
st_transform(mycrs) %>%
st_make_valid()-> city
从土地利用数据的介绍中可以看到，城镇区域的编号是 51，所以我们仅仅保留编号是 51 的栅格：
# 提取城镇建成区
rst[rst != 51] <- NA
rst %>%
raster::crop(city[1,]) %>%
raster::mask(city[1,]) %>%
raster::rasterToPolygons() -> rstpoly
rstpoly
#> class : SpatialPolygonsDataFrame
#> features : 1011
#> extent : 1503133, 1572133, 3381283, 3483283 (xmin, xmax, ymin, ymax)
#> crs : +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
#> variables : 1
#> names : w001001
#> min values : 51
#> max values : 51
# 转换成 sf 对象
rstpoly %>%
st_as_sf()
#> Simple feature collection with 1011 features and 1 field
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 1503133 ymin: 3381283 xmax: 1572133 ymax: 3483283
#> Projected CRS: +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
#> First 10 features:
#> w001001 geometry
#> 1 51 POLYGON ((1524133 3483283, ...
#> 2 51 POLYGON ((1526133 3483283, ...
#> 3 51 POLYGON ((1523133 3482283, ...
#> 4 51 POLYGON ((1524133 3482283, ...
#> 5 51 POLYGON ((1525133 3482283, ...
#> 6 51 POLYGON ((1529133 3482283, ...
#> 7 51 POLYGON ((1524133 3481283, ...
#> 8 51 POLYGON ((1525133 3481283, ...
#> 9 51 POLYGON ((1548133 3477283, ...
#> 10 51 POLYGON ((1547133 3476283, ...
然后我们可以使用 st_union() 函数把相邻近的 polygon “融合”成大的多边形，分别计算每个大多边形的面积，保留第一个（最大的）：
rstpoly %>%
st_as_sf() %>%
as_tibble() %>%
st_sf() %>%
st_union() %>%
st_cast("POLYGON") %>%
st_sf() %>%
mutate(area = st_area(.)) %>%
arrange(desc(area)) %>%
slice(1)
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 1523133 ymin: 3414283 xmax: 1557133 ymax: 3464283
#> Projected CRS: +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
#> area geometry
#> 1 759000000 [m^2] POLYGON ((1526133 3424283, ...
这样我们就构造好了第一个城市的，下面我们采用多线程的方法处理所有城市的：
library(parallel)
makeCluster(16) -> cl
clusterEvalQ(cl, ({
library(tidyverse)
library(sf)
library(geosphere)
library(raster)
})) -> tempres
clusterExport(cl, "city")
dir.create("cityres2")
parLapply(cl, 1:nrow(city), function(x){
lapply(fs::dir_ls("中国土地利用数据1980-2015",
regexp = "w001001.adf",
recurse = T),
function(y){
try({
str_extract(y, "lucc\\d{4}") -> newy
if(!file.exists(paste0("cityres2/", x, "_", newy, ".rds"))) {
raster(y) -> rst
rst[rst != 51] <- NA
rst %>%
raster::crop(city[x,]) %>%
raster::mask(city[x,]) %>%
raster::rasterToPolygons() -> rstpoly
if(nrow(rstpoly) >= 1) {
st_as_sf(rstpoly) %>%
as_tibble() %>%
st_sf() %>%
st_union() %>%
st_cast("POLYGON") %>%
st_sf() %>%
mutate(area = st_area(.)) %>%
arrange(desc(area)) %>%
slice(1) %>%
bind_cols(st_drop_geometry(city[x,])) %>%
mutate(file = y) %>%
write_rds(paste0("cityres2/", x, "_", newy, ".rds"))
}
}
})
})
}) -> tempres
由于土地利用数据中并不是每个城市都有“城镇”区域，所以上面的代码中增加了一个 if 判断。
然后就可以合并处理得到的 rds 文件了：
fs::dir_ls("cityres2") %>%
parLapply(cl, ., read_rds) %>%
bind_rows() -> maxcity2
maxcity2 %>%
mutate(file = str_extract(file, "lucc\\d{4}"),
file = str_remove_all(file, "lucc"),
file = as.numeric(file)) %>%
rename(year = file) %>%
as_tibble() %>%
st_sf() -> maxcity2
maxcity2 %>%
write_rds("maxcity2.rds")
maxcity2
#> Simple feature collection with 2035 features and 6 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: -1583867 ymin: 1883283 xmax: 2002578 ymax: 5628946
#> Projected CRS: +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
#> # A tibble: 2,035 × 7
#> area 省 省代码 市 市代码 year geometry
#> [m^2] <chr> <dbl> <chr> <dbl> <dbl> <POLYGON [m]>
#> 1 11000000 广东省 440000 梅州市 441400 1980 ((1127133 2600283, 1126133 260028…
#> 2 18000000 广东省 440000 梅州市 441400 1990 ((1125578 2599756, 1125578 260075…
#> 3 15000000 广东省 440000 梅州市 441400 1995 ((1127130 2599946, 1126130 259994…
#> 4 14000000 广东省 440000 梅州市 441400 2000 ((1127133 2600283, 1126133 260028…
#> 5 14000000 广东省 440000 梅州市 441400 2005 ((1127133 2600283, 1126133 260028…
#> 6 14000000 广东省 440000 梅州市 441400 2010 ((1126133 2600283, 1126133 260128…
#> 7 16000000 广东省 440000 梅州市 441400 2015 ((1126133 2601283, 1126133 260228…
#> 8 21000000 广东省 440000 汕头市 440500 1980 ((1196133 2504283, 1196133 250528…
#> 9 50000000 广东省 440000 汕头市 440500 1990 ((1195578 2504756, 1195578 250575…
#> 10 51000000 广东省 440000 汕头市 440500 1995 ((1196130 2501946, 1195130 250194…
#> # ℹ 2,025 more rows
平均区域内任意两点间的距离
计算点对的距离可以使用 geosphere 包：
library(geosphere)
编写一个计算区域内任意两点间平均距离的函数：
calculate_avg_distance <- function(province_polygon, n_points = 1000) {
random_points <- st_sample(province_polygon, size = n_points)
coords <- st_coordinates(random_points)
distances <- distm(coords)
# 计算平均距离（只取上三角矩阵，避免重复计算）
avg_distance <- mean(distances[upper.tri(distances)])
return(avg_distance)
}
然后就可以使用了：
clusterExport(cl, "calculate_avg_distance")
# calculate_avg_distance() 需要 4326 坐标系下的数据
maxcity2 %>%
st_transform(4326) -> maxcity2
maxcity2$avg_distance <- parallel::parSapply(cl, st_geometry(maxcity2), calculate_avg_distance)
# 显示结果
maxcity2 %>%
st_drop_geometry() %>%
dplyr::select(-contains("类型")) %>%
mutate(avg_distance = avg_distance / 1000) %>%
rename(avg_distance_km = avg_distance) -> maxcitydf2
maxcitydf2
#> # A tibble: 2,035 × 7
#> area 省 省代码 市 市代码 year avg_distance_km
#> [m^2] <chr> <dbl> <chr> <dbl> <dbl> <dbl>
#> 1 11000000 广东省 440000 梅州市 441400 1980 1.87
#> 2 18000000 广东省 440000 梅州市 441400 1990 2.29
#> 3 15000000 广东省 440000 梅州市 441400 1995 2.08
#> 4 14000000 广东省 440000 梅州市 441400 2000 2.13
#> 5 14000000 广东省 440000 梅州市 441400 2005 2.11
#> 6 14000000 广东省 440000 梅州市 441400 2010 2.10
#> 7 16000000 广东省 440000 梅州市 441400 2015 2.36
#> 8 21000000 广东省 440000 汕头市 440500 1980 3.11
#> 9 50000000 广东省 440000 汕头市 440500 1990 4.10
#> 10 51000000 广东省 440000 汕头市 440500 1995 4.07
#> # ℹ 2,025 more rows
这样我们就得到了最开始说的 dist1。
构建等面积圆
首先计算每个区域的面积：
maxcitydf2 %>%
left_join(st_transform(city, 4326)) %>%
st_sf() %>%
st_centroid() -> citycentroid
citycentroid
#> Simple feature collection with 2035 features and 7 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 85.19786 ymin: 18.36389 xmax: 131.3015 ymax: 50.21737
#> Geodetic CRS: WGS 84
#> # A tibble: 2,035 × 8
#> area 省 省代码 市 市代码 year avg_distance_km
#> * [m^2] <chr> <dbl> <chr> <dbl> <dbl> <dbl>
#> 1 11000000 广东省 440000 梅州市 441400 1980 1.87
#> 2 18000000 广东省 440000 梅州市 441400 1990 2.29
#> 3 15000000 广东省 440000 梅州市 441400 1995 2.08
#> 4 14000000 广东省 440000 梅州市 441400 2000 2.13
#> 5 14000000 广东省 440000 梅州市 441400 2005 2.11
#> 6 14000000 广东省 440000 梅州市 441400 2010 2.10
#> 7 16000000 广东省 440000 梅州市 441400 2015 2.36
#> 8 21000000 广东省 440000 汕头市 440500 1980 3.11
#> 9 50000000 广东省 440000 汕头市 440500 1990 4.10
#> 10 51000000 广东省 440000 汕头市 440500 1995 4.07
#> # ℹ 2,025 more rows
#> # ℹ 1 more variable: geometry <POINT [°]>
然后以每个多边形的质心为中心构造等面积圆、再计算圆内的平均距离：
dir.create("newcityres2")
clusterExport(cl, "citycentroid")
parLapply(cl, 1:nrow(citycentroid), function(x){
citycentroid %>%
slice(x) -> temppc
st_centroid(temppc) %>%
dplyr::select(-area) %>%
st_buffer(dist = units::set_units(sqrt(as.numeric(temppc$area) / pi), m)) -> equal_circle
calculate_avg_distance(equal_circle$geometry) -> ec_dist
citycentroid %>%
slice(x) %>%
st_drop_geometry() %>%
dplyr::select(-area) %>%
mutate(adk_ec = ec_dist / 1000) %>%
write_rds(paste0("newcityres2/", x, ".rds"))
}) -> tempres
合并 rds 文件：
fs::dir_ls("newcityres2") %>%
parLapply(cl, ., readr::read_rds) %>%
bind_rows() -> maxcitydf2a
这样就得到了 dist2 指标。
计算空间形态的紧凑度
使用上述的 dist1/dist2 即可：
maxcitydf2a %>%
mutate(spatial_compactness = avg_distance_km / adk_ec)
#> # A tibble: 2,035 × 8
#> 省 省代码 市 市代码 year avg_distance_km adk_ec spatial_compactness
#> <chr> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 广东省 440000 梅州市 441400 1980 1.87 1.71 1.09
#> 2 广东省 440000 汕头市 440500 1995 4.07 3.68 1.10
#> 3 广东省 440000 韶关市 440200 1990 3.56 2.60 1.37
#> 4 吉林省 220000 松原市 220700 2010 2.92 2.32 1.26
#> 5 吉林省 220000 松原市 220700 2015 2.99 2.36 1.27
#> 6 甘肃省 620000 白银市 620400 1980 2.61 2.26 1.16
#> 7 甘肃省 620000 白银市 620400 1990 2.74 2.19 1.26
#> 8 甘肃省 620000 白银市 620400 1995 3.35 2.52 1.33
#> 9 甘肃省 620000 白银市 620400 2000 2.87 2.35 1.22
#> 10 甘肃省 620000 白银市 620400 2005 2.99 2.59 1.15
#> # ℹ 2,025 more rows
这样就得到了各城市的空间形态紧凑度了，紧凑度越大，该城市的空间形态越“不紧凑”。感觉称为不紧凑度或许更合适。
绘图展示：
直播信息
该课程已经讲解完，感兴趣的小伙伴可以文末的阅读原文跳转到 RStata 平台观看录播～
直播地址：腾讯会议(需要报名 RStata 培训班参加) 讲义材料：需要购买 RStata 名师讲堂会员，详情可阅读：一起来学习 R 语言和 Stata 啦！学习过程中遇到的问题也可以随时提问！
更多关于 RStata 会员的更多信息可添加微信号 r_stata 咨询：
附件下载（点击文末的阅读原文即可跳转）：
https://rstata.duanshu.com/#/brief/course/5ebfe03c91774471b483104e9477e70f
