Coverage for src / CSET / operators / temperature.py: 100%
121 statements
« prev ^ index » next coverage.py v7.13.2, created at 2026-01-28 11:16 +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 References
72 ----------
73 .. [Bolton80] Bolton. D. (1980) "The computation of equivalent potential
74 temperature". Monthly Weather Review, vol. 108, 1046-1053,
75 doi: 10.1175/1520-0493(1980)108<1046:TCOEPT>2.0.CO;2
76 .. [Flack24] Flack, D.L.A. (2024) "Stratification of the vertical spread-skill
77 relation by radiosonde drift in a convective-scale ensemble."
78 Atmospheric Science Letters, vol. 25, e1194, doi: 10.1002/asl.1194
80 Examples
81 --------
82 >>> Td = temperature.dewpoint_temperature(T, RH)
83 """
84 Td = iris.cube.CubeList([])
85 for T, RH in zip(
86 iter_maybe(temperature), iter_maybe(relative_humidity), strict=True
87 ):
88 vp = vapour_pressure_from_relative_humidity(T, RH)
89 td = vp.copy()
90 td.data = ((243.5 * np.log(vp.core_data())) - 440.8) / (
91 19.48 - np.log(vp.core_data())
92 )
93 td.data[T.data - T0 < -35.0] = np.nan
94 td.data[T.data - T0 > 35.0] = np.nan
95 td.units = "Celsius"
96 td.rename("dewpoint_temperature")
97 td = convert_units(td, "K")
98 Td.append(td)
99 if len(Td) == 1:
100 return Td[0]
101 else:
102 return Td
105def virtual_temperature(
106 temperature: iris.cube.Cube | iris.cube.CubeList,
107 mixing_ratio: iris.cube.Cube | iris.cube.CubeList,
108) -> iris.cube.Cube | iris.cube.CubeList:
109 r"""Calculate the virtual temperature.
111 Arguments
112 ---------
113 temperature: iris.cube.Cube | iris.cube.CubeList
114 Cubes of temperature in Kelvin.
115 mixing_ratio: iris.cube.Cube | iris.cube.CubeList
116 Cubes of mixing ratio.
118 Returns
119 -------
120 iris.cube.Cube | iris.cube.CubeList
121 Calculated virtual temperature in Kelvin.
123 Notes
124 -----
125 The virtual temperature is that required for a dry parcel to have the
126 same density as a moist parcel at the same pressure. It is calculated
127 as
129 .. math:: T_v = T * \frac{w + \epsilon}{\epsilon (1 + w)}
131 for :math:`T_v` the virtual temperature, T the temperature, w the mixing
132 ratio, and :math:`\epsilon` the ratio between dry and moist air equating
133 to 0.622.
135 T is given in K and w is given in kg kg-1.
136 All cubes must be on the same grid.
138 Examples
139 --------
140 >>> Tv = temperature.virtual_temperature(T, w)
141 """
142 Tv = iris.cube.CubeList([])
143 for T, W in zip(iter_maybe(temperature), iter_maybe(mixing_ratio), strict=True):
144 virT = T * ((W + EPSILON) / (EPSILON * (1 + W)))
145 virT.rename("virtual_temperature")
146 Tv.append(virT)
147 if len(Tv) == 1:
148 return Tv[0]
149 else:
150 return Tv
153def wet_bulb_temperature(
154 temperature: iris.cube.Cube | iris.cube.CubeList,
155 relative_humidity: iris.cube.Cube | iris.cube.CubeList,
156) -> iris.cube.Cube | iris.cube.CubeList:
157 r"""Calculate the wet-bulb temperature.
159 Arguments
160 ---------
161 temperature: iris.cube.Cube | iris.cube.CubeList
162 Cubes of temperature.
163 relative_humidity: iris.cube.Cube | iris.cube.CubeList
164 Cubes of relative humidity.
166 Returns
167 -------
168 iris.cube.Cube | iris.cube.CubeList
169 Calculated wet-bulb temperature in Kelvin.
171 Notes
172 -----
173 The wet-bulb temperature is the temperature the air reaches when all the
174 water has been evaporated out of it. It can be calculated from temperature in Celsius
175 and relative humidity in percent following [Stull11]_
177 .. 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
179 for :math:`T_w` the wet-bulb temperature, T the temperature, and RH the relative
180 humidity.
182 The equation is valid for 5% < RH < 99% and -20 C < T < 50 C, and results are
183 only presented for these values.
185 The temperature and relative humidity unit conversions are applied in the
186 operator.
188 All cubes should be on the same grid.
190 References
191 ----------
192 .. [Stull11] Stull, R. (2011) "Wet-Bulb Temperature from Relative Humidity
193 and air temperature." Journal of Applied Meteorology and Climatology, vol. 50,
194 2267-2269, doi: 10.1175/JAMC-D-11-0143.1
196 Examples
197 --------
198 >>> Tw = temperature.wet_bulb_temperature(T, RH)
199 """
200 Tw = iris.cube.CubeList([])
201 for T, RH in zip(
202 iter_maybe(temperature), iter_maybe(relative_humidity), strict=True
203 ):
204 RH = convert_units(RH, "%")
205 T = convert_units(T, "Celsius")
206 wetT = (
207 T * np.arctan(0.151977 * (RH.core_data() + 8.313659) ** 0.5)
208 + np.arctan(T.core_data() + RH.core_data())
209 - np.arctan(RH.core_data() - 1.676331)
210 + 0.00391838
211 * (RH.core_data()) ** (3.0 / 2.0)
212 * np.arctan(0.023101 * RH.core_data())
213 - 4.686035
214 )
215 wetT.rename("wet_bulb_temperature")
216 wetT = convert_units(wetT, "K")
217 wetT.data[T.core_data() < -20.0] = np.nan
218 wetT.data[T.core_data() > 50.0] = np.nan
219 wetT.data[RH.core_data() < 5.0] = np.nan
220 wetT.data[RH.core_data() > 99.0] = np.nan
221 Tw.append(wetT)
222 if len(Tw) == 1:
223 return Tw[0]
224 else:
225 return Tw
228def potential_temperature(
229 temperature: iris.cube.Cube | iris.cube.CubeList,
230 pressure: iris.cube.Cube | iris.cube.CubeList,
231) -> iris.cube.Cube | iris.cube.CubeList:
232 r"""Calculate the potential temperature.
234 Arguments
235 ---------
236 temperature: iris.cube.Cube | iris.cube.CubeList
237 Cubes of temperature.
238 pressure: iris.cube.Cube | iris.cube.CuebList
239 Cubes of pressure.
241 Returns
242 -------
243 iris.cube.Cube | iris.cube.CubeList
244 Calculated potential temperature in Kelvin.
246 Notes
247 -----
248 The potential temperature is the temperature an air parcel would reach
249 if it was moved adiabatically to a reference pressure. Here we use 1000 hPa
250 as the reference pressure. It is calculated from
252 .. math:: \theta = \frac{T}{\Pi}
254 for :math:`\theta` the potential temperature, T the temperature, and :math:`\Pi`
255 the exner pressure. The exner pressure is calculated using `pressure.exner_pressure`.
257 Temperature must be in Kelvin.
259 All cubes must be on the same grid.
261 Examples
262 --------
263 >>> Theta = temperature.potential_temperature(T, P)
264 """
265 theta = iris.cube.CubeList([])
266 for T, P in zip(iter_maybe(temperature), iter_maybe(pressure), strict=True):
267 TH = T / exner_pressure(P)
268 TH.rename("potential_temperature")
269 theta.append(TH)
270 if len(theta) == 1:
271 return theta[0]
272 else:
273 return theta
276def virtual_potential_temperature(
277 temperature: iris.cube.Cube | iris.cube.CubeList,
278 mixing_ratio: iris.cube.Cube | iris.cube.CubeList,
279 pressure: iris.cube.Cube | iris.cube.CubeList,
280) -> iris.cube.Cube | iris.cube.CubeList:
281 r"""Calculate the virtual potential temperature.
283 Arguments
284 ---------
285 temperature: iris.cube.Cube | iris.cube.CubeList
286 Cubes of temperature.
287 mixing_ratio: iris.cube.Cube | iris.cube.CubeList
288 Cubes of mixing ratio.
289 pressure: iris.cube.Cube | iris.cube.CubeList
290 Cubes of pressure.
292 Returns
293 -------
294 iris.cube.Cube | iris.cube.CubeList
295 Calculated virtual potential temperature in Kelvin.
297 Notes
298 -----
299 The virtual potential temperature is mechanistically equivalent to the
300 potential temperature, except rather than using the (dry-bulb) temperature
301 the virtual temperature used. The virtual potential temperature is the potential
302 temperature a parcel would have if its temperature was replaced by its virtual temperature
303 and then the parcel is brought adiabatically to 1000 hPa. It is calculated as
305 .. math:: \theta_v = \frac{T_v}{\Pi}
307 for :math:`\theta_v` the virtual potential temperature, :math:`T_v` the
308 virtual temperature, and :math:`\Pi` the exner pressure. The exner pressure
309 is calculated using `pressure.exner_pressure`.
311 All cubes must be on the same grid.
313 Examples
314 --------
315 >>> Theta_v = temperature.virtual_potential_temperature(T, W, P)
316 """
317 theta_v = iris.cube.CubeList([])
318 for T, W, P in zip(
319 iter_maybe(temperature),
320 iter_maybe(mixing_ratio),
321 iter_maybe(pressure),
322 strict=True,
323 ):
324 TH_V = virtual_temperature(T, W) / exner_pressure(P)
325 TH_V.rename("virtual_potential_temperature")
326 theta_v.append(TH_V)
327 if len(theta_v) == 1:
328 return theta_v[0]
329 else:
330 return theta_v
333def equivalent_potential_temperature(
334 temperature: iris.cube.Cube | iris.cube.CubeList,
335 relative_humidity: iris.cube.Cube | iris.cube.CubeList,
336 pressure: iris.cube.Cube | iris.cube.CubeList,
337) -> iris.cube.Cube | iris.cube.CubeList:
338 r"""Calculate the equivalent potential temperature.
340 Arguments
341 ---------
342 temperature: iris.cube.Cube | iris.cube.CubeList
343 Cubes of temperature.
344 relative_humidity: iris.cube.Cube | iris.cube.CubeList
345 Cubes of relative humidity.
346 pressure: iris.cube.Cube | iris.cube.CubeList
347 Cubes of pressure.
349 Returns
350 -------
351 iris.cube.Cube | iris.cube.CubeList
352 Calculated equivalent potential temperature in Kelvin.
354 Notes
355 -----
356 The equivalent potential temperature is the temperature an air parcel
357 would have if it was raised until all of the water vapour in it was
358 condensed out and then brought adiabatically to the reference pressure,
359 here taken as 1000 hPa. It is calculated as
361 .. math:: \theta_e = \theta * RH^{- \frac{w R_v}{c_p}} * exp\left(\frac{L_v w}{c_p T} \right)
363 for :math:`\theta_e` the equivalent potential temperature, :math:`\theta` the
364 potential temperature, RH the relative humidity, w the mixing ratio,
365 :math:`R_v` the specific gas constant of water vapour (461
366 :math:`J kg^{-1} K^{-1}`), :math:`c_p` the specific heat capacity of dry air
367 (1005.7 :math:`J kg^{-1} K^{-1}`), :math:`L_v` the latent heat of vapourization
368 (2.5 x :math:`10^6 J kg^{-1}`), and T the temperature.
370 Potential temperature and temperature in K.
371 Relative humidity in percentage and will be converted to fraction.
372 Mixing ratio in kg kg-1 (dimensionless).
374 In this operator the mixing ratio is calculated from
375 `humidity.mixing_ratio_from_relative_humidity` and the potential temperature
376 is calculated from `temperature.potential_temperature`.
378 This equation is a simplification of [Paluch79]_ following [Emanuel94]_,
379 whilst still holding reasonable accuracy.
381 All cubes must be on the same grid.
383 References
384 ----------
385 .. [Emanuel94] Emanuel, K.A. (1994) "Atmospheric Convection" Oxford University
386 Press, 580 pp.
387 .. [Paluch79] Paluch, I.R., (1979) "The entrainment mechanism in Colorado
388 Cumuli" Journal of the Atmospheric Sciences, vol. 36, 2467-2478,
389 doi: 10.1175/1520-0469(1979)036<2467:TEMICC>2.0.CO;2
391 Examples
392 --------
393 >>> Theta_e = temperature.equivalent_potential_temperature(T, RH, P)
394 """
395 theta_e = iris.cube.CubeList([])
396 for T, RH, P in zip(
397 iter_maybe(temperature),
398 iter_maybe(relative_humidity),
399 iter_maybe(pressure),
400 strict=True,
401 ):
402 RH = convert_units(RH, "1")
403 theta = potential_temperature(T, P)
404 w = mixing_ratio_from_relative_humidity(T, P, RH)
405 second_term_power = -(w * RV) / CPD
406 second_term = RH.core_data() ** second_term_power.core_data()
407 third_term_power = LV * w / (CPD * T)
408 third_term = np.exp(third_term_power.core_data())
409 TH_E = theta * second_term * third_term
410 TH_E.rename("equivalent_potential_temperature")
411 TH_E.units = "K"
412 theta_e.append(TH_E)
413 if len(theta_e) == 1:
414 return theta_e[0]
415 else:
416 return theta_e
419def wet_bulb_potential_temperature(
420 temperature: iris.cube.Cube | iris.cube.CubeList,
421 relative_humidity: iris.cube.Cube | iris.cube.CubeList,
422 pressure: iris.cube.Cube | iris.cube.CubeList,
423) -> iris.cube.Cube | iris.cube.CubeList:
424 r"""Calculate wet-bulb potential temperature.
426 Arguments
427 ---------
428 temperature: iris.cube.Cube | iris.cube.CubeList
429 Cubes of temperature.
430 relative_humidity: iris.cube.Cube | iris.cube.CubeList
431 Cubes of relative humidity.
432 pressure: iris.cube.Cube | iris.cube.CubeList
433 Cubes of pressure.
435 Returns
436 -------
437 iris.cube.Cube | iris.cube.CubeList
438 Calculated wet-bulb potential temperature in Kelvin.
440 Notes
441 -----
442 The wet-bulb potential temperature represents the temperature an air parcel
443 would have if cooled adiabatically until saturation and then taken down to
444 a reference pressure (1000 hPa) whilst conserving the moisture (i.e. along a
445 psuedoadiabat). It can be calculated either through a series of iterations,
446 or through empirical relations. Here, we use the [DaviesJones08]_ formulation:
448 .. 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)
450 for :math:`\theta_w` the wet-bulb potential temperature, :math:`\theta_e` the
451 equivalent potential temperature, X = :math:`\theta_e` / 273.15 K, and a
452 series of constants :math:`a_0` = 7.101574, :math:`a_1` = -20.68208,
453 :math:`a_2` = 16.11182, :math:`a_3` = 2.574631, :math:`a_4` = -5.205688,
454 :math:`b_1` = -3.552497, :math:`b_2` = 3.781782, :math:`b_3` = -0.6899655, and
455 :math:`b_4` = -0.5929340.
457 The wet-bulb potential temperature calculated for this method is valid for
458 the range -30 C < :math:`\theta_w` < 50 C and is thus set to NaN outside of
459 this range.
461 All cubes must be on the same grid.
463 References
464 ----------
465 .. [DaviesJones08] Davies-Jones, R. (2008) "An Efficient and Accurate
466 Method for Computing the Wet-Bulb Temperature along Pseudoadiabats"
467 Monthly Weather Review, vol. 136, 2764-2785, doi: 10.1175/2007MWR2224.1
469 Examples
470 --------
471 >>> Theta_w = temperature.wet_bulb_potential_temperature(T, RH, P)
472 """
473 theta_w = iris.cube.CubeList([])
474 for T, RH, P in zip(
475 iter_maybe(temperature),
476 iter_maybe(relative_humidity),
477 iter_maybe(pressure),
478 strict=True,
479 ):
480 TH_E = equivalent_potential_temperature(T, RH, P)
481 X = TH_E / T0
482 X.units = "1"
483 A0 = 7.101574
484 A1 = -20.68208
485 A2 = 16.11182
486 A3 = 2.574631
487 A4 = -5.205688
488 B1 = -3.552497
489 B2 = 3.781782
490 B3 = -0.6899655
491 B4 = -0.5929340
492 exponent = (A0 + A1 * X + A2 * X**2 + A3 * X**3 + A4 * X**4) / (
493 1.0 + B1 * X + B2 * X**2 + B3 * X**3 + B4 * X**4
494 )
495 TH_W = TH_E.copy()
496 TH_W.data[:] = TH_E.core_data() - np.exp(exponent.core_data())
497 TH_W.rename("wet_bulb_potential_temperature")
498 TH_W.data[TH_W.data - T0 < -30.0] = np.nan
499 TH_W.data[TH_W.data - T0 > 50.0] = np.nan
500 theta_w.append(TH_W)
501 if len(theta_w) == 1:
502 return theta_w[0]
503 else:
504 return theta_w
507def saturation_equivalent_potential_temperature(
508 temperature: iris.cube.Cube | iris.cube.CubeList,
509 pressure: iris.cube.Cube | iris.cube.CubeList,
510) -> iris.cube.Cube | iris.cube.CubeList:
511 r"""Calculate the saturation equivalent potential temperature.
513 Arguments
514 ---------
515 temperature: iris.cube.Cube | iris.cube.CubeList
516 Cubes of temperature.
517 pressure: iris.cube.Cube | iris.cube.CubeList
518 Cubes of pressure.
520 Returns
521 -------
522 iris.cube.Cube | iris.cube.CubeList
523 Calculated saturation equivalent potential temperature in Kelvin.
525 Notes
526 -----
527 The saturation equivalent potential temperature, also referred to as the
528 saturation potential temperature is as the equivalent potential temperature
529 following a saturated process throughout. It is calculated as
531 .. math:: \theta_{es} = \theta * exp\left(\frac{L_v w}{c_p T} \right)
533 for :math:`\theta_{es}` the saturation equivalent potential temperature,
534 :math:`\theta` the potential temperature, w the mixing ratio,
535 :math:`c_p` the specific heat capacity of dry air (1005.7
536 :math:`J kg^{-1} K^{-1}`), :math:`L_v` the latent heat of vapourization
537 (2.5 x :math:`10^6 J kg^{-1}`), and T the temperature.
539 As a saturated process is assumed throughout the RH multiplier in the
540 equivalent potential temperature will always be a value of one, and is thus
541 omitted. In this operator the potential temperature is calculated from
542 `temperature.potential_temperature`.
544 All cubes must be on the same grid.
546 Examples
547 --------
548 >>> theta_es = temperature.saturation_equivalent_potential_temperature(T, P)
549 """
550 theta_es = iris.cube.CubeList([])
551 for T, P in zip(
552 iter_maybe(temperature),
553 iter_maybe(pressure),
554 strict=True,
555 ):
556 theta = potential_temperature(T, P)
557 ws = saturation_mixing_ratio(T, P)
558 second_term_power = LV * ws / (CPD * T)
559 second_term = np.exp(second_term_power.core_data())
560 TH_ES = theta * second_term
561 TH_ES.rename("saturation_equivalent_potential_temperature")
562 TH_ES.units = "K"
563 theta_es.append(TH_ES)
564 if len(theta_es) == 1:
565 return theta_es[0]
566 else:
567 return theta_es