简介：2022年7月版陕西省行政区划与路网SHP数据包，面向GIS专业学生、空间数据分析人员以及交通与城乡规划研究者，提供省、市、县三级行政边界和道路、铁路网络的矢量基础数据，可用于行政区划制图、交通可达性分析、物流路径规划及灾害风险评估等场景。压缩包共33个文件，约25.94MB，以shp几何数据、dbf属性表、prj坐标参考、shx索引及xml元数据等为主，配合sbn/sbx空间索引，便于在ArcGIS、QGIS等平台直接加载与叠加分析。目前已有2063人学习下载，说明其作为区域基础地理数据的参考价值较受关注。读者可获得较完整的陕西行政边界与路网图层，并借助属性表开展道路等级、铁路线路与行政区划关联分析，为人口分布、道路密度、城市间连接性等研究提供可直接使用的数据底图与练习素材。

1. 陕西这份 2022.07 的 SHP：边界与路网叠不上，问题多半不在数据

有位做货运调度的 A 同学，拿到一份 2022.07 时点的陕西 SHP 包：省、市、县三级行政边界面齐全，另带道路网和公路网两层线。他打开 QGIS 想把路网按地级市切分、统计高速和国省道里程，结果两层一叠，边界跑到赤道附近去了，路网还留在原地。他第一反应是数据错了，其实九成是坐标系没定义或定义错了。这套数据能解决区域统计、选址、路网密度评估、汇报底图这几类活，适合做城乡规划、物流网络、交通分析的从业者，也适合要交课程作业的学生。陕西行政边界、道路网、铁路网这几个词在 GIS 检索里经常一起出现，真正卡人的从来不是有没有数据，而是拿到手之后参数对不对。

2. 解包与图层盘点：用 ogrinfo 十分钟摸清字段和坐标系

先把压缩包解开，别急着往 QGIS 里拖。SHP 是组合格式，.shp存几何、.dbf存属性、.shx是索引、.prj是坐标系定义，四个文件缺一个都不完整——尤其是.prj，丢了它，软件只能靠猜，猜错就直接把边界甩到赤道。这一步的目标很明确：搞清楚有几个图层、每层多少要素、坐标系是什么、属性里有什么字段能拿来分组统计。

2.1 目录结构与命名规律

常见的组织方式是按用途分目录，行政边界、道路、公路各放一处，文件名里带省份缩写和层级后缀。层级一般用 prov / city / county 或者 省 / 市 / 县 区分，行政级别也可能直接写进属性字段。道路网与公路网这两层容易混：道路网通常覆盖面更宽，城市道路、乡镇道路、村道都收；公路网偏重国省县道这类带行政等级的线路，属性里有等级编码。两者要素数不一样，交叉在地图上看着重合，实际并不完全一样。

先用命令行把家底列出来，命令比在资源管理器里翻快得多：

# 列出目录下所有文件及大小，先确认 shp/dbf/shx/prj 是否成套 find ./shaanxi_2022_07 -type f -printf '%10s %p\n' | sort -k2 # 重点看每个 shp 旁边有没有同名 prj，缺了就得自己补坐标系 find ./shaanxi_2022_07 -name "*.shp" -exec sh -c 'test -f "${1%.shp}.prj" || echo "缺 prj: $1"' _ {} \;

-printf是 GNU find 的写法，macOS 上换成ls -lhR。第二条命令的逻辑是：对每个.shp，拼出同名.prj路径，存在就跳过，不存在就打印。先跑这两条，后面少踩一半坑。

如果一层里混了铁路和公路，属性表通常会有一个等级或类型字段把它们分开；没有这个字段，说明这层只覆盖公路，铁路要另找来源，别硬套。

2.2 用 ogrinfo 核对要素数与坐标系

GDAL 自带的ogrinfo是查 SHP 元信息最省事的工具，-so只输出摘要不打印全部要素，-al遍历所有图层：

# -so: summary only -al: all layers ogrinfo -so -al ./shaanxi_2022_07/adm/sx_adm_county.shp

输出里重点看四行：Geometry（面还是线）、Feature Count（要素数）、Extent（范围）、Layer SRS WKT（坐标系）。县级边界的要素数一般在一百上下，市级在十个左右，如果县级数量明显偏少，可能是只导出了部分区域，也可能是多部件几何被合并了——这层要单独用-geom=SUMMARY之外的选项看清楚。

范围这一行最直观：陕西跨东经 105° 到 111°、北纬 31° 到 39°，如果Extent显示成(0, 0) - (1, 1)或者(-20037508, ...)这种数量级，坐标系基本可以判定有问题。前者是没定义，后者是被误当成 Web 墨卡托了。

QGIS 里也有等价操作：右键图层 → 属性 → 信息，看“坐标参照系”那一栏是不是 EPSG 编码；如果是“未知”或者用户自定义，就得手动指定。批量看的话，Python 更快：

import glob, os import geopandas as gpd # recursive=True 让 glob 钻进子目录；读不了的文件直接跳过，不影响整体盘点 for f in sorted(glob.glob("shaanxi_2022_07/**/*.shp", recursive=True)): try: g = gpd.read_file(f) except Exception as e: print("读取失败", f, e) continue print(f"{os.path.basename(f):28s} 要素={len(g):6d} 类型={set(g.geom_type)} CRS={g.crs}")

geopandas.read_file底层就是 GDAL，.prj缺失时g.crs返回None，这时候别急着to_crs，先按来源把原始坐标系显式设上（g.set_crs(4490)），再做转换。

提示：read_file遇到中文属性乱码时加encoding="gbk"；实在读不出来，加engine="pyogrio"换条路试试。

2.3 属性字段决定你能算到什么粒度

拿到数据第一件事，是把字段列表打出来，看有没有能唯一标识行政单元的编码。国内行政边界数据常见的字段是六位行政区划代码：省级以 61 开头、后四位为 0，市级后两位为 0，县级六位全用。有这个字段，分组统计就稳；只有名称字段，遇到重名或者“市辖区”这类条目就得手工补。

import geopandas as gpd city = gpd.read_file("shaanxi_2022_07/adm/sx_adm_city.shp", encoding="gbk") print(city.columns.tolist()) # 看有哪些字段 print(city.dtypes) # 看编码是 int 还是 object print(city.head(3).T) # 转置看三行，比横向打印好读 # 检查编码位数是否统一，六位才是完整的县级/市级代码 if "adcode" in city.columns: print(city["adcode"].astype(str).str.len().value_counts())

dtypes里如果看到adcode是int64，要警惕前导零丢失——这种情况不会出现在省级（61 开头），但县级代码里如果某段以 0 打头就会出问题，统一astype(str).str.zfill(6)处理。

路网层的字段更关键。至少要有一个等级字段，能区分高速、国道、省道、县道、乡道；有的数据用文本，有的用数字编码加一张对照表。字段名常见为grade、type、road_class、DLDJ之类，不要靠猜，直接把唯一值打出来看：

road = gpd.read_file("shaanxi_2022_07/road/sx_road.shp") for col in road.columns: if road[col].dtype == object and road[col].nunique() < 30: print(col, road[col].value_counts().to_dict())

nunique() < 30是个经验阈值：唯一值少的字段大概率是分类字段，多的多半是名称或编号。把每个候选列的唯一值都扫一遍，等级体系就清楚了。顺带检查有没有length之类已经算好的里程字段——如果有，注意它的单位和计算时用的坐标系，投影坐标下算出来的米才可信，地理坐标下算出来的是度，不能直接乘系数。

这一步做完，你手里应该有一张清单：几个图层、每层多少要素、坐标系是什么、按哪个字段分组、等级怎么分。后面所有操作都建立在这张清单上。

3. 坐标与拓扑：4490、4544、3857 到底该选哪个

坐标系选错，后面所有统计都是错的，而且错得很隐蔽——图能画出来，数也在，就是不对。这一章把选型理由讲清楚，再给可执行的重投影与拓扑检查脚本。

3.1 三种坐标系的适用边界

边界和路网这类数据，绕不开三种坐标系，用错场景的后果完全不同：

陕西东西跨东经 105° 到 111°，正好落在 3 度带的三个带里：105°E、108°E、111°E，对应 EPSG:4544、4545、4546。全省级统计必须统一到一个带上，一般取 108°E（EPSG:4545），因为中央经线在省域中部，边缘变形相对可控。如果只做某一个县的分析，按它所在的带选更准。

Web 墨卡托（3857）只适合“看着对”的场景。它在中纬度会把东西方向拉伸，纬度越高拉伸越明显，用它算出来的陕西面积和实际差得远，路网里程同样偏大。图上好看和数上准确，是两件必须分开处理的事。

注意：.prj缺失时不要凭“看着像经纬度”就断言是 4490。先看数值范围，再看数据来源说明，实在没依据就跟提供方确认，这一步没有后悔药。

3.2 批量重投影：ogr2ogr 与 GeoPandas 两条路

数据量不大的话，用ogr2ogr一条命令转一层，还能顺手把 SHP 换成 GeoPackage——GPKG 单文件、字段名不截断到 10 字符、支持 UTF-8，比 SHP 舒服太多：

# -t_srs 目标坐标系 -nln 输出图层名 -lco ENCODING 属性表编码 ogr2ogr -f GPKG out/sx_4545.gpkg shaanxi_2022_07/adm/sx_adm_city.shp \ -t_srs EPSG:4545 \ -nln adm_city_4545 \ -lco ENCODING=UTF-8 \ -nlt MULTIPOLYGON ogr2ogr -f GPKG out/sx_4545.gpkg shaanxi_2022_07/road/sx_road.shp \ -t_srs EPSG:4545 -nln road_4545 -append -lco ENCODING=UTF-8

第一条建库并写入市级边界，-nlt MULTIPOLYGON强制统一几何类型，避免单部件和多部件混存导致后续操作报错。第二条用-append追加到同一个 GPKG 里，省得管理一堆散文件。两条命令的顺序不能反：-append的前提是库已存在。

用 Python 做批量转换更灵活，能顺手做类型清理：

import pathlib import geopandas as gpd TARGET = 4545 jobs = { "adm_prov": "shaanxi_2022_07/adm/sx_adm_prov.shp", "adm_city": "shaanxi_2022_07/adm/sx_adm_city.shp", "adm_county": "shaanxi_2022_07/adm/sx_adm_county.shp", "road": "shaanxi_2022_07/road/sx_road.shp", "highway": "shaanxi_2022_07/highway/sx_highway.shp", } pathlib.Path("out").mkdir(exist_ok=True) for layer, path in jobs.items(): g = gpd.read_file(path, encoding="gbk") if g.crs is None: # prj 缺失时按来源显式指定 g = g.set_crs(4490, allow_override=True) g = g.to_crs(TARGET) # 统一到 108E 中央经线 g.to_file("out/sx_4545.gpkg", layer=layer, driver="GPKG") print(layer, len(g), "→", g.crs.to_epsg())

关键是set_crs和to_crs的区别：前者只贴标签不改坐标数值，后者做真正的数学变换。很多人把没坐标系的图层直接to_crs(4545)，GeoPandas 会当成已经是 4490 处理，结果整体偏移几百公里——这就是开头 A 同学遇到的现象。判断依据很简单：.prj缺失、或者g.crs is None，就先set_crs。

3.3 拓扑检查：面要严丝合缝，线要通到底

行政边界的常见毛病是接缝：相邻县的共用边一个顶点位置有微小差异，导致叠在一起时有缝隙或碎屑面。直接后果是面积统计加起来比全省面积多或者少。检查办法是用shapely做有效性校验：

import geopandas as gpd from shapely.validation import make_valid county = gpd.read_file("out/sx_4545.gpkg", layer="adm_county") bad = ~county.geometry.is_valid print("无效几何：", bad.sum()) # 自相交、零面积环等用 make_valid 修复，别用 buffer(0)，后者会丢部件 county.loc[bad, "geometry"] = county.loc[bad, "geometry"].apply(make_valid) # 面积总和与省级面对比，差值超过千分之一就说明有缝隙或重叠 prov = gpd.read_file("out/sx_4545.gpkg", layer="adm_prov") diff = abs(county.area.sum() - prov.area.sum()) / prov.area.sum() print(f"面积偏差 {diff:.4%}")

偏差在千分之一以内属于正常（不同来源的边界精度差异），超过百分之一就要查是不是有重叠面。重叠用sjoin找交集非空的要素对，缝隙则比较麻烦，一般用省级面减去县集合的差集来看残余碎片。

路网层的拓扑问题不同：线要连通，不能断在半路。检查悬挂节点的思路是统计每个端点的出现次数，只出现一次的就是悬挂点。用shapely的unary_union加节点计数能粗略查出来；工程上更省事的做法是直接看gaps——把线按缓冲区融合成面，如果一条主干道中间断开，融合后会出现孤立的碎片面。对做路网密度统计来说，悬挂节点影响不大；对做路径分析或者连通性评估，这就是致命的。

还有一个容易被忽略的点：重复要素。同一段路在道路网和公路网里各出现一次，或者同一条线正反两个方向各存一条，里程统计直接翻倍。去重的办法是按几何 WKB 加等级字段做联合唯一：

road = gpd.read_file("out/sx_4545.gpkg", layer="road") before = len(road) road["_wkb"] = road.geometry.apply(lambda g: g.wkb_hex) road = road.drop_duplicates(subset=["_wkb"]).drop(columns=["_wkb"]) print(f"去重 {before} → {len(road)}")

注意这里只按几何去重，如果同名同几何但等级不同的两条记录确实代表不同信息（比如复合等级），就得改成按几何加等级两个字段一起去重。

4. 按市统计公路里程：overlay 与 sjoin 差出来的那 30%

这一章把数据用起来。最常见也最容易出错的诉求是“按市统计各类公路里程”，里面藏着一个方法论问题：跨市界的路，算给谁。两种做法给出两个数，差 30% 很平常。

4.1 面选线：sjoin 与 overlay 的本质区别

spatial join（sjoin）做的是要素归属判断：一条线只要和某个市的面相交，整条线都归到这个市。跨界的路会同时归给两个市，两边各算一遍全长。overlay（intersection）做的是几何切割：先把线在面边界处切断，再按段归属，各算各的那一段。

从统计口径上说，overlay是对的——里程应该按实际落在辖区内的长度算。但实际项目里两种都要会：做“本市涉及哪些道路”的清单，用 sjoin 更快；做里程汇总，必须用 overlay。

import geopandas as gpd import pandas as pd city = gpd.read_file("out/sx_4545.gpkg", layer="adm_city") road = gpd.read_file("out/sx_4545.gpkg", layer="road") # 只保留需要的字段，overlay 对字段多的大表很吃内存 city = city[["adcode", "name", "geometry"]] road = road[["grade", "name", "geometry"]] road["grade"] = road["grade"].fillna("未分级") # 按面切割线，跨界的线被切成多段，每段归属各自的面 seg = gpd.overlay(road, city, how="intersection", keep_geom_type=True) # 切割后重新算长度，单位是米（因为已经是 4545） seg["len_km"] = seg.length / 1000.0 stat = (seg.groupby(["adcode", "name", "grade"])["len_km"] .sum() .unstack(fill_value=0) .round(1)) print(stat.head()) stat.to_csv("out/road_len_by_city.csv", encoding="utf-8-sig")

几个参数值得说清楚。keep_geom_type=True保证输出还是线，不加的话面线相交可能生出点或者空几何，后面的length就没意义。how="intersection"是取交集，还有一种用法是how="difference"取差集，用来找“不在任何市界内的道路”——这招查数据错位特别好使。unstack(fill_value=0)把等级展开成列，方便直接做对比表。

len_km必须在切割之后算。如果先算长度再 overlay，切出来的段还带着原来整条线的长度字段，加总就重复了。这个坑很隐蔽，因为数字看着挺合理。

4.2 缓冲区与覆盖率：把米制单位当回事

路网覆盖率的常见算法是：以道路为中心做 N 公里缓冲区，算缓冲区面积占辖区面积的比例。核心前提是缓冲区半径和坐标单位一致，投影坐标系下半径写 2000 就是两公里，地理坐标系下写 2000 会被当成 2000 度，结果直接是空。

# 2000 米 = 2 公里，前提是 city 已在 EPSG:4545 city_buf = city.copy() city_buf["geometry"] = city.geometry.buffer(2000) # 先融合成一张网再算覆盖，避免重复覆盖导致比例超过 1 road_union = road.union_all() if hasattr(road, "union_all") else road.unary_union road_buf = gpd.GeoSeries([road_union.buffer(2000)], crs=road.crs) cover = gpd.overlay( gpd.GeoDataFrame(geometry=road_buf, crs=road.crs), city, how="intersection" ) cover["cov_km2"] = cover.length * 0 # 占位，下面用面积重算 cover["cov_km2"] = cover.area / 1e6 cover["city_km2"] = cover.apply(lambda r: 0, axis=1) res = cover.groupby("name")["cov_km2"].sum().to_frame() res["city_km2"] = city.set_index("name").area / 1e6 res["覆盖率"] = (res["cov_km2"] / res["city_km2"]).clip(upper=1).round(4) print(res.sort_values("覆盖率", ascending=False).head(10))

这里有个必须做的操作：先把所有道路融合成一张网再 buffer。逐条 buffer 会得到一堆重叠圆角，面积重复计算，覆盖率能算出大于 1 的荒唐值。union_all()是新版 GeoPandas 的方法，老版本用unary_union。clip(upper=1)是个兜底，正常情况下不应该触发，触发了说明融合没做对。

提示：缓冲区半径越大，圆角分段越多，内存涨得越快。半径超过 5 公里、要素超过十万条时，考虑先按市分组再算，别一次全上。

4.3 出图与样式：让边界和路网各司其职

出图这件事，QGIS 比代码省事。加载 GPKG 后按这个顺序配样式：行政边界用无填充加细描边，市级加粗、县级减细；道路网按等级分色，高速最粗最醒目，县道细线灰度处理；最后叠一层浅色底图压住视觉重心。

分级配色建议用色带而不是离散色，等级是有序的，色带的明度变化能直接传达这个信息。线宽按等级设成 0.2 / 0.4 / 0.8 / 1.2 毫米四档，超过四档人眼分不清。标注只开市级名称，放在面内部居中，县级名称在 1:50 万以下比例尺再出现。

导出时注意：如果目标是网页展示，先在 QGIS 里把 4545 转成 3857 再出图；如果目标是打印，保持 4545 直接用，别转。两种产物分开存，别指望一份文件两个用途。

5. 避坑与排查：五个把边界和路网用废的现场

下面这几条都是实际翻过车的，按“现象 → 原因 → 解决”写，遇到对得上的直接照做。

5.1 图层叠加整体偏移几百公里

现象：边界和路网加载后不在一个位置，边界跑到赤道附近或者经度 0 附近，路网留在原地。

原因：.prj缺失，或者存在但内容是空的、写成LOCAL_CS，GDAL 无法识别，读取后crs为None，QGIS 默认按未知坐标系处理，叠加时不做任何变换。

解决：先g.crs，为None就set_crs(4490)显式声明，再to_crs转换。同时检查.prj文件内容是不是一段完整的 WKT，几百字节以上才算正常，几十字节的多半是残文件。

5.2 属性表中文全是问号或乱码

现象：QGIS 打开后名称字段显示成????或者一堆生僻字，编码换了几次都不对。

原因：SHP 的.dbf默认用系统代码页，中文环境常见是 GBK，但文件里没有.cpg声明，软件按 UTF-8 读就乱了。

解决：read_file(..., encoding="gbk")指定编码；QGIS 里在图层属性的“数据源编码”手动选 GBK。彻底的办法是转成 GPKG 时用-lco ENCODING=UTF-8，一次性把编码钉死，后面不再折腾。

5.3 各市面积加起来比全省多

现象：按市汇总面积，总和比省级面大 3% 到 8%。

原因：相邻市界不重合，存在重叠带；或者市级面是从不同来源拼接的，边界精度不一致。

解决：用省级面做一次“权威裁剪”——以省面为掩膜对市级面做 intersection，强制让市级集合不超出省界；重叠部分用difference减去，或者按面积比例分摊。检查手段就是 3.3 节那个面积偏差公式，超过 1% 就得处理。

5.4 里程统计出现整数倍偏差

现象：某几个市的公路里程明显偏高，接近实际的两倍。

原因：重复要素。同一条路在道路网和公路网各存一份，或者数据导出时正反方向各生成一条。

解决：按几何 WKB 去重（见 3.3 节），再去比对两个图层之间的重合段。如果两层本就应该有交集，统计时只取其中一层，别叠加。

5.5 字段名被截断成 st_na_1

现象：转成 SHP 之后属性字段名全变成 10 个字符以内，road_name变成road_name（勉强），highway_grade变成highway_gr，两个字段截断后重名直接覆盖。

原因：DBF 格式的字段名上限就是 10 字节，中文占两字节，五个汉字就到顶了。

解决：中间产物一律用 GeoPackage，只在最终交付必须给 SHP 时才转，转之前把字段名改成短英文缩写，并附一份字段对照表。这个坑没有技术解法，只能靠流程规避。

6. 进阶：把这份 SHP 变成可复用的底图服务

一次性的分析做完就完了，可惜。这套数据真正的价值在于变成底图，后面每个项目都能直接调。我的做法是把它灌进 PostGIS，建成带空间索引的表，再用视图把常用的统计口径固化下来。

# 一次性把 GPKG 里的图层全部推进 PostGIS ogr2ogr -f PostgreSQL "PG:host=localhost dbname=gis user=gis" \ out/sx_4545.gpkg \ -lco GEOMETRY_NAME=geom \ -lco FID=gid \ -lco PRECISION=NO \ -nlt PROMOTE_TO_MULTI

GEOMETRY_NAME=geom让字段名统一，后面写 SQL 不用记每个表叫wkb_geometry还是the_geom；PRECISION=NO关掉坐标精度约束，避免导入时因为小数位超限报错；PROMOTE_TO_MULTI把单部件几何升级成多部件，混合类型入库不会中断。

进库之后建索引和视图：

-- 空间索引是必须的，没索引的 ST_Intersects 在大表上会跑到天荒地老 CREATE INDEX idx_city_geom ON adm_city_4545 USING GIST (geom); CREATE INDEX idx_road_geom ON road_4545 USING GIST (geom); ANALYZE adm_city_4545; ANALYZE road_4545; -- 把“按市按等级统计里程”固化成一个视图，以后直接 select CREATE OR REPLACE VIEW v_road_len_by_city AS SELECT c.adcode, c.name AS city_name, COALESCE(r.grade, '未分级') AS grade, ROUND(SUM(ST_Length(ST_Intersection(r.geom, c.geom)))::numeric / 1000, 2) AS len_km FROM adm_city_4545 c JOIN road_4545 r ON ST_Intersects(r.geom, c.geom) GROUP BY c.adcode, c.name, r.grade;

视图里用ST_Intersection而不是直接ST_Length(r.geom)，跟第 4 章的道理一样：先切再量。ST_Intersects只用来走索引筛候选，真正的计算交给ST_Intersection，这是 PostGIS 里最常用的组合写法。视图建好之后，前端或者报表工具直接查这个视图就行，不用每次重写 SQL。

还有一层可以补：把 4490 的原始版本也保留一份，别只存 4545。做 Web 展示、和第三方底图对接时，4490 更省事；做统计计算时用 4545。两套并存，靠命名区分（_4490/_4545后缀），谁也别覆盖谁。

血泪经验是：第一次入库千万别图省事跳过索引，我在一张十万级的道路表上试过一次全表ST_Intersects，跑了二十多分钟才出结果，加了 GIST 索引之后三秒。从那以后我每次建完表，第一件事就是确认索引和ANALYZE都跑过了，再动数据。希望帮到你。

本文还有配套的精品资源，点击获取