✍️ 田冬冬 , 姚家园 • 📅 2022-04-29

# 绘制地球内部主要界面

在利用地震波研究地球深部结构时，经常需要绘制震相在深度剖面下的射线路径， 同时也需要绘制地球内部的主要界面。

震相的射线路径可以用 TauP 提供的 `taup_path` 命令计算得到，然后在极坐标系 `-JP` 中绘制。难点在于如何绘制几个主要界面。

本示例将展示如何绘制震中距为30度的 PcP 和 PKiKP 震相的射线路径， 同时绘制地球内的 410、660界面以及 CMB 和 ICB。 脚本中调用 TauP 计算射线路径。 没有安装 TauP 的用户可以直接下载示例射线路径数据 `PKiKP.raypath.gmt` 和 `PcP.raypath.gmt`。

绘图脚本如下：

    #!/usr/bin/env bash
    #
    # 绘制地球内部主要界面
    #

    gmt begin earth-discontinuities
    gmt set MAP_GRID_PEN_PRIMARY 1p
    # 绘制地表
    gmt basemap -JP10c+a+t15 -R-10/40/0/6371 -Byg6371 -BS
    # 绘制 410 界面
    gmt basemap -Byg6371+5961 -BS
    # 绘制 660 界面
    gmt basemap -Byg6371+5711 -BS
    # 绘制 CMB
    gmt basemap -Byg6371+3480 -BS
    # 绘制 ICB
    gmt basemap -Byg6371+1221 -BS

    # 计算震相的射线路径，用户需安装 TauP，并取消以下注释
    # taup_path -mod prem -ph PcP -h 300 -deg 30 -o PcP.raypath
    # taup_path -mod prem -ph PKiKP -h 300 -deg 30 -o PKiKP.raypath

    # 绘制震相的射线路径
    gmt plot PcP.raypath.gmt -W1p,blue
    gmt plot PKiKP.raypath.gmt -W1p,red

    # 绘制震源和台站位置
    gmt plot -S -Gblack -N << EOF
    0 6071 0.4c a
    30 6471 0.4c i
    EOF

    # 添加标注
    gmt text -F+f11p+a -N << EOF
    -6 6170 20 Surface
    -6 3280 20 CMB
    -3 1020 20 ICB
    16 4100 0 @;blue;PcP@;;
    37 1600 0 @;red;PKiKP@;;
    EOF
    gmt end show

脚本使用了五次 `basemap` 命令，分别绘制五个界面：

- 第一次用于绘制地表：`-Byg6371` 表明要在 Y 方向（极坐标下即 R 方向）以6371为间隔绘制网格线， 由于地球半径是6371，所以理论上会在R=0和 R=6371 两处绘制网格线。 `-BS` 的作用是只绘制 R=6371 处的网格线，且只绘制网格线而不绘制刻度

- 第二次用于绘制410界面：与前一命令基本相同，唯一的区别是 `-Byg6371+5961` 中多了 `+5961`，其作用是定义网格线的起算点，即此时将在 R=5961（即410界面）处绘制网格线

- 其他同理