Coverage for src/CSET/operators/temperature.py: 100%
121 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-07 15:12 +0000
1# © Crown copyright, Met Office (2022-2026) and CSET contributors.
2#
3# Licensed under the Apache License, Version 2.0 (the "License");
4# you may not use this file except in compliance with the License.
5# You may obtain a copy of the License at
6#
7# http://www.apache.org/licenses/LICENSE-2.0
8#
9# Unless required by applicable law or agreed to in writing, software
10# distributed under the License is distributed on an "AS IS" BASIS,
11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12# See the License for the specific language governing permissions and
13# limitations under the License.
15"""Operators for temperature conversions."""
17import iris.cube
18import numpy as np
20from CSET._common import iter_maybe
21from CSET.operators._atmospheric_constants import CPD, EPSILON, LV, RV, T0
22from CSET.operators.humidity import (
23 mixing_ratio_from_relative_humidity,
24 saturation_mixing_ratio,
25)
26from CSET.operators.misc import convert_units
27from CSET.operators.pressure import (
28 exner_pressure,
29 vapour_pressure_from_relative_humidity,
30)
33def dewpoint_temperature(
34 temperature: iris.cube.Cube | iris.cube.CubeList,
35 relative_humidity: iris.cube.Cube | iris.cube.CubeList,
36) -> iris.cube.Cube | iris.cube.CubeList:
37 r"""Calculate the dewpoint temperature.
39 Arguments
40 ---------
41 temperature: iris.cube.Cube | iris.cube.CubeList
42 Cubes of temperature in Kelvin.
43 relative_humidity: iris.cube.Cube | iris.cube.CubeList
44 Cubes of relative humidity.
46 Returns
47 -------
48 iris.cube.Cube | iris.cueb.CubeList
49 Calculated dewpoint temperature in Kelvin.
51 Notes
52 -----
53 The dewpoint temperature provides the temperature at which the air must be
54 cooled for dew to form (i.e. the water vapour condenses to liquid water).
55 It is calculated based upon [Bolton80]_
57 .. math:: T_d = \frac{243.5 * ln(e) - 440.8}{19.48 - ln(e)}
59 for :math:`T_d` the dewpoint temperature and e the vapour pressure. Here the
60 vapour pressure is calculated from relative humidity using
61 `pressure.vapour_pressure_from_relative_humidity`.
63 The dewpoint temperature is presented when -35.0 C < T < 35.0 C as this is
64 when the calculation is the most accurate [Bolton80]_ this roughly equates to
65 exclusion of dewpoints on pressures of 400 hPa and above (e.g. [Flack24]_).
67 When :math:`T_d` is equivalent to T the relative humidity will be 100 %.
69 All cubes must be on the same grid.
71 Examples
72 --------
73 >>> Td = temperature.dewpoint_temperature(T, RH)
74 """
75 Td = iris.cube.CubeList([])
76 for T, RH in zip(
77 iter_maybe(temperature), iter_maybe(relative_humidity), strict=True
78 ):
79 vp = vapour_pressure_from_relative_humidity(T, RH)
80 td = vp.copy()
81 td.data = ((243.5 * np.log(vp.core_data())) - 440.8) / (
82 19.48 - np.log(vp.core_data())
83 )
84 td.data[T.data - T0 < -35.0] = np.nan
85 td.data[T.data - T0 > 35.0] = np.nan
86 td.units = "Celsius"
87 td.rename("dewpoint_temperature")
88 td = convert_units(td, "K")
89 Td.append(td)
90 if len(Td) == 1:
91 return Td[0]
92 else:
93 return Td
96def virtual_temperature(
97 temperature: iris.cube.Cube | iris.cube.CubeList,
98 mixing_ratio: iris.cube.Cube | iris.cube.CubeList,
99) -> iris.cube.Cube | iris.cube.CubeList:
100 r"""Calculate the virtual temperature.
102 Arguments
103 ---------
104 temperature: iris.cube.Cube | iris.cube.CubeList
105 Cubes of temperature in Kelvin.
106 mixing_ratio: iris.cube.Cube | iris.cube.CubeList
107 Cubes of mixing ratio.
109 Returns
110 -------
111 iris.cube.Cube | iris.cube.CubeList
112 Calculated virtual temperature in Kelvin.
114 Notes
115 -----
116 The virtual temperature is that required for a dry parcel to have the
117 same density as a moist parcel at the same pressure. It is calculated
118 as
120 .. math:: T_v = T * \frac{w + \epsilon}{\epsilon (1 + w)}
122 for :math:`T_v` the virtual temperature, T the temperature, w the mixing
123 ratio, and :math:`\epsilon` the ratio between dry and moist air equating
124 to 0.622.
126 T is given in K and w is given in kg kg-1.
127 All cubes must be on the same grid.
129 Examples
130 --------
131 >>> Tv = temperature.virtual_temperature(T, w)
132 """
133 Tv = iris.cube.CubeList([])
134 for T, W in zip(iter_maybe(temperature), iter_maybe(mixing_ratio), strict=True):
135 virT = T * ((W + EPSILON) / (EPSILON * (1 + W)))
136 virT.rename("virtual_temperature")
137 Tv.append(virT)
138 if len(Tv) == 1:
139 return Tv[0]
140 else:
141 return Tv
144def wet_bulb_temperature(
145 temperature: iris.cube.Cube | iris.cube.CubeList,
146 relative_humidity: iris.cube.Cube | iris.cube.CubeList,
147) -> iris.cube.Cube | iris.cube.CubeList:
148 r"""Calculate the wet-bulb temperature.
150 Arguments
151 ---------
152 temperature: iris.cube.Cube | iris.cube.CubeList
153 Cubes of temperature.
154 relative_humidity: iris.cube.Cube | iris.cube.CubeList
155 Cubes of relative humidity.
157 Returns
158 -------
159 iris.cube.Cube | iris.cube.CubeList
160 Calculated wet-bulb temperature in Kelvin.
162 Notes
163 -----
164 The wet-bulb temperature is the temperature the air cools to by
165 evaporating water into it. It can be calculated from temperature in Celsius
166 and relative humidity in percent following [Stull11]_
168 .. math:: T_w = T * arctan\left(0.151977*(RH + 8.313659)^{0.5}\right) + arctan(T + RH) - arctan(RH - 1.676331) + 0.00391838*RH^{\frac{3}{2}}*arctan(0.023101*RH) - 4.686035
170 for :math:`T_w` the wet-bulb temperature, T the temperature, and RH the relative
171 humidity.
173 The equation is valid for 5% < RH < 99% and -20 C < T < 50 C, and results are
174 only presented for these values.
176 The temperature and relative humidity unit conversions are applied in the
177 operator.
179 All cubes should be on the same grid.
181 Examples
182 --------
183 >>> Tw = temperature.wet_bulb_temperature(T, RH)
184 """
185 Tw = iris.cube.CubeList([])
186 for T, RH in zip(
187 iter_maybe(temperature), iter_maybe(relative_humidity), strict=True
188 ):
189 RH = convert_units(RH, "%")
190 T = convert_units(T, "Celsius")
191 wetT = (
192 T * np.arctan(0.151977 * (RH.core_data() + 8.313659) ** 0.5)
193 + np.arctan(T.core_data() + RH.core_data())
194 - np.arctan(RH.core_data() - 1.676331)
195 + 0.00391838
196 * (RH.core_data()) ** (3.0 / 2.0)
197 * np.arctan(0.023101 * RH.core_data())
198 - 4.686035
199 )
200 wetT.rename("wet_bulb_temperature")
201 wetT = convert_units(wetT, "K")
202 wetT.data[T.core_data() < -20.0] = np.nan
203 wetT.data[T.core_data() > 50.0] = np.nan
204 wetT.data[RH.core_data() < 5.0] = np.nan
205 wetT.data[RH.core_data() > 99.0] = np.nan
206 Tw.append(wetT)
207 if len(Tw) == 1:
208 return Tw[0]
209 else:
210 return Tw
213def potential_temperature(
214 temperature: iris.cube.Cube | iris.cube.CubeList,
215 pressure: iris.cube.Cube | iris.cube.CubeList,
216) -> iris.cube.Cube | iris.cube.CubeList:
217 r"""Calculate the potential temperature.
219 Arguments
220 ---------
221 temperature: iris.cube.Cube | iris.cube.CubeList
222 Cubes of temperature.
223 pressure: iris.cube.Cube | iris.cube.CuebList
224 Cubes of pressure.
226 Returns
227 -------
228 iris.cube.Cube | iris.cube.CubeList
229 Calculated potential temperature in Kelvin.
231 Notes
232 -----
233 The potential temperature is the temperature an air parcel would reach
234 if it was moved adiabatically to a reference pressure. Here we use 1000 hPa
235 as the reference pressure. It is calculated from
237 .. math:: \theta = \frac{T}{\Pi}
239 for :math:`\theta` the potential temperature, T the temperature, and :math:`\Pi`
240 the exner pressure. The exner pressure is calculated using `pressure.exner_pressure`.
242 Temperature must be in Kelvin.
244 All cubes must be on the same grid.
246 Examples
247 --------
248 >>> Theta = temperature.potential_temperature(T, P)
249 """
250 theta = iris.cube.CubeList([])
251 for T, P in zip(iter_maybe(temperature), iter_maybe(pressure), strict=True):
252 TH = T / exner_pressure(P)
253 TH.rename("potential_temperature")
254 theta.append(TH)
255 if len(theta) == 1:
256 return theta[0]
257 else:
258 return theta
261def virtual_potential_temperature(
262 temperature: iris.cube.Cube | iris.cube.CubeList,
263 mixing_ratio: iris.cube.Cube | iris.cube.CubeList,
264 pressure: iris.cube.Cube | iris.cube.CubeList,
265) -> iris.cube.Cube | iris.cube.CubeList:
266 r"""Calculate the virtual potential temperature.
268 Arguments
269 ---------
270 temperature: iris.cube.Cube | iris.cube.CubeList
271 Cubes of temperature.
272 mixing_ratio: iris.cube.Cube | iris.cube.CubeList
273 Cubes of mixing ratio.
274 pressure: iris.cube.Cube | iris.cube.CubeList
275 Cubes of pressure.
277 Returns
278 -------
279 iris.cube.Cube | iris.cube.CubeList
280 Calculated virtual potential temperature in Kelvin.
282 Notes
283 -----
284 The virtual potential temperature is mechanistically equivalent to the
285 potential temperature, except rather than using the (dry-bulb) temperature
286 the virtual temperature used. The virtual potential temperature is the potential
287 temperature a parcel would have if its temperature was replaced by its virtual temperature
288 and then the parcel is brought adiabatically to 1000 hPa. It is calculated as
290 .. math:: \theta_v = \frac{T_v}{\Pi}
292 for :math:`\theta_v` the virtual potential temperature, :math:`T_v` the
293 virtual temperature, and :math:`\Pi` the exner pressure. The exner pressure
294 is calculated using `pressure.exner_pressure`.
296 All cubes must be on the same grid.
298 Examples
299 --------
300 >>> Theta_v = temperature.virtual_potential_temperature(T, W, P)
301 """
302 theta_v = iris.cube.CubeList([])
303 for T, W, P in zip(
304 iter_maybe(temperature),
305 iter_maybe(mixing_ratio),
306 iter_maybe(pressure),
307 strict=True,
308 ):
309 TH_V = virtual_temperature(T, W) / exner_pressure(P)
310 TH_V.rename("virtual_potential_temperature")
311 theta_v.append(TH_V)
312 if len(theta_v) == 1:
313 return theta_v[0]
314 else:
315 return theta_v
318def equivalent_potential_temperature(
319 temperature: iris.cube.Cube | iris.cube.CubeList,
320 relative_humidity: iris.cube.Cube | iris.cube.CubeList,
321 pressure: iris.cube.Cube | iris.cube.CubeList,
322) -> iris.cube.Cube | iris.cube.CubeList:
323 r"""Calculate the equivalent potential temperature.
325 Arguments
326 ---------
327 temperature: iris.cube.Cube | iris.cube.CubeList
328 Cubes of temperature.
329 relative_humidity: iris.cube.Cube | iris.cube.CubeList
330 Cubes of relative humidity.
331 pressure: iris.cube.Cube | iris.cube.CubeList
332 Cubes of pressure.
334 Returns
335 -------
336 iris.cube.Cube | iris.cube.CubeList
337 Calculated equivalent potential temperature in Kelvin.
339 Notes
340 -----
341 The equivalent potential temperature is the temperature an air parcel
342 would have if it was raised until all of the water vapour in it was
343 condensed out and then brought adiabatically to the reference pressure,
344 here taken as 1000 hPa. It is calculated as
346 .. math:: \theta_e = \theta * RH^{- \frac{w R_v}{c_p}} * exp\left(\frac{L_v w}{c_p T} \right)
348 for :math:`\theta_e` the equivalent potential temperature, :math:`\theta` the
349 potential temperature, RH the relative humidity, w the mixing ratio,
350 :math:`R_v` the specific gas constant of water vapour (461
351 :math:`J kg^{-1} K^{-1}`), :math:`c_p` the specific heat capacity of dry air
352 (1005.7 :math:`J kg^{-1} K^{-1}`), :math:`L_v` the latent heat of vapourization
353 (2.5 x :math:`10^6 J kg^{-1}`), and T the temperature.
355 Potential temperature and temperature in K.
356 Relative humidity in percentage and will be converted to fraction.
357 Mixing ratio in kg kg-1 (dimensionless).
359 In this operator the mixing ratio is calculated from
360 `humidity.mixing_ratio_from_relative_humidity` and the potential temperature
361 is calculated from `temperature.potential_temperature`.
363 This equation is a simplification of [Paluch79]_ following [Emanuel94]_,
364 whilst still holding reasonable accuracy.
366 All cubes must be on the same grid.
368 Examples
369 --------
370 >>> Theta_e = temperature.equivalent_potential_temperature(T, RH, P)
371 """
372 theta_e = iris.cube.CubeList([])
373 for T, RH, P in zip(
374 iter_maybe(temperature),
375 iter_maybe(relative_humidity),
376 iter_maybe(pressure),
377 strict=True,
378 ):
379 RH = convert_units(RH, "1")
380 theta = potential_temperature(T, P)
381 w = mixing_ratio_from_relative_humidity(T, P, RH)
382 second_term_power = -(w * RV) / CPD
383 second_term = RH.core_data() ** second_term_power.core_data()
384 third_term_power = LV * w / (CPD * T)
385 third_term = np.exp(third_term_power.core_data())
386 TH_E = theta * second_term * third_term
387 TH_E.rename("equivalent_potential_temperature")
388 TH_E.units = "K"
389 theta_e.append(TH_E)
390 if len(theta_e) == 1:
391 return theta_e[0]
392 else:
393 return theta_e
396def wet_bulb_potential_temperature(
397 temperature: iris.cube.Cube | iris.cube.CubeList,
398 relative_humidity: iris.cube.Cube | iris.cube.CubeList,
399 pressure: iris.cube.Cube | iris.cube.CubeList,
400) -> iris.cube.Cube | iris.cube.CubeList:
401 r"""Calculate wet-bulb potential temperature.
403 Arguments
404 ---------
405 temperature: iris.cube.Cube | iris.cube.CubeList
406 Cubes of temperature.
407 relative_humidity: iris.cube.Cube | iris.cube.CubeList
408 Cubes of relative humidity.
409 pressure: iris.cube.Cube | iris.cube.CubeList
410 Cubes of pressure.
412 Returns
413 -------
414 iris.cube.Cube | iris.cube.CubeList
415 Calculated wet-bulb potential temperature in Kelvin.
417 Notes
418 -----
419 The wet-bulb potential temperature represents the temperature an air parcel
420 would have if cooled adiabatically until saturation and then taken down to
421 a reference pressure (1000 hPa) whilst conserving the moisture (i.e. along a
422 psuedoadiabat). It can be calculated either through a series of iterations,
423 or through empirical relations. Here, we use the [DaviesJones08]_ formulation:
425 .. math:: \theta_w = \theta_e - exp\left(\frac{a_0 + a_1 X + a_2 X^2 + a_3 X^3 + a_4 X^4}{1.0 + b_1 X + b_2 X^2 + b_3 X^3 + b_4 X^4} \right)
427 for :math:`\theta_w` the wet-bulb potential temperature, :math:`\theta_e` the
428 equivalent potential temperature, X = :math:`\theta_e` / 273.15 K, and a
429 series of constants :math:`a_0` = 7.101574, :math:`a_1` = -20.68208,
430 :math:`a_2` = 16.11182, :math:`a_3` = 2.574631, :math:`a_4` = -5.205688,
431 :math:`b_1` = -3.552497, :math:`b_2` = 3.781782, :math:`b_3` = -0.6899655, and
432 :math:`b_4` = -0.5929340.
434 The wet-bulb potential temperature calculated for this method is valid for
435 the range -30 C < :math:`\theta_w` < 50 C and is thus set to NaN outside of
436 this range.
438 All cubes must be on the same grid.
440 Examples
441 --------
442 >>> Theta_w = temperature.wet_bulb_potential_temperature(T, RH, P)
443 """
444 theta_w = iris.cube.CubeList([])
445 for T, RH, P in zip(
446 iter_maybe(temperature),
447 iter_maybe(relative_humidity),
448 iter_maybe(pressure),
449 strict=True,
450 ):
451 TH_E = equivalent_potential_temperature(T, RH, P)
452 X = TH_E / T0
453 X.units = "1"
454 A0 = 7.101574
455 A1 = -20.68208
456 A2 = 16.11182
457 A3 = 2.574631
458 A4 = -5.205688
459 B1 = -3.552497
460 B2 = 3.781782
461 B3 = -0.6899655
462 B4 = -0.5929340
463 exponent = (A0 + A1 * X + A2 * X**2 + A3 * X**3 + A4 * X**4) / (
464 1.0 + B1 * X + B2 * X**2 + B3 * X**3 + B4 * X**4
465 )
466 TH_W = TH_E.copy()
467 TH_W.data[:] = TH_E.core_data() - np.exp(exponent.core_data())
468 TH_W.rename("wet_bulb_potential_temperature")
469 TH_W.data[TH_W.data - T0 < -30.0] = np.nan
470 TH_W.data[TH_W.data - T0 > 50.0] = np.nan
471 theta_w.append(TH_W)
472 if len(theta_w) == 1:
473 return theta_w[0]
474 else:
475 return theta_w
478def saturation_equivalent_potential_temperature(
479 temperature: iris.cube.Cube | iris.cube.CubeList,
480 pressure: iris.cube.Cube | iris.cube.CubeList,
481) -> iris.cube.Cube | iris.cube.CubeList:
482 r"""Calculate the saturation equivalent potential temperature.
484 Arguments
485 ---------
486 temperature: iris.cube.Cube | iris.cube.CubeList
487 Cubes of temperature.
488 pressure: iris.cube.Cube | iris.cube.CubeList
489 Cubes of pressure.
491 Returns
492 -------
493 iris.cube.Cube | iris.cube.CubeList
494 Calculated saturation equivalent potential temperature in Kelvin.
496 Notes
497 -----
498 The saturation equivalent potential temperature, also referred to as the
499 saturation potential temperature is as the equivalent potential temperature
500 following a saturated process throughout. It is calculated as
502 .. math:: \theta_{es} = \theta * exp\left(\frac{L_v w}{c_p T} \right)
504 for :math:`\theta_{es}` the saturation equivalent potential temperature,
505 :math:`\theta` the potential temperature, w the mixing ratio,
506 :math:`c_p` the specific heat capacity of dry air (1005.7
507 :math:`J kg^{-1} K^{-1}`), :math:`L_v` the latent heat of vapourization
508 (2.5 x :math:`10^6 J kg^{-1}`), and T the temperature.
510 As a saturated process is assumed throughout the RH multiplier in the
511 equivalent potential temperature will always be a value of one, and is thus
512 omitted. In this operator the potential temperature is calculated from
513 `temperature.potential_temperature`.
515 All cubes must be on the same grid.
517 Examples
518 --------
519 >>> theta_es = temperature.saturation_equivalent_potential_temperature(T, P)
520 """
521 theta_es = iris.cube.CubeList([])
522 for T, P in zip(
523 iter_maybe(temperature),
524 iter_maybe(pressure),
525 strict=True,
526 ):
527 theta = potential_temperature(T, P)
528 ws = saturation_mixing_ratio(T, P)
529 second_term_power = LV * ws / (CPD * T)
530 second_term = np.exp(second_term_power.core_data())
531 TH_ES = theta * second_term
532 TH_ES.rename("saturation_equivalent_potential_temperature")
533 TH_ES.units = "K"
534 theta_es.append(TH_ES)
535 if len(theta_es) == 1:
536 return theta_es[0]
537 else:
538 return theta_es