SHP文件实战：山东省行政区划与路网数据GIS分析全流程

简介这份资源面向GIS从业者、城市规划人员及地理数据分析学习者提供2022年7月山东省的矢量地理数据涵盖省、市、县三级行政区划边界以及道路网与铁路网信息可用于空间查询、叠加分析、缓冲区分析及专题地图制作等场景。压缩包共33个文件约37.59MB以shp、dbf、prj、shx、xml等格式为主其中shp存储几何形状dbf保存道路名称、等级、行政代码等属性prj定义坐标系统shx与xml辅助索引和元数据描述各类型配合构成完整可用的GIS数据集。目前已有2018人学习下载适合需要快速获取山东基础地理底图、开展交通规划、灾害响应评估或商业选址分析的用户可直接在ArcGIS、QGIS等软件中加载使用省去繁琐的数据采集与整理环节。1. 拿到一份 SHP 文件之后山东省行政区划与路网数据能跑出什么很多做 GIS 分析的人都有过这种经历拿到一个压缩包解压出来一堆 .shp、.dbf、.shx、.prj双击打不开拖进 QGIS 能看但不知道坐标系对不对想按市县做统计又发现属性表里字段名全是拼音缩写。这份「山东省行政区划省、市、县道路网公路网」的 SHP 文件本质上就是一套典型的矢量地理数据集合覆盖省级、地级市、区县三级行政边界外加道路网和公路网两层线状要素。它能干的事很具体按行政区做数据聚合、做路网密度计算、做可达性分析、做专题图出图。适合谁做国土空间规划的技术人员、做物流选址的分析师、做城市研究的学生以及需要给业务系统灌底图数据的后端工程师。坐标系大概率是 CGCS2000 或 WGS84但别猜后面会讲怎么验。2. 先验坐标系和字段别让第一个图层就翻车2.1 用 ogrinfo 把元数据一次看穿拿到 SHP 第一件事不是打开 QGIS 画图而是用命令行把元数据读出来。GDAL 自带的 ogrinfo 是最稳的工具没有之一。# 查看 SHP 文件的图层摘要信息包括几何类型、要素数、坐标系 ogrinfo -so -al shandong_province.shp # 查看前 3 个要素的完整属性确认字段名和编码 ogrinfo -al -fid 1 -fid 2 -fid 3 shandong_city.shp-so表示 summary only只输出摘要不输出每个要素-al表示 all layers遍历所有图层。输出里重点看三行Geometry:告诉你这是 Polygon 还是 LineStringFeature Count:告诉你要素数量省级应该是 1地级市 16 左右区县 130 多Layer SRS WKT:里的AUTHORITY[EPSG,4490]就是坐标系代码4490 是 CGCS2000 地理坐标系4326 是 WGS84。如果 ogrinfo 输出里 SRS 显示为(unknown)说明 .prj 文件缺失或损坏。这时候不要硬猜用 QGIS 打开后手动指定坐标系再导出或者找原始数据提供方要 .prj 文件。我一般会先把所有图层的坐标系统一成 EPSG:4490因为国内行政区划数据用 CGCS2000 是主流做法后续做面积计算和距离量测时偏差最小。2.2 字段名乱码和编码问题SHP 的 .dbf 文件对中文支持一直是个老大难。用 ogrinfo 看到字段值显示为????或乱码基本就是编码问题。常见做法是设置SHAPE_ENCODING环境变量# 强制以 GBK 编码读取 dbf适用于国内老数据 SHAPE_ENCODINGGBK ogrinfo -al shandong_county.shp # 如果 GBK 不对试 UTF-8 SHAPE_ENCODINGUTF-8 ogrinfo -al shandong_county.shp参数说明SHAPE_ENCODING是 GDAL 的 shapefile 驱动专用环境变量只影响 .dbf 的文本读取不影响几何。如果两种编码都试过还是乱码那可能是数据生产时用了自定义编码页这时候用 Python 的 dbfread 库逐字节读出来看十六进制再决定怎么转。别在这上面耗太久实在不行就按 FID 顺序重新挂接一份干净的属性表。2.3 道路网和公路网的区别要先搞清楚标题里同时出现了「道路网」和「公路网」这两个不是一回事。道路网通常包含城市道路、乡村道路、匝道等所有可通行路径要素密度高、属性字段多公路网一般只包含国道、省道、县道等纳入公路管理体系的线路要素少但等级属性明确。拿到数据后先看两个图层的要素数量和字段结构图层典型要素数山东关键字段几何类型道路网数十万条道路名称、类型、宽度LineString公路网数万条路线编号、等级、路面材质LineString区县边界130区县代码、名称、面积Polygon如果道路网要素数明显偏少可能是只抽了主干道如果公路网里出现了「城市主干路」这种字段值说明数据生产方把两个概念混在一起了。这直接影响后续分析口径做路网密度时用哪个图层结论可能差一倍。3. 用 Python 把三级行政区划和路网串起来做分析3.1 环境准备与依赖安装# 推荐用 conda 建独立环境避免 GDAL 版本冲突 conda create -n shandong_gis python3.10 conda activate shandong_gis # 安装核心库geopandas 会自动带 GDAL 和 Fiona conda install -c conda-forge geopandas matplotlib contextily pip install shapely pyproj rtreegeopandas是操作矢量数据的核心库底层调 GDAL 做 I/O调 Shapely 做几何运算。contextily用来加底图出图时有用。rtree是空间索引库做空间连接时能快一个数量级。版本上 geopandas 0.14 以上对 pandas 2.x 兼容更好如果遇到AttributeError: GeometryArray object has no attribute crs这类报错八成是版本不匹配回退到 0.12 或升级到最新。3.2 读取三级行政区划并统一坐标系import geopandas as gpd # 读取三个层级的行政区划 province gpd.read_file(shandong_province.shp, encodingGBK) city gpd.read_file(shandong_city.shp, encodingGBK) county gpd.read_file(shandong_county.shp, encodingGBK) # 检查坐标系统一转为 EPSG:4490 for name, gdf in [(province, province), (city, city), (county, county)]: print(f{name}: CRS {gdf.crs}, 要素数 {len(gdf)}) if gdf.crs is None: gdf.set_crs(epsg4490, inplaceTrue) elif gdf.crs.to_epsg() ! 4490: gdf.to_crs(epsg4490, inplaceTrue) # 验证省级图层应该只有 1 个要素 assert len(province) 1, 省级图层要素数异常 # 验证区县图层要素数应大于地级市图层 assert len(county) len(city), 区县数应多于地级市数encodingGBK参数直接传给 Fiona 的打开选项解决中文乱码。set_crs和to_crs的区别要分清前者是「声明」坐标系用于原本没有 CRS 信息的情况后者是「转换」坐标系会实际重算坐标值。如果数据本身有 CRS 但你不确定对不对先用to_crs转成 4490 再看范围是否落在山东经纬度区间东经 114°–122°北纬 34°–38°偏得离谱就说明原始 CRS 标错了。3.3 按区县统计路网密度路网密度是路网分析里最常用的指标算法是「区域内道路总长度 ÷ 区域面积」。这里有个坑地理坐标系下直接算长度得到的是度数不是米。必须先投影到投影坐标系。# 读取路网数据 road gpd.read_file(shandong_road.shp, encodingGBK) # 投影到 CGCS2000 3 度带山东大致在 39 带中央经线 117°E # EPSG:4547 是 CGCS2000 / 3-degree Gauss-Kruger CM 117E road_proj road.to_crs(epsg4547) county_proj county.to_crs(epsg4547) # 计算每条道路的长度米 road_proj[length_m] road_proj.geometry.length # 空间连接把道路按位置归属到区县 road_joined gpd.sjoin(road_proj, county_proj[[区县代码, geometry]], howinner, predicateintersects) # 按区县汇总道路总长度 road_length_by_county road_joined.groupby(区县代码)[length_m].sum().reset_index() road_length_by_county.columns [区县代码, 道路总长_m] # 计算区县面积平方公里 county_proj[area_km2] county_proj.geometry.area / 1e6 # 合并计算路网密度 result county_proj.merge(road_length_by_county, on区县代码, howleft) result[路网密度_km_per_km2] (result[道路总长_m] / 1000) / result[area_km2] # 输出密度最高的 10 个区县 print(result.nlargest(10, 路网密度_km_per_km2)[[区县代码, 路网密度_km_per_km2]])关键参数说明epsg4547是 CGCS2000 高斯克吕格 3 度带投影中央经线 117°E适合山东大部分地区。如果数据覆盖范围偏东如威海、烟台东部可以用 4548中央经线 120°E减少投影变形。predicateintersects表示道路与区县边界相交即算归属如果一条路跨两个区县会被分别计入两个区县这是合理的——路网密度本来就是按行政单元统计的。howinner会丢弃落在区县边界外的道路如果发现汇总后总长度比原始数据少很多检查是否有道路坐标偏移到了省界外。3.4 用公路网等级字段做分层统计公路网数据通常带等级字段按等级分层统计能看出路网结构。# 假设公路网有 等级 字段值为 高速、国道、省道、县道 highway gpd.read_file(shandong_highway.shp, encodingGBK) highway_proj highway.to_crs(epsg4547) highway_proj[length_m] highway_proj.geometry.length # 按等级汇总 level_stats highway_proj.groupby(等级).agg( 总长度_km(length_m, lambda x: x.sum() / 1000), 路段数(length_m, count), 平均段长_km(length_m, lambda x: x.mean() / 1000) ).reset_index() print(level_stats.sort_values(总长度_km, ascendingFalse))如果字段名不是「等级」而是拼音如dengji或英文road_level先用highway.columns.tolist()打印所有字段名确认。agg里对同一列用了三个不同聚合函数这是 pandas 的命名聚合语法输出列名就是等号左边的名字。平均段长这个指标能反映数据采集的粒度如果高速的平均段长只有 0.5 公里说明数据被切得很碎可能是按匝道口分段了做路径分析时要先合并。4. 避坑与排查SHP 处理中最容易翻车的五个地方4.1 中文乱码反复出现改了编码还是不对现象用 geopandas 读取时设置了encodingGBK但字段值仍然显示为乱码或问号。原因SHP 的编码问题可能出在 .dbf 文件头部的 language driver ID 上GDAL 有时会忽略 encoding 参数而优先读文件头。解决先用ogrinfo确认乱码是否存在如果 ogrinfo 正常而 Python 不正常用gpd.read_file(..., encodingGBK)的同时加SHAPE_ENCODINGGBK环境变量如果都不行用 dbfread 单独读属性表再手动挂接。4.2 空间连接后要素数暴增现象sjoin之后行数比原始道路数多了好几倍。原因一条道路与多个区县相交时会被复制多份这是空间连接的正常行为但如果不做后续聚合就会导致重复计算。解决明确分析目标——如果要统计每个区县的道路长度sjoin 后必须 groupby 聚合如果要保留道路完整属性用predicatewithin只保留完全落在区县内的道路或者用dissolve先合并区县再连接。4.3 面积计算出来是负数或极小值现象geometry.area算出来是 0.00x 这种数量级。原因数据还在地理坐标系度下没有投影。解决任何涉及长度、面积的计算之前必须to_crs到投影坐标系。山东用 EPSG:4547 或 4548全国范围用 Albers 等积投影EPSG:9822 或自定义。记住一个原则地理坐标系只用来存和传投影坐标系才用来算。4.4 道路网和行政区划边界对不上现象道路看起来偏移到了海里或者省界外。原因两个图层的坐标系不一致或者其中一个的 CRS 标错了。解决分别用gdf.crs打印两个图层的 CRS确认 EPSG 代码一致如果不一致统一转成同一个。如果 CRS 一致但仍有偏移检查是否一个用了 CGCS2000 一个用了 WGS84——两者在山东地区差几十米视觉上不明显但做精确分析时会有影响。4.5 大文件读取内存爆掉现象读取道路网 SHP 时进程被 kill 或报 MemoryError。原因SHP 格式单文件上限 2GB但要素数几十万时 geopandas 全量加载仍然吃内存。解决用gpd.read_file(..., rowsslice(0, 50000))分批读或者用 Fiona 的迭代器逐要素处理。如果只是做统计不需要几何运算用gpd.read_file(..., ignore_geometryTrue)只读属性表内存占用能降一个数量级。5. 进阶技巧用区划数据做路网可达性插值前面讲的都是「统计」这一章讲「分析」。路网可达性是个典型场景给定若干节点比如物流园区、医院计算每个区县到最近节点的路网距离。这个分析比直线距离靠谱得多但实现上有个关键选择——用路网距离还是用欧氏距离做插值。我的做法是先用区县质心做欧氏距离的快速估算筛出候选区县再对这些区县做精确的路网距离计算。这样能把计算量从 O(n×m) 降到 O(k×m)k 是候选区县数通常只有总数的三分之一。import numpy as np from shapely.ops import nearest_points # 假设 facilities 是设施点图层county_proj 是投影后的区县 facilities_proj facilities.to_crs(epsg4547) # 计算每个区县质心到最近设施点的欧氏距离 county_proj[centroid] county_proj.geometry.centroid county_proj[euclid_dist_m] county_proj[centroid].apply( lambda pt: facilities_proj.distance(pt).min() ) # 取欧氏距离最小的前 30% 区县做精确路网分析 threshold county_proj[euclid_dist_m].quantile(0.3) candidates county_proj[county_proj[euclid_dist_m] threshold].copy() # 对候选区县用最近设施点的路网距离替代 # 这里简化处理实际应用需调 pgRouting 或 networkx 做路网最短路径 candidates[road_dist_m] candidates[euclid_dist_m] * 1.4 # 迂回系数经验值 # 合并结果 county_proj county_proj.merge( candidates[[区县代码, road_dist_m]], on区县代码, howleft ) county_proj[final_dist_m] county_proj[road_dist_m].fillna( county_proj[euclid_dist_m] )迂回系数 1.4 是平原路网的常见经验值山东大部分地区适用。山区或水网密集区可能要调到 1.6 以上。这个系数的校准方法选几个已知路网距离的样本点用实际值除以欧氏距离取平均。如果手头有路网拓扑数据直接用 networkx 构建图做 Dijkstra 更准但构建图的时间成本不低几十万条路段大概要几分钟。出图时用matplotlib配contextily加底图把可达性用分级设色渲染出来import matplotlib.pyplot as plt fig, ax plt.subplots(1, 1, figsize(12, 10)) county_proj.plot(columnfinal_dist_m, cmapRdYlGn_r, legendTrue, axax, edgecolorgray, linewidth0.3) ax.set_title(山东省区县路网可达性到最近设施点) ax.axis(off) plt.tight_layout() plt.savefig(accessibility.png, dpi200, bbox_inchestight)cmapRdYlGn_r是红黄绿反转色带距离越近越绿符合直觉。edgecolorgray给每个区县加细边界线出图时层次更清楚。dpi 设 200 以上打印或放报告里不糊。最后说个血泪教训做可达性分析时千万别直接用区县质心代表整个区县。面积大的区县比如沂蒙山区那些质心可能落在无人区算出来的距离没有实际意义。我一般会按人口或 GDP 做加权质心或者干脆用格网做分析单元。这个坑我踩过不止一次希望帮到你。本文还有配套的精品资源点击获取