Coverage for src / CSET / operators / temperature.py: 30%
117 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-08 16:49 +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 """Calculate the dewpoint temperature following Bolton (1980)."""
38 Td = iris.cube.CubeList([])
39 for T, RH in zip(
40 iter_maybe(temperature), iter_maybe(relative_humidity), strict=True
41 ):
42 vp = vapour_pressure_from_relative_humidity(T, RH)
43 td = vp.copy()
44 td.data = ((243.5 * np.log(vp.core_data())) - 440.8) / (
45 19.48 - np.log(vp.core_data())
46 )
47 td.data[T.data - T0 < -35.0] = np.nan
48 td.data[T.data - T0 > 35.0] = np.nan
49 td.units = "Celsius"
50 td.rename("dewpoint_temperature")
51 td = convert_units(td, "K")
52 Td.append(td)
53 if len(Td) == 1:
54 return Td[0]
55 else:
56 return Td
59def virtual_temperature(
60 temperature: iris.cube.Cube | iris.cube.CubeList,
61 mixing_ratio: iris.cube.Cube | iris.cube.CubeList,
62) -> iris.cube.Cube | iris.cube.CubeList:
63 """Calculate the virtual temperature."""
64 Tv = iris.cube.CubeList([])
65 for T, W in zip(iter_maybe(temperature), iter_maybe(mixing_ratio), strict=True):
66 virT = T * ((W + EPSILON) / (EPSILON * (1 + W)))
67 virT.rename("virtual_temperature")
68 Tv.append(virT)
69 if len(Tv) == 1:
70 return Tv[0]
71 else:
72 return Tv
75def wet_bulb_temperature(
76 temperature: iris.cube.Cube | iris.cube.CubeList,
77 relative_humidity: iris.cube.Cube | iris.cube.CubeList,
78) -> iris.cube.Cube | iris.cube.CubeList:
79 """Calculate the wet-bulb temperature following Stull (2011)."""
80 Tw = iris.cube.CubeList([])
81 for T, RH in zip(
82 iter_maybe(temperature), iter_maybe(relative_humidity), strict=True
83 ):
84 RH = convert_units(RH, "%")
85 T = convert_units(T, "Celsius")
86 wetT = (
87 T * np.arctan(0.151977 * (RH.core_data() + 8.313659) ** 0.5)
88 + np.arctan(T.core_data() + RH.core_data())
89 - np.arctan(RH.core_data() - 1.676331)
90 + 0.00391838
91 * (RH.core_data()) ** (3.0 / 2.0)
92 * np.arctan(0.023101 * RH.core_data())
93 - 4.686035
94 )
95 wetT.rename("wet_bulb_temperature")
96 wetT = convert_units(wetT, "K")
97 Tw.append(wetT)
98 if len(Tw) == 1:
99 return Tw[0]
100 else:
101 return Tw
104def potential_temperature(
105 temperature: iris.cube.Cube | iris.cube.CubeList,
106 pressure: iris.cube.Cube | iris.cube.CubeList,
107) -> iris.cube.Cube | iris.cube.CubeList:
108 """Calculate the potenital temperature."""
109 theta = iris.cube.CubeList([])
110 for T, P in zip(iter_maybe(temperature), iter_maybe(pressure), strict=True):
111 TH = T / exner_pressure(P)
112 TH.rename("potential_temperature")
113 theta.append(TH)
114 if len(theta) == 1:
115 return theta[0]
116 else:
117 return theta
120def virtual_potential_temperature(
121 temperature: iris.cube.Cube | iris.cube.CubeList,
122 mixing_ratio: iris.cube.Cube | iris.cube.CubeList,
123 pressure: iris.cube.Cube | iris.cube.CubeList,
124) -> iris.cube.Cube | iris.cube.CubeList:
125 """Calculate the virtual potential temperature."""
126 theta_v = iris.cube.CubeList([])
127 for T, W, P in zip(
128 iter_maybe(temperature),
129 iter_maybe(mixing_ratio),
130 iter_maybe(pressure),
131 strict=True,
132 ):
133 TH_V = virtual_temperature(T, W) / exner_pressure(P)
134 TH_V.rename("virtual_potential_temperature")
135 theta_v.append(TH_V)
136 if len(theta_v) == 1:
137 return theta_v[0]
138 else:
139 return theta_v
142def equivalent_potential_temperature(
143 temperature: iris.cube.Cube | iris.cube.CubeList,
144 relative_humidity: iris.cube.Cube | iris.cube.CubeList,
145 pressure: iris.cube.Cube | iris.cube.CubeList,
146) -> iris.cube.Cube | iris.cube.CubeList:
147 """Calculate the equivalent potenital temperature."""
148 theta_e = iris.cube.CubeList([])
149 for T, RH, P in zip(
150 iter_maybe(temperature),
151 iter_maybe(relative_humidity),
152 iter_maybe(pressure),
153 strict=True,
154 ):
155 RH = convert_units(RH, "1")
156 theta = potential_temperature(T, P)
157 w = mixing_ratio_from_relative_humidity(T, P, RH)
158 second_term_power = -(w * RV) / CPD
159 second_term = RH.core_data() ** second_term_power.core_data()
160 third_term_power = LV * w / (CPD * T)
161 third_term = np.exp(third_term_power.core_data())
162 TH_E = theta * second_term * third_term
163 TH_E.rename("equivalent_potential_temperature")
164 TH_E.units = "K"
165 theta_e.append(TH_E)
166 if len(theta_e) == 1:
167 return theta_e[0]
168 else:
169 return theta_e
172def wet_bulb_potential_temperature(
173 temperature: iris.cube.Cube | iris.cube.CubeList,
174 relative_humidity: iris.cube.Cube | iris.cube.CubeList,
175 pressure: iris.cube.Cube | iris.cube.CubeList,
176) -> iris.cube.Cube | iris.cube.CubeList:
177 """Calculate wet-bulb potential temperature after Davies-Jones (2008)."""
178 theta_w = iris.cube.CubeList([])
179 for T, RH, P in zip(
180 iter_maybe(temperature),
181 iter_maybe(relative_humidity),
182 iter_maybe(pressure),
183 strict=True,
184 ):
185 TH_E = equivalent_potential_temperature(T, P, RH)
186 X = TH_E / T0
187 X.units = "1"
188 A0 = 7.101574
189 A1 = -20.68208
190 A2 = 16.11182
191 A3 = 2.574631
192 A4 = -5.205688
193 B1 = -3.552497
194 B2 = 3.781782
195 B3 = -0.6899655
196 B4 = -0.5929340
197 exponent = (A0 + A1 * X + A2 * X**2 + A3 * X**3 + A4 * X**4) / (
198 1.0 + B1 * X + B2 * X**2 + B3 * X**3 + B4 * X**4
199 )
200 TH_W = TH_E.copy()
201 TH_W.data[:] = TH_E.core_data() - np.exp(exponent.core_data())
202 TH_W.rename("wet_bulb_potential_temperature")
203 TH_W.data[TH_W.data - T0 < -30.0] = np.nan
204 TH_W.data[TH_W.data - T0 > 50.0] = np.nan
205 theta_w.append(TH_W)
206 if len(theta_w) == 1:
207 return theta_w[0]
208 else:
209 return theta_w
212def saturation_equivalent_potential_temperature(
213 temperature: iris.cube.Cube | iris.cube.CubeList,
214 pressure: iris.cube.Cube | iris.cube.CubeList,
215) -> iris.cube.Cube | iris.cube.CubeList:
216 """Calculate the saturation equivalent potenital temperature."""
217 theta_s = iris.cube.CubeList([])
218 for T, P in zip(
219 iter_maybe(temperature),
220 iter_maybe(pressure),
221 strict=True,
222 ):
223 theta = potential_temperature(T, P)
224 ws = saturation_mixing_ratio(T, P)
225 second_term_power = LV * ws / (CPD * T)
226 second_term = np.exp(second_term_power.core_data())
227 TH_S = theta * second_term
228 TH_S.rename("saturation_equivalent_potential_temperature")
229 TH_S.units = "K"
230 theta_s.append(TH_S)
231 if len(theta_s) == 1:
232 return theta_s[0]
233 else:
234 return theta_s