2024行政区划数据处理实战：从SHP文件到洞察地图的Python全流程

最近在做一个区域商业分析的项目，客户扔过来一个最新的行政区划SHP文件包，说是2024年更新的数据。打开一看，里面包含了省、市、县三级的驻地点位信息，坐标是WGS1984。本以为用熟悉的工具简单处理一下就能出图，结果在实际操作中遇到了编码问题、几何校验失败、属性表字段混乱等一系列“坑”。这让我意识到，处理看似标准的GIS数据，远不是加载-显示那么简单，尤其是当数据时效性要求高、分析维度复杂时。

这篇文章就是基于那次实战经历整理出来的。如果你也是数据分析师、GIS开发者，或者任何需要处理最新行政区划点位数据的朋友，希望这份从数据清洗到可视化的完整Python方案能帮你少走弯路。我们将完全使用geopandas及相关生态库，专注于解决实际问题，而不是泛泛而谈库的功能。

1. 环境搭建与数据初探：避开第一个坑

工欲善其事，必先利其器。处理SHP文件，一个稳定且版本匹配的Python环境至关重要。我强烈建议使用Conda来管理你的GIS Python环境，因为它能很好地处理geopandas背后复杂的C库依赖（比如GDAL、Fiona）。

# 创建一个新的conda环境
conda create -n gis_analysis python=3.10
conda activate gis_analysis

# 安装核心地理数据处理库
conda install -c conda-forge geopandas
# 安装可视化相关库
conda install -c conda-forge matplotlib contextily folium
# 安装数据处理辅助库
conda install -c conda-forge pyproj rtree

注意：务必通过conda-forge频道安装geopandas。用pip直接安装geopandas极易因底层GDAL库版本冲突而失败，这是新手最常见的“开局雷击”。

环境准备好后，让我们加载数据。假设你的数据包解压后包含以下文件：

• province_points.shp (省驻地点位)
• city_points.shp (市驻地点位)
• county_points.shp (县驻地点位)

SHP格式实际由多个文件组成，.shp存储几何图形，.dbf存储属性数据，.prj存储坐标系统信息。用geopandas读取时，只需指定.shp文件路径即可。

import geopandas as gpd
import pandas as pd
import os

# 定义数据路径
data_dir = './2024_data/'
province_path = os.path.join(data_dir, 'province_points.shp')
city_path = os.path.join(data_dir, 'city_points.shp')
county_path = os.path.join(data_dir, 'county_points.shp')

# 读取数据
try:
    gdf_province = gpd.read_file(province_path, encoding='utf-8')
    gdf_city = gpd.read_file(city_path, encoding='utf-8')
    gdf_county = gpd.read_file(county_path, encoding='utf-8')
except UnicodeDecodeError:
    # 如果utf-8失败，尝试GBK或GB18030，这是中文GIS数据常见的编码
    gdf_province = gpd.read_file(province_path, encoding='gb18030')
    gdf_city = gpd.read_file(city_path, encoding='gb18030')
    gdf_county = gpd.read_file(county_path, encoding='gb18030')

print(f"省级数据记录数: {len(gdf_province)}")
print(f"市级数据记录数: {len(gdf_city)}")
print(f"县级数据记录数: {len(gdf_county)}")

初次读取后，别急着画图。先花几分钟做一次数据“体检”，查看数据结构、坐标系和基本信息：

# 查看省级数据的属性表前几行
print(gdf_province.head())
# 查看数据结构信息
print(gdf_province.info())
# 查看坐标系
print(gdf_province.crs)

一个健康的GeoDataFrame应该包含一个geometry列（存储点位几何信息）和若干属性列（如名称、代码、统计指标等）。crs属性如果显示EPSG:4326，则代表是WGS1984地理坐标系，这与我们输入信息一致。如果不是，后续的坐标转换步骤就必不可少。

2. 数据清洗与质量校验：构建可靠的分析基底

原始数据，尤其是从不同来源整合的最新数据，几乎总存在一些瑕疵。直接用于分析可能导致结果偏差或程序报错。我们的清洗工作主要围绕几何有效性、属性完整性和空间一致性展开。

首先处理几何问题。无效的几何图形（例如自相交、空几何）会导致许多空间操作失败。

# 检查是否存在无效几何
invalid_province = ~gdf_province.is_valid
if invalid_province.any():
    print(f"省级数据中发现 {invalid_province.sum()} 个无效几何体。")
    # 尝试修复无效几何（缓冲区为0是一个常用技巧）
    gdf_province.loc[invalid_province, 'geometry'] = gdf_province.loc[invalid_province, 'geometry'].buffer(0)

# 对市、县级数据执行同样的检查与修复

其次是属性数据清洗。常见问题包括字段名不统一、重要字段缺失、存在空值或异常值。

# 1. 统一关键字段名（假设我们需要‘名称’和‘行政代码’）
# 先查看所有列名，发现原始数据可能用‘NAME’, ‘名称’，‘xingzhengqu’等
print(gdf_city.columns.tolist())

# 假设我们发现市级数据中，名称字段为‘CITY_NAME’，代码字段为‘CODE’
# 我们将其标准化
gdf_city = gdf_city.rename(columns={'CITY_NAME': 'name', 'CODE': 'code'})

# 2. 处理缺失值
# 检查‘name’字段是否有空
missing_name = gdf_city['name'].isnull().sum()
if missing_name > 0:
    print(f"市级数据有 {missing_name} 条记录名称缺失。")
    # 根据实际情况处理：如果‘code’存在，可以用代码映射表补全；否则考虑删除或标记
    # 例如，删除名称缺失的记录（谨慎操作，需结合业务）
    # gdf_city = gdf_city.dropna(subset=['name'])

# 3. 去除完全重复的记录（几何和属性都相同）
gdf_county = gdf_county.drop_duplicates(subset=['geometry', 'name', 'code'], keep='first')

最后是空间参考系的确认与统一。虽然输入说明是WGS1984，但有时.prj文件可能丢失或错误。我们需要确保所有图层在同一个坐标系下进行分析。

# 检查并统一坐标系
target_crs = 'EPSG:4326' # WGS84

if gdf_province.crs is None:
    print("省级数据未定义坐标系，将设置为WGS84。")
    gdf_province.set_crs(target_crs, inplace=True)
elif gdf_province.crs != target_crs:
    print(f"省级数据坐标系为 {gdf_province.crs}，将转换至WGS84。")
    gdf_province = gdf_province.to_crs(target_crs)

# 对市、县级数据执行同样的操作
# 确保所有图层坐标系一致后，才能进行叠加分析或联合绘图

完成这些步骤，我们就得到了一个相对“干净”的数据集，为后续的深度分析和可视化打下了坚实基础。

3. 空间分析与关系构建：让数据产生连接

单纯的点位展示价值有限。将不同层级的行政区划数据关联起来，才能挖掘出更深层次的洞察。例如，我们可以计算每个省下有多少个市、每个市下有多少个县，或者计算县点到所属市点的距离。

由于我们目前只有点位数据，缺乏面状的行政区划边界，无法进行严格的空间包含分析。但我们可以利用行政代码（如国家标准行政区划代码）来建立层级关联。假设我们的数据中含有类似code（本级代码）和pcode（上级代码）的字段。

# 示例：建立省-市关联
# 假设省级数据有‘code’（省代码），市级数据有‘code’（市代码）和‘pcode’（所属省代码）
# 为每个市级点位找到所属的省级点位信息

# 首先将省级数据的代码和名称提取为字典，便于映射
province_dict = gdf_province.set_index('code')[['name']].to_dict()['name']

# 在市级数据中创建新列‘province_name’，通过‘pcode’从字典映射获取
gdf_city['province_name'] = gdf_city['pcode'].map(province_dict)

# 检查映射结果
print(gdf_city[['name', 'pcode', 'province_name']].head())

# 统计每个省有多少个市
city_count_by_province = gdf_city.groupby('province_name').size().sort_values(ascending=False)
print("各省下辖市数量统计（前10）：")
print(city_count_by_province.head(10))

如果数据中没有明确的层级代码字段，但有点位的经纬度，且我们知道行政区划在空间上是嵌套的，我们可以通过空间连接（Spatial Join）的近似方法来实现。但这需要面状的边界数据。这里引出一个高级技巧：如何利用公开的边界数据增强你的点位数据？

你可以从一些开源GIS数据平台获取到2024年（或近年）的省、市、县面状边界SHP文件。然后使用gpd.sjoin进行空间连接：

# 假设我们已经加载了省级面状边界GeoDataFrame: gdf_province_polygon
# 它包含几何列（多边形）和‘name’等属性

# 进行空间连接，将市级点位匹配到其所在的省面内
gdf_city_with_province = gpd.sjoin(gdf_city, gdf_province_polygon[['name', 'geometry']],
                                   how='left', predicate='within')
# ‘predicate=within’表示寻找完全在省面内的市点
# 结果中会多出一列‘name_right’，即所属省名称

除了层级关联，基础的空间度量也很有用，比如计算所有县点到其所属市点的球面距离（因为我们是地理坐标）：

from shapely.geometry import Point
import numpy as np

# 假设我们已经通过某种方式为每个县点匹配到了所属市点的几何对象‘city_geometry’
# 计算距离（单位：公里）
def calculate_distance(row):
    county_point = row['geometry']
    city_point = row['city_geometry']
    if pd.isna(city_point):
        return np.nan
    # 将Shapely几何对象转换为（经度，纬度）
    county_coords = (county_point.x, county_point.y)
    city_coords = (city_point.x, city_point.y)
    # 使用haversine公式计算大圆距离
    # 这里省略具体实现，可使用geopy库的`geodesic`函数更便捷准确
    # from geopy.distance import geodesic
    # distance_km = geodesic(county_coords, city_coords).km
    # return distance_km
    return None  # 此处为占位

# gdf_county['distance_to_city_km'] = gdf_county.apply(calculate_distance, axis=1)

这些分析结果可以作为新的属性字段添加到数据中，极大地丰富了数据的维度，为后续的可视化和统计分析提供了更多素材。

4. 多维度可视化：从静态地图到交互式探索

数据清洗和分析之后，可视化是将成果呈现给他人或自己探索的关键一步。我们将介绍三种不同侧重点的可视化方法：静态专题图、层级叠加图和交互式网页地图。

首先是静态专题图。我们可以用matplotlib配合geopandas快速绘制一张展示省级点位并标注名称的底图，同时用颜色或大小区分某些属性（假设我们有各省的‘GDP’字段）。

import matplotlib.pyplot as plt
import contextily as ctx

fig, ax = plt.subplots(1, 1, figsize=(15, 10))

# 绘制省级点位，点的大小和颜色可以映射到某个数值字段
# 假设gdf_province有‘gdp’字段
if 'gdp' in gdf_province.columns:
    # 归一化GDP用于点的大小
    norm_gdp = (gdf_province['gdp'] - gdf_province['gdp'].min()) / (gdf_province['gdp'].max() - gdf_province['gdp'].min())
    point_sizes = 50 + norm_gdp * 300  # 基础大小50，最大350
    scatter = gdf_province.plot(ax=ax, markersize=point_sizes, color='red', edgecolor='black', alpha=0.7, label='Province')
else:
    gdf_province.plot(ax=ax, markersize=50, color='red', edgecolor='black', label='Province')

# 添加标注（避免重叠，可以只标注重要的或使用工具优化）
for idx, row in gdf_province.iterrows():
    ax.annotate(text=row['name'], xy=(row.geometry.x, row.geometry.y),
                xytext=(3, 3), textcoords="offset points", fontsize=8, alpha=0.8)

# 添加底图（需要网络连接）
try:
    ctx.add_basemap(ax, crs=gdf_province.crs.to_string(), source=ctx.providers.CartoDB.Positron)
except:
    print("无法加载在线底图，将使用默认背景。")

ax.set_axis_off()
ax.set_title('2024年中国省级行政中心分布图', fontsize=16)
plt.legend()
plt.tight_layout()
# plt.savefig('province_points_map.png', dpi=300, bbox_inches='tight')
plt.show()

其次是层级叠加图。在一张图上同时展示省、市、县三级点位，并用不同的样式区分，可以直观感受行政层级和空间分布密度。

fig, ax = plt.subplots(1, 1, figsize=(18, 12))

# 绘制县级点位（最多，用最小最浅的点）
gdf_county.plot(ax=ax, markersize=2, color='grey', alpha=0.5, label='County (县)')
# 绘制市级点位
gdf_city.plot(ax=ax, markersize=15, color='blue', edgecolor='white', linewidth=0.5, label='City (市)')
# 绘制省级点位
gdf_province.plot(ax=ax, markersize=80, color='red', edgecolor='yellow', linewidth=1.5, label='Province (省)')

# 添加简洁的底图
try:
    ctx.add_basemap(ax, crs=gdf_province.crs.to_string(), source=ctx.providers.Stamen.TonerLite, alpha=0.6)
except:
    pass

ax.set_axis_off()
ax.set_title('2024年中国省-市-县三级行政中心分布', fontsize=18)
ax.legend(loc='upper left', fontsize=10)
plt.tight_layout()
# plt.savefig('hierarchy_points_map.png', dpi=300)
plt.show()

最后是交互式网页地图。使用folium库，我们可以生成一个HTML文件，在浏览器中实现缩放、点击查看属性等交互功能，非常适合汇报或探索。

import folium

# 以全国几何中心（近似）初始化地图
center_lat = gdf_province.geometry.y.mean()
center_lon = gdf_province.geometry.x.mean()
m = folium.Map(location=[center_lat, center_lon], zoom_start=4, tiles='CartoDB positron')

# 添加省级点位图层
for idx, row in gdf_province.iterrows():
    popup_text = f"<b>{row['name']}</b><br>Code: {row.get('code', 'N/A')}"
    folium.CircleMarker(
        location=[row.geometry.y, row.geometry.x],
        radius=6,
        color='red',
        fill=True,
        fillColor='red',
        fillOpacity=0.7,
        popup=folium.Popup(popup_text, max_width=250)
    ).add_to(m)

# 添加市级点位图层（可以控制显示层级，避免过于密集）
# 使用FeatureGroup来管理，方便控制显示/隐藏
city_fg = folium.FeatureGroup(name='Cities', show=False)
for idx, row in gdf_city.iterrows():
    popup_text = f"<b>{row['name']}</b><br>Province: {row.get('province_name', 'N/A')}"
    folium.CircleMarker(
        location=[row.geometry.y, row.geometry.x],
        radius=3,
        color='blue',
        fill=True,
        fillColor='blue',
        fillOpacity=0.5,
        popup=folium.Popup(popup_text, max_width=250)
    ).add_to(city_fg)
city_fg.add_to(m)

# 添加图层控制
folium.LayerControl().add_to(m)

# 保存为HTML文件
m.save('2024_admin_points_interactive_map.html')
print("交互式地图已生成，请在浏览器中打开 '2024_admin_points_interactive_map.html' 查看。")

这三种可视化方式各有优劣：静态图适合放入报告或论文；叠加图便于整体观察层级关系；交互式地图则提供了最强的探索性。你可以根据项目需求选择或组合使用。

5. 实战案例：构建区域经济密度热点图

让我们把前面所有的技巧串联起来，完成一个稍微复杂的实战任务：假设我们除了行政区划点位，还获得了一份附加了2024年预估GDP数据的属性表。我们的目标是生成一张区域经济密度热点图，这里我们用“省级单位GDP”和“市级点位空间密度”来模拟“经济热度”。

第一步，数据融合。将经济数据与空间点位数据关联。

# 假设我们有一个CSV文件‘2024_gdp_estimate.csv’，包含‘code’和‘gdp_estimate’两列
gdp_df = pd.read_csv('2024_gdp_estimate.csv')
# 将GDP数据合并到省级GeoDataFrame
gdf_province = gdf_province.merge(gdp_df, on='code', how='left')
print(gdf_province[['name', 'gdp_estimate']].head())

第二步，计算空间密度。我们使用核密度估计（Kernel Density Estimation, KDE）来可视化市级点位的聚集程度，这可以间接反映城市分布的密集区，通常也是经济活跃区。

import numpy as np
from scipy.stats import gaussian_kde

# 获取所有市级点位的坐标
city_coords = np.vstack([gdf_city.geometry.x, gdf_city.geometry.y])

# 计算核密度
kde = gaussian_kde(city_coords)
# 创建一个覆盖全国的网格
xmin, ymin, xmax, ymax = gdf_city.total_bounds
xi, yi = np.mgrid[xmin:xmax:100j, ymin:ymax:100j] # 100x100的网格
coords = np.vstack([xi.flatten(), yi.flatten()])
zi = kde(coords).reshape(xi.shape)

第三步，分层设色与叠加绘图。将省级GDP（用点的大小或颜色表示）与市级密度（用背景色图表示）叠加在一张图上。

fig, ax = plt.subplots(1, 1, figsize=(20, 15))

# 1. 绘制市级点位密度图（背景）
density_plot = ax.pcolormesh(xi, yi, zi, shading='auto', cmap='YlOrRd', alpha=0.6)
plt.colorbar(density_plot, ax=ax, label='City Point Density (KDE)')

# 2. 绘制省级点位，大小映射到GDP
if 'gdp_estimate' in gdf_province.columns:
    # 归一化GDP，用于点的大小和颜色
    gdp_norm = (gdf_province['gdp_estimate'] - gdf_province['gdp_estimate'].min()) / (gdf_province['gdp_estimate'].max() - gdf_province['gdp_estimate'].min())
    scatter_size = 100 + gdp_norm * 500
    scatter_color = plt.cm.viridis(gdp_norm) # 使用viridis配色表示GDP高低

    scatter = ax.scatter(gdf_province.geometry.x, gdf_province.geometry.y,
                         s=scatter_size, c=scatter_color, edgecolors='black', linewidth=1.5, zorder=5)
    # 为散点图添加图例（需要手动创建代理对象）
    from matplotlib.lines import Line2D
    legend_elements = [Line2D([0], [0], marker='o', color='w', label='High GDP',
                              markerfacecolor=plt.cm.viridis(0.9), markersize=15),
                       Line2D([0], [0], marker='o', color='w', label='Low GDP',
                              markerfacecolor=plt.cm.viridis(0.1), markersize=15)]
    ax.legend(handles=legend_elements, loc='upper right')
else:
    gdf_province.plot(ax=ax, markersize=100, color='red', edgecolor='black', zorder=5)

# 3. 添加省级标注
for idx, row in gdf_province.iterrows():
    ax.annotate(text=row['name'], xy=(row.geometry.x, row.geometry.y),
                xytext=(5, 5), textcoords="offset points", fontsize=9, weight='bold',
                bbox=dict(boxstyle="round,pad=0.3", facecolor="white", alpha=0.7, edgecolor='none'))

# 4. 添加简洁底图
try:
    ctx.add_basemap(ax, crs=gdf_province.crs.to_string(), source=ctx.providers.Stamen.TonerBackground, alpha=0.3)
except:
    pass

ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
ax.set_axis_off()
ax.set_title('2024年区域经济热度模拟图\n(背景：城市分布密度 | 气泡：省级经济规模)', fontsize=18, pad=20)
plt.tight_layout()
# plt.savefig('economic_hotspot_map.png', dpi=300, bbox_inches='tight')
plt.show()

这张合成图虽然是一个模拟，但它清晰地展示了如何将空间分布特征（密度）与属性特征（经济指标）结合起来，生成更有洞察力的可视化成果。在实际项目中，你可以替换为真实的人口密度、灯光指数、企业POI密度等数据，与行政区划点位进行叠加分析。

处理最新的行政区划数据，核心在于流程的严谨性和思维的灵活性。从数据读取的编码问题，到空间关系的构建，再到最终的可视化表达，每一步都可能遇到意想不到的情况。我自己的经验是，多写检查代码，多用print或logging输出中间结果，确保每一步都符合预期。当遇到坐标系转换、空间连接失败时，回头检查几何有效性往往是解决问题的关键。最后，别忘了可视化不仅是结果的展示，更是探索数据、发现新问题的过程。尝试不同的地图样式、叠加不同的数据层，或许会有意想不到的发现。