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

.. only:: html

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

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

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

.. _sphx_glr_gallery_isostatic_moho_airy.py:


Airy Isostatic Moho
===================

According to the Airy hypothesis of isostasy, topography above sea level is
supported by a thickening of the crust (a root) while oceanic basins are
supported by a thinning of the crust (an anti-root). Function
:func:`harmonica.isostatic_moho_airy` computes the depth to crust-mantle
interface (the Moho) according to Airy isostasy. The function takes the depth
to the crystalline basement and optionally any layers on top of it. Each layer
is defined by its thickness and its density. In addition, one must assume
a value for the reference thickness of the continental crust in order to
convert the root/anti-root thickness into Moho depth. The function contains
common default values for the reference thickness and crust, mantle
[TurcotteSchubert2014]_.

We'll use our sample topography data
(:func:`harmonica.datasets.fetch_topography_earth`) to calculate the Airy
isostatic Moho depth of Africa.

.. GENERATED FROM PYTHON SOURCE LINES 27-78



.. image-sg:: /gallery/images/sphx_glr_isostatic_moho_airy_001.png
   :alt: isostatic moho airy
   :srcset: /gallery/images/sphx_glr_isostatic_moho_airy_001.png
   :class: sphx-glr-single-img


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

 Out:

 .. code-block:: none

    Topography/bathymetry grid:
    <xarray.Dataset>
    Dimensions:     (longitude: 481, latitude: 511)
    Coordinates:
      * longitude   (longitude) float64 -20.0 -19.83 -19.67 ... 59.67 59.83 60.0
      * latitude    (latitude) float64 -40.0 -39.83 -39.67 ... 44.67 44.83 45.0
    Data variables:
        topography  (latitude, longitude) float32 -3.528e+03 -3.51e+03 ... 31.0 33.0

    Moho depth grid:
    <xarray.DataArray 'moho_depth' (latitude: 511, longitude: 481)>
    array([[17510.88  , 17574.602 , 17712.66  , ..., 11715.9   , 11602.619 ,
            11443.32  ],
           [17833.02  , 17942.76  , 18165.781 , ..., 11609.699 , 11581.381 ,
            11461.02  ],
           [17684.34  , 17797.62  , 18034.8   , ..., 11503.502 , 11570.762 ,
            11521.201 ],
           ...,
           [15167.399 , 15202.8   , 14937.301 , ..., 30151.2   , 30184.8   ,
            30184.8   ],
           [15347.939 , 15503.7   , 15266.5205, ..., 30145.6   , 30184.8   ,
            30184.8   ],
           [15362.101 , 15517.86  , 15252.359 , ..., 30151.2   , 30173.6   ,
            30184.8   ]], dtype=float32)
    Coordinates:
      * longitude  (longitude) float64 -20.0 -19.83 -19.67 ... 59.67 59.83 60.0
      * latitude   (latitude) float64 -40.0 -39.83 -39.67 -39.5 ... 44.67 44.83 45.0
    Attributes:
        isostasy:        Airy
        density_crust:   2800.0
        density_mantle:  3300.0
        density_water:   1030
    <IPython.core.display.Image object>






|

.. code-block:: default

    import ensaio
    import numpy as np
    import pygmt
    import xarray as xr

    import harmonica as hm

    # Load the elevation model and cut out the portion of the data corresponding to
    # Africa
    fname = ensaio.fetch_earth_topography(version=1)
    data = xr.load_dataset(fname)
    region = (-20, 60, -40, 45)
    data_africa = data.sel(latitude=slice(*region[2:]), longitude=slice(*region[:2]))
    print("Topography/bathymetry grid:")
    print(data_africa)

    # Calculate the water thickness
    oceans = np.array(data_africa.topography < 0)
    water_thickness = data_africa.topography * oceans * -1
    water_density = 1030

    # Calculate the isostatic Moho depth using the default values for densities and
    # reference Moho with water load. We neglect the effect of sediment here, so
    # basement elevation refers to topography.
    moho = hm.isostatic_moho_airy(
        basement=data_africa.topography,
        layers={"water": (water_thickness, water_density)},
    )
    print("\nMoho depth grid:")
    print(moho)

    # Draw the maps
    fig = pygmt.Figure()

    pygmt.grd2cpt(grid=moho, cmap="viridis", reverse=True, continuous=True)

    title = "Airy isostatic Moho depth of Africa"

    fig.grdimage(
        region=region,
        projection="Y20/0/10c",
        frame=["ag", f"+t{title}"],
        grid=moho,
        cmap=True,
    )

    fig.coast(shorelines="0.5p,black", resolution="crude")

    fig.colorbar(cmap=True, frame=["a10000f2500", "x+lmeters"])

    fig.show()


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

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


.. _sphx_glr_download_gallery_isostatic_moho_airy.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: isostatic_moho_airy.py <isostatic_moho_airy.py>`



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

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


.. only:: html

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

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