ARTICLE DETAIL

资讯详情

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

全球海底滑坡分布shp点位数据处理与空间分析全流程

全球海底滑坡分布shp点位数据处理与空间分析全流程 简介这份全球海底滑坡分布点文件以 shp 矢量格式收录了覆盖全球海域的 500 余处海底滑坡点位面向海洋地质、地质灾害评估、海底工程等方向的科研人员与 GIS 学习者可用于快速获取滑坡空间分布、辅助区域规律分析或教学示例。压缩包共 8 个文件除核心的 shp 点图层外还包含 dbf 属性表、prj 投影坐标系、shx 几何索引、sbn/sbx 空间索引及 cpg 编码等标准配套文件整体仅 9KB数据量轻简适合直接加载至 ArcGIS、QGIS 等平台使用。已有 57 人学习下载便于快速了解全球海底滑坡宏观格局。通过该 shp 文件读者可获得已配准 WGS1984 坐标的全球滑坡点位及其属性字段省去自行搜集、矢量化与坐标转换流程从而将更多精力投入滑坡分布成因、密度制图与风险分析等后续工作。1. 全球海底滑坡分布点shp500多个点位怎么变成论文底图做海底工程选址或者海洋灾害评估时最缺的不是专业方法而是能直接把分散文献里的滑坡事件汇总成一张点图的现成底图。这份压缩包解压后就是一个 shp 点文件全球 500 多个海底滑坡点位属性表里带经纬度坐标能直接拖进 ArcGIS 或 QGIS 显示也方便和海底地形、水深栅格做叠加分析。它适合刚接触海洋 GIS、需要快速出一张全球分布图的从业者也适合想核对自己编目有没有漏区域的老手。但要提醒一句它不是开箱即用的结论图坐标系、属性完整度、点位可信度都要自己先过一遍下面按实际操作顺序来拆。2. 打开前的三件事坐标系、属性表与投影变换2.1 加载shp第一步先查空间参考别直接画点shp 不是单文件而是由 .shp、.dbf、.shx、.prj 等一组文件组成的。拖动到 ArcGIS Pro 或 QGIS 时软件主要靠 .prj 判定坐标体系如果压缩包里缺 .prj图层会以“未知坐标”加载点位虽然能显示但后续一切投影、量算、叠加都可能错位。所以拿到 rar 解压后第一件事不是双击预览而是右键图层打开“属性 → 源”把空间参考看清楚。import arcpy arcpy.env.workspace rD:\landslide_data fc Global_Landslide_Points.shp desc arcpy.Describe(fc) print(坐标系:, desc.spatialReference.name) print(坐标单位:, desc.spatialReference.linearUnitName) print(范围:, desc.extent)这段代码是检查数据的第一道防线。spatialReference.name 会返回“GCS_WGS_1984”或“WGS_1984_Web_Mercator”这类名称linearUnitName 用来判断坐标单位是度还是米。范围同样关键如果坐标值落在经度 -180~180、纬度 -90~90 之内但坐标系标成了 Web Mercator说明数据在制作环节被错误定义。常见做法是把这种情况先记下来读完属性表再决定是补定义坐标系还是做真正的投影。如果发现空间参考缺失用 Define Projection 补上 WGS84。这个操作只是给数据贴一个正确标签不改变坐标值对后续分析最重要的一步就是这里import arcpy fc rD:\landslide_data\Global_Landslide_Points.shp sr arcpy.SpatialReference(4326) # WGS84 经纬度 arcpy.DefineProjection_management(fc, sr)注意经纬度数据用 4326不要顺手选成 3857。凡是从网上下载的旧点位文件缺 .prj 的概率很高这一步不做后面投影转换会全部白跑。2.2 属性表字段先判断这份shp能做什么分析右键打开属性表字段决定分析上限。这类全球滑坡点位 shp 常见的字段一般有编号ID、经纬度Lon/Lat、滑坡名称、水深、规模等级、来源文献等。但如果字段很少只有 FID 和一个坐标字段那只能做位置分布做不了规模统计。字段类型典型字段名用途位置类POINT_X / POINT_Y坐标转出、空间连接属性类名称 / 编号图例标注、检索环境类水深 / 坡度 / 区域与 DEM 叠加验证来源类文献 / 年份确定数据时间边界用 arcpy 快速打印全部字段和类型判断哪些能用import arcpy arcpy.env.workspace rD:\landslide_data fc Global_Landslide_Points.shp desc arcpy.Describe(fc) for f in desc.fields: print(f.name, f.type)输出里的 f.type 会返回 OID、Double、String 等类型这决定后续是用数值分级还是文本标注。如果水深字段是 String 类型做数值运算前要先转成 Double这是最常见的翻车点之一。经验是如果字段里同时有经纬度和来源文献这份数据适合做全球尺度编目对比如果只有坐标我一般会先做一步坐标转出转成 txt 或 GeoJSON方便交给不熟悉 ArcGIS 的同事核对转换方法在第 3 章写。2.3 投影变换全球展示与区域分析用两套坐标系底图的坐标体系决定你在什么尺度上做空间分析。全球分布展示用 WGS84 经纬度没问题但一旦做核密度、缓冲区、面积统计经纬度坐标会带来畸变——高纬度区域被拉大计算出的面积和距离不可信。常见做法是全球尺度的密度分析投影到 Web Mercator3857区域研究比如北大西洋、地中海用所在区域的 UTM 分带。import arcpy in_fc rD:\landslide_data\Global_Landslide_Points.shp out_fc rD:\landslide_data\Landslide_3857.shp arcpy.Project_management(in_fc, out_fc, arcpy.SpatialReference(3857))Project 和 DefineProjection 的区别在于Project 实际改变坐标值Define 只是贴标签。如果原数据是经纬度用 Project 转 3857 后坐标会从度变成米但空间关系不变。之后和海底地形栅格比如 GEBCO 转出的栅格叠加时要确保栅格也被投影到同一坐标系否则错位会被放大几十公里。QGIS 用户可以直接在图层上右键“导出 → 要素另存为”在 CRS 里选 EPSG:3857效果一样。要注意 QGIS 会按当前工程的 CRS 自动重投影显示可能掩盖底层坐标不一致问题所以文件命名最好带坐标系后缀比如写成 Landslide_4326.shp、Landslide_3857.shp免得后续接手的人去猜坐标。提示文件命名带上坐标系后缀是减少团队协作误判最简单的习惯。3. 从点位到结论核密度分析、水深叠加与成果导出3.1 核密度搜索半径与栅格像元的选值逻辑拿到 500 多个点第一反应通常是“画张点图”但点图只能看分布看不出相对密集程度。核密度分析能把点位的空间密度拟合成一张连续栅格在图上用暖色标出高密度区。这项分析的前提是点要素已经投影到米制坐标系所以先做第 2 章的投影转换再做密度。import arcpy from arcpy.sa import KernelDensity arcpy.env.workspace rD:\landslide_data arcpy.env.overwriteOutput True out_raster KernelDensity( Landslide_3857.shp, population_fieldNONE, cell_size1000, # 栅格分辨率1 km search_radius300000, # 搜索半径300 km area_unit_scale_factorSQUARE_KILOMETERS ) out_raster.save(kde_global.tif)cell_size 决定输出栅格分辨率1 km 适合全球尺度想细看区域细节可以调到 500 m但再小会出现大量空栅格。search_radius 是最难调的参数直接控制平滑程度300 km 的搜索结果能看出大陆坡区域的分片你也可以用“点位平均最近距离 × 3 到 5 倍”这个经验值反推。population_field 设为 NONE 表示每个点权重相等如果属性表里有规模等级字段可以改成该字段名让大滑坡在密度图上压过小滑坡。判断搜索半径是否合理看结果是不是“一片黑”或者“一堆碎斑块”。前者半径过大、平滑过度后者像元太小、噪声太多。这个调参过程没有绝对正确我一般会先按点集范围宽度的 2% 试一次再上下调整最终参数记录在地图文档里写报告时给评审交代参数来源。3.2 与海底地形叠加从DEM提取水深验证点位合理性海底滑坡的分布不是随机的集中在陆坡、海沟、海山边缘。如果你有全球海底地形栅格常见来源是 GEBCO 或 ETOPO1 分钟分辨率足够把点位的水深、坡度提取出来能直接判断这份 shp 的质量正常大陆坡滑坡水深多在几百米到数千米如果大量点位水深落在 0 米附近或者上万米说明坐标或高程基准有问题。import arcpy from arcpy.sa import ExtractMultiValuesToPoints arcpy.env.workspace rD:\landslide_data ExtractMultiValuesToPoints( Landslide_3857.shp, [[gebco_3857.tif, bathy], [slope_3857.tif, slope]], BILINEAR)BILINEAR 表示双线性内插值用于海底地形这种连续表面如果只取像元中心值可以改成 NONE。bathy 和 slope 是新建字段名提取后打开属性表就能看到每个滑坡点对应的水深和坡度。之后做一步简单统计把水深字段排序看中位数落在哪个区间整体是否符合陆坡环境。如果没有现成的海底地形栅格也可以从 DEM 里生成等深线 shp用 ArcGIS 空间分析工具箱里的 Contour等值线工具把 DEM 转成等深线再做空间连接把等深线属性挂到点图层。效果不如栅格提取精细但作为快速校验足够。要注意 DEM 和点图层的坐标系必须一致否则提取出的数值是错的。这一步和“从 DEM 提取 shp”的思路一样核心都是保证输入数据空间参考统一。3.3 成果导出点转txt、GeoJSON与3D Tiles发布分析做完成果往往不止一张 ArcGIS 工程图。报告协作时同事不一定装了 GIS 软件最常见的是要一份 txt 或表格。用 TableToTable 把 dbf 直接导成 txt保留经纬度和关键字段import arcpy arcpy.env.workspace rD:\landslide_data arcpy.TableToTable_conversion( Global_Landslide_Points.dbf, rD:\landslide_out, landslide_points.txt)TableToTable 会读取 .dbf 字段并输出成 txt 或 CSV分隔符默认是制表符适合直接进 Excel 或 Pandas。如果想要更通用的 GeoJSON用 ogr2ogr 更快一条命令搞定ogr2ogr -f GeoJSON landslide_points.geojson Global_Landslide_Points.shp注意如果原 shp 缺 .prjogr2ogr 会按 EPSG:4326 默认读取坐标值不变但标签可能被纠正输出后一定要再用 QGIS 打开检查一下属性字段有没有丢。如果想把点图层发布到三维地球或 Web 端shp 不能直接被前端解析常见流程是先转 GeoJSON再转成 3D Tiles 切片最后在 Cesium 或 Mapbox 里加载。ArcGIS Pro 自带“创建 3D Tiles”工具能直接输出QGIS 用户可以通过插件导出 glTF 后再切片。这个环节链路长最容易翻车的是坐标系3D Tiles 对坐标要求严格喂进去的数据最好明确标注 EPSG:3857千万别把投影坐标当经纬度传。导出前还有一个习惯先精简字段。用“字段 → 消除字段”把 FID、内部编号等技术字段删掉只留坐标和结论字段导出的 txt 和 GeoJSON 体积更小交给外部的人也更安全。海洋点位数据就是核心资产不要把无关内部字段一并带出去。4. 避坑手册这份shp最常见的五个翻车点这类点位 shp 我处理过不少多数坑跟坐标系、属性编码有关跟专业水平关系不大更多是数据来源太杂。下面五条是踩过的重灾区按现象到原因到解决写清楚你可以照着排查。4.1 加载后一片空白坐标系和视图没对上现象双击 shp 添加进 ArcGIS 后内容列表里有图层但画布空白。原因一是 shp 缺 .prj软件不知道点位在哪显示成“未知坐标系”二是图层没有被缩放到数据范围视图停留在一个完全无关的区域。解决先右键图层选“缩放至图层”排除视图问题若仍空白用 arcpy.Describe 打印坐标范围按 2.1 的方法补定义坐标系。遇到坐标数值正常但范围显示异常时基本就是 .prj 缺失别去动坐标值直接 Define Projection 为 4326。4.2 属性表中文乱码编码格式不一致现象名称类字段显示成问号或乱码。原因shp 的 dbf 表历史遗留编码是 GBK而 ArcGIS Pro 默认按 UTF-8 读取旧文件QGIS 有时能读但操作系统区域设置不同也会出错。解决在 QGIS 加载矢量图层时把编码选项改为 GBK 或 CP936。ArcGIS 里没有一次性转换入口常见做法是先转成 GeoPackage 或 File Geodatabase再重新加载。转完之后如果字段值还是乱的就检查源 dbf 的原始编码不要反复转。4.3 点位批量落到陆地上数值和坐标系对不上现象点图层和海岸线叠加后大量点位落在大陆内部或者 x、y 数值出现六位数。原因这份 shp 在制作时可能用经纬度保存但被贴了 Web Mercator 标签也可能数据源从 UTM 投影转回经纬度时少做了一步数值直接叠加。解决看坐标范围。经纬度应落在 -180~180、-90~90若 x 是几十万的量级说明数据被当作投影坐标使用这时候要用 Project 反向转回 WGS84而不是 DefineProjection。注意DefineProjection 不改坐标值Project 才改坐标值用错一步整批点位就偏了。4.4 核密度结果是一张空栅格参数组合不合适现象KernelDensity 跑完输出栅格几乎全为空值或只有一个极其孤立的高值点。原因search_radius 设置得比相邻点间距还小导致大多数像元周围没有点落入或 cell_size 设得过小数据点太少形不成连续分布。解决先用“平均最近邻”工具算出点平均间距或者用 arcpy.Describe 查看数据范围把搜索半径设为平均最近距的 3~5 倍输出像元设为比搜索半径小一到两个量级的数值。跑完后检查栅格最大值如果最大值大得离谱或者非空像元数量极少基本就是参数失衡。4.5 与海底地形错位几十公里基准面不一致现象点图层和 GEBCO 栅格叠加后点位整体朝一个方向偏移在海沟和陆坡折线处尤其明显。原因两个数据源的空间参考基准面不同点数据是 WGS84栅格或海岸线数据是其他基准面叠加时没有指定 datum transformation软件用了默认近似方式。解决在工程里先统一数据框坐标系再为图层间指定变换。ArcGIS 里做 Project 时把“地理变换”参数手动选择常见的 WGS84 与其他基准之间用带名称的变换方法。海洋数据错位问题肉眼不一定能看出来叠加等深线或海岸线后再下结论。5. 拿到手先做一遍点位校验重复点、离群点与字段缺失shp 能打开、能上图只是第一步点位本身可不可信是另一回事。我习惯拿到手先跑一遍重复点检查因为编目合并时两份文献引用同一个滑坡事件很容易被录两次。import arcpy fc rD:\landslide_data\Global_Landslide_Points.shp arcpy.AddGeometryAttributes_management(fc, POINT_X_Y_Z) count_before int(arcpy.GetCount_management(fc).getOutput(0)) freq_table rD:\landslide_data\freq.dbf arcpy.Statistics_analysis(fc, freq_table, [[POINT_X, FIRST], [POINT_Y, FIRST]], POINT_X;POINT_Y) dup_count 0 with arcpy.da.SearchCursor(freq_table, [POINT_X, POINT_Y, FREQUENCY]) as cur: for row in cur: if row[2] 1: dup_count 1 print(原始点数:, count_before) print(重复坐标组数:, dup_count)AddGeometryAttributes 会生成 POINT_X、POINT_Y 字段Statistics_analysis 按坐标分组统计FREQUENCY 大于 1 的坐标组就是重复点。重复点位多说明数据源合并时清洗不够做密度分析时权重会被放大。另一种校验是离群点。海底滑坡点位如果大量出现在深海平原或完全无地形变化的区域值得怀疑。用第 3 章提取出来的水深字段做极值检查超过 1 万米海沟极限深度或者为负值陆地高程的点位单独挑出来对照文献背景。字段缺失检查也顺手做掉。对每个字段统计空值数量缺失太多的字段直接放弃使用fields [f.name for f in arcpy.ListFields(fc) if f.type not in [OID, Geometry]] for f in fields: null_count 0 with arcpy.da.SearchCursor(fc, [f]) as cur: for row in cur: if row[0] is None: null_count 1 print(f, 空值数量:, null_count)注意 OID 和 Geometry 字段没有空值概念要跳过。空值率超过 40% 的字段后续分析里通常不值得信任除非你有原始文献能补齐。这套校验被我用成固定流程了。从那以后每次拿到陌生 shp不管来源写得多规范都先跑一遍检查再动手分析省下的都是后面的返工时间。希望帮到你。本文还有配套的精品资源点击获取
返回列表