Python 实现两步移动搜索法：基于 ArcPy 与 Pandas 的公共服务可达性计算
在城市规划与地理信息科学领域，公共服务设施的可达性分析是评估资源分配公平性的重要工具。两步移动搜索法（2SFCA）作为经典的空间可达性计算方法，能够量化居民获取公共服务的便利程度。本文将手把手带您用Python实现完整的2SFCA分析流程，结合ArcPy处理空间数据与Pandas进行高效计算，最终生成可直接用于决策的可视化结果。
1. 环境准备与数据预处理
1.1 工具库配置
实现2SFCA需要以下Python库支持：
import arcpy import pandas as pd import numpy as np from scipy.spatial import cKDTree import matplotlib.pyplot as plt关键工具说明：
arcpy：处理空间数据核心工具，需ArcGIS许可证pandas：表格数据处理引擎cKDTree：高效计算空间距离矩阵matplotlib：结果可视化
提示：若无法使用ArcPy，可用GeoPandas替代部分功能，但会损失部分ArcGIS专属工具
1.2 数据格式标准化
典型输入数据应包含两个要素类：
| 数据类型 | 必需字段 | 说明 |
|---|---|---|
| 需求点 | ID, Population, Geometry | 居民点/人口统计区 |
| 供给点 | ID, Capacity, Geometry | 医院/学校等设施 |
使用ArcPy进行数据检查：
def validate_data(input_fc): fields = [f.name for f in arcpy.ListFields(input_fc)] assert 'SHAPE' in fields, "缺少几何字段" if 'Population' not in fields and 'Capacity' not in fields: raise ValueError("缺失人口或容量字段")2. 核心算法实现
2.1 距离矩阵计算
采用空间索引加速计算，比传统循环效率提升50倍以上：
def create_distance_matrix(origins, destinations, max_dist): """ 构建稀疏距离矩阵 """ orig_coords = np.array([(pt.X, pt.Y) for pt in origins]) dest_coords = np.array([(pt.X, pt.Y) for pt in destinations]) tree = cKDTree(dest_coords) dists, indices = tree.query( orig_coords, k=len(dest_coords), distance_upper_bound=max_dist ) return dists, indices参数优化建议：
- 城市尺度：搜索半径建议1-5km
- 区域尺度：建议10-30km
- 医疗设施：建议使用分级半径（急诊15km，专科30km）
2.2 两步搜索实现
完整算法流程分四个阶段：
- 数据加载：将要素类转为Pandas DataFrame
- 第一步搜索：计算每个供给点的服务人口
- 第二步搜索：汇总需求点可达资源
- 结果标准化：按人口加权处理
def calculate_2sfca(demand_fc, supply_fc, threshold): # 转换为DataFrame demand_df = featureclass_to_dataframe(demand_fc) supply_df = featureclass_to_dataframe(supply_fc) # 第一步：供给点搜索 supply_df['Rj'] = supply_df['Capacity'] / demand_in_range( demand_df, supply_df, threshold ) # 第二步：需求点搜索 demand_df['Ai'] = supply_in_range( supply_df, demand_df, threshold ) return demand_df3. 可视化与结果分析
3.1 热力图生成
使用matplotlib绘制分级色彩图：
def plot_accessibility(result_df): fig, ax = plt.subplots(figsize=(12, 8)) sc = ax.scatter( result_df.geometry.x, result_df.geometry.y, c=result_df.Ai, cmap='YlOrRd', s=50, alpha=0.6 ) plt.colorbar(sc, label='可达性指数') plt.title('公共服务设施可达性分布')3.2 统计指标计算
关键评估指标：
| 指标 | 公式 | 解读 |
|---|---|---|
| 基尼系数 | $\frac{\sum\sum | A_i-A_j |
| 变异系数 | $\sigma/\mu$ | 区域差异度 |
| 首位比 | 最高值/次高值 | 资源集中度 |
实现代码：
def calculate_equity_metrics(access_series): gini = 0.5 * np.mean(np.abs(np.subtract.outer( access_series, access_series ))) / np.mean(access_series) cv = np.std(access_series) / np.mean(access_series) return {'Gini': gini, 'CV': cv}4. 性能优化技巧
4.1 并行计算加速
对于百万级数据点，采用Dask实现分布式计算：
import dask.dataframe as dd def parallel_2sfca(demand_df, supply_df): dask_demand = dd.from_pandas(demand_df, npartitions=4) results = dask_demand.map_partitions( calculate_partition, supply_df=supply_df, meta=('Ai', 'float64') ) return results.compute()4.2 内存管理
处理大型数据集时的优化策略：
- 分块处理：将研究区域划分为网格单元
- 稀疏矩阵：只存储阈值范围内的距离关系
- 数据压缩：使用category类型存储文本字段
# 内存优化示例 demand_df['Region'] = demand_df['Region'].astype('category') supply_df['Type'] = supply_df['Type'].astype('category')5. 实际应用案例
以某三甲医院选址分析为例，演示完整工作流：
数据准备：
- 居民点：街道人口统计数据
- 供给点：现有医院床位数量
参数设置：
config = { 'search_radius': 15000, # 15公里服务半径 'capacity_field': 'Beds', 'population_field': 'Pop2023' }结果解读：
- 发现西北城区存在明显服务盲区
- 基尼系数0.42显示资源分配不均
- 新院区选址建议优先考虑低值区
注意：实际应用中需考虑交通网络阻抗，而非直线距离
