✍️ 朱邓达 • 📅 2025-09-29

# 绘制三维速度模型的任意垂直切片

grdinterpolate 模块从三维网格中计算任意点位下方不同层处的值。 以地震学中的三维速度模型为例，三维网格水平方向为 (lon, lat)，深度方向为 depth ， 其中沿深度方向网格分为 nz 层，网格节点值以速度为例， grdinterpolate 模块采样流程为

1. 通过 -S 给定地表某点位坐标 (lon0, lat0) ，使用 grdtrack 模块采样每层二维网格中该点位的值，共得到 nz 个值。

2. 如果使用了 -T 自定义深度范围和间距，则使用 sample1d 模块将上一步得到的数据采样到自定义深度。

3. 使用 -E 定义测线则切片输出为网格文件， 使用 -S 定义点位则切片输出为ASCII表格文件。

由于目前难以导出通过 -E 选项生成的测线，为了更灵活地自定义测线， 这里使用 project 模块先生成测线，使用 grdinterpolate -S 计算切片， 再使用 xyz2grd 模块将结果转为网格文件。

    #!/usr/bin/env bash
    #
    # 绘制三维速度模型的任意垂直切片

    # 下载示例模型文件
    # wget "https://ds.iris.edu/files/products/emc/emc-files/S362ANI_percent.nc"

    gmt begin 3d_vert_slice
        # 构建底图
        gmt basemap -Rg -JQ10c -Baf
        # 绘制地形
        gmt grdimage @earth_relief_05m_g

        # 创建测线，点距为1度
        D1=0
        D2=180
        dD=1
        gmt project -C110/20 -A90 -L$D1/$D2 -G$dD > prof
        # 绘制测线
        gmt plot prof -W1.5p,red

        # 定义沿深度方向的插值
        T1=25
        T2=2890
        dT=10
        T="-T$T1/$T2/$dT"

        # 计算该测线下的垂直切片
        gmt grdinterpolate S362ANI_percent.nc?dvs -Sprof $T > samp
        # 在 GMT 6.5 以及更低版本中，grdinterpolate -S 和 -T 同时使用存在 bug。
        # 以上命令可暂时替换为下方命令来回避此 bug
        # gmt grdinterpolate S362ANI_percent.nc?dvs -Sprof | gmt sample1d $T -N2 > samp

        # 将表数据转为网格文件
        gmt xyz2grd samp -Gvs_slice.nc -i3,2,4 -R$D1/$D2/$T1/$T2 -I$dD+e/$dT+e

        # 创建色标
        gmt makecpt -Cseis -D -T-2/2

        # 构建底图，横坐标的距离为角度，故可使用极坐标投影
        gmt basemap -Rvs_slice.nc -JP10c+fp+a+t90 -Baf -Y-7c
        # 绘制切片
        gmt grdimage vs_slice.nc
        # 绘制色标棒
        gmt colorbar -Bxa1+l"dVs (%)"

    gmt end show