Added depth_from_pressure method required for the calculation of `v…

…ertical_offset` value (#1207)

* added depth_from_pressure method

* Function arguments can be scalars or sequences; and add more tests (#2)

* Reverse order of equation terms to match UNESCO 1983 source

* For depth_from_pressure, accept scalars, lists and arrays for all 3 arguments. Include consistency checks.

* Add a greater variety of depth_from_pressure tests.

* Update echopype/utils/misc.py

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

---------

Co-authored-by: Emilio Mayorga <[email protected]>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

echopype/tests/utils/test_utils_misc.py

@@ -0,0 +1,31 @@
+import pytest
+
+import numpy as np
+from echopype.utils.misc import depth_from_pressure
+
+
+def test_depth_from_pressure():
+    # A single pressure value and defaults for the other arguments
+    pressure = 100.0
+    depth = depth_from_pressure(pressure)
+    assert np.isclose(depth, 99.2954)
+
+    # Array of pressure and list of latitude values
+    pressure = np.array([100.0, 101.0, 101.0])
+    latitude = [0.0, 30.0, 50.0]
+    depth = depth_from_pressure(pressure, latitude)
+    assert np.allclose(depth, [99.4265, 100.2881, 100.1096])
+
+    # Scalars specified for all 3 arguments
+    pressure = 1000.0
+    latitude = 0.0
+    atm_pres_surf = 10.1325  # standard atm pressure at sea level
+    depth = depth_from_pressure(pressure, latitude, atm_pres_surf)
+    assert np.isclose(depth, 982.0882)
+
+    # ValueError triggered by argument arrays having different lengths
+    pressure = np.array([100.0, 101.0, 101.0])
+    latitude = [0.0, 30.0]
+    with pytest.raises(ValueError) as excinfo:
+        depth = depth_from_pressure(pressure, latitude)
+    assert str(excinfo.value) == "Sequence shape or size does not match pressure"

echopype/utils/misc.py

@@ -1,3 +1,11 @@
+from typing import List, Optional, Tuple, Union
+
+import numpy as np
+from numpy.typing import NDArray
+
+FloatSequence = Union[List[float], Tuple[float], NDArray[float]]
+
+
 def camelcase2snakecase(camel_case_str):
     """
     Convert string from CamelCase to snake_case
@@ -11,3 +19,69 @@ def camelcase2snakecase(camel_case_str):
             camel_case_str = camel_case_str[:i] + "_" + camel_case_str[i:]
 
     return camel_case_str.lower()
+
+
+def depth_from_pressure(
+    pressure: Union[float, FloatSequence],
+    latitude: Optional[Union[float, FloatSequence]] = 30.0,
+    atm_pres_surf: Optional[Union[float, FloatSequence]] = 0.0,
+) -> NDArray[float]:
+    """
+    Convert pressure to depth using UNESCO 1983 algorithm.
+
+    UNESCO. 1983. Algorithms for computation of fundamental properties of seawater (Pressure to
+    Depth conversion, pages 25-27). Prepared by Fofonoff, N.P. and Millard, R.C. UNESCO technical
+    papers in marine science, 44. http://unesdoc.unesco.org/images/0005/000598/059832eb.pdf
+
+    Parameters
+    ----------
+    pressure : Union[float, FloatSequence]
+        Pressure in dbar
+    latitude : Union[float, FloatSequence], default=30.0
+        Latitude in decimal degrees.
+    atm_pres_surf : Union[float, FloatSequence], default=0.0
+        Atmospheric pressure at the surface in dbar.
+        Use the default 0.0 value if pressure is corrected to be 0 at the surface.
+        Otherwise, enter a correction for pressure due to air, sea ice and any other
+        medium that may be present
+
+    Returns
+    -------
+    depth : NDArray[float]
+        Depth in meters
+    """
+
+    def _as_nparray_check(v, check_vs_pressure=False):
+        """
+        Convert to np.array if not already a np.array.
+        Ensure latitude and atm_pres_surf are of the same size and shape as
+        pressure if they are not scalar.
+        """
+        v_array = np.array(v) if not isinstance(v, np.ndarray) else v
+        if check_vs_pressure:
+            if v_array.size != 1:
+                if v_array.size != pressure.size or v_array.shape != pressure.shape:
+                    raise ValueError("Sequence shape or size does not match pressure")
+        return v_array
+
+    pressure = _as_nparray_check(pressure)
+    latitude = _as_nparray_check(latitude, check_vs_pressure=True)
+    atm_pres_surf = _as_nparray_check(atm_pres_surf, check_vs_pressure=True)
+
+    # Constants
+    g = 9.780318
+    c1 = 9.72659
+    c2 = -2.2512e-5
+    c3 = 2.279e-10
+    c4 = -1.82e-15
+    k1 = 5.2788e-3
+    k2 = 2.36e-5
+    k3 = 1.092e-6
+
+    # Calculate depth
+    pressure = pressure - atm_pres_surf
+    depth_w_g = c1 * pressure + c2 * pressure**2 + c3 * pressure**3 + c4 * pressure**4
+    x = np.sin(np.deg2rad(latitude))
+    gravity = g * (1.0 + k1 * x**2 + k2 * x**4) + k3 * pressure
+    depth = depth_w_g / gravity
+    return depth