尧图网页设计 
简介此套数据是覆盖河北省全域的30米分辨率数字高程模型数据源于ASTER GDEM V3采用GeoTiff格式并配以WGS84坐标系适合GIS使用者、地理分析人员及规划从业者用于地形分析、流域模拟、选址规划等专业场景也可服务于水利、地质、环境科学等领域的教学科研。包内共10个文件主文件HeBei_DEM_30m_ASTGTMV003.tif为高程栅格另含河北省域矢量边界shp及其属性表dbf、坐标投影prj、空间索引sbn/sbx/shx、栅格配准tfw等配套数据压缩包约229.09MB解压后即可载入ArcGIS、QGIS等软件开展分析目录结构清晰便于快速定位。目前已有698人学习下载。这份数据除提供30米高程栅格外还同步给出省界矢量文件配合使用可方便地进行坡度坡向提取、可视域分析、洪水淹没模拟、通信基站布点与工程选址为华北地区地学研究、工程实践及课堂教学提供了可靠的高精度地形底图。1. 河北省DEM30米分辨率不是一张图而是一条要自己拼完的数据链把河北省30米分辨率的DEM数据真正用起来的人多半不是想“看地图”。河北的地形从太行山、燕山一路跌到冀南平原省内高差超过两千米做坡度分析、流域划分、选址评估和土石方估算时30米栅格是性价比最好的起步分辨率比90米能看清沟谷山脊比12米又不需要动几十GB的数据。问题也在这里——网上没有一块现成的“河北DEM”拿到的都是按经度纬度切好的分幅文件要自己下载、检查、拼接、裁剪才能组成一套能进项目的河北省DEM数据。整条链路最容易让人翻车的不是下载那一关而是处理顺序和参数设置。2. 河北省30米DEM数据源选型四类来源怎么取舍、按用途挑哪个2.1 免费且稳定拿到河北省分幅DEM的四个来源做河北省全域分析常见的数据源就四种它们的原始分辨率都不小于30米但产品特性差异很大。我用过的免费来源主要是以下这些数据源名义分辨率覆盖与特点常见问题ASTER GDEM30米全球覆盖山区细节好可以覆盖到北纬83度平原噪声大城市区接近DSM分幅多SRTM 1弧秒30米覆盖60°N以南高程整体平稳接边好山谷和湿地有抬高森林区偏树冠顶ALOS AW3D3030米全球覆盖垂直精度在公开数据里口碑好分幅文件体积大下载麻烦国内镜像30米地理空间数据云等平台有上述数据镜像部分平台偶尔缺分幅要核对范围做“dem数据下载”时我一般先去国内镜像按分幅筛一遍优先选WGS84地理坐标的原始产品不要在下载阶段就先选了UTM投影版本。很多镜像平台上的文件名带经纬度网格号但不同平台命名规则不一样真正重要的不是文件名而是元数据里的分辨率、坐标系和值域。还有一种很容易踩的版本搜索“dem安装包”得到的东西别信。DEM数据没有“安装包”这种发布形式它永远是一堆分幅栅格文件谁给你一个“河北省DEM一键安装包”谁就是在拿预制数据骗你边界、投影和精度都没法追溯。2.2 按用途选数据源填洼、坡度、高程精度分别该信谁不同的应用对数据源的信任度完全不同这里没有哪个源能通吃。做水文分析和填洼优先SRTM。ASTER GDEM在冀南平原这种低起伏区域会有一堆伪洼地填洼时会填出大面积“锅底”最后河网提取出来的河道位置和历史河道对不上。SRTM整体更平滑伪洼地少填洼后的流向连续性明显更好。做坡度坡向分析优先ASTER。山区沟谷细节更锐利坡度分级结果在太行山、燕山区域和真实地形更贴。但ASTER在平原区的噪声无法忽视玉米地、农田边界会出现周期性的条纹坡度计算会把这种噪声放大所以平原区要配合平滑处理。做高程绝对值核对两个源都要怀疑。拿河北省已有的水准点检查SRTM整体偏差小但在河谷和湿地偏高ASTER在城市偏高。最好的做法是取两者逐点差值差值大的区域单独处理不要直接二选一。2.3 先算覆盖再下载河北省需要多少分幅、怎么批量挑河北省主体在36°N到43°N、113.5°E到120°E之间。按1°×1°分幅估算横跨7个经度带、7个纬度带最多五十个分幅以内还要剔除渤海里的纯海面分幅。下载前我会先用脚本把需要的分幅网格列出来再逐个检查平台里是否齐全。# 列出河北省1°×1°分幅的经纬度网格范围用于在数据平台挨个勾选 min_lat, max_lat 36, 43 min_lon, max_lon 113, 120 tiles [] for lat in range(min_lat, max_lat): for lon in range(min_lon, max_lon): tiles.append(fN{lat:02d}-E{lon:03d}) print(f分幅数量{len(tiles)}) print(tiles)这段代码只是把纬度带和经度带做一个笛卡尔积方便核对平台数据是否齐全。注意纬度范围从36到43是7个纬度带经度从113到120也是7个经度带总共49个网格其中渤海、天津和山西交界处的分幅要手动剔除。下载时建议把边缘多下一圈也就是河北省边界外围再扩展一个分幅避免后续省界缓冲区分析时缺数据。如果你同时在找山东省30米DEM思路完全一样。山东省经纬度范围大概是34.3°N到38.4°N、114.8°E到122.7°E分幅数量更少处理流程一点不用改只是边界文件和网格范围不同。3. 把分幅DEM整理成河北省30米DEM镶嵌、裁剪与投影一条路走完3.1 dem文件不是某个专有格式先检查栅格这颗“黑匣子”经常有人问“ArcGIS怎么把tiff转成dem文件”这个问法本身就有误区。DEM是数据内容不是文件后缀。GeoTIFF本身就是最常见的dem文件ArcGIS可以直接读你非要转成ESRI GRID格式用Copy Raster转一下就行但绝大多数场景没必要。下载完分幅数据第一步不是急着拼接而是做一次批量体检看每个分幅的像素类型、坐标系和值域是否正常。import arcpy arcpy.env.workspace rD:\hebei_dem\raw rasters arcpy.ListRasters(*) def get_prop(ras, prop): return arcpy.management.GetRasterProperties(ras, prop).getOutput(0) for ras in rasters: desc arcpy.Describe(ras) sr desc.spatialReference mn get_prop(ras, MINIMUM) mx get_prop(ras, MAXIMUM) cs get_prop(ras, CELLSIZE) print(ras, desc.pixelType, sr.name, mn, mx, cs)这段代码逐个分幅打印四个关键信息像素类型、坐标系、最小值、最大值和像元大小。正常高程分幅的像素类型一般是16位整型或32位浮点最小值是负几百米或0最大值是真实山峰高程。如果最小值是-9999或-32768这种“整齐的怪数”说明NoData值混进了有效统计后续处理会出问题。坐标系如果显示GCS_WGS_1984就先不要投影等镶嵌裁剪完成后统一做一次投影重采样减少插值次数。3.2 批量镶嵌接边不出现台阶的Mosaic参数组合分幅检查没问题后做镶嵌。用Mosaic To New Raster工具我通常把重叠区处理方式选BLEND而不是默认的LAST尤其是ASTER分幅接边有明显高差的时候BLEND能让重叠区渐变过渡避免出现一条亮边或暗边。inputs [rD:\hebei_dem\raw\{}.tif.format(r) for r in rasters] arcpy.management.MosaicToNewRaster( inputs, rD:\hebei_dem\out, hebei_mosaic.tif, arcpy.SpatialReference(WGS 1984), 32_BIT_FLOAT, 0.00027777778, 1, BLEND )这里几个参数必须说清楚。像素类型选32_BIT_FLOAT让高程值保留小数避免拼接后梯田化。CellSize传0.00027777778是30米在赤道附近的经纬度数值如果你下载的原始分幅就带着自己的像元大小最稳妥的办法是把原始CELLSIZE读取后原样填入。最后两个参数分别是波段数和镶嵌方法波段固定1方法用BLEND。如果你的分幅边缘有明确无效区域可以在镶嵌前先用Copy Raster统一NoData值为-9999否则不同来源的NoData值0、255、-9999会在拼接后混成一片。3.3 按河北省行政边界裁剪留一个像元的缓冲拼接完的栅格覆盖范围远大于河北省接下来裁剪。这里有个很多人忽视的细节直接用省界裁剪边界上的栅格容易丢半圈尤其是边界沿着山谷线走的时候裁剪出来的DEM边缘会出现锯齿状缺块。正确做法是先把省界向外Buffer一个像元再裁剪最后再把缓冲区多出来的区域用原始省界做一次掩膜处理。# 省界向外缓冲30米留出一个像元的余量 arcpy.analysis.Buffer( rD:\hebei_dem\hebei_boundary.shp, rD:\hebei_dem\hebei_buf.shp, 30 Meters ) # 用缓冲区面裁剪镶嵌栅格 arcpy.management.Clip( rD:\hebei_dem\out\hebei_mosaic.tif, rD:\hebei_dem\hebei_buf.shp, rD:\hebei_dem\hebei_clip.tif, #, -9999, ClippingGeometry )Clip参数里的-9999是写入NoData值要和镶嵌输出的NoData保持一致。ClippingGeometry表示按面要素的实际形状裁剪而不是按外接矩形。同时建议把环境设置里的Snap Raster指定为镶嵌后的原栅格让裁剪输出的格网和原栅格完全对齐避免重采样产生半个像元偏移。裁剪后第一时间检查属性表里的最小值和最大值如果突然出现一个很大的负数说明NoData没有正确写入。3.4 投影重采样为什么河北省不用UTM而用Albers等面积投影河北省横跨了UTM 50N和51N两个投影带用任何一个带的UTM都会让另一半扭曲。常见做法是选择WGS 1984 Albers等面积投影中央经线设在117°E双标准纬线取37°N和41°N附近这样整个省的面积变形都很小做坡度分级后的面积统计才有意义。arcpy.management.ProjectRaster( rD:\hebei_dem\hebei_clip.tif, rD:\hebei_dem\hebei_dem_30m_albers.tif, rD:\hebei_dem\wgs1984_albers.prj, BILINEAR, 30 )重采样方法这里选BILINEAR。高程是连续表面用NEAREST会让台阶化严重用CUBIC会产生过冲平地上容易造出小土包。BILINEAR是折中方案。投影输出的像元大小设为30注意投影后经纬度栅格里一个像元代表的地面大小是不均匀的只有在Albers这类投影下每个像元才真正对应地面30米乘30米。投影完成后河北省30米DEM才算在几何上可用了。4. 河北省DEM处理的五个坑现象、原因与解决办法4.1 太行山区出现条纹状负值掩膜没清干净现象保定、石家庄西侧太行山区域处理完的DEM上出现沿分幅边界走向的长条状负值看着像地形上被硬切了一道深沟实际地形根本没有这条沟。原因ASTER GDEM的分幅产品在云盖、雪盖和陡坎处有掩膜空洞下载后的分幅边缘经常残留一小块无效区。镶嵌时这些无效区如果没被正确识别为NoData就会以负值形式参与插值拼接后被放大成一条条“假深谷”。解决镶嵌前每个分幅单独做一次IsNull检查把无效像元统一设为NoData。已经拼出问题的用Focal Statistics窗口内中值填充这些区域然后和原始山峰边界做叠加检查确认没有把真实山谷填平。操作顺序是先修掩膜再镶嵌顺序反了后面要返工。4.2 石家庄周边高程整体偏高误把DSM当DEM用现象拿处理好的DEM和河北省已公开的高程控制点对比石家庄二环内普遍偏高1到2米城市外围很快就正常了。原因ASTER产品在城区保留了建筑顶和树冠高度本质上是数字表面模型DSM而不是纯地面高程。业内说的“dsm生成dem”就是用局部替换加邻域滤波把表面高度压回地面趋势这在局部小范围可行省域尺度直接套用山区沟谷会被填浅平原建筑区又压不到位。解决先查分幅元数据看产品描述里写的是DEM还是DSM。如果混入了DSM以SRTM为基准做逐点差值把高出SRTM超过阈值的城区像元替换为SRTM值再对替换边缘做一个轻量3×3平滑。不要对整个河北省做“一键DSM转DEM”那会在太行山造出一批错误地形。4.3 坝上高原坡向呈现条带纹理接边系统差被放大现象张家口张北、尚义一带坡向分析结果出现明显的南北向条带同一片平坦草原上坡向呈刷子状看着像人工拼图。原因分幅接边处存在0.1到0.3米的系统高程差。这个误差在高程上看几乎不算事但坡向是对邻域高差变化的敏感性极高接边处的微小错位被坡向算法放大成几十度的方向偏移正好沿接边形成条带。解决坡度小于0.5°的平坦区域坡向本身就没有物理意义直接设为NoData或单独归为一类。做法是先算坡度用SetNull把坡度过小的像元置空再做坡向分类。如果条带依然影响山区接边可以用3×3窗口对DEM做一次轻平滑后再重算坡向但代价是山脊线细节会被磨掉一点。4.4 河道呈现“脊形”假象水面高程被抬高现象在滹沱河、永定河河段做剖面线河道不但没有低于两岸反而中间凸起像一条土脊流向分析在河道里找不到连续的低谷。原因公开的光学雷达高程产品在水体上的反演信号不稳定宽河段和有植被的河滩尤其明显。SRTM在湿地和水面抬升是出了名的黑匣子ASTER也好不到哪去。这不是你处理错了是原始数据在河道里本身就不准。解决做水文分析前把已确认的河网矢量做缓冲区沿河道走向对缓冲区内的DEM做线性插值替换让河道恢复到“河谷低于两岸”的状态。更省事的办法是用水文专业站点发布的修正DEM比如HydroSHEDS之类的水文修正产品再把河北省范围内截取出来。千万别拿原始30米DEM的河道高程去做洪水淹没分析那会得出离谱的结论。4.5 裁剪后的边界出现黑色边条NoData值没统一现象按省界裁剪后河北省边界内侧出现一圈规则的黑条刚开始只看到一片黑线放大后发现和边界线完全重合。原因原始分幅的NoData值不是同一个数。有人裁剪时把NoData设成255有人设成0镶嵌时又没有统一导致边界处有些像元被当成有效高程裁剪时这些“伪有效值”被保留成了黑色区域。解决镶嵌前统一NoData。最稳的做法是在镶嵌后用栅格计算器执行Con(IsNull(dem), -9999, dem)把所有无效像元固定写成-9999。裁剪时参数里的NoData也填-9999裁剪后检查边界一圈像元的统计值如果最小值是-9999说明边界处理干净了。5. 把河北省30米DEM用起来坡度、坡向、地形湿度指数与河网提取5.1 坡度分级太行山与冀南平原要分开设阈值河北省DEM做坡度分析直接用统一分级会出问题。太行山和燕山的坡度峰段集中在15度到30度冀南平原的“坡度”基本是噪声很多地方明明是一马平川算出来却有3到5度。dem rD:\hebei_dem\hebei_dem_30m_albers.tif slope_deg arcpy.sa.Slope(dem, DEGREE) # 分级0-2°平地2-6°缓坡6-15°坡地15-25°陡坡25°山地 remap 0 2 0;2 6 1;6 15 2;15 25 3;25 90 4 slope_cls arcpy.sa.Reclassify(slope_deg, VALUE, remap, NODATA) slope_cls.save(rD:\hebei_dem\result\slope_class.tif)分级阈值不是固定的。平原农业区建议把第一级下限提到2度小于2度按平地处理山区项目可以把陡坡起点提到30度更贴合冀北山地的实际坡度分布。Reclassify的输入必须是原始坡度值域不要拿已经做过平滑的坡度去分级那会把真实陡坡削平。5.2 坡向计算的“平地区翻车”与掩膜办法坡向在平缓区是全靠噪声驱动的结果这点做河北全域坡向时特别明显。华北平原地势平缓坡向却四面八方乱跳这不是算法错了是坡向定义本身依赖邻域高差高差趋近于零时任何微小的DEM噪声都会主导方向。slope_deg arcpy.sa.Slope(dem, DEGREE) aspect_deg arcpy.sa.Aspect(dem) # 坡度小于0.5度的区域坡向没有物理意义直接置为NoData aspect_clean arcpy.sa.SetNull(slope_deg 0.5, aspect_deg) aspect_clean.save(rD:\hebei_dem\result\aspect_clean.tif)阈值0.5度是经验值。如果项目区域是张家口坝上那种相对平坦的高原阈值可以提到1度如果是太行山峡谷区0.3度都够。这样处理后的坡向分类图山区和坝上交界处才不会有那种“碎玻璃”噪声。5.3 地形湿度指数TWI的计算链TWI是很多土壤湿度、滑坡和生态分析的前置变量。它的公式是ln(单位宽度汇流面积除以坡度正切值)在ArcGIS里通过填洼、流向、汇流累积三步算出。dem_fill arcpy.sa.Fill(dem, 30) flow_dir arcpy.sa.FlowDirection(dem_fill) flow_acc arcpy.sa.FlowAccumulation(flow_dir) slope_rad arcpy.sa.Slope(dem_fill, RADIAN) tan_slope arcpy.sa.Tan(slope_rad) # 避免坡度为0时除零 tan_slope_safe arcpy.sa.Con(tan_slope 0.001, 0.001, tan_slope) # 简化TWI汇流面积用累积像元数×30米像元大小近似 twi arcpy.sa.Ln((flow_acc * 30.0 1.0) / tan_slope_safe) twi.save(rD:\hebei_dem\result\twi.tif)这里解释两点。第一FlowAccumulation先跑在填洼后的DEM上如果不填洼平原区的伪洼地会让水流方向断裂汇流累积完全乱掉。第二公式里flow_acc * 30用30米像元的单边长度做简化近似严格的地表汇流面积还要考虑坡度对实际面积的拉伸但省级分析尺度下这个简化足够稳定。5.4 提取河网阈值选择与省界缓冲区陷阱用汇流累积量提取河网是最常见的河网生成方式。阈值直接影响河网密度河北山区和平原要分别对待。# 先用1000像元阈值提取试算约等于0.9平方公里集水面积 stream_raw arcpy.sa.Con(flow_acc 1000, 1) stream_raw.save(rD:\hebei_dem\result\stream_threshold1000.tif) # 平原地区阈值过大河网会断头降到300到500再试 stream_plain arcpy.sa.Con(flow_acc 400, 1) stream_plain.save(rD:\hebei_dem\result\stream_threshold400.tif)阈值1000在太行山区能得到比较完整的水系骨架但到冀南平原很多支流在平原区汇流面积积累不够会断在半路。平原地区阈值降到400再看。更准的做法是拿已经公开的1比25万河网做对比把阈值调到提取河网和现有河网重合率最高为止。还有一个省级分析特有的陷阱只拿河北省边界裁剪后的DEM提取河网省界附近河流的上游汇水面积会被切掉河网在边界处明显变稀这不是算法问题是数据范围问题。正确的做法是前期把邻省数据拼进来或者下载数据时就刻意多覆盖周边几十公里等河网提取完再按河北省界裁剪结果。6. 投入前的最后一步填洼控制、差值验证与省界缓冲区习惯6.1 填洼别一上来就用默认参数ArcGIS的Fill工具默认把所有洼地填到没有低点为止这在华北平原非常危险。平原的伪洼地会被成片填平汇流累积变成一摊辫状水道提取的河网完全不能用。我会在填洼时限制最大填充深度用z_limit参数把填充量压到5米或10米以内只填掉噪声级的小洼地保留真实的地形洼地。# 限制最大填洼深度避免平原被成片填平 dem_fill_limit arcpy.sa.Fill(dem, 10)填完后做一次“填洼差值图”把填洼后的DEM减去原DEM如果某个区域的差值和真实地势完全不匹配而且连成大片说明参数设置过激需要降低z_limit或改用其他数据源。6.2 用差值验证法给最终DEM上保险处理完成不等于能用。我会拿SRTM和最终DEM做一次逐点差值检查均值和中位数是否异常。diff arcpy.sa.Minus(dem_fill_limit, srtm_mosaic) mean_diff arcpy.management.GetRasterProperties(diff, MEAN).getOutput(0) std_diff arcpy.management.GetRasterProperties(diff, STD).getOutput(0) print(平均差值, mean_diff, 标准差, std_diff)河北省全域范围内两者差值的均值应该落在正负1米以内标准差在山区可以放宽到5米左右。如果出现某个方向系统性偏差比如西部山区整体高5米东部平原整体低3米那就是镶嵌时某个分幅没对齐要回去查那个分幅的投影和数据版本。6.3 省级分析先问一句边界外带了吗我处理这类省级DEM的过程里返工最多的就是省界缓冲区问题。有一回做河北西部某流域的河网提取河北与山西交界处的河流全部断头查了三天才意识到不是算法问题而是裁剪时把山西境内的上游切掉了。那次之后我养成了一个习惯任何省级地形分析第一步先把省界往外扩展一百公里保证海河和滦河的完整上游都在数据处理范围内等所有分析做完再切回河北省界出图。30米分辨率的河北省DEM本身精度够做项目但能不能真正落地取决于你留给自己的验证余量。希望这篇笔记能帮你把下载、拼接、裁剪和验证的弯路一次走完。本文还有配套的精品资源点击获取 河北省30米DEM数据处理全流程：从分幅下载到镶嵌裁剪与验证
简介此套数据是覆盖河北省全域的30米分辨率数字高程模型数据源于ASTER GDEM V3采用GeoTiff格式并配以WGS84坐标系适合GIS使用者、地理分析人员及规划从业者用于地形分析、流域模拟、选址规划等专业场景也可服务于水利、地质、环境科学等领域的教学科研。包内共10个文件主文件HeBei_DEM_30m_ASTGTMV003.tif为高程栅格另含河北省域矢量边界shp及其属性表dbf、坐标投影prj、空间索引sbn/sbx/shx、栅格配准tfw等配套数据压缩包约229.09MB解压后即可载入ArcGIS、QGIS等软件开展分析目录结构清晰便于快速定位。目前已有698人学习下载。这份数据除提供30米高程栅格外还同步给出省界矢量文件配合使用可方便地进行坡度坡向提取、可视域分析、洪水淹没模拟、通信基站布点与工程选址为华北地区地学研究、工程实践及课堂教学提供了可靠的高精度地形底图。1. 河北省DEM30米分辨率不是一张图而是一条要自己拼完的数据链把河北省30米分辨率的DEM数据真正用起来的人多半不是想“看地图”。河北的地形从太行山、燕山一路跌到冀南平原省内高差超过两千米做坡度分析、流域划分、选址评估和土石方估算时30米栅格是性价比最好的起步分辨率比90米能看清沟谷山脊比12米又不需要动几十GB的数据。问题也在这里——网上没有一块现成的“河北DEM”拿到的都是按经度纬度切好的分幅文件要自己下载、检查、拼接、裁剪才能组成一套能进项目的河北省DEM数据。整条链路最容易让人翻车的不是下载那一关而是处理顺序和参数设置。2. 河北省30米DEM数据源选型四类来源怎么取舍、按用途挑哪个2.1 免费且稳定拿到河北省分幅DEM的四个来源做河北省全域分析常见的数据源就四种它们的原始分辨率都不小于30米但产品特性差异很大。我用过的免费来源主要是以下这些数据源名义分辨率覆盖与特点常见问题ASTER GDEM30米全球覆盖山区细节好可以覆盖到北纬83度平原噪声大城市区接近DSM分幅多SRTM 1弧秒30米覆盖60°N以南高程整体平稳接边好山谷和湿地有抬高森林区偏树冠顶ALOS AW3D3030米全球覆盖垂直精度在公开数据里口碑好分幅文件体积大下载麻烦国内镜像30米地理空间数据云等平台有上述数据镜像部分平台偶尔缺分幅要核对范围做“dem数据下载”时我一般先去国内镜像按分幅筛一遍优先选WGS84地理坐标的原始产品不要在下载阶段就先选了UTM投影版本。很多镜像平台上的文件名带经纬度网格号但不同平台命名规则不一样真正重要的不是文件名而是元数据里的分辨率、坐标系和值域。还有一种很容易踩的版本搜索“dem安装包”得到的东西别信。DEM数据没有“安装包”这种发布形式它永远是一堆分幅栅格文件谁给你一个“河北省DEM一键安装包”谁就是在拿预制数据骗你边界、投影和精度都没法追溯。2.2 按用途选数据源填洼、坡度、高程精度分别该信谁不同的应用对数据源的信任度完全不同这里没有哪个源能通吃。做水文分析和填洼优先SRTM。ASTER GDEM在冀南平原这种低起伏区域会有一堆伪洼地填洼时会填出大面积“锅底”最后河网提取出来的河道位置和历史河道对不上。SRTM整体更平滑伪洼地少填洼后的流向连续性明显更好。做坡度坡向分析优先ASTER。山区沟谷细节更锐利坡度分级结果在太行山、燕山区域和真实地形更贴。但ASTER在平原区的噪声无法忽视玉米地、农田边界会出现周期性的条纹坡度计算会把这种噪声放大所以平原区要配合平滑处理。做高程绝对值核对两个源都要怀疑。拿河北省已有的水准点检查SRTM整体偏差小但在河谷和湿地偏高ASTER在城市偏高。最好的做法是取两者逐点差值差值大的区域单独处理不要直接二选一。2.3 先算覆盖再下载河北省需要多少分幅、怎么批量挑河北省主体在36°N到43°N、113.5°E到120°E之间。按1°×1°分幅估算横跨7个经度带、7个纬度带最多五十个分幅以内还要剔除渤海里的纯海面分幅。下载前我会先用脚本把需要的分幅网格列出来再逐个检查平台里是否齐全。# 列出河北省1°×1°分幅的经纬度网格范围用于在数据平台挨个勾选 min_lat, max_lat 36, 43 min_lon, max_lon 113, 120 tiles [] for lat in range(min_lat, max_lat): for lon in range(min_lon, max_lon): tiles.append(fN{lat:02d}-E{lon:03d}) print(f分幅数量{len(tiles)}) print(tiles)这段代码只是把纬度带和经度带做一个笛卡尔积方便核对平台数据是否齐全。注意纬度范围从36到43是7个纬度带经度从113到120也是7个经度带总共49个网格其中渤海、天津和山西交界处的分幅要手动剔除。下载时建议把边缘多下一圈也就是河北省边界外围再扩展一个分幅避免后续省界缓冲区分析时缺数据。如果你同时在找山东省30米DEM思路完全一样。山东省经纬度范围大概是34.3°N到38.4°N、114.8°E到122.7°E分幅数量更少处理流程一点不用改只是边界文件和网格范围不同。3. 把分幅DEM整理成河北省30米DEM镶嵌、裁剪与投影一条路走完3.1 dem文件不是某个专有格式先检查栅格这颗“黑匣子”经常有人问“ArcGIS怎么把tiff转成dem文件”这个问法本身就有误区。DEM是数据内容不是文件后缀。GeoTIFF本身就是最常见的dem文件ArcGIS可以直接读你非要转成ESRI GRID格式用Copy Raster转一下就行但绝大多数场景没必要。下载完分幅数据第一步不是急着拼接而是做一次批量体检看每个分幅的像素类型、坐标系和值域是否正常。import arcpy arcpy.env.workspace rD:\hebei_dem\raw rasters arcpy.ListRasters(*) def get_prop(ras, prop): return arcpy.management.GetRasterProperties(ras, prop).getOutput(0) for ras in rasters: desc arcpy.Describe(ras) sr desc.spatialReference mn get_prop(ras, MINIMUM) mx get_prop(ras, MAXIMUM) cs get_prop(ras, CELLSIZE) print(ras, desc.pixelType, sr.name, mn, mx, cs)这段代码逐个分幅打印四个关键信息像素类型、坐标系、最小值、最大值和像元大小。正常高程分幅的像素类型一般是16位整型或32位浮点最小值是负几百米或0最大值是真实山峰高程。如果最小值是-9999或-32768这种“整齐的怪数”说明NoData值混进了有效统计后续处理会出问题。坐标系如果显示GCS_WGS_1984就先不要投影等镶嵌裁剪完成后统一做一次投影重采样减少插值次数。3.2 批量镶嵌接边不出现台阶的Mosaic参数组合分幅检查没问题后做镶嵌。用Mosaic To New Raster工具我通常把重叠区处理方式选BLEND而不是默认的LAST尤其是ASTER分幅接边有明显高差的时候BLEND能让重叠区渐变过渡避免出现一条亮边或暗边。inputs [rD:\hebei_dem\raw\{}.tif.format(r) for r in rasters] arcpy.management.MosaicToNewRaster( inputs, rD:\hebei_dem\out, hebei_mosaic.tif, arcpy.SpatialReference(WGS 1984), 32_BIT_FLOAT, 0.00027777778, 1, BLEND )这里几个参数必须说清楚。像素类型选32_BIT_FLOAT让高程值保留小数避免拼接后梯田化。CellSize传0.00027777778是30米在赤道附近的经纬度数值如果你下载的原始分幅就带着自己的像元大小最稳妥的办法是把原始CELLSIZE读取后原样填入。最后两个参数分别是波段数和镶嵌方法波段固定1方法用BLEND。如果你的分幅边缘有明确无效区域可以在镶嵌前先用Copy Raster统一NoData值为-9999否则不同来源的NoData值0、255、-9999会在拼接后混成一片。3.3 按河北省行政边界裁剪留一个像元的缓冲拼接完的栅格覆盖范围远大于河北省接下来裁剪。这里有个很多人忽视的细节直接用省界裁剪边界上的栅格容易丢半圈尤其是边界沿着山谷线走的时候裁剪出来的DEM边缘会出现锯齿状缺块。正确做法是先把省界向外Buffer一个像元再裁剪最后再把缓冲区多出来的区域用原始省界做一次掩膜处理。# 省界向外缓冲30米留出一个像元的余量 arcpy.analysis.Buffer( rD:\hebei_dem\hebei_boundary.shp, rD:\hebei_dem\hebei_buf.shp, 30 Meters ) # 用缓冲区面裁剪镶嵌栅格 arcpy.management.Clip( rD:\hebei_dem\out\hebei_mosaic.tif, rD:\hebei_dem\hebei_buf.shp, rD:\hebei_dem\hebei_clip.tif, #, -9999, ClippingGeometry )Clip参数里的-9999是写入NoData值要和镶嵌输出的NoData保持一致。ClippingGeometry表示按面要素的实际形状裁剪而不是按外接矩形。同时建议把环境设置里的Snap Raster指定为镶嵌后的原栅格让裁剪输出的格网和原栅格完全对齐避免重采样产生半个像元偏移。裁剪后第一时间检查属性表里的最小值和最大值如果突然出现一个很大的负数说明NoData没有正确写入。3.4 投影重采样为什么河北省不用UTM而用Albers等面积投影河北省横跨了UTM 50N和51N两个投影带用任何一个带的UTM都会让另一半扭曲。常见做法是选择WGS 1984 Albers等面积投影中央经线设在117°E双标准纬线取37°N和41°N附近这样整个省的面积变形都很小做坡度分级后的面积统计才有意义。arcpy.management.ProjectRaster( rD:\hebei_dem\hebei_clip.tif, rD:\hebei_dem\hebei_dem_30m_albers.tif, rD:\hebei_dem\wgs1984_albers.prj, BILINEAR, 30 )重采样方法这里选BILINEAR。高程是连续表面用NEAREST会让台阶化严重用CUBIC会产生过冲平地上容易造出小土包。BILINEAR是折中方案。投影输出的像元大小设为30注意投影后经纬度栅格里一个像元代表的地面大小是不均匀的只有在Albers这类投影下每个像元才真正对应地面30米乘30米。投影完成后河北省30米DEM才算在几何上可用了。4. 河北省DEM处理的五个坑现象、原因与解决办法4.1 太行山区出现条纹状负值掩膜没清干净现象保定、石家庄西侧太行山区域处理完的DEM上出现沿分幅边界走向的长条状负值看着像地形上被硬切了一道深沟实际地形根本没有这条沟。原因ASTER GDEM的分幅产品在云盖、雪盖和陡坎处有掩膜空洞下载后的分幅边缘经常残留一小块无效区。镶嵌时这些无效区如果没被正确识别为NoData就会以负值形式参与插值拼接后被放大成一条条“假深谷”。解决镶嵌前每个分幅单独做一次IsNull检查把无效像元统一设为NoData。已经拼出问题的用Focal Statistics窗口内中值填充这些区域然后和原始山峰边界做叠加检查确认没有把真实山谷填平。操作顺序是先修掩膜再镶嵌顺序反了后面要返工。4.2 石家庄周边高程整体偏高误把DSM当DEM用现象拿处理好的DEM和河北省已公开的高程控制点对比石家庄二环内普遍偏高1到2米城市外围很快就正常了。原因ASTER产品在城区保留了建筑顶和树冠高度本质上是数字表面模型DSM而不是纯地面高程。业内说的“dsm生成dem”就是用局部替换加邻域滤波把表面高度压回地面趋势这在局部小范围可行省域尺度直接套用山区沟谷会被填浅平原建筑区又压不到位。解决先查分幅元数据看产品描述里写的是DEM还是DSM。如果混入了DSM以SRTM为基准做逐点差值把高出SRTM超过阈值的城区像元替换为SRTM值再对替换边缘做一个轻量3×3平滑。不要对整个河北省做“一键DSM转DEM”那会在太行山造出一批错误地形。4.3 坝上高原坡向呈现条带纹理接边系统差被放大现象张家口张北、尚义一带坡向分析结果出现明显的南北向条带同一片平坦草原上坡向呈刷子状看着像人工拼图。原因分幅接边处存在0.1到0.3米的系统高程差。这个误差在高程上看几乎不算事但坡向是对邻域高差变化的敏感性极高接边处的微小错位被坡向算法放大成几十度的方向偏移正好沿接边形成条带。解决坡度小于0.5°的平坦区域坡向本身就没有物理意义直接设为NoData或单独归为一类。做法是先算坡度用SetNull把坡度过小的像元置空再做坡向分类。如果条带依然影响山区接边可以用3×3窗口对DEM做一次轻平滑后再重算坡向但代价是山脊线细节会被磨掉一点。4.4 河道呈现“脊形”假象水面高程被抬高现象在滹沱河、永定河河段做剖面线河道不但没有低于两岸反而中间凸起像一条土脊流向分析在河道里找不到连续的低谷。原因公开的光学雷达高程产品在水体上的反演信号不稳定宽河段和有植被的河滩尤其明显。SRTM在湿地和水面抬升是出了名的黑匣子ASTER也好不到哪去。这不是你处理错了是原始数据在河道里本身就不准。解决做水文分析前把已确认的河网矢量做缓冲区沿河道走向对缓冲区内的DEM做线性插值替换让河道恢复到“河谷低于两岸”的状态。更省事的办法是用水文专业站点发布的修正DEM比如HydroSHEDS之类的水文修正产品再把河北省范围内截取出来。千万别拿原始30米DEM的河道高程去做洪水淹没分析那会得出离谱的结论。4.5 裁剪后的边界出现黑色边条NoData值没统一现象按省界裁剪后河北省边界内侧出现一圈规则的黑条刚开始只看到一片黑线放大后发现和边界线完全重合。原因原始分幅的NoData值不是同一个数。有人裁剪时把NoData设成255有人设成0镶嵌时又没有统一导致边界处有些像元被当成有效高程裁剪时这些“伪有效值”被保留成了黑色区域。解决镶嵌前统一NoData。最稳的做法是在镶嵌后用栅格计算器执行Con(IsNull(dem), -9999, dem)把所有无效像元固定写成-9999。裁剪时参数里的NoData也填-9999裁剪后检查边界一圈像元的统计值如果最小值是-9999说明边界处理干净了。5. 把河北省30米DEM用起来坡度、坡向、地形湿度指数与河网提取5.1 坡度分级太行山与冀南平原要分开设阈值河北省DEM做坡度分析直接用统一分级会出问题。太行山和燕山的坡度峰段集中在15度到30度冀南平原的“坡度”基本是噪声很多地方明明是一马平川算出来却有3到5度。dem rD:\hebei_dem\hebei_dem_30m_albers.tif slope_deg arcpy.sa.Slope(dem, DEGREE) # 分级0-2°平地2-6°缓坡6-15°坡地15-25°陡坡25°山地 remap 0 2 0;2 6 1;6 15 2;15 25 3;25 90 4 slope_cls arcpy.sa.Reclassify(slope_deg, VALUE, remap, NODATA) slope_cls.save(rD:\hebei_dem\result\slope_class.tif)分级阈值不是固定的。平原农业区建议把第一级下限提到2度小于2度按平地处理山区项目可以把陡坡起点提到30度更贴合冀北山地的实际坡度分布。Reclassify的输入必须是原始坡度值域不要拿已经做过平滑的坡度去分级那会把真实陡坡削平。5.2 坡向计算的“平地区翻车”与掩膜办法坡向在平缓区是全靠噪声驱动的结果这点做河北全域坡向时特别明显。华北平原地势平缓坡向却四面八方乱跳这不是算法错了是坡向定义本身依赖邻域高差高差趋近于零时任何微小的DEM噪声都会主导方向。slope_deg arcpy.sa.Slope(dem, DEGREE) aspect_deg arcpy.sa.Aspect(dem) # 坡度小于0.5度的区域坡向没有物理意义直接置为NoData aspect_clean arcpy.sa.SetNull(slope_deg 0.5, aspect_deg) aspect_clean.save(rD:\hebei_dem\result\aspect_clean.tif)阈值0.5度是经验值。如果项目区域是张家口坝上那种相对平坦的高原阈值可以提到1度如果是太行山峡谷区0.3度都够。这样处理后的坡向分类图山区和坝上交界处才不会有那种“碎玻璃”噪声。5.3 地形湿度指数TWI的计算链TWI是很多土壤湿度、滑坡和生态分析的前置变量。它的公式是ln(单位宽度汇流面积除以坡度正切值)在ArcGIS里通过填洼、流向、汇流累积三步算出。dem_fill arcpy.sa.Fill(dem, 30) flow_dir arcpy.sa.FlowDirection(dem_fill) flow_acc arcpy.sa.FlowAccumulation(flow_dir) slope_rad arcpy.sa.Slope(dem_fill, RADIAN) tan_slope arcpy.sa.Tan(slope_rad) # 避免坡度为0时除零 tan_slope_safe arcpy.sa.Con(tan_slope 0.001, 0.001, tan_slope) # 简化TWI汇流面积用累积像元数×30米像元大小近似 twi arcpy.sa.Ln((flow_acc * 30.0 1.0) / tan_slope_safe) twi.save(rD:\hebei_dem\result\twi.tif)这里解释两点。第一FlowAccumulation先跑在填洼后的DEM上如果不填洼平原区的伪洼地会让水流方向断裂汇流累积完全乱掉。第二公式里flow_acc * 30用30米像元的单边长度做简化近似严格的地表汇流面积还要考虑坡度对实际面积的拉伸但省级分析尺度下这个简化足够稳定。5.4 提取河网阈值选择与省界缓冲区陷阱用汇流累积量提取河网是最常见的河网生成方式。阈值直接影响河网密度河北山区和平原要分别对待。# 先用1000像元阈值提取试算约等于0.9平方公里集水面积 stream_raw arcpy.sa.Con(flow_acc 1000, 1) stream_raw.save(rD:\hebei_dem\result\stream_threshold1000.tif) # 平原地区阈值过大河网会断头降到300到500再试 stream_plain arcpy.sa.Con(flow_acc 400, 1) stream_plain.save(rD:\hebei_dem\result\stream_threshold400.tif)阈值1000在太行山区能得到比较完整的水系骨架但到冀南平原很多支流在平原区汇流面积积累不够会断在半路。平原地区阈值降到400再看。更准的做法是拿已经公开的1比25万河网做对比把阈值调到提取河网和现有河网重合率最高为止。还有一个省级分析特有的陷阱只拿河北省边界裁剪后的DEM提取河网省界附近河流的上游汇水面积会被切掉河网在边界处明显变稀这不是算法问题是数据范围问题。正确的做法是前期把邻省数据拼进来或者下载数据时就刻意多覆盖周边几十公里等河网提取完再按河北省界裁剪结果。6. 投入前的最后一步填洼控制、差值验证与省界缓冲区习惯6.1 填洼别一上来就用默认参数ArcGIS的Fill工具默认把所有洼地填到没有低点为止这在华北平原非常危险。平原的伪洼地会被成片填平汇流累积变成一摊辫状水道提取的河网完全不能用。我会在填洼时限制最大填充深度用z_limit参数把填充量压到5米或10米以内只填掉噪声级的小洼地保留真实的地形洼地。# 限制最大填洼深度避免平原被成片填平 dem_fill_limit arcpy.sa.Fill(dem, 10)填完后做一次“填洼差值图”把填洼后的DEM减去原DEM如果某个区域的差值和真实地势完全不匹配而且连成大片说明参数设置过激需要降低z_limit或改用其他数据源。6.2 用差值验证法给最终DEM上保险处理完成不等于能用。我会拿SRTM和最终DEM做一次逐点差值检查均值和中位数是否异常。diff arcpy.sa.Minus(dem_fill_limit, srtm_mosaic) mean_diff arcpy.management.GetRasterProperties(diff, MEAN).getOutput(0) std_diff arcpy.management.GetRasterProperties(diff, STD).getOutput(0) print(平均差值, mean_diff, 标准差, std_diff)河北省全域范围内两者差值的均值应该落在正负1米以内标准差在山区可以放宽到5米左右。如果出现某个方向系统性偏差比如西部山区整体高5米东部平原整体低3米那就是镶嵌时某个分幅没对齐要回去查那个分幅的投影和数据版本。6.3 省级分析先问一句边界外带了吗我处理这类省级DEM的过程里返工最多的就是省界缓冲区问题。有一回做河北西部某流域的河网提取河北与山西交界处的河流全部断头查了三天才意识到不是算法问题而是裁剪时把山西境内的上游切掉了。那次之后我养成了一个习惯任何省级地形分析第一步先把省界往外扩展一百公里保证海河和滦河的完整上游都在数据处理范围内等所有分析做完再切回河北省界出图。30米分辨率的河北省DEM本身精度够做项目但能不能真正落地取决于你留给自己的验证余量。希望这篇笔记能帮你把下载、拼接、裁剪和验证的弯路一次走完。本文还有配套的精品资源点击获取