✍️ 朱邓达  •  📅 2025-05-31

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

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 png,pdf
    # 构建底图
    gmt basemap -Rg -JQ10c -Baf
    # 绘制地形
    gmt grdimage @earth_relief_05m

    # 创建测线，点距为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 6.5 版本下，grdinterpolate -S 和 -T 同时使用存在 bug，
    # gmt grdinterpolate S362ANI_percent.nc?dvs -Sprof $T > samp
    # 以上命令可暂时替换为下方命令来回避此 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