gasfca 可达性分析 —— 全流程技能

本技能覆盖 gasfca（Ga2SFCA）可达性分析的完整流程：数据准备 → 算法原理 → 参数设置 → 运行计算 → 结果验证与解读。

项目仓库：https://github.com/Yao2028908/gasfca

1. 全流程概览

数据准备(AOI/路网/边界)
  → 读取与清洗
  → 点/面自动适配
  → (可选)边界裁剪
  → 按类别拆分供给/需求
  → run_gasfca 参数设置与运行
  → 多格式输出(shp/gpkg/geojson/csv/xlsx/栅格/图)
  → 结果验证与解读

2. 数据源约定

| 数据 | 推荐格式 | 说明 | |---|---|---| | AOI（面）或供给/需求点 | shp / gpkg / geojson | 需含分类字段（如 类别）与权重字段（如 value） | | 路网 | gpkg（OSM 导出） | 常用图层如 gis_osm_roads_free | | 研究区边界 | gpkg / geojson | 可选，用于按研究区/核心区裁剪 |

3. 点/面自动适配（重要）

包内 _polygon_to_centroid 自动处理面数据：投影 → 计算面积（area 字段）→ 转质心
点数据直接使用；面数据无需手动转质心
日志中会打印 polygon features detected -> computing area and converting to centroids

4. 算法原理（Ga2SFCA）

两步浮动捕捉区法（Gaussian 2-Step Floating Catchment Area）：

第一步（供给端）：对每个供给点 j，在阈值 d0 范围内统计加权需求总量，计算供需比 Rj = Sj / Σ(Pk · G(dkj))，其中 Sj 为供给容量，Pk 为需求规模，G 为高斯衰减权重
第二步（需求端）：对每个需求点 i，在阈值 d0 范围内累加各供给点的供需比 Ai = Σ(Rj · G(dij))

高斯衰减函数： G(d) = [exp(-½(d/d0)²) - exp(-½)] / [1 - exp(-½)]（d ≤ d0），d > d0 时 G=0

含义：Ai 越大表示该点周边"人均可用供给"越多，可达性越好。

5. 标准处理流程

内置边界裁剪（推荐）：run_gasfca 支持 boundary_path / boundary_layer / boundary_buffer 参数，包内自动完成供给、需求（within 裁剪）与路网（缓冲裁剪）的边界裁剪，无需手动 sjoin/clip。以下手动画线裁剪仅作备用。

读取数据：gpd.read_file(path, engine="pyogrio")
（可选）裁剪到研究区：需限定范围时读边界 boundary，gpd.sjoin(aoi, boundary, predicate="within", how="inner")；全量分析可跳过
按类别拆分：
- 需求固定为住宅/需求大类，demand["value"] = 1.0
- 供给 = 除需求大类外的每个类别，逐类循环（每个类别一个案例）
（可选）路网裁剪：数据量大时对边界做小缓冲后 gpd.clip(roads, clip_box) 缩小计算规模；全量路网可直接使用
坐标系：target_crs 用投影坐标（如 EPSG:32650 UTM 或 3857），保证米制距离；boundary_buffer 单位即 target_crs 单位（米）
输出准备：统一 output_dir，case_name 传入类别名 → xlsx 自动合并为不同 Sheet

6. 参数详解

| 参数 | 说明 | 建议 | |---|---|---| | supply_field / demand_field | 供给容量 / 需求规模数值字段 | 面数据用 area（包内自动计算，按面积计容量）；点数据用 value，权重统一 1.0 时按个数计容量 | | d0_list | 距离阈值列表（米） | 社区级 1000–2000m；城市级 3000–5000m；多尺度对比用多个阈值 | | target_crs | 目标投影坐标系 | 必须为米制投影（UTM/3857），决定距离精度 | | cell_size | 栅格分辨率（米） | 与数据尺度匹配：街区级 50–200m | | output_formats | 输出格式 | 矢量 shp/gpkg/geojson；表格 csv/xlsx | | case_name | 案例名（xlsx Sheet 名） | 按类别传入以区分合并结果 | | boundary_path/layer/buffer | 边界裁剪 | 无边界时留 None | | auto_split_road | 路网交点打断 | 默认 True，修复拓扑 | | generate_plots | 分位数可视化 | 默认 True，出图 |

7. 典型脚本骨架

完整可运行示例随技能附带：examples/run_gasfca.py（修改顶部 CONFIG 配置后直接运行）。

import os
import geopandas as gpd
from shapely.geometry import box
from gasfca import run_gasfca

BASE = os.path.dirname(os.path.abspath(__file__))
aoi = gpd.read_file("aoi.shp", engine="pyogrio")

# (Optional) clip AOI to a study region / boundary
# boundary = gpd.read_file("boundary.gpkg", engine="pyogrio").to_crs(aoi.crs)
# aoi = gpd.sjoin(aoi, boundary, predicate="within", how="inner")

DEMAND_CAT = "住宅"
demand = aoi[aoi["类别"] == DEMAND_CAT].copy()
demand["value"] = 1.0                    # point data: count-based weight
demand.to_file("demand.geojson", driver="GeoJSON")

roads = gpd.read_file("road.gpkg", layer="gis_osm_roads_free", engine="pyogrio")
# (Optional) clip roads to control data size; skip for full-network analysis
# buf = 0.02
# b = boundary.total_bounds
# clip_box = box(b[0]-buf, b[1]-buf, b[2]+buf, b[3]+buf)
# roads = gpd.clip(roads, clip_box)
roads.to_file("roads.geojson", driver="GeoJSON")

for cat in [c for c in sorted(aoi["类别"].unique()) if c != DEMAND_CAT]:
    supply = aoi[aoi["类别"] == cat].copy()
    supply["value"] = 1.0
    supply.to_file("supply_tmp.geojson", driver="GeoJSON")
    run_gasfca(
        supply_path="supply_tmp.geojson",
        demand_path="demand.geojson",
        road_path="roads.geojson",
        supply_field="value", demand_field="value",   # polygon supply: use "area" instead of "value"
        d0_list=[2000, 5000],
        target_crs=32650,
        cell_size=200,
        output_dir="result",
        output_formats=("shp", "gpkg", "geojson", "csv", "xlsx"),
        case_name=cat,          # 每个类别一个 Sheet
        # Optional boundary clipping (built-in; replaces the manual steps above)
        # boundary_path="boundary.gpkg",
        # boundary_layer=None,
        # boundary_buffer=2000,   # meters (target_crs units)
        generate_plots=True,
    )

8. 输出结构与解读

result/
├── vector/accessibility.{shp,gpkg,geojson}   # 矢量结果（共享目录时仅保留最后一个类别）
├── table/accessibility.csv                    # 属性表
├── table/accessibility.xlsx                   # 多案例合并（每类别一个 Sheet）
├── raster/Ai_*.tif                            # IDW 栅格
└── raster/quantile_plots/*.png                # 分位数可视化

Ai_<d0> 字段：该阈值下的可达性指数，越大越好；跨需求点可比较
栅格：IDW 插值连续面，用于制图与空间分布判断
分位数图：每张图独立分位分级，适合偏态数据

9. 结果验证

无 NaN：需求点结果不应出现 NaN；若有，检查供给/需求字段与匹配
不可达点占比：Ai = 0 表示 d0 内无任何供给，可量化"服务盲区"
数值合理性：均值与分布符合直觉（如城市中心高、边缘低）；多个 d0 下趋势一致
多类别对比：读取 xlsx 各 Sheet 计算均值，比较不同供给类别的覆盖水平
图层回读：导出结果用 gpd.read_file 回读验证要素数与字段完整

10. 注意事项

面数据权重按面积计：包内为多边形自动计算 area 字段，supply_field 直接设为 "area" 即按面积计供给容量（面积已在 target_crs 下计算）；点数据才用 value=1.0（按个数计容量）
边界/路网裁剪均非必须：全量分析直接使用全部 AOI 与路网；仅当数据量大、需控制计算耗时或限定研究范围时才裁剪
首选包内置参数 boundary_path / boundary_layer / boundary_buffer 完成裁剪，日志会打印各要素裁剪前后数量
单类别运行耗时约 15~20s（中小数据规模），可接受
共享 output_dir 时矢量/栅格会被后运行的类别覆盖；如需全部保留，按类别分子目录
数据量过大时优先裁剪路网，路网规模决定最短路径计算耗时
供给点过少（如个位数）时可达性分布可能失真，需谨慎解读