04_PVGeo数据接口与可视化
内容摘要： 地球物理数据格式多种多样，有时候你要用其它软件的处理结果或者把结果导入到第三方软件可视化，都需要进行格式转换，这时候你就要有数据接口的概念。可视化是数据分析的基础，重要性不需要多说。本节课针对地球物理数据接口和可视化问题，给大家介绍一个软件包PVGeo。
1、PVGeo介绍
PVGeo是基于VTK的为地球科学中的数据可视化量身定制的地球物理数据可视化工具包，使用python语言编写（主要依赖PyVista的底层库开发），侧重于2D或3D时变网格等结构化数据集。
最新版本2.0.1(版本还没到1.0的包，意味着多半还在原型开发中，后期变动会很大，建议大家尽量不要用)，该工具包Github的项目地址是https://github.com/OpenGeoVis/PVGeo。官方网站https://pvgeo.org/index.html。
工具包还提供了与Paraview的插件接口（Paraview Macros），可以直接使用Paraview作为前端操作界面，实现多种地球物理数据的导入和格式转换等操作（Paraview需要特定版本支持，需要使用python2.7环境）。
下图是该工具包工作原理。
图1 PVGeo工作原理
2、函数用法
PVGeo主要实现了Grid网格类数据的读写操作，对ESRI和Surfer的网格格式支持较好，可是实现诸如用该类网格描述的地形和异常的读写操作。
第二类是GSLib的地质建模结果的三维可视化，包括点三维网格和Table的读取和可视化，在地球物理反演方面支持UBC数据格式的可视化，包括观测数据，模型和网格的可视化能力。
在地球物理反演方面，如果采用以上的网格和模型描述格式来保存数据，后续的三维可视化可以直接通过PVGeo工具包来实现，现在PVGeo包含的各Suites类数量统计如下：
| base | filter | macro | reader | source | writer |
|---|---|---|---|---|---|
| 12 | 33 | 8 | 17 | 5 | 6 |
除了读取数据可视化之外，PVGeo也能提供简单的建模能力，比如构建一些规则的网格，从地形抽取数据构建六面体网格等。PVGeo包含了各种的功能模块，统称为suites，实现了读写、过滤等数据操作。
PVGeo支持模型数据读入的基本数据类型包括三类：重力观测数据文件，网格化文件，模型网格和节点。这三种类型都支持print方法看实例信息。
- （1）观测数据文件：该文件类型为pyvista.core.pointset.PolyData，是用于存储三维点坐标的数据类，访问该类实例的方法和属性需要注意。
instance.points #三维xyz坐标值
instance.point_arrays[‘Grav’]
instance.point_arrays[‘Err’]
在PVGeo中的实例化方法：
grav = PVGeo.ubc.GravObsReader().apply(grav_file_name)
- （2）网格化文件：该文件类型为pyvista.core.grid.UniformGrid，是用于存储网格化后的数据类，成功读取数据后，要访问这些数组，通过以下方法：
instance.point_arrays['Data'] #每个网格节点的数据
instance.points #三列数组，对应每个点的x，y，z
surfer的grd为二维第三列全为零，也可以通过instance.x/y/z
在PVGeo中的实例化方法：
grd = PVGeo.grids.SurferGridReader().apply(grd_file_name)
- （3）模型网格和节点文件：该文件类型为pyvista.core.grid.RectilinearGrid，是用于表示三维网格模型和属性的数据类，根据模型定义数据文件和属性数据文件，可以实例化该类，在类实例对象中，数据访问方法：
instance.points #返回一个模型网格的全部节点数据
instance.point_data_to_cell_data(True) #查看模型参数信息
instance.x/y/z #返回每个模型节点的坐标
在PVGeo中的实例化方法：
reader = PVGeo.ubc.TensorMeshReader()
reader.set_mesh_filename(mesh_file)
reader.add_model_file_name(model_file)
mesh = reader.apply()
3、 实例计算
import PVGeo
import geoist.others.fetch_data as fetch_data
url = fetch_data._get_gitee_file_url('geodataset','demodata/LdM_grav_obs.grv')
grav_file = fetch_data._retrieve_file(url, 'ttt.grv')
print(grav_file)
grav = PVGeo.ubc.GravObsReader().apply(grav_file)
print(grav) #pyvista.core.pointset.PolyData
# np.array(grav.points) #XYZ坐标点
# grav.point_arrays['Grav'] #重力数据
# grav.point_arrays['Err'] #误差数据
grav.plot(render_points_as_spheres=True, point_size=10)
#读取Surfer格式文件
url2 = fetch_data._get_gitee_file_url('geodataset','demodata/surfer-grid.grd')
grd_file = fetch_data._retrieve_file(url2, 'test.grd')
grd = PVGeo.grids.SurferGridReader().apply(grd_file)
# grd是pyvist的Uniform Grid数据类型，常用的属性：dimension, origin, spaceing
# grd.X, .Y, .Z可以范围数据点坐标，point_arrays是字典型数据存储
# 通过['Data']键值可以得到读入的numpy数组
max(grd.point_arrays['Data']) # 访问网格数据
warped = grd.warp_by_scalar(scale_factor=300.)
warped.plot(cmap='terrain')
#读取ubc的msh和mod文件
url3 = fetch_data._get_gitee_file_url('geodataset','pvgeo/craig_chile.msh')
mesh_file = fetch_data._retrieve_file(url3, 'craig_chile.msh')
# gitee 大于1m文件需要登陆
url4 = fetch_data._get_coding_file_url('geodataset','geodataset','Lpout.mod')
print(url4)
model_file = fetch_data._retrieve_file(url4, 'Lpout.mod')
reader = PVGeo.ubc.TensorMeshReader()
reader.set_mesh_filename(mesh_file)
reader.add_model_file_name(model_file)
mesh = reader.apply()
print(mesh)
type(mesh) #pyvista.core.grid.RectilinearGrid
mesh.points #返回一个网格节点属性的数组 XYZ分别对应坐标
mesh.threshold().plot()
mesh.slice_orthogonal().plot()
mesh.threshold([-0.6, -0.3]).plot(clim=[-0.6, 0.3])
4、注意事项
Pyvista是一套继承自VTK数据格式的python工具包，由于同VTK数据对象兼容，所以，PVGeo为实现VTK支持的三维可视化在其基础上开发。
VTK是一套广泛使用的可视化库，在python环境中，可以通过pip install vtk来安装，VTK的核心是C++开发的，该包实现了python接口。
