1. 地形起伏图到底在画什么为什么值得用GMT折腾地形起伏图这东西乍一听像是地理专业的学生才需要碰的玩意儿但实际工作中它的出场频率远比想象中高。做区域规划的要拿它当底图写论文的要靠它撑起研究区概况那一节甚至做户外线路设计的也会用它来判断某段山脊的陡缓程度。说白了地形起伏图就是把一个区域的高程变化用颜色、阴影或者等值线的方式表达出来让人一眼就能看出哪里是山、哪里是谷、哪里平坦、哪里险峻。我最早接触这类图是用ArcGIS做的鼠标点一点、符号系统调一调确实上手快。但后来遇到几个问题一是出图风格千篇一律审稿人一看就知道是默认配色二是当研究区跨度大、投影变形明显的时候ArcGIS默认的渲染方式会让高纬度地区看起来“扁”得不对劲三是批量出图的时候手动操作效率太低。后来被同行安利了GMT全称Generic Mapping Tools中文一般叫通用制图工具才算是找到了一个既能精细控制又能脚本化批量处理的方案。GMT目前最新稳定版本是6.1.x系列相比5.x版本在语法上做了不少简化尤其是现代模式modern mode的引入让脚本写起来清爽了很多。这篇文章面向的是那些手里有DEM数据、想画出拿得出手的地形起伏图、但又不想在ArcGIS里反复调参数的人。不管你是地学专业的学生还是做环境评估的工程师只要跟着走一遍基本能掌握从数据获取到最终出图的完整链路。我会把踩过的坑、绕过的弯路都摊开来讲尽量让你少走冤枉路。2. 数据从哪来怎么选怎么处理2.1 DEM数据源的选择逻辑画地形起伏图核心输入就是DEM全称Digital Elevation Model中文叫数字高程模型。没有DEM后面的一切都是空谈。目前公开可获取的全球DEM数据主要有几个来源各有各的脾气。SRTM系列是最老牌的美国地质调查局发布覆盖全球南北纬60度之间的陆地分辨率有30米和90米两种。优点是获取方便、覆盖广、时间序列长缺点是高纬度地区没有而且30米版本在部分区域有数据空洞。ASTER GDEM是日本产的覆盖范围扩展到南北纬83度分辨率也是30米但它在植被茂密区和雪覆盖区的高程精度不如SRTM。ALOS World 3D是日本宇宙航空研究开发机构出的分辨率也是30米精度比前两者都好但下载流程稍微麻烦一些。再往上走就是商业数据了比如WorldDEM、LiDAR点云生成的DEM精度能到5米甚至更高但价格不菲。对于大多数科研和工程应用来说30米分辨率的公开数据已经够用了。如果你做的是城市尺度或者小流域分析可能需要考虑12.5米的TanDEM-X或者更高精度的数据。注意下载DEM之前一定要确认研究区的经纬度范围不同数据源的覆盖范围不一样别下完了才发现研究区有一半不在数据范围内。2.2 数据下载的实操路径以SRTM 30米数据为例最直接的获取方式是通过美国地质调查局的EarthExplorer平台。注册账号后在搜索框输入研究区名称或直接画矩形框选范围数据集选择“SRTM 1 Arc-Second Global”然后就能看到可下载的瓦片列表。每个瓦片覆盖1度×1度的范围命名规则是NxxExxx这样的格式比如N30E120表示北纬30度到31度、东经120度到121度。下载下来的是GeoTIFF格式可以直接用GMT读取。但如果你研究区跨了好几个瓦片就需要先做镶嵌。GMT本身有grdpaste命令可以拼接网格文件但更稳妥的做法是先用GDAL的gdal_merge.py做预处理因为GMT对GeoTIFF的投影信息读取有时候会出幺蛾子。gdal_merge.py -o merged_dem.tif -of GTiff N30E120.tif N30E121.tif N31E120.tif N31E121.tif拼接完之后用gdalinfo检查一下投影和范围是否正确。如果原始数据是地理坐标系经纬度而你想用投影坐标系出图还需要做重投影。这一步可以用gdalwarp完成gdalwarp -t_srs EPSG:32651 -r bilinear merged_dem.tif projected_dem.tif这里EPSG:32651对应的是UTM 51N带具体用哪个带号取决于你的研究区经度。重投影的时候重采样方法选双线性bilinear还是三次卷积cubic取决于你对高程精度的要求。一般来说地形起伏图用双线性就够了三次卷积虽然更平滑但计算量更大。2.3 数据预处理的几个关键动作拿到DEM之后别急着往GMT里塞先做几件事。第一是检查数据空洞SRTM数据在陡峭山区和城市区域可能有空洞表现为异常值通常是-32768或者NaN。可以用gdal_fillnodata.py做插值填补gdal_fillnodata.py -md 10 input_dem.tif filled_dem.tif第二是裁剪研究区如果下载的范围比实际需要的大很多用gdal_translate裁剪可以减小文件体积、加快后续处理速度gdal_translate -projwin 120.0 31.0 121.0 30.0 input_dem.tif clipped_dem.tif第三是确认高程单位大部分DEM的高程单位是米但也有个别数据源用分米或者英尺这个在元数据里会写清楚。如果单位不对后续的色标设置和等值线间隔都会跟着错。实操心得我习惯在预处理阶段就把DEM转成GMT的原生网格格式.grd用gdal_translate -of GMT就能完成。这样后续调用的时候不用每次都读GeoTIFF速度会快不少尤其是在做批量出图的时候。3. GMT6.1绘图的核心思路与脚本骨架3.1 现代模式与传统模式的取舍GMT6.1最大的变化就是引入了现代模式。传统模式下每个命令都是独立的你需要自己管理临时文件、自己控制图层叠加顺序脚本写长了容易乱。现代模式则把整个绘图过程当作一个会话session用gmt begin和gmt end包裹起来中间的命令自动按顺序叠加临时文件也由GMT自己管理。对于地形起伏图这种需要多层叠加底图阴影等值线标注的场景现代模式的优势非常明显。你不需要手动指定每个中间文件的名称也不需要担心图层顺序搞反。我现在的习惯是只要GMT版本在6.0以上一律用现代模式。gmt begin terrain_map png gmt grdimage dem.grd -Crelief -Id gmt grdcontour dem.grd -C100 -W0.5p,black gmt colorbar -DJBC -Baf gmt end上面这段就是一个最简化的地形起伏图脚本。grdimage负责渲染高程配色-Id表示自动计算阴影并叠加grdcontour画等值线colorbar加色标。三行命令一张基本的地形图就出来了。3.2 色标的选择与自定义GMT内置了不少色标比如relief、geo、topo、etopo1等都是为地形渲染设计的。relief偏暖色调低海拔偏绿、高海拔偏棕红geo偏冷色调适合表现海底地形topo是经典的绿-黄-棕-白过渡适合大多数陆地地形。但内置色标不一定符合你的审美或者期刊要求。GMT允许你用makecpt命令自定义色标。比如你想要一个从深绿到浅绿到黄色到棕色的过渡gmt makecpt -C0/2000/100 -Z -T0/5000/100 my.cpt这条命令的意思是在0到5000米范围内每100米一个色阶用-C指定颜色列表。更灵活的方式是直接写CPT文件格式是“高程值 红 绿 蓝 高程值 红 绿 蓝”这样一行一行定义。注意自定义色标的时候色阶的过渡要均匀否则会出现明显的色带banding。如果研究区高差不大比如只有几百米那色标范围就要相应缩小不然所有颜色都挤在一起看不出起伏变化。3.3 阴影效果的参数调校地形起伏图的立体感主要靠阴影hillshade来体现。GMT的grdimage命令通过-I选项控制阴影-Id是自动计算-Idirection/altitude可以手动指定光源方向和高度角。默认的光源方向是西北方向315度高度角45度。这个设置符合大多数人的视觉习惯因为自然界中阳光通常从偏北方向照射北半球。但如果研究区的地形走向比较特殊比如主要山脊是东西走向的那可能需要调整光源方向来突出山脊线。gmt grdimage dem.grd -Crelief -I315/45 -Q这里的-Q选项表示用阴影作为透明度调制而不是直接叠加。效果是阴影不会完全遮盖颜色而是让颜色在背光面变暗、在迎光面变亮看起来更自然。阴影的强度也可以通过-Id后面的参数微调比如-Id0.5表示阴影强度减半。如果觉得阴影太重、颜色被压得太暗就调小这个值反之则调大。4. 从零到一出图的完整实操流程4.1 环境准备与数据检查假设你已经装好了GMT6.1在终端输入gmt --version能正常显示版本号。如果还没装Windows用户可以直接下载安装包Linux用户用包管理器或者conda安装都行。conda的安装命令是conda install -c conda-forge gmt这个渠道的版本更新比较及时。数据方面假设你已经按照第2节的流程拿到了一个裁剪好的DEM文件命名为study_area.grd。先用gmt grdinfo看一眼基本信息gmt grdinfo study_area.grd输出会显示网格的行列数、经纬度范围、高程最小值和最大值。记下高程范围后面设置色标和等值线间隔的时候要用到。比如输出显示高程从200米到3200米那色标范围就可以设成200到3200等值线间隔可以设成200米或者250米。4.2 底图渲染与色标设置先做一个最基础的底图确认数据读取和渲染没有问题gmt begin base_map png gmt makecpt -Ctopo -T200/3200/200 -Z gmt grdimage study_area.grd -Id gmt colorbar -DJBCw10c/0.5c -Baf gmt endmakecpt的-T200/3200/200表示从200到3200米每200米一个色阶。-Z表示连续色标不加-Z则是离散色标。连续色标适合表现平滑的高程变化离散色标适合表现分级统计。colorbar的-DJBC表示色标放在底部居中w10c/0.5c指定色标的宽度和高度-Baf表示自动标注刻度。跑完这段脚本当前目录下会生成一个base_map.png。打开看看如果颜色过渡自然、阴影方向合理那底图就算成了。4.3 等值线叠加与标注底图之上叠加等值线可以让高程信息更精确。grdcontour命令的-C选项指定等值线间隔-W指定线宽和颜色-A控制标注。gmt grdcontour study_area.grd -C200 -W0.3p,gray40 -A400f8p,Helvetica,gray20这里-C200表示每200米画一条等值线-W0.3p,gray40表示线宽0.3磅、颜色为40%灰-A400表示每400米标注一次标注字体8磅、颜色20%灰。如果觉得等值线太密可以加大间隔如果觉得标注位置不理想可以用-An选项手动指定标注位置或者用-Af调整标注的避让策略。实操心得等值线的颜色不要用纯黑纯黑太抢眼会盖过底图的颜色变化。用灰色系gray30到gray50比较合适既能看清又不喧宾夺主。标注字体也不宜太大8到10磅就够了太大反而显得杂乱。4.4 地图边框与经纬度标注GMT的basemap命令负责画边框和经纬度标注。对于地形起伏图通常用-B选项指定标注间隔和样式gmt basemap -B10m/10m -BWSen -Lg120.5/30.2c0w100k-B10m/10m表示经纬度每隔10分标注一次-BWSen表示四边都画边框、标注放在西边和南边-L加比例尺位置在经度120.5、纬度30.2长度100公里。如果研究区范围比较小经纬度标注可能显得太密这时候可以改成-B5m/5m或者-B30s/30s。GMT支持度分秒的灵活组合d表示度m表示分s表示秒。4.5 输出格式与分辨率控制gmt begin的第二个参数指定输出格式可以是png、pdf、eps、jpg等。如果是要投稿的图建议用pdf或者eps矢量格式放大不糊。如果是网页展示或者PPT用png就够了。分辨率通过--DPI参数控制默认是300。如果要印刷级质量可以设成600gmt begin final_map pdf --DPI600 ... gmt end注意DPI设太高会导致文件体积急剧增大而且渲染时间也会变长。一般300 DPI已经能满足大多数期刊的要求600 DPI只在需要极高细节的时候才用。5. 常见问题与排查技巧实录5.1 数据读取报错与坐标系混乱最常见的问题就是GMT读不了GeoTIFF报错信息通常是“Unable to read file”或者“Unrecognized format”。这多半是因为GDAL库没有正确链接或者GeoTIFF的压缩格式GMT不支持。解决办法是用gdal_translate转成GMT原生格式gdal_translate -of GMT input.tif output.grd另一个常见问题是坐标系混乱。如果DEM是地理坐标系经纬度但你在grdimage里用了投影坐标系的参数出来的图会变形得亲妈都不认识。确认方法是用grdinfo看输出的“Projection”字段如果是“Geographic”就说明是经纬度。5.2 色标范围设置不当导致颜色失真色标范围设得太宽所有高程都挤在中间几个色阶里看起来一片糊设得太窄超出范围的高程会被截断成极值颜色看起来像贴了两块补丁。正确的做法是先看grdinfo输出的高程范围然后色标范围比实际范围略宽一点比如实际200到3200色标设成0到3500。如果研究区有负高程比如湖泊或者洼地色标范围要包含负值否则负高程会被渲染成最低色阶的颜色看起来像平地。5.3 阴影方向与地形走向不匹配默认的西北光源在大多数情况下没问题但如果研究区的主要山脊是南北走向的西北光源会让山脊的一侧过亮、另一侧过暗看起来不自然。这时候可以试试把光源方向调到东北45度或者正北0度看看哪个效果更好。gmt grdimage dem.grd -Crelief -I45/45多试几个方向选一个立体感最强、细节最清晰的。这个没有绝对标准以视觉效果为准。5.4 输出图片空白或只有边框这种情况通常是gmt begin和gmt end之间的命令没有正确执行或者图层叠加顺序有问题。排查方法是把脚本拆开一条一条命令单独跑看哪一步开始出问题。另外gmt begin之后的第一个绘图命令必须指定数据文件否则GMT不知道画什么。还有一种可能是输出格式不支持透明通道比如jpg不支持透明如果底图有透明区域就会显示成白色。换成png或者pdf就能解决。5.5 常见问题速查表问题现象可能原因解决方法数据读取报错GDAL链接问题或格式不支持用gdal_translate转GMT格式图变形严重坐标系不匹配用grdinfo确认投影类型颜色一片糊色标范围过宽缩小色标范围至高程实际范围阴影不自然光源方向与地形走向不匹配调整-I参数的光源方向输出空白命令执行失败或格式不支持透明逐条排查命令换png格式等值线太密间隔设置过小加大-C参数的间隔值标注重叠标注间隔过小加大-B参数的标注间隔实操心得GMT的报错信息有时候比较隐晦尤其是涉及投影转换的时候。我的习惯是在脚本开头加一行gmt set IO_N_HEADER_RECS 0这样可以避免一些因为头文件记录数不对导致的读取错误。另外gmt set FORMAT_GEO_MAP ddd:mm:ss可以统一经纬度标注格式避免出现度分秒混用的情况。6. 进阶技巧让地形图更专业、更耐看6.1 多图层叠加的透明度控制有时候需要在 terrain 底图之上叠加其他图层比如水系、道路、行政边界。如果直接叠加上层图层会完全遮盖下层。GMT的-t选项可以设置透明度比如-t50表示50%透明。gmt plot rivers.shp -W0.5p,blue -t30这样水系就会以半透明的方式叠加在地形之上既能看到水系走向又不影响地形的视觉表达。6.2 局部放大与插图如果研究区有一个重点区域需要放大展示可以用gmt inset命令在图上开一个插图窗口gmt inset begin -DjTRw4c/3co0.2c gmt grdimage dem.grd -Crelief -Id -R120.5/121.0/30.0/30.5 gmt basemap -B5m/5m -BWSen gmt inset end-DjTR表示插图放在右上角w4c/3c指定插图尺寸o0.2c指定偏移量。插图内部可以单独设置范围和标注间隔不受主图影响。6.3 批量出图的脚本化思路如果你需要为多个研究区出图手动改参数太累。可以把研究区名称、范围、色标范围写成变量用shell循环批量执行for region in area1 area2 area3; do gmt begin ${region}_map png gmt makecpt -Ctopo -T0/3000/200 -Z gmt grdimage ${region}.grd -Id gmt grdcontour ${region}.grd -C200 -W0.3p,gray40 gmt basemap -B10m/10m -BWSen gmt colorbar -DJBC gmt end done这样只要准备好每个研究区的DEM文件跑一遍脚本就能全部出图。如果色标范围需要根据每个区域的高程范围自动调整可以用gmt grdinfo提取高程极值再用awk或者sed动态生成makecpt的参数。6.4 色彩搭配的审美建议地形起伏图的色彩搭配没有绝对标准但有几个原则可以参考。第一是冷暖对比要适度低海拔用冷色绿、蓝高海拔用暖色黄、棕、红这样视觉上自然形成层次。第二是避免饱和度过高的颜色纯红纯绿纯蓝放在一起会显得刺眼用稍微灰一点的色调更耐看。第三是色标过渡要平滑不要出现明显的色带断裂。如果拿不准用什么配色可以参考GMT内置的relief、topo、geo这几个色标它们都是经过专业设计的适合大多数场景。如果期刊有特定要求比如要求灰度打印那就用gray色标或者自定义一个从浅灰到深灰的过渡。6.5 输出前的最终检查清单出图之前我习惯做一遍检查数据范围对不对、色标范围合不合理、阴影方向自不自然、等值线密不密、标注清不清楚、比例尺有没有、边框完不完整。这些检查花不了几分钟但能避免返工。另外如果是投稿用图还要确认字体是否嵌入、分辨率是否达标、文件格式是否符合期刊要求。有些期刊要求矢量图那就输出PDF或者EPS有些要求位图那就输出PNG或者TIFFDPI至少300。实操心得GMT的gmt begin命令支持--F参数指定输出文件名比如--Ffinal_map这样就不用每次改文件名了。另外gmt end之后可以用gmt clear清理临时文件避免工作目录里堆满中间文件。7. 我踩过的几个坑和最后的建议第一个坑是数据单位。有一次我拿到的DEM高程单位是分米但元数据里没写清楚我按米处理结果色标范围设成了0到3000实际高程只有0到300出来的图一片绿完全没有起伏感。后来用grdinfo看极值才发现问题。所以拿到数据第一件事就是确认单位。第二个坑是投影转换。有一次研究区跨了两个UTM带我图省事用了地理坐标系直接出图结果高纬度部分被拉伸得厉害山脊看起来比实际宽了一倍。后来老老实实做了重投影虽然多花了几分钟但图的几何精度对了。第三个坑是色标。我一开始喜欢用鲜艳的配色觉得好看。后来投期刊被审稿人提意见说颜色太饱和、打印出来效果差。后来改用偏灰的色调反而显得更专业。最后一个建议是GMT的文档虽然全但例子比较分散。我的习惯是建一个自己的脚本库把常用的命令组合、参数配置、色标文件都存下来下次遇到类似的任务直接改改就能用。这样积累下来出图效率会越来越高。地形起伏图这东西入门不难但要做到专业水准还是得在细节上花功夫。数据预处理、色标选择、阴影调校、标注布局每一个环节都有讲究。多画几张、多对比、多调整慢慢就能找到感觉。
