PVGeo-Examples 1.0 - Welcome

14 篇文章 2 订阅
import PVGeo
import numpy as np
import pandas as pd
import pyvista

print('NumPy Version: %s' % np.__version__)
print('PVGeo Version: %s' % PVGeo.__version__)
print('pyvista Version: %s' % pyvista.__version__)
NumPy Version: 1.18.5
PVGeo Version: 2.1.0
pyvista Version: 0.24.2

Welcome to PVGeo

Thanks for checking out this notebook! We hope this provides some insight on how you can get started using PVGeo in your Python (3) routines. Let’s get started!

At the top of this notebook, we import numpy and PVGeo and display the current environment’s version of thos packages. Since PVGeo is still in its infancy, feature development is rapid and new versions get deployed very often. With this in mind, be sure to keep the PVGeo in your environment up to date!

A simple way to update PVGeo from your Jupyter Notebook:

!pip install --upgrade PVGeo

We have placed some data files in the data/ directory for you to use. The remainder of this notebook will explore a few examples loading, filtering, and writing out data using a combination of PVGeo and vista

1. Introduction to PVGeo

What is PVGeo?

  • Python package for 3D/4D geovisualization.
  • Create compelling and integrated visualizations.
  • Built upon VTK, a scalable and well-maintained visualization library.
  • Extends geovisualization into ParaView, VTK.js, and Virtual Reality.
  • Open-source and automatically deployed.

在这里插入图片描述

Abstract

PVGeo is an open-source Python package for geoscientific visualization and analysis, harnessing an already powerful software platform: the Visualization Toolkit (VTK) and its front-end application, ParaView. The VTK software platform is well-maintained, contains an expansive set of native functionality, and provides a robust foundation for scientific visualization, yet the development of tools compatible for geoscience data and models has been limited. As a software extension package to VTK and ParaView, PVGeo addresses the lack of geoscientific compatibility by creating a framework for geovisualization. This framework is a set of tools for visually integrating geoscience data and models directly within ParaView’s graphical user interface, simplifying the required routines to make compelling visualizations of geoscientific datasets. PVGeo aims to make the process of importing data into ParaView simple and fluid for users while providing a guide for contributions avoiding the typical, ambitious programming endeavor of building ParaView plugins. The PVGeo package is available for download on PyPI (pip install PVGeo), documented online, and open-source on GitHub for community-driven development.

PVGeo Resources

Take aways

  • Join PVGeo on Slack
    • The slack workspace is for anyone using ParaView for geovisualization
  • Presentation at AGU in December 2018
  • ParaView natively extends into VR (dynamically linked)
  • VTK and ParaView are incredibly scalable
  • PVGeo is Python based and open-source

Filtering Data with PVGeo

Table to Points

Let’s go ahead and load a simple file that has XYZ coordinates and a boolean array for fault presence. This point cloud makes some sort of regular grid, but we have forgotten the deatials of the cell spacings and local coordinate rotations.

We will read in this data with pandas and send it to the PVGeo filter PointsToPolyData to create a vista.PolyData object (essentially a point cloud).

import pandas as pd
points = pd.read_csv('data/fault_points.csv')
points[0:2]
XYZFault
0326819.4974407450.6361287.50
1326834.3404407470.7531287.50
# Convert to vista.PolyData (assumes first this three columns are XYZ)
vtkpoints = PVGeo.points_to_poly_data(points)
vtkpoints
HeaderData Arrays
PolyDataInformation
N Cells499200
N Points499200
X Bounds3.268e+05, 3.302e+05
Y Bounds4.406e+06, 4.410e+06
Z Bounds1.250e+01, 1.288e+03
N Arrays1
NameFieldTypeN CompMinMax
FaultPointsfloat6410.000e+001.000e+00

Note that we have a vista.PolyData object now which allows us to do all types of immediate plotting of our data. First, lets threshold our points as the point cloud has a bunch of zeros and ones throughout the dataspace to describe the presence of a fault.

To threshold the points, we call the threshold filter directly on our data object and pass the thresholding value.

We can then plot the result by calling the plot function. (Note: change the notebook parameter to False for an interactive window)

vtkpoints.plot(clim=[0,1])

在这里插入图片描述

Points to Voxelized Volume

The above figure is pretty cools! But its a point cloud which means out filtering options are pretty limited. Fortunately, we know that the point cloud represents some sort of regularlized gridded volume of data and PVGeo has a filter to recover that volume. This will allow further volumetric operations can be performed with other PVGeo or VTK filters.

Remember that these points are rotated and we do not know the cell sizes… this is okay! The VoxelizePoints filter from PVGeo will handle the recovory of the coordinate rotation and grid our data without running an interpolation scheme. The VoxelizePoints filter assumes that the points are structure on some rotated XY-plane with regular cell spacings and does the rest on its own! Check out VoxelizePoints code docs for more details.

# The full pipeline method
print('Voxelizing... ', end='')
voxelizer = PVGeo.filters.VoxelizePoints()
#print(dir(voxelizer))
grid = voxelizer.apply(vtkpoints)
print('done.')

# Output the results
print('Recovered Angle (deg.): %.3f' % voxelizer.get_recovered_angle())
print('Recovered Cell Sizes: (%.2f, %.2f, %.2f)' % voxelizer.get_spacing())
grid
Voxelizing... done.
Recovered Angle (deg.): 53.550
Recovered Cell Sizes: (25.00, 25.00, 25.00)
HeaderData Arrays
UnstructuredGridInformation
N Cells499200
N Points524064
X Bounds3.268e+05, 3.302e+05
Y Bounds4.406e+06, 4.410e+06
Z Bounds0.000e+00, 1.300e+03
N Arrays3
NameFieldTypeN CompMinMax
FaultCellsfloat6410.000e+001.000e+00
Recovered Angle (Deg.)Fieldsfloat6415.355e+015.355e+01
Recovered Cell SizesFieldsfloat6432.500e+012.500e+01
# A simpler method to voxelize
# grid = PVGeo.filters.VoxelizePoints().Apply(vtkpoints)
# grid

Now we have a volumetric dataset in the form of a vista.UnstructuredGrid! This means we can perform volumetric operations like slicing, clipping, and the works!

Slice Volumetric Data

Now lets use one of vista's filters to create slices of the thresholded dataset. Specifically, we are using the slice_orthogonal filter that will create 3 orthogonal slices through a data volume.

slices = grid.slice_orthogonal()
slices
InformationBlocks
MultiBlockValues
N Blocks3
X Bounds326804.336, 330185.260
Y Bounds4406253.954, 4409862.606
Z Bounds0.000, 1300.000
IndexNameType
0YZPolyData
1XZPolyData
2XYPolyData

Integrated Visualization

Up to this point, we have filtered a single dataset and plotted the result by itself; this is usefult, but what if we have all kinds of data we want to throw into one rendering environment? Its pretty easy:

For a simple case, see below. Otherwise, move on to our next notebooks that run through bigger datasets and show how to combine data feature like a 3D model with topography surfaces, well trajectories, and more!

clip = grid.clip(normal='x').clip(normal='-y').threshold(0.5)
p = pyvista.Plotter()
p.add_mesh(slices)
p.add_mesh(clip)
p.show_grid()
p.show()

在这里插入图片描述

Reference

PVGeo-Examples

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值