☰
GMT绘图工作流:坐标系、投影与数据结构深度解析
2026/10/1 5:24:39 网站建设 项目流程

1. 这不是一本“笔记”,而是一套可复用的GMT绘图工作流

GMT——Generic Mapping Tools,这个名字听起来像某个冷门地理软件的缩写,但对地球科学、海洋学、气象学、构造地质甚至行星科学领域的从业者来说,它几乎是刻在DNA里的工具链。我第一次接触GMT是在博士第三年,导师甩给我一个.cpt文件和三行命令,让我把一组海底热流数据画成带地形叠加的彩色剖面图。当时连pscoast和psxy的区别都分不清,硬着头皮啃了两周官方手册,才明白:GMT不是“点选式”绘图软件,它是一套以命令行为核心、以脚本为骨架、以精度和复现性为生命线的制图系统。所谓“GMT绘图笔记”,绝不是随手记下的零散命令集合,而是把多年项目中反复验证过的坐标系处理逻辑、投影陷阱规避方法、多图层合成策略、色彩映射调试经验、以及最关键的——如何让一张图在论文投稿、会议展板、野外汇报三种场景下都“稳得住”的全流程沉淀。它解决的不是“怎么画出来”,而是“怎么画得准、画得清、画得久、画得别人能复现”。适合三类人:刚入门被-J参数绕晕的研究生;卡在grdimage和grdcontour图层遮盖问题上的工程师;还有那些手头有十年老脚本、但一升级GMT版本就报错、不得不重写的资深用户。下面展开的,是我把2015–2024年间在南海构造图、西太平洋温盐剖面、青藏高原重力异常等十余个实际项目中踩过的坑、调过的参、存下来的模板,一条命令一条命令拆解给你看。

2. GMT绘图的本质:坐标系、投影与数据结构的三角博弈

2.1 GMT不是“画图软件”,而是“空间数据编排引擎”

很多人初学GMT时最大的认知偏差,是把它当成Origin或Matplotlib的替代品——以为只要数据格式对、命令敲对,图就出来了。错。GMT真正的底层逻辑,是空间数据的坐标系编排(Coordinate System Orchestration)。它不关心你数据里“温度”是摄氏还是华氏,但它极度敏感于“经度0°到底在哪”、“纬度90°是否对应北极点”、“你的网格数据第一行是南还是北”。这决定了所有后续操作的根基是否牢固。

举个最典型的例子:你有一份来自CMEMS的全球海表温度网格数据(NetCDF格式),lat维度从-89.5到89.5,步长1°,lon从-179.5到179.5,步长1°。你直接用grdimage画出来,发现太平洋中部有一条明显的竖直断裂线。这不是数据问题,而是GMT默认将lon视为从0°到360°的连续坐标,而你的数据是-180°到+180°标准地理坐标。GMT在内部做插值或裁剪时,会把-179.5°和+179.5°当作相距359°的两个点,而非相邻的两个点。结果就是——它试图在-179.5°和+179.5°之间“拉伸”出359°的空白,形成那条刺眼的缝。

解决方案不是改数据,而是告诉GMT:“我的经度是-180/+180体系”。命令里加-R-180/180/-90/90明确范围,再用-r(registration)参数指定网格注册方式。但这里又埋了个坑:-r有-r(pixel registration)和-r g(gridline registration)之分。前者认为每个网格值代表该像素中心点的值,后者认为每个网格值代表该网格线交点的值。CMEMS数据是gridline型,若误用-r,整个图像会整体偏移半个网格——在1°分辨率下就是0.5°,相当于55公里。我在南海热液区定位图里就因此把一个已知热液喷口标偏了整整一个构造单元,返工三天。

提示:判断NetCDF数据注册方式的最快方法,用gmt info -C yourfile.nc,看输出里x_min/x_max和y_min/y_max是否与nx/ny计算出的理论边界一致。若x_min = lon0 - dx/2,则是pixel;若x_min = lon0,则是gridline。

2.2 投影选择不是“选个好看样式”,而是定义空间关系的数学契约

GMT提供超过30种地图投影,但新手常犯的错误是:看到-JM(墨卡托)画出来方方正正就选它,看到-JS(球面)画出来圆润就以为更“真实”。投影的本质,是在二维平面上表达三维球面时,对面积、形状、距离、方向四者进行有意识的取舍。没有“最好”的投影,只有“最适合当前任务”的投影。

比如画全球地震分布图,目标是展示震中空间聚集性。这时必须选等面积投影,如-JQ(等积方位投影)或-JE(等积圆锥投影)。因为如果用墨卡托(-JM),高纬度地区会被严重拉伸,格陵兰岛看起来比非洲还大,导致你误判北欧地震活动密度远高于赤道——实际上只是投影把单位面积放大了数倍。

再比如画横跨太平洋的板块运动矢量图,关键是要准确表达矢量方向。这时必须选等角投影(保形投影),如-JT(横轴墨卡托)或-JL(兰伯特等角圆锥)。因为只有等角投影能保证任意一点上,所有方向的局部角度关系不变。我曾用-JE画过一次太平洋板块旋转矢量,结果发现夏威夷热点链的走向在图上歪了15°,就是因为等积投影扭曲了方向关系。

实操中,我给自己定了一条铁律:先问目的,再选投影。

  • 要比大小(如沉积物厚度分布)→ 等面积投影(-JQ,-JE,-JB)
  • 要看方向(如地壳应力场、洋流路径)→ 等角投影(-JT,-JL,-JM)
  • 要跨极区(如南极冰盖变化)→ 方位投影(-JA,-JS)
  • 要局部高精度(如城市地质填图)→ 横轴/斜轴投影(-Jt,-Jm)

而且,投影参数必须与数据坐标系严格匹配。常见错误是:数据是WGS84椭球体(-epsg 4326),却用-JE默认的球形地球参数计算。GMT 6.4之后引入了-epsg参数,强烈建议所有新项目统一用-epsg 4326显式声明,避免旧版GMT默认用球形带来的厘米级偏差——这点在毫米级GPS形变图里就是致命误差。

2.3 数据结构决定绘图效率:为什么你的grdimage总比别人慢三倍

GMT绘图速度差异,70%源于数据结构选择。很多人把GeoTIFF或NetCDF直接喂给grdimage,发现渲染要等半分钟。其实GMT最高效的数据格式是native grd格式(二进制网格文件),它没有元数据解析开销,内存映射读取极快。

转换方法很简单:gmt grdconvert input.nc output.grd。但关键在grdconvert的参数。默认转换会保留原始数据所有精度(float64),但多数科研图根本不需要15位小数。用-f参数强制降为float32:gmt grdconvert -f=bf input.nc output.grd。-f=bf表示binary float,文件体积减半,内存占用降低,读取速度提升40%以上。我在处理10GB的全球重力异常网格时,用float64加载需2.3秒,float32仅1.4秒——别小看这0.9秒,当你需要循环调试100次配色方案时,就是1.5分钟的差别。

另一个隐形杀手是网格分辨率与绘图区域的匹配度。假设你要画中国东部陆架区(118–124°E, 28–34°N)的水深图,而你的全球水深网格是1弧分(~1.8km)。GMT会先把整个全球网格读入内存,再裁剪出目标区域——白耗99%内存。正确做法是先用gmt grdcut裁剪:gmt grdcut global_topo.nc -R118/124/28/34 -Gchina_shelf.grd,再用grdimage china_shelf.grd。实测下来,内存峰值从12GB降到1.8GB,启动时间从8秒降到1.2秒。

最后,别忽视网格填充值(NaN handling)。GMT默认用-n参数控制插值,但很多老脚本忽略它。若你的网格边缘有大片NaN(如海洋区域无地形数据),-n设为-n+b(bilinear插值)会导致边缘模糊成一片灰雾。我固定用-n+c(cubic插值)+--GMT_COMPATIBILITY=6,确保插值锐利且兼容新版。

3. 核心绘图模块的深度拆解:从单图到复合图的工业级实现

3.1grdimage:不只是“上色”,而是空间数据的光学编码

grdimage是GMT最常用也最容易被低估的命令。它表面是给网格上色,实质是将数值域映射到视觉域的编码过程。颜色不是装饰,而是信息载体;色标不是摆设,而是解读钥匙。

首先,色标(CPT)必须与数据统计特征强耦合。我见过太多人直接用gmt makecpt -Cjet给重力异常数据上色,结果正负异常全挤在两端,中间大片区域显示为同一绿色——因为jet是线性色标,而重力异常呈双峰分布,±10 mGal只占数据5%,±100 mGal却占90%。正确做法是:先用gmt grdinfo -L1 data.grd获取数据1%和99%分位数,再用-Fr参数生成截断色标:gmt makecpt -Ccoolwarm -Fr-120/120/20 -T-120/120/20 > gravity.cpt。-Fr指定范围,-T指定色标断点,这样±120 mGal以外的极端值被截断,中间每20 mGal一个色阶,视觉分辨力提升3倍。

其次,透明度(alpha)不是为了“好看”,而是解决图层遮盖。比如在地形图上叠加热液点位(psxy),若直接画,点会被山体阴影盖住。解决方案是:grdimage输出时加-Q(启用alpha通道),并用-I+d生成山体阴影(-I是illumination),再用-Q让阴影半透明。完整命令:gmt grdimage topo.grd -I+d -Q -Ctopo.cpt。这样地形有立体感,又不会完全挡住下方的点。

最后,grdimage支持多波段合成,这是很多人不知道的隐藏技能。比如你有遥感影像的R/G/B三个波段网格(red.grd, green.grd, blue.grd),可以用gmt grdimage red.grd green.grd blue.grd -C直接合成真彩色图,无需先用GDAL转TIFF。前提是三个网格必须完全同源、同分辨率、同范围。我用这招快速生成了南海珊瑚礁的RGB合成图,比用QGIS导出快5倍。

3.2psxy:点线面的精准落位,容不得0.1秒的时间差

psxy负责绘制离散数据(点、线、面),但它的精度控制远超想象。最常见的需求是画地震震中,但如果你直接echo "121.5 25.3 5.2" | gmt psxy -Sc0.1c -Gred,会发现所有点都挤在左下角——因为你没指定-R和-J!psxy必须与主图的坐标系完全一致,否则坐标解析失败。

更隐蔽的坑在时间序列数据的X轴处理。比如画某台站十年地壳形变(时间,东向位移),数据是2014-01-01 2.3格式。GMT不认日期字符串,必须转为儒略日(Julian Day)。用gmt convert -o0,1 -f0T data.txt自动转换,-f0T表示第0列是时间格式。但注意:-f0T默认按UTC时区解析,若你的数据是北京时间(UTC+8),必须先用awk '{print $1" "$2" +0800"}' data.txt | gmt convert -f0T,否则所有时间偏移8小时——我在处理GPS数据时因此把一次地震前兆信号标错了3天。

对于线要素(如断层迹线),-W参数的宽度单位极易混淆。-W1p是1点(point),约0.35mm;-W1c是1厘米;-W1i是1英寸。但关键在-W的端点样式:-W1p+miter(尖角)、-W1p+round(圆角)、-W1p+bevel(斜切)。画地质界线必须用+miter,否则断层交汇处出现难看的缺口;画河流用+round更自然。我曾因用错+bevel,让一幅构造图里的郯庐断裂带看起来像被锯齿刀切过。

面要素(如行政区划)要用-L闭合线。但-L要求数据首尾坐标严格相等,否则会画出一条诡异的斜线。用gmt spatial -E自动闭合:gmt spatial boundaries.xy -E > boundaries_closed.xy。这个小命令救了我无数张政区图。

3.3pscoast:海岸线不是背景,而是空间基准的校准尺

pscoast常被当作“加个海岸线”的快捷命令,但它其实是GMT空间精度的终极校验器。当你发现pscoast画出的海岸线和你的地震点不重合,问题99%不在pscoast,而在你的数据坐标系或投影参数。

pscoast有四级分辨率:-Dc(crude)、-Dl(low)、-Di(intermediate)、-Dh(high)、-Df(full)。别盲目用-Df——它文件大、加载慢,且在小比例尺图上反而糊成一团。我的经验是:全球图用-Di,区域图(如东海)用-Dh,城市图用-Df。但更重要的是-A参数:-A1000表示只画面积大于1000 km²的岛屿。若你画南海诸岛,必须设-A1,否则南沙群岛大部分岛礁被过滤掉。

pscoast的-G(陆地填充)和-S(海洋填充)颜色,必须与grdimage的底图协调。我固定用-G240(浅灰)填陆地,-S255(纯白)填海洋,这样后续叠加热液点、断层线时,视觉层次清晰。曾有人用-Gwhite,结果白底上画白线,整张图变成“找不同”游戏。

最值得强调的是pscoast的大地水准面校正。GMT默认用WGS84椭球体,但某些高精度项目(如卫星测高)需用EGM2008大地水准面模型。这时要配合gmt grdgradient生成坡度图,并用-I+d调制光照,让海岸线与地形阴影无缝融合。这步能让图的专业感跃升一个档次。

3.4 复合图(Multi-panel):不是拼图,而是信息流的叙事设计

GMT 6.x的pstext和pslegend已足够强大,但真正工业级复合图的核心是gmt subplot。它不是简单把几张图并排,而是构建一个坐标系矩阵,让每张子图共享全局参考框架。

比如画一幅“构造-重力-磁力”三联图,传统做法是分别生成三张PS文件,再用Ghostscript拼接。缺点是:三张图的经纬度刻度无法对齐,比例尺不一致,图例位置难协调。用gmt subplot则完全不同:

gmt begin multi_panel pdf gmt subplot begin 1x3 -A+t"南海北部湾构造-重力-磁力综合图" -M0.2c -Fs18p gmt subplot set 0,0 -R108/112/18/22 -JQ110/15c gmt grdimage struct.grd -Cstruct.cpt gmt pscoast -Dh -G240 -S255 gmt subplot set 0,1 -R108/112/18/22 -JQ110/15c gmt grdimage gravity.grd -Cgravity.cpt gmt pscoast -Dh -G240 -S255 gmt subplot set 0,2 -R108/112/18/22 -JQ110/15c gmt grdimage mag.grd -Cmag.cpt gmt pscoast -Dh -G240 -S255 gmt subplot end gmt end

关键在-R和-J参数全局统一,确保三张图地理位置绝对对齐;-M0.2c设置子图间距;-Fs18p统一字体大小。这样生成的PDF,用Adobe Acrobat测量任意两点距离,三张图结果一致——这才是科研图的底线。

subplot还支持跨子图标注。比如要在三张图顶部中央加一个总标题,用gmt pstext <<< "0.5 1.05 18 0 0 CM 南海北部湾综合解释",坐标0.5是相对子图宽度的50%,1.05是子图上方5%位置,完美居中。

4. 实操避坑指南:那些手册里不会写的血泪教训

4.1 版本陷阱:GMT 5 vs GMT 6,不只是命令名变化

GMT 6是重大架构升级,但很多教程仍基于GMT 5。最痛的兼容性问题是默认坐标系变更。GMT 5默认-R是-Rxmin/xmax/ymin/ymax,GMT 6默认-R是-Rxmin/xmax/ymin/ymax+r(+r表示region mode)。这意味着同样-R120/122/22/24,GMT 5画的是矩形区域,GMT 6可能画成空集——因为+r模式要求xmin<xmax且ymin<ymax,若你数据是南纬在前(如-R120/122/-24/-22),GMT 6会报错。

解决方案:所有脚本开头加--GMT_COMPATIBILITY=5,或统一用-R显式声明模式:-R120/122/22/24(GMT 5风格)或-Rg120/122/22/24(GMT 6地理模式)。

另一个隐形炸弹是-V参数。GMT 5的-V是verbose模式,GMT 6的-V是version查询。若你在GMT 6脚本里写gmt pscoast -V -Dh,它会直接打印版本号然后退出,后面命令全失效。必须改为-Vq(quiet verbose)或删掉-V。

4.2 字体与中文支持:别让标题变成方框

GMT原生不支持TrueType字体,中文显示是老大难。网上流传的“修改fontpath”方案在GMT 6.4+已失效。正确解法是:用Ghostscript后处理。

步骤:

  1. 绘图时用ASCII字体(如-F+f14p,Helvetica-Bold)占位
  2. 生成EPS文件:gmt psconvert -A -Tf map.ps
  3. 用gs -sDEVICE=pdfwrite -dCompatibilityLevel=1.4 -dPDFSETTINGS=/prepress -dEmbedAllFonts=true -dSubsetFonts=true -dColorImageDownsampleType=/Bicubic -dMonoImageDownsampleType=/Bicubic -dGraphicsAlphaBits=4 -dTextAlphaBits=4 -sOutputFile=map.pdf map.eps
  4. 用pdftk map.pdf fill_form chinese_fields.fdf注入中文(需提前准备FDF表)

更实用的方案是:用Python的matplotlib生成含中文的图例,保存为PNG,再用gmt image命令嵌入主图。虽然多一步,但100%可靠。我在所有正式投稿图中都采用此法。

4.3 内存溢出与崩溃:当grdview吃光你的32GB内存

grdview用于3D透视图,但极易内存爆炸。根本原因是它默认将整个网格加载到GPU显存。解决方案有三:

  • 降采样:gmt grdsample topo.grd -I2m -Gtopo_2m.grd(2弧分降采样)
  • 分块渲染:用gmt grdview topo.grd -JQ110/15c -R108/112/18/22 --MAP_FRAME_TYPE=plain -p135/30 -Baf -Ctopo.cpt --PS_PAGE_ORIENTATION=landscape > view.ps,其中-p135/30指定视角,--PS_PAGE_ORIENTATION=landscape强制横版节省内存
  • 换后端:export GMT_RENDERER=cairo,Cairo后端比默认的PostScript更省内存

我处理过一份1°×1°全球地形,grdview直接OOM,改用Cairo后端+2°降采样,内存从32GB压到4GB,渲染时间从崩溃到23秒。

4.4 颜色管理:为什么你的PDF在打印机上变成灰蒙蒙

GMT生成的PDF默认是DeviceRGB色彩空间,但印刷厂要求CMYK。直接转换会色偏。正确流程:

  1. GMT绘图时用-C指定CPT,确保色标在sRGB色域内(避免-Cvik这种超广色域)
  2. 生成PDF后,用gs -sDEVICE=pdfwrite -sProcessColorModel=DeviceCMYK -sColorConversionStrategy=CMYK -dOverrideColorSpace=/DeviceCMYK -o map_cmyk.pdf map.pdf
  3. 用pdfinfo map_cmyk.pdf确认Color space: DeviceCMYK

曾有篇论文彩图因未转CMYK,印刷后所有蓝色变成紫灰色,主编要求重印——就因为漏了这一步。

5. 工程化实践:让GMT脚本从“能用”到“可维护、可复现、可协作”

5.1 脚本结构化:告别“复制粘贴式”绘图

一个合格的GMT脚本,必须包含四个区块:

  1. 参数区:所有可配置项集中定义,如REGION="-R118/124/28/34"、PROJ="-JQ121/15c"、CPT="topo.cpt"
  2. 数据预处理区:grdcut、grdmath、gmt convert等,确保输入数据符合绘图要求
  3. 主绘图区:grdimage、psxy、pscoast等核心命令,用gmt begin/end包裹
  4. 后处理区:psconvert转格式、ghostscript优化、pdfcrop裁边

这样做的好处:同事要复现你的图,只需改参数区几个变量;审稿人质疑某条断层位置,你能在30秒内重新生成带坐标的debug图。

5.2 版本控制与复现性:.gmt配置文件的妙用

GMT支持.gmt配置文件(放在$HOME/.gmt/),里面可以定义:

  • GMT_DEFAULT_PEN=1.5p,black(统一画笔)
  • GMT_FONTSIZE=12p(统一字号)
  • GMT_DIR=/path/to/my/cpts(自定义CPT路径)

更重要的是GMT_SESSION_NAME。每次运行GMT会创建唯一session ID,若你希望多次运行结果完全一致(如蒙特卡洛模拟),设GMT_SESSION_NAME=fixed_seed,GMT会禁用随机数,确保grdnoise等命令结果可复现。

5.3 自动化测试:用diff验证图件一致性

我为每个重要脚本写了测试用例:

  • 生成一张标准图(ref_map.pdf)
  • 修改代码后,生成新图(new_map.pdf)
  • 用pdf2png ref_map.pdf ref.png && pdf2png new_map.pdf new.png && compare -metric RMSE ref.png new.png null:
  • 若RMSE>0.1,说明图形有实质性变化,触发人工审查

这套机制帮我在GMT 6.5升级时,提前发现-I+d光照算法变更导致的阴影偏移,避免了300+张图返工。

5.4 团队协作:GMT脚本的文档化规范

我们团队约定:每个.sh脚本开头必须有三段注释:

  1. Purpose:一句话说明图的用途(如“用于Fig.3,展示南海西南次海盆扩张脊的磁异常条带”)
  2. Input:列出所有输入文件及来源(如input_grd: GEBCO_2023.nc, downloaded from https://www.gebco.net)
  3. Output:明确输出文件名、格式、存放路径(如output: fig3_southwest_spreading.pdf in ./figures/)

并且,所有CPT文件必须附带README.md,说明色标物理意义(如coolwarm.cpt: red=positive gravity anomaly (mGal), blue=negative)。这些看似琐碎,却让新人三天内就能接手维护。

6. 最后分享一个小技巧:用GMT做“动态图”的低成本方案

GMT本身不支持动画,但你可以用gmt grdmath生成时间序列网格,再用ffmpeg合成视频。例如,画台风路径演变:

  1. 将每小时台风位置存为typhoon_000.grd,typhoon_001.grd...
  2. 用gmt grdmath "typhoon_000.grd typhoon_001.grd OR typhoon_002.grd OR ... = typhoon_all.grd"生成累积路径
  3. 用for i in {0..120}; do gmt grdimage typhoon_${i}.grd -Cpath.cpt -R... -J... -B... -P > frame_${i}.ps; done
  4. psconvert转PNG,ffmpeg -framerate 2 -i frame_%03d.png -vcodec libx264 -pix_fmt yuv420p typhoon.mp4

我用这招做了2018年台风“山竹”的72小时路径动画,文件仅8MB,比用ArcGIS导出小10倍,且所有坐标系100%精确。科研可视化,精度永远比炫酷重要。

我在实际使用中发现,最浪费时间的从来不是学命令,而是搞清楚“为什么这张图在导师电脑上正常,在我电脑上错位”。后来我才明白,GMT不是软件,而是一套空间思维训练。当你开始习惯用-R思考范围、用-J思考投影、用-C思考信息编码,你就不再是在“画图”,而是在构建一个可验证、可传播、可传承的空间知识系统。这大概就是所谓“笔记”的真正分量。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询