Better relative humidity - allow invalid_values (#2254)

<!--Please ensure the PR fulfills the following requirements! -->
<!-- If this is your first PR, make sure to add your details to the
AUTHORS.rst! -->
### Pull Request Checklist:
- [ ] This PR addresses an already opened issue (for bug fixes /
features)
    - This PR fixes #xyz
- [x] Tests for the changes have been added (for bug fixes / features)
- [x] (If applicable) Documentation has been added / updated (for bug
fixes / features)
- [x] CHANGELOG.rst has been updated (with summary of main changes)
- [x] Link to issue (:issue:`number`) and pull request (:pull:`number`)
has been added

### What kind of change does this PR introduce?
The previous way of computing `hurs` from `huss`, `tas` and `ps` was
taking the ratio of the mixing ratios. It seems that is not how it
should be done... MetPy does something similar, but with an added term.

I decided to instead use the base definition : hurs is the ratio of
vapor pressure to saturation vapor pressure. It corresponds to a WMO
equation to, so that's better. I can't find where I took the previous
equation. It's been 5 years and we were not adding references are
thoroughly back then.

I harmonized the $$\epsilon$$ constant (molecular weight ratio of vapor
to dry air) to 0.62198, the value the WMO uses.

Also, I exposed the `invalid_values` to the `relative_humidity`
indicator. So instead of masking hurs > 100, one can now either leave
them be (None) or clip to 100 ('clip').

### Does this PR introduce a breaking change?
Yes as the output has changed. The change is small, but large enough for
some use cases.

### Other information:

$w$ is the mixing ratio, $w_{sat}$ the saturated mixing ratio, $q$ the
specific humidity, $p$ the pressure, $T$ the temperature. $P_w$ the
vapor pressure, $P_{wsat}$ the saturation vapor pressure. $\epsilon =
0.62198$ the molecular weight ratio of vapor to dry air, a constant.

```math
w = \frac{q}{1-q},    w_{sat} = \epsilon\frac{P_{wsat}}{p - P_{wsat}},  P_w = \frac{pq}{\epsilon + (1 - \epsilon)q}
```
Previous : $$RH = 100\frac{w}{w_{sat}}$$
MetPy (WMO eq. 4.A.16) : $$RH = 100 \frac{w}{w_{sat}}\frac{w_{sat} +
\epsilon}{w + \epsilon}$$
New (WMO eq. 4.A.15) : $$RH = 100 \frac{P_w}{P_{wsat}}$$

Author: aulemahal

CHANGELOG.rst

@@ -12,6 +12,7 @@ New indicators and features
 
 Breaking changes
 ^^^^^^^^^^^^^^^^
+* The relative humidity computations from specific humidity, pressure and temperature (``vapor_pressure`` and ``relative_humidity``) were modified to use the fraction of vapour pressure to saturation vapour pressure instead of an incomplete equation with the mixing ratios. Changes are small, but sometimes not negligible. (:pull:`2254`).
 * `black` and `blackdoc` are no longer required for development. `ruff` is now exclusively used for code and code-block formatting. (:pull:`2249`).
 * Python HDF5 libraries now have lower pins to ensure modern versions are preferably installed (`h5netcdf >=1.5.0` and `h5py >=3.12.1`) (:pull:`2253`).
 

src/xclim/indicators/convert/_conversion.py

@@ -246,7 +246,7 @@ def tg(*args, **kwargs):  # numpydoc ignore=GL08
         "tdps": None,
         "huss": {"kind": InputKind.VARIABLE},
         "ps": {"kind": InputKind.VARIABLE},
-        "invalid_values": "mask",
+        "invalid_values": {"default": "mask"},
     },
 )
 

src/xclim/indices/converters.py

@@ -628,13 +628,13 @@ def vapor_pressure(huss: xr.DataArray, ps: xr.DataArray):
 
     .. math::
 
-        e = \frac{pq}{\epsilon\left(1 + q\left(\frac{1}{\epsilon} - 1\right)\right)}
+        e = \frac{pq}{\epsilon + (1 - \epsilon)q}
 
     Where :math:`p` is the pressure, :math:`q` is the specific humidity and :math:`\epsilon` us the ratio of the dry air
-    gas constant to the water vapor gas constant : :math:`\frac{R_{dry}}{R_{vapor}} = 0.621981`.
+    gas constant to the water vapor gas constant : :math:`\frac{R_{dry}}{R_{vapor}} = 0.62198`.
     """
-    eps = 0.621981
-    e = cast(xr.DataArray, ps * huss / (eps * (1 + huss * (1 / eps - 1))))
+    eps = 0.62198
+    e = cast(xr.DataArray, ps * huss / (eps + (1 - eps) * huss))
     return e.assign_attrs(units=ps.attrs["units"])
 
 
@@ -772,16 +772,19 @@ def relative_humidity(
 
        RH = 100\frac{e_{sat}(T_d)}{e_{sat}(T)}
 
-    Otherwise, the specific humidity and the air pressure must be given so relative humidity can be computed as:
+    Otherwise, the specific humidity and the air pressure must be given so relative humidity
+    can be computed as the ratio of actual vapor pressure to saturation vapor pressure:
 
     .. math::
 
-       RH = 100\frac{w}{w_{sat}}
-       w = \frac{q}{1-q}
-       w_{sat} = 0.622\frac{e_{sat}}{P - e_{sat}}
+       RH = 100\frac{P_w}{P_{wsat}}
+       P_w = \frac{pq}{\epsilon\left(1 + q\left(\frac{1}{\epsilon} - 1\right)\right)}
+       \epsilon = 0.62198
 
-    The methods differ by how :math:`e_{sat}` is computed.
-    See the doc of :py:func:`xclim.core.utils.saturation_vapor_pressure`.
+    The methods differ by how :math:`P_{wsat}` is computed.
+    See the doc of :py:func:`saturation_vapor_pressure` and :py:func:`vapor_pressure`.
+    This equation for RH is the same as eq. 4.A.15 of :cite:p:`world_meteorological_organization_guide_2008`
+    and differs very slightly from MetPy which uses 4.A.16 by computing the mixing ratios first.
 
     References
     ----------
@@ -820,15 +823,13 @@ def relative_humidity(
     elif huss is not None and ps is not None:
         ps = convert_units_to(ps, "Pa")
         huss = convert_units_to(huss, "")
-        tas = convert_units_to(tas, "K")
 
-        e_sat = saturation_vapor_pressure(
+        P_w = vapor_pressure(huss=huss, ps=ps)
+        P_wsat = saturation_vapor_pressure(
             tas=tas, ice_thresh=ice_thresh, method=method, interp_power=interp_power, water_thresh=water_thresh
         )
 
-        w = huss / (1 - huss)
-        w_sat = 0.62198 * e_sat / (ps - e_sat)  # type: ignore
-        hurs = 100 * w / w_sat
+        hurs = 100 * P_w / P_wsat
     else:
         raise ValueError("`huss` and `ps` must be provided if `tdps` is not given.")
 
@@ -933,7 +934,7 @@ def specific_humidity(
         tas=tas, ice_thresh=ice_thresh, method=method, interp_power=interp_power, water_thresh=water_thresh
     )
 
-    w_sat = 0.621981 * e_sat / (ps - e_sat)  # type: ignore
+    w_sat = 0.62198 * e_sat / (ps - e_sat)  # type: ignore
     w = w_sat * hurs
     q: xr.DataArray = w / (1 + w)
 
@@ -1009,7 +1010,7 @@ def specific_humidity_from_dewpoint(
     ...     method="wmo08",
     ... )
     """
-    EPSILON = 0.621981  # molar weight of water vs dry air []
+    EPSILON = 0.62198  # molar weight of water vs dry air []
     e = saturation_vapor_pressure(
         tas=tdps, method=method, ice_thresh=ice_thresh, interp_power=interp_power, water_thresh=water_thresh
     )  # vapour pressure [Pa]

tests/test_converters.py

@@ -151,7 +151,7 @@ def test_saturation_vapor_pressure(tas_series):
 
 def test_relative_humidity(tas_series, hurs_series, huss_series, ps_series):
     tas = tas_series(np.array([-10, -10, 10, 20, 35, 50, 75, 95]) + K2C)
-    hurs_exp = hurs_series([np.nan, 63.0, 66.0, 34.0, 14.0, 6.0, 1.0, 0.0])
+    hurs_exp = hurs_series([np.nan, 62.5, 66.0, 35.0, 14.5, 6.5, 2.0, 1.0])
     ps = ps_series([101325] * 8)
     huss = huss_series([0.003, 0.001] + [0.005] * 7)
 

tests/test_indices.py

@@ -3250,8 +3250,8 @@ def test_vapor_pressure_deficit(tas_series, hurs_series, method):
 def test_relative_humidity(tas_series, hurs_series, huss_series, ps_series, method, invalid_values, exp0):
     tas = tas_series(np.array([-10, -10, 10, 20, 35, 50, 75, 95]) + K2C)
 
-    # Expected values obtained with the Sonntag90 method
-    hurs_exp = hurs_series([exp0, 63.0, 66.0, 34.0, 14.0, 6.0, 1.0, 0.0])
+    # Expected values obtained with the Sonntag90 method, rounded to the nearest half unit
+    hurs_exp = hurs_series([exp0, 62.5, 66.0, 35.0, 14.5, 6.5, 2.0, 1.0])
     ps = ps_series([101325] * 8)
     huss = huss_series([0.003, 0.001] + [0.005] * 7)