大概半年前我接到一个需求把某个城区的公交站点按500米服务半径做覆盖分析评估覆盖了多少面积、覆盖了多少常住人口。这类问题放到ArcGIS里点几下就能出图但这次卡着一个要求——结果要可复现、半径参数要随时能改、后续还要接其他交通模型。我转了一圈最后用geopandas写了一套覆盖分析脚本从清洗数据到输出覆盖率报表一次性打通。上周又有同事问起这套逻辑我盘了一下决定写成这篇东西把覆盖率的业务口径、坐标系和数据准备的坑、核心计算流程以及我实际跑数据时踩过的几个错误完整讲一遍。1. 公交站点覆盖率为什么值得用geopandas算一遍1.1 覆盖率的业务口径它到底在衡量什么公交站点覆盖率简单说就是以公交站点为圆心按一个步行可达距离画出服务圈再看这些服务圈在研究区里覆盖了多少面积、多少人口。它回答的问题是居民从家里走到公交站的这段路是不是大部分人能在可接受的时间内走完。在公共交通规划里最常用的服务半径是300米、500米和800米。300米约等于步行3到4分钟通常用来衡量高密度城区的高频站点服务500米大约步行6到8分钟是一般公交线网规划的基准值800米常被拿来评估远郊或者低密度片区的覆盖下限。这意味着做覆盖率分析时不应该只跑一个半径至少要按300/500/800三档都计算一遍才能看出不同服务强度下的差异。核心口径其实就两个面积覆盖率 公交站点服务范围与研究区交集的面积 / 研究区总面积人口覆盖率 公交站点服务范围内的人口数 / 研究区总人口这两个指标单独看都有盲区后面第4章会专门讲怎么组合解读。另外这套分析逻辑不只是公交能用把站点数据换成公园入口、菜市场、幼儿园同样的代码就可以做公园绿地覆盖率、生活服务设施覆盖率。所以它是空间统计里非常通用的一套框架。1.2 为什么最终选geopandas而不是ArcGIS或QGIS工具选型这事我纠结过一阵。ArcGIS功能全网络分析尤其强但授权成本高脚本化又依赖ArcPy在服务器环境里部署很麻烦。QGIS免费可视化很好模型构建器也能做批处理可一旦要嵌入数据管道、频繁调参、周期性重跑还是不如直接写代码灵活。维度ArcGIS / ArcPyQGISgeopandas授权与部署商业授权部署较重免费窗口环境为主免费可无界面运行批量与参数化一般依赖ArcPy环境中等需模型构建器或Python API天然脚本化与Python生态衔接较弱中等极好pandas语法直接上手适合场景复杂网络分析一次性可视化和探索可复现数据管道、批量测算我最后的结论是如果只是单独做一次分析、快速出图用QGIS可能更顺手但如果你要的是参数化、批量、可持续更新的数据流程geopandas的优势非常明显。它本质上是pandas加上shapely和pyproj的封装把矢量空间操作变成DataFrame里的一列操作学过pandas的人很快就能适应。而且它能在没有任何图形界面的服务器上运行方便接到自动化调度里。1.3 分析的整体路径覆盖率分析跑通一次不难难的是别漏步骤。我建议所有项目都按下面这条路径走收集公交站点、行政区边界和人口统计单元数据统一坐标系并清洗站点坐标按服务半径生成缓冲区并合并重叠区域与研究区求交计算面积覆盖率和人口覆盖率按街道拆分统计出图出表。下面每个环节我都会展开讲并给出可以直接跑的代码。尤其是第2步坐标系问题如果处理不好后面的缓冲区形状和面积数字全是错的而且这个坑特别隐蔽。2. 数据准备坐标系、数据格式、清洗流程2.1 三类基础数据做一次完整的覆盖率分析至少需要三类数据公交站点点数据每个站点有名称和经纬度坐标。可以从公交线路表里整理出来也可以从开放数据平台下载。格式无所谓CSV、GeoJSON、Shapefile都可以geopandas都能读。研究区边界面数据最基础的是区县边界或街道边界。它有两个作用一是作为覆盖率的分母二是用于裁切覆盖范围避免把研究区外的站点服务圈算进来。人口数据如果没有人口数据就只能输出面积覆盖率也能说明一定问题。但真正做规划评估时人口覆盖率才是更贴近实际的指标。人口数据通常是带人口字段的统计单元多边形比如社区边界、普查区划。如果只需要快速出面积覆盖率其实数据就两类。但我的建议是尽量把人口数据也准备好因为面积覆盖率在很多场景下会骗人后面会细说。2.2 坐标系是覆盖率分析的第一道关这是我认为整个分析里最重要、也最容易被新手踩坑的地方。通俗地讲地理坐标系比如WGS84EPSG:4326存储的是经纬度单位是度投影坐标系比如UTM投影带EPSG:326xx把地球表面投影到一个平面上单位是米。覆盖率要算面积还要生成500米这种带真实距离的缓冲区就必须在投影坐标系下进行。直接在EPSG:4326的数据上执行buffer(500)表面上看是半径500实际上执行的是在经度和纬度方向各偏移500度在真实世界里完全不是半径500米。可以这样理解把地图想象成橘子皮经纬度是在球面上的位置想要按米数去量距离就得先把球面摊平到桌面上这个摊平的方式就是投影。在纬度较高的地方经度1度对应的地面距离远小于纬度1度所以不投影直接buffer出来的形状还会变成一个横向和纵向比例完全失真的椭圆。选择投影时一般选研究区所在的UTM投影带或者直接使用研究区所在城市的地方坐标系。核心原则只有一个所有参与分析的图层必须统一到同一个投影坐标系再做缓冲区和面积计算。2.3 站点数据清洗三个必然遇到的细节站点数据通常比想象中脏我每次都会遇到以下三类问题。第一个是重复。公交数据经常按线路存储同一个站名、同一个位置在数据里会出现很多次。做覆盖率分析前必须去重不然同一个位置会生成多个重叠缓冲区后面如果做站点密度类统计数字会虚高。我常用的做法是按站名和坐标双重判断把坐标四舍五入到小数点后约4位约10米级别相同的归并为一个站点。第二个是异常值。经纬度为0、明显超出研究区范围的点都要在进入分析前剔除。这里建议保留清洗后的中间文件后面多次跑不同半径时就不需要重复清洗。第三个是坐标基准偏移。如果站点坐标来自国内某些在线地图平台坐标可能存在偏移而有些公开底图是WGS84基准两者叠加时站点可能会偏离真实位置几百米。这里提醒一句geopandas的to_crs只能做坐标系数学变换不能解决坐标基准偏移。所以尽量让所有数据来自同一基准或者先统一纠偏再做分析否则缓冲区会盖在错误的位置上覆盖率数字自然也不对。3. 核心实现从坐标点到覆盖率结果3.1 读取数据并统一坐标系先通过CSV读入站点数据再组装成GeoDataFrame。这里有一个重要细节CSV文件本身不携带坐标系信息所以读进来之后要手动指定crsEPSG:4326否则后面调用to_crs时geopandas根本不知道数据在哪个坐标参考系里。import pandas as pd import geopandas as gpd from shapely.geometry import Point # 1. 站点CSV转GeoDataFrame stations_df pd.read_csv(bus_stations.csv) stations gpd.GeoDataFrame( stations_df, geometrygpd.points_from_xy( stations_df[lng], stations_df[lat] ), crsEPSG:4326 ) # 2. 读取边界与人口数据 district gpd.read_file(district.geojson) blocks gpd.read_file(population_blocks.geojson) # 3. 统一投影坐标系 LOCAL_CRS EPSG:32650 # 按研究区所在UTM带调整 stations stations.to_crs(LOCAL_CRS) district district.to_crs(LOCAL_CRS) blocks blocks.to_crs(LOCAL_CRS)研究区如果横跨多个UTM带就要谨慎使用单一带号了这个放到第5章讲。现阶段只要保证所有图层都用同一个投影坐标分析逻辑就成立。3.2 生成缓冲区并合并重叠接下来生成每个站点的服务范围。站点数量几千个时buffer运算非常快不用太担心性能。# 生成500米服务范围 coverage_geom stations.buffer(500).unary_union()buffer(500)会在每个站点周围生成半径500米的圆形多边形unary_union()则会把所有重叠的圆合并成一个整体的覆盖范围。这一步的意义在于站点之间服务圈大量重叠如果不合并直接把每块面积相加覆盖率会远超100%完全没有意义。这里有个参数可以提一下buffer方法有个quad_segs参数控制圆弧被分割成多少段默认值是16对于常规分析足够用。如果对边界光滑度要求极高可以调大到32或更高但生成的多边形顶点会更多文件体积也会变大。3.3 面积覆盖率与人口覆盖率计算覆盖范围生成后先用它与研究区求交再做面积占比这是面积覆盖率的正解。# 覆盖率-面积口径 cover_intersect district.geometry.intersection(coverage_geom) total_area district.geometry.area.sum() covered_area cover_intersect.area.sum() area_coverage covered_area / total_area print(f面积覆盖率: {area_coverage:.2%})人口覆盖率有两种算法。第一种粗暴但快只要人口统计单元与覆盖范围相交就把整个单元的人口计入第二种是面积加权更精确。# 覆盖率-人口口径方法一单元全部计入偏粗 covered_blocks blocks[blocks.geometry.intersects(coverage_geom)] pop_covered covered_blocks[population].sum() pop_total blocks[population].sum() pop_coverage_method1 pop_covered / pop_total # 覆盖率-人口口径方法二面积加权更精确 overlap_ratio blocks.geometry.intersection( coverage_geom ).area / blocks.geometry.area blocks[covered_pop_est] blocks[population] * overlap_ratio pop_coverage_method2 blocks[covered_pop_est].sum() / pop_total方法二的基本假设是人口在统计单元内均匀分布。如果统计单元划分得细两种方法差异很小但如果统计单元面积很大比如一个街道只有一个人口总数而它的边缘恰好被覆盖了一小块方法一会严重高估覆盖人口方法二就合理得多。在正式报告里我一般使用方法二的数据。3.4 多半径批量计算与街道级统计实际项目很少只跑一个半径。我通常会一次跑三档用一个循环搞定radius_results {} for radius in [300, 500, 800]: geom stations.buffer(radius).unary_union() covered district.geometry.intersection(geom).area.sum() radius_results[radius] covered / total_area如果需要按街道拆分统计别去对每个街道单独跑一遍缓冲区后再求交。更好的做法是先生成整体覆盖范围再与每个街道的多边形求交这样能天然处理站点在研究区外但服务圈伸进研究区的情况street_stats [] for _, row in district.iterrows(): inter_geom row.geometry.intersection(coverage_geom) street_area row.geometry.area street_stats.append({ name: row[name], covered_area: inter_geom.area, area_ratio: inter_geom.area / street_area }) street_df pd.DataFrame(street_stats) street_df.to_csv(street_coverage.csv, indexFalse)这个处理顺序是我踩过坑后才固定的。如果反过来先判断哪个街道离站点近再逐街道单独求覆盖很容易漏掉边界外站点的贡献导致街道覆盖率整体偏低。3.5 完整代码骨架把上面的片段串成一个函数方便后续调整参数复用import pandas as pd import geopandas as gpd def coverage_analysis(bus_stations_csv, boundary_path, pop_block_pathNone, radii(300, 500, 800)): stations_df pd.read_csv(bus_stations_csv) stations gpd.GeoDataFrame( stations_df, geometrygpd.points_from_xy(stations_df[lng], stations_df[lat]), crsEPSG:4326 ).to_crs(LOCAL_CRS) district gpd.read_file(boundary_path).to_crs(LOCAL_CRS) total_area district.geometry.area.sum() results {} for radius in radii: coverage_geom stations.buffer(radius).unary_union() covered_area district.geometry.intersection( coverage_geom ).area.sum() results[radius] round(covered_area / total_area, 4) # 可选人口覆盖率 if pop_block_path: blocks gpd.read_file(pop_block_path).to_crs(LOCAL_CRS) pop_total blocks[population].sum() for radius in radii: coverage_geom stations.buffer(radius).unary_union() overlap_ratio blocks.geometry.intersection( coverage_geom ).area / blocks.geometry.area pop_covered_est (blocks[population] * overlap_ratio).sum() results[fpop_{radius}] round(pop_covered_est / pop_total, 4) return results这个函数的设计思路是半径参数留空调用时传什么就跑什么人口数据可选传不传就只输出面积覆盖率。在实际项目里这种一个函数管多档参数的方式比对着界面点来点去高效太多。4. 结果解读与可视化数字背后的问题4.1 两个指标分开看面积口径容易被没人住的地方拉低算出了覆盖率数字不等于能回答业务问题。我发现很多人把面积覆盖率当成唯一结论这其实不够。面积覆盖率是非常保守的口径。研究区里如果有大片的公园、水域、山体或者待开发空地这些区域虽然可以被站点服务圈覆盖但几乎没有常住人口把它们计入覆盖率会拉低整体数字。所以做面积覆盖率时建议先从研究区里扣除水域和不可进入的用地否则不同城市、不同城区的对比会失公平。人口覆盖率则相反它更贴近居民是否方便走到公交站这个真实问题。两组数字组合起来解读时大致可以有三种判断面积覆盖率与人口覆盖率接近说明站点空间分布均匀但没有明显突出人口集聚区。面积覆盖率明显低于人口覆盖率说明站点相对集中在居住区人口密度高的地方覆盖得好但全域服务并不均衡。面积覆盖率明显高于人口覆盖率说明很多站点铺在了人口稀疏的区域比如工业园区、公园周边对常住人口的效率不高。覆盖率不是越高越好关键还要看覆盖到的人是谁、在哪。4.2 模拟数据结果解读示例这里演示一组模拟数据的解读逻辑不代表任何真实城市。假设研究区面积110平方公里常住人口76万公交站点420个。500米服务半径下面积覆盖率47.8%人口覆盖率63.5%。两个指标相差约15个百分点说明站点分布相对偏向居住集聚区服务效率尚可但全域空间均衡性一般。再把街道拆开看人口密度最高的几个街道覆盖率超过75%而两个低密度街道只有22%到30%。这类结果在报告里怎么用通常可以支撑两类结论一是中心城区站点加密已接近饱和下一步应关注中低密度片区的线网延伸二是面积覆盖率偏低的区域需要结合用地性质判断新增站点的优先级而不是简单按面积补盲。4.3 出图与报告呈现的细节覆盖率分析最终总要落到图上。我的标配出图方式是服务圈用半透明蓝色面站点用红色点边界用灰色线。import matplotlib.pyplot as plt fig, ax plt.subplots(figsize(10, 12)) district.boundary.plot(axax, colorgray, linewidth1) gpd.GeoSeries(coverage_geom).plot( axax, color#74c7b8, alpha0.55 ) stations.geometry.plot( axax, color#e74c3c, markersize6, alpha0.8 ) plt.axis(off) plt.tight_layout() plt.savefig(coverage_map.png, dpi300)报告里除了总图我还会同时给出两张表。一张是不同半径下的总体覆盖率对照表另一张是分街道覆盖率排行表。表格比图更容易让业务方快速定位问题区域尤其当街道数量很多时。这里再强调一次导出的统计结果一定要带CSV或Excel源文件方便后续管线复用而不只是给一张截图。5. 实战中的坑与优化建议5.1 椭圆缓冲区不是算法bug是坐标系没统一我第一次跑通脚本后盯着出图发现缓冲区整体像被压扁的椭圆直觉告诉我这不是正常的圆形服务范围。排查了半天才发现站点数据读取时没有手动指定CRSGeoDataFrame的crs属性是空的我直接对空坐标系的几何执行了buffergeopandas没有报错但得到的几何体跟投影后的边界叠加时完全对不上。旧事重提是想强调排查顺序遇到覆盖率明显不合理时第一件事不是检查算法而是打印所有图层的crs确认它们是否一致。读取CSV时立刻指定crsEPSG:4326这个习惯能帮你省掉大量排查时间。另外to_crs只有在GeoDataFrame已经带CRS信息时才有效。如果crs是None转换不会发生甚至可能报错。所以数据加载阶段就要把CRS问题解决掉。5.2 投影带对整个区域真的准确吗UTM投影在经度跨度不大的城市没问题但如果研究区范围很大比如跨了多个UTM带单一投影带在边缘区域的面积和距离误差就会明显放大。覆盖率是分子分母同源的比值误差会被部分抵消但如果报告里要写覆盖面积是多少平方公里这种绝对数字就要留意。还有一种场景是研究区不在标准UTM带中央比如城市正好卡在带边界这时可以选用适合该区域的地方投影或者使用自定义中央经线的投影参数。原则很简单所有图层统一到同一个投影且这个投影对研究区造成的变形尽量小。5.3 重叠缓冲区与重复计算这是新手最容易算错的地方。多个站点的服务圈只要彼此重叠就不能简单把每个站点的缓冲区面积加起来当总覆盖面积。正确做法永远是把缓冲区合并成一个整体几何再求面积。类似地判断某个小区是否在覆盖范围内时用小区中心点代表整个小区也会漏掉边缘被覆盖的情况尤其当小区多边形很大时中心点法会产生肉眼可见的偏差。我实际项目里宁可多做一步面积加权也尽量避免中心点判断。多花的时间很少但结果的可靠性完全不一样。5.4 大数据量下的性能优化站点数量几千个时buffer和unary_union几秒钟就能完成不需要特殊优化。但做人口覆盖率时如果人口统计单元有几万个逐条做intersects就会明显变慢。这种情况下可以用空间连接geopandas的sjoin内部会走R树空间索引比全量相交快一个数量级joined gpd.sjoin( blocks[[geometry, population]], gpd.GeoDataFrame(geometry[coverage_geom], crsLOCAL_CRS), howinner, predicateintersects )如果是循环里反复对同一个大几何做判断还可以用shapely.prepared.prep把几何预处理一下批量点查询时提速明显from shapely.prepared import prep coverage_prep prep(coverage_geom) blocks[covered_flag] blocks.geometry.apply( lambda g: coverage_prep.intersects(g) )内存方面读入人口数据后可以先只保留必要的字段中间过程及时清理临时GeoDataFrame。如果数据太大跑不动按街道分片计算再汇总也是可行的效率虽然低一点但至少不会把内存撑爆。做覆盖分析到现在我慢慢形成了固定的习惯动手之前先确认两个零——坐标系是否是投影坐标、重复站点是否归零然后直接跑300/500/800三档半径同时输出面积口径和人口口径给业务方两张表、两张图。这套流程用geopandas维护成本很低半径参数改了重跑一遍就行也方便把结果递交接续的线路优化模型。如果你也正准备做类似的空间覆盖率分析我建议直接按这个思路起步应该能少走不少弯路。