API reference

The thermal-index functions are available at the top level of the package, for example thermofeel.calculate_utci(...) or thermofeel.calculate_excess_heat_factor(...). The supporting functions for the excess heat/cold factors (daily mean temperature, significance and acclimatisation indices, heatwave severity) live in the thermofeel.excess_heat submodule. This page is generated from the source docstrings.

All functions are vectorised over NumPy and take SI inputs; pass arrays and wrap scalars in an array (e.g. np.array([300.0])). See Calling convention for details.

Thermal indices and supporting quantities

thermofeel is a library to calculate human thermal comfort indexes.

Currently calculates the thermal indexes: * Universal Thermal Climate Index * Apparent Temperature * Apparent Temperature (radiation) * Heat Index Adjusted * Heat Index Simplified * Humidex * Discomfort Index * Summer Simmer Index * Normal Effective Temperature * Relative Strain Index * Wet Bulb Globe Temperature * Wet Bulb Globe Temperature Simple * Wet Bulb Globe Temperature (Liljegren method) * Heat Force (KNMI 0-10 heat-stress scale) * Excess Heat Factor and Excess Cold Factor * Predicted Mean Vote and Predicted Percentage of Dissatisfied (PMV/PPD) * Wind Chill

In support of the above indexes, it also calculates: * Globe Temperature * Mean Radiant Temperature * Mean Radiant Temperature from Globe Temperature * Relative Humidity Percentage * Saturation vapour pressure * Wet Bulb Temperature

The thermofeel.approximations namespace additionally provides clearly-labelled estimators (currently the direct solar radiation fdir from global radiation) for datasets that lack an input that a radiation index needs.

To calculate the cos of the solar zenith angle, we suggest to use the earthkit-meteo library (github.com:ecmwf/earthkit-meteo.git)

calculate_relative_humidity_percent

calculate_relative_humidity_percent(
    t2_k: ArrayLike, td_k: ArrayLike
) -> np.ndarray

Relative Humidity in percent :param t2_k: (float array) 2m temperature [K] :param td_k: (float array) dew point temperature [K] returns relative humidity [%]

Uses the Magnus-Tetens saturation vapour pressure (coefficients over water) - a different empirical form to calculate_saturation_vapour_pressure (Hardy 1998). The result is not clamped: when td_k > t2_k (supersaturation) it exceeds 100%.

Reference: Tetens (1930); coefficients per Murray (1967) https://doi.org/10.1175/1520-0450(1967)006<0203:OTCOSV>2.0.CO;2

Source code in thermofeel/thermofeel.py

def calculate_relative_humidity_percent(t2_k: ArrayLike, td_k: ArrayLike) -> np.ndarray:
    """
    Relative Humidity in percent
        :param t2_k: (float array) 2m temperature [K]
        :param td_k: (float array) dew point temperature [K]
        returns relative humidity [%]

    Uses the Magnus-Tetens saturation vapour pressure (coefficients over water) -
    a different empirical form to ``calculate_saturation_vapour_pressure``
    (Hardy 1998). The result is not clamped: when ``td_k > t2_k`` (supersaturation)
    it exceeds 100%.

    Reference: Tetens (1930); coefficients per Murray (1967)
    https://doi.org/10.1175/1520-0450(1967)006<0203:OTCOSV>2.0.CO;2
    """

    t2_c = kelvin_to_celsius(t2_k)
    td_c = kelvin_to_celsius(td_k)
    # saturated vapour pressure
    es = 6.11 * 10.0 ** (7.5 * t2_c / (237.3 + t2_c))
    # vapour pressure
    e = 6.11 * 10.0 ** (7.5 * td_c / (237.3 + td_c))
    rh = (e / es) * 100
    return rh

calculate_saturation_vapour_pressure

calculate_saturation_vapour_pressure(
    t2_k: ArrayLike,
) -> np.ndarray

Saturation vapour pressure over water :param t2_k: (float array) 2m temperature [K] returns saturation vapor pressure over water in the pure phase [hPa] == [mBar] Reference: Hardy (1998) https://www.decatur.de/javascript/dew/resources/its90formulas.pdf

Source code in thermofeel/thermofeel.py

def calculate_saturation_vapour_pressure(t2_k: ArrayLike) -> np.ndarray:
    """
    Saturation vapour pressure over water
        :param t2_k: (float array) 2m temperature [K]
        returns saturation vapor pressure over water in the pure phase [hPa] == [mBar]
    Reference: Hardy (1998)
    https://www.decatur.de/javascript/dew/resources/its90formulas.pdf
    """

    g = [
        -2.8365744e3,
        -6.028076559e3,
        1.954263612e1,
        -2.737830188e-2,
        1.6261698e-5,
        7.0229056e-10,
        -1.8680009e-13,
        2.7150305,
    ]
    ess = g[7] * np.log(t2_k)
    for i in range(7):
        ess = ess + g[i] * np.power(t2_k, (i - 2))

    ess = np.exp(ess) * 0.01  # hPa

    return ess

calculate_saturation_vapour_pressure_multiphase

calculate_saturation_vapour_pressure_multiphase(
    t2_k: ArrayLike, phase: ArrayLike
) -> np.ndarray

Saturation vapour pressure over liquid water and ice :param t2_k: (float array) 2m temperature [K] :param phase: (int array) 0 over liquid water, 1 over ice (same shape as t2_k) returns pressure of water vapor over a surface of liquid water or ice [hPa] == [mBar]

t2_k and phase are array-like (wrap a scalar in an array). Only elements with phase equal to 0 (liquid) or 1 (ice) are computed; any other phase value leaves that element at 0 hPa.

Reference: ECMWF IFS Documentation CY45R1 - Part IV : Physical processes (2018) pp. 116 https://doi.org/10.21957/4whwo8jw0 https://metview.readthedocs.io/en/latest/api/functions/saturation_vapour_pressure.html

Source code in thermofeel/thermofeel.py

def calculate_saturation_vapour_pressure_multiphase(
    t2_k: ArrayLike, phase: ArrayLike
) -> np.ndarray:
    """
    Saturation vapour pressure over liquid water and ice
        :param t2_k: (float array) 2m temperature [K]
        :param phase: (int array) 0 over liquid water, 1 over ice (same shape as
            t2_k)
        returns pressure of water vapor over a surface of liquid water or ice [hPa] == [mBar]

    ``t2_k`` and ``phase`` are array-like (wrap a scalar in an array). Only
    elements with ``phase`` equal to 0 (liquid) or 1 (ice) are computed; any
    other ``phase`` value leaves that element at 0 hPa.

    Reference: ECMWF IFS Documentation CY45R1 - Part IV : Physical processes (2018) pp. 116
    https://doi.org/10.21957/4whwo8jw0
    https://metview.readthedocs.io/en/latest/api/functions/saturation_vapour_pressure.html
    """
    T0 = 273.16  # triple point of water 273.16 K (0.01 °C) at 611.73 Pa
    # float dtype so integer-typed temperature input is not truncated on assignment
    es = np.zeros_like(t2_k, dtype=float)
    y = (t2_k - T0) / (t2_k - 32.19)  # over liquid water
    es[phase == 0] = 6.1121 * np.exp(17.502 * y[phase == 0])
    y = (t2_k - T0) / (t2_k + 0.7)  # over ice
    es[phase == 1] = 6.1121 * np.exp(22.587 * y[phase == 1])

    return es

calculate_nonsaturation_vapour_pressure

calculate_nonsaturation_vapour_pressure(
    t2_k: ArrayLike, rh: ArrayLike
) -> np.ndarray

Non saturated vapour pressure :param t2_k: (float array) 2m temperature [K] :param rh: (float array) relative humidity percentage [%] returns non saturated vapor pressure [hPa] == [mBar] Reference: Bureau of Meteorology (2010) http://www.bom.gov.au/info/thermal_stress/#approximation

Source code in thermofeel/thermofeel.py

def calculate_nonsaturation_vapour_pressure(
    t2_k: ArrayLike, rh: ArrayLike
) -> np.ndarray:
    """
    Non saturated vapour pressure
        :param t2_k: (float array) 2m temperature [K]
        :param rh: (float array) relative humidity percentage [%]
        returns non saturated vapor pressure [hPa] == [mBar]
    Reference: Bureau of Meteorology (2010)
    http://www.bom.gov.au/info/thermal_stress/#approximation
    """

    t2_c = kelvin_to_celsius(t2_k)
    ens = rh / 100 * 6.105 * np.exp(17.27 * t2_c / (237.7 + t2_c))

    return ens

scale_windspeed

scale_windspeed(va: ArrayLike, h: ArrayLike) -> np.ndarray

Scaling wind speed from 10 metres to height h :param va: (float array) 10m wind speed [m/s] :param h: (float array) height at which wind speed needs to be scaled [m] returns wind speed at height h Reference: Bröde et al. (2012) https://doi.org/10.1007/s00484-011-0454-1

Source code in thermofeel/thermofeel.py

def scale_windspeed(va: ArrayLike, h: ArrayLike) -> np.ndarray:
    """
    Scaling wind speed from 10 metres to height h
        :param va: (float array) 10m wind speed [m/s]
        :param h: (float array) height at which wind speed needs to be scaled [m]
        returns wind speed at height h
    Reference: Bröde et al. (2012)
    https://doi.org/10.1007/s00484-011-0454-1
    """
    # log-law scaling factor 1 / log10(z_ref / z0), with the 10 m reference
    # height and a 0.01 m roughness length (evaluates to 1/3).
    c = 1 / np.log10(10 / 0.01)
    vh = va * np.log10(h / 0.01) * c

    return vh

approximate_dsrp

approximate_dsrp(
    fdir: ArrayLike,
    cossza: ArrayLike,
    threshold: float = 0.1,
) -> np.ndarray

Helper function to approximate dsrp from fdir and cossza

By geometry the direct solar radiation perpendicular to the beam is dsrp = fdir / cos(sza) = fdir / cossza; this is applied where cossza > threshold and left as fdir below it. The approximation introduces large errors as cossza approaches zero, so use it only if dsrp is not available in your dataset. fdir and cossza are array-like (wrap a scalar in an array). To compute cossza consider using earthkit-meteo.solar.cos_solar_zenith_angle. :param fdir: (float array) total sky direct solar radiation at surface [W m-2] :param cossza: (float array) cosine of solar zenith angle [dimentionless] :param threshold: (float) minimum cossza for which fdir is divided (default 0.1) returns direct radiation from the Sun [W m-2]

Source code in thermofeel/thermofeel.py

def approximate_dsrp(
    fdir: ArrayLike, cossza: ArrayLike, threshold: float = 0.1
) -> np.ndarray:
    """
    Helper function to approximate dsrp from fdir and cossza

    By geometry the direct solar radiation perpendicular to the beam is
    ``dsrp = fdir / cos(sza) = fdir / cossza``; this is applied where
    ``cossza > threshold`` and left as ``fdir`` below it. The approximation
    introduces large errors as cossza approaches zero, so use it only if dsrp is
    not available in your dataset. ``fdir`` and ``cossza`` are array-like (wrap a
    scalar in an array). To compute cossza consider using
    earthkit-meteo.solar.cos_solar_zenith_angle.
        :param fdir: (float array) total sky direct solar radiation at surface [W m-2]
        :param cossza: (float array) cosine of solar zenith angle [dimentionless]
        :param threshold: (float) minimum cossza for which fdir is divided (default 0.1)
        returns direct radiation from the Sun [W m-2]
    """
    # filter statement for solar zenith angle to avoid division by zero.
    csza_filter1 = np.where((cossza > threshold))
    # float copy so integer-typed fdir input is not truncated on assignment below
    dsrp = np.array(fdir, dtype=float)  # leave dsrp = fdir where cossza <= threshold
    dsrp[csza_filter1] = dsrp[csza_filter1] / cossza[csza_filter1]
    return dsrp

calculate_dew_point_from_relative_humidity

calculate_dew_point_from_relative_humidity(
    rh: ArrayLike, t2_k: ArrayLike
) -> np.ndarray

Dew point temperature at 2m from relative humidity in percent :param rh: (float array) relative humidity [%] :param t2_k: (float array) 2m temperature [K] returns dew point temperature [K] Reference: Alduchov and Eskridge (1996) https://doi.org/10.1175/1520-0450(1996)035<0601:IMFAOS>2.0.CO;2

Source code in thermofeel/thermofeel.py

def calculate_dew_point_from_relative_humidity(
    rh: ArrayLike, t2_k: ArrayLike
) -> np.ndarray:
    """
    Dew point temperature at 2m from relative humidity in percent
        :param rh: (float array) relative humidity [%]
        :param t2_k: (float array) 2m temperature [K]
        returns dew point temperature [K]
    Reference: Alduchov and Eskridge (1996)
    https://doi.org/10.1175/1520-0450(1996)035<0601:IMFAOS>2.0.CO;2
    """
    t2_c = kelvin_to_celsius(t2_k)
    td_c = (
        243.04
        * (np.log(rh / 100) + ((17.625 * t2_c) / (243.04 + t2_c)))
        / (17.625 - np.log(rh / 100) - ((17.625 * t2_c) / (243.04 + t2_c)))
    )
    td_k = celsius_to_kelvin(td_c)
    return td_k

calculate_mean_radiant_temperature

calculate_mean_radiant_temperature(
    ssrd: ArrayLike,
    ssr: ArrayLike,
    dsrp: ArrayLike,
    strd: ArrayLike,
    fdir: ArrayLike,
    strr: ArrayLike,
    cossza: ArrayLike,
) -> np.ndarray

MRT - Mean Radiant Temperature To compute cossza consider using earhkit-meteo.solar.calculate_cos_solar_zenith_angle :param ssrd: (float array) surface solar radiation downwards [W m-2] :param ssr: (float array) surface net solar radiation [W m-2] :param dsrp: (float array) direct solar radiation [W m-2] :param strd: (float array) surface thermal radiation downwards [W m-2] :param fdir: (float array) total sky direct solar radiation at surface [W m-2] :param strr: (float array) surface net thermal radiation [W m-2] :param cossza: (float array) cosine of solar zenith angle [dimentionless] returns mean radiant temperature [K] Reference: Di Napoli et al. (2020) https://link.springer.com/article/10.1007/s00484-020-01900-5

Source code in thermofeel/thermofeel.py

def calculate_mean_radiant_temperature(
    ssrd: ArrayLike,
    ssr: ArrayLike,
    dsrp: ArrayLike,
    strd: ArrayLike,
    fdir: ArrayLike,
    strr: ArrayLike,
    cossza: ArrayLike,
) -> np.ndarray:
    """
    MRT - Mean Radiant Temperature
    To compute cossza consider using earhkit-meteo.solar.calculate_cos_solar_zenith_angle
        :param ssrd: (float array) surface solar radiation downwards [W m-2]
        :param ssr: (float array) surface net solar radiation [W m-2]
        :param dsrp: (float array) direct solar radiation [W m-2]
        :param strd: (float array) surface thermal radiation downwards [W m-2]
        :param fdir: (float array) total sky direct solar radiation at surface [W m-2]
        :param strr: (float array) surface net thermal radiation [W m-2]
        :param cossza: (float array) cosine of solar zenith angle [dimentionless]
        returns mean radiant temperature [K]
    Reference: Di Napoli et al. (2020)
    https://link.springer.com/article/10.1007/s00484-020-01900-5
    """

    dsw = ssrd - fdir
    rsw = ssrd - ssr
    lur = strd - strr
    # Istar = dsrp

    # calculate fp projected factor area
    gamma = np.arcsin(cossza) * 180 / np.pi
    fp = 0.308 * np.cos(to_radians * gamma * (0.998 - gamma * gamma / 50000))

    # calculate mean radiant temperature
    mrt = np.power(
        (
            (1 / 0.0000000567)
            * (
                0.5 * strd
                + 0.5 * lur
                + (0.7 / 0.97) * (0.5 * dsw + 0.5 * rsw + fp * dsrp)
            )
        ),
        0.25,
    )

    return mrt

calculate_utci_polynomial

calculate_utci_polynomial(
    t2m: ArrayLike,
    mrt: ArrayLike,
    va: ArrayLike,
    wvp: ArrayLike,
) -> np.ndarray

Helper function to calculate the UTCI polynomial approximation :param t2m: (float array) is 2m temperature [C] :param mrt: (float array) is mean radiant temperature [C] :param va: (float array) is wind speed at 10 meters [m/s] :param wvp: (float array) is water vapour pressure [kPa] returns UTCI [K] Reference: Brode et al. (2012) https://doi.org/10.1007/s00484-011-0454-1

Source code in thermofeel/thermofeel.py

def calculate_utci_polynomial(
    t2m: ArrayLike, mrt: ArrayLike, va: ArrayLike, wvp: ArrayLike
) -> np.ndarray:
    """
    Helper function to calculate the UTCI polynomial approximation
        :param t2m: (float array) is 2m temperature [C]
        :param mrt: (float array) is mean radiant temperature [C]
        :param va: (float array) is wind speed at 10 meters [m/s]
        :param wvp: (float array) is water vapour pressure [kPa]
    returns UTCI [K]
    Reference: Brode et al. (2012)
    https://doi.org/10.1007/s00484-011-0454-1
    """
    e_mrt = np.subtract(mrt, t2m)

    t2m2 = t2m * t2m
    t2m3 = t2m2 * t2m
    t2m4 = t2m3 * t2m
    t2m5 = t2m4 * t2m
    t2m6 = t2m5 * t2m

    va2 = va * va
    va3 = va2 * va
    va4 = va3 * va
    va5 = va4 * va
    va6 = va5 * va

    e_mrt2 = e_mrt * e_mrt
    e_mrt3 = e_mrt2 * e_mrt
    e_mrt4 = e_mrt3 * e_mrt
    e_mrt5 = e_mrt4 * e_mrt
    e_mrt6 = e_mrt5 * e_mrt

    wvp2 = wvp * wvp
    wvp3 = wvp2 * wvp
    wvp4 = wvp3 * wvp
    wvp5 = wvp4 * wvp
    wvp6 = wvp5 * wvp

    varh2 = va * wvp2
    va2_rh = va2 * wvp
    va2_e_mrt = va2 * e_mrt
    e_mrt_rh = e_mrt * wvp
    e_mrt_rh2 = e_mrt * wvp2
    e_mrt2_rh = e_mrt2 * wvp
    e_mrt2_rh2 = e_mrt2 * wvp2
    e_mrt_rh3 = e_mrt * wvp3
    va_e_mrt = va * e_mrt
    va_e_mrt2 = va * e_mrt2
    va_rh = va * wvp
    t2m_va = t2m * va
    e_mrt3_rh = e_mrt3 * wvp
    e_mrt4_rh = e_mrt4 * wvp

    utci = (
        t2m
        + 6.07562052e-01
        + -2.27712343e-02 * t2m
        + 8.06470249e-04 * t2m2
        + -1.54271372e-04 * t2m3
        + -3.24651735e-06 * t2m4
        + 7.32602852e-08 * t2m5
        + 1.35959073e-09 * t2m6
        + -2.25836520e00 * va
        + 8.80326035e-02 * t2m * va
        + 2.16844454e-03 * t2m2 * va
        + -1.53347087e-05 * t2m3 * va
        + -5.72983704e-07 * t2m4 * va
        + -2.55090145e-09 * t2m5 * va
        + -7.51269505e-01 * va2
        + -4.08350271e-03 * t2m * va2
        + -5.21670675e-05 * t2m2 * va2
        + 1.94544667e-06 * t2m3 * va2
        + 1.14099531e-08 * t2m4 * va2
        + 1.58137256e-01 * va3
        + -6.57263143e-05 * t2m * va3
        + 2.22697524e-07 * t2m2 * va3
        + -4.16117031e-08 * t2m3 * va3
        + -1.27762753e-02 * va4
        + 9.66891875e-06 * t2m * va4
        + 2.52785852e-09 * t2m2 * va4
        + 4.56306672e-04 * va5
        + -1.74202546e-07 * t2m * va5
        + -5.91491269e-06 * va6
        + 3.98374029e-01 * e_mrt
        + 1.83945314e-04 * t2m * e_mrt
        + -1.73754510e-04 * t2m2 * e_mrt
        + -7.60781159e-07 * t2m3 * e_mrt
        + 3.77830287e-08 * t2m4 * e_mrt
        + 5.43079673e-10 * t2m5 * e_mrt
        + -2.00518269e-02 * va_e_mrt
        + 8.92859837e-04 * t2m * va_e_mrt
        + 3.45433048e-06 * t2m2 * va_e_mrt
        + -3.77925774e-07 * t2m3 * va_e_mrt
        + -1.69699377e-09 * t2m4 * va_e_mrt
        + 1.69992415e-04 * va2_e_mrt
        + -4.99204314e-05 * t2m * va2_e_mrt
        + 2.47417178e-07 * t2m2 * va2_e_mrt
        + 1.07596466e-08 * t2m3 * va2_e_mrt
        + 8.49242932e-05 * va3 * e_mrt
        + 1.35191328e-06 * t2m * va3 * e_mrt
        + -6.21531254e-09 * t2m2 * va3 * e_mrt
        + -4.99410301e-06 * va4 * e_mrt
        + -1.89489258e-08 * t2m * va4 * e_mrt
        + 8.15300114e-08 * va5 * e_mrt
        + 7.55043090e-04 * e_mrt2
        + -5.65095215e-05 * t2m * e_mrt2
        + -4.52166564e-07 * t2m2 * e_mrt2
        + 2.46688878e-08 * t2m3 * e_mrt2
        + 2.42674348e-10 * t2m4 * e_mrt2
        + 1.54547250e-04 * va_e_mrt2
        + 5.24110970e-06 * t2m * va_e_mrt2
        + -8.75874982e-08 * t2m2 * va_e_mrt2
        + -1.50743064e-09 * t2m3 * va_e_mrt2
        + -1.56236307e-05 * va2 * e_mrt2
        + -1.33895614e-07 * t2m * va2 * e_mrt2
        + 2.49709824e-09 * t2m2 * va2 * e_mrt2
        + 6.51711721e-07 * va3 * e_mrt2
        + 1.94960053e-09 * t2m * va3 * e_mrt2
        + -1.00361113e-08 * va4 * e_mrt2
        + -1.21206673e-05 * e_mrt3
        + -2.18203660e-07 * t2m * e_mrt3
        + 7.51269482e-09 * t2m2 * e_mrt3
        + 9.79063848e-11 * t2m3 * e_mrt3
        + 1.25006734e-06 * va * e_mrt3
        + -1.81584736e-09 * t2m_va * e_mrt3
        + -3.52197671e-10 * t2m2 * va * e_mrt3
        + -3.36514630e-08 * va2 * e_mrt3
        + 1.35908359e-10 * t2m * va2 * e_mrt3
        + 4.17032620e-10 * va3 * e_mrt3
        + -1.30369025e-09 * e_mrt4
        + 4.13908461e-10 * t2m * e_mrt4
        + 9.22652254e-12 * t2m2 * e_mrt4
        + -5.08220384e-09 * va * e_mrt4
        + -2.24730961e-11 * t2m_va * e_mrt4
        + 1.17139133e-10 * va2 * e_mrt4
        + 6.62154879e-10 * e_mrt5
        + 4.03863260e-13 * t2m * e_mrt5
        + 1.95087203e-12 * va * e_mrt5
        + -4.73602469e-12 * e_mrt6
        + 5.12733497e00 * wvp
        + -3.12788561e-01 * t2m * wvp
        + -1.96701861e-02 * t2m2 * wvp
        + 9.99690870e-04 * t2m3 * wvp
        + 9.51738512e-06 * t2m4 * wvp
        + -4.66426341e-07 * t2m5 * wvp
        + 5.48050612e-01 * va_rh
        + -3.30552823e-03 * t2m * va_rh
        + -1.64119440e-03 * t2m2 * va_rh
        + -5.16670694e-06 * t2m3 * va_rh
        + 9.52692432e-07 * t2m4 * va_rh
        + -4.29223622e-02 * va2_rh
        + 5.00845667e-03 * t2m * va2_rh
        + 1.00601257e-06 * t2m2 * va2_rh
        + -1.81748644e-06 * t2m3 * va2_rh
        + -1.25813502e-03 * va3 * wvp
        + -1.79330391e-04 * t2m * va3 * wvp
        + 2.34994441e-06 * t2m2 * va3 * wvp
        + 1.29735808e-04 * va4 * wvp
        + 1.29064870e-06 * t2m * va4 * wvp
        + -2.28558686e-06 * va5 * wvp
        + -3.69476348e-02 * e_mrt_rh
        + 1.62325322e-03 * t2m * e_mrt_rh
        + -3.14279680e-05 * t2m2 * e_mrt_rh
        + 2.59835559e-06 * t2m3 * e_mrt_rh
        + -4.77136523e-08 * t2m4 * e_mrt_rh
        + 8.64203390e-03 * va * e_mrt_rh
        + -6.87405181e-04 * t2m_va * e_mrt_rh
        + -9.13863872e-06 * t2m2 * va * e_mrt_rh
        + 5.15916806e-07 * t2m3 * va * e_mrt_rh
        + -3.59217476e-05 * va2 * e_mrt_rh
        + 3.28696511e-05 * t2m * va2 * e_mrt_rh
        + -7.10542454e-07 * t2m2 * va2 * e_mrt_rh
        + -1.24382300e-05 * va3 * e_mrt_rh
        + -7.38584400e-09 * t2m * va3 * e_mrt_rh
        + 2.20609296e-07 * va4 * e_mrt_rh
        + -7.32469180e-04 * e_mrt2_rh
        + -1.87381964e-05 * t2m * e_mrt2_rh
        + 4.80925239e-06 * t2m2 * e_mrt2_rh
        + -8.75492040e-08 * t2m3 * e_mrt2_rh
        + 2.77862930e-05 * va * e_mrt2_rh
        + -5.06004592e-06 * t2m_va * e_mrt2_rh
        + 1.14325367e-07 * t2m2 * va * e_mrt2_rh
        + 2.53016723e-06 * va2 * e_mrt2_rh
        + -1.72857035e-08 * t2m * va2 * e_mrt2_rh
        + -3.95079398e-08 * va3 * e_mrt2_rh
        + -3.59413173e-07 * e_mrt3_rh
        + 7.04388046e-07 * t2m * e_mrt3_rh
        + -1.89309167e-08 * t2m2 * e_mrt3_rh
        + -4.79768731e-07 * va * e_mrt3_rh
        + 7.96079978e-09 * t2m_va * e_mrt3_rh
        + 1.62897058e-09 * va2 * e_mrt3_rh
        + 3.94367674e-08 * e_mrt4_rh
        + -1.18566247e-09 * t2m * e_mrt4_rh
        + 3.34678041e-10 * va * e_mrt4_rh
        + -1.15606447e-10 * e_mrt5 * wvp
        + -2.80626406e00 * wvp2
        + 5.48712484e-01 * t2m * wvp2
        + -3.99428410e-03 * t2m2 * wvp2
        + -9.54009191e-04 * t2m3 * wvp2
        + 1.93090978e-05 * t2m4 * wvp2
        + -3.08806365e-01 * varh2
        + 1.16952364e-02 * t2m * varh2
        + 4.95271903e-04 * t2m2 * varh2
        + -1.90710882e-05 * t2m3 * varh2
        + 2.10787756e-03 * va2 * wvp2
        + -6.98445738e-04 * t2m * va2 * wvp2
        + 2.30109073e-05 * t2m2 * va2 * wvp2
        + 4.17856590e-04 * va3 * wvp2
        + -1.27043871e-05 * t2m * va3 * wvp2
        + -3.04620472e-06 * va4 * wvp2
        + 5.14507424e-02 * e_mrt_rh2
        + -4.32510997e-03 * t2m * e_mrt_rh2
        + 8.99281156e-05 * t2m2 * e_mrt_rh2
        + -7.14663943e-07 * t2m3 * e_mrt_rh2
        + -2.66016305e-04 * va * e_mrt_rh2
        + 2.63789586e-04 * t2m_va * e_mrt_rh2
        + -7.01199003e-06 * t2m2 * va * e_mrt_rh2
        + -1.06823306e-04 * va2 * e_mrt_rh2
        + 3.61341136e-06 * t2m * va2 * e_mrt_rh2
        + 2.29748967e-07 * va3 * e_mrt_rh2
        + 3.04788893e-04 * e_mrt2_rh2
        + -6.42070836e-05 * t2m * e_mrt2_rh2
        + 1.16257971e-06 * t2m2 * e_mrt2_rh2
        + 7.68023384e-06 * va * e_mrt2_rh2
        + -5.47446896e-07 * t2m_va * e_mrt2_rh2
        + -3.59937910e-08 * va2 * e_mrt2_rh2
        + -4.36497725e-06 * e_mrt3 * wvp2
        + 1.68737969e-07 * t2m * e_mrt3 * wvp2
        + 2.67489271e-08 * va * e_mrt3 * wvp2
        + 3.23926897e-09 * e_mrt4 * wvp2
        + -3.53874123e-02 * wvp3
        + -2.21201190e-01 * t2m * wvp3
        + 1.55126038e-02 * t2m2 * wvp3
        + -2.63917279e-04 * t2m3 * wvp3
        + 4.53433455e-02 * va * wvp3
        + -4.32943862e-03 * t2m_va * wvp3
        + 1.45389826e-04 * t2m2 * va * wvp3
        + 2.17508610e-04 * va2 * wvp3
        + -6.66724702e-05 * t2m * va2 * wvp3
        + 3.33217140e-05 * va3 * wvp3
        + -2.26921615e-03 * e_mrt_rh3
        + 3.80261982e-04 * t2m * e_mrt_rh3
        + -5.45314314e-09 * t2m2 * e_mrt_rh3
        + -7.96355448e-04 * va * e_mrt_rh3
        + 2.53458034e-05 * t2m_va * e_mrt_rh3
        + -6.31223658e-06 * va2 * e_mrt_rh3
        + 3.02122035e-04 * e_mrt2 * wvp3
        + -4.77403547e-06 * t2m * e_mrt2 * wvp3
        + 1.73825715e-06 * va * e_mrt2 * wvp3
        + -4.09087898e-07 * e_mrt3 * wvp3
        + 6.14155345e-01 * wvp4
        + -6.16755931e-02 * t2m * wvp4
        + 1.33374846e-03 * t2m2 * wvp4
        + 3.55375387e-03 * va * wvp4
        + -5.13027851e-04 * t2m_va * wvp4
        + 1.02449757e-04 * va2 * wvp4
        + -1.48526421e-03 * e_mrt * wvp4
        + -4.11469183e-05 * t2m * e_mrt * wvp4
        + -6.80434415e-06 * va * e_mrt * wvp4
        + -9.77675906e-06 * e_mrt2 * wvp4
        + 8.82773108e-02 * wvp5
        + -3.01859306e-03 * t2m * wvp5
        + 1.04452989e-03 * va * wvp5
        + 2.47090539e-04 * e_mrt * wvp5
        + 1.48348065e-03 * wvp6
    )

    return utci

calculate_utci

calculate_utci(
    t2_k: ArrayLike,
    va: ArrayLike,
    mrt: ArrayLike,
    td_k: ArrayLike | None = None,
    ehPa: ArrayLike | None = None,
) -> np.ndarray

UTCI - Universal Thermal Climate Index :param t2_k: (float array) is 2m temperature [K] :param va: (float array) is wind speed at 10 meters [m/s] :param mrt: (float array) is mean radiant temperature [K] :param td_k: (float array) is 2m dew point temperature [K] :param ehPa: (float array) is water vapour pressure [hPa] returns UTCI [K]

Validity: the polynomial is fitted for air temperature -50...+50 degC, 10 m wind speed 0.5...17 m/s, mean radiant temperature from 30 degC below to 70 degC above air temperature, and water vapour pressure up to 50 hPa. Inputs are not clamped; outside this range the approximation extrapolates.

Reference: Brode et al. (2012) https://doi.org/10.1007/s00484-011-0454-1

Source code in thermofeel/thermofeel.py

def calculate_utci(
    t2_k: ArrayLike,
    va: ArrayLike,
    mrt: ArrayLike,
    td_k: ArrayLike | None = None,
    ehPa: ArrayLike | None = None,
) -> np.ndarray:
    """
    UTCI - Universal Thermal Climate Index
        :param t2_k: (float array) is 2m temperature [K]
        :param va: (float array) is wind speed at 10 meters [m/s]
        :param mrt: (float array) is mean radiant temperature [K]
        :param td_k: (float array) is 2m dew point temperature [K]
        :param ehPa: (float array) is water vapour pressure [hPa]
    returns UTCI [K]

    Validity: the polynomial is fitted for air temperature -50...+50 degC, 10 m
    wind speed 0.5...17 m/s, mean radiant temperature from 30 degC below to
    70 degC above air temperature, and water vapour pressure up to 50 hPa.
    Inputs are not clamped; outside this range the approximation extrapolates.

    Reference: Brode et al. (2012)
    https://doi.org/10.1007/s00484-011-0454-1
    """

    if ehPa is not None:
        vapour = ehPa
    elif td_k is not None:
        rh_pc = calculate_relative_humidity_percent(t2_k, td_k)
        vapour = calculate_saturation_vapour_pressure(t2_k) * rh_pc / 100.0
    else:
        raise ValueError("Missing input ehPa or td_k")
    wvp = vapour / 10.0  # water vapour pressure in kPa

    t2_c = kelvin_to_celsius(t2_k)  # polynomial approx. is in Celsius
    mrt_c = kelvin_to_celsius(mrt)  # polynomial approx. is in Celsius

    utci = calculate_utci_polynomial(t2_c, mrt_c, va, wvp)
    utci_k = celsius_to_kelvin(utci)

    return utci_k

calculate_wbgt_simple

calculate_wbgt_simple(
    t2_k: ArrayLike, rh: ArrayLike
) -> np.ndarray

WBGT - Wet Bulb Globe Temperature computed by a the simpler algorithm :param t2_k: (float array) 2m temperature [K] :param rh: (float array) relative humidity percentage [%] returns Wet Bulb Globe Temperature [K]

Validity: an empirical screening estimate (no radiation/wind term), intended for moderate-to-warm outdoor conditions; it is not a substitute for the physically based calculate_wbgt_liljegren where accuracy matters.

Reference: ACSM (1984) https://doi.org/10.1080/00913847.1984.11701899 See also: http://www.bom.gov.au/info/thermal_stress/#approximation https://www.jstage.jst.go.jp/article/indhealth/50/4/50_MS1352/_pdf

Source code in thermofeel/thermofeel.py

def calculate_wbgt_simple(t2_k: ArrayLike, rh: ArrayLike) -> np.ndarray:
    """
    WBGT - Wet Bulb Globe Temperature computed by a the simpler algorithm
        :param t2_k: (float array) 2m temperature [K]
        :param rh: (float array) relative humidity percentage [%]
        returns Wet Bulb Globe Temperature [K]

    Validity: an empirical screening estimate (no radiation/wind term), intended
    for moderate-to-warm outdoor conditions; it is not a substitute for the
    physically based ``calculate_wbgt_liljegren`` where accuracy matters.

    Reference: ACSM (1984)
    https://doi.org/10.1080/00913847.1984.11701899
    See also: http://www.bom.gov.au/info/thermal_stress/#approximation
    https://www.jstage.jst.go.jp/article/indhealth/50/4/50_MS1352/_pdf
    """
    t2_c = kelvin_to_celsius(t2_k)
    e = calculate_nonsaturation_vapour_pressure(t2_k, rh)
    wbgt = 0.567 * t2_c + 0.393 * e + 3.94
    wbgt_k = celsius_to_kelvin(wbgt)

    return wbgt_k

calculate_wbt

calculate_wbt(t2_k: ArrayLike, rh: ArrayLike) -> np.ndarray

Wet Bulb Temperature :param t2_k: (float array) 2m temperature [K] :param rh: (float array) relative humidity percentage [%] returns wet bulb temperature [K]

Validity: Stull's regression is fitted for relative humidity 5...99% and air temperature -20...+50 degC at standard sea-level pressure (1013.25 hPa).

Reference: Stull (2011) https://doi.org/10.1175/JAMC-D-11-0143.1

Source code in thermofeel/thermofeel.py

def calculate_wbt(t2_k: ArrayLike, rh: ArrayLike) -> np.ndarray:
    """
    Wet Bulb Temperature
        :param t2_k: (float array) 2m temperature [K]
        :param rh: (float array) relative humidity percentage [%]
        returns wet bulb temperature [K]

    Validity: Stull's regression is fitted for relative humidity 5...99% and air
    temperature -20...+50 degC at standard sea-level pressure (1013.25 hPa).

    Reference: Stull (2011)
    https://doi.org/10.1175/JAMC-D-11-0143.1
    """
    t2_c = kelvin_to_celsius(t2_k)
    tw = (
        t2_c * np.arctan(0.151977 * np.sqrt(rh + 8.313659))
        + np.arctan(t2_c + rh)
        - np.arctan(rh - 1.676331)
        + 0.00391838 * (rh) ** (3 / 2) * np.arctan(0.023101 * rh)
        - 4.686035
    )
    tw_k = celsius_to_kelvin(tw)

    return tw_k

calculate_bgt

calculate_bgt(
    t2_k: ArrayLike, mrt: ArrayLike, va: ArrayLike
) -> np.ndarray

Globe temperature :param t2_k: (float array) 2m temperature [K] :param mrt: (float array) mean radiant temperature [K] :param va: (float array) wind speed at 10 meters [m/s] returns globe temperature [K]

Solves the globe energy balance with a closed-form quartic root.

Edge case — calm air: the closed form is a 0/0 indeterminate at exactly zero wind speed. The limit as va -> 0 is the mean radiant temperature (with no convection the globe reaches radiative equilibrium, bgt -> mrt), so mrt is returned where va == 0. Invalid negative wind still yields NaN.

Reference: Guo et al. 2018 https://doi.org/10.1016/j.enbuild.2018.08.029

Source code in thermofeel/thermofeel.py

def calculate_bgt(t2_k: ArrayLike, mrt: ArrayLike, va: ArrayLike) -> np.ndarray:
    """
    Globe temperature
        :param t2_k: (float array) 2m temperature [K]
        :param mrt: (float array) mean radiant temperature [K]
        :param va: (float array) wind speed at 10 meters [m/s]
        returns globe temperature [K]

    Solves the globe energy balance with a closed-form quartic root.

    Edge case — calm air: the closed form is a ``0/0`` indeterminate at exactly
    zero wind speed. The limit as ``va -> 0`` is the mean radiant temperature
    (with no convection the globe reaches radiative equilibrium, ``bgt -> mrt``),
    so ``mrt`` is returned where ``va == 0``. Invalid negative wind still yields
    ``NaN``.

    Reference: Guo et al. 2018
    https://doi.org/10.1016/j.enbuild.2018.08.029
    """
    v = scale_windspeed(
        va, 1.1
    )  # formula requires wind speed at 1.1m (i.e., at the level of the globe)

    # a = 1
    d = (1.1e8 * v**0.6) / (0.95 * 0.15**0.4)
    e = -(mrt**4) - d * t2_k

    q = 12 * e
    s = 27 * (d**2)
    # At v == 0 the closed form evaluates 0/0; guard the warnings and substitute
    # the analytic calm-air limit (bgt -> mrt) below.
    with np.errstate(invalid="ignore", divide="ignore"):
        delta = ((s + np.sqrt(s**2 - 4 * (q**3))) / 2) ** (1 / 3)
        Q = 0.5 * np.sqrt((1 / 3) * (delta + q / delta))
        bgt = -Q + 0.5 * np.sqrt(-4 * (Q**2) + d / Q)

    # Calm-air limit: no convection -> globe at radiative equilibrium (= mrt).
    bgt = np.where(v == 0.0, np.asarray(mrt, dtype=float), bgt)

    return bgt

calculate_wbgt

calculate_wbgt(
    t2_k: ArrayLike,
    mrt: ArrayLike,
    va: ArrayLike,
    td_k: ArrayLike,
) -> np.ndarray

WBGT - Wet Bulb Globe Temperature :param t2_k: (float array) 2m temperature [K] :param mrt: (float array) mean radiant temperature [K] :param va: (float array) wind speed at 10 meters [m/s] :param td_k: (float array) dew point temperature [K] returns wet bulb globe temperature [K] Reference: Stull (2011) https://doi.org/10.1175/JAMC-D-11-0143.1 See also: http://www.bom.gov.au/info/thermal_stress/

Source code in thermofeel/thermofeel.py

def calculate_wbgt(
    t2_k: ArrayLike, mrt: ArrayLike, va: ArrayLike, td_k: ArrayLike
) -> np.ndarray:
    """
    WBGT - Wet Bulb Globe Temperature
        :param t2_k: (float array) 2m temperature [K]
        :param mrt: (float array) mean radiant temperature [K]
        :param va: (float array) wind speed at 10 meters [m/s]
        :param td_k: (float array) dew point temperature [K]
        returns wet bulb globe temperature [K]
    Reference: Stull (2011)
    https://doi.org/10.1175/JAMC-D-11-0143.1
    See also: http://www.bom.gov.au/info/thermal_stress/
    """

    bgt_k = calculate_bgt(t2_k, mrt, va)
    bgt_c = kelvin_to_celsius(bgt_k)

    rh = calculate_relative_humidity_percent(t2_k, td_k)
    t2_c = kelvin_to_celsius(t2_k)
    tw_k = calculate_wbt(t2_k, rh)
    tw_c = kelvin_to_celsius(tw_k)

    wbgt = 0.7 * tw_c + 0.2 * bgt_c + 0.1 * t2_c
    wbgt_k = celsius_to_kelvin(wbgt)

    return wbgt_k

calculate_wind_speed_2m_liljegren

calculate_wind_speed_2m_liljegren(
    va: ArrayLike, cossza: ArrayLike, ssrd: ArrayLike
) -> np.ndarray

Scale 10m wind speed to 2m using the Liljegren stability-dependent profile :param va: (float array) wind speed at 10 metres [m/s] :param cossza: (float array) cosine of the solar zenith angle [dimensionless] :param ssrd: (float array) instantaneous downward shortwave radiation at the surface (SSRD) [W/m2] returns wind speed at 2 metres [m/s]

This is the 10m -> 2m conversion used operationally by KNMI within the Liljegren WBGT (va * (2/10)**p), where the power-law exponent p comes from a Pasquill-Gifford stability class derived from the solar elevation, incoming radiation and wind speed. The result is floored at 0.13 m/s.

Reference: Liljegren et al. (2008) https://doi.org/10.1080/15459620802310770 See also: Kong and Huber (2022) https://doi.org/10.1029/2021EF002334

Source code in thermofeel/thermofeel.py

def calculate_wind_speed_2m_liljegren(
    va: ArrayLike, cossza: ArrayLike, ssrd: ArrayLike
) -> np.ndarray:
    """
    Scale 10m wind speed to 2m using the Liljegren stability-dependent profile
        :param va: (float array) wind speed at 10 metres [m/s]
        :param cossza: (float array) cosine of the solar zenith angle
            [dimensionless]
        :param ssrd: (float array) instantaneous downward shortwave
            radiation at the surface (SSRD) [W/m2]
        returns wind speed at 2 metres [m/s]

    This is the 10m -> 2m conversion used operationally by KNMI within the
    Liljegren WBGT (``va * (2/10)**p``), where the power-law exponent ``p`` comes
    from a Pasquill-Gifford stability class derived from the solar elevation,
    incoming radiation and wind speed. The result is floored at 0.13 m/s.

    Reference: Liljegren et al. (2008)
    https://doi.org/10.1080/15459620802310770
    See also: Kong and Huber (2022) https://doi.org/10.1029/2021EF002334
    """
    return _liljegren_wind_speed_2m(va, cossza, ssrd)

calculate_wbgt_liljegren

calculate_wbgt_liljegren(
    t2_k: ArrayLike,
    rh: ArrayLike,
    pressure: ArrayLike,
    va: ArrayLike,
    ssrd: ArrayLike,
    fdir: ArrayLike,
    cossza: ArrayLike,
    wind_scaling: str = "liljegren",
) -> np.ndarray

WBGT - Wet Bulb Globe Temperature, Liljegren method :param t2_k: (float array) 2m temperature [K] :param rh: (float array) relative humidity [%] :param pressure: (float array) surface air pressure [hPa] :param va: (float array) wind speed at 10 metres [m/s] :param ssrd: (float array) instantaneous downward shortwave radiation at the surface (SSRD) [W/m2] :param fdir: (float array) fraction of ssrd that is direct beam [dimensionless, 0-1] :param cossza: (float array) cosine of the solar zenith angle [dimensionless] :param wind_scaling: (str) how to convert the 10m wind to the 2m sensor height: "liljegren" (default, the KNMI stability-dependent profile, see calculate_wind_speed_2m_liljegren) or "brode" (the generic scale_windspeed log profile) returns wet bulb globe temperature [K]

Physically based WBGT after Liljegren et al. (2008), the method used operationally by KNMI. The globe and natural-wet-bulb temperatures are each solved from their steady-state energy balance by fixed-point iteration and combined with the air temperature as WBGT = 0.7*Tnw + 0.2*Tg + 0.1*Ta.

The KNMI operational guards are applied: wind below 0.62 m/s at 10 m is raised to that floor (it is then scaled to 2 m for the sensor model); the direct-beam fraction is clamped to [0, 0.9] and set to 0 when the sun is at or below 89.5 degrees zenith. Instantaneous ssrd and cossza are supplied by the caller (e.g. via earthkit-meteo). NaN is returned where the iteration does not converge.

The Liljegren energy-balance solvers and physical constants live in the thermofeel.liljegren submodule.

Reference: Liljegren et al. (2008) https://doi.org/10.1080/15459620802310770 See also: Kong and Huber (2022) https://doi.org/10.1029/2021EF002334

Source code in thermofeel/thermofeel.py

def calculate_wbgt_liljegren(
    t2_k: ArrayLike,
    rh: ArrayLike,
    pressure: ArrayLike,
    va: ArrayLike,
    ssrd: ArrayLike,
    fdir: ArrayLike,
    cossza: ArrayLike,
    wind_scaling: str = "liljegren",
) -> np.ndarray:
    """
    WBGT - Wet Bulb Globe Temperature, Liljegren method
        :param t2_k: (float array) 2m temperature [K]
        :param rh: (float array) relative humidity [%]
        :param pressure: (float array) surface air pressure [hPa]
        :param va: (float array) wind speed at 10 metres [m/s]
        :param ssrd: (float array) instantaneous downward shortwave
            radiation at the surface (SSRD) [W/m2]
        :param fdir: (float array) fraction of ssrd that is direct beam
            [dimensionless, 0-1]
        :param cossza: (float array) cosine of the solar zenith angle
            [dimensionless]
        :param wind_scaling: (str) how to convert the 10m wind to the 2m sensor
            height: "liljegren" (default, the KNMI stability-dependent profile,
            see calculate_wind_speed_2m_liljegren) or "brode" (the generic
            scale_windspeed log profile)
        returns wet bulb globe temperature [K]

    Physically based WBGT after Liljegren et al. (2008), the method used
    operationally by KNMI. The globe and natural-wet-bulb temperatures are each
    solved from their steady-state energy balance by fixed-point iteration and
    combined with the air temperature as ``WBGT = 0.7*Tnw + 0.2*Tg + 0.1*Ta``.

    The KNMI operational guards are applied: wind below 0.62 m/s at 10 m is
    raised to that floor (it is then scaled to 2 m for the sensor model); the
    direct-beam fraction is clamped to [0, 0.9] and set to 0 when the sun is at
    or below 89.5 degrees zenith. Instantaneous ``ssrd`` and ``cossza`` are
    supplied by the caller (e.g. via earthkit-meteo). NaN is returned where the
    iteration does not converge.

    The Liljegren energy-balance solvers and physical constants live in the
    ``thermofeel.liljegren`` submodule.

    Reference: Liljegren et al. (2008)
    https://doi.org/10.1080/15459620802310770
    See also: Kong and Huber (2022) https://doi.org/10.1029/2021EF002334
    """
    va = np.asarray(va, dtype=float)

    # KNMI minimum wind floor at 10 m, then scale to the 2 m sensor height.
    va = np.maximum(va, _LILJEGREN_MIN_WIND_10M)
    if wind_scaling == "liljegren":
        speed = _liljegren_wind_speed_2m(va, cossza, ssrd)
    elif wind_scaling == "brode":
        speed = scale_windspeed(va, 2.0)
    else:
        raise ValueError(
            f"wind_scaling must be 'liljegren' or 'brode', got {wind_scaling!r}"
        )

    return _liljegren_wbgt(t2_k, rh, pressure, speed, ssrd, fdir, cossza)

calculate_heat_force

calculate_heat_force(wbgt_k: ArrayLike) -> np.ndarray

Heat Force - KNMI 0-10 heat-stress communication scale (hittekracht) :param wbgt_k: (float array) wet bulb globe temperature [K] returns heat force on a 0-10 scale [dimensionless]

Translates WBGT onto the integer 0-10 "heat force" scale used by KNMI for public communication of heat stress, analogous to wind force and the UV index. The bands are fixed 2 degC intervals of WBGT (lower-closed): heat force 0 is WBGT < 14 degC, 1 is [14, 16) degC, ..., 9 is [30, 32) degC, and 10 is WBGT >= 32 degC. Values are returned as whole numbers (as a float array, so NaN inputs propagate).

Reference: KNMI Technical Report TR-26-04 (2026), Table 1.

Source code in thermofeel/thermofeel.py

def calculate_heat_force(wbgt_k: ArrayLike) -> np.ndarray:
    """
    Heat Force - KNMI 0-10 heat-stress communication scale (hittekracht)
        :param wbgt_k: (float array) wet bulb globe temperature [K]
        returns heat force on a 0-10 scale [dimensionless]

    Translates WBGT onto the integer 0-10 "heat force" scale used by KNMI for
    public communication of heat stress, analogous to wind force and the UV
    index. The bands are fixed 2 degC intervals of WBGT (lower-closed): heat
    force 0 is WBGT < 14 degC, 1 is [14, 16) degC, ..., 9 is [30, 32) degC, and
    10 is WBGT >= 32 degC. Values are returned as whole numbers (as a float
    array, so NaN inputs propagate).

    Reference: KNMI Technical Report TR-26-04 (2026), Table 1.
    """
    wbgt_c = kelvin_to_celsius(np.asarray(wbgt_k, dtype=float))
    return np.clip(np.floor((wbgt_c - 14.0) / 2.0) + 1.0, 0.0, 10.0)

calculate_excess_heat_factor

calculate_excess_heat_factor(
    ehi_sig: ArrayLike,
    ehi_accl: ArrayLike,
    clip: bool = False,
) -> np.ndarray

EXHF - Excess Heat Factor :param ehi_sig: (float array) significance index, e.g. relative to the 95th percentile of daily mean temperature over a reference period. :param ehi_accl: (float array) acclimatisation index. :param clip: (bool) clip the significance index at zero (default False). returns excess heat factor [input unit squared, e.g. K^2]

Top-level convenience wrapper for thermofeel.excess_heat.excess_heat_factor. The supporting daily mean temperature, significance and acclimatisation indices, and heatwave severity, live in the thermofeel.excess_heat submodule.

Reference: Nairn and Fawcett (2014), equation (3) https://doi.org/10.3390/ijerph120100227

Source code in thermofeel/thermofeel.py

def calculate_excess_heat_factor(
    ehi_sig: ArrayLike, ehi_accl: ArrayLike, clip: bool = False
) -> np.ndarray:
    """
    EXHF - Excess Heat Factor
        :param ehi_sig: (float array) significance index, e.g. relative to the
            95th percentile of daily mean temperature over a reference period.
        :param ehi_accl: (float array) acclimatisation index.
        :param clip: (bool) clip the significance index at zero (default False).
        returns excess heat factor [input unit squared, e.g. K^2]

    Top-level convenience wrapper for
    ``thermofeel.excess_heat.excess_heat_factor``. The supporting daily mean
    temperature, significance and acclimatisation indices, and heatwave
    severity, live in the ``thermofeel.excess_heat`` submodule.

    Reference: Nairn and Fawcett (2014), equation (3)
    https://doi.org/10.3390/ijerph120100227
    """
    return _excess_heat_factor(ehi_sig, ehi_accl, clip=clip)

calculate_excess_cold_factor

calculate_excess_cold_factor(
    ehi_sig: ArrayLike,
    ehi_accl: ArrayLike,
    clip: bool = False,
) -> np.ndarray

EXCF - Excess Cold Factor :param ehi_sig: (float array) significance index, e.g. relative to the 5th percentile of daily mean temperature over a reference period. :param ehi_accl: (float array) acclimatisation index. :param clip: (bool) clip the significance index at zero (default False). returns excess cold factor [input unit squared, e.g. K^2]

Top-level convenience wrapper for thermofeel.excess_heat.excess_cold_factor. The supporting daily mean temperature, significance and acclimatisation indices live in the thermofeel.excess_heat submodule.

Reference: Nairn (2013), equation (7) https://www.cawcr.gov.au/technical-reports/CTR_060.pdf

Source code in thermofeel/thermofeel.py

def calculate_excess_cold_factor(
    ehi_sig: ArrayLike, ehi_accl: ArrayLike, clip: bool = False
) -> np.ndarray:
    """
    EXCF - Excess Cold Factor
        :param ehi_sig: (float array) significance index, e.g. relative to the
            5th percentile of daily mean temperature over a reference period.
        :param ehi_accl: (float array) acclimatisation index.
        :param clip: (bool) clip the significance index at zero (default False).
        returns excess cold factor [input unit squared, e.g. K^2]

    Top-level convenience wrapper for
    ``thermofeel.excess_heat.excess_cold_factor``. The supporting daily mean
    temperature, significance and acclimatisation indices live in the
    ``thermofeel.excess_heat`` submodule.

    Reference: Nairn (2013), equation (7)
    https://www.cawcr.gov.au/technical-reports/CTR_060.pdf
    """
    return _excess_cold_factor(ehi_sig, ehi_accl, clip=clip)

calculate_mrt_from_bgt

calculate_mrt_from_bgt(
    t2_k: ArrayLike, bgt_k: ArrayLike, va: ArrayLike
) -> np.ndarray

Mean radiant temperature from globe temperature :param t2_k: (float array) 2m temperature [K] :param bgt_k: (float array) globe temperature [K] :param va: (float array) wind speed at 10 meters [m/s] returns mean radiant temperature [K] Reference: Brimicombe et al. (2023) https://doi.org/10.1029/2022GH000701

Source code in thermofeel/thermofeel.py

def calculate_mrt_from_bgt(
    t2_k: ArrayLike, bgt_k: ArrayLike, va: ArrayLike
) -> np.ndarray:
    """
    Mean radiant temperature from globe temperature
        :param t2_k: (float array) 2m temperature [K]
        :param bgt_k: (float array) globe temperature [K]
        :param va: (float array) wind speed at 10 meters [m/s]
        returns mean radiant temperature [K]
    Reference: Brimicombe et al. (2023)
    https://doi.org/10.1029/2022GH000701
    """
    v = scale_windspeed(
        va, 1.1
    )  # formula requires wind speed at 1.1m (i.e., at the level of the globe)
    f = (1.1e8 * v**0.6) / (0.95 * 0.15**0.4)
    bgt4 = bgt_k**4
    mrtc = bgt4 + f * (bgt_k - t2_k)
    mrtc2 = np.sqrt(np.sqrt(mrtc))

    return mrtc2

calculate_humidex

calculate_humidex(
    t2_k: ArrayLike, td_k: ArrayLike
) -> np.ndarray

Humidex
:param t2_k: (float array) 2m temperature [K]
:param td_k: (float array) dew point temperature [K]
returns humidex [K]

Validity: Environment Canada's index; it is most meaningful in warm, humid conditions and is not reported by Environment Canada below ~20 degC.

Reference: Environment Canada https://climate.weather.gc.ca/glossary_e.html#humidex

Source code in thermofeel/thermofeel.py

def calculate_humidex(t2_k: ArrayLike, td_k: ArrayLike) -> np.ndarray:
    """
        Humidex
        :param t2_k: (float array) 2m temperature [K]
        :param td_k: (float array) dew point temperature [K]
        returns humidex [K]

    Validity: Environment Canada's index; it is most meaningful in warm, humid
    conditions and is not reported by Environment Canada below ~20 degC.

    Reference: Environment Canada
    https://climate.weather.gc.ca/glossary_e.html#humidex
    """
    vp = 6.11 * np.exp(5417.7530 * ((1 / 273.15) - (1 / td_k)))  # vapour pressure [hPa]
    h = 0.5555 * (vp - 10.0)
    humidex = t2_k + h

    return humidex

calculate_discomfort_index

calculate_discomfort_index(
    t2_k: ArrayLike, rh: ArrayLike
) -> np.ndarray

Discomfort Index (Thom) :param t2_k: (float array) 2m temperature [K] :param rh: (float array) relative humidity [%] returns discomfort index [K]

Thom's Discomfort Index (also called the Temperature-Humidity Index) estimates heat discomfort from air temperature and relative humidity. This is the Celsius/relative-humidity formulation of Thom's index given by Giles et al. (1990), DI = T - 0.55 (1 - 0.01 RH)(T - 14.5), with T in degC and RH in %; at 100% RH the index equals the air temperature, and in warm conditions (above 14.5 degC) drier air lowers it below the air temperature. It is a warm-season heat-stress indicator and is not clamped - out-of-range inputs return the raw value (the caller masks).

Reference (this temperature/relative-humidity formulation): Giles et al. (1990) https://doi.org/10.1007/BF01093455 Index origin (Thom's Discomfort Index / Temperature-Humidity Index): Thom (1959) https://doi.org/10.1080/00431672.1959.9926960 See also (modern review and discomfort categories; a different wet-bulb formulation): Epstein and Moran (2006) https://doi.org/10.2486/indhealth.44.388

Source code in thermofeel/thermofeel.py

def calculate_discomfort_index(t2_k: ArrayLike, rh: ArrayLike) -> np.ndarray:
    """
    Discomfort Index (Thom)
        :param t2_k: (float array) 2m temperature [K]
        :param rh: (float array) relative humidity [%]
        returns discomfort index [K]

    Thom's Discomfort Index (also called the Temperature-Humidity Index)
    estimates heat discomfort from air temperature and relative humidity. This is
    the Celsius/relative-humidity formulation of Thom's index given by Giles et
    al. (1990), DI = T - 0.55 (1 - 0.01 RH)(T - 14.5), with T in degC and RH in
    %; at 100% RH the index equals the air temperature, and in warm conditions
    (above 14.5 degC) drier air lowers it below the air temperature. It is a
    warm-season heat-stress indicator and is not clamped - out-of-range inputs
    return the raw value (the caller masks).

    Reference (this temperature/relative-humidity formulation): Giles et al.
    (1990) https://doi.org/10.1007/BF01093455
    Index origin (Thom's Discomfort Index / Temperature-Humidity Index):
    Thom (1959) https://doi.org/10.1080/00431672.1959.9926960
    See also (modern review and discomfort categories; a different wet-bulb
    formulation): Epstein and Moran (2006)
    https://doi.org/10.2486/indhealth.44.388
    """
    t2_c = kelvin_to_celsius(t2_k)
    di = t2_c - 0.55 * (1.0 - 0.01 * rh) * (t2_c - 14.5)
    di_k = celsius_to_kelvin(di)

    return di_k

calculate_summer_simmer_index

calculate_summer_simmer_index(
    t2_k: ArrayLike, rh: ArrayLike
) -> np.ndarray

Summer Simmer Index (Pepi) :param t2_k: (float array) 2m temperature [K] :param rh: (float array) relative humidity [%] returns summer simmer index [K]

The Summer Simmer Index (SSI) estimates warm-season heat discomfort from air temperature and relative humidity. This is the common 1987 closed form, evaluated in Fahrenheit, SSI = 1.98 (Tf - (0.55 - 0.0055 RH)(Tf - 58)) - 56.83, with Tf the air temperature in degF and RH in %, the result converted back to Kelvin. The inner bracket is Thom's Temperature-Humidity Index in Fahrenheit, so the SSI is an affine image of that index (SSI_F = 1.98 THI_F - 56.83), the near-exact Fahrenheit sibling of calculate_discomfort_index (which uses 14.5 degC where 58 degF = 14.44... degC). It is a warm-season heat-stress indicator and is not clamped - out-of-range inputs return the raw value (the caller masks).

Reference: Pepi, J.W. (1987) The Summer Simmer Index, Weatherwise 40(3): 143-145 https://doi.org/10.1080/00431672.1987.9933356

Provenance caveat: the 1987 article is not openly available, so the equation above is reproduced from secondary sources. It is the common 1987 closed form - an affine transform of Thom's Fahrenheit Temperature-Humidity Index - NOT the author's later tabulated "New Summer Simmer Index", which is a different relationship and is not implemented here.

Source code in thermofeel/thermofeel.py

def calculate_summer_simmer_index(t2_k: ArrayLike, rh: ArrayLike) -> np.ndarray:
    """
    Summer Simmer Index (Pepi)
        :param t2_k: (float array) 2m temperature [K]
        :param rh: (float array) relative humidity [%]
        returns summer simmer index [K]

    The Summer Simmer Index (SSI) estimates warm-season heat discomfort from air
    temperature and relative humidity. This is the common 1987 closed form,
    evaluated in Fahrenheit,
    SSI = 1.98 (Tf - (0.55 - 0.0055 RH)(Tf - 58)) - 56.83, with Tf the air
    temperature in degF and RH in %, the result converted back to Kelvin. The
    inner bracket is Thom's Temperature-Humidity Index in Fahrenheit, so the SSI
    is an affine image of that index (SSI_F = 1.98 THI_F - 56.83), the near-exact
    Fahrenheit sibling of calculate_discomfort_index (which uses 14.5 degC where
    58 degF = 14.44... degC). It is a warm-season heat-stress indicator and is
    not clamped - out-of-range inputs return the raw value (the caller masks).

    Reference: Pepi, J.W. (1987) The Summer Simmer Index, Weatherwise 40(3):
    143-145 https://doi.org/10.1080/00431672.1987.9933356

    Provenance caveat: the 1987 article is not openly available, so the equation
    above is reproduced from secondary sources. It is the *common* 1987 closed
    form - an affine transform of Thom's Fahrenheit Temperature-Humidity Index -
    NOT the author's later tabulated "New Summer Simmer Index", which is a
    different relationship and is not implemented here.
    """
    t2_f = kelvin_to_fahrenheit(t2_k)
    thi_f = t2_f - (0.55 - 0.0055 * rh) * (t2_f - 58.0)
    ssi_f = 1.98 * thi_f - 56.83
    ssi_k = fahrenheit_to_kelvin(ssi_f)

    return ssi_k

calculate_normal_effective_temperature

calculate_normal_effective_temperature(
    t2_k: ArrayLike, va: ArrayLike, rh: ArrayLike
) -> np.ndarray

NET - Normal Effective Temperature :param t2_k: (float array) 2m temperature [K] :param va: (float array) wind speed at 10 meters [m/s] :param rh: (float array) relative humidity percentage [%] returns normal effective temperature [K]

Validity: an empirical index derived for subtropical (Hong Kong) conditions (Li and Chan 2000); it has no sharply defined input range.

Reference: Li and Chan (2006) https://doi.org/10.1017/S1350482700001602

Source code in thermofeel/thermofeel.py

def calculate_normal_effective_temperature(
    t2_k: ArrayLike, va: ArrayLike, rh: ArrayLike
) -> np.ndarray:
    """
    NET - Normal Effective Temperature
        :param t2_k: (float array) 2m temperature [K]
        :param va: (float array) wind speed at 10 meters [m/s]
        :param rh: (float array) relative humidity percentage [%]
        returns normal effective temperature [K]

    Validity: an empirical index derived for subtropical (Hong Kong) conditions
    (Li and Chan 2000); it has no sharply defined input range.

    Reference: Li and Chan (2006)
    https://doi.org/10.1017/S1350482700001602
    """
    t2_c = kelvin_to_celsius(t2_k)
    v = scale_windspeed(va, 1.2)  # formula requires wind speed at 1.2m
    ditermeq = 1 / (1.76 + 1.4 * v**0.75)
    net = (
        37
        - ((37 - t2_c) / (0.68 - 0.0014 * rh + ditermeq))
        - 0.29 * t2_c * (1 - 0.01 * rh)
    )
    net_k = celsius_to_kelvin(net)

    return net_k

calculate_relative_strain_index

calculate_relative_strain_index(
    t2_k: ArrayLike, rh: ArrayLike
) -> np.ndarray

Relative Strain Index (RSI) :param t2_k: (float array) 2m temperature [K] :param rh: (float array) relative humidity percentage [%] returns relative strain index [dimensionless]

The Relative Strain Index expresses warm-environment heat strain on a young, healthy adult from air temperature and ambient water-vapour pressure. This is the peer-reviewed hectopascal closed form stated with units by Asghari et al. (2020), RSI = (Ta - 21) / (58 - e), with Ta the air temperature in degC and e the ambient water-vapour pressure in hPa (from calculate_nonsaturation_vapour_pressure). It is dimensionless and not clamped - out-of-range inputs return the raw value (the caller masks).

Domain edge: the denominator 58 - e shrinks as e approaches 58 hPa (near saturation around 35-36 degC), so RSI grows without bound (diverging to +/-inf), and for e > 58 hPa (very hot and near-saturated, beyond the ~35 degC validity range) the denominator turns negative and RSI becomes negative - a spurious value for a heat-strain index. These out-of-domain elements are returned raw and are deliberately not clamped; the caller masks inputs above the validity range.

Variant caveat: a different literature form, (10.7 + 0.74 (Ta - 35)) / (44 - Pa), appears in secondary sources with Pa in other units (likely mmHg). It could not be verified against the primary text and is deliberately NOT implemented here; only this hPa closed form is provided.

Reference (this hPa closed form and the five-level assessment bands): Asghari et al. (2020) https://doi.org/10.2174/1874213002013010011 Index origin (Relative Strain Index): Lee and Henschel (1966) https://doi.org/10.1111/j.1749-6632.1966.tb43059.x

Source code in thermofeel/thermofeel.py

def calculate_relative_strain_index(t2_k: ArrayLike, rh: ArrayLike) -> np.ndarray:
    """
    Relative Strain Index (RSI)
        :param t2_k: (float array) 2m temperature [K]
        :param rh: (float array) relative humidity percentage [%]
        returns relative strain index [dimensionless]

    The Relative Strain Index expresses warm-environment heat strain on a young,
    healthy adult from air temperature and ambient water-vapour pressure. This is
    the peer-reviewed hectopascal closed form stated with units by Asghari et al.
    (2020), RSI = (Ta - 21) / (58 - e), with Ta the air temperature in degC and e
    the ambient water-vapour pressure in hPa (from
    ``calculate_nonsaturation_vapour_pressure``). It is dimensionless and not
    clamped - out-of-range inputs return the raw value (the caller masks).

    Domain edge: the denominator ``58 - e`` shrinks as e approaches 58 hPa (near
    saturation around 35-36 degC), so RSI grows without bound (diverging to
    +/-inf), and for e > 58 hPa (very hot and near-saturated, beyond the ~35 degC
    validity range) the denominator turns negative and RSI becomes negative - a
    spurious value for a heat-strain index. These out-of-domain elements are
    returned raw and are deliberately not clamped; the caller masks inputs above
    the validity range.

    Variant caveat: a different literature form, (10.7 + 0.74 (Ta - 35)) / (44 -
    Pa), appears in secondary sources with Pa in other units (likely mmHg). It
    could not be verified against the primary text and is deliberately NOT
    implemented here; only this hPa closed form is provided.

    Reference (this hPa closed form and the five-level assessment bands): Asghari
    et al. (2020) https://doi.org/10.2174/1874213002013010011
    Index origin (Relative Strain Index): Lee and Henschel (1966)
    https://doi.org/10.1111/j.1749-6632.1966.tb43059.x
    """
    t2_c = kelvin_to_celsius(t2_k)
    e = calculate_nonsaturation_vapour_pressure(t2_k, rh)
    rsi = (t2_c - 21.0) / (58.0 - e)

    return rsi

calculate_apparent_temperature

calculate_apparent_temperature(
    t2_k: ArrayLike, va: ArrayLike, rh: ArrayLike
) -> np.ndarray

Apparent Temperature - version without radiation :param t2_k: (float array) 2m temperature [K] :param va: (float array) wind speed at 10 meters [m/s] :param rh: (float array) relative humidity percentage [%] returns apparent temperature [K]

Validity: the Bureau of Meteorology non-radiation form of Steadman's apparent temperature; an empirical estimate without sharply defined input bounds.

Reference: Steadman (1984) https://doi.org/10.1175/1520-0450(1984)023%3C1674:AUSOAT%3E2.0.CO;2 See also: http://www.bom.gov.au/info/thermal_stress/#atapproximation

Source code in thermofeel/thermofeel.py

def calculate_apparent_temperature(
    t2_k: ArrayLike, va: ArrayLike, rh: ArrayLike
) -> np.ndarray:
    """
    Apparent Temperature - version without radiation
        :param t2_k: (float array) 2m temperature [K]
        :param va: (float array) wind speed at 10 meters [m/s]
        :param rh: (float array) relative humidity percentage [%]
        returns apparent temperature [K]

    Validity: the Bureau of Meteorology non-radiation form of Steadman's apparent
    temperature; an empirical estimate without sharply defined input bounds.

    Reference: Steadman (1984)
    https://doi.org/10.1175/1520-0450(1984)023%3C1674:AUSOAT%3E2.0.CO;2
    See also: http://www.bom.gov.au/info/thermal_stress/#atapproximation
    """
    t2_c = kelvin_to_celsius(t2_k)
    e = calculate_nonsaturation_vapour_pressure(t2_k, rh)
    at = t2_c + 0.33 * e - 0.7 * va - 4
    at_k = celsius_to_kelvin(at)

    return at_k

calculate_apparent_temperature_radiation

calculate_apparent_temperature_radiation(
    t2_k: ArrayLike,
    va: ArrayLike,
    rh: ArrayLike,
    q: ArrayLike,
) -> np.ndarray

Apparent Temperature - version including radiation :param t2_k: (float array) 2m temperature [K] :param va: (float array) wind speed at 10 meters [m/s] :param rh: (float array) relative humidity percentage [%] :param q: (float array) net radiation absorbed per unit body-surface area [W m-2] returns apparent temperature [K]

The radiation-inclusive form of Steadman's apparent temperature, as published operationally by the Australian Bureau of Meteorology: AT = Ta + 0.348*e - 0.70*va + 0.70*q/(va + 10) - 4.25 (Ta in degC), with e the ambient water-vapour pressure in hPa. It shares the vapour path of calculate_apparent_temperature via calculate_nonsaturation_vapour_pressure (whose constants match BoM's e = (rh/100)*6.105*exp(17.27*Ta/(237.7+Ta))).

q is the net radiation absorbed per unit area of body surface. It is a caller-supplied input (like cossza), NOT an NWP surface flux and NOT the mean radiant temperature; callers must supply their own q.

Validity: an empirical estimate without sharply defined input bounds; the result is not clamped.

Reference: Steadman, R.G. (1994) Norms of apparent temperature in Australia, Aust. Met. Mag. 43:1-16. https://doi.org/10.1071/es94001 See also: http://www.bom.gov.au/info/thermal_stress/#atapproximation

Source code in thermofeel/thermofeel.py

def calculate_apparent_temperature_radiation(
    t2_k: ArrayLike, va: ArrayLike, rh: ArrayLike, q: ArrayLike
) -> np.ndarray:
    """
    Apparent Temperature - version including radiation
        :param t2_k: (float array) 2m temperature [K]
        :param va: (float array) wind speed at 10 meters [m/s]
        :param rh: (float array) relative humidity percentage [%]
        :param q: (float array) net radiation absorbed per unit body-surface
            area [W m-2]
        returns apparent temperature [K]

    The radiation-inclusive form of Steadman's apparent temperature, as
    published operationally by the Australian Bureau of Meteorology:
    ``AT = Ta + 0.348*e - 0.70*va + 0.70*q/(va + 10) - 4.25`` (Ta in degC),
    with ``e`` the ambient water-vapour pressure in hPa. It shares the vapour
    path of ``calculate_apparent_temperature`` via
    ``calculate_nonsaturation_vapour_pressure`` (whose constants match BoM's
    ``e = (rh/100)*6.105*exp(17.27*Ta/(237.7+Ta))``).

    ``q`` is the net radiation absorbed per unit area of body surface. It is a
    caller-supplied input (like ``cossza``), NOT an NWP surface flux and NOT the
    mean radiant temperature; callers must supply their own ``q``.

    Validity: an empirical estimate without sharply defined input bounds; the
    result is not clamped.

    Reference: Steadman, R.G. (1994) Norms of apparent temperature in Australia,
    Aust. Met. Mag. 43:1-16.
    https://doi.org/10.1071/es94001
    See also: http://www.bom.gov.au/info/thermal_stress/#atapproximation
    """
    t2_c = kelvin_to_celsius(t2_k)
    e = calculate_nonsaturation_vapour_pressure(t2_k, rh)
    at = t2_c + 0.348 * e - 0.70 * va + 0.70 * q / (va + 10) - 4.25
    at_k = celsius_to_kelvin(at)

    return at_k

calculate_wind_chill

calculate_wind_chill(
    t2_k: ArrayLike, va: ArrayLike
) -> np.ndarray

Wind Chill :param t2_k: (float array) 2m Temperature [K] :param va: (float array) wind speed at 10 meters [m/s] returns wind chill [K] Computation is only valid for temperatures between -50°C and 5°C and wind speeds between 5km/h and 80km/h. For input values outside those ranges, computed results not be considered valid. Reference: Blazejczyk et al. (2012) https://doi.org/10.1007/s00484-011-0453-2 See also: https://web.archive.org/web/20130627223738/http://climate.weatheroffice.gc.ca/prods_servs/normals_documentation_e.html # noqa

Source code in thermofeel/thermofeel.py

def calculate_wind_chill(t2_k: ArrayLike, va: ArrayLike) -> np.ndarray:
    """
    Wind Chill
        :param t2_k: (float array) 2m Temperature [K]
        :param va: (float array) wind speed at 10 meters [m/s]
        returns wind chill [K]
        Computation is only valid for temperatures between -50°C and 5°C and wind speeds between 5km/h and 80km/h.
        For input values outside those ranges, computed results not be considered valid.
    Reference: Blazejczyk et al. (2012)
    https://doi.org/10.1007/s00484-011-0453-2
    See also: https://web.archive.org/web/20130627223738/http://climate.weatheroffice.gc.ca/prods_servs/normals_documentation_e.html  # noqa
    """
    t2_c = kelvin_to_celsius(t2_k)
    v = va * 3.6  # convert to kilometers per hour
    windchill = 13.12 + 0.6215 * t2_c - 11.37 * v**0.16 + 0.3965 * t2_c * v**0.16
    windchill_k = celsius_to_kelvin(windchill)

    return windchill_k

calculate_heat_index_simplified

calculate_heat_index_simplified(
    t2_k: ArrayLike, rh: ArrayLike
) -> np.ndarray

Heat Index :param t2_k: (float array) 2m temperature [K] :param rh: (float array) relative humidity [%] returns heat index [K]

The regression is applied only where air temperature exceeds 20 degC (below that the air temperature is returned unchanged). Inputs are array-like (wrap a scalar in an array), as the function masks elementwise internally.

Reference: Blazejczyk et al. (2012) https://doi.org/10.1007/s00484-011-0453-2

Source code in thermofeel/thermofeel.py

def calculate_heat_index_simplified(t2_k: ArrayLike, rh: ArrayLike) -> np.ndarray:
    """
    Heat Index
        :param t2_k: (float array) 2m temperature [K]
        :param rh: (float array) relative humidity [%]
        returns heat index [K]

    The regression is applied only where air temperature exceeds 20 degC (below
    that the air temperature is returned unchanged). Inputs are array-like (wrap
    a scalar in an array), as the function masks elementwise internally.

    Reference: Blazejczyk et al. (2012)
    https://doi.org/10.1007/s00484-011-0453-2
    """

    t2_c = kelvin_to_celsius(t2_k)

    hiarray = [
        8.784695,
        1.61139411,
        2.338549,
        0.14611605,
        1.2308094e-2,
        1.6424828e-2,
        2.211732e-3,
        7.2546e-4,
        3.582e-6,
    ]
    hi = np.copy(t2_c)

    hi_filter1 = np.where(t2_c > 20)

    hi[hi_filter1] = (
        -hiarray[0]
        + hiarray[1] * t2_c[hi_filter1]
        + hiarray[2] * rh[hi_filter1]
        - hiarray[3] * t2_c[hi_filter1] * rh[hi_filter1]
        - hiarray[4] * t2_c[hi_filter1] ** 2
        - hiarray[5] * rh[hi_filter1] ** 2
        + hiarray[6] * t2_c[hi_filter1] ** 2 * rh[hi_filter1]
        + hiarray[7] * t2_c[hi_filter1] * rh[hi_filter1] ** 2
        - hiarray[8] * t2_c[hi_filter1] ** 2 * rh[hi_filter1] ** 2
    )

    hi_k = celsius_to_kelvin(hi)

    return hi_k

calculate_heat_index_adjusted

calculate_heat_index_adjusted(
    t2_k: ArrayLike, td_k: ArrayLike
) -> np.ndarray

Heat Index adjusted :param t2_k: (float array) 2m temperature [K] :param td_k: (float array) 2m dewpoint temperature [K] returns heat index [K] Reference: https://www.wpc.ncep.noaa.gov/html/heatindex_equation.shtml

Source code in thermofeel/thermofeel.py

def calculate_heat_index_adjusted(t2_k: ArrayLike, td_k: ArrayLike) -> np.ndarray:
    """
    Heat Index adjusted
      :param t2_k: (float array) 2m temperature [K]
      :param td_k: (float array) 2m dewpoint temperature  [K]
      returns heat index [K]
    Reference: https://www.wpc.ncep.noaa.gov/html/heatindex_equation.shtml
    """

    rh = calculate_relative_humidity_percent(t2_k, td_k)
    t2_f = kelvin_to_fahrenheit(t2_k)

    # Simple formula
    A = 0.5 * (t2_f + 61.0 + ((t2_f - 68.0) * 1.2) + (rh * 0.094))

    # Use regression only when the average between A and air temperature is or exceeds 80°F
    warm_mask = ((A + t2_f) / 2) >= 80

    hi = A  # (T+A)/2 < 80F (reverse of warm_mask)

    # Regression formula
    rh2 = rh**2
    t2f2 = t2_f**2
    B = (
        -42.379
        + 2.04901523 * t2_f
        + 10.14333127 * rh
        - 0.22475541 * t2_f * rh
        - 6.83783e-3 * t2f2
        - 5.481717e-2 * rh2
        + 1.22874e-3 * t2f2 * rh
        + 8.5282e-4 * t2_f * rh2
        - 1.99e-6 * t2f2 * rh2
    )

    # Replace with regression where applicable
    hi[warm_mask] = B[warm_mask]

    # Adjustments — only for warm_mask points
    # 1) low humidity and high temperatures
    mask = warm_mask & (t2_f >= 80.0) & (t2_f <= 112.0) & (rh < 13.0)
    lowH_highT = ((13.0 - rh) / 4.0) * np.sqrt(
        np.maximum((17.0 - np.abs(t2_f - 95.0)), 0.0) / 17.0
    )
    hi[mask] -= lowH_highT[mask]

    # 2) High humidity and moderate temperatures
    mask = warm_mask & (t2_f >= 80.0) & (t2_f <= 87.0) & (rh > 85.0)
    # highH_modT = ((rh - 85.0) / 10.0) * ((87.0 - t2_f) / 5.0)
    highH_modT = 0.02 * (rh - 85.0) * (87.0 - t2_f)
    hi[mask] += highH_modT[mask]

    # Cold case (t ≤ 40°F)
    cold_mask = t2_f <= 40.0
    hi[cold_mask] = t2_f[cold_mask]

    # Convert to Kelvin
    hi_k = fahrenheit_to_kelvin(hi)

    return hi_k

calculate_pmv

calculate_pmv(
    t2_k: ArrayLike,
    mrt_k: ArrayLike,
    var: ArrayLike,
    rh: ArrayLike | None = None,
    vapour_pressure_hpa: ArrayLike | None = None,
    met: ArrayLike = 1.2,
    clo: ArrayLike = 0.5,
    wme: ArrayLike = 0.0,
) -> np.ndarray

PMV - Predicted Mean Vote (Fanger's thermal comfort model) :param t2_k: (float array) 2m air temperature [K] :param mrt_k: (float array) mean radiant temperature [K] :param var: (float array) relative air velocity at the body [m/s]. This is the air movement felt at the body (the still-air velocity plus any movement the occupant's own activity induces), NOT the 10 m meteorological wind speed - do not pass the model wind field. :param rh: (float array) relative humidity [%]. Provide exactly one of rh or vapour_pressure_hpa. :param vapour_pressure_hpa: (float array) water vapour pressure [hPa]. Provide exactly one of rh or vapour_pressure_hpa. :param met: (float array) metabolic rate [met] (1 met = 58.15 W m-2); default 1.2 met (ISO 8996 sedentary/office activity). :param clo: (float array) clothing insulation [clo] (1 clo = 0.155 m2 K W-1); default 0.5 clo (light indoor clothing). :param wme: (float array) external (mechanical) work [met]; default 0.0. returns predicted mean vote [dimensionless]. The thermal-sensation scale runs from -3 (cold) through 0 (neutral) to +3 (hot); the index is not clamped, so inputs beyond the ISO 7730 validity range can yield values outside +/-3 (ISO recommends interpreting PMV only within about +/-2).

Fanger's steady-state heat-balance comfort equation as standardised in ISO 7730:2005. The clothing-surface temperature is obtained by the Annex D fixed-point iteration (tolerance 0.00015, capped at 150 sweeps). The whole array is iterated together (mirroring the Liljegren energy-balance solvers) and any element that has not converged within the cap is returned as NaN.

met and clo are exposed as parameters because PMV is only defined for a stated activity and clothing level; the defaults (1.2 met, 0.5 clo) are common office conventions rather than fixed constants of the standard.

Reference: ISO 7730:2005 (Annex D); Fanger, P.O. (1970) Thermal Comfort: Analysis and Applications in Environmental Engineering, McGraw-Hill.

Source code in thermofeel/thermofeel.py

def calculate_pmv(
    t2_k: ArrayLike,
    mrt_k: ArrayLike,
    var: ArrayLike,
    rh: ArrayLike | None = None,
    vapour_pressure_hpa: ArrayLike | None = None,
    met: ArrayLike = 1.2,
    clo: ArrayLike = 0.5,
    wme: ArrayLike = 0.0,
) -> np.ndarray:
    """
    PMV - Predicted Mean Vote (Fanger's thermal comfort model)
        :param t2_k: (float array) 2m air temperature [K]
        :param mrt_k: (float array) mean radiant temperature [K]
        :param var: (float array) relative air velocity at the body [m/s].
            This is the air movement felt at the body (the still-air velocity
            plus any movement the occupant's own activity induces), NOT the
            10 m meteorological wind speed - do not pass the model wind field.
        :param rh: (float array) relative humidity [%]. Provide exactly one of
            ``rh`` or ``vapour_pressure_hpa``.
        :param vapour_pressure_hpa: (float array) water vapour pressure [hPa].
            Provide exactly one of ``rh`` or ``vapour_pressure_hpa``.
        :param met: (float array) metabolic rate [met] (1 met = 58.15 W m-2);
            default 1.2 met (ISO 8996 sedentary/office activity).
        :param clo: (float array) clothing insulation [clo] (1 clo = 0.155
            m2 K W-1); default 0.5 clo (light indoor clothing).
        :param wme: (float array) external (mechanical) work [met]; default 0.0.
        returns predicted mean vote [dimensionless]. The thermal-sensation
        scale runs from -3 (cold) through 0 (neutral) to +3 (hot); the index is
        not clamped, so inputs beyond the ISO 7730 validity range can yield
        values outside +/-3 (ISO recommends interpreting PMV only within
        about +/-2).

    Fanger's steady-state heat-balance comfort equation as standardised in
    ISO 7730:2005. The clothing-surface temperature is obtained by the Annex D
    fixed-point iteration (tolerance 0.00015, capped at 150 sweeps). The whole
    array is iterated together (mirroring the Liljegren energy-balance solvers)
    and any element that has not converged within the cap is returned as NaN.

    ``met`` and ``clo`` are exposed as parameters because PMV is only defined for
    a stated activity and clothing level; the defaults (1.2 met, 0.5 clo) are
    common office conventions rather than fixed constants of the standard.

    Reference: ISO 7730:2005 (Annex D); Fanger, P.O. (1970) Thermal Comfort:
    Analysis and Applications in Environmental Engineering, McGraw-Hill.
    """
    ta = kelvin_to_celsius(np.asarray(t2_k, dtype=float))
    tr = kelvin_to_celsius(np.asarray(mrt_k, dtype=float))
    var = np.asarray(var, dtype=float)
    met = np.asarray(met, dtype=float)
    clo = np.asarray(clo, dtype=float)
    wme = np.asarray(wme, dtype=float)

    # Water vapour partial pressure [Pa]: require exactly one humidity input.
    if rh is not None and vapour_pressure_hpa is None:
        pa = (
            np.asarray(rh, dtype=float)
            * 10.0
            * np.exp(16.6536 - 4030.183 / (ta + 235.0))
        )
    elif vapour_pressure_hpa is not None and rh is None:
        pa = np.asarray(vapour_pressure_hpa, dtype=float) * 100.0
    else:
        raise ValueError("Provide exactly one of rh (%) or vapour_pressure_hpa (hPa)")

    icl = 0.155 * clo  # clothing insulation [m2 K W-1]
    m = 58.15 * met  # metabolic rate [W m-2]
    w = 58.15 * wme  # external work [W m-2]
    mw = m - w  # internal heat production [W m-2]

    fcl = np.where(icl <= 0.078, 1.0 + 1.29 * icl, 1.05 + 0.645 * icl)
    hcf = 12.1 * np.sqrt(var)  # forced convective heat transfer coeff.

    taa = ta + 273.0
    tra = tr + 273.0
    tcla = taa + (35.5 - ta) / (3.5 * icl + 0.1)

    p1 = icl * fcl
    p2 = p1 * 3.96
    p3 = p1 * 100.0
    p4 = p1 * taa
    p5 = 308.7 - 0.028 * mw + p2 * (tra / 100.0) ** 4

    # Vectorised fixed-point iteration for the clothing-surface temperature,
    # mirroring the Liljegren globe/wet-bulb solvers: the whole array is swept
    # together and each element is frozen once it meets the ISO tolerance. The
    # convergence mask replaces the standard's per-point "fail if N > 150".
    eps = 0.00015
    bshape = np.broadcast_shapes(
        np.shape(taa),
        np.shape(hcf),
        np.shape(p2),
        np.shape(p3),
        np.shape(p4),
        np.shape(p5),
        np.shape(tcla),
    )
    xn = np.broadcast_to(tcla / 100.0, bshape)
    xf = xn
    hc = np.zeros(bshape, dtype=float)
    converged = np.zeros(bshape, dtype=bool)
    for _ in range(150):
        xf_new = (xf + xn) / 2.0
        hcn = 2.38 * np.abs(100.0 * xf_new - taa) ** 0.25
        hc_new = np.maximum(hcf, hcn)
        xn_new = (p5 + p4 * hc_new - p2 * xf_new**4) / (100.0 + p3 * hc_new)
        update = ~converged
        xf = np.where(update, xf_new, xf)
        xn = np.where(update, xn_new, xn)
        hc = np.where(update, hc_new, hc)
        converged = converged | (np.abs(xn_new - xf_new) <= eps)
        if converged.all():
            break

    tcl = 100.0 * xn - 273.0  # clothing-surface temperature [degC]

    # Heat-loss components [W m-2]
    hl1 = 3.05e-3 * (5733.0 - 6.99 * mw - pa)  # skin diffusion
    hl2 = np.where(mw > 58.15, 0.42 * (mw - 58.15), 0.0)  # sweating
    hl3 = 1.7e-5 * m * (5867.0 - pa)  # latent respiration
    hl4 = 0.0014 * m * (34.0 - ta)  # dry respiration
    hl5 = 3.96e-8 * fcl * ((tcl + 273.0) ** 4 - (tr + 273.0) ** 4)  # radiation
    hl6 = fcl * hc * (tcl - ta)  # convection

    ts = 0.303 * np.exp(-0.036 * m) + 0.028  # thermal-sensation transfer coeff.
    pmv = ts * (mw - hl1 - hl2 - hl3 - hl4 - hl5 - hl6)

    return np.where(converged, pmv, np.nan)

calculate_ppd

calculate_ppd(pmv: ArrayLike) -> np.ndarray

PPD - Predicted Percentage of Dissatisfied :param pmv: (float array) predicted mean vote [dimensionless] returns predicted percentage of dissatisfied [%]

ISO 7730:2005 / Fanger relation mapping the Predicted Mean Vote onto the percentage of occupants likely to be thermally dissatisfied. It is symmetric in PMV (warm and cold deviations of equal magnitude give the same PPD) with a minimum of 5% at PMV = 0 (thermal neutrality). A NaN PMV (e.g. a non-converged calculate_pmv element) propagates to NaN.

Reference: ISO 7730:2005 (Annex D); Fanger, P.O. (1970) Thermal Comfort: Analysis and Applications in Environmental Engineering, McGraw-Hill.

Source code in thermofeel/thermofeel.py

def calculate_ppd(pmv: ArrayLike) -> np.ndarray:
    """
    PPD - Predicted Percentage of Dissatisfied
        :param pmv: (float array) predicted mean vote [dimensionless]
        returns predicted percentage of dissatisfied [%]

    ISO 7730:2005 / Fanger relation mapping the Predicted Mean Vote onto the
    percentage of occupants likely to be thermally dissatisfied. It is symmetric
    in PMV (warm and cold deviations of equal magnitude give the same PPD) with a
    minimum of 5% at PMV = 0 (thermal neutrality). A NaN PMV (e.g. a non-converged
    ``calculate_pmv`` element) propagates to NaN.

    Reference: ISO 7730:2005 (Annex D); Fanger, P.O. (1970) Thermal Comfort:
    Analysis and Applications in Environmental Engineering, McGraw-Hill.
    """
    pmv = np.asarray(pmv, dtype=float)
    return 100.0 - 95.0 * np.exp(-0.03353 * pmv**4 - 0.2179 * pmv**2)

Excess heat and cold factors

The Excess Heat Factor and Excess Cold Factor are available at the top level as thermofeel.calculate_excess_heat_factor(...) and thermofeel.calculate_excess_cold_factor(...) (listed above). The thermofeel.excess_heat submodule below additionally provides the supporting daily mean temperature, significance and acclimatisation indices, and heatwave severity.

daily_mean_temperature

daily_mean_temperature(
    t2_min: ArrayLike, t2_max: ArrayLike
) -> np.ndarray

Daily mean temperature computed from min and max

:param t2_min: 2-metre daily minimum temperature.
:param t2_max: 2-metre daily maximum temperature.

Reference: Nairn and Fawcett (2014) https://doi.org/10.3390/ijerph120100227

Source code in thermofeel/excess_heat.py

def daily_mean_temperature(t2_min: ArrayLike, t2_max: ArrayLike) -> np.ndarray:
    """
    Daily mean temperature computed from min and max

        :param t2_min: 2-metre daily minimum temperature.
        :param t2_max: 2-metre daily maximum temperature.

    Reference: Nairn and Fawcett (2014)
    https://doi.org/10.3390/ijerph120100227
    """
    return 0.5 * (t2_min + t2_max)

significance_index

significance_index(
    dmt: ArrayLike, threshold: ArrayLike
) -> np.ndarray

Significance index

:param dmt: 3-day running mean daily mean temperature (day i, i+1, i+2).
:param threshold: climatological threshold, e.g., 95th (heat extremes)
    or 5th (cold extremes) percentile of daily mean temperature over a
    reference period.

Reference: Nairn and Fawcett (2014), equation (1) https://doi.org/10.3390/ijerph120100227

Source code in thermofeel/excess_heat.py

def significance_index(dmt: ArrayLike, threshold: ArrayLike) -> np.ndarray:
    """
    Significance index

        :param dmt: 3-day running mean daily mean temperature (day i, i+1, i+2).
        :param threshold: climatological threshold, e.g., 95th (heat extremes)
            or 5th (cold extremes) percentile of daily mean temperature over a
            reference period.

    Reference: Nairn and Fawcett (2014), equation (1)
    https://doi.org/10.3390/ijerph120100227
    """
    return dmt - threshold

acclimatisation_index

acclimatisation_index(
    dmt: ArrayLike, threshold: ArrayLike
) -> np.ndarray

Acclimatisation index

:param dmt: 3-day running mean daily mean temperature (days i, i+1, i+2).
:param threshold: 30-day running mean daily mean temperature (days i-30, ..., i-1).

Reference: Nairn and Fawcett (2014), equation (2) https://doi.org/10.3390/ijerph120100227

Source code in thermofeel/excess_heat.py

def acclimatisation_index(dmt: ArrayLike, threshold: ArrayLike) -> np.ndarray:
    """
    Acclimatisation index

        :param dmt: 3-day running mean daily mean temperature (days i, i+1, i+2).
        :param threshold: 30-day running mean daily mean temperature (days i-30, ..., i-1).

    Reference: Nairn and Fawcett (2014), equation (2)
    https://doi.org/10.3390/ijerph120100227
    """
    return dmt - threshold

excess_heat_factor

excess_heat_factor(
    ehi_sig: ArrayLike,
    ehi_accl: ArrayLike,
    clip: bool = False,
) -> np.ndarray

Excess heat factor

:param ehi_sig: Significance index, e.g., computed with respect to 95th
    percentile of daily mean temperature over a reference period.
:param ehi_accl: Acclimatisation index.
:param clip: Whether to clip the lower value range at zero. Disabled by default.

Reference: Nairn and Fawcett (2014), equation (3) https://doi.org/10.3390/ijerph120100227

Source code in thermofeel/excess_heat.py

def excess_heat_factor(
    ehi_sig: ArrayLike, ehi_accl: ArrayLike, clip: bool = False
) -> np.ndarray:
    """
    Excess heat factor

        :param ehi_sig: Significance index, e.g., computed with respect to 95th
            percentile of daily mean temperature over a reference period.
        :param ehi_accl: Acclimatisation index.
        :param clip: Whether to clip the lower value range at zero. Disabled by default.

    Reference: Nairn and Fawcett (2014), equation (3)
    https://doi.org/10.3390/ijerph120100227
    """
    if clip:
        ehi_sig = np.maximum(0, ehi_sig)
    return ehi_sig * np.maximum(1.0, ehi_accl)

heatwave_severity

heatwave_severity(
    exhf: ArrayLike, threshold: ArrayLike
) -> np.ndarray

Heatwave severity index

:param exhf: Excess heat factor.
:param threshold: Excess heat factor threshold, the 85th percentile
    of all positive excess heat factor values in a reference period.

Reference: Nairn (2018), equation (4) https://doi.org/10.3390/ijerph15112494

Source code in thermofeel/excess_heat.py

def heatwave_severity(exhf: ArrayLike, threshold: ArrayLike) -> np.ndarray:
    """
    Heatwave severity index

        :param exhf: Excess heat factor.
        :param threshold: Excess heat factor threshold, the 85th percentile
            of all positive excess heat factor values in a reference period.

    Reference: Nairn (2018), equation (4)
    https://doi.org/10.3390/ijerph15112494
    """
    return exhf / threshold

excess_cold_factor

excess_cold_factor(
    ehi_sig: ArrayLike,
    ehi_accl: ArrayLike,
    clip: bool = False,
) -> np.ndarray

Excess cold factor

:param ehi_sig: Significance index, e.g., computed with respect to 5th
    percentile of daily mean temperature over a reference period..
:param ehi_accl: Acclimatisation index.
:param clip: Whether to clip the upper value range at zero. Disabled by default.

Reference: Nairn (2013), equation (7) https://www.cawcr.gov.au/technical-reports/CTR_060.pdf

Source code in thermofeel/excess_heat.py

def excess_cold_factor(
    ehi_sig: ArrayLike, ehi_accl: ArrayLike, clip: bool = False
) -> np.ndarray:
    """
    Excess cold factor

        :param ehi_sig: Significance index, e.g., computed with respect to 5th
            percentile of daily mean temperature over a reference period..
        :param ehi_accl: Acclimatisation index.
        :param clip: Whether to clip the upper value range at zero. Disabled by default.

    Reference: Nairn (2013), equation (7)
    https://www.cawcr.gov.au/technical-reports/CTR_060.pdf
    """
    if clip:
        ehi_sig = np.minimum(0, ehi_sig)
    return -ehi_sig * np.minimum(-1.0, ehi_accl)

Approximations

The thermofeel.approximations submodule provides clearly-labelled estimators (not validated index formulas), kept out of the top-level namespace so calling code always shows the .approximations. marker. It currently estimates the direct solar radiation fdir (needed for the mean radiant temperature, and hence UTCI / WBGT / PMV) from the global horizontal radiation ssrd and solar geometry, for datasets that lack fdir. These are demonstration/estimation grade — see the Approximations guide.

Approximations used to estimate an input that a thermal index needs but that a dataset does not provide.

These are clearly-labelled empirical estimates, deliberately kept separate from the validated index formulas in thermofeel.thermofeel and NOT re-exported at the top level: calling code carries the .approximations. marker so it is always visible that an approximation is in play. Their outputs are demonstration / estimation grade, not observations; use a dataset that provides the real field for quantitative work.

Currently provides estimators of the direct (beam) solar radiation on a horizontal surface -- ECMWF fdir -- from the global horizontal radiation ssrd and the solar geometry, for datasets (e.g. ECMWF open data) that carry ssrd but not fdir. fdir is required by calculate_mean_radiant_temperature and hence by UTCI / WBGT / PMV. The direct/diffuse split cannot be reconstructed exactly from ssrd (two skies with the same ssrd can have very different beam fractions), so these use an empirical decomposition driven by the clearness index; expect errors of order 0.1-0.25 of the global radiation, larger in broken cloud and at low sun.

References

• Erbs, D.G., Klein, S.A. & Duffie, J.A. (1982). Estimation of the diffuse radiation fraction for hourly, daily and monthly-average global radiation. Solar Energy 28(4), 293-302. https://doi.org/10.1016/0038-092X(82)90302-4
• Maxwell, E.L. (1987). A Quasi-Physical Model for Converting Hourly Global Horizontal to Direct Normal Insolation. SERI/TR-215-3087, Solar Energy Research Institute, Golden, CO.
• Spencer, J.W. (1971). Fourier series representation of the position of the Sun. Search 2(5), 172. (Earth-Sun distance factor.)
• Kasten, F. (1966). A new table and approximation formula for the relative optical air mass. Arch. Meteorol. Geophys. Bioklimatol. B 14, 206-223.

earth_sun_distance_factor

earth_sun_distance_factor(doy: ArrayLike) -> np.ndarray

Earth-Sun distance (eccentricity) factor (mean-distance normal irradiance scaling), via the Spencer (1971) Fourier series.

:param doy: (float array) day of the year (1..365/366)
returns the dimensionless factor (a0/a)^2, ranging ~0.967 (aphelion,
early July) to ~1.035 (perihelion, early January)

Multiply the solar constant by this factor to obtain the actual extraterrestrial normal irradiance for the day.

Reference: Spencer (1971); see also Iqbal, An Introduction to Solar Radiation (1983).

Source code in thermofeel/approximations.py

def earth_sun_distance_factor(doy: ArrayLike) -> np.ndarray:
    """
    Earth-Sun distance (eccentricity) factor (mean-distance normal irradiance
    scaling), via the Spencer (1971) Fourier series.

        :param doy: (float array) day of the year (1..365/366)
        returns the dimensionless factor (a0/a)^2, ranging ~0.967 (aphelion,
        early July) to ~1.035 (perihelion, early January)

    Multiply the solar constant by this factor to obtain the actual
    extraterrestrial normal irradiance for the day.

    Reference: Spencer (1971); see also Iqbal, An Introduction to Solar
    Radiation (1983).
    """
    doy = np.asarray(doy, dtype=float)
    b = 2.0 * np.pi * (doy - 1.0) / 365.0
    return (
        1.00011
        + 0.034221 * np.cos(b)
        + 0.00128 * np.sin(b)
        + 0.000719 * np.cos(2.0 * b)
        + 7.7e-05 * np.sin(2.0 * b)
    )

approximate_fdir_erbs

approximate_fdir_erbs(
    ssrd: ArrayLike,
    cossza: ArrayLike,
    *,
    doy: ArrayLike | None = None,
    solar_constant: float = SOLAR_CONSTANT,
    min_cossza: float = 0.065
) -> np.ndarray

Estimate direct horizontal solar radiation (ECMWF fdir) from global horizontal ssrd and cos(solar zenith angle) via the Erbs et al. (1982) diffuse-fraction decomposition.

:param ssrd: (float array) surface global horizontal solar radiation
    (ECMWF ssrd), a mean flux [W m-2] (or any energy unit; the output
    follows the input unit)
:param cossza: (float array) cosine of the solar zenith angle
:param doy: (optional, float array) day of the year; if given, the
    extraterrestrial irradiance is corrected for the Earth-Sun distance
    (``earth_sun_distance_factor``); if omitted the mean-distance solar
    constant is used
:param solar_constant: (float) solar constant [W m-2] (default 1361)
:param min_cossza: (float) below this cos(zenith) the Sun is treated as
    below the horizon; ``fdir`` returns 0 there (default 0.065)
returns the estimated direct horizontal radiation, in the same unit as
``ssrd``, clipped to [0, ssrd] (and to 0 where ``ssrd`` is negative)

The Erbs correlation splits global horizontal radiation into its direct and diffuse parts from the clearness index kt = ssrd / TOA. It is an empirical hourly fit; applying it to longer step-mean fluxes is a further approximation. NaN inputs propagate to NaN. Estimation grade only -- use a real fdir field for quantitative work.

Reference: Erbs, Klein & Duffie (1982), Solar Energy 28(4):293-302, https://doi.org/10.1016/0038-092X(82)90302-4

Source code in thermofeel/approximations.py

def approximate_fdir_erbs(
    ssrd: ArrayLike,
    cossza: ArrayLike,
    *,
    doy: ArrayLike | None = None,
    solar_constant: float = SOLAR_CONSTANT,
    min_cossza: float = 0.065,
) -> np.ndarray:
    """
    Estimate direct horizontal solar radiation (ECMWF ``fdir``) from global
    horizontal ``ssrd`` and cos(solar zenith angle) via the Erbs et al. (1982)
    diffuse-fraction decomposition.

        :param ssrd: (float array) surface global horizontal solar radiation
            (ECMWF ssrd), a mean flux [W m-2] (or any energy unit; the output
            follows the input unit)
        :param cossza: (float array) cosine of the solar zenith angle
        :param doy: (optional, float array) day of the year; if given, the
            extraterrestrial irradiance is corrected for the Earth-Sun distance
            (``earth_sun_distance_factor``); if omitted the mean-distance solar
            constant is used
        :param solar_constant: (float) solar constant [W m-2] (default 1361)
        :param min_cossza: (float) below this cos(zenith) the Sun is treated as
            below the horizon; ``fdir`` returns 0 there (default 0.065)
        returns the estimated direct horizontal radiation, in the same unit as
        ``ssrd``, clipped to [0, ssrd] (and to 0 where ``ssrd`` is negative)

    The Erbs correlation splits global horizontal radiation into its direct and
    diffuse parts from the clearness index ``kt = ssrd / TOA``. It is an
    empirical hourly fit; applying it to longer step-mean fluxes is a further
    approximation. NaN inputs propagate to NaN. Estimation grade only -- use a
    real ``fdir`` field for quantitative work.

    Reference: Erbs, Klein & Duffie (1982), Solar Energy 28(4):293-302,
    https://doi.org/10.1016/0038-092X(82)90302-4
    """
    ssrd = np.asarray(ssrd, dtype=float)
    cossza = np.asarray(cossza, dtype=float)
    e0 = earth_sun_distance_factor(doy) if doy is not None else 1.0
    ssrd, cossza, e0 = np.broadcast_arrays(ssrd, cossza, np.asarray(e0, dtype=float))
    finite = np.isfinite(ssrd) & np.isfinite(cossza) & np.isfinite(e0)

    toa = solar_constant * e0 * np.maximum(cossza, min_cossza)
    day = (cossza > min_cossza) & (ssrd > 0)
    with np.errstate(invalid="ignore", divide="ignore"):
        kt = np.clip(np.where(day, ssrd / toa, 0.0), 0.0, 1.0)
    kd = np.where(
        kt <= 0.22,
        1.0 - 0.09 * kt,
        np.where(
            kt <= 0.80,
            0.9511 - 0.1604 * kt + 4.388 * kt**2 - 16.638 * kt**3 + 12.336 * kt**4,
            0.165,
        ),
    )
    fdir = np.clip(np.where(day, ssrd * (1.0 - kd), 0.0), 0.0, np.maximum(ssrd, 0.0))
    return np.where(finite, fdir, np.nan)

approximate_fdir_disc

approximate_fdir_disc(
    ssrd: ArrayLike,
    cossza: ArrayLike,
    doy: ArrayLike,
    *,
    pressure_hpa: ArrayLike = 1013.25,
    solar_constant: float = SOLAR_CONSTANT,
    min_cossza: float = 0.065,
    max_airmass: float = 12.0
) -> np.ndarray

Estimate direct horizontal solar radiation (ECMWF fdir) from global horizontal ssrd and the solar geometry via the DISC model (Maxwell, 1987) -- a quasi-physical direct-normal estimator that generally improves on the simpler Erbs split.

:param ssrd: (float array) surface global horizontal solar radiation
    (ECMWF ssrd), a mean flux [W m-2]
:param cossza: (float array) cosine of the solar zenith angle
:param doy: (float array) day of the year (required; sets the
    Earth-Sun distance correction of the extraterrestrial irradiance)
:param pressure_hpa: (float array) surface pressure [hPa] for the
    pressure-corrected air mass (default 1013.25)
:param solar_constant: (float) solar constant [W m-2] (default 1361).
    Maxwell derived the DISC coefficients with a solar constant of
    ~1367-1370 W m-2, so pass ``solar_constant=1370`` to reproduce the
    original model (and pvlib) exactly.
:param min_cossza: (float) below this cos(zenith) the Sun is treated as
    below the horizon; ``fdir`` returns 0 there (default 0.065)
:param max_airmass: (float) air mass is capped at this value in the Kn
    fit, per the original model (default 12)
returns the estimated direct horizontal radiation [W m-2], clipped to
[0, ssrd] (and to 0 where ``ssrd`` is negative)

DISC converts the global clearness index to a direct-beam clearness index Kn using empirical polynomials in kt and the (pressure-corrected Kasten 1966) air mass, then dni = Kn * I0 and fdir = dni * cossza. NaN inputs propagate to NaN. Estimation grade only.

Reference: Maxwell, E.L. (1987), A Quasi-Physical Model for Converting Hourly Global Horizontal to Direct Normal Insolation, SERI/TR-215-3087. The coefficients are those of the DISC model as implemented in pvlib (pvlib.irradiance.disc).

Source code in thermofeel/approximations.py

def approximate_fdir_disc(
    ssrd: ArrayLike,
    cossza: ArrayLike,
    doy: ArrayLike,
    *,
    pressure_hpa: ArrayLike = 1013.25,
    solar_constant: float = SOLAR_CONSTANT,
    min_cossza: float = 0.065,
    max_airmass: float = 12.0,
) -> np.ndarray:
    """
    Estimate direct horizontal solar radiation (ECMWF ``fdir``) from global
    horizontal ``ssrd`` and the solar geometry via the DISC model (Maxwell,
    1987) -- a quasi-physical direct-normal estimator that generally improves on
    the simpler Erbs split.

        :param ssrd: (float array) surface global horizontal solar radiation
            (ECMWF ssrd), a mean flux [W m-2]
        :param cossza: (float array) cosine of the solar zenith angle
        :param doy: (float array) day of the year (required; sets the
            Earth-Sun distance correction of the extraterrestrial irradiance)
        :param pressure_hpa: (float array) surface pressure [hPa] for the
            pressure-corrected air mass (default 1013.25)
        :param solar_constant: (float) solar constant [W m-2] (default 1361).
            Maxwell derived the DISC coefficients with a solar constant of
            ~1367-1370 W m-2, so pass ``solar_constant=1370`` to reproduce the
            original model (and pvlib) exactly.
        :param min_cossza: (float) below this cos(zenith) the Sun is treated as
            below the horizon; ``fdir`` returns 0 there (default 0.065)
        :param max_airmass: (float) air mass is capped at this value in the Kn
            fit, per the original model (default 12)
        returns the estimated direct horizontal radiation [W m-2], clipped to
        [0, ssrd] (and to 0 where ``ssrd`` is negative)

    DISC converts the global clearness index to a direct-beam clearness index
    ``Kn`` using empirical polynomials in ``kt`` and the (pressure-corrected
    Kasten 1966) air mass, then ``dni = Kn * I0`` and ``fdir = dni * cossza``.
    NaN inputs propagate to NaN. Estimation grade only.

    Reference: Maxwell, E.L. (1987), A Quasi-Physical Model for Converting Hourly
    Global Horizontal to Direct Normal Insolation, SERI/TR-215-3087. The
    coefficients are those of the DISC model as implemented in pvlib
    (``pvlib.irradiance.disc``).
    """
    ssrd = np.asarray(ssrd, dtype=float)
    cossza = np.asarray(cossza, dtype=float)
    pressure_hpa = np.asarray(pressure_hpa, dtype=float)
    doy = np.asarray(doy, dtype=float)
    ssrd, cossza, pressure_hpa, doy = np.broadcast_arrays(
        ssrd, cossza, pressure_hpa, doy
    )
    finite = (
        np.isfinite(ssrd)
        & np.isfinite(cossza)
        & np.isfinite(pressure_hpa)
        & np.isfinite(doy)
    )

    i0 = solar_constant * earth_sun_distance_factor(doy)  # extraterrestrial normal
    day = (cossza > min_cossza) & (ssrd > 0)
    with np.errstate(invalid="ignore", divide="ignore"):
        kt = np.clip(
            np.where(day, ssrd / (i0 * np.maximum(cossza, min_cossza)), 0.0), 0.0, 1.0
        )
        am = _relative_airmass_kasten1966(cossza) * pressure_hpa / 1013.25
    am = np.minimum(am, max_airmass)

    cloudy = kt <= 0.6
    a = np.where(
        cloudy,
        0.512 + kt * (-1.56 + kt * (2.286 - 2.222 * kt)),
        -5.743 + kt * (21.77 + kt * (-27.49 + 11.56 * kt)),
    )
    b = np.where(
        cloudy,
        0.37 + 0.962 * kt,
        41.4 + kt * (-118.5 + kt * (66.05 + 31.9 * kt)),
    )
    c = np.where(
        cloudy,
        -0.28 + kt * (0.932 - 2.048 * kt),
        -47.01 + kt * (184.2 + kt * (-222.0 + 73.81 * kt)),
    )
    with np.errstate(over="ignore"):
        delta_kn = a + b * np.exp(c * am)
    knc = 0.866 + am * (-0.122 + am * (0.0121 + am * (-0.000653 + 1.4e-05 * am)))
    kn = knc - delta_kn
    dni = np.where(day & (kn > 0), kn * i0, 0.0)
    fdir = np.clip(dni * cossza, 0.0, np.maximum(ssrd, 0.0))
    return np.where(finite, fdir, np.nan)

Unit converters