全球90米地形数据实战从零获取MERIT DEM到专业级地图制作最近在做一个全球尺度的水文模拟项目找地形数据时真是费了不少功夫。SRTM 90米数据虽然经典但海洋区域的数据空洞和山地地区的条带噪声一直是个头疼的问题。后来发现了东京大学发布的MERIT DEM同样是90米分辨率但在数据质量上做了显著的改进特别是通过融合多种数据源有效减少了误差。不过这个数据的下载、处理流程和官方文档里写得不太一样尤其是对于不熟悉地理信息处理流程的朋友可能会在数据转换、坐标系统这些环节卡住。今天我就把自己从注册账号、下载数据到完成预处理、符号化渲染的完整操作链梳理出来希望能帮你绕过那些我踩过的坑。这篇文章面向的是需要使用全球地形数据进行科研分析、环境建模或地图制作的GIS初学者和研究人员。我们不会停留在单纯介绍数据特性而是聚焦于实战操作确保你拿到手的是一套可以直接用于分析的、干净的地形数据。整个过程会涉及一些命令行操作和GIS软件的基本使用但别担心每一步我都会解释清楚背后的原理。1. 数据获取解密MERIT DEM的存储结构与高效下载策略MERIT DEM的全称是“Multi-Error-Removed Improved-Terrain DEM”由东京大学的研究团队开发。它的核心价值在于通过智能算法融合了SRTM、AW3D等多套数据并移除了多种系统性误差最终生成了这套目前全球范围内质量最高的90米分辨率公开地形数据集。数据版本目前稳定在v1.0.3覆盖范围是南纬60度到北纬90度南极洲不在其内。数据以5度×5度的瓦片形式组织每个瓦片是一个6000×6000像素的GeoTIFF文件。这个设计非常巧妙既保证了单文件不会过大又方便用户按需下载特定区域。文件命名规则是理解下载的关键。例如文件n30w120_dem.tif代表了一个地理范围其左下角像素中心位于北纬30度、西经120度因此这个文件实际覆盖了N30° 到 N35°W120° 到 W115°的区域。这里有一个细节需要注意由于像素中心定位实际边界会有约0.0004度的微小偏移在大多数宏观分析中可忽略但在需要极高精度的边界拼接时需留意。提示在规划下载区域时建议先用GIS软件如QGIS绘制一个研究区的矢量边界框计算出其经纬度的最大最小值再根据5度的整数倍去反推需要下载哪些瓦片这样可以避免遗漏或下载过多无用数据。数据提供方为了便于分发将原始的5度瓦片进一步打包成了更大的30度×30度压缩包。例如dem_tif_n30w120.tar这个压缩包就包含了从N30到N60、W120到W090这个巨大区域内的所有5度瓦片。对于需要大范围数据如整个大陆的用户直接下载这种大包效率更高。下载实战步骤访问与注册首先打开数据官网。页面设计比较简洁找到注册入口通常只需提供邮箱和设置密码即可完成注册。注册后可能需要等待一封激活邮件有时在垃圾箱激活后即可登录。定位数据文件登录后页面会列出不同版本和数据格式的目录。我们找到MERIT_DEM/v1.0.3/目录。里面通常有按不同压缩格式如tif格式、mrf格式分类的文件夹。对于大多数用户选择dem_tif/这个目录里面就是按30度区块打包的.tar文件。选择与下载确定你的研究区域所对应的30度区块编号。例如你需要中国东部地区的数据可能涉及n20e100,n30e100等区块。直接点击对应的.tar文件链接即可开始下载。由于文件较大每个tar包约1-2GB建议使用具有断点续传功能的下载工具如aria2c,wget -c并确保网络稳定。# 使用wget命令进行下载示例需替换为实际文件名和cookie信息 # 注意由于需要登录直接wget链接可能无效。更可靠的方法是使用浏览器插件获取下载链接和Cookie后用以下格式 wget --headerCookie: 你的登录Cookie信息 -c 完整的.tar文件下载链接 -O dem_tif_n30w120.tar解压与整理下载完成后使用解压工具如7-Zip、Bandizip或命令行tar解压tar包。解压后会得到一堆5度瓦片的.tif文件。建议立即建立清晰的文件夹结构进行管理例如按大洲或项目名称分类存放。区块打包文件名示例覆盖经纬度范围包含的5度瓦片数量预估解压后大小dem_tif_n30w120.tarN30-N60, W120-W09036个约 4 GBdem_tif_s10e140.tarS10-N20, E140-E17036个约 4 GBdem_tif_n00w060.tarN00-N30, W060-W03036个约 4 GB2. 核心预处理高程基准转换与数据镶嵌下载得到的原始DEM数据其高程基准是EGM96大地水准面也就是我们常说的“海拔高”。然而在许多空间分析、特别是涉及GPS观测值基于WGS84椭球或进行精确三维建模的场景中我们需要的是基于WGS84椭球的“椭球高”。这两者之间的差值就是“大地水准面起伏”在某些地区可达数十米甚至百米不进行转换会引入显著误差。为什么必须转换想象一下你要计算两个点之间的真实坡度或者将地形数据与卫星激光测高数据如ICESat-2进行融合分析。如果你的地形是海拔高而其他数据是椭球高这就好比用米尺和市尺去测量同一个物体然后直接比较结果必然失真。因此将EGM96正高转换为椭球高是让MERIT DEM数据能够与其他现代地理空间数据集“无缝对话”的关键一步。转换实战使用GDAL进行高效批量处理GDALGeospatial Data Abstraction Library是处理栅格数据的瑞士军刀。我们将使用其命令行工具gdalwarp来完成基准转换。这里需要一个关键文件EGM96大地水准面起伏格网文件通常为.gtx格式。你可以从NASA等机构下载或者使用一些GIS软件自带的如QGIS安装目录下可能包含。假设你已经将EGM96的格网文件egm96_15.gtx放在工作目录下以下是一个批处理脚本示例用于处理单个目录下的所有.tif文件#!/bin/bash # 批量EGM96转椭球高脚本 # 将脚本保存为 convert_height.sh并在终端中运行bash convert_height.sh INPUT_DIR./raw_dem_tiles # 存放原始瓦片的文件夹 OUTPUT_DIR./converted_dem # 输出文件夹 EGM96_GRID./egm96_15.gtx # EGM96大地水准面模型文件路径 mkdir -p $OUTPUT_DIR for tif_file in $INPUT_DIR/*_dem.tif; do if [ -f $tif_file ]; then filename$(basename $tif_file) output_file$OUTPUT_DIR/${filename%.tif}_ellipsoidal.tif echo 正在处理: $filename - ${filename%.tif}_ellipsoidal.tif # 使用gdalwarp进行转换 # -s_srs: 源坐标系WGS84 EGM96高程 # -t_srs: 目标坐标系WGS84椭球 # -srcnodata/-dstnodata: 处理无数据值 gdalwarp -s_srs projlonglat datumWGS84 no_defs geoidgrids$EGM96_GRID \ -t_srs projlonglat datumWGS84 no_defs \ -srcnodata -9999 -dstnodata -9999 \ -r bilinear -of GTiff \ $tif_file $output_file echo 完成: $output_file fi done echo 所有文件转换完成数据镶嵌Mosaicking转换完成后你得到的是一个个独立的5度瓦片。对于区域分析我们需要将它们拼接成一张完整的地图。这里推荐使用gdal_merge.py或gdalwarp进行镶嵌并注意处理瓦片边缘可能存在的微小值差异。# 使用gdal_merge.py进行简单镶嵌 # 该命令会将所有匹配的tif文件合并成一个大的mosaic.tif gdal_merge.py -o ./study_area_mosaic.tif -of GTiff -n -9999 -a_nodata -9999 ./converted_dem/*_ellipsoidal.tif # 更推荐使用gdalwarp功能更强能进行重投影和更好的重采样 gdalwarp -of GTiff -t_srs EPSG:4326 \ -srcnodata -9999 -dstnodata -9999 \ -r bilinear \ ./converted_dem/*_ellipsoidal.tif \ ./study_area_mosaic_final.tif注意镶嵌时务必设置正确的-srcnodata和-dstnodata参数MERIT DEM的无数据值通常是-9999否则无数据区域可能会被插值成异常值在后续分析中产生“悬崖”效应。3. 在QGIS中的高级可视化与符号化技巧将镶嵌好的DEM数据加载到QGIS中默认的灰度渲染可能无法有效展示地形特征。一套精心设计的色彩方案能极大提升数据的表现力和信息传递效率。QGIS提供了多种灵活的渲染器我们可以根据分析目的进行选择。1. 单波段伪彩色拉伸渲染这是最常用也最有效的地形展示方式。它将连续的高程值映射到一个渐变的色彩条带上。在图层属性中切换到“符号化”选项卡。将渲染类型从“单波段灰度”改为“单波段伪彩色”。在“渐变”下拉菜单中选择适合地形的色带如“地形”、“等离子”、“Spectral”等。“地形”色带从绿色到棕色到白色是最直观的。关键步骤是调整“最小值/最大值”和“裁剪”方式。不要直接使用“累积计数裁剪”这可能导致色彩对比度不足。建议在“加载”下拉框中选择“最小最大值”然后点击“加载”按钮读取实际数据范围。手动微调根据你的研究区可以手动设置一个更合理的范围。例如如果研究区是平原可以设置0-500米如果是山区可以设置0-4000米。这样能突出区域内的地形起伏。勾选“分类化”并设置一个合适的分类数如15-20这会使色带呈现为离散的色阶增强不同高程带的区分度。2. 山体阴影Hillshade叠加单纯的颜色渲染缺乏立体感。创建山体阴影图层能极大地增强地形的三维视觉效果。在QGIS菜单栏选择“栅格” - “分析” - “山体阴影”。输入图层选择你的DEM。设置“Z因子”垂直夸大系数。对于全球90米数据由于水平尺度大适当增大Z因子如2.0或3.0可以使地形起伏更明显。生成的山体阴影图层是一个灰度栅格。将其置于DEM图层之上并将混合模式改为“叠加”、“柔光”或“线性减淡”。调整其透明度通常在40%-70%之间使其与下层的彩色DEM完美融合。3. 等高线生成对于需要精确读图或进行地形特征提取如流域分析的场景等高线必不可少。在菜单栏选择“栅格” - “提取” - “等高线”。设置“间隔米”根据你的比例尺和需求来定如100米、250米。在“高级”选项中可以勾选“属性名称”将高程值写入生成的等高线矢量图层的属性表中。生成的等高线可以设置为细的、半透明的灰色或棕色线条叠加在彩色DEM上既能提供精确高程信息又不会过于喧宾夺主。一个专业地形图的图层顺序与样式示例底层彩色渲染的DEM单波段伪彩色地形色带。中层山体阴影图层混合模式叠加透明度50%。上层等高线图层线宽0.1mm颜色深灰色透明度30%。最上层可选重要的点、线矢量要素如河流、站点。通过这样的组合你得到的不再是原始数据而是一幅信息丰富、视觉上专业的地形图。4. 在ArcGIS Pro中的处理流程与对比对于ArcGIS Pro用户流程在逻辑上相似但操作界面和工具名称有所不同。这里简要对比关键步骤方便熟悉ArcGIS生态的用户快速上手。高程基准转换在ArcGIS Pro中可以使用“投影栅格”工具来完成基准转换但其设置更为隐蔽。在“分析”工具箱中找到“投影栅格”工具。输入栅格选择你的原始DEM。在“输出坐标系”中不能直接选择WGS84。你需要点击旁边的“地球”图标打开“坐标系”窗口。选择“地理坐标系” - “World” - “WGS 1984”。关键步骤在右侧的“XY坐标系详细信息”下方找到“垂直坐标系”选项卡。点击“转换”在弹出窗口中选择“从EGM96高度到WGS84椭球体高度的转换”。这样工具就会在转换水平坐标的同时应用垂直基准的改正。符号化与渲染ArcGIS Pro的符号化窗格功能非常强大。在“外观”选项卡下点击“符号系统”。在“主符号系统”下拉框中选择“拉伸”。在“色带”中选择一个合适的渐变色如“Elevation #1”。点击“直方图”按钮打开“拉伸”对话框。这里可以手动拖动滑块调整显示的最小最大值或者使用“类型”下拉菜单中的“最小最大值”、“标准差”等统计方式进行拉伸以优化视觉效果。要创建山体阴影可以使用“影像”选项卡下的“山体阴影”函数或者使用“3D Analyst”工具箱中的“山体阴影”工具生成独立图层再通过图层的混合模式进行叠加。数据管理对比表操作环节QGIS (开源方案)ArcGIS Pro (商业方案)核心要点与建议高程基准转换使用GDAL命令行gdalwarp灵活精准可批量脚本化。使用“投影栅格”工具需在垂直坐标系设置中指定转换。QGIS方案更适合自动化处理大批量数据ArcGIS Pro方案图形化操作更友好但批量处理需借助ModelBuilder。数据镶嵌使用gdal_merge.py或gdalwarp命令行或Processing工具箱中的“合并”工具。使用“镶嵌至新栅格”或“镶嵌数据集”工具。对于简单拼接两者皆可。对于需要复杂逻辑如像素优先级的大规模镶嵌ArcGIS的镶嵌数据集管理能力更强。山体阴影生成“栅格分析”-“山体阴影”工具或gdaldem hillshade命令行。“影像”选项卡-“山体阴影”函数或“3D Analyst”工具箱。QGIS的gdaldem命令行参数更丰富ArcGIS Pro的“函数”链可以实时渲染不生成中间文件。符号化灵活性通过“符号化”面板手动调整色带、拉伸、分类支持混合模式。“符号系统”窗格功能类似预置色带和模板丰富与整个平台风格统一。QGIS在自定义色带和与开源样式库如cpt-city结合上更有优势ArcGIS Pro在出图的整体美观度和与专题图模板的集成上更胜一筹。5. 进阶应用从DEM到水文分析与地形指数计算得到高质量的地形数据后它的价值才真正开始显现。MERIT DEM的一个主要设计目的就是支持高精度的水文分析。这里介绍两个最常用的衍生分析流向计算和地形湿度指数TWI计算。水文校正填洼与流向计算原始DEM可能存在凹陷点导致水流无法流出形成“内流盆地”。在进行水文分析前必须进行填洼处理。工具在QGIS中可以使用“Processing Toolbox”中的“GRASS” - “r.fill.dir”工具或者“SAGA GIS” - “填洼Wang Liu”工具。在ArcGIS Pro中使用“水文分析”工具箱中的“填洼”工具。输出填洼后的DEM以及流向栅格每个像元的水流方向常用D8算法用1,2,4,8,16,32,64,128八个值表示八个方向。计算地形湿度指数TWITWI是描述地形控制下土壤水分空间分布的重要指数公式为TWI ln(α / tanβ)其中α为上坡汇流面积β为坡度。流程使用填洼后的DEM计算坡度Slope栅格。使用流向栅格计算上坡汇流面积Flow Accumulation栅格。这个值代表流入每个像元的上游像元总数。利用栅格计算器套用TWI公式。在QGIS的Raster Calculator或ArcGIS Pro的“栅格函数”中都可以实现。# 以下是一个在Python环境中使用richdem库计算TWI的简化示例代码 # 这展示了核心算法的实现逻辑实际在GIS软件中通常有现成工具 import richdem as rd import numpy as np # 加载填洼后的DEM dem_filled rd.LoadGDAL(filled_dem.tif) # 计算坡度弧度制 slope rd.TerrainAttribute(dem_filled, attribslope_radians) # 计算D8流向 flow_dir rd.FlowDirection(dem_filled, methodD8) # 计算汇流面积 flow_acc rd.FlowAccumulation(flow_dir) # 计算TWI (注意处理坡度为0的情况通常加一个极小值) twi np.log((flow_acc 1) / (np.tan(slope) 0.001)) # 保存结果 rd.SaveGDAL(twi_output.tif, twi)注意实际计算中需要对汇流面积进行单位换算乘以像元面积并处理坡度为0的情况避免除零错误。不同的软件工具如SAGA GIS的“Topographic Wetness Index”工具、Whitebox Tools的“WetnessIndex”工具已经内置了这些处理逻辑通常比手动计算更稳健。通过这些衍生分析MERIT DEM就从静态的地形数据变成了动态的水文过程模拟的基础可以用于洪水风险区划、土壤侵蚀评估、生境适宜性分析等多个领域。我上次做一个流域的非点源污染模拟就是基于这套处理好的DEM生成了流向、汇流网络和TWI图作为模型的关键输入参数效果非常可靠。