Modern GMT Series：Slice in 3D View (三维切片图)

2022-04-14 00:59:23

Fancy的版本见九天学者的个人博客，关注文集博士干点啥或者微信公众号九天学者及时获取连载更新。

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

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

提出问题

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

Modern GMT Series：Slice in 3D View (三维切片图)

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)，或者直接跟方位角/仰角。方位角的范围是 $Modern GMT Series：Slice in 3D View (三维切片图)$ ，仰角范围是 $Modern GMT Series：Slice in 3D View (三维切片图)$ ，这个很重要，比如方位角或者仰角为0，则gmt会报错。

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

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

绘图实例

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

代码

fig_name=test_basemap_p
fig_fmt=png
gmt 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参数参见官网说明

参数解释

图像纵轴高度 与-JX参数类似，后面直接跟此轴方向的图像大小，单位可以是c厘米或者i英尺。
坐标范围 gmt绘图必须给定坐标范围即-R参数（当然了这是针对第一个绘图命令而言的），除了新版(6.0版)的-Re和-Ra之外，对于二维绘图必须给出四个数字(-Rxmin/xmax/ymin/ymax)，而对于三维绘图（也就是有-JZ参数出现的时候）必须给出六个数字(-Rxmin/xmax/ymin/ymax/zmin/zmax)
坐标轴 坐标轴的label格式用类似命令-Bza250f50g250(主刻度间隔250，副刻度间隔50，显示网格间隔250)设置；如果想显示纵轴，则用-BWseNZ类似的命令设置，大写字母表示要显示哪个轴，如果要显示立方体框(cubic outline)，只要在后面跟一个+b即可(-BWseNZ+b)。

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

不同视角下的-JZ参数绘图实例

代码

fig_name=test_basemap_JZ
fig_fmt=png
gmt 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} -X19c
gmt end

open ${fig_name}.${fig_fmt}

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

三维切片绘图实例

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

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

#!/bin/bash
gmt set GMT_COMPATIBILITY=5 MAP_FRAME_TYPE=plain
# 1. 数据和绘图坐标参数设置
fig_width_x=14
fig_width_z=8
# 数据范围
lon_min=-75
lon_max=-60
lat_min=-50
lat_max=-40
zmin=0
zmax=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 $lat0
range_lon_lat=$lon_min/$lon_max/$lat_min/$lat_max
range_z=$zmin/$zmax
range=$range_lon_lat/$range_z
# 投影参数
proj_xy=M$lon0/$lat0/$fig_width_x  #x-y平面的投影参数
# 数据文件名
data_xy=base.nc
data_xy_cart=base_cart.nc
data_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.5
lat_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_slice
fig_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_z
proj_xz=X$fig_width_x/$fig_width_z
rm gmt.history gmt.conf
# 8.1 将数据显示在平面内
fig_name=vslice_gmt
fig_fmt=png
gmt 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 end

open ${fig_name}.${fig_fmt}

# 8.2 显示为3维切片形式
# 方位角设置
azimuth=135
elevation=45
angle_view=$azimuth/$elevation
fig_name=vslice_gmt_3D
fig_fmt=png
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 -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 end

open ${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重新计算后，得到了正确的结果，如下图

Debug之后的结果，位置正确

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

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

如果此文对您有启发，感谢您支持原创: 有钱的动手打个赏，没钱的动手点个赞
您的鼓励就是我原创的动力。个人水平有限，若有问题，可在下方留言讨论。

码农公寓