Repository: Ouranosinc/xclim

Issue:
Humidex indicator
### Description
<!--Describe what you were trying to get done.
Tell us what happened, what went wrong, and what you expected to happen.-->

Add humidex indicator. 

https://www.canada.ca/en/environment-climate-change/services/seasonal-weather-hazards/warm-season-weather-hazards.html#toc7

https://en.wikipedia.org/wiki/Humidex


Discussion:

Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index f812dc8a..bfb2f18d 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -34,6 +34,7 @@ New indicators
 ~~~~~~~~~~~~~~
 * `atmos.corn_heat_units` computes the daily temperature-based index for corn growth.
 * New indices and indicators for `tx_days_below`, `tg_days_above`, `tg_days_below`, and `tn_days_above`.
+* `atmos.humidex` returns the Canadian *humidex*, an indicator of perceived temperature account for relative humidity.
 * `multiday_temperature_swing` indice for returning general statistics based on spells of doubly-thresholded temperatures (Tmin < T1, Tmax > T2).
 * New indicators `atmos.freezethaw_frequency`, `atmos.freezethaw_spell_mean_length`, `atmos.freezethaw_spell_max_length` for statistics of Tmin < 0 degC and Tmax > 0 deg C days now available (wrapped from `multiday_temperature_swing`).
 * `atmos.wind_chill_index` computes the daily wind chill index. The default is similar to what Environment and Climate Change Canada does, options are tunable to get the version of the National Weather Service.
diff --git a/xclim/data/fr.json b/xclim/data/fr.json
index e0b161cd..064b4d7e 100644
--- a/xclim/data/fr.json
+++ b/xclim/data/fr.json
@@ -860,5 +860,11 @@
     "abstract": "Le facteur éolien est un indice qui équivaut à la sensation du froid par un individu moyen. Il est calculé à partir de la température et de la vitesse du vent à 10 m. Tel que définit par Environnement et Changement Climatique Canada, une seconde formule est utilisé pour les vents faibles. La formule standard est sinon la même qu'utilisée aux États-Unis.",
     "description": "Refroiddissement éolien, équivalent à la sensation du froid due au vent par in individu moyen.",
     "long_name": "Facteur éolien."
+  },
+  "HUMIDEX": {
+    "title": "Indice humidex",
+    "abstract": "L'indice humidex décrit la température ressentie par une personne lorsqu'on tient compte de l'humidité relative. Il peut être interpreté comme la température équivalente ressentie lorsque l'air est sec.",
+    "description": "Indice humidex décrivant la température ressentie par une personne en réponse à l'humidité relative.",
+    "long_name": "indice humidex"
   }
 }
diff --git a/xclim/indicators/atmos/_conversion.py b/xclim/indicators/atmos/_conversion.py
index 59b58535..bd090e21 100644
--- a/xclim/indicators/atmos/_conversion.py
+++ b/xclim/indicators/atmos/_conversion.py
@@ -6,6 +6,7 @@ from xclim.core.indicator import Indicator
 from xclim.core.utils import wrapped_partial
 
 __all__ = [
+    "humidex",
     "tg",
     "wind_speed_from_vector",
     "wind_vector_from_speed",
@@ -25,6 +26,18 @@ class Converter(Indicator):
     missing = "skip"
 
 
+humidex = Converter(
+    identifier="humidex",
+    nvar=2,
+    units="C",
+    standard_name="air_temperature",
+    long_name="humidex index",
+    description="Humidex index describing the temperature felt by the average person in response to relative humidity.",
+    cell_methods="",
+    compute=indices.humidex,
+)
+
+
 tg = Converter(
     identifier="tg",
     nvar=2,
diff --git a/xclim/indices/_conversion.py b/xclim/indices/_conversion.py
index 9c6cba25..5b9789cf 100644
--- a/xclim/indices/_conversion.py
+++ b/xclim/indices/_conversion.py
@@ -4,9 +4,10 @@ from typing import Tuple
 import numpy as np
 import xarray as xr
 
-from xclim.core.units import convert_units_to, declare_units
+from xclim.core.units import convert_units_to, declare_units, units2pint
 
 __all__ = [
+    "humidex",
     "tas",
     "uas_vas_2_sfcwind",
     "sfcwind_2_uas_vas",
@@ -19,6 +20,75 @@ __all__ = [
 ]
 
 
+@declare_units(tas="[temperature]", tdps="[temperature]")
+def humidex(
+    tas: xr.DataArray,
+    tdps: xr.DataArray,
+) -> xr.DataArray:
+    r"""Humidex index.
+
+    The humidex indicates how hot the air feels to an average person, accounting for the effect of humidity. It
+    can be loosely interpreted as the equivalent perceived temperature when the air is dry.
+
+    Parameters
+    ----------
+    tas : xarray.DataArray
+      Air temperature.
+    tdps : xarray.DataArray,
+      Dewpoint temperature.
+
+    Returns
+    -------
+    xarray.DataArray, [temperature]
+      The humidex index.
+
+    Notes
+    -----
+    The humidex is usually computed using hourly observations of dry bulb and dewpoint temperatures. It is computed
+    using the formula based on [masterton79]_:
+
+    .. math::
+
+       T_{\text{dry bulb}}+{\frac {5}{9}}\left[6.11\times \exp(5417.7530\left({\frac {1}{273.16}}-{\frac {1}{T_{\text{
+       dewpoint}}}}\right)-10\right]
+
+    where :math:`T_{dry bulb}` is the dry bulb air temperature (°C), and :math:`T_{dewpoint}` the dewpoint
+    temperature (°K). The constant 5417.753 reflects the molecular weight of water, latent heat of vaporization,
+    and the universal gas constant ([mekis15]_).
+
+    The humidex *comfort scale* ([eccc]_) can be interpreted as follows:
+
+    - 20 to 29 : no discomfort;
+    - 30 to 39 : some discomfort;
+    - 40 to 45 : great discomfort, avoid exertion;
+    - 46 and over : dangerous, possible heat stroke;
+
+    References
+    ----------
+    .. [masterton79] Masterton, J. M., & Richardson, F. A. (1979). HUMIDEX, A method of quantifying human discomfort due to excessive heat and humidity, CLI 1-79. Downsview, Ontario: Environment Canada, Atmospheric Environment Service.
+    .. [mekis15] Éva Mekis, Lucie A. Vincent, Mark W. Shephard & Xuebin Zhang (2015) Observed Trends in Severe Weather Conditions Based on Humidex, Wind Chill, and Heavy Rainfall Events in Canada for 1953–2012, Atmosphere-Ocean, 53:4, 383-397, DOI: 10.1080/07055900.2015.1086970
+    .. [eccc] https://climate.weather.gc.ca/glossary_e.html
+    """
+    # Convert dewpoint temperature to Kelvins
+    tdps = convert_units_to(tdps, "kelvin")
+
+    # Vapour pressure in hPa
+    e = 6.11 * np.exp(5417.7530 * (1 / 273.16 - 1.0 / tdps))
+
+    # Temperature delta due to humidity in delta_degC
+    h = 5 / 9 * (e - 10)
+    h.attrs["units"] = "delta_degree_Celsius"
+
+    # Get delta_units for output
+    du = (1 * units2pint(tas) - 0 * units2pint(tas)).units
+    h = convert_units_to(h, du)
+
+    # Add the delta to the input temperature
+    out = h + tas
+    out.attrs["units"] = tas.units
+    return out
+
+
 @declare_units(tasmin="[temperature]", tasmax="[temperature]")
 def tas(tasmin: xr.DataArray, tasmax: xr.DataArray) -> xr.DataArray:
     """Average temperature from minimum and maximum temperatures.


Test patch:
diff --git a/xclim/testing/tests/test_atmos.py b/xclim/testing/tests/test_atmos.py
index eb0148df..272383f3 100644
--- a/xclim/testing/tests/test_atmos.py
+++ b/xclim/testing/tests/test_atmos.py
@@ -69,6 +69,20 @@ def test_relative_humidity_dewpoint(tas_series, hurs_series):
     )
 
 
+def test_humidex(tas_series):
+
+    tas = tas_series([15, 25, 35, 40])
+    tas.attrs["units"] = "C"
+
+    dtas = tas_series([10, 15, 25, 25])
+    dtas.attrs["units"] = "C"
+
+    # expected values from https://en.wikipedia.org/wiki/Humidex
+    h = atmos.humidex(tas, dtas)
+    np.testing.assert_array_almost_equal(h, [16, 29, 47, 52], 0)
+    assert h.name == "humidex"
+
+
 def test_saturation_vapor_pressure(tas_series):
     tas = tas_series(np.array([-20, -10, -1, 10, 20, 25, 30, 40, 60]) + K2C)
     e_sat_exp = [103, 260, 563, 1228, 2339, 3169, 4247, 7385, 19947]
diff --git a/xclim/testing/tests/test_indices.py b/xclim/testing/tests/test_indices.py
index ad6ff594..e953d23c 100644
--- a/xclim/testing/tests/test_indices.py
+++ b/xclim/testing/tests/test_indices.py
@@ -24,7 +24,7 @@ import xarray as xr
 from xclim import indices as xci
 from xclim.core.calendar import percentile_doy
 from xclim.core.options import set_options
-from xclim.core.units import ValidationError, convert_units_to
+from xclim.core.units import ValidationError, convert_units_to, units
 from xclim.testing import open_dataset
 
 K2C = 273.15
@@ -1827,6 +1827,30 @@ def test_winter_storm(snd_series):
     np.testing.assert_array_equal(out, [3])
 
 
+def test_humidex(tas_series):
+
+    tas = tas_series([15, 25, 35, 40])
+    tas.attrs["units"] = "C"
+
+    dtas = tas_series([10, 15, 25, 25])
+    dtas.attrs["units"] = "C"
+
+    # expected values from https://en.wikipedia.org/wiki/Humidex
+    expected = np.array([16, 29, 47, 52]) * units.degC
+
+    # Celcius
+    hc = xci.humidex(tas, dtas)
+    np.testing.assert_array_almost_equal(hc, expected, 0)
+
+    # Kelvin
+    hk = xci.humidex(convert_units_to(tas, "K"), dtas)
+    np.testing.assert_array_almost_equal(hk, expected.to("K"), 0)
+
+    # Fahrenheit
+    hf = xci.humidex(convert_units_to(tas, "fahrenheit"), dtas)
+    np.testing.assert_array_almost_equal(hf, expected.to("fahrenheit"), 0)
+
+
 @pytest.mark.parametrize(
     "op,exp", [("max", 11), ("sum", 21), ("count", 3), ("mean", 7)]
 )