免费获取学习方案
ARTICLE DETAIL

资讯详情

深耕编程基础知识与建站技术分享的一线实战洞察。

GMT 6.1 地形起伏图绘制全流程:从 DEM 数据到出图避坑指南

GMT 6.1 地形起伏图绘制全流程:从 DEM 数据到出图避坑指南 地形起伏图这东西说简单也简单说坑也多。GMTGeneric Mapping Tools作为地学圈里画图的老牌工具6.1版本在配色、投影和渲染上已经相当成熟但真正上手画一张能直接放进论文或报告里的地形起伏图从数据下载到最终出图中间踩的坑能写满一页纸。这篇内容就是把我自己反复折腾出来的完整流程摊开讲——从DEM数据去哪儿拿、怎么裁、怎么配色、怎么叠光照到出图时那些让人抓狂的报错和偏移问题。适合已经装好GMT、想画一张正经地形图但被各种参数劝退的人也适合之前用其他工具画图、想转到GMT命令行工作流的朋友。整篇不绕弯子直接给命令、给参数、给理由能抄的地方直接抄该解释的地方把原理说透。1. 先把数据这件事说清楚DEM从哪来、选哪个、怎么裁1.1 常见DEM数据源的取舍逻辑画地形起伏图第一步永远是搞到高程数据。GMT本身不带数据它只是个渲染引擎你得喂给它栅格。目前圈子里用得最多的几套数据各有各的脾气。SRTM系列是最经典的选择。30米分辨率的SRTM1和90米的SRTM3覆盖全球南北纬60度之间数据稳定、格式规整很多公开镜像都能直接下载。它的优点是拿来就能用缺点是高纬度地区没有覆盖而且城区和植被茂密区的高程会有偏差——因为SRTM测的是地表反射面不是裸地面。ASTER GDEM覆盖范围更广能到南北纬83度分辨率也是30米。但它的版本差异比较大早期V2版本噪声明显V3好了不少。如果你画的是大范围区域图ASTER的拼接痕迹有时候会在图上显出来需要额外做平滑。ALOS World 3D是后起之秀30米分辨率精度口碑比前两者都好尤其在山地和高差大的区域。缺点是下载流程相对麻烦需要注册账号逐个瓦片拿。TanDEM-X的90米全球版是商业数据里开放程度较高的精度极佳但免费版分辨率有限适合大尺度区域图。选哪个核心看三件事区域范围、纬度位置、精度要求。画中国东部SRTM1足够画青藏高原ALOS或TanDEM-X更稳画极地只能上ASTER或专门的极地DEM。提示不管用哪套数据下载后先用gdalinfo看一眼投影和NoData值。很多坑都是因为没检查元数据后面裁切和配色全乱套。1.2 下载与格式转换的实操细节以SRTM1为例公开镜像通常按经纬度瓦片组织命名规则是NxxExxx这种格式每个瓦片覆盖1度×1度。下载下来是.hgt格式GMT不能直接读需要转成.nc或.grd。转换用GDAL最省事gdal_translate -of GMT N30E120.hgt N30E120.grd如果区域跨多个瓦片先合并再转换gdalbuildvrt merged.vrt N30E120.hgt N30E121.hgt N31E120.hgt N31E121.hgt gdal_translate -of GMT merged.vrt merged.grd这里有个容易忽略的点.hgt是big-endian字节序GDAL能自动识别但如果你用其他工具处理过字节序错了高程值会变成天文数字。转换完用grdinfo确认一下数值范围正常陆地高程在-500到9000之间超出这个范围基本就是出问题了。1.3 裁切区域为什么不能直接画整幅很多人图省事下载完直接整幅渲染结果图上一大半是海或者无关区域出图比例失调。正确做法是先裁到目标范围。GMT6.1里裁切有两种方式。一种是用grdcut按经纬度矩形裁grdcut merged.grd -R115/125/30/40 -Gtarget.grd另一种是如果研究区是不规则形状比如某个流域可以用grdmask配合多边形文件做掩膜。这一步在画流域地形图时特别关键否则边界外的数据会干扰配色范围。裁切时要注意留一点缓冲。比如你最终要画115到125度裁的时候可以裁114.5到125.5给投影变换和边缘平滑留余量。GMT在投影边缘会有轻微重采样裁得太紧容易出现边缘锯齿。2. GMT6.1的配色与光照地形图好不好看全在这一步2.1 色标选择不是越鲜艳越好GMT6.1内置了大量CPTColor Palette Tablegmt show-cpt能列出全部。地形图最常用的是geo、topo、relief这几套。但内置色标直接用在专业图里往往太艳尤其是rainbow类颜色过渡生硬打印出来层次全糊。我的做法是基于内置色标改。比如以geo为基础把低海拔的绿色调暗一点高海拔的白色不要纯白改成略带灰的250/250/248这样在白色背景上不会飘。生成自定义CPTgmt makecpt -Cgeo -T-500/8000/100 -Z mytopo.cpt-T指定范围和步长-Z让颜色连续过渡。步长的选择有讲究如果区域高差大步长可以粗一点100米一档如果画的是丘陵地带高差只有几百米步长要细到20米甚至10米否则颜色分层太明显像等高线填充图而不是地形起伏图。2.2 光照渲染gradient和hillshade的区别地形起伏图的立体感来自光照。GMT里做光照有两条路grdgradient算梯度然后grdimage用-I参数叠加或者直接用grdimage的-I选项。grdgradient的经典用法gmt grdgradient target.grd -A45 -Ne0.6 -Gshadow.grd-A45是光照方位角45度意味着光从西北方向打过来——这是地图制图的惯例因为人眼习惯左上角光源看起来最自然。-Ne0.6是归一化参数0.6是个比较稳的值太高阴影太重太低立体感不足。然后渲染gmt grdimage target.grd -Ishadow.grd -Cmytopo.cpt -Jm -R115/125/30/40 -pdf topo这里有个关键坑-I叠加光照时如果CPT范围和数据范围不匹配会出现灰蒙蒙一层。原因是光照强度被错误地映射到了颜色上。解决办法是确保grdgradient输出的shadow文件是归一化的并且grdimage里不要额外加-E参数去调整亮度。2.3 水体处理海面不能跟着地形一起渲染如果区域包含海洋或湖泊直接用DEM渲染会把海底地形也画出来视觉上很乱。正确做法是先用水深或海平面掩膜。简单办法是用grdclip把海平面以下的值统一设为0grdclip target.grd -Sb0/0 -Gtarget_clip.grd这样海面就是平的再用CPT里0值对应的颜色通常是浅蓝填充。更精细的做法是单独准备一份水体多边形用pscoast叠加。pscoast在GMT6.1里的用法gmt pscoast -R115/125/30/40 -Jm -Df -S200/230/255 -W0.5 -O -K topo.ps-S填海面颜色-W画海岸线。注意-Df是全分辨率数据量大时渲染慢可以降到-Dh或-Di。3. 从命令行到成图完整脚本拆解与参数逐条解释3.1 一个可复用的完整脚本框架把前面所有步骤串起来下面是一个可以直接改区域和文件名就用的脚本#!/bin/bash # 地形起伏图绘制脚本 # 输入target.grd已裁切的高程数据 # 输出topo.pdf R115/125/30/40 JM6i CPTmytopo.cpt # 1. 生成色标 gmt makecpt -Cgeo -T-500/8000/100 -Z $CPT # 2. 计算光照 gmt grdgradient target.grd -A45 -Ne0.6 -Gshadow.grd # 3. 渲染地形 gmt grdimage target.grd -Ishadow.grd -C$CPT -J$J -R$R -K -P topo.ps # 4. 叠加海岸线和水体 gmt pscoast -R$R -J$J -Df -S200/230/255 -W0.3 -O -K topo.ps # 5. 加经纬网格 gmt psbasemap -R$R -J$J -Ba2f1 -O -K topo.ps # 6. 加色标 gmt psscale -C$CPT -Dx8c/2cw10c/0.5ch -B2000 -O topo.ps # 7. 转PDF gmt psconvert topo.ps -A -Tf这个框架里每一步的-K和-O是GMT的图层叠加机制-K表示还没画完后面还有-O表示这是叠加层不是新图。漏掉任何一个要么图出不来要么后一层把前一层覆盖掉。3.2 投影选择为什么Mercator不是默认最优脚本里用了-Jm即Mercator投影。它在低纬度地区变形小但高纬度会严重拉伸。画中国全图Mercator会让东北和新疆看起来比实际大很多。更合适的选择是Albers等面积投影或Lambert等角圆锥投影。GMT里Albers的写法-JB116/35/25/47/6i参数依次是中央经线、中央纬线、第一标准纬线、第二标准纬线、图幅宽度。标准纬线的选择有讲究一般取区域的1/6和5/6纬度位置这样变形最均匀。如果只是画小区域比如一个省Mercator和UTM差别不大用哪个都行。但如果是发表级图件投影说明要写清楚审稿人会看。3.3 网格和标注的细节控制psbasemap的-Ba2f1表示标注每2度一个网格线每1度一个。这个比例不是固定的要看图幅大小。图幅宽10度、标注间隔2度大概5个标注视觉上比较舒服。如果图幅只有3度宽标注间隔要降到0.5度。标注字体大小用-B的扩展参数控制-Ba2f1::.:WSen这里WSen表示只在西边和南边标注东边和北边不标。这是专业地图的惯例避免四周都是数字显得乱。网格线样式可以单独设gmt set MAP_GRID_PEN_PRIMARY 0.1p,gray0.1p是线宽gray是颜色。网格线太粗会抢地形的主体太细又看不清0.1到0.2磅之间比较合适。4. 出图避坑那些让我重画三次的问题4.1 数据范围与CPT范围不匹配导致的灰图这是最常遇到的坑。现象是地形图整体发灰颜色层次出不来像蒙了一层雾。根因是grdimage在渲染时如果CPT的范围-T参数和数据实际范围不一致超出CPT范围的值会被映射到CPT两端的颜色而光照叠加时这些溢出区域会异常。排查方法先跑grdinfo target.grd看数据的最小最大值再确认makecpt的-T范围覆盖了它。如果数据有负值比如海底CPT的-T下限要低于数据最小值。修复重新生成CPT范围比数据范围各留10%余量。4.2 光照方向反了山看起来像谷grdgradient的-A参数是方位角从正北顺时针算。-A45是光从东北方向来。但有些人习惯用-A315西北方向这两个出来的效果完全相反。判断标准很简单看山脊线。如果山脊亮、山谷暗说明光照方向对如果反过来把-A加180度。另外-Ne参数控制阴影强度。默认是1但实际用0.5到0.7效果最好。太高会让背光面全黑丢失细节。4.3 边缘白边裁切和投影的边界问题出图后图幅边缘有一条白边或者地形数据没有填满整个图框。这通常是两个原因一是裁切范围小于绘图范围。grdcut裁的是115到125但psbasemap的-R也是115到125理论上应该对齐但投影变换时边缘会有半个像素的偏移。解决办法是裁切时各方向多裁0.1度。二是NoData值被渲染成白色。DEM数据边缘常有NoDataGMT默认用白色填充。可以在grdimage里加-Q让NoData透明或者用grdclip把NoData设成背景色。4.4 色标位置和大小psscale的参数陷阱psscale的-D参数格式是x位置/y位置宽度/高度对齐方式。常见错误是位置单位搞混x8c/2c表示距左边8厘米、距底边2厘米但如果图幅本身只有10厘米宽色标就跑到图外了。更稳的做法是用相对位置-DjBCw10c/0.5co0c/1cjBC表示底部居中o是偏移量。这样不管图幅多大色标都在底部中间。色标的标注间隔用-B控制-B2000表示每2000米一个标注。如果高程范围是0到8000就是5个标注比较合适。范围小的时候要相应调小。4.5 字体和中文支持GMT6.1默认字体对中文支持不好如果图上要标中文地名需要设置字体gmt set FONT_ANNOT_PRIMARY SimHei gmt set FONT_LABEL SimHei但SimHei在Linux下不一定有可以用fc-list查系统里可用的中文字体。如果没有装一个开源的思源黑体然后指定字体文件路径。注意中文标注在psconvert转PDF时容易出问题建议先转成PNG确认效果再转PDF。5. 进阶玩法让地形图更有信息量5.1 叠加水系和道路纯地形图信息量有限叠加水系能立刻提升可读性。如果有矢量水系数据比如从OSM提取的河流用psxy叠加gmt psxy rivers.gmt -R$R -J$J -W0.5p,blue -O -K topo.ps河流线宽要细颜色用深蓝不要用亮蓝否则会抢地形。道路类似但通常用更细的线和灰色。5.2 等高线叠加在地形起伏图上叠等高线适合需要精确读数的场景。用grdcontourgmt grdcontour target.grd -R$R -J$J -C200 -W0.1p,gray -O -K topo.ps-C200表示每200米一条等高线。线要细要淡否则地形本身就看不清了。5.3 剖面线在地形图上画一条线旁边附剖面这是发表级图件常用的手法。先用psxy画一条剖面线然后用grdtrack提取沿线高程再用psxy画剖面图。echo 116 32 line.txt echo 124 38 line.txt gmt grdtrack line.txt -Gtarget.grd profile.txt然后在另一个面板里画profile.txt。GMT6.1支持多面板布局用-X和-Y参数控制每个面板的位置。5.4 批量出图改区域不改脚本如果要做多个区域的地形图把区域参数提取成变量用循环批量跑for region in 115/125/30/40 105/115/25/35; do R$region # 后续命令 done注意每次循环要改输出文件名否则会覆盖。可以用区域字符串拼文件名outfiletopo_$(echo $R | tr / _).pdf6. 我踩过的那些坑和最后的经验说几个文档里不会写、但实际一定会遇到的问题。第一个是内存。GMT渲染大区域高分辨率DEM时如果数据超过几千万个格点会直接吃满内存然后被系统杀掉。解决办法是先用grdsample降采样或者分块渲染再拼接。我画青藏高原全图时原始30米数据直接崩了三次降到250米才跑通。第二个是psconvert的-A参数。这个参数是自动裁剪白边但有时候会把图例或色标裁掉。如果发现输出PDF缺东西去掉-A试试。第三个是颜色在不同设备上的差异。屏幕上看着刚好的配色打印出来可能偏暗。建议出图前用gmt psconvert -Tg转一张高分辨率PNG在手机和电脑上都看一眼确认对比度够。第四个是版本兼容。GMT6.1和6.0的语法有细微差别比如grdimage的-I参数在6.1里支持扩展6.0不支持。如果脚本是从旧版本抄来的先查官方迁移文档。最后分享一个提高效率的习惯把常用参数写成GMT的配置文件。在~/.gmt/gmt.conf里设好默认字体、网格线样式、色标位置脚本里就不用每次都写。比如FONT_ANNOT_PRIMARY 10p,SimHei MAP_FRAME_TYPE plain MAP_GRID_PEN_PRIMARY 0.1p,gray这样每次跑脚本基础样式都是统一的只需要关注区域和配色这些变量。我现在的脚本从最早的80多行压缩到30行左右大部分重复参数都进了配置文件。地形起伏图这东西工具只是手段核心还是对数据的理解和对自己需求的清楚认知。想清楚图给谁看、要传达什么信息再去调参数比盲目试配色高效得多。
返回列表