我们通过计算获得了地壳三维速度模型、Moho面起伏，为例对比壳内低速区与Moho面起伏的关系，可以将这些结果用Paraview进行三维可视化。由于数据格式的不同，我们需要借助GMT和PyVista辅助绘制。

准备工作#

原材料包括三维速度模型、Moho面起伏数据以及地形数据、断层数据等，我们需要将这些不同格式的文件转换成Paraview识别的.vtk文件，这里主要用到的工具是PyVista，这是一个Python与Visualization Toolkit (VTK)的接口，通过这个模块我们可以将Python常用的np.ndarray数据转换成VTK文件格式。

定义一个类#

我们先定义一个Python类，初始化一些每种数据共用的参数，如经纬度范围，深度范围等。

1
import numpy as np
2
from scipy.interpolate import griddata
3
from scipy import ndimage
4
import scipy.signal as signal
5
import pyvista as pv
6
from netCDF4 import Dataset
7

8
class XYZ2VTK():
9
    def __init__(self, lat1=25, lon1=99, lat2=33, lon2=106, val=0.1, maxdep=80):
10
        self.lat1 = lat1
11
        self.lat2 = lat2
12
        self.lon1 = lon1
13
        self.lon2 = lon2
14
        self.lats = np.arange(lat1, lat2, 0.1)
15
        self.lons = np.arange(lon1, lon2, 0.1)
16
        self.real_dep = np.arange(-maxdep, 2, 2)
17
        self.dep = self.real_dep / 10
18
        self.vs_file = 'ETvs.xyz'
19
        self.moho_file = '../ccp3d/good_moho_ET2020_repick.dat'
20
        self.topo_file = '../figures/topo_resample.grd'

对于深度范围需要做一些说明，我们需要的深度范围是0-80 km，但水平范围分别是7（lon2-lon1）和8（lat2-lat1）度，两者单位不匹配，但是如果将经纬度换算成大地坐标，深度范围又将远小于水平范围导致图像纵向很扁，难以达到我们想要的效果，这时有两种选择

第一种是将经纬度换算成大地坐标，然后将深度范围人为放大，这时水平方向单位是km。
第二种是直接将深度范围缩小10倍，这样水平范围是7和8，深度范围也是8，比例较合适。

这里经纬度更符合我们的习惯，因此选择类第二种方法。self.real_dep是实际的深度范围，self.dep是缩小10倍后绘图用的深度范围。

三维速度模型转VTK网格#

原始的三维速度模型是.xyz数据表文件，包括4列：经度、纬度、深度和速度大小，我们现将其转换为VTK文件，共分为以下几步：

将数据表文件网格化为Numpy矩阵#

1
    def read_xyz(self):
2
        lat, lon, dep, vs = np.loadtxt(self.vs_file, comments='>', unpack=True)
3
        dep /= -1
4
        new_lon, new_lat, new_dep = np.meshgrid(self.lons, self.lats, self.real_dep, indexing='ij')
5
        grid_data = griddata((lon, lat, dep), vs, (new_lon, new_lat, new_dep))
6
        return grid_data

lat, lon, dep, vs是数据文件每列的内容。
grid_data是网格后的的三维矩阵。

生成三维速度模型的VTK文件#

这一步将速度模型的三维矩阵转换为VTK文件并保存。

1
    def create_vtk(self):
2
        grid = pv.UniformGrid()
3
        grid_vs = self.read_xyz()
4
        grid.dimensions = grid_vs.shape
5
        grid.origin = (self.lon1, self.lat1, self.dep[0])
6
        grid.spacing = (0.1, 0.1, 0.2)
7
        grid.point_arrays["values"] = grid_vs.flatten(order="F")
8
        grid.save('ETvs.vtk')

首先通过pv.UniformGrid()定义一个用于可视化的三维网格。
grid_vs.shape、grid.origin和grid.spacing分别定义了网格的个数，起始点坐标和网格间隔。
最后将速度值赋给grid.point_arrays["values"]并保存至文件ETvs.vtk。

将Moho面起伏转换为VTK多边形#

Moho文件也是.xyz文件，包括经度、纬度和深度三列，Paraview中可以通过多边形数据展示界面起伏。

将Moho面起伏网格化为二维矩阵#

1
    def read_moho(self):
2
        lat, lon, moho, _ = np.loadtxt(self.moho_file, comments='>', unpack=True)
3
        idx = np.where(np.logical_not(np.isnan(moho)))[0]
4
        lat = lat[idx]
5
        lon = lon[idx]
6
        moho = -moho[idx]
7
        new_lon, new_lat = np.meshgrid(self.lons, self.lats)
8
        grid_data = griddata((lon, lat), moho, (new_lon, new_lat))
9
        grid_data = ndimage.gaussian_filter(grid_data, sigma=2)
10
        return new_lon, new_lat, grid_data

与速度模型类似，先从.xyz文件中读取Moho面深度数据，再将其网格话为二维矩阵。为了更好地可视化，我们将数据进行了高斯滤波，高斯系数为2。

将Moho面二维矩阵转换为多边形网格#

1
    def create_moho_vtk(self):
2
        xx, yy, moho = self.read_moho()
3
        points = np.empty([0, 3])
4
        for i, ix in enumerate(self.lons):
5
            for j, jy in enumerate(self.lats):
6
                points = np.vstack((points, [ix, jy, moho[j, i]]))
7
        grid = pv.PolyData(points)
8
        grid['values'] = points[:, 2]
9
        surf = grid.delaunay_2d()
10
        surf.save('moho_real.vtk')

首先将二维矩阵的Moho面起伏转写为三列数据点格式，分别为经度、纬度和Moho面深度。然后用pv.PolyData(points)定义一个包含这些数据点的多边形实例，再将数值定义为Moho面深度，这时多边形数据只在网格点中有数据，而网格点间的多边形内没有数据，这时我们用PyVista自带的三角剖分算法delaunay_2d，对多边形进行数据填充，最后保存该VTK文件。

制作地形起伏数据#

通过GMT生成地形起伏数据#

首先用GMT中grdcut截取研究区域的地形起伏。
由于通常地形数据精度较高，反映大范围的地形变化时可能不直观，因此我们用grdsample命令对生成的网格文件进行减采样。最后保存为netCDF格式文件。

将地形数据转换为多边形网格#

1
    def create_topo_vtk(self):
2
        ds = Dataset(self.topo_file)
3
        z = ds.variables['z'][:].data / 10
4
        points = np.empty([0, 3])
5
        for i, ix in enumerate(self.lons):
6
            for j, jy in enumerate(self.lats):
7
                points = np.vstack((points, [ix, jy, z[j, i]]))
8
        grid = pv.PolyData(points)
9
        grid['values'] = points[:, 2]
10
        surf = grid.delaunay_2d()
11
        surf.save('topo.vtk')

首先用netcdf4模块读取减采样后的地形数据文件，然后同Moho面起伏，将其转换为多边形网格并保存。

将以上三种数据的VTK文件导入Paraview，即可绘制他们之间的关系。

TODO#

用Pavaview取速度等值面
用Pavaview将断层投影至地形
用Pavaview绘制真实深度刻度