本文所说的 3D 数据是指 XYZD 数据。常见的情况是数据中包含四列：经度、纬度、深度、值。

GMT 几乎是不可以处理 3D 数据的，比如与网格相关的命令只支持 2D netCDF 文件。另外，GMT 绘制 3D 图的效果也很一般，因而通常情况下都需要从 3D 数据中截取剖面，绘制剖面图。

本文试着介绍如何从 3D 数据中截取自己想要的剖面。

数据准备

因为我手头没有真实的数据，所以这里虚构一个 3D 数据。就示例而言，虚构的数据也许比真实的数据的演示效果更好。

虚构数据如下：

  • X 的范围为 0 到 9，间隔为 1；
  • Y 的范围为 10 到 19，间隔为 1；
  • Z 的范围为 20 到 29，间隔为 1；
  • D 的值由 X、Y、Z 的值唯一决定；

生成虚构数据的 C 代码如下：

#include <stdio.h>

int main()
{
    int i, j, k;
    for (i=0; i<10; i++)
        for (j = 10; j < 20; ++j)
            for (k = 20; k < 30; ++k)
                printf("%d %d %d %d\n", i, j, k, i*100+(j-10)*10+(k-20) );
    return 0;
}

截取垂直于坐标轴的剖面

最简单的情况是从 3D 速度结构中截取某个深度的速度结构，比如 Z=24 这个垂直于 Z 轴的平面。此时直接用 gawk 将 3D 数据中 Z=24 的点提取出来:

gawk '$3==24 {print $1,$2,$4}' data.xyzd > Z24.xyd

如果数据的间隔比较小，用 $3==4 去判断可能有些危险，更安全点的方法是:

gawk '$3>23.999 && $3<24.001 {print $1,$2,$4}' data.xyzd > Z24.xyd

截取了这个剖面之后，下面就是 xyz2grd 和 grdimage 了，不多说。

截取垂直于坐标平面的平面

最常见的情况是，地表有一个测线，需要做一个沿着这个测线的垂直于地表的剖面。

这里用 project 命令虚构起点在 (2,12)，终点为 (5,16)，间距为 1 的剖面:

project -C2/12 -E5/16 -G1 > track.data

track.data 内容如下，三列数据分别为经度、纬度以及离起点的距离（大圆距离）:

2   12  0
2.59831503664   12.8115096006   1
3.20049364771   13.6216595304   2
3.80680955458   14.4303509263   3
4.41754162203   15.2374825707   4
5   16  4.94662175502

思路其实很简单，如上面那个例子，先用 gawk 取深度 Z=20 处的平面，然后将其转换为网格文件，并用 grdtrack 取该深度在测线上的值；然后取下一个深度，依次处理下去。

下面给出了一个示例脚本，该脚本可以解决本文的问题，但灵活性不够，所以使用者需要根据自己的情况改写或重写，理解原理最重要。

#!/bin/bash

# 用于清除上一次执行该脚本时生成的文件
if [-e track.profile]; then
    rm track.profile
fi

# 对深度进行循环
for ((i = 20; i < 30; i++)); do
    # 取某个深度的剖面
    gawk -v dep=$i '$3==dep {print $1,$2,$4}' data.xyzd > $i.xyd
    # 将数据转换为 netCDF 格式
    xyz2grd $i.xyd -R0/9/10/19 -I1/1 -G$i.nc -V
    # 在该深度处做 grdtrack，最终结果重定向到 track.profile 中
    grdtrack track.data -G$i.nc | gawk -v dep=$i '{print $3, dep, $4}' >> track.profile

    rm $i.xyd $i.nc
done;

几点说明:

  1. 使用 gawk 时，涉及到要向 gawk 传递一个 bash 中的变量，需要使用 -v 选项；
  2. grdtrack 命令的输出有四列，分别是经度、纬度、离起点的距离以及 Z；
  3. 绘制剖面时，一般需要离起点的距离，以及深度这两个信息，所以用 gawk 处理了一下；
  4. 所有的循环的输出都重定向到文件 track.profile，最终该文件中有三列，离起点的距离、深度以及值；

任意一个倾斜剖面

这个说难不难，说简单也不简单。首先，如何用一些参数去定义这个剖面就是个问题。假设定义好了，就可以通过某些计算，得到这个剖面与地表的交线以及剖面与各个深度水平面的交线。

有了这些交线，就可以用与上面的例子类似的方法对每个深度处做 grdtrack，然后把结果收集起来即可。

思路有了，具体的就不写啦。