GMT在绘制辽宁地形图中的应用

孔祥雪，王万宁

（辽宁省地震局，辽宁 沈阳 110034）

0 引言

GMT，全称Generic Mapping Tools，一般翻译成“通用制图工具”。其是拥有80多种命令行工具的开源绘图软件，支持30多种地图投影，自带模块包括海岸线、河流、国界等数据信息，用于绘制不同类型的专题地图[1-2]。GMT是纯命令行的输入方式，所有的绘图操作都是通过命令行代码执行的。命令行相对于图形界面的优势在于内存占用更低；精确控制图形的显示，例如线条宽度、圆的大小；方便写成脚本，具有可批量处理、可重复、可自动化的特点。GMT将绘图及数据处理功能划分到不同的模块，各个模块之间相互独立且代码量少，易于更新和维护；绘图步骤之间以及涉及的数据类型独立性强，因而应用领域较广；可在脚本语句中调用一系列程序，进而绘制更为复杂的图片。GMT输出的图件格式为矢量图片格式，具有可以任意放大缩小而不失真的特点。GMT充分利用这种特性，可以生成高质量的矢量图件，并可以很容易地转换为其他图片格式。基于上述优点，GMT在地球科学领域得到了广泛的应用[3-4]。

目前GMT还提供了与MATLAB兼容的程序接口，可以在MATLAB中直接调用GMT命令；而GMT命令生成的结果，如格网数据（GRD）、表格数据 （TABLE）、颜色表 （CPT）、文本文件（TXT）、图片（PSFILE）等，都可以作为MATLAB变量调用运算；此外，MATLAB中的矩阵变量可以直接作为GMT的输入数据[5]。上述特性将MATLAB数学运算与GMT成果绘制的功能优势有效地融合在一起，进一步促进了GMT在地球科学研究中的应用。

目前，辽宁省地形图通常使用MapInfo和ArcGIS等软件进行绘制。其中，MapInfo的含义是“Mapping+Information（地图+信息）”即：地图对象+属性数据。MapInfo是一种操作简便的桌面地图信息系统，它具有图形的输入与编辑、图形的查询与显示、数据库操作、空间分析和图形的输出等基本功能，通过菜单条命令或工具条按钮进行操作，提供地图窗口、浏览窗口、统计窗口以及帮助输出设计的布局窗口等查看表窗口，并将结果方便地输出到打印机或绘图仪。总体来说，MapInfo软件更加侧重于地图信息的属性查询，在地图绘制功能方面相对较弱。ArcGIS for Desktop是对地理信息进行编辑、创建以及分析的一种GIS软件，提供了一系列的工具用于数据采集和管理、可视化、空间建模分析以及高级制图。ArcGIS不仅支持单用户和多用户的编辑，还可以进行复杂的自动化工作流程，同时便捷地实现栅格图像数字化输入、地图的编辑与制图输出、数据的编辑与管理，以及空间分析等功能，是专题图绘制的较好选择[6]。虽然MapInfo和ArcGIS软件操作简单，使用方便，但是在地形绘制和等值线绘制方面存在一定的局限性。以辽宁省地形图绘制为例，辽宁省基本地貌可概括为“六山一水三分田”，地势大致呈现自北向南、由东西两侧向中部倾斜的趋势，山地丘陵分列辽东、辽西两侧，辽河平原位于中部，呈马蹄形指向辽西渤海沿岸狭长的海滨平原[7]，境内山脉分别位于东西两侧。考虑到辽宁省内平原与山地丘陵地区比例均衡，且境内山地丘陵海拔相对较低，与平原地区海拔差较小，利用MapInfo和ArcGIS软件等现有手段绘制的二维地形图存在山体阴影效果不明显的现象，且难以实现三维地形图的绘制。为了更好地展现辽宁省地形地貌特征，本文拟采用GMT软件实现辽宁省二维及三维地形图的绘制。

1 GMT数据类型介绍

GMT作为一种绘图工具，同时作为一种数据处理工具，模块涉及多种输入和输出数据，但概括来讲，GMT软件的输入输出数据类型可划分为以下几类。

第一类为表数据，英文称为table data，也称为列数据或者多列数据，常用于表示点和线。表数据中有N个记录，每个记录都有M个字段。表数据可以有三种形式：ASCII表、二进制表和netCDF表。最常用的表数据形式是ASCII表，可以使用编辑器直接编辑。ASCII表中有N行M列，每行代表一个记录，每列代表一个字段。一个记录内的字段之间可以用空格、制表符、逗号或分号分隔。

第二类为网格文件，GMT可以绘制2D网格数据。通常，2D网格文件的X方向代表经度、Y方向代表纬度，Z值可以表示高程、重力值、温度、速度等。也可以将XYZ格式的1D表数据网格化得到2D网格数据。GMT可以使用的2D网格数据有三种格式：netCDF格式、Sun光栅文件和自定义的二进制数据格式。其中最常见也较为推荐的网格数据格式是netCDF格式，一般以.nc或.grd作为文件后缀。

第三类为CPT文件，CPT全称是color palette table，也称为调色板或色标文件。CPT文件可以在grdimage、psxy、psxyz等命令中使用。一般情况下，不同的数据范围颜色要求也不尽相同，所以可以使用makecpt或grd2cpt命令对已有的CPT（例如GMT内置的CPT）文件进行自定义编辑。当然，也可以根据不同程序中的颜色配比需求手写CPT文件，或者使用awk、perl之类的文本处理工具自动生成CPT文件。

第四类为PS文件，GMT生成的图片为PS文件格式，全称是PostScript。PostScript是一种用于描述矢量图形的页面描述语言。PS是矢量图形格式，即用点、线段或多边形等数学方程的几何元素来表示图像。因而可以任意旋转与缩放而不会出现图像失真的情况。GMT提供了psconvert命令，可以很方便地将PS文件以任意精度转换为jpeg、eps、png、pdf等图片格式，以满足不同情形下的需求。

2 地形图的绘制

GMT绘图命令较为复杂，从该软件网站可以获取较为详尽的使用手册和模块手册介绍，本文以辽宁省及周边作为绘图区域，利用手册总结出二维和三维地形图的绘制方法。

2.1 2D地形图的绘制

利用GMT绘制2D地形图的步骤概括如下：

（1）绘制地形底图，加载1′×1′格网分辨率的earth_relief_01m.grd全球地形模型文件。

gmt grdcut earth_relief_01m.grd-R＄R-G%grd%

gmt grdgradient%grd%-Ne0.7-A50-G%grd%

grdcut命令指从大区域地形中裁剪出自定义范围的数据（本例中经度为118°～126°E，纬度38°～44°N），该操作可降低后续梯度计算的复杂度，提高绘图效率。grdgradient命令为制作梯度信息的网格文件，得到某一个方向的方向导数。

（2）生成并编辑调试CPT文件。

gmt grd2cpt%grd%-Cglobe-S-10000/10000/200-Z-D＞%cpt%

grd2cpt命令为制作颜色表文件。-C选项指定生产的cpt文件类型，-S选项设置颜色表最小值和最大值和间隔，-Z选项生产连续颜色的CPT文件。

（3）绘制2D地形图。

gmt grdimage%grd%-R-J-I%grd%-C%cpt%-Q-O-K ＞＞＄PS

grdimage命令为在地图上绘制网格数据文件，-I选项接grdgradient计算得到的梯度文件，-C选项接颜色表文件。

（4）加载行政区边界等信息。

gmt psxy%CN-border-La% -J-R＄R-W0.5p-O-K ＞＞＄PS

gmt psxy%prov%-J-R＄R-W1.5p/red-O-K ＞＞＄PS

gmt psxy%city%-J＄J-R＄R-Sc0.08c-Gblack-K-O＞＞＄PS

gmt psxy%city%-J＄J-R＄R-Sc0.15c-W0.2p,black-K-O ＞＞＄PS

gmt pstext%city%-J＄J-R＄R-Dj0.15c/0.15c-P-K-O ＞＞＄PS

psxy为绘制线段，符号命令%CN-border-La%，%prov%，%city%分别为国界、省界、市政府所在地，-W选项设置线段或符号轮廓的画笔属性，-G选项设置符号或多边形的填充色，-S选项表示要绘制符号类型。pstext在图上写文本，-D选项表示文本在指定坐标的基础上相对x轴和y轴的偏移。pstext命令中加载的市政府所在地标注信息，其数据格式是经度、纬度、字符大小、旋转角度、字体、文本。

（5）编辑图像的标题、图例等信息。

gmt psscale-DjCB+w18c/0.3c+o0/-1.7c+h-R-J-C%cpt%-BWSEN-Bxa100f50-G0/400-O-K＞＞＄PS

psscale命令为绘制色标，-D选项表示色标在图中的位置，+w18c/0.3c表示色标长18厘米，宽0.3厘米，-C选项接颜色表文件，-B选项表示色标间隔，-G选项表示截断范围为0～400 m。

（6）将绘制的地形图输出为可移植网络图形格式（PNG）文件。

gmt psconvert＄PS-A-Tg-P-Z

-A选项表示对PS文件进行裁剪，仅保留其中有绘图的部分。-T选项指定要转换的图片格式。-Z选项转换完成后删除输入的PS文件。

基于上述步骤，利用GMT软件绘制的辽宁省2D地形图如图1（左）所示。

作为对比，利用ArcGIS软件绘制2D地形图的步骤如下：

（1）使用ArcToolbox→Spatial Analyst工具→表面分析→山体阴影工具，输入dem图层，生成新的栅格数据图层；

（2）将栅格数据图层透明度设置为35%；

（3）将栅格数据图层移至dem图层上方，选择dem图层中符号系统选项中适合的色带。

根据上述步骤，利用ArcGIS软件绘制的辽宁省2D地形图如图1（右）所示。

图1 基于GMT（左）和ArcGIS（右）绘制的辽宁省2D地形图Fig.1 The 2D topographic map of Liaoning Province by using GMT （left） and ArcGIS （right） respectively

通过根据上述命令制作的地形图，可以大致上了解辽宁地形地貌，但从图1（左）中可以看出在地势比较平缓的中部地区颜色的区分差别不大，东西两侧山体阴影效果不是十分显著。通过左右两幅图对比，利用ArcGIS绘制的地形图山体阴影效果相对较弱，仅仅在颜色上面区分高程较为明显。

业务部门在使用ArcGIS软件绘制专题图时通常只绘制二维图形，在三维展示模块具有一定的局限性。因此GMT的灵活、开源等优势更加得到凸显，下面将进一步讨论GMT的三维成图效果，并与图1进行对比。

2.2 3D地形图的绘制

利用GMT绘制3D地形图的步骤概括如下：

（1）网格化数据，生成grd文件。

surface%dat%-Gtmp.grd-R118/126/38/44-I2m

surface命令是数据网格化的一种方法。%dat%选项为网格数据，-G选项为网格化后的数据文件名，-R选项为网格化经纬度范围，-I选项为指定x方向和y方向的增量。

（2）生成并编辑调试CPT文件。

grd2cpt tmp.grd-Cglobe-S-10000/10000/200-Z-V＞%cpt%

grd2cpt命令制作颜色表文件。同二维地形图。

（3）制作梯度信息的网格文件。

grdgradient tmp.grd-A270-Gtmp.grad grdhisteq tmp.grad-Gtmp.his-N grdmath tmp.his 0.4 x=tmp.int

grdgradient命令是制作梯度信息的网格文件，得到某一个方向的方向导数，输入文件为tmp.grd，-A选项设置方向导数的方向为270°，输出的文件为tmp.grad。grdhisteq命令表示网格文件柱状图均化。grdmath命令是对grd文件进行数字计算。

（4）绘制3D地形图。

grdview-K-O-S-P tmp.grd-JM10-JZ2.5-Ba1/a1/800:"Altitude(m)":SEwnZ-Qs-N-1/235-Itmp.int-X2-C%cpt%-R-E150/45-Y13＞＞＄psf

grdcontour tmp.grd-J-P-O-K-E230/30-RC100-A1000+p+g255+s5-Wa3/0-X-1-Y-11.5＞＞＄psf

grdview-O-K-S-P tmp.grd-J-JZ2.5i-Ba1/a1/800SeWnz-Qs-N-100-Itmp.int-C%cpt%-RE230/30＞＞＄psf

grdview命令为基于二维网格数据文件绘制三维灰度或彩色图像。-B选项设置地图的边框注记以及刻度间隔，-Qs选项表示采用表面绘制，-I选项表示给出一个含有照明强度的文件名，-E选项设置观察点的方位角和高度角。

grdcontour命令为绘制等值线。-C选项表示指定等值线绘制方式，-A选项表示标注标签信息，-W选项表示指定画笔属性。

（5）编辑图例等信息。

psscale-E-D15/5/7/0.3-C%cpt%-Ba1000g 500/:"(m)":-O-P-V＞＞＄psf

（6）将绘制的地形图输出为可移植网络图形格式（JPEG）文件。

ps2raster-A-Tj＄psf

通过上述步骤，绘制辽宁省3D地形图，如图2所示。

图2 辽宁省3D地形图Fig.2 3D topographic map of Liaoning Province

图2上图为俯视的辽宁省三维地形图，下图为含有等高线的三维地形图。图2相比于图1来说，在展现山体阴影方面有了很大的提高，同时在中部地区的平地也可以看出地形的凹凸起伏。相比二维地形图，三维地形图在地形地貌的表达上具有较大的优势。

3 结论

通过上述例子可以看出，GMT可以绘制出精美直观的二维和三维的地形图，并且在此基础上可添加行政边界和政府所在地等信息，结果表明GMT在绘制地形图上具有操作方便、图形美观等优势。此外，GMT软件同样实现了地形图由平面到立体，由抽象到直观，打破了二维平面展示和二维可视的限制，更好的展现了辽宁省地形，为后续基于地形图绘制地震应急专题图等方面提供积极性的参考。