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

.. only:: html

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

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

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

.. _sphx_glr_gallery_forward_prisms_topo_gravity.py:


Gravitational effect of topography
==================================

One possible application of the :func:`harmonica.prism_layer` function is to
create a model of the terrain and compute its gravity effect. Here we will use
a regular grid of topographic and bathymetric heights for South Africa to
create a prisms layer that model the terrain with a density of 2670 kg/m^3 and
the ocean with a density contrast of -1900 kg/m^3 obtained as the difference
between the density of water (1000 kg/m^3) and the normal density of upper
crust (2900 kg/m^3). Then we will use :func:`harmonica.prism_gravity` to
compute the gravity effect of the model on a regular grid of observation
points.

.. GENERATED FROM PYTHON SOURCE LINES 21-79



.. image:: /gallery/forward/images/sphx_glr_prisms_topo_gravity_001.png
    :alt: Gravitational acceleration of the topography
    :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none

    /usr/share/miniconda3/envs/testing/lib/python3.8/site-packages/cartopy/mpl/geoaxes.py:1797: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
      result = matplotlib.axes.Axes.pcolormesh(self, *args, **kwargs)






|

.. code-block:: default

    import pyproj
    import numpy as np
    import verde as vd
    import harmonica as hm
    import matplotlib.pyplot as plt
    import cartopy.crs as ccrs


    # Read South Africa topography
    south_africa_topo = hm.datasets.fetch_south_africa_topography()

    # 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},
    )

    # Compute gravity field on a regular grid located at 4000m above the ellipsoid
    coordinates = vd.grid_coordinates(
        region=(12, 33, -35, -18), spacing=0.2, extra_coords=4000
    )
    easting, northing = projection(*coordinates[:2])
    coordinates_projected = (easting, northing, coordinates[-1])
    prisms_gravity = prisms.prism_layer.gravity(coordinates_projected, field="g_z")

    # Make a plot of the computed gravity
    plt.figure(figsize=(8, 8))
    ax = plt.axes(projection=ccrs.Mercator())
    maxabs = vd.maxabs(prisms_gravity)
    tmp = ax.pcolormesh(
        *coordinates[:2],
        prisms_gravity,
        vmin=-maxabs,
        vmax=maxabs,
        cmap="RdBu_r",
        transform=ccrs.PlateCarree()
    )
    ax.set_extent(vd.get_region(coordinates), crs=ccrs.PlateCarree())
    plt.title("Gravitational acceleration of the topography")
    plt.colorbar(
        tmp, label="mGal", orientation="horizontal", shrink=0.93, pad=0.01, aspect=50
    )
    plt.show()


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

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


.. _sphx_glr_download_gallery_forward_prisms_topo_gravity.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: prisms_topo_gravity.py <prisms_topo_gravity.py>`



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

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


.. only:: html

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

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