尧图网页设计 
简介河北省30米分辨率DEM数据包是一份面向GIS从业者、测绘与地理分析人员的地理空间数据产品。数据源自ASTER GDEM V3版本以GeoTiff栅格格式存储海拔信息并采用WGS84坐标系可直接用于地形坡度、坡向分析洪水模拟、规划选址及三维景观渲染等场景。压缩包共10个文件总大小约229.09MB除核心的高程栅格文件外还附带河北省行政边界shp矢量文件及配套的dbf、prj、sbn等属性与投影辅助文件便于在ArcGIS等平台中与边界叠加使用省去自行配准与裁剪的步骤。发布至今已有698人学习下载数据版本较新时效性有保障。对于需要快速获取河北全境高分辨率地形数据并进行流域分析、交通选线等工作的用户这份打包好的DEM与边界文件组合可直接投入使用减少数据预处理工作量。1. 河北省DEM 30米数据先想清楚要拿它做什么拿到这份河北省DEM 30米分辨率数据第一件事不是急着拖进ArcGIS渲染成彩色而是先想清楚要拿它做什么。河北从坝上高原、燕山、太行山一路过渡到华北平原30米像元对这种省级尺度的分析刚刚好全省坡度分级、区域汇水分析、灾害易发性初筛都能扛住但要拿来算单边坡精细坡度或做乡镇级土方设计分辨率和精度都不够。这份数据适合GIS从业者、规划师、水利与测绘背景的工程师直接当底图和分析输入新手也能用它练完整的地形分析流程。动手之前先花五分钟确认投影、值域、NoData、空洞分布四个参数它们决定后面整个工作流要不要返工。2. 选型与元数据体检这份DEM凭什么能用2.1 30米分辨率在河北能扛住的场景和扛不住的任务先说一个直观的数字30米像元意味着每个栅格代表 900 平方米的地面。河北全省面积约 18.88 万平方公里一张覆盖全省的30米栅格大约有 2.1 亿个像元单波段Int16裸数据大概 400MB 出头加上金字塔和压缩后也就几百兆。如果换成10米数据像元数直接涨到 19 亿数据量接近10倍ArcGIS里做一次填洼和流向计算的时间也差不多翻倍。对区域级项目来说这个性价比非常关键。实用场景里30米DEM适合做的是省级国土空间规划里的坡度分级水利上的汇水区划分和河网初提地质灾害易发性的区域性初筛以及林业、农业的立地条件分类。这些任务的共同特点是“看趋势、看格局”单个像元的误差不会颠覆结论。反过来它扛不住的任务也很明确单边坡开挖设计、村庄尺度的排水管线分析、矿区沉降监测、工程土方量精算这些场合需要的是厘米级到米级数据30米DEM连“参考”都很勉强。具体适合不适合可以对照下面这张表来判断任务类型是否适合原因省级坡度分级适合30米能识别出大尺度地形趋势坡度分级误差可控流域汇水区划适合汇水面积通常以km²计像元尺度影响有限地质灾害区域初筛适合关注的是地形组合关系不是单点高程单边坡精细设计不适合900m²一个像元边坡边界完全无法表达村庄级排水分析不适合排水管线精度要求远高于30米土方量精算不适合地形细节被平滑挖填方误差会放大2.2 主流30米数据源对比SRTM、ASTER、ALOS与Copernicus市面上30米DEM的全球公开数据源就那么几个特性差异不小。如果你拿到手的数据没标注生产源先按这个表做个预判断数据源标称分辨率公开精度指标已知软肋河北适用性SRTM v31弧秒30米相对垂直精度约±6m高纬度覆盖到60°N为止河北没问题但在复杂山地有空洞平原区表现平滑适合ASTER GDEM v330米官方称±7~14m山区条带噪声、云遮蔽空洞明显太行山、燕山要检查空洞ALOS AW3D30 v330米由5米重采样相对精度约5m量级云覆盖区有缺失但比ASTER干净综合质量好山地表现优于ASTERCopernicus GLO-3030米垂直精度优于SRTM来源于TanDEM-X平面位置精度高噪声均匀适合做分析和出图经验上讲在河北这种地貌组合复杂的地方SRTM的平原高程误差小ASTER在山地容易出现细长条带状空洞ALOS整体干净但偶尔有云补丁Copernicus是后起之秀均匀性好。我一般拿到资源后不纠结于它属于哪个源直接跑一遍体检脚本看空洞分布和噪声形态就能猜个八九不离十。另外要提醒一句如果数据表面在建筑区和茂密林区明显隆起那你拿到的可能是DSM数字表面模型而不是DEM。DSM带屋顶和树冠直接做水文分析会出现径流在房顶和树梢上流的怪事这是很多新手把河网提得歪七扭八的隐性原因。2.3 数据到手先查四个参数投影、值域、NoData、空洞打开数据之前先写个检查脚本比在ArcGIS里肉眼刷屏高效得多。我用Python的rasterio来做几分钟能跑完import rasterio import numpy as np path heb_dem_30m.tif with rasterio.open(path) as ds: print(投影:, ds.crs) print(范围:, ds.bounds) print(像元尺寸:, ds.res) arr ds.read(1) nodata ds.nodata valid arr[arr ! nodata] total arr.size nodata_count total - valid.size print(有效高程min/max:, valid.min(), valid.max()) print(NoData占比: {:.2f}%.format(nodata_count / total * 100))这段脚本的输出要对着四个检查点看第一投影如果是WGS84地理坐标系后面做面积计算前必须转投影第二像元尺寸必须接近30×30实际出现29.99或30.00都正常但出现1或2开头的数说明数据被重采样过第三有效高程最小值如果是0或-32768基本可以断定边界或空洞被填了0值要回炉处理第四NoData占比超过2%就得把NoData像元分布图导出来看确认空洞是不是集中在太行山、燕山的阴坡。这四个参数过关再往下做坡度水文分析才踏实。空洞位置怎么看简单做法是把NoData重分类成1、有效值设成0用ArcGIS的Spatial Analyst画出来。如果空洞呈细长条带、走向和山脊一致这是典型的ASTER天窗特征如果呈块状可能是云或原始数据缺失。记住这个判断后面第5章补洞时能少走弯路。3. 把分幅数据变成河北省一张图镶嵌、裁剪与投影选择3.1 分幅TIFF镶嵌成全省栅格如果这份河北省DEM是按标准分幅交付的第一步是把几十个分幅拼成一张全省图。命令行环境里我习惯先用GDAL把分幅“虚拟拼接”成VRT再做一次实际转换这样中途出错不用重新读几十个文件# 先生成VRT虚拟栅格不实际复制数据 gdalbuildvrt heb_dem_merged.vrt *.tif # 再转成正式TIFF统一像元类型和NoData gdal_translate -ot Int16 -a_nodata -9999 -co COMPRESSDEFLATE heb_dem_merged.vrt heb_dem_merged.tif关键参数在第二行-ot Int16把栅格统一成16位有符号整型防止某个分幅是浮点型导致整张图体积膨胀-a_nodata -9999强制把NoData设为-9999这是后续避免边界0值巨坑的第一道防线-co COMPRESSDEFLATE用无损压缩30米全省数据压完通常能缩到200MB以内。如果分幅文件分散在多个子目录把文件列表写进一个txt用gdalbuildvrt -input_file_list list.txt指定避免通配符匹配到无关文件。在ArcGIS里对应的是Mosaic To New Raster工具。注意把“Pixel Type”选为16 Bit Signed“NoData”填-9999Mosaic Operator选First或Blend都行关键在于每个分幅的边缘值要一致。分幅重叠区如果有高程差可以用Mosaic Tool底下的Blend做羽化但30米省级数据我一般不用缝隙小到对分析没影响。3.2 按河北省界裁剪别让边缘NoData污染分析裁剪看起来是最简单的操作但水文分析对边缘极其敏感。直接用省界shp裁剪后边界外的NoData会在坡度坡向计算时参与邻域运算导致边界一圈出现坡度假高值甚至把流向引向栅格外部。我常用的做法是裁剪前先判断后续做什么。如果只是出图直接剪到省界没问题如果要做汇水分析或流向计算先用Buffer把省界向外扩10个像元300米用这个缓冲后的范围裁剪分析完成后再把结果裁回省界。命令行裁剪gdalwarp -cutline hebei_buf.shp -crop_to_cutline \ -tr 30 30 -r near -dstnodata -9999 \ heb_dem_merged.tif heb_dem_bufclip.tif这里-tr 30 30强制像元尺寸不变形避免GDAL自动重采样把30米悄悄变成几十米-r near对整型高程用最近邻重采样保持原始高程值不插值。裁剪完之后用2.3节的脚本再跑一次值域统计重点确认两件事一是裁剪后范围的最外层一圈是否还有-9999二是有效高程最小值和原始数据是否一致。不一致就说明裁剪时发生了重采样回去查cutline有没有问题。3.3 投影坐标系为什么省级分析我更推荐Albers等积圆锥很多新手拿到DEM直接就用WGS84地理坐标系做坡度分析这在河北这种中纬度地区不是不能用但算面积、算坡度后处理时会有畸变。河北地处北纬36°到43°之间经度跨度约6.5°一个3度带覆盖不了全省所以省级尺度的分析我更推荐Albers等积圆锥投影中央经线选117°E双标准纬线选37°N和41°N这套参数能把河北省绝大部分区域的面积变形控制在可忽略范围。ArcGIS和GDAL里都可以自定义这套投影。要注意的是转换前用gdalwarp -r bilinear还是-r near对DEM来说我坚持用Near坡度分析里高程不能靠插值“发明”出中间值。投影转换后重新检查像元尺寸Albers投影下30米像元不会正好是30.00整出现30.00x29.98之类很正常只要偏差小于1%不用重新采样。如果要发布Web地图最后再转成Web墨卡托EPSG:3857供切片显示分析工作流里不要用Web墨卡托算面积和坡度不然河北北部的面积会被放大得离谱。4. 地貌分析实战坡度、坡向、填洼与河网提取的参数怎么定4.1 坡度坡向提取像元邻域、Z因子和单位ArcGIS的Slope工具默认用3×3像元邻域计算每个像元的坡度这个窗口对30米DEM意味着每个坡度值由周围270米范围内的地形决定局部微地形被天然平滑。执行前要注意三个参数输出单位选degree还是percent水文领域习惯用度Z factor当水平单位和高程单位都是米时必须设为1如果高程是英尺而水平是米Z factor要填0.3048方向反了会得到灾难性的错误结果第三个是输出栅格类型选浮点坡度不要存整型否则1.7°会被取整成2°平原区分析直接失掉一多半信息。坡向Aspect同样基于3×3邻域输出的是0到360度的方位角其中-1表示平地。坡向是环形数据后续做统计时不能直接求平均东北坡的平均值不是45和90取均值得到67.5这种算法要用三角函数转换后再平均。这个细节在山地生态分析里经常被忽略但在河北省的坡度分级和光照分析中很常见。提取完之后用栅格计算器把坡度小于某个阈值的像元挑出来看分布比直接看渲染图靠谱。河北平原区坡度大量集中在0°到1°之间这不代表数据坏了而是30米像元对微地貌的平滑结果。4.2 填洼不是无脑填先量化“填了多深”水文分析前必须做填洼Fill但填洼是最容易翻车的一步。ArcGIS的Fill工具默认会把所有检测到的洼地填平Z limit参数可以限制最大填深但默认行为在不同版本间语义有差异我不建议新手一上来就设某个经验值而是先跑一次默认填洼然后对比填洼前后差异import rasterio import numpy as np with rasterio.open(dem_filled.tif) as ds: filled ds.read(1) with rasterio.open(heb_dem_bufclip.tif) as ds: orig ds.read(1) depth filled.astype(np.float32) - orig.astype(np.float32) depth[orig -9999] np.nan print(最大填深(m):, np.nanmax(depth)) print(平均填深(m):, np.nanmean(depth)) print(填深超过5m的像元数:, np.sum(depth 5))这段脚本的价值在于把“填了多少”量化出来。如果最大填深超过20米、平均填深超过5米就要警觉要么数据源在山区有真实洼地喀斯特漏斗、采煤塌陷区要么原始DEM噪声太大。河北省太行山区有采空塌陷区平原区有古河道洼地真实洼地不该被填平。处理办法是分两步第一步用Z limit10限制最大填深只填浅层伪洼地第二步把填洼结果和高程影像叠加目视检查那些被填掉的深洼地是不是地形真相。4.3 河网提取全流程与阈值试错法ArcGIS/QGIS对照河网提取的标准链路是填洼DEM → 流向Flow Direction → 累积流量Flow Accumulation → 阈值提取 → 河网矢量。ArcGIS里每一步都有对应工具QGIS里则通过SAGA。流程和参数对照如下步骤ArcGIS工具关键参数QGISSAGA对应填洼FillZ limit视探查结果定Fill SinksWang Liu流向Flow DirectionD8算法默认即可SAGA Catchment Area累积流量Flow Accumulation输出浮点型Flow AccumulationParallelizable阈值提取栅格计算器Con(acc≥阈值, 1)SAGA Channel Definition转矢量Stream to Feature简化、平滑按需SAGA Channels阈值怎么定是河网提取里最像“玄学”的部分。我一般用试错法先按总像元数的万分之一到十万分之一作为初次阈值。河北全省约2.1亿像元万分之一对应2.1万像元适用于提取主要河流十万分之一对应2100像元能提出来大量沟壑。然后把这个结果和真实水系叠加调整阈值直到干流不破碎、支流不泛滥。平原区阈值要调高因为平坦地区累积流量随机性大太行山、燕山区阈值可以调低地形约束更强提取的河网可信度更高。还有一个容易被忽略的操作流向计算前如果DEM边缘还有NoData流向会朝无数据区域倾斜导致边界河网扭曲。所以第3章的外扩裁剪在这里就体现出价值了。QGIS用户要注意SAGA工具里的“Fill Sinks”默认会连真实洼地一起填处理逻辑和ArcGIS不完全一致对比结果时别拿两套结果直接PK要看相对真实水系的一致性。5. 避坑指南河北DEM处理中容易翻车的五个现场5.1 平原区坡度渲染成“一块饼”看着像坏数据现象把河北省DEM放进ArcGIS算坡度渲染出来华北平原一大片全是同一个颜色坡度像被“拍平”了甚至怀疑数据源出了问题。原因30米像元在平原区对微地形的分辨率只有零点几米大片区域坡度确实在0°到1°之间默认色带把1°和0.5°都映射成相近颜色视觉上会糊成一片。解决别改坡度数据改渲染方式。用山体阴影做底图叠加或者改用局部地形起伏指数TPI表达平缓地区的微地形如果一定用坡度做分级把分级边界设在0.5°、1°、2°、5°而不是等间距。再不行就用等高线辅助判读30米DEM在平原区提取5米等高线完全可行。从那以后我做河北平原的分析默认不再用单一色带渲染坡度。5.2 省界外出现一圈0m巨坑现象裁剪后的DEM在省界外出现一圈黑色或白色的“0米”像元统计高程最小值直接变成0整个栅格的山体阴影在边界处塌陷。原因原始TIFF里NoData值没被正确识别或者裁剪工具把掩膜外的值直接写成了0。0被当作有效高程参与计算后坡度、坡向在水系边界处全部失真。解决裁剪前先确认原始NoData值裁剪后立刻用2.3节的脚本查值域。如果已经出现0值用栅格计算器把高程等于0的像元重分类为NoData然后重新做一次裁剪。GDAL里加-dstnodata -9999能避免大部分情况ArcGIS的Extract by Mask工具里要手动指定NoData值。这个坑我在刚接触DEM那年踩过之后每次裁剪都强制跑一遍值域统计。5.3 TIFF转DEM格式后高程值全乱套现象把GeoTIFF用ArcGIS右键Export Data导出成DEM/Grid格式后高程从三百多变成几万或者全变成127、128这类奇怪的数看起来像被“加了盐”。原因TIFF里存的是浮点型高程导出成Esri Grid整数格式时发生了一次不恰当的整型转换或者导出过程中NoData值被重新赋值导致原本的-9999被当成有效高程读出来。解决不要用右键导出用ArcToolbox里的Copy Raster工具在Environment设置里把“Pixel Type”显式指定为16 Bit SignedNoData指定为-9999。转完以后再随机抽查三五个已知高程点对比真实值别只看最小最大值就放行。这个坑在把DEM交给只认旧格式的合作方时特别常见转换后多花两分钟做一次抽查能省掉后面一整轮返工。5.4 平原河网提取出来打转、拐直角现象用填洼后的DEM提取河网平原区河道不走自然弯曲而是沿着道路、田埂拐直角甚至一段上游流向和下游流向完全相反成“打转”状态。原因30米DEM在平原根本区分不了车辙、田埂、人工沟渠和自然微地形。道路高出地面几十厘米在30米采样后就是一个“伪山脊”径流就沿着路网走了农田的畦埂形成方格状的“伪地形”流向被带着拐直角。解决提取河网前先对平原区DEM做一次低通滤波把短波长噪声压掉。用Focal Statistics做一个5×5均值滤波然后再填洼、算流向。另一个办法是提取后做“矢量修剪”拿真实水系线对提取的河网做缓冲匹配只保留距离真实水系一定范围内的分支。做完这两个操作平原河网还打转就把阈值调高提取更“干流化”的河网。5.5 太行山、燕山出现条带状黑缝现象山体阴影渲染后太行山和燕山的阴坡出现细长条带状的黑色缝隙这些区域没有高程值像被刀划了一道。原因这是ASTER GDEM系数据的典型特征陡峭山谷的立体像对形成天窗加上云遮挡在阴坡会留下连续的空洞。如果数据源是ASTER衍生品几乎绕不开这个问题。解决不要试图用插值硬填大缝隙。先用ALOS AW3D30或Copernicus GLO-30作为参考数据把同区域的高程贴进来再用gdal_fillnodata做小范围补洞gdal_fillnodata.py -md 10 -si 1 heb_dem_ast.tif heb_dem_fill.tif-md 10表示最大搜索半径为10个像元超过10个像元的空洞不要补——补出来也是“发明地形”-si 1是平滑插值。补完以后把填充区域单独做个掩膜后续在坡度分析中把这部分像元标记为“低置信度”比装作没有空洞要诚实得多。6. 进阶给DEM做一次精度体检和一个出图习惯6.1 用实测高程点算RMSE判断数据够不够格把DEM当黑匣子用之前我建议花半小时做一次精度体检。方法不复杂拿一批已知高程的控制点水准点、GPS实测点、工程勘察点都可以在ArcGIS里用Extract Multi Values to Points把DEM高程提取到点上然后算均方根误差RMSE sqrt(mean((z_dem - z_ref)^2))。对30米产品RMSE落在5~10米以内算正常超过15米就要怀疑配准或来源问题。体检结果在项目报告里写一句“基于某批次控制点DEM中误差约X米”比空口说“用30米DEM够用”有说服力得多。6.2 山体阴影叠加分层设色一个能直接用的出图组合河北地形出图我常用的组合是DEM分层设色在下层山体阴影在上层透明度40%到50%。山体阴影参数固定为方位角315°、太阳高度角45°、Z因子1这组参数对河北地貌出图不容易出现“半边黑、半边亮”的失衡。叠加后如果觉得阴坡过重把透明度往上调到60%或者改成Multiply混合模式视觉效果立刻立体。这个组合同时适用于ArcGIS和QGIS打印输出记得把分辨率设到300dpi以上。6.3 一个值得长期保留的处理习惯把第2.3节的体检脚本存成一个固定工具每次接手任何DEM数据第一件事就是跑一遍投影、值域、NoData、空洞占比然后把输出截图存进项目文件夹。这套流程救过我很多次最典型的一次是换了一份看起来没问题的新数据结果空洞占比3.2%分布全在蓄水区范围里差点把整个库区的汇水分析带偏。从那以后我要求自己无论多急DEM进项目先体检绝不跳步。希望这个习惯对你也有用。本文还有配套的精品资源点击获取 河北省30米DEM数据处理实战：从精度评估到河网提取
简介河北省30米分辨率DEM数据包是一份面向GIS从业者、测绘与地理分析人员的地理空间数据产品。数据源自ASTER GDEM V3版本以GeoTiff栅格格式存储海拔信息并采用WGS84坐标系可直接用于地形坡度、坡向分析洪水模拟、规划选址及三维景观渲染等场景。压缩包共10个文件总大小约229.09MB除核心的高程栅格文件外还附带河北省行政边界shp矢量文件及配套的dbf、prj、sbn等属性与投影辅助文件便于在ArcGIS等平台中与边界叠加使用省去自行配准与裁剪的步骤。发布至今已有698人学习下载数据版本较新时效性有保障。对于需要快速获取河北全境高分辨率地形数据并进行流域分析、交通选线等工作的用户这份打包好的DEM与边界文件组合可直接投入使用减少数据预处理工作量。1. 河北省DEM 30米数据先想清楚要拿它做什么拿到这份河北省DEM 30米分辨率数据第一件事不是急着拖进ArcGIS渲染成彩色而是先想清楚要拿它做什么。河北从坝上高原、燕山、太行山一路过渡到华北平原30米像元对这种省级尺度的分析刚刚好全省坡度分级、区域汇水分析、灾害易发性初筛都能扛住但要拿来算单边坡精细坡度或做乡镇级土方设计分辨率和精度都不够。这份数据适合GIS从业者、规划师、水利与测绘背景的工程师直接当底图和分析输入新手也能用它练完整的地形分析流程。动手之前先花五分钟确认投影、值域、NoData、空洞分布四个参数它们决定后面整个工作流要不要返工。2. 选型与元数据体检这份DEM凭什么能用2.1 30米分辨率在河北能扛住的场景和扛不住的任务先说一个直观的数字30米像元意味着每个栅格代表 900 平方米的地面。河北全省面积约 18.88 万平方公里一张覆盖全省的30米栅格大约有 2.1 亿个像元单波段Int16裸数据大概 400MB 出头加上金字塔和压缩后也就几百兆。如果换成10米数据像元数直接涨到 19 亿数据量接近10倍ArcGIS里做一次填洼和流向计算的时间也差不多翻倍。对区域级项目来说这个性价比非常关键。实用场景里30米DEM适合做的是省级国土空间规划里的坡度分级水利上的汇水区划分和河网初提地质灾害易发性的区域性初筛以及林业、农业的立地条件分类。这些任务的共同特点是“看趋势、看格局”单个像元的误差不会颠覆结论。反过来它扛不住的任务也很明确单边坡开挖设计、村庄尺度的排水管线分析、矿区沉降监测、工程土方量精算这些场合需要的是厘米级到米级数据30米DEM连“参考”都很勉强。具体适合不适合可以对照下面这张表来判断任务类型是否适合原因省级坡度分级适合30米能识别出大尺度地形趋势坡度分级误差可控流域汇水区划适合汇水面积通常以km²计像元尺度影响有限地质灾害区域初筛适合关注的是地形组合关系不是单点高程单边坡精细设计不适合900m²一个像元边坡边界完全无法表达村庄级排水分析不适合排水管线精度要求远高于30米土方量精算不适合地形细节被平滑挖填方误差会放大2.2 主流30米数据源对比SRTM、ASTER、ALOS与Copernicus市面上30米DEM的全球公开数据源就那么几个特性差异不小。如果你拿到手的数据没标注生产源先按这个表做个预判断数据源标称分辨率公开精度指标已知软肋河北适用性SRTM v31弧秒30米相对垂直精度约±6m高纬度覆盖到60°N为止河北没问题但在复杂山地有空洞平原区表现平滑适合ASTER GDEM v330米官方称±7~14m山区条带噪声、云遮蔽空洞明显太行山、燕山要检查空洞ALOS AW3D30 v330米由5米重采样相对精度约5m量级云覆盖区有缺失但比ASTER干净综合质量好山地表现优于ASTERCopernicus GLO-3030米垂直精度优于SRTM来源于TanDEM-X平面位置精度高噪声均匀适合做分析和出图经验上讲在河北这种地貌组合复杂的地方SRTM的平原高程误差小ASTER在山地容易出现细长条带状空洞ALOS整体干净但偶尔有云补丁Copernicus是后起之秀均匀性好。我一般拿到资源后不纠结于它属于哪个源直接跑一遍体检脚本看空洞分布和噪声形态就能猜个八九不离十。另外要提醒一句如果数据表面在建筑区和茂密林区明显隆起那你拿到的可能是DSM数字表面模型而不是DEM。DSM带屋顶和树冠直接做水文分析会出现径流在房顶和树梢上流的怪事这是很多新手把河网提得歪七扭八的隐性原因。2.3 数据到手先查四个参数投影、值域、NoData、空洞打开数据之前先写个检查脚本比在ArcGIS里肉眼刷屏高效得多。我用Python的rasterio来做几分钟能跑完import rasterio import numpy as np path heb_dem_30m.tif with rasterio.open(path) as ds: print(投影:, ds.crs) print(范围:, ds.bounds) print(像元尺寸:, ds.res) arr ds.read(1) nodata ds.nodata valid arr[arr ! nodata] total arr.size nodata_count total - valid.size print(有效高程min/max:, valid.min(), valid.max()) print(NoData占比: {:.2f}%.format(nodata_count / total * 100))这段脚本的输出要对着四个检查点看第一投影如果是WGS84地理坐标系后面做面积计算前必须转投影第二像元尺寸必须接近30×30实际出现29.99或30.00都正常但出现1或2开头的数说明数据被重采样过第三有效高程最小值如果是0或-32768基本可以断定边界或空洞被填了0值要回炉处理第四NoData占比超过2%就得把NoData像元分布图导出来看确认空洞是不是集中在太行山、燕山的阴坡。这四个参数过关再往下做坡度水文分析才踏实。空洞位置怎么看简单做法是把NoData重分类成1、有效值设成0用ArcGIS的Spatial Analyst画出来。如果空洞呈细长条带、走向和山脊一致这是典型的ASTER天窗特征如果呈块状可能是云或原始数据缺失。记住这个判断后面第5章补洞时能少走弯路。3. 把分幅数据变成河北省一张图镶嵌、裁剪与投影选择3.1 分幅TIFF镶嵌成全省栅格如果这份河北省DEM是按标准分幅交付的第一步是把几十个分幅拼成一张全省图。命令行环境里我习惯先用GDAL把分幅“虚拟拼接”成VRT再做一次实际转换这样中途出错不用重新读几十个文件# 先生成VRT虚拟栅格不实际复制数据 gdalbuildvrt heb_dem_merged.vrt *.tif # 再转成正式TIFF统一像元类型和NoData gdal_translate -ot Int16 -a_nodata -9999 -co COMPRESSDEFLATE heb_dem_merged.vrt heb_dem_merged.tif关键参数在第二行-ot Int16把栅格统一成16位有符号整型防止某个分幅是浮点型导致整张图体积膨胀-a_nodata -9999强制把NoData设为-9999这是后续避免边界0值巨坑的第一道防线-co COMPRESSDEFLATE用无损压缩30米全省数据压完通常能缩到200MB以内。如果分幅文件分散在多个子目录把文件列表写进一个txt用gdalbuildvrt -input_file_list list.txt指定避免通配符匹配到无关文件。在ArcGIS里对应的是Mosaic To New Raster工具。注意把“Pixel Type”选为16 Bit Signed“NoData”填-9999Mosaic Operator选First或Blend都行关键在于每个分幅的边缘值要一致。分幅重叠区如果有高程差可以用Mosaic Tool底下的Blend做羽化但30米省级数据我一般不用缝隙小到对分析没影响。3.2 按河北省界裁剪别让边缘NoData污染分析裁剪看起来是最简单的操作但水文分析对边缘极其敏感。直接用省界shp裁剪后边界外的NoData会在坡度坡向计算时参与邻域运算导致边界一圈出现坡度假高值甚至把流向引向栅格外部。我常用的做法是裁剪前先判断后续做什么。如果只是出图直接剪到省界没问题如果要做汇水分析或流向计算先用Buffer把省界向外扩10个像元300米用这个缓冲后的范围裁剪分析完成后再把结果裁回省界。命令行裁剪gdalwarp -cutline hebei_buf.shp -crop_to_cutline \ -tr 30 30 -r near -dstnodata -9999 \ heb_dem_merged.tif heb_dem_bufclip.tif这里-tr 30 30强制像元尺寸不变形避免GDAL自动重采样把30米悄悄变成几十米-r near对整型高程用最近邻重采样保持原始高程值不插值。裁剪完之后用2.3节的脚本再跑一次值域统计重点确认两件事一是裁剪后范围的最外层一圈是否还有-9999二是有效高程最小值和原始数据是否一致。不一致就说明裁剪时发生了重采样回去查cutline有没有问题。3.3 投影坐标系为什么省级分析我更推荐Albers等积圆锥很多新手拿到DEM直接就用WGS84地理坐标系做坡度分析这在河北这种中纬度地区不是不能用但算面积、算坡度后处理时会有畸变。河北地处北纬36°到43°之间经度跨度约6.5°一个3度带覆盖不了全省所以省级尺度的分析我更推荐Albers等积圆锥投影中央经线选117°E双标准纬线选37°N和41°N这套参数能把河北省绝大部分区域的面积变形控制在可忽略范围。ArcGIS和GDAL里都可以自定义这套投影。要注意的是转换前用gdalwarp -r bilinear还是-r near对DEM来说我坚持用Near坡度分析里高程不能靠插值“发明”出中间值。投影转换后重新检查像元尺寸Albers投影下30米像元不会正好是30.00整出现30.00x29.98之类很正常只要偏差小于1%不用重新采样。如果要发布Web地图最后再转成Web墨卡托EPSG:3857供切片显示分析工作流里不要用Web墨卡托算面积和坡度不然河北北部的面积会被放大得离谱。4. 地貌分析实战坡度、坡向、填洼与河网提取的参数怎么定4.1 坡度坡向提取像元邻域、Z因子和单位ArcGIS的Slope工具默认用3×3像元邻域计算每个像元的坡度这个窗口对30米DEM意味着每个坡度值由周围270米范围内的地形决定局部微地形被天然平滑。执行前要注意三个参数输出单位选degree还是percent水文领域习惯用度Z factor当水平单位和高程单位都是米时必须设为1如果高程是英尺而水平是米Z factor要填0.3048方向反了会得到灾难性的错误结果第三个是输出栅格类型选浮点坡度不要存整型否则1.7°会被取整成2°平原区分析直接失掉一多半信息。坡向Aspect同样基于3×3邻域输出的是0到360度的方位角其中-1表示平地。坡向是环形数据后续做统计时不能直接求平均东北坡的平均值不是45和90取均值得到67.5这种算法要用三角函数转换后再平均。这个细节在山地生态分析里经常被忽略但在河北省的坡度分级和光照分析中很常见。提取完之后用栅格计算器把坡度小于某个阈值的像元挑出来看分布比直接看渲染图靠谱。河北平原区坡度大量集中在0°到1°之间这不代表数据坏了而是30米像元对微地貌的平滑结果。4.2 填洼不是无脑填先量化“填了多深”水文分析前必须做填洼Fill但填洼是最容易翻车的一步。ArcGIS的Fill工具默认会把所有检测到的洼地填平Z limit参数可以限制最大填深但默认行为在不同版本间语义有差异我不建议新手一上来就设某个经验值而是先跑一次默认填洼然后对比填洼前后差异import rasterio import numpy as np with rasterio.open(dem_filled.tif) as ds: filled ds.read(1) with rasterio.open(heb_dem_bufclip.tif) as ds: orig ds.read(1) depth filled.astype(np.float32) - orig.astype(np.float32) depth[orig -9999] np.nan print(最大填深(m):, np.nanmax(depth)) print(平均填深(m):, np.nanmean(depth)) print(填深超过5m的像元数:, np.sum(depth 5))这段脚本的价值在于把“填了多少”量化出来。如果最大填深超过20米、平均填深超过5米就要警觉要么数据源在山区有真实洼地喀斯特漏斗、采煤塌陷区要么原始DEM噪声太大。河北省太行山区有采空塌陷区平原区有古河道洼地真实洼地不该被填平。处理办法是分两步第一步用Z limit10限制最大填深只填浅层伪洼地第二步把填洼结果和高程影像叠加目视检查那些被填掉的深洼地是不是地形真相。4.3 河网提取全流程与阈值试错法ArcGIS/QGIS对照河网提取的标准链路是填洼DEM → 流向Flow Direction → 累积流量Flow Accumulation → 阈值提取 → 河网矢量。ArcGIS里每一步都有对应工具QGIS里则通过SAGA。流程和参数对照如下步骤ArcGIS工具关键参数QGISSAGA对应填洼FillZ limit视探查结果定Fill SinksWang Liu流向Flow DirectionD8算法默认即可SAGA Catchment Area累积流量Flow Accumulation输出浮点型Flow AccumulationParallelizable阈值提取栅格计算器Con(acc≥阈值, 1)SAGA Channel Definition转矢量Stream to Feature简化、平滑按需SAGA Channels阈值怎么定是河网提取里最像“玄学”的部分。我一般用试错法先按总像元数的万分之一到十万分之一作为初次阈值。河北全省约2.1亿像元万分之一对应2.1万像元适用于提取主要河流十万分之一对应2100像元能提出来大量沟壑。然后把这个结果和真实水系叠加调整阈值直到干流不破碎、支流不泛滥。平原区阈值要调高因为平坦地区累积流量随机性大太行山、燕山区阈值可以调低地形约束更强提取的河网可信度更高。还有一个容易被忽略的操作流向计算前如果DEM边缘还有NoData流向会朝无数据区域倾斜导致边界河网扭曲。所以第3章的外扩裁剪在这里就体现出价值了。QGIS用户要注意SAGA工具里的“Fill Sinks”默认会连真实洼地一起填处理逻辑和ArcGIS不完全一致对比结果时别拿两套结果直接PK要看相对真实水系的一致性。5. 避坑指南河北DEM处理中容易翻车的五个现场5.1 平原区坡度渲染成“一块饼”看着像坏数据现象把河北省DEM放进ArcGIS算坡度渲染出来华北平原一大片全是同一个颜色坡度像被“拍平”了甚至怀疑数据源出了问题。原因30米像元在平原区对微地形的分辨率只有零点几米大片区域坡度确实在0°到1°之间默认色带把1°和0.5°都映射成相近颜色视觉上会糊成一片。解决别改坡度数据改渲染方式。用山体阴影做底图叠加或者改用局部地形起伏指数TPI表达平缓地区的微地形如果一定用坡度做分级把分级边界设在0.5°、1°、2°、5°而不是等间距。再不行就用等高线辅助判读30米DEM在平原区提取5米等高线完全可行。从那以后我做河北平原的分析默认不再用单一色带渲染坡度。5.2 省界外出现一圈0m巨坑现象裁剪后的DEM在省界外出现一圈黑色或白色的“0米”像元统计高程最小值直接变成0整个栅格的山体阴影在边界处塌陷。原因原始TIFF里NoData值没被正确识别或者裁剪工具把掩膜外的值直接写成了0。0被当作有效高程参与计算后坡度、坡向在水系边界处全部失真。解决裁剪前先确认原始NoData值裁剪后立刻用2.3节的脚本查值域。如果已经出现0值用栅格计算器把高程等于0的像元重分类为NoData然后重新做一次裁剪。GDAL里加-dstnodata -9999能避免大部分情况ArcGIS的Extract by Mask工具里要手动指定NoData值。这个坑我在刚接触DEM那年踩过之后每次裁剪都强制跑一遍值域统计。5.3 TIFF转DEM格式后高程值全乱套现象把GeoTIFF用ArcGIS右键Export Data导出成DEM/Grid格式后高程从三百多变成几万或者全变成127、128这类奇怪的数看起来像被“加了盐”。原因TIFF里存的是浮点型高程导出成Esri Grid整数格式时发生了一次不恰当的整型转换或者导出过程中NoData值被重新赋值导致原本的-9999被当成有效高程读出来。解决不要用右键导出用ArcToolbox里的Copy Raster工具在Environment设置里把“Pixel Type”显式指定为16 Bit SignedNoData指定为-9999。转完以后再随机抽查三五个已知高程点对比真实值别只看最小最大值就放行。这个坑在把DEM交给只认旧格式的合作方时特别常见转换后多花两分钟做一次抽查能省掉后面一整轮返工。5.4 平原河网提取出来打转、拐直角现象用填洼后的DEM提取河网平原区河道不走自然弯曲而是沿着道路、田埂拐直角甚至一段上游流向和下游流向完全相反成“打转”状态。原因30米DEM在平原根本区分不了车辙、田埂、人工沟渠和自然微地形。道路高出地面几十厘米在30米采样后就是一个“伪山脊”径流就沿着路网走了农田的畦埂形成方格状的“伪地形”流向被带着拐直角。解决提取河网前先对平原区DEM做一次低通滤波把短波长噪声压掉。用Focal Statistics做一个5×5均值滤波然后再填洼、算流向。另一个办法是提取后做“矢量修剪”拿真实水系线对提取的河网做缓冲匹配只保留距离真实水系一定范围内的分支。做完这两个操作平原河网还打转就把阈值调高提取更“干流化”的河网。5.5 太行山、燕山出现条带状黑缝现象山体阴影渲染后太行山和燕山的阴坡出现细长条带状的黑色缝隙这些区域没有高程值像被刀划了一道。原因这是ASTER GDEM系数据的典型特征陡峭山谷的立体像对形成天窗加上云遮挡在阴坡会留下连续的空洞。如果数据源是ASTER衍生品几乎绕不开这个问题。解决不要试图用插值硬填大缝隙。先用ALOS AW3D30或Copernicus GLO-30作为参考数据把同区域的高程贴进来再用gdal_fillnodata做小范围补洞gdal_fillnodata.py -md 10 -si 1 heb_dem_ast.tif heb_dem_fill.tif-md 10表示最大搜索半径为10个像元超过10个像元的空洞不要补——补出来也是“发明地形”-si 1是平滑插值。补完以后把填充区域单独做个掩膜后续在坡度分析中把这部分像元标记为“低置信度”比装作没有空洞要诚实得多。6. 进阶给DEM做一次精度体检和一个出图习惯6.1 用实测高程点算RMSE判断数据够不够格把DEM当黑匣子用之前我建议花半小时做一次精度体检。方法不复杂拿一批已知高程的控制点水准点、GPS实测点、工程勘察点都可以在ArcGIS里用Extract Multi Values to Points把DEM高程提取到点上然后算均方根误差RMSE sqrt(mean((z_dem - z_ref)^2))。对30米产品RMSE落在5~10米以内算正常超过15米就要怀疑配准或来源问题。体检结果在项目报告里写一句“基于某批次控制点DEM中误差约X米”比空口说“用30米DEM够用”有说服力得多。6.2 山体阴影叠加分层设色一个能直接用的出图组合河北地形出图我常用的组合是DEM分层设色在下层山体阴影在上层透明度40%到50%。山体阴影参数固定为方位角315°、太阳高度角45°、Z因子1这组参数对河北地貌出图不容易出现“半边黑、半边亮”的失衡。叠加后如果觉得阴坡过重把透明度往上调到60%或者改成Multiply混合模式视觉效果立刻立体。这个组合同时适用于ArcGIS和QGIS打印输出记得把分辨率设到300dpi以上。6.3 一个值得长期保留的处理习惯把第2.3节的体检脚本存成一个固定工具每次接手任何DEM数据第一件事就是跑一遍投影、值域、NoData、空洞占比然后把输出截图存进项目文件夹。这套流程救过我很多次最典型的一次是换了一份看起来没问题的新数据结果空洞占比3.2%分布全在蓄水区范围里差点把整个库区的汇水分析带偏。从那以后我要求自己无论多急DEM进项目先体检绝不跳步。希望这个习惯对你也有用。本文还有配套的精品资源点击获取