简介青海省30米分辨率DEM数据包基于ASTER GDEM V3全球高程数据制作面向GIS从业者、地理科研人员及环境规划相关师生可用于地形分析、流域研究、灾害评估与生态制图等场景。压缩包共10个文件核心为GeoTIFF格式30米高程栅格及其tfw坐标参考、xml元数据另含青海省边界Shapefileshp、dbf、prj、sbn等方便在ArcGIS、QGIS中直接叠加裁剪或限定研究范围。包体约944.59MB已有494人浏览学习。数据采用WGS84坐标系全球通用、空间匹配可靠加载后可进行高程提取、坡度坡向计算、地形剖面与可视域分析结合边界矢量数据可快速研究青海高山、湖泊、草原等复杂地貌对植被分布、河流走向及气候变化的响应为区域地理研究、地质灾害评估与国土空间规划提供扎实的基础底图。1. 青海省30米DEM三件决定成败的准备工作“青海省DEM30米分辨率”乍看是个数据文件名实际做下来是一条完整的数据链下载、选源、镶嵌、投影、裁剪、补空洞、质检。我接过不少青海项目——水文分析、地质灾害排查、光伏选址最终都落在同一份30米DEM上而每回绕不开的问题几乎一样用哪个数据源上百幅分幅怎么拼省界裁剪后的白边怎么处理祁连山雪线附近的空洞拿什么补。这篇就把这套流程讲透。适用人群是那些要做省级尺度地形分析、又不想在数据预处理上耗掉整个工期的人。跟着复现两三个小时可以拿到一份能直接入库的青海省30米DEM。2. 数据源怎么选ASTER GDEM v3、SRTM与AW3D30的取舍与镶嵌2.1 市面上能拿到的省级30米DEM不止一个很多人以为青海省30米DEM有现成整包下载实际不是。省级产品通常是从全球或全国尺度数据里按省界裁剪出来的而能落到“30米”这个档位的公开免费数据源主要有三个。数据源官方分辨率覆盖范围质量特点适合场景ASTER GDEM v31弧秒约30米全球83°N–83°Sv3比v2少很多伪坑但雪线、裸岩区仍有局部空洞省级地形骨架、坡度坡向、插值底图SRTM v3SRTMGL31弧秒约30米北纬60°–南纬56°2000年采集空洞少但已插值填平细节偏老与ASTER互检、稳定性要求高的批量处理ALOS AW3D30约30米网格全球由5米DSM抽稀而来山区纹理好、空洞少但属于DSM河谷、沟谷提取以及作为空洞替换的替补源选型上我一般主用ASTER GDEM v3辅以AW3D30做补洞和交叉验证。原因很简单ASTER在青海的覆盖完整、下载渠道多官方就是按30米发布的做省级分析足够AW3D30的原始数据分辨率更高山区细节更可信但它是DSM树冠和房顶会让高程略微偏高。青海整体植被稀疏这个偏差影响很小但在湟水河谷的灌丛带和城区周边要留个心眼。2.2 下载前的参数确认分辨率、分幅与坐标系数据下载之前先把青海的范围框清楚。青海大致位于89.4°E–103.1°E、31.4°N–39.2°N按1度乘1度的标准分幅全境约占满14列乘8行的矩形框也就是110多幅图。单幅30米GeoTIFF文件大小在20MB上下整批下载总量约2–3GB。这个体量不要用浏览器一个一个点常见的做法是在地理空间数据云、USGS EarthExplorer这类平台上按范围框选后批量加入下载队列再用下载工具批量拉取。拿到压缩包后不要急着解压拼接先抽样检查三件事坐标系是不是WGS84经纬度EPSG:4326像元尺寸是不是0.0002778度左右位深是不是16位整型。这三项不一致的分幅混在一起后面merge出来的就是一张废图。尤其是部分镜像站会顺手把数据重投影到UTM同一批数据里混着两种坐标系拼出来会出现几十米的错位。2.3 用rasterio把上百幅tif一次拼起来分幅检查完用Python的rasterio做镶嵌是最高效的办法。下面这个脚本把指定目录下所有tif按文件名排序后合并适用于ASTER和AW3D30。import glob import rasterio from rasterio.merge import merge # 按文件名排序避免merge结果不稳定 tiff_list sorted(glob.glob(/data/qinghai_dem/tiles/*.tif)) # 打开所有分幅文件注意这里假设所有文件已是EPSG:4326且分辨率一致 src_files [rasterio.open(p) for p in tiff_list] # methodlast 表示重叠区用后打开的那幅覆盖前一幅 mosaic, out_transform merge(src_files, methodlast, nodata0) # 沿用第一幅的元数据更新尺寸和仿射变换 profile src_files[0].profile.copy() profile.update( heightmosaic.shape[0], widthmosaic.shape[1], transformout_transform, compressdeflate ) with rasterio.open(/data/qinghai_dem/qinghai_mosaic.tif, w, **profile) as dst: dst.write(mosaic, 1) # 释放文件句柄 for f in src_files: f.close()这段脚本的核心是merge(src_files, methodlast, nodata0)。method参数决定重叠区像素的取值策略last用后读入的文件first用先读入的min和max取重叠区的最小或最大高程值。我建议先用last拼一版接着统计空洞占比如果某个区域恰好是两幅文件的空洞重叠再改用min或max试试往往能救回一部分像元。nodata0是把0值统一识别为无效值因为ASTER分幅里常把背景区域填0或-9999拼之前最好先确认所有文件的实际NoData值。2.4 空洞修补fillnodata与跨源替换拼完的第一版数据马上要查空洞比例。青海的高原雪线、祁连山裸岩区以及可可西里的冻土带都是ASTER立体匹配容易失败的典型区域屏幕上表现为一块块黑色斑块。用一行Python就能统计整体空洞占比import rasterio with rasterio.open(/data/qinghai_dem/qinghai_mosaic.tif) as src: arr src.read(1) nodata src.nodata print(nodata占比: {:.2f}%.format((arr nodata).mean() * 100))空洞占比低于1%直接用rasterio的fillnodata做邻域内插即可超过5%建议用AW3D30做跨源替换因为大面积插值会把地形抹成平坦的“补丁”后续水文分析会翻车。跨源替换的常规做法是把两种数据重采样到同一网格然后以ASTER为主、AW3D30补洞gdal_calc.py -A qinghai_mosaic.tif -B aw3d30_resampled.tif \ --outfileqinghai_mosaic_filled.tif \ --calcwhere((AA_NoData), B, A) --NoDataValue-9999这条命令的--calc表达式意思是AASTER分幅拼接结果里凡是NoData的像元用B重采样后的AW3D30对应位置替换其余保留A。注意两个输入的坐标系和像元尺寸必须一致不一致时先用gdalwarp -tr 0.0002778 0.0002778 -r bilinear把B重采样到A的网格。替换完再跑一遍空洞统计目标是把NoData占比压到0.1%以下。3. 投影与裁剪如何把全球分幅变成青海省界内的可用DEM3.1 为什么不能抱着WGS84直接算坡度很多人做完镶嵌就直接用EPSG:4326计算坡度和坡长结果算出来的坡度值明显偏小、坡向也乱。原因不复杂在经纬度坐标系下X方向一个像元是约0.0002778度经度而青海纬度在31°N到39°N之间同样0.0002778度经度对应的地面距离只有约74米纬度方向却是约30.8米像元根本不是正方形。坡度是基于水平距离和垂直高差的比值算的水平距离都算错了结果自然全错。省级分析我一般用Albers等积投影青海全境落在中央经线96°E、双标准纬线32°N和37°N这一组参数内东西方向的形变控制得比较均衡。如果项目只涉及某个小区域比如海西州一个县直接按UTM分带处理更方便但青海东西方向跨了约14度经度UTM要切到45N、46N、47N三个带全省一张图时拼接边界的麻烦远大于收益。3.2 用缓冲裁剪一次解决白边裁剪到青海省界这件事最容易出的问题是“贴边白边”。直接拿省界矢量去裁DEM边界内侧往往留有一圈NoData因为这些像元在原始的1度分幅里就落在边缘重采样后没有被赋值。常见做法是先把省界向外缓冲一段距离再裁剪裁完再填一次洞从根本上避开白边。import rasterio from rasterio.mask import mask from rasterio.fill import fillnodata import geopandas as gpd # 读取省界确保shp与DEM使用同一地理坐标系 border gpd.read_file(/data/qinghai_dem/qinghai_boundary.shp) # 向外缓冲约1公里0.01度在青海约0.9公里够用 border_buffered border.buffer(0.01) with rasterio.open(/data/qinghai_dem/qinghai_mosaic_filled.tif) as src: out_image, out_transform mask( src, border_buffered.geometry, cropTrue, nodatasrc.nodata ) profile src.profile.copy() profile.update( heightout_image.shape[0], widthout_image.shape[1], transformout_transform, compressdeflate ) # max_search_distance单位是像元10像元约300米只修补边界附近的细小空洞 filled fillnodata(out_image, max_search_distance10) with rasterio.open(/data/qinghai_dem/qinghai_clip.tif, w, **profile) as dst: dst.write(filled, 1)这段里有三个参数值得说。border.buffer(0.01)的0.01是经纬度单位在青海纬度上约等于0.9公里只做保险作用不要设得太大否则把邻省地形也裁进来了。max_search_distance10的单位是像元而不是米它限制挖洞填补的最大搜索半径设太大会把空洞区填成一块过度平滑的“锅盖”设太小则补不干净。nodatasrc.nodata必须显式传递如果省略mask函数会用默认值容易把负值或0值误判成有效高程。3.3 重采样方式与水文填洼的准备裁剪完成的DEM还差最后一道预处理确认输出数据类型。坡度、坡向这类参数建议在浮点型DEM上计算整型数据在高差小的地区会出现大量相同坡度值的“台阶”影响后续分级统计。如果之前保存的是整型用gdal_translate加-ot Float32转一下即可。重采样到其他网格时凡是用于地形分析的都用bilinear或cubic不要用nearest——后者会把山脊线切出锯齿河谷提取时容易出现平行伪河道。另外一个容易忽略的点水文分析前要先填洼。青海内陆河流域多、盐湖周边地势平坦DEM里普遍存在伪凹陷直接提取河网会在盆地中央断头或绕圈。填洼是独立的处理步骤常见做法是交给Whitebox Tools或TauDEM这类专门工具处理不要把fillnodata当成填洼用它只是修补NoData不会处理真实地形里的凹陷。4. 避坑清单青海DEM处理中容易被坑的四个环节这一步说的都是我自己在青海项目里碰到过的真实问题每条按“现象→原因→解决”来写处理完这一轮后面的坡度、坡向、水文分析才能睡得着觉。4.1 “tiff转dem文件”不是格式魔术现象不少朋友搜“tiff转dem文件”以为DEM是一种需要特殊转换才能得到的专属格式拿到GeoTIFF之后到处找转换工具。原因这是把“DEM产品”和“文件格式”两件事混在一起了。DEM描述的是数据内容——每个像元存高程值而GeoTIFF、Esri Grid、SRTM的.hgt都只是承载它的文件容器。ASTER和SRTM分发的.tif本身就是一份完整的DEM文件不存在“不转就不能用”的问题。解决先确认你真正要的是什么。ArcGIS里常见的需求是把整型高程转成浮点型用于坡度计算或者把DEM导出成ASCII用于外部程序这些在Data Management Tools里都能完成本质是改数据类型或交换格式不是给DEM做“激活”。如果你只是想要一份能直接用的青海省30米DEM把第2、3章的产物保存成GeoTIFF就够了不需要再做任何格式转换。4.2 祁连山和可可西里上空的黑斑现象镶嵌后的DEM上祁连山雪线附近、可可西里冻土区出现一片片黑色斑块这些位置高程值为NoData坡度计算后形成夸张的深坑。原因ASTER GDEM v3虽然比v2改善明显但立体像对匹配在积雪、强反照率裸岩和纹理贫乏的平坦冻土上仍然会失败。青海正好把这些失败场景占全了。解决先按2.4的做法统计空洞占比和空间分布。如果空洞集中在少数几个区域用AW3D30局部替换如果分布零散直接用fillnodata内插。替换后一定要做目视检查重点看空洞边界是否出现“补丁棱”——也就是插值区域和原始地形之间有肉眼可见的折痕有的话用3x3低通滤波在边界上抹一下。4.3 省界裁剪后的贴边白边现象直接用青海省界shp裁剪得到的DEM在省界内侧有一圈白色NoData带宽度从几十米到几百米不等河网提取到这里全部断头。原因裁剪是严格的几何切割分幅镶嵌时边缘像元没有足够的外扩数据省界线经过的位置恰好落在NoData像元上。这个现象在ASTER分幅边缘尤其明显。解决第3章的缓冲裁剪是预防手段。如果已经裁出白边补救方法是把白边像元先转化为NoData再用fillnodata补一遍最后重新按省界精确裁剪。不要试图用ArcGIS的“边界平滑”功能去抹白边它处理的是矢量的几何形状对栅格NoData无效。4.4 高程里的负值和异常尖峰现象统计DEM最小值时发现大量负值比如-200、-1500位置集中在柴达木盆地的盐湖和水体周边最大值区域则出现明显高于周围地形的“尖塔”。原因ASTER在生产过程中水体、盐碱地和低反照率区域的立体匹配经常产生异常高程官方文档里也承认存在残余的伪值。负值不是真实海拔尖峰则是匹配错误导致的飞点。解决处理流程固定为“先统计再截断后内插”。先用gdalinfo -stats或Python看min/max把高程合理范围之外的像元设为NoData再对NoData做邻域内插。青海全省的高程范围大约在1500米到6860米之间低于1000米或高于7000米的值基本可以判定为异常。截断操作不要直接设为0那会让水面变成平地影响后续的水文分析。4.5 DSM当DEM用的偏差现象有人直接拿ALOS AW3D30当DEM做坡度分析和汇水提取发现河谷地带的坡度和预期明显不符山坡像元高差不自然。原因AW3D30的本质是DSM也就是地表模型它记录的是树冠、屋顶、植被表面的高度而不是裸地高程。青海虽然整体植被稀疏但河谷灌木带和城镇建成区的DSM比真实地形高出数米到十几米这在30米尺度上足以干扰坡度分级统计。解决如果你需要的是“DSM生成DEM”核心是滤波去掉地物高度不是做格式转换。常见做法是对DSM做形态学开运算或低通滤波把比地形更粗糙的地物细节磨平。青海这种植被稀疏区域可以先做差值对比把AW3D30与ASTER相减差值中系统偏高且聚集的区域就是DSM地物影响区有针对性地修正即可。5. 进阶验证坡度分级统计与多源DEM互检5.1 用gdaldem批量生成坡度并分级统计DEM入库后第一个常规用途就是坡度分级。用GDAL自带的DEMProcessing生成坡度再叠加像元面积做分级统计整个过程不需要打开ArcGIS。from osgeo import gdal import numpy as np import rasterio # 生成坡度栅格单位为度输入必须是投影坐标系下的DEM gdal.DEMProcessing( qhs_slope.tif, qinghai_clip.tif, slope, formatGTiff, slopeFormatdegree ) with rasterio.open(qhs_slope.tif) as src: slope src.read(1) transform src.transform # 投影后像元是正方形像元面积直接用仿射参数计算 pixel_area abs(transform.a * transform.e) # 常见坡度分级0-5、5-8、8-15、15-25、25-90 for lo, hi in [(0, 5), (5, 8), (8, 15), (15, 25), (25, 90)]: count np.sum((slope lo) (slope hi)) area_km2 count * pixel_area / 1e6 print(f{lo}-{hi}度: {area_km2:.0f} km², 占比 {count / slope.size * 100:.1f}%)需要注意DEMProcessing的输入必须是米制投影的DEM直接喂EPSG:4326会得到一堆看似合理实则错误的坡度值。分级标准可以根据项目调整但统计时要确认NoData没有参与任何一档面积累计。5.2 多源互差质检拿到成品DEM我习惯马上与另一份30米数据做互差这一步能快速暴露系统性偏移和局部异常。import numpy as np import rasterio with rasterio.open(qinghai_clip.tif) as a, \ rasterio.open(qinghai_aw3d30_clip.tif) as b: da a.read(1).astype(np.float32) db b.read(1).astype(np.float32) valid (da ! a.nodata) (db ! b.nodata) diff da[valid] - db[valid] print(f平均差值: {diff.mean():.2f} m) print(f标准差: {diff.std():.2f} m) print(f95%绝对差值: {np.percentile(np.abs(diff), 95):.2f} m)平均差值在±5米以内说明两个数据源整体一致标准差在山区达到15–20米是正常的因为两套数据的采集时间和匹配算法不同。如果某个区域差值超过30米用县界做分区统计定位重点检查是不是空洞修补区或者盐湖边缘的异常带。我现在的习惯是任何一份DEM不管来源多权威先跑一遍多源互差再决定要不要进入生产流程。数据源本身是个黑匣子用统计学尺子量一量比肉眼看来得可靠。希望帮到你。本文还有配套的精品资源点击获取