资讯详情

资讯详情

江苏GIS SHP数据实战:行政区划路网坐标系与密度出图

简介面向从事GIS空间分析与地图制图的用户这份江苏省基础地理数据包整合了省级、地市级、区县级三级行政区划边界以及道路网与铁路网线状要素适合用于区域分析、专题制图、交通可达性评估和空间统计等场景。压缩包共33个文件约41.72MB以SHP矢量格式为主体配套SHX空间索引、DBF属性表、PRJ坐标系定义及XML元数据等文件便于在ArcGIS、QGIS等常用软件中直接加载与二次处理。已有1973人学习下载说明其在江苏相关研究和课程实践中具有一定参考价值。数据按省、市、县及道路、铁路分图层组织边界与路网可叠加使用既能支撑行政区划统计、路网密度与连通性分析也方便读者快速搭建研究底图、核对空间范围并开展制图表达。1. 拿到这份江苏 GIS 数据包先别急着拖进软件三级面加两类线能干的活比想象中多上周帮一个做物流选址的同行看数据他手上有一份 2022.07 版的江苏行政区划 SHP省、市、县三级面数据外加道路网和公路网两套线数据结果卡在第一步算某个县的公路密度面积跑出来是 0.00x 平方度除下来的数完全没法看。这不是数据坏了是口径没对齐——SHP 里的几何默认是地理坐标单位是度直接拿去算面积和长度等于白算。这份资源真正的价值在于面线配套齐全三级行政单元能当统计口径用道路网和公路网能算密度、做可达性、叠底图出图省掉了自己拼边界、找路网的功夫。适合做区域分析、规划制图、GIS 二次开发的从业者也适合刚上手 QGIS、想找一份结构规整的数据练手的人。下面按我拆包的顺序写清楚。2. 拆包先看伴生文件.shp/.shx/.dbf/.prj/.cpg 少一个都可能白干解压之后第一件事不是双击打开是把目录列一遍。SHP 这个格式年头的确久久到很多人以为它就是一个文件实际上它是一组同名文件的组合体缺了谁都会在后面的某一步给你脸色看。血泪经验是从压缩包里挑着复制 .shp 出来的操作十次有八次会在投影或中文上翻车。2.1 五个必须凑齐的伴生文件各自负责什么文件作用缺失后的典型表现.shp存几何坐标主体直接打不开.shx几何索引按要素序号定位部分软件能读但属性与几何会错位或读不全.dbf属性表dBase 格式只剩下图形没有任何名称和代码.prj坐标系定义WKT 文本CRS 显示为 Unknown投影转换全错.cpg属性表字符编码声明中文变问号或乱码除了这五个还可能出现 .sbn/.sbx空间索引、.qixQGIS 空间索引、.lyr/.qmd样式文件、.xml元数据。这些属于可选项缺了不影响读取只影响打开速度或者样式还原。真正不能缺的是上表五个其中 .cpg 最容易被忽略因为它通常只有一行文本内容就是UTF-8或GBK也有写成936的肉眼看上去毫无存在感但它决定了属性表里的县名、道路名能不能正常显示。检查方式不用开图形界面命令行一条就够# 列目录重点确认 .prj 和 .cpg 是否在同级 ls -lh js_admin_2022/ # -so 只输出摘要不遍历要素-al 输出所有图层的完整信息 ogrinfo -so -al js_admin_2022/js_county.shp-so这个参数很关键。县级面数据要素数一般上百直接ogrinfo不带-so会把每个要素的几何都打印一遍终端滚半天。加上-so之后输出里会有 Feature Count、Geometry、Extent、Layer SRS WKT 和字段定义一眼就能判断投影在不在、字段有没有被截断。如果输出里Layer SRS WKT那一段是空的说明 .prj 没跟过来先别往下走。2.2 属性表怎么读行政代码、名称、道路等级字段的实际口径三级面数据的字段结构通常比较规整省、市、县各一层字段名大同小异。需要留意的是行政代码的位数省级两位、地级市四位、县级六位这是 GB/T 2260 的口径。很多分析场景下你要用这个代码做关联键而不是用名称——名称可能有市辖区城区这类重复值代码不会。道路网的字段就杂一些。常见的字段包括道路名称、等级或类型、长度等级字段名可能是ROAD_TYPE、GRADE、TYPE、CLASS中的任意一个取值可能是中文高速公路、国道、省道、县道、乡道也可能是数字编码或者字母前缀G、S。这两套口径混在一份数据里也不奇怪清洗之前先看一眼值分布再动手。还有一个必须提前知道的限制DBF 的字段名上限是 10 个字节。COUNTY_NAME是 11 个字符存进去会被截成COUNTY_NAM如果有两个字段前 10 个字符相同就会出现重名读进来之后 GeoPandas 会给后者加后缀。看字段列表时发现名字怪怪的八成就是这个原因不是数据制作者的锅。import os import geopandas as gpd # 在 read_file 之前设置GDAL/OGR 系读取时按这个编码解析 DBF os.environ[SHAPE_ENCODING] GBK county gpd.read_file(js_admin_2022/js_county.shp) road gpd.read_file(js_admin_2022/js_road.shp) print(len(county), county.geom_type.unique()) # 要素数、几何类型 print(county.columns.tolist()) # 字段名重点看有没有被截断 print(county.head(3).drop(columnsgeometry)) # 只打属性避免刷屏 print(county.crs) # None 表示 .prj 没生效 print(road[ROAD_TYPE].value_counts(dropnaFalse)) # 等级字段的值分布SHAPE_ENCODING这个环境变量对 fiona 和 pyogrio 两种读取引擎都生效比在read_file里单独传encoding更稳因为有些版本的参数透传会丢。value_counts(dropnaFalse)里的dropnaFalse别省道路数据里等级为空是常态空值占多大比例直接决定你清洗时要不要给它一个未分类的兜底值。2.3 第一次加载的三条验证路径新手容易犯的错是能打开就算成功。能打开只是第一步几何和属性对得上、坐标系认得出、要素数没丢才算真的加载成功。我一般会走三条路交叉验证任意一条出问题就停下来查。第一条是 QGIS 直接拖入看右下角坐标显示的单位是度还是米看属性表里中文是否正常看图层信息里的要素数与ogrinfo输出是否一致。第二条是上文的 GeoPandas 脚本重点看crs和columns。第三条是ogrinfo的输出它是最接近原始文件的视角不受任何软件默认设置影响。检查项期望结果不通过说明要素数与 ogrinfo 的 Feature Count 一致读取引擎或索引文件有问题坐标系显示 CGCS2000 或 WGS 84 地理坐标.prj 缺失或被改过属性表中文正常显示县名、道路名编码不对查 .cpg几何类型面层为 MultiPolygon线层为 MultiLineString混合几何需要先统一几何类型这一项容易被跳过。面数据里混进单部件 Polygon 和多部件 MultiPolygon 在多数场景下不影响显示但做overlay或者入库建表时PostGIS 的几何类型约束会直接报错。发现混合就统一转成 Multi 类型面层用MultiPolygon线层用MultiLineString后面省事。3. 坐标系与面积口径投影选错后面算出来的数全是废的地理坐标下的面积单位是平方度这个量纲没有任何业务意义而且不同纬度上一度对应的实际距离差得远。省内南北跨度摆在那同一份数据里南北两端的县用平方度算出来的面积和用平方米算出来的面积比例关系完全对不上。所以只要涉及面积、长度、密度这类量算第一个动作永远是转投影坐标。3.1 地理坐标算面积的坑平方度到底错在哪经度一度的地面距离随纬度收窄纬度 30 度附近约 96 公里到 35 度附近只剩 91 公里左右纬度一度倒是稳定在 111 公里上下。这意味着一个南北长、东西窄的县用平方度算出来的面积会被系统性低估而且越往北低估得越多。更麻烦的是这种误差不是常数比例你没法用一个系数统一修正。正确的做法是转到等积投影。省内单点分析用高斯-克吕格三度带中央经线取 117°E 或 120°ECGCS2000 基准精度最好但江苏东西跨了两个三度带做全省统计时会在带边界出现接缝要素被切开。所以全省口径的统计我一般用 Albers 等积圆锥投影一套参数覆盖全省牺牲一点点精度换一致性跨省对比也方便。import geopandas as gpd county gpd.read_file(js_admin_2022/js_county.shp, encodingGBK) # 地理坐标下直接算面积单位是平方度只能用来做相对比较 county[area_deg] county.geometry.area print(county[area_deg].head()) # 0.0x 量级别拿去做密度 # Albers 等积投影单位米ellpsGRS80 与 CGCS2000 椭球一致 ALBERS (projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 ellpsGRS80 unitsm no_defs) county_m county.to_crs(ALBERS) county_m[area_km2] county_m.geometry.area / 1e6 print(county_m[area_km2].describe())参数里lat_1和lat_2是两条标准纬线取 25 和 47 是中国区域常用的组合两条标准纬线之间的面积变形最小lon_0是中央经线取 105 是为了全国范围对称。如果只做省内、对面积精度要求到千分位以内把lon_0改成 119 附近、标准纬线收窄到省内纬度范围误差还能再降一档。to_crs返回的是新对象原county不会变别指望它原地修改。3.2 用 GeoPandas 做批量投影与面积核算一份数据包里通常有四五个 SHP逐个手动转会漏写个循环一次处理完。这里有两个细节值得单独说一是判断几何类型决定算面积还是算长度二是输出格式换成 GPKG 而不是继续存 SHP。from pathlib import Path import geopandas as gpd ALBERS (projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 ellpsGRS80 unitsm no_defs) for shp in Path(js_admin_2022).glob(*.shp): gdf gpd.read_file(shp, encodingGBK) if gdf.crs is None: # 没有 .prj 就不要硬转先补定义再谈转换 print(f{shp.name} 缺少坐标系定义跳过) continue gdf_m gdf.to_crs(ALBERS) geom_type gdf_m.geom_type.iloc[0] if geom_type in (Polygon, MultiPolygon): gdf_m[area_km2] gdf_m.geometry.area / 1e6 else: gdf_m[length_km] gdf_m.geometry.length / 1e3 # 输出 GPKG不受 DBF 字段名 10 字节限制也没有单文件 2GB 上限 out shp.with_name(shp.stem _albers.gpkg) gdf_m.to_file(out, layershp.stem, driverGPKG) print(shp.name, len(gdf_m), round(gdf_m.geometry.area.sum() / 1e6, 1))geom_type.iloc[0]只看第一个要素前提是前面已经确认几何类型统一如果没确认改成gdf_m.geom_type.unique()打印出来人工判断更保险。输出用layershp.stem是为了 GPKG 里每层有个能对上的名字一个 GPKG 可以装多个图层管理起来比一堆散文件清爽。最后那行打印总面积是给面数据做一次总量校验——全省县级面面积加总应该和省级面的面积对得上差太多说明有碎面或者重叠这个校验在后面第 5 章还会用到。3.3 叠在线底图上偏移几百米加密坐标的锅数据本身是 CGCS2000 或 WGS 84 的经纬度叠加到国内主流在线底图上会看到整体偏移几百米方向还挺规律。这不是数据错了是底图用了 GCJ-02 加密坐标两套坐标之间有一个非线性偏移。判断方法很简单找一条明显的线性地物比如一条高速和影像上的路面量一下偏移方向和距离如果全省范围内偏移量大致相当、方向一致基本可以确认是坐标系基准不一致。处理思路分场景。做分析、做统计、入库一律用 CGCS2000别碰加密坐标因为那套坐标是为显示服务的逆向转换有精度损失转来转去最后自己都说不清数据是什么基准。做制图出图换成不加密的地形底图或者影像底图问题自然消失。如果非要叠加密底图就在出图前临时转一份显示用的副本原数据保持不动这个习惯能省掉很多到底哪份是对的扯皮。注意坐标系转换分两种操作定义投影只是给数据贴个标签不改坐标数值投影转换才会重算坐标。前者用错了数据会跑到几万公里外而且很难一眼看出来。4. 道路网与公路网落地分级、打断、密度统计一条链面数据相对好处理线数据的坑集中在几何本身。道路网和公路网两套线放在一起先想清楚是合并还是分开用做路网密度统计通常分开算两套口径做可达性分析往往要合并成一张连通图。这一章按清洗—修复—统计的顺序走一遍。4.1 等级字段清洗与分级符号化等级字段的原值分布是清洗的唯一依据别凭经验猜。跑一遍value_counts之后常见的情况是全角半角混用、带前后空格、高速公路和高速并存、部分记录为空。清洗不是把值改漂亮是把口径统一到能分组统计的程度。import pandas as pd GRADE_MAP { 高速: 高速, 高速公路: 高速, 高速路: 高速, 国道: 国道, G: 国道, 普通国道: 国道, 省道: 省道, S: 省道, 普通省道: 省道, 县道: 县道, 乡道: 乡道, 城市道路: 城市道路, } def normalize(v): if pd.isna(v): return 未分类 key str(v).strip().replace( , ).replace( , ) return GRADE_MAP.get(key, key) road[GRADE] road[ROAD_TYPE].map(normalize) print(road[GRADE].value_counts(dropnaFalse)) print(road.groupby(GRADE)[length_km].sum().sort_values(ascendingFalse))replace( , )那个全角空格是特意加的中文数据里全角空格肉眼和半角几乎没区别但strip()不会处理它分组时会莫名其妙多出一个类别。GRADE_MAP.get(key, key)的第二个参数是兜底映射表里没有的原值直接保留这样清洗后还能看到哪些值没覆盖到比一律归到其他要好排查。分级符号化建议用字典驱动而不是在 QGIS 里手点。字段值清洗干净之后QGIS 图层属性里选分类渲染、字段选GRADE再按等级调线宽高速 1.2mm、国道 0.9mm、省道 0.7mm、县乡道 0.4mm。这个层级关系是有讲究的主干路在视觉上必须压得住支路否则出图一片糊。4.2 线要素打断与拓扑检查道路数据从测绘成果导出来线要素往往在交叉口处并没有真正打断两条路交叉但几何上各走各的。显示上完全看不出来一做路径分析就露馅算法认为这两条路不连通。同样的问题还有悬挂点——线的端点悬在半空没有和其他线相接形成断头路。打断操作在 QGIS 里有现成算法处理工具箱里搜在交叉处分割线Split lines at intersections输入线图层跑完输出一份在交点处全部打断的新图层。代价是要素数会涨不少属性会按分割段复制做长度统计没问题做属性唯一性校验要重新处理。PostGIS 里对应的是ST_Node但它只保留几何不保留属性属性要事后按位置回挂麻烦一些-- 在交点处节点化输出的几何可直接用于连通性分析 CREATE TABLE road_noded AS SELECT (ST_Dump(ST_Node(ST_Collect(geom)))).geom AS geom FROM js_road; -- 建空间索引后续按位置回挂属性全靠它 CREATE INDEX idx_road_noded_geom ON road_noded USING GIST (geom);ST_Collect把所有线合成一个集合ST_Node在集合内所有交点处插入节点ST_Dump再把结果拆回单条线。这套组合适合做连通图不适合直接出图因为属性丢了。如果只是想做一次长度统计没必要走这一步直接在当前几何上算长度就行打断是用来解决连通性问题的不是万能的预处理步骤。重复线是另一个高频问题多份数据拼接或者双向道路画成两条都会产生。检查方法很直接-- 几何完全相同的重复线常见于多版本数据叠加 SELECT COUNT(*) AS dup_groups FROM ( SELECT geom FROM js_road GROUP BY geom HAVING COUNT(*) 1 ) t;这个查询在大表上会慢因为GROUP BY geom走的是几何相等比较。可以先按道路名称和起终点坐标粗筛再对候选集做精确比对。数量不多的话直接用 QGIS 的按几何删除重复项更省事。4.3 按区县汇总路网密度空间连接与面积口径对齐密度统计的完整链路是道路和县界都转到同一个等积投影道路在县界上求交切开按县代码分组累加长度再除以县面积。这里最容易错的是口径——到底算哪些等级的道路。只算高速加国道省道得到的数叫干线密度把县乡道全算进去叫路网总密度。两个数差好几倍报告里不写清楚口径读者会以为你算错了。import geopandas as gpd ALBERS (projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 ellpsGRS80 unitsm no_defs) county gpd.read_file(js_admin_2022/js_county.shp, encodingGBK).to_crs(ALBERS) road gpd.read_file(js_admin_2022/js_road.shp, encodingGBK).to_crs(ALBERS) # 干线口径高速 国道 省道 main road[road[GRADE].isin([高速, 国道, 省道])].copy() # 求交把每条道路切到县界内同时带上县代码 seg gpd.overlay(main[[GRADE, geometry]], county[[CODE, geometry]], howintersection) seg[len_km] seg.geometry.length / 1e3 stat seg.groupby(CODE)[len_km].sum().rename(main_km) area (county.set_index(CODE).geometry.area / 1e6).rename(area_km2) result county.set_index(CODE).join(stat).join(area) result[density_km_per_km2] result[main_km] / result[area_km2] print(result[density_km_per_km2].describe())overlay的howintersection只保留两个图层重叠的部分道路落在县界外的部分自动丢弃这正是密度统计要的效果。传入的子集只取GRADE和geometry是为了避免两个表里同名字段比如都有NAME在结果里变成NAME_1、NAME_2后面还得改回来。join用索引对齐前提是CODE在两个表里都唯一县级行政代码满足这个条件。一个容易忽略的点切在县界上的道路长度会按县界精确分配但道路本身的宽度信息是被忽略的也就是按中心线长度算。这对密度统计没影响但如果要做道路面积占地分析得另外用缓冲区。还有overlay对几万条线乘上百个面会跑几分钟别以为卡死了中途强退只会得到半份数据。5. 避坑清单中文乱码、.prj 丢失与面自相交的五个现场这一章是我这几年拆同类数据包时反复遇到的坑按现象—原因—解决记下来遇到时可以直接对号入座。每条都不是理论推演是实际处理过的现场。5.1 属性表中文变问号或方块现象QGIS 打开属性表县名一列显示成????或者½ËÕò这类字符但几何正常地图上也能画出来。原因DBF 是按字节存的本身不自带编码信息全靠 .cpg 文件声明。压缩包里 .cpg 丢了或者内容是空的时候GDAL 会按 UTF-8 去解析 GBK 编码的字节一个汉字两个字节被拆开解释自然就乱了。解决先看目录里有没有 .cpg。有的话打开确认内容是不是GBK或936没有就手工建一个内容写GBK重新加载。如果不想改文件GeoPandas 读取时设SHAPE_ENCODINGGBKQGIS 在图层属性 → 源 → 数据源编码里选 GBK 后重新加载。已经读进内存的图层改编码没用必须重新读一次。5.2 .prj 丢失或坐标系标错现象gdf.crs返回None或者显示的是一个明显不对的坐标系转投影之后数据跑到非洲附近去了。原因拷贝时只复制了 .shp 和 .dbf或者原始数据被人用定义投影改过标签数值没动标签和实际不符。后者更隐蔽因为文件看起来是完整的。解决从同批数据的其他图层里找一个 .prj 复制过来前提是确认整批数据同一基准。确认不了就做实测校验取两个已知地物的坐标跟公开的地理位置对一下量级经纬度应该在 116 到 122、30 到 36 这个范围内落在别处就是坐标系标错了。修正用定义投影重新贴标签不要用投影转换两者方向完全相反。5.3 县级面自相交与狭长碎面现象面积算出负数或者 0overlay抛TopologyException报某个坐标点 self-intersection。原因数字化时节点交叉、边界重叠或者从更大范围裁剪时切出了极窄的碎片。这类问题在省级边界和海岸线附近特别多。解决用 Shapely 2.x 的make_valid修复几何比老办法buffer(0)稳后者在某些自相交形态下会丢掉整个部件。from shapely.validation import make_valid bad ~county.geometry.is_valid print(f无效几何 {bad.sum()} 个) county.loc[bad, geometry] county.loc[bad, geometry].apply(make_valid) print(county.geometry.is_valid.all())is_valid是逐要素检查要素多的时候会慢几秒。修完之后一定要再跑一次面积汇总和省级面的面积对一下差超过千分之一就说明修复过程改变了拓扑得回头看是哪个要素出的问题不能直接拿去用。5.4 道路重复线与悬挂点现象路网总长度比预期多出近一倍或者路径分析跑出来一堆断头路。原因重复线来自多份数据叠加同一条路存了两遍长度自然翻倍悬挂点来自数据拼接时端点没对齐差几米也算不连通肉眼在小比例尺下完全看不出来。解决先按几何去重再按道路名称 起终点坐标取整去重两轮下来能清掉大部分。悬挂点用 QGIS 拓扑检查器规则选线不能有悬挂点容差设 1 到 5 米跑完会列出所有端点人工判断哪些是真的断头比如省界处本来就该断、哪些是数据问题。这一步别偷懒用自动吸附容差设大了会把平行的两条路吸到一起产生的错误比修掉的还多。5.5 字段名被截断与 DBF 体积上限现象COUNTY_NAME读进来变成COUNTY_NAM或者两个字段名撞车数据量再大一点写入时报字段超长或者文件超过 2GB。原因DBF 格式的字段名上限 10 字节、单文件 2GB、单字段字符数上限 254这些都是格式本身的历史包袱跟你用什么软件无关。解决分析链路里尽早换 GPKG 或 GeoPackage字段名长度、文件体积、字符编码都不受这些限制。确实要给下游交 SHP就在导出前做一次字段名映射表把长名字压到 10 字节以内并保证唯一同时把超长文本字段截断或者拆表。这个映射表建议存成 CSV 放在数据包旁边下次别人接手能直接看懂字段对应关系。6. 把数据包变成可交付成果批量出图与入库的一个固定套路数据清洗完之后落地方式基本就两条入库供查询分析出图供汇报展示。我的固定套路是先把清洗结果统一存成 GPKG再按需分发到 PostGIS 和制图脚本中间不做二次转换避免每一步都引入新的编码或精度问题。入库这一步用 ogr2ogr 最省事一条命令搞定编码、字段名、几何类型和空间索引ogr2ogr -f PostgreSQL PG:hostlocalhost port5432 dbnamegisdb usergis \ js_admin_2022/js_road.shp \ -nln js_road \ -lco GEOMETRY_NAMEgeom \ -lco FIDgid \ -lco SPATIAL_INDEXGIST \ -nlt MULTILINESTRING \ -t_srs EPSG:4490 \ --config SHAPE_ENCODING GBK-lco GEOMETRY_NAMEgeom把几何列名固定成geom后面写 SQL 不用每次查列名-lco FIDgid是给主键起个明确的名字默认的ogc_fid在跨表关联时容易和别的表撞-t_srs EPSG:4490是入库时统一到 CGCS2000 地理坐标量算再在 SQL 里用ST_Transform转到投影这样库里存的是标准基准不锁死某一种投影。SPATIAL_INDEXGIST别忘了几万条线没有空间索引一个范围查询能跑到你怀疑人生。出图这一端我用 matplotlib 直接读 GPKG比在 QGIS 里手调样式更容易复现import matplotlib.pyplot as plt import geopandas as gpd county gpd.read_file(js_county_albers.gpkg) road gpd.read_file(js_road_albers.gpkg) fig, ax plt.subplots(figsize(7, 9)) county.plot(axax, facecolor#f7f7f7, edgecolor#999999, linewidth0.5) road[road[GRADE] 高速].plot(axax, color#8e2b2b, linewidth0.9) road[road[GRADE].isin([国道, 省道])].plot(axax, color#c0703a, linewidth0.4) ax.set_axis_off() fig.savefig(js_main_road.png, dpi300, bbox_inchestight)分级画的时候注意顺序主干路后画才能压住支路bbox_inchestight去掉白边出图直接能贴进文档。中文标注要在开头设plt.rcParams[font.sans-serif]字体名按本机有的填不然标题会变成一串方框——这个坑我踩过不止一次每次换机器都要重新确认。从那以后我每次拿到新的 SHP第一件事都是跑一遍ogrinfo -so看 SRS 和字段再决定后面怎么走绝不先拖进软件。看着多花两分钟实际省掉的是后面几小时的返工。希望帮到你。本文还有配套的精品资源点击获取
觉得有用,分享给同行:

为您的企业打造数字门面

稳重轻奢商务风格,端正雅致视觉,长效耐看不易过时。

立即咨询 →