相信理工科方面的科研人员经常会遇到做三维切片图的时候，以地球物理为例，比如重力、磁异常三维反演，地震速度剖面+地形图，热液循环三维数值模拟（CFD范畴），地球动力学三维模拟等。当然了做三维切片图也有商业软件(比如Tecplot 360)和开源软件（比如Paraview），前者缺点是收费收费收费，后者导出的三维场景图效果很差（不是矢量图，colorbar位置错乱等问题），其优点是做现场演示。还有我最喜欢的编程语言python也可以做这个工作，但是我觉得它对坐标轴的label处理的不是那么完美，但也是个不错的选择。我的原则就是用开源的程序或者软件做出最漂亮最简洁并且高质量(本人追求矢量的eps、pdf、ps格式)的学术论文图。经过我四年的不断折腾和筛选，发现到目前为止，GMT 6.0最新版（我是直接从github克隆下来的开发版）对三维图支持已经非常棒了。不仅可以做三维地形图（上面标注文字，数据点，线，legend等），还可以做切片的等值线图或者image。

一贯的风格：先来一段回顾和概述，然后提出问题，再解决问题，最后点评一下存在的问题或可选方案

提出问题

这种图在什么情况下使用？直接上图感受一下

Olive, et al., 2015, Science

Schmid F. and Schlindwein V., 2016, G3

Zhao, et al., 2013, G3

GMT如何实现三维切片图

gmt新版推出的三维成图主要通过-p和-JZ实现的。前者进行旋转，后者给定垂向轴相关信息。gmt官网上有个用户给出的例子，但是没有讲清楚这个切片图到底如何准确实现，但是还是有借鉴意义的。因为本人的博士课题有部分数据涉密，不方面放上来，就以上面提到的这个例子里面的函数造一个数据进行成图。

-p参数

这个参数是新版gmt里面一个很通用的参数，在basemap, grdimage, grdcontour, grdview等都有这个参数，下面以basemap为例介绍-p参数的用法。具体请参见官网英文解释

• 参数解释

-p后面可以跟x,y,z（这三个字母表示你的此图的纵轴是哪个），然后再跟方位角(azimuth)/仰角(elevation)，或者直接跟方位角/仰角。方位角的范围是(0, 360]，仰角范围是(0, 90]，这个很重要，比如方位角或者仰角为0，则gmt会报错。 1. 方位角 表示的是绕纵轴(由-p后面跟的字母决定，可以是x,y,z如果没给出则默认为z)旋转多少度，方位角=180的情况相当于没有旋转，小于180°表示顺时针旋转；大于180°表示顺时针旋转。 2. 仰角 过视角点做一个与此图平面P₁垂直的面P₂。在P₂面内，视角点与参考点的直线假设为L₀，那么L₀与P₁面之间的夹角就是仰角 （挺复杂，姑且这么理解，慢慢体会）。所以，如果仰角=90表示俯视图，也就是下图(b)所示的只是在(a)图的基础上绕z轴顺时针旋转了180^° − 135^° = 45^°；(c)图在(b)图的基础上上仰了45°。 3. 切片位置 -p后面除了方位角和仰角，还可以再跟第三个数字，类似这样的格式:-p方位角/仰角/在纵轴的位置 即-pazimuth/elevation/pos，这个pos设定了此图在纵轴的什么位置处。比如-pz135/60/400表示此图平面是x-y平面（与纸张面平行），位于纵轴(也就是z轴）的z = 400处；-px135/60/-50表示此图的平面与y-z平面平行，位于x = −50处，这是设置切片位置的关键。这个功能必须与-JZ参数配个才生效，同时要求-R后面必须跟六个参数表示坐标范围（见下面的-JZ解释）。

注意：下面的所有操作都是在Mac系统下进行的，而且用到了awk程序（配合gmt做一些文字处理和简单计算），如何安装和使用gmt开发版和如何配置gmt请阅读gmt第一课

• 绘图实例

-p参数在底图中的使用。（a）没有使用-p参数；（b）方位角135°/仰角90°；（c）仰角45°

• 代码

fig_name=test_basemap_pfig_fmt=pnggmt begin $fig_name $fig_fmt    gmt gmtset MAP_FRAME_TYPE=fancy+    fig_width=14    lon_min=-75    lon_max=-60    lat_min=-50    lat_max=-40    # 方位角    azimuth=135    # 仰角    elevation=90    range_lon_lat=$lon_min/$lon_max/$lat_min/$lat_max    gmt basemap -JM${fig_width}c -R$range_lon_lat -Ba5f1 -BWSen+t"(a) without -p" -TdjBL+o0c/0c+w1.5c+f+l,E,,N 
    gmt basemap -B -BWSen+t"(b) -p"$azimuth"/"$elevation -TdjBL+o0c/0c+w1.5c+f+l,E,,N -X16c -p${azimuth}/${elevation}    # 更改仰角    elevation=45    gmt basemap -B -BWSen+t"(c) -p"$azimuth"/"$elevation -TdjBL+o0c/0c+w1.5c+f+l,E,,N -X19c -p${azimuth}/${elevation}gmt end# 打开文件查看open ${fig_name}.${fig_fmt}

• 关于代码

这个代码风格是新版gmt推出的modern模式，带来了很大方便。 (1) 不用每个绘图命令后面都要跟一个-O,-K,>>xxx.ps等 (2) 不用再用psconvert将ps格式转换为png,pdf,svg了，直接在开头给出文件名和格式即可 (3) MAP_FRAME_TYPE=fancy+ fancy表示坐标轴显示为火车道样式；+表示拐角显示为圆角(不信你就上翻在看一眼图哟) (4) 如果投影参数-J, 范围参数-R与前一个绘图完全一致，则可以省略这两个参数。但是-B不能省略，不然将不显示坐标轴，虽然坐标轴设置参数也完全一致

-JZ参数

上面讲的-p参数只是对图做了方位旋转，还不是真正的三维作图，因为所有的绘图元素都在一个平面内。要实现三维成图，需要附加-JZ参数，一看-J就知道肯定是与投影有关了，此参数设定了纵轴方向的长度。更多-J参数参见官网说明

• 参数解释

1. 图像纵轴高度 与-JX参数类似，后面直接跟此轴方向的图像大小，单位可以是c厘米或者i英尺。

2. 坐标范围 gmt绘图必须给定坐标范围即-R参数（当然了这是针对第一个绘图命令而言的），除了新版(6.0版)的-Re和-Ra之外，对于二维绘图必须给出四个数字(-Rxmin/xmax/ymin/ymax)，而对于三维绘图（也就是有-JZ参数出现的时候）必须给出六个数字(-Rxmin/xmax/ymin/ymax/zmin/zmax)

3. 坐标轴 坐标轴的label格式用类似命令-Bza250f50g250(主刻度间隔250，副刻度间隔50，显示网格间隔250)设置；如果想显示纵轴，则用-BWseNZ类似的命令设置，大写字母表示要显示哪个轴，如果要显示立方体框(cubic outline)，只要在后面跟一个+b即可(-BWseNZ+b)。

• 绘图实例

体会不同旋转角度的视角差异，注意（a）和（b）中的玫瑰方向标的变化。

• 代码

fig_name=test_basemap_JZfig_fmt=pnggmt begin $fig_name $fig_fmt    gmt gmtset MAP_FRAME_TYPE=fancy+    fig_width=14    # 增加图像高度：及纵轴方向的长度    fig_height=8    lon_min=-75    lon_max=-60    lat_min=-50    lat_max=-40    # 增加纵轴方向的范围    z_min=0    z_max=999    # 方位角    azimuth=135    # 仰角    elevation=45    range_lon_lat=$lon_min/$lon_max/$lat_min/$lat_max    range_z=$z_min/$z_max    range=$range_lon_lat/$range_z    gmt basemap -R$range -JM$fig_width -JZ$fig_height -Ba5f1 -Bza250f50g250+l"km" -BWSenZ+b+t"3D basemap: "$azimuth"/"$elevation -TdjBL+o0c/0c+w1.5c+f+l,E,,N -pz${azimuth}/${elevation} 
    azimuth=45    gmt basemap -JZ -B -Bza250f50g250+l"km" -BWSenZ+b+t"3D basemap: "$azimuth"/"$elevation -TdjBL+o0c/0c+w1.5c+f+l,E,,N -pz${azimuth}/${elevation} -X19c    elevation=30    gmt basemap -JZ -B -Bza250f50g250+l"km" -BWSenZ+b+t"3D basemap: "$azimuth"/"$elevation -TdjBL+o0c/0c+w1.5c+f+l,E,,N -pz${azimuth}/${elevation} -X19cgmt endopen ${fig_name}.${fig_fmt}

• 关于代码

（1）如果使用了-JZ参数，必须设定六位数的-R参数 （2）-J参数可以省略（继承上一个绘图命令的-J参数），但是-JZ不能省略 （3）-Bz也不能省略

三维切片绘图实例

利用GMT绘制的两个切片:平面模式

利用GMT绘制的两个切片：三维立体模式

#!/bin/bashgmt set GMT_COMPATIBILITY=5 MAP_FRAME_TYPE=plain# 1. 数据和绘图坐标参数设置fig_width_x=14fig_width_z=8# 数据范围lon_min=-75lon_max=-60lat_min=-50lat_max=-40zmin=0zmax=999# 数据中心点坐标lon0=`echo $lon_min $lon_max | awk '{print $1+($2-$1)/2.0}'`lat0=`echo $lat_min $lat_max | awk '{print $1+($2-$1)/2.0}'`echo $lon0 $lat0range_lon_lat=$lon_min/$lon_max/$lat_min/$lat_maxrange_z=$zmin/$zmaxrange=$range_lon_lat/$range_z# 投影参数proj_xy=M$lon0/$lat0/$fig_width_x  #x-y平面的投影参数# 数据文件名data_xy=base.ncdata_xy_cart=base_cart.ncdata_yz=slice_cut.nc# 2. x-y平面内的网格数据# gmt grdmath -R$range_lon_lat -I0.005 X D2R Y D2R ADD STO@xySum SIN @xySum 3 MUL NEG EXP MUL = $data_xy# 3. 沿纬度方向的垂直切片lat_min_slice=-47.5lat_max_slice=-42.5# gmt grdmath -R${lat_min_slice}/${lat_max_slice}/${zmin}/${zmax} -I0.005/0.5 X D2R -67.5 D2R ADD STO@xySum SIN @xySum 3 Y 1E4 DIV SUB MUL NEG EXP MUL = slice.nc# 4. 将数据范围投影到图像空间的大小:x对应经度方向；y对应纬度方向xmin=`echo $lon_min $lat_min | gmt mapproject -R$range_lon_lat -J$proj_xy | awk '{print $1}'`ymin=`echo $lon_min $lat_min | gmt mapproject -R$range_lon_lat -J$proj_xy | awk '{print $2}'`xmax=`echo $lon_max $lat_max | gmt mapproject -R$range_lon_lat -J$proj_xy | awk '{print $1}'`ymax=`echo $lon_max $lat_max | gmt mapproject -R$range_lon_lat -J$proj_xy | awk '{print $2}'`ymin_slice=`echo $lon0 $lat_min_slice | gmt mapproject -R$range_lon_lat -J$proj_xy | awk '{print $2}'`ymax_slice=`echo $lon0 $lat_max_slice | gmt mapproject -R$range_lon_lat -J$proj_xy | awk '{print $2}'`echo $xmin $xmax $ymin $ymax $ymin_slice $ymax_slicefig_width_y=$ymax# 5. 将切片数据投影到图像大小的空间：最好利用中心经度，也就是后面的-JM参数中的中心经度值# mapproject输入经纬度值（读取前两列），投影参数与后面作图的时候一致，输出的前两列是x,y值，位于[0,fig_width_x]和[0,fig_width_y]范围内,fig_width_y目前还不知道，需要计算# gmt grd2xyz slice.nc | awk '{print "'$lon0'",$0}' | gmt mapproject -R$range_lon_lat -J$proj_xy |awk '{print $2,$3,$4}'> points.txt# gmt grdproject $data_xy -R$range_lon_lat -J$proj_xy -r -E300 -G$data_xy_cart# 6. 重新网格化切片# gmt surface points.txt -G$data_yz -R$ymin_slice/$ymax_slice/$zmin/$zmax -I1500+/2000+ -C0.1 -T0.25# # 开始绘图proj_yz=X$fig_width_y/$fig_width_zproj_xz=X$fig_width_x/$fig_width_zrm gmt.history gmt.conf# 8.1 将数据显示在平面内fig_name=vslice_gmtfig_fmt=pnggmt begin $fig_name $fig_fmt    # x-y平面的数据显示    gmt basemap -J$proj_xy -R$range_lon_lat -Ba5f1 -BWSen+t"Data on x-y plane" 
    gmt grdimage $data_xy -Ccpt-city/oc/sst    gmt coast -Df -A0/0/1 -N1/0.5p,black,-..- -W0.5p,black    # y-z切片数据显示    move_x=16    gmt basemap -J$proj_yz -R$ymin/$ymax/$zmin/$zmax -Bxa5f1+l"y (latitude)" -Bya200f40+l"z (km)" -BWSen+t"Slice data on y-z plane" -X$move_x -p90/90    gmt grdimage $data_yz -Crainbow -p    # x-z切片数据显示    gmt basemap -J$proj_yz -R$xmin/$xmax/$zmin/$zmax -Bxa5f1+l"x (longitude)" -Bya200f40+l"z (km)" -BWSen+t"Slice data on x-z plane" -X-$move_x -Y-10 
    gmt grdimage $data_yz -Chot 
gmt endopen ${fig_name}.${fig_fmt}# 8.2 显示为3维切片形式# 方位角设置azimuth=135elevation=45angle_view=$azimuth/$elevationfig_name=vslice_gmt_3Dfig_fmt=pnggmt begin $fig_name $fig_fmt    echo "figwidth "$fig_width_x $fig_width_y $fig_width_z    echo "range " $xmin $xmax $ymin $ymax $zmin $zmax    # 三维框架    gmt basemap -J$proj_xy -JZ$fig_width_z -R$range -pz$angle_view -Ba5f1 -Bza200f40g200+l"Z (km)" -BwSEnZ+b 
    # 贴x-y平面数据图    gmt grdimage $data_xy -J$proj_xy -JZ$fig_width_z -R$range -pz$angle_view/0 -Ccpt-city/oc/sst --MAP_FRAME_AXES=''    gmt coast -J$proj_xy -JZ$fig_width_z -R$range -pz$angle_view/0 -Df -A0/0/1 -N1/0.5p,black,-..- -W0.5p,black 
    # # 贴y-z平面的切片数据    gmt grdimage $data_yz -JX$fig_width_y/$fig_width_z -JZ$fig_width_x -R$ymin/$ymax/$zmin/$zmax/$xmin/$xmax  -px$angle_view/7 -Crainbow 
    gmt grdimage $data_yz -JX$fig_width_x/$fig_width_z -JZ$fig_width_y -R$xmin/$xmax/$zmin/$zmax/$ymin/$ymax -py$angle_view/0 -Chot 
gmt endopen ${fig_name}.${fig_fmt}

存在BUG

错误展示

旋转之后切片位置不能很好的归位，因为源代码里面有问题，比如下图

两个切片分别位于xmin和ymin位置：显示位置错乱

查找问题

找到三维透视图绘图的源代码gmt_plot.c->gmt_plane_perspective函数（位于7429行）

double a, b, c, d, e, f;struct PSL_CTRL *PSL= GMT->PSL;
....//省略一些代码....switch (plane % 3) {            case GMT_X: /* Constant x, Convert y,z to x',y' */
                a = GMT->current.proj.z_project.sin_az;
                b = -GMT->current.proj.z_project.cos_az * GMT->current.proj.z_project.sin_el;
                c = 0.0;
                d = GMT->current.proj.z_project.cos_el;
                e = GMT->current.proj.z_project.x_off - level * GMT->current.proj.z_project.cos_az;
                f = GMT->current.proj.z_project.y_off - level * GMT->current.proj.z_project.sin_az * GMT->current.proj.z_project.sin_el;
。。。//省略PSL_command (PSL, "%s [%g %g %g %g %g %g] concat\n",
            (GMT->current.proj.z_project.plane >= 0) ? "PSL_GPP setmatrix" : "/PSL_GPP matrix currentmatrix def",
            a, b, c, d, e * PSL->internal.x2ix, f * PSL->internal.y2iy);

根据坐标轴范围(-R参数)和方位角,仰角,level值（也即是-p后面的第三个参数）计算a,b,c,d,e,f，然后将这六个数字写入ps文件里面。这里面的计算公式是有问题的，我已将其修改，具体修改方法见下回分解.

修正后

将e和f重新计算后，得到了正确的结果，如下图

上面两个图的绘图代码均为

gmt begin $fig_name $fig_fmt    echo "figwidth "$fig_width_x $fig_width_y $fig_width_z    echo "range " $xmin $xmax $ymin $ymax $zmin $zmax    # 三维框架    gmt basemap -J$proj_xy -JZ$fig_width_z -R$range -pz$angle_view -Ba5f1 -Bza200f40+l"Z (km)" -BwSEnZ+b+t"After debug" 
    # 贴x-y平面数据图    gmt grdimage $data_xy -J$proj_xy -JZ$fig_width_z -R$range -pz$angle_view/0 -Ccpt-city/oc/sst --MAP_FRAME_AXES=''    gmt coast -J$proj_xy -JZ$fig_width_z -R$range -pz$angle_view/0 -Df -A0/0/1 -N1/0.5p,black,-..- -W0.5p,black 
    # 贴y-z平面的切片数据    gmt grdimage $data_yz -JX$fig_width_y/$fig_width_z -JZ$fig_width_x -R$ymin/$ymax/$zmin/$zmax/$xmin/$xmax  -px$angle_view/0 -Crainbow 
    # X-Z平面    gmt grdimage $data_yz -JX$fig_width_x/$fig_width_z -JZ$fig_width_y -R$xmin/$xmax/$zmin/$zmax/$ymin/$ymax -py$angle_view/0 -Chot 
gmt end