.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "gallery/visualization/prism_layer_pyvista.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        Click :ref:`here <sphx_glr_download_gallery_visualization_prism_layer_pyvista.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_gallery_visualization_prism_layer_pyvista.py:


Plot a prism layer in 3D
========================

The :func:`harmonica.prism_layer` allows to easily create a layer of prisms
whose top and bottom boundaries might drape a certain surface, like topography,
bathymetry or the Moho discontinuity. It returns a :class:`xarray.Dataset`
object with the horizontal coordinates of the center of each prism, their top
and bottom boundaries and their physical properties, like their density.
Through the ``prism_layer`` accessor (see
:class:`harmonica.DatasetAccessorPrismLayer`) we can call some methods for our
prism layer. For example, the
:meth:`harmonica.DatasetAccessorPrismLayer.gravity` method lets us compute the
gravitational fields of the layer on any set of observation points.
Another interesting method is the
:meth:`harmonica.DatasetAccessorPrismLayer.to_pyvista`: it converts the prism
layer into a :class:`pyvista.UnstructuredGrid` that could be easily plotted
through ``pyvista``.

In this example we will show how we can build a prism layer out of a topography
and bathymetry grid for South Africa and how we can visualize the layer as a 3D
plot using ``pyvista``.

.. GENERATED FROM PYTHON SOURCE LINES 31-68

.. code-block:: default

    import ensaio
    import pyproj
    import pyvista as pv
    import verde as vd
    import xarray as xr

    import harmonica as hm

    # Read South Africa topography
    fname = ensaio.fetch_southern_africa_topography(version=1)
    south_africa_topo = xr.load_dataset(fname)

    # Project the grid
    projection = pyproj.Proj(proj="merc", lat_ts=south_africa_topo.latitude.values.mean())
    south_africa_topo = vd.project_grid(south_africa_topo.topography, projection=projection)

    # Create a 2d array with the density of the prisms Points above the geoid will
    # have a density of 2670 kg/m^3 Points below the geoid will have a density
    # contrast equal to the difference between the density of the ocean and the
    # density of the upper crust: # 1000 kg/m^3 - 2900 kg/m^3
    density = south_africa_topo.copy()  # copy topography to a new xr.DataArray
    density.values[:] = 2670.0  # replace every value for the density of the topography
    # Change density values of ocean points
    density = density.where(south_africa_topo >= 0, 1000 - 2900)

    # Create layer of prisms
    prisms = hm.prism_layer(
        (south_africa_topo.easting, south_africa_topo.northing),
        surface=south_africa_topo,
        reference=0,
        properties={"density": density},
    )

    # Create a pyvista UnstructuredGrid from the prism layer
    pv_grid = prisms.prism_layer.to_pyvista()
    pv_grid






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <table><tr><th>Header</th><th>Data Arrays</th></tr><tr><td>
    <table>
    <tr><th>UnstructuredGrid</th><th>Information</th></tr>
    <tr><td>N Cells</td><td>1620522</td></tr>
    <tr><td>N Points</td><td>12964176</td></tr>
    <tr><td>X Bounds</td><td>1.091e+06, 3.375e+06</td></tr>
    <tr><td>Y Bounds</td><td>-3.841e+06, -1.642e+06</td></tr>
    <tr><td>Z Bounds</td><td>-5.634e+03, 3.358e+03</td></tr>
    <tr><td>N Arrays</td><td>1</td></tr>
    </table>

    </td><td>
    <table>
    <tr><th>Name</th><th>Field</th><th>Type</th><th>N Comp</th><th>Min</th><th>Max</th></tr>
    <tr><td><b>density</b></td><td>Cells</td><td>float64</td><td>1</td><td>-1.900e+03</td><td>2.670e+03</td></tr>
    </table>

    </td></tr> </table>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 69-93

.. code-block:: default


    # Plot with pyvista
    plotter = pv.Plotter(lighting="three_lights", window_size=(1000, 800))
    plotter.add_mesh(pv_grid, scalars="density")
    plotter.set_scale(zscale=75)  # exaggerate the vertical coordinate
    plotter.camera_position = "xz"
    plotter.camera.elevation = 20
    plotter.camera.azimuth = 35
    plotter.camera.zoom(1.2)

    # Add a ceiling light
    west, east, south, north = vd.get_region((prisms.easting, prisms.northing))
    easting_center, northing_center = (east + west) / 2, (north + south) / 2
    light = pv.Light(
        position=(easting_center, northing_center, 10e3),
        focal_point=(easting_center, northing_center, 0),
        intensity=0.3,
        light_type="scene light",  # the light doesn't move with the camera
        positional=False,  # the light comes from infinity
    )
    plotter.add_light(light)

    plotter.show_axes()
    plotter.show()



.. image-sg:: /gallery/visualization/images/sphx_glr_prism_layer_pyvista_001.png
   :alt: prism layer pyvista
   :srcset: /gallery/visualization/images/sphx_glr_prism_layer_pyvista_001.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** ( 1 minutes  22.671 seconds)


.. _sphx_glr_download_gallery_visualization_prism_layer_pyvista.py:


.. only :: html

 .. container:: sphx-glr-footer
    :class: sphx-glr-footer-example



  .. container:: sphx-glr-download sphx-glr-download-python

     :download:`Download Python source code: prism_layer_pyvista.py <prism_layer_pyvista.py>`



  .. container:: sphx-glr-download sphx-glr-download-jupyter

     :download:`Download Jupyter notebook: prism_layer_pyvista.ipynb <prism_layer_pyvista.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_