Skip to content

Adds windspeed below the MAUL #2140

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

Open
David Flack (daflack) wants to merge 9 commits into
base: wind_speed_direction_calculator
Choose a base branch
from
Open

Adds windspeed below the MAUL #2140

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
0befe8e
Adds windspeed below the MAUL
daflack May 18, 2026
3231dec
Adds initial implementation of wind_below_maul
daflack May 18, 2026
2ae246f
Adds additional data for tests
daflack May 18, 2026
993cf74
Corrects test data
daflack May 18, 2026
0f670b5
Update looping to use itertools.product
daflack May 18, 2026
d433322
Fixes tests
daflack May 18, 2026
7739590
Add tests, fixes operator, adds test data
daflack May 19, 2026
9308b86
Update recipes and add new recipe for wind below MAUL
daflack May 19, 2026
b0a6337
Adds plumbing of recipe into the workflow
daflack May 19, 2026
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
20 changes: 17 additions & 3 deletions src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1823,7 +1823,7 @@ title= Number of moist absolutely unstable layers (MAULs) in a column
description=PROCESS-BASED DIAGNOSTIC.
Determines the number of moist absolutely unstable layers (MAUL) present
in a grid column regardless of depth.
Requires "air_temperature", "air_pressure", and "specific_humidity" on model levels.
Requires "air_temperature", "air_pressure", "specific_humidity", "x_wind" and "y_wind" on model levels.
help=This diagnostic identifies the numbers of MAUL that are present within the column.
It does not provide any information about the properties of the MAUL. A MAUL
is defined as an unstable layer (equivalent potential temperature gradient
Expand All @@ -1840,7 +1840,7 @@ title= Depth of moist absolutely unstable layer (MAUL) in a column
description=PROCESS-BASED DIAGNOSTIC.
Determines the depth of the deepest moist absolutely unstable layer (MAUL)
in a grid column.
Requires "air_temperature", "air_pressure", and "specific_humidity" on model levels.
Requires "air_temperature", "air_pressure", "specific_humidity", "x_wind" and "y_wind" on model levels.
help=This diagnostic identifies only the deepest MAUL within the column. A MAUL
is defined as an unstable layer (equivalent potential temperature gradient
with height less than zero) that has a relatively humidity above '90%' following
Expand All @@ -1856,7 +1856,7 @@ title= Base (AGL) of moist absolutely unstable layer (MAUL) in a column
description=PROCESS-BASED DIAGNOSTIC.
Determines the base height (AGL) of the deepest moist absolutely unstable layer (MAUL)
in a grid column.
Requires "air_temperature", "air_pressure", and "specific_humidity" on model levels.
Requires "air_temperature", "air_pressure", "specific_humidity", "x_wind" and "y_wind" on model levels.
help=This diagnostic identifies the base of the deepest MAUL within the column. A MAUL
is defined as an unstable layer (equivalent potential temperature gradient
with height less than zero) that has a relatively humidity above '90%' following
Expand All @@ -1866,6 +1866,20 @@ type=python_boolean
compulsory=true
sort-key=xp-maulp3

[tepmplate variable=AVERAGE_WIND_BELOW_MAUL]
ns = Diagnostics/Derived/XPPN
title= Average windspeed below the deepest moist absolutely unstable layer (MAUL) in a column
description=PROCESS-BASED DIAGNOSTIC.
Determines the average windspeed below the deepest moist absolutely unstable layer (MAUL)

@jwarner8 James Warner (jwarner8) May 28, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

worth being clear that this is the horizontal component of wind, not the true windspeed in the model (which would include a vertical component which is not used here).

in a grid column.
Requires "air_temperature", "air_pressure", "specific_humidity", "x_wind" and "y_wind" on model levels.
help=This diagnostic identifies the base of the deepest MAUL within the column and
averages the windspeed below that base. A MAUL is defined as an unstable layer
(equivalent potential temperature gradient with height less than zero) that has
a relatively humidity above '90%' following Takemi and Unuma (2020).
type=python_boolean
compulsory=true
sort-key=xp-maulp4
###################################
# Ensembles
[Diagnostics/Ensembles]
Expand Down
1 change: 1 addition & 0 deletions src/CSET/cset_workflow/rose-suite.conf.example
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@ ANALYSIS_LENGTH=""
!!AOA_CYCLIC=False
AOA_DIAG=False
!!AOA_PLEV=[]
AVERAGE_WIND_BELOW_MAUL=False
AVIATION_COLOUR_STATE=False
AVIATION_COLOUR_STATE_CLOUD_BASE=False
AVIATION_COLOUR_STATE_VISIBILITY=False
Expand Down
16 changes: 16 additions & 0 deletions src/CSET/loaders/spatial_field.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -651,6 +651,22 @@ def load(conf: Config):
aggregation=False,
)

# Wind Below Moist Absolutely Unstable Layer (of deepest)
if conf.AVERAGE_WIND_BELOW_MAUL:
for model in models:
yield RawRecipe(
recipe="average_windspeed_below_maul_spatial_plot.yaml",
variables={
"MODEL_NAME": model["name"],
"SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None,
"SUBAREA_EXTENT": conf.SUBAREA_EXTENT
if conf.SELECT_SUBAREA
else None,
},
model_ids=model["id"],
aggregation=False,
)

# Screen-level temperature probabilities
for model, condition, threshold in itertools.product(
models,
Expand Down
183 changes: 169 additions & 14 deletions src/CSET/operators/precipitation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,11 +22,14 @@
from skimage.measure import label

from CSET._common import iter_maybe
from CSET.operators.wind import calculate_vector_wind


def MAUL_properties(
cubes: iris.cube.Cube | iris.cube.CubeList,
output: Literal["number", "base", "depth"],
u_cubes: iris.cube.Cube | iris.cube.CubeList,
v_cubes: iris.cube.Cube | iris.cube.CubeList,
output: Literal["number", "base", "depth", "wind_below"],
) -> iris.cube.Cube | iris.cube.CubeList:
"""Identify properties of Moist Absolutely Unstable Layers.

Expand All @@ -35,11 +38,16 @@ def MAUL_properties(
cubes: iris.cube.Cube | iris.cube.CubeList
A cube or cubelist of a mask(s) as to whether a MAUL exists.
This input must be a binary field.
output: Literal["number", "base", "depth"]
u_cubes: iris.cube.Cube | iris.cube.CubeList
A cube or cubelist of the wind in the u direction.
v_cubes: iris.cube.Cube | iris.cube.CubeList
A cube or cubelist of the wind in the v direction.
output: Literal["number", "base", "depth", "wind_below"]
The output is the desired property required. It can be
number, base, depth for the number of MAULs, base height
of the deepest MAUL, or the depth of the deepest MAUL,
respectively.
of the deepest MAUL, the depth of the deepest MAUL, or the
average windspeed below the MAUL, respectively.


Returns
-------
Expand All @@ -50,7 +58,7 @@ def MAUL_properties(
------
ValueError: Data contains values that are not 0 or 1, only masked data should be used.
This error is raised when a mask field is not provided to the operator.
ValueError: Unexpected value for output. Expected number, depth or base. Got {output}.
ValueError: Unexpected value for output. Expected number, base, depth or wind_below. Got {output}.
This error is raised when the wrong output string is specified.

Notes
Expand All @@ -60,9 +68,10 @@ def MAUL_properties(
set out in a recipe. The operator applies image processing to the mask
to each point in the latitude/longitude coordinates. It uses the image
processing to identify continuous layers (1s), and labels them.
It identifies the number of layesr by identifying the maximum label number,
and then finds the top and base of each layer. Depending on the output
desired it will output information for the deepest MAUL.
It identifies the number of layers by identifying the maximum label number,
and then finds the top and base of each layer. It will also find the average
windspeed below the MAUL for indications of presence of low-level jets.
Depending on the output desired it will output information for the deepest MAUL.

When a MAUL is not present the output will be set to NaN for depth and base.
If number of MAULs is the desired output it will be set to zero.
Expand All @@ -72,13 +81,16 @@ def MAUL_properties(
num_MAULs = iris.cube.CubeList([])
maul_d = iris.cube.CubeList([])
maul_b = iris.cube.CubeList([])
windspeed_below_MAUL = iris.cube.CubeList([])

if output not in ("number", "base", "depth"):
if output not in ("number", "base", "depth", "wind_below"):
raise ValueError(
f"""Unexpected value for output. Expected number, depth or base. Got {output}."""
f"""Unexpected value for output. Expected number, base, depth or wind_below. Got {output}."""
)

for cube in iter_maybe(cubes):
for cube, u, v in zip(
iter_maybe(cubes), iter_maybe(u_cubes), iter_maybe(v_cubes), strict=True
):
# Check for binary fields.
if not np.array_equal(cube.data, cube.data.astype(bool)):
raise ValueError(
Expand All @@ -90,6 +102,12 @@ def MAUL_properties(
number_of_MAULs.data[:] = 0.0
maul_depth = number_of_MAULs.copy()
maul_base = number_of_MAULs.copy()
wind_below_maul = number_of_MAULs.copy()
# Calculate windspeed and direction.
windspeed_and_direction = calculate_vector_wind(u, v)
# Select windspeed, hard coded as always in same position from output
# of calculate_vector_wind.
windspeed = windspeed_and_direction[0]
# Loop over realization.
for mem_number, member in enumerate(cube.slices_over("realization")):
# Loop over time.
Expand Down Expand Up @@ -129,7 +147,7 @@ def MAUL_properties(
)
else:
number_of_MAULs.data[lat_point, lon_point] = np.max(labels)
if output != "number":
if output not in ("number", "wind_below"):
# Find the base, top, and depth for each object
# using cube metadata.
maul_start = []
Expand Down Expand Up @@ -247,6 +265,135 @@ def MAUL_properties(
else:
maul_depth.data[lat_point, lon_point] = np.nan
maul_base.data[lat_point, lon_point] = np.nan
# Separate loop for calculating wind properties.
elif output not in ("number"):
# Find the base, top, and depth for each object
# using cube metadata.
maul_start = []
maul_end = []
maul_dep = []
# Loop over the number of MAULs (plus one to ensure
# the case for only one MAUL being present).
for maul in range(1, np.max(labels) + 1):

@jwarner8 James Warner (jwarner8) May 28, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

not sure I follow the range starting from one here, even with the comment

# Find all vertical indices belonging to a MAUL.
maul_range = np.where(labels == maul)
# Find the height at the base of the MAUL
# (lowest level).
maul_start_point = lon.coord("level_height").points[
maul_range[0][0]
]
# Find the height at the top of the MAUL
# (highest level).
maul_end_point = lon.coord("level_height").points[
maul_range[0][-1]
]
# Calculate the MAUL depth, and store
# base and top heights.
maul_dep.append(maul_end_point - maul_start_point)
maul_start.append(maul_start_point)
maul_end.append(maul_end_point)
try:
# Idendtify where the deepest MAUL is.

@jwarner8 James Warner (jwarner8) May 28, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# Idendtify where the deepest MAUL is.
# Identify where the deepest MAUL is.

index = int(

@jwarner8 James Warner (jwarner8) May 28, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

argmax? I.e. index = np.argmax(maul_dep)

np.where(maul_dep == np.max(maul_dep))[0][0]
)
maul_base_value = maul_start[index]
height_index = np.abs(
lon.coord("level_height").points - maul_base_value
).argmin()
# As with number the code checks for whether
# there are multiple realization and/or time
# points for correct indexing of the output data
# and applies accordingly.
if (
len(number_of_MAULs.coord("realization").points)
!= 1
and len(number_of_MAULs.coord("time").points) != 1
):
# Store and calculate the windspeed below the
# deepest MAUL.
wind_below_maul.data[
mem_number, time_point, lat_point, lon_point
] = np.mean(
windspeed[
mem_number,
time_point,
0:height_index,
lat_point,
lon_point,
].data
)
elif (
len(number_of_MAULs.coord("realization").points)
!= 1
and len(number_of_MAULs.coord("time").points) == 1
):
wind_below_maul.data[
mem_number, lat_point, lon_point
] = np.mean(
windspeed[
mem_number,
0:height_index,
lat_point,
lon_point,
].data
)
elif (
len(number_of_MAULs.coord("time").points) != 1
and len(number_of_MAULs.coord("realization").points)
== 1
):
wind_below_maul.data[
time_point, lat_point, lon_point
] = np.mean(
windspeed[
time_point,
0:height_index,
lat_point,
lon_point,
].data
)
else:
wind_below_maul.data[lat_point, lon_point] = (
np.mean(
windspeed[
0:height_index, lat_point, lon_point
].data
)
)
# Here a ValueError is raised if a MAUL is not found, or an
# IndexError if the MAUL starts at the surface and so there
# is no wind below the MAUL however these are a valid answers,
# and so output data is set to NaN.
# The dimensionality logic for output data is identical
# to that used previously.
except (ValueError, IndexError):
if (
len(number_of_MAULs.coord("realization").points)
!= 1
and len(number_of_MAULs.coord("time").points) != 1
):
wind_below_maul.data[
mem_number, time_point, lat_point, lon_point
] = np.nan
elif (
len(number_of_MAULs.coord("realization").points)
!= 1
and len(number_of_MAULs.coord("time").points) == 1
):
wind_below_maul.data[
mem_number, lat_point, lon_point
] = np.nan
elif (
len(number_of_MAULs.coord("time").points) != 1
and len(number_of_MAULs.coord("realization").points)
== 1
):
wind_below_maul.data[
time_point, lat_point, lon_point
] = np.nan
else:
wind_below_maul.data[lat_point, lon_point] = np.nan

# Units and renaming for number, depth and base (the other case).
match output:
Expand All @@ -258,10 +405,14 @@ def MAUL_properties(
maul_depth.units = "m"
maul_depth.rename("MAUL_depth")
maul_d.append(maul_depth)
case _:
case "base":
maul_base.units = "m"
maul_base.rename("MAUL_base_height")
maul_b.append(maul_base)
case _:
wind_below_maul.units = "m s^-1"
wind_below_maul.rename("windspeed_below_MAUL")
windspeed_below_MAUL.append(wind_below_maul)

# Output data.
match output:
Expand All @@ -275,5 +426,9 @@ def MAUL_properties(
return maul_d
case "base" if len(maul_b) == 1:
return maul_b[0]
case _:
case "base":
return maul_b
case "wind_below" if len(windspeed_below_MAUL) == 1:
return windspeed_below_MAUL[0]
case _:
return windspeed_below_MAUL
Loading