ARTICLE DETAIL

资讯详情

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

河南省30米DEM数据处理全流程:从拼接裁剪到地形分析

河南省30米DEM数据处理全流程:从拼接裁剪到地形分析 简介河南省三十米分辨率数字高程模型DEM数据包面向地理信息系统、遥感与地理信息分析人员基于ASTER GDEM第三版拼接生成覆盖河南省全域。数据采用GeoTiff格式和WGS84地理坐标系每个像元代表实地三十米×三十米范围可用于地形起伏分析、坡度坡向提取、流域特征识别及土地利用规划等场景也常用于水文分析和工程选址。压缩包共十个文件大小约一百七十六兆字节除核心高程栅格文件外还包含河南省行政边界矢量文件以及地理配准、投影坐标系和属性表等辅助文件可无缝导入常用GIS平台使用。已有五百九十人学习下载说明这类高程数据在区域地形研究中有一定需求。对需要快速获取河南省级DEM底图的用户而言这份数据省去了自行下载多幅ASTER GDEM数据再拼接的步骤配合配套边界矢量可直接用于裁剪、叠加、出图等操作同时基于WGS84参考系可与其他地理数据灵活集成适合作为区域地理研究的统一底图。1. 河南省DEM30米分辨率这张栅格图里面装的不是颜色而是高度拿到一版河南省DEM30米分辨率的GeoTIFF第一反应是“这不就是一张地形图吗”真正在GIS里放大到某个县、某个山沟时才会意识到每个像素记录的是该位置的绝对高程而不是遥感影像里的光谱反射值。30米分辨率意味着每个栅格单元对应地面上大约30米乘30米的面积河南省约16.7万平方公里的范围在这个尺度下大约会展开成1.8亿个有效像元。做水文分析、洪水淹没模拟、道路选线、光伏选址、或者把DEM和DOM、DSM叠在一起做地物分类时这套数据解决的核心问题是如何在省级尺度上拿到一套空间位置一致、高程基准统一、能被ArcGIS/QGIS和GDAL直接读取的连续高程场。适合用它的人包括自然资源行业的实施工程师、做国土空间规划的GIS分析师、以及需要把地形因子作为输入特征做机器学习建模的研究人员。需要先说清楚的是河南省DEM并不存在唯一官方版本常见的数据源是SRTM、ALOS AW3D30、ASTER GDEM等公开全球数据以及测绘部门发布的成果数据不同来源在局部区域的精度和空值分布有差异。这篇内容就按“数据判断、拼接、裁剪、地形分析、质量验收”这条线把河南省DEM从原始瓦片到可用成果的完整流程讲清楚。2. 河南省DEM数据源选型与格式基础先搞清楚30米从哪里来2.1 不同来源的河南省DEM精度差异比想象中大30米分辨率并不是某一个数据源的专有名称。目前业界常用的公开高程数据里SRTM 1弧秒约30米覆盖全球北纬60°到南纬56°河南省全境都在覆盖范围内数据特点是在中低纬度地区水平精度相对稳定但局部山谷地区容易出现噪声ALOS AW3D30同样是1弧秒网格基于PRISM立体像对提取在垂直精度上通常比SRTM更平滑尤其在地形起伏大的区域表现更好ASTER GDEM标称也是30米但经常存在残余的云遮挡伪影和条纹噪声用之前必须先看空值分布。测绘部门发布的省级DEM成果通常来源于更高分辨率的立体航摄影像或LiDAR点云经过滤波、编辑、接边处理后重采样到30米网格质量最好但获取门槛和许可限制也最高。做河南省全域分析时我一般建议优先选ALOS AW3D30或测绘成果而不是默认下载SRTM。数据源标称分辨率河南省覆盖垂直精度参考优点典型问题SRTM 1弧秒30米全覆盖常见说法约6至10米LEM90数据完整度高接边成熟软件兼容好局部山谷噪声水域可能异常ALOS AW3D30约30米全覆盖常见说法RMSE约5米垂直一致性较好平滑度好下载分幅多处理量大ASTER GDEM V330米全覆盖变化较大覆盖完整条纹、云伪影、空洞测绘省级成果30米重采样全覆盖按国家规范分级精度可靠有质量报告数据许可受限2.2 DEM文件格式与坐标系为什么我不建议直接存成IMG先说一个常见误区有人把河南省DEM从原始来源转成IMG或IMG压缩格式理由是“省空间”。实际工作中ESRI GRID、IMG这类格式在处理大范围连续栅格时随机读取性能明显不如GeoTIFF加LZW/Deflate压缩而且GDAL之外的很多Python库如rasterio、xarray对GeoTIFF的支持最完整。完整的dem文件内部必须包含有效的Nodata值定义、地理变换参数六参数和坐标系描述缺任意一个后续裁剪和坡度分析都会出问题。gdalinfo henan_aw3d30_mosaic.tif输出里需要重点看这几行Size is 15440, 12340Origin (110.350000000000000,36.180000000000000)Pixel Size (0.000277777777778,-0.000277777777778)Coordinate Reference SystemNoData Value-9999。如果NoData Value没有显示后续计算坡度时会把空值当作0高程参与运算生成的结果就是一张布满黑色沟壑的坡度图。处理这类问题我会用一条命令统一修正Nodata和数据类型gdal_translate -a_nodata -9999 -ot Float32 -co COMPRESSDEFLATE -co TILEDYES \ henan_raw.tif henan_dem_clean.tif-a_nodata -9999把没有高程值的像元统一标记为-9999-ot Float32保证高程不因整数截断损失精度-co TILEDYES把内部存储切成块ArcGIS和QGIS读取大文件时速度会快很多。另一个值得注意的点是坐标系如果原始数据是WGS84经纬度像素尺寸显示为0.0002777778度约30米这个状态下直接算坡度需要指定比例尺否则软件默认按度计算结果完全不可用。常规做法是先投影到省级或区域统一坐标系后再做地形分析省界矢量边界用河南省级行政区划的Shapefile或GeoPackage文件。3. 用GDAL把河南省DEM从原始瓦片拼成分幅成果再按省界精确裁剪3.1 先拼VRT再做物理裁剪省掉大量中间文件公开数据源下载的河南省DEM通常不是一个整幅tif而是按经度和纬度切割成若干1乘1度的分幅文件。把这些dem文件合起来最直接的办法是写一个循环逐幅镶嵌但这样做会需要巨大的磁盘临时空间。正确做法是用GDAL的VRT虚拟栅格机制先用一个XML文件把所有分幅文件逻辑组织起来不复制实际像素然后再按河南省范围裁剪。gdalbuildvrt henan_raw.vrt \ N34E110.hgt N34E111.hgt N34E112.hgt \ N35E110.hgt N35E111.hgt N35E112.hgt \ N36E110.hgt N36E111.hgt N36E112.hgt这里假设每幅文件的命名采用类似SRTM的经纬度分幅规则实际文件名以你下载结果为准也可以用通配符写成gdalbuildvrt henan_raw.vrt N3*E11*.hgt。gdalbuildvrt不会检查所有瓦片是否严格对齐但会记录每幅图的范围和分辨率所以后续读取VRT时会自动做一个“最近邻”的像素匹配。这个阶段不要急着做重投影保持经纬度网格等裁剪时一并处理。3.2 用矢量边界裁剪河南省DEM并把坐标系转到投影坐标做完VRT之后下一步就是把河南省DEM裁剪成业务可用的成果文件。裁剪dem数据的常规做法是用行政区划边界作为-cutline这比用矩形范围-te更符合应用需求。边界文件用GeoPackage或Shapefile都行关键是边界图层本身和DEM必须在空间上有交集且坐标系统可以不同GDAL会自动处理投影转换。gdalwarp -cutline henan_province.gpkg -crop_to_cutline \ -dstnodata -9999 -r cubic -of GTiff \ -co COMPRESSDEFLATE -co TILEDYES \ henan_raw.vrt henan_dem_cut.tif这条命令里-cutline指定省界矢量-crop_to_cutline表示把输出栅格的范围严格对齐到省界外接矩形边界外的像元全部设为Nodata-r cubic用三次卷积重采样在地形连续变化区域比双线性更平滑比最近邻更不容易产生锯齿。这里有两个参数需要特别留意第一-dstnodata -9999最好和源数据Nodata保持一致后续做坡度或填洼时才不会引入“假高程”第二如果需要在投影坐标系下输出需要加上-t_srs EPSG:xxxx河南跨多个投影带小区域用UTM 50N或者CGCS2000 3度带全省范围则常用Albers等积投影具体选择取决于你的分析目标是算面积还是量距离。3.3 如果只需要河南省某个市或某个县的DEM就要先分割再分别处理全省范围的数据不是所有项目都用得上有些场景只要洛阳、南阳或者某个流域。把大文件拆成若干小文件可以避免在ArcGIS里每次加载都卡顿。分割dem数据有两种思路一种是在GIS里用“按掩膜提取”工具逐县提取操作简单但重复另一种是用gdalwarp批量循环提取脚本化后更高效。按县界分割的脚本逻辑是这样的for f in *.gpkg; do name$(basename $f .gpkg) gdalwarp -cutline $f -crop_to_cutline -dstnodata -9999 \ -r near -of GTiff henan_dem_cut.tif dem_${name}.tif done-r near在裁剪县级数据时已经够用因为只涉及裁切不涉及分辨率变化重采样方式对结果影响很小。县级分割出来的dem文件通常尺寸在几千乘几千像素后续无论是计算坡度、提取等高线还是做可视域分析响应速度都会有明显提升。这里要补充一个常见坑如果县界矢量属性表里的字段名包含中文或特殊字符basename或循环变量解析可能失败建议先把待处理矢量统一重命名为拼音或英文名再跑脚本。4. 河南省DEM的地形分析实操坡度、坡向、山体阴影与等高线生成4.1 投影之后再算坡度别拿经纬度数据直接开算30米分辨率的河南省DEM做完裁剪之后就能直接出坡度了吗实际操作中如果栅格还是经纬度坐标系gdaldem虽然可以运行但默认按度为单位计算梯度结果数值会严重偏小或异常。常见做法是先投影到Albers等积投影或UTM使X和Y方向单位统一为米再开始地形分析。投影这一步放在裁剪前或裁剪后都可以我习惯放在裁剪后因为裁完的数据量更小重投影速度更快。gdalwarp -t_srs EPSG:3857 -r cubic -tr 30 30 -dstnodata -9999 \ henan_dem_cut.tif henan_dem_proj.tifEPSG:3857是Web墨卡托如果分析场景偏WebGIS可视化用这个可以但如果是专业水文分析建议换成EPSG:4528或EPSG:4547这类CGCS2000 3度分带投影长度变形更小。-tr 30 30强制重采样到30米网格避免投影转换过程中像素尺寸变成29.999米或30.001米这种微小不一致在后续拼接多个分幅时会被放大。4.2 用gdaldem批量生成坡度、坡向与山体阴影参数按省级地形特点来调高程数据处理最常用的一组命令是gdaldem下的坡度、坡向和山体阴影。河南省地貌类型比较杂西部有伏牛山、太行山余脉东部是平原因此省级坡度图如果按照全国统一标准分级平原地区会大面积落在0到2度无法体现微地形。算法上可以用-p选项计算百分比坡度也可以保留默认的度数值关键是要把输出类型设为Float32坡度图精度比原始DEM降一档会丢失细节。gdaldem slope henan_dem_proj.tif henan_slope_deg.tif \ -s 1.0 -p -co COMPRESSDEFLATE gdaldem aspect henan_dem_proj.tif henan_aspect.tif \ -zero_for_flat -co COMPRESSDEFLATE gdaldem hillshade henan_dem_proj.tif henan_hs.tif \ -z 2.0 -az 315 -alt 45 -co COMPRESSDEFLATE-s 1.0在这里是垂直比例因子针对已投影的30米数据一般设为1如果是经纬度数据想近似换算需要用-s 111120这类值-p改变的是输出单位不改变算法-az 315是光源方位角对河南省多山地区而言西北方向打光能更好显示山脊线走向-alt 45是太阳高度角抬高会让阴影变短降低会拉长阴影。hillshade输出的是一张灰度图本身不是高程派生因子但拿它做DOM或遥感影像的叠加坡度图对判断地形与植被分布的关系很有用。做土地利用分析时我通常把坡度小于2度的区域单独提取出来作为平原候选区再用gdal_calc.py做条件判断。gdal_calc.py -A henan_slope_deg.tif --outfileplain_2deg.tif \ --calcA2 --NoDataValue0 --typeByte--calcA2逐个像元判断坡度是否小于等于2度结果为1表示满足0表示不满足。输出用Byte类型存储文件只有几十MB加载和统计速度比浮点栅格快很多。如果你拿到的成果里已经包含DSM地表模型而你需要的是DEM地形模型操作上可以用点云滤波或形态学开运算去除地表建筑物和植被冠层把DSM生成DEM的过程放到LiDAR数据处理软件里完成GDAL命令行无法直接完成这一步。4.3 等高线生成在QGIS里做矢量化的常用路线DEM转等高线是很多规划项目的刚需。QGIS里做这个操作不复杂但参数设置决定生成的等高线是“能看”还是“能用”。在工具箱里搜索“等高线”对应算法是“GDAL 等高线”输入河南省DEM设置间隔为100米可以快速生成全省骨架等高线区域分析时通常加密到10米或5米输出为GeoPackage格式比Shapefile更能避免字段名长度限制。等高距的选择没有绝对标准常见做法是参考地形图比例尺1比5万图用20米或50米1比1万图用5米或10米。生成之后必须检查等高线与原始DEM的套合关系把等高线叠加到山体阴影图上如果出现等高线穿越河谷但高程值明显偏高或偏低的情况多半是原始DEM边缘的空值没有被正确处理。5. 质量验收的关键一步拿一个Python脚本把河南省DEM的无效值、范围和坐标系检查完数据做完不等于数据可用。交付或入库之前我会跑一套固定检查流程避免项目后期才发现某个县的DEM是空白或范围没对齐。用GDAL的Python绑定直接读取tif元数据和像元块比反复打开QGIS工程更高效。脚本不长核心逻辑是并行读取统计信息和有效数据范围最后按省界边界做一次空间套合判断。from osgeo import gdal import numpy as np def check_dem(path, expected_crsEPSG:4528): ds gdal.Open(path) if ds is None: return {error: 文件无法打开} band ds.GetRasterBand(1) nodata band.GetNoDataValue() arr band.ReadAsArray() if nodata is not None: invalid (arr nodata) | ~np.isfinite(arr) else: invalid ~np.isfinite(arr) total arr.size invalid_ratio float(invalid.sum()) / total gt ds.GetGeoTransform() xmin gt[0] ymax gt[3] xmax xmin gt[1] * ds.RasterXSize ymin ymax gt[5] * ds.RasterYSize crs ds.GetProjection() esri_crs_ok 4528 in crs or CGCS2000 in crs print(f无效值占比: {invalid_ratio:.4%}) print(f范围: {xmin:.3f}, {ymin:.3f}, {xmax:.3f}, {ymax:.3f})这里~np.isfinite(arr)用于捕获NaN和无穷值避免只检查Nodata导致漏检gt[5]通常是负值所以ymin ymax gt[5] * RasterYSize得到的是南边界坐标。广东省或河南省这类跨带省域 DEM如果打印出东西跨度超过500公里说明投影选取或范围检查有误需要重新确认是否按省界裁剪。把无效值占比控制在0.1%以内是比较理想的水平超过这个数值就说明原始数据存在漏洞或拼接时没有设置好Nodata。最后的建议是在处理河南省DEM这类大范围数据时每一步都保留一个中间VRT文件它本身不占空间但能让你随时回到某个处理节点重新调整参数而不必把裁剪、投影、填洼整体重跑一遍。gdaldem hillshade这一步生成的山体阴影文件也可以作为质检底图肉眼检查时重点看黄河滩区是否有明显的条带伪影以及伏牛山和太行山交界处是否存在接边高差。如果发现接边处出现“台阶状”错位可以先在原始瓦片层面重新拼接再回到裁剪流程。所有检查项通过之后这个河南省DEM30米分辨率文件就能进入后续的坡度分级、汇水区提取或三维地形可视化环节了。本文还有配套的精品资源点击获取
返回列表