三维数据体生成与切片

我们在计算得到的数据通常是非等间隔网格，为了用GMT进行2D可视化，需要将数据网格化后进行插值。

1. 将文本文件转换为三维数据体#

我们拿到的原始文件是包含经度、纬度、深度、和每个数据点的数据（速度、温度、密度等）组成的文本格式文件，以这次的数据为例，我将Excel文件直接复制到文本中形成了以空格为分隔符的文本文件，以下列出部分行：

1
10.55467676 108.4399258 175.207137 0.01217283 -0.015367995 -9.20E-05 1681.16272 3291.964844 5848311808 1.00E+21
2
10.55113868 107.7219827 175.2071289 0.012363088 -0.01525117 -0.00010885 1675.885254 3292.487305 5847836672 1.00E+21
3
10.25584183 108.1032219 146.0060139 0.012991064 -0.014229504 -9.73E-05 1679.978882 3292.082031 4886236672 1.00E+21
4
9.998090145 108.4349265 146.030284 0.01289901 -0.014223722 -8.98E-05 1681.883859 3291.893512 4886096558 1.00E+21
5
9.998090145 107.8137658 146.0285942 0.013020781 -0.014152945 -0.000117022 1675.698183 3292.505842 4885533575 1.00E+21
6
10.55468488 108.439848 146.0060673 0.012965109 -0.014302315 -7.66E-05 1684.43457 3291.640869 4887052800 1.00E+21
7
10.5511462 107.7219104 146.0062614 0.013087932 -0.014204981 -0.000102527 1677.521362 3292.325439 4886397440 1.00E+21
8
10.25585045 108.1031426 116.804997 0.013968474 -0.01280891 -1.20E-05 1596.835083 3300.313232 3924829952 1.00E+21
9
9.998090145 108.4348529 116.8293289 0.013916044 -0.012845078 -5.31E-06 1600.383762 3299.961929 3924717946 1.00E+21

网格化说明#

由于这样的数据点是非等间隔的，没有办法进行切片时所需要的插值操作，所以我先将数据进行了三维网格化。这里我选用了Python语言和它的科学计算模块Numpy、Scipy进行网格化。Matlab也可以实现类似功能，核心函数也叫griddata。

1
import numpy as np
2
from scipy.interpolate import griddata
3

4
def  read_dat(fname):
5
   data = np.loadtxt(fname)
6
   data[:, 3] *=  100
7
   data[:, 4] *=  100
8
   data[:, 5] *=  100
9
   return data
10

11
def dat2npz(fname='data', latmin=10, latmax=60, lonmin=70, lonmax=140, depmin=0, depmax=700):
12
   # 在三个纬度上设置研究区域和网格间隔，网格间隔根据需要选取。
13
   lat = np.arange(latmin, latmax+0.1, 0.1)
14
   lon = np.arange(lonmin, lonmax+0.1, 0.1)
15
   dep = np.arange(depmin, depmax+5, 5)
16
   new_dep, new_lat, new_lon = np.meshgrid(dep, lat, lon, indexing='ij')
17
   # 读取数据列表
18
   data =  read_dat(fname)
19
   # 设置原始数据点位置
20
   points = np.array([data[:, 2], data[:, 0], data[:, 1]]).T
21
   # 对每一种数据进行网格化
22
   vew =  griddata(points, data[:, 3], (new_dep, new_lat, new_lon), method='linear')
23
   vns =  griddata(points, data[:, 4], (new_dep, new_lat, new_lon), method='linear')
24
   vrad =  griddata(points, data[:, 5], (new_dep, new_lat, new_lon), method='linear')
25
   t =  griddata(points, data[:, 6], (new_dep, new_lat, new_lon), method='linear')
26
   dens =  griddata(points, data[:, 7], (new_dep, new_lat, new_lon), method='linear')
27
   # 将网格化后每个维度坐标和数据保存在一个.npz文件中
28
   np.savez('data_0.1.npz', dep=dep, lat=lat, lon=lon, vew=vew, vns=vns, vrad=vrad, t=t, dens=dens)

2. 平面和剖面的切片#

用GMT绘制二维图像需要从上一步网格好的数据文件中对某个平面或剖面进行切片。对平面的切片简单一些，对剖面的切片则需要先对大圆弧路径进行插值，生成数据点坐标，再用这些坐标对三维数据体进行插值。

2.1 对数据的平面切片 (`cut_plane.py`)#

1
import numpy as np
2
from scipy.interpolate import interpn
3

4
# 在指定的区域内生成数据点坐标，理论上说这些点也可以是不等间隔的，但是GMT绘图时依旧需要网格化，所以这里还是使用等间隔插值。
5
def init_plane(lat1=15, lon1=65, lat2=55, lon2=140, val=2, depth=200):
6
    lat = np.arange(lat1, lat2, val)
7
    lon = np.arange(lon1, lon2, val)
8
    points = []
9
    for la in lat:
10
        for lo in lon:
11
            points.append([depth, la, lo])
12
    return np.array(points)
13

14
# 生成的数据点对三维数据体插值
15
def interall(data, points):
16
    inter_data = np.zeros([points.shape[0], 2])
17
    inter_data[:, 0] = interpn((data['dep'], data['lat'], data['lon']), data['vew'], points, bounds_error=False, fill_value=None)
18
    inter_data[:, 1] = interpn((data['dep'], data['lat'], data['lon']), data['vns'], points, bounds_error=False, fill_value=None)
19
    # 返回时同时将数据点坐标那三列于插值得到的数据拼贴后一起输出
20
    return np.hstack((points[:, 1:3], inter_data))
21

22

23
def cut_plane(fname='data_0.1.npz'):
24
    data=np.load(fname)
25
    points = init_plane()
26
    inter_data = interall(data, points)
27
    # 最后输出到一个文本文本文件
28
    np.savetxt('vel_200.dat', inter_data)

那么这时输出的数据包括了4列，依次为经度，纬度，南北方向的速度和东西方向的速度。

2.2 对数据的剖面切片 (`cut_sec.py`)#

虽热都是对数据体的切片，但是剖面的切片要更复杂一些。主要的难点是：

对大圆弧路径的等距离划分。其中要考虑球面三角函数和椭球矫正（有现成代码）。但是这里为了方便直接调用了GMT里的project命令，直接生成横向网格，更加方便。
对原先的南北和东西分量进行坐标旋转，旋转角度是剖面两个端点之间的方位角（请参考IRIS 的一个教程）。

那么先导入函数库

1
import numpy as np
2
from scipy.interpolate import interpn
3
import subprocess
4
import re
5
import sys
6
from seispy_distaz import *

生成数据点#

以下是生成据点的函数。需要说明的是gmt project是GMT内一个用于投影点的命令，这里直接在Python中调用了，这条命令的具体说明请参考GMT中文手册。所以如果系统中没有GMT，就只能尝试用别的方法了。

1
def gen_points(lat1, lon1, lat2, lon2, val=10, depmax=700):
2
    # 生成GMT命令字符串
3
    cmd = 'gmt project -C{}/{} -E{}/{} -G{} -Q'.format(lon1, lat1, lon2, lat2, val)
4
    # 开启一个命令行子进程，并将cmd传入
5
    p = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE)
6
    # 读取网格点经纬度和相对第一个点的距离
7
    line_pos = np.loadtxt(p.stdout)
8
    points = []
9
    # 对深度进行网格有depth和val确定最大深度和间隔
10
    for dep in np.arange(0, depmax+val, val):
11
        for row in line_pos:
12
            # 生成数据列表，这里/111是将千米转化为度
13
            points.append([dep, row[1], row[0], row[2]/111])
14
    return np.array(points)

对矢量数据进行坐标系旋转#

大部分都是标准公式。需要一提的是distaz来自seispy_distaz.py，一个在原版distaz基础上改进过的版本，说不定以后做其他处理还能用的上。

1
def rotate(vew, vns, baz):
2
    angle = np.mod(baz+180, 360)
3
    vr = vns*cosd(angle) + vew*sind(angle)
4
    vt = vew*cosd(angle) - vns*sind(angle)
5
    return vr, vt
6

7
def velproj(vew, vns, lat1, lon1, lat2, lon2):
8
    baz = distaz(lat1, lon1, lat2, lon2).baz
9
    vr, vt = rotate(vew, vns, baz)
10
    return vr, vt

数据插值#

插值和平面切片是类似的。只有以下几点不同：

速度为坐标系转化后的速度。
需要对垂直方向速度，温度和密度等进行插值。

1
def interall(data, vr, vt, points):
2
    inter_data = np.zeros([points.shape[0], 5])
3
    inter_data[:, 0] = interpn((data['dep'], data['lat'], data['lon']), vr, points[:, 0:3], bounds_error=False, fill_value=None)
4
    inter_data[:, 1] = interpn((data['dep'], data['lat'], data['lon']), vt, points[:, 0:3], bounds_error=False, fill_value=None)
5
    inter_data[:, 2] = interpn((data['dep'], data['lat'], data['lon']), data['vrad'], points[:, 0:3], bounds_error=False, fill_value=None)
6
    inter_data[:, 3] = interpn((data['dep'], data['lat'], data['lon']), data['t'], points[:, 0:3], bounds_error=False, fill_value=None)
7
    inter_data[:, 4] = interpn((data['dep'], data['lat'], data['lon']), data['dens'], points[:, 0:3], bounds_error=False, fill_value=None)
8
    return np.hstack((points, inter_data))

这个脚本中一些其他小技巧#

由于剖面不止一条，所以我按照我在地震学上的处理方式做了一个剖面列表，给每个剖面编上序号，切剖面和画剖面时只需要改动剖面编号。剖面列表lines.lst内容如下：

1
25 80 45 87 a
2
24 85 25 113 b
3
25 80 40 118 c

最后最好空一行，否则有的语言读取可能会有bug。

所以在剖面切片的脚本中就加入了输入剖面编号，剖面切片的功能。以下是输入编号输出对应剖面经纬度的函数：

1
def load_line(line_name, line_file='lines.lst'):
2
    with open(line_file) as f:
3
        line_content = f.readlines()
4
    for line_str in line_content:
5
        line_str = line_str.strip()
6
        if line_str[-1] == line_name:
7
            return [float(value) for value in line_str.split()[:-1]]