From 3cc13faf357961557aa2d7f6c35b67ca467747d6 Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Tue, 12 May 2026 15:49:08 +0100 Subject: [PATCH 01/21] Add function to convert w'T' covariance into sensible heat flux W m-2 --- src/CSET/operators/misc.py | 73 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 73 insertions(+) diff --git a/src/CSET/operators/misc.py b/src/CSET/operators/misc.py index 34e59aefb..897696fae 100644 --- a/src/CSET/operators/misc.py +++ b/src/CSET/operators/misc.py @@ -615,3 +615,76 @@ def differentiate( return new_cubelist[0] else: return new_cubelist + + +def sensible_heat_units(cubes, **kwargs): + """ + Compute surface upward sensible heat flux from covariances. + + EXPECTS (exactly one of each, preselected upstream): + - wt_covariance_2m + - air_temperature_rtd_1p2m + - barometric_pressure + UM cubes are passed through unchanged. + HEIGHT is treated as a *nominal* height (for UM selection / labelling), + not as a measurement height for Cardington. + """ + from cf_units import Unit + + Cp = 1004.67 # J kg-1 K-1 + Rd = 287.05 # J kg-1 K-1 + cubes = ( + iris.cube.CubeList(cubes) + if not isinstance(cubes, iris.cube.CubeList) + else cubes + ) + + # --- Extract Cardington inputs explicitly listed upstream --- + if "CARDINGTON_VARNAMES" not in kwargs: + raise ValueError("sensible_heat_units requires CARDINGTON_VARNAMES") + + wanted = set(kwargs["CARDINGTON_VARNAMES"].split(",")) + selected = {c.var_name: c for c in cubes if c.var_name in wanted} + missing = wanted - set(selected) + if missing: + raise ValueError( + f"sensible_heat_units missing Cardington inputs: {sorted(missing)}" + ) + + wT = next(v for k, v in selected.items() if k.startswith("wt_covariance_")) + temp = next(v for k, v in selected.items() if k.startswith("air_temperature_rtd_")) + pressure = selected["pressure_barometric"] + + # --- Unit handling --- + temp_K = temp.copy() + if temp_K.units is None or temp_K.units.is_unknown(): + temp_K.units = Unit("degC") + temp_K.convert_units("K") + + pres_Pa = pressure.copy() + if pres_Pa.units is None or pres_Pa.units.is_unknown(): + pres_Pa.units = Unit("hPa") + pres_Pa.convert_units("Pa") + + # --- Compute sensible heat flux --- + rho_air = pres_Pa.data / (Rd * temp_K.data) + shf = wT.copy() + shf.data = Cp * rho_air * wT.data + shf.units = "W m-2" + shf.rename("surface_upward_sensible_heat_flux_cardington") + shf.var_name = "surface_upward_sensible_heat_flux_cardington" + + # --- Metadata: be explicit about mixed heights --- + shf.attributes["model_name"] = wT.attributes.get("model_name") + shf.attributes["cardington_measurement_heights"] = { + "wt_covariance": wT.var_name.split("_")[-1], + "air_temperature": temp.var_name.split("_")[-1], + "air_pressure": "1.2 m", + } + if "HEIGHT" in kwargs: + shf.attributes["nominal_height"] = f"{kwargs['HEIGHT']} m" + + # --- Return: passthrough everything except Cardington inputs, plus SHF --- + out = iris.cube.CubeList(c for c in cubes if c.var_name not in wanted) + out.append(shf) + return out if len(out) > 1 else out[0] From cf0435b9f11017dabe33e9a34f72c80465eb263e Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Fri, 29 May 2026 09:41:28 +0100 Subject: [PATCH 02/21] Add series of tests for sensible heat flux operator --- tests/operators/test_misc.py | 106 +++++++++++++++++++++++++++++++++++ 1 file changed, 106 insertions(+) diff --git a/tests/operators/test_misc.py b/tests/operators/test_misc.py index 41cfc7266..475ddeac6 100644 --- a/tests/operators/test_misc.py +++ b/tests/operators/test_misc.py @@ -23,6 +23,7 @@ import iris.exceptions import numpy as np import pytest +from cf_units import Unit from CSET.operators import misc, read @@ -541,3 +542,108 @@ def test_extract_common_points_nocommonpoints(vertical_profile_cube): misc.extract_common_points( cubes=iris.cube.CubeList([cube1, cube2]), coordinate="pressure" ) + + +def _make_scalar_cube( + value, + var_name, + units=None, + model_name="Cardington", +): + """Make tiny 1x1 cube with var_name, units, and model_name attribute.""" + data = np.array([[value]], dtype=np.float64) + lat = iris.coords.DimCoord([0.0], standard_name="latitude", units="degrees") + lon = iris.coords.DimCoord([0.0], standard_name="longitude", units="degrees") + + cube = iris.cube.Cube( + data, + dim_coords_and_dims=[(lat, 0), (lon, 1)], + units=units if units is not None else Unit("unknown"), + ) + cube.var_name = var_name + cube.attributes["model_name"] = model_name + return cube + + +def test_sensible_heat_units_requires_cardington_varnames(): + """Must provide CARDINGTON_VARNAMES kwarg.""" + wT = _make_scalar_cube(0.1, "wt_covariance_2m", units=Unit("K m s-1")) + t = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m", units=Unit("degC")) + p = _make_scalar_cube(1000.0, "pressure_barometric", units=Unit("hPa")) + with pytest.raises(ValueError, match="requires CARDINGTON_VARNAMES"): + misc.sensible_heat_units([wT, t, p]) + + +def test_sensible_heat_units_missing_required_inputs(): + """If any of the specified Cardington inputs are missing, raise.""" + wT = _make_scalar_cube(0.1, "wt_covariance_2m", units=Unit("K m s-1")) + t = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m", units=Unit("degC")) + # pressure missing + kwargs = { + "CARDINGTON_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric" + } + with pytest.raises(ValueError, match="missing Cardington inputs"): + misc.sensible_heat_units([wT, t], **kwargs) + + +def test_sensible_heat_units(): + """Test that core calculation works.""" + wT = _make_scalar_cube(0.1, "wt_covariance_2m", units=Unit("K m s-1")) + temp = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m", units=Unit("degC")) + press = _make_scalar_cube(1000.0, "pressure_barometric", units=Unit("hPa")) + + kwargs = { + "CARDINGTON_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric", + } + out = misc.sensible_heat_units([wT, temp, press], **kwargs) + arr = out.data + if hasattr(arr, "compute"): + arr = arr.compute() + + Cp = 1004.67 + Rd = 287.05 + T = 20.0 + 273.15 + pPa = 1000.0 * 100.0 + rho = pPa / (Rd * T) + expected = Cp * rho * 0.1 + assert np.isclose(arr[0, 0], expected) + + +def test_sensible_heat_units_returns_cube_when_only_shf_remains(): + """If no cubes remain, function should return single Cube.""" + wT = _make_scalar_cube(0.1, "wt_covariance_2m", units=Unit("K m s-1")) + temp = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m", units=Unit("degC")) + press = _make_scalar_cube(1000.0, "pressure_barometric", units=Unit("hPa")) + + kwargs = { + "CARDINGTON_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric" + } + + out = misc.sensible_heat_units([wT, temp, press], **kwargs) + assert isinstance(out, iris.cube.Cube) + assert out.var_name == "surface_upward_sensible_heat_flux_cardington" + assert out.units == "W m-2" + + +def test_sensible_heat_units_passthrough_and_filtering(): + """Ensure input cubes are removed whilst others are kept.""" + wT = _make_scalar_cube(0.1, "wt_covariance_2m") + temp = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m") + press = _make_scalar_cube(1000.0, "pressure_barometric") + extra = _make_scalar_cube(5.0, "other_var") + kwargs = { + "CARDINGTON_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric", + } + + out = misc.sensible_heat_units([wT, temp, press, extra], **kwargs) + + assert isinstance(out, iris.cube.CubeList) + varnames = [c.var_name for c in out] + # input variables removed + assert "wt_covariance_2m" not in varnames + assert "air_temperature_rtd_1p2m" not in varnames + assert "pressure_barometric" not in varnames + # extra retained + assert "other_var" in varnames + # SHF added + assert "surface_upward_sensible_heat_flux_cardington" in varnames From 5e7fb26d998d498c6f8865807ddfc1ad6e4cee7c Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Fri, 12 Jun 2026 10:39:17 +0100 Subject: [PATCH 03/21] Move sensible heat function into new fluxes.py and associated tests into new test_fluxes.py --- src/CSET/operators/fluxes.py | 133 +++++++++++++++++++++++++++++++++ src/CSET/operators/misc.py | 73 ------------------ tests/operators/test_fluxes.py | 129 ++++++++++++++++++++++++++++++++ tests/operators/test_misc.py | 106 -------------------------- 4 files changed, 262 insertions(+), 179 deletions(-) create mode 100644 src/CSET/operators/fluxes.py create mode 100644 tests/operators/test_fluxes.py diff --git a/src/CSET/operators/fluxes.py b/src/CSET/operators/fluxes.py new file mode 100644 index 000000000..2312a5c1d --- /dev/null +++ b/src/CSET/operators/fluxes.py @@ -0,0 +1,133 @@ +# © Crown copyright, Met Office (2022-2025) 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 to calculate kinematic heat fluxes from covariances.""" + +import iris +import iris.cube + + +def sensible_heat_units(cubes, **kwargs): + """ + Convert covariance measurements into surface upward sensible heat flux. + + This operator computes sensible heat flux (SHF) from turbulent covariance + measurements using: + SHF = ρ * Cp * (w'T') + where: + ρ = air density (derived from pressure and temperature) + Cp = specific heat capacity of dry air + w'T' = vertical wind–temperature covariance + + Parameters + ---------- + cubes : iris.cube.Cube or iterable of Cube + Collection of input cubes. Must include exactly one cube for each of: + - wt_covariance_* + - air_temperature_rtd_* + - pressure_barometric + Other cubes are passed through unchanged. + + kwargs : dict + Must contain: + WT_VARNAMES : str + Comma-separated list of variable names to extract from `cubes`. + Expected to include the three required inputs. + Optional: + HEIGHT : float or int + Nominal height used for metadata tagging only (not used in computations). + + Returns + ------- + iris.cube.Cube or iris.cube.CubeList + Output cubes with: + - All non-WT related cubes passed through unchanged + - One additional cube: + surface_upward_sensible_heat_flux[W m-2] + + Returns a single Cube if input was a single Cube, otherwise a CubeList. + + Assumptions and Notes + --------------------- + - Input cubes must be preselected so that exactly one of each required variable + is present. No disambiguation is performed within this function. + - Measurement heights of inputs may differ; these are recorded in output metadata. + - The `HEIGHT` kwarg is treated as a nominal reporting height only and does not + affect the calculation. + - Units are normalised internally: + - Temperature → K (defaults to degC if unknown) + - Pressure → Pa (defaults to hPa if unknown) + + Raises + ------ + ValueError + If required inputs are missing or WT_VARNAMES is not provided. + """ + from cf_units import Unit + + Cp = 1004.67 # J kg-1 K-1 + Rd = 287.05 # J kg-1 K-1 + cubes = ( + iris.cube.CubeList(cubes) + if not isinstance(cubes, iris.cube.CubeList) + else cubes + ) + + # --- Extract inputs explicitly listed upstream --- + if "WT_VARNAMES" not in kwargs: + raise ValueError("sensible_heat_units requires WT_VARNAMES") + + wanted = set(kwargs["WT_VARNAMES"].split(",")) + selected = {c.var_name: c for c in cubes if c.var_name in wanted} + missing = wanted - set(selected) + if missing: + raise ValueError(f"sensible_heat_units missing inputs: {sorted(missing)}") + + wT = next(v for k, v in selected.items() if k.startswith("wt_covariance_")) + temp = next(v for k, v in selected.items() if k.startswith("air_temperature_rtd_")) + pressure = selected["pressure_barometric"] + + # --- Unit handling --- + temp_K = temp.copy() + if temp_K.units is None or temp_K.units.is_unknown(): + temp_K.units = Unit("degC") + temp_K.convert_units("K") + + pres_Pa = pressure.copy() + if pres_Pa.units is None or pres_Pa.units.is_unknown(): + pres_Pa.units = Unit("hPa") + pres_Pa.convert_units("Pa") + + # --- Compute sensible heat flux --- + rho_air = pres_Pa.data / (Rd * temp_K.data) + shf = wT.copy() + shf.data = Cp * rho_air * wT.data + shf.units = "W m-2" + shf.rename("surface_upward_sensible_heat_flux") + shf.var_name = "surface_upward_sensible_heat_flux" + + # --- Metadata: be explicit about mixed heights --- + shf.attributes["model_name"] = wT.attributes.get("model_name") + shf.attributes["measurement_heights"] = { + "wt_covariance": wT.var_name.split("_")[-1], + "air_temperature": temp.var_name.split("_")[-1], + "air_pressure": "1.2 m", + } + if "HEIGHT" in kwargs: + shf.attributes["nominal_height"] = f"{kwargs['HEIGHT']} m" + + # --- Return: passthrough everything except inputs, plus SHF --- + out = iris.cube.CubeList(c for c in cubes if c.var_name not in wanted) + out.append(shf) + return out if len(out) > 1 else out[0] diff --git a/src/CSET/operators/misc.py b/src/CSET/operators/misc.py index 897696fae..34e59aefb 100644 --- a/src/CSET/operators/misc.py +++ b/src/CSET/operators/misc.py @@ -615,76 +615,3 @@ def differentiate( return new_cubelist[0] else: return new_cubelist - - -def sensible_heat_units(cubes, **kwargs): - """ - Compute surface upward sensible heat flux from covariances. - - EXPECTS (exactly one of each, preselected upstream): - - wt_covariance_2m - - air_temperature_rtd_1p2m - - barometric_pressure - UM cubes are passed through unchanged. - HEIGHT is treated as a *nominal* height (for UM selection / labelling), - not as a measurement height for Cardington. - """ - from cf_units import Unit - - Cp = 1004.67 # J kg-1 K-1 - Rd = 287.05 # J kg-1 K-1 - cubes = ( - iris.cube.CubeList(cubes) - if not isinstance(cubes, iris.cube.CubeList) - else cubes - ) - - # --- Extract Cardington inputs explicitly listed upstream --- - if "CARDINGTON_VARNAMES" not in kwargs: - raise ValueError("sensible_heat_units requires CARDINGTON_VARNAMES") - - wanted = set(kwargs["CARDINGTON_VARNAMES"].split(",")) - selected = {c.var_name: c for c in cubes if c.var_name in wanted} - missing = wanted - set(selected) - if missing: - raise ValueError( - f"sensible_heat_units missing Cardington inputs: {sorted(missing)}" - ) - - wT = next(v for k, v in selected.items() if k.startswith("wt_covariance_")) - temp = next(v for k, v in selected.items() if k.startswith("air_temperature_rtd_")) - pressure = selected["pressure_barometric"] - - # --- Unit handling --- - temp_K = temp.copy() - if temp_K.units is None or temp_K.units.is_unknown(): - temp_K.units = Unit("degC") - temp_K.convert_units("K") - - pres_Pa = pressure.copy() - if pres_Pa.units is None or pres_Pa.units.is_unknown(): - pres_Pa.units = Unit("hPa") - pres_Pa.convert_units("Pa") - - # --- Compute sensible heat flux --- - rho_air = pres_Pa.data / (Rd * temp_K.data) - shf = wT.copy() - shf.data = Cp * rho_air * wT.data - shf.units = "W m-2" - shf.rename("surface_upward_sensible_heat_flux_cardington") - shf.var_name = "surface_upward_sensible_heat_flux_cardington" - - # --- Metadata: be explicit about mixed heights --- - shf.attributes["model_name"] = wT.attributes.get("model_name") - shf.attributes["cardington_measurement_heights"] = { - "wt_covariance": wT.var_name.split("_")[-1], - "air_temperature": temp.var_name.split("_")[-1], - "air_pressure": "1.2 m", - } - if "HEIGHT" in kwargs: - shf.attributes["nominal_height"] = f"{kwargs['HEIGHT']} m" - - # --- Return: passthrough everything except Cardington inputs, plus SHF --- - out = iris.cube.CubeList(c for c in cubes if c.var_name not in wanted) - out.append(shf) - return out if len(out) > 1 else out[0] diff --git a/tests/operators/test_fluxes.py b/tests/operators/test_fluxes.py new file mode 100644 index 000000000..286442eda --- /dev/null +++ b/tests/operators/test_fluxes.py @@ -0,0 +1,129 @@ +# © Crown copyright, Met Office (2022-2025) 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 turbulent flux operators.""" + +import iris +import iris.coords +import iris.cube +import numpy as np +import pytest +from cf_units import Unit + +from CSET.operators import fluxes + + +def _make_scalar_cube( + value, + var_name, + units=None, + model_name="Cardington", +): + """Make tiny 1x1 cube with var_name, units, and model_name attribute.""" + data = np.array([[value]], dtype=np.float64) + lat = iris.coords.DimCoord([0.0], standard_name="latitude", units="degrees") + lon = iris.coords.DimCoord([0.0], standard_name="longitude", units="degrees") + + cube = iris.cube.Cube( + data, + dim_coords_and_dims=[(lat, 0), (lon, 1)], + units=units if units is not None else Unit("unknown"), + ) + cube.var_name = var_name + cube.attributes["model_name"] = model_name + return cube + + +def test_sensible_heat_units_requires_cardington_varnames(): + """Must provide CARDINGTON_VARNAMES kwarg.""" + wT = _make_scalar_cube(0.1, "wt_covariance_2m", units=Unit("K m s-1")) + t = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m", units=Unit("degC")) + p = _make_scalar_cube(1000.0, "pressure_barometric", units=Unit("hPa")) + with pytest.raises(ValueError, match="requires CARDINGTON_VARNAMES"): + fluxes.sensible_heat_units([wT, t, p]) + + +def test_sensible_heat_units_missing_required_inputs(): + """If any of the specified Cardington inputs are missing, raise.""" + wT = _make_scalar_cube(0.1, "wt_covariance_2m", units=Unit("K m s-1")) + t = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m", units=Unit("degC")) + # pressure missing + kwargs = { + "CARDINGTON_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric" + } + with pytest.raises(ValueError, match="missing Cardington inputs"): + fluxes.sensible_heat_units([wT, t], **kwargs) + + +def test_sensible_heat_units(): + """Test that core calculation works.""" + wT = _make_scalar_cube(0.1, "wt_covariance_2m", units=Unit("K m s-1")) + temp = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m", units=Unit("degC")) + press = _make_scalar_cube(1000.0, "pressure_barometric", units=Unit("hPa")) + + kwargs = { + "CARDINGTON_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric", + } + out = fluxes.sensible_heat_units([wT, temp, press], **kwargs) + arr = out.data + if hasattr(arr, "compute"): + arr = arr.compute() + + Cp = 1004.67 + Rd = 287.05 + T = 20.0 + 273.15 + pPa = 1000.0 * 100.0 + rho = pPa / (Rd * T) + expected = Cp * rho * 0.1 + assert np.isclose(arr[0, 0], expected) + + +def test_sensible_heat_units_returns_cube_when_only_shf_remains(): + """If no cubes remain, function should return single Cube.""" + wT = _make_scalar_cube(0.1, "wt_covariance_2m", units=Unit("K m s-1")) + temp = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m", units=Unit("degC")) + press = _make_scalar_cube(1000.0, "pressure_barometric", units=Unit("hPa")) + + kwargs = { + "CARDINGTON_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric" + } + + out = fluxes.sensible_heat_units([wT, temp, press], **kwargs) + assert isinstance(out, iris.cube.Cube) + assert out.var_name == "surface_upward_sensible_heat_flux_cardington" + assert out.units == "W m-2" + + +def test_sensible_heat_units_passthrough_and_filtering(): + """Ensure input cubes are removed whilst others are kept.""" + wT = _make_scalar_cube(0.1, "wt_covariance_2m") + temp = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m") + press = _make_scalar_cube(1000.0, "pressure_barometric") + extra = _make_scalar_cube(5.0, "other_var") + kwargs = { + "CARDINGTON_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric", + } + + out = fluxes.sensible_heat_units([wT, temp, press, extra], **kwargs) + + assert isinstance(out, iris.cube.CubeList) + varnames = [c.var_name for c in out] + # input variables removed + assert "wt_covariance_2m" not in varnames + assert "air_temperature_rtd_1p2m" not in varnames + assert "pressure_barometric" not in varnames + # extra retained + assert "other_var" in varnames + # SHF added + assert "surface_upward_sensible_heat_flux_cardington" in varnames diff --git a/tests/operators/test_misc.py b/tests/operators/test_misc.py index 475ddeac6..41cfc7266 100644 --- a/tests/operators/test_misc.py +++ b/tests/operators/test_misc.py @@ -23,7 +23,6 @@ import iris.exceptions import numpy as np import pytest -from cf_units import Unit from CSET.operators import misc, read @@ -542,108 +541,3 @@ def test_extract_common_points_nocommonpoints(vertical_profile_cube): misc.extract_common_points( cubes=iris.cube.CubeList([cube1, cube2]), coordinate="pressure" ) - - -def _make_scalar_cube( - value, - var_name, - units=None, - model_name="Cardington", -): - """Make tiny 1x1 cube with var_name, units, and model_name attribute.""" - data = np.array([[value]], dtype=np.float64) - lat = iris.coords.DimCoord([0.0], standard_name="latitude", units="degrees") - lon = iris.coords.DimCoord([0.0], standard_name="longitude", units="degrees") - - cube = iris.cube.Cube( - data, - dim_coords_and_dims=[(lat, 0), (lon, 1)], - units=units if units is not None else Unit("unknown"), - ) - cube.var_name = var_name - cube.attributes["model_name"] = model_name - return cube - - -def test_sensible_heat_units_requires_cardington_varnames(): - """Must provide CARDINGTON_VARNAMES kwarg.""" - wT = _make_scalar_cube(0.1, "wt_covariance_2m", units=Unit("K m s-1")) - t = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m", units=Unit("degC")) - p = _make_scalar_cube(1000.0, "pressure_barometric", units=Unit("hPa")) - with pytest.raises(ValueError, match="requires CARDINGTON_VARNAMES"): - misc.sensible_heat_units([wT, t, p]) - - -def test_sensible_heat_units_missing_required_inputs(): - """If any of the specified Cardington inputs are missing, raise.""" - wT = _make_scalar_cube(0.1, "wt_covariance_2m", units=Unit("K m s-1")) - t = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m", units=Unit("degC")) - # pressure missing - kwargs = { - "CARDINGTON_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric" - } - with pytest.raises(ValueError, match="missing Cardington inputs"): - misc.sensible_heat_units([wT, t], **kwargs) - - -def test_sensible_heat_units(): - """Test that core calculation works.""" - wT = _make_scalar_cube(0.1, "wt_covariance_2m", units=Unit("K m s-1")) - temp = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m", units=Unit("degC")) - press = _make_scalar_cube(1000.0, "pressure_barometric", units=Unit("hPa")) - - kwargs = { - "CARDINGTON_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric", - } - out = misc.sensible_heat_units([wT, temp, press], **kwargs) - arr = out.data - if hasattr(arr, "compute"): - arr = arr.compute() - - Cp = 1004.67 - Rd = 287.05 - T = 20.0 + 273.15 - pPa = 1000.0 * 100.0 - rho = pPa / (Rd * T) - expected = Cp * rho * 0.1 - assert np.isclose(arr[0, 0], expected) - - -def test_sensible_heat_units_returns_cube_when_only_shf_remains(): - """If no cubes remain, function should return single Cube.""" - wT = _make_scalar_cube(0.1, "wt_covariance_2m", units=Unit("K m s-1")) - temp = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m", units=Unit("degC")) - press = _make_scalar_cube(1000.0, "pressure_barometric", units=Unit("hPa")) - - kwargs = { - "CARDINGTON_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric" - } - - out = misc.sensible_heat_units([wT, temp, press], **kwargs) - assert isinstance(out, iris.cube.Cube) - assert out.var_name == "surface_upward_sensible_heat_flux_cardington" - assert out.units == "W m-2" - - -def test_sensible_heat_units_passthrough_and_filtering(): - """Ensure input cubes are removed whilst others are kept.""" - wT = _make_scalar_cube(0.1, "wt_covariance_2m") - temp = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m") - press = _make_scalar_cube(1000.0, "pressure_barometric") - extra = _make_scalar_cube(5.0, "other_var") - kwargs = { - "CARDINGTON_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric", - } - - out = misc.sensible_heat_units([wT, temp, press, extra], **kwargs) - - assert isinstance(out, iris.cube.CubeList) - varnames = [c.var_name for c in out] - # input variables removed - assert "wt_covariance_2m" not in varnames - assert "air_temperature_rtd_1p2m" not in varnames - assert "pressure_barometric" not in varnames - # extra retained - assert "other_var" in varnames - # SHF added - assert "surface_upward_sensible_heat_flux_cardington" in varnames From 501d2ac75ec2f9229df11235106c979c9a0c3ab6 Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Fri, 12 Jun 2026 10:54:02 +0100 Subject: [PATCH 04/21] Remove tiny test cube name of 'Cardington' --- tests/operators/test_fluxes.py | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/tests/operators/test_fluxes.py b/tests/operators/test_fluxes.py index 286442eda..17f2f8bb4 100644 --- a/tests/operators/test_fluxes.py +++ b/tests/operators/test_fluxes.py @@ -28,9 +28,8 @@ def _make_scalar_cube( value, var_name, units=None, - model_name="Cardington", ): - """Make tiny 1x1 cube with var_name, units, and model_name attribute.""" + """Make tiny 1x1 cube with var_name and units.""" data = np.array([[value]], dtype=np.float64) lat = iris.coords.DimCoord([0.0], standard_name="latitude", units="degrees") lon = iris.coords.DimCoord([0.0], standard_name="longitude", units="degrees") @@ -41,28 +40,27 @@ def _make_scalar_cube( units=units if units is not None else Unit("unknown"), ) cube.var_name = var_name - cube.attributes["model_name"] = model_name return cube -def test_sensible_heat_units_requires_cardington_varnames(): - """Must provide CARDINGTON_VARNAMES kwarg.""" +def test_sensible_heat_units_requires_WT_Varnames(): + """Must provide WT_VARNAMES kwarg.""" wT = _make_scalar_cube(0.1, "wt_covariance_2m", units=Unit("K m s-1")) t = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m", units=Unit("degC")) p = _make_scalar_cube(1000.0, "pressure_barometric", units=Unit("hPa")) - with pytest.raises(ValueError, match="requires CARDINGTON_VARNAMES"): + with pytest.raises(ValueError, match="requires WT_VARNAMES"): fluxes.sensible_heat_units([wT, t, p]) def test_sensible_heat_units_missing_required_inputs(): - """If any of the specified Cardington inputs are missing, raise.""" + """If any of the specified inputs are missing, raise.""" wT = _make_scalar_cube(0.1, "wt_covariance_2m", units=Unit("K m s-1")) t = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m", units=Unit("degC")) # pressure missing kwargs = { - "CARDINGTON_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric" + "WT_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric" } - with pytest.raises(ValueError, match="missing Cardington inputs"): + with pytest.raises(ValueError, match="missing inputs"): fluxes.sensible_heat_units([wT, t], **kwargs) @@ -73,7 +71,7 @@ def test_sensible_heat_units(): press = _make_scalar_cube(1000.0, "pressure_barometric", units=Unit("hPa")) kwargs = { - "CARDINGTON_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric", + "WT_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric", } out = fluxes.sensible_heat_units([wT, temp, press], **kwargs) arr = out.data @@ -96,12 +94,12 @@ def test_sensible_heat_units_returns_cube_when_only_shf_remains(): press = _make_scalar_cube(1000.0, "pressure_barometric", units=Unit("hPa")) kwargs = { - "CARDINGTON_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric" + "WT_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric" } out = fluxes.sensible_heat_units([wT, temp, press], **kwargs) assert isinstance(out, iris.cube.Cube) - assert out.var_name == "surface_upward_sensible_heat_flux_cardington" + assert out.var_name == "surface_upward_sensible_heat_flux" assert out.units == "W m-2" @@ -112,7 +110,7 @@ def test_sensible_heat_units_passthrough_and_filtering(): press = _make_scalar_cube(1000.0, "pressure_barometric") extra = _make_scalar_cube(5.0, "other_var") kwargs = { - "CARDINGTON_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric", + "WT_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric", } out = fluxes.sensible_heat_units([wT, temp, press, extra], **kwargs) @@ -126,4 +124,4 @@ def test_sensible_heat_units_passthrough_and_filtering(): # extra retained assert "other_var" in varnames # SHF added - assert "surface_upward_sensible_heat_flux_cardington" in varnames + assert "surface_upward_sensible_heat_flux" in varnames From 002297ce81865f33dcbefc2d82dacc26c6f228e2 Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Fri, 12 Jun 2026 13:10:37 +0100 Subject: [PATCH 05/21] Move latent heat operator to fluxes.py and associated tests to test_fluxes.py --- src/CSET/operators/__init__.py | 4 ++ src/CSET/operators/fluxes.py | 70 ++++++++++++++++++++++++++++++++++ tests/operators/test_fluxes.py | 42 ++++++++++++++++++++ 3 files changed, 116 insertions(+) diff --git a/src/CSET/operators/__init__.py b/src/CSET/operators/__init__.py index 7274e39ca..5f4fd7752 100644 --- a/src/CSET/operators/__init__.py +++ b/src/CSET/operators/__init__.py @@ -29,11 +29,13 @@ ageofair, aggregate, aviation, + cardington, collapse, constraints, convection, ensembles, filters, + fluxes, humidity, imageprocessing, mesoscale, @@ -54,12 +56,14 @@ "ageofair", "aggregate", "aviation", + "cardington", "collapse", "constraints", "convection", "ensembles", "execute_recipe", "filters", + "fluxes", "humidity", "get_operator", "imageprocessing", diff --git a/src/CSET/operators/fluxes.py b/src/CSET/operators/fluxes.py index 2312a5c1d..936b141d8 100644 --- a/src/CSET/operators/fluxes.py +++ b/src/CSET/operators/fluxes.py @@ -16,6 +16,10 @@ import iris import iris.cube +from cf_units import Unit +from iris.cube import Cube, CubeList + +from CSET._common import iter_maybe def sensible_heat_units(cubes, **kwargs): @@ -131,3 +135,69 @@ def sensible_heat_units(cubes, **kwargs): out = iris.cube.CubeList(c for c in cubes if c.var_name not in wanted) out.append(shf) return out if len(out) > 1 else out[0] + + +def latent_heat_units( + cubes: Cube | CubeList, + **kwargs, +) -> Cube | CubeList: + """ + Convert covariance into latent heat flux units. + + This operator converts the covariance of vertical wind and + specific humidity (w'q') from mass flux units to latent heat flux (W m-2). + + This function operates on one or more Iris cubes. Any cube with + units convertible to mass flux (kg m-2 s-1) is multiplied by a + constant latent heat of vaporisation to produce a latent heat flux. + Cubes with incompatible, missing, or unknown units are passed through + unchanged. + + Parameters + ---------- + cubes : Cube or CubeList + Input cube(s), typically containing w'q' covariance or other flux-like + quantities. + + **kwargs : dict + Unused; accepted for interface consistency with other operators. + + Returns + ------- + Cube or CubeList + Output cube(s) where: + - Cubes with units convertible to kg m-2 s-1 are converted to W m-2. + - All other cubes are returned unchanged. + - The return type matches the input type (single Cube or CubeList). + + Notes + ----- + - The conversion uses a fixed latent heat of vaporisation: + Lc = 2.45 × 10^6 J kg-1 + - In reality, Lc varies with temperature (~5% variation between -20 °C + and +40 °C). This dependency is currently neglected but could be + included in future improvements. + - This function does not attempt to identify specific variables; it relies + solely on unit convertibility to determine applicability. + """ + REQUIRED_UNITS = Unit("kg m-2 s-1") + OUTPUT_UNITS = Unit("W m-2") + Lc = 2.45e6 # J kg-1 + + out = iris.cube.CubeList() + for cube in iter_maybe(cubes): + # ACT ON MASS FLUXES + if cube.units is None or cube.units.is_unknown(): + out.append(cube) + continue + if not cube.units.is_convertible(REQUIRED_UNITS): + # e.g. if UM LE or some other diagnostic — leave untouched + out.append(cube) + continue + + cube_a = cube.copy() + cube_a = cube_a * Lc + cube_a.units = OUTPUT_UNITS + out.append(cube_a) + + return out[0] if len(out) == 1 else out diff --git a/tests/operators/test_fluxes.py b/tests/operators/test_fluxes.py index 17f2f8bb4..4f95407ab 100644 --- a/tests/operators/test_fluxes.py +++ b/tests/operators/test_fluxes.py @@ -125,3 +125,45 @@ def test_sensible_heat_units_passthrough_and_filtering(): assert "other_var" in varnames # SHF added assert "surface_upward_sensible_heat_flux" in varnames + + +def test_latent_heat_units_conversion(): + """Test unit conversion of latent heat units.""" + wq = _make_scalar_cube(0.001, "wq_covariance", units=Unit("kg m-2 s-1")) + + out = fluxes.latent_heat_units(wq) + arr = out.data + if hasattr(arr, "compute"): + arr = arr.compute() + + expected = 2.45e6 * 0.001 + assert np.isclose(arr[0, 0], expected) + assert out.units == "W m-2" + + +def test_latent_heat_units_passthrough_non_convertible(): + """Test operator returns unchanged cube if not convertible.""" + cube = _make_scalar_cube(5.0, "not_flux", units=Unit("K")) + out = fluxes.latent_heat_units(cube) + # should be unchanged + assert out is cube + + +def test_latent_heat_units_cubelist_mixed(): + """Test operator with mixed cubelist.""" + wq = _make_scalar_cube(0.001, "wq_covariance", units=Unit("kg m-2 s-1")) + other = _make_scalar_cube(10.0, "temperature", units=Unit("K")) + out = fluxes.latent_heat_units([wq, other]) + + assert isinstance(out, iris.cube.CubeList) + assert len(out) == 2 + # check conversion happened only for first cube + converted = next(c for c in out if c.units == "W m-2") + arr = converted.data + if hasattr(arr, "compute"): + arr = arr.compute() + + assert np.isclose(arr[0, 0], 2.45e6 * 0.001) + assert converted.units == "W m-2" + # check passthrough + assert any(c is other for c in out) From 6ead89fce01367722daedbc9716deff48ab512c0 Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Fri, 12 Jun 2026 13:15:41 +0100 Subject: [PATCH 06/21] small mod to __init__.py --- src/CSET/operators/__init__.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/CSET/operators/__init__.py b/src/CSET/operators/__init__.py index 5f4fd7752..5347880c0 100644 --- a/src/CSET/operators/__init__.py +++ b/src/CSET/operators/__init__.py @@ -29,7 +29,6 @@ ageofair, aggregate, aviation, - cardington, collapse, constraints, convection, @@ -56,7 +55,6 @@ "ageofair", "aggregate", "aviation", - "cardington", "collapse", "constraints", "convection", From 84215c3c71108815d6b0233a5fadbe2b06ff063c Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Fri, 12 Jun 2026 13:54:48 +0100 Subject: [PATCH 07/21] Improve heat flux tests for latent heat units --- src/CSET/operators/fluxes.py | 13 ++++++++++--- tests/operators/test_fluxes.py | 16 +++++++++++++--- 2 files changed, 23 insertions(+), 6 deletions(-) diff --git a/src/CSET/operators/fluxes.py b/src/CSET/operators/fluxes.py index 936b141d8..cb6570d81 100644 --- a/src/CSET/operators/fluxes.py +++ b/src/CSET/operators/fluxes.py @@ -144,8 +144,14 @@ def latent_heat_units( """ Convert covariance into latent heat flux units. - This operator converts the covariance of vertical wind and - specific humidity (w'q') from mass flux units to latent heat flux (W m-2). + This operator converts any cube with units convertible to kg m-2 s-1 + (i.e. water mass flux) into latent heat flux (W m-2) by multiplying + by a constant latent heat of vaporisation. + + No attempt is made to distinguish between turbulent fluxes (e.g. w'q') + and other water mass fluxes. This generalisation seems reasonable + given that interpreting rainfall or dewfall, for example, as an + equivalent heat flux is physically meaningful. This function operates on one or more Iris cubes. Any cube with units convertible to mass flux (kg m-2 s-1) is multiplied by a @@ -197,7 +203,8 @@ def latent_heat_units( cube_a = cube.copy() cube_a = cube_a * Lc - cube_a.units = OUTPUT_UNITS + cube_a.units = cube.units * Unit("J kg-1") + cube_a.convert_units(OUTPUT_UNITS) out.append(cube_a) return out[0] if len(out) == 1 else out diff --git a/tests/operators/test_fluxes.py b/tests/operators/test_fluxes.py index 4f95407ab..128b2e5c7 100644 --- a/tests/operators/test_fluxes.py +++ b/tests/operators/test_fluxes.py @@ -138,7 +138,7 @@ def test_latent_heat_units_conversion(): expected = 2.45e6 * 0.001 assert np.isclose(arr[0, 0], expected) - assert out.units == "W m-2" + assert out.units == Unit("W m-2") def test_latent_heat_units_passthrough_non_convertible(): @@ -146,7 +146,8 @@ def test_latent_heat_units_passthrough_non_convertible(): cube = _make_scalar_cube(5.0, "not_flux", units=Unit("K")) out = fluxes.latent_heat_units(cube) # should be unchanged - assert out is cube + assert out == cube + assert out.units == cube.units def test_latent_heat_units_cubelist_mixed(): @@ -158,7 +159,9 @@ def test_latent_heat_units_cubelist_mixed(): assert isinstance(out, iris.cube.CubeList) assert len(out) == 2 # check conversion happened only for first cube - converted = next(c for c in out if c.units == "W m-2") + converted = [c for c in out if c.units.is_convertible("W m-2")] + assert len(converted) == 1 + converted = converted[0] arr = converted.data if hasattr(arr, "compute"): arr = arr.compute() @@ -167,3 +170,10 @@ def test_latent_heat_units_cubelist_mixed(): assert converted.units == "W m-2" # check passthrough assert any(c is other for c in out) + + +def test_latent_heat_units_unknown_units_passthrough(): + """Test operator if units unknown.""" + cube = _make_scalar_cube(1.0, "unknown", units=None) + out = fluxes.latent_heat_units(cube) + assert out is cube From 1e01e82828e470f3aab2e2520dffaed5d4781e26 Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Fri, 12 Jun 2026 15:31:32 +0100 Subject: [PATCH 08/21] Refactored sensible heat flux operator, and associated tests, to make it universal --- src/CSET/operators/fluxes.py | 209 +++++++++++++++++----------- tests/operators/test_fluxes.py | 240 +++++++++++++++++++++++++-------- 2 files changed, 310 insertions(+), 139 deletions(-) diff --git a/src/CSET/operators/fluxes.py b/src/CSET/operators/fluxes.py index cb6570d81..b7dea5c2d 100644 --- a/src/CSET/operators/fluxes.py +++ b/src/CSET/operators/fluxes.py @@ -22,119 +22,164 @@ from CSET._common import iter_maybe +def _exactly_one(matches, role): + if len(matches) == 0: + raise ValueError(f"sensible_heat_units could not identify a unique {role} cube") + if len(matches) > 1: + names = [getattr(c, "var_name", None) or c.name() for c in matches] + raise ValueError( + f"sensible_heat_units found multiple possible {role} cubes: {names}" + ) + return matches[0] + + +def _is_p_cube(cube): + return ( + cube.units is not None + and not cube.units.is_unknown() + and cube.units.is_convertible(Unit("Pa")) + ) + + +def _is_T_cube(cube): + if cube.units is None or cube.units.is_unknown(): + return False + return cube.units.is_convertible(Unit("K")) or cube.units.is_convertible( + Unit("degC") + ) + + +def _is_wt_covar_cube(cube): + if cube.units is None or cube.units.is_unknown(): + return False + # turbulence covariance may be recorded either as K m s-1 + # or degC m s-1; these are equivalent + return cube.units.is_convertible(Unit("K m s-1")) or cube.units.is_convertible( + Unit("degC m s-1") + ) + + +def _score_T_cube(cube): + score = 0 + if cube.standard_name == "air_temperature": + score += 3 + if cube.var_name and "temp" in cube.var_name.lower(): + score += 2 + if "temperature" in cube.name().lower(): + score += 1 + return score + + +def _score_p_cube(cube): + score = 0 + if cube.standard_name == "air_pressure": + score += 3 + if cube.var_name: + v = cube.var_name.lower() + if v == "barometric_pressure": + score += 3 + elif "press" in v: + score += 2 + if cube.name() and "pressure" in cube.name().lower(): + score += 1 + + return score + + +def _score_covar_cube(cube): + score = 0 + text = " ".join( + str(x).lower() + for x in [cube.var_name, cube.standard_name, cube.long_name, cube.name()] + if x + ) + if "cov" in text: + score += 2 + if "wt" in text or "w't" in text: + score += 2 + return score + + def sensible_heat_units(cubes, **kwargs): """ - Convert covariance measurements into surface upward sensible heat flux. + Convert turbulent temperature covariance into sensible heat flux. - This operator computes sensible heat flux (SHF) from turbulent covariance - measurements using: - SHF = ρ * Cp * (w'T') - where: - ρ = air density (derived from pressure and temperature) - Cp = specific heat capacity of dry air - w'T' = vertical wind–temperature covariance + This operator identifies: + - one covariance cube with units compatible with temperature × velocity + - one air temperature cube + - one pressure cube - Parameters - ---------- - cubes : iris.cube.Cube or iterable of Cube - Collection of input cubes. Must include exactly one cube for each of: - - wt_covariance_* - - air_temperature_rtd_* - - pressure_barometric - Other cubes are passed through unchanged. - - kwargs : dict - Must contain: - WT_VARNAMES : str - Comma-separated list of variable names to extract from `cubes`. - Expected to include the three required inputs. - Optional: - HEIGHT : float or int - Nominal height used for metadata tagging only (not used in computations). + It then computes: + SHF = rho * Cp * (w'T') - Returns - ------- - iris.cube.Cube or iris.cube.CubeList - Output cubes with: - - All non-WT related cubes passed through unchanged - - One additional cube: - surface_upward_sensible_heat_flux[W m-2] - - Returns a single Cube if input was a single Cube, otherwise a CubeList. - - Assumptions and Notes - --------------------- - - Input cubes must be preselected so that exactly one of each required variable - is present. No disambiguation is performed within this function. - - Measurement heights of inputs may differ; these are recorded in output metadata. - - The `HEIGHT` kwarg is treated as a nominal reporting height only and does not - affect the calculation. - - Units are normalised internally: - - Temperature → K (defaults to degC if unknown) - - Pressure → Pa (defaults to hPa if unknown) - - Raises - ------ - ValueError - If required inputs are missing or WT_VARNAMES is not provided. + using air density derived from pressure and temperature. + + Cubes not used in the calculation are passed through unchanged. """ from cf_units import Unit Cp = 1004.67 # J kg-1 K-1 Rd = 287.05 # J kg-1 K-1 + cubes = ( iris.cube.CubeList(cubes) if not isinstance(cubes, iris.cube.CubeList) else cubes ) - # --- Extract inputs explicitly listed upstream --- - if "WT_VARNAMES" not in kwargs: - raise ValueError("sensible_heat_units requires WT_VARNAMES") + p_cand = [c for c in cubes if _is_p_cube(c)] + T_cand = [c for c in cubes if _is_T_cube(c)] + covar_cand = [c for c in cubes if _is_wt_covar_cube(c)] + + # Optional: use scoring if more than one candidate exists + if len(p_cand) > 1: + p_cand = sorted(p_cand, key=_score_p_cube, reverse=True) + if len(p_cand) > 1 and _score_p_cube(p_cand[0]) == _score_p_cube(p_cand[1]): + raise ValueError("Multiple plausible pressure cubes found") - wanted = set(kwargs["WT_VARNAMES"].split(",")) - selected = {c.var_name: c for c in cubes if c.var_name in wanted} - missing = wanted - set(selected) - if missing: - raise ValueError(f"sensible_heat_units missing inputs: {sorted(missing)}") + if len(T_cand) > 1: + T_cand = sorted(T_cand, key=_score_T_cube, reverse=True) + if len(T_cand) > 1 and _score_T_cube(T_cand[0]) == _score_T_cube(T_cand[1]): + raise ValueError("Multiple plausible temperature cubes found") - wT = next(v for k, v in selected.items() if k.startswith("wt_covariance_")) - temp = next(v for k, v in selected.items() if k.startswith("air_temperature_rtd_")) - pressure = selected["pressure_barometric"] + if len(covar_cand) > 1: + covar_cand = sorted(covar_cand, key=_score_covar_cube, reverse=True) + if len(covar_cand) > 1 and _score_covar_cube( + covar_cand[0] + ) == _score_covar_cube(covar_cand[1]): + raise ValueError("Multiple plausible covariance cubes found") + + pressure = _exactly_one(p_cand[:1], "pressure") + temp = _exactly_one(T_cand[:1], "temperature") + wT = _exactly_one(covar_cand[:1], "w'T' covariance") - # --- Unit handling --- temp_K = temp.copy() - if temp_K.units is None or temp_K.units.is_unknown(): - temp_K.units = Unit("degC") - temp_K.convert_units("K") + if temp_K.units.is_convertible(Unit("degC")): + temp_K.convert_units("K") pres_Pa = pressure.copy() - if pres_Pa.units is None or pres_Pa.units.is_unknown(): - pres_Pa.units = Unit("hPa") pres_Pa.convert_units("Pa") - # --- Compute sensible heat flux --- + # Treat degC covariance numerically as K covariance for fluctuations + wT_cov = wT.copy() + if str(wT_cov.units) == "degC m s-1": + wT_cov.units = Unit("K m s-1") + rho_air = pres_Pa.data / (Rd * temp_K.data) - shf = wT.copy() - shf.data = Cp * rho_air * wT.data - shf.units = "W m-2" + + shf = wT_cov.copy() + shf.data = Cp * rho_air * wT_cov.data + shf.units = Unit("W m-2") shf.rename("surface_upward_sensible_heat_flux") shf.var_name = "surface_upward_sensible_heat_flux" - - # --- Metadata: be explicit about mixed heights --- - shf.attributes["model_name"] = wT.attributes.get("model_name") - shf.attributes["measurement_heights"] = { - "wt_covariance": wT.var_name.split("_")[-1], - "air_temperature": temp.var_name.split("_")[-1], - "air_pressure": "1.2 m", - } if "HEIGHT" in kwargs: shf.attributes["nominal_height"] = f"{kwargs['HEIGHT']} m" - # --- Return: passthrough everything except inputs, plus SHF --- - out = iris.cube.CubeList(c for c in cubes if c.var_name not in wanted) + used_ids = {id(wT), id(temp), id(pressure)} + out = iris.cube.CubeList(c for c in cubes if id(c) not in used_ids) out.append(shf) - return out if len(out) > 1 else out[0] + + return out[0] if len(out) == 1 else out def latent_heat_units( diff --git a/tests/operators/test_fluxes.py b/tests/operators/test_fluxes.py index 128b2e5c7..bdd8ce36f 100644 --- a/tests/operators/test_fluxes.py +++ b/tests/operators/test_fluxes.py @@ -28,8 +28,10 @@ def _make_scalar_cube( value, var_name, units=None, + standard_name=None, + long_name=None, ): - """Make tiny 1x1 cube with var_name and units.""" + """Make tiny 1x1 cube with metadata and units.""" data = np.array([[value]], dtype=np.float64) lat = iris.coords.DimCoord([0.0], standard_name="latitude", units="degrees") lon = iris.coords.DimCoord([0.0], standard_name="longitude", units="degrees") @@ -38,42 +40,54 @@ def _make_scalar_cube( data, dim_coords_and_dims=[(lat, 0), (lon, 1)], units=units if units is not None else Unit("unknown"), + standard_name=standard_name, + long_name=long_name, ) cube.var_name = var_name return cube -def test_sensible_heat_units_requires_WT_Varnames(): - """Must provide WT_VARNAMES kwarg.""" - wT = _make_scalar_cube(0.1, "wt_covariance_2m", units=Unit("K m s-1")) - t = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m", units=Unit("degC")) - p = _make_scalar_cube(1000.0, "pressure_barometric", units=Unit("hPa")) - with pytest.raises(ValueError, match="requires WT_VARNAMES"): - fluxes.sensible_heat_units([wT, t, p]) - - def test_sensible_heat_units_missing_required_inputs(): - """If any of the specified inputs are missing, raise.""" - wT = _make_scalar_cube(0.1, "wt_covariance_2m", units=Unit("K m s-1")) - t = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m", units=Unit("degC")) + """Raise if one of the required physical inputs cannot be identified.""" + wT = _make_scalar_cube( + 0.1, + "wt_covariance_2m", + units=Unit("K m s-1"), + long_name="vertical wind-temperature covariance", + ) + temp = _make_scalar_cube( + 20.0, + "air_temperature_rtd_1p2m", + units=Unit("degC"), + standard_name="air_temperature", + ) # pressure missing - kwargs = { - "WT_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric" - } - with pytest.raises(ValueError, match="missing inputs"): - fluxes.sensible_heat_units([wT, t], **kwargs) - - -def test_sensible_heat_units(): - """Test that core calculation works.""" - wT = _make_scalar_cube(0.1, "wt_covariance_2m", units=Unit("K m s-1")) - temp = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m", units=Unit("degC")) - press = _make_scalar_cube(1000.0, "pressure_barometric", units=Unit("hPa")) - - kwargs = { - "WT_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric", - } - out = fluxes.sensible_heat_units([wT, temp, press], **kwargs) + with pytest.raises(ValueError, match="pressure"): + fluxes.sensible_heat_units([wT, temp]) + + +def test_sensible_heat_units_core_calculation(): + """Compute sensible heat flux from covariance, temperature, and pressure.""" + wT = _make_scalar_cube( + 0.1, + "wt_covariance_2m", + units=Unit("K m s-1"), + long_name="vertical wind-temperature covariance", + ) + temp = _make_scalar_cube( + 20.0, + "air_temperature_rtd_1p2m", + units=Unit("degC"), + standard_name="air_temperature", + ) + press = _make_scalar_cube( + 1000.0, + "barometric_pressure", + units=Unit("hPa"), + long_name="barometric pressure", + ) + + out = fluxes.sensible_heat_units([wT, temp, press]) arr = out.data if hasattr(arr, "compute"): arr = arr.compute() @@ -85,51 +99,163 @@ def test_sensible_heat_units(): rho = pPa / (Rd * T) expected = Cp * rho * 0.1 assert np.isclose(arr[0, 0], expected) + assert out.var_name == "surface_upward_sensible_heat_flux" + assert out.units == Unit("W m-2") def test_sensible_heat_units_returns_cube_when_only_shf_remains(): - """If no cubes remain, function should return single Cube.""" - wT = _make_scalar_cube(0.1, "wt_covariance_2m", units=Unit("K m s-1")) - temp = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m", units=Unit("degC")) - press = _make_scalar_cube(1000.0, "pressure_barometric", units=Unit("hPa")) - - kwargs = { - "WT_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric" - } + """Return a single Cube if only the derived SHF cube remains.""" + wT = _make_scalar_cube( + 0.1, + "wt_covariance_2m", + units=Unit("K m s-1"), + long_name="vertical wind-temperature covariance", + ) + temp = _make_scalar_cube( + 20.0, + "air_temperature_rtd_1p2m", + units=Unit("degC"), + standard_name="air_temperature", + ) + press = _make_scalar_cube( + 1000.0, + "barometric_pressure", + units=Unit("hPa"), + long_name="barometric pressure", + ) - out = fluxes.sensible_heat_units([wT, temp, press], **kwargs) + out = fluxes.sensible_heat_units([wT, temp, press]) assert isinstance(out, iris.cube.Cube) assert out.var_name == "surface_upward_sensible_heat_flux" - assert out.units == "W m-2" + assert out.units == Unit("W m-2") def test_sensible_heat_units_passthrough_and_filtering(): - """Ensure input cubes are removed whilst others are kept.""" - wT = _make_scalar_cube(0.1, "wt_covariance_2m") - temp = _make_scalar_cube(20.0, "air_temperature_rtd_1p2m") - press = _make_scalar_cube(1000.0, "pressure_barometric") - extra = _make_scalar_cube(5.0, "other_var") - kwargs = { - "WT_VARNAMES": "wt_covariance_2m,air_temperature_rtd_1p2m,pressure_barometric", - } - - out = fluxes.sensible_heat_units([wT, temp, press, extra], **kwargs) + """Remove only the cubes used in SHF calculation and retain unrelated cubes.""" + wT = _make_scalar_cube( + 0.1, + "wt_covariance_2m", + units=Unit("K m s-1"), + long_name="vertical wind-temperature covariance", + ) + temp = _make_scalar_cube( + 20.0, + "air_temperature_rtd_1p2m", + units=Unit("degC"), + standard_name="air_temperature", + ) + press = _make_scalar_cube( + 1000.0, + "barometric_pressure", + units=Unit("hPa"), + long_name="barometric pressure", + ) + extra = _make_scalar_cube( + 5.0, + "other_var", + units=Unit("1"), + ) + out = fluxes.sensible_heat_units([wT, temp, press, extra]) assert isinstance(out, iris.cube.CubeList) varnames = [c.var_name for c in out] # input variables removed assert "wt_covariance_2m" not in varnames assert "air_temperature_rtd_1p2m" not in varnames - assert "pressure_barometric" not in varnames - # extra retained + assert "barometric_pressure" not in varnames + # unrelated cube retained assert "other_var" in varnames # SHF added assert "surface_upward_sensible_heat_flux" in varnames +def test_sensible_heat_units_adds_nominal_height_metadata(): + """Attach nominal_height metadata when HEIGHT is provided.""" + wT = _make_scalar_cube( + 0.1, + "wt_covariance_2m", + units=Unit("K m s-1"), + long_name="vertical wind-temperature covariance", + ) + temp = _make_scalar_cube( + 20.0, + "air_temperature_rtd_1p2m", + units=Unit("degC"), + standard_name="air_temperature", + ) + press = _make_scalar_cube( + 1000.0, + "barometric_pressure", + units=Unit("hPa"), + long_name="barometric pressure", + ) + out = fluxes.sensible_heat_units([wT, temp, press], HEIGHT=2) + assert out.attributes["nominal_height"] == "2 m" + + +def test_sensible_heat_units_prefers_barometric_pressure(): + """Prefer the strongest-scoring pressure cube.""" + wT = _make_scalar_cube( + 0.1, + "wt_covariance_2m", + units=Unit("K m s-1"), + ) + temp = _make_scalar_cube( + 20.0, + "air_temperature_rtd_1p2m", + units=Unit("degC"), + standard_name="air_temperature", + ) + press1 = _make_scalar_cube( + 1000.0, + "barometric_pressure", + units=Unit("hPa"), + ) + press2 = _make_scalar_cube( + 1002.0, + "surface_pressure", + units=Unit("hPa"), + ) + out = fluxes.sensible_heat_units([wT, temp, press1, press2]) + + # extract SHF cube + shf = next(c for c in out if c.var_name == "surface_upward_sensible_heat_flux") + arr = shf.data + if hasattr(arr, "compute"): + arr = arr.compute() + + Cp = 1004.67 + Rd = 287.05 + T = 20.0 + 273.15 + pPa = 1000.0 * 100.0 # should use barometric_pressure + rho = pPa / (Rd * T) + expected = Cp * rho * 0.1 + + assert np.isclose(arr[0, 0], expected) + + +def _make_scalar_latent_cube( + value, + var_name, + units=None, +): + """Make tiny 1x1 cube with var_name and units.""" + data = np.array([[value]], dtype=np.float64) + lat = iris.coords.DimCoord([0.0], standard_name="latitude", units="degrees") + lon = iris.coords.DimCoord([0.0], standard_name="longitude", units="degrees") + + cube = iris.cube.Cube( + data, + dim_coords_and_dims=[(lat, 0), (lon, 1)], + units=units if units is not None else Unit("unknown"), + ) + cube.var_name = var_name + return cube + + def test_latent_heat_units_conversion(): """Test unit conversion of latent heat units.""" - wq = _make_scalar_cube(0.001, "wq_covariance", units=Unit("kg m-2 s-1")) + wq = _make_scalar_latent_cube(0.001, "wq_covariance", units=Unit("kg m-2 s-1")) out = fluxes.latent_heat_units(wq) arr = out.data @@ -143,7 +269,7 @@ def test_latent_heat_units_conversion(): def test_latent_heat_units_passthrough_non_convertible(): """Test operator returns unchanged cube if not convertible.""" - cube = _make_scalar_cube(5.0, "not_flux", units=Unit("K")) + cube = _make_scalar_latent_cube(5.0, "not_flux", units=Unit("K")) out = fluxes.latent_heat_units(cube) # should be unchanged assert out == cube @@ -152,8 +278,8 @@ def test_latent_heat_units_passthrough_non_convertible(): def test_latent_heat_units_cubelist_mixed(): """Test operator with mixed cubelist.""" - wq = _make_scalar_cube(0.001, "wq_covariance", units=Unit("kg m-2 s-1")) - other = _make_scalar_cube(10.0, "temperature", units=Unit("K")) + wq = _make_scalar_latent_cube(0.001, "wq_covariance", units=Unit("kg m-2 s-1")) + other = _make_scalar_latent_cube(10.0, "temperature", units=Unit("K")) out = fluxes.latent_heat_units([wq, other]) assert isinstance(out, iris.cube.CubeList) @@ -174,6 +300,6 @@ def test_latent_heat_units_cubelist_mixed(): def test_latent_heat_units_unknown_units_passthrough(): """Test operator if units unknown.""" - cube = _make_scalar_cube(1.0, "unknown", units=None) + cube = _make_scalar_latent_cube(1.0, "unknown", units=None) out = fluxes.latent_heat_units(cube) assert out is cube From 7c9cca65373af44996b646b92b1f08348aa329fd Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Fri, 12 Jun 2026 15:53:22 +0100 Subject: [PATCH 09/21] Improve doc string in senisble heat flux operator --- src/CSET/operators/fluxes.py | 73 ++++++++++++++++++++++++++++++---- tests/operators/test_fluxes.py | 12 +++--- 2 files changed, 71 insertions(+), 14 deletions(-) diff --git a/src/CSET/operators/fluxes.py b/src/CSET/operators/fluxes.py index b7dea5c2d..a965f737d 100644 --- a/src/CSET/operators/fluxes.py +++ b/src/CSET/operators/fluxes.py @@ -100,21 +100,78 @@ def _score_covar_cube(cube): return score -def sensible_heat_units(cubes, **kwargs): +def sensible_heat_flux_from_covariance(cubes, **kwargs): """ Convert turbulent temperature covariance into sensible heat flux. - This operator identifies: + This operator computes surface upward sensible heat flux (SHF) from + turbulent temperature covariance using: + SHF = ρ * Cp * (w'T') + where: + ρ = air density (derived from pressure and temperature) + Cp = specific heat capacity of dry air + w'T' = covariance between vertical wind and temperature fluctuations + The operator identifies the required input cubes from a collection + based on their physical units, with optional use of metadata (e.g. + standard names, variable names, long names) to resolve ambiguities. + + Specifically, it attempts to identify: - one covariance cube with units compatible with temperature × velocity - - one air temperature cube - - one pressure cube + (e.g. "K m s-1" or "degC m s-1") + - one air temperature cube (units convertible to "K" or "degC") + - one pressure cube (units convertible to "Pa") - It then computes: - SHF = rho * Cp * (w'T') + If multiple plausible candidates exist for a given role, a scoring + system based on metadata is used to select the most likely match. + If ambiguity remains after scoring, a ValueError is raised. - using air density derived from pressure and temperature. + Parameters + ---------- + cubes : iris.cube.Cube or iterable of Cube + Collection of input cubes. Must contain exactly one physically + plausible cube for each of: + - temperature covariance (w'T') + - air temperature + - air pressure + Additional cubes are passed through unchanged. + + **kwargs : dict, optional + Optional keyword arguments. + HEIGHT : float or int + Nominal measurement height (in metres), used only for + metadata annotation of the output cube. - Cubes not used in the calculation are passed through unchanged. + Returns + ------- + iris.cube.Cube or iris.cube.CubeList + Output cubes where: + - Input cubes used in the SHF calculation are removed + - A new cube is added: + surface_upward_sensible_heat_flux [W m-2] + - All other input cubes are passed through unchanged + + If only the SHF cube remains, a single Cube is returned. + Otherwise, a CubeList is returned. + + Assumptions and Notes + --------------------- + - Input cubes are assumed to be on a consistent grid; no spatial + alignment or regridding is performed. + - Temperature is internally converted to Kelvin; pressure is converted + to Pascals. + - Temperature covariance may be expressed in either "degC m s-1" or + "K m s-1". These are treated as numerically equivalent for turbulent + fluctuations (i.e. offsets cancel in fluctuations). + - No attempt is made to distinguish between resolved and turbulent + flux components beyond unit consistency. + - Identification of cubes is primarily unit-based; metadata is used + only to resolve ambiguity between multiple candidates. + + Raises + ------ + ValueError + If required inputs (pressure, temperature, or covariance) cannot + be uniquely identified from the provided cubes. """ from cf_units import Unit diff --git a/tests/operators/test_fluxes.py b/tests/operators/test_fluxes.py index bdd8ce36f..c35fb506b 100644 --- a/tests/operators/test_fluxes.py +++ b/tests/operators/test_fluxes.py @@ -63,7 +63,7 @@ def test_sensible_heat_units_missing_required_inputs(): ) # pressure missing with pytest.raises(ValueError, match="pressure"): - fluxes.sensible_heat_units([wT, temp]) + fluxes.sensible_heat_flux_from_covariance([wT, temp]) def test_sensible_heat_units_core_calculation(): @@ -87,7 +87,7 @@ def test_sensible_heat_units_core_calculation(): long_name="barometric pressure", ) - out = fluxes.sensible_heat_units([wT, temp, press]) + out = fluxes.sensible_heat_flux_from_covariance([wT, temp, press]) arr = out.data if hasattr(arr, "compute"): arr = arr.compute() @@ -124,7 +124,7 @@ def test_sensible_heat_units_returns_cube_when_only_shf_remains(): long_name="barometric pressure", ) - out = fluxes.sensible_heat_units([wT, temp, press]) + out = fluxes.sensible_heat_flux_from_covariance([wT, temp, press]) assert isinstance(out, iris.cube.Cube) assert out.var_name == "surface_upward_sensible_heat_flux" assert out.units == Unit("W m-2") @@ -156,7 +156,7 @@ def test_sensible_heat_units_passthrough_and_filtering(): units=Unit("1"), ) - out = fluxes.sensible_heat_units([wT, temp, press, extra]) + out = fluxes.sensible_heat_flux_from_covariance([wT, temp, press, extra]) assert isinstance(out, iris.cube.CubeList) varnames = [c.var_name for c in out] # input variables removed @@ -189,7 +189,7 @@ def test_sensible_heat_units_adds_nominal_height_metadata(): units=Unit("hPa"), long_name="barometric pressure", ) - out = fluxes.sensible_heat_units([wT, temp, press], HEIGHT=2) + out = fluxes.sensible_heat_flux_from_covariance([wT, temp, press], HEIGHT=2) assert out.attributes["nominal_height"] == "2 m" @@ -216,7 +216,7 @@ def test_sensible_heat_units_prefers_barometric_pressure(): "surface_pressure", units=Unit("hPa"), ) - out = fluxes.sensible_heat_units([wT, temp, press1, press2]) + out = fluxes.sensible_heat_flux_from_covariance([wT, temp, press1, press2]) # extract SHF cube shf = next(c for c in out if c.var_name == "surface_upward_sensible_heat_flux") From 4e08e84621569717aa809fcb8c64d35462aff291 Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Fri, 12 Jun 2026 16:20:36 +0100 Subject: [PATCH 10/21] Use library scientific constants --- src/CSET/operators/fluxes.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/CSET/operators/fluxes.py b/src/CSET/operators/fluxes.py index a965f737d..31d3d7f72 100644 --- a/src/CSET/operators/fluxes.py +++ b/src/CSET/operators/fluxes.py @@ -175,8 +175,8 @@ def sensible_heat_flux_from_covariance(cubes, **kwargs): """ from cf_units import Unit - Cp = 1004.67 # J kg-1 K-1 - Rd = 287.05 # J kg-1 K-1 + from CSET.operators._atmospheric_constants import CPD as Cp + from CSET.operators._atmospheric_constants import RD as Rd cubes = ( iris.cube.CubeList(cubes) @@ -281,7 +281,7 @@ def latent_heat_units( Notes ----- - The conversion uses a fixed latent heat of vaporisation: - Lc = 2.45 × 10^6 J kg-1 + Lc = 2.5 × 10^6 J kg-1 - In reality, Lc varies with temperature (~5% variation between -20 °C and +40 °C). This dependency is currently neglected but could be included in future improvements. @@ -290,7 +290,7 @@ def latent_heat_units( """ REQUIRED_UNITS = Unit("kg m-2 s-1") OUTPUT_UNITS = Unit("W m-2") - Lc = 2.45e6 # J kg-1 + from CSET.operators._atmospheric_constants import LV as Lc out = iris.cube.CubeList() for cube in iter_maybe(cubes): From 6343367bbad195e8251ab5ee078064a7d2e19273 Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Fri, 12 Jun 2026 16:37:49 +0100 Subject: [PATCH 11/21] Use CSET naming consistently for physical constants --- src/CSET/operators/fluxes.py | 19 ++++++---------- tests/operators/test_fluxes.py | 41 ++++++---------------------------- 2 files changed, 14 insertions(+), 46 deletions(-) diff --git a/src/CSET/operators/fluxes.py b/src/CSET/operators/fluxes.py index 31d3d7f72..974663c56 100644 --- a/src/CSET/operators/fluxes.py +++ b/src/CSET/operators/fluxes.py @@ -20,6 +20,7 @@ from iris.cube import Cube, CubeList from CSET._common import iter_maybe +from CSET.operators._atmospheric_constants import CPD, LV, RD def _exactly_one(matches, role): @@ -106,10 +107,10 @@ def sensible_heat_flux_from_covariance(cubes, **kwargs): This operator computes surface upward sensible heat flux (SHF) from turbulent temperature covariance using: - SHF = ρ * Cp * (w'T') + SHF = ρ * CPD * (w'T') where: ρ = air density (derived from pressure and temperature) - Cp = specific heat capacity of dry air + CPD = specific heat capacity of dry air w'T' = covariance between vertical wind and temperature fluctuations The operator identifies the required input cubes from a collection based on their physical units, with optional use of metadata (e.g. @@ -175,9 +176,6 @@ def sensible_heat_flux_from_covariance(cubes, **kwargs): """ from cf_units import Unit - from CSET.operators._atmospheric_constants import CPD as Cp - from CSET.operators._atmospheric_constants import RD as Rd - cubes = ( iris.cube.CubeList(cubes) if not isinstance(cubes, iris.cube.CubeList) @@ -222,15 +220,13 @@ def sensible_heat_flux_from_covariance(cubes, **kwargs): if str(wT_cov.units) == "degC m s-1": wT_cov.units = Unit("K m s-1") - rho_air = pres_Pa.data / (Rd * temp_K.data) + rho_air = pres_Pa.data / (RD * temp_K.data) shf = wT_cov.copy() - shf.data = Cp * rho_air * wT_cov.data + shf.data = CPD * rho_air * wT_cov.data shf.units = Unit("W m-2") shf.rename("surface_upward_sensible_heat_flux") shf.var_name = "surface_upward_sensible_heat_flux" - if "HEIGHT" in kwargs: - shf.attributes["nominal_height"] = f"{kwargs['HEIGHT']} m" used_ids = {id(wT), id(temp), id(pressure)} out = iris.cube.CubeList(c for c in cubes if id(c) not in used_ids) @@ -281,7 +277,7 @@ def latent_heat_units( Notes ----- - The conversion uses a fixed latent heat of vaporisation: - Lc = 2.5 × 10^6 J kg-1 + LV = 2.5 × 10^6 J kg-1 - In reality, Lc varies with temperature (~5% variation between -20 °C and +40 °C). This dependency is currently neglected but could be included in future improvements. @@ -290,7 +286,6 @@ def latent_heat_units( """ REQUIRED_UNITS = Unit("kg m-2 s-1") OUTPUT_UNITS = Unit("W m-2") - from CSET.operators._atmospheric_constants import LV as Lc out = iris.cube.CubeList() for cube in iter_maybe(cubes): @@ -304,7 +299,7 @@ def latent_heat_units( continue cube_a = cube.copy() - cube_a = cube_a * Lc + cube_a = cube_a * LV cube_a.units = cube.units * Unit("J kg-1") cube_a.convert_units(OUTPUT_UNITS) out.append(cube_a) diff --git a/tests/operators/test_fluxes.py b/tests/operators/test_fluxes.py index c35fb506b..ad8c05fdf 100644 --- a/tests/operators/test_fluxes.py +++ b/tests/operators/test_fluxes.py @@ -22,6 +22,7 @@ from cf_units import Unit from CSET.operators import fluxes +from CSET.operators._atmospheric_constants import CPD, LV, RD def _make_scalar_cube( @@ -92,12 +93,10 @@ def test_sensible_heat_units_core_calculation(): if hasattr(arr, "compute"): arr = arr.compute() - Cp = 1004.67 - Rd = 287.05 T = 20.0 + 273.15 pPa = 1000.0 * 100.0 - rho = pPa / (Rd * T) - expected = Cp * rho * 0.1 + rho = pPa / (RD * T) + expected = CPD * rho * 0.1 assert np.isclose(arr[0, 0], expected) assert out.var_name == "surface_upward_sensible_heat_flux" assert out.units == Unit("W m-2") @@ -169,30 +168,6 @@ def test_sensible_heat_units_passthrough_and_filtering(): assert "surface_upward_sensible_heat_flux" in varnames -def test_sensible_heat_units_adds_nominal_height_metadata(): - """Attach nominal_height metadata when HEIGHT is provided.""" - wT = _make_scalar_cube( - 0.1, - "wt_covariance_2m", - units=Unit("K m s-1"), - long_name="vertical wind-temperature covariance", - ) - temp = _make_scalar_cube( - 20.0, - "air_temperature_rtd_1p2m", - units=Unit("degC"), - standard_name="air_temperature", - ) - press = _make_scalar_cube( - 1000.0, - "barometric_pressure", - units=Unit("hPa"), - long_name="barometric pressure", - ) - out = fluxes.sensible_heat_flux_from_covariance([wT, temp, press], HEIGHT=2) - assert out.attributes["nominal_height"] == "2 m" - - def test_sensible_heat_units_prefers_barometric_pressure(): """Prefer the strongest-scoring pressure cube.""" wT = _make_scalar_cube( @@ -224,12 +199,10 @@ def test_sensible_heat_units_prefers_barometric_pressure(): if hasattr(arr, "compute"): arr = arr.compute() - Cp = 1004.67 - Rd = 287.05 T = 20.0 + 273.15 pPa = 1000.0 * 100.0 # should use barometric_pressure - rho = pPa / (Rd * T) - expected = Cp * rho * 0.1 + rho = pPa / (RD * T) + expected = CPD * rho * 0.1 assert np.isclose(arr[0, 0], expected) @@ -262,7 +235,7 @@ def test_latent_heat_units_conversion(): if hasattr(arr, "compute"): arr = arr.compute() - expected = 2.45e6 * 0.001 + expected = LV * 0.001 assert np.isclose(arr[0, 0], expected) assert out.units == Unit("W m-2") @@ -292,7 +265,7 @@ def test_latent_heat_units_cubelist_mixed(): if hasattr(arr, "compute"): arr = arr.compute() - assert np.isclose(arr[0, 0], 2.45e6 * 0.001) + assert np.isclose(arr[0, 0], LV * 0.001) assert converted.units == "W m-2" # check passthrough assert any(c is other for c in out) From 0ba821ac4d35e5675d0177d8d989defcdb823ac2 Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Tue, 16 Jun 2026 06:59:31 +0100 Subject: [PATCH 12/21] Remove realization in SH calculation --- src/CSET/operators/fluxes.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/CSET/operators/fluxes.py b/src/CSET/operators/fluxes.py index 974663c56..48e70873e 100644 --- a/src/CSET/operators/fluxes.py +++ b/src/CSET/operators/fluxes.py @@ -220,10 +220,10 @@ def sensible_heat_flux_from_covariance(cubes, **kwargs): if str(wT_cov.units) == "degC m s-1": wT_cov.units = Unit("K m s-1") - rho_air = pres_Pa.data / (RD * temp_K.data) + rho_air = pres_Pa / (RD * temp_K) shf = wT_cov.copy() - shf.data = CPD * rho_air * wT_cov.data + shf = CPD * rho_air * wT_cov shf.units = Unit("W m-2") shf.rename("surface_upward_sensible_heat_flux") shf.var_name = "surface_upward_sensible_heat_flux" From 94d6358390dd28b7968143a00da77d39867ceb53 Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Tue, 16 Jun 2026 07:26:27 +0100 Subject: [PATCH 13/21] Fix air density calculation for data attributes --- src/CSET/operators/fluxes.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/CSET/operators/fluxes.py b/src/CSET/operators/fluxes.py index 48e70873e..974663c56 100644 --- a/src/CSET/operators/fluxes.py +++ b/src/CSET/operators/fluxes.py @@ -220,10 +220,10 @@ def sensible_heat_flux_from_covariance(cubes, **kwargs): if str(wT_cov.units) == "degC m s-1": wT_cov.units = Unit("K m s-1") - rho_air = pres_Pa / (RD * temp_K) + rho_air = pres_Pa.data / (RD * temp_K.data) shf = wT_cov.copy() - shf = CPD * rho_air * wT_cov + shf.data = CPD * rho_air * wT_cov.data shf.units = Unit("W m-2") shf.rename("surface_upward_sensible_heat_flux") shf.var_name = "surface_upward_sensible_heat_flux" From 1aa940d37e3097a33804e0e105b2f5c85751f99f Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Mon, 6 Jul 2026 11:55:42 +0100 Subject: [PATCH 14/21] Refactor sensible_heat_flux_from_covariance function Refactor and improve documentation for sensible heat flux calculation. Remove _score_ internal functions and do selection of variable names within the sensible heat function. Latent heat function unchanged. --- src/CSET/operators/fluxes.py | 238 ++++++++++++++++------------------- 1 file changed, 107 insertions(+), 131 deletions(-) diff --git a/src/CSET/operators/fluxes.py b/src/CSET/operators/fluxes.py index 974663c56..1ff86d167 100644 --- a/src/CSET/operators/fluxes.py +++ b/src/CSET/operators/fluxes.py @@ -60,120 +60,65 @@ def _is_wt_covar_cube(cube): ) -def _score_T_cube(cube): - score = 0 - if cube.standard_name == "air_temperature": - score += 3 - if cube.var_name and "temp" in cube.var_name.lower(): - score += 2 - if "temperature" in cube.name().lower(): - score += 1 - return score - - -def _score_p_cube(cube): - score = 0 - if cube.standard_name == "air_pressure": - score += 3 - if cube.var_name: - v = cube.var_name.lower() - if v == "barometric_pressure": - score += 3 - elif "press" in v: - score += 2 - if cube.name() and "pressure" in cube.name().lower(): - score += 1 - - return score - - -def _score_covar_cube(cube): - score = 0 - text = " ".join( - str(x).lower() - for x in [cube.var_name, cube.standard_name, cube.long_name, cube.name()] - if x - ) - if "cov" in text: - score += 2 - if "wt" in text or "w't" in text: - score += 2 - return score - - def sensible_heat_flux_from_covariance(cubes, **kwargs): """ - Convert turbulent temperature covariance into sensible heat flux. - - This operator computes surface upward sensible heat flux (SHF) from - turbulent temperature covariance using: - SHF = ρ * CPD * (w'T') - where: - ρ = air density (derived from pressure and temperature) - CPD = specific heat capacity of dry air - w'T' = covariance between vertical wind and temperature fluctuations - The operator identifies the required input cubes from a collection - based on their physical units, with optional use of metadata (e.g. - standard names, variable names, long names) to resolve ambiguities. - - Specifically, it attempts to identify: - - one covariance cube with units compatible with temperature × velocity - (e.g. "K m s-1" or "degC m s-1") - - one air temperature cube (units convertible to "K" or "degC") - - one pressure cube (units convertible to "Pa") - - If multiple plausible candidates exist for a given role, a scoring - system based on metadata is used to select the most likely match. - If ambiguity remains after scoring, a ValueError is raised. +Convert turbulent temperature covariance into sensible heat flux. - Parameters - ---------- - cubes : iris.cube.Cube or iterable of Cube - Collection of input cubes. Must contain exactly one physically - plausible cube for each of: - - temperature covariance (w'T') - - air temperature - - air pressure - Additional cubes are passed through unchanged. - - **kwargs : dict, optional - Optional keyword arguments. - HEIGHT : float or int - Nominal measurement height (in metres), used only for - metadata annotation of the output cube. +This operator computes surface upward sensible heat flux (SHF) from +temperature covariance using: - Returns - ------- - iris.cube.Cube or iris.cube.CubeList - Output cubes where: - - Input cubes used in the SHF calculation are removed - - A new cube is added: - surface_upward_sensible_heat_flux [W m-2] - - All other input cubes are passed through unchanged - - If only the SHF cube remains, a single Cube is returned. - Otherwise, a CubeList is returned. - - Assumptions and Notes - --------------------- - - Input cubes are assumed to be on a consistent grid; no spatial - alignment or regridding is performed. - - Temperature is internally converted to Kelvin; pressure is converted - to Pascals. - - Temperature covariance may be expressed in either "degC m s-1" or - "K m s-1". These are treated as numerically equivalent for turbulent - fluctuations (i.e. offsets cancel in fluctuations). - - No attempt is made to distinguish between resolved and turbulent - flux components beyond unit consistency. - - Identification of cubes is primarily unit-based; metadata is used - only to resolve ambiguity between multiple candidates. - - Raises - ------ - ValueError - If required inputs (pressure, temperature, or covariance) cannot - be uniquely identified from the provided cubes. - """ + SHF = ρ * CPD * (w'T') + +where air density is calculated from pressure and temperature via the +ideal gas law. + +The required input cubes are identified primarily from their physical +units: + + - temperature covariance (e.g. K m s-1 or degC m s-1) + - air temperature (convertible to K or degC) + - air pressure (convertible to Pa) + +If multiple physically plausible candidates are found, CF metadata +(e.g. standard names) are used as a secondary disambiguation step. +A ValueError is raised if the required cubes cannot be uniquely +identified. + +Parameters +---------- +cubes : Cube or CubeList + Input cube(s) containing exactly one identifiable covariance, + temperature and pressure cube. Additional cubes are passed through + unchanged. + +**kwargs : dict, optional + Additional keyword arguments. + +Returns +------- +Cube or CubeList + Input cubes with the pressure, temperature and covariance cubes + removed and a new + ``surface_upward_sensible_heat_flux`` cube added. Unrelated cubes + are passed through unchanged. + +Notes +----- +- Pressure is converted internally to Pa and temperature to K. +- Covariance units of ``degC m s-1`` are treated as numerically + equivalent to ``K m s-1`` because temperature offsets cancel when + forming fluctuations. +- Input cubes are assumed to be physically compatible; no regridding or + coordinate alignment is performed. +- Identification is unit-based, with metadata used only to resolve + ambiguities. + +Raises +------ +ValueError + If suitable pressure, temperature or covariance cubes cannot be + uniquely identified. +""" from cf_units import Unit cubes = ( @@ -182,32 +127,58 @@ def sensible_heat_flux_from_covariance(cubes, **kwargs): else cubes ) + # Pressure cube p_cand = [c for c in cubes if _is_p_cube(c)] - T_cand = [c for c in cubes if _is_T_cube(c)] - covar_cand = [c for c in cubes if _is_wt_covar_cube(c)] - - # Optional: use scoring if more than one candidate exists if len(p_cand) > 1: - p_cand = sorted(p_cand, key=_score_p_cube, reverse=True) - if len(p_cand) > 1 and _score_p_cube(p_cand[0]) == _score_p_cube(p_cand[1]): - raise ValueError("Multiple plausible pressure cubes found") + preferred = [ + c for c in p_cand + if c.standard_name == "air_pressure" + ] + if len(preferred) == 1: + p_cand = preferred + + pressure = _exactly_one(p_cand, "pressure") + + # Temperature cube + T_cand = [c for c in cubes if _is_T_cube(c)] if len(T_cand) > 1: - T_cand = sorted(T_cand, key=_score_T_cube, reverse=True) - if len(T_cand) > 1 and _score_T_cube(T_cand[0]) == _score_T_cube(T_cand[1]): - raise ValueError("Multiple plausible temperature cubes found") + preferred = [ + c for c in T_cand + if c.standard_name == "air_temperature" + ] - if len(covar_cand) > 1: - covar_cand = sorted(covar_cand, key=_score_covar_cube, reverse=True) - if len(covar_cand) > 1 and _score_covar_cube( - covar_cand[0] - ) == _score_covar_cube(covar_cand[1]): - raise ValueError("Multiple plausible covariance cubes found") + if len(preferred) == 1: + T_cand = preferred - pressure = _exactly_one(p_cand[:1], "pressure") - temp = _exactly_one(T_cand[:1], "temperature") - wT = _exactly_one(covar_cand[:1], "w'T' covariance") + temp = _exactly_one(T_cand, "temperature") + # Covariance cube + covar_cand = [c for c in cubes if _is_wt_covar_cube(c)] + if len(covar_cand) > 1: + preferred = [] + for cube in covar_cand: + text = " ".join( + str(x).lower() + for x in ( + cube.standard_name, + cube.var_name, + cube.long_name, + cube.name(), + ) + if x + ) + if "wt" in text or "w't" in text: + preferred.append(cube) + + if len(preferred) == 1: + covar_cand = preferred + + wT = _exactly_one(covar_cand, "w'T' covariance") + + # + # Unit conversions + # temp_K = temp.copy() if temp_K.units.is_convertible(Unit("degC")): temp_K.convert_units("K") @@ -215,7 +186,7 @@ def sensible_heat_flux_from_covariance(cubes, **kwargs): pres_Pa = pressure.copy() pres_Pa.convert_units("Pa") - # Treat degC covariance numerically as K covariance for fluctuations + # Treat degC covariance numerically as K covariance wT_cov = wT.copy() if str(wT_cov.units) == "degC m s-1": wT_cov.units = Unit("K m s-1") @@ -229,7 +200,12 @@ def sensible_heat_flux_from_covariance(cubes, **kwargs): shf.var_name = "surface_upward_sensible_heat_flux" used_ids = {id(wT), id(temp), id(pressure)} - out = iris.cube.CubeList(c for c in cubes if id(c) not in used_ids) + + out = iris.cube.CubeList( + c for c in cubes + if id(c) not in used_ids + ) + out.append(shf) return out[0] if len(out) == 1 else out From 0c373ee371d663fd12bb835fb19c6d368c1049ec Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Mon, 6 Jul 2026 13:46:10 +0100 Subject: [PATCH 15/21] Reconfigured pressure naming logic within operator and expanded tests --- src/CSET/operators/fluxes.py | 93 ++++++++++++++++++---------------- tests/operators/test_fluxes.py | 23 +++++++++ 2 files changed, 71 insertions(+), 45 deletions(-) diff --git a/src/CSET/operators/fluxes.py b/src/CSET/operators/fluxes.py index 1ff86d167..2e2eb43de 100644 --- a/src/CSET/operators/fluxes.py +++ b/src/CSET/operators/fluxes.py @@ -62,63 +62,63 @@ def _is_wt_covar_cube(cube): def sensible_heat_flux_from_covariance(cubes, **kwargs): """ -Convert turbulent temperature covariance into sensible heat flux. + Convert turbulent temperature covariance into sensible heat flux. -This operator computes surface upward sensible heat flux (SHF) from -temperature covariance using: + This operator computes surface upward sensible heat flux (SHF) from + temperature covariance using: SHF = ρ * CPD * (w'T') -where air density is calculated from pressure and temperature via the -ideal gas law. + where air density is calculated from pressure and temperature via the + ideal gas law. -The required input cubes are identified primarily from their physical -units: + The required input cubes are identified primarily from their physical + units: - temperature covariance (e.g. K m s-1 or degC m s-1) - air temperature (convertible to K or degC) - air pressure (convertible to Pa) -If multiple physically plausible candidates are found, CF metadata -(e.g. standard names) are used as a secondary disambiguation step. -A ValueError is raised if the required cubes cannot be uniquely -identified. + If multiple physically plausible candidates are found, CF metadata + (e.g. standard names) are used as a secondary disambiguation step. + A ValueError is raised if the required cubes cannot be uniquely + identified. -Parameters ----------- -cubes : Cube or CubeList + Parameters + ---------- + cubes : Cube or CubeList Input cube(s) containing exactly one identifiable covariance, temperature and pressure cube. Additional cubes are passed through unchanged. -**kwargs : dict, optional + **kwargs : dict, optional Additional keyword arguments. -Returns -------- -Cube or CubeList + Returns + ------- + Cube or CubeList Input cubes with the pressure, temperature and covariance cubes removed and a new ``surface_upward_sensible_heat_flux`` cube added. Unrelated cubes are passed through unchanged. -Notes ------ -- Pressure is converted internally to Pa and temperature to K. -- Covariance units of ``degC m s-1`` are treated as numerically - equivalent to ``K m s-1`` because temperature offsets cancel when - forming fluctuations. -- Input cubes are assumed to be physically compatible; no regridding or - coordinate alignment is performed. -- Identification is unit-based, with metadata used only to resolve - ambiguities. - -Raises ------- -ValueError + Notes + ----- + - Pressure is converted internally to Pa and temperature to K. + - Covariance units of ``degC m s-1`` are treated as numerically + equivalent to ``K m s-1`` because temperature offsets cancel when + forming fluctuations. + - Input cubes are assumed to be physically compatible; no regridding or + coordinate alignment is performed. + - Identification is unit-based, with metadata used only to resolve + ambiguities. + + Raises + ------ + ValueError If suitable pressure, temperature or covariance cubes cannot be uniquely identified. -""" + """ from cf_units import Unit cubes = ( @@ -130,23 +130,29 @@ def sensible_heat_flux_from_covariance(cubes, **kwargs): # Pressure cube p_cand = [c for c in cubes if _is_p_cube(c)] if len(p_cand) > 1: - preferred = [ - c for c in p_cand - if c.standard_name == "air_pressure" - ] + preferred = [c for c in p_cand if c.name() == "barometric_pressure"] if len(preferred) == 1: p_cand = preferred + print("=== pressure candidates ===") + for c in p_cand: + print( + "name=", + c.name(), + "var_name=", + c.var_name, + "standard_name=", + c.standard_name, + "long_name=", + c.long_name, + ) pressure = _exactly_one(p_cand, "pressure") # Temperature cube T_cand = [c for c in cubes if _is_T_cube(c)] if len(T_cand) > 1: - preferred = [ - c for c in T_cand - if c.standard_name == "air_temperature" - ] + preferred = [c for c in T_cand if c.standard_name == "air_temperature"] if len(preferred) == 1: T_cand = preferred @@ -201,10 +207,7 @@ def sensible_heat_flux_from_covariance(cubes, **kwargs): used_ids = {id(wT), id(temp), id(pressure)} - out = iris.cube.CubeList( - c for c in cubes - if id(c) not in used_ids - ) + out = iris.cube.CubeList(c for c in cubes if id(c) not in used_ids) out.append(shf) diff --git a/tests/operators/test_fluxes.py b/tests/operators/test_fluxes.py index ad8c05fdf..9ef1d5b12 100644 --- a/tests/operators/test_fluxes.py +++ b/tests/operators/test_fluxes.py @@ -48,6 +48,29 @@ def _make_scalar_cube( return cube +def test_sensible_heat_units_accepts_surface_pressure(): + """Verify fallback path for pressure name.""" + wT = _make_scalar_cube( + 0.1, + "wt_covariance_2m", + units=Unit("K m s-1"), + ) + temp = _make_scalar_cube( + 20.0, + "air_temperature_rtd_1p2m", + units=Unit("degC"), + standard_name="air_temperature", + ) + press = _make_scalar_cube( + 1000.0, + "surface_pressure", + units=Unit("hPa"), + ) + out = fluxes.sensible_heat_flux_from_covariance([wT, temp, press]) + assert isinstance(out, iris.cube.Cube) + assert out.var_name == "surface_upward_sensible_heat_flux" + + def test_sensible_heat_units_missing_required_inputs(): """Raise if one of the required physical inputs cannot be identified.""" wT = _make_scalar_cube( From 3690eee0b14d3236ce078440aa73abe3a939c167 Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Wed, 8 Jul 2026 11:35:26 +0100 Subject: [PATCH 16/21] Refactor sensible_heat_flux_from_covariance parameters Refactor sensible_heat_flux_from_covariance function to use explicit parameters for wt_flux, air_temperature, and pressure. Update documentation and improve unit conversion handling. --- src/CSET/operators/fluxes.py | 219 +++++++---------------------------- 1 file changed, 45 insertions(+), 174 deletions(-) diff --git a/src/CSET/operators/fluxes.py b/src/CSET/operators/fluxes.py index 2e2eb43de..e7b4a3401 100644 --- a/src/CSET/operators/fluxes.py +++ b/src/CSET/operators/fluxes.py @@ -23,195 +23,66 @@ from CSET.operators._atmospheric_constants import CPD, LV, RD -def _exactly_one(matches, role): - if len(matches) == 0: - raise ValueError(f"sensible_heat_units could not identify a unique {role} cube") - if len(matches) > 1: - names = [getattr(c, "var_name", None) or c.name() for c in matches] - raise ValueError( - f"sensible_heat_units found multiple possible {role} cubes: {names}" - ) - return matches[0] - - -def _is_p_cube(cube): - return ( - cube.units is not None - and not cube.units.is_unknown() - and cube.units.is_convertible(Unit("Pa")) - ) - - -def _is_T_cube(cube): - if cube.units is None or cube.units.is_unknown(): - return False - return cube.units.is_convertible(Unit("K")) or cube.units.is_convertible( - Unit("degC") - ) - - -def _is_wt_covar_cube(cube): - if cube.units is None or cube.units.is_unknown(): - return False - # turbulence covariance may be recorded either as K m s-1 - # or degC m s-1; these are equivalent - return cube.units.is_convertible(Unit("K m s-1")) or cube.units.is_convertible( - Unit("degC m s-1") - ) - - -def sensible_heat_flux_from_covariance(cubes, **kwargs): +def sensible_heat_flux_from_covariance( + wt_flux: iris.cube.Cube | iris.cube.CubeList, + air_temperature: iris.cube.Cube | iris.cube.CubeList, + pressure: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.CubeList: """ Convert turbulent temperature covariance into sensible heat flux. - This operator computes surface upward sensible heat flux (SHF) from - temperature covariance using: - - SHF = ρ * CPD * (w'T') - - where air density is calculated from pressure and temperature via the - ideal gas law. + Computes: + SHF = rho * CPD * w'T' - The required input cubes are identified primarily from their physical - units: - - - temperature covariance (e.g. K m s-1 or degC m s-1) - - air temperature (convertible to K or degC) - - air pressure (convertible to Pa) - - If multiple physically plausible candidates are found, CF metadata - (e.g. standard names) are used as a secondary disambiguation step. - A ValueError is raised if the required cubes cannot be uniquely - identified. + where: + rho = pressure / (RD * temperature) Parameters ---------- - cubes : Cube or CubeList - Input cube(s) containing exactly one identifiable covariance, - temperature and pressure cube. Additional cubes are passed through - unchanged. + wt_flux : iris.cube.Cube or iris.cube.CubeList + Turbulent temperature covariance (w'T'). + + air_temperature : iris.cube.Cube or iris.cube.CubeList + Air temperature. - **kwargs : dict, optional - Additional keyword arguments. + pressure : iris.cube.Cube or iris.cube.CubeList + Air pressure. Returns ------- - Cube or CubeList - Input cubes with the pressure, temperature and covariance cubes - removed and a new - ``surface_upward_sensible_heat_flux`` cube added. Unrelated cubes - are passed through unchanged. - - Notes - ----- - - Pressure is converted internally to Pa and temperature to K. - - Covariance units of ``degC m s-1`` are treated as numerically - equivalent to ``K m s-1`` because temperature offsets cancel when - forming fluctuations. - - Input cubes are assumed to be physically compatible; no regridding or - coordinate alignment is performed. - - Identification is unit-based, with metadata used only to resolve - ambiguities. - - Raises - ------ - ValueError - If suitable pressure, temperature or covariance cubes cannot be - uniquely identified. + iris.cube.CubeList + Surface upward sensible heat flux cube(s). """ from cf_units import Unit - cubes = ( - iris.cube.CubeList(cubes) - if not isinstance(cubes, iris.cube.CubeList) - else cubes - ) - - # Pressure cube - p_cand = [c for c in cubes if _is_p_cube(c)] - if len(p_cand) > 1: - preferred = [c for c in p_cand if c.name() == "barometric_pressure"] - - if len(preferred) == 1: - p_cand = preferred - - print("=== pressure candidates ===") - for c in p_cand: - print( - "name=", - c.name(), - "var_name=", - c.var_name, - "standard_name=", - c.standard_name, - "long_name=", - c.long_name, - ) - pressure = _exactly_one(p_cand, "pressure") - - # Temperature cube - T_cand = [c for c in cubes if _is_T_cube(c)] - if len(T_cand) > 1: - preferred = [c for c in T_cand if c.standard_name == "air_temperature"] - - if len(preferred) == 1: - T_cand = preferred - - temp = _exactly_one(T_cand, "temperature") - - # Covariance cube - covar_cand = [c for c in cubes if _is_wt_covar_cube(c)] - if len(covar_cand) > 1: - preferred = [] - for cube in covar_cand: - text = " ".join( - str(x).lower() - for x in ( - cube.standard_name, - cube.var_name, - cube.long_name, - cube.name(), - ) - if x - ) - if "wt" in text or "w't" in text: - preferred.append(cube) - - if len(preferred) == 1: - covar_cand = preferred - - wT = _exactly_one(covar_cand, "w'T' covariance") - - # - # Unit conversions - # - temp_K = temp.copy() - if temp_K.units.is_convertible(Unit("degC")): + out = iris.cube.CubeList() + for wt_cube, temp_cube, pres_cube in zip( + iter_maybe(wt_flux), + iter_maybe(air_temperature), + iter_maybe(pressure), + strict=True, + ): + + # Unit conversions + temp_K = temp_cube.copy() temp_K.convert_units("K") - - pres_Pa = pressure.copy() - pres_Pa.convert_units("Pa") - - # Treat degC covariance numerically as K covariance - wT_cov = wT.copy() - if str(wT_cov.units) == "degC m s-1": - wT_cov.units = Unit("K m s-1") - - rho_air = pres_Pa.data / (RD * temp_K.data) - - shf = wT_cov.copy() - shf.data = CPD * rho_air * wT_cov.data - shf.units = Unit("W m-2") - shf.rename("surface_upward_sensible_heat_flux") - shf.var_name = "surface_upward_sensible_heat_flux" - - used_ids = {id(wT), id(temp), id(pressure)} - - out = iris.cube.CubeList(c for c in cubes if id(c) not in used_ids) - - out.append(shf) - - return out[0] if len(out) == 1 else out + pres_Pa = pres_cube.copy() + pres_Pa.convert_units("Pa") + wt_cov = wt_cube.copy() + + # Treat degC covariance numerically as K covariance + if str(wt_cov.units) == "degC m s-1": + wt_cov.units = Unit("K m s-1") + + # Density and SHF + rho_air = pres_Pa / (RD * temp_K) + shf = CPD * rho_air * wt_cov + + shf.units = Unit("W m-2") + shf.rename("surface_upward_sensible_heat_flux") + out.append(shf) + + return out def latent_heat_units( From 04d848c149ffa132e0e6c56b03b28c4edc87fe40 Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Wed, 8 Jul 2026 11:41:10 +0100 Subject: [PATCH 17/21] Update sensible heat flux tests for clarity Only two tests required now for a simplified operator --- tests/operators/test_fluxes.py | 158 +++++---------------------------- 1 file changed, 20 insertions(+), 138 deletions(-) diff --git a/tests/operators/test_fluxes.py b/tests/operators/test_fluxes.py index 9ef1d5b12..d096e41a6 100644 --- a/tests/operators/test_fluxes.py +++ b/tests/operators/test_fluxes.py @@ -48,8 +48,9 @@ def _make_scalar_cube( return cube -def test_sensible_heat_units_accepts_surface_pressure(): - """Verify fallback path for pressure name.""" +def test_sensible_heat_flux_core_calculation(): + """Compute sensible heat flux from covariance, temperature and pressure.""" + wT = _make_scalar_cube( 0.1, "wt_covariance_2m", @@ -61,57 +62,16 @@ def test_sensible_heat_units_accepts_surface_pressure(): units=Unit("degC"), standard_name="air_temperature", ) - press = _make_scalar_cube( + pressure = _make_scalar_cube( 1000.0, "surface_pressure", units=Unit("hPa"), ) - out = fluxes.sensible_heat_flux_from_covariance([wT, temp, press]) - assert isinstance(out, iris.cube.Cube) - assert out.var_name == "surface_upward_sensible_heat_flux" - - -def test_sensible_heat_units_missing_required_inputs(): - """Raise if one of the required physical inputs cannot be identified.""" - wT = _make_scalar_cube( - 0.1, - "wt_covariance_2m", - units=Unit("K m s-1"), - long_name="vertical wind-temperature covariance", + out = fluxes.sensible_heat_flux_from_covariance( + wt_flux=wT, + air_temperature=temp, + pressure=pressure, ) - temp = _make_scalar_cube( - 20.0, - "air_temperature_rtd_1p2m", - units=Unit("degC"), - standard_name="air_temperature", - ) - # pressure missing - with pytest.raises(ValueError, match="pressure"): - fluxes.sensible_heat_flux_from_covariance([wT, temp]) - - -def test_sensible_heat_units_core_calculation(): - """Compute sensible heat flux from covariance, temperature, and pressure.""" - wT = _make_scalar_cube( - 0.1, - "wt_covariance_2m", - units=Unit("K m s-1"), - long_name="vertical wind-temperature covariance", - ) - temp = _make_scalar_cube( - 20.0, - "air_temperature_rtd_1p2m", - units=Unit("degC"), - standard_name="air_temperature", - ) - press = _make_scalar_cube( - 1000.0, - "barometric_pressure", - units=Unit("hPa"), - long_name="barometric pressure", - ) - - out = fluxes.sensible_heat_flux_from_covariance([wT, temp, press]) arr = out.data if hasattr(arr, "compute"): arr = arr.compute() @@ -119,115 +79,37 @@ def test_sensible_heat_units_core_calculation(): T = 20.0 + 273.15 pPa = 1000.0 * 100.0 rho = pPa / (RD * T) + expected = CPD * rho * 0.1 assert np.isclose(arr[0, 0], expected) - assert out.var_name == "surface_upward_sensible_heat_flux" - assert out.units == Unit("W m-2") - - -def test_sensible_heat_units_returns_cube_when_only_shf_remains(): - """Return a single Cube if only the derived SHF cube remains.""" - wT = _make_scalar_cube( - 0.1, - "wt_covariance_2m", - units=Unit("K m s-1"), - long_name="vertical wind-temperature covariance", - ) - temp = _make_scalar_cube( - 20.0, - "air_temperature_rtd_1p2m", - units=Unit("degC"), - standard_name="air_temperature", - ) - press = _make_scalar_cube( - 1000.0, - "barometric_pressure", - units=Unit("hPa"), - long_name="barometric pressure", - ) - - out = fluxes.sensible_heat_flux_from_covariance([wT, temp, press]) - assert isinstance(out, iris.cube.Cube) - assert out.var_name == "surface_upward_sensible_heat_flux" + assert out.name() == "surface_upward_sensible_heat_flux" assert out.units == Unit("W m-2") -def test_sensible_heat_units_passthrough_and_filtering(): - """Remove only the cubes used in SHF calculation and retain unrelated cubes.""" - wT = _make_scalar_cube( - 0.1, - "wt_covariance_2m", - units=Unit("K m s-1"), - long_name="vertical wind-temperature covariance", - ) - temp = _make_scalar_cube( - 20.0, - "air_temperature_rtd_1p2m", - units=Unit("degC"), - standard_name="air_temperature", - ) - press = _make_scalar_cube( - 1000.0, - "barometric_pressure", - units=Unit("hPa"), - long_name="barometric pressure", - ) - extra = _make_scalar_cube( - 5.0, - "other_var", - units=Unit("1"), - ) - - out = fluxes.sensible_heat_flux_from_covariance([wT, temp, press, extra]) - assert isinstance(out, iris.cube.CubeList) - varnames = [c.var_name for c in out] - # input variables removed - assert "wt_covariance_2m" not in varnames - assert "air_temperature_rtd_1p2m" not in varnames - assert "barometric_pressure" not in varnames - # unrelated cube retained - assert "other_var" in varnames - # SHF added - assert "surface_upward_sensible_heat_flux" in varnames - +def test_sensible_heat_flux_accepts_degc_covariance(): + """degC m s-1 covariance should be treated as K m s-1.""" -def test_sensible_heat_units_prefers_barometric_pressure(): - """Prefer the strongest-scoring pressure cube.""" wT = _make_scalar_cube( 0.1, "wt_covariance_2m", - units=Unit("K m s-1"), + units=Unit("degC m s-1"), ) temp = _make_scalar_cube( 20.0, "air_temperature_rtd_1p2m", units=Unit("degC"), - standard_name="air_temperature", ) - press1 = _make_scalar_cube( + pressure = _make_scalar_cube( 1000.0, - "barometric_pressure", - units=Unit("hPa"), - ) - press2 = _make_scalar_cube( - 1002.0, "surface_pressure", units=Unit("hPa"), ) - out = fluxes.sensible_heat_flux_from_covariance([wT, temp, press1, press2]) - - # extract SHF cube - shf = next(c for c in out if c.var_name == "surface_upward_sensible_heat_flux") - arr = shf.data - if hasattr(arr, "compute"): - arr = arr.compute() - - T = 20.0 + 273.15 - pPa = 1000.0 * 100.0 # should use barometric_pressure - rho = pPa / (RD * T) - expected = CPD * rho * 0.1 - - assert np.isclose(arr[0, 0], expected) + out = fluxes.sensible_heat_flux_from_covariance( + wt_flux=wT, + air_temperature=temp, + pressure=pressure, + ) + assert out.units == Unit("W m-2") def _make_scalar_latent_cube( From a954be95a1a9ff48086985f760a909a92df66b44 Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Wed, 8 Jul 2026 11:43:15 +0100 Subject: [PATCH 18/21] Fix typo in test_sensible_heat_flux_accepts_degc_covariance --- tests/operators/test_fluxes.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/operators/test_fluxes.py b/tests/operators/test_fluxes.py index d096e41a6..51da197dc 100644 --- a/tests/operators/test_fluxes.py +++ b/tests/operators/test_fluxes.py @@ -87,7 +87,7 @@ def test_sensible_heat_flux_core_calculation(): def test_sensible_heat_flux_accepts_degc_covariance(): - """degC m s-1 covariance should be treated as K m s-1.""" + """DegC m s-1 covariance should be treated as K m s-1.""" wT = _make_scalar_cube( 0.1, From a3a7dc493534ad2eb126452dfb080486560b9e28 Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Wed, 8 Jul 2026 11:52:55 +0100 Subject: [PATCH 19/21] ruff format fixes --- src/CSET/operators/fluxes.py | 3 +-- tests/operators/test_fluxes.py | 3 --- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/src/CSET/operators/fluxes.py b/src/CSET/operators/fluxes.py index e7b4a3401..13c105340 100644 --- a/src/CSET/operators/fluxes.py +++ b/src/CSET/operators/fluxes.py @@ -62,7 +62,6 @@ def sensible_heat_flux_from_covariance( iter_maybe(pressure), strict=True, ): - # Unit conversions temp_K = temp_cube.copy() temp_K.convert_units("K") @@ -73,7 +72,7 @@ def sensible_heat_flux_from_covariance( # Treat degC covariance numerically as K covariance if str(wt_cov.units) == "degC m s-1": wt_cov.units = Unit("K m s-1") - + # Density and SHF rho_air = pres_Pa / (RD * temp_K) shf = CPD * rho_air * wt_cov diff --git a/tests/operators/test_fluxes.py b/tests/operators/test_fluxes.py index 51da197dc..f7f81880a 100644 --- a/tests/operators/test_fluxes.py +++ b/tests/operators/test_fluxes.py @@ -18,7 +18,6 @@ import iris.coords import iris.cube import numpy as np -import pytest from cf_units import Unit from CSET.operators import fluxes @@ -50,7 +49,6 @@ def _make_scalar_cube( def test_sensible_heat_flux_core_calculation(): """Compute sensible heat flux from covariance, temperature and pressure.""" - wT = _make_scalar_cube( 0.1, "wt_covariance_2m", @@ -88,7 +86,6 @@ def test_sensible_heat_flux_core_calculation(): def test_sensible_heat_flux_accepts_degc_covariance(): """DegC m s-1 covariance should be treated as K m s-1.""" - wT = _make_scalar_cube( 0.1, "wt_covariance_2m", From 4412b9c0c1638e0ff91976b044c20ee83eb4c6ae Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Wed, 8 Jul 2026 12:07:25 +0100 Subject: [PATCH 20/21] Rename variables for consistency in tests --- tests/operators/test_fluxes.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/operators/test_fluxes.py b/tests/operators/test_fluxes.py index f7f81880a..ae19c73f1 100644 --- a/tests/operators/test_fluxes.py +++ b/tests/operators/test_fluxes.py @@ -51,18 +51,18 @@ def test_sensible_heat_flux_core_calculation(): """Compute sensible heat flux from covariance, temperature and pressure.""" wT = _make_scalar_cube( 0.1, - "wt_covariance_2m", + "wt_covariance", units=Unit("K m s-1"), ) temp = _make_scalar_cube( 20.0, - "air_temperature_rtd_1p2m", + "air_temperature", units=Unit("degC"), standard_name="air_temperature", ) pressure = _make_scalar_cube( 1000.0, - "surface_pressure", + "surface_air_pressure", units=Unit("hPa"), ) out = fluxes.sensible_heat_flux_from_covariance( From 896096a5087d08d165fedf23febc3638899e3470 Mon Sep 17 00:00:00 2001 From: Simon Osborne Date: Wed, 8 Jul 2026 13:09:51 +0100 Subject: [PATCH 21/21] return cube rather than cubelist from sensible heat flux function --- src/CSET/operators/fluxes.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/CSET/operators/fluxes.py b/src/CSET/operators/fluxes.py index 13c105340..65fcd5ccc 100644 --- a/src/CSET/operators/fluxes.py +++ b/src/CSET/operators/fluxes.py @@ -79,9 +79,12 @@ def sensible_heat_flux_from_covariance( shf.units = Unit("W m-2") shf.rename("surface_upward_sensible_heat_flux") + # Keep compatibility with existing tests/output + shf.var_name = "surface_upward_sensible_heat_flux" + out.append(shf) - return out + return out[0] if len(out) == 1 else out def latent_heat_units(