zhikuiguo的个人博客分享 http://blog.sciencenet.cn/u/zhikuiguo

博文

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

已有 6613 次阅读 2018-11-6 01:09 |个人分类:读博那点事儿|系统分类:科研笔记| GMT三维切片图

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

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

提出问题

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

Olive, et al., 2015, Science

Olive, et al., 2015, Science

Schmid F. and Schlindwein V., 2016, G3

Schmid F. and Schlindwein V., 2016, G3

Zhao, et al., 2013, 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. 仰角 过视角点做一个与此图平面P1垂直的面P2。在P2面内,视角点与参考点的直线假设为L0,那么L0P1面之间的夹角就是仰角挺复杂,姑且这么理解,慢慢体会)。所以,如果仰角=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°

-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) 不用再用psconvertps格式转换为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)中的玫瑰方向标的变化。 \"不同视角下的-JZ参数绘图实例\"/

  • 代码

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绘制的两个切片:平面模式

\"利用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位置:显示位置错乱

两个切片分别位于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文件里面。这里面的计算公式是有问题的,我已将其修改,具体修改方法见下回分解.

修正后

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


https://wap.sciencenet.cn/blog-621290-1144740.html

上一篇:地球上不止有谷歌学术
下一篇:Ubuntu下源码编译安装GMT
收藏 IP: 134.245.223.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-29 04:26

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部