简介：青海省30米分辨率DEM数据包，基于ASTER GDEM V3全球高程数据制作，面向GIS从业者、地理科研人员及环境规划相关师生，可用于地形分析、流域研究、灾害评估与生态制图等场景。压缩包共10个文件，核心为GeoTIFF格式30米高程栅格及其tfw坐标参考、xml元数据，另含青海省边界Shapefile（shp、dbf、prj、sbn等），方便在ArcGIS、QGIS中直接叠加裁剪或限定研究范围。包体约944.59MB，已有494人浏览学习。数据采用WGS84坐标系，全球通用、空间匹配可靠；加载后可进行高程提取、坡度坡向计算、地形剖面与可视域分析，结合边界矢量数据可快速研究青海高山、湖泊、草原等复杂地貌对植被分布、河流走向及气候变化的响应，为区域地理研究、地质灾害评估与国土空间规划提供扎实的基础底图。
1. 青海省30米DEM：三件决定成败的准备工作
“青海省DEM（30米分辨率）”乍看是个数据文件名，实际做下来是一条完整的数据链：下载、选源、镶嵌、投影、裁剪、补空洞、质检。我接过不少青海项目——水文分析、地质灾害排查、光伏选址，最终都落在同一份30米DEM上，而每回绕不开的问题几乎一样：用哪个数据源，上百幅分幅怎么拼，省界裁剪后的白边怎么处理，祁连山雪线附近的空洞拿什么补。这篇就把这套流程讲透。适用人群是那些要做省级尺度地形分析、又不想在数据预处理上耗掉整个工期的人。跟着复现，两三个小时可以拿到一份能直接入库的青海省30米DEM。
2. 数据源怎么选：ASTER GDEM v3、SRTM与AW3D30的取舍与镶嵌
2.1 市面上能拿到的省级30米DEM不止一个
很多人以为青海省30米DEM有现成整包下载，实际不是。省级产品通常是从全球或全国尺度数据里按省界裁剪出来的，而能落到“30米”这个档位的公开免费数据源，主要有三个。
| 数据源 | 官方分辨率 | 覆盖范围 | 质量特点 | 适合场景 |
|---|---|---|---|---|
| ASTER GDEM v3 | 1弧秒（约30米） | 全球83°N–83°S | v3比v2少很多伪坑，但雪线、裸岩区仍有局部空洞 | 省级地形骨架、坡度坡向、插值底图 |
| SRTM v3（SRTMGL3） | 1弧秒（约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] # method='last' 表示重叠区用后打开的那幅覆盖前一幅 mosaic, out_transform = merge(src_files, method='last', nodata=0) # 沿用第一幅的元数据，更新尺寸和仿射变换 profile = src_files[0].profile.copy() profile.update( height=mosaic.shape[0], width=mosaic.shape[1], transform=out_transform, compress='deflate' ) 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, method='last', nodata=0)。method参数决定重叠区像素的取值策略：'last'用后读入的文件，'first'用先读入的，'min'和'max'取重叠区的最小或最大高程值。我建议先用'last'拼一版，接着统计空洞占比，如果某个区域恰好是两幅文件的空洞重叠，再改用'min'或'max'试试，往往能救回一部分像元。nodata=0是把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 \ --outfile=qinghai_mosaic_filled.tif \ --calc="where((A==A_NoData), B, A)" --NoDataValue=-9999这条命令的--calc表达式意思是：A（ASTER分幅拼接结果）里凡是NoData的像元，用B（重采样后的AW3D30）对应位置替换；其余保留A。注意两个输入的坐标系和像元尺寸必须一致，不一致时先用gdalwarp -tr 0.0002778 0.0002778 -r bilinear把B重采样到A的网格。替换完再跑一遍空洞统计，目标是把NoData占比压到0.1%以下。
3. 投影与裁剪：如何把全球分幅变成青海省界内的可用DEM
3.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, crop=True, nodata=src.nodata ) profile = src.profile.copy() profile.update( height=out_image.shape[0], width=out_image.shape[1], transform=out_transform, compress='deflate' ) # max_search_distance单位是像元，10像元约300米，只修补边界附近的细小空洞 filled = fillnodata(out_image, max_search_distance=10) 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_distance=10的单位是像元而不是米，它限制挖洞填补的最大搜索半径，设太大会把空洞区填成一块过度平滑的“锅盖”，设太小则补不干净。nodata=src.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', format='GTiff', slopeFormat='degree' ) 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(f'95%绝对差值: {np.percentile(np.abs(diff), 95):.2f} m')平均差值在±5米以内说明两个数据源整体一致；标准差在山区达到15–20米是正常的，因为两套数据的采集时间和匹配算法不同。如果某个区域差值超过30米，用县界做分区统计定位，重点检查是不是空洞修补区或者盐湖边缘的异常带。我现在的习惯是：任何一份DEM，不管来源多权威，先跑一遍多源互差再决定要不要进入生产流程。数据源本身是个黑匣子，用统计学尺子量一量，比肉眼看来得可靠。希望帮到你。
本文还有配套的精品资源，点击获取
