简介全国地质灾害点空间矢量shp分布数据涵盖崩塌、塌陷、泥石流、地面沉降、地裂缝、滑坡、不稳定斜坡七大类隐患点面向GIS数据处理、灾害防治规划与科研分析人员可直接用于区域灾害监测预警、风险评估、防治区划及灾后重建等实际工作。压缩包共含十八个文件以shp空间矢量文件为主体配套dbf属性表、prj坐标投影、shx索引与xml元数据整体约40.49MB文件体系完整便于在主流GIS平台中直接调用。目前已有267人浏览学习。借助这些数据可快速还原各地质灾害点的空间位置、属性信息与分布态势为地灾成因分析、历史灾情重建、未来趋势预测以及城乡规划和工程建设选址提供基础数据支持。对于地灾教学培训而言也能作为直观的案例素材帮助提升相关人员对灾害类型的识别与应对能力实用价值显著。1. 全国地质灾害点空间矢量shp分布数据这套底数到底能干什么做地质灾害易发性评价时最耗时间的往往不是算法而是先凑一套能用的底数。全国地质灾害点空间矢量shp分布数据把崩塌、塌陷、泥石流、地面沉降、地裂缝、滑坡、不稳定斜坡这7类灾害点位整理成一份shp打开就能叠到影像、行政区划和地形图上直接看。这份数据的价值在“空间矢量”四个字上分类清楚、有点位坐标、带属性字段能做缓冲区、密度、叠加统计这类正经空间分析而不是只躺在Excel表格里。适合做地灾评估、国土空间规划双评价、应急管理隐患核查和GIS数据开发的人使用。下面我不讲地质灾害的机理只讲怎么把它用好、少踩坑。2. shp数据内容拆解7类地质灾害在矢量里到底存了什么2.1 崩塌、塌陷、泥石流……7类灾害怎么区分点位落在哪儿拿到shp先别急着打开先搞明白每一类点在空间上代表什么。这7类灾害的形态差异很大点的空间含义完全不同。崩塌点一般标在陡峭山体的崩落源区也就是“从哪里掉下来”的位置而不是掉落后的堆积区。野外调查时GPS定位通常打在危岩体下方所以一个崩塌点的周边几百米都可能受威胁。塌陷点标的是地表陷落坑的中心常见于岩溶发育区或地下采空区点位和实际塌陷范围往往存在偏移需要结合坑的直径判断影响半径。泥石流的点大多标在沟口就是流通区和堆积区交界的位置这一点很关键——泥石流的影响范围在纵向上可以延伸出几公里单看点位会严重低估威胁范围。地面沉降点不一样它是面的缓变过程shp里的点更多是监测井或代表观测点不能当作灾害发生点去理解。地裂缝点则可能是某条裂缝线的端点或中心位置严格说应该用线要素表达但点数据会在后续分析中造成范围误判。滑坡点标在滑坡后缘或主滑方向轴线附近常用于圈定滑坡边界也需要结合DEM去看。不稳定斜坡的点是“还没滑但趋势明显”的潜在隐患点这类在评估里权重往往更高。一句话每个点的“空间含义”决定了后面做分析时的缓冲半径怎么设。如果一份数据的属性字段里带“灾害类型”或“类型编号”先用这个字段把7类分开检查一遍分布比直接全量做分析稳妥得多。2.2 一个完整的shp文件包不止.shp还有.shx、.dbf、.prj很多新手只把“.shp”一个文件拷走结果发给同事后打开报错。shapefile虽然叫“shp”但它是一个复合格式至少需要三个文件配套才能正常显示文件后缀作用.shp要素几何主体存放每个点的坐标和形状.shx几何索引文件没有它很多软件会提示缺失索引.dbf属性表存放灾害类型、编号、地名等字段.prj坐标系描述没有它就失去空间基准.cpg属性表字符编码声明经常被忽略.sbn / .sbx空间索引部分软件打开大数据量时用到我一般会在拷贝时把.shp、.shx、.dbf、.prj、.cpg五个文件压缩成一个zip发出去缺.prj和.cpg后面一定会来问你“为什么坐标不对”“为什么中文乱码”。更重要的是做任何空间分析前先打开.prj文件看一眼坐标系统哪怕用记事本看也可以这一眼能省掉后面大半的坐标系翻车。2.3 用Python打开地质灾害点shp坐标系统、字段与点数一次查清在ArcGIS/QGIS里打开当然可以但批量和自动化场景下用Python检查要素几何类型、字段构成和坐标系更快。常见做法是用pyshp读基础信息配合GDAL的ogr读坐标系import shapefile # 读取地灾点shpcpg声明为UTF-8时指定encoding with shapefile.Reader(disaster_points.shp, encodingutf-8) as sf: print(几何类型编号, sf.shapeType) # 1点 3线 5面能确定这份数据到底是点还是线 print(要素数量, len(sf)) print(字段列表, [f[0] for f in sf.fields[1:]]) # 打印前3条属性记录看灾害类型字段的实际值 for i, rec in enumerate(sf.records()): if i 3: break print(rec)这段代码解决三个问题shapeType确认是不是点数据len(sf)确认数据量级字段列表确认属性表结构。sf.fields[1:]是为了跳过pyshp默认的DeletionFlag字段。再看坐标系用ogr能直接导出prj的文本描述from osgeo import ogr ds ogr.Open(disaster_points.shp, 0) layer ds.GetLayer(0) sr layer.GetSpatialRef() if sr: print(坐标系名称, sr.GetName()) print(Proj4定义, sr.ExportToProj4()) else: print(未定义坐标系prj文件缺失)如果输出是“未定义坐标系”后面所有叠加分析都可能出错这是这份数据能不能直接用的分水岭。多数公开的地质灾害shp使用CGCS2000经纬度坐标系也就是EPSG:4490少数是WGS84EPSG:4326两者在低纬度地区平面误差小于一米但混用时还是会偏。3. 拿到数据先做三步坐标系统一、编码修复与字段规整3.1 坐标系的“玄学”先确认prj再谈投影转换地质灾害shp最容易踩的就是坐标系。常见来源里自然资源系统提供的数据多是CGCS2000早期历史数据可能是西安80。CGCS2000和WGS84在平面坐标下差异不大但西安80与CGCS2000的偏移是几十米到上百米量级跟在线底图叠加时一眼就能看出“点不在路上”。所以第一步永远是确认prj。prj丢失时只能用已知的控制点反推或者直接找数据来源方要坐标信息。最稳妥的做法是拿到一份带完整prj的shp作为基准把其它来源的点位向它对齐。但多数情况下你的地质灾害点shp已经带prj真正需要做的是投影转换。如果原数据是CGCS2000经纬度EPSG:4490需要把它转到适合测距的投影坐标系才能做缓冲区等操作。用一个完整脚本示例import os from osgeo import ogr, osr src_path disaster_points_4490.shp dst_path disaster_points_4527.shp # 源坐标系CGCS2000 经纬度 src_sr osr.SpatialReference() src_sr.ImportFromEPSG(4490) # 目标坐标系CGCS2000 / 3-degree Gauss-Kruger zone 39117°E # 具体带号要根据数据主经度换算这里以117°E为例 dst_sr osr.SpatialReference() dst_sr.ImportFromEPSG(4527) transform osr.CoordinateTransformation(src_sr, dst_sr) # 打开源文件 src_ds ogr.Open(src_path, 0) src_layer src_ds.GetLayer(0) src_defn src_layer.GetLayerDefn() # 输出文件存在则先删除避免追加导致重复要素 driver ogr.GetDriverByName(ESRI Shapefile) if os.path.exists(dst_path): driver.DeleteDataSource(dst_path) dst_ds driver.CreateDataSource(dst_path) dst_layer dst_ds.CreateLayer(disaster_points_4527, dst_sr, ogr.wkbPoint) # 复制原属性字段到新图层 for i in range(src_defn.GetFieldCount()): fd src_defn.GetFieldDefn(i) dst_layer.CreateField(fd) # 逐要素读取几何并做坐标转换再写入新图层 for src_feat in src_layer: geom src_feat.GetGeometryRef() geom.Transform(transform) dst_feat ogr.Feature(dst_layer.GetLayerDefn()) dst_feat.SetGeometry(geom) for i in range(src_defn.GetFieldCount()): dst_feat.SetField(i, src_feat.GetField(i)) dst_layer.CreateFeature(dst_feat) dst_feat None src_ds dst_ds None print(转换完成)这个脚本的核心逻辑是通过osr.SpatialReference定义源与目标两套坐标系用CoordinateTransformation做转换对象然后逐要素读取几何执行geom.Transform(transform)。属性字段必须手动逐个复制因为CreateLayer只建几何和默认字段不会自动搬属性。EPSG:4527只是CGCS2000三度带第39带117度中央经线如果数据覆盖范围很大要按经度选择对应的带号或者改用Albers等积投影用于全国尺度分析。其实日常工作中最常用的是命令行的ogr2ogr比上面的Python脚本更容易写也更不容易漏步骤ogr2ogr -s_srs EPSG:4490 -t_srs EPSG:4527 disaster_points_4527.shp disaster_points_4490.shp注意-s_srs指定源坐标系、-t_srs指定目标坐标系。如果原文件的prj本身正确连-s_srs都可以省略但如果prj缺失必须手动指定否则转换结果就是错的。这一步是“坐标系是玄学”这句玩笑话的核心——不是坐标系复杂是无数人跳过了prj检查。3.2 中文乱码GBK与UTF-8的dbf编码修复打开属性表发现全是一串问号或乱码是地灾shp非常常见的症状。原因是dbf属性表有内部编码.cpg文件声明了它但很多历史数据用的是GBK生成时没有写cpg或者cpg声明成ISO-8859-1。在QGIS里可以直接改图层编码右键图层 → 属性 → 数据源 → 编码逐个试GBK和UTF-8。但批量处理时用ogr2ogr重编码更高效。先把原文件编码改成UTF-8避免在其它软件里再乱码ogr2ogr -lco ENCODINGUTF-8 disaster_points_utf8.shp disaster_points_gbk.shp-lco是图层创建选项ENCODINGUTF-8让新写入的dbf声明为UTF-8字段里的中文内容也会被正确转换。这个操作不会丢几何属于可逆修复可以放心使用。还有一个小技巧转换完用上一章的pyshp脚本重新读一次记录确认中文字段能正常print出来再继续分析。乱码数据就算强做分析后面出图也是一个一个乱码标签返工成本很高。3.3 给7类灾害建type_code把零散字段归到一个标准不同来源的shp灾害类型字段名称五花八门有的叫“灾害类型”有的叫“类型名称”“LEIBIE”值是中文缩写还有带括号备注的比如“崩塌岩质”。要统计就必须先统一编码。常见做法是新建一个整数字段type_code把字符串映射成固定数字from osgeo import ogr # 以可写方式打开直接修改原shp建议先用一份副本操作 ds ogr.Open(disaster_points_clean.shp, 1) layer ds.GetLayer(0) if layer.FindFieldIndex(type_code, 0) 0: field ogr.FieldDefn(type_code, ogr.OFTInteger) field.SetWidth(2) layer.CreateField(field) # 灾害类型字符串 - 编码映射 type_map { 崩塌: 1, 塌陷: 2, 泥石流: 3, 地面沉降: 4, 地裂缝: 5, 滑坡: 6, 不稳定斜坡: 7 } for feat in layer: raw feat.GetField(灾害类型) if raw is None: feat.SetField(type_code, 0) else: code type_map.get(raw.strip(), 0) feat.SetField(type_code, code) layer.SetFeature(feat) ds None print(type_code字段生成完毕)这里用findFieldIndex先检查字段是否已存在避免重复运行时报错。type_map.get(raw.strip(), 0)是安全写法字段值里带空格、带括号备注时都归到0方便后续排查有没有没映射上的类型。我一般会再打印一下code为0的记录数量如果占比超过1%说明原数据的类型写法有额外变化需要补映射关系。这一步做完后面不管是按类型过滤、图例配色还是统计出表都会轻松很多。没有统一编码之前直接做符号化图例上会出现十几个颜色根本没法看。4. 地灾点shp的空间分析缓冲区、点密度与县域叠加统计4.1 缓冲区分析评估道路与居民区周边的灾害威胁范围地灾点是点位但实际威胁是面状的。做风险评估时最常用的操作是缓冲区分析比如“距灾害点500米内的道路有哪些”“居民区是否在滑坡威胁范围内”。缓冲区最核心的坑是单位。如果shp还是经纬度坐标系缓冲区工具会把距离当作“度”来处理500米被当成500度整个图层被放大到无法直视。所以做缓冲区之前必须先把数据投影到以“米”为单位的坐标系EPSG:4527这类高斯投影或UTM都可以。用ogr2ogr配合SQLite方言可以一条命令完成ogr2ogr buffer_result.shp disaster_points_4527.shp -dialect sqlite -sql SELECT ST_Buffer(geometry, 500) AS geometry, type_code, 灾害类型 FROM disaster_points_4527ST_Buffer(geometry, 500)里geometry是投影坐标系下的几何对象500的单位跟随数据坐标系也就是米。-dialect sqlite启用GDAL内嵌的SQLite空间函数可以直接在SQL里做空间计算。字段部分只挑选后面分析需要的列能减小输出文件体积。缓冲区半径怎么设不是技术问题是业务问题。一般参考值崩塌和滑坡取300—500米泥石流取1000—2000米地面沉降作为面状缓变过程不适合用点缓冲区描述。如果不知道该怎么设先对7类灾害分别做几个不同半径的缓冲区叠到遥感影像上看对比哪个半径能覆盖实际灾点周围的地形变化特征这一步是值得多试几次的。4.2 点密度分析用热力图识别灾害高发区密度分析解决的是“哪些区域灾害点最密集”的问题服务易发性分区和隐患排查。QGIS里用热力图Heatmap工具ArcGIS里叫核密度分析Kernel Density原理差不多以每个点为圆心按一个搜索半径生成加权密度曲面。热力图的两个参数直接影响结果形态搜索半径Radius和像元大小Pixel size。半径设太小图上到处是孤立的“火山口”看不出区域趋势太大则一片糊边界信息丢失。我一般按数据的平均点间距来试比如点位平均相距3公里就从10公里半径起步分别生成5公里、10公里、20公里三版对比。输出栅格的像元大小用500米或1公里即可像元过小只是文件变大对结论影响有限。还有人会做格网统计先把研究区用渔网分割shp切成规则格网再统计每个格网内的灾害点数。这个思路适合做县级或网格化的易发性底图但要注意渔网必须与灾害点层用同一套投影坐标系且研究区边缘的半格网很容易漏统。性价比最高的还是热力图调参快、出图直观。4.3 与县域行政区划边界shp叠加按县统计各类灾害数量业务报告里最常出现的表格就是“xx县有多少处滑坡、多少处泥石流”。这需要把灾害点shp与县域行政区划边界shp做空间连接。用ogr2ogr的SQLite方言可以一条命令完成ogr2ogr county_stats.shp county_boundary.shp -dialect sqlite -sql SELECT b.name AS county_name, count(p.type_code) AS cnt FROM county_boundary AS b LEFT JOIN disaster_points AS p ON ST_Intersects(p.geometry, b.geometry) GROUP BY b.name关键点在LEFT JOIN。它保证没有灾害点的县也会出现在结果表里count值为0而不是直接消失。如果用INNER JOIN那些“零灾点”的县会被过滤掉统计表就不完整领导问起来很难解释。ST_Intersects是空间谓词表示两个几何相交即命中count(p.type_code)按县分组计数如果只想统计单一灾害类型在SQL后面加WHERE p.type_code 6即可。执行完这条命令后原来的面图层会带上灾害数量属性可以直接按“数量等级”做 choropleth 配色。但这里也有一个隐藏bug落在县界上的点会被相邻两个县同时统计导致总数比原始点数多这个问题我放在下一章专门讲。5. 地灾点shp使用避坑5个高频翻车现场与处理办法5.1 翻车一只拷了一个.shp图层打不开现象同事发来一个.shp文件ArcGIS图层加载时提示“无法打开”。原因shapefile是复合文件主文件记录几何、dbf记录属性、shx做索引缺一不可。只传.shp等于只传了一半数据。解决拷贝时连同.shx、.dbf、.prj、.cpg一起打包确认收到的是zip压缩包再解压。我自己的习惯是交接数据时列一个文件清单逐个打钩后再发避免来回折腾。5.2 翻车二点和遥感影像错位几百米坐标系对不上现象灾害点叠加到在线影像底图上点位偏离实际位置几百米甚至更远。原因大概率是prj文件缺失或声明坐标系与实际不符。很多历史数据用西安80或北京54却被误标成WGS84。解决先看.prj里写的坐标系名称没有prj时通过影像上几个已知地物点反推。真正定位到坐标系后再按第3章的代码重投影统一到目标坐标系。几百米的偏移做宏观统计看不出问题一旦做现场核查导航到了错误位置才是真正的麻烦。5.3 翻车三属性表全是问号中文编码翻车现象QGIS打开属性表灾害类型字段显示为“”或乱码。原因dbf内部是GBK编码但打开时被当成UTF-8解析或者.cpg文件缺失导致软件猜错编码。解决先试QGIS图层属性里的“数据源编码”切到GBK或GB2312如果还乱就说明数据本身编码可能损坏更稳妥的是用ogr2ogr -lco ENCODINGUTF-8重写一份。重写前留好原始文件这条是后悔药别一上来就原地覆盖。5.4 翻车四数据没有时间注记拿旧底数做现状评价现象用一套2010年的地灾点数据写当前现状评估报告被人质疑“这些点大部分已经治理销号了”。原因地灾点数据有极强的时效性隐患点经过降险、治理或自然稳定后会动态销号历史shp只代表调查时点的状态。解决拿到数据先看元数据或XML里有没有调查年份字段没有就向数据提供方确认写报告时注明数据时间基准并对已经在治理名录中销号的点位做剔除或降权。这一条不是技术问题却是所有分析里最致命的问题。5.5 翻车五点在县界上被重复统计数量虚高现象按县统计的各类型灾害数量相加后比原始总点数多了几十个。原因点要素落在县界线上ST_Intersects同时命中了相邻两个县。解决按县统计前先对点做微量负缓冲比如ST_Buffer(geometry, -0.000001)个位数的百万分之一度相当于把点内缩几厘米足以让贴近边界的点只在一边命中。也可以先做空间连接再按点的唯一ID去重保留行政区划面积占比最大的那个归属。用前者最简单一条SQL就能解决也是我常用的做法。6. 把地灾点shp用起来shp转wkt/txt对接业务、转3dtiles上三维6.1 shp转wkt/txt给后端系统提供文本坐标业务系统不一定能直接解析shp后端开发经常要求提供文本格式坐标。把灾害点shp导出成WKT或者纯文本是Shp转TXT这类需求的常见落地方式尤其适合导入PostGIS或MySQL空间字段。用shapely可以非常简洁地完成import shapefile from shapely.geometry import shape # 读入shp逐要素转WKT sf shapefile.Reader(disaster_points_4490.shp, encodingutf-8) with open(disaster_points.wkt, w, encodingutf-8) as f: f.write(type_code,wkt\n) for s, rec in zip(sf.shapes(), sf.records()): geom shape(s.__geo_interface__) # 把记录的第二个字段视为灾害类型编码按实际结构调整 code rec[1] if len(rec) 1 else f.write(f{code},{geom.wkt}\n)__geo_interface__是pyshp转换给shapely的标准接口一行就能把点变成shapely几何对象再调用geom.wkt得到WKT文本。写出的文本可以直接用于PostGIS的ST_GeomFromText(POINT(...))。如果后端只需要x、y坐标改成f.write(f{geom.x},{geom.y}\n)即可。注意导出前把坐标系统一成WGS84经纬度业务系统拿到4326坐标后按高德或天地图的要求再做投影偏移。6.2 shp转3dtiles把灾害点搬到三维地图上在Cesium、UE这类三维场景里做地灾态势研判需要把shp转成3dtiles。过去我用的路径是先把shp转成GeoJSON再借助桌面端切片工具生成3dtilesogr2ogr -f GeoJSON disaster_points.geojson disaster_points_4490.shpGeoJSON是Web三维引擎的基础交换格式绝大多数切片工具都能直接消费。转换时注意属性字段要保留type_code设置点样式时才能按灾害类型分别染色。三维场景对坐标系的敏感度比二维低一些CGCS2000经纬度和WGS84经纬度在Cesium上叠加误差肉眼基本看不出来但严谨起见还是在转换前统一成WGS84免得日后做量测时被问到“这个距离准不准”。这里有一个经验顺便给到准备上三维的同行切片工具的坐标系和原点设置一定要和原始shp一致别在转换链路上混用不同坐标系。我有一回就是shp是CGCS2000中途用WGS84生成过一次GeoJSON切出来的3dtiles在Cesium里偏移了大概十几米。排查了半天最后还是回到shp重新转。地面沉降这类点位做三维标注有一定效果但完整的沉降趋势还是靠面状插值更可靠。说回这套全国地质灾害点shp本身。它适合做宏观尺度上的底数核查、易发性初筛和规划参考但不适合直接替代野外的隐患点台账。不管后续是转wkt、转3dtiles还是做密度分析我都会先确认两件事坐标系有没有prj数据是什么年份的。这两个确认动作是我处理这类数据不翻车的习惯少一个都会在某个节点回来补课。希望帮到你。本文还有配套的精品资源点击获取