Control zero humidity in thermo computations

Author: sandorkertesz

src/earthkit/meteo/thermo/array/thermo.py

@@ -13,6 +13,10 @@
 from earthkit.meteo import constants
 
 
+def _valid_number(x):
+    return x is not None and not np.isnan(x)
+
+
 def celsius_to_kelvin(t):
     """Converts temperature values from Celsius to Kelvin.
 
@@ -259,9 +263,21 @@ def compute_slope(self, t, phase):
         elif phase == "ice":
             return self._es_ice_slope(t)
 
-    def t_from_es(self, es):
-        v = np.log(es / self.c1)
-        return (v * self.c4w - self.c3w * self.t0) / (v - self.c3w)
+    def t_from_es(self, es, eps=1e-8, out=None):
+        def _comp(x):
+            x = x / self.c1
+            x = np.log(x)
+            return (x * self.c4w - self.c3w * self.t0) / (x - self.c3w)
+
+        if out is not None:
+            v = np.asarray(es)
+            z_mask = v < eps
+            v_mask = ~z_mask
+            v[v_mask] = _comp(v[v_mask])
+            v[z_mask] = out
+            return v
+        else:
+            return _comp(es)
 
     def _es_water(self, t):
         return self.c1 * np.exp(self.c3w * (t - self.t0) / (t - self.c4w))
@@ -556,13 +572,21 @@ def saturation_specific_humidity_slope(t, p, es=None, es_slope=None, phase="mixe
     return constants.epsilon * es_slope * p / v
 
 
-def temperature_from_saturation_vapour_pressure(es):
-    r"""Computes the temperature from saturation vapour pressure.
+def temperature_from_saturation_vapour_pressure(es, eps=1e-8, out=None):
+    r"""Compute the temperature from saturation vapour pressure.
 
     Parameters
     ----------
     es: ndarray
         :func:`saturation_vapour_pressure` (Pa)
+    eps: number
+        If ``out`` is not None, return ``out`` when ``es`` < ``eps``. If out
+        is None, ``eps`` is ignored and return np.nan for ``es`` values very
+        close to zero.
+    out: number or None
+        If not None, return ``out`` when ``es`` < ``eps``. If None, ``eps`` is
+        ignored and return np.nan for ``es`` values very close to zero.
+
 
     Returns
     -------
@@ -575,7 +599,7 @@ def temperature_from_saturation_vapour_pressure(es):
     phase ``es`` was computed to.
 
     """
-    return _EsComp().t_from_es(es)
+    return _EsComp().t_from_es(es, eps=eps, out=out)
 
 
 def relative_humidity_from_dewpoint(t, td):
@@ -750,20 +774,27 @@ def specific_humidity_from_relative_humidity(t, r, p):
     return specific_humidity_from_vapour_pressure(e, p)
 
 
-def dewpoint_from_relative_humidity(t, r):
-    r"""Computes the dewpoint temperature from relative humidity.
+def dewpoint_from_relative_humidity(t, r, eps=1e-8, out=None):
+    r"""Compute the dewpoint temperature from relative humidity.
 
     Parameters
     ----------
     t: ndarray
         Temperature (K)
     r: ndarray
         Relative humidity (%)
+    eps: number
+        If ``out`` is not None, return ``out`` when ``r`` < ``eps``.
+        If out is None, ``eps`` is ignored and return np.nan for ``q``
+        values very close to zero.
+    out: number or None
+        If not None, return ``out`` when ``r`` < ``eps``. If None, ``eps`` is
+        ignored and return np.nan for ``r`` values very close to zero.
 
     Returns
     -------
     ndarray
-        Dewpoint (K)
+        Dewpoint temperature (K)
 
 
     The computation starts with determining the the saturation vapour pressure over
@@ -782,11 +813,23 @@ def dewpoint_from_relative_humidity(t, r):
     equations used in :func:`saturation_vapour_pressure`.
 
     """
+    # by masking upfront we avoid RuntimeWarning when calling log() in
+    # the computation of td when r is very small
+    if out is not None:
+        r = np.asarray(r)
+        mask = r < eps
+        r[mask] = np.nan
+
     es = saturation_vapour_pressure(t, phase="water") * r / 100.0
-    return temperature_from_saturation_vapour_pressure(es)
+    td = temperature_from_saturation_vapour_pressure(es)
 
+    if out is not None and not np.isnan(out):
+        td[mask] = out
 
-def dewpoint_from_specific_humidity(q, p):
+    return td
+
+
+def dewpoint_from_specific_humidity(q, p, eps=1e-8, out=None):
     r"""Computes the dewpoint temperature from specific humidity.
 
     Parameters
@@ -795,11 +838,20 @@ def dewpoint_from_specific_humidity(q, p):
         Specific humidity (kg/kg)
     p: ndarray
         Pressure (Pa)
+    eps: number
+        As an intermediate result the saturation vapour pressure
+        (``es``) is computed (see details below). If ``out`` is not
+        None, return ``out`` when ``es`` < ``eps``. If out is None,
+        ``eps`` is ignored and return np.nan for ``es`` values very
+        close to zero.
+    out: number or None
+        If not None, return ``out`` when ``es`` < ``eps``. If None, ``eps`` is
+        ignored and return np.nan for ``es`` values very close to zero.
 
     Returns
     -------
     ndarray
-        Dewpoint (K)
+        Dewpoint temperature (K)
 
 
     The computation starts with determining the the saturation vapour pressure over
@@ -819,7 +871,20 @@ def dewpoint_from_specific_humidity(q, p):
     used in :func:`saturation_vapour_pressure`.
 
     """
-    return temperature_from_saturation_vapour_pressure(vapour_pressure_from_specific_humidity(q, p))
+    # by masking upfront we avoid RuntimeWarning when calling log() in
+    # the computation of td when q is very small
+    if out is not None:
+        q = np.asarray(q)
+        mask = q < eps
+        q[mask] = np.nan
+
+    es = vapour_pressure_from_specific_humidity(q, p)
+    td = temperature_from_saturation_vapour_pressure(es)
+
+    if out is not None and not np.isnan(out):
+        td[mask] = out
+
+    return td
 
 
 def virtual_temperature(t, q):
@@ -828,7 +893,7 @@ def virtual_temperature(t, q):
     Parameters
     ----------
     t: number or ndarray
-        Temperature (K)
+        Temperature (K)s
     q: number or ndarray
         Specific humidity (kg/kg)
 

tests/thermo/test_thermo_array.py

@@ -10,6 +10,7 @@
 import os
 
 import numpy as np
+import pytest
 
 from earthkit.meteo import thermo
 
@@ -320,7 +321,7 @@ def test_saturation_specific_humidity_slope():
         np.testing.assert_allclose(svp, v_ref[i])
 
 
-def test_temperature_from_saturation_vapour_pressure():
+def test_temperature_from_saturation_vapour_pressure_1():
     ref_file = "sat_vp.csv"
     d = np.genfromtxt(
         data_file(ref_file),
@@ -329,7 +330,32 @@ def test_temperature_from_saturation_vapour_pressure():
     )
 
     t = thermo.array.temperature_from_saturation_vapour_pressure(d["water"])
-    np.testing.assert_allclose(t, d["t"])
+    assert np.allclose(t, d["t"], equal_nan=True)
+
+
+@pytest.mark.parametrize(
+    "es,kwargs,expected_values",
+    [
+        (4.2, {}, 219.7796336743947),
+        ([4.2, 0, 0.001, np.nan], {"eps": 1e-2, "out": 100}, [219.7796336743947, 100, 100, np.nan]),
+        ([4.2, 0, 0.001, np.nan], {"eps": 1e-2, "out": np.nan}, [219.7796336743947, np.nan, np.nan, np.nan]),
+        (0, {}, np.nan),
+        (0.001, {"eps": 1e-2, "out": 100}, 100.0),
+        (0.001, {"eps": 1e-2, "out": np.nan}, np.nan),
+    ],
+)
+def test_temperature_from_saturation_vapour_pressure_2(es, kwargs, expected_values):
+
+    multi = isinstance(es, list)
+    if multi:
+        es = np.array(es)
+        expected_values = np.array(expected_values)
+
+    t = thermo.array.temperature_from_saturation_vapour_pressure(es, **kwargs)
+    if multi:
+        np.testing.assert_allclose(t, expected_values, equal_nan=True)
+    else:
+        assert np.isclose(t, expected_values, equal_nan=True)
 
 
 def test_relative_humidity_from_dewpoint():
@@ -421,43 +447,112 @@ def test_specific_humidity_from_relative_humidity():
     np.testing.assert_allclose(q, v_ref)
 
 
-def test_dewpoint_from_relative_humidity():
-    t = thermo.array.celsius_to_kelvin(np.array([20.0, 20, 0, 35, 5, -15, 25]))
-    r = np.array(
-        [
-            100.0000000000,
-            52.5224541378,
-            46.8714823296,
-            84.5391163313,
-            21.9244774232,
-            46.1081101229,
-            15.4779832381,
-        ]
-    )
-    v_ref = thermo.array.celsius_to_kelvin(np.array([20.0, 10, -10, 32, -15, -24, -3]))
+@pytest.mark.parametrize(
+    "t,r,kwargs,expected_values",
+    [
+        (
+            [20.0, 20, 0, 35, 5, -15, 25, 25],
+            [
+                100.0000000000,
+                52.5224541378,
+                46.8714823296,
+                84.5391163313,
+                21.9244774232,
+                46.1081101229,
+                15.4779832381,
+                0,
+            ],
+            {},
+            [20.0, 10, -10, 32, -15, -24, -3, np.nan],
+        ),
+        (
+            [20.0, 20.0, 20.0],
+            [
+                52.5224541378,
+                0.0,
+                0.000001,
+            ],
+            {"eps": 1e-3, "out": thermo.array.celsius_to_kelvin(100)},
+            [10, 100, 100],
+        ),
+        (
+            [20.0, 20.0, 20.0],
+            [
+                52.5224541378,
+                0.0,
+                0.000001,
+            ],
+            {"eps": 1e-3, "out": np.nan},
+            [10, np.nan, np.nan],
+        ),
+    ],
+)
+def test_dewpoint_from_relative_humidity(t, r, kwargs, expected_values):
     # reference was tested with an online relhum calculator at:
     # https://bmcnoldy.rsmas.miami.edu/Humidity.html
-    td = thermo.array.dewpoint_from_relative_humidity(t, r)
-    np.testing.assert_allclose(td, v_ref)
-
 
-def test_dewpoint_from_specific_humidity():
-    p = np.array([967.5085, 936.3775, 872.248, 756.1647, 649.157, 422.4207]) * 100
-    q = np.array(
-        [
-            0.0169461501,
-            0.0155840075,
-            0.0134912382,
-            0.0083409720,
-            0.0057268584,
-            0.0025150791,
-        ]
-    )
-    v_ref = thermo.array.celsius_to_kelvin(
-        np.array([21.78907, 19.90885, 16.50236, 7.104064, -0.3548709, -16.37916])
-    )
-    td = thermo.array.dewpoint_from_specific_humidity(q, p)
-    np.testing.assert_allclose(td, v_ref)
+    multi = isinstance(t, list)
+    if multi:
+        t = np.array(t)
+        r = np.array(r)
+        expected_values = np.array(expected_values)
+
+    t = thermo.array.celsius_to_kelvin(t)
+    v_ref = thermo.array.celsius_to_kelvin(expected_values)
+
+    td = thermo.array.dewpoint_from_relative_humidity(t, r, **kwargs)
+    if multi:
+        assert np.allclose(td, v_ref, equal_nan=True)
+    else:
+        assert np.isclose(td, v_ref, equal_nan=True)
+
+
+@pytest.mark.parametrize(
+    "q,p,kwargs,expected_values",
+    [
+        (
+            [0.0169461501, 0.0155840075, 0.0134912382, 0.0083409720, 0.0057268584, 0.0025150791, 0],
+            [967.5085, 936.3775, 872.248, 756.1647, 649.157, 422.4207, 422.4207],
+            {},
+            [21.78907, 19.90885, 16.50236, 7.104064, -0.3548709, -16.37916, np.nan],
+        ),
+        (
+            [
+                0.0169461501,
+                0.0,
+                0.000001,
+            ],
+            [967.5085, 967.5085, 967.5085],
+            {"eps": 1e-3, "out": thermo.array.celsius_to_kelvin(100)},
+            [21.78907, 100, 100],
+        ),
+        (
+            [
+                0.0169461501,
+                0.0,
+                0.000001,
+            ],
+            [967.5085, 967.5085, 967.5085],
+            {"eps": 1e-3, "out": np.nan},
+            [21.78907, np.nan, np.nan],
+        ),
+    ],
+)
+def test_dewpoint_from_specific_humidity(q, p, kwargs, expected_values):
+    multi = isinstance(q, list)
+    if multi:
+        q = np.array(q)
+        p = np.array(p)
+        expected_values = np.array(expected_values)
+
+    p = p * 100.0
+    v_ref = thermo.array.celsius_to_kelvin(expected_values)
+
+    td = thermo.array.dewpoint_from_specific_humidity(q, p, **kwargs)
+    if multi:
+        assert np.allclose(td, v_ref, equal_nan=True)
+    else:
+        assert np.isclose(td, v_ref, equal_nan=True)
 
 
 def test_virtual_temperature():