省域SHP包实战：从坐标系到空间统计的完整流程-育师

简介：面向GIS从业者、城市规划与交通研究人员及地理数据学习者，这份资源提供2022年7月山东省的矢量地理底图，可用于行政区划制图、路网可达性分析、物流选址与灾害应急路线评估等场景。压缩包共33个文件、约37.59MB，以shp几何文件、dbf属性表、prj坐标系统定义、shx索引及xml元数据为主，另有sbn、sbx空间索引文件，覆盖省、市、县三级行政边界以及道路网、铁路网图层，属性表记录名称、等级、里程等信息，prj文件保证投影一致、定位准确。已有2028人学习下载，说明数据在同类资源中具备一定参考价值。读者可直接在ArcGIS、QGIS等软件中加载，完成空间查询、叠加与缓冲区分析，并快速出图展示山东行政边界与交通网络，省去自行采集与整理底图的时间。

1. 拿到一份 2022.07 的省域 SHP 包，先想清楚它能干什么

一份标注了 2022.07 的省域行政区划加道路网、公路网 SHP 包，落到手里第一件事不是双击打开看颜色，而是判断它能不能撑起你要做的那个分析。省、市、县三级面图层加上线状路网，这套组合最常见的用途是空间叠加统计：把点位、栅格、业务表格挂到县界上算密度，或者用路网做可达性、缓冲区、最近设施匹配。它适合做区域规划、选址评估、物流路径粗算、专题制图这类工作，不适合做导航级路径规划，也不适合做需要实时路况的动态分析。很多人拿到 SHP 就急着往 GIS 软件里拖，结果坐标系一错、拓扑一乱，后面所有统计全是玄学。这一章先把这份数据的边界讲清楚，后面几章再一步步落到可复现的处理流程上。

2. 省域三级面图层与路网图层的结构拆解：字段、坐标系与拓扑先摸清

2.1 一个 SHP 包实际包含哪些文件，缺一个就打不开

SHP 从来不是单个文件，而是一组同名文件的集合。拿到压缩包解压后，先按扩展名清点，缺任何一个都会在读取时报错或静默丢数据。

扩展名	作用	缺失后果
.shp	存储几何坐标	完全无法读取
.shx	几何索引	部分库报错，或读取极慢
.dbf	属性表	几何能读，字段全空
.prj	坐标系定义	坐标被当成无投影的经纬度或平面值
.cpg	属性表编码	中文县名变乱码

我一般会先跑一段脚本做体检，而不是靠肉眼在文件夹里数。下面这段用geopandas把每个图层的要素数、几何类型、坐标系、字段名一次性打出来，省得后面反复翻。

import geopandas as gpd import os base = "./shandong_2022_07" # 解压后的目录 for fname in sorted(os.listdir(base)): if fname.endswith(".shp"): path = os.path.join(base, fname) gdf = gpd.read_file(path) print("=" * 50) print("图层:", fname) print("要素数:", len(gdf)) print("几何类型:", gdf.geom_type.unique()) print("坐标系:", gdf.crs) print("字段:", list(gdf.columns)) # 检查几何有效性，面图层尤其重要 invalid = (~gdf.geometry.is_valid).sum() print("无效几何数:", invalid)

这段代码的关键不在打印，而在最后一行。面图层里常见的自相交、悬挂节点，在is_valid为 False 时不会阻止你读取，但会在做叠加分析时让结果凭空少几块。参数上没什么可调的，gpd.read_file默认按.prj识别坐标系，如果.prj缺失，gdf.crs会是 None，这时候不要猜，先去找原始来源确认，猜错坐标系是后面所有面积、距离算错的根源。

2.2 坐标系怎么定：地理坐标系与投影坐标系别混用

省域数据最常见的坑是坐标系混用。行政区划面图层如果是 CGCS2000 地理坐标系（单位是度），而路网图层被转成了某个投影坐标系（单位是米），两者直接叠加，路网会飞到几千公里外。判断方法很简单：看gdf.crs.is_geographic，为 True 就是经纬度。

做面积统计、缓冲区、路网长度这类需要真实距离的计算，必须统一到投影坐标系。省域尺度我一般用 Albers 等积投影或高斯-克吕格对应带号，具体带号按数据覆盖的经度范围选。下面是把三级面图层和路网统一投影的写法：

import geopandas as gpd from pyproj import CRS # 省域尺度常用 Albers 等积投影，参数按覆盖范围调整 albers = CRS.from_proj4( "+proj=aea +lat_1=34 +lat_2=38 +lat_0=36 +lon_0=118 +x_0=0 +y_0=0 +datum=CGCS2000 +units=m +no_defs" ) for name in ["province", "city", "county", "road", "highway"]: gdf = gpd.read_file(f"./shandong_2022_07/{name}.shp") if gdf.crs is None: raise ValueError(f"{name} 缺少坐标系定义，先确认来源再继续") gdf_proj = gdf.to_crs(albers) gdf_proj.to_file(f"./out/{name}_aea.shp", encoding="utf-8")

lat_1、lat_2是标准纬线，取覆盖区域南北边界附近能让变形更均匀；lon_0是中央经线，取区域中心经度。这几个参数不是随便填的，填偏了面积误差会累积到统计结果里。转换后统一写回 SHP 时务必带encoding="utf-8"，否则中文属性在部分环境下会变成问号。

2.3 路网与行政区划的层级关系：别把高速和乡道当一回事

道路网和公路网两个线图层，字段结构往往不同。公路网通常带等级字段（高速、国道、省道、县道），道路网可能更细，含城市道路。做可达性分析时，如果把所有线一视同仁算平均速度，结果会严重失真。我一般先看属性表里有没有grade、type、name这类字段，没有就按图层来源分开处理。

road = gpd.read_file("./out/road_aea.shp") print(road.columns.tolist()) # 假设存在等级字段，按等级统计总长度 if "grade" in road.columns: road["length_km"] = road.geometry.length / 1000 print(road.groupby("grade")["length_km"].sum().sort_values(ascending=False))

geometry.length在投影坐标系下单位是米，除以 1000 得公里。这一步的意义在于先建立对数据的量级认知：如果某等级总长度明显异常（比如县道比乡道还短），大概率是字段值不规范或存在重复线。后面做网络分析前，这个统计能帮你快速发现数据质量问题。

3. 用 GeoPandas 把三级行政区划和路网跑通：叠加统计的最小闭环

3.1 从点位到县界：空间连接算密度

最常见的需求是：手上一批点位（门店、监测站、事件记录），想知道每个县有多少个、密度多少。核心操作是sjoin，把点挂到面上。

import geopandas as gpd import pandas as pd county = gpd.read_file("./out/county_aea.shp") points = gpd.read_file("./data/points.shp").to_crs(county.crs) # 空间连接，谓词用 within 表示点落在面内 joined = gpd.sjoin(points, county, how="left", predicate="within") # 按县统计点数和密度 county["pt_count"] = county.index.map(joined.groupby("index_right").size()).fillna(0).astype(int) county["area_km2"] = county.geometry.area / 1e6 county["density"] = county["pt_count"] / county["area_km2"] print(county[["name", "pt_count", "area_km2", "density"]].head())

predicate="within"表示点完全落在面内，落在边界上的点可能被漏掉，需要时改用intersects。how="left"保证没有点的县也保留，计数补 0，否则统计表会缺行。index_right是 sjoin 自动加的目标面索引，用它做分组映射比循环快得多。面积除以 1e6 是把平方米换成平方公里，这一步在投影坐标系下才成立，地理坐标系下算出来的“面积”是平方度，没有物理意义。

3.2 路网缓冲区覆盖县界比例：一个能直接复用的指标

判断某县路网密度，除了线长度除以面积，还可以算缓冲区覆盖率：路网两侧各缓冲一定距离，看覆盖了县面积的多少。这个指标对选址和通达性评估很直观。

county = gpd.read_file("./out/county_aea.shp") road = gpd.read_file("./out/road_aea.shp") # 路网两侧各缓冲 500 米，融合成整体 road_buf = road.geometry.buffer(500).unary_union results = [] for idx, row in county.iterrows(): inter = row.geometry.intersection(road_buf) ratio = inter.area / row.geometry.area if row.geometry.area > 0 else 0 results.append({"name": row["name"], "cover_ratio": round(ratio, 4)}) print(pd.DataFrame(results).sort_values("cover_ratio", ascending=False).head(10))

buffer(500)的单位是米，因为已经投影。unary_union把相邻缓冲区合并，避免重叠部分重复计算。逐县求交在县数量多时会慢，可以先用gpd.overlay做批量相交再分组，但逐县写法逻辑最清楚，适合先跑通再优化。覆盖率接近 1 说明路网密集，接近 0 说明该县大部分区域远离道路，这个结论比单纯看总长度更有解释力。

3.3 把结果落成可交付的表和图层

分析跑完，交付物一般两类：一张带统计字段的面图层，和一份纯表格。前者给 GIS 用，后者给业务方看。

# 写回带统计字段的面图层 county.to_file("./out/county_stats.shp", encoding="utf-8") # 导出纯表格，去掉几何列 county.drop(columns="geometry").to_csv("./out/county_stats.csv", index=False, encoding="utf-8-sig")

encoding="utf-8-sig"是为了让表格软件打开 CSV 时中文不乱码，这个细节踩过一次就记住了。SHP 字段名有长度限制（常见 10 字符），统计字段名尽量短，比如pt_count、density，太长会被截断导致后续读取对不上。

4. 避坑与排查：SHP 处理里最容易翻车的五件事

4.1 中文属性变乱码：现象是县名显示为问号或方块

现象：打开属性表，中文全是乱码。原因：.dbf的编码没有声明，或.cpg文件缺失，读取库按默认编码解析。解决：读取时显式指定encoding="gbk"或"utf-8"试一次，确认后补一个.cpg文件写明编码，再统一转存为 UTF-8。我一般会在第一次读取时把encoding参数试两遍，看哪个能正确显示县名，再固定下来。

4.2 面积算出来大得离谱：现象是县面积是几万平方公里

现象：一个县算出几万平方公里，明显不对。原因：在地理坐标系下直接调geometry.area，得到的是平方度。解决：先to_crs到投影坐标系再算面积。这个坑几乎每个新手都会踩一次，判断方法是看数值量级，正常县面积在几百到几千平方公里，出现几万、几十万就要怀疑坐标系。

4.3 空间连接后点数对不上：现象是总数比原始点位少

现象：sjoin 之后统计的总点数少于原始点位数。原因：有点落在县界外，或落在边界线上被within排除，或点位坐标系与面图层不一致。解决：先检查points.crs == county.crs，不一致先to_crs；再用predicate="intersects"对比结果差异；最后统计joined["index_right"].isna().sum()看有多少点没匹配上，逐个排查。

4.4 路网长度统计重复：现象是总里程明显偏大

现象：路网总长度比预期大很多。原因：线图层存在重复要素，或一条路被拆成多段且属性重复。解决：先按几何去重，road.drop_duplicates(subset="geometry")，再按名称和等级检查是否有重叠段。去重前先备份原始文件，去重是不可逆操作。

4.5 写回 SHP 字段被截断：现象是统计字段读出来是空或错位

现象：写回的 SHP 里，长字段名变成截断后的名字，或者字段值错位。原因：SHP 格式对字段名长度有限制，超出部分被截。解决：字段名控制在 10 字符以内，写回前用gdf.columns检查一遍，必要时先重命名。需要保留长字段名时，改用 GeoPackage 或 GeoJSON 格式交付。

5. 进阶：把这份 2022.07 数据做成可复用的分析底座

5.1 建立统一的投影与字段规范

一份数据用一次就扔，是浪费。我习惯在第一次处理时就把规范定下来：所有图层统一到同一个 Albers 投影，字段名统一小写下划线风格，面图层保留name、code、level三个基础字段，路网保留grade、length_km。规范定好后写成一个配置字典，后续脚本都从配置读参数，换一批数据只改配置不改代码。

CONFIG = { "crs": "+proj=aea +lat_1=34 +lat_2=38 +lat_0=36 +lon_0=118 +datum=CGCS2000 +units=m +no_defs", "layers": { "county": {"name_field": "name", "code_field": "code"}, "road": {"grade_field": "grade"}, }, "buffer_m": 500, }

这样做的价值在于可复现：半年后换一份新时点的数据，跑同一套脚本就能得到可比结果，而不是每次重新拍脑袋定参数。

5.2 用分区统计做时间对比的基线

2022.07 这个时点本身就是一个基线。把县界的面积、路网长度、覆盖率、点密度存成一张宽表，后续拿到新时点数据，按code字段做左连接对比，就能算出变化量。对比时注意行政区划调整：如果新数据里某些县被合并或拆分，直接按名称连接会错位，必须用稳定的行政区划代码做键。这一步没有捷径，只能先核对两期的代码表。

对比项	基线字段	新时点字段	变化计算
面积	area_km2	area_km2_new	差值，异常时查边界调整
路网长度	road_km	road_km_new	差值，注意等级口径一致
覆盖率	cover_ratio	cover_ratio_new	差值，缓冲区半径须相同
点密度	density	density_new	比值，注意点位口径