From 0b52c83ab3a7af8b67b4b3087fea4d0243ed7eee Mon Sep 17 00:00:00 2001 From: Bernard Claxton Date: Thu, 4 Jun 2026 09:34:52 +0100 Subject: [PATCH 01/13] 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 000000000..82683d7e3 --- /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 02/13] 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 ae74148a0..5f53462d6 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 03/13] 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 634506dd6..94637911c 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 000000000..4350099f9 --- /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 04/13] 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 000000000..dd1f31f2d --- /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 05/13] 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 5f53462d6..f46f26e24 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 06/13] 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 6a6e3800c..f3444ac7c 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 513f9d0f7..2a8e792ae 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 07/13] 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 dd1f31f2d..83184b886 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 08/13] 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 4350099f9..2d4a4c3d0 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 09/13] 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 2d4a4c3d0..5268b176d 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 10/13] 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 f46f26e24..bd330303f 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 11/13] 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 a8a72a51e..501176f3d 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 12/13] 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 f3444ac7c..dd042043e 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 13/13] 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 5268b176d..8f5473417 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])