diff --git a/docs/source/conf.py b/docs/source/conf.py index e104f4220..1834f0222 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -8,7 +8,7 @@ # https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information project = "CSET" -copyright = "Crown copyright, Met Office (2022-2025) and CSET contributors." +copyright = "Crown copyright, Met Office (2022-2026) and CSET contributors." author = "Met Office and Partners" # -- General configuration --------------------------------------------------- diff --git a/docs/source/index.rst b/docs/source/index.rst index 1d103a465..53c35565c 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -104,7 +104,7 @@ Conduct`_. Licence ------- -© Crown copyright, Met Office (2022-2025) and CSET contributors. +© Crown copyright, Met Office (2022-2026) and CSET contributors. Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a diff --git a/docs/source/reference/operators.rst b/docs/source/reference/operators.rst index 259bf9815..db6ab6aa1 100644 --- a/docs/source/reference/operators.rst +++ b/docs/source/reference/operators.rst @@ -95,6 +95,15 @@ CSET.operators.ensembles .. automodule:: CSET.operators.ensembles :members: +Humidity Operators +~~~~~~~~~~~~~~~~~~ + +CSET.operators.humidity +----------------------- + +.. automodule:: CSET.operators.humidity + :members: + Image Processing Operators ~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -113,6 +122,24 @@ CSET.operators.mesoscale .. automodule:: CSET.operators.mesoscale :members: +Pressure Operators +~~~~~~~~~~~~~~~~~~ + +CSET.operators.pressure +----------------------- + +.. automodule:: CSET.operators.pressure + :members: + +Temperature Operators +~~~~~~~~~~~~~~~~~~~~~ + +CSET.operators.temperature +-------------------------- + +.. automodule:: CSET.operators.temperature + :members: + Wind Operators ~~~~~~~~~~~~~~~~~~~ diff --git a/src/CSET/operators/__init__.py b/src/CSET/operators/__init__.py index ef40fdf22..cf583c30d 100644 --- a/src/CSET/operators/__init__.py +++ b/src/CSET/operators/__init__.py @@ -1,4 +1,4 @@ -# © Crown copyright, Met Office (2022-2024) and CSET contributors. +# © Crown copyright, Met Office (2022-2026) and CSET contributors. # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. @@ -34,12 +34,15 @@ convection, ensembles, filters, + humidity, imageprocessing, mesoscale, misc, plot, + pressure, read, regrid, + temperature, transect, wind, write, @@ -56,13 +59,16 @@ "ensembles", "execute_recipe", "filters", + "humidity", "get_operator", "imageprocessing", "mesoscale", "misc", "plot", + "pressure", "read", "regrid", + "temperature", "transect", "wind", "write", diff --git a/src/CSET/operators/_atmospheric_constants.py b/src/CSET/operators/_atmospheric_constants.py new file mode 100644 index 000000000..41e55c358 --- /dev/null +++ b/src/CSET/operators/_atmospheric_constants.py @@ -0,0 +1,41 @@ +# © Crown copyright, Met Office (2022-2026) and CSET contributors. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +"""Constants for the atmosphere.""" + +# Reference pressure. +P0 = 1000.0 # hPa. + +# Specific gas constant for dry air. +RD = 287.0 # J/kg/K. + +# Specific gas constant for water vapour. +RV = 461.0 # J/kg/K. + +# Specific heat capacity for dry air. +CPD = 1005.7 # J/kg/K. + +# Latent heat of vaporization. +LV = 2.501e6 # J/kg/K. + +# Reference vapour pressure. +E0 = 6.1078 # hPa. + +# Reference temperature. +T0 = 273.15 # K. + +# Ratio between mixing ratio of dry and moist air. +EPSILON = 0.622 + +# Ratio between specific gas constant and specific heat capacity. +KAPPA = RD / CPD diff --git a/src/CSET/operators/humidity.py b/src/CSET/operators/humidity.py new file mode 100644 index 000000000..85856f14a --- /dev/null +++ b/src/CSET/operators/humidity.py @@ -0,0 +1,445 @@ +# © Crown copyright, Met Office (2022-2026) and CSET contributors. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Operators for humidity conversions.""" + +import iris.cube + +from CSET._common import iter_maybe +from CSET.operators._atmospheric_constants import EPSILON +from CSET.operators.misc import convert_units +from CSET.operators.pressure import vapour_pressure + + +def mixing_ratio_from_specific_humidity( + specific_humidity: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Convert specific humidity to mixing ratio. + + Arguments + --------- + specific_humidity: iris.cube.Cube | iris.cube.CubeList + Cubes of specific humidity to be converted to mixing ratio. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + Converted mixing ratio. + + Notes + ----- + Atmospheric water vapour can be described by multiple quantities. Here, + we convert the specific humidity to the mixing ratio using the following + relation + + .. math:: w = \frac{q}{1 - q} + + with w the mixing ratio and q the specific humidity. + + Larger mixing ratios imply more moisture in the atmosphere. The mixing + ratio will have the same units as the specific humidity (kg/kg). + + + Examples + -------- + >>> w = humidity.mixing_ratio_from_specific_humidity(specific_humidity) + """ + w = iris.cube.CubeList([]) + for q in iter_maybe(specific_humidity): + mr = q.copy() + mr = q / (1 - q) + mr.rename("mixing_ratio") + w.append(mr) + if len(w) == 1: + return w[0] + else: + return w + + +def specific_humidity_from_mixing_ratio( + mixing_ratio: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Convert mixing ratio to specific humidity. + + Arguments + --------- + mixing_ratio: iris.cube.Cube | iris.Cube.CubeList + Cubes of mixing ratio to be converted to specific humidity. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + Converted specific humidity. + + Notes + ----- + Here, we invert the relationship from `humidity.mixing_ratio_from_specific_humidity` + for the following relation + + .. math:: q = \frac{w}{1 + w} + + with q the specific humidity and w the mixing ratio. + + A larger specific humidity implies a more moist atmosphere. The specific + humidity will have the same units as the mixing ratio (kg/kg). + + Examples + -------- + >>> q = humidity.specific_humidity_from_mixing_ratio(mixing_ratio) + """ + q = iris.cube.CubeList([]) + for w in iter_maybe(mixing_ratio): + sh = w.copy() + sh = w / (1 + w) + sh.rename("specific_humidity") + q.append(sh) + if len(q) == 1: + return q[0] + else: + return q + + +def saturation_mixing_ratio( + temperature: iris.cube.Cube | iris.cube.CubeList, + pressure: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Calculate saturation mixing ratio. + + Arguments + --------- + temperature: iris.cube.Cube | iris.cube.CubeList + Cubes of temperature in Kelvin. + pressure: iris.cube.Cube | iris.cube.CubeList + Cubes of pressure. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + Saturation mixing ratio in kg/kg. + + Notes + ----- + The saturation mixing ratio is required to help calculate the relative + humidity and other diagnostics with respect to the mixing ratio. It can + be calculated from + + .. math:: w = \epsilon \frac{e}{P - e} + + for w the mixing ratio, :math:`\epsilon` the ratio between the mixing ratio + of dry and moist air equating to 0.622, P the pressure and e the vapour + pressure. To ensure that the saturation mixing ratio (:math:`w_s`) is + calculated the vapour pressure calculated should be with the (dry-bulb) + temperature to ensure it is the saturated vapour pressure. + + All cubes need to be on the same grid. + + Examples + -------- + >>> w_s = humidity.saturation_mixing_ratio(temperature, pressure) + """ + w = iris.cube.CubeList([]) + for T, P in zip(iter_maybe(temperature), iter_maybe(pressure), strict=True): + P = convert_units(P, "hPa") + mr = (EPSILON * vapour_pressure(T)) / (P - vapour_pressure(T)) + mr.units = "kg/kg" + mr.rename("saturation_mixing_ratio") + w.append(mr) + if len(w) == 1: + return w[0] + else: + return w + + +def saturation_specific_humidity( + temperature: iris.cube.Cube | iris.cube.CubeList, + pressure: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Calculate saturation specific humidity. + + Arguments + --------- + temperature: iris.cube.Cube | iris.cube.CubeList + Cubes of temperature in Kelvin. + pressure: iris.cube.Cube | iris.cube.CubeList + Cubes of pressure. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + Saturation specific humidity in kg/kg. + + Notes + ----- + The saturation specific humidity is required to help calculate the relative + humidity and other diagnostics with respect to the mixing ratio. It can + be calculated from + + .. math:: q = \epsilon \frac{e}{P} + + for q the specific humidity, :math:`\epsilon` the ratio between the mixing ratio + of dry and moist air equating to 0.622, P the pressure and e the vapour + pressure. To ensure that the saturation specific humidity (:math:`q_{sat}`) is + calculated the vapour pressure calculated should be with the (dry-bulb) + temperature to ensure it is the saturated vapour pressure. + + All cubes need to be on the same grid. + + Examples + -------- + >>> q_sat = humidity.saturation_specific_humidity(temperature, pressure) + """ + q = iris.cube.CubeList([]) + for T, P in zip(iter_maybe(temperature), iter_maybe(pressure), strict=True): + P = convert_units(P, "hPa") + sh = (EPSILON * vapour_pressure(T)) / P + sh.units = "kg/kg" + sh.rename("saturation_specific_humidity") + q.append(sh) + if len(q) == 1: + return q[0] + else: + return q + + +def mixing_ratio_from_relative_humidity( + temperature: iris.cube.Cube | iris.cube.CubeList, + pressure: iris.cube.Cube | iris.cube.CubeList, + relative_humidity: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Calculate the mixing ratio from RH. + + Arguments + --------- + temperature: iris.cube.Cube | iris.cube.CubeList + Cubes of temperature in Kelvin. + pressure: iris.cube.Cube | iris.cube.CubeList + Cubes of pressure. + relative_humidity: iris.cube.Cube | iris.cube.CubeList + Cubes of relative humidity. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + Calculated mixing ratio from relative humidity in kg/kg. + + Notes + ----- + The mixing ratio can be calculated from temperature, pressure, and + relative humidity using the following relation + + .. math:: w = RH * w_s + + for w the mixing ratio, :math:`w_s` the saturation mixing ratio, and + RH the relative humidity. RH is converted to dimensionless fraction rather + than percentage. + + The operator uses `humidity.saturation_mixing_ratio` to calculate the + saturation mixing ratio from the temperature and pressure. The relative + humidity is converted into a decimal before the multiplication occurs. + + All cubes need to be on the same grid. + + Examples + -------- + >>> w = humidity.mixing_ratio_from_relative_humidity(T, P, RH) + """ + w = iris.cube.CubeList([]) + for T, P, RH in zip( + iter_maybe(temperature), + iter_maybe(pressure), + iter_maybe(relative_humidity), + strict=True, + ): + RH = convert_units(RH, "1") + mr = saturation_mixing_ratio(T, P) * RH + mr.rename("mixing_ratio") + mr.units = "kg/kg" + w.append(mr) + if len(w) == 1: + return w[0] + else: + return w + + +def specific_humidity_from_relative_humidity( + temperature: iris.cube.Cube | iris.cube.CubeList, + pressure: iris.cube.Cube | iris.cube.CubeList, + relative_humidity: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Calculate the specific humidity from relative humidity. + + Arguments + --------- + temperature: iris.cube.Cube | iris.cube.CubeList + Cubes of temperature in Kelvin. + pressure: iris.cube.Cube | iris.cube.CubeList + Cubes of pressure. + relative_humidity: iris.cube.Cube | iris.cube.CubeList + Cubes of relative humidity. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + Calculated specific humidity from relative humidity in kg/kg. + + Notes + ----- + The specific humidity can be calculated from temperature, pressure, and + relative humidity using the following relation + + .. math:: q = RH * q_{sat} + + for q the specific humidity, :math:`q_{sat}` the saturation specific + humidity, and RH the relative humidity. + + The operator uses `humidity.saturation_specific_humidity` to calculate the + saturation specific humidity from the temperature and pressure. The relative + humidity is converted into a decimal before the multiplication occurs. + + All cubes need to be on the same grid. + + Examples + -------- + >>> q = humidity.specific_humidity_from_relative_humidity(T, P, RH) + """ + q = iris.cube.CubeList([]) + for T, P, RH in zip( + iter_maybe(temperature), + iter_maybe(pressure), + iter_maybe(relative_humidity), + strict=True, + ): + RH = convert_units(RH, "1") + sh = saturation_specific_humidity(T, P) * RH + sh.rename("specific_humidity") + sh.units = "kg/kg" + q.append(sh) + if len(q) == 1: + return q[0] + else: + return q + + +def relative_humidity_from_mixing_ratio( + mixing_ratio: iris.cube.Cube | iris.cube.CubeList, + temperature: iris.cube.Cube | iris.cube.CubeList, + pressure: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Convert mixing ratio to relative humidity. + + Arguments + --------- + mixing_ratio: iris.cube.Cube | iris.cube.CubeList + Cubes of mixing ratio. + temperature: iris.cube.Cube | iris.cube.CuebList + Cubes of temperature in Kelvin. + pressure: iris.cube.Cube | iris.cube.CubeList + Cubes of pressure. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + Relative humidity calculated from mixing ratio. + + Notes + ----- + The relative humidity can be calculated from the mixing ratio following + + .. math:: RH = \frac{w}{w_s} + + for RH the relative humidity, w the mixing ratio, and :math:`w_s` the + saturation mixing ratio. The saturation mixing ratio is calculated using + `humidity.saturation_mixing_ratio`. + + The RH varies predominatly between zero (completely dry) and one (saturated). + Values larger than one are possible and imply supersaturation. + + All cubes must be on the same grid. + + Examples + -------- + >>> RH = humidity.relative_humidity_from_mixing_ratio(w, T, P) + """ + RH = iris.cube.CubeList([]) + for W, T, P in zip( + iter_maybe(mixing_ratio), + iter_maybe(temperature), + iter_maybe(pressure), + strict=True, + ): + rel_h = W / saturation_mixing_ratio(T, P) + rel_h.rename("relative_humidity") + rel_h = convert_units(rel_h, "%") + RH.append(rel_h) + if len(RH) == 1: + return RH[0] + else: + return RH + + +def relative_humidity_from_specific_humidity( + specific_humidity: iris.cube.Cube | iris.cube.CubeList, + temperature: iris.cube.Cube | iris.cube.CubeList, + pressure: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Convert specific humidity to relative humidity. + + Arguments + --------- + specific_humidity: iris.cube.Cube | iris.cube.CubeList + Cubes of specific humidity. + temperature: iris.cube.Cube | iris.cube.CuebList + Cubes of temperature in Kelvin. + pressure: iris.cube.Cube | iris.cube.CubeList + Cubes of pressure. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + Relative humidity calculated from specific humidity. + + Notes + ----- + The relative humidity can be calculated from the specific humidity following + + .. math:: RH = \frac{q}{q_{sat}} + + for RH the relative humidity, q the specific humidity, and :math:`q_{sat}` the + saturation specific humidity. The saturation specific humidity is calculated using + `humidity.saturation_specific_humidity`. + + The RH varies predominatly between zero (completely dry) and one (saturated). + Values larger than one are possible and imply supersaturation. + + All cubes must be on the same grid. + + Examples + -------- + >>> RH = humidity.relative_humidity_from_specific_humidity(q, T, P) + """ + RH = iris.cube.CubeList([]) + for Q, T, P in zip( + iter_maybe(specific_humidity), + iter_maybe(temperature), + iter_maybe(pressure), + strict=True, + ): + rel_h = Q / saturation_specific_humidity(T, P) + rel_h.rename("relative_humidity") + rel_h = convert_units(rel_h, "%") + RH.append(rel_h) + if len(RH) == 1: + return RH[0] + else: + return RH diff --git a/src/CSET/operators/pressure.py b/src/CSET/operators/pressure.py new file mode 100644 index 000000000..544ea19c0 --- /dev/null +++ b/src/CSET/operators/pressure.py @@ -0,0 +1,191 @@ +# © Crown copyright, Met Office (2022-2026) and CSET contributors. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Operators for pressure conversions.""" + +import iris.cube +import numpy as np + +from CSET._common import iter_maybe +from CSET.operators._atmospheric_constants import E0, KAPPA, P0 +from CSET.operators.misc import convert_units + + +def vapour_pressure( + temperature: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Calculate the vapour pressure of the atmosphere. + + Arguments + --------- + temperature: iris.cube.Cube | iris.cube.CubeList + Cubes of temperature to be converted into vapour pressure. + Temperature must be provided in Kelvin. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + Vapour pressure in hPa. + + Notes + ----- + The vapour pressure represents the pressure exerted by water vapour on the + atmosphere. It is related to atmospheric temperature through the + Clausius-Clapeyron equation. There are several different formulations based + on empirical relations that are used to calculate the vapour pressure. Here, + we calculate the vapour pressure based on [Tetens30]_ equation: + + .. math:: e = e_0 * exp\left(17.27\frac{T - 273.16}{T - 35.86}\right) + + for e the vapour pressure, :math:`e_0` the saturation vapour pressure at + a reference temperature of 273.15 K with a value of 6.1078 hPa, and T the + temperature. + + The temperature is provided in Kelvin, and the vapour pressure is returned + in hPa. + + If the (dry-bulb) temperature is used the value given by the vapour pressure + will be the saturation vapour pressure. On the other hand, if the dewpoint + temperature is used the vapour pressure will be calculated. + + References + ---------- + .. [Tetens30] Tetens, O. (1930) "Über einige meteorologische Begriffe". + Z. Geophys 6: 297-309. + + Examples + -------- + >>> vapour_pressure = pressure.vapour_pressure(temperature) + """ + v_pressure = iris.cube.CubeList([]) + for T in iter_maybe(temperature): + es = T.copy() + exponent = 17.27 * (T - 273.16) / (T - 35.86) + es.data = E0 * np.exp(exponent.core_data()) + es.units = "hPa" + es.rename("vapour_pressure") + v_pressure.append(es) + if len(v_pressure) == 1: + return v_pressure[0] + else: + return v_pressure + + +def vapour_pressure_from_relative_humidity( + temperature: iris.cube.Cube | iris.cube.CubeList, + relative_humidity: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Calculate the vapour pressure using relative humidity RH. + + Arguments + --------- + temperature: iris.cube.Cube | iris.cube.CubeList + Cubes of temperature to be converted into saturation vapour pressure. + Temperature must be provided in Kelvin. + relative_humidity: iris.cube.Cube | iris.cube.CubeList + Cubes of relative humidity to be converted into vapour pressure. + Relative humidity must be provided in percentage and + will be converted to a decimal by the operator. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + Vapour pressure in hPa. + + Notes + ----- + The vapour pressure can be derived from the relative humidity and + temperature based on the following relation + + .. math:: e = RH * e_s + + for e the vapour pressure, :math:`e_s` the saturation vapour pressure, + and RH the relative humidity. + + The relative humidity is converted to a decimal. The saturation vapour + pressure is calculated using `pressure.vapour_pressure`. + + Examples + -------- + >>> vapour_pressure = pressure.vapour_pressure_from_relative_humidity(temperature, relative_humidity) + """ + v_pressure = iris.cube.CubeList([]) + for T, RH in zip( + iter_maybe(temperature), iter_maybe(relative_humidity), strict=True + ): + RH = convert_units(RH, "1") + vp = vapour_pressure(T) * RH + vp.units = "hPa" + vp.rename("vapour_pressure") + v_pressure.append(vp) + if len(v_pressure) == 1: + return v_pressure[0] + else: + return v_pressure + + +def exner_pressure( + pressure: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Calculate the exner pressure. + + Arguments + --------- + pressure: iris.cube.Cube | iris.cube.CubeList + Cubes of pressure to be converted. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + Exner pressure. + + Notes + ----- + The Exner pressure is also referred to as the Exner function. It is a + dimensionless parameter that can be used either as a vertical coordinate or + more frequently as a means to simplifying conversions between different + parameters (e.g. Temperature and Potential Temperature; [Holton13]_). + It is calculated as + + .. math:: \Pi = \left(\frac{P}{P_0}\right)^{\frac{R_d}{c_p}} + + for :math:`\Pi` the Exner Pressure, P the pressure in hPa, :math:`P_0` a reference pressure + of 1000 hPa, :math:`R_d` the specific gas constant of dry air taken as 297 :math:`J kg^{-1} K^{-1}`, + :math:`c_p` the specific heat capacity at constant pressure taken as 1005.7 :math:`J kg^{-1} K^{-1}`. + + A value below one implies the pressure is higher than the reference pressure; + values above one implies the pressure is lower than the reference pressure; a + value of one implies the pressure is equal to the reference pressure. + + References + ---------- + .. [Holton13] Holton, J.R. and Hakim G.J. (2013) "An Introduction to Dynamic + Meteorology." 5th Edition, Burlington, MA, Elsevier Academic Press, 532 pp. + + Examples + -------- + >>> Exner_pressure = pressure.exner_pressure(pressure) + """ + pi = iris.cube.CubeList([]) + for P in iter_maybe(pressure): + PI = P.copy() + P = convert_units(P, "hPa") + PI.data = (P.core_data() / P0) ** KAPPA + PI.rename("exner_pressure") + PI.units = "1" + pi.append(PI) + if len(pi) == 1: + return pi[0] + else: + return pi diff --git a/src/CSET/operators/temperature.py b/src/CSET/operators/temperature.py new file mode 100644 index 000000000..1ecbcadc7 --- /dev/null +++ b/src/CSET/operators/temperature.py @@ -0,0 +1,567 @@ +# © Crown copyright, Met Office (2022-2026) and CSET contributors. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Operators for temperature conversions.""" + +import iris.cube +import numpy as np + +from CSET._common import iter_maybe +from CSET.operators._atmospheric_constants import CPD, EPSILON, LV, RV, T0 +from CSET.operators.humidity import ( + mixing_ratio_from_relative_humidity, + saturation_mixing_ratio, +) +from CSET.operators.misc import convert_units +from CSET.operators.pressure import ( + exner_pressure, + vapour_pressure_from_relative_humidity, +) + + +def dewpoint_temperature( + temperature: iris.cube.Cube | iris.cube.CubeList, + relative_humidity: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Calculate the dewpoint temperature. + + Arguments + --------- + temperature: iris.cube.Cube | iris.cube.CubeList + Cubes of temperature in Kelvin. + relative_humidity: iris.cube.Cube | iris.cube.CubeList + Cubes of relative humidity. + + Returns + ------- + iris.cube.Cube | iris.cueb.CubeList + Calculated dewpoint temperature in Kelvin. + + Notes + ----- + The dewpoint temperature provides the temperature at which the air must be + cooled for dew to form (i.e. the water vapour condenses to liquid water). + It is calculated based upon [Bolton80]_ + + .. math:: T_d = \frac{243.5 * ln(e) - 440.8}{19.48 - ln(e)} + + for :math:`T_d` the dewpoint temperature and e the vapour pressure. Here the + vapour pressure is calculated from relative humidity using + `pressure.vapour_pressure_from_relative_humidity`. + + The dewpoint temperature is presented when -35.0 C < T < 35.0 C as this is + when the calculation is the most accurate [Bolton80]_ this roughly equates to + exclusion of dewpoints on pressures of 400 hPa and above (e.g. [Flack24]_). + + When :math:`T_d` is equivalent to T the relative humidity will be 100 %. + + All cubes must be on the same grid. + + References + ---------- + .. [Bolton80] Bolton. D. (1980) "The computation of equivalent potential + temperature". Monthly Weather Review, vol. 108, 1046-1053, + doi: 10.1175/1520-0493(1980)108<1046:TCOEPT>2.0.CO;2 + .. [Flack24] Flack, D.L.A. (2024) "Stratification of the vertical spread-skill + relation by radiosonde drift in a convective-scale ensemble." + Atmospheric Science Letters, vol. 25, e1194, doi: 10.1002/asl.1194 + + Examples + -------- + >>> Td = temperature.dewpoint_temperature(T, RH) + """ + Td = iris.cube.CubeList([]) + for T, RH in zip( + iter_maybe(temperature), iter_maybe(relative_humidity), strict=True + ): + vp = vapour_pressure_from_relative_humidity(T, RH) + td = vp.copy() + td.data = ((243.5 * np.log(vp.core_data())) - 440.8) / ( + 19.48 - np.log(vp.core_data()) + ) + td.data[T.data - T0 < -35.0] = np.nan + td.data[T.data - T0 > 35.0] = np.nan + td.units = "Celsius" + td.rename("dewpoint_temperature") + td = convert_units(td, "K") + Td.append(td) + if len(Td) == 1: + return Td[0] + else: + return Td + + +def virtual_temperature( + temperature: iris.cube.Cube | iris.cube.CubeList, + mixing_ratio: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Calculate the virtual temperature. + + Arguments + --------- + temperature: iris.cube.Cube | iris.cube.CubeList + Cubes of temperature in Kelvin. + mixing_ratio: iris.cube.Cube | iris.cube.CubeList + Cubes of mixing ratio. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + Calculated virtual temperature in Kelvin. + + Notes + ----- + The virtual temperature is that required for a dry parcel to have the + same density as a moist parcel at the same pressure. It is calculated + as + + .. math:: T_v = T * \frac{w + \epsilon}{\epsilon (1 + w)} + + for :math:`T_v` the virtual temperature, T the temperature, w the mixing + ratio, and :math:`\epsilon` the ratio between dry and moist air equating + to 0.622. + + T is given in K and w is given in kg kg-1. + All cubes must be on the same grid. + + Examples + -------- + >>> Tv = temperature.virtual_temperature(T, w) + """ + Tv = iris.cube.CubeList([]) + for T, W in zip(iter_maybe(temperature), iter_maybe(mixing_ratio), strict=True): + virT = T * ((W + EPSILON) / (EPSILON * (1 + W))) + virT.rename("virtual_temperature") + Tv.append(virT) + if len(Tv) == 1: + return Tv[0] + else: + return Tv + + +def wet_bulb_temperature( + temperature: iris.cube.Cube | iris.cube.CubeList, + relative_humidity: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Calculate the wet-bulb temperature. + + Arguments + --------- + temperature: iris.cube.Cube | iris.cube.CubeList + Cubes of temperature. + relative_humidity: iris.cube.Cube | iris.cube.CubeList + Cubes of relative humidity. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + Calculated wet-bulb temperature in Kelvin. + + Notes + ----- + The wet-bulb temperature is the temperature the air reaches when all the + water has been evaporated out of it. It can be calculated from temperature in Celsius + and relative humidity in percent following [Stull11]_ + + .. 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 + + for :math:`T_w` the wet-bulb temperature, T the temperature, and RH the relative + humidity. + + The equation is valid for 5% < RH < 99% and -20 C < T < 50 C, and results are + only presented for these values. + + The temperature and relative humidity unit conversions are applied in the + operator. + + All cubes should be on the same grid. + + References + ---------- + .. [Stull11] Stull, R. (2011) "Wet-Bulb Temperature from Relative Humidity + and air temperature." Journal of Applied Meteorology and Climatology, vol. 50, + 2267-2269, doi: 10.1175/JAMC-D-11-0143.1 + + Examples + -------- + >>> Tw = temperature.wet_bulb_temperature(T, RH) + """ + Tw = iris.cube.CubeList([]) + for T, RH in zip( + iter_maybe(temperature), iter_maybe(relative_humidity), strict=True + ): + RH = convert_units(RH, "%") + T = convert_units(T, "Celsius") + wetT = ( + T * np.arctan(0.151977 * (RH.core_data() + 8.313659) ** 0.5) + + np.arctan(T.core_data() + RH.core_data()) + - np.arctan(RH.core_data() - 1.676331) + + 0.00391838 + * (RH.core_data()) ** (3.0 / 2.0) + * np.arctan(0.023101 * RH.core_data()) + - 4.686035 + ) + wetT.rename("wet_bulb_temperature") + wetT = convert_units(wetT, "K") + wetT.data[T.core_data() < -20.0] = np.nan + wetT.data[T.core_data() > 50.0] = np.nan + wetT.data[RH.core_data() < 5.0] = np.nan + wetT.data[RH.core_data() > 99.0] = np.nan + Tw.append(wetT) + if len(Tw) == 1: + return Tw[0] + else: + return Tw + + +def potential_temperature( + temperature: iris.cube.Cube | iris.cube.CubeList, + pressure: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Calculate the potenital temperature. + + Arguments + --------- + temperature: iris.cube.Cube | iris.cube.CubeList + Cubes of temperature. + pressure: iris.cube.Cube | iris.cube.CuebList + Cubes of pressure. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + Calculated potential temperature in Kelvin. + + Notes + ----- + The potential temperature is the temperature an air parcel would reach + if it was moved adiabatically to a reference pressure. Here we use 1000 hPa + as the reference pressure. It is calculated from + + .. math:: \theta = \frac{T}{\Pi} + + for :math:`\theta` the potential temperature, T the temperature, and :math:`\Pi` + the exner pressure. The exner pressure is calculated using `pressure.exner_pressure`. + + Temperature must be in Kelvin. + + All cubes must be on the same grid. + + Examples + -------- + >>> Theta = temperature.potential_temperature(T, P) + """ + theta = iris.cube.CubeList([]) + for T, P in zip(iter_maybe(temperature), iter_maybe(pressure), strict=True): + TH = T / exner_pressure(P) + TH.rename("potential_temperature") + theta.append(TH) + if len(theta) == 1: + return theta[0] + else: + return theta + + +def virtual_potential_temperature( + temperature: iris.cube.Cube | iris.cube.CubeList, + mixing_ratio: iris.cube.Cube | iris.cube.CubeList, + pressure: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Calculate the virtual potential temperature. + + Arguments + --------- + temperature: iris.cube.Cube | iris.cube.CubeList + Cubes of temperature. + mixing_ratio: iris.cube.Cube | iris.cube.CubeList + Cubes of mixing ratio. + pressure: iris.cube.Cube | iris.cube.CubeList + Cubes of pressure. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + Calculated virtual potential temperature in Kelvin. + + Notes + ----- + The virtual potential temperature is mechanistically equivalent to the + potential temperature, except rather than using the (dry-bulb) temperature + the virtual temperature used. The virtual potential temperature is the potential + temperature a parcel would have if its temperature was replaced by its virtual temperature + and then the parcel is brought adiabatically to 1000 hPa. It is calculated as + + .. math:: \theta_v = \frac{T_v}{\Pi} + + for :math:`\theta_v` the virtual potenital temperature, :math:`T_v` the + virtual temperature, and :math:`\Pi` the exner pressure. The exner pressure + is calculated using `pressure.exner_pressure`. + + All cubes must be on the same grid. + + Examples + -------- + >>> Theta_v = temperature.virtual_potenital_temperature(T, W, P) + """ + theta_v = iris.cube.CubeList([]) + for T, W, P in zip( + iter_maybe(temperature), + iter_maybe(mixing_ratio), + iter_maybe(pressure), + strict=True, + ): + TH_V = virtual_temperature(T, W) / exner_pressure(P) + TH_V.rename("virtual_potential_temperature") + theta_v.append(TH_V) + if len(theta_v) == 1: + return theta_v[0] + else: + return theta_v + + +def equivalent_potential_temperature( + temperature: iris.cube.Cube | iris.cube.CubeList, + relative_humidity: iris.cube.Cube | iris.cube.CubeList, + pressure: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Calculate the equivalent potenital temperature. + + Arguments + --------- + temperature: iris.cube.Cube | iris.cube.CubeList + Cubes of temperature. + relative_humidity: iris.cube.Cube | iris.cube.CubeList + Cubes of relative humidity. + pressure: iris.cube.Cube | iris.cube.CubeList + Cubes of pressure. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + Calculated equivalent potential temperature in Kelvin. + + Notes + ----- + The equivalent potential temperature is the temperature an air parcel + would have if it was raised until all of the water vapour in it was + condensed out and then brought adiabatically to the reference pressure, + here taken as 1000 hPa. It is calculated as + + .. math:: \theta_e = \theta * RH^{- \frac{w R_v}{c_p}} * exp\left(\frac{L_v w}{c_p T} \right) + + for :math:`\theta_e` the equivalent potenital temperature, :math:`\theta` the + potential temperature, RH the relative humidity, w the mixing ratio, + :math:`R_v` the specific gas constant of water vapour (461 + :math:`J kg^{-1} K^{-1}`), :math:`c_p` the specific heat capacity of dry air + (1005.7 :math:`J kg^{-1} K^{-1}`), :math:`L_v` the latent heat of vapourization + (2.5 x :math:`10^6 J kg^{-1}`), and T the temperature. + + Potential temperature and temperature in K. + Relative humidity in percentage and will be converted to fraction. + Mixing ratio in kg kg-1 (dimensionless). + + In this operator the mixing ratio is calculated from + `humidity.mixing_ratio_from_relative_humidity` and the potential temperature + is calculated from `temperature.potential_temperature`. + + This equation is a simplification of [Paluch79]_ following [Emanuel94]_, + whilst still holding reasonable accuracy. + + All cubes must be on the same grid. + + References + ---------- + .. [Emanuel94] Emanuel, K.A. (1994) "Atmospheric Convection" Oxford University + Press, 580 pp. + .. [Paluch79] Paluch, I.R., (1979) "The entrainment mechanism in Colorado + Cumuli" Journal of the Atmospheric Sciences, vol. 36, 2467-2478, + doi: 10.1175/1520-0469(1979)036<2467:TEMICC>2.0.CO;2 + + Examples + -------- + >>> Theta_e = temperature.equivalent_potential_temperature(T, RH, P) + """ + theta_e = iris.cube.CubeList([]) + for T, RH, P in zip( + iter_maybe(temperature), + iter_maybe(relative_humidity), + iter_maybe(pressure), + strict=True, + ): + RH = convert_units(RH, "1") + theta = potential_temperature(T, P) + w = mixing_ratio_from_relative_humidity(T, P, RH) + second_term_power = -(w * RV) / CPD + second_term = RH.core_data() ** second_term_power.core_data() + third_term_power = LV * w / (CPD * T) + third_term = np.exp(third_term_power.core_data()) + TH_E = theta * second_term * third_term + TH_E.rename("equivalent_potential_temperature") + TH_E.units = "K" + theta_e.append(TH_E) + if len(theta_e) == 1: + return theta_e[0] + else: + return theta_e + + +def wet_bulb_potential_temperature( + temperature: iris.cube.Cube | iris.cube.CubeList, + relative_humidity: iris.cube.Cube | iris.cube.CubeList, + pressure: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Calculate wet-bulb potential temperature. + + Arguments + --------- + temperature: iris.cube.Cube | iris.cube.CubeList + Cubes of temperature. + relative_humidity: iris.cube.Cube | iris.cube.CubeList + Cubes of relative humidity. + pressure: iris.cube.Cube | iris.cube.CubeList + Cubes of pressure. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + Calculated wet-bulb potential temperature in Kelvin. + + Notes + ----- + The wet-bulb potential temperature represents the temperature an air parcel + would have if cooled adiabatically until saturation and then taken down to + a reference pressure (1000 hPa) whilst conserving the moisture (i.e. along a + psuedoadiabat). It can be calculated either through a series of iterations, + or through empirical relations. Here, we use the [DaviesJones08]_ formulation: + + .. 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) + + for :math:`\theta_w` the wet-bulb potential temperature, :math:`\theta_e` the + equivalent potential temperature, X = :math:`\theta_e` / 273.15 K, and a + series of constants :math:`a_0` = 7.101574, :math:`a_1` = -20.68208, + :math:`a_2` = 16.11182, :math:`a_3` = 2.574631, :math:`a_4` = -5.205688, + :math:`b_1` = -3.552497, :math:`b_2` = 3.781782, :math:`b_3` = -0.6899655, and + :math:`b_4` = -0.5929340. + + The wet-bulb potential temperature calculated for this method is valid for + the range -30 C < :math:`\theta_w` < 50 C and is thus set to NaN outside of + this range. + + All cubes must be on the same grid. + + References + ---------- + .. [DaviesJones08] Davies-Jones, R. (2008) "An Efficient and Accurate + Method for Computing the Wet-Bulb Temperature along Pseudoadiabats" + Monthly Weather Review, vol. 136, 2764-2785, doi: 10.1175/2007MWR2224.1 + + Examples + -------- + >>> Theta_w = temperature.wet_bulb_potential_temperature(T, RH, P) + """ + theta_w = iris.cube.CubeList([]) + for T, RH, P in zip( + iter_maybe(temperature), + iter_maybe(relative_humidity), + iter_maybe(pressure), + strict=True, + ): + TH_E = equivalent_potential_temperature(T, RH, P) + X = TH_E / T0 + X.units = "1" + A0 = 7.101574 + A1 = -20.68208 + A2 = 16.11182 + A3 = 2.574631 + A4 = -5.205688 + B1 = -3.552497 + B2 = 3.781782 + B3 = -0.6899655 + B4 = -0.5929340 + exponent = (A0 + A1 * X + A2 * X**2 + A3 * X**3 + A4 * X**4) / ( + 1.0 + B1 * X + B2 * X**2 + B3 * X**3 + B4 * X**4 + ) + TH_W = TH_E.copy() + TH_W.data[:] = TH_E.core_data() - np.exp(exponent.core_data()) + TH_W.rename("wet_bulb_potential_temperature") + TH_W.data[TH_W.data - T0 < -30.0] = np.nan + TH_W.data[TH_W.data - T0 > 50.0] = np.nan + theta_w.append(TH_W) + if len(theta_w) == 1: + return theta_w[0] + else: + return theta_w + + +def saturation_equivalent_potential_temperature( + temperature: iris.cube.Cube | iris.cube.CubeList, + pressure: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + r"""Calculate the saturation equivalent potential temperature. + + Arguments + --------- + temperature: iris.cube.Cube | iris.cube.CubeList + Cubes of temperature. + pressure: iris.cube.Cube | iris.cube.CubeList + Cubes of pressure. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + Calculated saturation equivalent potenital temperature in Kelvin. + + Notes + ----- + The saturation equivalent potential temperature, also referred to as the + saturation potential temperature is as the equivalent potential temperature + following a saturated process throughout. It is calculated as + + .. math:: \theta_{es} = \theta * exp\left(\frac{L_v w}{c_p T} \right) + + for :math:`\theta_{es}` the saturation equivalent potential temperature, + :math:`\theta` the potential temperature, w the mixing ratio, + :math:`c_p` the specific heat capacity of dry air (1005.7 + :math:`J kg^{-1} K^{-1}`), :math:`L_v` the latent heat of vapourization + (2.5 x :math:`10^6 J kg^{-1}`), and T the temperature. + + As a saturated process is assumed throughout the RH multiplier in the + equivalent potenital temperature will always be a value of one, and is thus + omitted. In this operator the potential temperature is calculated from + `temperature.potential_temperature`. + + All cubes must be on the same grid. + + Examples + -------- + >>> theta_es = temperature.saturation_equivalent_potential_temperature(T, P) + """ + theta_es = iris.cube.CubeList([]) + for T, P in zip( + iter_maybe(temperature), + iter_maybe(pressure), + strict=True, + ): + theta = potential_temperature(T, P) + ws = saturation_mixing_ratio(T, P) + second_term_power = LV * ws / (CPD * T) + second_term = np.exp(second_term_power.core_data()) + TH_ES = theta * second_term + TH_ES.rename("saturation_equivalent_potential_temperature") + TH_ES.units = "K" + theta_es.append(TH_ES) + if len(theta_es) == 1: + return theta_es[0] + else: + return theta_es diff --git a/src/CSET/operators/wind.py b/src/CSET/operators/wind.py index 4ea767f63..a9faf2fe6 100644 --- a/src/CSET/operators/wind.py +++ b/src/CSET/operators/wind.py @@ -42,7 +42,7 @@ def convert_to_beaufort_scale( The relationship used to convert the windspeed from m/s to the Beaufort Scale is an empirical relationship (e.g., [Beer96]_): - .. math:: F = (\frac{v}{0.836})^{2/3} + .. math:: F = \left(\frac{v}{0.836}\right)^{2/3} for F the Beaufort Force, and v the windspeed at 10 m in m/s. diff --git a/tests/conftest.py b/tests/conftest.py index 7bc1c40bb..a04207353 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,4 +1,4 @@ -# © Crown copyright, Met Office (2022-2024) and CSET contributors. +# © Crown copyright, Met Office (2022-2026) and CSET contributors. # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. @@ -323,3 +323,69 @@ def orography_4D_cube_read_only(): def orography_4D_cube(orography_4D_cube_read_only): """Get 4D orography cube to run tests on. It is safe to modify.""" return orography_4D_cube_read_only.copy() + + +@pytest.fixture() +def temperature_for_conversions_cube_read_only(): + """Get temperature cube for conversions to run tests on. It is NOT safe to modify.""" + return read.read_cube("tests/test_data/pressure/air_temperature.nc") + + +@pytest.fixture() +def temperature_for_conversions_cube(temperature_for_conversions_cube_read_only): + """Get temperature cube for conversions to run tests on. It is safe to modify.""" + return temperature_for_conversions_cube_read_only.copy() + + +@pytest.fixture() +def pressure_for_conversions_cube_read_only(): + """Get pressure cube for conversions to run tests on. It is NOT safe to modify.""" + return read.read_cube("tests/test_data/pressure/pressure.nc") + + +@pytest.fixture() +def pressure_for_conversions_cube(pressure_for_conversions_cube_read_only): + """Get pressure cube for conversions to run tests on. It is safe to modify.""" + return pressure_for_conversions_cube_read_only.copy() + + +@pytest.fixture() +def relative_humidity_for_conversions_cube_read_only(): + """Get relative humidity cube for conversions to run tests on. It is NOT safe to modify.""" + return read.read_cube("tests/test_data/pressure/relative_humidity.nc") + + +@pytest.fixture() +def relative_humidity_for_conversions_cube( + relative_humidity_for_conversions_cube_read_only, +): + """Get relative humidity cube for conversions to run tests on. It is safe to modify.""" + return relative_humidity_for_conversions_cube_read_only.copy() + + +@pytest.fixture() +def specific_humidity_for_conversions_cube_read_only(): + """Get specific humidity cube for conversions to run tests on. It is NOT safe to modify.""" + return read.read_cube("tests/test_data/humidity/specific_humidity.nc") + + +@pytest.fixture() +def specific_humidity_for_conversions_cube( + specific_humidity_for_conversions_cube_read_only, +): + """Get specific humidity cube for conversions to run tests on. It is safe to modify.""" + return specific_humidity_for_conversions_cube_read_only.copy() + + +@pytest.fixture() +def mixing_ratio_for_conversions_cube_read_only(): + """Get mixing ratio cube for conversions to run tests on. It is NOT safe to modify.""" + return read.read_cube("tests/test_data/humidity/mixing_ratio.nc") + + +@pytest.fixture() +def mixing_ratio_for_conversions_cube( + mixing_ratio_for_conversions_cube_read_only, +): + """Get mixing ratio cube for conversions to run tests on. It is safe to modify.""" + return mixing_ratio_for_conversions_cube_read_only.copy() diff --git a/tests/operators/test_humidity.py b/tests/operators/test_humidity.py new file mode 100644 index 000000000..9ff94d14d --- /dev/null +++ b/tests/operators/test_humidity.py @@ -0,0 +1,724 @@ +# © Crown copyright, Met Office (2022-2026) and CSET contributors. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Test humidity operators.""" + +import cf_units +import iris.cube +import numpy as np + +from CSET.operators import _atmospheric_constants, humidity, misc, pressure + + +def test_mixing_ratio_from_specific_humidity(specific_humidity_for_conversions_cube): + """Test calculation of mixing ratio from specific humidity.""" + expected_data = specific_humidity_for_conversions_cube / ( + 1 - specific_humidity_for_conversions_cube + ) + assert np.allclose( + expected_data.data, + humidity.mixing_ratio_from_specific_humidity( + specific_humidity_for_conversions_cube + ).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_mixing_ratio_from_specific_humidity_name( + specific_humidity_for_conversions_cube, +): + """Test naming of mixing ratio cube.""" + expected_name = "mixing_ratio" + assert ( + expected_name + == humidity.mixing_ratio_from_specific_humidity( + specific_humidity_for_conversions_cube + ).name() + ) + + +def test_mixing_ratio_from_specific_humidity_unit( + specific_humidity_for_conversions_cube, +): + """Test units for mixing ratio.""" + expected_units = cf_units.Unit("kg/kg") + assert ( + expected_units + == humidity.mixing_ratio_from_specific_humidity( + specific_humidity_for_conversions_cube + ).units + ) + + +def test_mixing_ratio_from_specific_humidity_cubelist( + specific_humidity_for_conversions_cube, +): + """Test calculation of mixing ratio from specific humidity for a CubeList.""" + expected_data = specific_humidity_for_conversions_cube / ( + 1 - specific_humidity_for_conversions_cube + ) + expected_list = iris.cube.CubeList([expected_data, expected_data]) + input_list = iris.cube.CubeList( + [specific_humidity_for_conversions_cube, specific_humidity_for_conversions_cube] + ) + actual_cubelist = humidity.mixing_ratio_from_specific_humidity(input_list) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) + + +def test_mixing_ratio_from_specific_humidity_using_calculated( + mixing_ratio_for_conversions_cube, +): + """Test mixing ratio calculation using calculated q from mixing ratio.""" + calculated_w = humidity.mixing_ratio_from_specific_humidity( + humidity.specific_humidity_from_mixing_ratio(mixing_ratio_for_conversions_cube) + ) + assert np.allclose( + mixing_ratio_for_conversions_cube.data, calculated_w.data, rtol=1e-6, atol=1e-2 + ) + + +def test_specific_humidity_from_mixing_ratio(mixing_ratio_for_conversions_cube): + """Test specific humidity calculation from mixing ratio.""" + expected_data = mixing_ratio_for_conversions_cube / ( + 1 + mixing_ratio_for_conversions_cube + ) + assert np.allclose( + expected_data.data, + humidity.specific_humidity_from_mixing_ratio( + mixing_ratio_for_conversions_cube + ).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_specific_humidity_from_mixing_ratio_name(mixing_ratio_for_conversions_cube): + """Test specific humidity cube is named.""" + expected_name = "specific_humidity" + assert ( + expected_name + == humidity.specific_humidity_from_mixing_ratio( + mixing_ratio_for_conversions_cube + ).name() + ) + + +def test_specific_humidity_from_mixing_ratio_units(mixing_ratio_for_conversions_cube): + """Test units of specific humidity.""" + expected_units = cf_units.Unit("kg/kg") + assert ( + expected_units + == humidity.specific_humidity_from_mixing_ratio( + mixing_ratio_for_conversions_cube + ).units + ) + + +def test_specific_humidity_from_mixing_ratio_cubelist( + mixing_ratio_for_conversions_cube, +): + """Test specific humidity calculation from mixing ratio with a CuebList.""" + expected_data = mixing_ratio_for_conversions_cube / ( + 1 + mixing_ratio_for_conversions_cube + ) + expected_list = iris.cube.CubeList([expected_data, expected_data]) + input_list = iris.cube.CubeList( + [mixing_ratio_for_conversions_cube, mixing_ratio_for_conversions_cube] + ) + actual_cubelist = humidity.specific_humidity_from_mixing_ratio(input_list) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) + + +def test_specific_humidity_from_mixing_ratio_using_calculated( + specific_humidity_for_conversions_cube, +): + """Test specific humidity calculation using calculated w from specific humidity.""" + calculated_q = humidity.specific_humidity_from_mixing_ratio( + humidity.mixing_ratio_from_specific_humidity( + specific_humidity_for_conversions_cube + ) + ) + assert np.allclose( + specific_humidity_for_conversions_cube.data, + calculated_q.data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_saturation_mixing_ratio( + temperature_for_conversions_cube, pressure_for_conversions_cube +): + """Test calculation of saturation mixing ratio.""" + expected_data = ( + _atmospheric_constants.EPSILON + * pressure.vapour_pressure(temperature_for_conversions_cube) + ) / ( + (misc.convert_units(pressure_for_conversions_cube, "hPa")) + - pressure.vapour_pressure(temperature_for_conversions_cube) + ) + assert np.allclose( + expected_data.data, + humidity.saturation_mixing_ratio( + temperature_for_conversions_cube, pressure_for_conversions_cube + ).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_saturation_mixing_ratio_name( + temperature_for_conversions_cube, pressure_for_conversions_cube +): + """Test naming of saturation mixing ratio cube.""" + expected_name = "saturation_mixing_ratio" + assert ( + expected_name + == humidity.saturation_mixing_ratio( + temperature_for_conversions_cube, pressure_for_conversions_cube + ).name() + ) + + +def test_saturation_mixing_ratio_units( + temperature_for_conversions_cube, pressure_for_conversions_cube +): + """Test units of saturation mixing ratio cube.""" + expected_units = cf_units.Unit("kg/kg") + assert ( + expected_units + == humidity.saturation_mixing_ratio( + temperature_for_conversions_cube, pressure_for_conversions_cube + ).units + ) + + +def test_saturation_mixing_ratio_cubelist( + temperature_for_conversions_cube, pressure_for_conversions_cube +): + """Test calculation of saturation mixing ratio as CubeList.""" + expected_data = ( + _atmospheric_constants.EPSILON + * pressure.vapour_pressure(temperature_for_conversions_cube) + ) / ( + (misc.convert_units(pressure_for_conversions_cube, "hPa")) + - pressure.vapour_pressure(temperature_for_conversions_cube) + ) + expected_list = iris.cube.CubeList([expected_data, expected_data]) + temperature_list = iris.cube.CubeList( + [temperature_for_conversions_cube, temperature_for_conversions_cube] + ) + pressure_list = iris.cube.CubeList( + [pressure_for_conversions_cube, pressure_for_conversions_cube] + ) + actual_cubelist = humidity.saturation_mixing_ratio(temperature_list, pressure_list) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) + + +def test_saturation_specific_humidity( + temperature_for_conversions_cube, pressure_for_conversions_cube +): + """Test calculation of saturation specific humidity.""" + expected_data = ( + _atmospheric_constants.EPSILON + * pressure.vapour_pressure(temperature_for_conversions_cube) + ) / misc.convert_units(pressure_for_conversions_cube, "hPa") + assert np.allclose( + expected_data.data, + humidity.saturation_specific_humidity( + temperature_for_conversions_cube, pressure_for_conversions_cube + ).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_saturation_specific_humidity_name( + temperature_for_conversions_cube, pressure_for_conversions_cube +): + """Test naming of saturation specific humidity cube.""" + expected_name = "saturation_specific_humidity" + assert ( + expected_name + == humidity.saturation_specific_humidity( + temperature_for_conversions_cube, pressure_for_conversions_cube + ).name() + ) + + +def test_saturation_specific_humidity_units( + temperature_for_conversions_cube, pressure_for_conversions_cube +): + """Test units of saturation specific humidity.""" + expected_units = cf_units.Unit("kg/kg") + assert ( + expected_units + == humidity.saturation_specific_humidity( + temperature_for_conversions_cube, pressure_for_conversions_cube + ).units + ) + + +def test_saturation_specific_humidity_cubelist( + temperature_for_conversions_cube, pressure_for_conversions_cube +): + """Test calculation of saturation specific humidity for a CubeList.""" + expected_data = ( + _atmospheric_constants.EPSILON + * pressure.vapour_pressure(temperature_for_conversions_cube) + ) / misc.convert_units(pressure_for_conversions_cube, "hPa") + expected_list = iris.cube.CubeList([expected_data, expected_data]) + temperature_list = iris.cube.CubeList( + [temperature_for_conversions_cube, temperature_for_conversions_cube] + ) + pressure_list = iris.cube.CubeList( + [pressure_for_conversions_cube, pressure_for_conversions_cube] + ) + actual_cubelist = humidity.saturation_specific_humidity( + temperature_list, pressure_list + ) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) + + +def test_mixing_ratio_from_relative_humidity( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + relative_humidity_for_conversions_cube, +): + """Test calculation of mixing ratio from relative humidity.""" + expected_data = humidity.saturation_mixing_ratio( + temperature_for_conversions_cube, pressure_for_conversions_cube + ) * misc.convert_units(relative_humidity_for_conversions_cube, "1") + assert np.allclose( + expected_data.data, + humidity.mixing_ratio_from_relative_humidity( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + relative_humidity_for_conversions_cube, + ).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_mixing_ratio_from_relative_humidity_name( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + relative_humidity_for_conversions_cube, +): + """Test naming of mixing ratio cube.""" + expected_name = "mixing_ratio" + assert ( + expected_name + == humidity.mixing_ratio_from_relative_humidity( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + relative_humidity_for_conversions_cube, + ).name() + ) + + +def test_mixing_ratio_from_relative_humidity_units( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + relative_humidity_for_conversions_cube, +): + """Test units of mixing ratio cube.""" + expected_units = cf_units.Unit("kg/kg") + assert ( + expected_units + == humidity.mixing_ratio_from_relative_humidity( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + relative_humidity_for_conversions_cube, + ).units + ) + + +def test_mixing_ratio_from_relative_humidity_cubelist( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + relative_humidity_for_conversions_cube, +): + """Test calculation of mixing ratio from relative humidity for a CubeList.""" + expected_data = humidity.saturation_mixing_ratio( + temperature_for_conversions_cube, pressure_for_conversions_cube + ) * misc.convert_units(relative_humidity_for_conversions_cube, "1") + expected_list = iris.cube.CubeList([expected_data, expected_data]) + temperature_list = iris.cube.CubeList( + [temperature_for_conversions_cube, temperature_for_conversions_cube] + ) + pressure_list = iris.cube.CubeList( + [pressure_for_conversions_cube, pressure_for_conversions_cube] + ) + rh_list = iris.cube.CubeList( + [relative_humidity_for_conversions_cube, relative_humidity_for_conversions_cube] + ) + actual_cubelist = humidity.mixing_ratio_from_relative_humidity( + temperature_list, pressure_list, rh_list + ) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) + + +def test_mixing_ratio_from_relative_humidity_calculated( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + mixing_ratio_for_conversions_cube, +): + """Test mixing ratio from relative humidity using calculated relative humidity.""" + calculated_w = humidity.mixing_ratio_from_relative_humidity( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + humidity.relative_humidity_from_mixing_ratio( + mixing_ratio_for_conversions_cube, + temperature_for_conversions_cube, + pressure_for_conversions_cube, + ), + ) + assert np.allclose( + calculated_w.data, + mixing_ratio_for_conversions_cube.data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_specific_humidity_from_relative_humidity( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + relative_humidity_for_conversions_cube, +): + """Test calculation of specific humidity from relative humidity.""" + expected_data = humidity.saturation_specific_humidity( + temperature_for_conversions_cube, pressure_for_conversions_cube + ) * misc.convert_units(relative_humidity_for_conversions_cube, "1") + assert np.allclose( + expected_data.data, + humidity.specific_humidity_from_relative_humidity( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + relative_humidity_for_conversions_cube, + ).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_specific_humidity_from_relative_humidity_name( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + relative_humidity_for_conversions_cube, +): + """Test naming of specific humidity cube.""" + expected_name = "specific_humidity" + assert ( + expected_name + == humidity.specific_humidity_from_relative_humidity( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + relative_humidity_for_conversions_cube, + ).name() + ) + + +def test_specific_humidity_from_relative_humidity_units( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + relative_humidity_for_conversions_cube, +): + """Test units of specific humidity cube.""" + expected_units = cf_units.Unit("kg/kg") + assert ( + expected_units + == humidity.specific_humidity_from_relative_humidity( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + relative_humidity_for_conversions_cube, + ).units + ) + + +def test_specific_humidity_from_relative_humidity_cubelist( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + relative_humidity_for_conversions_cube, +): + """Test calculation of specific humidity from relative humidity for a CubeList.""" + expected_data = humidity.saturation_specific_humidity( + temperature_for_conversions_cube, pressure_for_conversions_cube + ) * misc.convert_units(relative_humidity_for_conversions_cube, "1") + expected_list = iris.cube.CubeList([expected_data, expected_data]) + temperature_list = iris.cube.CubeList( + [temperature_for_conversions_cube, temperature_for_conversions_cube] + ) + pressure_list = iris.cube.CubeList( + [pressure_for_conversions_cube, pressure_for_conversions_cube] + ) + rh_list = iris.cube.CubeList( + [relative_humidity_for_conversions_cube, relative_humidity_for_conversions_cube] + ) + actual_cubelist = humidity.specific_humidity_from_relative_humidity( + temperature_list, pressure_list, rh_list + ) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) + + +def test_specific_humidity_from_relative_humidity_calculated( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + specific_humidity_for_conversions_cube, +): + """Test specific humidity from relative humidity using calculated relative humidity.""" + calculated_q = humidity.specific_humidity_from_relative_humidity( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + humidity.relative_humidity_from_specific_humidity( + specific_humidity_for_conversions_cube, + temperature_for_conversions_cube, + pressure_for_conversions_cube, + ), + ) + assert np.allclose( + calculated_q.data, + specific_humidity_for_conversions_cube.data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_relative_humidity_from_mixing_ratio( + mixing_ratio_for_conversions_cube, + temperature_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test calculation of relative humidity from mixing ratio.""" + expected_data = 100.0 * ( + mixing_ratio_for_conversions_cube + / humidity.saturation_mixing_ratio( + temperature_for_conversions_cube, pressure_for_conversions_cube + ) + ) + assert np.allclose( + expected_data.data, + humidity.relative_humidity_from_mixing_ratio( + mixing_ratio_for_conversions_cube, + temperature_for_conversions_cube, + pressure_for_conversions_cube, + ).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_relative_humidity_from_mixing_ratio_calculated( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + relative_humidity_for_conversions_cube, +): + """Test relative humidity from mixing ratio using calculated mixing ratio.""" + calculated_rh = humidity.relative_humidity_from_mixing_ratio( + humidity.mixing_ratio_from_relative_humidity( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + relative_humidity_for_conversions_cube, + ), + temperature_for_conversions_cube, + pressure_for_conversions_cube, + ) + assert np.allclose( + calculated_rh.data, + relative_humidity_for_conversions_cube.data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_relative_humidity_from_mixing_ratio_name( + mixing_ratio_for_conversions_cube, + temperature_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test name of relative humidity cube.""" + expected_name = "relative_humidity" + assert ( + expected_name + == humidity.relative_humidity_from_mixing_ratio( + mixing_ratio_for_conversions_cube, + temperature_for_conversions_cube, + pressure_for_conversions_cube, + ).name() + ) + + +def test_relative_humidity_from_mixing_ratio_units( + mixing_ratio_for_conversions_cube, + temperature_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test units of relative humidity cube.""" + expected_units = cf_units.Unit("%") + assert ( + expected_units + == humidity.relative_humidity_from_mixing_ratio( + mixing_ratio_for_conversions_cube, + temperature_for_conversions_cube, + pressure_for_conversions_cube, + ).units + ) + + +def test_relative_humidity_from_mixing_ratio_cubelist( + mixing_ratio_for_conversions_cube, + temperature_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test calculation of relative humidity from mixing ratio for a CubeList.""" + expected_data = 100.0 * ( + mixing_ratio_for_conversions_cube + / humidity.saturation_mixing_ratio( + temperature_for_conversions_cube, pressure_for_conversions_cube + ) + ) + expected_list = iris.cube.CubeList([expected_data, expected_data]) + mixing_ratio_list = iris.cube.CubeList( + [mixing_ratio_for_conversions_cube, mixing_ratio_for_conversions_cube] + ) + temperature_list = iris.cube.CubeList( + [temperature_for_conversions_cube, temperature_for_conversions_cube] + ) + pressure_list = iris.cube.CubeList( + [pressure_for_conversions_cube, pressure_for_conversions_cube] + ) + actual_cubelist = humidity.relative_humidity_from_mixing_ratio( + mixing_ratio_list, temperature_list, pressure_list + ) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) + + +def test_relative_humidity_from_specific_humidity( + specific_humidity_for_conversions_cube, + temperature_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test calculation of relative humidity from specific humidity.""" + expected_data = 100.0 * ( + specific_humidity_for_conversions_cube + / humidity.saturation_specific_humidity( + temperature_for_conversions_cube, pressure_for_conversions_cube + ) + ) + assert np.allclose( + expected_data.data, + humidity.relative_humidity_from_specific_humidity( + specific_humidity_for_conversions_cube, + temperature_for_conversions_cube, + pressure_for_conversions_cube, + ).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_relative_humidity_from_specific_humidity_calculated( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + relative_humidity_for_conversions_cube, +): + """Test relative humidity from specific humidity using calculated specific humidity.""" + calculated_rh = humidity.relative_humidity_from_specific_humidity( + humidity.specific_humidity_from_relative_humidity( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + relative_humidity_for_conversions_cube, + ), + temperature_for_conversions_cube, + pressure_for_conversions_cube, + ) + assert np.allclose( + calculated_rh.data, + relative_humidity_for_conversions_cube.data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_relative_humidity_from_specific_humidity_name( + specific_humidity_for_conversions_cube, + temperature_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test name of relative humidity cube.""" + expected_name = "relative_humidity" + assert ( + expected_name + == humidity.relative_humidity_from_specific_humidity( + specific_humidity_for_conversions_cube, + temperature_for_conversions_cube, + pressure_for_conversions_cube, + ).name() + ) + + +def test_relative_humidity_from_specific_humidity_units( + specific_humidity_for_conversions_cube, + temperature_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test units of relative humidity cube.""" + expected_units = cf_units.Unit("%") + assert ( + expected_units + == humidity.relative_humidity_from_specific_humidity( + specific_humidity_for_conversions_cube, + temperature_for_conversions_cube, + pressure_for_conversions_cube, + ).units + ) + + +def test_relative_humidity_from_specific_humidity_cubelist( + specific_humidity_for_conversions_cube, + temperature_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test calculation of relative humidity from specific humidity for a CubeList.""" + expected_data = 100.0 * ( + specific_humidity_for_conversions_cube + / humidity.saturation_specific_humidity( + temperature_for_conversions_cube, pressure_for_conversions_cube + ) + ) + expected_list = iris.cube.CubeList([expected_data, expected_data]) + specific_humidity_list = iris.cube.CubeList( + [specific_humidity_for_conversions_cube, specific_humidity_for_conversions_cube] + ) + temperature_list = iris.cube.CubeList( + [temperature_for_conversions_cube, temperature_for_conversions_cube] + ) + pressure_list = iris.cube.CubeList( + [pressure_for_conversions_cube, pressure_for_conversions_cube] + ) + actual_cubelist = humidity.relative_humidity_from_specific_humidity( + specific_humidity_list, temperature_list, pressure_list + ) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) diff --git a/tests/operators/test_pressure.py b/tests/operators/test_pressure.py new file mode 100644 index 000000000..42e7d78ff --- /dev/null +++ b/tests/operators/test_pressure.py @@ -0,0 +1,183 @@ +# © Crown copyright, Met Office (2022-2026) and CSET contributors. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Test pressure operators.""" + +import cf_units +import iris.cube +import numpy as np + +from CSET.operators import _atmospheric_constants, pressure + + +def test_vapour_pressure(temperature_for_conversions_cube): + """Test calculation of vapour pressure for a cube.""" + expected_data = temperature_for_conversions_cube.copy() + exponent = ( + 17.27 + * (temperature_for_conversions_cube - 273.16) + / (temperature_for_conversions_cube - 35.86) + ) + expected_data.data = _atmospheric_constants.E0 * np.exp(exponent.core_data()) + assert np.allclose( + expected_data.data, + pressure.vapour_pressure(temperature_for_conversions_cube).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_vapour_pressure_name(temperature_for_conversions_cube): + """Test renaming of vapour pressure cube.""" + expected_name = "vapour_pressure" + assert ( + expected_name + == pressure.vapour_pressure(temperature_for_conversions_cube).name() + ) + + +def test_vapour_pressure_units(temperature_for_conversions_cube): + """Test units for vapour pressure cube.""" + expected_units = cf_units.Unit("hPa") + assert ( + expected_units + == pressure.vapour_pressure(temperature_for_conversions_cube).units + ) + + +def test_vapour_pressure_cube_list(temperature_for_conversions_cube): + """Test calculation of vapour pressure for a CubeList.""" + expected_data = temperature_for_conversions_cube.copy() + exponent = ( + 17.27 + * (temperature_for_conversions_cube - 273.16) + / (temperature_for_conversions_cube - 35.86) + ) + expected_data.data = _atmospheric_constants.E0 * np.exp(exponent.core_data()) + expected_list = iris.cube.CubeList([expected_data, expected_data]) + input_cubelist = iris.cube.CubeList( + [temperature_for_conversions_cube, temperature_for_conversions_cube] + ) + actual_cubelist = pressure.vapour_pressure(input_cubelist) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) + + +def test_vapour_pressure_from_relative_humidity( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube +): + """Test calculation of vapour pressure from relative humidity.""" + expected_data = pressure.vapour_pressure(temperature_for_conversions_cube) * ( + relative_humidity_for_conversions_cube / 100.0 + ) + assert np.allclose( + expected_data.data, + pressure.vapour_pressure_from_relative_humidity( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube + ).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_vapour_pressure_from_relative_humidity_name( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube +): + """Test naming of vapour pressure cube.""" + expected_name = "vapour_pressure" + assert ( + expected_name + == pressure.vapour_pressure_from_relative_humidity( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube + ).name() + ) + + +def test_vapour_pressure_from_relative_humidity_units( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube +): + """Test units of vapour pressure cube.""" + expected_units = cf_units.Unit("hPa") + assert ( + expected_units + == pressure.vapour_pressure_from_relative_humidity( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube + ).units + ) + + +def test_vapour_pressure_from_relative_humidity_cubelist( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube +): + """Test calculation of vapour pressure from relative humidity as a CubeList.""" + expected_data = pressure.vapour_pressure(temperature_for_conversions_cube) * ( + relative_humidity_for_conversions_cube / 100.0 + ) + expected_list = iris.cube.CubeList([expected_data, expected_data]) + temperature_input_list = iris.cube.CubeList( + [temperature_for_conversions_cube, temperature_for_conversions_cube] + ) + relative_humidity_input_list = iris.cube.CubeList( + [relative_humidity_for_conversions_cube, relative_humidity_for_conversions_cube] + ) + actual_cubelist = pressure.vapour_pressure_from_relative_humidity( + temperature_input_list, relative_humidity_input_list + ) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) + + +def test_exner_pressure(pressure_for_conversions_cube): + """Test calculation of exner pressure.""" + expected_data = pressure_for_conversions_cube.copy() + expected_data.data = ( + (pressure_for_conversions_cube / 100).core_data() / _atmospheric_constants.P0 + ) ** _atmospheric_constants.KAPPA + assert np.allclose( + expected_data.data, + pressure.exner_pressure(pressure_for_conversions_cube).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_exner_pressure_name(pressure_for_conversions_cube): + """Test naming of exner pressure.""" + expected_name = "exner_pressure" + assert ( + expected_name == pressure.exner_pressure(pressure_for_conversions_cube).name() + ) + + +def test_exner_pressure_units(pressure_for_conversions_cube): + """Test units of exner pressure.""" + expected_units = cf_units.Unit("1") + assert ( + expected_units == pressure.exner_pressure(pressure_for_conversions_cube).units + ) + + +def test_exner_pressure_cubelist(pressure_for_conversions_cube): + """Test calculation of exner pressure for a CubeList.""" + expected_data = pressure_for_conversions_cube.copy() + expected_data.data = ( + (pressure_for_conversions_cube / 100).core_data() / _atmospheric_constants.P0 + ) ** _atmospheric_constants.KAPPA + expected_list = iris.cube.CubeList([expected_data, expected_data]) + input_list = iris.cube.CubeList( + [pressure_for_conversions_cube, pressure_for_conversions_cube] + ) + actual_cubelist = pressure.exner_pressure(input_list) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) diff --git a/tests/operators/test_temperature.py b/tests/operators/test_temperature.py new file mode 100644 index 000000000..1af6adcae --- /dev/null +++ b/tests/operators/test_temperature.py @@ -0,0 +1,723 @@ +# © Crown copyright, Met Office (2022-2026) and CSET contributors. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Test temperature operators.""" + +import cf_units +import iris.cube +import numpy as np + +from CSET.operators import _atmospheric_constants, humidity, misc, pressure, temperature + + +def test_dewpoint_temperature( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube +): + """Test calculation of dewpoint temperature.""" + vp = pressure.vapour_pressure_from_relative_humidity( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube + ) + expected_data = vp.copy() + expected_data.data = (243.5 * np.log(vp.core_data()) - 440.8) / ( + 19.48 - np.log(vp.core_data()) + ) + expected_data.data[ + temperature_for_conversions_cube.data - _atmospheric_constants.T0 < -35.0 + ] = np.nan + expected_data.data[ + temperature_for_conversions_cube.data - _atmospheric_constants.T0 > 35.0 + ] = np.nan + expected_data.units = "Celsius" + expected_data = misc.convert_units(expected_data, "K") + assert np.allclose( + expected_data.data, + temperature.dewpoint_temperature( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube + ).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_dewpoint_temperature_name( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube +): + """Test name of dewpoint temperature cube.""" + expected_name = "dewpoint_temperature" + assert ( + expected_name + == temperature.dewpoint_temperature( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube + ).name() + ) + + +def test_dewpoint_temperature_unit( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube +): + """Test units for dewpoint temperature cube.""" + expected_units = cf_units.Unit("K") + assert ( + expected_units + == temperature.dewpoint_temperature( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube + ).units + ) + + +def test_dewpoint_temperature_cubelist( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube +): + """Test calculation of dewpoint temperature for a CubeList.""" + vp = pressure.vapour_pressure_from_relative_humidity( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube + ) + expected_data = vp.copy() + expected_data.data = (243.5 * np.log(vp.core_data()) - 440.8) / ( + 19.48 - np.log(vp.core_data()) + ) + expected_data.data[ + temperature_for_conversions_cube.data - _atmospheric_constants.T0 < -35.0 + ] = np.nan + expected_data.data[ + temperature_for_conversions_cube.data - _atmospheric_constants.T0 > 35.0 + ] = np.nan + expected_data.units = "Celsius" + expected_data = misc.convert_units(expected_data, "K") + expected_list = iris.cube.CubeList([expected_data, expected_data]) + temperature_list = iris.cube.CubeList( + [temperature_for_conversions_cube, temperature_for_conversions_cube] + ) + rh_list = iris.cube.CubeList( + [relative_humidity_for_conversions_cube, relative_humidity_for_conversions_cube] + ) + actual_cubelist = temperature.dewpoint_temperature(temperature_list, rh_list) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) + + +def test_virtual_temperature( + temperature_for_conversions_cube, mixing_ratio_for_conversions_cube +): + """Test to calculate virtual temperature.""" + expected_data = temperature_for_conversions_cube * ( + (mixing_ratio_for_conversions_cube + _atmospheric_constants.EPSILON) + / (_atmospheric_constants.EPSILON * (1 + mixing_ratio_for_conversions_cube)) + ) + assert np.allclose( + expected_data.data, + temperature.virtual_temperature( + temperature_for_conversions_cube, mixing_ratio_for_conversions_cube + ).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_virtual_temperature_name( + temperature_for_conversions_cube, mixing_ratio_for_conversions_cube +): + """Test name of virtual temperature cube.""" + expected_name = "virtual_temperature" + assert ( + expected_name + == temperature.virtual_temperature( + temperature_for_conversions_cube, mixing_ratio_for_conversions_cube + ).name() + ) + + +def test_virtual_temperature_units( + temperature_for_conversions_cube, mixing_ratio_for_conversions_cube +): + """Test units of virtual temperature cube.""" + expected_units = cf_units.Unit("K") + assert ( + expected_units + == temperature.virtual_temperature( + temperature_for_conversions_cube, mixing_ratio_for_conversions_cube + ).units + ) + + +def test_virtual_temperature_cubelist( + temperature_for_conversions_cube, mixing_ratio_for_conversions_cube +): + """Test to calculate virtual temperature for a CubeList.""" + expected_data = temperature_for_conversions_cube * ( + (mixing_ratio_for_conversions_cube + _atmospheric_constants.EPSILON) + / (_atmospheric_constants.EPSILON * (1 + mixing_ratio_for_conversions_cube)) + ) + expected_list = iris.cube.CubeList([expected_data, expected_data]) + temperature_list = iris.cube.CubeList( + [temperature_for_conversions_cube, temperature_for_conversions_cube] + ) + mixing_ratio_list = iris.cube.CubeList( + [mixing_ratio_for_conversions_cube, mixing_ratio_for_conversions_cube] + ) + actual_cubelist = temperature.virtual_temperature( + temperature_list, mixing_ratio_list + ) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) + + +def test_wet_bulb_temperature( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube +): + """Test to calculate wet-bulb temperature.""" + T_units = misc.convert_units(temperature_for_conversions_cube, "Celsius") + RH_units = misc.convert_units(relative_humidity_for_conversions_cube, "%") + TW = ( + T_units * np.arctan(0.151977 * (RH_units.core_data() + 8.313659) ** 0.5) + + np.arctan(T_units.core_data() + RH_units.core_data()) + - np.arctan(RH_units.core_data() - 1.676331) + + 0.00391838 + * (RH_units.core_data()) ** (3.0 / 2.0) + * np.arctan(0.023101 * RH_units.core_data()) + - 4.686035 + ) + TW = misc.convert_units(TW, "K") + TW.data[T_units.data < -20.0] = np.nan + TW.data[T_units.data > 50.0] = np.nan + TW.data[RH_units.data < 5.0] = np.nan + TW.data[RH_units.data > 99.0] = np.nan + assert np.allclose( + TW.data, + temperature.wet_bulb_temperature( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube + ).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_wet_bulb_temperature_units( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube +): + """Test units for wet bulb temperature cube.""" + expected_units = cf_units.Unit("K") + assert ( + expected_units + == temperature.wet_bulb_temperature( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube + ).units + ) + + +def test_wet_bulb_temperature_name( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube +): + """Test name for wet bulb temperature cube.""" + expected_name = "wet_bulb_temperature" + assert ( + expected_name + == temperature.wet_bulb_temperature( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube + ).name() + ) + + +def test_wet_bulb_temperature_cubelist( + temperature_for_conversions_cube, relative_humidity_for_conversions_cube +): + """Test to calculate wet-bulb temperature for a CubeList.""" + T_units = misc.convert_units(temperature_for_conversions_cube, "Celsius") + RH_units = misc.convert_units(relative_humidity_for_conversions_cube, "%") + TW = ( + T_units * np.arctan(0.151977 * (RH_units.core_data() + 8.313659) ** 0.5) + + np.arctan(T_units.core_data() + RH_units.core_data()) + - np.arctan(RH_units.core_data() - 1.676331) + + 0.00391838 + * (RH_units.core_data()) ** (3.0 / 2.0) + * np.arctan(0.023101 * RH_units.core_data()) + - 4.686035 + ) + TW = misc.convert_units(TW, "K") + TW.data[T_units.data < -20.0] = np.nan + TW.data[T_units.data > 50.0] = np.nan + TW.data[RH_units.data < 5.0] = np.nan + TW.data[RH_units.data > 99.0] = np.nan + expected_list = iris.cube.CubeList([TW, TW]) + temperature_list = iris.cube.CubeList( + [temperature_for_conversions_cube, temperature_for_conversions_cube] + ) + rh_list = iris.cube.CubeList( + [relative_humidity_for_conversions_cube, relative_humidity_for_conversions_cube] + ) + actual_cubelist = temperature.wet_bulb_temperature(temperature_list, rh_list) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) + + +def test_potential_temperature( + temperature_for_conversions_cube, pressure_for_conversions_cube +): + """Test to calculate potenital temperature.""" + expected_data = temperature_for_conversions_cube / pressure.exner_pressure( + pressure_for_conversions_cube + ) + assert np.allclose( + expected_data.data, + temperature.potential_temperature( + temperature_for_conversions_cube, pressure_for_conversions_cube + ).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_potenital_temperature_name( + temperature_for_conversions_cube, pressure_for_conversions_cube +): + """Test name of potential temperature cube.""" + expected_name = "potential_temperature" + assert ( + expected_name + == temperature.potential_temperature( + temperature_for_conversions_cube, pressure_for_conversions_cube + ).name() + ) + + +def test_potential_temperature_units( + temperature_for_conversions_cube, pressure_for_conversions_cube +): + """Test units of potential temperature cueb.""" + expected_units = cf_units.Unit("K") + assert ( + expected_units + == temperature.potential_temperature( + temperature_for_conversions_cube, pressure_for_conversions_cube + ).units + ) + + +def test_potenital_temperature_cubelist( + temperature_for_conversions_cube, pressure_for_conversions_cube +): + """Test calculation of potenital temperature for a CubeList.""" + expected_data = temperature_for_conversions_cube / pressure.exner_pressure( + pressure_for_conversions_cube + ) + expected_list = iris.cube.CubeList([expected_data, expected_data]) + temperature_list = iris.cube.CubeList( + [temperature_for_conversions_cube, temperature_for_conversions_cube] + ) + pressure_list = iris.cube.CubeList( + [pressure_for_conversions_cube, pressure_for_conversions_cube] + ) + actual_cubelist = temperature.potential_temperature(temperature_list, pressure_list) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) + + +def test_virtual_potential_temperature( + temperature_for_conversions_cube, + mixing_ratio_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test calculation of virtual potential temperature.""" + expected_data = temperature.virtual_temperature( + temperature_for_conversions_cube, mixing_ratio_for_conversions_cube + ) / pressure.exner_pressure(pressure_for_conversions_cube) + assert np.allclose( + expected_data.data, + temperature.virtual_potential_temperature( + temperature_for_conversions_cube, + mixing_ratio_for_conversions_cube, + pressure_for_conversions_cube, + ).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_virtual_potential_temperature_name( + temperature_for_conversions_cube, + mixing_ratio_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test name of virtual potential temperature cube.""" + expected_name = "virtual_potential_temperature" + assert ( + expected_name + == temperature.virtual_potential_temperature( + temperature_for_conversions_cube, + mixing_ratio_for_conversions_cube, + pressure_for_conversions_cube, + ).name() + ) + + +def test_virtual_potential_temperature_units( + temperature_for_conversions_cube, + mixing_ratio_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test units of virtual potenital temperature cube.""" + expected_units = cf_units.Unit("K") + assert ( + expected_units + == temperature.virtual_potential_temperature( + temperature_for_conversions_cube, + mixing_ratio_for_conversions_cube, + pressure_for_conversions_cube, + ).units + ) + + +def test_virtual_potential_temperature_cubelist( + temperature_for_conversions_cube, + mixing_ratio_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test calculation of virtual potenital temperature for a CubeList.""" + expected_data = temperature.virtual_temperature( + temperature_for_conversions_cube, mixing_ratio_for_conversions_cube + ) / pressure.exner_pressure(pressure_for_conversions_cube) + expected_list = iris.cube.CubeList([expected_data, expected_data]) + temperature_list = iris.cube.CubeList( + [temperature_for_conversions_cube, temperature_for_conversions_cube] + ) + mixing_ratio_list = iris.cube.CubeList( + [mixing_ratio_for_conversions_cube, mixing_ratio_for_conversions_cube] + ) + pressure_list = iris.cube.CubeList( + [pressure_for_conversions_cube, pressure_for_conversions_cube] + ) + actual_cubelist = temperature.virtual_potential_temperature( + temperature_list, mixing_ratio_list, pressure_list + ) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) + + +def test_equivalent_potential_temperature( + temperature_for_conversions_cube, + relative_humidity_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test calculation of equivalent potential temperature.""" + RH = misc.convert_units(relative_humidity_for_conversions_cube, "1") + theta = temperature.potential_temperature( + temperature_for_conversions_cube, pressure_for_conversions_cube + ) + w = humidity.mixing_ratio_from_relative_humidity( + temperature_for_conversions_cube, pressure_for_conversions_cube, RH + ) + second_term_power = -(w * _atmospheric_constants.RV) / _atmospheric_constants.CPD + second_term = RH.core_data() ** second_term_power.core_data() + third_term_power = ( + _atmospheric_constants.LV + * w + / (_atmospheric_constants.CPD * temperature_for_conversions_cube) + ) + third_term = np.exp(third_term_power.core_data()) + expected_data = theta * second_term * third_term + assert np.allclose( + expected_data.data, + temperature.equivalent_potential_temperature( + temperature_for_conversions_cube, + relative_humidity_for_conversions_cube, + pressure_for_conversions_cube, + ).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_equivalent_potenital_temperature_name( + temperature_for_conversions_cube, + relative_humidity_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test name of equivalent potential temperature cube.""" + expected_name = "equivalent_potential_temperature" + assert ( + expected_name + == temperature.equivalent_potential_temperature( + temperature_for_conversions_cube, + relative_humidity_for_conversions_cube, + pressure_for_conversions_cube, + ).name() + ) + + +def test_equivalent_potenital_temperature_units( + temperature_for_conversions_cube, + relative_humidity_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test units of equivalent potential temperature cube.""" + expected_units = cf_units.Unit("K") + assert ( + expected_units + == temperature.equivalent_potential_temperature( + temperature_for_conversions_cube, + relative_humidity_for_conversions_cube, + pressure_for_conversions_cube, + ).units + ) + + +def test_equivalent_potential_temperature_cubelist( + temperature_for_conversions_cube, + relative_humidity_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test calculation of equivalent potential temperature for a CubeList.""" + RH = misc.convert_units(relative_humidity_for_conversions_cube, "1") + theta = temperature.potential_temperature( + temperature_for_conversions_cube, pressure_for_conversions_cube + ) + w = humidity.mixing_ratio_from_relative_humidity( + temperature_for_conversions_cube, pressure_for_conversions_cube, RH + ) + second_term_power = -(w * _atmospheric_constants.RV) / _atmospheric_constants.CPD + second_term = RH.core_data() ** second_term_power.core_data() + third_term_power = ( + _atmospheric_constants.LV + * w + / (_atmospheric_constants.CPD * temperature_for_conversions_cube) + ) + third_term = np.exp(third_term_power.core_data()) + expected_data = theta * second_term * third_term + expected_list = iris.cube.CubeList([expected_data, expected_data]) + temperature_list = iris.cube.CubeList( + [temperature_for_conversions_cube, temperature_for_conversions_cube] + ) + rh_list = iris.cube.CubeList( + [relative_humidity_for_conversions_cube, relative_humidity_for_conversions_cube] + ) + pressure_list = iris.cube.CubeList( + [pressure_for_conversions_cube, pressure_for_conversions_cube] + ) + actual_cubelist = temperature.equivalent_potential_temperature( + temperature_list, rh_list, pressure_list + ) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) + + +def test_saturation_equivalent_potential_temperature( + temperature_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test calculation of saturation equivalent potential temperature.""" + theta = temperature.potential_temperature( + temperature_for_conversions_cube, pressure_for_conversions_cube + ) + ws = humidity.saturation_mixing_ratio( + temperature_for_conversions_cube, pressure_for_conversions_cube + ) + third_term_power = ( + _atmospheric_constants.LV + * ws + / (_atmospheric_constants.CPD * temperature_for_conversions_cube) + ) + third_term = np.exp(third_term_power.core_data()) + expected_data = theta * third_term + assert np.allclose( + expected_data.data, + temperature.saturation_equivalent_potential_temperature( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + ).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_saturation_equivalent_potenital_temperature_name( + temperature_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test name of saturation potential temperature cube.""" + expected_name = "saturation_equivalent_potential_temperature" + assert ( + expected_name + == temperature.saturation_equivalent_potential_temperature( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + ).name() + ) + + +def test_saturation_equivalent_potenital_temperature_units( + temperature_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test units of saturation equivalent potential temperature cube.""" + expected_units = cf_units.Unit("K") + assert ( + expected_units + == temperature.saturation_equivalent_potential_temperature( + temperature_for_conversions_cube, + pressure_for_conversions_cube, + ).units + ) + + +def test_saturation_equivalent_potential_temperature_cubelist( + temperature_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test calculation of saturation equivalent potential temperature for a CubeList.""" + theta = temperature.potential_temperature( + temperature_for_conversions_cube, pressure_for_conversions_cube + ) + ws = humidity.saturation_mixing_ratio( + temperature_for_conversions_cube, pressure_for_conversions_cube + ) + third_term_power = ( + _atmospheric_constants.LV + * ws + / (_atmospheric_constants.CPD * temperature_for_conversions_cube) + ) + third_term = np.exp(third_term_power.core_data()) + expected_data = theta * third_term + expected_list = iris.cube.CubeList([expected_data, expected_data]) + temperature_list = iris.cube.CubeList( + [temperature_for_conversions_cube, temperature_for_conversions_cube] + ) + pressure_list = iris.cube.CubeList( + [pressure_for_conversions_cube, pressure_for_conversions_cube] + ) + actual_cubelist = temperature.saturation_equivalent_potential_temperature( + temperature_list, pressure_list + ) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) + + +def test_wet_bulb_potential_temperature( + temperature_for_conversions_cube, + relative_humidity_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test calculation for wet-bulb potential temperature.""" + TH_E = temperature.equivalent_potential_temperature( + temperature_for_conversions_cube, + relative_humidity_for_conversions_cube, + pressure_for_conversions_cube, + ) + X = TH_E / _atmospheric_constants.T0 + X.units = "1" + A0 = 7.101574 + A1 = -20.68208 + A2 = 16.11182 + A3 = 2.574631 + A4 = -5.205688 + B1 = -3.552497 + B2 = 3.781782 + B3 = -0.6899655 + B4 = -0.5929340 + exponent = (A0 + A1 * X + A2 * X**2 + A3 * X**3 + A4 * X**4) / ( + 1.0 + B1 * X + B2 * X**2 + B3 * X**3 + B4 * X**4 + ) + expected_data = TH_E.copy() + expected_data.data[:] = TH_E.core_data() - np.exp(exponent.core_data()) + expected_data.data[expected_data.data - _atmospheric_constants.T0 < -30.0] = np.nan + expected_data.data[expected_data.data - _atmospheric_constants.T0 > 50.0] = np.nan + assert np.allclose( + expected_data.data, + temperature.wet_bulb_potential_temperature( + temperature_for_conversions_cube, + relative_humidity_for_conversions_cube, + pressure_for_conversions_cube, + ).data, + rtol=1e-6, + atol=1e-2, + ) + + +def test_wet_bulb_potential_temperature_name( + temperature_for_conversions_cube, + relative_humidity_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test name of wet-bulb potenital temperature cube.""" + expected_name = "wet_bulb_potential_temperature" + assert ( + expected_name + == temperature.wet_bulb_potential_temperature( + temperature_for_conversions_cube, + relative_humidity_for_conversions_cube, + pressure_for_conversions_cube, + ).name() + ) + + +def test_wet_bulb_potential_temperature_units( + temperature_for_conversions_cube, + relative_humidity_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test units of wet-bulb potenital temperature cube.""" + expected_units = cf_units.Unit("K") + assert ( + expected_units + == temperature.wet_bulb_potential_temperature( + temperature_for_conversions_cube, + relative_humidity_for_conversions_cube, + pressure_for_conversions_cube, + ).units + ) + + +def test_wet_bulb_potential_temperature_cubelist( + temperature_for_conversions_cube, + relative_humidity_for_conversions_cube, + pressure_for_conversions_cube, +): + """Test calculation for wet-bulb potential temperature.""" + TH_E = temperature.equivalent_potential_temperature( + temperature_for_conversions_cube, + relative_humidity_for_conversions_cube, + pressure_for_conversions_cube, + ) + X = TH_E / _atmospheric_constants.T0 + X.units = "1" + A0 = 7.101574 + A1 = -20.68208 + A2 = 16.11182 + A3 = 2.574631 + A4 = -5.205688 + B1 = -3.552497 + B2 = 3.781782 + B3 = -0.6899655 + B4 = -0.5929340 + exponent = (A0 + A1 * X + A2 * X**2 + A3 * X**3 + A4 * X**4) / ( + 1.0 + B1 * X + B2 * X**2 + B3 * X**3 + B4 * X**4 + ) + expected_data = TH_E.copy() + expected_data.data[:] = TH_E.core_data() - np.exp(exponent.core_data()) + expected_data.data[expected_data.data - _atmospheric_constants.T0 < -30.0] = np.nan + expected_data.data[expected_data.data - _atmospheric_constants.T0 > 50.0] = np.nan + expected_list = iris.cube.CubeList([expected_data, expected_data]) + temperature_list = iris.cube.CubeList( + [temperature_for_conversions_cube, temperature_for_conversions_cube] + ) + rh_list = iris.cube.CubeList( + [relative_humidity_for_conversions_cube, relative_humidity_for_conversions_cube] + ) + pressure_list = iris.cube.CubeList( + [pressure_for_conversions_cube, pressure_for_conversions_cube] + ) + actual_cubelist = temperature.wet_bulb_potential_temperature( + temperature_list, rh_list, pressure_list + ) + for cube_a, cube_b in zip(expected_list, actual_cubelist, strict=True): + assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2) diff --git a/tests/test_data/humidity/mixing_ratio.nc b/tests/test_data/humidity/mixing_ratio.nc new file mode 100644 index 000000000..a30b3751b Binary files /dev/null and b/tests/test_data/humidity/mixing_ratio.nc differ diff --git a/tests/test_data/humidity/specific_humidity.nc b/tests/test_data/humidity/specific_humidity.nc new file mode 100644 index 000000000..e85049faf Binary files /dev/null and b/tests/test_data/humidity/specific_humidity.nc differ diff --git a/tests/test_data/pressure/air_temperature.nc b/tests/test_data/pressure/air_temperature.nc new file mode 100644 index 000000000..c49d20d9d Binary files /dev/null and b/tests/test_data/pressure/air_temperature.nc differ diff --git a/tests/test_data/pressure/pressure.nc b/tests/test_data/pressure/pressure.nc new file mode 100644 index 000000000..9d3f1864c Binary files /dev/null and b/tests/test_data/pressure/pressure.nc differ diff --git a/tests/test_data/pressure/relative_humidity.nc b/tests/test_data/pressure/relative_humidity.nc new file mode 100644 index 000000000..908b6e32c Binary files /dev/null and b/tests/test_data/pressure/relative_humidity.nc differ