ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

全国地质灾害点shp数据全流程解析:坐标系转换与空间分析实战

全国地质灾害点shp数据全流程解析:坐标系转换与空间分析实战 简介涵盖崩塌、塌陷、泥石流、地面沉降、地裂缝、滑坡和不稳定斜坡等七类灾害的全国地质灾害点空间矢量shp数据面向GIS分析师、地学研究者和城乡规划人员可直接用于灾害空间分布制图与风险早期识别为灾害防治与应急管理提供基础底图。压缩包共18个文件包含shp要素文件、dbf属性表、prj坐标系、shx几何索引及xml元数据其中dbf记录灾害类型、位置等属性prj确保地理坐标精度xml提供元数据说明整体40.49MB结构清晰适配ArcGIS、QGIS等主流平台使用。目前已有259人学习/下载过该数据包。借助这套矢量数据可快速掌握各类隐患点的位置、属性与空间关系服务于地质灾害易发性评价、监测点布设、防灾规划编制亦可用于高校地学相关课程的教学演示与科研分析为人身财产安全保障提供基础数据支撑。1. 全国地质灾害点 shp 数据一份能直接落进 GIS 的灾害分布底图做地质灾害评估、国土空间规划或者环评前期踏勘的人最头疼的往往不是算法而是底图。你手上可能握着几万条灾害点台账但它们是 Excel 表格、是分散的 Word 文档甚至只是纸面上的坐标记录——要把这些点落到地图上跟地形、流域、行政区叠加做分析第一步就得把它们变成空间矢量数据。这份全国地质灾害点 shp 分布数据包含崩塁、塌陷、泥石流、地面沉降、地裂缝、滑坡、不稳定斜坡七类灾害点直接把灾害分布做成了标准 Shapefile 空间图层省掉从零开始建库、配坐标、录属性的过程。适合做区域灾害易发性评价、危险性分区、边坡与矿山勘查前期筛选以及 GIS 教学里最缺的“真实灾害数据”素材。2. 数据内部结构七类灾害点怎么分类、怎么挂属性2.1 灾害类型字段从符号化到统计的关键依据这份数据的核心价值在于“分类”。拿到手后先打开属性表重点看灾害类型字段——它记录每一个灾害点的具体类型也就是崩塁、塌陷、泥石流、地面沉降、地裂缝、滑坡、不稳定斜坡这七类。这个字段是你后续做筛选、做符号化、做统计的起点。打开属性表的常见方式是在 ArcGIS Pro 或 ArcMap 里右键图层打开属性表或者用 QGIS 直接查看。我一般会先按灾害类型字段做一次分组统计确认每一类有多少个点、字段值是否有空值或者错别字。这个动作花不了两分钟但能避免后面做专题图时出现符号错配或者漏统计。import geopandas as gpd gdf gpd.read_file(地质灾害点.shp, encodingutf-8) print(gdf.head()) # 查看字段名找到灾害类型字段 print(gdf.columns.tolist()) # 按灾害类型统计数量 type_field 类型 # 实际字段名以属性表为准 count gdf[type_field].value_counts() print(count)统计完成后把数量跟数据包自带的说明文档做一次比对能初步判断数据是否完整。如果没有说明文档就按七类标准灾害类型去核对发现异常值再排查原因。2.2 坐标信息文件里藏着“哪个坐标系”的答案很多人拿到 shp 文件后第一步就叠加底图发现点全部跑到海里或者飞到天上大概率就是坐标系不匹配。打开图层的属性查看坐标系定义这一步不能省。常见的坐标系有三种WGS84 地理坐标系、CGCS2000 地理坐标系、以及高斯投影平面坐标系。这份数据如果是全国范围的点分布大概率是地理坐标系如果是分省或分县数据可能是投影坐标系。你需要做的是把图层的坐标系定义记下来后面叠加底图或者做面积、距离计算时再决定是否投影。# 查看坐标系 print(gdf.crs) # 查看经纬度范围 print(gdf.total_bounds)这里有个判断技巧total_bounds输出的范围如果接近[73, 3, 135, 53]就是全国范围且是地理坐标系。如果数值特别大比如[20000000, 4000000, 50000000, 6000000]多半是投影坐标系。知道这个区别后再做任何转换都心里有底。2.3 图形要素结构点要素的 Geometry 怎么读这份数据是点要素不是面也不是线。点要素的 Geometry 比较简单就是每个灾害点的 XY 坐标。但要注意的是点要素的精度取决于原始采集方式有的点是野外 GPS 实测有的点是从遥感影像上判读后人工标注还有的点可能是纸质图件数字化出来的。精度不同后续做缓冲区和叠置分析的可靠性就不同。# 检查几何类型 print(gdf.geom_type.unique()) # 检查是否有空几何 print(gdf.is_valid.sum())一般不会出现空几何但如果出现可以用gdf.dropna(subset[geometry])把空值行删掉避免后续分析报错。2.4 字段到底有哪些一份实用的列字段清单数据包里具体有哪些字段以实物为准但一份标准的灾害点 shp 通常会包含以下五类字段编号字段用于唯一标识每个灾害点灾害类型字段记录崩塁、滑坡等分类地理位置字段记录省、市、县、乡镇、村等行政信息规模或险情等级字段按小、中、大、特大为标准划分坐标字段记录经度、纬度或者投影坐标。拿到数据后最值得做的第一步就是把关键字段整理成一张表确认哪些字段能直接用、哪些字段需要清洗。地质灾害点的属性清洗是一个脏活比如同一个乡镇的写法可能是“XX乡”也可能是“XX乡镇”用文字做关联时容易出问题。我一般会先把字段裁剪出来只保留分析会用到的列导出成一个精简版 shp后续操作会快很多。# 保留关键字段 cols [编号, 类型, 省, 市, 县, 规模等级, 经度, 纬度, geometry] gdf_clean gdf[cols].copy() gdf_clean gdf_clean.drop_duplicates() gdf_clean.to_file(地灾点_精简.shp, encodingutf-8, driverESRI Shapefile)到这一步数据已经能用了。接下来最常做的事就是把数据转成下游需要的格式。3. shp 转 kml、json、txt从 ArcGIS 到在线地图的转换实录3.1 在线地图分享shp 转 kml 的两个选择把灾害点数据放到 Google Earth 或者手机上查看最常见的方式是转成 kml/kmz。ArcGIS Pro 自带的图层转 KML 工具就可以但需要注意单位问题。转 KML 时如果图层是地理坐标系直接转换即可如果是投影坐标系需要先做投影转换否则生成的 kml 位置可能会偏。# ArcGIS Pro 中操作路径 # 工具箱 - 转换工具 - KML - 图层转KML # 输入图层: 地灾点_精简 # 输出文件: 地灾点.kmz如果是 QGIS 用户也可以用右键图层另存为格式选 Keyhole Markup Language [KML]再设置坐标系为 WGS84。这里的关键参数是Output file name和CRSCRS 务必选择 WGS84。用 Python 转换时最简单的方式是用geopandas配合simplekml或者直接用fiona读取后手动构造 KML 节点。后者写起来代码较长建议直接用 GDAL 的命令行ogr2ogr -f KML 地灾点.kml 地灾点.shp -t_srs EPSG:4326这行命令会把 shp 转成 kml-t_srs EPSG:4326是强制转换为 WGS84 经纬度防止目标坐标系不一致导致位置偏移。KML 文件用文本编辑器就能打开是纯 XML 结构如果出现乱码或者节点信息丢失可以检查 XML 的编码声明和字段映射。3.2 给 WebGIS 用shp 转 GeoJSON 的参数细节目前很多在线可视化平台和自研 WebGIS 都直接吃 GeoJSON转换时要注意两点坐标系要转成 WGS84 经纬度编码必须是 UTF-8。以前踩过坑用 ArcGIS 的“要素转 JSON”直接转出来可能是 GCS_WGS_1984 但字段编码不是 UTF-8导致网页上显示中文问号。常见的做法是直接用 QGIS 另存为 GeoJSONEncoding 选择 UTF-8CRS 选择 WGS84。如果使用geopandas一行代码就能搞定gdf_wgs84 gdf.to_crs(epsg4326) gdf_wgs84.to_file(地灾点.geojson, driverGeoJSON, encodingutf-8)to_crs(epsg4326)是精度损失最小的做法比先转投影再转回来靠谱。GeoJSON 的字段名不要用中文部分下游工具解析中文键名会报错建议转换前先把中文列名改成英文。3.3 回表格shp 转 txt/csv 的两种路径做统计分析和归档管理时经常需要把属性表导成 txt 或 csv。最简单的办法是右键图层——打开属性表——表选项——导出这个操作在 ArcGIS 里人人会做。但如果是批量导出且每个 shp 的字段不一致脚本化的方式更有优势。import csv with open(地灾点.csv, w, newline, encodingutf-8-sig) as f: writer csv.writer(f) writer.writerow(gdf_clean.columns.tolist()) for idx, row in gdf_clean.iterrows(): values [str(row[col]) for col in gdf_clean.columns if col ! geometry] writer.writerow(values)utf-8-sig的编码方式是为了兼容 Excel 直接打开 CSV 不乱码这一点比 utf-8 更实用。如果表格里经纬度字段已经有现成的坐标值直接选对应列导出即可如果只有 geometry则需要用gdf.geometry.x和gdf.geometry.y把坐标取出来再写入。3.4 dwg 转 shp、json 转 shp反向转换的对称坑做工程的人经常手上有 CAD 的地形图或者灾害点测图需要从 dwg 转成 shp。跟 shp 转 kml 相反dwg 转 shp 最大的坑是比例尺和坐标系。CAD 图形往往没有定义坐标系只有图上坐标转 shp 之前必须找到 CAD 图上的控制点坐标做一次空间校正否则结果就是一张不能跟真实地物对齐的空图。常见的做法是先把 dwg 中需要的图层单独导出为 dxf再在 ArcGIS Pro 中用“CAD 至地理数据库”工具导入。这里要注意Import CAD工具里的空间参考设置必须手动给定坐标系。如果是点状灾害标记CAD 里的TEXT和INSERT图元会分别变成注记和点要素需要额外关联处理。json 转 shp 的情况类似GeoJSON 一般自带坐标系声明直接用geopandas.read_file读取再另存为 shp 即可。唯一需要注意的是 GeoJSON 里如果存在MultiPoint或GeometryCollection写入 shp 时可能报错因为 Shapefile 对几何类型的支持有限必须先explode()成单点再导出。# GeoJSON 转 shp json_gdf gpd.read_file(地灾点.geojson) json_gdf json_gdf.explode(index_partsFalse) # 多部件要素拆分为单部件 json_gdf.to_file(地灾点_from_json.shp, encodingutf-8, driverESRI Shapefile)explode()是处理多部件要素的关键不做这一步shp 写入时容易遇到几何类型不匹配的报错。4. 坐标系与属性问题排查避坑记录与修复方法4.1 现象叠加后点整体偏移 50 到 200 米方向固定原因源数据坐标系是 CGCS2000 或西安 80叠加的底图是 WGS84或者反过来。两者在局部区域的偏差可以到几十米甚至上百米肉眼看不出来但叠加道路或地块边界时很明显。解决明确源数据的原始坐标系。右键图层属性查看坐标系定义如果不是自己想要的坐标系用“投影”工具做转换而不是直接在图层属性里改。改图层属性里的坐标系只是改“声明”不会改变实际坐标数值属于自欺欺人。4.2 现象属性表中文变成问号或乱码导入 geopandas 后字段名一堆奇怪字符原因shp 文件的 .dbf 属性表编码跟读取工具预期的不一致。老数据常用 GBK 或 GB2312 编码新工具默认按 UTF-8 读取导致乱码。解决读取时显式指认编码。如果是 QGIS在添加矢量数据时会询问编码如果是 geopandas指定encodinggbk或encodingutf-8。试错型的处理方式是用chardet先检测一下 dbf 文件的编码再重新读取。import chardet with open(地灾点.dbf, rb) as f: raw f.read() result chardet.detect(raw) print(result)检测结果如果显示gb2312或gbk重新读取时指认对应编码。编码错的成本很高属性表全部乱套白做半天分析。4.3 现象输出 shp 后属性表字段名被截断为 10 个字符原因Shapefile 格式的 dbf 字段名最多支持 10 个字符超长部分会被截断。这是格式的硬限制不是操作问题。解决输出 shp 前把所有字段名手动改成 10 字符以内的英文字母例如disaster_type改为d_type。如果字段很多写个小脚本统一改名更高效gdf.columns [col[:10] for col in gdf.columns] gdf.to_file(地灾点_截断.shp, encodingutf-8)col[:10]是简单的截断处理实际项目里建议直接重新命名为有意义的短字段名比截断更可控。4.4 现象同一个灾害点在两张图中重复出现属性略有不同原因数据来自不同普查阶段或不同来源合并时没有做去重。全国尺度的灾害点数据往往由多个省份、多次调查拼接而来点重叠是高频问题。解决先按类型、坐标、行政字段做联合去重。我可以接受两个点相距小于 5 米且类型相同视为重复点这比精确定位去重更合理# 生成空间坐标字符串 gdf[geom_key] gdf.geometry.apply(lambda g: f{round(g.x, 5)}_{round(g.y, 5)}) gdf[type_key] gdf[type_field].astype(str) gdf[dup_key] gdf[geom_key] _ gdf[type_key] gdf gdf.drop_duplicates(subsetdup_key, keepfirst) gdf gdf.drop(columns[geom_key, type_key, dup_key])保留五位小数相当于约 1 米的去重精度这个参数要根据数据精度调整。野外 GPS 精度一般就是 3 到 5 米所以五位到六位小数是不错的选择。4.5 现象用 select by attributes 筛选崩塁时结果包含滑坡原因灾害类型字段存在同词异写比如“崩塁”在部分属性里写成“崩蹋”或“崩塌”。不同普查团队的术语不统一这种坑在合并数据里非常常见。解决筛选前先做字段值归一化。查看字段的所有唯一值建立映射表统一替换mapping { 崩塁: 崩塌, 崩蹋: 崩塌, 塌陷: 塌陷, 不稳定斜坡: 不稳定边坡 } gdf[类型] gdf[类型].replace(mapping)做地质灾害分析时类型字段直接影响统计口径和后续分区评价权重这一步不能跳过。替换完后再做一次value_counts()检查各类别数量是否合理。5. 进阶实战灾害点密度分析与风险分区叠加拿到这份全国灾害点数据后真正有实操价值的用法是把灾害点跟行政区、流域、地形数据叠加做密度分析。地质灾害易发性评价的基础工作就是计算单位面积内灾害点的数量、频率和规模再叠加坡度、岩性、降雨等因子做加权。以下步骤可以复现一份最简单的灾害点密度分区图。先把数据投影到适合做面积计算的坐标系。全国尺度用 Albers 等积投影省级尺度用高斯克吕格投影不要直接在经纬度下做面积统计结果会失真。# 投影到 Albers 等积投影 albers ESRI:102025 # 适合中国区域的 Albers 投影 gdf_proj gdf_wgs84.to_crs(albers) print(gdf_proj.crs)投影后把灾害点跟行政区面要素做空间连接统计每个行政区的灾害点数量再除以行政区面积得到点密度。admin gpd.read_file(市级行政区.shp, encodingutf-8) admin_proj admin.to_crs(albers) # 空间连接把灾害点属性落到行政面 joined gpd.sjoin(gdf_proj, admin_proj, howinner, predicateintersects) # 分行政区统计灾害点数量 counts joined.groupby(市代码).size().reset_index(name灾害点数) admin_density admin_proj.merge(counts, left_on市代码, right_on市代码, howleft) admin_density[灾害点数] admin_density[灾害点数].fillna(0) admin_density[点密度] admin_density[灾害点数] / (admin_density.geometry.area / 1e6)geometry.area是投影坐标系下的面积单位为平方米除以 1e6 转成平方公里。得到的点密度字段按自然间断点分级设色就能产出一张直观的灾害点分布密度图。这里提醒一点行政面数据和灾害点数据要统一代码字段比如都用区划代码前六位避免名称匹配时因为“区”和“县”写法的差异对不上。如果想进一步做缓冲区分析可以按灾害点规模生成不同半径的缓冲区再做合并和叠加得到灾害影响范围面。# 按规模等级差异设置缓冲区半径 buf_dict {巨型: 3000, 大型: 1500, 中型: 800, 小型: 300} gdf_proj[buf_r] gdf_proj[规模等级].map(buf_dict).fillna(300) buf gdf_proj.buffer(gdf_proj[buf_r])缓冲区半径要结合数据本身的精度和用途设定地质灾害影响范围并不完全等于圆缓冲区但作为前期初筛手段效率很高。做生态环境评价或者道路选线回避灾害易发区的时候这个缓冲区图层能直接拿去跟路线方案做相交统计。最后建议做一次完整性校验。比对灾害点图层和原始统计表中的总点数、各类型点数是否一致同时用 Python 跑一遍几何有效性检查确认没有自相交或空几何。数据基本可靠之后再开始做论文、报告里的图件。从那以后我每次拿到灾害点数据都会强制走一遍坐标系检查、字段值归一化、去重三步流程再做任何叠加分析。这套习惯救了我很多次希望帮到你。本文还有配套的精品资源点击获取
返回列表