diff --git a/doc/_static/figures/airy-isostasy.svg b/doc/_static/figures/airy-isostasy.svg deleted file mode 100644 index 69a958f12..000000000 --- a/doc/_static/figures/airy-isostasy.svg +++ /dev/null @@ -1,915 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - image/svg+xml - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/doc/api/index.rst b/doc/api/index.rst index 2bde2afce..9347a2775 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -83,13 +83,6 @@ Isostatic Moho :toctree: generated/ isostatic_moho_airy - isostasy_airy (**DEPRECATED**) - -.. warning:: - - The :func:`harmonica.isostasy_airy` function will be deprecated in - Harmonica v0.6. Please use :func:`harmonica.isostatic_moho_airy` - instead. Input and Output ---------------- diff --git a/harmonica/__init__.py b/harmonica/__init__.py index 32c8e70a5..dd9203ba2 100644 --- a/harmonica/__init__.py +++ b/harmonica/__init__.py @@ -19,7 +19,7 @@ from .forward.tesseroid import tesseroid_gravity from .forward.tesseroid_layer import DatasetAccessorTesseroidLayer, tesseroid_layer from .gravity_corrections import bouguer_correction -from .isostasy import isostasy_airy, isostatic_moho_airy +from .isostasy import isostatic_moho_airy from .transformations import ( derivative_easting, derivative_northing, diff --git a/harmonica/isostasy.py b/harmonica/isostasy.py index 6565bc702..0a2adbea0 100644 --- a/harmonica/isostasy.py +++ b/harmonica/isostasy.py @@ -8,9 +8,6 @@ Function to calculate the thickness of the roots and antiroots assuming the Airy isostatic hypothesis. """ -import warnings - -import numpy as np import xarray as xr @@ -162,100 +159,3 @@ def isostatic_moho_airy( for name, (_, density) in layers.items(): moho.attrs[f"density_{name}"] = density return moho - - -def isostasy_airy( - topography, - density_crust=2.8e3, - density_mantle=3.3e3, - density_water=1e3, - reference_depth=30e3, -): - r""" - Calculate the isostatic Moho depth from topography using Airy's hypothesis. - - .. warning:: - - The :func:`harmonica.isostasy_airy` function will be deprecated in - Harmonica v0.6. Please use :func:`harmonica.isostatic_moho_airy` - instead. - - 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). This assumption is - usually - - .. figure:: ../../_static/figures/airy-isostasy.svg - :align: center - :width: 400px - - *Schematic of isostatic compensation following the Airy hypothesis.* - - The relationship between the topographic/bathymetric heights (:math:`h`) - and the root thickness (:math:`r`) is governed by mass balance relations - and can be found in classic textbooks like [TurcotteSchubert2014]_ and - [Hofmann-WellenhofMoritz2006]_. - - On the continents (positive topographic heights): - - .. math :: - - r = \frac{\rho_{c}}{\rho_m - \rho_{c}} h - - while on the oceans (negative topographic heights): - - .. math :: - r = \frac{\rho_{c} - \rho_w}{\rho_m - \rho_{c}} h - - in which :math:`h` is the topography/bathymetry, :math:`\rho_m` is the - density of the mantle, :math:`\rho_w` is the density of the water, and - :math:`\rho_{c}` is the density of the crust. - - The computed root thicknesses will be added to the given reference Moho - depth (:math:`H`) to arrive at the isostatic Moho depth. Use - ``reference_depth=0`` if you want the values of the root thicknesses - instead. - - Parameters - ---------- - topography : array or :class:`xarray.DataArray` - Topography height and bathymetry depth in meters. It is usually prudent - to use floating point values instead of integers to avoid integer - division errors. - density_crust : float - Density of the crust in :math:`kg/m^3`. - density_mantle : float - Mantle density in :math:`kg/m^3`. - density_water : float - Water density in :math:`kg/m^3`. - reference_depth : float - The reference Moho depth (:math:`H`) in meters. - - Returns - ------- - moho_depth : array or :class:`xarray.DataArray` - The isostatic Moho depth in meters. - - """ - warnings.warn( - "The harmonica.isostasy_airy function will be deprecated in " - + "Harmonica v0.6. Please use harmonica.isostatic_moho_airy " - + "instead.", - FutureWarning, - ) - # Need to cast to array to make sure numpy indexing works as expected for - # 1D DataArray topography - oceans = np.array(topography < 0) - continent = np.logical_not(oceans) - scale = np.full(topography.shape, np.nan, dtype="float") - scale[continent] = density_crust / (density_mantle - density_crust) - scale[oceans] = (density_crust - density_water) / (density_mantle - density_crust) - moho = topography * scale + reference_depth - if isinstance(moho, xr.DataArray): - moho.name = "moho_depth" - moho.attrs["isostasy"] = "Airy" - moho.attrs["density_crust"] = str(density_crust) - moho.attrs["density_mantle"] = str(density_mantle) - moho.attrs["density_water"] = str(density_water) - moho.attrs["reference_depth"] = str(reference_depth) - return moho diff --git a/harmonica/tests/test_isostasy_old.py b/harmonica/tests/test_isostasy_old.py deleted file mode 100644 index f97c86ccc..000000000 --- a/harmonica/tests/test_isostasy_old.py +++ /dev/null @@ -1,69 +0,0 @@ -# Copyright (c) 2018 The Harmonica Developers. -# Distributed under the terms of the BSD 3-Clause License. -# SPDX-License-Identifier: BSD-3-Clause -# -# This code is part of the Fatiando a Terra project (https://www.fatiando.org) -# -""" -Testing isostasy_airy function calculation (soon to be deprecated) -""" -import numpy as np -import numpy.testing as npt -import pytest -import xarray as xr - -from ..isostasy import isostasy_airy - - -def test_deprecation_warning(): - """ - Test if isostasy_airy raises a FutureWarning - """ - topography = np.zeros(20, dtype=np.float64) - with pytest.warns( - FutureWarning, match="The harmonica.isostasy_airy function will be deprecated" - ): - isostasy_airy(topography) - - -def test_isostasy_airy_zero_topography(): - "Root should be zero for zero topography" - topography = np.zeros(20, dtype=np.float64) - npt.assert_equal(isostasy_airy(topography, reference_depth=0), 0) - npt.assert_equal(isostasy_airy(topography, reference_depth=30e3), 30e3) - # Check that the shape of the topography is preserved - topography = np.zeros((20, 31), dtype=np.float64) - assert isostasy_airy(topography).shape == topography.shape - npt.assert_equal(isostasy_airy(topography, reference_depth=0), 0) - npt.assert_equal(isostasy_airy(topography, reference_depth=30e3), 30e3) - - -def test_isostasy_airy(): - "Use a simple integer topography to check the calculations" - topography = np.array([-2, -1, 0, 1, 2, 3]) - true_root = np.array([-0.5, -0.25, 0, 0.5, 1, 1.5]) - root = isostasy_airy( - topography, - density_crust=1, - density_mantle=3, - density_water=0.5, - reference_depth=0, - ) - npt.assert_equal(root, true_root) - - -def test_isostasy_airy_dataarray(): - "Pass in a DataArray and make sure things work" - topography = xr.DataArray( - np.array([-2, -1, 0, 1, 2, 3]), coords=(np.arange(6),), dims=["something"] - ) - true_root = np.array([-0.5, -0.25, 0, 0.5, 1, 1.5]) - root = isostasy_airy( - topography, - density_crust=1, - density_mantle=3, - density_water=0.5, - reference_depth=0, - ) - assert isinstance(root, xr.DataArray) - npt.assert_equal(root.values, true_root)