绘制三维有限断层模型

示例贡献者:

菲林可乐（作者）、李水平（作者）、陈箫翰（修订）

使用 GMT 绘制有限断层模型的关键，是利用 plot3d -C -L 命令选项绘制一个个小多边形并填充颜色。 这就要求用户在绘图之前必须以特定的格式准备数据。以绘制三维有限断层模型为例，用户必须准备好如下格式的多段数据:

# 数据头段
>  -Z1.800000e-01
 # 经度 纬度 深度
 99.3785  34.5324  -0.0010
 99.3574  34.5279  -0.0010
 99.3569  34.5297  -1.9900
 99.3779  34.5342  -1.9900
>  -Z5.000000e-02
 99.3574  34.5279  -0.0010
 99.3363  34.5233  -0.0010
 99.3358  34.5252  -1.9900
 99.3569  34.5297  -1.9900
...

不同小多边形的数据以 > 符号分开，每段数据给出这个小多边形的所有顶点坐标，按顺时针或逆时针顺序排列。 > 符号后面以 -Z 选项给出这个小多边形对应的 Z 值。绘图时将根据这个 Z 值对多边形填充不同的颜色。 这种多段数据格式的详细说明请参考 table_ascii_id4 。

有些三维有限断层反演程序给出的不是以上格式的数据，而是给出每一个小矩形中心坐标、矩形长宽、走向倾向倾角。例如:

     lat_deg     lon_deg    depth_km   length_km    width_km  slp_strk_m  slp_ddip_m    slp_am_m  strike_deg     dip_deg    rake_deg sig_stk_MPa sig_ddi_MPa sig_nrm_MPa
     34.7982     97.5406      0.8874      1.9546      1.8890      0.2690     -0.3206      0.4186    102.8566     69.9500     50.0000      1.7759     -1.5403      0.2174
     34.7925     97.5389      2.6611      1.9548      1.8890      0.2156     -0.2569      0.3354    102.8566     69.8445     50.0000      0.8002     -0.7488      0.0076
     34.7868     97.5372      4.4338      1.9549      1.8890      0.2255     -0.2688      0.3508    102.8567     69.7391     50.0000      0.7526     -0.8695     -0.0233
...

对于这种数据可以参考 plotting_3D_FFM 进行转换。

下面的示例展示了如何基于以上方法绘制三维有限断层模型。

断层示例数据下载： slip_3dgrid.gmt

#!/usr/bin/env bash
#
# 绘制三维有限断层模型

gmt begin 3d_faultslip
    gmt set MAP_FRAME_TYPE plain #设置边框类型
    gmt set MAP_GRID_PEN 0p,gray,- #设置网格线
    
    gmt basemap -R97.3/99.5/34.3/35/-25/0 -JX15c/10c -JZ5c -Bxa0.25fg -Bya0.25fg -Bza10fg+l"Depth (km)" -BwSEnZ -p160/20
    # 为滑移位错量制作 CPT。滑移位错量最小值为 0， 最大值为 5.5
    gmt makecpt -Cwhite,lightblue,yellow,red -T0/5.5/0.2 -Z
    # 绘制有限断层模型
    gmt plot3d slip_3dgrid.gmt -C -L -W0p,gray -p
    
    # 绘制 colorbar
    gmt set FONT_ANNOT_PRIMARY 7p FONT_LABEL 8p
    gmt set MAP_ANNOT_OFFSET_PRIMARY 2p MAP_FRAME_PEN thinner MAP_TICK_LENGTH 2p
    gmt colorbar -C -DjBL+w2.5c/0.3c+h+o1.3c/2c+ml -Ba1 -Bx+l"slip(cm)"
gmt end show
../../_images/ca530adf4c88d02f61decf4a4b296529.png