从 3D 数据中截取剖面

本文所说的 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 代码如下:

1
2
3
4
5
6
7
8
9
10
11
#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

截取了这个剖面之后,下面就是 xyz2grdgrdimage 了,不多说。

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

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

这里用 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 取该深度在测线上的值;然后取下一个深度,依次处理下去。

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#!/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,然后把结果收集起来即可。

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