Skip to content

Commit

Permalink
Control zero humidity in thermo computations (#23)
Browse files Browse the repository at this point in the history 
* Control zero humidity in thermo computations
  • Loading branch information
@sandorkertesz
sandorkertesz authored Nov 4, 2024
1 parent d6175a8 commit 3b404c0
Show file tree
Hide file tree
Showing 2 changed files with 220 additions and 50 deletions.
93 changes: 79 additions & 14 deletions src/earthkit/meteo/thermo/array/thermo.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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
-------
Expand All @@ -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):
Expand Down Expand Up @@ -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 ``r``
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
Expand All @@ -782,11 +813,24 @@ 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 = np.asarray(td)
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
Expand All @@ -795,11 +839,18 @@ def dewpoint_from_specific_humidity(q, p):
Specific humidity (kg/kg)
p: ndarray
Pressure (Pa)
eps: number
If ``out`` is not None, return ``out`` when ``q`` < ``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 ``q`` < ``eps``. If None, ``eps`` is
ignored and return np.nan for ``q`` values very close to zero.
Returns
-------
ndarray
Dewpoint (K)
Dewpoint temperature (K)
The computation starts with determining the the saturation vapour pressure over
Expand All @@ -819,7 +870,21 @@ 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 = np.asarray(td)
td[mask] = out

return td


def virtual_temperature(t, q):
Expand All @@ -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)
Expand Down
177 changes: 141 additions & 36 deletions tests/thermo/test_thermo_array.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
import os

import numpy as np
import pytest

from earthkit.meteo import thermo

Expand Down Expand Up @@ -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),
Expand All @@ -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():
Expand Down Expand Up @@ -421,43 +447,122 @@ 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],
),
(20.0, 52.5224541378, {}, 10.0),
(20.0, 0.0, {"eps": 1e-3, "out": thermo.array.celsius_to_kelvin(100)}, 100),
(20.0, 0.000001, {"eps": 1e-3, "out": thermo.array.celsius_to_kelvin(100)}, 100),
(20.0, 0.0, {"eps": 1e-3, "out": np.nan}, np.nan),
(20, 0.000001, {"eps": 1e-3, "out": 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],
),
(0.0169461501, 967.508, {}, 21.78907),
(0.0, 967.5085, {"eps": 1e-3, "out": thermo.array.celsius_to_kelvin(100)}, 100),
(0.000001, 967.5085, {"eps": 1e-3, "out": thermo.array.celsius_to_kelvin(100)}, 100),
(0.0, 967.5085, {"eps": 1e-3, "out": np.nan}, np.nan),
(0.000001, 967.5085, {"eps": 1e-3, "out": 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():
Expand Down

0 comments on commit 3b404c0

Please sign in to comment.