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

.. only:: html

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

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

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

.. _sphx_glr_gallery_forward_tesseroid_layer.py:


Layer of tesseroids
===================

.. GENERATED FROM PYTHON SOURCE LINES 11-68



.. image-sg:: /gallery/forward/images/sphx_glr_tesseroid_layer_001.png
   :alt: tesseroid layer
   :srcset: /gallery/forward/images/sphx_glr_tesseroid_layer_001.png
   :class: sphx-glr-single-img


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

 Out:

 .. code-block:: none

    <IPython.core.display.Image object>






|

.. code-block:: default

    import boule as bl
    import ensaio
    import numpy as np
    import pygmt
    import verde as vd
    import xarray as xr

    import harmonica as hm

    fname = ensaio.fetch_earth_topography(version=1)
    topo = xr.load_dataarray(fname)

    region = (-78, -53, -57, -20)
    topo = topo.sel(latitude=slice(*region[2:]), longitude=slice(*region[:2]))

    ellipsoid = bl.WGS84

    longitude, latitude = np.meshgrid(topo.longitude, topo.latitude)
    reference = ellipsoid.geocentric_radius(latitude)
    surface = topo + reference
    density = xr.where(topo > 0, 2670.0, 1040.0 - 2670.0)

    tesseroids = hm.tesseroid_layer(
        coordinates=(topo.longitude, topo.latitude),
        surface=surface,
        reference=reference,
        properties={"density": density},
    )

    # Create a regular grid of computation points located at 10km above reference
    grid_longitude, grid_latitude = vd.grid_coordinates(region=region, spacing=0.5)
    grid_radius = ellipsoid.geocentric_radius(grid_latitude) + 10e3
    grid_coords = (grid_longitude, grid_latitude, grid_radius)

    # Compute gravity field of tesseroids on a regular grid of observation points
    gravity = tesseroids.tesseroid_layer.gravity(grid_coords, field="g_z")
    gravity = vd.make_xarray_grid(
        grid_coords,
        gravity,
        data_names="g_z",
        dims=("latitude", "longitude"),
        extra_coords_names="radius",
    )

    # Plot gravity field
    fig = pygmt.Figure()
    maxabs = vd.maxabs(gravity.g_z)
    pygmt.makecpt(cmap="polar", series=(-maxabs, maxabs))
    fig.grdimage(
        gravity.g_z,
        projection="M15c",
        nan_transparent=True,
    )
    fig.basemap(frame=True)
    fig.colorbar(frame='af+l"Gravity [mGal]"', position="JCR")
    fig.coast(shorelines="0.5p,black", borders=["1/0.5p,black"])
    fig.show()


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

   **Total running time of the script:** ( 0 minutes  55.301 seconds)


.. _sphx_glr_download_gallery_forward_tesseroid_layer.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: tesseroid_layer.py <tesseroid_layer.py>`



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

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


.. only:: html

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

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