GMT三維繪圖有BUG? 修復(fù)它！

前期一篇文章:Modern GMT Series：Slice in 3D View (三維切片圖)中提到了在科研作圖中會(huì)經(jīng)常遇到三維作圖的問(wèn)題，而且GMT可以做三維圖且導(dǎo)出的圖片質(zhì)量非常高！但是也提到了當(dāng)前GMT版本中的三維繪圖存在漏洞。主要兩個(gè)問(wèn)題：（1）切片位置錯(cuò)亂；（2）三維文字鏡像。就這兩個(gè)問(wèn)題極大的限制了GMT三維作圖的應(yīng)用。因?yàn)樽罱撐睦锩嬉鋈S立體圖(地理坐標(biāo)+笛卡爾坐標(biāo)+海底地形+剖面切片等)，找了找也沒(méi)找到更好的解決方案，干脆自己動(dòng)手豐衣足食——修改GMT的源代碼。經(jīng)過(guò)一天半的努力，已經(jīng)比較完美的修復(fù)了這兩個(gè)漏洞。目前官方版本中并沒(méi)有修復(fù)此bug。 下面就介紹一下修復(fù)漏洞的過(guò)程和最終效果！

在線視頻講解

先上一張BUG修復(fù)后的效果圖鑒賞一下

三維地形圖+剖面切片+線+文字+方向標(biāo)+ModernGMTlogo

注： 由于我的論文正在審稿中，不便上傳較為復(fù)雜的繪圖結(jié)果，暫用一個(gè)相對(duì)簡(jiǎn)單的圖來(lái)代替。

簡(jiǎn)單例子

問(wèn)題描述：繪制一個(gè)三維立方體并在每個(gè)面上標(biāo)注文字以及設(shè)置每個(gè)面具有不同顏色且半透明；水平旋轉(zhuǎn)45度，然后x軸和y軸標(biāo)注坐標(biāo)標(biāo)簽。

GMT官方版本運(yùn)行的結(jié)果

GMT6.0三維繪圖結(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上肯定是沒(méi)有相關(guān)的文件和函數(shù)。直接復(fù)制上面的代碼應(yīng)該是可以直接出圖的。注意這是Mac系統(tǒng)下的腳本，也可以在Linux系統(tǒng)下使用，但是windows系統(tǒng)下需要將變量引用從$xxx改為%xxx%。

都有什么BUG？

從上面的第二個(gè)圖中可以看出，官方版本的三維繪圖結(jié)果存在一下幾個(gè)問(wèn)題：

  1. 切片位置錯(cuò)亂：需要向y軸正方向平移一個(gè)量

  2. 文字出現(xiàn)鏡像：x軸和y軸的標(biāo)注文字都變成了鏡像且有一個(gè)翻轉(zhuǎn)角度

  3. 普通文字鏡像：除了坐標(biāo)軸的標(biāo)注文字外，剖面切片上的文字的角度也有錯(cuò)亂和鏡像現(xiàn)象

如何修復(fù)BUG？

簡(jiǎn)單介紹以上三個(gè)BUG是如何修復(fù)的，具體見(jiàn)修復(fù)后的源代碼(我的資源庫(kù)中有介紹如何獲取源代碼和可用軟件)。

切片位置錯(cuò)亂問(wèn)題

GMT繪圖是基于PostScript的，三維繪圖的坐標(biāo)旋轉(zhuǎn)理論見(jiàn)GMT三維坐標(biāo)旋轉(zhuǎn)理論。經(jīng)測(cè)試發(fā)現(xiàn)問(wè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)矩陣寫(xiě)入ps文件就能得到三維的效果。

** GMT_Z的情況都是OK的，正如上圖所示三維坐標(biāo)軸主框架沒(méi)錯(cuò)。但是 GMT_X和 GMT_Y**的情況，也就是分別于X軸和Y軸垂直的切片的情況，出現(xiàn)了問(wèn)題。經(jīng)測(cè)試發(fā)現(xiàn)，問(wèn)題在于e, f計(jì)算錯(cuò)誤。三種情況的e, f值需要一致才不會(huì)發(fā)生未知平移錯(cuò)亂。

解決方案

這里有一個(gè)很容易的解決方案就是在case GMT_Z計(jì)算得到e,f值的時(shí)候，打開(kāi)一個(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é)果如下：

修復(fù)剖面位置問(wèn)題后的結(jié)果

坐標(biāo)軸標(biāo)注文字鏡像問(wèn)題

上面一步修復(fù)了剖面位置問(wèn)題，剖面坐標(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命令真是不容易，參見(jiàn)一個(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除了鏡像問(wèn)題還存在一個(gè)180度的旋轉(zhuǎn)角度問(wè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é)果

坐標(biāo)軸label和tick label鏡像問(wèn)題修復(fù)后的效果

剖面文字的鏡像和旋轉(zhuǎn)角度問(wèn)題

第二個(gè)BUG修復(fù)之后，坐標(biāo)軸上的文字正確歸為了。但是剖面上的文字還存在問(wèn)題：文字鏡像以及某些情況下存在一個(gè)角度偏轉(zhuǎn)。找到繪制文字的代碼位置：pstext.c文件中的GMT_pstext函數(shù)，代碼太長(zhǎng)，這里只貼關(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è)問(wèn)題的修復(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è)問(wèn)題已全部修復(fù)完畢，修復(fù)bug之后的GMT姑且稱之為ModernGMT(GMT私房菜??)。運(yùn)行結(jié)果如下。

ModernGMT運(yùn)行結(jié)果

本人修復(fù)BUG之后的版本稱之為ModernGMT，這是一個(gè)持續(xù)更新的私房版本，不僅修復(fù)各種bug還會(huì)添加一些有用的地球物理重磁位場(chǎng)相關(guān)的程序。而且我每個(gè)星期都會(huì)同步官方更新，所以ModernGMT總是比官方版本新一些且功能多一些。但是由于我時(shí)間有限不能暫時(shí)還沒(méi)有把自己做的更新推送給官方團(tuán)隊(duì)，以后會(huì)考慮！

ModernGMT三維繪圖結(jié)果

后記

經(jīng)過(guò)上面的修復(fù)之后，目前繪制正常的三維圖應(yīng)該是沒(méi)有什么大的問(wèn)題了。如果有人嘗試一些奇怪的想法，比如旋轉(zhuǎn)角度很畸形或者其他的參數(shù)不和常理，不保證不出問(wèn)題。

依然存在的問(wèn)題

當(dāng)使用-JM投影方式的時(shí)候，如果數(shù)據(jù)范圍很大，比如100度X80度的跨度，由于投影計(jì)算過(guò)程中的誤差等原因，可能會(huì)出現(xiàn)切片不貼合的現(xiàn)象。如下圖所示：

大范圍三維繪圖切片長(zhǎng)度計(jì)算誤差

解決這個(gè)問(wèn)題的技巧就是直接在GMT繪圖代碼里面調(diào)整一下切片的長(zhǎng)度即可，比如增加一個(gè)小的倍數(shù)，稍微調(diào)整然后看著沒(méi)問(wèn)題就行。以后有空了會(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

小范圍是沒(méi)有問(wèn)題的，比如下面這張圖

西南印度洋脊龍旗熱液區(qū)附近的ETOPO1水深數(shù)據(jù)

兩個(gè)剖面上的數(shù)據(jù)是隨便用程序生成的高斯分布數(shù)據(jù)，不代表任何意義，只是為了說(shuō)明在剖面切片中繪制grdimage的問(wèn)題。

GMT更新程序獲取方法

由于時(shí)間和精力有限，暫不為伸手黨直接開(kāi)放源代碼和程序。源代碼以及編譯之后的安裝文件我會(huì)上傳到云端，請(qǐng)到我的資源庫(kù)中查看代碼和程序獲取方法。

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書(shū)系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容