基于GeoPandas的空间数据分析实战：历史区划与现代经济数据融合

最近在技术社区看到一个很有意思的讨论：如果抛开现实行政约束，单纯从数据分析和地理信息系统的角度，去模拟一个“假如河南恢复清朝区划”的GDP变化，会得到什么结果？这听起来像是一个历史或地理话题，但对于我们开发者而言，其内核是一个 典型的多源数据融合、空间分析与经济指标重算的实战项目 。

很多开发者学了Python、Pandas、GeoPandas，也了解一些GIS（地理信息系统）概念，但往往停留在处理单一数据源或完成教程案例上。当面对“将历史行政区划映射到现代经济数据”这类复杂、跨域、数据脏乱的问题时，就不知从何下手。这背后涉及 历史地图矢量化、空间连接（Spatial Join）、数据聚合与归一化、缺失值处理 等一系列关键技术点。

本文将以“河南清朝区划GDP模拟”为引子，拆解一个完整的数据工程项目。我不会空谈理论，而是带你从零开始，用代码一步步解决：如何获取并清洗清朝地图数据？如何与当代的县级GDP数据做空间匹配？匹配不上怎么办？人均GDP又该如何科学地重新计算？最终，你将获得一套可复用的 空间数据分析方法论 和 完整的Python代码实现 ，能够处理类似的“古今对照”、“区域重组”数据分析需求。

1. 核心问题拆解：这不是历史题，而是数据工程题

首先，我们必须明确技术边界。我们不是在讨论行政区划调整的可行性或合理性，而是将其视为一个 纯粹的数据处理与模拟问题 。这能让我们聚焦于技术实现。

要实现“恢复清朝区划并计算GDP”，我们需要解决几个核心的数据挑战：

1. 空间数据对齐 ：清朝的“府”、“州”、“县”边界与今天的“市”、“县”、“区”边界绝大多数都不重合。一个清朝的府可能覆盖今天几个市的部分区域。我们需要将今天以区县为单位的GDP数据，“分配”到历史的行政区划单元上。
2. 数据粒度转换 ：我们拥有的现代经济数据（如GDP）通常以县级行政区为最小统计单元。而我们的目标单元是清朝的府级行政区。这涉及到数据的 聚合（Aggregation） 。
3. 权重分配难题 ：如何将某个现代区县的GDP“分配”给覆盖它的多个清朝府？按面积平分？这显然不合理，因为经济产出不是均匀分布的。我们需要更合理的代理指标，例如 夜间灯光数据 或 人口分布数据 作为权重。
4. 人均GDP计算 ：人均GDP = GDP / 人口。清朝的人口数据与当代人口数据天差地别，直接使用当代人口除以历史区划没有意义。因此，这个模拟通常 只模拟GDP总量在不同历史区划下的重新排列 ，而“人均GDP”则需要基于模拟后的GDP总量和 某种估算的历史人口 来计算，这引入了更大的不确定性。在本文中，我们将重点放在GDP总量的空间重算上。

所以，本文的真正目标是： 教会你如何使用GeoPandas等工具，完成跨时空、非对齐边界下的经济数据空间重分配，并输出可视化结果。

2. 技术栈与核心概念

在开始编码前，快速理解几个关键概念和我们将要使用的工具：

• GeoPandas ：Python中处理地理空间数据的“神级”库。它扩展了Pandas，使得DataFrame可以存储几何列（点、线、面），并能进行空间运算。它是本项目的核心。
• Shapely ：用于操作和分析平面几何对象的库。GeoPandas的几何运算依赖于它。
• Fiona ：用于读写空间数据文件（如Shapefile）的库，通常由GeoPandas间接调用。
• 空间连接（Spatial Join） ：GIS中的核心操作。它根据两个空间数据集的地理位置关系（如相交、包含、Within）来连接它们的属性表。在本项目中，我们将用现代区县面（Polygon）去连接清朝府的面。
• Shapefile ：一种常见的空间矢量数据格式，由 .shp （几何）、 .dbf （属性）、 .shx （索引）等文件组成。
• CRS（坐标参考系统） ：定义几何对象如何与地球表面关联的坐标系。在进行空间运算前， 必须确保所有数据层处于同一CRS ，否则结果毫无意义。常用的是 EPSG:4326 （WGS84，经纬度）和 EPSG:3857 （Web墨卡托投影）。

工作流程概览 ：

1. 获取并加载清朝河南行政区划矢量数据（面数据）。
2. 获取并加载当代河南区县级行政区划矢量数据及附带的GDP属性数据。
3. 统一两者的坐标参考系统（CRS）。
4. 进行空间连接，找出每个当代区县与哪些清朝府相交。
5. 设计算法，将区县的GDP按一定权重（如相交面积占区县面积的比例）分配给它所在的各个清朝府。
6. 按清朝府聚合分配后的GDP，得到模拟结果。
7. 可视化对比。

3. 环境准备与数据获取

3.1 Python环境与库安装

建议使用Conda或venv创建独立的Python环境。

# 使用pip安装核心库
pip install geopandas pandas numpy matplotlib contextily jupyterlab

• geopandas : 核心空间数据处理库。
• pandas , numpy : 数据处理基础。
• matplotlib : 绘图。
• contextily : 为地图添加底图。
• jupyterlab : 推荐在Notebook中交互式运行代码，便于调试和可视化。

3.2 关键数据获取

这是项目的难点和起点。数据质量直接决定结果的可信度。

1. 清朝河南行政区划矢量数据（ qing_henan.shp ）

• 来源 ：学术机构或历史GIS项目。例如，复旦大学历史地理研究中心发布的“中国历史地理信息系统（CHGIS）”数据可能包含清代矢量数据。本项目假设你已通过合法渠道获得了一份简化版的清朝河南府级Shapefile数据。
• 数据内容 ：每个面要素代表一个清朝的“府”或“直隶州”，属性表至少包含 name （府名）字段。

2. 当代河南区县矢量数据与GDP数据（ modern_henan.shp 及 gdp_data.csv ）

• 矢量数据 ：可以从资源环境科学与数据中心、GADM等网站获取中国县级行政区划Shapefile。属性表应包含 code （行政区划代码）、 name （区县名）、 city （所属地市）等字段。
• GDP数据 ：需要从《河南统计年鉴》或相关统计网站收集，整理成CSV格式，至少包含 code （与矢量数据对应）、 gdp_total （GDP总量，单位：亿元）字段。 注意：务必使用同一年份的数据以保证一致性。

假设的数据文件结构 ：

你的项目目录/
├── data/
│   ├── qing_henan/          # 清朝数据文件夹
│   │   ├── qing_henan.shp
│   │   ├── qing_henan.dbf
│   │   └── ...
│   ├── modern_henan/        # 现代数据文件夹
│   │   ├── modern_henan.shp
│   │   ├── modern_henan.dbf
│   │   └── ...
│   └── gdp_2023.csv         # 2023年河南区县GDP数据
└── qing_gdp_simulation.ipynb # 你的Jupyter Notebook

4. 核心流程拆解与代码实现

让我们开始真正的代码实战。以下代码块需要按顺序在Jupyter Notebook中执行。

4.1 导入库与加载数据

import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import contextily as ctx
from shapely.geometry import Polygon, MultiPolygon
# 确保绘图在Notebook内显示
%matplotlib inline

# 1. 加载清朝河南行政区划数据
qing_gdf = gpd.read_file('./data/qing_henan/qing_henan.shp')
print("清朝数据概览：")
print(qing_gdf.info())
print(f"清朝府级数量：{len(qing_gdf)}")
print(qing_gdf.head())

# 2. 加载当代河南区县行政区划数据
modern_gdf = gpd.read_file('./data/modern_henan/modern_henan.shp')
print("\n当代数据概览：")
print(modern_gdf.info())
print(f"当代区县数量：{len(modern_gdf)}")
print(modern_gdf.head())

# 3. 加载当代区县GDP数据
gdp_df = pd.read_csv('./data/gdp_2023.csv', dtype={'code': str}) # code列读为字符串
print("\nGDP数据概览：")
print(gdp_df.info())
print(gdp_df.head())

关键检查点 ：

• 查看 qing_gdf.crs 和 modern_gdf.crs ，确认它们的坐标参考系统。如果不同，必须转换。
• 检查 modern_gdf 和 gdp_df 能否通过 code 字段正确关联。

4.2 数据清洗与预处理

# 1. 统一CRS (假设现代数据是EPSG:4326，清朝数据也需要转换到此)
target_crs = 'EPSG:4326'
if qing_gdf.crs != target_crs:
    qing_gdf = qing_gdf.to_crs(target_crs)
if modern_gdf.crs != target_crs:
    modern_gdf = modern_gdf.to_crs(target_crs)
print(f"清朝数据CRS：{qing_gdf.crs}")
print(f"当代数据CRS：{modern_gdf.crs}")

# 2. 将GDP数据合并到当代区县GeoDataFrame中
# 确保modern_gdf中也有字符串类型的code字段用于合并
modern_gdf['code'] = modern_gdf['code'].astype(str)
modern_gdf_with_gdp = modern_gdf.merge(gdp_df, on='code', how='left')
# 检查是否有区县缺失GDP数据
missing_gdp = modern_gdf_with_gdp['gdp_total'].isnull().sum()
print(f"缺失GDP数据的区县数量：{missing_gdp}")
if missing_gdp > 0:
    print(modern_gdf_with_gdp[modern_gdf_with_gdp['gdp_total'].isnull()][['name', 'code']])
    # 处理缺失值：这里简单用0填充，实际项目需根据情况处理（如用上一级市均值、或剔除）
    modern_gdf_with_gdp['gdp_total'] = modern_gdf_with_gdp['gdp_total'].fillna(0)

# 3. 计算每个当代区县的面积（平方公里）
# 注意：EPSG:4326是经纬度，计算面积需转换为投影坐标系（如EPSG:3857或等面积投影）
modern_gdf_with_gdp['area_km2'] = modern_gdf_with_gdp.to_crs('EPSG:3857').geometry.area / 10**6
print(f"当代区县平均面积：{modern_gdf_with_gdp['area_km2'].mean():.2f} km²")

4.3 空间连接与GDP分配（核心算法）

这是最核心的一步。我们采用“面积权重法”作为基础模型：假设一个区县内的GDP在其辖区内是均匀分布的，那么分配给某个相交清朝府的GDP，比例等于相交面积占该区县总面积的比例。

# 1. 执行空间连接，找出每个当代区县与哪些清朝府相交
# 使用`how='inner'`只保留有相交的部分
intersection_gdf = gpd.overlay(modern_gdf_with_gdp, qing_gdf, how='intersection')
print(f"空间连接后生成的相交单元数量：{len(intersection_gdf)}")
print(intersection_gdf.head())

# 2. 计算每个相交单元的面积
intersection_gdf['intersection_area_km2'] = intersection_gdf.to_crs('EPSG:3857').geometry.area / 10**6

# 3. 将每个相交单元的面积除以其所属原始区县的面积，得到面积比例
# 首先需要将区县面积信息合并到相交单元中
area_lookup = modern_gdf_with_gdp[['code', 'area_km2']].set_index('code')
intersection_gdf['parent_area_km2'] = intersection_gdf['code'].map(area_lookup['area_km2'])
intersection_gdf['area_ratio'] = intersection_gdf['intersection_area_km2'] / intersection_gdf['parent_area_km2']

# 4. 根据面积比例，分配GDP
intersection_gdf['allocated_gdp'] = intersection_gdf['gdp_total'] * intersection_gdf['area_ratio']

# 检查分配是否合理：一个区县分配出去的总GDP应约等于其原始GDP（由于计算误差和非常小的碎片，可能略有出入）
check_allocation = intersection_gdf.groupby('code')['allocated_gdp'].sum()
original_gdp = modern_gdf_with_gdp.set_index('code')['gdp_total']
comparison = pd.DataFrame({'allocated': check_allocation, 'original': original_gdp})
comparison['diff_pct'] = (comparison['allocated'] - comparison['original']) / comparison['original'] * 100
print("\nGDP分配一致性检查（前5个区县）：")
print(comparison.head())
print(f"\n分配误差绝对值平均：{comparison['diff_pct'].abs().mean():.4f}%")

4.4 按清朝府聚合GDP并分析结果

# 1. 按清朝府聚合分配到的GDP
qing_gdp_result = intersection_gdf.groupby('name_2')['allocated_gdp'].sum().reset_index() # 假设清朝府名列在‘name_2’
qing_gdp_result.rename(columns={'name_2': 'qing_fu', 'allocated_gdp': 'simulated_gdp_total'}, inplace=True)
qing_gdp_result['simulated_gdp_total'] = qing_gdp_result['simulated_gdp_total'].round(2) # 保留两位小数

# 2. 排序查看结果
qing_gdp_result_sorted = qing_gdp_result.sort_values('simulated_gdp_total', ascending=False)
print("模拟的清朝河南各府GDP总量排名（前10）：")
print(qing_gdp_result_sorted.head(10))

# 3. 计算总量并与当代河南全省GDP对比
simulated_total_gdp = qing_gdp_result['simulated_gdp_total'].sum()
modern_total_gdp = modern_gdf_with_gdp['gdp_total'].sum()
print(f"\n模拟的清朝区划下河南总GDP：{simulated_total_gdp:.2f} 亿元")
print(f"当代区划下河南总GDP：{modern_total_gdp:.2f} 亿元")
print(f"两者差异（模拟-当代）：{simulated_total_gdp - modern_total_gdp:.2f} 亿元")
print(f"差异百分比：{(simulated_total_gdp - modern_total_gdp)/modern_total_gdp*100:.4f}%")
# 理论上，由于我们只是重新分配，总量应该几乎相等。微小差异来源于计算误差和可能的几何碎片。

4.5 结果可视化

# 1. 将模拟结果合并回清朝GeoDataFrame用于绘图
qing_gdf_with_result = qing_gdf.merge(qing_gdp_result, left_on='name', right_on='qing_fu', how='left')
qing_gdf_with_result['simulated_gdp_total'] = qing_gdf_with_result['simulated_gdp_total'].fillna(0) # 处理可能没有匹配到的府

# 2. 绘制模拟GDP分布图
fig, ax = plt.subplots(1, 2, figsize=(20, 8))

# 子图1：模拟GDP分级色彩图
qing_gdf_with_result.plot(column='simulated_gdp_total',
                          ax=ax[0],
                          legend=True,
                          cmap='OrRd', # 橙红色系
                          legend_kwds={'label': "模拟GDP总量（亿元）", 'orientation': "horizontal"},
                          edgecolor='black',
                          linewidth=0.3)
ax[0].set_title('模拟：清朝河南区划下的GDP分布')
ctx.add_basemap(ax[0], crs=qing_gdf_with_result.crs.to_string(), source=ctx.providers.CartoDB.Positron) # 添加底图

# 子图2：现代河南区县GDP分布图（作为对比）
modern_gdf_with_gdp.plot(column='gdp_total',
                         ax=ax[1],
                         legend=True,
                         cmap='OrRd',
                         legend_kwds={'label': "实际GDP总量（亿元）", 'orientation': "horizontal"},
                         edgecolor='black',
                         linewidth=0.2)
ax[1].set_title('实际：当代河南区划下的GDP分布')
ctx.add_basemap(ax[1], crs=modern_gdf_with_gdp.crs.to_string(), source=ctx.providers.CartoDB.Positron)

plt.tight_layout()
plt.show()

# 3. 绘制GDP排名柱状图
plt.figure(figsize=(12, 6))
top_n = 15
top_data = qing_gdp_result_sorted.head(top_n)
bars = plt.barh(range(top_n), top_data['simulated_gdp_total'].values[::-1]) # 反转使最高的在上方
plt.yticks(range(top_n), top_data['qing_fu'].values[::-1])
plt.xlabel('模拟GDP总量（亿元）')
plt.title(f'模拟GDP总量排名（前{top_n}）')
# 在柱子上添加数值
for i, (bar, v) in enumerate(zip(bars, top_data['simulated_gdp_total'].values[::-1])):
    plt.text(v + max(top_data['simulated_gdp_total'])*0.01, bar.get_y() + bar.get_height()/2,
             f'{v:.0f}', va='center')
plt.tight_layout()
plt.show()

5. 运行结果与解读

运行上述代码后，你将得到：

1. 控制台输出 ：包含数据基本信息、GDP分配误差检查、模拟GDP排名及总量对比。理想情况下，模拟总量与实际总量差异应小于0.1%，这验证了分配算法的基本一致性。
2. 两张对比地图 ：左图为模拟的清朝区划下GDP分布，右图为实际的当代区划下GDP分布。你可以直观地看到经济重心（如郑州、洛阳周边）在两种区划下的形态差异。例如，“开封府”可能因包含当代郑州部分区域而GDP大增，“河南府”（洛阳）的格局也可能发生变化。
3. 一张排名柱状图 ：清晰展示在模拟情境下，哪些清朝府会成为经济重镇。

核心发现示例（基于模拟结果，非真实结论） ：

• “开封府”可能成为绝对核心 ：因为其历史辖区覆盖了当代郑州、开封等经济强市的大部分区域。
• “河南府”（洛阳）地位依旧突出 。
• 一些边缘府的经济总量被稀释 ，因为现代经济高度集中在少数城市群。
• 总量不变，但分布格局重塑 ：这直观地展示了行政区划如何影响我们对区域经济实力的“观感”。

6. 常见问题与排查思路

7. 模型优化与最佳实践

基础的“面积权重法”模型假设经济密度均匀，这显然粗糙。在实际研究中，我们可以引入更精细的权重：

1. 夜间灯光数据权重 ：使用DMSP/OLS或VIIRS夜间灯光数据作为经济活动的代理指标。将灯光亮度值聚合到区县，然后按“相交区域内灯光总值占区县灯光总值的比例”来分配GDP。这比单纯按面积分配更合理。

   # 伪代码思路
   # 假设已有每个区县的灯光总值 `light_total`
   # 计算相交区域的灯光值（需要灯光栅格数据与多边形做zonal statistics）
   intersection_gdf['light_ratio'] = intersection_gdf['intersection_light'] / intersection_gdf['parent_light_total']
   intersection_gdf['allocated_gdp'] = intersection_gdf['gdp_total'] * intersection_gdf['light_ratio']
   

2. 人口分布权重 ：如果有所需年份的网格化人口数据，原理同上，按人口比例分配。

3. 土地利用类型权重 ：将建设用地、耕地的经济产出系数设高，林地、水域设低，进行加权分配。

工程实践建议 ：

• 版本控制 ：将数据处理流程（如 .ipynb 或 .py 脚本）和关键中间数据（如清洗后的CSV）纳入Git管理。
• 模块化 ：将数据加载、清洗、空间连接、分配计算、可视化分别写成函数，提高代码可读性和复用性。
• 参数化 ：将关键参数（如GDP年份、权重方法、面积过滤阈值）放在脚本开头，便于调整实验。
• 数据验证 ：在每一步关键操作后，都进行数据完整性检查（如非空检查、总和校验）。
• 文档化 ：在代码中清晰注释数据来源、假设、算法选择理由。

8. 总结与拓展方向

通过这个项目，我们完成了一次完整的 空间数据分析和经济指标重算 的实战。核心收获不在于“河南GDP会怎么变”这个具体答案，而在于掌握了一套处理 时空非对齐数据 的方法论：

1. 空间思维 ：将抽象的经济数据与具体的地理空间绑定。
2. 工具链 ：熟练使用GeoPandas进行数据操作、空间连接和可视化。
3. 算法设计 ：理解并实现了基于空间关系的指标分配算法。
4. 数据批判性思维 ：认识到模型假设（如均匀分布）的局限性，并知道如何引入更合理的代理变量进行优化。

你可以将这套方法轻松迁移到其他有趣的问题上 ：

• 历史分析 ：模拟唐朝“道”、元朝“行省”下的现代经济格局。
• 规划评估 ：评估某个新规划的行政区（如新区、都市圈）对现有经济统计指标的影响。
• 资源分配 ：将医院、学校等点状公共服务设施的服务能力，按人口权重分配到各个居住小区。
• 环境研究 ：计算不同历史时期湖泊、森林的面积变化，并关联气候变化数据。

记住，技术是骨架，数据是血肉，而提出一个好问题才是灵魂。希望本文提供的代码和思路，能成为你探索更多空间数据奥秘的起点。建议收藏本文，在遇到类似的多源数据融合与空间分析任务时，随时回来查阅这套标准化的解决流程。