#!/usr/bin/env python
# coding: utf-8

# # PyVista Damavand Volcano
# 
# [PyVista](https://www.pyvista.org) is leveraged to quickly visualize 3D models of Damavand Volcano, Alborz, Iran.
# 
# Find this notebook [on GitHub](https://github.com/banesullivan/damavand-volcano).
# 
# Creator of this notebook: [Bane Sullivan](https://banesullivan.com)
# 
# This is an adaption of [Alexey Pechnikov](https://orcid.org/0000-0001-9626-8615) and [A.V.Durandin](https://orcid.org/0000-0001-6468-9757)'s [ParaView-MoshaFault](https://github.com/mobigroup/ParaView-MoshaFault)
# 
# See LinkedIn posts for more details:
# 
# - [The slices of the 3D model of the density on the Mosha fault area, North Iran](https://www.linkedin.com/posts/activity-6610080454911631360-97-V/)
# 
# - [Comparing Magnetic and Gravity Data to the Mosha Fault Area](https://www.linkedin.com/posts/activity-6609736436344201216-Kxls/)
# 
# - [North Iran, Mosha fault](https://www.linkedin.com/posts/activity-6609681862937853952-2BPG/)
# 
# - [North Iran](https://www.linkedin.com/posts/activity-6609486793676996608-ZF-J/)
# 
# ## Installation
# 
# Simply install PyVista (see [installation](https://docs.pyvista.org/getting-started/installation.html) for tips):
# 
# ```bash
# pip install -U pyvista panel
# ```
# 
# or 
# 
# ```bash
# conda install -c conda-forge pyvista panel
# ```
# 
# The `panel` installation is optional. `panel` enables embeddable renderings within the notebook (it is an experimental feature).

# In[1]:


import pyvista as pv
import numpy as np


# In[2]:


gebco = pv.read("data/GEBCO_2019.subset.32639.0.5km_b_gamma15.0km,20.0km.anomaly.vtk")
gebco_a = pv.read("data/gebco7510_49cl.stl")
gebco_b = pv.read("data/gebco7510_55cl.stl")
aoi = pv.read("data/AOI.Damavand.32639.vtp")


# In[3]:


def extract_subset(volume, voi):
    """Helper for filter not yet in PyVista.
    
    Coming in https://github.com/pyvista/pyvista/pull/569
    """
    import vtk
    alg = vtk.vtkExtractVOI()
    alg.SetVOI(voi)
    alg.SetInputDataObject(volume)
    alg.Update()
    return pv.wrap(alg.GetOutput())

gebco_clipped = extract_subset(gebco, [175, 200, 105, 132, 98, 170])


# In[4]:


contours = gebco_clipped.contour(np.arange(5,55,5))
contours


# In[5]:


contours.plot(cmap="nipy_spectral", opacity=0.15, use_panel=True)


# In[6]:


roi = [*gebco_clipped.bounds[0:4],
       *aoi.bounds[4:6]]
aoi_clipped = aoi.clip_box(roi, invert=False)
pv.plot([aoi, pv.Box(roi).outline()], cpos="xy", use_panel=True)


# In[7]:


p = pv.Plotter(window_size=np.array([1024, 768])*2)

# Add all the data we want to see
p.add_mesh(contours, cmap="nipy_spectral", opacity=0.15)
p.add_mesh(gebco_a, color="#ff0000")
p.add_mesh(gebco_b, color="#ff0000")
p.add_mesh(aoi_clipped, cmap="coolwarm", opacity=0.7)

# Add a title
p.add_text("Vent and Magma Chambers\nDamavand Volcano, Alborz")

# A nice perspective
p.camera_position = [(544065.5831913119, 3924518.576093113, 24324.3096344195),
                     (597885.1732914157, 3982998.0900773173, -12587.537450058662),
                     (0.33162714740718435, 0.26609487244915314, 0.9051060456978746)]
p.show(use_panel=True, screenshot="volcano.png")


# In[ ]: