Damavand Volcano

Visualize 3D models of Damavand Volcano, Alborz, Iran.

This is an adaption of Alexey Pechnikov and A.V.Durandin’s ParaView-MoshaFault.

See LinkedIn posts for more details:

• The slices of the 3D model of the density on the Mosha fault area, North Iran

• Comparing Magnetic and Gravity Data to the Mosha Fault Area

• North Iran, Mosha fault

• North Iran

Originally posted: banesullivan/damavand-volcano

# sphinx_gallery_thumbnail_number = 6
import numpy as np
import pooch
import pyvista as pv
from pyvista import examples

base_url = r"https://raw.githubusercontent.com/pyvista/vtk-data/master/Data/{}"
a = pooch.retrieve(url=base_url.format("gebco7510_49cl.stl"), known_hash=None)
b = pooch.retrieve(url=base_url.format("gebco7510_55cl.stl"), known_hash=None)
c = pooch.retrieve(url=base_url.format("AOI.Damavand.32639.vtp"), known_hash=None)

gebco = examples.download_damavand_volcano()
gebco_a = pv.read(a)
gebco_b = pv.read(b)
aoi = pv.read(c)

opacity = [0, 0.75, 0, 0.75, 1.0]
clim = [0, 100]

p = pv.Plotter()
p.add_volume(
    gebco,
    cmap="magma",
    clim=clim,
    opacity=opacity,
    opacity_unit_distance=6000,
)
p.show()

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

p = pv.Plotter()
p.add_mesh(gebco.outline(), color="k")
p.add_mesh(voi, cmap="magma")
p.show()

p = pv.Plotter()
p.add_volume(voi, cmap="magma", clim=clim, opacity=opacity, opacity_unit_distance=2000)
p.camera_position = [
    (531554.5542909054, 3944331.800171338, 26563.04809259223),
    (599088.1433822059, 3982089.287834022, -11965.14728669936),
    (0.3738545892415734, 0.244312810377319, 0.8947312427698892),
]
p.show()

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

Header	Data Arrays
PolyData Information N Cells 49181 N Points 25586 N Strips 0 X Bounds 5.880e+05, 6.130e+05 Y Bounds 3.969e+06, 3.996e+06 Z Bounds -2.015e+04, 4.949e+03 N Arrays 1	Name Field Type N Comp Min Max data Points float32 1 5.000e+00 5.000e+01

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

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

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()