Introduction of Icing Multivariate Regression

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",

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

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

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

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)