简介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坐标系。县级边界的要素数一般在一百上下市级在十个左右如果县级数量明显偏少可能是只导出了部分区域也可能是多部件几何被合并了——这层要单独用-geomSUMMARY之外的选项看清楚。范围这一行最直观陕西跨东经 105° 到 111°、北纬 31° 到 39°如果Extent显示成(0, 0) - (1, 1)或者(-20037508, ...)这种数量级坐标系基本可以判定有问题。前者是没定义后者是被误当成 Web 墨卡托了。QGIS 里也有等价操作右键图层 → 属性 → 信息看“坐标参照系”那一栏是不是 EPSG 编码如果是“未知”或者用户自定义就得手动指定。批量看的话Python 更快import glob, os import geopandas as gpd # recursiveTrue 让 glob 钻进子目录读不了的文件直接跳过不影响整体盘点 for f in sorted(glob.glob(shaanxi_2022_07/**/*.shp, recursiveTrue)): 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遇到中文属性乱码时加encodinggbk实在读不出来加enginepyogrio换条路试试。2.3 属性字段决定你能算到什么粒度拿到数据第一件事是把字段列表打出来看有没有能唯一标识行政单元的编码。国内行政边界数据常见的字段是六位行政区划代码省级以 61 开头、后四位为 0市级后两位为 0县级六位全用。有这个字段分组统计就稳只有名称字段遇到重名或者“市辖区”这类条目就得手工补。import geopandas as gpd city gpd.read_file(shaanxi_2022_07/adm/sx_adm_city.shp, encodinggbk) 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 三种坐标系的适用边界边界和路网这类数据绕不开三种坐标系用错场景的后果完全不同坐标系EPSG单位适用场景常见误用CGCS2000 地理坐标4490度数据交换、入库、叠加在线底图前直接算长度和面积CGCS2000 3 度带高斯投影4544 / 4545 / 4546米省内长度、面积、缓冲区分析跨带混用、中央经线选错Web 墨卡托3857米切片出图、网页展示拿来做面积统计陕西东西跨东经 105° 到 111°正好落在 3 度带的三个带里105°E、108°E、111°E对应 EPSG:4544、4545、4546。全省级统计必须统一到一个带上一般取 108°EEPSG: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 ENCODINGUTF-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 ENCODINGUTF-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_okTrue) for layer, path in jobs.items(): g gpd.read_file(path, encodinggbk) if g.crs is None: # prj 缺失时按来源显式指定 g g.set_crs(4490, allow_overrideTrue) g g.to_crs(TARGET) # 统一到 108E 中央经线 g.to_file(out/sx_4545.gpkg, layerlayer, driverGPKG) 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, layeradm_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, layeradm_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, layerroad) 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 joinsjoin做的是要素归属判断一条线只要和某个市的面相交整条线都归到这个市。跨界的路会同时归给两个市两边各算一遍全长。overlayintersection做的是几何切割先把线在面边界处切断再按段归属各算各的那一段。从统计口径上说overlay是对的——里程应该按实际落在辖区内的长度算。但实际项目里两种都要会做“本市涉及哪些道路”的清单用 sjoin 更快做里程汇总必须用 overlay。import geopandas as gpd import pandas as pd city gpd.read_file(out/sx_4545.gpkg, layeradm_city) road gpd.read_file(out/sx_4545.gpkg, layerroad) # 只保留需要的字段overlay 对字段多的大表很吃内存 city city[[adcode, name, geometry]] road road[[grade, name, geometry]] road[grade] road[grade].fillna(未分级) # 按面切割线跨界的线被切成多段每段归属各自的面 seg gpd.overlay(road, city, howintersection, keep_geom_typeTrue) # 切割后重新算长度单位是米因为已经是 4545 seg[len_km] seg.length / 1000.0 stat (seg.groupby([adcode, name, grade])[len_km] .sum() .unstack(fill_value0) .round(1)) print(stat.head()) stat.to_csv(out/road_len_by_city.csv, encodingutf-8-sig)几个参数值得说清楚。keep_geom_typeTrue保证输出还是线不加的话面线相交可能生出点或者空几何后面的length就没意义。howintersection是取交集还有一种用法是howdifference取差集用来找“不在任何市界内的道路”——这招查数据错位特别好使。unstack(fill_value0)把等级展开成列方便直接做对比表。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)], crsroad.crs) cover gpd.overlay( gpd.GeoDataFrame(geometryroad_buf, crsroad.crs), city, howintersection ) cover[cov_km2] cover.length * 0 # 占位下面用面积重算 cover[cov_km2] cover.area / 1e6 cover[city_km2] cover.apply(lambda r: 0, axis1) 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(upper1).round(4) print(res.sort_values(覆盖率, ascendingFalse).head(10))这里有个必须做的操作先把所有道路融合成一张网再 buffer。逐条 buffer 会得到一堆重叠圆角面积重复计算覆盖率能算出大于 1 的荒唐值。union_all()是新版 GeoPandas 的方法老版本用unary_union。clip(upper1)是个兜底正常情况下不应该触发触发了说明融合没做对。提示缓冲区半径越大圆角分段越多内存涨得越快。半径超过 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_CSGDAL 无法识别读取后crs为NoneQGIS 默认按未知坐标系处理叠加时不做任何变换。解决先g.crs为None就set_crs(4490)显式声明再to_crs转换。同时检查.prj文件内容是不是一段完整的 WKT几百字节以上才算正常几十字节的多半是残文件。5.2 属性表中文全是问号或乱码现象QGIS 打开后名称字段显示成????或者一堆生僻字编码换了几次都不对。原因SHP 的.dbf默认用系统代码页中文环境常见是 GBK但文件里没有.cpg声明软件按 UTF-8 读就乱了。解决read_file(..., encodinggbk)指定编码QGIS 里在图层属性的“数据源编码”手动选 GBK。彻底的办法是转成 GPKG 时用-lco ENCODINGUTF-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:hostlocalhost dbnamegis usergis \ out/sx_4545.gpkg \ -lco GEOMETRY_NAMEgeom \ -lco FIDgid \ -lco PRECISIONNO \ -nlt PROMOTE_TO_MULTIGEOMETRY_NAMEgeom让字段名统一后面写 SQL 不用记每个表叫wkb_geometry还是the_geomPRECISIONNO关掉坐标精度约束避免导入时因为小数位超限报错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都跑过了再动数据。希望帮到你。本文还有配套的精品资源点击获取