diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index ae74148a0..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) @@ -95,6 +98,98 @@ 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" + ] + + # # 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 @@ -135,6 +230,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. + 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 "", @@ -161,6 +290,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( 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/_colormaps.py b/src/CSET/operators/_colormaps.py index 6a6e3800c..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 @@ -544,6 +548,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 +561,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/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), diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py new file mode 100644 index 000000000..8f5473417 --- /dev/null +++ b/src/CSET/operators/radar_filter.py @@ -0,0 +1,259 @@ +# © 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, + # 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 + + +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 = radar_apply_mask( surface_microphysical_rainfall_rate, Nimrod2km) + """ + # 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): + masked_field = F.copy() + + # 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") + masked_field = masked_field.regrid(M, scheme) + + # 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: + # return masked_fields + return masked_fields.merge() 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 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..83184b886 --- /dev/null +++ b/src/CSET/recipes/surface_fields/radar_masked_mean_time_series.yaml @@ -0,0 +1,43 @@ +category: Histogram +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 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)