名师讲堂|使用 Python 提取城市共同边界、边界乡镇识别及距离计算

由于借助 AI 工具学习编程已经变得非常容易了，因此之后的课程就不再默认进行视频讲解了，如果特别需要视频讲解也可以联系李老师预约讲解～讲义材料学习过程中遇到的问题也可以及时与李老师联系。

购买 RStata 名师讲堂会员即可参加该课程啦（之前的和未来的都可以参加）！

价格：2800/年 或者 4800/长期

购买会员可以从这里下单：https://rstata.duanshu.com/#/card/list/

名师讲堂会员权益：

1. 参加每个月 3～4 次的名师讲堂课程；
2. 参加平台上的其他 R 语言和 Stata 的课程；
3. 以会员折扣价购买我们分享的数据资料（10 元/份）；
4. 课程内外的提问解答服务（课程外的尽量帮忙解决）。

* 如果发票可添加小编微信 r_stata2 （RStata 李老师）开具。如需数据资料，购买后可添加小编微信免费领取数据折扣卡。

更多关于 RStata 会员的更多信息可添加微信号 r_stata2 咨询：

课程主页（点击文末的阅读原文即可跳转）：https://rstata.duanshu.com/#/brief/course/fd6b33d0630545daafa61ecd607a55f5

今天给大家分享使用 Python 提取城市共同边界、识别边界相交乡镇并计算各乡镇到边界距离的完整方法。该指标可用于研究城市交界地区的经济发展、人口流动、环境溢出等空间边界效应，是空间计量与城市经济学研究中的重要变量。

处理过程基于 2021 年行政区划矢量数据，识别出全国 955 对相邻城市，提取其共同边界，并计算出各边界乡镇到共同边界的距离。最终产出包括：城市共同边界矢量文件、边界相交乡镇列表及距离矩阵，可直接用于 Stata 或 Python 的实证分析。

指标介绍

城市共同边界：指两个相邻城市行政区划多边形相交形成的线段（shared border）。例如成都市与德阳市相邻，两者的行政边界相交形成一条共同边界线。

边界相交乡镇：指行政边界多边形与城市共同边界线相交的乡镇（街道）。这些乡镇位于城市交界处，是研究"边界效应"（border effect）的核心样本。

距边界距离：指乡镇行政中心（质心）到最近的城市共同边界线的直线距离（米），用于构造连续的空间边界距离变量。

乡镇对-城市对距离矩阵：对于每一对相邻城市，分别列出两侧所有边界乡镇及其到共同边界的距离，形成"乡镇对"级别的分析数据。

整个处理流程分为九个步骤，下面逐一介绍。

使用 reticulate 创建与管理 Python 虚拟环境

在 R 中通过 reticulate 包来调用 Python，最好的实践是为项目创建一个专属的 Python 虚拟环境，将所需依赖隔离到独立空间，避免与系统 Python（如 Anaconda）发生版本冲突。

重要说明（避免"已初始化"报错）：reticulate 在 R 会话中只能绑定一次 Python——一旦某个 {python} 代码块运行，Python 解释器就被锁定，之后再调用 use_virtualenv() 会报错：

ERROR: The requested version of Python cannot be used, as another version has already been initialized.

因此，虚拟环境的激活必须在所有 {python} 代码块之前完成。本文档的解决方案是在 setup chunk 中通过 Sys.setenv(RETICULATE_PYTHON = ...) 提前锁定 Python 路径，这是 reticulate 选取 Python 的最高优先级入口。

安装 reticulate（仅首次）

虚拟环境初始化原理（已在 setup chunk 中完成）

本文档的 setup chunk（隐藏运行）包含如下逻辑：

这样做的关键在于：knitr 在处理第一个 {python} chunk 时，reticulate 已经通过 RETICULATE_PYTHON 环境变量知道要使用 .venv，不会再去碰 Anaconda。

在虚拟环境中安装 Python 包（仅首次）

验证激活状态

查看已安装的包

虚拟环境管理常用命令

数据准备：将 town.rds 转为 GeoPackage

Python 无法直接读取 R 的 .rds 格式，需在 R 中先将乡镇矢量数据转出为通用格式（GeoPackage）：

第 1 步：读取数据

首先加载 Python 包，读取 2021 年城市级和乡镇级行政区划矢量数据，并统一投影坐标系（Albers 等面积投影）：

关键说明：

• town.gpkg 为从 town.rds 转换而来的乡镇级矢量数据（含 43,366 个乡镇多边形），已使用 Albers 等面积投影；
• 城市数据同步转换到相同 CRS，确保后续空间运算准确无误；
• MAX_TOWNS_PER_SIDE = 500 控制每侧最多保留的乡镇数量；
• GRID_SIZE = 400000 为格网大小（400 km），用于后续空间索引加速。

第 2 步：预计算乡镇质心坐标

为提高后续大规模距离计算的速度，先一次性计算所有乡镇的质心坐标并存入数据框：

关键说明：

• geometry.centroid 计算每个乡镇多边形的质心（几何中心）；
• .x 和 .y 属性直接提取质心坐标；
• 将几何信息转为普通数据框列（cx、cy），为后续格网索引和空间匹配做准备。

第 3 步：判断乡镇所属城市

乡镇矢量数据本身不包含"所属城市"字段，用空间连接（sjoin within）将乡镇质心与城市多边形匹配：

第 4 步：构建格网空间索引（加速查询）

将整个研究区域划分为 400 km × 400 km 的网格，每个乡镇根据其质心归入对应格网，大幅缩短候选乡镇的筛查时间：

第 5 步：识别相邻城市对

使用 Shapely 的 .touches() 方法识别哪些城市多边形在空间上相邻（存在公共边界）：

关键说明：

• .touches() 判断两个几何对象是否有公共边界但不重叠，等价于 R 的 st_touches()；
• range(i+1, n_cities) 确保每个城市对只记录一次；
• 最终得到 955 对相邻城市。

第 6 步：提取城市共同边界

对每一对相邻城市，用 intersection() 计算两个城市多边形的几何交集，并从中提取线要素：

关键说明：

• intersection() 计算两个几何对象的交集，等价于 R 的 st_intersection()；
• 从 GeometryCollection 中提取线要素，等价于 R 的 st_collection_extract(..., "LINESTRING")；
• unary_union() 合并多段边界，等价于 R 的 st_union()。

第 6b 步：保存城市对共同边界矢量数据

将所有城市对的共同边界几何数据导出为 shp 矢量文件：

关键说明：

• GeoDataFrame() 构造包含属性和几何数据的空间数据框，等价于 R 的 st_sf()；
• city_pair 字段为"市代码1-市代码2"格式，可直接与其他数据关联；
• shp 格式字段名不超过 10 字符，因此使用短英文名。

输出文件说明：

第 7 步：构建乡镇空间数据（含所属城市）

将乡镇矢量数据与城市匹配结果合并，为后续相交判断做准备：

第 8 步：查找边界相交乡镇并计算距离

这是整个流程的核心步骤。对每对城市，用格网索引快速筛选候选乡镇，再用 .intersects() 精确判断，最后计算各乡镇质心到共同边界的距离：

找到相交乡镇后，用 .distance() 计算每个相交乡镇质心到共同边界的最短距离，然后将两侧乡镇交叉配对，形成"乡镇对"级别的分析数据。最终得到 82,497 对乡镇对，涉及 14,789 个边界相交乡镇。

第 9 步：输出结果文件

所有结果均输出为 Stata .dta 格式（version=118 支持 UTF-8 中文），并添加了详细的变量标签：

输出文件说明：

专题地图展示

全国城市共同边界及边界乡镇分布图

代码文件：china_border_map_touching.py

该图使用 draw_china_map.py 的辅助函数绘制专业级中国地图，包含南海小地图、省名标注、九段线/海岸线、比例尺和指北针等标准地图元素。乡镇多边形按距边界距离用 Lipari 渐变色填色。

输出图件：

成都-德阳城市对专题放大图

代码文件：chengdu_deyang_map.py

输出图件：

如何参加课程？

购买 RStata 名师讲堂会员即可参加该课程啦（之前的和未来的都可以参加）！

价格：2800/年 或者 4800/长期

购买会员可以从这里下单：https://rstata.duanshu.com/#/card/list/

名师讲堂会员权益：

1. 参加每个月 3～4 次的名师讲堂课程；
2. 参加平台上的其他 R 语言和 Stata 的课程；
3. 以会员折扣价购买我们分享的数据资料（10 元/份）；
4. 课程内外的提问解答服务（课程外的尽量帮忙解决）。

* 如果发票可添加小编微信 r_stata2 （RStata 李老师）开具。如需数据资料，购买后可添加小编微信免费领取数据折扣卡。

更多关于 RStata 会员的更多信息可添加微信号 r_stata2 咨询：

课程主页（点击文末的阅读原文即可跳转）：https://rstata.duanshu.com/#/brief/course/fd6b33d0630545daafa61ecd607a55f5