2.0 - PVGeo+Discretize.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PVGeo-Discretize\n",
    "\n",
    "This notebook demonstrates how to pair `PVGeo` and `discretize` for simple processing routines. \n",
    "\n",
    "This notebook is outlined into four sections:\n",
    "\n",
    "1. Introduction to PVGeo\n",
    "2. Overview of new VTK interface in `discretize`\n",
    "3. Pairing PVGeo and `discretize`\n",
    "4. Examples of PVGeo in ParaView"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import discretize\n",
    "import PVGeo\n",
    "import numpy as np\n",
    "import vista\n",
    "import vtk\n",
    "\n",
    "print('NumPy Version: %s' % np.__version__)\n",
    "print('PVGeo Version: %s' % PVGeo.__version__)\n",
    "print('vista Version: %s' % vista.__version__)\n",
    "print('vtk Version: %s' % vtk.VTK_VERSION)\n",
    "\n",
    "vista.set_plot_theme('document')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Learn about PVGeo\n",
    "\n",
    "To learn more about PVGeo, please refer to the first notebook: `1.0 - Welcome`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Discretize VTK Mixin\n",
    "\n",
    "This section demonstrates the VTK interface in `discretize` outlined in [SimPEG/discretize#114](https://github.com/simpeg/discretize/pull/114).\n",
    "\n",
    "Let's check out how the new VTK interface can be used on simple mesh objects to create a VTK data object ready for VTK and/or PVGeo processing routines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a simple TensorMesh\n",
    "h1 = np.linspace(.1, .5, 3)\n",
    "h2 = np.linspace(.1, .5, 5)\n",
    "h3 = np.linspace(.1, .5, 3)\n",
    "mesh = discretize.TensorMesh([h1, h2, h3])\n",
    "mesh.plotGrid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have a `TensorMesh` object, we can call the `toVTK()` method to yield the proper VTK data object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get a VTK data object\n",
    "grid = mesh.toVTK()\n",
    "grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An additional feature added in [SimPEG/discretize#114](https://github.com/simpeg/discretize/pull/114) is the ability to specify rotated reference frames for any given mesh object. Let's rotate our reference frame and then convert the mesh to a VTK data object. Note that we no longer have a `vtkRectilinearGrid` but a `vtkStructuredGrid` due to having that `TensorMesh` rotated off of the traditional reference frame (traditional being <1,0,0>, <0,1,0>, <0,0,1>)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defined a rotated reference frame\n",
    "mesh.axis_u = (1,-1,0)\n",
    "mesh.axis_v = (-1,-1,0)\n",
    "mesh.axis_w = (0,0,1)\n",
    "\n",
    "# Check that the referenc fram is valid\n",
    "mesh._validate_orientation()\n",
    "\n",
    "# At this time, the grid code in discretize is not updated to plot the rotated grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Yield the rotated vtkStructuredGrid\n",
    "grid_r = mesh.toVTK()\n",
    "grid_r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a rendering of these two meshes to demonstrate the rotation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note: you could also use vista.BackgroundPlotter() for interactivity\n",
    "p = vista.Plotter(notebook=True)\n",
    "p.add_mesh(grid, color='green', show_edges=True)\n",
    "p.add_mesh(grid_r, color='orange', show_edges=True)\n",
    "p.show_grid()\n",
    "p.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Pairing discretize+PVGeo\n",
    "\n",
    "In this example, we load a `discretize` 3D model and a topography surface that generally covers that model in space. We then use PVGeo to provide a boolean array to describe whether any given cell is above/below the topography surface for the mesh.\n",
    "\n",
    "This is a fairly simple example... we want you to focus less on the specific task of extracting the topography and more on the idea that `discretize` and `PVGeo` are able to talk to eachother and share their processing results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Data\n",
    "\n",
    "Here we load in some data we'd like to process: the [**Laguna del Maule Bouguer Gravity**](http://docs.simpeg.xyz/content/examples/04-grav/plot_laguna_del_maule_inversion.html#sphx-glr-content-examples-04-grav-plot-laguna-del-maule-inversion-py) example from the SimPEG docs.\n",
    "\n",
    "This data scene is was produced from the [Laguna del Maule Bouguer Gravity](http://docs.simpeg.xyz/content/examples/04-grav/plot_laguna_del_maule_inversion.html#sphx-glr-content-examples-04-grav-plot-laguna-del-maule-inversion-py) example provided by [Craig Miller](https://github.com/craigmillernz) (Maule volcanic field, Chile. Refer to Miller et al 2016 EPSL for full details.)\n",
    "\n",
    "> Miller, C. A., Williams-Jones, G., Fournier, D., & Witter, J. (2017). 3D gravity inversion and thermodynamic modelling reveal properties of shallow silicic magma reservoir beneath Laguna del Maule, Chile. Earth and Planetary Science Letters, 459, 14–27. https://doi.org/10.1016/j.epsl.2016.11.007\n",
    "\n",
    "The rendering below shows several data sets and a model integrated together:\n",
    "\n",
    "* Point Data: the Bouguer gravity anomalies\n",
    "* Topography Surface\n",
    "* Inverted Model: The model has been both sliced and thresholded for low values\n",
    "\n",
    "\n",
    "This rendering was created in ParaView using file I/O methods in PVGeo for UBC fomrats and general VTK filters available in ParaView. A ParaView state file is included in the data directory to recreate this scene.\n",
    "\n",
    "![scene](./images/craig-example.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!head data/Craig-Chile/LdM_topo.topo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!head data/Craig-Chile/craig_chile.msh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's load the data files using a mixture of `discretize` and `PVGeo` for now to demo how PVGeo and discretize can talk to eachother."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the TensorMesh and some already processed model data \n",
    "mesh = discretize.TensorMesh.readUBC('craig_chile.msh', directory='data/Craig-Chile')\n",
    "models = {'lpout': mesh.readModelUBC(fileName='Lpout.mod', directory='data/Craig-Chile')}\n",
    "mesh.plot_3d_slicer(v=models['lpout'], zslice=2350)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh.nC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Process the Data\n",
    "\n",
    "Now let's use the mesh from `discretize` and the topo surface in `PVGeo` to create a model array that describes whether or not any given cell in the model space is above/below the topogrpahy surface.\n",
    "\n",
    "Also, let's ignore the fact the given model data is already accounts for topography (NaN values). Let's suppose for a moment that you are designing/inspecting your model space: Simply load the topography into a `vtkPolyData` object in `PVGeo` and feed it to the `ExtractTopography` algorithm.\n",
    "\n",
    "Since the given topography file is in the [3D GIF Topography](https://giftoolscookbook.readthedocs.io/en/latest/content/fileFormats/topoGIF3Dfile.html) format, we can use the `PVGeo.ubc.TopoReader` file reader to read the file and automatically construct the `vtkPolyData`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load topography data using PVGeo\n",
    "topo = PVGeo.ubc.TopoReader().Apply('data/Craig-Chile/LdM_topo.topo')\n",
    "topo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Call the ExtractTopography algorithm and have it apply on the \n",
    "#    discretize mesh and the topography\n",
    "extracted = PVGeo.grids.ExtractTopography().Apply(mesh.toVTK(models), topo)\n",
    "extracted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "active = extracted.get_scalar('Extracted')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "models['active'] = active"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh.plot_3d_slicer(v=models['active'], zslice=2350)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What about plotting those orthographic slices in 3D and showing the topography and creating a whole integrated scene?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slices = extracted.slice_orthogonal()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# And then we can add more features to the scene\n",
    "p = vista.Plotter()\n",
    "p.add_mesh(topo, opacity=0.5, psize=1.0, \n",
    "           cmap='gist_earth', clim=[1.7e+03, 3.104e+03])\n",
    "p.add_mesh(slices, cmap='coolwarm')\n",
    "p.show_grid()\n",
    "p.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What about using a VTK algorthm?\n",
    "\n",
    "Easy! Simply pass the VTK object to the VTK algorithm and use PVGeo's top level functions to yield the output in NumPy or Pandas friendly data structures!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate your algorithm\n",
    "alg = vtk.vtkCellSizeFilter()\n",
    "# Set the inputs\n",
    "alg.SetInputDataObject(mesh.toVTK(models=models))\n",
    "# Run the algorithm\n",
    "alg.Update()\n",
    "# Yield the output on the 0th port\n",
    "out = alg.GetOutputDataObject(0)\n",
    "\n",
    "# Get the Volme array via PVGeo\n",
    "counts = PVGeo.getArray(out, 'Volume')\n",
    "# Use that arry to plot up the results!\n",
    "mesh.plot_3d_slicer(counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or better yet, utilize `vista`'s streamlined interface to VTK:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_sizes = mesh.toVTK(models=models).compute_cell_sizes()\n",
    "cell_sizes.plot(scalars='Volume')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Now What?\n",
    "\n",
    "There are tons of awesome algorithms in ParaView/VTK and PVGeo that you could use to integrate datasets and produce meaningful visualizations. At this time, using those algorithms in a standard Python environment doesn't really make sense... They should be used directly in ParaView or `pvpython` to create visualizations until we have a stable toolset for interactive VTK visualizations in Jupyter Notebooks.\n",
    "\n",
    "Note that the new VTK interface in `discretize` enables PVGeo to build an interface for `discretize` directly in [ParaView](https://paraview.org).\n",
    "\n",
    "Use the UBC suite in PVGeo which has I/O functionality for `discretize` to create compelling 3D visualizations of all your data!\n",
    "\n",
    "Some algorithms that might be of interest to you:\n",
    "\n",
    "- [Extract Topography](http://pvgeo.org/examples/grids/extract-topography/): Use a topography surface to add an active cells field to an input dataset\n",
    "- [Create Rectilinear Grid](http://pvgeo.org/examples/model-building/create-rectilinear-grid/) : Create a rectilinear grid / tensor mesh (`vtkRectilinearGrid`)\n",
    "- [The UBC suite in PVGeo](http://pvgeo.org/examples/contents/#ubc-mesh-tools): Contains file I/O for TensorMeshes, TreeMeshes, and time series model data as well as file readers for general data formats like Grav/Mag observations or topography surfaces. The status of the UBC suite can be found on [OpenGeoVis/PVGeo#28](https://github.com/OpenGeoVis/PVGeo/issues/28).\n",
    "- [Many Slices Along Axis](http://pvgeo.org/examples/filters-general/many-slices-along-axis/): Generate N slices of a dataset along a specified axis\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:dev]",
   "language": "python",
   "name": "dev"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}