简介中国范围数字高程模型栅格数据是一份面向地理信息分析、课程设计与科研应用的全国尺度高程数据。原始30米分辨率的SRTM高程数据经ArcGIS镶嵌拼接与重采样后得到1公里分辨率的完整全国栅格采用WGS-84全球坐标系以标准TIFF栅格格式存储可直接在ArcGIS等GIS软件中加载使用。压缩包共7个文件以TIFF主栅格为核心另含tfw坐标配准信息、ovr金字塔加速显示、xml元数据等辅助文件整体仅26.51MB便于下载与分发。目前已有695人学习或下载适合需要快速获取全国尺度高程数据开展地形分析、栅格计算、空间分析等工作的GIS使用者。数据经过质量检查完整正确包含完整的高程信息覆盖全国范围省去自行拼接与重采样的繁琐步骤可让读者将精力集中在后续分析与应用上。1. DEM栅格数据不是图片中国范围落地的第一道坎拿到一份标注“中国范围”的DEM栅格数据很多人第一反应是把它拖进GIS当图片看结果要么是黑乎乎一片要么是花屏色带折腾半天才发现问题出在不理解栅格的组织形式——DEM的每个像元存的是高程数值不是颜色值。真正让新手翻车的往往是分辨率不一致、投影坐标系混乱、NoData被当成0参与计算这三件事处理不好后面拼接、裁剪、算坡度全跟着错。这篇文章想解决的就是“DEM栅格数据怎么在中国范围内落地成能用的高程底图”从数据源选型、拼接裁剪参数到黑边接缝排错一条线走完你照着做至少能避开八成常见坑。2. 读懂DEM栅格的组织形式分辨率、像元类型与NoData的坑2.1 栅格数据的三种组织形式BSQ、BIL与BIP栅格数据落盘时的组织形式常见有三种BSQ波段顺序、BIL波段按行交错、BIP像元按点交错。DEM一般是单波段单文件BSQ最简单也最容易读绝大多数情况你拿到的.tif或.img都是BSQ组织但遇到.img的旧格式或者ENVI标准格式时BIL和BIP也有出现。三者的差别在于数据在磁盘上的排列顺序。BIP把同一个像元的高程值连续排在一起做逐像元计算时缓存命中率最高BIL在波段数少的场景下读写均衡BSQ则是最直观——文件里先存完整个波段的所有像元再存下一个波段。写处理脚本前先看一眼数据格式类型否则用numpy直接读二进制时会读出一堆乱码一样的数。# Python里快速判断Img格式的组织方式不需要打开桌面软件 import struct def peek_img_header(img_path): with open(img_path, rb) as f: header f.read(256) # IMAGINE .img文件头部通常会在第80字节附近标识数据类型 # 0BSQ, 1BIL, 2BIPENVI标准下不同但IMAGINE适用 print([hex(b) for b in header[80:96]])这段代码只做头部探测。实际工作中我更建议直接用GDAL读取因为它已经帮你处理了格式差异gdalinfo可以打印出完整的组织形式信息不需要自己去解析二进制。记住一个原则优先用GDAL统一读取只有做底层性能优化时才需要关心BSQ、BIL、BIP的物理排布否则不值得在这上面花时间。2.2 像元类型与NoData高程为什么要用浮点存DEM的高程值在数据文件里通常用整型或浮点型保存。SRTM原始数据是16位有符号整型单位是米但ASTER GDEM和ALOS AW3D30则常用浮点型因为浮点能表达亚米级的小数高程差。这里埋着一个大坑整型DEM做坡度计算时NoData常常被记为-9999或-32768如果你把整个栅格转成浮点后忘记处理NoData最小值和最大值的统计会被这些哨兵值污染拉伸显示时整幅图变成灰蒙蒙一片。我的做法是拿到数据先跑一遍GDAL统计确认NoData值到底是什么再统一设置成相同的NoData标记。尤其在中国范围内分幅下载的数据不同分幅的NoData可能不完全一致——有的写-9999有的写-32768拼接前必须统一。# 用gdalinfo查看NoData定义这一步跑完再开始拼接 gdalinfo SRTM_56_04.tif | grep NoData如果输出里没有NoData行说明这副数据没有设置NoData需要自己补上。常见的脏数据表现为边缘有一圈-9999或0的像元后面做坡度或水文分析时会出现离谱的负值或零值洼地。设置的方法很简单在ArcGIS里用Copy Raster勾选NoData值或GDAL命令行加一句-a_nodata -9999。在选像元类型时能用整数就用整数存储但出分析结果坡度、坡向、累积流量时务必转回浮点否则截断误差会让你在验证时对不上数。这一点没有玄学纯粹是数值精度问题。2.3 SRTM、ASTER与ALOS混拼前先检查高程基准中国范围的DEM主要涉及三套公开数据SRTM覆盖全球分辨率30米或90米、ASTER GDEM V330米、ALOS AW3D3030米标称实际部分地区12.5米效果。它们的椭球和高程基准并不完全一致SRTM使用EGM96大地水准面改正ASTER GDEM V3也做了类似改正但二者在山区地形上存在数米到十几米的差异。拿到不同源数据拼接前最稳妥的做法是在重叠区域抽样对比一下高程差差异大就老老实实只用一个数据源。理论上这些数据都已基于WGS84椭球但“基于WGS84”不代表高程基准一致。EGM96与CGCS2000的高程异常在中国境内通常有几米偏差这不会影响做相对分析坡度、坡向、流域但如果你要把DEM用于洪水淹没模拟或工程勘察必须用一个已知水准点把高程系统校准到1985国家高程基准。# 用gdalwarp统一重投影到CGCS2000 / 高斯-克吕格投影并重采样到30米 gdalwarp -t_srs EPSG:4534 -tr 30 30 -r bilinear \ -of GTiff SRTM_56_04.tif SRTM_56_04_cgcs2000.tif参数说明-t_srs EPSG:4534是CGCS2000 / 3-degree Gauss-Kruger zone 26的坐标系适合中国中东部区域-tr 30 30强制输出为30米像元-r bilinear做双线性重采样适合高程这类连续表面。如果你的研究区跨两个投影带通常改用Albers等积投影而不是高斯投影后面第4章会说。3. 下载中国范围DEM自由30米与12.5米的三个主力源3.1 ASTER GDEM V3地理空间数据云的30米亚洲区域数据对于中国范围做分析我最常用的是从地理空间数据云下载ASTER GDEM V3它覆盖了北纬83度到南纬83度的绝大部分陆地区域中国全境都能拿到。这个源的优点是在国内下载速度快、不需要额外折腾整幅数据按1度分幅每幅约3601行乘3601列30米分辨率文件格式是GeoTIFF。下载前先规划好目标分幅。中国经度跨度约73度到135度纬度约18度到54度按1度分幅大概要六七十幅一次性全下载不太现实。我的习惯是先画个面状边界然后在平台里按边界筛选瓦片或者直接用行列号规律批量记下需要的分幅号。常见的一个误解是“分辨率是30米就一定能看清山脊”但ASTER GDEM在云覆盖多的区域可能有空洞和条带噪声这些坑换数据源解决不了。对比之下ALOS AW3D30在垂直精度上通常优于ASTER但下载渠道需要从日本的JAXA官网走大文件下载速度不稳定适合做小范围的精细化分析。3.2 ALOS AW3D3012.5米细节与条带噪声的取舍ALOS AW3D30是JAXA发布的全球30米网格DEM但它的原始数据来自ALOS卫星的立体像对实际地形细节明显好于ASTER在一些区域甚至可以提取出类似12.5米的地形纹理。做山体阴影可视化或提取等高线时AW3D30的地形细节优势非常明显。不过AW3D30有一个比较烦的特性数据存在南北向的条带噪声。在低纬度和地形平坦地区尤其明显表现为山体阴影图上一道一道近乎平行的明暗条纹。追究起来是轨道方向相关的处理误差残留。要消除这种条纹常见做法是做一个低通滤波或者在坡度计算后用中值滤波抹掉条带方向的高频。# GDAL 3.x用gdal_fillnodata填补数据空洞再用gdal_translate转成tif gdal_fillnodata -md 30 -si 5 AW3D30_tile.tif AW3D30_filled.tif gdal_translate -ot Float32 -co TILEDYES -co COMPRESSDEFLATE \ AW3D30_filled.tif AW3D30_ready.tif这段命令里的-md 30意思是搜索半径30个像元即把直径30像元范围内的空洞用周围有效值插值填补-si 5是平滑迭代5次。值得注意的是gdal_fillnodata只对栅格空洞有效对条带噪声需要另一套办法。条带消除我更倾向于在ArcGIS里用Focal Statistics做3×3或5×5的均值滤波滤波半径大了会抹平真实地形半径小了条带压不掉参数要盯着山体阴影图反复试。3.3 GLO-30与SRTM全球数据的中国分幅特点Copernicus GLO-30是欧洲空间局发布的30米全球DEM在中国范围内同样覆盖完整垂直精度普遍评价不错。相对ASTERGLO-30的洞和条带少很多但它的数据分幅按经纬度不规则切下载方式需要按AWS S3 bucket的目录结构来找对不熟悉命令行的用户来说门槛偏高。SRTM虽然名声大但它的原始30米数据只在北纬60度到南纬56度之间中国范围内西南边境地区覆盖有缺失。国内经常能下载到的SRTM 90米数据USGS版本更适合做大尺度地形分析比如全国尺度的坡度分区不太适合做县级尺度的高精度工程分析。90米像元在地形起伏大的区域会明显抹掉山谷细节做汇水面积分析时河道位置会偏移几十米到上百米。我在做全国尺度项目时通常用GLO-30或ASTER GDEM V3统一拼接输出成Albers等积投影的30米镶嵌结果。做省级或流域精细化分析时才换ALOS因为它对微小地形的描述更好。3.4 下载前的瓦片清单规划按经纬度网格批量落位无论从哪个平台下载第一步肯定是确定需要哪些瓦片。中国范围跨30多个经度、30多个纬度如果你想省事直接下载全中国所有分幅50多GB的原始文件不是问题问题是你拼接时内存会爆掉。更务实的做法是按目标区域框选瓦片只下载覆盖范围内的分幅。# 用Python根据经纬度边界生成ASTER GDEM的瓦片行列号清单 def dem_tile_list(lon_min, lon_max, lat_min, lat_max): tiles [] for lat in range(int(lat_min), int(lat_max) 1): for lon in range(int(lon_min), int(lon_max) 1): # ASTER GDEM V3瓦片命名ASTGTMV3_XX_YYDEM.tif tiles.append(fASTGTMV3_{lat:02d}_{lon:03d}DEM.tif) return tiles tiles dem_tile_list(110, 115, 30, 35) # 山东省中部某区域示例 print(len(tiles), tiles[:5])这段代码逻辑很简单按整经纬度度生成瓦片文件名。ASTER瓦片以纬度带为行、经度为列命名负值区域要做偏移处理但中国区域都在北半球东半球行列号相对好算。下载前先按这个清单核对平台上的文件是否存在能少走弯路。实际下载时把清单列表存成文本用浏览器的多线程下载工具或平台的批量下载功能比一页页手动点快得多。4. 把分幅DEM拼成中国范围拼接裁剪与坐标系的落地操作4.1 Mosaic to New Raster拼接中国全境的参数组合拿到若干分幅DEM后最直接的做法是在ArcGIS里用Mosaic to New Raster工具。第一步把格式、像元类型、波段数、NoData值统一定下来第二步设置Mosaic Method和Blend Width。# 使用arcpy执行Mosaic to New RasterPython窗口内运行 import arcpy arcpy.env.workspace rD:\dem_tiles arcpy.env.outputCoordinateSystem arcpy.SpatialReference(CGCS2000 Albers) arcpy.MosaicToNewRaster_management( input_rastersASTGTMV3_30_110.tif;ASTGTMV3_31_110.tif;ASTGTMV3_30_111.tif, output_locationrD:\dem_mosaic, raster_dataset_name_with_extensionchina_dem_30m.tif, coordinate_system_for_the_raster#, cellsize30, pixel_type16_BIT_SIGNED, number_of_bands1, mosaic_methodBLEND, mosaic_colormap_modeFIRST )参数说明pixel_type选择16_BIT_SIGNED能保留-9999这类NoData值但如果你的源数据是浮点型这里改成32_BIT_FLOAT更稳妥BLEND会在瓦片重叠区做平滑过渡避免生硬的接缝线cellsize必须统一为30否则不同分辨率瓦片混拼后输出网格会错位。拼接完成后立刻检查栅格统计值最小值和最大值如果出现-9999或0说明NoData被混了进来需要在后续处理步骤前重设。4.2 用中国范围矢量边界裁剪掩膜外的NoData才是黑边元凶拼接完成后下一步通常是按中国国界或者省界做裁剪。常见的错误是只用矢量边界做“外矩形裁剪”然后手工去抠边界外的区域结果裁剪结果外圈一圈白边或黑边怎么调都难看。这是NoData和背景0值混淆导致的问题。正确做法是使用带掩膜的Extract by Mask并且确认矢量边界和栅格投影一致。如果你的矢量边界是CGCS2000地理坐标栅格已经转成了Albers投影务必先把矢量做投影转换再裁否则边界会偏移几十米。在ArcGIS里用Project工具对矢量执行投影投影参数与栅格的Albers参数保持一致就避免了错位。# 用gdalwarp做带掩膜的裁剪一步完成投影和裁剪 gdalwarp -cutline china_province.shp -crop_to_cutline \ -t_srs EPSG:4529 -tr 30 30 -r bilinear \ china_dem_30m.tif shandong_dem_30m.tif注意这里EPSG:4529是CGCS2000 / 3-degree Gauss-Kruger CM 117E适合山东省这种位于带内的区域。如果你的目标区域跨度超过3度建议改用EPSG:102025中国双标准纬线Albers这类等积投影避免边缘拉伸变形影响面积计算。裁剪完成后要检查黑色边缘把NoData再次统一设置确保边缘外是透明而不是0。4.3 投影转换到CGCS2000 Albers面积计算不变形在全国尺度拼接中国范围DEM时Albers等积投影是最稳妥的选择。它的好处是面积不变形坡度、坡向在中等纬度区域变形也可接受。国内很多成果规范要求使用Albers或高斯投影前者用于区域综合分析后者用于大比例尺工程制图。# 从WGS84地理坐标转换到CGCS2000 Albers等积投影 gdalwarp -overwrite -s_srs EPSG:4326 -t_srs EPSG:102025 \ -tr 30 30 -r cubic -of GTiff \ china_dem_30m_wgs84.tif china_dem_30m_albers.tif这里的-s_srs EPSG:4326指定源数据坐标-t_srs EPSG:102025指定目标投影-r cubic使用三次卷积重采样对高程这种连续表面来说比双线性更平滑。重采样方法会影响高程值最近邻会保留原始值但产生锯齿边缘双线性和三次卷积会平滑地形但也可能让峰谷值略微钝化。做水文分析推荐用双线性做可视化用三次卷积效果好没有绝对最优。4.4 提取像元到表格把DEM栅格转成Excel可读的高程表许多从业者需要把栅格高程导出到Excel配合采样点做统计。ArcGIS里可以用Raster to Point把栅格转成点要素再用Table to Excel导出但数据量大时这个流程很慢更高效的方式是直接生成ASCII再读入表格。# 先转成ASCII网格再用Python直接写CSV from osgeo import gdal import pandas as pd ds gdal.Open(shandong_dem_30m_albers.tif) band ds.GetRasterBand(1) array band.ReadAsArray() rows, cols array.shape lon, dx, _, lat, _, dy ds.GetGeoTransform() xs [lon i * dx for i in range(cols)] ys [lat j * dy for j in range(rows)] grid pd.DataFrame(array) grid.to_csv(dem_export.csv, index_labelrow)这段代码把栅格从GDAL读入numpy数组然后借助GeoTransform计算每个像元的经纬度坐标最后输出CSVExcel直接打开就是一张二维高程表。要注意的是band.ReadAsArray()在大范围数据上容易占满内存更稳妥的做法是分块读取一次只读几十行循环写入CSV。采样点提取则建议直接用ArcGIS的Extract Multi Values to Points输出带高程字段的属性表再导出dbf或Excel速度比Raster to Point快一个量级。5. DEM拼接与裁剪的5条踩坑记录黑边、接缝与坐标偏移5.1 拼接后两幅数据之间出现明显接缝现象用Mosaic to New Raster拼完两幅相邻DEM重叠区域出现一条横向或纵向的高程突变带看起来像台阶但两幅图单独看都没问题。原因常见原因是两幅DEM高程基准或NoData值不一致比如一个用ASTER源、一个用SRTM源二者在山区重叠区高程相差数米也可能是相邻瓦片的分辨率不一致重采样到同一网格时产生了系统偏移。解决把项目统一为单一数据源不要混源拼接。如果必须混源先在重叠区计算两份数据的平均差用一个常数修正后重拼。在Python里可以用差值统计后叠加修正值但这属于后期校正能不做尽量不做。5.2 裁剪后边界外是纯黑而不是透明现象裁剪后的DEM外框是黑色或Z值异常的区域在ArcMap里用拉伸显示时一片漆黑做坡度分析时边界一圈出现离谱的大角度值。原因这是因为矢量边界外的栅格被赋予了一个非NoData的值常见为0而0米在海洋区域看起来合理但在山区就被拉低整体色带导致陆地部分一片黑。严格说是NoData设置被Extract by Mask覆盖掉了。解决在裁剪工具参数里指定NoData值为-9999或者裁剪后用Copy Raster重设NoData再检查统计值。裁剪完的栅格用gdalinfo -stats看一眼最小最大值发现最小值是0且不该有0的区域就要警惕。5.3 坡度与坡向图出现横条纹或网格状纹理现象输出坡度图时平缓区域出现规律性的横条纹、竖条纹或网格状纹理看起来像印刷网点比例尺拉大后纹理更明显。原因原始DEM存在系统条带噪声或地板量化误差ASTER GDEM和ALOS常见的条带噪声、SRTM的网格状噪声在坡度计算时被一阶差分放大。解决对DEM做3×3均值滤波后再算坡度。如果你担心滤波损失地形细节可以用Focal Statistics的MEDIAN类型替代均值保留边缘的同时压制冲激噪声。滤波半径参数建议在3×3和7×7之间试先用山体阴影图目测再定量比较滤波前后的坡度均值变化。5.4 裁剪结果和高分影像明显错位半个像元现象把DEM生成的等高线叠加到高分影像上等高线与山脊线整体偏移距离接近15米或30米的半像元左右。原因裁剪时矢量边界和栅格没有对齐可能是投影转换方法用了最小公分母或者栅格本身的角点坐标在重采样时被取整到相邻网格。解决先用gdalinfo查看裁剪前后栅格的原点坐标确认是否为30米的整数倍。如果不是用gdal_translate加-a_ullr参数手动校正角点坐标。更常见的做法是所有中间步骤都保留原始像元对齐只在最后一步输出时做重采样减少多重采样的叠加误差。5.5 大批量拼接时内存爆掉或程序闪退现象在中国全境的拼接过程中ArcGIS直接无响应或者Python进程在读取第5幅瓦片时内存溢出退出。原因大范围高分辨率栅格数据占用的内存很容易超过16GB。中国全境30米DEM的阵列尺寸大约是22000行乘20000列按16位整型算约880MB看起来不大但GDAL的缓存机制、显示刷新、金字塔构建叠加之后内存占用会飙到好几GB各种临时文件和程序一起吃满资源。解决断掉ArcGIS的自动金字塔构建使用GDAL分块拼接函数gdal.BuildVRT先建虚拟栅格再转成单一GeoTIFF转的时候设置-co TILEDYES并按512×512块写入。另外处理时限制GDAL_CACHEMAX为512MB避免缓存无限膨胀。# 使用BuildVRT合并瓦片避免一次性ReadAsArray gdalbuildvrt china_dem.vrt ASTGTMV3_*.tif gdal_translate -co TILEDYES -co BIGTIFFYES \ -co COMPRESSDEFLATE china_dem.vrt china_dem_30m.tif参数说明BIGTIFFYES让输出超过4GB时自动升级为BigTIFF格式避免文件大小上限报错COMPRESSDEFLATE压缩率高适合高程栅格但读取时会消耗一点CPU作为存储格式比较理想。这种方法在把中国全境30米DEM拼成单一文件时非常利索内存占用稳定在2GB以内。6. 用已知高程点验证DEM质量顺势解决DSM转DEM的粗差做完拼接裁剪和投影转换最不该省的一步是用已知高程点验证。我在山东省做过一次30米DEM质检把无人机LiDAR点云抽稀成地面控制点和ASTER GDEM V3对比中误差在平原有3到5米山地达到8到12米个别植被茂密的山沟出现20米以上的异常差。如果你没有LiDAR也可以用国家测绘地理信息发布的水准点成果或者Google Earth里筛选的高精度地标点原则是选取地形平缓、没有建筑物遮挡的位置。具体操作上我在ArcGIS里用Extract Multi Values to Points把DEM高程提取到测量点上算差值的均方根误差和中误差再看有没有系统偏差。差值均值如果始终为正或始终为负说明DEM存在整体高程偏置可能是高程基准不一致修正方法是整体加减一个常数。差值标准差过大则说明局部地形失真需要检查是否有条带噪声或空洞填补过度。剖面线验证是我惯用的第二步沿山脊线画一条Profile看高程曲线是否平滑、有无锯齿跳动。锯齿多的地方对应原始瓦片接缝或填补空洞区域直接在剖面图上就能定位。顺手还能发现DSM混入DEM的问题——如果剖面线穿过树林边缘原本平滑的地形却出现一个突兀的小包子很可能原始数据是DSM而非DEM。说到DSM转DEM很多热词检索里都在问“DSM生成DEM”。这是一个滤波问题DSM包含地表建筑和植被的高程DEM只保留裸地面。常见做法是对DSM做形态学开运算或渐进式数学形态学滤波把比周围明显突出的像元削平。在OpenCV里做底帽变换能提取出植被和建筑的“局部突起”然后从DSM里减去这部分就是近似的地面高程。# 用OpenCV形态学重建从DSM中提取地物并生成近似DEM import cv2 import numpy as np from osgeo import gdal ds gdal.Open(dsm_30m.tif) dem ds.GetRasterBand(1).ReadAsArray().astype(np.float32) dem[dem -9999] np.nan kernel cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (5, 5)) opening cv2.morphologyEx(dem, cv2.MORPH_OPEN, kernel) ground np.where(np.isnan(dem), np.nan, opening) ds_out gdal.GetDriverByName(GTiff).Create(dem_est.tif, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Float32) ds_out.GetRasterBand(1).WriteArray(ground)这段代码做了3×3的椭圆核开运算把比核尺寸小的凸起物去掉用于估计地面。MORPH_OPEN先腐蚀后膨胀形态学上可以消除小于结构元素的亮色细节——对应树冠、房屋等局部凸起。核越大滤掉的建筑物尺度越大但也会抹平真实山脊。这个方法的缺点是遇到陡峭地形时会把山谷填起来和真实DEM出现系统性偏差。如果你只追求精度别自己滤波直接下载成品DEM更靠谱只有拿不到理想数据源时才值得走这一步。最后的习惯是每次交付DEM成果前我会在ArcMap里叠加等高线和高分影像做一遍目视检查重点看河流谷地是否连贯、山脊线是否圆滑、城市周边有没有异常坑洼。这一步玄学成分确实有点大但往往能发现数值检验发现不了的问题。处理DEM这条路的坑很密集但只要数据源、投影参数、NoData这三个关口守住后面的大多数分析都不会翻车。希望帮到你。本文还有配套的精品资源点击获取