Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

MOBT-494: Phase probability plugin without percentile generation #1903

Merged
merged 8 commits into from
Jun 5, 2023
Merged

MOBT-494: Phase probability plugin without percentile generation #1903

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
4bc3df1
Starting to rewrite tests.
bayliffe May 16, 2023
659c9fb
Unit tests updated.
bayliffe May 16, 2023
671f927
Tweaks to unit tests and new acceptance tests.
bayliffe May 17, 2023
dfb2d27
Further tests added. Comparator for snow changed to ge to ensure that…
bayliffe May 18, 2023
77ea6c3
Checksum update and style fixes.
bayliffe May 18, 2023
c528b5b
Add extra test.
bayliffe May 23, 2023
e590b26
Add value error to doc-string.
bayliffe May 24, 2023
9051f20
Review changes.
bayliffe Jun 5, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
45 changes: 30 additions & 15 deletions improver/cli/phase_probability.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,31 +35,46 @@

@cli.clizefy
@cli.with_output
def process(*cubes: cli.inputcube, radius: float = None):
def process(*cubes: cli.inputcube):
"""
Converts a phase-change-level cube into the
probability of a specific precipitation phase being found at the surface.
Converts a phase-change-level cube into the probability of a specific
precipitation phase being found at the surface or at a site's altitude.

If provided with a phase-change level cube without a percentile coordinate,
the phase-change level is compared directly to the orographic / site height
for each grid cell / site. A binary probability of the phase is returned for
each grid cell / site.

If the phase-change level cube has a percentile coordinate specific
percentiles will be used for the snow, rain or hail falling-levels. If a
snow-falling-level diagnostic is provided, the 80th percentile altitude will
be used. If a rain-falling-level or hail-falling-level diagnostic is provided,
the 20th percentile altitude will be used.

Args:
cubes (iris.cube.CubeList or list):
Contains cubes of the altitude of the phase-change level (this
can be snow->sleet, hail->rain or sleet->rain) and the altitude of the
orography. The name of the phase-change level cube must be
either "altitude_of_snow_falling_level", "altitude_of_rain_from_hail_falling_level or
"altitude_of_rain_falling_level". The name of the orography
cube must be "surface_altitude".
radius (float or None):
Neighbourhood radius from which the 20th or 80th percentile is found (m).
Restrictions in the neighbourhooding code limit the use of a
neighbourhood to data on equal area projections.
If set to None percentiles are not generated, and the input
altitudes are used directly as the phase discriminators.
can be snow->sleet, hail->rain or sleet->rain), and the altitude of the
orography at grid cells, or the altitude of sites at which the phase
probability should be returned.

The name of the phase-change level cube must be either:

- "altitude_of_snow_falling_level"
- "altitude_of_rain_from_hail_falling_level"
- "altitude_of_rain_falling_level"

If the phase-change level cube contains percentiles, these must include
the 80th percentile for the snow-falling-level, and the 20th percentile
for any other phase.

The name of the orography cube must be "surface_altitude".
The name of the site ancillary most be "grid_neighbours".
"""
from iris.cube import CubeList

from improver.psychrometric_calculations.precip_phase_probability import (
PrecipPhaseProbability,
)

return PrecipPhaseProbability(radius=radius)(CubeList(cubes))
return PrecipPhaseProbability()(CubeList(cubes))
185 changes: 90 additions & 95 deletions improver/psychrometric_calculations/precip_phase_probability.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@
"""Module for calculating the probability of specific precipitation phases."""

import operator
from typing import List, Optional, Union
from typing import List, Union

import iris
import numpy as np
Expand All @@ -43,7 +43,6 @@
create_new_diagnostic_cube,
generate_mandatory_attributes,
)
from improver.nbhood.nbhood import GeneratePercentilesFromANeighbourhood
from improver.utilities.cube_checker import spatial_coords_match


Expand All @@ -52,56 +51,55 @@ class PrecipPhaseProbability(BasePlugin):
This plugin converts a falling-phase-change-level cube into the
probability of a specific precipitation phase being found at the surface.

For snow; the 80th percentile is taken from a neighbourhood around
each point and is compared with the orography. Where the orography is
higher, the returned probability-of-snow is 1, else 0.
For hail; the 20th percentile is also taken from a neighbourhood around
each point and is compared with the orography. Where the orography is
higher, the returned probability-of-rain-from-hail is 1, else 0.
For rain, the above method is modified to get the 20th percentile
and where the orography is lower than the percentile value, the returned
probability-of-rain is 1, else 0.
Consider the case in which a snow-falling-level diagnostic is given as
input. For a deterministic snow-falling-level diagnostic (i.e. no percentile
coordinate), these falling levels are compared directly to the altitudes
of interest to determine the binary probability (0 or 1) of snow at these
altitudes. A 0 is returned if the altitude is below the snow-falling-level
and a 1 if the altitude if above the snow-falling-level. If a probabilistic
snow-falling-level is provided, this plugin will seek to use the 80th
percentile falling-level to compare with the altitudes. This provides a
greater level of certainty that the phase is snow, and in conjunction with
the percentile chosen for rain (see below), provides for a greater depth
of sleet (mixed phase precipitation).

The plugin behaves in the same way for a rain-falling-level, but in the
case of a probabalistic input, the 20th percentile is used. This lowers
the falling-level within the atmosphere, providing for greater certainty
that the phase is entirely rain.

The altitudes used in these comparisons may be provided as gridded
orographies or as spot forecast ancillaries that contain an altitude
coordinate. In the former case the phase probability will be determined at
each grid point at the orographic altitude. In the latter case, the phase
probability will be determined at each site's specific altitude.
"""

def __init__(self, radius: Optional[float] = None) -> None:
"""
Initialise plugin

Args:
radius:
Neighbourhood radius from which the 20th or 80th percentile is found (m).
Restrictions in the neighbourhooding code limit the use of a
neighbourhood to data on equal area projections.
If set to None percentiles are not generated, and the input
altitudes are used directly as the phase discriminators.
"""
self.percentile_plugin = GeneratePercentilesFromANeighbourhood
self.radius = radius

def _extract_input_cubes(self, cubes: Union[CubeList, List[Cube]]) -> None:
"""
Separates the input list into the required cubes for this plugin,
detects whether snow, rain from hail or rain are required from the input phase-level
cube name and appropriately initialises the percentile_plugin, sets
the appropriate comparator operator for comparing with orography and
the unique part of the output cube name.
detects whether snow, rain from hail or rain are required from the input
phase-level cube name, and as required, appropriately extracts the relevant
percentile.

Converts units of falling_level_cube to that of orography_cube if
necessary. Sets flag for snow, rain from hail or rain depending on name of
falling_level_cube.
Converts units of falling_level_cube to that of orography_cube / site
altitudes if necessary. Sets flag for snow, rain from hail or rain depending
on name of falling_level_cube.

Args:
cubes:
Contains cubes of the altitude of the phase-change level (this
can be snow->sleet, hail->rain or sleet->rain) and the altitude of the
orography. The name of the phase-change level cube must be
orography or sites. The name of the phase-change level cube must be
"altitude_of_snow_falling_level", "altitude_of_rain_from_hail_falling_level" or
"altitude_of_rain_falling_level". The name of the orography
cube must be "surface_altitude".
"altitude_of_rain_falling_level". If a gridded orography is provided it must
be named "surface_altitude". If a spot forecast ancillary is provided it
must be named "grid_neighbours".

Raises:
ValueError: If cubes with the expected names cannot be extracted.
ValueError: If cubes does not have the expected length of 2.
ValueError: If a percentile cube does not contain the expected percentiles.
ValueError: If the extracted cubes do not have matching spatial
coordinates.
"""
Expand All @@ -114,92 +112,89 @@ def _extract_input_cubes(self, cubes: Union[CubeList, List[Cube]]) -> None:
raise ValueError(
"Spatial coords mismatch between " f"{cubes[0]} and " f"{cubes[1]}"
)
extracted_cube = cubes.extract("altitude_of_snow_falling_level")
if extracted_cube:
(self.falling_level_cube,) = extracted_cube
self.param = "snow"
self.comparator = operator.gt
if self.radius is not None:
self.get_discriminating_percentile = self.percentile_plugin(
self.radius, percentiles=[80.0]
)
elif cubes.extract("altitude_of_rain_falling_level"):
extracted_cube = cubes.extract("altitude_of_rain_falling_level")
(self.falling_level_cube,) = extracted_cube
self.param = "rain"
self.comparator = operator.lt
# We want rain that has come from sleet at or above the surface, so inverse of 80th
# centile is the 20th centile.
if self.radius is not None:
self.get_discriminating_percentile = self.percentile_plugin(
self.radius, percentiles=[20.0]
)
else:
extracted_cube = cubes.extract("altitude_of_rain_from_hail_falling_level")
if not extracted_cube:

definitions = {
"snow": {"comparator": operator.ge, "percentile": 80},
"rain": {"comparator": operator.lt, "percentile": 20},
"rain_from_hail": {"comparator": operator.lt, "percentile": 20},
}

for diagnostic, definition in definitions.items():
extracted_cube = cubes.extract(f"altitude_of_{diagnostic}_falling_level")
if extracted_cube:
# Once a falling-level cube has been extracted, exit this loop.
break
if not extracted_cube:
raise ValueError(
"Could not extract a rain, rain from hail or snow falling-level "
f"cube from {', '.join([cube.name() for cube in cubes])}"
)
(self.falling_level_cube,) = extracted_cube
self.param = diagnostic
self.comparator = definition["comparator"]
if self.falling_level_cube.coords("percentile"):
constraint = iris.Constraint(percentile=definition["percentile"])
required_percentile = self.falling_level_cube.extract(constraint)
if not required_percentile:
raise ValueError(
"Could not extract a rain, rain from hail or snow falling-level "
f"cube from {cubes}"
)
(self.falling_level_cube,) = extracted_cube
self.param = "rain_from_hail"
self.comparator = operator.lt
if self.radius is not None:
self.get_discriminating_percentile = self.percentile_plugin(
self.radius, percentiles=[20.0]
f"Cube {self.falling_level_cube.name()} does not "
"contain the required percentile "
f"{definition['percentile']}."
)
self.falling_level_cube = required_percentile
self.falling_level_cube.remove_coord("percentile")

orography_name = "surface_altitude"
extracted_cube = cubes.extract(orography_name)
if extracted_cube:
(self.orography_cube,) = extracted_cube
self.altitudes = extracted_cube[0].data
altitude_units = extracted_cube[0].units
elif cubes.extract("grid_neighbours"):
(extracted_cube,) = cubes.extract("grid_neighbours")
self.altitudes = extracted_cube.coord("altitude").points
altitude_units = extracted_cube.coord("altitude").units
else:
raise ValueError(
f"Could not extract {orography_name} cube from " f"{cubes}"
f"Could not extract {orography_name} cube from "
f"cube from {', '.join([cube.name() for cube in cubes])}"
)

if self.falling_level_cube.units != self.orography_cube.units:
if self.falling_level_cube.units != altitude_units:
self.falling_level_cube = self.falling_level_cube.copy()
self.falling_level_cube.convert_units(self.orography_cube.units)
self.falling_level_cube.convert_units(altitude_units)

def process(self, cubes: Union[CubeList, List[Cube]]) -> Cube:
"""
Derives the probability of a precipitation phase at the surface. If
the snow-sleet falling-level is supplied, this is the probability of
snow at (or below) the surface. If the sleet-rain falling-level is
supplied, this is the probability of rain at (or above) the surface.
If the hail-rain falling-level is supplied, this is the probability
of rain from hail at (or above) the surface.
Derives the probability of a precipitation phase at the surface /
site altitude. If the snow-sleet falling-level is supplied, this is
the probability of snow at the surface / site altitude. If the sleet-rain
falling-level is supplied, this is the probability of rain at the surface
/ site altitude. If the hail-rain falling-level is supplied, this is the
probability of rain from hail at the surface / site altitude.

Args:
cubes:
Contains cubes of the altitude of the phase-change level (this
can be snow->sleet, hail->rain or sleet->rain) and the altitude
of the orography.
of the orography or a spot forecast ancillary that defines site
altitudes.

Returns:
Cube containing the probability of a specific precipitation phase
reaching the surface orography. If the falling_level_cube was
snow->sleet, then this will be the probability of snow at the
surface. If the falling_level_cube was sleet->rain, then this
will be the probability of rain from sleet at the surface.
If the falling_level_cube was hail->rain, then this
will be the probability of rain from hail at the surface.
The probabilities are categorical (1 or 0) allowing
reaching the surface orography or site altitude. If the
falling_level_cube was snow->sleet, then this will be the probability
of snow at the surface. If the falling_level_cube was sleet->rain,
then this will be the probability of rain from sleet at the surface
or site altitude. If the falling_level_cube was hail->rain, then this
will be the probability of rain from hail at the surface or site
altitude. The probabilities are categorical (1 or 0) allowing
precipitation to be divided uniquely between snow, sleet and
rain phases.
"""
self._extract_input_cubes(cubes)
if self.radius is not None:
processed_falling_level = iris.util.squeeze(
self.get_discriminating_percentile(self.falling_level_cube)
)
else:
processed_falling_level = self.falling_level_cube

result_data = np.where(
self.comparator(self.orography_cube.data, processed_falling_level.data),
1,
0,
self.comparator(self.altitudes, self.falling_level_cube.data), 1, 0,
).astype(np.int8)
mandatory_attributes = generate_mandatory_attributes([self.falling_level_cube])

Expand Down
Loading