# GIS数据实战：如何用Python快速处理全国乡镇级Shapefile与GeoJSON数据 最近在做一个智慧社区的项目，需要分析全国乡镇层级的服务设施覆盖情况。拿到数据的那一刻，我有点懵——一个包含全国近四万个乡镇的Shapefile文件，解压后足足有十几个G，用传统的桌面GIS软件打开都卡得不行，更别说进行复杂的空间分析和属性筛选了。相信不少GIS开发者和数据分析师都遇到过类似的困境：数据体量巨大，处理流程繁琐，传统工具在效率和灵活性上捉襟见肘。 这正是Python生态大显身手的地方。通过`geopandas`、`shapely`、`fiona`等库构建的轻量级处理流水线，我们完全可以在Jupyter Notebook里，用几行清晰的代码完成从数据加载、清洗、转换到分析的全过程，将原本需要数小时的手动操作压缩到几分钟内。这篇文章，我就结合自己处理全国乡镇级行政区划数据（包括Shapefile和GeoJSON格式）的实际经验，分享一套高效、可复用的Python实战技巧。无论你是要筛选特定省份的乡镇、进行坐标系统一，还是做空间合并与聚合，这里都有现成的代码块和思路供你参考。 ## 1. 环境搭建与数据高效加载 工欲善其事，必先利其器。处理大规模地理数据，第一步不是急着写代码，而是搭建一个稳定且资源充足的环境。我强烈建议使用**Anaconda**来管理Python环境，它能很好地解决地理空间库复杂的依赖问题。 首先，我们创建一个专属的conda环境。打开终端（或Anaconda Prompt），执行以下命令： ```bash conda create -n gis_processing python=3.9 conda activate gis_processing ``` 接下来，安装核心的“地理数据科学三件套”。使用conda安装是首选，因为它能自动处理GDAL、PROJ等底层C库的依赖，比pip安装省心很多。 ```bash conda install -c conda-forge geopandas ``` 这一条命令会同时安装好`geopandas`、`pandas`、`shapely`、`fiona`、`pyproj`等所有必需库。为了后续的数据可视化，我们还可以安装`matplotlib`和`contextily`（用于添加在线底图）。 ```bash conda install -c conda-forge matplotlib contextily ``` 环境准备好后，启动Jupyter Notebook，我们就可以开始数据加载了。面对一个巨大的全国乡镇Shapefile文件，直接使用`gpd.read_file()`可能会耗尽内存。这里有几个策略： **策略一：利用过滤条件仅加载所需数据。** 如果我们只需要处理“浙江省”的数据，可以在读取时通过`bbox`（地理范围框）或`where`（属性过滤）参数进行筛选。但Shapefile的过滤有时不如数据库高效。 **策略二（推荐）：先转换为高效格式。** 我的做法是，先将庞大的Shapefile数据加载到**GeoPackage**或**Parquet**格式中。这两种格式支持更快的查询和选择性列读取。下面是一个转换示例： ```python import geopandas as gpd # 1. 首次加载巨型Shapefile（可能需要一些时间） print("正在加载原始Shapefile...") national_townships = gpd.read_file('path/to/全国乡镇.shp') print(f"数据加载完成，共{len(national_townships)}个要素。") # 2. 保存为GeoPackage格式 national_townships.to_file('national_townships.gpkg', driver='GPKG') print("已保存为GeoPackage格式。") # 3. 后续使用中，可以按需读取特定区域 import sqlite3 # 例如，只读取‘province’字段为‘浙江省’的要素 zhejiang_towns = gpd.read_file('national_townships.gpkg', layer='national_townships', where="province='浙江省'") ``` 对于GeoJSON数据，如果文件过大，同样可以考虑先转换为GeoPackage。此外，`geopandas`从`0.12`版本开始加强了对**Parquet**格式的支持，这是一种列式存储格式，对于只需要部分属性的查询速度极快。 > 注意：处理全国性数据时，务必留意你的内存大小。如果数据实在太大（比如超过内存的50%），考虑使用Dask-geopandas进行分块处理，或者将数据导入PostGIS数据库中进行操作。 ## 2. 核心数据处理操作：坐标转换、清洗与筛选 数据加载进来后，通常不能直接使用，需要经过一系列的清洗和转换。全国公开的行政区划数据，常见的坐标系统是**WGS84**（EPSG:4326）。但在进行面积计算、距离测量或与某些在线地图服务叠加时，我们可能需要将其转换为**投影坐标系**。 ### 2.1 坐标系统一与转换 使用`geopandas`转换坐标系统非常简单，一行代码即可。例如，转换为适用于中国区域的**兰伯特等角圆锥投影**（EPSG:4526）或**Web墨卡托**（EPSG:3857，常用于网络地图）。 ```python # 查看当前坐标参考系 print(zhejiang_towns.crs) # 输出应类似 EPSG:4326 # 转换为Web墨卡托投影（用于与在线瓦片地图叠加） zhejiang_towns_web = zhejiang_towns.to_crs(epsg=3857) # 转换为CGCS2000高斯-克吕格投影（例如3度带，中央经线120度，EPSG:4547） # 注意：全国数据需要分带处理，这里以120度带为例 zhejiang_towns_gk = zhejiang_towns.to_crs(epsg=4547) ``` 转换后，几何对象的坐标单位会从度变为米，此时我们就可以正确计算面积和周长了。 ```python # 在投影坐标系下计算每个乡镇的面积（平方米） zhejiang_towns_web['area_sq_m'] = zhejiang_towns_web.geometry.area # 转换为平方公里 zhejiang_towns_web['area_sq_km'] = zhejiang_towns_web['area_sq_m'] / 1_000_000 ``` ### 2.2 属性数据的清洗与筛选 原始数据中，属性字段可能存在命名不规范、值缺失或格式不一致的问题。我们可以像使用`pandas`处理普通表格一样来处理GeoDataFrame的属性表。 * **字段重命名**：将晦涩的字段名改为易懂的英文或拼音。 * **处理缺失值**：检查关键字段（如行政区划代码、名称）是否有空值。 * **值标准化**：例如，确保“类型”字段中的“镇”、“乡镇”统一为“镇”。 ```python # 假设原始属性表有字段‘NAME’、‘CODE’ # 1. 重命名 zhejiang_towns = zhejiang_towns.rename(columns={'NAME': 'town_name', 'CODE': 'town_code'}) # 2. 检查并处理缺失值 missing_name = zhejiang_towns['town_name'].isnull().sum() print(f"乡镇名称缺失的数量：{missing_name}") # 如果缺失，可以用上级名称+序号填充（根据实际情况） # zhejiang_towns['town_name'] = zhejiang_towns['town_name'].fillna(zhejiang_towns['county_name'] + '_未知') # 3. 筛选数据：例如，筛选出所有‘街道’类型的行政区 streets = zhejiang_towns[zhejiang_towns['type'] == '街道'] # 4. 复杂筛选：筛选出面积大于50平方公里的镇 large_towns = zhejiang_towns_web[(zhejiang_towns_web['type'] == '镇') & (zhejiang_towns_web['area_sq_km'] > 50)] ``` 为了更直观地对比不同数据处理操作，可以参考下表： | 操作类型 | 常用方法/函数 | 关键参数/说明 | 应用场景举例 | | :--- | :--- | :--- | :--- | | **坐标转换** | `GeoDataFrame.to_crs()` | `epsg` (目标坐标系代码) | 将WGS84数据转为投影坐标以计算面积。 | | **属性筛选** | 布尔索引 `df[condition]` | 基于字段条件创建布尔序列 | 提取特定省份、特定类型的乡镇。 | | **空间筛选** | `gpd.sjoin()` | `how` (‘inner’, ‘left’), `predicate` | 筛选出落在某个城市边界内的所有乡镇。 | | **几何简化** | `GeoSeries.simplify()` | `tolerance` (简化容差) | 减少数据量，用于快速预览或网络传输。 | | **几何修复** | `GeoSeries.buffer(0)` | 距离设为0 | 修复无效的几何图形（如自相交）。 | ## 3. 空间关系分析与数据融合 处理行政区划数据时，单要素的分析往往不够，我们经常需要分析乡镇之间的空间关系，或者将乡镇数据与其他图层（如河流、道路、兴趣点）进行叠加分析。 ### 3.1 空间连接：为乡镇附加外部属性 假设我们有一个全国学校点的GeoJSON文件，想要统计每个乡镇内有多少所学校。`geopandas`的`spatial join`功能可以轻松实现。 ```python # 加载学校点数据 schools = gpd.read_file('path/to/schools.geojson') # 确保学校数据与乡镇数据为同一坐标系 schools = schools.to_crs(zhejiang_towns.crs) # 进行空间连接，使用‘within’谓词（点在学校内） towns_with_schools = gpd.sjoin(zhejiang_towns, schools, how='left', predicate='contains') # 分组统计每个乡镇的学校数量 school_count = towns_with_schools.groupby('town_code').size().reset_index(name='school_count') # 将统计结果合并回原乡镇数据 zhejiang_towns = zhejiang_towns.merge(school_count, on='town_code', how='left') zhejiang_towns['school_count'] = zhejiang_towns['school_count'].fillna(0) # 填充NaN为0 ``` ### 3.2 数据融合：乡镇合并为区县 有时，我们需要将乡镇级的统计数据聚合到区县级。这涉及到两个步骤：**空间融合**和**属性聚合**。 * **空间融合**：使用`dissolve`方法，根据某个字段（如区县代码）将相邻的乡镇面合并成大的区县面。 * **属性聚合**：在融合过程中，对乡镇的其他属性（如人口、学校数量）进行求和、求平均等操作。 ```python # 假设数据中有‘county_code’（区县代码）和‘population’（人口）字段 counties = zhejiang_towns.dissolve(by='county_code', aggfunc={ 'town_name': 'count', # 统计该区县下有多少个乡镇 'population': 'sum', # 将人口求和 'area_sq_km': 'sum', # 将面积求和 'school_count': 'sum' # 将学校数量求和 }).reset_index() # 重命名聚合后的列 counties = counties.rename(columns={'town_name': 'town_count'}) print(counties[['county_code', 'town_count', 'population', 'school_count']].head()) ``` `dissolve`操作非常强大，但要注意几何合并可能产生复杂多边形，有时需要后续的几何简化（`simplify`）来优化。 ## 4. 性能优化与实战工作流 当数据量达到全国乡镇级别时，性能就成为不可忽视的因素。除了之前提到的使用高效存储格式（GeoPackage/Parquet）外，还有以下优化技巧： - **使用空间索引**：`geopandas`在执行空间连接（`sjoin`）时，会自动利用`rtree`空间索引来加速查询。确保你的库版本支持此功能。 - **选择性读取列**：使用`GeoPackage`或`Parquet`时，可以通过`columns`参数只读取需要的属性字段，减少内存占用。 - **分块处理**：对于全国数据，可以按省或大区进行循环处理，每次只处理一小部分。 下面是一个将上述所有步骤串联起来的**实战工作流示例**，目标是为每个乡镇计算其到所属地级市市政府的距离： ```python import geopandas as gpd from shapely.geometry import Point # 1. 加载数据（假设已转换为GeoPackage） towns_gdf = gpd.read_file('national_townships.gpkg', columns=['code', 'name', 'city_code', 'geometry']) cities_gdf = gpd.read_file('city_centroids.gpkg') # 假设这是地级市市政府的点数据 # 2. 确保坐标系一致，并转换为投影坐标系以计算米制距离 target_crs = 'EPSG:4526' # 例如使用CGCS2000投影 towns_proj = towns_gdf.to_crs(target_crs) cities_proj = cities_gdf.to_crs(target_crs) # 3. 为每个乡镇找到其所属地级市的市政府点 # 首先通过属性连接（如果数据有city_code字段） merged = towns_proj.merge(cities_proj[['city_code', 'geometry']], left_on='city_code', right_on='city_code', suffixes=('_town', '_city')) # 重命名几何列 merged = merged.set_geometry('geometry_town') merged['city_center'] = merged['geometry_city'] # 4. 计算距离 merged['dist_to_city_center_m'] = merged.apply(lambda row: row.geometry_town.distance(row.city_center), axis=1) # 5. 保存结果 merged[['code', 'name', 'dist_to_city_center_m']].to_file('towns_with_distance.gpkg', driver='GPKG') print("处理完成，结果已保存。") ``` 这个工作流涵盖了数据读取、坐标转换、属性连接、空间计算和结果输出，是一个典型的端到端处理案例。在实际项目中，你可能还需要加入更多的数据验证和异常处理逻辑。 处理完数据，最后别忘了可视化检查一下。用`geopandas`的`.plot()`方法快速查看一下某个省的乡镇分布，或者用`contextily`叠加一个在线底图，可以直观地验证数据处理结果是否正确。比如，检查合并后的区县边界是否光滑连续，筛选出的“街道”是否确实位于城市核心区域。地图是检验空间数据处理质量最直接的“试金石”。