diff --git a/improver/api/__init__.py b/improver/api/__init__.py index 73688a91f8..b9ab29abaa 100644 --- a/improver/api/__init__.py +++ b/improver/api/__init__.py @@ -71,6 +71,7 @@ "HailSize": "improver.psychrometric_calculations.hail_size", "height_of_maximum": "improver.utilities.cube_manipulation", "HumidityMixingRatio": "improver.psychrometric_calculations.psychrometric_calculations", + "IcingSeverityMultivariateRegression_USAF2024": "improver.icing", "Integration": "improver.utilities.mathematical_operations", "InterpolateUsingDifference": "improver.utilities.interpolation", "LapseRate": "improver.temperature.lapse_rate", diff --git a/improver/cli/icing_severity_multivariate_regression_usaf2024.py b/improver/cli/icing_severity_multivariate_regression_usaf2024.py new file mode 100755 index 0000000000..66e6193b01 --- /dev/null +++ b/improver/cli/icing_severity_multivariate_regression_usaf2024.py @@ -0,0 +1,71 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# ----------------------------------------------------------------------------- +# (C) British Crown copyright. The Met Office. +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# * Redistributions of source code must retain the above copyright notice, this +# list of conditions and the following disclaimer. +# +# * Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# * Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from +# this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +"""Script to create aircraft icing severidy indicies from multi-parameter datasets.""" + +from improver import cli + + +@cli.clizefy +@cli.with_output +def process( + *cubes: cli.inputcube, model_id_attr: str = None, +): + """ + From the supplied following cubes: + Air Temperature (t in K), + Relative Humidity (rh in %), + calculate a Aircraft Icing Severity Index using regression equation. + + The cubes for T and RH must match spatial and temporal coordinates. + + Does not collapse a realization coordinate. + + Args: + cubes (list of iris.cube.Cube): + Cubes to be processed. + model_id_attr (str): + Name of the attribute used to identify the source model for + blending. + + Returns: + iris.cube.Cube: + Cube of Aircraft Icing Severity Index + """ + from iris.cube import CubeList + + from improver.icing import IcingSeverityMultivariateRegression_USAF2024 + + result = IcingSeverityMultivariateRegression_USAF2024()( + CubeList(cubes), model_id_attr=model_id_attr + ) + + return result diff --git a/improver/icing.py b/improver/icing.py new file mode 100755 index 0000000000..584bbfaaeb --- /dev/null +++ b/improver/icing.py @@ -0,0 +1,179 @@ +# -*- coding: utf-8 -*- +# ----------------------------------------------------------------------------- +# (C) British Crown copyright. The Met Office. +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# * Redistributions of source code must retain the above copyright notice, this +# list of conditions and the following disclaimer. +# +# * Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# * Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from +# this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +"""Module containing aviation icing classes.""" +from datetime import timedelta +from typing import Tuple + +import iris +import numpy as np +from iris.cube import Cube, CubeList + +from improver import PostProcessingPlugin +from improver.metadata.constants import FLOAT_DTYPE +from improver.metadata.utilities import ( + create_new_diagnostic_cube, + generate_mandatory_attributes, +) +from improver.utilities.cube_checker import spatial_coords_match + +class IcingSeverityMultivariateRegression_USAF2024(PostProcessingPlugin): + """ + The algorithm outputs the unitless aircraft icing severity index. + This index can be converted directly to categorical icing severity level + using the category definitions below. Alternatively, the probability of + reaching or exceeding these categorical icing severity levels can be + calculated in a downstream thresholding operation. + + Inputs: + Temperature (T) in units of K + Relative Humidity (RH) in units of % + + Outputs: + Aircraft icing severity index (AISI) unitless + + Description of the algorithm: + + IF RH is greater than 70% and T is between 250.0K and 273.15K THEN: + AISI=100*TANH(0.06*RH-4.0)[TANH(0.1(T-247.0)] + + Categorical icing severity levels are defined as + AISI < 58 : "No Icing" + 58 <= AISI < 85 : "Light Icing" + 85 <= AISI < 92 : "Moderate Icing" + 92 <= AISI : "Severe Icing" + + """ + + @staticmethod + def _extract_input(cubes: CubeList, cube_name: str) -> Cube: + """Extract the relevant cube based on the cube name. + + Args: + cubes: Cubes from which to extract required input. + cube_name: Name of cube to extract. + + Returns: + The extracted cube. + """ + try: + cube = cubes.extract_cube(iris.Constraint(cube_name)) + except iris.exceptions.ConstraintMismatchError: + raise ValueError(f"No cube named {cube_name} found in {cubes}") + return cube + + def _get_inputs(self, cubes: CubeList) -> Tuple[Cube, Cube]: + """ + Separates T and RH cubes and checks that the following match: + forecast_reference_time, spatial coords, and time. + """ + + output_cubes = iris.cube.CubeList() + input_names = { + "air_temperature": ["K"], + "relative_humidity": ["%"], + } + + for input_name, units in input_names.items(): + output_cubes.append(self._extract_input(cubes, input_name)) + if not output_cubes[-1].units in units: + expected_unit_string = " or ".join(map(str, units)) + received_unit_string = str(output_cubes[-1].units) + raise ValueError( + f"The {output_cubes[-1].name()} units are incorrect, expected " + f"units as {expected_unit_string} but received {received_unit_string})." + ) + + t,rh = output_cubes + + if t.coord("forecast_reference_time") != rh.coord("forecast_reference_time"): + raise ValueError( + f"{t.name()} and {rh.name()} do not have the same forecast reference time" + ) + + if not spatial_coords_match([t, rh]): + raise ValueError( + f"{t.name()} and {rh.name()} do not have the same spatial " + f"coordinates" + ) + + if t.coord("time") != rh.coord("time"): + raise ValueError( + f"{t.name()} and {rh.name()} do not have the same valid time" + ) + + return t, rh + + def process(self, cubes: CubeList, model_id_attr: str = None) -> Cube: + """ + From the supplied Air Temperature and Relative Humidity cubes, calculate the Aircraft + Icing Severity Index. + + Args: + cubes: + Cubes of Air Temperature and Relative Humidity. + model_id_attr: + The name of the dataset attribute to be used to identify the source + model when blending data from different models. + + Returns: + Cube of Aircraft Icing Severity Index + + Raises: + ValueError: + If one of the cubes is not found, doesn't match the others, or has incorrect units + """ + t, rh = self._get_inputs(cubes) + + # Regression equations require math on cubes with incompatible units, so strip data + template = t.copy() + t = t.data + rh = rh.data + + + # Regression equation if RH is greater than 70% and T is between 250.0K and 273.15K + aisi = 100*np.tanh(0.06*rh-4.0)*(np.tanh(0.1*(t-247.0))) + aisi[np.where(rh<70)] = 0 + aisi[np.where(t<250.0)] = 0 + aisi[np.where(t>273.15)] = 0 + + cube = create_new_diagnostic_cube( + name=( + "aircraft_icing_severity_index" + ), + units="1", + template_cube=template, + data=aisi.astype(FLOAT_DTYPE), + mandatory_attributes=generate_mandatory_attributes( + cubes, model_id_attr=model_id_attr + ), + ) + + return cube diff --git a/improver_tests/acceptance/SHA256SUMS b/improver_tests/acceptance/SHA256SUMS index 91e7aa1ac8..259fac262e 100644 --- a/improver_tests/acceptance/SHA256SUMS +++ b/improver_tests/acceptance/SHA256SUMS @@ -572,6 +572,9 @@ ce09aafc535c43611df16bb8b0bf4a43c49b2a8302cc34018b72d9373841508b ./hail-size/wi 7e315c8eae130125ce2eace27cac08b89233765f3f4fc55c6f1ce30b7da77a80 ./height-of-max-vertical-velocity/kgo.nc 90ac17c415ba2b0249de3f304bf2f511a9a294710e0577dac9231b6ab822660d ./height-of-max-vertical-velocity/max_vertical_velocity.nc 929f98fa947ca8b635d76f6d3509b368fe7780019af76172dddbae4fef21822d ./height-of-max-vertical-velocity/vertical_velocity_on_height_levels.nc +80cf3edc4bcbd78998313df1f0c007e5610828cee7e006c0cadabf1e0e670e57 ./icing-severity-multivariate-regression-usaf2024/kgo.nc +d0833dac2ef303ddfe33a525a255e71037283c1615a5fc1500dc0f270303db7d ./icing-severity-multivariate-regression-usaf2024/rh.nc +63687402be8df9ccbfdcf6cfd5a843066ef8f6ce3517aa1c2885ca4fb9663bc0 ./icing-severity-multivariate-regression-usaf2024/t.nc e4002b78026bf59b8a2274872dd87d17b4c6f54085ba75573f0b6099e3f62ae6 ./integrate-time-bounds/basic/kgo.nc edc20b73a66f29159ee676b98eae8eed9b8d5b2a1d7b7b906415d3e452cdb195 ./integrate-time-bounds/basic/kgo_renamed.nc 5aaa03199faf9df5fda699936b33df862b071b3790b04791fef31d8cc0fd074a ./integrate-time-bounds/basic/lightning_frequency.nc diff --git a/improver_tests/acceptance/test_icing_severity_multivariate_regression_usaf2024.py b/improver_tests/acceptance/test_icing_severity_multivariate_regression_usaf2024.py new file mode 100644 index 0000000000..d2ce1a0f4a --- /dev/null +++ b/improver_tests/acceptance/test_icing_severity_multivariate_regression_usaf2024.py @@ -0,0 +1,33 @@ +# (C) Crown Copyright, Met Office. All rights reserved. +# +# This file is part of 'IMPROVER' and is released under the BSD 3-Clause license. +# See LICENSE in the root of the repository for full licensing details. +"""Tests for the lightning_usaf""" + +import pytest + +from . import acceptance as acc + +#pytestmark = [pytest.mark.acc, acc.skip_if_kgo_missing] +CLI = acc.cli_name_with_dashes(__file__) +run_cli = acc.run_cli(CLI) + + +@pytest.mark.parametrize("with_model_attr", (True, False)) +def test_basic(tmp_path, with_model_attr): + """Test basic invocation""" + kgo_dir = acc.kgo_root() / "icing-severity-multivariate-regression-usaf2024" + kgo_path = kgo_dir / "kgo.nc" + output_path = tmp_path / "output.nc" + + args = [ + kgo_dir / "rh.nc", + kgo_dir / "t.nc", + "--output", + f"{output_path}", + ] + + + + run_cli(args) + acc.compare(output_path, kgo_path) diff --git a/improver_tests/icing/__init__.py b/improver_tests/icing/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/improver_tests/icing/test_IcingSeverityMultivariateRegression_USAF2024.py b/improver_tests/icing/test_IcingSeverityMultivariateRegression_USAF2024.py new file mode 100644 index 0000000000..c9f80e3c50 --- /dev/null +++ b/improver_tests/icing/test_IcingSeverityMultivariateRegression_USAF2024.py @@ -0,0 +1,185 @@ +# -*- coding: utf-8 -*- +# ----------------------------------------------------------------------------- +# (C) British Crown copyright. The Met Office. +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# * Redistributions of source code must retain the above copyright notice, this +# list of conditions and the following disclaimer. +# +# * Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# * Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from +# this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +"""Test methods in icing.IcingSeverityMultivariateRegression_USAF2024""" + +from datetime import datetime + +import numpy as np +import pytest +from iris.cube import Cube, CubeList + +from improver.icing import IcingSeverityMultivariateRegression_USAF2024 +from improver.metadata.constants.attributes import MANDATORY_ATTRIBUTE_DEFAULTS +from improver.synthetic_data.set_up_test_cubes import ( + add_coordinate, + set_up_variable_cube, +) +from improver.utilities.spatial import create_vicinity_coord + + +@pytest.fixture(name="temp_cube") +def temp_cube_fixture() -> Cube: + """ + Set up an air temperature cube for use in tests over a variety of conditions. + """ + + data = np.array([[270, 272], [266, 280]], dtype=np.float32) + cube = set_up_variable_cube( + data, + name="air_temperature", + units="K", + time=datetime(2017, 11, 10, 3, 0), + time_bounds=None, + attributes=MANDATORY_ATTRIBUTE_DEFAULTS, + standard_grid_metadata="gl_ens", + domain_corner=(-20, 0), + x_grid_spacing=20, + y_grid_spacing=20, + ) + return cube + + +@pytest.fixture(name="rh_cube") +def rh_cube_fixture() -> Cube: + """ + Set up a relative humidity cube for use in tests over a variety of conditions. + """ + + data = np.array([[65, 75], [90, 80]], dtype=np.float32) + cube = set_up_variable_cube( + data, + name="relative_humidity", + units="%", + time=datetime(2017, 11, 10, 3, 0), + time_bounds=None, + attributes=MANDATORY_ATTRIBUTE_DEFAULTS, + standard_grid_metadata="gl_ens", + domain_corner=(-20, 0), + x_grid_spacing=20, + y_grid_spacing=20, + ) + return cube + + + + +@pytest.fixture(name="expected_cube") +def expected_cube_fixture() -> Cube: + """ + Set up the Icing cube that we expect to get from the plugin. + """ + + data = np.array([[0, 45.593143], [84.660645, 0 ]], dtype=np.float32) + cube = set_up_variable_cube( + data, + name="aircraft_icing_severity_index", + units="1", + time=datetime(2017, 11, 10, 3, 0), + attributes=MANDATORY_ATTRIBUTE_DEFAULTS, + domain_corner=(-20, 0), + x_grid_spacing=20, + y_grid_spacing=20, + ) + + + return cube + + +def test_basic(temp_cube, rh_cube, expected_cube): + """Run the plugin and check the result cube matches the expected_cube""" + result = IcingSeverityMultivariateRegression_USAF2024()( + CubeList([temp_cube, rh_cube]), None + ) + assert result.xml().splitlines(keepends=True) == expected_cube.xml().splitlines( + keepends=True + ) + assert np.allclose(result.data, expected_cube.data) + + +def break_temp_name(temp_cube, rh_cube): + """Modifies temp_cube name to be changed to appear missing and + returns the error message this will trigger""" + temp_cube.rename("TEMP") + return r"No cube named air_temperature found in .*" + + +def break_rh_name(temp_cube, rh_cube): + """Modifies precip_cube name to be changed to appear missing and + returns the error message this will trigger""" + rh_cube.rename("RH") + return r"No cube named relative_humidity found in .*" + +def break_reference_time(temp_cube, rh_cube): + """Modifies temp_cube time points to be incremented by 1 second and + returns the error message this will trigger""" + temp_cube.coord("forecast_reference_time").points = ( + temp_cube.coord("forecast_reference_time").points + 1 + ) + return r".* and .* do not have the same forecast reference time" + + +def break_units(temp_cube, rh_cube): + """Modifies rh_cube units to be incorrect as K and + returns the error message this will trigger""" + rh_cube.units = "K" + return r"The .* units are incorrect, expected units as .* but received .*" + + +def break_coordinates(temp_cube, rh_cube): + """Modifies the first latitude point on the temp_cube (adds one degree) + and returns the error message this will trigger""" + points = list(temp_cube.coord("latitude").points) + points[0] = points[0] + 1 + temp_cube.coord("latitude").points = points + return ".* and .* do not have the same spatial coordinates" + +@pytest.mark.parametrize( + "breaking_function", + ( + break_temp_name, + break_rh_name, + break_reference_time, + break_units, + break_coordinates, + ), +) +def test_exceptions( + temp_cube, rh_cube, breaking_function +): + """Tests that a suitable exception is raised when the cube meta-data does + not match what is expected""" + error_msg = breaking_function( + temp_cube, rh_cube + ) + with pytest.raises(ValueError, match=error_msg): + IcingSeverityMultivariateRegression_USAF2024()( + CubeList([temp_cube, rh_cube]) + )