From 0b52c83ab3a7af8b67b4b3087fea4d0243ed7eee Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 4 Jun 2026 09:34:52 +0100 Subject: [PATCH 001/105] Recipe for histograms of masked Nimrod accumulations #2113 --- .../recipes/surface_fields/radar_dev3.yaml | 49 +++++++++++++++++++ 1 file changed, 49 insertions(+) create mode 100644 src/CSET/recipes/surface_fields/radar_dev3.yaml diff --git a/src/CSET/recipes/surface_fields/radar_dev3.yaml b/src/CSET/recipes/surface_fields/radar_dev3.yaml new file mode 100644 index 0000000000..82683d7e36 --- /dev/null +++ b/src/CSET/recipes/surface_fields/radar_dev3.yaml @@ -0,0 +1,49 @@ +category: Histogram +title: "Histogram masked Nimrod $VARNAME\n$RADAR_NAME" +description: Histograms of Nimrod radar observations +# Extracts and plots the probability density of surface `$VARNAME`. It uses +# [`plt.hist`](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.hist.html). +# +# The default method generates a probability density so that the area under the histogram +# normalized to 1. +# Histograms of rainfall or snowfall rate are plotted using a logarithmic scale. +# Histograms of rainfall or snowfall amount are based on the +# [Klingaman et al. 2017](https://gmd.copernicus.org/articles/10/57/2017/gmd-10-57-2017.html) +# ASoP method, where histograms show the fractional contributions from each precipitation bin +# to the total precipitation. The area under the histogram shows the total precipitation. + +steps: + - operator: read.read_cubes + file_paths: $INPUT_PATHS +# model_names: [$RADAR_NAME, $WEIGHTS_NAME] + model_names: $ALL_NAME + subarea_type: $SUBAREA_TYPE + subarea_extent: $SUBAREA_EXTENT + constraint: + operator: constraints.combine_constraints +# variable_constraint: +# operator: constraints.generate_var_constraint +# varname: "Hourly rain accumulation" +# variable_constraint: +# operator: constraints.generate_var_constraint +# varname: "Hourly wts accumulation" + cell_methods_constraint: + operator: constraints.generate_cell_methods_constraint + cell_methods: [] + varname: $VARNAME + pressure_level_constraint: + operator: constraints.generate_level_constraint + coordinate: pressure + levels: [] +# constraint: ["Hourly rain accumulation", "Hourly wts accumulation"] + + + - operator: radar_filter.mask_by_weights + model_names: $RADAR_NAME + weights_names: $WEIGHTS_NAME + + - operator: write.write_cube_to_nc + overwrite: True + + - operator: plot.plot_histogram_series + sequence_coordinate: $SEQUENCE From 391c43305bc959624163ff20a492e402ce4c39e8 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 4 Jun 2026 09:44:01 +0100 Subject: [PATCH 002/105] Loader for Nimrod masked histograms #2113 --- src/CSET/loaders/radar.py | 41 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index ae74148a0a..5f53462d6f 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -95,6 +95,13 @@ def load(conf: Config): if radar["varname"] == "Hourly rain accumulation" ] + # Form the list of accumulated hourly weights for Nimrod radar sources. + wts_radars = [ + radar + for radar in radar_sources + if radar["varname"] == "Hourly wts accumulation" + ] + # Surface (2D) fields for Nimrod radar rainfall. # # The different sources of Nimrod rainfall accumulation have @@ -135,6 +142,40 @@ def load(conf: Config): "SEQUENCE": "time" if conf.HISTOGRAM_SURFACE_FIELD_SEQUENCE else "realization", + "SUBAREA_NAME": "", + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, + }, + aggregation=False, + ) + + # Histograms for surface (2D) Nimrod radar hourly accumulated rainfall. + # + # The histograms are produced after the rainfall obs have been masked using + # the associated Nimrod weights file. + # + # To get multiple radar sources plotted on the histogram the + # recipe must be done by passing lists of the radar_ids and + # the radar_names. As this is a multiline plot, all radar sources + # share the same radar variable name. + radar_obs_ids = [radar["id"] for radar in accum_radars] + radar_wts_ids = [radar["id"] for radar in wts_radars] + combined_ids = radar_obs_ids + radar_wts_ids + print(" combined_ids: ", combined_ids) + if conf.HISTOGRAM_SURFACE_FIELD: + yield RawRecipe( + recipe="radar_dev3.yaml", + # model_ids -> Becomes $INPUT_PATHS + # model_ids=[ radar_obs_ids, radar_wts_ids], + model_ids=combined_ids, + variables={ + "VARNAME": next(radar["varname"] for radar in accum_radars), + "ALL_NAME": combined_ids, + "RADAR_NAME": [radar["id"] for radar in accum_radars], + "WEIGHTS_NAME": [radar["id"] for radar in wts_radars], + "SEQUENCE": "time" + if conf.HISTOGRAM_SURFACE_FIELD_SEQUENCE + else "realization", "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", From 8a091ece9f6837b3afdc73f9e7890e27f2bd321f Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 4 Jun 2026 09:59:33 +0100 Subject: [PATCH 003/105] Function to produce Nimrod masked cubes #2113 --- src/CSET/operators/__init__.py | 2 + src/CSET/operators/radar_filter.py | 196 +++++++++++++++++++++++++++++ 2 files changed, 198 insertions(+) create mode 100644 src/CSET/operators/radar_filter.py diff --git a/src/CSET/operators/__init__.py b/src/CSET/operators/__init__.py index 634506dd6a..94637911c7 100644 --- a/src/CSET/operators/__init__.py +++ b/src/CSET/operators/__init__.py @@ -44,6 +44,7 @@ power_spectrum, precipitation, pressure, + radar_filter, read, regrid, scoreswrappers, @@ -75,6 +76,7 @@ "power_spectrum", "precipitation", "pressure", + "radar_filter", "read", "regrid", "temperature", diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py new file mode 100644 index 0000000000..4350099f92 --- /dev/null +++ b/src/CSET/operators/radar_filter.py @@ -0,0 +1,196 @@ +# © Crown copyright, Met Office (2022-2026) 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 perform various kind of filtering.""" + +import iris +import iris.cube +import iris.exceptions + +from CSET._common import iter_maybe +from CSET.operators.filters import apply_mask, generate_mask + + +def mask_list(model_names: list[str]) -> list[str]: + """Determine the Nimrod weights files to use. + + Parameters + ---------- + model_names: list[str] + A list of model and Nimrod hourly rainfall accumulation files. + + Returns + ------- + list[str] + A list of the Nimrod weights files to use with each of the input + model / observations files. + + Notes + ----- + At lest one of the entries in the input list must be a Nimrod hourly + rainfall accumulation file. + + If just one Nimrod file is specified, then then the weights file associated + with this field is used. + + If more than one Nimrod file is in the input list, then each of the these + Nimrod files is associated with its own weights file e.g. if the input list + contains ["Nimrod1km", "Nimrod2km"] then the weights files for these will + be ["Nimrod1km_weights", "Nimrod2km_weights"]. Any model fields in the input + list will be allocated a weights file according to the order of preference + specified in the list nimrod_preference e.g. if the input list is + ["UM_model", "Nimrod1km", "Nimrod2km"] then the output weights files list will + be ["Nimrod2km_weights", "Nimrod1km_weights", "Nimrod2km_weights"] as the Nimrod + weights for 2km data are preferred over those for 1km. + + Examples + -------- + >>> list_weights = make_list( ["UM_model", "Nimrod1km", "Nimrod2km"] ) + >>> print(list_weights) + ["Nimrod2km_weights", "Nimrod1km_weights", "Nimrod2km_weights"] + + """ + # Set the preference order for choosing a Nimrod radar weights source + # in order of most to least preferred. + nimrod_preference = [ + "Nimrod2km", + "Nimrod_2km", + "Nimrodxkm", + "Nimrod_xkm", + "Nimrod1km", + "Nimrod_1km", + ] + + # Define the string that helps form a Nimrod weights file. + wei = "_weights" + + # Determine the preferred Nimrod mask to use. + empty_string = "" + preferred_nimrod = empty_string + for prefer in reversed(nimrod_preference): + if any(prefer in model for model in model_names): + preferred_nimrod = prefer + + # Create the list of the required Nimrod masks. + mask_names_list = [] + if preferred_nimrod != empty_string: + # Loop over the input model_names. + for model in model_names: + if any(model in nimrod for nimrod in nimrod_preference): + nimrod_mask = model + wei + else: + nimrod_mask = preferred_nimrod + wei + mask_names_list.append(nimrod_mask) + + return mask_names_list + + +def mask_by_weights( + cubes: iris.cube.CubeList, + model_names: list[str], + weights_names: list[str], + **kwargs, +) -> iris.cube.CubeList: + """Filter a field using a second field as a mask. + + Parameters + ---------- + cubes: iris.cube.CubeList + Two cubes containing the radar observations and their weights. + + Returns + ------- + Cube + + Raises + ------ + ValueError, iris.exceptions.NotYetImplementedError + When the cubes are not compatible. + + Notes + ----- + This is a simple operator designed for combination of diagnostics or + creating new diagnostics by using recipes. + + Examples + -------- + >>> field_filtered = mask_by_weights(cubelist, model_names) + + """ + print("model_names are: ", model_names) + print("weights_names", weights_names) + + for cube in cubes: + print(" cube.var_name ", cube.var_name) + print(" cube.name ", cube.name) + print(" cube: ") + print(cube) + print(" cube.attributes.model_name ", cube.attributes["model_name"]) + + # Check the input unfiltered cubes and the mask cubes are both cubelists + # with the same number of cubes. If not, then add extra mask cubes. + if len(model_names) != len(weights_names): + weights_names = mask_list(model_names) + + # Create an empty cubelist to hold the filtered fields. + filtered_list = iris.cube.CubeList([]) + + # Loop over the fields to filter. + var_constraint = iris.NameConstraint(var_name="hourly_rain_accumulation") + mask_var_constraint = iris.NameConstraint(var_name="hourly_wts_accumulation") + for model, mask in zip( + iter_maybe(model_names), iter_maybe(weights_names), strict=True + ): + print(" model, mask ", model, mask) + + # grab the field to filter + model_constraint = iris.AttributeConstraint(model_name=model) + unfiltered_field = cubes.extract_cube(var_constraint & model_constraint) + + # Select the field to use as the mask. + # Nice to do - put in support for a static mask. + mask_constraint = iris.AttributeConstraint(model_name=mask) + mask_field = cubes.extract_cube(mask_var_constraint & mask_constraint) + + # Create the mask - note that the condition e.g. "ge" can be set by a loader + # as can the threshold value. + mask_radar_wts = generate_mask(mask_field, "ge", 11) + + # print(" This is cube radar_obs: ", radar_obs) + # print(" This is cube radar_weights: ", radar_wts) + # print(" This is cube unfiltered: ", unfiltered_field) + + # check the coords of the unfiltered field and the mask field. + # If these do not match, then regrid the unfiltered field onto + # the grid used for the mask field. + # For radar weights fields can use the function regrid_onto_xyspacing in regrid.py, + # but then might have to extract a subarea to match the mask grid. + # Might have to consider serval cases for regridding: + # (1) model_field(lat, lon) to radar_weights_field(x, y) + # (2) model_field(lat, lon) to other_model_field(lat, lon) + # (3) Nimrod_field(x, y) to radar_weights_field(x, y) + # (4) Nimrod_field(x, y) to model_field(lat,lon) ? + # + + # Apply the mask. + masked_radar_obs = apply_mask(unfiltered_field, mask_radar_wts) + + # Put the filtered cube into the list of filtered cubes. + filtered_list.append(masked_radar_obs) + + # Preserve returning a cube if only a cube has been supplied to filter. + if len(filtered_list) == 1: + return filtered_list[0] + else: + return filtered_list From a1fc5753c35686447941daa66248aa5bd1869a6e Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 4 Jun 2026 11:38:29 +0100 Subject: [PATCH 004/105] Recipe to produce mean time series of masked Nimrod accumulations #2113 --- .../radar_masked_mean_time_series.yaml | 44 +++++++++++++++++++ 1 file changed, 44 insertions(+) create mode 100644 src/CSET/recipes/surface_fields/radar_masked_mean_time_series.yaml diff --git a/src/CSET/recipes/surface_fields/radar_masked_mean_time_series.yaml b/src/CSET/recipes/surface_fields/radar_masked_mean_time_series.yaml new file mode 100644 index 0000000000..dd1f31f2d3 --- /dev/null +++ b/src/CSET/recipes/surface_fields/radar_masked_mean_time_series.yaml @@ -0,0 +1,44 @@ +category: Histogram +#title: "Histogram masked Nimrod $VARNAME\n$RADAR_NAME" +title: "Masked Nimrod mean times series $VARNAME\n$RADAR_NAME" +description: Timeseries of masked Nimrod radar observations. + +steps: + - operator: read.read_cubes + file_paths: $INPUT_PATHS +# model_names: [$RADAR_NAME, $WEIGHTS_NAME] + model_names: $ALL_NAME + subarea_type: $SUBAREA_TYPE + subarea_extent: $SUBAREA_EXTENT + constraint: + operator: constraints.combine_constraints +# variable_constraint: +# operator: constraints.generate_var_constraint +# varname: "Hourly rain accumulation" +# variable_constraint: +# operator: constraints.generate_var_constraint +# varname: "Hourly wts accumulation" + cell_methods_constraint: + operator: constraints.generate_cell_methods_constraint + cell_methods: [] + varname: $VARNAME + pressure_level_constraint: + operator: constraints.generate_level_constraint + coordinate: pressure + levels: [] +# constraint: ["Hourly rain accumulation", "Hourly wts accumulation"] + + + - operator: radar_filter.mask_by_weights + model_names: $RADAR_NAME + weights_names: $WEIGHTS_NAME + + - operator: collapse.collapse + coordinate: [projection_y_coordinate, projection_x_coordinate] + method: MEAN + + - operator: write.write_cube_to_nc + overwrite: True + + - operator: plot.plot_line_series + sequence_coordinate: $SEQUENCE From 4cfbdfd78b00d6181835237e522ad316ee3a1d2f Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 4 Jun 2026 11:39:30 +0100 Subject: [PATCH 005/105] Loader for recipe to produce mean time series of masked Nimrod accumulations #2113 --- src/CSET/loaders/radar.py | 42 +++++++++++++++++++++++++++++++++++---- 1 file changed, 38 insertions(+), 4 deletions(-) diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index 5f53462d6f..f46f26e246 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -158,11 +158,11 @@ def load(conf: Config): # recipe must be done by passing lists of the radar_ids and # the radar_names. As this is a multiline plot, all radar sources # share the same radar variable name. - radar_obs_ids = [radar["id"] for radar in accum_radars] - radar_wts_ids = [radar["id"] for radar in wts_radars] - combined_ids = radar_obs_ids + radar_wts_ids - print(" combined_ids: ", combined_ids) if conf.HISTOGRAM_SURFACE_FIELD: + radar_obs_ids = [radar["id"] for radar in accum_radars] + radar_wts_ids = [radar["id"] for radar in wts_radars] + combined_ids = radar_obs_ids + radar_wts_ids + print(" combined_ids: ", combined_ids) yield RawRecipe( recipe="radar_dev3.yaml", # model_ids -> Becomes $INPUT_PATHS @@ -202,6 +202,40 @@ def load(conf: Config): aggregation=False, ) + # Timeseries for surface (2D) Nimrod radar hourly accumulated rainfall. + # + # The timeseries are produced after the rainfall obs have been masked using + # the associated Nimrod weights file. + # + # To get multiple radar sources plotted on the histogram the + # recipe must be done by passing lists of the radar_ids and + # the radar_names. As this is a multiline plot, all radar sources + # share the same radar variable name. + if conf.TIMESERIES_SURFACE_FIELD: + radar_obs_ids = [radar["id"] for radar in accum_radars] + radar_wts_ids = [radar["id"] for radar in wts_radars] + combined_ids = radar_obs_ids + radar_wts_ids + print(" combined_ids: ", combined_ids) + yield RawRecipe( + recipe="radar_masked_mean_time_series.yaml", + # model_ids -> Becomes $INPUT_PATHS + # model_ids=[ radar_obs_ids, radar_wts_ids], + model_ids=combined_ids, + variables={ + "VARNAME": next(radar["varname"] for radar in accum_radars), + "ALL_NAME": combined_ids, + "RADAR_NAME": [radar["id"] for radar in accum_radars], + "WEIGHTS_NAME": [radar["id"] for radar in wts_radars], + "SEQUENCE": "realisation", + # "SEQUENCE": "time" + # if conf.HISTOGRAM_SURFACE_FIELD_SEQUENCE + # else "realization", + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, + }, + aggregation=False, + ) + # Timeseries plot of Nimrod hourly surface rainfall accumulation. if conf.TIMESERIES_SURFACE_FIELD and accum_radars: yield RawRecipe( From fa0f7cd099778dde581cda69468e0acb2ea64faf Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 4 Jun 2026 16:16:50 +0100 Subject: [PATCH 006/105] Nimrod weights 2-d plots to use gray for out of radar range domain #2113 --- src/CSET/operators/_colormaps.py | 2 +- tests/operators/test_colormaps.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/CSET/operators/_colormaps.py b/src/CSET/operators/_colormaps.py index 6a6e3800ce..f3444ac7c4 100644 --- a/src/CSET/operators/_colormaps.py +++ b/src/CSET/operators/_colormaps.py @@ -544,6 +544,7 @@ def custom_colourmap_nimrod_weights(cube: iris.cube.Cube, cmap, levels, norm): ] norm = mcolors.BoundaryNorm(levels, cmap.N) colours = [ + "#dcdcdc", "#d10000", "purple", "#8f00d6", @@ -556,7 +557,6 @@ def custom_colourmap_nimrod_weights(cube: iris.cube.Cube, cmap, levels, norm): "#37a648", "#8edc64", "#c5ffc5", - "#dcdcdc", "#ffffff", ] # Create a custom colormap. diff --git a/tests/operators/test_colormaps.py b/tests/operators/test_colormaps.py index 513f9d0f72..2a8e792ae6 100644 --- a/tests/operators/test_colormaps.py +++ b/tests/operators/test_colormaps.py @@ -404,6 +404,7 @@ def test_colorbar_map_nimrod_wts(cube): cube.rename("Hourly wts accumulation") expected_levels = np.arange(-0.5, 14.5, 1.0) expected_colors = [ + "#dcdcdc", "#d10000", "purple", "#8f00d6", @@ -416,7 +417,6 @@ def test_colorbar_map_nimrod_wts(cube): "#37a648", "#8edc64", "#c5ffc5", - "#dcdcdc", "#ffffff", ] expected_cmap = mpl.colors.ListedColormap(expected_colors) From d6d441bc445ad2ca2ced5342b2ce8b380519b0c5 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 18 Jun 2026 11:23:40 +0100 Subject: [PATCH 007/105] jaml recipe for plotting masked model rainfall #2113 --- .../recipes/surface_fields/radar_masked_mean_time_series.yaml | 1 - 1 file changed, 1 deletion(-) diff --git a/src/CSET/recipes/surface_fields/radar_masked_mean_time_series.yaml b/src/CSET/recipes/surface_fields/radar_masked_mean_time_series.yaml index dd1f31f2d3..83184b8861 100644 --- a/src/CSET/recipes/surface_fields/radar_masked_mean_time_series.yaml +++ b/src/CSET/recipes/surface_fields/radar_masked_mean_time_series.yaml @@ -1,5 +1,4 @@ category: Histogram -#title: "Histogram masked Nimrod $VARNAME\n$RADAR_NAME" title: "Masked Nimrod mean times series $VARNAME\n$RADAR_NAME" description: Timeseries of masked Nimrod radar observations. From 93825e25ab6e3fa1f4c74245ebaf3109e7f2f101 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 18 Jun 2026 11:34:10 +0100 Subject: [PATCH 008/105] function radar_apply_mask in radar_filter.py #2113 --- src/CSET/operators/radar_filter.py | 152 ++++++++++++++++++++++++++++- 1 file changed, 151 insertions(+), 1 deletion(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index 4350099f92..2d4a4c3d04 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -14,9 +14,12 @@ """Operators to perform various kind of filtering.""" +import logging + import iris import iris.cube import iris.exceptions +import numpy as np from CSET._common import iter_maybe from CSET.operators.filters import apply_mask, generate_mask @@ -150,7 +153,10 @@ def mask_by_weights( var_constraint = iris.NameConstraint(var_name="hourly_rain_accumulation") mask_var_constraint = iris.NameConstraint(var_name="hourly_wts_accumulation") for model, mask in zip( - iter_maybe(model_names), iter_maybe(weights_names), strict=True + iter_maybe(model_names), + iter_maybe(weights_names), + strict=True, + # iter_maybe(model_names), iter_maybe(weights_names), strict=True ): print(" model, mask ", model, mask) @@ -194,3 +200,147 @@ def mask_by_weights( return filtered_list[0] else: return filtered_list + + +def radar_apply_mask( + original_field: iris.cube.Cube | iris.cube.CubeList, + mask: iris.cube.Cube | iris.cube.CubeList, +) -> iris.cube.Cube | iris.cube.CubeList: + """Apply a mask to given data as a masked array. + + Parameters + ---------- + original_field: iris.cube.Cube | iris.cube.CubeList + The field(s) to be masked. + mask: iris.cube.Cube | iris.cube.CubeList + The mask(s) being applied to the original field(s). + + Returns + ------- + masked_field: iris.cube.Cube | iris.cube.CubeList + A cube or cubelist of the masked field(s). + + Notes + ----- + The mask is first converted to 1s and NaNs before multiplication with + the original data. + + As discussed in generate_mask, you can combine multiple masks in a + recipe using other functions before applying the mask to the data. + + Examples + -------- + >>> land_points_only = apply_mask(temperature, land_mask) + """ + print("original_field cube: ", original_field) + print("mask cube", mask) + + # Create an empty cubelist to hold the filtered fields. + masked_fields = iris.cube.CubeList([]) + + # Loop over the input mask and field cubes. + for M, F in zip(iter_maybe(mask), iter_maybe(original_field), strict=True): + # Ensure mask data are floats and only 1s or NaNs. + print(" type(M.data)", type(M.data)) + M.data = np.float64(M.data) + M.data[M.data == 0.0] = np.nan + M.data[~np.isnan(M.data)] = 1.0 + print(" type(M.data) v2", type(M.data)) + print("---> M[0][300][:].data ", M[0][300][:].data) + print( + "---> max(M) ", max(M[0][300][:].data), " min(M) ", min(M[0][300][:].data) + ) + logging.info( + "Mask set to 1 or 0s, if addition of multiple masks results" + "in values > 1 these are set to 1." + ) + + print("") + print("---> here is F[0]") + print(F[0]) + print("") + + print("") + print("---> here is F[1]") + print(F[1]) + print("") + + # print("") + # print("---> here is F[24]") + # print(F[24]) + # print("") + + # Apply the mask + masked_field = F.copy() + print(" M.shape is ", M.shape) + print(" masked_field.shape is ", masked_field.shape) + print(" M[0].shape is ", M[0].shape) + print(" masked_field[0].shape is ", masked_field[0].shape) + + # If the field and mask on on different grids, then regrid the field. + if M[0].shape != masked_field[0].shape: + regridded_cube = masked_field.regrid(M, iris.analysis.Linear()) + masked_field = regridded_cube + + print(" ---> Have regridded the field") + + print("") + print("---> here is masked_field[0]") + print(masked_field[0]) + print("") + + print(" M.shape is ", M.shape) + print(" masked_field.shape is ", masked_field.shape) + print(" M[0].shape is ", M[0].shape) + print(" masked_field[0].shape is ", masked_field[0].shape) + + print("") + print("---> check how many time instances there are of both field and mask") + print("---> masked_field.shape[0] ", masked_field.shape[0]) + print("---> M.shape[0] ", M.shape[0]) + + print("") + # print("---> max(M)", max(list(M.data))) + + loop_i = min(M.shape[0], masked_field.shape[0]) + for i in range(loop_i): + print(" i is ", i) + # print("---> max(M[i]) ", max(M[i].data), " min(M[i]) ", min(M[i].data)) + # print("---> max(M[i]) ", max(list(M[i].data)), " min(M[i]) ", min(list(M[i].data))) + # masked_field[i].data *= M[i].data + # masked_field[i] = apply_mask(masked_field[i], M[i]) + # mask_array = M[i][:][:].data.filled(np.array) + mask_array = M[i].data.filled(np.array) + print(" type(M.data)", type(M.data)) + print(" type(mask_array)", type(mask_array)) + masked_field[i][:][:].data *= mask_array + # masked_field[i][:][:].data *= M[i][:][:].data + # print("---> max(mask_array) ", max(mask_array), " min(mask_array ) ", min(mask_array) ) + print( + "---> max(M[i][300]) ", + max(M[i][300].data), + " min(M[i][300]) ", + min(M[i][300].data), + ) + print( + "---> max(F[i][300]) ", + max(masked_field[i][300].data), + " min(F[i][300]) ", + min(masked_field[i][300].data), + ) + # print(" check 1 ") + masked_field[i].attributes["mask"] = f"mask_of_{F.name()}" + # print(" check 2 ") + masked_fields.append(masked_field[i]) + + # if M[0].shape == F[0].shape: + # masked_field.data *= M.data + # masked_field.attributes["mask"] = f"mask_of_{F.name()}" + # masked_fields.append(masked_field) + + print(" len(masked_fields) ", len(masked_fields)) + if len(masked_fields) == 1: + return masked_fields[0] + else: + # return masked_fields + return masked_fields.merge() From 99cc5370d4524855e5abb9c2b9d22a37c5cc52e4 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 18 Jun 2026 15:57:57 +0100 Subject: [PATCH 009/105] radar_apply_mask now uses function apply_mask in filters.py #2113 --- src/CSET/operators/radar_filter.py | 118 ++++------------------------- 1 file changed, 15 insertions(+), 103 deletions(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index 2d4a4c3d04..5268b176d3 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -14,12 +14,9 @@ """Operators to perform various kind of filtering.""" -import logging - import iris import iris.cube import iris.exceptions -import numpy as np from CSET._common import iter_maybe from CSET.operators.filters import apply_mask, generate_mask @@ -230,115 +227,30 @@ def radar_apply_mask( Examples -------- - >>> land_points_only = apply_mask(temperature, land_mask) + >>> land_points_only = radar_apply_mask( surface_microphysical_rainfall_rate, Nimrod2km) """ - print("original_field cube: ", original_field) - print("mask cube", mask) - # Create an empty cubelist to hold the filtered fields. masked_fields = iris.cube.CubeList([]) # Loop over the input mask and field cubes. for M, F in zip(iter_maybe(mask), iter_maybe(original_field), strict=True): - # Ensure mask data are floats and only 1s or NaNs. - print(" type(M.data)", type(M.data)) - M.data = np.float64(M.data) - M.data[M.data == 0.0] = np.nan - M.data[~np.isnan(M.data)] = 1.0 - print(" type(M.data) v2", type(M.data)) - print("---> M[0][300][:].data ", M[0][300][:].data) - print( - "---> max(M) ", max(M[0][300][:].data), " min(M) ", min(M[0][300][:].data) - ) - logging.info( - "Mask set to 1 or 0s, if addition of multiple masks results" - "in values > 1 these are set to 1." - ) - - print("") - print("---> here is F[0]") - print(F[0]) - print("") - - print("") - print("---> here is F[1]") - print(F[1]) - print("") - - # print("") - # print("---> here is F[24]") - # print(F[24]) - # print("") - - # Apply the mask masked_field = F.copy() - print(" M.shape is ", M.shape) - print(" masked_field.shape is ", masked_field.shape) - print(" M[0].shape is ", M[0].shape) - print(" masked_field[0].shape is ", masked_field[0].shape) - # If the field and mask on on different grids, then regrid the field. + # If the field and mask are on different grids, then regrid the field. if M[0].shape != masked_field[0].shape: - regridded_cube = masked_field.regrid(M, iris.analysis.Linear()) - masked_field = regridded_cube - - print(" ---> Have regridded the field") - - print("") - print("---> here is masked_field[0]") - print(masked_field[0]) - print("") - - print(" M.shape is ", M.shape) - print(" masked_field.shape is ", masked_field.shape) - print(" M[0].shape is ", M[0].shape) - print(" masked_field[0].shape is ", masked_field[0].shape) - - print("") - print("---> check how many time instances there are of both field and mask") - print("---> masked_field.shape[0] ", masked_field.shape[0]) - print("---> M.shape[0] ", M.shape[0]) - - print("") - # print("---> max(M)", max(list(M.data))) - - loop_i = min(M.shape[0], masked_field.shape[0]) - for i in range(loop_i): - print(" i is ", i) - # print("---> max(M[i]) ", max(M[i].data), " min(M[i]) ", min(M[i].data)) - # print("---> max(M[i]) ", max(list(M[i].data)), " min(M[i]) ", min(list(M[i].data))) - # masked_field[i].data *= M[i].data - # masked_field[i] = apply_mask(masked_field[i], M[i]) - # mask_array = M[i][:][:].data.filled(np.array) - mask_array = M[i].data.filled(np.array) - print(" type(M.data)", type(M.data)) - print(" type(mask_array)", type(mask_array)) - masked_field[i][:][:].data *= mask_array - # masked_field[i][:][:].data *= M[i][:][:].data - # print("---> max(mask_array) ", max(mask_array), " min(mask_array ) ", min(mask_array) ) - print( - "---> max(M[i][300]) ", - max(M[i][300].data), - " min(M[i][300]) ", - min(M[i][300].data), - ) - print( - "---> max(F[i][300]) ", - max(masked_field[i][300].data), - " min(F[i][300]) ", - min(masked_field[i][300].data), - ) - # print(" check 1 ") - masked_field[i].attributes["mask"] = f"mask_of_{F.name()}" - # print(" check 2 ") - masked_fields.append(masked_field[i]) - - # if M[0].shape == F[0].shape: - # masked_field.data *= M.data - # masked_field.attributes["mask"] = f"mask_of_{F.name()}" - # masked_fields.append(masked_field) - - print(" len(masked_fields) ", len(masked_fields)) + masked_field = masked_field.regrid(M, iris.analysis.Linear()) + + # Apply the mask. + min_timesteps = min(M.shape[0], masked_field.shape[0]) + masked_field = apply_mask(masked_field[0:min_timesteps], M[0:min_timesteps]) + + # Attach and attribute to the masked field detailing the mask used. + masked_field.attributes["mask"] = f"mask_of_{F.name()}" + + # Append the masked field to the output list of masked fields. + masked_fields.append(masked_field) + + # Return either a single cube or a cubelist. if len(masked_fields) == 1: return masked_fields[0] else: From 326873644dabb83c0c50dc4d2ebbac5db4eae752 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 18 Jun 2026 15:59:59 +0100 Subject: [PATCH 010/105] loaders for radar_mask_model.yaml #2113 --- src/CSET/loaders/radar.py | 88 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 88 insertions(+) diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index f46f26e246..bd330303fc 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -85,6 +85,9 @@ def get_radar_sources(conf) -> list[dict]: def load(conf: Config): """Yield recipes from the given workflow configuration.""" + # Load a list of model detail dictionaries. + # models = get_models(conf.asdict()) + # Load the required radar observation sources. radar_sources = get_radar_sources(conf) @@ -102,6 +105,91 @@ def load(conf: Config): if radar["varname"] == "Hourly wts accumulation" ] + # # Radar masking based on sea mask. + # if conf.SPATIAL_SURFACE_FIELD: + # for field in conf.SURFACE_FIELDS: + # yield RawRecipe( + # recipe="sea_mask_for_surface_domain_mean_time_series.yaml", + # variables={ + # "VARNAME": field, + # "MODEL_NAME": [model["name"] for model in models], + # "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + # "SUBAREA_EXTENT": conf.SUBAREA_EXTENT + # if conf.SELECT_SUBAREA + # else None, + # "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", + # }, + # model_ids=[model["id"] for model in models], + # aggregation=False, + # ) + + # Radar masking of radar obs based on sea mask. + if conf.SPATIAL_SURFACE_FIELD: + field = "Hourly rain accumulation" + yield RawRecipe( + recipe="radar_mask_model.yaml", + variables={ + "VARNAME": field, + "MODEL_LABEL": "Nimrod2km", + "MASK_LABEL": "Nimrod2km", + "METHOD": "SEQ", + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, + "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", + }, + model_ids=["Nimrod2km", "Nimrod2km_weights"], + aggregation=False, + ) + + # Radar masking of model rainfall based on sea mask. + if conf.SPATIAL_SURFACE_FIELD: + field = "surface_microphysical_rainfall_rate" + yield RawRecipe( + recipe="radar_mask_model.yaml", + variables={ + "VARNAME": field, + "MODEL_LABEL": "ModelA", + "MASK_LABEL": "Nimrod2km", + "METHOD": "SEQ", + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, + "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", + }, + model_ids=["1", "Nimrod2km_weights"], + aggregation=False, + ) + + # Surface (2D) fields for model rainfall masked by Nimrod radar. + # + # The different sources of Nimrod rainfall accumulation have + # different spatial grids. So each source requires its own + # recipe to prevent incompatible cubes being created. + # if conf.SPATIAL_SURFACE_FIELD: + # radar_source = ["Nimrod_2km"] + # for radar in radar_source: + # model_labels = [model["id"] for model in models] + # radar_label = ["Nimrod2km_weights"] + # combined_ids = [model_labels[0]] + radar_label + # print("Combined ids is: ", combined_ids) + # yield RawRecipe( + # recipe="radar_plot_sequence_rainfall.yaml", + ## model_ids=radar["id"], # -> Becomes $INPUT_PATHS + # model_ids=combined_ids, + # variables={ + ## "VARNAME": radar["varname"], + ## "RADAR_NAME": radar["name"], + # "RADAR_NAME": "Nimrod_2km_weights", + ## "MODEL_NAME": [model["name"] for model in models], + # "MODEL_NAME": "ModelA", + # "METHOD": "SEQ", + # "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + # "SUBAREA_EXTENT": conf.SUBAREA_EXTENT + # if conf.SELECT_SUBAREA + # else None, + # }, + # aggregation=False, + # ) + # Surface (2D) fields for Nimrod radar rainfall. # # The different sources of Nimrod rainfall accumulation have From 33a20f2f5cf0a1e48c5f449b6296b30b122461b6 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 18 Jun 2026 16:01:38 +0100 Subject: [PATCH 011/105] for 2-d surface plots replaced np.min, np.max, np.mean with np.nanmin, np.nanmax, np.nanmean #2113 --- src/CSET/operators/plot.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index a8a72a51e4..501176f3d8 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -717,7 +717,7 @@ def _plot_and_save_spatial_plot( # Add watermark with min/max/mean. Currently not user togglable. # In the bbox dictionary, fc and ec are hex colour codes for grey shade. axes.annotate( - f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}", + f"Min: {np.nanmin(cube.data):.3g} Max: {np.nanmax(cube.data):.3g} Mean: {np.nanmean(cube.data):.3g}", xy=(0.025, yinfopad), xycoords="axes fraction", xytext=(-5, 5), From 4b81b55f129f453b61cc2aebcbb241ab5a93f2f2 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Fri, 19 Jun 2026 10:24:00 +0100 Subject: [PATCH 012/105] Grey out domain beyond radar range for rainfall masked spatial plots #2113 --- src/CSET/operators/_colormaps.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/CSET/operators/_colormaps.py b/src/CSET/operators/_colormaps.py index f3444ac7c4..dd042043ef 100644 --- a/src/CSET/operators/_colormaps.py +++ b/src/CSET/operators/_colormaps.py @@ -512,6 +512,10 @@ def custom_colormap_precipitation(cube: iris.cube.Cube, cmap, levels, norm): # Normalize the levels norm = mcolors.BoundaryNorm(levels, cmap.N) logging.info("Using custom rainfall colourmap.") + + # Set any Nan values to be plotted a light grey. + cmap.set_bad("#dcdcdc") + return cmap, levels, norm From 493ddfa4c205fcf32ec0c38e6cdd7fbb431b276a Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Mon, 20 Jul 2026 14:53:26 +0100 Subject: [PATCH 013/105] When regridding model using nimrod wts field, set extrapolated values to NaN #2113 --- src/CSET/operators/radar_filter.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index 5268b176d3..8f5473417a 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -238,7 +238,8 @@ def radar_apply_mask( # If the field and mask are on different grids, then regrid the field. if M[0].shape != masked_field[0].shape: - masked_field = masked_field.regrid(M, iris.analysis.Linear()) + scheme = iris.analysis.Linear(extrapolation_mode="nan") + masked_field = masked_field.regrid(M, scheme) # Apply the mask. min_timesteps = min(M.shape[0], masked_field.shape[0]) From 79d947636a831df83e63b8ec6358e28648a0a934 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 27 Aug 2026 09:38:33 +0100 Subject: [PATCH 014/105] Radar operators in radar_filter.py #2113 --- src/CSET/operators/radar_filter.py | 222 ++++++++++++++++++++++++++++- 1 file changed, 221 insertions(+), 1 deletion(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index 8f5473417a..e4b0019fcb 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -17,6 +17,7 @@ import iris import iris.cube import iris.exceptions +import numpy as np from CSET._common import iter_maybe from CSET.operators.filters import apply_mask, generate_mask @@ -56,7 +57,7 @@ def mask_list(model_names: list[str]) -> list[str]: Examples -------- - >>> list_weights = make_list( ["UM_model", "Nimrod1km", "Nimrod2km"] ) + >>> list_weights = mask_list( ["UM_model", "Nimrod1km", "Nimrod2km"] ) >>> print(list_weights) ["Nimrod2km_weights", "Nimrod1km_weights", "Nimrod2km_weights"] @@ -202,6 +203,7 @@ def mask_by_weights( def radar_apply_mask( original_field: iris.cube.Cube | iris.cube.CubeList, mask: iris.cube.Cube | iris.cube.CubeList, + boundary_margin: int = 8, ) -> iris.cube.Cube | iris.cube.CubeList: """Apply a mask to given data as a masked array. @@ -211,6 +213,9 @@ def radar_apply_mask( The field(s) to be masked. mask: iris.cube.Cube | iris.cube.CubeList The mask(s) being applied to the original field(s). + boundary_margin: int, optional + Number of grid points from the domain boundary considered "unreliable". + Defaults to 8. Returns ------- @@ -236,6 +241,28 @@ def radar_apply_mask( for M, F in zip(iter_maybe(mask), iter_maybe(original_field), strict=True): masked_field = F.copy() + # TODO set the model perimeter to NaN as these gridpoints contain no data. + # c.f. boundary_margin in regrid.py + # Get bounds + # Get axis + # lat, lon = M.coord(y_coord), M.coord(x_coord) + # lat_min, lon_min = lat.points.min(), lon.points.min() + # lat_max, lon_max = lat.points.max(), lon.points.max() + margin_width = boundary_margin + if margin_width > 0: + masked_field.data[:, -margin_width - 1 :, :] = np.nan + masked_field.data[:, :, -margin_width - 1 :] = np.nan + masked_field.data[:, :, 0:margin_width] = np.nan + masked_field.data[:, 0:margin_width, :] = np.nan + # masked_field.data[:, -margin_width-1:-1, :] = 3.0 + # masked_field.data[:, :, -margin_width-1:-1] = 1.5 + # masked_field.data[:, :, 0:margin_width] = 50.0 + # masked_field.data[:, 0:margin_width, :] = 30.0 + # masked_field.data[ masked_field.data < 0.01 ] = 200.0 + # masked_field.data[ masked_field.data == np.nan ] = 20.0 + + # TODO must mask the radar obs using the model domain + # If the field and mask are on different grids, then regrid the field. if M[0].shape != masked_field[0].shape: scheme = iris.analysis.Linear(extrapolation_mode="nan") @@ -257,3 +284,196 @@ def radar_apply_mask( else: # return masked_fields return masked_fields.merge() + + +def radar_mask( + model_field: iris.cube.Cube | iris.cube.CubeList, + nimrod_field: iris.cube.Cube | iris.cube.CubeList, + nimrod_mask: iris.cube.Cube | iris.cube.CubeList, + boundary_margin: int = 8, + outputs: str = "Nimrod", +) -> iris.cube.Cube | iris.cube.CubeList: + """Apply a mask to given data as a masked array.""" + # Create an empty cubelist to hold the filtered fields. + filtered_fields = iris.cube.CubeList([]) + filtered_radar = iris.cube.CubeList([]) + filtered_model = iris.cube.CubeList([]) + + # ensure the three inputs to this function are model_field, nimrod_field, nimrod_mask + + for M, F, N in zip( + iter_maybe(nimrod_mask), + iter_maybe(model_field), + iter_maybe(nimrod_field), + strict=True, + ): + # masked_field = F.copy() + + # apply the function radar_apply_mask to generate the re-gridded and masked model_field + # i.e. generate masked_model_field + masked_model_field = radar_apply_mask(F, M, boundary_margin=boundary_margin) + + # use the model_field as the mask for the nimrod_field - note: no re-gridding required + # i.e. generate masked_nimrod_field + + min_timesteps = min(N.shape[0], masked_model_field.shape[0]) + + temp_mask = masked_model_field[0:min_timesteps].copy() + # temp_mask.data[ temp_mask.data != np.nan ] = 1.0 + temp_mask.data[~np.isnan(temp_mask.data)] = 1.0 + + masked_nimrod_field = N[0:min_timesteps].copy() + # masked_nimrod_field.data = N[0:min_timesteps].data * temp_mask[0:min_timesteps].data + masked_nimrod_field.data *= temp_mask.data + + # Append the masked field to the output list of masked fields. + filtered_model.append(masked_model_field) + filtered_radar.append(masked_nimrod_field) + + # return masked_model_field and masked_nimrod_field --> these can then be passed to either + # the histogram or time series plotting operators. + + print(" test point bmc1") + if outputs == "Nimrod": + filtered_fields.append(filtered_radar.merge_cube()) + print(" test point bmc2") + if outputs == "model": + filtered_fields.append(filtered_model.merge_cube()) + print(" test point bmc3") + if outputs == "all": + print(" test point bmc4") + filtered_fields.append(filtered_model.merge_cube()) + print(" test point bmc5") + filtered_fields.append(filtered_radar.merge_cube()) + + print("+++++++++++++++++++++++++++++") + print(filtered_fields) + print("+++++++++++++++++++++++++++++") + # Return either a single cube or a cubelist. + print(" test point bmc6") + if len(filtered_fields) == 1: + print(" test point bmc7a") + return filtered_fields[0] + else: + # return masked_fields + # return filtered_fields.merge() + print(" test point bmc7") + return filtered_fields + + +def radar_mask_loop( + model_field: iris.cube.Cube | iris.cube.CubeList, + nimrod_field: iris.cube.Cube | iris.cube.CubeList, + nimrod_mask: iris.cube.Cube | iris.cube.CubeList, + boundary_margin: int = 8, + outputs: str = "Nimrod", +) -> iris.cube.Cube | iris.cube.CubeList: + """Find common domains between a list of models and radar observations.""" + # Create an empty cubelist to hold the filtered fields. + filtered_cubes = iris.cube.CubeList([]) + + # if len(nimrod_field) > 1: + # use_nimrod_field = nimrod_field[0] + # use_nimrod_mask = nimrod_mask[0] + # else: + print("-----> len(nimrod_field) ", len(nimrod_field)) + use_nimrod_field = nimrod_field + use_nimrod_mask = nimrod_mask + + # Loop over the models. + for model in model_field: + print("-------> using model ", model) + filtered_model = radar_mask( + model, + use_nimrod_field, + use_nimrod_mask, + boundary_margin=boundary_margin, + outputs="model", + ) + filtered_cubes.append(filtered_model) + + # Grab the filtered radar observations. + print("-----> len(model_field) ", len(model_field)) + if len(model_field) == 1: + filtered_radar = radar_mask( + model_field[0], + use_nimrod_field, + use_nimrod_mask, + boundary_margin=boundary_margin, + outputs="Nimrod", + ) + else: + filtered_radar = radar_mask( + model_field[0], + use_nimrod_field, + use_nimrod_mask, + boundary_margin=boundary_margin, + outputs="Nimrod", + ) + + filtered_cubes.append(filtered_radar) + + print("-----------------> returning from radar_mask_loop") + return filtered_cubes + + +def match_varname_and_units(cubes: iris.cube.Cube | iris.cube.CubeList): + """Match the varname and units of a cube list. + + Arguments + --------- + cubes: iris.cube.Cube | iris.cube.CubeList + A Cube or CubeList of a field to be matched. + + Returns + ------- + iris.cube.Cube | iris.cube.CubeList + The matched cubes. + + + Notes + ----- + This function converts the names and units of a cube list to match + the first cube in the list. If just one cube is input, then this is + returned. + """ + # If just one cube, then no need to match so return. + if len(cubes) == 1: + return cubes + + print("@@@@@@@@@@@@@@@@@@@@@ cubes[0] @@@@@@@@@@@@") + print(cubes[0]) + print("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@") + + # Initialise the list of matched cubes. + new_cubelist = iris.cube.CubeList([]) + new_cubelist.append(cubes[0]) + + # Loop the cubes to match to the first cube. + base_cube = cubes[0] + print("@@@@@@@@@@@@@@@@@@@@@ base_cube @@@@@@@@@@@@") + print("base_cube.name : ", base_cube.name) + print("base_cube.standard_name : ", base_cube.standard_name) + print("base_cube.long_name : ", base_cube.long_name) + print("base_cube.var_name : ", base_cube.var_name) + print("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@") + for cube in cubes[1:]: + print("@@@@@@@@@@@@@@@ raw @@@@@@@@@@@@@@@@@@") + print(cube) + print("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@") + new_cube = cube.copy() + new_cube.rename(base_cube.long_name) + new_cube.long_name = base_cube.long_name + new_cube.units = base_cube.units + print("@@@@@@@@@@@@@@@ matched @@@@@@@@@@@@@@@@@@") + print("new_cube.name : ", new_cube.name) + print("new_cube.standard_name : ", new_cube.standard_name) + print("new_cube.long_name : ", new_cube.long_name) + print("new_cube.var_name : ", new_cube.var_name) + print(new_cube) + print("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@") + + # Append the matched cube to the output cube list. + new_cubelist.append(new_cube) + + return new_cubelist From da7576d9154738beea90882f4333d4f9763b8788 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 27 Aug 2026 09:40:24 +0100 Subject: [PATCH 015/105] Radar loaders #2113 --- src/CSET/loaders/radar.py | 180 +++++++++++++++++++++++++++++++++++++- 1 file changed, 178 insertions(+), 2 deletions(-) diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index bd330303fc..175260b78e 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -14,7 +14,20 @@ """Load radar observation recipes.""" -from CSET.recipes import Config, RawRecipe +from CSET.recipes import Config, RawRecipe, get_models + + +def select_radar_source(input_source_list): + """Select the preferred radar source from a list of sources.""" + # Preference list for Nimrod sources, with most preferred heading the list. + nimrod_preference = ["Nimrod2km", "Nimrodxkm", "Nimrod1km"] + + empty_string = "" + preferred_nimrod = empty_string + for prefer in reversed(nimrod_preference): + if any(prefer in source for source in input_source_list): + preferred_nimrod = prefer + return preferred_nimrod def get_radar_sources(conf) -> list[dict]: @@ -86,7 +99,7 @@ def get_radar_sources(conf) -> list[dict]: def load(conf: Config): """Yield recipes from the given workflow configuration.""" # Load a list of model detail dictionaries. - # models = get_models(conf.asdict()) + models = get_models(conf.asdict()) # Load the required radar observation sources. radar_sources = get_radar_sources(conf) @@ -123,8 +136,171 @@ def load(conf: Config): # aggregation=False, # ) + # Common domain between model and radar observations. + # if conf.SPATIAL_SURFACE_FIELD: + if False: + # field = "Hourly rain accumulation" + field = "surface_microphysical_rainfall_rate" + yield RawRecipe( + recipe="radar_mask.yaml", + variables={ + "VARNAME": field, + "MODEL_LABEL": "ModelA", + "NIMROD_LABEL": "Nimrod2km", + "MASK_LABEL": "Nimrod2km_wts", + "SEQUENCE": "time", + "OUTPUTS": "all", + "METHOD": "SEQ", + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, + "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", + }, + model_ids=["1", "Nimrod2km", "Nimrod2km_weights"], + aggregation=False, + ) + + # Histogram sequence using common domain between + # model rainfall and radar observations. + if conf.HISTOGRAM_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: + # Select the radar source to use. + radar_source = select_radar_source([source["id"] for source in accum_radars]) + radar_obs_ids = [radar_source] + radar_wts_ids = [radar_source + "_weights"] + model_names_list = [model["name"] for model in models] + model_ids_list = [model["id"] for model in models] + combined_names = model_names_list + radar_obs_ids + radar_wts_ids + combined_ids = model_ids_list + radar_obs_ids + radar_wts_ids + yield RawRecipe( + recipe="radar_common_domain_histogram_sequence.yaml", + variables={ + "VARNAME": "surface_microphysical_rainfall_rate", + "ALL_LABEL": combined_names, + "MODEL_LABEL": model_names_list, + "NIMROD_LABEL": radar_source, + "SEQUENCE": "time", + "OUTPUTS": "all", + "METHOD": "SEQ", + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, + "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", + }, + model_ids=combined_ids, + aggregation=False, + ) + + # Histogram domain mean using common domain between + # model rainfall and radar observations. + # if conf.HISTOGRAM_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: + if False: + # Select the radar source to use. + radar_source = select_radar_source([source["id"] for source in accum_radars]) + radar_obs_ids = [radar_source] + radar_wts_ids = [radar_source + "_weights"] + model_names_list = [model["name"] for model in models] + model_ids_list = [model["id"] for model in models] + combined_names = model_names_list + radar_obs_ids + radar_wts_ids + combined_ids = model_ids_list + radar_obs_ids + radar_wts_ids + yield RawRecipe( + recipe="radar_common_domain_histogram_domain_mean.yaml", + variables={ + "VARNAME": "surface_microphysical_rainfall_rate", + "ALL_LABEL": combined_names, + "MODEL_LABEL": model_names_list, + "NIMROD_LABEL": radar_source, + "SEQUENCE": "time", + "OUTPUTS": "all", + "METHOD": "SEQ", + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, + "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", + }, + model_ids=combined_ids, + aggregation=False, + ) + + # Domain mean timeseries using common domain between + # model rainfall and radar observations. + if False: + # if conf.HISTOGRAM_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: + # Select the radar source to use. + radar_source = select_radar_source([source["id"] for source in accum_radars]) + radar_obs_ids = [radar_source] + radar_wts_ids = [radar_source + "_weights"] + model_names_list = [model["name"] for model in models] + model_ids_list = [model["id"] for model in models] + combined_names = model_names_list + radar_obs_ids + radar_wts_ids + combined_ids = model_ids_list + radar_obs_ids + radar_wts_ids + yield RawRecipe( + recipe="radar_common_domain_timeseries_domain_mean.yaml", + variables={ + "VARNAME": "surface_microphysical_rainfall_rate", + "ALL_LABEL": combined_names, + "MODEL_LABEL": model_names_list, + "NIMROD_LABEL": radar_source, + "SEQUENCE": "time", + "OUTPUTS": "all", + "METHOD": "SEQ", + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, + "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", + }, + model_ids=combined_ids, + aggregation=False, + ) + + # Domain timeseries sequence using common domain between + # model rainfall and radar observations. + # if False: + if conf.HISTOGRAM_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: + # Select the radar source to use. + radar_source = select_radar_source([source["id"] for source in accum_radars]) + radar_obs_ids = [radar_source] + radar_wts_ids = [radar_source + "_weights"] + model_names_list = [model["name"] for model in models] + model_ids_list = [model["id"] for model in models] + combined_names = model_names_list + radar_obs_ids + radar_wts_ids + combined_ids = model_ids_list + radar_obs_ids + radar_wts_ids + yield RawRecipe( + recipe="radar_common_domain_timeseries_sequence.yaml", + variables={ + "VARNAME": "surface_microphysical_rainfall_rate", + "ALL_LABEL": combined_names, + "MODEL_LABEL": model_names_list, + "NIMROD_LABEL": radar_source, + "SEQUENCE": "time", + "OUTPUTS": "all", + "METHOD": "SEQ", + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, + "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", + }, + model_ids=combined_ids, + aggregation=False, + ) + # Radar masking of radar obs based on sea mask. if conf.SPATIAL_SURFACE_FIELD: + # field = "Hourly rain accumulation" + field = "surface_microphysical_rainfall_rate" + yield RawRecipe( + recipe="radar_mask_model_obs.yaml", + variables={ + "VARNAME": field, + "MODEL_LABEL": "ModelA", + "NIMROD_LABEL": "Nimrod2km", + "MASK_LABEL": "Nimrod2km_wts", + "METHOD": "SEQ", + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, + "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", + }, + model_ids=["1", "Nimrod2km", "Nimrod2km_weights"], + aggregation=False, + ) + + # Radar masking of radar obs based on sea mask. + # if conf.SPATIAL_SURFACE_FIELD: + if False: field = "Hourly rain accumulation" yield RawRecipe( recipe="radar_mask_model.yaml", From 2d973f2dc4395b6e2d455b9eab5357850e6e3cab Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 27 Aug 2026 09:43:28 +0100 Subject: [PATCH 016/105] Radar recipe for domain mean histogram #2113 --- ...r_common_domain_histogram_domain_mean.yaml | 63 +++++++++++++++++++ 1 file changed, 63 insertions(+) create mode 100644 src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml diff --git a/src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml b/src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml new file mode 100644 index 0000000000..d06388f99f --- /dev/null +++ b/src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml @@ -0,0 +1,63 @@ +category: Surface Spatial Plot +title: "bmc9 Domain mean histogram $VARNAME\n common domain model and radar $NIMROD_LABEL" +description: | + Creates domain mean histogram of the model rainfall rate and + the observed radar rainfall rate. The model field is regridded onto + the radar grid. Both regridded model and radar grids are then masked + so that only data within a domain common to both are used. + +steps: + - operator: read.read_cubes + file_paths: $INPUT_PATHS + model_names: $ALL_LABEL + subarea_type: $SUBAREA_TYPE + subarea_extent: $SUBAREA_EXTENT + constraint: [$VARNAME, "Hourly rain accumulation", "Hourly wts accumulation"] + + - operator: radar_filter.radar_mask_loop + model_field: + operator: filters.filter_multiple_cubes + constraint: + operator: constraints.combine_constraints + varname_constraint: + operator: constraints.generate_var_constraint + varname: $VARNAME + cell_methods_constraint: + operator: constraints.generate_cell_methods_constraint + cell_methods: [] + varname: $VARNAME + pressure_level_constraint: + operator: constraints.generate_level_constraint + coordinate: "pressure" + levels: [] + + nimrod_field: + operator: filters.filter_multiple_cubes + constraint: + operator: constraints.generate_var_constraint + varname: "Hourly rain accumulation" + + nimrod_mask: + operator: filters.generate_mask + mask_field: + operator: filters.filter_multiple_cubes + constraint: + operator: constraints.generate_var_constraint + varname: "Hourly wts accumulation" + condition: 'ge' + value: 11 + + boundary_margin: 16 + outputs: $OUTPUTS + + - operator: radar_filter.match_varname_and_units + + - operator: collapse.collapse + coordinate: [projection_y_coordinate, projection_x_coordinate] + method: MEAN + + - operator: write.write_cube_to_nc + overwrite: True + + - operator: plot.plot_histogram_series + sequence_coordinate: $SEQUENCE From a9e0d1303c174af82fd9f77884d868034dc620fe Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 27 Aug 2026 09:44:47 +0100 Subject: [PATCH 017/105] Radar recipe for histogram sequences #2113 --- ...adar_common_domain_histogram_sequence.yaml | 59 +++++++++++++++++++ 1 file changed, 59 insertions(+) create mode 100644 src/CSET/recipes/surface_fields/radar_common_domain_histogram_sequence.yaml diff --git a/src/CSET/recipes/surface_fields/radar_common_domain_histogram_sequence.yaml b/src/CSET/recipes/surface_fields/radar_common_domain_histogram_sequence.yaml new file mode 100644 index 0000000000..b03359916e --- /dev/null +++ b/src/CSET/recipes/surface_fields/radar_common_domain_histogram_sequence.yaml @@ -0,0 +1,59 @@ +category: Surface Spatial Plot +title: "bmc8b Histogram $VARNAME\n common domain model and radar $NIMROD_LABEL" +description: | + Creates time sequence histograms of the model rainfall rate and + the observed radar rainfall rate. The model field is regridded onto + the radar grid. Both regridded model and radar grids are then masked + so that only data within a domain common to both are used. + +steps: + - operator: read.read_cubes + file_paths: $INPUT_PATHS + model_names: $ALL_LABEL + subarea_type: $SUBAREA_TYPE + subarea_extent: $SUBAREA_EXTENT + constraint: [$VARNAME, "Hourly rain accumulation", "Hourly wts accumulation"] + + - operator: radar_filter.radar_mask_loop + model_field: + operator: filters.filter_multiple_cubes + constraint: + operator: constraints.combine_constraints + varname_constraint: + operator: constraints.generate_var_constraint + varname: $VARNAME + cell_methods_constraint: + operator: constraints.generate_cell_methods_constraint + cell_methods: [] + varname: $VARNAME + pressure_level_constraint: + operator: constraints.generate_level_constraint + coordinate: "pressure" + levels: [] + + nimrod_field: + operator: filters.filter_multiple_cubes + constraint: + operator: constraints.generate_var_constraint + varname: "Hourly rain accumulation" + + nimrod_mask: + operator: filters.generate_mask + mask_field: + operator: filters.filter_multiple_cubes + constraint: + operator: constraints.generate_var_constraint + varname: "Hourly wts accumulation" + condition: 'ge' + value: 11 + + boundary_margin: 16 + outputs: $OUTPUTS + + - operator: radar_filter.match_varname_and_units + + - operator: write.write_cube_to_nc + overwrite: True + + - operator: plot.plot_histogram_series + sequence_coordinate: $SEQUENCE From 6073ff92e8415ffe77ab4a638f06fb8a7f61e173 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 27 Aug 2026 11:18:14 +0100 Subject: [PATCH 018/105] Radar loaders #2113 --- src/CSET/loaders/radar.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index 175260b78e..57b7fd0d67 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -190,8 +190,8 @@ def load(conf: Config): # Histogram domain mean using common domain between # model rainfall and radar observations. - # if conf.HISTOGRAM_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: - if False: + if conf.HISTOGRAM_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: + # if False: # Select the radar source to use. radar_source = select_radar_source([source["id"] for source in accum_radars]) radar_obs_ids = [radar_source] @@ -207,9 +207,10 @@ def load(conf: Config): "ALL_LABEL": combined_names, "MODEL_LABEL": model_names_list, "NIMROD_LABEL": radar_source, - "SEQUENCE": "time", + # "SEQUENCE": "time", "OUTPUTS": "all", - "METHOD": "SEQ", + # "METHOD": "SEQ", + "METHOD": "MEAN", "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", @@ -220,8 +221,8 @@ def load(conf: Config): # Domain mean timeseries using common domain between # model rainfall and radar observations. + # if conf.HISTOGRAM_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: if False: - # if conf.HISTOGRAM_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: # Select the radar source to use. radar_source = select_radar_source([source["id"] for source in accum_radars]) radar_obs_ids = [radar_source] @@ -250,8 +251,8 @@ def load(conf: Config): # Domain timeseries sequence using common domain between # model rainfall and radar observations. - # if False: - if conf.HISTOGRAM_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: + # if conf.HISTOGRAM_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: + if False: # Select the radar source to use. radar_source = select_radar_source([source["id"] for source in accum_radars]) radar_obs_ids = [radar_source] From 7acc815e5f3d953e011c253d22766708759f9097 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 27 Aug 2026 11:23:32 +0100 Subject: [PATCH 019/105] recipe histogram mean #2113 --- .../radar_common_domain_histogram_domain_mean.yaml | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml b/src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml index d06388f99f..f94e701d47 100644 --- a/src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml +++ b/src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml @@ -53,11 +53,12 @@ steps: - operator: radar_filter.match_varname_and_units - operator: collapse.collapse - coordinate: [projection_y_coordinate, projection_x_coordinate] + #coordinate: [projection_y_coordinate, projection_x_coordinate] + coordinate: [time] method: MEAN - operator: write.write_cube_to_nc overwrite: True - operator: plot.plot_histogram_series - sequence_coordinate: $SEQUENCE + #sequence_coordinate: $SEQUENCE From 36dd595435320c20735ee44505bb357c65f180db Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 27 Aug 2026 11:25:40 +0100 Subject: [PATCH 020/105] recipe 2-d radar plots #2113 --- .../generic_surface_spatial_plot_sequence_radar_rainfall.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml b/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml index 66c8606d8e..e6524b0249 100644 --- a/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml +++ b/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml @@ -1,5 +1,5 @@ category: Surface Spatial Plot -title: $RADAR_NAME $VARNAME $SUBAREA_NAME +title: "bmc3 $RADAR_NAME $VARNAME $SUBAREA_NAME" description: Extracts and plots radar rainfall data. steps: From ae5ba6ca2ee6fb056aeb52949fcf1e3f86f1fe0e Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 27 Aug 2026 11:26:24 +0100 Subject: [PATCH 021/105] recipe radar_dev3 plots #2113 --- src/CSET/recipes/surface_fields/radar_dev3.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/recipes/surface_fields/radar_dev3.yaml b/src/CSET/recipes/surface_fields/radar_dev3.yaml index 82683d7e36..4d62b52abf 100644 --- a/src/CSET/recipes/surface_fields/radar_dev3.yaml +++ b/src/CSET/recipes/surface_fields/radar_dev3.yaml @@ -1,5 +1,5 @@ category: Histogram -title: "Histogram masked Nimrod $VARNAME\n$RADAR_NAME" +title: "dev3 Histogram masked Nimrod $VARNAME\n$RADAR_NAME" description: Histograms of Nimrod radar observations # Extracts and plots the probability density of surface `$VARNAME`. It uses # [`plt.hist`](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.hist.html). From 596cace613dfe1e7b7843602e9b63eb9fb6eb37b Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 27 Aug 2026 11:27:15 +0100 Subject: [PATCH 022/105] recipe radar_masked_mean_time_series #2113 --- .../recipes/surface_fields/radar_masked_mean_time_series.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/recipes/surface_fields/radar_masked_mean_time_series.yaml b/src/CSET/recipes/surface_fields/radar_masked_mean_time_series.yaml index 83184b8861..f723f53033 100644 --- a/src/CSET/recipes/surface_fields/radar_masked_mean_time_series.yaml +++ b/src/CSET/recipes/surface_fields/radar_masked_mean_time_series.yaml @@ -1,5 +1,5 @@ category: Histogram -title: "Masked Nimrod mean times series $VARNAME\n$RADAR_NAME" +title: "bmc5 Masked Nimrod mean times series $VARNAME\n$RADAR_NAME" description: Timeseries of masked Nimrod radar observations. steps: From b4390e1f445e021d20a171f82c8dc64fb00a0195 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 27 Aug 2026 11:27:54 +0100 Subject: [PATCH 023/105] recipe radar_masked_time_series #2113 --- src/CSET/recipes/surface_fields/radar_mean_time_series.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/recipes/surface_fields/radar_mean_time_series.yaml b/src/CSET/recipes/surface_fields/radar_mean_time_series.yaml index aaf417f4eb..293dc8c9e8 100644 --- a/src/CSET/recipes/surface_fields/radar_mean_time_series.yaml +++ b/src/CSET/recipes/surface_fields/radar_mean_time_series.yaml @@ -1,5 +1,5 @@ category: Surface Time Series -title: Domain mean $VARNAME time series $SUBAREA_NAME +title: "bmc6 Domain mean $VARNAME time series $SUBAREA_NAME" description: Plots a time series of the domain horizontal mean $VARNAME. steps: From 0b096aeb5d76b0acc1ecb0924ed98daea2c1d56b Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 27 Aug 2026 11:29:57 +0100 Subject: [PATCH 024/105] radar rose-meta.conf #2113 --- src/CSET/cset_workflow/meta/observations/rose-meta.conf | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/CSET/cset_workflow/meta/observations/rose-meta.conf b/src/CSET/cset_workflow/meta/observations/rose-meta.conf index b77b90b7b9..b5367960b8 100644 --- a/src/CSET/cset_workflow/meta/observations/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/observations/rose-meta.conf @@ -159,8 +159,8 @@ trigger=template variables=NIMROD_COMP_XKM: True; [template variables=NIMROD_COMP_XKM] ns=Observations/Radar -title=UK Nimrod radar accumulated rainfall -description=Default UK Nimrod rainfall accumulation +title=UK Nimrod radar accumulated rainfall 5km gridding +description=5km gridded UK Nimrod rainfall accumulation help=Switch to select default gridded UK Nimrod rainfall accumulation. compulsory=true type=python_boolean From 5323901f513edf1b735f28ec2d8914783a082953 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 3 Sep 2026 16:09:00 +0100 Subject: [PATCH 025/105] modified collapse operator #2113 --- src/CSET/operators/collapse.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/src/CSET/operators/collapse.py b/src/CSET/operators/collapse.py index dcd2994844..796309a09e 100644 --- a/src/CSET/operators/collapse.py +++ b/src/CSET/operators/collapse.py @@ -74,6 +74,30 @@ def collapse( if method == "PERCENTILE" and additional_percent is None: raise ValueError("Must specify additional_percent") + print("bmcline77 starting collapse.collapse and value of method is ", method) + + cube_number = 0 + for new_cube in cubes: + print( + "@@@@ collapse.collapse @@@@@@@ cube number ", + cube_number, + "@@@@@@@@@@@@@@@@@@", + ) + print("new_cube.name : ", new_cube.name) + print("new_cube.standard_name : ", new_cube.standard_name) + print("new_cube.long_name : ", new_cube.long_name) + print("new_cube.var_name : ", new_cube.var_name) + print("new_cube.coords(time) : ", new_cube.coord("time")) + print( + "new_cube.coords(forecast_reference_time : ", + new_cube.coord("forecast_reference_time"), + ) + print("new_cube.coords(forecast_period) : ", new_cube.coord("forecast_period")) + print(new_cube) + print("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@") + cube_number = cube_number + 1 + # print(bmc_halt_this) + # Retain only common time points between different models if multiple model inputs. if isinstance(cubes, iris.cube.CubeList) and len(cubes) > 1: logging.debug( @@ -88,6 +112,8 @@ def collapse( if len(cubes) == 0: raise ValueError("No overlapping times detected in input cubes.") + print("bmcline92 collapse.collapse after common time points section") + collapsed_cubes = iris.cube.CubeList([]) with warnings.catch_warnings(): warnings.filterwarnings( @@ -116,6 +142,7 @@ def collapse( collapsed_cubes.append( cube.collapsed(coordinate, getattr(iris.analysis, method)) ) + print("bmcline123 exiting collapse.collapse and value of method is ", method) if len(collapsed_cubes) == 1: return collapsed_cubes[0] else: From 282e4166612490341368e3d0d330a6e5b0f9b37b Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 3 Sep 2026 16:09:48 +0100 Subject: [PATCH 026/105] modified radar_filter operator #2113 --- src/CSET/operators/radar_filter.py | 57 +++++++++++++++++++++++++++++- 1 file changed, 56 insertions(+), 1 deletion(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index e4b0019fcb..41735f1a6f 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -441,13 +441,16 @@ def match_varname_and_units(cubes: iris.cube.Cube | iris.cube.CubeList): if len(cubes) == 1: return cubes - print("@@@@@@@@@@@@@@@@@@@@@ cubes[0] @@@@@@@@@@@@") + print("@@@@@@@@@@@@@@@@@@@@@ base cubes[0] @@@@@@@@@@@@") print(cubes[0]) print("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@") + # TODO when testing completed remove the chop [2:25] [3:46] + # Initialise the list of matched cubes. new_cubelist = iris.cube.CubeList([]) new_cubelist.append(cubes[0]) + # new_cubelist.append(cubes[0][2:45]) # Loop the cubes to match to the first cube. base_cube = cubes[0] @@ -456,24 +459,76 @@ def match_varname_and_units(cubes: iris.cube.Cube | iris.cube.CubeList): print("base_cube.standard_name : ", base_cube.standard_name) print("base_cube.long_name : ", base_cube.long_name) print("base_cube.var_name : ", base_cube.var_name) + print("base_cube.coords(time) : ", base_cube.coord("time")) + print("base_cube.coords(forecast_period) : ", base_cube.coord("forecast_period")) + # print("------- base cube below -----------") + # print(base_cube) print("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@") for cube in cubes[1:]: print("@@@@@@@@@@@@@@@ raw @@@@@@@@@@@@@@@@@@") print(cube) print("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@") new_cube = cube.copy() + # if new_cube.attributes["model_name"] == "Nimrod2km": + # new_cube = new_cube[3:46] + # else: + # new_cube = new_cube[2:45] + + # Match the cube varname. new_cube.rename(base_cube.long_name) new_cube.long_name = base_cube.long_name + new_cube.var_name = base_cube.var_name + + # Match the cube units. new_cube.units = base_cube.units + + # Match the cube forecast_reference_time, a single value e.g. 2023-11-12 00:00:00. + # TODO move this to the read function for radar + forecast_reference_time_use = base_cube.coord("forecast_reference_time") + new_cube.remove_coord("forecast_reference_time") + new_cube.add_aux_coord(forecast_reference_time_use) + # new_cube.replace_coord(base_cube.coord('forecast_reference_time')) + + # Match the cube forecast_period. + # forecast_period_use = base_cube.coord('forecast_period') + # new_cube.remove_coord("forecast_period") + # new_cube.replace_coord(forecast_period_use, data_dims=1) + + # coord_dims = base_cube.coord_dims('forecast_period') + # print("---> coord_dims", coord_dims) + # dim_index = coord_dims[0] if coord_dims else 0 + # print("---> dim_index", dim_index) + # new_coord = base_cube.coord('forecast_period') + # print("---> new_coord", new_coord) + # new_cube.remove_coord("forecast_period") + # new_cube.add_dim_coord(new_coord, dim_index) + + # Match the cube forecast_period - an array points: [ 1., 2., ..., 47., 48.]. + # TODO move this to the read function for radar + coord_dims = base_cube.coord_dims("forecast_period") + print("---> coord_dims", coord_dims) + new_coord = base_cube.coord("forecast_period") + print("---> new_coord", new_coord) + if new_cube.coords("forecast_period"): + new_cube.remove_coord("forecast_period") + new_cube.add_aux_coord(new_coord, coord_dims) + print("@@@@@@@@@@@@@@@ matched @@@@@@@@@@@@@@@@@@") print("new_cube.name : ", new_cube.name) print("new_cube.standard_name : ", new_cube.standard_name) print("new_cube.long_name : ", new_cube.long_name) print("new_cube.var_name : ", new_cube.var_name) + print("new_cube.coords(time) : ", new_cube.coord("time")) + print("new_cube.coords(forecast_period) : ", new_cube.coord("forecast_period")) print(new_cube) print("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@") # Append the matched cube to the output cube list. new_cubelist.append(new_cube) + print( + "---------------->bmc exiting radar_filter.match_varname_and_units len(cubes) is ", + len(cubes), + ) + # print(bmc_halt_this) return new_cubelist From d2aff8b0135202022e9b706a4adfbc5baa8e89e3 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Wed, 9 Sep 2026 09:15:02 +0100 Subject: [PATCH 027/105] Added NIMROD_RADAR_OBS guards in radar loaders #2113 --- src/CSET/loaders/radar.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index 57b7fd0d67..e16f494f6f 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -280,7 +280,7 @@ def load(conf: Config): ) # Radar masking of radar obs based on sea mask. - if conf.SPATIAL_SURFACE_FIELD: + if conf.SPATIAL_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: # field = "Hourly rain accumulation" field = "surface_microphysical_rainfall_rate" yield RawRecipe( @@ -300,7 +300,7 @@ def load(conf: Config): ) # Radar masking of radar obs based on sea mask. - # if conf.SPATIAL_SURFACE_FIELD: + # if conf.SPATIAL_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: if False: field = "Hourly rain accumulation" yield RawRecipe( From 2c7731dc3586e07e1ff3747dda06ca5b35d8ae65 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Wed, 9 Sep 2026 14:21:38 +0100 Subject: [PATCH 028/105] fixed a radar recipe to produce mean histograms over a case study #2113 --- .../radar_common_domain_histogram_domain_mean.yaml | 1 + 1 file changed, 1 insertion(+) diff --git a/src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml b/src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml index f94e701d47..644f5a5554 100644 --- a/src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml +++ b/src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml @@ -61,4 +61,5 @@ steps: overwrite: True - operator: plot.plot_histogram_series + sequence_coordinate: "realization" #sequence_coordinate: $SEQUENCE From 5c0407f2fb1c83ad1ae967da542f02a4184f8a15 Mon Sep 17 00:00:00 2001 From: Sylviabohnenstengel Date: Wed, 9 Sep 2026 16:41:16 +0100 Subject: [PATCH 029/105] wire SEQUENCE into radar histogram recipe. --- src/CSET/loaders/radar.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index e16f494f6f..9223a9c1d6 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -207,7 +207,9 @@ def load(conf: Config): "ALL_LABEL": combined_names, "MODEL_LABEL": model_names_list, "NIMROD_LABEL": radar_source, - # "SEQUENCE": "time", + "SEQUENCE": "time" + if conf.HISTOGRAM_SURFACE_FIELD_SEQUENCE + else "realization", "OUTPUTS": "all", # "METHOD": "SEQ", "METHOD": "MEAN", From 39dcac25d8db816ef491c5a75458ba8006d458df Mon Sep 17 00:00:00 2001 From: Sylviabohnenstengel Date: Wed, 9 Sep 2026 17:49:24 +0100 Subject: [PATCH 030/105] remove time mean of spatial fields for calculation of radar histograms. --- .../radar_common_domain_histogram_domain_mean.yaml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml b/src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml index 644f5a5554..2229807878 100644 --- a/src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml +++ b/src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml @@ -52,10 +52,10 @@ steps: - operator: radar_filter.match_varname_and_units - - operator: collapse.collapse - #coordinate: [projection_y_coordinate, projection_x_coordinate] - coordinate: [time] - method: MEAN + #- operator: collapse.collapse + # #coordinate: [projection_y_coordinate, projection_x_coordinate] + # coordinate: [time] + # method: MEAN - operator: write.write_cube_to_nc overwrite: True From 2d3e72d382c49188ae732abc719d43fe87746429 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 10 Sep 2026 08:15:18 +0100 Subject: [PATCH 031/105] update collapse operator #2113 --- src/CSET/operators/collapse.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/operators/collapse.py b/src/CSET/operators/collapse.py index 796309a09e..534d44081c 100644 --- a/src/CSET/operators/collapse.py +++ b/src/CSET/operators/collapse.py @@ -112,7 +112,7 @@ def collapse( if len(cubes) == 0: raise ValueError("No overlapping times detected in input cubes.") - print("bmcline92 collapse.collapse after common time points section") + # print("bmcline92 collapse.collapse after common time points section") collapsed_cubes = iris.cube.CubeList([]) with warnings.catch_warnings(): From 0901dd3ab3a3955c5e211991fb8e44ebc667faf1 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 10 Sep 2026 08:15:54 +0100 Subject: [PATCH 032/105] update radar_filter operator #2113 --- src/CSET/operators/radar_filter.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index 41735f1a6f..75bca16543 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -466,6 +466,12 @@ def match_varname_and_units(cubes: iris.cube.Cube | iris.cube.CubeList): print("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@") for cube in cubes[1:]: print("@@@@@@@@@@@@@@@ raw @@@@@@@@@@@@@@@@@@") + print("cube.name : ", cube.name) + print("cube.standard_name : ", cube.standard_name) + print("cube.long_name : ", cube.long_name) + print("cube.var_name : ", cube.var_name) + print("cube.coords(time) : ", cube.coord("time")) + print("cube.coords(forecast_period) : ", cube.coord("forecast_period")) print(cube) print("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@") new_cube = cube.copy() From 9870de728402be703d762d41be1599e69f7a71db Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 10 Sep 2026 09:54:55 +0100 Subject: [PATCH 033/105] Added +1 hour offset to the start time of the radar files to fetch #2113 --- .../app/fetch_nimrod/bin/fetch_nimrod.py | 22 +++++++++++++------ 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py b/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py index 59ecc765fe..717885faf1 100755 --- a/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py +++ b/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py @@ -5,7 +5,7 @@ import json import logging import os -from datetime import datetime +from datetime import datetime, timedelta from pathlib import Path import iris @@ -99,10 +99,6 @@ def retrieve_nimrod(): with open(nimrod_met_office, "rt") as fp: # pragma: no cover nimrod_dict = json.load(fp) - # Form the Nimrod start and end dates. - date_start = v["data_time"] - date_end = v["data_time"] + v["forecast_length"] - # Loop over the required Nimrod fields, i.e. 1km 2km or xkm rainfall # accumulation composites or the 5 minute rainfall rate composites. for nimrod_field in v["field"]: @@ -122,9 +118,21 @@ def retrieve_nimrod(): os.makedirs(nimrod_dir_wei, exist_ok=True) logging.info("Cylc-run Nimrod weights directory: %s", nimrod_dir_wei) + # Put +1 hour offset for accumulation radar files as + # the time stamps for these files mark the end of the + # accumulation period rather than the beginning. + radar_offset = timedelta(hours=0.0) + if nimrod_field == "Nimrod_comp_xkm": + radar_offset = timedelta(hours=1.0) + if nimrod_field == "Nimrod_comp_1km": + radar_offset = timedelta(hours=1.0) + if nimrod_field == "Nimrod_comp_2km": + radar_offset = timedelta(hours=1.0) + date_start_offset = v["data_time"] + radar_offset + # Process Nimrod data between the start and end dates. - date_use = date_start - while date_use <= date_end: + date_use = date_start_offset + while date_use <= date_start_offset + v["forecast_length"]: # Load the Nimrod data into an Iris cube. nimrod_obs_exist = "False" nimrod_obs = ( From 20f101d3fe6b87d1c9a2906b9c8186f881a85cf3 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 10 Sep 2026 11:33:28 +0100 Subject: [PATCH 034/105] New button in Rose GUI to select radar histogram outputs #2113 --- .../cset_workflow/meta/observations/rose-meta.conf | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/CSET/cset_workflow/meta/observations/rose-meta.conf b/src/CSET/cset_workflow/meta/observations/rose-meta.conf index b5367960b8..882a743e4c 100644 --- a/src/CSET/cset_workflow/meta/observations/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/observations/rose-meta.conf @@ -205,6 +205,20 @@ compulsory=true type=python_boolean sort-key=nimrod-radar-7 +[template variables=PROCESS_RADAR_HISTOGRAMS] +ns=Observations/Radar +title=Produce Nimrod radar histograms +description=Produce histograms incorporating radar observations. +help=Switch to produce histograms utilising Nimrod radar observations. + The histograms will combine both model rainfall and radar observations. + The model rainfall variable is selected in in the Rose tab Diagnostics/Surface (2D) fields. +compulsory=true +type=python_boolean +sort-key=nimrod-radar-8 + + + + ################################################################ # Cardington UK Research Observations [Observations/Cardington] From 99478b692d6e22e26ce102699844f2b274bbda7c Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 10 Sep 2026 12:39:44 +0100 Subject: [PATCH 035/105] New recipe and loaders for Nimrod radar histograms #2113 --- src/CSET/loaders/radar.py | 46 +++++++------- .../radar_common_domain_histogram.yaml | 61 +++++++++++++++++++ 2 files changed, 86 insertions(+), 21 deletions(-) create mode 100644 src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index 9223a9c1d6..9ed868664a 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -159,9 +159,13 @@ def load(conf: Config): aggregation=False, ) - # Histogram sequence using common domain between - # model rainfall and radar observations. - if conf.HISTOGRAM_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: + # Histogram sequence rainfall using common domain between + # model and Nimrod radar observations. + if ( + conf.HISTOGRAM_SURFACE_FIELD + and conf.NIMROD_RADAR_OBS + and conf.PROCESS_RADAR_HISTOGRAMS + ): # Select the radar source to use. radar_source = select_radar_source([source["id"] for source in accum_radars]) radar_obs_ids = [radar_source] @@ -171,15 +175,15 @@ def load(conf: Config): combined_names = model_names_list + radar_obs_ids + radar_wts_ids combined_ids = model_ids_list + radar_obs_ids + radar_wts_ids yield RawRecipe( - recipe="radar_common_domain_histogram_sequence.yaml", + recipe="radar_common_domain_histogram.yaml", variables={ - "VARNAME": "surface_microphysical_rainfall_rate", + "MODEL_VARNAME": "surface_microphysical_rainfall_rate", + "RADAR_VARNAME": "Hourly rain accumulation", + "RADAR_WTS_VARNAME": "Hourly wts accumulation", "ALL_LABEL": combined_names, - "MODEL_LABEL": model_names_list, - "NIMROD_LABEL": radar_source, "SEQUENCE": "time", "OUTPUTS": "all", - "METHOD": "SEQ", + "TITLE_STRING": "sequence plots", "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", @@ -188,10 +192,13 @@ def load(conf: Config): aggregation=False, ) - # Histogram domain mean using common domain between - # model rainfall and radar observations. - if conf.HISTOGRAM_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: - # if False: + # Histogram case study rainfall using common domain between + # model and Nimrod radar observations. + if ( + conf.HISTOGRAM_SURFACE_FIELD + and conf.NIMROD_RADAR_OBS + and conf.PROCESS_RADAR_HISTOGRAMS + ): # Select the radar source to use. radar_source = select_radar_source([source["id"] for source in accum_radars]) radar_obs_ids = [radar_source] @@ -201,18 +208,15 @@ def load(conf: Config): combined_names = model_names_list + radar_obs_ids + radar_wts_ids combined_ids = model_ids_list + radar_obs_ids + radar_wts_ids yield RawRecipe( - recipe="radar_common_domain_histogram_domain_mean.yaml", + recipe="radar_common_domain_histogram.yaml", variables={ - "VARNAME": "surface_microphysical_rainfall_rate", + "MODEL_VARNAME": "surface_microphysical_rainfall_rate", + "RADAR_VARNAME": "Hourly rain accumulation", + "RADAR_WTS_VARNAME": "Hourly wts accumulation", "ALL_LABEL": combined_names, - "MODEL_LABEL": model_names_list, - "NIMROD_LABEL": radar_source, - "SEQUENCE": "time" - if conf.HISTOGRAM_SURFACE_FIELD_SEQUENCE - else "realization", + "SEQUENCE": "realization", "OUTPUTS": "all", - # "METHOD": "SEQ", - "METHOD": "MEAN", + "TITLE_STRING": "case study", "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", diff --git a/src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml b/src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml new file mode 100644 index 0000000000..129fcae38b --- /dev/null +++ b/src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml @@ -0,0 +1,61 @@ +category: Surface Spatial Plot +title: "bmc9b Histogram $MODEL_VARNAME\n $TITLE_STRING common domain model and Nimrod radar" +description: | + Creates histograms of the model rainfall rate and + the observed radar rainfall rate. The model field is regridded onto + the radar grid. Both regridded model and radar grids are then masked + so that only data within a domain common to both are used. + If $SEQUENCE is set to "time" then time sequence plots are produced. + If $SEQUENCE is set to "realization" then a single case-study mean plot is produced. + +steps: + - operator: read.read_cubes + file_paths: $INPUT_PATHS + model_names: $ALL_LABEL + subarea_type: $SUBAREA_TYPE + subarea_extent: $SUBAREA_EXTENT + constraint: [$MODEL_VARNAME, $RADAR_VARNAME, $RADAR_WTS_VARNAME] + + - operator: radar_filter.radar_mask_loop + model_field: + operator: filters.filter_multiple_cubes + constraint: + operator: constraints.combine_constraints + varname_constraint: + operator: constraints.generate_var_constraint + varname: $MODEL_VARNAME + cell_methods_constraint: + operator: constraints.generate_cell_methods_constraint + cell_methods: [] + varname: $MODEL_VARNAME + pressure_level_constraint: + operator: constraints.generate_level_constraint + coordinate: "pressure" + levels: [] + + nimrod_field: + operator: filters.filter_multiple_cubes + constraint: + operator: constraints.generate_var_constraint + varname: $RADAR_VARNAME + + nimrod_mask: + operator: filters.generate_mask + mask_field: + operator: filters.filter_multiple_cubes + constraint: + operator: constraints.generate_var_constraint + varname: $RADAR_WTS_VARNAME + condition: 'ge' + value: 11 + + boundary_margin: 16 + outputs: $OUTPUTS + + - operator: radar_filter.match_varname_and_units + + - operator: write.write_cube_to_nc + overwrite: True + + - operator: plot.plot_histogram_series + sequence_coordinate: $SEQUENCE From 6544638564ea814a4a003ae1ca41e89295dd4af7 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 10 Sep 2026 15:19:51 +0100 Subject: [PATCH 036/105] Tidied radar histogram recipe #2113 --- .../recipes/surface_fields/radar_common_domain_histogram.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml b/src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml index 129fcae38b..b2c98a1728 100644 --- a/src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml +++ b/src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml @@ -1,5 +1,5 @@ category: Surface Spatial Plot -title: "bmc9b Histogram $MODEL_VARNAME\n $TITLE_STRING common domain model and Nimrod radar" +title: "Histogram $MODEL_VARNAME\n $TITLE_STRING common domain model and Nimrod radar" description: | Creates histograms of the model rainfall rate and the observed radar rainfall rate. The model field is regridded onto From bac16f577c3ba086537875111e0ee1e6e4786641 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 10 Sep 2026 15:23:21 +0100 Subject: [PATCH 037/105] New radar Rose GUI button for PROCESS_RADAR_2D #2113 --- .../cset_workflow/meta/observations/rose-meta.conf | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/CSET/cset_workflow/meta/observations/rose-meta.conf b/src/CSET/cset_workflow/meta/observations/rose-meta.conf index 882a743e4c..99315cdd7a 100644 --- a/src/CSET/cset_workflow/meta/observations/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/observations/rose-meta.conf @@ -205,6 +205,17 @@ compulsory=true type=python_boolean sort-key=nimrod-radar-7 +[template variables=PROCESS_RADAR_2D] +ns=Observations/Radar +title=Produce Nimrod 2D surface plots +description=Produce 2D surface plots of radar observations. +help=Switch to produce 2D surface plots utilising Nimrod radar observations. + Plots are produced of the 2D surface rainfall for each of the selected Nimrod sources. + Plots of the associated Nimrod weights fields are also produced. +compulsory=true +type=python_boolean +sort-key=nimrod-radar-8 + [template variables=PROCESS_RADAR_HISTOGRAMS] ns=Observations/Radar title=Produce Nimrod radar histograms @@ -214,7 +225,7 @@ help=Switch to produce histograms utilising Nimrod radar observations. The model rainfall variable is selected in in the Rose tab Diagnostics/Surface (2D) fields. compulsory=true type=python_boolean -sort-key=nimrod-radar-8 +sort-key=nimrod-radar-9 From 5a6e22849bc3dd07b818c9cd62018cf3c3eb3a8c Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 10 Sep 2026 15:24:35 +0100 Subject: [PATCH 038/105] Tidied recipe to produce 2D surface plots of radar rainfall and weights #2113 --- .../generic_surface_spatial_plot_sequence_radar_rainfall.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml b/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml index e6524b0249..740eb2841e 100644 --- a/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml +++ b/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml @@ -1,5 +1,5 @@ category: Surface Spatial Plot -title: "bmc3 $RADAR_NAME $VARNAME $SUBAREA_NAME" +title: "$RADAR_NAME $VARNAME $SUBAREA_NAME" description: Extracts and plots radar rainfall data. steps: From 5e0aa20ad159dad1ed0724cdd1079e4caeb0deb8 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 10 Sep 2026 15:26:02 +0100 Subject: [PATCH 039/105] Loader to produce 2D surface plots of radar rainfall and weights #2113 --- src/CSET/loaders/radar.py | 157 +++++++++----------------------------- 1 file changed, 37 insertions(+), 120 deletions(-) diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index 9ed868664a..549ea4872a 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -118,54 +118,35 @@ def load(conf: Config): if radar["varname"] == "Hourly wts accumulation" ] - # # Radar masking based on sea mask. - # if conf.SPATIAL_SURFACE_FIELD: - # for field in conf.SURFACE_FIELDS: - # yield RawRecipe( - # recipe="sea_mask_for_surface_domain_mean_time_series.yaml", - # variables={ - # "VARNAME": field, - # "MODEL_NAME": [model["name"] for model in models], - # "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - # "SUBAREA_EXTENT": conf.SUBAREA_EXTENT - # if conf.SELECT_SUBAREA - # else None, - # "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", - # }, - # model_ids=[model["id"] for model in models], - # aggregation=False, - # ) - - # Common domain between model and radar observations. - # if conf.SPATIAL_SURFACE_FIELD: - if False: - # field = "Hourly rain accumulation" - field = "surface_microphysical_rainfall_rate" - yield RawRecipe( - recipe="radar_mask.yaml", - variables={ - "VARNAME": field, - "MODEL_LABEL": "ModelA", - "NIMROD_LABEL": "Nimrod2km", - "MASK_LABEL": "Nimrod2km_wts", - "SEQUENCE": "time", - "OUTPUTS": "all", - "METHOD": "SEQ", - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, - "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", - }, - model_ids=["1", "Nimrod2km", "Nimrod2km_weights"], - aggregation=False, - ) + # Surface (2D) fields for Nimrod radar rainfall. + # + # These loaders produce 2D plots of both the Nimrod + # surface rainfall and the Nimrod weights field. + # + # The different sources of Nimrod rainfall accumulation have + # different spatial grids. So each source requires its own + # recipe to prevent incompatible cubes being created. + if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D: + for radar in radar_sources: + yield RawRecipe( + recipe="generic_surface_spatial_plot_sequence_radar_rainfall.yaml", + model_ids=radar["id"], # -> Becomes $INPUT_PATHS + variables={ + "VARNAME": radar["varname"], + "RADAR_NAME": radar["name"], + "METHOD": "SEQ", + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT + if conf.SELECT_SUBAREA + else None, + "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", + }, + aggregation=False, + ) # Histogram sequence rainfall using common domain between # model and Nimrod radar observations. - if ( - conf.HISTOGRAM_SURFACE_FIELD - and conf.NIMROD_RADAR_OBS - and conf.PROCESS_RADAR_HISTOGRAMS - ): + if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_HISTOGRAMS: # Select the radar source to use. radar_source = select_radar_source([source["id"] for source in accum_radars]) radar_obs_ids = [radar_source] @@ -194,11 +175,7 @@ def load(conf: Config): # Histogram case study rainfall using common domain between # model and Nimrod radar observations. - if ( - conf.HISTOGRAM_SURFACE_FIELD - and conf.NIMROD_RADAR_OBS - and conf.PROCESS_RADAR_HISTOGRAMS - ): + if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_HISTOGRAMS: # Select the radar source to use. radar_source = select_radar_source([source["id"] for source in accum_radars]) radar_obs_ids = [radar_source] @@ -225,55 +202,18 @@ def load(conf: Config): aggregation=False, ) - # Domain mean timeseries using common domain between - # model rainfall and radar observations. - # if conf.HISTOGRAM_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: - if False: - # Select the radar source to use. - radar_source = select_radar_source([source["id"] for source in accum_radars]) - radar_obs_ids = [radar_source] - radar_wts_ids = [radar_source + "_weights"] - model_names_list = [model["name"] for model in models] - model_ids_list = [model["id"] for model in models] - combined_names = model_names_list + radar_obs_ids + radar_wts_ids - combined_ids = model_ids_list + radar_obs_ids + radar_wts_ids - yield RawRecipe( - recipe="radar_common_domain_timeseries_domain_mean.yaml", - variables={ - "VARNAME": "surface_microphysical_rainfall_rate", - "ALL_LABEL": combined_names, - "MODEL_LABEL": model_names_list, - "NIMROD_LABEL": radar_source, - "SEQUENCE": "time", - "OUTPUTS": "all", - "METHOD": "SEQ", - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, - "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", - }, - model_ids=combined_ids, - aggregation=False, - ) - - # Domain timeseries sequence using common domain between - # model rainfall and radar observations. - # if conf.HISTOGRAM_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: + # Common domain between model and radar observations. + # if conf.SPATIAL_SURFACE_FIELD: if False: - # Select the radar source to use. - radar_source = select_radar_source([source["id"] for source in accum_radars]) - radar_obs_ids = [radar_source] - radar_wts_ids = [radar_source + "_weights"] - model_names_list = [model["name"] for model in models] - model_ids_list = [model["id"] for model in models] - combined_names = model_names_list + radar_obs_ids + radar_wts_ids - combined_ids = model_ids_list + radar_obs_ids + radar_wts_ids + # field = "Hourly rain accumulation" + field = "surface_microphysical_rainfall_rate" yield RawRecipe( - recipe="radar_common_domain_timeseries_sequence.yaml", + recipe="radar_mask.yaml", variables={ - "VARNAME": "surface_microphysical_rainfall_rate", - "ALL_LABEL": combined_names, - "MODEL_LABEL": model_names_list, - "NIMROD_LABEL": radar_source, + "VARNAME": field, + "MODEL_LABEL": "ModelA", + "NIMROD_LABEL": "Nimrod2km", + "MASK_LABEL": "Nimrod2km_wts", "SEQUENCE": "time", "OUTPUTS": "all", "METHOD": "SEQ", @@ -281,7 +221,7 @@ def load(conf: Config): "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", }, - model_ids=combined_ids, + model_ids=["1", "Nimrod2km", "Nimrod2km_weights"], aggregation=False, ) @@ -373,29 +313,6 @@ def load(conf: Config): # aggregation=False, # ) - # Surface (2D) fields for Nimrod radar rainfall. - # - # The different sources of Nimrod rainfall accumulation have - # different spatial grids. So each source requires its own - # recipe to prevent incompatible cubes being created. - if conf.SPATIAL_SURFACE_FIELD: - for radar in radar_sources: - yield RawRecipe( - recipe="generic_surface_spatial_plot_sequence_radar_rainfall.yaml", - model_ids=radar["id"], # -> Becomes $INPUT_PATHS - variables={ - "VARNAME": radar["varname"], - "RADAR_NAME": radar["name"], - "METHOD": "SEQ", - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT - if conf.SELECT_SUBAREA - else None, - "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", - }, - aggregation=False, - ) - # Histograms for surface (2D) Nimrod radar hourly accumulated rainfall. # # To get multiple radar sources plotted on the histogram the From b3271357e8871135d68c44c51a3f30b00c259918 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 10 Sep 2026 19:24:37 +0100 Subject: [PATCH 040/105] New recipe for common domain 2D radar plots #2113 --- .../radar_common_domain_2d_plot.yaml | 61 +++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 src/CSET/recipes/surface_fields/radar_common_domain_2d_plot.yaml diff --git a/src/CSET/recipes/surface_fields/radar_common_domain_2d_plot.yaml b/src/CSET/recipes/surface_fields/radar_common_domain_2d_plot.yaml new file mode 100644 index 0000000000..360da20a1a --- /dev/null +++ b/src/CSET/recipes/surface_fields/radar_common_domain_2d_plot.yaml @@ -0,0 +1,61 @@ +category: Surface Spatial Plot +title: "bmc2c $NIMROD_LABEL $RADAR_VARNAME\n common domain $MODEL_LABEL and radar" +description: | + Creates a 2-d surface plot of the model field using the + grid and masking appropriate to the Nimrod file $MASK_LABEL. An + example would be to produce a plot of the LFRic UKV rainfall rate + re-gridded and masked by the Nimrod 2km data. + +steps: + - operator: read.read_cubes + file_paths: $INPUT_PATHS + model_names: [$MODEL_LABEL, $NIMROD_LABEL, $MASK_LABEL] + subarea_type: $SUBAREA_TYPE + subarea_extent: $SUBAREA_EXTENT + constraint: [$MODEL_VARNAME, "Hourly rain accumulation", "Hourly wts accumulation"] + + - operator: radar_filter.radar_mask + model_field: + operator: filters.filter_multiple_cubes + constraint: + operator: constraints.combine_constraints + varname_constraint: + operator: constraints.generate_var_constraint + varname: $MODEL_VARNAME + cell_methods_constraint: + operator: constraints.generate_cell_methods_constraint + cell_methods: [] + varname: $MODEL_VARNAME + pressure_level_constraint: + operator: constraints.generate_level_constraint + coordinate: "pressure" + levels: [] + + nimrod_field: + operator: filters.filter_multiple_cubes + constraint: + operator: constraints.generate_var_constraint + varname: "Hourly rain accumulation" + + nimrod_mask: + operator: filters.generate_mask + mask_field: + operator: filters.filter_multiple_cubes + constraint: + operator: constraints.generate_var_constraint + varname: "Hourly wts accumulation" + condition: 'ge' + value: 11 + + boundary_margin: 16 + outputs: $OUTPUTS + + - operator: collapse.collapse + coordinate: [time] + method: $METHOD + + - operator: write.write_cube_to_nc + overwrite: True + + - operator: plot.spatial_pcolormesh_plot + sequence_coordinate: time From 232795290621d153984f222cd7a6c9eac2c99fee Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Fri, 11 Sep 2026 16:15:17 +0100 Subject: [PATCH 041/105] Recipe for 2D plots of model rainfall over domain common to both model and radar #2113 --- .../model_radar_common_domain_2d_model.yaml | 61 +++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml diff --git a/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml b/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml new file mode 100644 index 0000000000..45d3bb76d0 --- /dev/null +++ b/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml @@ -0,0 +1,61 @@ +category: Surface Spatial Plot +title: "bmc2f $MODEL_LABEL $MODEL_VARNAME\n common domain $MODEL_LABEL and $RADAR_LABEL" +description: | + Creates a 2-d surface plot of the model rainfall field over the + domain common to both the model and the radar observation network. + To achieve this, the model field is regridded onto the radar grid, + before being masked by the common domain. + +steps: + - operator: read.read_cubes + file_paths: $INPUT_PATHS + model_names: [$MODEL_LABEL, $RADAR_LABEL, $MASK_LABEL] + subarea_type: $SUBAREA_TYPE + subarea_extent: $SUBAREA_EXTENT + constraint: [$MODEL_VARNAME, "Hourly rain accumulation", "Hourly wts accumulation"] + + - operator: radar_filter.radar_mask + model_field: + operator: filters.filter_multiple_cubes + constraint: + operator: constraints.combine_constraints + varname_constraint: + operator: constraints.generate_var_constraint + varname: $MODEL_VARNAME + cell_methods_constraint: + operator: constraints.generate_cell_methods_constraint + cell_methods: [] + varname: $MODEL_VARNAME + pressure_level_constraint: + operator: constraints.generate_level_constraint + coordinate: "pressure" + levels: [] + + nimrod_field: + operator: filters.filter_multiple_cubes + constraint: + operator: constraints.generate_var_constraint + varname: "Hourly rain accumulation" + + nimrod_mask: + operator: filters.generate_mask + mask_field: + operator: filters.filter_multiple_cubes + constraint: + operator: constraints.generate_var_constraint + varname: "Hourly wts accumulation" + condition: 'ge' + value: 11 + + boundary_margin: 16 + outputs: "model" + + - operator: collapse.collapse + coordinate: [time] + method: $METHOD + + - operator: write.write_cube_to_nc + overwrite: True + + - operator: plot.spatial_pcolormesh_plot + sequence_coordinate: time From ae2380d86303bf791a0f90e98a05c59e83e5abce Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Fri, 11 Sep 2026 16:19:02 +0100 Subject: [PATCH 042/105] Recipe for 2D plots of radar rainfall over domain common to both model and radar #2113 --- .../model_radar_common_domain_2d_radar.yaml | 61 +++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml diff --git a/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml b/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml new file mode 100644 index 0000000000..70d9d0c736 --- /dev/null +++ b/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml @@ -0,0 +1,61 @@ +category: Surface Spatial Plot +title: "bmc2c $RADAR_LABEL $RADAR_VARNAME\n common domain $MODEL_LABEL and $RADAR_LABEL" +description: | + Creates a 2-d surface plot of the observed radar rain accumualation + field over the domain common to both the model and the radar observation network. + To achieve this, the model field is regridded onto the radar grid, enabling the + mapping of the common domain that is then used to mask the radar field. + +steps: + - operator: read.read_cubes + file_paths: $INPUT_PATHS + model_names: [$MODEL_LABEL, $RADAR_LABEL, $MASK_LABEL] + subarea_type: $SUBAREA_TYPE + subarea_extent: $SUBAREA_EXTENT + constraint: [$MODEL_VARNAME, "Hourly rain accumulation", "Hourly wts accumulation"] + + - operator: radar_filter.radar_mask + model_field: + operator: filters.filter_multiple_cubes + constraint: + operator: constraints.combine_constraints + varname_constraint: + operator: constraints.generate_var_constraint + varname: $MODEL_VARNAME + cell_methods_constraint: + operator: constraints.generate_cell_methods_constraint + cell_methods: [] + varname: $MODEL_VARNAME + pressure_level_constraint: + operator: constraints.generate_level_constraint + coordinate: "pressure" + levels: [] + + nimrod_field: + operator: filters.filter_multiple_cubes + constraint: + operator: constraints.generate_var_constraint + varname: "Hourly rain accumulation" + + nimrod_mask: + operator: filters.generate_mask + mask_field: + operator: filters.filter_multiple_cubes + constraint: + operator: constraints.generate_var_constraint + varname: "Hourly wts accumulation" + condition: 'ge' + value: 11 + + boundary_margin: 16 + outputs: "radar" + + - operator: collapse.collapse + coordinate: [time] + method: $METHOD + + - operator: write.write_cube_to_nc + overwrite: True + + - operator: plot.spatial_pcolormesh_plot + sequence_coordinate: time From 24701dc012252ba77bdbdbfcc2d357f268df7025 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Mon, 14 Sep 2026 14:56:47 +0100 Subject: [PATCH 043/105] fetch_nimrod.py now retrieves the correct number of radar files #2113 --- src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py b/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py index ea14d3d987..8386664348 100755 --- a/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py +++ b/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py @@ -134,7 +134,7 @@ def retrieve_nimrod(): # Process Nimrod data between the start and end dates. date_use = date_start_offset - while date_use <= date_start_offset + v["forecast_length"]: + while date_use < date_start_offset + v["forecast_length"]: # Load the Nimrod data into an Iris cube. nimrod_obs_exist = "False" nimrod_obs = ( From e388bcc98b9216c39a347a0049adf80a19695455 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Mon, 14 Sep 2026 16:29:51 +0100 Subject: [PATCH 044/105] Put strict_title into ploy.py to turn off prepending plot tiltes with model_name #2113 --- src/CSET/operators/plot.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index 80eb88a832..6497881e99 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -1871,6 +1871,7 @@ def _spatial_plot( filename: str | None, sequence_coordinate: str, stamp_coordinate: str, + strict_title: bool = False, overlay_cube: iris.cube.Cube | None = None, contour_cube: iris.cube.Cube | None = None, point_cube: iris.cube.Cube | None = None, @@ -1966,6 +1967,8 @@ def _spatial_plot( model_name = cube.attributes["model_name"] else: model_name = None + if strict_title: + model_name = None plot_title, plot_filename = _set_title_and_filename( seq_coord, nplot, recipe_title, filename, model_name=model_name @@ -2054,6 +2057,7 @@ def spatial_pcolormesh_plot( filename: str | None = None, sequence_coordinate: str = "time", stamp_coordinate: str = "realization", + strict_title: bool = False, **kwargs, ) -> iris.cube.Cube: """Plot a spatial variable onto a map from a 2D, 3D, or 4D cube. @@ -2100,6 +2104,7 @@ def spatial_pcolormesh_plot( filename, sequence_coordinate, stamp_coordinate, + strict_title, **kwargs, ) elif isinstance(cubes, iris.cube.Cube): @@ -2109,6 +2114,7 @@ def spatial_pcolormesh_plot( filename, sequence_coordinate, stamp_coordinate, + strict_title, **kwargs, ) return cubes From 141c2463967343d82ff0473cccb7d5a745beca40 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Mon, 14 Sep 2026 16:30:49 +0100 Subject: [PATCH 045/105] Put strict_title into a recipe to turn off prepending plot title with model_name #2113 --- .../surface_fields/model_radar_common_domain_2d_model.yaml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml b/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml index 45d3bb76d0..63e26d8873 100644 --- a/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml +++ b/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml @@ -1,5 +1,5 @@ category: Surface Spatial Plot -title: "bmc2f $MODEL_LABEL $MODEL_VARNAME\n common domain $MODEL_LABEL and $RADAR_LABEL" +title: "$MODEL_LABEL $MODEL_VARNAME\n common domain $MODEL_LABEL and $RADAR_LABEL" description: | Creates a 2-d surface plot of the model rainfall field over the domain common to both the model and the radar observation network. @@ -58,4 +58,5 @@ steps: overwrite: True - operator: plot.spatial_pcolormesh_plot + strict_title: True sequence_coordinate: time From 145aab7b4fd8119b95bcb1dbc325a7cdfb662445 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Mon, 14 Sep 2026 16:31:18 +0100 Subject: [PATCH 046/105] Put strict_title into a recipe to turn off prepending plot title with model_name #2113 --- .../surface_fields/model_radar_common_domain_2d_radar.yaml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml b/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml index 70d9d0c736..fa5291075e 100644 --- a/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml +++ b/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml @@ -1,5 +1,5 @@ category: Surface Spatial Plot -title: "bmc2c $RADAR_LABEL $RADAR_VARNAME\n common domain $MODEL_LABEL and $RADAR_LABEL" +title: "$RADAR_LABEL $RADAR_VARNAME\n common domain $MODEL_LABEL and $RADAR_LABEL" description: | Creates a 2-d surface plot of the observed radar rain accumualation field over the domain common to both the model and the radar observation network. @@ -58,4 +58,5 @@ steps: overwrite: True - operator: plot.spatial_pcolormesh_plot + strict_title: True sequence_coordinate: time From 175101bb6ec3e71707be17e06bf130cabf0815b6 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Tue, 15 Sep 2026 15:20:38 +0100 Subject: [PATCH 047/105] Put strict_title into doc strings in ploy.py #2113 --- src/CSET/operators/plot.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index 6497881e99..ed07d736da 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -1905,6 +1905,9 @@ def _spatial_plot( stamp_coordinate: str Coordinate about which to plot postage stamp plots. Defaults to ``"realization"``. + strict_title: bool, optional + Logical switch that if set to True will ensure that the MODEL_NAME + string is not prepended to the plot title. The default is False. overlay_cube: Cube | None, optional Optional 2 dimensional (lat and lon) Cube of data to overplot on top of base cube contour_cube: Cube | None, optional @@ -2085,6 +2088,9 @@ def spatial_pcolormesh_plot( stamp_coordinate: str, optional Coordinate about which to plot postage stamp plots. Defaults to ``"realization"``. + strict_title: bool, optional + Logical switch that if set to True will ensure that the MODEL_NAME + string is not prepended to the plot title. The default is False. Returns ------- From b5b7b806500a937405feafe0f00213d78eadc65f Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 09:16:20 +0100 Subject: [PATCH 048/105] Removed redundant code from radar_filter.radar_apply_mask #2113 --- src/CSET/operators/radar_filter.py | 27 ++++++++------------------- 1 file changed, 8 insertions(+), 19 deletions(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index 75bca16543..a0866fad41 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -241,27 +241,14 @@ def radar_apply_mask( for M, F in zip(iter_maybe(mask), iter_maybe(original_field), strict=True): masked_field = F.copy() - # TODO set the model perimeter to NaN as these gridpoints contain no data. + # Set the model perimeter to NaN as these gridpoints contain no useful data. # c.f. boundary_margin in regrid.py - # Get bounds - # Get axis - # lat, lon = M.coord(y_coord), M.coord(x_coord) - # lat_min, lon_min = lat.points.min(), lon.points.min() - # lat_max, lon_max = lat.points.max(), lon.points.max() margin_width = boundary_margin if margin_width > 0: masked_field.data[:, -margin_width - 1 :, :] = np.nan masked_field.data[:, :, -margin_width - 1 :] = np.nan masked_field.data[:, :, 0:margin_width] = np.nan masked_field.data[:, 0:margin_width, :] = np.nan - # masked_field.data[:, -margin_width-1:-1, :] = 3.0 - # masked_field.data[:, :, -margin_width-1:-1] = 1.5 - # masked_field.data[:, :, 0:margin_width] = 50.0 - # masked_field.data[:, 0:margin_width, :] = 30.0 - # masked_field.data[ masked_field.data < 0.01 ] = 200.0 - # masked_field.data[ masked_field.data == np.nan ] = 20.0 - - # TODO must mask the radar obs using the model domain # If the field and mask are on different grids, then regrid the field. if M[0].shape != masked_field[0].shape: @@ -291,7 +278,7 @@ def radar_mask( nimrod_field: iris.cube.Cube | iris.cube.CubeList, nimrod_mask: iris.cube.Cube | iris.cube.CubeList, boundary_margin: int = 8, - outputs: str = "Nimrod", + outputs: str = "radar", ) -> iris.cube.Cube | iris.cube.CubeList: """Apply a mask to given data as a masked array.""" # Create an empty cubelist to hold the filtered fields. @@ -334,10 +321,12 @@ def radar_mask( # the histogram or time series plotting operators. print(" test point bmc1") - if outputs == "Nimrod": + if outputs == "radar": + print(" test point bmc1a") filtered_fields.append(filtered_radar.merge_cube()) print(" test point bmc2") if outputs == "model": + print(" test point bmc2a") filtered_fields.append(filtered_model.merge_cube()) print(" test point bmc3") if outputs == "all": @@ -366,7 +355,7 @@ def radar_mask_loop( nimrod_field: iris.cube.Cube | iris.cube.CubeList, nimrod_mask: iris.cube.Cube | iris.cube.CubeList, boundary_margin: int = 8, - outputs: str = "Nimrod", + outputs: str = "radar", ) -> iris.cube.Cube | iris.cube.CubeList: """Find common domains between a list of models and radar observations.""" # Create an empty cubelist to hold the filtered fields. @@ -400,7 +389,7 @@ def radar_mask_loop( use_nimrod_field, use_nimrod_mask, boundary_margin=boundary_margin, - outputs="Nimrod", + outputs="radar", ) else: filtered_radar = radar_mask( @@ -408,7 +397,7 @@ def radar_mask_loop( use_nimrod_field, use_nimrod_mask, boundary_margin=boundary_margin, - outputs="Nimrod", + outputs="radar", ) filtered_cubes.append(filtered_radar) From f04834810402df328cafa016dfe7e683e94a0f9f Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 10:29:56 +0100 Subject: [PATCH 049/105] Radar loaders updated #2113 --- src/CSET/loaders/radar.py | 74 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 71 insertions(+), 3 deletions(-) diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index 549ea4872a..394b113034 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -144,6 +144,70 @@ def load(conf: Config): aggregation=False, ) + # Radar 2D plots using a common domain between model and radar network. + if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D: + # field = "Hourly rain accumulation" + field = "surface_microphysical_rainfall_rate" + radar_source = select_radar_source([source["id"] for source in accum_radars]) + radar_obs_ids = [radar_source] + radar_wts_ids = [radar_source + "_weights"] + model_names_list = [model["name"] for model in models] + model_ids_list = [model["id"] for model in models] + yield RawRecipe( + recipe="model_radar_common_domain_2d_radar.yaml", + variables={ + "MODEL_VARNAME": "surface_microphysical_rainfall_rate", + "RADAR_VARNAME": "Hourly rain accumulation", + "MODEL_LABEL": model_names_list[0], + # "MODEL_LABEL": "ModelA", + # "RADAR_LABEL": "Nimrod2km", + # "MASK_LABEL": "Nimrod2km_wts", + "RADAR_LABEL": radar_obs_ids[0], + "MASK_LABEL": radar_wts_ids[0], + "METHOD": "SEQ", + # "OUTPUTS": "radar", + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, + "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", + }, + model_ids=[model_ids_list[0], "Nimrod2km", "Nimrod2km_weights"], + # model_ids=["1", "Nimrod2km", "Nimrod2km_weights"], + # model_ids=["ModelA", "Nimrod2km", "Nimrod2km_weights"], + aggregation=False, + ) + + # Radar 2D plots using a common domain between model and radar network. + if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D: + # field = "Hourly rain accumulation" + field = "surface_microphysical_rainfall_rate" + radar_source = select_radar_source([source["id"] for source in accum_radars]) + radar_obs_ids = [radar_source] + radar_wts_ids = [radar_source + "_weights"] + model_names_list = [model["name"] for model in models] + model_ids_list = [model["id"] for model in models] + yield RawRecipe( + recipe="model_radar_common_domain_2d_model.yaml", + variables={ + "MODEL_VARNAME": "surface_microphysical_rainfall_rate", + "RADAR_VARNAME": "Hourly rain accumulation", + "MODEL_LABEL": model_names_list[0], + # "MODEL_LABEL": "ModelA", + # "RADAR_LABEL": "Nimrod2km", + # "MASK_LABEL": "Nimrod2km_wts", + "RADAR_LABEL": radar_obs_ids[0], + "MASK_LABEL": radar_wts_ids[0], + "METHOD": "SEQ", + # "OUTPUTS": "model", + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, + "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", + }, + model_ids=[model_ids_list[0], "Nimrod2km", "Nimrod2km_weights"], + # model_ids=["1", "Nimrod2km", "Nimrod2km_weights"], + # model_ids=["ModelA", "Nimrod2km", "Nimrod2km_weights"], + aggregation=False, + ) + # Histogram sequence rainfall using common domain between # model and Nimrod radar observations. if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_HISTOGRAMS: @@ -202,8 +266,10 @@ def load(conf: Config): aggregation=False, ) + ################################################################################# + # Common domain between model and radar observations. - # if conf.SPATIAL_SURFACE_FIELD: + # if conf.SPATIAL_SURFACE_FIELD and False: if False: # field = "Hourly rain accumulation" field = "surface_microphysical_rainfall_rate" @@ -226,7 +292,8 @@ def load(conf: Config): ) # Radar masking of radar obs based on sea mask. - if conf.SPATIAL_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: + # if conf.SPATIAL_SURFACE_FIELD and conf.NIMROD_RADAR_OBS and False: + if False: # field = "Hourly rain accumulation" field = "surface_microphysical_rainfall_rate" yield RawRecipe( @@ -265,7 +332,8 @@ def load(conf: Config): ) # Radar masking of model rainfall based on sea mask. - if conf.SPATIAL_SURFACE_FIELD: + # if conf.SPATIAL_SURFACE_FIELD and False: + if False: field = "surface_microphysical_rainfall_rate" yield RawRecipe( recipe="radar_mask_model.yaml", From 26bb073d13645008215cfd2c7b5a87ce03b30222 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 10:31:21 +0100 Subject: [PATCH 050/105] Spare recipe #2113 --- .../recipes/surface_fields/radar_common_domain_2d_plot.yaml | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/CSET/recipes/surface_fields/radar_common_domain_2d_plot.yaml b/src/CSET/recipes/surface_fields/radar_common_domain_2d_plot.yaml index 360da20a1a..7cb809bfab 100644 --- a/src/CSET/recipes/surface_fields/radar_common_domain_2d_plot.yaml +++ b/src/CSET/recipes/surface_fields/radar_common_domain_2d_plot.yaml @@ -1,5 +1,6 @@ category: Surface Spatial Plot -title: "bmc2c $NIMROD_LABEL $RADAR_VARNAME\n common domain $MODEL_LABEL and radar" +title: "bmc2d" +#title: "bmc2c $RADAR_VARNAME\n common domain $MODEL_LABEL and radar" description: | Creates a 2-d surface plot of the model field using the grid and masking appropriate to the Nimrod file $MASK_LABEL. An @@ -9,7 +10,7 @@ description: | steps: - operator: read.read_cubes file_paths: $INPUT_PATHS - model_names: [$MODEL_LABEL, $NIMROD_LABEL, $MASK_LABEL] + model_names: [$MODEL_LABEL, $RADAR_LABEL, $MASK_LABEL] subarea_type: $SUBAREA_TYPE subarea_extent: $SUBAREA_EXTENT constraint: [$MODEL_VARNAME, "Hourly rain accumulation", "Hourly wts accumulation"] From d73552eaf565d1cc3246c0edef322e0932a5359d Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 13:44:21 +0100 Subject: [PATCH 051/105] Updated recipe description for generic_surface_spatial_plot_sequence_radar_rainfall.yaml #2113 --- ...ric_surface_spatial_plot_sequence_radar_rainfall.yaml | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml b/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml index 740eb2841e..77deba33f5 100644 --- a/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml +++ b/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml @@ -1,6 +1,13 @@ category: Surface Spatial Plot title: "$RADAR_NAME $VARNAME $SUBAREA_NAME" -description: Extracts and plots radar rainfall data. +description: | + Extracts and plots radar rainfall data. + ATTENTION: the radar data can be 1 hour accumulations. Because showers + typically advect, on a 2D plot of accumulated rainfall, showers show + as streaks i.e. lower intensity rainfall over a greater area than a + plot of instantaneous rainfall rate. This results in accumulated rainfall + being damped in intensity compared to rainfall rates so caution is + required when comparing model rainfall rate with observed rainfall accumulation. steps: - operator: read.read_cube From 0bcacf81eb203bf5eaf6ed80c9e4cbad44f3fea4 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 13:45:09 +0100 Subject: [PATCH 052/105] Updated recipe description for model_radar_common_domain_2d_model.yaml #2113 --- .../surface_fields/model_radar_common_domain_2d_model.yaml | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml b/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml index 63e26d8873..3def2fe945 100644 --- a/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml +++ b/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml @@ -5,6 +5,12 @@ description: | domain common to both the model and the radar observation network. To achieve this, the model field is regridded onto the radar grid, before being masked by the common domain. + ATTENTION: the radar data are 1 hour accumulations. Because showers + typically advect, on a 2D plot of accumulated rainfall, showers show + as streaks i.e. lower intensity rainfall over a greater area than a + plot of instantaneous rainfall rate. This results in accumulated rainfall + being damped in intensity compared to rainfall rates so caution is + required when comparing model rainfall rate with observed rainfall accumulation. steps: - operator: read.read_cubes From d0c00fe416a59a2e80e47afec45ccad05fbe3634 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 13:45:24 +0100 Subject: [PATCH 053/105] Updated recipe description for model_radar_common_domain_2d_radar.yaml #2113 --- .../surface_fields/model_radar_common_domain_2d_radar.yaml | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml b/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml index fa5291075e..c5a29c3311 100644 --- a/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml +++ b/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml @@ -5,6 +5,12 @@ description: | field over the domain common to both the model and the radar observation network. To achieve this, the model field is regridded onto the radar grid, enabling the mapping of the common domain that is then used to mask the radar field. + ATTENTION: the radar data are 1 hour accumulations. Because showers + typically advect, on a 2D plot of accumulated rainfall, showers show + as streaks i.e. lower intensity rainfall over a greater area than a + plot of instantaneous rainfall rate. This results in accumulated rainfall + being damped in intensity compared to rainfall rates so caution is + required when comparing model rainfall rate with observed rainfall accumulation. steps: - operator: read.read_cubes From c5683f36fb1f4912d323ad16e30da44356b9f957 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 13:46:13 +0100 Subject: [PATCH 054/105] Updated recipe description for radar_common_domain_histogram.yaml #2113 --- .../radar_common_domain_histogram.yaml | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml b/src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml index b2c98a1728..eb997e4b9b 100644 --- a/src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml +++ b/src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml @@ -2,11 +2,18 @@ category: Surface Spatial Plot title: "Histogram $MODEL_VARNAME\n $TITLE_STRING common domain model and Nimrod radar" description: | Creates histograms of the model rainfall rate and - the observed radar rainfall rate. The model field is regridded onto - the radar grid. Both regridded model and radar grids are then masked - so that only data within a domain common to both are used. + the observed radar rainfall accumulation. + The model field is regridded onto the radar grid. Both regridded model + and radar grids are then masked so that only data within a domain common + to both are used. If $SEQUENCE is set to "time" then time sequence plots are produced. If $SEQUENCE is set to "realization" then a single case-study mean plot is produced. + ATTENTION: the radar data are 1 hour accumulations. Because showers + typically advect, on a 2D plot of accumulated rainfall, showers show + as streaks i.e. lower intensity rainfall over a greater area than a + plot of instantaneous rainfall rate. This results in accumulated rainfall + being damped in intensity compared to rainfall rates so caution is + required when comparing model rainfall rate with observed rainfall accumulation. steps: - operator: read.read_cubes From 0e33a0cfacf302eaeafc69fa141462add82b143b Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 13:50:36 +0100 Subject: [PATCH 055/105] Updated radar loaders #2113 --- src/CSET/loaders/radar.py | 285 ++------------------------------------ 1 file changed, 13 insertions(+), 272 deletions(-) diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index 394b113034..9b08952e09 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -111,13 +111,6 @@ def load(conf: Config): if radar["varname"] == "Hourly rain accumulation" ] - # Form the list of accumulated hourly weights for Nimrod radar sources. - wts_radars = [ - radar - for radar in radar_sources - if radar["varname"] == "Hourly wts accumulation" - ] - # Surface (2D) fields for Nimrod radar rainfall. # # These loaders produce 2D plots of both the Nimrod @@ -126,7 +119,7 @@ def load(conf: Config): # The different sources of Nimrod rainfall accumulation have # different spatial grids. So each source requires its own # recipe to prevent incompatible cubes being created. - if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D: + if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D and (len(accum_radars) > 0): for radar in radar_sources: yield RawRecipe( recipe="generic_surface_spatial_plot_sequence_radar_rainfall.yaml", @@ -145,9 +138,8 @@ def load(conf: Config): ) # Radar 2D plots using a common domain between model and radar network. - if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D: + if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D and (len(accum_radars) > 0): # field = "Hourly rain accumulation" - field = "surface_microphysical_rainfall_rate" radar_source = select_radar_source([source["id"] for source in accum_radars]) radar_obs_ids = [radar_source] radar_wts_ids = [radar_source + "_weights"] @@ -177,9 +169,8 @@ def load(conf: Config): ) # Radar 2D plots using a common domain between model and radar network. - if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D: + if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D and (len(accum_radars) > 0): # field = "Hourly rain accumulation" - field = "surface_microphysical_rainfall_rate" radar_source = select_radar_source([source["id"] for source in accum_radars]) radar_obs_ids = [radar_source] radar_wts_ids = [radar_source + "_weights"] @@ -210,7 +201,11 @@ def load(conf: Config): # Histogram sequence rainfall using common domain between # model and Nimrod radar observations. - if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_HISTOGRAMS: + if ( + conf.NIMROD_RADAR_OBS + and conf.PROCESS_RADAR_HISTOGRAMS + and (len(accum_radars) > 0) + ): # Select the radar source to use. radar_source = select_radar_source([source["id"] for source in accum_radars]) radar_obs_ids = [radar_source] @@ -239,7 +234,11 @@ def load(conf: Config): # Histogram case study rainfall using common domain between # model and Nimrod radar observations. - if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_HISTOGRAMS: + if ( + conf.NIMROD_RADAR_OBS + and conf.PROCESS_RADAR_HISTOGRAMS + and (len(accum_radars) > 0) + ): # Select the radar source to use. radar_source = select_radar_source([source["id"] for source in accum_radars]) radar_obs_ids = [radar_source] @@ -265,261 +264,3 @@ def load(conf: Config): model_ids=combined_ids, aggregation=False, ) - - ################################################################################# - - # Common domain between model and radar observations. - # if conf.SPATIAL_SURFACE_FIELD and False: - if False: - # field = "Hourly rain accumulation" - field = "surface_microphysical_rainfall_rate" - yield RawRecipe( - recipe="radar_mask.yaml", - variables={ - "VARNAME": field, - "MODEL_LABEL": "ModelA", - "NIMROD_LABEL": "Nimrod2km", - "MASK_LABEL": "Nimrod2km_wts", - "SEQUENCE": "time", - "OUTPUTS": "all", - "METHOD": "SEQ", - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, - "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", - }, - model_ids=["1", "Nimrod2km", "Nimrod2km_weights"], - aggregation=False, - ) - - # Radar masking of radar obs based on sea mask. - # if conf.SPATIAL_SURFACE_FIELD and conf.NIMROD_RADAR_OBS and False: - if False: - # field = "Hourly rain accumulation" - field = "surface_microphysical_rainfall_rate" - yield RawRecipe( - recipe="radar_mask_model_obs.yaml", - variables={ - "VARNAME": field, - "MODEL_LABEL": "ModelA", - "NIMROD_LABEL": "Nimrod2km", - "MASK_LABEL": "Nimrod2km_wts", - "METHOD": "SEQ", - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, - "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", - }, - model_ids=["1", "Nimrod2km", "Nimrod2km_weights"], - aggregation=False, - ) - - # Radar masking of radar obs based on sea mask. - # if conf.SPATIAL_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: - if False: - field = "Hourly rain accumulation" - yield RawRecipe( - recipe="radar_mask_model.yaml", - variables={ - "VARNAME": field, - "MODEL_LABEL": "Nimrod2km", - "MASK_LABEL": "Nimrod2km", - "METHOD": "SEQ", - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, - "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", - }, - model_ids=["Nimrod2km", "Nimrod2km_weights"], - aggregation=False, - ) - - # Radar masking of model rainfall based on sea mask. - # if conf.SPATIAL_SURFACE_FIELD and False: - if False: - field = "surface_microphysical_rainfall_rate" - yield RawRecipe( - recipe="radar_mask_model.yaml", - variables={ - "VARNAME": field, - "MODEL_LABEL": "ModelA", - "MASK_LABEL": "Nimrod2km", - "METHOD": "SEQ", - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, - "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", - }, - model_ids=["1", "Nimrod2km_weights"], - aggregation=False, - ) - - # Surface (2D) fields for model rainfall masked by Nimrod radar. - # - # The different sources of Nimrod rainfall accumulation have - # different spatial grids. So each source requires its own - # recipe to prevent incompatible cubes being created. - # if conf.SPATIAL_SURFACE_FIELD: - # radar_source = ["Nimrod_2km"] - # for radar in radar_source: - # model_labels = [model["id"] for model in models] - # radar_label = ["Nimrod2km_weights"] - # combined_ids = [model_labels[0]] + radar_label - # print("Combined ids is: ", combined_ids) - # yield RawRecipe( - # recipe="radar_plot_sequence_rainfall.yaml", - ## model_ids=radar["id"], # -> Becomes $INPUT_PATHS - # model_ids=combined_ids, - # variables={ - ## "VARNAME": radar["varname"], - ## "RADAR_NAME": radar["name"], - # "RADAR_NAME": "Nimrod_2km_weights", - ## "MODEL_NAME": [model["name"] for model in models], - # "MODEL_NAME": "ModelA", - # "METHOD": "SEQ", - # "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - # "SUBAREA_EXTENT": conf.SUBAREA_EXTENT - # if conf.SELECT_SUBAREA - # else None, - # }, - # aggregation=False, - # ) - - # Histograms for surface (2D) Nimrod radar hourly accumulated rainfall. - # - # To get multiple radar sources plotted on the histogram the - # recipe must be done by passing lists of the radar_ids and - # the radar_names. As this is a multiline plot, all radar sources - # share the same radar variable name. - if conf.HISTOGRAM_SURFACE_FIELD and accum_radars: - yield RawRecipe( - recipe="generic_surface_histogram_series.yaml", - # model_ids -> Becomes $INPUT_PATHS - model_ids=[radar["id"] for radar in accum_radars], - variables={ - "VARNAME": next(radar["varname"] for radar in accum_radars), - "MODEL_NAME": [radar["name"] for radar in accum_radars], - "SEQUENCE": "time" - if conf.HISTOGRAM_SURFACE_FIELD_SEQUENCE - else "realization", - "SUBAREA_NAME": "", - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, - }, - aggregation=False, - ) - - # Histograms for surface (2D) Nimrod radar hourly accumulated rainfall. - # - # The histograms are produced after the rainfall obs have been masked using - # the associated Nimrod weights file. - # - # To get multiple radar sources plotted on the histogram the - # recipe must be done by passing lists of the radar_ids and - # the radar_names. As this is a multiline plot, all radar sources - # share the same radar variable name. - if conf.HISTOGRAM_SURFACE_FIELD: - radar_obs_ids = [radar["id"] for radar in accum_radars] - radar_wts_ids = [radar["id"] for radar in wts_radars] - combined_ids = radar_obs_ids + radar_wts_ids - print(" combined_ids: ", combined_ids) - yield RawRecipe( - recipe="radar_dev3.yaml", - # model_ids -> Becomes $INPUT_PATHS - # model_ids=[ radar_obs_ids, radar_wts_ids], - model_ids=combined_ids, - variables={ - "VARNAME": next(radar["varname"] for radar in accum_radars), - "ALL_NAME": combined_ids, - "RADAR_NAME": [radar["id"] for radar in accum_radars], - "WEIGHTS_NAME": [radar["id"] for radar in wts_radars], - "SEQUENCE": "time" - if conf.HISTOGRAM_SURFACE_FIELD_SEQUENCE - else "realization", - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, - "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", - }, - aggregation=False, - ) - - # Histograms for surface (2D) Nimrod radar 5 minute rainfall rate. - if conf.HISTOGRAM_SURFACE_FIELD and conf.NIMROD_COMP_5MIN: - yield RawRecipe( - recipe="generic_surface_histogram_series.yaml", - # model_ids -> Becomes $INPUT_PATHS - model_ids="Nimrod5min", - variables={ - "VARNAME": "Rainfall rate Composite", - "MODEL_NAME": "Nimrod_5min", - "SEQUENCE": "time" - if conf.HISTOGRAM_SURFACE_FIELD_SEQUENCE - else "realization", - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, - "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", - }, - aggregation=False, - ) - - # Timeseries for surface (2D) Nimrod radar hourly accumulated rainfall. - # - # The timeseries are produced after the rainfall obs have been masked using - # the associated Nimrod weights file. - # - # To get multiple radar sources plotted on the histogram the - # recipe must be done by passing lists of the radar_ids and - # the radar_names. As this is a multiline plot, all radar sources - # share the same radar variable name. - if conf.TIMESERIES_SURFACE_FIELD: - radar_obs_ids = [radar["id"] for radar in accum_radars] - radar_wts_ids = [radar["id"] for radar in wts_radars] - combined_ids = radar_obs_ids + radar_wts_ids - print(" combined_ids: ", combined_ids) - yield RawRecipe( - recipe="radar_masked_mean_time_series.yaml", - # model_ids -> Becomes $INPUT_PATHS - # model_ids=[ radar_obs_ids, radar_wts_ids], - model_ids=combined_ids, - variables={ - "VARNAME": next(radar["varname"] for radar in accum_radars), - "ALL_NAME": combined_ids, - "RADAR_NAME": [radar["id"] for radar in accum_radars], - "WEIGHTS_NAME": [radar["id"] for radar in wts_radars], - "SEQUENCE": "realisation", - # "SEQUENCE": "time" - # if conf.HISTOGRAM_SURFACE_FIELD_SEQUENCE - # else "realization", - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, - }, - aggregation=False, - ) - - # Timeseries plot of Nimrod hourly surface rainfall accumulation. - if conf.TIMESERIES_SURFACE_FIELD and accum_radars: - yield RawRecipe( - recipe="radar_mean_time_series.yaml", - # model_ids -> Becomes $INPUT_PATHS - model_ids=[radar["id"] for radar in accum_radars], - variables={ - "VARNAME": next(radar["varname"] for radar in accum_radars), - "MODEL_NAME": [radar["name"] for radar in accum_radars], - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, - "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", - }, - aggregation=False, - ) - - # Timeseries plot of Nimrod 5 minute rainfall rate. - if conf.TIMESERIES_SURFACE_FIELD and conf.NIMROD_COMP_5MIN: - yield RawRecipe( - recipe="radar_mean_time_series.yaml", - # model_ids -> Becomes $INPUT_PATHS - model_ids="Nimrod5min", - variables={ - "VARNAME": "Rainfall rate Composite", - "MODEL_NAME": "Nimrod_5min", - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, - "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", - }, - aggregation=False, - ) From a2433b4d109935fcf75b99b0ec9dea2cb6460ade Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 13:52:18 +0100 Subject: [PATCH 056/105] Typos in rose-meta.conf for Rose GUI #2113 --- .../meta/observations/rose-meta.conf | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/src/CSET/cset_workflow/meta/observations/rose-meta.conf b/src/CSET/cset_workflow/meta/observations/rose-meta.conf index e21174c3d8..68316fdccf 100644 --- a/src/CSET/cset_workflow/meta/observations/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/observations/rose-meta.conf @@ -286,10 +286,21 @@ title=Produce Nimrod radar histograms description=Produce histograms incorporating radar observations. help=Switch to produce histograms utilising Nimrod radar observations. The histograms will combine both model rainfall and radar observations. - The model rainfall variable is selected in in the Rose tab Diagnostics/Surface (2D) fields. + The model rainfall variable is selected in the Rose tab Diagnostics/Surface (2D) fields. compulsory=true type=python_boolean -sort-key=nimrod-radar-9 +sort-key=nimrod-radar-9a + +#[template variables=PROCESS_RADAR_TIMESERIES] +#ns=Observations/Radar +#title=Produce Nimrod radar timeseries +#description=Produce timeseries incorporating radar observations. +#help=Switch to produce timeseries utilising Nimrod radar observations. +# The timeseries will combine both model rainfall and radar observations. +# The model rainfall variable is selected in the Rose tab Diagnostics/Surface (2D) fields. +#compulsory=true +#type=python_boolean +#sort-key=nimrod-radar-9b From 771adf7f1b40f51235b6901c3f0b932a34cf46b8 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 13:57:21 +0100 Subject: [PATCH 057/105] Tidied operator radar_filter.mask_list #2113 --- src/CSET/operators/radar_filter.py | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index a0866fad41..953b0ece97 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -29,24 +29,29 @@ def mask_list(model_names: list[str]) -> list[str]: Parameters ---------- model_names: list[str] - A list of model and Nimrod hourly rainfall accumulation files. + A list containing model names and at least one Nimrod hourly + rainfall accumulation source. + Possible radar sources are: + "Nimrod2km", "Nimrod_2km". + "Nimrodxkm", "Nimrod_xkm". + "Nimrod1km", "Nimrod_1km". Returns ------- list[str] A list of the Nimrod weights files to use with each of the input - model / observations files. + models and radar sources. Notes ----- - At lest one of the entries in the input list must be a Nimrod hourly - rainfall accumulation file. + At least one of the entries in the input list must be a Nimrod hourly + rainfall accumulation source. - If just one Nimrod file is specified, then then the weights file associated - with this field is used. + If just one Nimrod source is specified, then the weights file associated + with this source is used. - If more than one Nimrod file is in the input list, then each of the these - Nimrod files is associated with its own weights file e.g. if the input list + If more than one Nimrod source is in the input list, then each of the these + Nimrod sources is associated with its own weights file e.g. if the input list contains ["Nimrod1km", "Nimrod2km"] then the weights files for these will be ["Nimrod1km_weights", "Nimrod2km_weights"]. Any model fields in the input list will be allocated a weights file according to the order of preference From aa56edf041439008a9b84cc26bc76adac6e6958b Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 14:12:50 +0100 Subject: [PATCH 058/105] Tidied operator radar_filter.mask_by_weights #2113 --- src/CSET/operators/radar_filter.py | 53 ++++++------------------------ 1 file changed, 10 insertions(+), 43 deletions(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index 953b0ece97..a9e32bc4fe 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -108,42 +108,28 @@ def mask_by_weights( weights_names: list[str], **kwargs, ) -> iris.cube.CubeList: - """Filter a field using a second field as a mask. + """Filter a field using a radar weights field as a mask. Parameters ---------- cubes: iris.cube.CubeList - Two cubes containing the radar observations and their weights. + CubeList containing fields to mask and radar weights fields to use as the masks. + model_names: list[str] + A list of the model_names or radar sources to mask. + weights_names: list[str] + A list of radar weights sources to use as masks. There should be an entry + in weights_names to correspond with every entry in model_names. Returns ------- - Cube - - Raises - ------ - ValueError, iris.exceptions.NotYetImplementedError - When the cubes are not compatible. - - Notes - ----- - This is a simple operator designed for combination of diagnostics or - creating new diagnostics by using recipes. + CubeList: + A CubeList of masked fields. Examples -------- >>> field_filtered = mask_by_weights(cubelist, model_names) """ - print("model_names are: ", model_names) - print("weights_names", weights_names) - - for cube in cubes: - print(" cube.var_name ", cube.var_name) - print(" cube.name ", cube.name) - print(" cube: ") - print(cube) - print(" cube.attributes.model_name ", cube.attributes["model_name"]) - # Check the input unfiltered cubes and the mask cubes are both cubelists # with the same number of cubes. If not, then add extra mask cubes. if len(model_names) != len(weights_names): @@ -159,11 +145,8 @@ def mask_by_weights( iter_maybe(model_names), iter_maybe(weights_names), strict=True, - # iter_maybe(model_names), iter_maybe(weights_names), strict=True ): - print(" model, mask ", model, mask) - - # grab the field to filter + # Grab the field to filter. model_constraint = iris.AttributeConstraint(model_name=model) unfiltered_field = cubes.extract_cube(var_constraint & model_constraint) @@ -176,22 +159,6 @@ def mask_by_weights( # as can the threshold value. mask_radar_wts = generate_mask(mask_field, "ge", 11) - # print(" This is cube radar_obs: ", radar_obs) - # print(" This is cube radar_weights: ", radar_wts) - # print(" This is cube unfiltered: ", unfiltered_field) - - # check the coords of the unfiltered field and the mask field. - # If these do not match, then regrid the unfiltered field onto - # the grid used for the mask field. - # For radar weights fields can use the function regrid_onto_xyspacing in regrid.py, - # but then might have to extract a subarea to match the mask grid. - # Might have to consider serval cases for regridding: - # (1) model_field(lat, lon) to radar_weights_field(x, y) - # (2) model_field(lat, lon) to other_model_field(lat, lon) - # (3) Nimrod_field(x, y) to radar_weights_field(x, y) - # (4) Nimrod_field(x, y) to model_field(lat,lon) ? - # - # Apply the mask. masked_radar_obs = apply_mask(unfiltered_field, mask_radar_wts) From d45db95a2e25db271d5f152b745d1460b98c857a Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 14:19:22 +0100 Subject: [PATCH 059/105] Tidied operator radar_filter.radar_apply_mask #2113 --- src/CSET/operators/radar_filter.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index a9e32bc4fe..05c2800d0e 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -177,7 +177,7 @@ def radar_apply_mask( mask: iris.cube.Cube | iris.cube.CubeList, boundary_margin: int = 8, ) -> iris.cube.Cube | iris.cube.CubeList: - """Apply a mask to given data as a masked array. + """Apply a mask to given field as a masked array. Parameters ---------- @@ -192,19 +192,19 @@ def radar_apply_mask( Returns ------- masked_field: iris.cube.Cube | iris.cube.CubeList - A cube or cubelist of the masked field(s). + A cube or CubeList of the masked field(s). Notes ----- The mask is first converted to 1s and NaNs before multiplication with the original data. - As discussed in generate_mask, you can combine multiple masks in a + As discussed in filters.generate_mask, you can combine multiple masks in a recipe using other functions before applying the mask to the data. Examples -------- - >>> land_points_only = radar_apply_mask( surface_microphysical_rainfall_rate, Nimrod2km) + >>> radar_domain_only = radar_apply_mask( surface_microphysical_rainfall_rate, Nimrod2km_wts) """ # Create an empty cubelist to hold the filtered fields. masked_fields = iris.cube.CubeList([]) From 6c0a7adae21c782fbb1c068b0b37b978f671200b Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 14:30:45 +0100 Subject: [PATCH 060/105] Tidied operator radar_filter.radar_mask #2113 --- src/CSET/operators/radar_filter.py | 35 +++++------------------------- 1 file changed, 6 insertions(+), 29 deletions(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index 05c2800d0e..ec43aea463 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -258,67 +258,44 @@ def radar_mask( filtered_radar = iris.cube.CubeList([]) filtered_model = iris.cube.CubeList([]) - # ensure the three inputs to this function are model_field, nimrod_field, nimrod_mask - + # Loop over nimrod_field, model_field and nimrod_mask. for M, F, N in zip( iter_maybe(nimrod_mask), iter_maybe(model_field), iter_maybe(nimrod_field), strict=True, ): - # masked_field = F.copy() - - # apply the function radar_apply_mask to generate the re-gridded and masked model_field - # i.e. generate masked_model_field + # Apply the function radar_apply_mask to generate the re-gridded + # and masked model field. masked_model_field = radar_apply_mask(F, M, boundary_margin=boundary_margin) - # use the model_field as the mask for the nimrod_field - note: no re-gridding required - # i.e. generate masked_nimrod_field - + # Use the masked model field as the mask for the Nimrod field. + # Note: no re-gridding required. min_timesteps = min(N.shape[0], masked_model_field.shape[0]) temp_mask = masked_model_field[0:min_timesteps].copy() - # temp_mask.data[ temp_mask.data != np.nan ] = 1.0 temp_mask.data[~np.isnan(temp_mask.data)] = 1.0 masked_nimrod_field = N[0:min_timesteps].copy() - # masked_nimrod_field.data = N[0:min_timesteps].data * temp_mask[0:min_timesteps].data masked_nimrod_field.data *= temp_mask.data # Append the masked field to the output list of masked fields. filtered_model.append(masked_model_field) filtered_radar.append(masked_nimrod_field) - # return masked_model_field and masked_nimrod_field --> these can then be passed to either - # the histogram or time series plotting operators. - - print(" test point bmc1") + # Return either the masked model or Nimrod fields, or both. if outputs == "radar": - print(" test point bmc1a") filtered_fields.append(filtered_radar.merge_cube()) - print(" test point bmc2") if outputs == "model": - print(" test point bmc2a") filtered_fields.append(filtered_model.merge_cube()) - print(" test point bmc3") if outputs == "all": - print(" test point bmc4") filtered_fields.append(filtered_model.merge_cube()) - print(" test point bmc5") filtered_fields.append(filtered_radar.merge_cube()) - print("+++++++++++++++++++++++++++++") - print(filtered_fields) - print("+++++++++++++++++++++++++++++") # Return either a single cube or a cubelist. - print(" test point bmc6") if len(filtered_fields) == 1: - print(" test point bmc7a") return filtered_fields[0] else: - # return masked_fields - # return filtered_fields.merge() - print(" test point bmc7") return filtered_fields From 48192abcdecd5d3dd5c5bb6503e23c0127461f85 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 14:40:53 +0100 Subject: [PATCH 061/105] Tidied and added to docstring operator radar_filter.radar_mask #2113 --- src/CSET/operators/radar_filter.py | 31 +++++++++++++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index ec43aea463..2f7668ad00 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -252,7 +252,36 @@ def radar_mask( boundary_margin: int = 8, outputs: str = "radar", ) -> iris.cube.Cube | iris.cube.CubeList: - """Apply a mask to given data as a masked array.""" + """Apply a mask to given fields using a masked array. + + Parameters + ---------- + model_field: iris.cube.Cube | iris.cube.CubeList + The model field(s) to be masked. + nimrod_field: iris.cube.Cube | iris.cube.CubeList + The Nimrod field(s) to be masked. + nimrod_mask: iris.cube.Cube | iris.cube.CubeList + The Nimrod mask(s) to use. These are normally Nimrod wts fields. + boundary_margin: int, optional + Number of grid points from the domain boundary considered "unreliable". + Defaults to 8. + outputs: str, optional + Specifies which outputs are required: + "radar" outputs masked Nimrod rainfall field(s) + "model" outputs masked model field(s) + "all" outputs both masked model and masked Nimrod field(s). + + Returns + ------- + masked_field: iris.cube.Cube | iris.cube.CubeList + A cube or CubeList of the masked field(s). + + Examples + -------- + To mask both model and Nimrod fields using the Nimrod weights in nimrod_wts: + >>> masked_fields = radar_mask( model_fields, nimrod_fields, nimrod_wts, outputs="all") + + """ # Create an empty cubelist to hold the filtered fields. filtered_fields = iris.cube.CubeList([]) filtered_radar = iris.cube.CubeList([]) From 5bfa342ad46d92399bcc9699dc9ef80d069b59c0 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 14:48:40 +0100 Subject: [PATCH 062/105] Tidied and added to docstring operator radar_filter.radar_mask_loop #2113 --- src/CSET/operators/radar_filter.py | 41 ++++++++++++++++++++++-------- 1 file changed, 31 insertions(+), 10 deletions(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index 2f7668ad00..5ffd773bf3 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -335,15 +335,39 @@ def radar_mask_loop( boundary_margin: int = 8, outputs: str = "radar", ) -> iris.cube.Cube | iris.cube.CubeList: - """Find common domains between a list of models and radar observations.""" + """Find common domains between a list of models and radar sources. + + Parameters + ---------- + model_field: iris.cube.Cube | iris.cube.CubeList + The model field(s) to be masked. + nimrod_field: iris.cube.Cube | iris.cube.CubeList + The Nimrod field(s) to be masked. + nimrod_mask: iris.cube.Cube | iris.cube.CubeList + The Nimrod mask(s) to use. These are normally Nimrod wts fields. + boundary_margin: int, optional + Number of grid points from the domain boundary considered "unreliable". + Defaults to 8. + outputs: str, optional + Specifies which outputs are required: + "radar" outputs masked Nimrod rainfall field(s) + "model" outputs masked model field(s) + "all" outputs both masked model and masked Nimrod field(s). + + Returns + ------- + masked_field: iris.cube.Cube | iris.cube.CubeList + A cube or CubeList of the masked field(s). + + Examples + -------- + To mask both model and Nimrod fields using the Nimrod weights in nimrod_wts: + >>> masked_fields = radar_mask_loop( model_fields, nimrod_fields, nimrod_wts, outputs="all") + + """ # Create an empty cubelist to hold the filtered fields. filtered_cubes = iris.cube.CubeList([]) - # if len(nimrod_field) > 1: - # use_nimrod_field = nimrod_field[0] - # use_nimrod_mask = nimrod_mask[0] - # else: - print("-----> len(nimrod_field) ", len(nimrod_field)) use_nimrod_field = nimrod_field use_nimrod_mask = nimrod_mask @@ -359,8 +383,7 @@ def radar_mask_loop( ) filtered_cubes.append(filtered_model) - # Grab the filtered radar observations. - print("-----> len(model_field) ", len(model_field)) + # Filter the radar observations. if len(model_field) == 1: filtered_radar = radar_mask( model_field[0], @@ -377,10 +400,8 @@ def radar_mask_loop( boundary_margin=boundary_margin, outputs="radar", ) - filtered_cubes.append(filtered_radar) - print("-----------------> returning from radar_mask_loop") return filtered_cubes From d8b2c68af12e886ec4d4530a7f38646ea1c491a5 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 14:55:40 +0100 Subject: [PATCH 063/105] Tidied and added to docstring operator radar_filter.match_varname_and_units #2113 --- src/CSET/operators/radar_filter.py | 69 ++---------------------------- 1 file changed, 4 insertions(+), 65 deletions(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index 5ffd773bf3..56f0569f5b 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -429,44 +429,16 @@ def match_varname_and_units(cubes: iris.cube.Cube | iris.cube.CubeList): if len(cubes) == 1: return cubes - print("@@@@@@@@@@@@@@@@@@@@@ base cubes[0] @@@@@@@@@@@@") - print(cubes[0]) - print("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@") - - # TODO when testing completed remove the chop [2:25] [3:46] - # Initialise the list of matched cubes. new_cubelist = iris.cube.CubeList([]) - new_cubelist.append(cubes[0]) - # new_cubelist.append(cubes[0][2:45]) - # Loop the cubes to match to the first cube. + # Use the first cube in the CubeList as the base cube. base_cube = cubes[0] - print("@@@@@@@@@@@@@@@@@@@@@ base_cube @@@@@@@@@@@@") - print("base_cube.name : ", base_cube.name) - print("base_cube.standard_name : ", base_cube.standard_name) - print("base_cube.long_name : ", base_cube.long_name) - print("base_cube.var_name : ", base_cube.var_name) - print("base_cube.coords(time) : ", base_cube.coord("time")) - print("base_cube.coords(forecast_period) : ", base_cube.coord("forecast_period")) - # print("------- base cube below -----------") - # print(base_cube) - print("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@") + new_cubelist.append(base_cube) + + # Loop over the cubes matching each to the base cube. for cube in cubes[1:]: - print("@@@@@@@@@@@@@@@ raw @@@@@@@@@@@@@@@@@@") - print("cube.name : ", cube.name) - print("cube.standard_name : ", cube.standard_name) - print("cube.long_name : ", cube.long_name) - print("cube.var_name : ", cube.var_name) - print("cube.coords(time) : ", cube.coord("time")) - print("cube.coords(forecast_period) : ", cube.coord("forecast_period")) - print(cube) - print("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@") new_cube = cube.copy() - # if new_cube.attributes["model_name"] == "Nimrod2km": - # new_cube = new_cube[3:46] - # else: - # new_cube = new_cube[2:45] # Match the cube varname. new_cube.rename(base_cube.long_name) @@ -477,52 +449,19 @@ def match_varname_and_units(cubes: iris.cube.Cube | iris.cube.CubeList): new_cube.units = base_cube.units # Match the cube forecast_reference_time, a single value e.g. 2023-11-12 00:00:00. - # TODO move this to the read function for radar forecast_reference_time_use = base_cube.coord("forecast_reference_time") new_cube.remove_coord("forecast_reference_time") new_cube.add_aux_coord(forecast_reference_time_use) # new_cube.replace_coord(base_cube.coord('forecast_reference_time')) - # Match the cube forecast_period. - # forecast_period_use = base_cube.coord('forecast_period') - # new_cube.remove_coord("forecast_period") - # new_cube.replace_coord(forecast_period_use, data_dims=1) - - # coord_dims = base_cube.coord_dims('forecast_period') - # print("---> coord_dims", coord_dims) - # dim_index = coord_dims[0] if coord_dims else 0 - # print("---> dim_index", dim_index) - # new_coord = base_cube.coord('forecast_period') - # print("---> new_coord", new_coord) - # new_cube.remove_coord("forecast_period") - # new_cube.add_dim_coord(new_coord, dim_index) - # Match the cube forecast_period - an array points: [ 1., 2., ..., 47., 48.]. - # TODO move this to the read function for radar coord_dims = base_cube.coord_dims("forecast_period") - print("---> coord_dims", coord_dims) new_coord = base_cube.coord("forecast_period") - print("---> new_coord", new_coord) if new_cube.coords("forecast_period"): new_cube.remove_coord("forecast_period") new_cube.add_aux_coord(new_coord, coord_dims) - print("@@@@@@@@@@@@@@@ matched @@@@@@@@@@@@@@@@@@") - print("new_cube.name : ", new_cube.name) - print("new_cube.standard_name : ", new_cube.standard_name) - print("new_cube.long_name : ", new_cube.long_name) - print("new_cube.var_name : ", new_cube.var_name) - print("new_cube.coords(time) : ", new_cube.coord("time")) - print("new_cube.coords(forecast_period) : ", new_cube.coord("forecast_period")) - print(new_cube) - print("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@") - # Append the matched cube to the output cube list. new_cubelist.append(new_cube) - print( - "---------------->bmc exiting radar_filter.match_varname_and_units len(cubes) is ", - len(cubes), - ) - # print(bmc_halt_this) return new_cubelist From f40b3a64c27a50cf0acc001ffd14e9dd2a0c07ba Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 15:08:56 +0100 Subject: [PATCH 064/105] Model boundary margin is now templated in model_radar_common_domain_2d_model.yaml #2113 --- .../surface_fields/model_radar_common_domain_2d_model.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml b/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml index 3def2fe945..92ceb9eefd 100644 --- a/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml +++ b/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml @@ -53,7 +53,7 @@ steps: condition: 'ge' value: 11 - boundary_margin: 16 + boundary_margin: $BOUNDARY_MARGIN outputs: "model" - operator: collapse.collapse From 0775861c6d5e132c04363c12bcca4c6d65a72ec5 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 15:09:12 +0100 Subject: [PATCH 065/105] Model boundary margin is now templated in model_radar_common_domain_2d_radar.yaml #2113 --- .../surface_fields/model_radar_common_domain_2d_radar.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml b/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml index c5a29c3311..f73b4067d9 100644 --- a/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml +++ b/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml @@ -53,7 +53,7 @@ steps: condition: 'ge' value: 11 - boundary_margin: 16 + boundary_margin: $BOUNDARY_MARGIN outputs: "radar" - operator: collapse.collapse From ad49d0061f171812ed89d989cd2577814675d43a Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 15:09:59 +0100 Subject: [PATCH 066/105] Model boundary margin is now templated in radar_common_domain_histogram.yaml #2113 --- .../recipes/surface_fields/radar_common_domain_histogram.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml b/src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml index eb997e4b9b..00e2614961 100644 --- a/src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml +++ b/src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml @@ -56,7 +56,7 @@ steps: condition: 'ge' value: 11 - boundary_margin: 16 + boundary_margin: $BOUNDARY_MARGIN outputs: $OUTPUTS - operator: radar_filter.match_varname_and_units From ae32363f14693cafe79fbc3652b3b7a5238f90f3 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 15:10:27 +0100 Subject: [PATCH 067/105] Model boundary margin is now templated in radar loaders #2113 --- src/CSET/loaders/radar.py | 23 ++++++++--------------- 1 file changed, 8 insertions(+), 15 deletions(-) diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index 9b08952e09..cb0059d9b5 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -111,6 +111,9 @@ def load(conf: Config): if radar["varname"] == "Hourly rain accumulation" ] + # Specify the boundary margin width to use when trimming model fields. + boundary_margin = 16 + # Surface (2D) fields for Nimrod radar rainfall. # # These loaders produce 2D plots of both the Nimrod @@ -139,7 +142,6 @@ def load(conf: Config): # Radar 2D plots using a common domain between model and radar network. if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D and (len(accum_radars) > 0): - # field = "Hourly rain accumulation" radar_source = select_radar_source([source["id"] for source in accum_radars]) radar_obs_ids = [radar_source] radar_wts_ids = [radar_source + "_weights"] @@ -151,26 +153,20 @@ def load(conf: Config): "MODEL_VARNAME": "surface_microphysical_rainfall_rate", "RADAR_VARNAME": "Hourly rain accumulation", "MODEL_LABEL": model_names_list[0], - # "MODEL_LABEL": "ModelA", - # "RADAR_LABEL": "Nimrod2km", - # "MASK_LABEL": "Nimrod2km_wts", "RADAR_LABEL": radar_obs_ids[0], "MASK_LABEL": radar_wts_ids[0], "METHOD": "SEQ", - # "OUTPUTS": "radar", + "BOUNDARY_MARGIN": boundary_margin, "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", }, model_ids=[model_ids_list[0], "Nimrod2km", "Nimrod2km_weights"], - # model_ids=["1", "Nimrod2km", "Nimrod2km_weights"], - # model_ids=["ModelA", "Nimrod2km", "Nimrod2km_weights"], aggregation=False, ) - # Radar 2D plots using a common domain between model and radar network. + # Model rainfall 2D plots using a common domain between model and radar network. if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D and (len(accum_radars) > 0): - # field = "Hourly rain accumulation" radar_source = select_radar_source([source["id"] for source in accum_radars]) radar_obs_ids = [radar_source] radar_wts_ids = [radar_source + "_weights"] @@ -182,20 +178,15 @@ def load(conf: Config): "MODEL_VARNAME": "surface_microphysical_rainfall_rate", "RADAR_VARNAME": "Hourly rain accumulation", "MODEL_LABEL": model_names_list[0], - # "MODEL_LABEL": "ModelA", - # "RADAR_LABEL": "Nimrod2km", - # "MASK_LABEL": "Nimrod2km_wts", "RADAR_LABEL": radar_obs_ids[0], "MASK_LABEL": radar_wts_ids[0], "METHOD": "SEQ", - # "OUTPUTS": "model", + "BOUNDARY_MARGIN": boundary_margin, "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", }, model_ids=[model_ids_list[0], "Nimrod2km", "Nimrod2km_weights"], - # model_ids=["1", "Nimrod2km", "Nimrod2km_weights"], - # model_ids=["ModelA", "Nimrod2km", "Nimrod2km_weights"], aggregation=False, ) @@ -223,6 +214,7 @@ def load(conf: Config): "ALL_LABEL": combined_names, "SEQUENCE": "time", "OUTPUTS": "all", + "BOUNDARY_MARGIN": boundary_margin, "TITLE_STRING": "sequence plots", "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, @@ -256,6 +248,7 @@ def load(conf: Config): "ALL_LABEL": combined_names, "SEQUENCE": "realization", "OUTPUTS": "all", + "BOUNDARY_MARGIN": boundary_margin, "TITLE_STRING": "case study", "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, From 6d06605c1b6acab4b6b4e4cfafba789a5b910106 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 17 Sep 2026 16:20:47 +0100 Subject: [PATCH 068/105] Radar loader for model common domain loops over models #2113 --- src/CSET/loaders/radar.py | 44 ++++++++++++++++++++++----------------- 1 file changed, 25 insertions(+), 19 deletions(-) diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index cb0059d9b5..17e26b853a 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -122,7 +122,7 @@ def load(conf: Config): # The different sources of Nimrod rainfall accumulation have # different spatial grids. So each source requires its own # recipe to prevent incompatible cubes being created. - if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D and (len(accum_radars) > 0): + if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D and (len(radar_sources) > 0): for radar in radar_sources: yield RawRecipe( recipe="generic_surface_spatial_plot_sequence_radar_rainfall.yaml", @@ -161,7 +161,8 @@ def load(conf: Config): "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", }, - model_ids=[model_ids_list[0], "Nimrod2km", "Nimrod2km_weights"], + # model_ids=[model_ids_list[0], "Nimrod2km", "Nimrod2km_weights"], + model_ids=[model_ids_list[0], radar_obs_ids[0], radar_wts_ids[0]], aggregation=False, ) @@ -172,23 +173,28 @@ def load(conf: Config): radar_wts_ids = [radar_source + "_weights"] model_names_list = [model["name"] for model in models] model_ids_list = [model["id"] for model in models] - yield RawRecipe( - recipe="model_radar_common_domain_2d_model.yaml", - variables={ - "MODEL_VARNAME": "surface_microphysical_rainfall_rate", - "RADAR_VARNAME": "Hourly rain accumulation", - "MODEL_LABEL": model_names_list[0], - "RADAR_LABEL": radar_obs_ids[0], - "MASK_LABEL": radar_wts_ids[0], - "METHOD": "SEQ", - "BOUNDARY_MARGIN": boundary_margin, - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, - "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", - }, - model_ids=[model_ids_list[0], "Nimrod2km", "Nimrod2km_weights"], - aggregation=False, - ) + # Loop over the model names. + for idx, model_use in enumerate(model_names_list): + yield RawRecipe( + recipe="model_radar_common_domain_2d_model.yaml", + variables={ + "MODEL_VARNAME": "surface_microphysical_rainfall_rate", + "RADAR_VARNAME": "Hourly rain accumulation", + "MODEL_LABEL": model_use, + "RADAR_LABEL": radar_obs_ids[0], + "MASK_LABEL": radar_wts_ids[0], + "METHOD": "SEQ", + "BOUNDARY_MARGIN": boundary_margin, + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT + if conf.SELECT_SUBAREA + else None, + "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", + }, + # model_ids=[model_ids_list[0], "Nimrod2km", "Nimrod2km_weights"], + model_ids=[model_ids_list[idx], radar_obs_ids[0], radar_wts_ids[0]], + aggregation=False, + ) # Histogram sequence rainfall using common domain between # model and Nimrod radar observations. From 5ebe82741ccac609721a5f2cd45a7cbdb5cfc307 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Mon, 21 Sep 2026 09:44:53 +0100 Subject: [PATCH 069/105] Tidied radar loaders #2113 --- src/CSET/loaders/radar.py | 81 +++++++++++++++++++++++---------------- 1 file changed, 47 insertions(+), 34 deletions(-) diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index 17e26b853a..7a74a1cc08 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -96,6 +96,29 @@ def get_radar_sources(conf) -> list[dict]: return radar_sources +def _model_and_radar_names_ids(models, radar_source, accum_radars): + """Select the model and radar names and ids.""" + # Model and radar names and ids. + model_names_list = [model["name"] for model in models] + model_ids_list = [model["id"] for model in models] + radar_source = select_radar_source([source["id"] for source in accum_radars]) + radar_obs_ids = [radar_source] + radar_wts_ids = [radar_source + "_weights"] + + # Construct the dictionary to return. + model_and_radar = { + "radar_source": radar_source, + "radar_obs_ids": radar_obs_ids, + "radar_wts_ids": radar_wts_ids, + "model_names_list": model_names_list, + "model_ids_list": model_ids_list, + "combined_names": model_names_list + radar_obs_ids + radar_wts_ids, + "combined_ids": model_ids_list + radar_obs_ids + radar_wts_ids, + } + + return model_and_radar + + def load(conf: Config): """Yield recipes from the given workflow configuration.""" # Load a list of model detail dictionaries. @@ -143,46 +166,43 @@ def load(conf: Config): # Radar 2D plots using a common domain between model and radar network. if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D and (len(accum_radars) > 0): radar_source = select_radar_source([source["id"] for source in accum_radars]) - radar_obs_ids = [radar_source] - radar_wts_ids = [radar_source + "_weights"] - model_names_list = [model["name"] for model in models] - model_ids_list = [model["id"] for model in models] + names_ids = _model_and_radar_names_ids(models, radar_source, accum_radars) yield RawRecipe( recipe="model_radar_common_domain_2d_radar.yaml", variables={ "MODEL_VARNAME": "surface_microphysical_rainfall_rate", "RADAR_VARNAME": "Hourly rain accumulation", - "MODEL_LABEL": model_names_list[0], - "RADAR_LABEL": radar_obs_ids[0], - "MASK_LABEL": radar_wts_ids[0], + "MODEL_LABEL": names_ids["model_names_list"][0], + "RADAR_LABEL": names_ids["radar_obs_ids"][0], + "MASK_LABEL": names_ids["radar_wts_ids"][0], "METHOD": "SEQ", "BOUNDARY_MARGIN": boundary_margin, "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", }, - # model_ids=[model_ids_list[0], "Nimrod2km", "Nimrod2km_weights"], - model_ids=[model_ids_list[0], radar_obs_ids[0], radar_wts_ids[0]], + model_ids=[ + names_ids["model_ids_list"][0], + names_ids["radar_obs_ids"][0], + names_ids["radar_wts_ids"][0], + ], aggregation=False, ) # Model rainfall 2D plots using a common domain between model and radar network. if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D and (len(accum_radars) > 0): radar_source = select_radar_source([source["id"] for source in accum_radars]) - radar_obs_ids = [radar_source] - radar_wts_ids = [radar_source + "_weights"] - model_names_list = [model["name"] for model in models] - model_ids_list = [model["id"] for model in models] + names_ids = _model_and_radar_names_ids(models, radar_source, accum_radars) # Loop over the model names. - for idx, model_use in enumerate(model_names_list): + for idx, model_use in enumerate(names_ids["model_names_list"]): yield RawRecipe( recipe="model_radar_common_domain_2d_model.yaml", variables={ "MODEL_VARNAME": "surface_microphysical_rainfall_rate", "RADAR_VARNAME": "Hourly rain accumulation", "MODEL_LABEL": model_use, - "RADAR_LABEL": radar_obs_ids[0], - "MASK_LABEL": radar_wts_ids[0], + "RADAR_LABEL": names_ids["radar_obs_ids"][0], + "MASK_LABEL": names_ids["radar_wts_ids"][0], "METHOD": "SEQ", "BOUNDARY_MARGIN": boundary_margin, "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, @@ -191,8 +211,11 @@ def load(conf: Config): else None, "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", }, - # model_ids=[model_ids_list[0], "Nimrod2km", "Nimrod2km_weights"], - model_ids=[model_ids_list[idx], radar_obs_ids[0], radar_wts_ids[0]], + model_ids=[ + names_ids["model_ids_list"][idx], + names_ids["radar_obs_ids"][0], + names_ids["radar_wts_ids"][0], + ], aggregation=False, ) @@ -205,19 +228,14 @@ def load(conf: Config): ): # Select the radar source to use. radar_source = select_radar_source([source["id"] for source in accum_radars]) - radar_obs_ids = [radar_source] - radar_wts_ids = [radar_source + "_weights"] - model_names_list = [model["name"] for model in models] - model_ids_list = [model["id"] for model in models] - combined_names = model_names_list + radar_obs_ids + radar_wts_ids - combined_ids = model_ids_list + radar_obs_ids + radar_wts_ids + names_ids = _model_and_radar_names_ids(models, radar_source, accum_radars) yield RawRecipe( recipe="radar_common_domain_histogram.yaml", variables={ "MODEL_VARNAME": "surface_microphysical_rainfall_rate", "RADAR_VARNAME": "Hourly rain accumulation", "RADAR_WTS_VARNAME": "Hourly wts accumulation", - "ALL_LABEL": combined_names, + "ALL_LABEL": names_ids["combined_names"], "SEQUENCE": "time", "OUTPUTS": "all", "BOUNDARY_MARGIN": boundary_margin, @@ -226,7 +244,7 @@ def load(conf: Config): "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", }, - model_ids=combined_ids, + model_ids=names_ids["combined_ids"], aggregation=False, ) @@ -239,19 +257,14 @@ def load(conf: Config): ): # Select the radar source to use. radar_source = select_radar_source([source["id"] for source in accum_radars]) - radar_obs_ids = [radar_source] - radar_wts_ids = [radar_source + "_weights"] - model_names_list = [model["name"] for model in models] - model_ids_list = [model["id"] for model in models] - combined_names = model_names_list + radar_obs_ids + radar_wts_ids - combined_ids = model_ids_list + radar_obs_ids + radar_wts_ids + names_ids = _model_and_radar_names_ids(models, radar_source, accum_radars) yield RawRecipe( recipe="radar_common_domain_histogram.yaml", variables={ "MODEL_VARNAME": "surface_microphysical_rainfall_rate", "RADAR_VARNAME": "Hourly rain accumulation", "RADAR_WTS_VARNAME": "Hourly wts accumulation", - "ALL_LABEL": combined_names, + "ALL_LABEL": names_ids["combined_names"], "SEQUENCE": "realization", "OUTPUTS": "all", "BOUNDARY_MARGIN": boundary_margin, @@ -260,6 +273,6 @@ def load(conf: Config): "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", }, - model_ids=combined_ids, + model_ids=names_ids["combined_ids"], aggregation=False, ) From 0f66a8967dc72dad70f1f8424e5faa53f01cdd93 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Mon, 21 Sep 2026 11:10:27 +0100 Subject: [PATCH 070/105] Test for radar_filter.mask_list #2113 --- tests/operators/test_radar_filter.py | 29 ++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) create mode 100644 tests/operators/test_radar_filter.py diff --git a/tests/operators/test_radar_filter.py b/tests/operators/test_radar_filter.py new file mode 100644 index 0000000000..29b61d35b0 --- /dev/null +++ b/tests/operators/test_radar_filter.py @@ -0,0 +1,29 @@ +# © Crown copyright, Met Office (2022-2026) 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 filter operators.""" + +from CSET.operators import constraints, radar_filter + +constraint_single = constraints.combine_constraints( + constraints.generate_stash_constraint("m01s03i236"), + a=constraints.generate_cell_methods_constraint([]), +) + + +def test_mask_list(): + """Test the mask_list function.""" + list_weights = radar_filter.mask_list(["UM_model", "Nimrod1km", "Nimrod2km"]) + expected_list = ["Nimrod2km_weights", "Nimrod1km_weights", "Nimrod2km_weights"] + assert list_weights == expected_list From caf8c62b53831268b634eab080aed69dd1f28ec1 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Mon, 21 Sep 2026 11:22:28 +0100 Subject: [PATCH 071/105] Changed L537 in _colormaps.py from logging.info to logger.info #2113 --- src/CSET/operators/_colormaps.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/operators/_colormaps.py b/src/CSET/operators/_colormaps.py index 3061a3ca74..9e19025615 100644 --- a/src/CSET/operators/_colormaps.py +++ b/src/CSET/operators/_colormaps.py @@ -534,7 +534,7 @@ def custom_colormap_precipitation(cube: iris.cube.Cube, cmap, levels, norm): cmap = mcolors.ListedColormap(colors) # Normalize the levels norm = mcolors.BoundaryNorm(levels, cmap.N) - logging.info("Using custom rainfall colourmap.") + logger.info("Using custom rainfall colourmap.") # Set any Nan values to be plotted a light grey. cmap.set_bad("#dcdcdc") From 103ada3b3f20c53862ef68ef48e7af6559b7a0ce Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Mon, 21 Sep 2026 11:31:55 +0100 Subject: [PATCH 072/105] Updated rose-suite.conf.example with PROCESS_RADAR_2D and PROCESSS_RADAR_HISTOGRAMS #2113 --- src/CSET/cset_workflow/rose-suite.conf.example | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/CSET/cset_workflow/rose-suite.conf.example b/src/CSET/cset_workflow/rose-suite.conf.example index af76fd0172..a9095bc8e5 100644 --- a/src/CSET/cset_workflow/rose-suite.conf.example +++ b/src/CSET/cset_workflow/rose-suite.conf.example @@ -104,6 +104,8 @@ MODERATE_RAIN_PRESENCE_SPATIAL_PLOT=False !!NIMROD_COMP_2KM=False !!NIMROD_COMP_5MIN=False !!NIMROD_COMP_XKM=False +!!PROCESS_RADAR_2D=False +!!PROCESS_RADAR_HISTOGRAMS=False NIMROD_RADAR_OBS=False !!NIMROD_WEIGHTS=False !!PLEVEL_TRANSECT_AGGREGATION=False,False,False,False From 8af3a65a5d01f303ab603532ad7f0a66ee2926d7 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Mon, 21 Sep 2026 13:23:01 +0100 Subject: [PATCH 073/105] Cleaned collapse.py - removed test prints #2113 --- src/CSET/operators/collapse.py | 27 --------------------------- 1 file changed, 27 deletions(-) diff --git a/src/CSET/operators/collapse.py b/src/CSET/operators/collapse.py index 99ee68c7cb..44e4b48af5 100644 --- a/src/CSET/operators/collapse.py +++ b/src/CSET/operators/collapse.py @@ -76,30 +76,6 @@ def collapse( if method == "PERCENTILE" and additional_percent is None: raise ValueError("Must specify additional_percent") - print("bmcline77 starting collapse.collapse and value of method is ", method) - - cube_number = 0 - for new_cube in cubes: - print( - "@@@@ collapse.collapse @@@@@@@ cube number ", - cube_number, - "@@@@@@@@@@@@@@@@@@", - ) - print("new_cube.name : ", new_cube.name) - print("new_cube.standard_name : ", new_cube.standard_name) - print("new_cube.long_name : ", new_cube.long_name) - print("new_cube.var_name : ", new_cube.var_name) - print("new_cube.coords(time) : ", new_cube.coord("time")) - print( - "new_cube.coords(forecast_reference_time : ", - new_cube.coord("forecast_reference_time"), - ) - print("new_cube.coords(forecast_period) : ", new_cube.coord("forecast_period")) - print(new_cube) - print("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@") - cube_number = cube_number + 1 - # print(bmc_halt_this) - # Retain only common time points between different models if multiple model inputs. if isinstance(cubes, iris.cube.CubeList) and len(cubes) > 1: logger.debug("Extracting common time points as multiple model inputs detected.") @@ -112,8 +88,6 @@ def collapse( if len(cubes) == 0: raise ValueError("No overlapping times detected in input cubes.") - # print("bmcline92 collapse.collapse after common time points section") - collapsed_cubes = iris.cube.CubeList([]) with warnings.catch_warnings(): warnings.filterwarnings( @@ -142,7 +116,6 @@ def collapse( collapsed_cubes.append( cube.collapsed(coordinate, getattr(iris.analysis, method)) ) - print("bmcline123 exiting collapse.collapse and value of method is ", method) if len(collapsed_cubes) == 1: return collapsed_cubes[0] else: From 4e42c4dd8285de1752cba8d8b081e7f8e3afcde2 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Mon, 21 Sep 2026 13:48:00 +0100 Subject: [PATCH 074/105] Removed commented out lines from meta/observations/rose-meta.conf #2113 --- .../cset_workflow/meta/observations/rose-meta.conf | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/src/CSET/cset_workflow/meta/observations/rose-meta.conf b/src/CSET/cset_workflow/meta/observations/rose-meta.conf index 68316fdccf..9601a12a89 100644 --- a/src/CSET/cset_workflow/meta/observations/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/observations/rose-meta.conf @@ -291,20 +291,6 @@ compulsory=true type=python_boolean sort-key=nimrod-radar-9a -#[template variables=PROCESS_RADAR_TIMESERIES] -#ns=Observations/Radar -#title=Produce Nimrod radar timeseries -#description=Produce timeseries incorporating radar observations. -#help=Switch to produce timeseries utilising Nimrod radar observations. -# The timeseries will combine both model rainfall and radar observations. -# The model rainfall variable is selected in the Rose tab Diagnostics/Surface (2D) fields. -#compulsory=true -#type=python_boolean -#sort-key=nimrod-radar-9b - - - - ################################################################ # Cardington UK Research Observations [ns=Observations/Cardington] From 7af61398f06c6acd4b31b8043708574207df6458 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Mon, 21 Sep 2026 13:53:42 +0100 Subject: [PATCH 075/105] Added the radar processing to trigger in meta/observations/rose-meta.conf #2113 --- src/CSET/cset_workflow/meta/observations/rose-meta.conf | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/CSET/cset_workflow/meta/observations/rose-meta.conf b/src/CSET/cset_workflow/meta/observations/rose-meta.conf index 9601a12a89..674ee63c96 100644 --- a/src/CSET/cset_workflow/meta/observations/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/observations/rose-meta.conf @@ -111,6 +111,8 @@ trigger=template variables=NIMROD_COMP_XKM: True; template variables=NIMROD_COMP_2KM: True; template variables=NIMROD_WEIGHTS: True; template variables=NIMROD_COMP_5MIN: True; + template variables=PROCESS_RADAR_2D: True; + template variables=PROCESS_RADAR_HISTOGRAMS: True; [template variables=CARDINGTON_OBS] ns=Observations From a61ebfcd6f8dbc2eb1724063a2a49a670812755a Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Mon, 21 Sep 2026 14:32:58 +0100 Subject: [PATCH 076/105] Updated test_fetch_nimrod to account for the +1hr for radar accumulations #2113 --- tests/test_fetch_nimrod.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_fetch_nimrod.py b/tests/test_fetch_nimrod.py index 7c2f2cb52b..5a4ec2fc34 100644 --- a/tests/test_fetch_nimrod.py +++ b/tests/test_fetch_nimrod.py @@ -105,8 +105,8 @@ def test_retrieve_nimrod_saves_observation_and_weights(tmp_path, monkeypatch): monkeypatch.setenv("NIMROD_COMP_2KM", "False") monkeypatch.setenv("NIMROD_COMP_5MIN", "False") monkeypatch.setenv("NIMROD_WEIGHTS", "True") - monkeypatch.setenv("CYLC_TASK_CYCLE_POINT", "2026-05-07T12:00:00") - monkeypatch.setenv("ANALYSIS_LENGTH", "PT00H") + monkeypatch.setenv("CYLC_TASK_CYCLE_POINT", "2026-04-07T12:00:00") + monkeypatch.setenv("ANALYSIS_LENGTH", "PT01H") rose_datac = tmp_path / "rose_datac" monkeypatch.setenv("ROSE_DATAC", str(rose_datac)) @@ -144,13 +144,13 @@ def save_side_effect(cube, filename): fetch_nimrod.apply_radar_weights.assert_called_once_with(obs_cube, weights_cube) expected_obs_file = str( - Path(rose_datac) / "data" / "nimrod_xkm" / "202605071200_nimrod_xkm.nc" + Path(rose_datac) / "data" / "nimrod_xkm" / "202604071300_nimrod_xkm.nc" ) expected_weight_file = str( Path(rose_datac) / "data" / "nimrod_xkm_weights" - / "202605071200_nimrod_xkm_weights.nc" + / "202604071300_nimrod_xkm_weights.nc" ) saved_paths = [path for _, path in save_calls] From 4875f8e2a794b138dc4576f2d6a5caae2c4cb9af Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Mon, 21 Sep 2026 15:41:18 +0100 Subject: [PATCH 077/105] Added test for radar_filter.mask_by_weights #2113 --- tests/operators/test_radar_filter.py | 35 ++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/tests/operators/test_radar_filter.py b/tests/operators/test_radar_filter.py index 29b61d35b0..453347f33a 100644 --- a/tests/operators/test_radar_filter.py +++ b/tests/operators/test_radar_filter.py @@ -14,6 +14,10 @@ """Test filter operators.""" +import iris +import numpy as np +from iris.cube import Cube + from CSET.operators import constraints, radar_filter constraint_single = constraints.combine_constraints( @@ -27,3 +31,34 @@ def test_mask_list(): list_weights = radar_filter.mask_list(["UM_model", "Nimrod1km", "Nimrod2km"]) expected_list = ["Nimrod2km_weights", "Nimrod1km_weights", "Nimrod2km_weights"] assert list_weights == expected_list + + +def test_mask_by_weights(): + """Test the mask_by_weights function.""" + obs = Cube([[1.0, 2.0], [3.0, 4.0]]) + wts = Cube([[13, 0], [10, 11]]) + out = Cube([[1.0, np.nan], [np.nan, 4.0]]) + print("----> 1 obs ", obs) + + var_name_obs = "hourly_rain_accumulation" + var_name_wts = "hourly_wts_accumulation" + + obs_name = "Radar" + wts_name = "Radar_wts" + obs.attributes["model_name"] = obs_name + wts.attributes["model_name"] = wts_name + + obs.rename(var_name_obs) + obs.long_name = var_name_obs + obs.var_name = var_name_obs + + wts.rename(var_name_wts) + wts.long_name = var_name_wts + wts.var_name = var_name_wts + + cube_list = iris.cube.CubeList([]) + cube_list.append(obs) + cube_list.append(wts) + masked = radar_filter.mask_by_weights(cube_list, [obs_name], [wts_name]) + + assert np.array_equal(masked.data, out.data, equal_nan=True) From d0c202964ceab40e13de2075a7143290e967b147 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Mon, 21 Sep 2026 15:46:47 +0100 Subject: [PATCH 078/105] Added test for radar_filter.mask_by_weights #2113 --- tests/operators/test_radar_filter.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/operators/test_radar_filter.py b/tests/operators/test_radar_filter.py index 453347f33a..b6cf749526 100644 --- a/tests/operators/test_radar_filter.py +++ b/tests/operators/test_radar_filter.py @@ -60,5 +60,4 @@ def test_mask_by_weights(): cube_list.append(obs) cube_list.append(wts) masked = radar_filter.mask_by_weights(cube_list, [obs_name], [wts_name]) - assert np.array_equal(masked.data, out.data, equal_nan=True) From f691244a4374a2e6af259550cc9efac3aa442ce3 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Tue, 22 Sep 2026 08:07:45 +0100 Subject: [PATCH 079/105] Minor edit to a radar test #2113 --- tests/operators/test_radar_filter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/operators/test_radar_filter.py b/tests/operators/test_radar_filter.py index b6cf749526..895ef28a4e 100644 --- a/tests/operators/test_radar_filter.py +++ b/tests/operators/test_radar_filter.py @@ -38,7 +38,6 @@ def test_mask_by_weights(): obs = Cube([[1.0, 2.0], [3.0, 4.0]]) wts = Cube([[13, 0], [10, 11]]) out = Cube([[1.0, np.nan], [np.nan, 4.0]]) - print("----> 1 obs ", obs) var_name_obs = "hourly_rain_accumulation" var_name_wts = "hourly_wts_accumulation" @@ -59,5 +58,6 @@ def test_mask_by_weights(): cube_list = iris.cube.CubeList([]) cube_list.append(obs) cube_list.append(wts) + masked = radar_filter.mask_by_weights(cube_list, [obs_name], [wts_name]) assert np.array_equal(masked.data, out.data, equal_nan=True) From 83ed54dc86846e2ed2331770911507c2597a8044 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Tue, 22 Sep 2026 11:20:56 +0100 Subject: [PATCH 080/105] Minor typo corrected #2113 --- src/CSET/operators/radar_filter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index 56f0569f5b..e2efbeaaab 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -231,7 +231,7 @@ def radar_apply_mask( min_timesteps = min(M.shape[0], masked_field.shape[0]) masked_field = apply_mask(masked_field[0:min_timesteps], M[0:min_timesteps]) - # Attach and attribute to the masked field detailing the mask used. + # Attach an attribute to the masked field detailing the mask used. masked_field.attributes["mask"] = f"mask_of_{F.name()}" # Append the masked field to the output list of masked fields. From d52c766cbf6a0ef4f3581c30a2606999bf2f12c8 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Wed, 23 Sep 2026 16:11:40 +0100 Subject: [PATCH 081/105] Test for radar_filter.match_varname_and_units #2113 --- tests/operators/test_radar_filter.py | 154 +++++++++++++++++++++++++++ 1 file changed, 154 insertions(+) diff --git a/tests/operators/test_radar_filter.py b/tests/operators/test_radar_filter.py index 895ef28a4e..2d8b4fb964 100644 --- a/tests/operators/test_radar_filter.py +++ b/tests/operators/test_radar_filter.py @@ -61,3 +61,157 @@ def test_mask_by_weights(): masked = radar_filter.mask_by_weights(cube_list, [obs_name], [wts_name]) assert np.array_equal(masked.data, out.data, equal_nan=True) + + +def test_match_varname_and_units(): + """Test the match_varname_and_units function.""" + var_name_cube1 = "surface_microphysical_rainfall_rate" + var_name_cube3 = "Hourly rain accumulation" + + units_cube1 = "mm hr-1" + units_cube3 = "mm" + + # Construct test_cube1. + shape = (5, 4, 4) + data = np.ones(shape, dtype=np.float32) + + lat_coord = iris.coords.DimCoord( + np.linspace(-90, 90, 4), standard_name="latitude", units="degrees" + ) + lon_coord = iris.coords.DimCoord( + np.linspace(-180, 180, 4), standard_name="longitude", units="degrees" + ) + + time_points1 = np.array([0, 6, 12, 18, 24], dtype=np.float64) + time_coord1 = iris.coords.DimCoord( + time_points1, + standard_name="time", + units="hours since 2026-09-23 00:00:00", + ) + + fp_coord1 = iris.coords.DimCoord( + time_points1, standard_name="forecast_period", units="hours" + ) + + test_cube1 = iris.cube.Cube( + data, + long_name=var_name_cube1, + units=units_cube1, + dim_coords_and_dims=[ + (time_coord1, 0), + (lat_coord, 1), + (lon_coord, 2), + ], + ) + + test_cube1.add_aux_coord(fp_coord1, data_dims=0) + + frt_coord1 = iris.coords.AuxCoord( + [0.0], + standard_name="forecast_reference_time", + units="hours since 2026-09-23 00:00:00", + ) + test_cube1.add_aux_coord(frt_coord1) + + # Construct test_cube2. + test_cube2 = test_cube1.copy() + + # Construct test_cube3. + time_points3 = np.array([3, 9, 15, 21, 27], dtype=np.float64) + + time_coord3 = iris.coords.DimCoord( + time_points3, + standard_name="time", + units="hours since 2026-09-23 00:00:00", + ) + + fp_coord3 = iris.coords.DimCoord( + time_points3, standard_name="forecast_period", units="hours" + ) + + test_cube3 = iris.cube.Cube( + data, + long_name=var_name_cube3, + units=units_cube3, + dim_coords_and_dims=[ + (time_coord3, 0), + (lat_coord, 1), + (lon_coord, 2), + ], + ) + + test_cube3.add_aux_coord(fp_coord3, data_dims=0) + + frt_coord3 = iris.coords.AuxCoord( + [1.0], + standard_name="forecast_reference_time", + units="hours since 2026-09-23 00:00:00", + ) + test_cube3.add_aux_coord(frt_coord3) + + print(" ") + print("-----> test_cube1 ", test_cube1) + print(" ") + print("-----> test_cube3 ", test_cube3) + + ######## + + # Construct the CubeList to test. + cube_list = iris.cube.CubeList([]) + cube_list.append(test_cube1) + cube_list.append(test_cube2) + cube_list.append(test_cube3) + + # Push the test CubeList through the function match_varname_and_units. + cubes_matched = radar_filter.match_varname_and_units(cube_list) + + assert cubes_matched[0].var_name == cubes_matched[1].var_name + assert cubes_matched[0].var_name == cubes_matched[2].var_name + + assert cubes_matched[0].long_name == cubes_matched[1].long_name + assert cubes_matched[0].long_name == cubes_matched[2].long_name + + assert cubes_matched[0].units == cubes_matched[1].units + assert cubes_matched[0].units == cubes_matched[2].units + + assert cubes_matched[0].coord("forecast_reference_time") == cubes_matched[1].coord( + "forecast_reference_time" + ) + assert cubes_matched[0].coord("forecast_reference_time") == cubes_matched[2].coord( + "forecast_reference_time" + ) + + assert cubes_matched[0].coord("forecast_period") == cubes_matched[1].coord( + "forecast_period" + ) + assert cubes_matched[0].coord("forecast_period") == cubes_matched[2].coord( + "forecast_period" + ) + + print(" ") + print("-----> cubes_matched[2] ", cubes_matched[2]) + + +# def test_radar_apply_mask(): +# """Test the radar_apply_mask function.""" +# obs = Cube([0,[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]]) +# #wts = Cube([[13, 0, 13], [9, 12, 11], [13, 11, 10]]) +# msk = Cube([0,[[1, np.nan, 1], [np.nan, 1, 1], [1, 1, 1]]]) +# +# var_name_obs = "hourly_rain_accumulation" +# var_name_msk = "hourly_wts_accumulation" +# +# obs_name = "Radar" +# msk_name = "Radar_msk" +# obs.attributes["model_name"] = obs_name +# msk.attributes["model_name"] = msk_name +# +# obs.rename(var_name_obs) +# obs.long_name = var_name_obs +# obs.var_name = var_name_obs +# +# msk.rename(var_name_msk) +# msk.long_name = var_name_msk +# msk.var_name = var_name_msk +# +# masked = radar_filter.radar_apply_mask(obs, msk, boundary_margin = 1) From 6315d478cd154afab26bfa7f759d04545c910d45 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 24 Sep 2026 16:00:23 +0100 Subject: [PATCH 082/105] Tidied radar_filter.py #2113 --- src/CSET/operators/radar_filter.py | 29 +++++++++-------------------- 1 file changed, 9 insertions(+), 20 deletions(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index e2efbeaaab..1473f0ef41 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -139,8 +139,6 @@ def mask_by_weights( filtered_list = iris.cube.CubeList([]) # Loop over the fields to filter. - var_constraint = iris.NameConstraint(var_name="hourly_rain_accumulation") - mask_var_constraint = iris.NameConstraint(var_name="hourly_wts_accumulation") for model, mask in zip( iter_maybe(model_names), iter_maybe(weights_names), @@ -148,12 +146,12 @@ def mask_by_weights( ): # Grab the field to filter. model_constraint = iris.AttributeConstraint(model_name=model) - unfiltered_field = cubes.extract_cube(var_constraint & model_constraint) + unfiltered_field = cubes.extract_cube(model_constraint) # Select the field to use as the mask. # Nice to do - put in support for a static mask. mask_constraint = iris.AttributeConstraint(model_name=mask) - mask_field = cubes.extract_cube(mask_var_constraint & mask_constraint) + mask_field = cubes.extract_cube(mask_constraint) # Create the mask - note that the condition e.g. "ge" can be set by a loader # as can the threshold value. @@ -426,18 +424,22 @@ def match_varname_and_units(cubes: iris.cube.Cube | iris.cube.CubeList): returned. """ # If just one cube, then no need to match so return. - if len(cubes) == 1: + if isinstance(cubes, iris.cube.Cube): + cubes_in = iris.cube.CubeList([cubes]) + else: + cubes_in = cubes + if len(cubes_in) == 1: return cubes # Initialise the list of matched cubes. new_cubelist = iris.cube.CubeList([]) # Use the first cube in the CubeList as the base cube. - base_cube = cubes[0] + base_cube = cubes_in[0] new_cubelist.append(base_cube) # Loop over the cubes matching each to the base cube. - for cube in cubes[1:]: + for cube in cubes_in[1:]: new_cube = cube.copy() # Match the cube varname. @@ -448,19 +450,6 @@ def match_varname_and_units(cubes: iris.cube.Cube | iris.cube.CubeList): # Match the cube units. new_cube.units = base_cube.units - # Match the cube forecast_reference_time, a single value e.g. 2023-11-12 00:00:00. - forecast_reference_time_use = base_cube.coord("forecast_reference_time") - new_cube.remove_coord("forecast_reference_time") - new_cube.add_aux_coord(forecast_reference_time_use) - # new_cube.replace_coord(base_cube.coord('forecast_reference_time')) - - # Match the cube forecast_period - an array points: [ 1., 2., ..., 47., 48.]. - coord_dims = base_cube.coord_dims("forecast_period") - new_coord = base_cube.coord("forecast_period") - if new_cube.coords("forecast_period"): - new_cube.remove_coord("forecast_period") - new_cube.add_aux_coord(new_coord, coord_dims) - # Append the matched cube to the output cube list. new_cubelist.append(new_cube) From cdee4888855cb8561d90394d779f862d0a3513ee Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 24 Sep 2026 16:01:37 +0100 Subject: [PATCH 083/105] Added some radar tests #2113 --- tests/operators/test_radar_filter.py | 248 +++++++++++++++------------ 1 file changed, 135 insertions(+), 113 deletions(-) diff --git a/tests/operators/test_radar_filter.py b/tests/operators/test_radar_filter.py index 2d8b4fb964..ae6c496305 100644 --- a/tests/operators/test_radar_filter.py +++ b/tests/operators/test_radar_filter.py @@ -16,7 +16,7 @@ import iris import numpy as np -from iris.cube import Cube +import pytest from CSET.operators import constraints, radar_filter @@ -26,141 +26,152 @@ ) -def test_mask_list(): - """Test the mask_list function.""" - list_weights = radar_filter.mask_list(["UM_model", "Nimrod1km", "Nimrod2km"]) - expected_list = ["Nimrod2km_weights", "Nimrod1km_weights", "Nimrod2km_weights"] - assert list_weights == expected_list - - -def test_mask_by_weights(): - """Test the mask_by_weights function.""" - obs = Cube([[1.0, 2.0], [3.0, 4.0]]) - wts = Cube([[13, 0], [10, 11]]) - out = Cube([[1.0, np.nan], [np.nan, 4.0]]) - - var_name_obs = "hourly_rain_accumulation" - var_name_wts = "hourly_wts_accumulation" - - obs_name = "Radar" - wts_name = "Radar_wts" - obs.attributes["model_name"] = obs_name - wts.attributes["model_name"] = wts_name - - obs.rename(var_name_obs) - obs.long_name = var_name_obs - obs.var_name = var_name_obs - - wts.rename(var_name_wts) - wts.long_name = var_name_wts - wts.var_name = var_name_wts - - cube_list = iris.cube.CubeList([]) - cube_list.append(obs) - cube_list.append(wts) - - masked = radar_filter.mask_by_weights(cube_list, [obs_name], [wts_name]) - assert np.array_equal(masked.data, out.data, equal_nan=True) - - -def test_match_varname_and_units(): - """Test the match_varname_and_units function.""" - var_name_cube1 = "surface_microphysical_rainfall_rate" - var_name_cube3 = "Hourly rain accumulation" - - units_cube1 = "mm hr-1" - units_cube3 = "mm" - - # Construct test_cube1. - shape = (5, 4, 4) +# Session scope fixtures, so the test data only has to be loaded once. +@pytest.fixture(scope="session") +def cube_radar() -> iris.cube.Cube: + """Construct a radar rainfall accumulation cube.""" + shape = (1, 3, 3) data = np.ones(shape, dtype=np.float32) lat_coord = iris.coords.DimCoord( - np.linspace(-90, 90, 4), standard_name="latitude", units="degrees" + np.linspace(-90, 90, 3), standard_name="latitude", units="degrees" ) lon_coord = iris.coords.DimCoord( - np.linspace(-180, 180, 4), standard_name="longitude", units="degrees" + np.linspace(-180, 180, 3), standard_name="longitude", units="degrees" ) - time_points1 = np.array([0, 6, 12, 18, 24], dtype=np.float64) - time_coord1 = iris.coords.DimCoord( - time_points1, + time_points = np.array([3], dtype=np.float64) + time_coord = iris.coords.DimCoord( + time_points, standard_name="time", units="hours since 2026-09-23 00:00:00", ) - fp_coord1 = iris.coords.DimCoord( - time_points1, standard_name="forecast_period", units="hours" + fp_coord = iris.coords.DimCoord( + time_points, standard_name="forecast_period", units="hours" ) - test_cube1 = iris.cube.Cube( + test_cube = iris.cube.Cube( data, - long_name=var_name_cube1, - units=units_cube1, + long_name="Hourly rain accumulation", + units="mm", dim_coords_and_dims=[ - (time_coord1, 0), + (time_coord, 0), (lat_coord, 1), (lon_coord, 2), ], ) - test_cube1.add_aux_coord(fp_coord1, data_dims=0) + test_cube.add_aux_coord(fp_coord, data_dims=0) - frt_coord1 = iris.coords.AuxCoord( + frt_coord = iris.coords.AuxCoord( [0.0], standard_name="forecast_reference_time", units="hours since 2026-09-23 00:00:00", ) - test_cube1.add_aux_coord(frt_coord1) + test_cube.add_aux_coord(frt_coord) + + test_cube.attributes["model_name"] = "RadarA" + + new_data = np.array([[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]]) + test_cube.data = new_data + + return test_cube + + +@pytest.fixture(scope="session") +def cube_radar_wts(cube_radar) -> iris.cube.Cube: + """Construct a radar rainfall accumulation weights cube.""" + test_cube = cube_radar.copy() + new_data = np.array([[[13, 0, 13], [9, 12, 11], [13, 11, 10]]]) + test_cube.data = new_data + test_cube.attributes["model_name"] = "RadarA_wts" + test_cube.long_name = "Hourly wts accumulation" + test_cube.units = "1" + return test_cube + + +@pytest.fixture(scope="session") +def cube_radar_masked_by_wts(cube_radar) -> iris.cube.Cube: + """Construct a radar rainfall accumulation cube that has been masked by the weights.""" + test_cube = cube_radar.copy() + new_data = np.array([[[1.0, np.nan, 3.0], [np.nan, 5.0, 6.0], [7.0, 8.0, np.nan]]]) + test_cube.data = new_data + test_cube.attributes["mask"] = "mask_of_" + test_cube.long_name + return test_cube + + +@pytest.fixture(scope="session") +def cube_model(cube_radar) -> iris.cube.Cube: + """Construct a model rainfall rate cube.""" + test_cube = cube_radar.copy() + new_data = np.array([[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]]) + test_cube.data = new_data + test_cube.attributes["model_name"] = "ModelA" + test_cube.long_name = "surface_microphysical_rainfall_rate" + test_cube.units = "mm hr-1" + return test_cube + + +def test_print_cube_radar(cube_radar): + """Print the radar cube.""" + print("----> cube_radar ", cube_radar) + print() + print("----> cube_radar.data ", cube_radar.data) + print() + print( + '----> cube_radar.attributes["model_name"] ', + cube_radar.attributes["model_name"], + ) - # Construct test_cube2. - test_cube2 = test_cube1.copy() - # Construct test_cube3. - time_points3 = np.array([3, 9, 15, 21, 27], dtype=np.float64) +def test_print_cube_radar_mask(cube_radar_wts): + """Print the radar wts cube.""" + print("----> cube_radar_wts ", cube_radar_wts) + print() + print("----> cube_radar_wts.data ", cube_radar_wts.data) - time_coord3 = iris.coords.DimCoord( - time_points3, - standard_name="time", - units="hours since 2026-09-23 00:00:00", - ) - fp_coord3 = iris.coords.DimCoord( - time_points3, standard_name="forecast_period", units="hours" - ) +def test_print_cube_radar_masked_by_wts(cube_radar_masked_by_wts): + """Print the radar masked by wts cube.""" + print("----> cube_radar_masked_by_wts ", cube_radar_masked_by_wts) + print() + print("----> cube_radar_masked_by_wts.data ", cube_radar_masked_by_wts.data) - test_cube3 = iris.cube.Cube( - data, - long_name=var_name_cube3, - units=units_cube3, - dim_coords_and_dims=[ - (time_coord3, 0), - (lat_coord, 1), - (lon_coord, 2), - ], - ) - test_cube3.add_aux_coord(fp_coord3, data_dims=0) +def test_mask_list(): + """Test the mask_list function.""" + list_weights = radar_filter.mask_list(["UM_model", "Nimrod1km", "Nimrod2km"]) + expected_list = ["Nimrod2km_weights", "Nimrod1km_weights", "Nimrod2km_weights"] + assert list_weights == expected_list - frt_coord3 = iris.coords.AuxCoord( - [1.0], - standard_name="forecast_reference_time", - units="hours since 2026-09-23 00:00:00", - ) - test_cube3.add_aux_coord(frt_coord3) - print(" ") - print("-----> test_cube1 ", test_cube1) - print(" ") - print("-----> test_cube3 ", test_cube3) +def test_mask_by_weights(cube_radar, cube_radar_wts, cube_radar_masked_by_wts): + """Test the mask_by_weights function.""" + obs_name = cube_radar.attributes["model_name"] + wts_name = cube_radar_wts.attributes["model_name"] + cube_list = iris.cube.CubeList([]) + cube_list.append(cube_radar) + cube_list.append(cube_radar_wts) + masked = radar_filter.mask_by_weights(cube_list, [obs_name], [wts_name]) + print() + print("---->masked ", masked) + print() + print("----> masked.data ", masked.data) + assert np.array_equal(masked.data, cube_radar_masked_by_wts.data, equal_nan=True) - ######## + +def test_match_varname_and_units(cube_model, cube_radar): + """Test the match_varname_and_units function.""" + cube1 = cube_model.copy() + cube2 = cube_model.copy() + cube3 = cube_radar.copy() # Construct the CubeList to test. cube_list = iris.cube.CubeList([]) - cube_list.append(test_cube1) - cube_list.append(test_cube2) - cube_list.append(test_cube3) + cube_list.append(cube1) + cube_list.append(cube2) + cube_list.append(cube3) # Push the test CubeList through the function match_varname_and_units. cubes_matched = radar_filter.match_varname_and_units(cube_list) @@ -174,22 +185,33 @@ def test_match_varname_and_units(): assert cubes_matched[0].units == cubes_matched[1].units assert cubes_matched[0].units == cubes_matched[2].units - assert cubes_matched[0].coord("forecast_reference_time") == cubes_matched[1].coord( - "forecast_reference_time" - ) - assert cubes_matched[0].coord("forecast_reference_time") == cubes_matched[2].coord( - "forecast_reference_time" - ) - assert cubes_matched[0].coord("forecast_period") == cubes_matched[1].coord( - "forecast_period" - ) - assert cubes_matched[0].coord("forecast_period") == cubes_matched[2].coord( - "forecast_period" - ) +def test_match_varname_and_units_single_cube(cube_model): + """Test the match_varname_and_units function.""" + cube1 = cube_model.copy() + + # Push the test Cube through the function match_varname_and_units. + cubes_matched = radar_filter.match_varname_and_units(cube1) + + assert cubes_matched.var_name == cube1.var_name + assert cubes_matched.long_name == cube1.long_name + assert cubes_matched.units == cube1.units + + +def test_match_varname_and_units_single_cube_in_cubelist(cube_model): + """Test the match_varname_and_units function.""" + cube1 = cube_model.copy() + + # Construct the CubeList to test. + cube_list = iris.cube.CubeList([]) + cube_list.append(cube1) + + # Push the test CubeList through the function match_varname_and_units. + cubes_matched = radar_filter.match_varname_and_units(cube_list) - print(" ") - print("-----> cubes_matched[2] ", cubes_matched[2]) + assert cubes_matched[0].var_name == cube_model.var_name + assert cubes_matched[0].long_name == cube_model.long_name + assert cubes_matched[0].units == cube_model.units # def test_radar_apply_mask(): From ea30b526b0f32048106ee84894ddb4428b332526 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 1 Oct 2026 10:44:16 +0100 Subject: [PATCH 084/105] Added more tests to test_radar_filter.py #2113 --- tests/operators/test_radar_filter.py | 141 ++++++++++++++++++++++++++- 1 file changed, 140 insertions(+), 1 deletion(-) diff --git a/tests/operators/test_radar_filter.py b/tests/operators/test_radar_filter.py index ae6c496305..d69bcd97c5 100644 --- a/tests/operators/test_radar_filter.py +++ b/tests/operators/test_radar_filter.py @@ -26,6 +26,134 @@ ) +def make_cube(data, name, model_name=None): + """Create a cube with a model_name attribute when requested.""" + cube = iris.cube.Cube(np.asarray(data, dtype=float), long_name=name) + cube.rename(name) + if model_name is not None: + cube.attributes["model_name"] = model_name + return cube + + +class TestMaskList: + """Tests for the Nimrod mask selection helper.""" + + def test_returns_empty_list_when_no_nimrod_source_present(self): + """No Nimrod source should produce no mask names.""" + assert radar_filter.mask_list(["UM_model", "UKV_model"]) == [] + + def test_prefers_highest_resolution_nimrod_source(self): + """The preferred Nimrod weights should match the highest resolution source.""" + result = radar_filter.mask_list(["UM_model", "Nimrod1km", "Nimrod2km"]) + assert result == [ + "Nimrod2km_weights", + "Nimrod1km_weights", + "Nimrod2km_weights", + ] + + +class TestMaskByWeights: + """Tests for applying a radar weight mask to a field.""" + + def test_applies_mask_using_weight_threshold(self): + """Only cells with a valid Nimrod weight should remain in the field.""" + field = make_cube(np.arange(12, dtype=float).reshape(1, 3, 4), "rain_rate") + field.attributes["model_name"] = "UM_model" + + mask_values = np.array( + [ + [ + [0.0, 11.0, 12.0, 4.0], + [13.0, 0.0, 5.0, 14.0], + [10.0, 11.0, 0.0, 21.0], + ] + ], + dtype=float, + ) + weights = make_cube( + mask_values, "Nimrod2km_weights", model_name="Nimrod2km_weights" + ) + + result = radar_filter.mask_by_weights( + iris.cube.CubeList([field, weights]), + ["UM_model"], + ["Nimrod2km_weights"], + ) + + expected = np.full((1, 3, 4), np.nan, dtype=float) + valid = np.array( + [ + [ + [False, True, True, False], + [True, False, False, True], + [False, True, False, True], + ] + ], + dtype=bool, + ) + expected[valid] = field.data[valid] + + np.testing.assert_allclose(result.data, expected, equal_nan=True) + + +class TestRadarApplyMask: + """Tests for the boundary-aware radar mask application.""" + + def test_masks_boundary_points_and_applies_valid_mask(self): + """Boundary cells and mask zeros should both become NaN.""" + field = make_cube(np.arange(120, dtype=float).reshape(1, 10, 12), "rain_rate") + mask = make_cube(np.ones((1, 10, 12), dtype=float), "mask") + + result = radar_filter.radar_apply_mask(field, mask, boundary_margin=2) + + assert np.isnan(result.data[0, 0, 0]) + assert np.isnan(result.data[0, 9, 11]) + assert result.data[0, 2, 2] == 26.0 + assert result.data[0, 5, 5] == 65.0 + + +class TestRadarMask: + """Tests for masking radar and model data together.""" + + def test_returns_masked_radar_cube(self): + """Valid model points should propagate to the radar field output.""" + model_field = make_cube( + np.arange(120, dtype=float).reshape(1, 10, 12), "model_field" + ) + radar_field = make_cube( + np.arange(120, 240, dtype=float).reshape(1, 10, 12), "radar_field" + ) + mask = make_cube(np.ones((1, 10, 12), dtype=float), "nimrod_mask") + + result = radar_filter.radar_mask( + model_field, radar_field, mask, boundary_margin=2, outputs="radar" + ) + + assert isinstance(result, iris.cube.Cube) + assert np.isnan(result.data[0, 0, 0]) + assert np.isnan(result.data[0, 9, 11]) + assert result.data[0, 2, 2] == 146.0 + + def test_returns_both_model_and_radar_cubes_for_all_output(self): + """The all-output mode should return both model and radar cubes.""" + model_field = make_cube( + np.arange(120, dtype=float).reshape(1, 10, 12), "model_field" + ) + radar_field = make_cube( + np.arange(120, 240, dtype=float).reshape(1, 10, 12), "radar_field" + ) + mask = make_cube(np.ones((1, 10, 12), dtype=float), "nimrod_mask") + + result = radar_filter.radar_mask( + model_field, radar_field, mask, boundary_margin=2, outputs="all" + ) + + assert isinstance(result, iris.cube.CubeList) + assert len(result) == 2 + assert result[0].name() == "model_field" + assert result[1].name() == "radar_field" + + # Session scope fixtures, so the test data only has to be loaded once. @pytest.fixture(scope="session") def cube_radar() -> iris.cube.Cube: @@ -105,7 +233,7 @@ def cube_radar_masked_by_wts(cube_radar) -> iris.cube.Cube: def cube_model(cube_radar) -> iris.cube.Cube: """Construct a model rainfall rate cube.""" test_cube = cube_radar.copy() - new_data = np.array([[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]]) + new_data = np.array([[[11.0, 12.0, 13.0], [14.0, 15.0, 16.0], [17.0, 18.0, 19.0]]]) test_cube.data = new_data test_cube.attributes["model_name"] = "ModelA" test_cube.long_name = "surface_microphysical_rainfall_rate" @@ -161,6 +289,17 @@ def test_mask_by_weights(cube_radar, cube_radar_wts, cube_radar_masked_by_wts): assert np.array_equal(masked.data, cube_radar_masked_by_wts.data, equal_nan=True) +def test_radar_apply_mask(cube_model, cube_radar_wts): + """Test the mask_by_weights function.""" + masked = radar_filter.radar_apply_mask( + cube_model, cube_radar_wts, boundary_margin=0 + ) + print() + print("---->masked from test_radar_apply_mask", masked) + print() + print("----> masked.data ", masked.data) + + def test_match_varname_and_units(cube_model, cube_radar): """Test the match_varname_and_units function.""" cube1 = cube_model.copy() From 33a3ba05072259fdac8967238657013ce96af478 Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 1 Oct 2026 11:05:14 +0100 Subject: [PATCH 085/105] Added yet more tests to test_radar_filter.py #2113 --- tests/operators/test_radar_filter.py | 72 ++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) diff --git a/tests/operators/test_radar_filter.py b/tests/operators/test_radar_filter.py index d69bcd97c5..a57dcf8929 100644 --- a/tests/operators/test_radar_filter.py +++ b/tests/operators/test_radar_filter.py @@ -154,6 +154,78 @@ def test_returns_both_model_and_radar_cubes_for_all_output(self): assert result[1].name() == "radar_field" +def test_handles_multiple_field_and_mask_cubes(): + """Test the loop over both the mask and rainfall fields.""" + field1 = make_cube(np.ones((1, 3, 3), dtype=float), "field_a") + field2 = make_cube(np.full((1, 3, 3), 2.0, dtype=float), "field_b") + + mask1 = make_cube(np.ones((1, 3, 3), dtype=float), "mask_a") + mask2 = make_cube(np.ones((1, 3, 3), dtype=float), "mask_b") + + result = radar_filter.radar_apply_mask( + iris.cube.CubeList([field1, field2]), + iris.cube.CubeList([mask1, mask2]), + boundary_margin=0, + ) + + assert isinstance(result, iris.cube.CubeList) + assert len(result) == 2 + assert np.allclose(result[0].data, 1.0) + assert np.allclose(result[1].data, 2.0) + + +def test_returns_masked_model_cube_when_output_is_model(): + """Test the branch that executes if just model output is required.""" + model_field = make_cube( + np.arange(120, dtype=float).reshape(1, 10, 12), "model_field" + ) + radar_field = make_cube( + np.arange(120, 240, dtype=float).reshape(1, 10, 12), "radar_field" + ) + mask = make_cube(np.ones((1, 10, 12), dtype=float), "nimrod_mask") + + result = radar_filter.radar_mask( + model_field, radar_field, mask, boundary_margin=2, outputs="model" + ) + + assert isinstance(result, iris.cube.Cube) + assert result.name() == "model_field" + assert np.isnan(result.data[0, 0, 0]) + + +def test_radar_mask_loop_returns_model_and_radar_outputs(): + """Test the loop over all models.""" + model_field = iris.cube.CubeList( + [ + make_cube(np.ones((1, 3, 3), dtype=float), "model_1"), + make_cube(np.full((1, 3, 3), 2.0, dtype=float), "model_2"), + ] + ) + radar_field = make_cube(np.full((1, 3, 3), 9.0, dtype=float), "radar_field") + mask = make_cube(np.ones((1, 3, 3), dtype=float), "nimrod_mask") + + result = radar_filter.radar_mask_loop( + model_field, + radar_field, + mask, + boundary_margin=0, + outputs="radar", + ) + + assert isinstance(result, iris.cube.CubeList) + assert len(result) >= 2 + + +def test_boundary_margin_zero_keeps_all_inner_points(): + """Test zero width boundary on model field.""" + field = make_cube(np.arange(9, dtype=float).reshape(1, 3, 3), "rain_rate") + mask = make_cube(np.ones((1, 3, 3), dtype=float), "mask") + + result = radar_filter.radar_apply_mask(field, mask, boundary_margin=0) + + assert np.all(np.isfinite(result.data)) + + # Session scope fixtures, so the test data only has to be loaded once. @pytest.fixture(scope="session") def cube_radar() -> iris.cube.Cube: From 86d2784fd2c146747c16e0a5ef16346d0d26a5d0 Mon Sep 17 00:00:00 2001 From: James Ross Date: Fri, 16 Oct 2026 10:05:21 +0100 Subject: [PATCH 086/105] Apply suggestion from @Sylviabohnenstengel Co-authored-by: Sylvia Bohnenstengel <62748926+Sylviabohnenstengel@users.noreply.github.com> --- .../generic_surface_spatial_plot_sequence_radar_rainfall.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml b/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml index 77deba33f5..e3fbb35abc 100644 --- a/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml +++ b/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml @@ -2,7 +2,7 @@ category: Surface Spatial Plot title: "$RADAR_NAME $VARNAME $SUBAREA_NAME" description: | Extracts and plots radar rainfall data. - ATTENTION: the radar data can be 1 hour accumulations. Because showers + ATTENTION: the radar data can be 1 hour accumulation amounts in [mm]. Because showers typically advect, on a 2D plot of accumulated rainfall, showers show as streaks i.e. lower intensity rainfall over a greater area than a plot of instantaneous rainfall rate. This results in accumulated rainfall From 51fb232da5e5a1eb1e2f4e2dbc7adb3b2f2b5a85 Mon Sep 17 00:00:00 2001 From: James Ross Date: Fri, 16 Oct 2026 10:06:20 +0100 Subject: [PATCH 087/105] Update src/CSET/operators/radar_filter.py --- src/CSET/operators/radar_filter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index 1473f0ef41..d6818983cc 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -94,7 +94,7 @@ def mask_list(model_names: list[str]) -> list[str]: # Loop over the input model_names. for model in model_names: if any(model in nimrod for nimrod in nimrod_preference): - nimrod_mask = model + wei + nimrod_mask = f"{model}{wei}" else: nimrod_mask = preferred_nimrod + wei mask_names_list.append(nimrod_mask) From 9b4dcaaf5981fd5fc24d624a8d09af231b34e0de Mon Sep 17 00:00:00 2001 From: James Ross Date: Fri, 16 Oct 2026 10:07:11 +0100 Subject: [PATCH 088/105] Update src/CSET/loaders/radar.py Co-authored-by: Sylvia Bohnenstengel <62748926+Sylviabohnenstengel@users.noreply.github.com> --- src/CSET/loaders/radar.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index 7a74a1cc08..5c937dc245 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -140,7 +140,7 @@ def load(conf: Config): # Surface (2D) fields for Nimrod radar rainfall. # # These loaders produce 2D plots of both the Nimrod - # surface rainfall and the Nimrod weights field. + # surface rainfall amount in [mm] and the Nimrod weights field, which gives an indication of data quality. # # The different sources of Nimrod rainfall accumulation have # different spatial grids. So each source requires its own From b7535056036bbc009f6c56f4679f1b0435498948 Mon Sep 17 00:00:00 2001 From: James Ross Date: Fri, 16 Oct 2026 10:07:34 +0100 Subject: [PATCH 089/105] Update src/CSET/operators/radar_filter.py --- src/CSET/operators/radar_filter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index d6818983cc..d4595f23cc 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -96,7 +96,7 @@ def mask_list(model_names: list[str]) -> list[str]: if any(model in nimrod for nimrod in nimrod_preference): nimrod_mask = f"{model}{wei}" else: - nimrod_mask = preferred_nimrod + wei + nimrod_mask = f"{preferred_nimrod}{wei}" mask_names_list.append(nimrod_mask) return mask_names_list From d7b621c266cf4198d7ad33e3eee4079f145c1224 Mon Sep 17 00:00:00 2001 From: James Ross Date: Fri, 16 Oct 2026 10:09:09 +0100 Subject: [PATCH 090/105] type hinting --- src/CSET/loaders/radar.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index 5c937dc245..033ddb097a 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -17,7 +17,7 @@ from CSET.recipes import Config, RawRecipe, get_models -def select_radar_source(input_source_list): +def select_radar_source(input_source_list: list[str]) -> str: """Select the preferred radar source from a list of sources.""" # Preference list for Nimrod sources, with most preferred heading the list. nimrod_preference = ["Nimrod2km", "Nimrodxkm", "Nimrod1km"] From 3bfbe9c7af7dd31415e30e0d489a1d79e4c7f5ea Mon Sep 17 00:00:00 2001 From: James Ross Date: Fri, 16 Oct 2026 10:20:04 +0100 Subject: [PATCH 091/105] outputs -> required_outputs. match statement instead of if statements --- src/CSET/operators/radar_filter.py | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index d4595f23cc..16d9a012ff 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -248,7 +248,7 @@ def radar_mask( nimrod_field: iris.cube.Cube | iris.cube.CubeList, nimrod_mask: iris.cube.Cube | iris.cube.CubeList, boundary_margin: int = 8, - outputs: str = "radar", + required_outputs: str = "radar", ) -> iris.cube.Cube | iris.cube.CubeList: """Apply a mask to given fields using a masked array. @@ -263,7 +263,7 @@ def radar_mask( boundary_margin: int, optional Number of grid points from the domain boundary considered "unreliable". Defaults to 8. - outputs: str, optional + required_outputs: str, optional Specifies which outputs are required: "radar" outputs masked Nimrod rainfall field(s) "model" outputs masked model field(s) @@ -277,7 +277,7 @@ def radar_mask( Examples -------- To mask both model and Nimrod fields using the Nimrod weights in nimrod_wts: - >>> masked_fields = radar_mask( model_fields, nimrod_fields, nimrod_wts, outputs="all") + >>> masked_fields = radar_mask( model_fields, nimrod_fields, nimrod_wts, required_outputs="all") """ # Create an empty cubelist to hold the filtered fields. @@ -311,13 +311,14 @@ def radar_mask( filtered_radar.append(masked_nimrod_field) # Return either the masked model or Nimrod fields, or both. - if outputs == "radar": - filtered_fields.append(filtered_radar.merge_cube()) - if outputs == "model": - filtered_fields.append(filtered_model.merge_cube()) - if outputs == "all": - filtered_fields.append(filtered_model.merge_cube()) - filtered_fields.append(filtered_radar.merge_cube()) + match required_outputs: + case "radar": + filtered_fields.append(filtered_radar.merge_cube()) + case "model": + filtered_fields.append(filtered_model.merge_cube()) + case "all": + filtered_fields.append(filtered_model.merge_cube()) + filtered_fields.append(filtered_radar.merge_cube()) # Return either a single cube or a cubelist. if len(filtered_fields) == 1: @@ -377,7 +378,7 @@ def radar_mask_loop( use_nimrod_field, use_nimrod_mask, boundary_margin=boundary_margin, - outputs="model", + required_outputs="model", ) filtered_cubes.append(filtered_model) @@ -388,7 +389,7 @@ def radar_mask_loop( use_nimrod_field, use_nimrod_mask, boundary_margin=boundary_margin, - outputs="radar", + required_outputs="radar", ) else: filtered_radar = radar_mask( @@ -396,7 +397,7 @@ def radar_mask_loop( use_nimrod_field, use_nimrod_mask, boundary_margin=boundary_margin, - outputs="radar", + required_outputs="radar", ) filtered_cubes.append(filtered_radar) From 1a0bcb4b946fb659cb053f5faf8bd6c4eda70627 Mon Sep 17 00:00:00 2001 From: James Ross Date: Fri, 16 Oct 2026 10:40:02 +0100 Subject: [PATCH 092/105] simiplify radar_mask_loop --- src/CSET/operators/radar_filter.py | 30 ++++++------------------------ 1 file changed, 6 insertions(+), 24 deletions(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index 16d9a012ff..1bad2b9bad 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -332,7 +332,7 @@ def radar_mask_loop( nimrod_field: iris.cube.Cube | iris.cube.CubeList, nimrod_mask: iris.cube.Cube | iris.cube.CubeList, boundary_margin: int = 8, - outputs: str = "radar", + required_outputs: str = "all", ) -> iris.cube.Cube | iris.cube.CubeList: """Find common domains between a list of models and radar sources. @@ -347,7 +347,7 @@ def radar_mask_loop( boundary_margin: int, optional Number of grid points from the domain boundary considered "unreliable". Defaults to 8. - outputs: str, optional + required_outputs: str, optional Specifies which outputs are required: "radar" outputs masked Nimrod rainfall field(s) "model" outputs masked model field(s) @@ -361,7 +361,7 @@ def radar_mask_loop( Examples -------- To mask both model and Nimrod fields using the Nimrod weights in nimrod_wts: - >>> masked_fields = radar_mask_loop( model_fields, nimrod_fields, nimrod_wts, outputs="all") + >>> masked_fields = radar_mask_loop( model_fields, nimrod_fields, nimrod_wts, required_outputs="all") """ # Create an empty cubelist to hold the filtered fields. @@ -371,36 +371,18 @@ def radar_mask_loop( use_nimrod_mask = nimrod_mask # Loop over the models. + # using required_outputs="all" returns both + # masked model and radar fields for model in model_field: - print("-------> using model ", model) filtered_model = radar_mask( model, use_nimrod_field, use_nimrod_mask, boundary_margin=boundary_margin, - required_outputs="model", + required_outputs=required_outputs, ) filtered_cubes.append(filtered_model) - # Filter the radar observations. - if len(model_field) == 1: - filtered_radar = radar_mask( - model_field[0], - use_nimrod_field, - use_nimrod_mask, - boundary_margin=boundary_margin, - required_outputs="radar", - ) - else: - filtered_radar = radar_mask( - model_field[0], - use_nimrod_field, - use_nimrod_mask, - boundary_margin=boundary_margin, - required_outputs="radar", - ) - filtered_cubes.append(filtered_radar) - return filtered_cubes From d381a8e3f2bdbf3429b3252859b7feb9e0f472ac Mon Sep 17 00:00:00 2001 From: James Ross Date: Fri, 16 Oct 2026 10:44:05 +0100 Subject: [PATCH 093/105] fix unit tests --- tests/operators/test_radar_filter.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/operators/test_radar_filter.py b/tests/operators/test_radar_filter.py index a57dcf8929..4efa9f3c49 100644 --- a/tests/operators/test_radar_filter.py +++ b/tests/operators/test_radar_filter.py @@ -126,7 +126,7 @@ def test_returns_masked_radar_cube(self): mask = make_cube(np.ones((1, 10, 12), dtype=float), "nimrod_mask") result = radar_filter.radar_mask( - model_field, radar_field, mask, boundary_margin=2, outputs="radar" + model_field, radar_field, mask, boundary_margin=2, required_outputs="radar" ) assert isinstance(result, iris.cube.Cube) @@ -145,7 +145,7 @@ def test_returns_both_model_and_radar_cubes_for_all_output(self): mask = make_cube(np.ones((1, 10, 12), dtype=float), "nimrod_mask") result = radar_filter.radar_mask( - model_field, radar_field, mask, boundary_margin=2, outputs="all" + model_field, radar_field, mask, boundary_margin=2, required_outputs="all" ) assert isinstance(result, iris.cube.CubeList) @@ -185,7 +185,7 @@ def test_returns_masked_model_cube_when_output_is_model(): mask = make_cube(np.ones((1, 10, 12), dtype=float), "nimrod_mask") result = radar_filter.radar_mask( - model_field, radar_field, mask, boundary_margin=2, outputs="model" + model_field, radar_field, mask, boundary_margin=2, required_outputs="model" ) assert isinstance(result, iris.cube.Cube) @@ -209,7 +209,7 @@ def test_radar_mask_loop_returns_model_and_radar_outputs(): radar_field, mask, boundary_margin=0, - outputs="radar", + required_outputs="radar", ) assert isinstance(result, iris.cube.CubeList) From 465fc5b6b953780f88a2537f6c204ebd7e97259a Mon Sep 17 00:00:00 2001 From: James Ross Date: Fri, 16 Oct 2026 11:14:34 +0100 Subject: [PATCH 094/105] Update src/CSET/loaders/radar.py Co-authored-by: Sylvia Bohnenstengel <62748926+Sylviabohnenstengel@users.noreply.github.com> --- src/CSET/loaders/radar.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index 033ddb097a..01f4ef40d5 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -134,7 +134,7 @@ def load(conf: Config): if radar["varname"] == "Hourly rain accumulation" ] - # Specify the boundary margin width to use when trimming model fields. + # Specify the lateral boundary margin width to use when trimming model fields. boundary_margin = 16 # Surface (2D) fields for Nimrod radar rainfall. From d4911069c72a53609d62cd39abc9b2b7f28241ba Mon Sep 17 00:00:00 2001 From: James Ross Date: Fri, 16 Oct 2026 11:31:48 +0100 Subject: [PATCH 095/105] combine conditional --- src/CSET/loaders/radar.py | 71 +++++++++++++-------------------------- 1 file changed, 24 insertions(+), 47 deletions(-) diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index 01f4ef40d5..32d1417f86 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -229,50 +229,27 @@ def load(conf: Config): # Select the radar source to use. radar_source = select_radar_source([source["id"] for source in accum_radars]) names_ids = _model_and_radar_names_ids(models, radar_source, accum_radars) - yield RawRecipe( - recipe="radar_common_domain_histogram.yaml", - variables={ - "MODEL_VARNAME": "surface_microphysical_rainfall_rate", - "RADAR_VARNAME": "Hourly rain accumulation", - "RADAR_WTS_VARNAME": "Hourly wts accumulation", - "ALL_LABEL": names_ids["combined_names"], - "SEQUENCE": "time", - "OUTPUTS": "all", - "BOUNDARY_MARGIN": boundary_margin, - "TITLE_STRING": "sequence plots", - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, - "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", - }, - model_ids=names_ids["combined_ids"], - aggregation=False, - ) - - # Histogram case study rainfall using common domain between - # model and Nimrod radar observations. - if ( - conf.NIMROD_RADAR_OBS - and conf.PROCESS_RADAR_HISTOGRAMS - and (len(accum_radars) > 0) - ): - # Select the radar source to use. - radar_source = select_radar_source([source["id"] for source in accum_radars]) - names_ids = _model_and_radar_names_ids(models, radar_source, accum_radars) - yield RawRecipe( - recipe="radar_common_domain_histogram.yaml", - variables={ - "MODEL_VARNAME": "surface_microphysical_rainfall_rate", - "RADAR_VARNAME": "Hourly rain accumulation", - "RADAR_WTS_VARNAME": "Hourly wts accumulation", - "ALL_LABEL": names_ids["combined_names"], - "SEQUENCE": "realization", - "OUTPUTS": "all", - "BOUNDARY_MARGIN": boundary_margin, - "TITLE_STRING": "case study", - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, - "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", - }, - model_ids=names_ids["combined_ids"], - aggregation=False, - ) + for sequence, title in ( + ("time", "sequence plots"), + ("realization", "case study"), + ): + yield RawRecipe( + recipe="radar_common_domain_histogram.yaml", + variables={ + "MODEL_VARNAME": "surface_microphysical_rainfall_rate", + "RADAR_VARNAME": "Hourly rain accumulation", + "RADAR_WTS_VARNAME": "Hourly wts accumulation", + "ALL_LABEL": names_ids["combined_names"], + "SEQUENCE": sequence, + "OUTPUTS": "all", + "BOUNDARY_MARGIN": boundary_margin, + "TITLE_STRING": title, + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT + if conf.SELECT_SUBAREA + else None, + "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", + }, + model_ids=names_ids["combined_ids"], + aggregation=False, + ) From 660b7766668ab8c4d4fa4597f839acc197fb82b0 Mon Sep 17 00:00:00 2001 From: James Ross Date: Fri, 16 Oct 2026 11:59:02 +0100 Subject: [PATCH 096/105] fix uneven boudnary --- src/CSET/operators/radar_filter.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index 1bad2b9bad..bb175aab64 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -213,12 +213,12 @@ def radar_apply_mask( # Set the model perimeter to NaN as these gridpoints contain no useful data. # c.f. boundary_margin in regrid.py - margin_width = boundary_margin - if margin_width > 0: - masked_field.data[:, -margin_width - 1 :, :] = np.nan - masked_field.data[:, :, -margin_width - 1 :] = np.nan - masked_field.data[:, :, 0:margin_width] = np.nan - masked_field.data[:, 0:margin_width, :] = np.nan + if boundary_margin > 0: + # Nan a strip of width boundary_margin around all four edges + masked_field.data[:, -boundary_margin:, :] = np.nan + masked_field.data[:, :, :boundary_margin] = np.nan + masked_field.data[:, :, -boundary_margin:] = np.nan + masked_field.data[:, :boundary_margin, :] = np.nan # If the field and mask are on different grids, then regrid the field. if M[0].shape != masked_field[0].shape: From e597460a70abcb0065edcdb31db5d121d491e24c Mon Sep 17 00:00:00 2001 From: James Ross Date: Fri, 16 Oct 2026 12:32:37 +0100 Subject: [PATCH 097/105] make mask boundary seperate function and unit test it --- src/CSET/operators/radar_filter.py | 29 +++++++++++++++----- tests/operators/test_radar_filter.py | 40 ++++++++++++++++++++++++++++ 2 files changed, 63 insertions(+), 6 deletions(-) diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py index bb175aab64..cdfcd4188a 100644 --- a/src/CSET/operators/radar_filter.py +++ b/src/CSET/operators/radar_filter.py @@ -213,12 +213,7 @@ def radar_apply_mask( # Set the model perimeter to NaN as these gridpoints contain no useful data. # c.f. boundary_margin in regrid.py - if boundary_margin > 0: - # Nan a strip of width boundary_margin around all four edges - masked_field.data[:, -boundary_margin:, :] = np.nan - masked_field.data[:, :, :boundary_margin] = np.nan - masked_field.data[:, :, -boundary_margin:] = np.nan - masked_field.data[:, :boundary_margin, :] = np.nan + _mask_boundary(masked_field.data, boundary_margin) # If the field and mask are on different grids, then regrid the field. if M[0].shape != masked_field[0].shape: @@ -243,6 +238,28 @@ def radar_apply_mask( return masked_fields.merge() +def _mask_boundary(data: np.ndarray, boundary_margin: int) -> None: + """Put a NaN strip along margin of width boundary_margin. + + Parameters + ---------- + data: array + The data array to be masked. + boundary_margin + Number of points to set as NaN at each edge. + """ + if boundary_margin < 0: + raise ValueError(f"boundary_margin must be >= 0, got {boundary_margin}") + if boundary_margin == 0: + return + + m = boundary_margin + data[..., :m, :] = np.nan # first rows + data[..., -m:, :] = np.nan # last rows + data[..., :m] = np.nan # first columns + data[..., -m:] = np.nan # last columns + + def radar_mask( model_field: iris.cube.Cube | iris.cube.CubeList, nimrod_field: iris.cube.Cube | iris.cube.CubeList, diff --git a/tests/operators/test_radar_filter.py b/tests/operators/test_radar_filter.py index 4efa9f3c49..e6d9c55bd2 100644 --- a/tests/operators/test_radar_filter.py +++ b/tests/operators/test_radar_filter.py @@ -425,6 +425,46 @@ def test_match_varname_and_units_single_cube_in_cubelist(cube_model): assert cubes_matched[0].units == cube_model.units +def test_mask_boundary_zero_margin_leaves_data_unchanged(): + """Test the mask_boundary function with boundary_margin set to 0.""" + data = np.zeros((2, 5, 6)) + radar_filter._mask_boundary(data, 0) + assert not np.isnan(data).any() + + +def test_mask_boundary_blanks_edges_only(): + """Test the mask_boundary function with boundary_margin set to 1.""" + data = np.zeros((5, 6)) + radar_filter._mask_boundary(data, 1) + expected_nan = np.ones((5, 6), dtype=bool) + expected_nan[1:-1, 1:-1] = False + np.testing.assert_array_equal(np.isnan(data), expected_nan) + + +def test_mask_boundary_is_symmetric(): + """Test the mask_boundary function with boundary_margin set to 3.""" + data = np.zeros((11, 11)) + radar_filter._mask_boundary(data, 3) + nan_rows = np.isnan(data[:, 5]) + np.testing.assert_array_equal(np.flatnonzero(nan_rows), [0, 1, 2, 8, 9, 10]) + + +def test_mask_boundary_applies_to_all_leading_dims(): + """Test the mask_boundary function with boundary_margin set to 1 and with time + realisation coords.""" + data = np.zeros((3, 2, 5, 5)) # e.g. (time, realization, y, x) + radar_filter._mask_boundary(data, 1) + for t in range(3): + for r in range(2): + assert np.isnan(data[t, r, 0, :]).all() + assert not np.isnan(data[t, r, 2, 2]) + + +def test_mask_boundary_negative_margin_raises(): + """Test bad boundary margin raises.""" + with pytest.raises(ValueError): + radar_filter._mask_boundary(np.zeros((4, 4)), -1) + + # def test_radar_apply_mask(): # """Test the radar_apply_mask function.""" # obs = Cube([0,[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]]) From 65ed6d7873fd32a55a32a900b74ee2ac3321a694 Mon Sep 17 00:00:00 2001 From: James Ross Date: Fri, 16 Oct 2026 13:07:45 +0100 Subject: [PATCH 098/105] move radar recipes to own folder --- .../generic_surface_spatial_plot_sequence_radar_rainfall.yaml | 0 .../{ => nimrod}/model_radar_common_domain_2d_model.yaml | 0 .../{ => nimrod}/model_radar_common_domain_2d_radar.yaml | 0 .../surface_fields/{ => nimrod}/radar_common_domain_2d_plot.yaml | 0 .../{ => nimrod}/radar_common_domain_histogram.yaml | 0 .../{ => nimrod}/radar_common_domain_histogram_domain_mean.yaml | 0 .../{ => nimrod}/radar_common_domain_histogram_sequence.yaml | 0 src/CSET/recipes/surface_fields/{ => nimrod}/radar_dev3.yaml | 0 .../{ => nimrod}/radar_masked_mean_time_series.yaml | 0 .../surface_fields/{ => nimrod}/radar_mean_time_series.yaml | 0 10 files changed, 0 insertions(+), 0 deletions(-) rename src/CSET/recipes/surface_fields/{ => nimrod}/generic_surface_spatial_plot_sequence_radar_rainfall.yaml (100%) rename src/CSET/recipes/surface_fields/{ => nimrod}/model_radar_common_domain_2d_model.yaml (100%) rename src/CSET/recipes/surface_fields/{ => nimrod}/model_radar_common_domain_2d_radar.yaml (100%) rename src/CSET/recipes/surface_fields/{ => nimrod}/radar_common_domain_2d_plot.yaml (100%) rename src/CSET/recipes/surface_fields/{ => nimrod}/radar_common_domain_histogram.yaml (100%) rename src/CSET/recipes/surface_fields/{ => nimrod}/radar_common_domain_histogram_domain_mean.yaml (100%) rename src/CSET/recipes/surface_fields/{ => nimrod}/radar_common_domain_histogram_sequence.yaml (100%) rename src/CSET/recipes/surface_fields/{ => nimrod}/radar_dev3.yaml (100%) rename src/CSET/recipes/surface_fields/{ => nimrod}/radar_masked_mean_time_series.yaml (100%) rename src/CSET/recipes/surface_fields/{ => nimrod}/radar_mean_time_series.yaml (100%) diff --git a/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml b/src/CSET/recipes/surface_fields/nimrod/generic_surface_spatial_plot_sequence_radar_rainfall.yaml similarity index 100% rename from src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml rename to src/CSET/recipes/surface_fields/nimrod/generic_surface_spatial_plot_sequence_radar_rainfall.yaml diff --git a/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml b/src/CSET/recipes/surface_fields/nimrod/model_radar_common_domain_2d_model.yaml similarity index 100% rename from src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml rename to src/CSET/recipes/surface_fields/nimrod/model_radar_common_domain_2d_model.yaml diff --git a/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml b/src/CSET/recipes/surface_fields/nimrod/model_radar_common_domain_2d_radar.yaml similarity index 100% rename from src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml rename to src/CSET/recipes/surface_fields/nimrod/model_radar_common_domain_2d_radar.yaml diff --git a/src/CSET/recipes/surface_fields/radar_common_domain_2d_plot.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_2d_plot.yaml similarity index 100% rename from src/CSET/recipes/surface_fields/radar_common_domain_2d_plot.yaml rename to src/CSET/recipes/surface_fields/nimrod/radar_common_domain_2d_plot.yaml diff --git a/src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram.yaml similarity index 100% rename from src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml rename to src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram.yaml diff --git a/src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_domain_mean.yaml similarity index 100% rename from src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml rename to src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_domain_mean.yaml diff --git a/src/CSET/recipes/surface_fields/radar_common_domain_histogram_sequence.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_sequence.yaml similarity index 100% rename from src/CSET/recipes/surface_fields/radar_common_domain_histogram_sequence.yaml rename to src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_sequence.yaml diff --git a/src/CSET/recipes/surface_fields/radar_dev3.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_dev3.yaml similarity index 100% rename from src/CSET/recipes/surface_fields/radar_dev3.yaml rename to src/CSET/recipes/surface_fields/nimrod/radar_dev3.yaml diff --git a/src/CSET/recipes/surface_fields/radar_masked_mean_time_series.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_masked_mean_time_series.yaml similarity index 100% rename from src/CSET/recipes/surface_fields/radar_masked_mean_time_series.yaml rename to src/CSET/recipes/surface_fields/nimrod/radar_masked_mean_time_series.yaml diff --git a/src/CSET/recipes/surface_fields/radar_mean_time_series.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_mean_time_series.yaml similarity index 100% rename from src/CSET/recipes/surface_fields/radar_mean_time_series.yaml rename to src/CSET/recipes/surface_fields/nimrod/radar_mean_time_series.yaml From b2642c3f1e5b9f79af9d1b8bba62e5a02dc2dbd7 Mon Sep 17 00:00:00 2001 From: James Ross Date: Fri, 16 Oct 2026 13:38:35 +0100 Subject: [PATCH 099/105] removing some redundant code --- src/CSET/loaders/radar.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index 32d1417f86..2478a9636c 100644 --- a/src/CSET/loaders/radar.py +++ b/src/CSET/loaders/radar.py @@ -96,7 +96,7 @@ def get_radar_sources(conf) -> list[dict]: return radar_sources -def _model_and_radar_names_ids(models, radar_source, accum_radars): +def _model_and_radar_names_ids(models, accum_radars): """Select the model and radar names and ids.""" # Model and radar names and ids. model_names_list = [model["name"] for model in models] @@ -165,8 +165,7 @@ def load(conf: Config): # Radar 2D plots using a common domain between model and radar network. if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D and (len(accum_radars) > 0): - radar_source = select_radar_source([source["id"] for source in accum_radars]) - names_ids = _model_and_radar_names_ids(models, radar_source, accum_radars) + names_ids = _model_and_radar_names_ids(models, accum_radars) yield RawRecipe( recipe="model_radar_common_domain_2d_radar.yaml", variables={ @@ -191,8 +190,7 @@ def load(conf: Config): # Model rainfall 2D plots using a common domain between model and radar network. if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D and (len(accum_radars) > 0): - radar_source = select_radar_source([source["id"] for source in accum_radars]) - names_ids = _model_and_radar_names_ids(models, radar_source, accum_radars) + names_ids = _model_and_radar_names_ids(models, accum_radars) # Loop over the model names. for idx, model_use in enumerate(names_ids["model_names_list"]): yield RawRecipe( @@ -226,9 +224,7 @@ def load(conf: Config): and conf.PROCESS_RADAR_HISTOGRAMS and (len(accum_radars) > 0) ): - # Select the radar source to use. - radar_source = select_radar_source([source["id"] for source in accum_radars]) - names_ids = _model_and_radar_names_ids(models, radar_source, accum_radars) + names_ids = _model_and_radar_names_ids(models, accum_radars) for sequence, title in ( ("time", "sequence plots"), ("realization", "case study"), From df5aace08cdfa3db7a1a7fcef4b185de6cf5fab8 Mon Sep 17 00:00:00 2001 From: James Ross Date: Fri, 16 Oct 2026 15:20:34 +0100 Subject: [PATCH 100/105] Update src/CSET/recipes/surface_fields/nimrod/radar_common_domain_2d_plot.yaml Co-authored-by: Carol Halliwell <150140850+cehalliwell@users.noreply.github.com> --- .../surface_fields/nimrod/radar_common_domain_2d_plot.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_2d_plot.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_2d_plot.yaml index 7cb809bfab..c89a06ffca 100644 --- a/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_2d_plot.yaml +++ b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_2d_plot.yaml @@ -48,7 +48,7 @@ steps: condition: 'ge' value: 11 - boundary_margin: 16 + boundary_margin: $BOUNDARY_MARGIN outputs: $OUTPUTS - operator: collapse.collapse From fc74417aadb24f864d9255cf9a471e3dd30ae350 Mon Sep 17 00:00:00 2001 From: James Ross Date: Fri, 16 Oct 2026 15:20:48 +0100 Subject: [PATCH 101/105] Update src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_domain_mean.yaml Co-authored-by: Carol Halliwell <150140850+cehalliwell@users.noreply.github.com> --- .../nimrod/radar_common_domain_histogram_domain_mean.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_domain_mean.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_domain_mean.yaml index 2229807878..6bf45ecb77 100644 --- a/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_domain_mean.yaml +++ b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_domain_mean.yaml @@ -47,7 +47,7 @@ steps: condition: 'ge' value: 11 - boundary_margin: 16 + boundary_margin: $BOUNDARY_MARGIN outputs: $OUTPUTS - operator: radar_filter.match_varname_and_units From 245f35637be37afabf358e3477f943080d94d32e Mon Sep 17 00:00:00 2001 From: James Ross Date: Fri, 16 Oct 2026 15:21:44 +0100 Subject: [PATCH 102/105] Update src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_sequence.yaml Co-authored-by: Carol Halliwell <150140850+cehalliwell@users.noreply.github.com> --- .../nimrod/radar_common_domain_histogram_sequence.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_sequence.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_sequence.yaml index b03359916e..b4ebbfa3f6 100644 --- a/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_sequence.yaml +++ b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_sequence.yaml @@ -47,7 +47,7 @@ steps: condition: 'ge' value: 11 - boundary_margin: 16 + boundary_margin: $BOUNDARY_MARGIN outputs: $OUTPUTS - operator: radar_filter.match_varname_and_units From 59c954626c8803623ed5e07a83505a04e53cc8b5 Mon Sep 17 00:00:00 2001 From: James Ross Date: Tue, 20 Oct 2026 14:30:40 +0100 Subject: [PATCH 103/105] code cleanup --- .../app/fetch_nimrod/bin/fetch_nimrod.py | 71 +++++++++---------- 1 file changed, 35 insertions(+), 36 deletions(-) diff --git a/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py b/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py index 8386664348..9e3af2c0f8 100755 --- a/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py +++ b/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py @@ -85,6 +85,14 @@ def apply_radar_weights(cube_obs: iris.cube.Cube, cube_wei: iris.cube.Cube): return cube_obs_weighted, cube_wei +def _load_nimrod(path: Path) -> iris.cube.Cube | None: + """Load a Nimrod file, or return None if it's missing.""" + if not path.exists(): + logger.warning("Nimrod file not found: %s", path) + return None + return iris.load_cube(path) + + def retrieve_nimrod(): """Fetch the observations corresponding to a model run. @@ -106,17 +114,15 @@ def retrieve_nimrod(): for nimrod_field in v["field"]: if nimrod_field: logger.info("Processing Nimrod field: %s", nimrod_field) - + nimrod_cfg = nimrod_dict[nimrod_field] # Prepare the output directory for the Nimrod field. - nimrod_dir = f"{v['rose_datac']}/data/{nimrod_dict[nimrod_field]['obs_id']}" + nimrod_dir = Path(v["rose_datac"]) / "data" / nimrod_cfg["obs_id"] os.makedirs(nimrod_dir, exist_ok=True) logger.info("Cylc-run Nimrod directory: %s", nimrod_dir) # Prepare the output directory for the Nimrod weights field. if nimrod_dict[nimrod_field]["weights_fname"]: - nimrod_dir_wei = ( - f"{v['rose_datac']}/data/{nimrod_dict[nimrod_field]['wei_id']}" - ) + nimrod_dir_wei = Path(v["rose_datac"]) / "data" / nimrod_cfg["wei_id"] os.makedirs(nimrod_dir_wei, exist_ok=True) logger.info("Cylc-run Nimrod weights directory: %s", nimrod_dir_wei) @@ -136,48 +142,41 @@ def retrieve_nimrod(): date_use = date_start_offset while date_use < date_start_offset + v["forecast_length"]: # Load the Nimrod data into an Iris cube. - nimrod_obs_exist = "False" - nimrod_obs = ( - f"{nimrod_dict[nimrod_field]['basedir']}/" - f"{nimrod_dict[nimrod_field]['obs_dir']}/" - f"{date_use.year}/{date_use.strftime('%Y%m%d%H%M')}" - f"{nimrod_dict[nimrod_field]['obs_fname']}" + + nimrod_obs_path = ( + Path(nimrod_cfg["basedir"]) + / nimrod_cfg["obs_dir"] + / str(date_use.year) + / f"{date_use:%Y%m%d%H%M}{nimrod_cfg['obs_fname']}" ) - try: - nimrod_cube_obs = iris.load_cube(nimrod_obs) - nimrod_obs_exist = "True" - except OSError: - logger.warning("Iris cannot find Nimrod file %s", nimrod_obs) + + nimrod_cube_obs = _load_nimrod(nimrod_obs_path) # Load the Nimrod weights into an Iris cube. - nimrod_weights_exist = "False" - if nimrod_dict[nimrod_field]["weights_fname"] != "": - nimrod_weights = ( - f"{nimrod_dict[nimrod_field]['basedir']}/" - f"{nimrod_dict[nimrod_field]['weights_dir']}/" - f"{date_use.year}/{date_use.strftime('%Y%m%d%H%M')}" - f"{nimrod_dict[nimrod_field]['weights_fname']}" + nimrod_cube_weights = None + if nimrod_cfg["weights_fname"]: + nimrod_weights_path = ( + Path(nimrod_cfg["basedir"]) + / nimrod_cfg["weights_dir"] + / str(date_use.year) + / f"{date_use:%Y%m%d%H%M}{nimrod_cfg['weights_fname']}" ) - try: - nimrod_cube_weights = iris.load_cube(nimrod_weights) - nimrod_weights_exist = "True" - except OSError: - logger.warning( - "Iris cannot find Nimrod weights file %s", nimrod_weights - ) - # Ensure that the weights are unitless. - nimrod_cube_weights.units = "1" + nimrod_cube_weights = _load_nimrod(nimrod_weights_path) + + if nimrod_cube_weights is not None: + # Ensure that the weights are unitless. + nimrod_cube_weights.units = "1" # QC the the Nimrod observations using the weights field. - if nimrod_obs_exist == "True" and nimrod_weights_exist == "True": + if nimrod_cube_obs is not None and nimrod_cube_weights is not None: nimrod_cube_obs, nimrod_cube_weights = apply_radar_weights( nimrod_cube_obs, nimrod_cube_weights ) # Now that the accumulation weights have been applied if required, # write the Nimrod obs to a NetCDF file. - if nimrod_obs_exist == "True": + if nimrod_cube_obs is not None: filename_obs_nc = ( f"{nimrod_dir}/{date_use.strftime('%Y%m%d%H%M')}" f"_{nimrod_dict[nimrod_field]['obs_id']}" @@ -187,7 +186,7 @@ def retrieve_nimrod(): # Write the Nimrod weights to a NetCDF file if the switch # for this is set in v["weights"]. - if v["weights"] == "True" and nimrod_weights_exist == "True": + if v["weights"] == "True" and nimrod_cube_weights is not None: filename_wei_nc = ( f"{nimrod_dir_wei}/{date_use.strftime('%Y%m%d%H%M')}" f"_{nimrod_dict[nimrod_field]['obs_id']}_weights" @@ -202,6 +201,6 @@ def retrieve_nimrod(): ) -# Run the function that fetches the NImrod radar obs. +# Run the function that fetches the Nimrod radar obs. if __name__ == "__main__": retrieve_nimrod() From 66e37dfe2c4506ca7025f3c5d0c82d4a524a53bb Mon Sep 17 00:00:00 2001 From: James Ross Date: Tue, 20 Oct 2026 15:15:38 +0100 Subject: [PATCH 104/105] fix test --- .../app/fetch_nimrod/bin/fetch_nimrod.py | 16 ++++++++-------- tests/test_fetch_nimrod.py | 10 ++++++++-- 2 files changed, 16 insertions(+), 10 deletions(-) diff --git a/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py b/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py index 9e3af2c0f8..d3b338d6fb 100755 --- a/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py +++ b/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py @@ -85,14 +85,6 @@ def apply_radar_weights(cube_obs: iris.cube.Cube, cube_wei: iris.cube.Cube): return cube_obs_weighted, cube_wei -def _load_nimrod(path: Path) -> iris.cube.Cube | None: - """Load a Nimrod file, or return None if it's missing.""" - if not path.exists(): - logger.warning("Nimrod file not found: %s", path) - return None - return iris.load_cube(path) - - def retrieve_nimrod(): """Fetch the observations corresponding to a model run. @@ -102,7 +94,15 @@ def retrieve_nimrod(): * ANALYSIS_LENGTH * ROSE_DATAC """ + # Grab the environment variables required for handling Nimrod data. + def _load_nimrod(path: Path) -> iris.cube.Cube | None: + """Load a Nimrod file, or return None if it's missing.""" + if not path.exists(): + logger.warning("Nimrod file not found: %s", path) + return None + return iris.load_cube(path) + v = _get_needed_environment_variables_nimrod() # Grab the Nimrod handling dictionary. diff --git a/tests/test_fetch_nimrod.py b/tests/test_fetch_nimrod.py index 5a4ec2fc34..d0337a5e8b 100644 --- a/tests/test_fetch_nimrod.py +++ b/tests/test_fetch_nimrod.py @@ -84,11 +84,12 @@ def test_get_needed_environment_variables_nimrod(monkeypatch): def test_retrieve_nimrod_saves_observation_and_weights(tmp_path, monkeypatch): """Test retrieve_nimrod saves both observation and weight files when available.""" # Create a temporary nimrod dictionary file. + basedir = tmp_path / "nimrod" nimrod_dict = { "Nimrod_comp_xkm": { "obs_id": "nimrod_xkm", "wei_id": "nimrod_xkm_weights", - "basedir": "/data/nimrod", + "basedir": str(basedir), "obs_dir": "obs", "obs_fname": "_obs.nc", "weights_dir": "weights", @@ -96,9 +97,14 @@ def test_retrieve_nimrod_saves_observation_and_weights(tmp_path, monkeypatch): "freq": "PT01H", } } + + for subdir, fname in [("obs", "_obs.nc"), ("weights", "_weights.nc")]: + f = basedir / subdir / "2026" / f"202604071300{fname}" + f.parent.mkdir(parents=True) + f.touch() + config_path = tmp_path / "restricted_nimrod_met_office.json" config_path.write_text(json.dumps(nimrod_dict)) - # Set the required environment variables. monkeypatch.setenv("NIMROD_COMP_XKM", "True") monkeypatch.setenv("NIMROD_COMP_1KM", "False") From bce263e29796d4ad38600a2c4f670e8ac5cbc074 Mon Sep 17 00:00:00 2001 From: James Ross Date: Thu, 22 Oct 2026 08:58:52 +0100 Subject: [PATCH 105/105] outputs -> required_outputs --- .../surface_fields/nimrod/radar_common_domain_2d_plot.yaml | 2 +- .../surface_fields/nimrod/radar_common_domain_histogram.yaml | 2 +- .../nimrod/radar_common_domain_histogram_domain_mean.yaml | 2 +- .../nimrod/radar_common_domain_histogram_sequence.yaml | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_2d_plot.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_2d_plot.yaml index c89a06ffca..4aab5a7768 100644 --- a/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_2d_plot.yaml +++ b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_2d_plot.yaml @@ -49,7 +49,7 @@ steps: value: 11 boundary_margin: $BOUNDARY_MARGIN - outputs: $OUTPUTS + required_outputs: $OUTPUTS - operator: collapse.collapse coordinate: [time] diff --git a/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram.yaml index 00e2614961..b6e4ed506e 100644 --- a/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram.yaml +++ b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram.yaml @@ -57,7 +57,7 @@ steps: value: 11 boundary_margin: $BOUNDARY_MARGIN - outputs: $OUTPUTS + required_outputs: $OUTPUTS - operator: radar_filter.match_varname_and_units diff --git a/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_domain_mean.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_domain_mean.yaml index 6bf45ecb77..5102312e06 100644 --- a/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_domain_mean.yaml +++ b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_domain_mean.yaml @@ -48,7 +48,7 @@ steps: value: 11 boundary_margin: $BOUNDARY_MARGIN - outputs: $OUTPUTS + required_outputs: $OUTPUTS - operator: radar_filter.match_varname_and_units diff --git a/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_sequence.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_sequence.yaml index b4ebbfa3f6..3c9116b9b9 100644 --- a/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_sequence.yaml +++ b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_sequence.yaml @@ -48,7 +48,7 @@ steps: value: 11 boundary_margin: $BOUNDARY_MARGIN - outputs: $OUTPUTS + required_outputs: $OUTPUTS - operator: radar_filter.match_varname_and_units