Fix delta_pressure to handle the case when pressure levels are greate…

…r than the surface pressure (#571)

* formatting

* add to release notes

* fix type problem

* make relative pressure warning an error

* update release notes and change delta pressure warning to error

* add additional delta pressure test coverage

* formatting

docs/release-notes.rst

@@ -5,7 +5,7 @@
 Release Notes
 =============
 
-v2024.03.0 (Unreleased)
+v2024.04.0 (Unreleased)
 -----------------------
 
 Internal Changes
@@ -15,6 +15,11 @@ Internal Changes
 * Add M1 runners to CI by `Katelyn FitzGerald`_ in (:pr:`581`)
 * Reorganize dask compatibility tests by `Anissa Zacharias`_ in (:pr:`568`)
 
+Bug Fixes
+^^^^^^^^^
+* Fix ``delta_pressure`` to handle the case where pressure level(s) are greater than surface pressure by `Katelyn FitzGerald`_ in (:pr:`571`)
+
+
 v2024.02.0 (February 28, 2024)
 ------------------------------
 This release switches the package to use an implicit namespace and adds support

geocat/comp/meteorology.py

@@ -1401,13 +1401,17 @@ def _delta_pressure1D(pressure_lev, surface_pressure):
         The pressure layer thickness array. Shares dimensions and units of
         `pressure_lev`.
     """
+
+    # Leaving this here to lay groundwork for including it as a separate argument
     pressure_top = min(pressure_lev)
 
     # Safety checks
-    if pressure_top <= 0:
-        warnings.warn("'pressure_lev` values must all be positive.")
+    if pressure_top < 0:
+        raise ValueError(
+            "`pressure_lev` values must all be greater than or equal to 0.")
+
     if pressure_top > surface_pressure:
-        warnings.warn(
+        raise ValueError(
             "`surface_pressure` must be greater than minimum `pressure_lev` value."
         )
 
@@ -1417,15 +1421,25 @@ def _delta_pressure1D(pressure_lev, surface_pressure):
         pressure_lev = np.flip(pressure_lev)
 
     # Calculate delta pressure
-    delta_pressure = np.empty_like(pressure_lev)
+    delta_pressure = np.full_like(pressure_lev, np.nan)
+
+    [indices] = np.nonzero(np.array(pressure_lev) <= surface_pressure)
 
-    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])
+    start_level = min(indices)
+    end_level = max(indices)
+
+    delta_pressure[start_level] = (pressure_lev[start_level] + pressure_lev[
+        start_level + 1]) / 2 - pressure_top  # top level
+
+    delta_pressure[start_level + 1:end_level] = [
+        (a - b) / 2
+        for a, b in zip(pressure_lev[start_level + 2:end_level +
+                                     1], pressure_lev[start_level:end_level])
     ]
-    delta_pressure[-1] = surface_pressure - (
-        pressure_lev[-1] + pressure_lev[-2]) / 2  # bottom level
+
+    delta_pressure[end_level] = surface_pressure - (
+        pressure_lev[end_level] +
+        pressure_lev[end_level - 1]) / 2  # bottom level
 
     # Return delta_pressure to original order
     if is_pressuredecreasing:

test/test_meteorology.py

@@ -515,15 +515,29 @@ def test_delta_pressure1D(self) -> None:
         delta_p = delta_pressure(pressure_lev, self.surface_pressure_scalar)
         assert sum(delta_p) == (self.surface_pressure_scalar - pressure_top)
 
-    def test_negative_pressure_warning(self) -> None:
+    def test_delta_pressure_level_below_surface(self) -> None:
+        pressure_lev = [float(i) for i in self.pressure_lev]
+        surface_pressure_adjusted = 900.0
+        delta_p = delta_pressure(pressure_lev, surface_pressure_adjusted)
+        assert np.nansum(delta_p) == (surface_pressure_adjusted -
+                                      min(pressure_lev))
+
+    def test_delta_pressure_levels_below_surface(self) -> None:
+        pressure_lev = [float(i) for i in self.pressure_lev]
+        surface_pressure_adjusted = 50.0
+        delta_p = delta_pressure(pressure_lev, surface_pressure_adjusted)
+        assert np.nansum(delta_p) == (surface_pressure_adjusted -
+                                      min(pressure_lev))
+
+    def test_negative_pressure_error(self) -> None:
         pressure_lev_negative = self.pressure_lev.copy()
         pressure_lev_negative[0] = -5
-        with pytest.warns(UserWarning):
+        with pytest.raises(ValueError):
             delta_pressure(pressure_lev_negative, self.surface_pressure_scalar)
 
-    def test_relative_pressure_warning(self) -> None:
+    def test_relative_pressure_error(self) -> None:
         surface_pressure_low = 0.5
-        with pytest.warns(UserWarning):
+        with pytest.raises(ValueError):
             delta_pressure(self.pressure_lev, surface_pressure_low)
 
     def test_output_type(self) -> None: