Adding method to calculate delta pressure (#338)

* first pass at calculating delta pressure

* linting

* much better! flattened arrays, warnings, and docstrings

* fix warnings (had reverse logic), preserve Xarray attributes

* Create test_deltapressure

* Update calc_deltapressure.py

* Update and rename test_deltapressure to test_deltapressure.py

* write tests

* pre-commit

* pre-commit

* fix import

* fix imports again

* test class inherits from TestCase

* wrong py name

* change plev to array of floats

* fix or remove some tests

* Update __init__.py

* missing sys import

* move to meteorology

* tests passing

exit()
:q
''

* add new section and PR functions to documentation

* Update src/geocat/comp/meteorology.py

Co-authored-by: Anissa Zacharias <[email protected]>

* Update src/geocat/comp/meteorology.py

Co-authored-by: Anissa Zacharias <[email protected]>

* rename methods and fix docstring autolinks

* remove conversion to numpy

* pre-commit

* attempt at docstring for dpres_plev wrapper

* dpres_plev docs

* attempt at fixing links and docstring

* Quick formatting fox

* make doc string variable hidden

* add line break to wrapper docstring

* Update src/geocat/comp/meteorology.py

Co-authored-by: Anissa Zacharias <[email protected]>

* slash going wrong way

* annotation in wrapper

* remove unnecessary checks for nan pressure or dims

* rm check for dims and nan

* Update test/test_meteorology.py

Co-authored-by: Heather Craker <[email protected]>

* Update src/geocat/comp/__init__.py

Co-authored-by: Heather Craker <[email protected]>

* Update test/test_meteorology.py

* rm test for 4d

---------

Co-authored-by: anissa111 <[email protected]>
Co-authored-by: Heather Craker <[email protected]>

docs/internal_api/index.rst

@@ -43,6 +43,8 @@ Internal API
 
    geocat.comp.interpolation._post_interp_multidim
 
+   geocat.comp.meteorology._delta_pressure1D
+
    geocat.comp.meteorology._dewtemp
 
    geocat.comp.meteorology._heat_index

docs/user_api/index.rst

@@ -58,6 +58,7 @@ Meteorology
    :toctree: ./generated/
 
    actual_saturation_vapor_pressure
+   delta_pressure
    dewtemp
    heat_index
    max_daylight
@@ -98,7 +99,6 @@ GeoCAT-comp routines from GeoCAT-f2py
    :nosignatures:
    :toctree: ./generated/
 
-   dpres_plevel
    grid_to_triple
    linint1
    linint2
@@ -109,6 +109,16 @@ GeoCAT-comp routines from GeoCAT-f2py
    rgrid2rcm
    triple_to_grid
 
+NCL Function Name Wrappers
+--------------------------
+.. currentmodule:: geocat.comp
+.. autosummary::
+   :nosignatures:
+   :toctree: ./generated/
+
+   meteorology.dpres_plev
+
+
 Deprecated Functions
 --------------------
 Climatologies

src/geocat/comp/__init__.py

@@ -8,7 +8,7 @@
 from .meteorology import (dewtemp, heat_index, relhum, relhum_ice, relhum_water,
                           actual_saturation_vapor_pressure, max_daylight,
                           psychrometric_constant, saturation_vapor_pressure,
-                          saturation_vapor_pressure_slope)
+                          saturation_vapor_pressure_slope, delta_pressure)
 from .spherical import decomposition, recomposition, scale_voronoi
 from .stats import eofunc, eofunc_eofs, eofunc_pcs, eofunc_ts, pearson_r
 # bring all functions from geocat.f2py into the geocat.comp namespace

src/geocat/comp/meteorology.py

@@ -1379,3 +1379,142 @@ def saturation_vapor_pressure_slope(
             tfill)
 
     return svp_slope
+
+
+def _delta_pressure1D(pressure_lev, surface_pressure):
+    """Helper function for `delta_pressure`. Calculates the pressure layer
+    thickness (delta pressure) of a one-dimensional pressure level array.
+
+    Returns an array of length matching `pressure_lev`.
+
+    Parameters
+    ----------
+    pressure_lev : :class:`numpy.ndarray`
+        The pressure level array. May be in ascending or descending order.
+        Must have the same units as `surface_pressure`.
+
+    surface_pressure : :class:`float`
+        The scalar surface pressure. Must have the same units as
+        `pressure_lev`.
+
+    Returns
+    -------
+    delta_pressure : :class:`numpy.ndarray`
+        The pressure layer thickness array. Shares dimensions and units of
+        `pressure_lev`.
+    """
+    pressure_top = min(pressure_lev)
+
+    # Safety checks
+    if pressure_top <= 0:
+        warnings.warn("'pressure_lev` values must all be positive.")
+    if pressure_top > surface_pressure:
+        warnings.warn(
+            "`surface_pressure` must be greater than minimum `pressure_lev` value."
+        )
+
+    # Sort so pressure increases (array goes from top of atmosphere to bottom)
+    is_pressuredecreasing = pressure_lev[1] < pressure_lev[0]
+    if is_pressuredecreasing:
+        pressure_lev = np.flip(pressure_lev)
+
+    # Calculate delta pressure
+    delta_pressure = np.empty_like(pressure_lev)
+
+    delta_pressure[0] = (pressure_lev[0] +
+                         pressure_lev[1]) / 2 - pressure_top  # top level
+    delta_pressure[1:-1] = [
+        (a - b) / 2 for a, b in zip(pressure_lev[2:], pressure_lev[:-1])
+    ]
+    delta_pressure[-1] = surface_pressure - (
+        pressure_lev[-1] + pressure_lev[-2]) / 2  # bottom level
+
+    # Return delta_pressure to original order
+    if is_pressuredecreasing:
+        delta_pressure = np.flip(delta_pressure)
+
+    return delta_pressure
+
+
+def delta_pressure(pressure_lev, surface_pressure):
+    """Calculates the pressure layer thickness (delta pressure) of a constant
+    pressure level coordinate system.
+
+    Returns an array of shape matching (`surface_pressure`, `pressure_lev`).
+
+    Parameters
+    ----------
+    pressure_lev : :class:`numpy.ndarray`, :class:'xarray.DataArray`
+        The pressure level array. May be in ascending or descending order.
+        Must have the same units as `surface_pressure`.
+    surface_pressure : :class:`np.Array`, :class:'xr.DataArray`
+        The scalar or N-dimensional surface pressure array. Must have the same
+        units as `pressure_lev`.
+
+    Returns
+    -------
+    delta_pressure : :class:`numpy.ndarray`, :class:'xarray.DataArray`
+        The pressure layer thickness array. Shares units with `pressure_lev`.
+        If `surface_pressure` is scalar, shares dimensions with
+        `pressure_level`. If `surface_pressure` is an array than the returned
+        array will have an additional dimension [e.g. (lat, lon, time) becomes
+        (lat, lon, time, lev)].
+
+    See Also
+    --------
+    Related NCL Functions:
+    `dpres_plev <https://www.ncl.ucar.edu/Document/Functions/Built-in/dpres_plevel.shtml>`__
+    """
+    # Get original array types
+    type_surface_pressure = type(
+        surface_pressure
+    )  # save type for delta_pressure to same type as surface_pressure at end
+    type_pressure_level = type(pressure_lev)
+
+    # Preserve attributes for Xarray
+    if type_surface_pressure == xr.DataArray:
+        da_coords = dict(surface_pressure.coords)
+        da_attrs = dict(surface_pressure.attrs)
+        da_dims = surface_pressure.dims
+    if type_pressure_level == xr.DataArray:
+        da_attrs = dict(
+            pressure_lev.attrs)  # Overwrite attributes to match pressure_lev
+
+    # Calculate delta pressure
+    if np.isscalar(surface_pressure):  # scalar case
+        delta_pressure = _delta_pressure1D(pressure_lev, surface_pressure)
+    else:  # multi-dimensional cases
+        shape = surface_pressure.shape
+        delta_pressure_shape = shape + (len(pressure_lev),
+                                       )  # preserve shape for reshaping
+
+        surface_pressure_flattened = np.ravel(
+            surface_pressure)  # flatten to avoid nested for loops
+        delta_pressure = [
+            _delta_pressure1D(pressure_lev, e)
+            for e in surface_pressure_flattened
+        ]
+
+        delta_pressure = np.array(delta_pressure).reshape(delta_pressure_shape)
+
+    # If passed in an Xarray array, return an Xarray array
+    # Change this to return a dataset that has both surface pressure and delta pressure?
+    if type_surface_pressure == xr.DataArray:
+        da_coords['lev'] = pressure_lev.values
+        da_dims = da_dims + ("lev",)
+        da_attrs.update({"long name": "pressure layer thickness"})
+        delta_pressure = xr.DataArray(delta_pressure,
+                                      coords=da_coords,
+                                      dims=da_dims,
+                                      attrs=da_attrs,
+                                      name="delta pressure")
+
+    return delta_pressure
+
+
+def dpres_plev(pressure_lev, surface_pressure):
+    return delta_pressure(pressure_lev, surface_pressure)
+
+
+_dpres_plev_doc_str = f".. attention:: This method is a wrapper for `delta_pressure <https://geocat-comp.readthedocs.io/en/stable/user_api/generated/geocat.comp.meteorology.delta_pressure.html>`_.\n\n    {delta_pressure.__doc__}"
+setattr(dpres_plev, '__doc__', _dpres_plev_doc_str)

test/test_meteorology.py

@@ -16,14 +16,14 @@
     from src.geocat.comp.meteorology import (
         dewtemp, heat_index, relhum, relhum_ice, relhum_water,
         actual_saturation_vapor_pressure, max_daylight, psychrometric_constant,
-        saturation_vapor_pressure, saturation_vapor_pressure_slope)
+        saturation_vapor_pressure, saturation_vapor_pressure_slope,
+        delta_pressure)
 else:
-    from geocat.comp.meteorology import (dewtemp, heat_index, relhum,
-                                         relhum_ice, relhum_water,
-                                         actual_saturation_vapor_pressure,
-                                         max_daylight, psychrometric_constant,
-                                         saturation_vapor_pressure,
-                                         saturation_vapor_pressure_slope)
+    from geocat.comp.meteorology import (
+        dewtemp, heat_index, relhum, relhum_ice, relhum_water,
+        actual_saturation_vapor_pressure, max_daylight, psychrometric_constant,
+        saturation_vapor_pressure, saturation_vapor_pressure_slope,
+        delta_pressure)
 
 
 class Test_dewtemp(unittest.TestCase):
@@ -645,3 +645,89 @@ def test_dask_lazy(self):
 
         assert isinstance((saturation_vapor_pressure_slope(tempf)).data,
                           dask.array.Array)
+
+
+class TestDeltaPressure(unittest.TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        cls.pressure_lev = np.array([1, 5, 100, 1000])
+        cls.pressure_lev_da = xr.DataArray(cls.pressure_lev)
+        cls.pressure_lev_da.attrs = {
+            "long name": "pressure level",
+            "units": "hPa",
+            "direction": "descending"
+        }
+
+        cls.surface_pressure_scalar = 1018
+        cls.surface_pressure_1D = np.array([1018, 1019])
+        cls.surface_pressure_2D = np.array([[1018, 1019], [1017, 1019.5]])
+        cls.surface_pressure_3D = np.array([[[1018, 1019], [1017, 1019.5]],
+                                            [[1019, 1020], [1018, 1020.5]]])
+
+        coords = {'time': [1, 2], 'lat': [3, 4], 'lon': [5, 6]}
+        dims = ["time", "lat", "lon"]
+        attrs = {"long name": "surface pressure", "units": "hPa"}
+        cls.surface_pressure_3D_da = xr.DataArray(cls.surface_pressure_3D,
+                                                  coords=coords,
+                                                  dims=dims,
+                                                  attrs=attrs)
+
+    def test_delta_pressure1D(self):
+        pressure_lev = [float(i) for i in self.pressure_lev]
+        pressure_top = min(pressure_lev)
+        delta_p = delta_pressure(pressure_lev, self.surface_pressure_scalar)
+        self.assertEqual(sum(delta_p),
+                         self.surface_pressure_scalar - pressure_top)
+
+    def test_negative_pressure_warning(self):
+        pressure_lev_negative = self.pressure_lev.copy()
+        pressure_lev_negative[0] = -5
+        with self.assertWarns(Warning):
+            delta_p = delta_pressure(pressure_lev_negative,
+                                     self.surface_pressure_scalar)
+
+    def test_relative_pressure_warning(self):
+        surface_pressure_low = 0.5
+        with self.assertWarns(Warning):
+            delta_p = delta_pressure(self.pressure_lev, surface_pressure_low)
+
+    def test_output_type(self):
+        delta_pressure_da = delta_pressure(self.pressure_lev_da,
+                                           self.surface_pressure_3D_da)
+        self.assertIsInstance(delta_pressure_da, xr.DataArray)
+
+        delta_pressure_np = delta_pressure(self.pressure_lev,
+                                           self.surface_pressure_3D)
+        self.assertIsInstance(delta_pressure_np, np.ndarray)
+
+    def test_output_dimensions(self):
+        delta_pressure_scalar = delta_pressure(self.pressure_lev,
+                                               self.surface_pressure_scalar)
+        self.assertEqual(delta_pressure_scalar.shape, (4,))
+
+        delta_pressure_1D = delta_pressure(self.pressure_lev,
+                                           self.surface_pressure_1D)
+        self.assertEqual(delta_pressure_1D.shape, (2, 4))
+
+        delta_pressure_2D = delta_pressure(self.pressure_lev,
+                                           self.surface_pressure_2D)
+        self.assertEqual(delta_pressure_2D.shape, (2, 2, 4))
+
+        delta_pressure_3D = delta_pressure(self.pressure_lev,
+                                           self.surface_pressure_3D)
+        self.assertEqual(delta_pressure_3D.shape, (2, 2, 2, 4))
+
+    def test_output_attrs(self):
+        delta_pressure_da = delta_pressure(self.pressure_lev_da,
+                                           self.surface_pressure_3D_da)
+        for item in self.pressure_lev_da.attrs:
+            self.assertIn(item, delta_pressure_da.attrs)
+
+    def test_output_coords(self):
+        delta_pressure_da = delta_pressure(self.pressure_lev_da,
+                                           self.surface_pressure_3D_da)
+        for item in self.surface_pressure_3D_da.coords:
+            self.assertIn(item, delta_pressure_da.coords)
+        for item in self.pressure_lev_da.coords:
+            self.assertIn(item, delta_pressure_da.coords)