今天给大家分享城市共同边界及边界乡镇识别与距离计算的完整处理方法。该指标可用于研究城市交界地区的经济发展、人口流动、环境溢出等空间边界效应，是空间计量与城市经济学研究中的重要变量。
处理过程基于 2021 年行政区划矢量数据，识别出全国 955 对相邻城市，提取其共同边界，并计算出各边界乡镇到共同边界的距离。最终产出包括：城市共同边界矢量文件、边界相交乡镇列表及距离矩阵，可直接用于 Stata 或 R 的实证分析。
指标介绍
城市共同边界：指两个相邻城市行政区划多边形相交形成的线段（shared border）。例如成都市与德阳市相邻，两者的行政边界相交形成一条共同边界线。
边界相交乡镇：指行政边界多边形与城市共同边界线相交的乡镇（街道）。这些乡镇位于城市交界处，是研究”边界效应”（border effect）的核心样本。
距边界距离：指乡镇行政中心（质心）到最近的城市共同边界线的直线距离（米），用于构造连续的空间边界距离变量。
乡镇对-城市对距离矩阵：对于每一对相邻城市，分别列出两侧所有边界乡镇及其到共同边界的距离，形成”乡镇对”级别的分析数据（例如成都侧某乡镇与德阳侧某乡镇各为一行，配对形成乡镇对）。
下图展示了全国城市共同边界及边界相交乡镇的分布情况：
计算方法
整个处理流程分为九个步骤，完整代码来自 border_distance_touching.R，下面逐一介绍。
第 1 步：读取行政区划数据
首先加载所需 R 包，读取 2021 年城市级和乡镇级行政区划矢量数据，并统一投影坐标系（Albers 等面积投影）：
suppressPackageStartupMessages({
library(tidyverse)
library(sf)
library(haven)
})
if (!dir.exists("res")) dir.create("res")
MAX_TOWNS_PER_SIDE <- 500
GRID_SIZE <- 400000
message("[1/9] 读取数据...")
city_raw <- st_read("2021行政区划/市.shp", quiet = TRUE)
town_raw <- readRDS("town.rds")
proj_crs <- st_crs(town_raw)
message(sprintf(" CRS: %s", proj_crs$input))
message(sprintf(" 城市数：%d，乡镇数：%d", nrow(city_raw), nrow(town_raw)))
city <- city_raw %>%
select(市代码, 市名称 = 市, 省名称 = 省) %>%
st_transform(crs = proj_crs)
|
关键说明：
town.rds 为预处理后的乡镇级矢量数据（含 43,366 个乡镇多边形），已使用 Albers 等面积投影；- 城市数据同步转换到相同 CRS，确保后续空间运算准确无误；
MAX_TOWNS_PER_SIDE <- 500 控制每侧最多保留的乡镇数量，防止数据量过大；GRID_SIZE <- 400000 为格网大小（400 km），用于后续空间索引加速。
第 2 步：预计算乡镇质心坐标
为了提高后续大规模距离计算的速度，先一次性计算所有乡镇的质心坐标并存入数据框（避免反复调用 st_centroid()）：
message("\n[2/9] 预计算乡镇质心坐标...")
town_cent <- st_centroid(town_raw)
town_coords <- st_coordinates(town_cent)
town_base <- town_raw %>%
st_drop_geometry() %>%
mutate(cx = town_coords[, 1], cy = town_coords[, 2])
|
关键说明：
st_centroid() 计算每个乡镇多边形的质心（几何中心）；st_coordinates() 提取质心的 x、y 坐标；- 将几何信息转为普通数据框列（
cx、cy），为后续格网索引和空间匹配做准备。
第 3 步：判断乡镇所属城市
乡镇矢量数据本身不包含”所属城市”字段，需要用空间分析方法进行匹配。方法为：先将质心转回 sf 对象，再用 st_within() 判断每个质心落在哪个城市多边形内部。
# ── 3. 用质心判断乡镇所属城市
# 此处代码需下载讲义材料查看～
|
少数乡镇质心可能因精度问题落在所有城市之外（返回 NA），此时用”最近城市”法则补充：
unmatched <- which(is.na(matched_city_idx))
if (length(unmatched) > 0) {
message(sprintf(" %d 个乡镇未落入城市，补充最近城市匹配...", length(unmatched)))
unmatched_coords <- st_coordinates(town_cent_sf)[unmatched, , drop = FALSE]
city_coords <- st_coordinates(st_centroid(city))
for (k in seq_along(unmatched)) {
dists <- sqrt(rowSums((city_coords - rep(unmatched_coords[k, ], each = nrow(city_coords)))^2))
matched_city_idx[unmatched[k]] <- which.min(dists)[1]
}
}
town_tbl <- town_base %>%
mutate(
城市索引 = matched_city_idx,
市代码 = city_codes[matched_city_idx],
市名称 = city_names[matched_city_idx],
省名称 = city_provs[matched_city_idx]
) %>%
select(乡镇代码 = code, 乡镇名称 = Name, 市代码, 市名称, 省名称, cx, cy)
town_tbl
message(sprintf(" 匹配完成：%d 个乡镇", nrow(town_tbl)))
|
第 4 步：构建格网空间索引（加速查询）
如果对每个城市对都遍历全部 43,366 个乡镇，计算量将非常大。为此引入格网空间索引：将整个研究区域划分为 400 km × 400 km 的网格，每个乡镇根据其质心坐标归入对应格网。
message("\n[4/9] 构建格网空间索引...")
town_tbl <- town_tbl %>%
mutate(
gx = as.integer(floor(cx / GRID_SIZE)),
gy = as.integer(floor(cy / GRID_SIZE))
)
grid_list <- split(town_tbl$乡镇代码, paste(town_tbl$gx, town_tbl$gy, sep = "_"))
|
对于任意一个城市对，只需提取其共同边界外接矩形（bounding box）覆盖的格网内的乡镇作为候选集，再在这些候选乡镇中进行精确的空间相交判断：
此处代码需下载讲义材料查看～
这一步能将每次判断的候选乡镇数量从 4 万多个降低到数百个，大幅缩短运行时间。
第 5 步：识别相邻城市对
使用 sf 包的 st_touches() 函数识别哪些城市多边形在空间上相邻（即存在公共边界）：
message("\n[5/9] 识别相邻城市对...")
adj <- st_touches(city, city, sparse = TRUE)
adj_df <- tibble(
idx1 = rep(seq_len(nrow(city)), lengths(adj)),
idx2 = unlist(adj)
) %>%
filter(idx1 < idx2) %>%
mutate(
市代码1 = city$市代码[idx1], 市名称1 = city$市名称[idx1],
市代码2 = city$市代码[idx2], 市名称2 = city$市名称[idx2],
城市对 = paste0(市代码1, "-", 市代码2)
)
adj_df
message(sprintf(" 相邻城市对：%d", nrow(adj_df)))
|
关键说明：
st_touches() 返回”接触”关系（多边形有公共边界但不重叠）；sparse = TRUE 返回列表格式，每个城市对应一个相邻城市索引向量；filter(idx1 < idx2) 确保每个城市对只记录一次；- 最终得到 955 对相邻城市。
第 6 步：提取城市共同边界
对每一对相邻城市，用 st_intersection() 计算两个城市多边形的几何交集。如果交集包含线要素（LINESTRING），即为两城市的共同边界：
# ── 6. 提取共同边界
# 此处代码需下载讲义材料查看～
|
关键说明：
st_intersection(A, B) 返回两个几何对象的交集；st_collection_extract(..., "LINESTRING") 从交集的 GEOMETRYCOLLECTION 中提取线段部分；st_union() 将可能断裂的多段边界合并为一条完整的共同边界线；- 使用
tryCatch() 包裹，避免个别异常城市对导致整个循环中断；
- 最终有效共同边界 955 条。
第 7 步：保存城市对共同边界矢量数据
将所有城市对的共同边界几何数据导出为 shp 矢量文件，方便后续直接调用或叠加到地图中使用：
message("\n[6b/9] 保存城市对共同边界矢量数据...")
border_sf <- st_sf(
city_pair = adj_df$城市对,
citycode1 = adj_df$市代码1,
cityname1 = adj_df$市名称1,
citycode2 = adj_df$市代码2,
cityname2 = adj_df$市名称2,
geometry = st_sfc(border_list, crs = proj_crs)
)
border_sf
st_write(border_sf, "res/城市对共同边界.shp", delete_dsn = TRUE, quiet = TRUE)
message(sprintf(" ✓ res/城市对共同边界.shp（%d 条）", nrow(border_sf)))
|
关键说明：
- st_sf() 用于构造包含属性数据和几何数据的空间数据框；
- 城市对 字段为”市代码1-市代码2”格式，可直接与其他数据关联；
- geometry 列存储共同边界的线要素（LINESTRING），坐标系与输入数据一致；
- 使用 delete_dsn = TRUE 自动覆盖同名文件，避免重复追加。
输出文件说明：
| 文件名 |
内容 |
行数 |
城市对共同边界.shp |
所有城市对的共同边界线要素 |
955 条 |
第 8 步：查找边界相交乡镇并计算距离
这是整个流程的核心步骤。对每一对城市，先通过其共同边界的 bbox 利用格网索引快速筛选出候选乡镇，再用 st_intersects() 精确判断哪些乡镇与共同边界相交：
# ── 8. 查找边界相交乡镇
# 此处代码需下载讲义材料查看～
|
找到相交乡镇后，用 st_distance() 计算每个相交乡镇质心到共同边界的最短距离，然后将两侧乡镇交叉配对，形成”乡镇对”级别的分析数据。为控制数据规模，每一侧最多保留距离最近的 500 个乡镇。
最终得到 82,497 对乡镇对，涉及 14,789 个边界相交乡镇。
第 9 步：输出结果文件
所有结果均输出为 Stata .dta 格式，并添加了详细的变量标签，可直接用于 Stata 实证分析：
message("\n[9/9] 写出 DTA...")
add_lab <- function(df, lm) {
for (nm in names(lm))
if (nm %in% names(df)) attr(df[[nm]], "label") <- lm[[nm]]
df
}
valid_results <- result_list[sapply(result_list, function(x) !is.null(x) && nrow(x) > 0)]
if (length(valid_results) == 0) {
result_df <- tibble(
乡镇对 = character(), 城市对 = character(),
乡镇代码1 = character(), 乡镇名称1 = character(),
所在城市代码1 = character(), 所在城市代码2 = character(),
所在城市名称1 = character(), 所在城市名称2 = character(),
距离1_米 = numeric(), 距离2_米 = numeric()
)
} else {
result_df <- bind_rows(valid_results) %>%
mutate(乡镇对 = paste0(乡镇代码1, "-", 乡镇代码2)) %>%
select(乡镇对, 城市对, 乡镇代码1, 乡镇名称1, 乡镇代码2, 乡镇名称2,
所在城市代码1, 所在城市代码2,
所在城市名称1, 所在城市名称2,
距离1_米, 距离2_米)
}
message(sprintf("\n 总乡镇对：%d", nrow(result_df)))
result_out <- result_df %>%
add_lab(list(
乡镇对="乡镇对（乡镇代码1-乡镇代码2）",
城市对="所在城市对（市代码1-市代码2）",
乡镇代码1="乡镇1行政区划代码（12位）",
乡镇名称1="乡镇1名称",
乡镇代码2="乡镇2行政区划代码（12位）",
乡镇名称2="乡镇2名称",
所在城市代码1="乡镇1所在城市代码（6位）",
所在城市代码2="乡镇2所在城市代码（6位）",
所在城市名称1="乡镇1所在城市名称",
所在城市名称2="乡镇2所在城市名称",
距离1_米="乡镇1质心到城市对共同边界距离（米）",
距离2_米="乡镇2质心到城市对共同边界距离（米）"
))
result_out
write_dta(result_out, "res/乡镇对边界距离.dta", label = "数据处理：微信公众号 RStata")
message(sprintf(" ✓ res/乡镇对边界距离.dta（%d 行）", nrow(result_out)))
|
输出文件说明：
| 文件名 |
内容 |
行数 |
乡镇对边界距离.dta |
乡镇对-城市对距离矩阵 |
82,497 行 |
边界相交乡镇.dta |
所有边界相交乡镇及距离 |
14,789 行 |
乡镇代码名称.dta |
乡镇代码-名称对照表 |
43,366 行 |
城市代码名称.dta |
城市代码-名称-省份对照表 |
371 行 |
专题地图展示
除了数据处理代码，还提供了两组地图绘制代码，分别展示全国概览和局部放大效果。
全国城市共同边界及边界乡镇分布图
代码文件：china_border_map_touching.R
该图使用 Albers 等面积投影，底图为全国城市行政区划，乡镇多边形按照”距边界距离”用 scico 色彩渐变填色，颜色越深表示距离边界越远。
完整绘图代码如下：
附件中提供了 r-china-map.zip，把代码中的 r-china-map 路径换成自己电脑上的即可。
source("/Users/ac/.Rprofile")
skill_path <- "/Users/ac/.workbuddy/skills/r-china-map"
data_path <- file.path(skill_path, "data")
suppressPackageStartupMessages({
library(tidyverse)
library(sf)
library(haven)
library(ggspatial)
library(scico)
})
set.seed(42)
mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
message("[1/3] 读取数据...")
dist_df <- read_dta("res/乡镇对边界距离.dta")
city_raw <- st_read("2021行政区划/市.shp", quiet = TRUE)
town_raw <- readRDS("town.rds")
message("[2/3] 读取地图...")
plotbbox <- st_bbox(
c(xmin = -2725586, xmax = 2982768, ymax = 6000000, ymin = 1800655),
crs = st_crs(mycrs)
) %>% st_as_sfc()
citymap <- read_sf(file.path(data_path, "chinacity2021mini/chinacity2021mini.shp")) %>%
filter(!is.na(省代码))
cityline <- read_sf(file.path(data_path, "chinacity2021mini/chinacity2021mini_line.shp")) %>%
filter(class %in% c("九段线", "海岸线", "小地图框格", "省界")) %>%
select(class)
china_neigh <- read_sf(file.path(data_path, "china_neighboring/china_neighboring.shp"))
message("[3/3] 处理乡镇数据...")
town_touching <- bind_rows(
dist_df %>% select(乡镇代码 = 乡镇代码1, 所在城市代码 = 所在城市代码1, 距离_米 = 距离1_米),
dist_df %>% select(乡镇代码 = 乡镇代码2, 所在城市代码 = 所在城市代码2, 距离_米 = 距离2_米)
) %>% distinct() %>%
mutate(乡镇代码 = as.character(乡镇代码))
town_sf <- town_raw %>%
rename(乡镇代码 = code, 乡镇名称 = Name) %>%
mutate(乡镇代码 = as.character(乡镇代码)) %>%
filter(乡镇代码 %in% town_touching$乡镇代码) %>%
st_transform(crs = mycrs)
valid_idx <- st_is_valid(town_sf) & !st_is_empty(town_sf)
town_sf_valid <- town_sf[which(valid_idx), ]
town_poly <- town_sf_valid %>%
mutate(乡镇代码 = as.character(乡镇代码)) %>%
left_join(town_touching %>% mutate(乡镇代码 = as.character(乡镇代码)),
by = "乡镇代码") %>%
mutate(距离_km = as.numeric(距离_米) / 1000) %>%
filter(!is.na(距离_km))
ggsave("res/全国边界乡镇分布.png", p_main,
width = 12, height = 9, device = png, dpi = 400)
|
输出图件：
成都-德阳城市对专题放大图
代码文件：chengdu_deyang_map.R
针对具体城市对（以成都-德阳为例），可以将共同边界区域放大展示。该图分别用不同颜色填充两个城市范围，用粗蓝线标注共同边界，并在乡镇上标注名称。
完整绘图代码如下：
source("/Users/ac/.Rprofile")
suppressPackageStartupMessages({
library(tidyverse)
library(sf)
library(haven)
library(ggspatial)
})
set.seed(42)
mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
message("[1/4] 读取数据...")
dist_df <- read_dta("res/乡镇对边界距离.dta")
city_raw <- st_read("2021行政区划/市.shp", quiet = TRUE)
town_raw <- readRDS("town.rds")
message("\n[2/4] 筛选成都-德阳城市对...")
chengdu_deyang <- dist_df %>%
filter((所在城市名称1 == "成都市" & 所在城市名称2 == "德阳市") |
(所在城市名称1 == "德阳市" & 所在城市名称2 == "成都市"))
chengdu_towns <- chengdu_deyang %>%
select(乡镇代码 = 乡镇代码1, 乡镇名称 = 乡镇名称1,
所在城市名称 = 所在城市名称1, 距离_米 = 距离1_米) %>%
distinct()
deyang_towns <- chengdu_deyang %>%
select(乡镇代码 = 乡镇代码2, 乡镇名称 = 乡镇名称2,
所在城市名称 = 所在城市名称2, 距离_米 = 距离2_米) %>%
distinct()
border_towns <- bind_rows(chengdu_towns, deyang_towns) %>% distinct()
message("\n[3/4] 处理空间数据...")
city_prj <- city_raw %>%
select(市代码, 市名称 = 市, 省名称 = 省) %>%
st_transform(crs = mycrs)
city_chengdu <- city_prj %>% filter(市名称 == "成都市")
city_deyang <- city_prj %>% filter(市名称 == "德阳市")
common_border <- st_intersection(st_geometry(city_chengdu), st_geometry(city_deyang))
common_border_line <- st_collection_extract(common_border, "LINESTRING") %>%
st_union() %>% st_sfc(crs = mycrs)
town_prj <- town_raw %>%
rename(乡镇代码 = code, 乡镇名称 = Name) %>%
mutate(乡镇代码 = as.character(乡镇代码)) %>%
st_transform(crs = mycrs)
town_chengdu_deyang <- town_prj %>%
filter(乡镇代码 %in% border_towns$乡镇代码) %>%
left_join(border_towns %>%
mutate(乡镇代码 = as.character(乡镇代码)) %>%
select(乡镇代码, 所在城市名称, 距离_米),
by = "乡镇代码") %>%
mutate(距离_km = as.numeric(距离_米) / 1000)
valid_idx <- st_is_valid(town_chengdu_deyang) & !st_is_empty(town_chengdu_deyang)
town_valid <- town_chengdu_deyang[which(valid_idx), ]
ggsave("res/成都德阳边界地图.png", p,
width = 12, height = 9, device = png, dpi = 400)
|
输出图件：
同时生成了成都-德阳边界乡镇列表（res/成都德阳边界乡镇列表.csv），方便查看具体乡镇名称和距离。
点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 R 语言提取城市共同边界、边界乡镇识别及距离计算
评论