城市边界和农村边界这两个概念，在地理信息处理里经常被混着用，但真正落到数据层面，它们完全是两套逻辑。我最近在做一个区域统计分析的小项目，需要把某片城区的周边农村范围圈出来，手头正好有一份全球城市边界数据集，原本以为直接拿城市多边形往外扩一圈就完事了，实测下来发现坑不少。这篇就聊聊我是怎么用 Python 把“城市边界”转换成“城市附近的农村边界”的，包括数据理解、空间运算的选型、具体实现步骤，以及踩过的那些坑。适合有一定 Python 基础、做过地理数据处理、或者正在做城乡结合部相关分析的朋友参考。哪怕你之前没碰过空间数据，跟着思路走也能明白每一步在干什么。
1. 先搞清楚数据长什么样再动手
1.1 全球城市边界数据的典型结构
这类全球城市边界数据集，通常是以矢量多边形（Polygon 或 MultiPolygon）的形式组织的，每个要素代表一个城市建成区的外轮廓。属性表里一般会有城市名称、所属区域、面积、人口之类的字段，几何字段则是那个不规则的多边形。我第一次打开数据的时候，习惯性地先看坐标系——这一步千万别省。很多全球尺度的数据集用的是地理坐标系（比如 WGS84，单位是度），而你要做的缓冲区、距离计算如果直接按度来算，结果会离谱到没法看。因为一度经度在不同纬度对应的实际距离差别巨大，赤道附近一度约 111 公里，到了高纬度可能只剩几十公里。
所以拿到数据的第一件事，我一般会做三件事：看坐标系、看几何类型、看属性字段。用geopandas读进来之后，gdf.crs看坐标系，gdf.geom_type看几何类型，gdf.columns看字段。如果坐标系是地理坐标，后面所有涉及距离的运算之前，必须先投影到合适的投影坐标系。这一步是后面所有操作的地基，地基歪了，上面盖什么都是歪的。
提示：不要假设数据集一定是干净的。全球尺度的数据经常存在自相交多边形、空几何、重复要素，这些在后续空间运算里会直接报错或者产生错误结果，读进来之后先做一次几何有效性检查。
1.2 “农村边界”到底该怎么定义
这是整个项目里最需要想清楚的问题，也是标题里那个“对应城市附近的农村边界”最容易含糊的地方。严格来说，农村边界并没有一个全球统一的、现成的定义。城市边界数据集只告诉你“哪里是城市建成区”，它不会告诉你“哪里是农村”。所以我们要做的，本质上是基于城市边界去推导一个“城市周边、但不在城市内部”的范围，把它当作农村区域的近似。
这里有两种常见思路。第一种是“缓冲区法”：以城市多边形为中心，向外扩一个固定距离（比如 10 公里、20 公里），得到一个大范围，再把这个大范围减去城市本身，剩下的环状区域就近似为城市附近的农村。第二种是“邻近行政区法”：如果有行政区划数据，取与城市相交或相邻的行政单元，再减去城市建成区。第一种方法不依赖额外数据，纯靠几何运算，通用性强，适合快速出结果；第二种更贴近真实的行政边界，但需要额外的行政区划数据源，数据获取和匹配的成本高不少。
我这次选的是缓冲区法，原因很直接：手头只有城市边界数据，没有可靠的行政区划数据，而且项目需求只是要一个“城市附近农村”的近似范围用于统计，不要求精确到行政村级别。选型的时候一定要想清楚你的下游用途——如果是要做精确的行政统计，缓冲区法出来的结果是不能直接用的；如果只是做空间分布的可视化或者粗略的面积估算，缓冲区法完全够用。
1.3 为什么不用简单的“外扩一圈”
有人可能会想，直接把城市多边形往外 offset 一圈不就行了？理论上缓冲区就是这么干的，但实际操作里有个关键细节：城市多边形往往是不规则的，有的地方凹进去，有的地方凸出来，简单外扩之后，凹进去的部分可能会被“填平”，导致缓冲区把一些本来不属于城市周边的区域也包进来。更麻烦的是，多个城市挨得近的时候，各自的缓冲区会重叠，重叠区域如果不去重，面积统计就会重复计算。
所以我的处理流程里，缓冲区生成之后一定要做两件事：一是把缓冲区结果做一次融合（dissolve），把重叠部分合并掉；二是用缓冲区减去原始城市多边形，得到纯粹的“城市外、缓冲区内的环状区域”。这个环状区域，才是我们想要的“城市附近农村边界”的几何表达。听起来简单，但每一步都有坑，下面细说。
2. 核心实现：从城市多边形到农村环状区域
2.1 环境准备与依赖选型
这套流程用到的核心库就三个：geopandas负责矢量数据的读写和属性操作，shapely负责几何运算（缓冲区、差集这些），pyproj负责坐标系转换。如果你还要做可视化，可以加matplotlib。安装的时候建议用 conda 或者 pip 装，geopandas依赖 GDAL 那套底层库，有时候 pip 直接装会出问题，conda 装省心一些。
conda install -c conda-forge geopandas shapely pyproj matplotlib选型上我为什么不用fiona单独读数据？因为geopandas已经把fiona封装好了，gpd.read_file()一行就能读，没必要自己再调底层。为什么不用rasterio？因为我们的数据是矢量，不是栅格，用不上。工具选型的原则就是：能用高层封装就不碰底层，除非你有特殊的性能或者格式需求。
2.2 读取数据与坐标系处理
读数据这一步，我习惯先打印基本信息，确认数据没问题再往下走。
import geopandas as gpd # 读取城市边界数据 cities = gpd.read_file("city_boundaries.shp") print(cities.crs) print(cities.geom_type.value_counts()) print(cities.head())如果cities.crs显示的是类似EPSG:4326的地理坐标系，那就要投影。投影坐标系的选择有个原则：选和你研究区域纬度匹配的等距投影，或者直接用该区域常用的投影带。全球数据的话，可以用EPSG:3857（Web 墨卡托），但它在高纬度地区面积变形很大；更稳妥的做法是按城市所在纬度选 UTM 带。不过为了简化，如果只是做近似分析，EPSG:3857也能凑合，但你要清楚它的变形问题。
# 投影到米制坐标系，方便按米计算缓冲区 cities_proj = cities.to_crs("EPSG:3857")这里有个实操心得：投影之后，几何的坐标单位就变成米了，后面buffer(10000)就是 10 公里，直观很多。如果忘了投影直接 buffer，buffer(0.1)在度制下大概是 11 公里（赤道附近），但到了高纬度就完全不是这个数，结果不可控。
2.3 缓冲区距离的参数选择与计算
缓冲区距离选多大，取决于你对“城市附近”的定义。我这次用的是 10 公里，理由是：对于大多数中小城市，10 公里基本能覆盖到城郊结合部和近郊农村；对于大城市，10 公里可能还在城区范围内，这时候可能需要 20 公里甚至更大。这个参数没有标准答案，你得结合研究目的和城市规模来定。
计算缓冲区之前，先做几何有效性修复，这一步能避免后面大量莫名其妙的报错。
from shapely.validation import make_valid # 修复无效几何 cities_proj["geometry"] = cities_proj["geometry"].apply( lambda geom: make_valid(geom) if not geom.is_valid else geom ) # 生成缓冲区 buffer_dist = 10000 # 单位：米 cities_proj["buffer_geom"] = cities_proj["geometry"].buffer(buffer_dist)buffer操作本身不复杂，但要注意：如果某个城市多边形是 MultiPolygon 且各部分离得很远，buffer 之后可能产生多个分离的环，这在后面做差集的时候要留意。另外，buffer 的默认分辨率是 16，对于大范围分析够用了，如果你要更平滑的边界可以调高，但计算量会上去。
2.4 融合缓冲区并减去城市本体
这一步是整个流程的核心。先把所有城市的缓冲区融合成一个大几何，去掉重叠部分，然后再减去所有城市多边形的并集。
# 融合所有缓冲区 buffer_union = cities_proj["buffer_geom"].unary_union # 融合所有城市本体 city_union = cities_proj["geometry"].unary_union # 差集：缓冲区减去城市本体，得到农村环状区域 rural_ring = buffer_union.difference(city_union)unary_union这个操作会把所有几何合并成一个，重叠部分自动去重。difference就是几何差集，A 减去 B，得到 A 中不属于 B 的部分。这两步做完，rural_ring就是我们要的“城市附近农村边界”的几何表达。
这里有个容易忽略的点：如果两个城市离得很近，它们的缓冲区融合后会连成一片，差集之后得到的农村环状区域也会是连通的。这在几何上没问题，但如果你后续要按城市分别统计农村面积，就需要在差集之前先按城市分组处理，而不是一股脑融合。我这次的需求是整体范围，所以融合处理没问题；如果你要分城市统计，记得改成按城市逐个 buffer、逐个 difference。
2.5 结果导出与可视化验证
算完之后，把结果存成矢量文件，方便后续在 GIS 软件里查看或者做进一步分析。
# 构造结果 GeoDataFrame rural_gdf = gpd.GeoDataFrame( {"type": ["rural_ring"]}, geometry=[rural_ring], crs=cities_proj.crs ) # 导出 rural_gdf.to_file("rural_boundary.geojson", driver="GeoJSON")导出之后，我强烈建议做一次可视化验证。把城市多边形和农村环状区域叠在一起画出来，肉眼看一下环状区域是不是合理——有没有把明显不该包进来的区域包进来，有没有因为缓冲区太小导致环状区域断掉。可视化用matplotlib几行就能搞定。
import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(10, 10)) cities_proj.plot(ax=ax, color="red", alpha=0.5, label="City") rural_gdf.plot(ax=ax, color="green", alpha=0.3, label="Rural Ring") plt.legend() plt.savefig("validation.png", dpi=150)这一步别偷懒。我见过太多人算完直接导出，结果打开一看环状区域缺了一大块，回头排查半天才发现是某个城市几何无效导致 buffer 失败。可视化是最快的验证手段。
3. 实操中必须注意的细节与参数
3.1 几何有效性与空几何的处理
全球尺度的城市边界数据，几何无效是常态。自相交、重复点、环方向错误，这些都会让buffer和difference产生意外结果。make_valid能修复大部分问题，但它不是万能的，修复后可能把 Polygon 变成 MultiPolygon 或者 GeometryCollection，后续处理要能兼容这些类型。
空几何也要单独处理。有些城市的几何字段可能是空的，buffer对空几何会返回空，unary_union遇到空几何一般会跳过，但如果你在循环里逐个处理，空几何会导致报错。我的做法是读进来之后先过滤掉空几何：
cities_proj = cities_proj[~cities_proj["geometry"].is_empty]3.2 缓冲区距离的单位陷阱
前面提过，缓冲区距离的单位取决于坐标系。投影坐标系下是米，地理坐标系下是度。这个坑我踩过不止一次：有一次数据没投影，直接buffer(0.1)，心想 0.1 度大概 11 公里，结果那个城市在高纬度，实际缓冲区只有几公里，环状区域窄得几乎看不见。所以记住一条铁律：做任何距离相关的运算之前，先确认坐标系单位。
3.3 多城市重叠缓冲区的去重逻辑
多个城市缓冲区重叠时，unary_union会自动去重，这是它的优势。但如果你先逐个城市算差集再合并，重叠区域会被重复计算两次，面积统计就偏大了。所以顺序很重要：先融合缓冲区，再融合城市本体，最后做一次差集。这个顺序能保证重叠区域只算一次。
3.4 大范围数据的性能优化
全球数据动辄几千上万个城市多边形，逐个 buffer 再 union，计算量不小。如果跑起来太慢，有两个优化方向：一是先用空间索引（geopandas的sindex）筛出研究范围内的城市，减少参与运算的要素数量；二是把 buffer 的分辨率调低（默认 16，调到 8 能明显提速，精度损失在粗略分析里可以接受）。我这次数据量不大，没做优化，但如果你要处理全球全量数据，这两招能省不少时间。
4. 常见问题排查速查表
实际跑这套流程，遇到的问题基本集中在几何运算和坐标系上。我把踩过的坑整理成表，方便对照排查。
| 问题现象 | 可能原因 | 排查与解决 |
|---|---|---|
| buffer 结果范围明显偏小或偏大 | 坐标系是地理坐标，单位是度 | 检查gdf.crs，投影到米制坐标系后再 buffer |
| difference 报拓扑错误 | 几何自相交或无效 | 先make_valid修复，再运算 |
| 环状区域出现断裂或缺口 | 缓冲区距离太小，城市间空隙未被覆盖 | 增大 buffer 距离，或检查是否有城市几何缺失 |
| 面积统计偏大 | 重叠缓冲区被重复计算 | 确保先 union 再 difference，不要逐个算再合并 |
| 导出文件在 GIS 软件里打不开 | 坐标系丢失或格式不兼容 | 导出时显式指定crs，优先用 GeoJSON 或 GPKG |
| 运算速度极慢 | 要素数量多，逐个运算 | 用空间索引预筛选，或降低 buffer 分辨率 |
除了表里这些，还有一个隐蔽的问题：unary_union在处理大量几何时可能返回 GeometryCollection 而不是 MultiPolygon，如果你的下游代码只接受 MultiPolygon，需要额外做一次类型转换。我一般会在差集之后检查一下几何类型，必要时用geopandas的explode或者手动提取多边形部分。
5. 这套方法能用到哪些场景
缓冲区法生成农村边界的思路，适用范围比想象中广。城乡结合部的土地利用变化监测是一个典型场景——把城市边界和农村环状区域叠起来，看环状区域内的建设用地扩张速度，能反映城市向外蔓延的趋势。另一个场景是公共服务设施的覆盖分析，比如分析农村环状区域内学校、卫生室的分布密度，判断城郊农村的公共服务是否跟得上。
再往大了说，这套几何运算的套路可以迁移到很多“基于一个区域推导周边区域”的问题上。比如给定一个保护区边界，推导周边的缓冲带用于生态影响评估；给定一个污染源点位，推导周边的影响范围。核心逻辑都是一样的：缓冲区生成、融合去重、差集提取。把这一套吃透，换个数据换个场景，改改参数就能复用。
我个人在实际操作中的体会是，这类空间推导的活儿，难点从来不在代码本身，代码就那么几行，难的是想清楚“你要的边界到底怎么定义”以及“数据本身靠不靠谱”。定义想清楚了，坐标系确认了，几何修复做了，剩下的就是按部就班地跑。反倒是那些急着写代码、跳过数据检查的，最后往往在排查几何错误上花掉更多时间。所以我的习惯是，动手之前先花十分钟把数据摸清楚，这十分钟能省掉后面一小时的 debug。
排查指南)
