Skip to content

Commit

Permalink
Merge pull request #39 from ecmwf/feature/rename-utci-variable
Browse files Browse the repository at this point in the history 
Rename utci variable
  • Loading branch information
@cladinapoli
cladinapoli authored Nov 29, 2024
2 parents 2b3ef66 + 0f8c204 commit 620e949
Showing 1 changed file with 86 additions and 76 deletions.
162 changes: 86 additions & 76 deletions thermofeel/thermofeel.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -220,7 +220,17 @@ def calculate_mean_radiant_temperature(ssrd, ssr, dsrp, strd, fdir, strr, cossza
return mrt


def calculate_utci_polynomial(t2m, mrt, va, rh):
def calculate_utci_polynomial(t2m, mrt, va, wvp):
"""
Helper function to calculate the UTCI polynomial approximation
:param t2_k: (float array) is 2m temperature [K]
:param mrt: (float array) is mean radiant temperature [K]
: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
Expand All @@ -241,26 +251,26 @@ def calculate_utci_polynomial(t2m, mrt, va, rh):
e_mrt5 = e_mrt4 * e_mrt
e_mrt6 = e_mrt5 * e_mrt

rh2 = rh * rh
rh3 = rh2 * rh
rh4 = rh3 * rh
rh5 = rh4 * rh
rh6 = rh5 * rh
wvp2 = wvp * wvp
wvp3 = wvp2 * wvp
wvp4 = wvp3 * wvp
wvp5 = wvp4 * wvp
wvp6 = wvp5 * wvp

varh2 = va * rh2
va2_rh = va2 * rh
varh2 = va * wvp2
va2_rh = va2 * wvp
va2_e_mrt = va2 * e_mrt
e_mrt_rh = e_mrt * rh
e_mrt_rh2 = e_mrt * rh2
e_mrt2_rh = e_mrt2 * rh
e_mrt2_rh2 = e_mrt2 * rh2
e_mrt_rh3 = e_mrt * rh3
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 * rh
va_rh = va * wvp
t2m_va = t2m * va
e_mrt3_rh = e_mrt3 * rh
e_mrt4_rh = e_mrt4 * rh
e_mrt3_rh = e_mrt3 * wvp
e_mrt4_rh = e_mrt4 * wvp

utci = (
t2m
Expand Down Expand Up @@ -348,12 +358,12 @@ def calculate_utci_polynomial(t2m, mrt, va, rh):
+ 4.03863260e-13 * t2m * e_mrt5
+ 1.95087203e-12 * va * e_mrt5
+ -4.73602469e-12 * e_mrt6
+ 5.12733497e00 * rh
+ -3.12788561e-01 * t2m * rh
+ -1.96701861e-02 * t2m2 * rh
+ 9.99690870e-04 * t2m3 * rh
+ 9.51738512e-06 * t2m4 * rh
+ -4.66426341e-07 * t2m5 * rh
+ 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
Expand All @@ -363,12 +373,12 @@ def calculate_utci_polynomial(t2m, mrt, va, rh):
+ 5.00845667e-03 * t2m * va2_rh
+ 1.00601257e-06 * t2m2 * va2_rh
+ -1.81748644e-06 * t2m3 * va2_rh
+ -1.25813502e-03 * va3 * rh
+ -1.79330391e-04 * t2m * va3 * rh
+ 2.34994441e-06 * t2m2 * va3 * rh
+ 1.29735808e-04 * va4 * rh
+ 1.29064870e-06 * t2m * va4 * rh
+ -2.28558686e-06 * va5 * 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
Expand Down Expand Up @@ -403,22 +413,22 @@ def calculate_utci_polynomial(t2m, mrt, va, 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 * rh
+ -2.80626406e00 * rh2
+ 5.48712484e-01 * t2m * rh2
+ -3.99428410e-03 * t2m2 * rh2
+ -9.54009191e-04 * t2m3 * rh2
+ 1.93090978e-05 * t2m4 * rh2
+ -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 * rh2
+ -6.98445738e-04 * t2m * va2 * rh2
+ 2.30109073e-05 * t2m2 * va2 * rh2
+ 4.17856590e-04 * va3 * rh2
+ -1.27043871e-05 * t2m * va3 * rh2
+ -3.04620472e-06 * va4 * rh2
+ 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
Expand All @@ -435,45 +445,45 @@ def calculate_utci_polynomial(t2m, mrt, va, rh):
+ 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 * rh2
+ 1.68737969e-07 * t2m * e_mrt3 * rh2
+ 2.67489271e-08 * va * e_mrt3 * rh2
+ 3.23926897e-09 * e_mrt4 * rh2
+ -3.53874123e-02 * rh3
+ -2.21201190e-01 * t2m * rh3
+ 1.55126038e-02 * t2m2 * rh3
+ -2.63917279e-04 * t2m3 * rh3
+ 4.53433455e-02 * va * rh3
+ -4.32943862e-03 * t2m_va * rh3
+ 1.45389826e-04 * t2m2 * va * rh3
+ 2.17508610e-04 * va2 * rh3
+ -6.66724702e-05 * t2m * va2 * rh3
+ 3.33217140e-05 * va3 * rh3
+ -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 * rh3
+ -4.77403547e-06 * t2m * e_mrt2 * rh3
+ 1.73825715e-06 * va * e_mrt2 * rh3
+ -4.09087898e-07 * e_mrt3 * rh3
+ 6.14155345e-01 * rh4
+ -6.16755931e-02 * t2m * rh4
+ 1.33374846e-03 * t2m2 * rh4
+ 3.55375387e-03 * va * rh4
+ -5.13027851e-04 * t2m_va * rh4
+ 1.02449757e-04 * va2 * rh4
+ -1.48526421e-03 * e_mrt * rh4
+ -4.11469183e-05 * t2m * e_mrt * rh4
+ -6.80434415e-06 * va * e_mrt * rh4
+ -9.77675906e-06 * e_mrt2 * rh4
+ 8.82773108e-02 * rh5
+ -3.01859306e-03 * t2m * rh5
+ 1.04452989e-03 * va * rh5
+ 2.47090539e-04 * e_mrt * rh5
+ 1.48348065e-03 * rh6
+ 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
Expand All @@ -493,19 +503,19 @@ def calculate_utci(t2_k, va, mrt, td_k=None, ehPa=None):
"""

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

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, kPa)
utci = calculate_utci_polynomial(t2_c, mrt_c, va, wvp)
utci_k = celsius_to_kelvin(utci)

return utci_k
Expand Down

0 comments on commit 620e949

Please sign in to comment.