1. 为什么选择GMT6.1来画地形起伏图地形起伏图这东西说白了就是把DEM数据里的高程值用颜色、阴影和光照效果“翻译”成一张能让人一眼看出哪里是山、哪里是谷、哪里是平原的图。我最早接触这类图的时候用的是ArcGIS那一套后来转到QGIS再后来因为要批量出图、要精确控制每一个像素的颜色和光照角度才彻底投奔了GMT。GMT全称Generic Mapping Tools在海洋地球物理和地学圈子里是公认的“出图神器”到了6.1版本它的现代模式Modern Mode已经非常成熟一条命令就能完成从数据到成图的全流程不用像老版本那样手动管理一堆中间文件。你可能会问现在Python的matplotlib、cartopy也能画地形图为什么还要折腾GMT我的体会是GMT在栅格数据的渲染质量和投影精度上依然是独一档的存在。尤其是画地形起伏图时GMT的grdimage配合grdgradient做光照阴影出来的效果非常自然山体的立体感很强不会像某些工具那样把山脊线画得跟刀切的一样生硬。而且GMT支持超过30种地图投影从常见的墨卡托到兰勃特等角圆锥再到各种方位投影基本上你能想到的它都有。这篇文章面向的读者是那些手里已经有了DEM数据、或者知道去哪里下载DEM数据但面对GMT那一堆参数不知道从何下手的人。我会从数据下载开始讲一直讲到出图时怎么避开那些让人抓狂的坑。整个过程我会尽量用“人话”来解释每个参数背后的逻辑而不是甩给你一堆命令让你自己猜。提示GMT6.1的安装方式有很多种Windows下推荐用官方提供的安装包Linux下可以用conda或者aptmacOS用Homebrew。安装完成后在终端输入gmt --version如果显示6.1.x就说明装好了。2. DEM数据下载与预处理2.1 主流DEM数据源对比与选择画地形起伏图第一步永远是搞到DEM数据。DEM就是数字高程模型简单理解就是一张巨大的表格每个格子记录了一个高程值。目前公开可下载的全球DEM数据主要有这么几种数据源分辨率覆盖范围获取方式适用场景SRTM30m/90m全球60°N~56°S美国地质调查局网站中大比例尺地形图ASTER GDEM30m全球日本METI/NASA区域地形分析ALOS World 3D30m全球日本JAXA高精度地形Copernicus DEM30m/90m全球欧洲空间局全球任意区域NASADEM30m全球60°N~56°SNASA LP DAACSRTM改进版我平时用得最多的是Copernicus DEM的30米数据原因是它的覆盖范围真正做到了全球无死角连极地地区都有而且数据质量比早期的SRTM干净不少空洞和噪声少。下载地址在欧空局的Copernicus Data Space Ecosystem上注册一个账号就能免费下载。如果你只是做一个小区域的教学演示SRTM 90米数据也够用文件小、下载快。下载的时候有个细节要注意不要一次性下载太大范围的数据。我见过有人直接下载整个亚洲的30米DEM结果文件好几个GBGMT读进去之后内存直接爆掉。正确的做法是先确定你的研究区域范围用经纬度框定一个矩形只下载这个矩形内的数据。比如你要画青藏高原东缘的地形大概范围是东经95°到105°、北纬28°到35°那就只下载这个范围。2.2 数据格式转换与裁剪下载下来的DEM数据通常是GeoTIFF格式GMT6.1可以直接读GeoTIFF但为了后续处理方便我习惯先把它转成GMT自己的网格格式NetCDF。转换命令很简单gmt grdconvert input.tif output.grd -V-V是打开详细输出让你看到转换进度。如果数据量比较大这一步可能需要等几分钟。转换完成后用gmt grdinfo output.grd查看一下网格的基本信息包括行列数、经纬度范围、高程最大最小值。这个信息很重要后面设置颜色表的时候要用到。有时候下载的DEM范围比你实际需要的要大这时候需要裁剪。GMT里裁剪网格用gmt grdcutgmt grdcut input.grd -R95/105/28/35 -Gcut.grd-R后面跟的是西经/东经/南纬/北纬注意GMT里经度范围是-180到180如果你研究的是西半球经度要写成负数。裁剪完之后再grdinfo看一眼确认范围对了再往下走。注意有些DEM数据的高程值是整数有些是浮点数。如果你发现出图后地形看起来“台阶感”很重多半是数据精度不够。这时候可以考虑用gmt grdsample做一次重采样把网格加密但要注意重采样不会增加真实信息只是让视觉上更平滑。2.3 数据空洞处理与异常值剔除DEM数据里经常会有空洞尤其是SRTM数据在陡峭山区和水体区域。空洞在GMT里表现为NaNNot a Number出图时这些区域会显示为白色或者透明。如果你不想看到这些空洞可以用gmt grdfill来填充gmt grdfill input.grd -An -Gfilled.grd-An表示用周围有效值的平均值来填充空洞。还有一种情况是数据里存在异常值比如某些像素的高程是-9999或者32767这些是标记值不是真实高程。处理方法是先用grdclip把超出合理范围的值截断gmt grdclip input.grd -Sa-500/0 -Sb9000/9000 -Gclipped.grd这条命令的意思是小于-500的值设为0大于9000的值设为9000。为什么是这两个阈值因为地球表面最低点死海大约在-430米最高点珠峰8848米留一点余量就够了。3. 地形起伏图的核心原理与参数解析3.1 颜色表地形的“翻译词典”地形起伏图好看不好看颜色表Color Palette TableCPT至少占一半功劳。GMT自带了很多经典的颜色表比如geo、topo、etopo1、relief这些都在GMT的share/cpt目录下。你可以用gmt makecpt命令基于这些内置颜色表生成自己的CPT文件gmt makecpt -Cgeo -T-8000/9000/500 -Z topo.cpt-Cgeo指定用geo颜色表-T-8000/9000/500表示高程范围从-8000到9000每500米一个色阶-Z表示连续渐变而不是分块。这里的高程范围要根据你实际数据的最大最小值来定不要照搬。比如你画的是某个小流域高程范围可能只有200到3000米那就写成-T200/3000/100。我个人的经验是色阶间隔不要太小。有人为了追求“平滑”把间隔设成10米结果颜色表里几千个色阶出图时颜色过渡反而显得脏。一般来说间隔设在100到500米之间比较合适具体看你的高程跨度。如果跨度是5000米间隔500米就是10个色阶视觉上层次分明如果跨度只有1000米间隔100米就是10个色阶同样合理。还有一个技巧是用-I参数反转颜色表。默认情况下低海拔是绿色、高海拔是白色或红色但有时候为了配合特定的视觉风格你可能想要反过来。gmt makecpt -Cgeo -I就能实现反转。3.2 光照阴影让山“立”起来光有颜色还不够地形图要让人看出立体感必须加光照阴影。GMT里做阴影用gmt grdgradientgmt grdgradient input.grd -A45 -Ne0.6 -Gshadow.grd这三个参数是核心-A45光照方位角45度表示光源从西北方向打过来。为什么是西北因为这是地图制图的传统惯例人眼已经习惯了这种光照方向如果改成东南方向山体会看起来像凹陷的坑而不是凸起的山。-Ne0.6-N表示用梯度模的指数归一化e表示指数0.6是归一化因子。这个参数控制阴影的强度值越大阴影越重。0.6是一个比较温和的设置适合大多数场景。-Gshadow.grd输出阴影网格文件。有时候你会看到别人用-Nt0.5t表示用梯度模的平方根归一化。这两种归一化方式的区别在于e是指数归一化t是平方根归一化。实测下来e方式在陡峭山区表现更好t方式在平缓地区更自然。你可以两种都试试看哪个顺眼就用哪个。提示如果你觉得阴影太重山体看起来像被墨泼过一样可以把归一化因子从0.6降到0.4或者0.3。反之如果觉得太平淡就往上加。这个参数没有绝对标准多试几次找到你满意的效果。3.3 投影选择不同场景用不同投影GMT支持的地图投影非常多画地形起伏图常用的有这几种墨卡托投影-JM适合低纬度地区高纬度变形严重。如果你画的是赤道附近区域用这个没问题。兰勃特等角圆锥投影-JL适合中纬度东西延伸的区域比如中国全图。标准纬线一般设在区域南北边界的1/6和5/6处。阿尔伯斯等面积圆锥投影-JB适合需要保持面积准确的场景比如计算某个区域的面积。通用横轴墨卡托-JU适合南北延伸的区域比如南美洲西海岸。选投影的核心原则是你的研究区域是什么形状就选什么投影。东西长南北短用圆锥投影南北长东西短用横轴墨卡托接近方形用等距圆柱或者墨卡托都行。不要小看投影选择选错了会让你的地图看起来“歪”得很别扭。4. 完整出图流程与命令详解4.1 现代模式下的脚本结构GMT6.1的现代模式让出图流程变得非常清晰。一个完整的脚本通常长这样#!/bin/bash gmt begin topo_map png,pdf gmt set FONT_ANNOT_PRIMARY 10p,Helvetica gmt set MAP_FRAME_TYPE plain gmt makecpt -Cgeo -T-500/6000/200 -Z topo.cpt gmt grdgradient dem.grd -A45 -Ne0.6 -Gshadow.grd gmt grdimage dem.grd -Itop shadow.grd -Ctopo.cpt -JM15c -R95/105/28/35 -Bafg gmt colorbar -DJBCw10c/0.5ch -BaflElevation (m) gmt endgmt begin和gmt end之间的所有命令会自动共享状态不需要像老版本那样手动指定输出文件。gmt begin topo_map png,pdf表示同时输出PNG和PDF两种格式PNG用于快速预览PDF用于印刷或投稿。4.2 关键参数逐个拆解gmt grdimage是出图的核心命令它的参数比较多我逐个解释dem.grd输入的DEM网格文件。-Ishadow.grd指定阴影网格。注意这里写的是-I后面直接跟文件名不是-Ishadow.grd。有些教程会写错导致阴影加载不上。-Ctopo.cpt指定颜色表文件。-JM15c墨卡托投影图幅宽度15厘米。15c表示15厘米你也可以用15i表示15英寸。-R95/105/28/35地图范围西经95度到东经105度南纬28度到北纬35度。-Bafg边框和网格线设置。a表示标注间隔自动f表示边框g表示网格线。更精细的写法是-B5/5:.Title:表示经纬度标注间隔5度标题写在底部中间。gmt colorbar是色标命令-DJBC位置在底部居中Bottom CenterJ表示图外。w10c/0.5c色标宽度10厘米高度0.5厘米。h水平放置。-Baf标注间隔自动。-lElevation (m)色标标签。4.3 出图后的检查与调整第一版图出来之后不要急着收工。我通常会检查这几个方面第一颜色过渡是否自然。如果发现某个高程段颜色跳变太厉害说明色阶间隔设得不合理回去调整makecpt的-T参数。第二阴影方向是否一致。有时候因为数据范围跨了多个投影带阴影方向会看起来不统一。这时候可以考虑用grdgradient的-D参数指定方向网格而不是用固定的方位角。第三标注是否清晰。如果经纬度标注太密或者太疏调整-B参数里的间隔值。标注字体大小用gmt set FONT_ANNOT_PRIMARY来改。第四色标是否遮挡地图内容。如果色标压住了地图上的重要区域用-DJBCo0c/1c把色标往下移1厘米。注意GMT的现代模式里所有gmt set命令必须在gmt begin之后、第一个绘图命令之前执行否则不生效。这个坑我踩过好几次明明设了字体大小但出图还是默认值后来才发现是gmt set写在了gmt begin前面。5. 常见问题与排查技巧实录5.1 数据读不进去怎么办最常见的问题是grdimage报错说“无法读取网格文件”。排查顺序是这样的先用gmt grdinfo dem.grd确认文件本身没问题。如果grdinfo都读不了说明文件损坏或者格式不对。检查文件路径。GMT对相对路径和绝对路径都支持但如果你在脚本里用了~表示家目录有时候会解析失败。建议统一用绝对路径。检查文件权限。Linux下如果文件没有读权限GMT会报错但错误信息可能很模糊。用ls -l看一眼权限位。如果是从GeoTIFF转过来的确认转换过程中没有报错。有时候GeoTIFF的坐标系信息不完整转换出来的grd文件范围会不对。5.2 出图后一片空白或者全黑这种情况通常是颜色表的问题。如果CPT文件里的高程范围和你数据的高程范围完全不重叠GMT就不知道该给每个像素上什么颜色结果就是一片空白或者全黑。解决办法是用gmt grdinfo看一下数据的实际最大最小值然后重新生成CPTgmt grdinfo dem.grd # 假设输出显示最小200最大4500 gmt makecpt -Cgeo -T200/4500/100 -Z topo.cpt还有一种可能是阴影网格的归一化因子设得太极端导致整个图要么全白要么全黑。把-Ne后面的值调到0.5左右试试。5.3 阴影效果不自然阴影看起来“假”通常有几个原因光照角度不对。默认的45度方位角在大多数情况下没问题但如果你画的是南北走向的山脉45度光照会让山脊线看起来很奇怪。这时候可以试试-A315让光源从东北方向打过来。归一化因子太大。-Ne0.8会让阴影非常重山体看起来像被墨汁泼过。降到0.4到0.5之间会自然很多。DEM分辨率太低。90米的数据画小区域地形阴影会显得很“糊”。如果条件允许尽量用30米数据。5.4 色标和地图对不齐色标和地图对不齐通常是因为色标的宽度和地图宽度不一致。gmt colorbar的w参数控制色标宽度gmt grdimage的-J参数控制地图宽度。把这两个值设成一样色标就会和地图等宽。比如地图是-JM15c色标就写w15c/0.5c。5.5 输出文件太大PDF输出如果包含高分辨率栅格文件可能会非常大。解决办法是在gmt begin里指定PDF的压缩级别gmt begin topo_map png,pdf gmt set PS_MEDIA A4 gmt set PS_PAGE_ORIENTATION landscape # ... 绘图命令 ... gmt end或者在gmt grdimage里用-Q参数降低栅格采样精度。但要注意降低精度会影响出图质量只建议在文件大小实在无法接受时使用。6. 进阶技巧与效率提升6.1 批量出图脚本模板如果你需要为多个区域出图手动改参数太慢了。我通常写一个bash脚本用循环遍历区域列表#!/bin/bash regions(95/105/28/35 100/110/25/32 90/100/30/38) names(east_tibet southwest central_tibet) for i in ${!regions[]}; do gmt begin ${names[$i]} png gmt makecpt -Cgeo -T-500/6000/200 -Z topo.cpt gmt grdgradient dem.grd -A45 -Ne0.6 -Gshadow.grd gmt grdimage dem.grd -Ishadow.grd -Ctopo.cpt -JM15c -R${regions[$i]} -Bafg gmt colorbar -DJBCw10c/0.5ch -BaflElevation (m) gmt end done这个模板的关键是把区域范围和输出文件名做成数组循环的时候用索引对应。这样一次就能出好几张图效率提升非常明显。6.2 用grdview画三维地形除了二维平面图GMT还能画三维地形。gmt grdview命令可以生成带透视效果的三维地形图gmt grdview dem.grd -JM15c -R95/105/28/35 -JZ5c -Ctopo.cpt -Ishadow.grd -Qm -Bafg -p135/30-JZ5c表示Z轴高度5厘米-Qm表示用网格线而不是面来渲染-p135/30表示视角方位角135度、仰角30度。三维图在展示地形起伏的宏观特征时非常直观但细节不如二维图清晰适合放在报告的概览部分。6.3 叠加水系和断层线地形起伏图如果只有DEM信息量还是有限。我经常会在上面叠加水系和断层线。水系数据可以从公开的HydroSHEDS下载断层线可以从全球活动断层数据库获取。叠加的方法是用gmt plotgmt plot rivers.shp -W0.5p,blue -R95/105/28/35 -JM15c gmt plot faults.shp -W1p,red,- -R95/105/28/35 -JM15c-W控制线宽和颜色0.5p表示0.5磅blue是颜色,-表示虚线。叠加的时候要注意图层顺序先画DEM再画水系最后画断层这样断层线会压在水系上面视觉层次更清晰。6.4 自定义颜色表GMT内置的颜色表虽然好用但有时候你需要根据特定的视觉风格自定义。CPT文件其实就是文本文件格式是“高程值 红 绿 蓝 高程值 红 绿 蓝”。你可以用文本编辑器直接改也可以用gmt makecpt生成后再微调。我常用的一个技巧是在低海拔用深绿色中海拔用浅绿色到黄色高海拔用棕色到白色这种配色方案在展示植被-地形关系时特别直观。7. 我踩过的那些坑第一个坑是忘记设置gmt set的位置。前面提过gmt set必须在gmt begin之后。我一开始不知道把字体设置写在了脚本最前面结果出图字体一直是默认的12磅改了跟没改一样。第二个坑是阴影网格和DEM网格范围不一致。有一次我用裁剪后的DEM出图但阴影网格还是用原始DEM生成的结果阴影和地形对不上山脊线错位了。后来养成习惯每次裁剪DEM之后重新生成阴影网格。第三个坑是色标标签里的单位没加。-lElevation和-lElevation (m)看起来差别不大但后者明显更专业。审稿人或者读者看到没有单位的色标会觉得你不够严谨。第四个坑是输出格式选错。PNG适合屏幕预览但分辨率有限PDF是矢量格式放大不糊但文件大。如果是投稿用建议同时输出PDF和PNGPDF用于印刷PNG用于在线预览。第五个坑是没有检查数据的坐标系。有些DEM数据用的是地理坐标系经纬度有些用的是投影坐标系米。GMT默认按经纬度处理如果你拿到的数据是投影坐标需要先用gmt grdproject转回地理坐标否则出图范围会完全不对。提示GMT的官方文档非常详细但全是英文而且命令参数太多初学者容易迷失。我的建议是先照着能跑通的脚本改改一个参数看一个效果慢慢就摸清每个参数的作用了。不要试图一次搞懂所有参数那不现实。8. 数据下载的替代方案如果你在Copernicus或者USGS下载数据时遇到网络问题还有一些替代方案。比如OpenTopography网站提供了多个DEM数据集的在线裁剪和下载服务你只需要框选区域它就会帮你裁好并打包下载。另外一些大学和科研机构也会在GitHub上分享处理好的区域DEM数据搜索“区域名DEMdownload”往往能找到惊喜。对于国内用户清华大学地球系统科学系维护了一个不透水面数据平台虽然主要是城市不透水面数据但他们的数据下载页面也提供了一些基础地理数据的链接。另外国家青藏高原科学数据中心提供了青藏高原区域的DEM数据下载注册后即可获取。下载数据时要注意数据许可协议。大多数科研DEM数据允许免费用于非商业用途但如果你要用于商业项目需要仔细阅读许可条款。有些数据要求你在发表成果时引用特定的论文这个引用格式一般在数据下载页面会有说明。9. 出图效率的优化建议最后分享几个提升出图效率的实操经验。第一把常用的参数写成GMT的配置文件。GMT支持在~/.gmt.conf里设置默认参数比如字体、边框样式、网格线颜色等。这样你就不用每次都在脚本里写一堆gmt set了。第二用gmt begin的-V参数控制输出详细程度。默认情况下GMT会输出很多中间信息如果你只想看错误信息用-Vqquiet模式。调试的时候用-Vddebug模式能看到每个命令的详细执行过程。第三把DEM数据和阴影网格缓存起来。如果你需要反复出图每次重新生成阴影网格很浪费时间。可以在第一次生成后把shadow.grd保存好后续出图直接加载。第四用GMT的-c参数指定多面板布局。如果你要在一张图上放多个子图比如不同区域的地形对比用gmt subplot可以自动排列不用手动计算每个子图的位置。第五输出前用gmt psconvert做最终转换。虽然gmt begin已经自动处理了格式转换但如果你需要更精细的控制比如设置DPI、裁剪空白边缘gmt psconvert提供了更多选项。比如gmt psconvert -A -Tg -E300表示自动裁剪边缘、输出PNG、分辨率300 DPI。我在实际使用中发现GMT6.1的现代模式虽然方便但有时候自动生成的中间文件会留在当前目录里时间长了目录会很乱。建议在脚本开头加一行rm -f gmt.*清理旧的临时文件或者在gmt begin里用-C参数指定一个专门的缓存目录。这个细节虽然小但能让你的工作目录保持整洁找文件的时候不至于抓狂。