前期一篇文章:Modern GMT Series：Slice in 3D View (三維切片圖)中提到了在科研作圖中會(huì)經(jīng)常遇到三維作圖的問題源譬，而且GMT可以做三維圖且導(dǎo)出的圖片質(zhì)量非常高集惋！但是也提到了當(dāng)前GMT版本中的三維繪圖存在漏洞。主要兩個(gè)問題：（1）切片位置錯(cuò)亂踩娘；（2）三維文字鏡像刮刑。就這兩個(gè)問題極大的限制了GMT三維作圖的應(yīng)用。因?yàn)樽罱撐睦锩嬉鋈S立體圖(地理坐標(biāo)+笛卡爾坐標(biāo)+海底地形+剖面切片等)养渴，找了找也沒找到更好的解決方案雷绢，干脆自己動(dòng)手豐衣足食——修改GMT的源代碼。經(jīng)過一天半的努力理卑，已經(jīng)比較完美的修復(fù)了這兩個(gè)漏洞翘紊。目前官方版本中并沒有修復(fù)此bug。 下面就介紹一下修復(fù)漏洞的過程和最終效果傻工！
先上一張BUG修復(fù)后的效果圖鑒賞一下
注： 由于我的論文正在審稿中霞溪，不便上傳較為復(fù)雜的繪圖結(jié)果，暫用一個(gè)相對(duì)簡單的圖來代替中捆。
簡單例子
問題描述：繪制一個(gè)三維立方體并在每個(gè)面上標(biāo)注文字以及設(shè)置每個(gè)面具有不同顏色且半透明鸯匹；水平旋轉(zhuǎn)45度，然后x軸和y軸標(biāo)注坐標(biāo)標(biāo)簽泄伪。
GMT官方版本運(yùn)行的結(jié)果
繪圖代碼
# Test 3D View of GMT: simple example
# Zhikui Guo, 2018-12-09, GEOMAR, Germany
function preset()
{
figname=Test_3d
angle_view=-45/25
width_fig_x=10
width_fig_y=10
width_fig_z=7
# 透明度
alpha_profile=50
# range
xmin=0
xmax=100
ymin=0
ymax=50
zmin=-20
zmax=0
xc=`echo $xmin $xmax | awk '{print $1+($2-$1)/2}'`
yc=`echo $ymin $ymax | awk '{print $1+($2-$1)/2}'`
zc=`echo $zmin $zmax | awk '{print $1+($2-$1)/2}'`
len_z=`echo $zmin $zmax | awk '{print ($2-$1)}'`
}
# 1. apply theme
# . styles.sh
# MonokaiTheme
# 常用代碼頭文件:繪制logo
# . stdafx.sh
# 范圍和顏色等設(shè)置
preset
# plot
gmt begin $figname pdf,png
gmt basemap -JX$width_fig_x/$width_fig_y -JZ$width_fig_z -R$xmin/$xmax/$ymin/$ymax/$zmin/$zmax -Ba -Bza+l"Z(m)" -BwsenZ+gred -pz$angle_view/$zmin
echo "$xc $yc zmin plane" | gmt pstext -JZ -p -F+f20p,Helvetica-Bold,blue=thinner,white+jCM+a190 -Dj0c/0c
gmt basemap -JX$width_fig_x/$width_fig_y -JZ$width_fig_z -R$xmin/$xmax/$ymin/$ymax/$zmin/$zmax -Bwsenz+ggray -pz$angle_view/$zmax
echo "$xc $yc zmax plane" | gmt pstext -JZ -p -F+f20p,Helvetica-Bold,blue=thinner,white+jCM+a190 -Dj0c/0c
gmt basemap -JX$width_fig_y/$width_fig_z -JZ$width_fig_x -R$ymin/$ymax/$zmin/$zmax/$xmin/$xmax -Ba -BwSenz+glightblue@$alpha_profile -Bx+l"Y(m)" --MAP_FRAME_PEN=1,black@0 --MAP_ANNOT_OFFSET=-0.2 -px$angle_view/$xminn
echo "$yc $zc ymin plane" | gmt pstext -JZ -p -F+f20p,Helvetica-Bold,blue=thinner,white+jCM+a45 -Dj0c/0c
gmt basemap -JX$width_fig_x/$width_fig_z -JZ$width_fig_y -R$xmin/$xmax/$zmin/$zmax/$ymin/$ymax -Ba -BwSenz+glightgreen@$alpha_profile -Bx+l"X(m)" --MAP_FRAME_PEN=1,black@0 --MAP_ANNOT_OFFSET=-0.2 -py$angle_view/$ymax
echo "$xc $zc xminn plane" | gmt pstext -JZ -p -F+f20p,Helvetica-Bold,blue=thinner,white+jCM+a-45 -Dj0c/0c
# 使用logo
# add_logo 0 0.5 moderngmt -grid=false
gmt end
# open $figname.pdf
# rm tmp* gmt.conf gmt.history
這段代碼不涉及什么數(shù)據(jù)殴蓬，直接可以復(fù)制運(yùn)行的。這里我將主題和logo注釋了，這是我自己定義的染厅，因?yàn)槟愕碾娔X上肯定是沒有相關(guān)的文件和函數(shù)痘绎。直接復(fù)制上面的代碼應(yīng)該是可以直接出圖的。注意這是Mac系統(tǒng)下的腳本肖粮，也可以在Linux系統(tǒng)下使用孤页，但是windows系統(tǒng)下需要將變量引用從
$xxx改為%xxx%。
都有什么BUG涩馆？
從上面的第二個(gè)圖中可以看出行施，官方版本的三維繪圖結(jié)果存在一下幾個(gè)問題：
切片位置錯(cuò)亂：需要向y軸正方向平移一個(gè)量
文字出現(xiàn)鏡像：x軸和y軸的標(biāo)注文字都變成了鏡像且有一個(gè)翻轉(zhuǎn)角度
普通文字鏡像：除了坐標(biāo)軸的標(biāo)注文字外，剖面切片上的文字的角度也有錯(cuò)亂和鏡像現(xiàn)象
如何修復(fù)BUG魂那？
簡單介紹以上三個(gè)BUG是如何修復(fù)的蛾号，具體見修復(fù)后的源代碼(我的資源庫中有介紹如何獲取源代碼和可用軟件)。
切片位置錯(cuò)亂問題
GMT繪圖是基于PostScript的涯雅，三維繪圖的坐標(biāo)旋轉(zhuǎn)理論見GMT三維坐標(biāo)旋轉(zhuǎn)理論鲜结。經(jīng)測試發(fā)現(xiàn)問題出現(xiàn)在gmt_plot.c源文件里面的gmt_plane_perspective函數(shù)
a = b = c = d = e = f = 0.0;
if (plane < 0) /* Reset to original matrix */
PSL_command (PSL, "PSL_GPP setmatrix\n");
else { /* New perspective plane: compute all derivatives and use full matrix */
if (plane >= GMT_ZW) level = gmt_z_to_zz (GMT, level); /* First convert world z coordinate to projected z coordinate */
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;
break;
case GMT_Y: /* Constant y. Convert x,z to x',y' */
a = -GMT->current.proj.z_project.cos_az;
b = -GMT->current.proj.z_project.sin_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.sin_az;
f = GMT->current.proj.z_project.y_off - level * GMT->current.proj.z_project.cos_az * GMT->current.proj.z_project.sin_el;
break;
case GMT_Z: /* Constant z. Convert x,y to x',y' */
a = -GMT->current.proj.z_project.cos_az;
b = -GMT->current.proj.z_project.sin_az * GMT->current.proj.z_project.sin_el;
c = GMT->current.proj.z_project.sin_az;
d = -GMT->current.proj.z_project.cos_az * GMT->current.proj.z_project.sin_el;
e = GMT->current.proj.z_project.x_off;
f = GMT->current.proj.z_project.y_off + level * GMT->current.proj.z_project.cos_el;
break;
}
// printf("ix: %f iy: %f\n",PSL->internal.x2ix,PSL->internal.y2iy);
/* First restore the old matrix or save the old one when that was not done before */
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);
}
正如GMT三維坐標(biāo)旋轉(zhuǎn)理論所提到的公式，gmt_plane_perspective函數(shù)的這段代碼目的就是根據(jù)MGT繪圖命令中的-R, -J, -JZ, -p參數(shù)計(jì)算坐標(biāo)旋轉(zhuǎn)矩陣參數(shù)a,b,c,d,e,f活逆，最后使用PSL_command將坐標(biāo)旋轉(zhuǎn)矩陣寫入ps文件就能得到三維的效果精刷。
** GMT_Z的情況都是OK的，正如上圖所示三維坐標(biāo)軸主框架沒錯(cuò)蔗候。但是 GMT_X和 GMT_Y**的情況贬养，也就是分別于X軸和Y軸垂直的切片的情況，出現(xiàn)了問題琴庵。經(jīng)測試發(fā)現(xiàn)，問題在于e, f計(jì)算錯(cuò)誤仰美。三種情況的e, f值需要一致才不會(huì)發(fā)生未知平移錯(cuò)亂迷殿。
解決方案
這里有一個(gè)很容易的解決方案就是在case GMT_Z計(jì)算得到e,f值的時(shí)候，打開一個(gè)臨時(shí)文件將這兩個(gè)值保存到臨時(shí)文件里面咖杂，然后如果執(zhí)行到case GMT_X或者case GMT_Y的時(shí)候?qū)⑦@個(gè)臨時(shí)文件里面的e,f值讀入覆蓋掉程序計(jì)算的結(jié)果庆寺。這雖然不是一種完美的解決方案，但是可以實(shí)實(shí)在在的修復(fù)這個(gè)bug诉字。以后有時(shí)間了仔細(xì)研究GMT內(nèi)部是如何計(jì)算這些參數(shù)然后找到計(jì)算錯(cuò)誤的位置懦尝，再做完美修改方案。修改源代碼后壤圃，重新編譯得到更新的gmt程序陵霉，然后再運(yùn)行這個(gè)例子，得到的結(jié)果如下：
坐標(biāo)軸標(biāo)注文字鏡像問題
上面一步修復(fù)了剖面位置問題伍绳，剖面坐標(biāo)軸標(biāo)注文字依然存在鏡像翻轉(zhuǎn)的現(xiàn)象踊挠。PS的鏡像翻轉(zhuǎn)命令是-1 1 scale，如果相對(duì)某個(gè)文字進(jìn)行鏡像翻轉(zhuǎn)冲杀，只需要在ps文件里面找到這個(gè)文字對(duì)應(yīng)的代碼效床，然后在前面加上這個(gè)命令即可實(shí)現(xiàn)鏡像翻轉(zhuǎn)睹酌。找到這個(gè)magic命令真是不容易，參見一個(gè)論壇剩檀。
依照這個(gè)邏輯憋沿，在源代碼中跟蹤找到繪制坐標(biāo)軸label的位置：gmt_plot.c文件中的gmt_xy_axis函數(shù)，關(guān)鍵代碼如下：
/* Finally do axis label */
if (A->label[0] && annotate && !gmt_M_axis_is_geo_strict (GMT, axis)) {
unsigned int far_ = !below;
char *this_label = (far_ && A->secondary_label[0]) ? A->secondary_label : A->label; /* Get primary or secondary axis label */
if (!MM_set) PSL_command (PSL, "/MM {%s%sM} def\n", neg ? "neg " : "", (axis != GMT_X) ? "exch " : "");
form = gmt_setfont (GMT, &GMT->current.setting.font_label);
PSL_command (PSL, "/PSL_LH ");
PSL_deftextdim (PSL, "-h", GMT->current.setting.font_label.size, "M");
PSL_command (PSL, "def\n");
PSL_command (PSL, "/PSL_L_y PSL_A0_y PSL_A1_y mx %d add %sdef\n", PSL_IZ (PSL, GMT->current.setting.map_label_offset), (neg == horizontal) ? "PSL_LH add " : "");
/* Move to new anchor point */
PSL_command (PSL, "%d PSL_L_y MM\n", PSL_IZ (PSL, 0.5 * length));
if (axis == GMT_Y && A->label_mode) {
i = (below) ? PSL_MR : PSL_ML;
PSL_plottext (PSL, 0.0, 0.0, -GMT->current.setting.font_label.size, this_label, 0.0, i, form);
}
else
PSL_plottext (PSL, 0.0, 0.0, -GMT->current.setting.font_label.size, this_label, horizontal ? 0.0 : 90.0, PSL_BC, form);
}
代碼段里面調(diào)用PSL_plottext繪制坐標(biāo)軸label(也就是圖中的X(m)和Y(m))沪猴，為了方便管理和代碼重用辐啄，我在gmt_plot.c里面自定義了兩個(gè)函數(shù): AxisLable_Flip_GMT_X_Y和AxisTickLabel_Flip_GMT_X_Y分別對(duì)axis label和axis tick label進(jìn)行操作。注意：tick label除了鏡像問題還存在一個(gè)180度的旋轉(zhuǎn)角度問題字币，這些都需要根據(jù)-p參數(shù)進(jìn)行重新計(jì)算和判斷然后修復(fù)则披。
加入自定義函數(shù)后的代碼段如下：
/* Move to new anchor point */
PSL_command (PSL, "%d PSL_L_y MM\n", PSL_IZ (PSL, 0.5 * length));
AxisLable_Flip_GMT_X_Y(GMT);// X labbel: 官方版本中的對(duì)于x和y剖面切片的文字出現(xiàn)了鏡像，所以這里只對(duì)x和y剖面進(jìn)行翻轉(zhuǎn)
if (axis == GMT_Y && A->label_mode) {
i = (below) ? PSL_MR : PSL_ML;
PSL_plottext (PSL, 0.0, 0.0, -GMT->current.setting.font_label.size, this_label, 0.0, i, form);
}
else
PSL_plottext (PSL, 0.0, 0.0, -GMT->current.setting.font_label.size, this_label, horizontal ? 0.0 : 90.0, PSL_BC, form);
AxisLable_Flip_GMT_X_Y(GMT);//恢復(fù)翻轉(zhuǎn)
ticklabel 和 axis label都在同一個(gè)函數(shù)里面洗出，調(diào)用自定義函數(shù)AxisTickLabel_Flip_GMT_X_Y修復(fù)之后的代碼片段如下：
for (i = 0; i < nx1; i++) {
if (gmtlib_annot_pos (GMT, val0, val1, T, &knots[i], &t_use)) continue; /* Outside range */
if (axis == GMT_Z && fabs (knots[i] - GMT->current.proj.z_level) < GMT_CONV8_LIMIT) continue; /* Skip z annotation coinciding with z-level plane */
if (GMT->current.setting.map_frame_type & GMT_IS_INSIDE && (fabs (knots[i] - val0) < GMT_CONV8_LIMIT || fabs (knots[i] - val1) < GMT_CONV8_LIMIT)) continue; /* Skip annotation on edges when MAP_FRAME_TYPE = inside */
if (!is_interval && plot_skip_second_annot (k, knots[i], knots_p, np, primary)) continue; /* Secondary annotation skipped when coinciding with primary annotation */
x = (*xyz_fwd) (GMT, t_use); /* Convert to inches on the page */
/* Move to new anchor point */
PSL_command (PSL, "%d PSL_A%d_y MM\n", PSL_IZ (PSL, x), annot_pos);
if (label_c && label_c[i] && label_c[i][0])
strncpy (string, label_c[i], GMT_LEN256-1);
else
gmtlib_get_coordinate_label (GMT, string, &GMT->current.plot.calclock, format, T, knots[i]); /* Get annotation string */
double angle_text2=AxisTickLabel_Flip_GMT_X_Y(GMT,text_angle); //坐標(biāo)軸標(biāo)注文字翻轉(zhuǎn)
PSL_plottext (PSL, 0.0, 0.0, -font.size, string, angle_text2, justify, form);
AxisTickLabel_Flip_GMT_X_Y(GMT,text_angle); //恢復(fù)：坐標(biāo)軸標(biāo)注文字翻轉(zhuǎn)
}
if (!faro) PSL_command (PSL, "/PSL_A%d_y PSL_A%d_y PSL_AH%d add def\n", annot_pos, annot_pos, annot_pos);
這里省略細(xì)節(jié)士复，感興趣的可以直接看我的源代碼
第二個(gè)BUG修復(fù)之后的結(jié)果
剖面文字的鏡像和旋轉(zhuǎn)角度問題
第二個(gè)BUG修復(fù)之后，坐標(biāo)軸上的文字正確歸為了翩活。但是剖面上的文字還存在問題：文字鏡像以及某些情況下存在一個(gè)角度偏轉(zhuǎn)阱洪。找到繪制文字的代碼位置：pstext.c文件中的GMT_pstext函數(shù)，代碼太長菠镇，這里只貼關(guān)鍵代碼附近的幾行：
c_txt[m] = strdup (curr_txt);
c_x[m] = plot_x;
c_y[m] = plot_y;
c_just[m] = T.block_justify;
c_font[m] = T.font;
m++;
}
else {
PSL_plottext (PSL, plot_x, plot_y, T.font.size, curr_txt, T.paragraph_angle, T.block_justify, fmode);
}
if (Ctrl->A.active) T.paragraph_angle = save_angle; /* Restore original angle */
}
就是在這里調(diào)用PSL_plottext函數(shù)進(jìn)行文字繪制冗荸，而PSL_plottext函數(shù)定義在postscriptlight.h里面，函數(shù)實(shí)現(xiàn)在postscriptlight.c里面蚌本。進(jìn)入這兩個(gè)文件里面發(fā)現(xiàn)不太容易傳遞GMT_CTRL *GMT參數(shù)(GMT這個(gè)數(shù)據(jù)結(jié)構(gòu)中包含了很多很多信息)，為了兼容性轴猎。重新自定義一個(gè)int PSL_plottext_mirror (int isMirrow,struct PSL_CTRL *PSL, .....)的函數(shù)，在這個(gè)函數(shù)中增加一個(gè)參數(shù)isMirrow锐峭。用PSL_plottext_mirror替換上面代碼段中的PSL_plottext沿癞。與前面的幾個(gè)問題的修復(fù)方法類似抛寝，在pstext.c中新增函數(shù)int Text_Flip_GMT_X_Y(struct GMT_CTRL *GMT)判斷何時(shí)需要對(duì)文字進(jìn)行鏡像翻轉(zhuǎn)盗舰，然后將這個(gè)參數(shù)傳遞給新定義的文字繪制函數(shù)PSL_plottext_mirror去執(zhí)行相應(yīng)的操作川陆。
c_txt[m] = strdup (curr_txt);
c_x[m] = plot_x;
c_y[m] = plot_y;
c_just[m] = T.block_justify;
c_font[m] = T.font;
m++;
}
else {
// PSL_command(PSL,"%%繪制文字主控函數(shù)\n");
// Text_Flip_GMT_X_Y(GMT); //文字鏡像翻轉(zhuǎn)
PSL_plottext_mirror (Text_Flip_GMT_X_Y(GMT),PSL, plot_x, plot_y, T.font.size, curr_txt, T.paragraph_angle, T.block_justify, fmode);
// Text_Flip_GMT_X_Y(GMT); //恢復(fù)
}
if (Ctrl->A.active) T.paragraph_angle = save_angle; /* Restore original angle */
}
使用自定義的文字繪制函數(shù)PSL_plottext_mirror。
至此失仁，上面三個(gè)問題已全部修復(fù)完畢控轿，修復(fù)bug之后的GMT姑且稱之為ModernGMT(GMT私房菜??)。運(yùn)行結(jié)果如下在抛。
ModernGMT運(yùn)行結(jié)果
本人修復(fù)BUG之后的版本稱之為ModernGMT刚梭，這是一個(gè)持續(xù)更新的私房版本望浩，不僅修復(fù)各種bug還會(huì)添加一些有用的地球物理重磁位場相關(guān)的程序缘回。而且我每個(gè)星期都會(huì)同步官方更新酥宴，所以ModernGMT總是比官方版本新一些且功能多一些授滓。但是由于我時(shí)間有限不能暫時(shí)還沒有把自己做的更新推送給官方團(tuán)隊(duì)在孝，以后會(huì)考慮私沮！
后記
經(jīng)過上面的修復(fù)之后，目前繪制正常的三維圖應(yīng)該是沒有什么大的問題了晰搀。如果有人嘗試一些奇怪的想法，比如旋轉(zhuǎn)角度很畸形或者其他的參數(shù)不和常理，不保證不出問題。
依然存在的問題
當(dāng)使用-JM投影方式的時(shí)候，如果數(shù)據(jù)范圍很大，比如100度X80度的跨度店溢，由于投影計(jì)算過程中的誤差等原因，可能會(huì)出現(xiàn)切片不貼合的現(xiàn)象壕吹。如下圖所示：
解決這個(gè)問題的技巧就是直接在GMT繪圖代碼里面調(diào)整一下切片的長度即可，比如增加一個(gè)小的倍數(shù)耳贬，稍微調(diào)整然后看著沒問題就行泳姐。以后有空了會(huì)從代碼里面直接解決
# # profile along lat
width_fig_y0=`echo $width_fig_y | awk '{print $1+0.3}'`
gmt psbasemap -R$range_LatZLon -JX$width_fig_y/$width_fig_z -JZ$width_fig_x -BS -Ba --MAP_ANNOT_OFFSET=$offset_ticklabel_lon -px$angle_view/$pos_profile_lon
小范圍是沒有問題的，比如下面這張圖
兩個(gè)剖面上的數(shù)據(jù)是隨便用程序生成的高斯分布數(shù)據(jù)摇肌，不代表任何意義，只是為了說明在剖面切片中繪制grdimage的問題围小。
GMT更新程序獲取方法
由于時(shí)間和精力有限变秦，暫不為伸手黨直接開放源代碼和程序。源代碼以及編譯之后的安裝文件我會(huì)上傳到云端框舔，請(qǐng)到我的資源庫中查看代碼和程序獲取方法蹦玫。
