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 7dc01377b0..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 @@ -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 @@ -101,10 +101,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"]: @@ -124,9 +120,21 @@ def retrieve_nimrod(): os.makedirs(nimrod_dir_wei, exist_ok=True) logger.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 = ( diff --git a/src/CSET/cset_workflow/meta/observations/rose-meta.conf b/src/CSET/cset_workflow/meta/observations/rose-meta.conf index b786240ac2..e21174c3d8 100644 --- a/src/CSET/cset_workflow/meta/observations/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/observations/rose-meta.conf @@ -223,8 +223,8 @@ title=Radar observations [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 @@ -269,6 +269,31 @@ 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 +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-9 + + + + ################################################################ # Cardington UK Research Observations [ns=Observations/Cardington] diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index ae74148a0a..549ea4872a 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]: @@ -85,6 +98,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,12 +111,22 @@ 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 + # 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.SPATIAL_SURFACE_FIELD: + 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", @@ -118,6 +144,175 @@ def load(conf: Config): aggregation=False, ) + # Histogram sequence rainfall using common domain between + # model and Nimrod radar observations. + 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] + 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.yaml", + variables={ + "MODEL_VARNAME": "surface_microphysical_rainfall_rate", + "RADAR_VARNAME": "Hourly rain accumulation", + "RADAR_WTS_VARNAME": "Hourly wts accumulation", + "ALL_LABEL": combined_names, + "SEQUENCE": "time", + "OUTPUTS": "all", + "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=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: + # 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.yaml", + variables={ + "MODEL_VARNAME": "surface_microphysical_rainfall_rate", + "RADAR_VARNAME": "Hourly rain accumulation", + "RADAR_WTS_VARNAME": "Hourly wts accumulation", + "ALL_LABEL": combined_names, + "SEQUENCE": "realization", + "OUTPUTS": "all", + "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=combined_ids, + 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, + ) + + # Radar masking of radar obs based on sea mask. + if conf.SPATIAL_SURFACE_FIELD and conf.NIMROD_RADAR_OBS: + # 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: + 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 @@ -135,6 +330,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 +390,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 01ccb8e575..2b9e4971cf 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", "scoreswrappers", diff --git a/src/CSET/operators/_colormaps.py b/src/CSET/operators/_colormaps.py index 9d37fc5db1..3061a3ca74 100644 --- a/src/CSET/operators/_colormaps.py +++ b/src/CSET/operators/_colormaps.py @@ -534,7 +534,11 @@ def custom_colormap_precipitation(cube: iris.cube.Cube, cmap, levels, norm): cmap = mcolors.ListedColormap(colors) # Normalize the levels norm = mcolors.BoundaryNorm(levels, cmap.N) - logger.info("Using custom rainfall colourmap.") + logging.info("Using custom rainfall colourmap.") + + # Set any Nan values to be plotted a light grey. + cmap.set_bad("#dcdcdc") + return cmap, levels, norm @@ -567,6 +571,7 @@ def custom_colourmap_nimrod_weights(cube: iris.cube.Cube, cmap, levels, norm): ] norm = mcolors.BoundaryNorm(levels, cmap.N) colours = [ + "#dcdcdc", "#d10000", "purple", "#8f00d6", @@ -579,7 +584,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/collapse.py b/src/CSET/operators/collapse.py index 44e4b48af5..99ee68c7cb 100644 --- a/src/CSET/operators/collapse.py +++ b/src/CSET/operators/collapse.py @@ -76,6 +76,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: logger.debug("Extracting common time points as multiple model inputs detected.") @@ -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: diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index 6138512eaa..ed07d736da 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -752,7 +752,7 @@ def _plot_and_save_spatial_plot( # In the bbox dictionary, fc and ec are hex colour codes for grey shade. cube_min, cube_max, cube_mean = calc_array_stats(cube.data) axes.annotate( - f"Min: {cube_min:.3g} Max: {cube_max:.3g} Mean: {cube_mean:.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), @@ -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, @@ -1904,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 @@ -1966,6 +1970,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 +2060,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. @@ -2081,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 ------- @@ -2100,6 +2110,7 @@ def spatial_pcolormesh_plot( filename, sequence_coordinate, stamp_coordinate, + strict_title, **kwargs, ) elif isinstance(cubes, iris.cube.Cube): @@ -2109,6 +2120,7 @@ def spatial_pcolormesh_plot( filename, sequence_coordinate, stamp_coordinate, + strict_title, **kwargs, ) return cubes diff --git a/src/CSET/operators/radar_filter.py b/src/CSET/operators/radar_filter.py new file mode 100644 index 0000000000..75bca16543 --- /dev/null +++ b/src/CSET/operators/radar_filter.py @@ -0,0 +1,540 @@ +# © 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 +import numpy as np + +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 = mask_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, + boundary_margin: int = 8, +) -> 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). + boundary_margin: int, optional + Number of grid points from the domain boundary considered "unreliable". + Defaults to 8. + + 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() + + # 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") + 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() + + +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("@@@@@@@@@@@@@@@@@@@@@ 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] + 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("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@") + 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) + 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 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..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: $RADAR_NAME $VARNAME $SUBAREA_NAME +title: "$RADAR_NAME $VARNAME $SUBAREA_NAME" description: Extracts and plots radar rainfall data. steps: 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..63e26d8873 --- /dev/null +++ b/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_model.yaml @@ -0,0 +1,62 @@ +category: Surface Spatial Plot +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. + 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 + strict_title: True + sequence_coordinate: time 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..fa5291075e --- /dev/null +++ b/src/CSET/recipes/surface_fields/model_radar_common_domain_2d_radar.yaml @@ -0,0 +1,62 @@ +category: Surface Spatial Plot +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. + 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 + strict_title: True + sequence_coordinate: time 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 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..b2c98a1728 --- /dev/null +++ b/src/CSET/recipes/surface_fields/radar_common_domain_histogram.yaml @@ -0,0 +1,61 @@ +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. + 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 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..2229807878 --- /dev/null +++ b/src/CSET/recipes/surface_fields/radar_common_domain_histogram_domain_mean.yaml @@ -0,0 +1,65 @@ +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] + # coordinate: [time] + # method: MEAN + + - operator: write.write_cube_to_nc + overwrite: True + + - operator: plot.plot_histogram_series + sequence_coordinate: "realization" + #sequence_coordinate: $SEQUENCE 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 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..4d62b52abf --- /dev/null +++ b/src/CSET/recipes/surface_fields/radar_dev3.yaml @@ -0,0 +1,49 @@ +category: Histogram +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). +# +# 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 0000000000..f723f53033 --- /dev/null +++ b/src/CSET/recipes/surface_fields/radar_masked_mean_time_series.yaml @@ -0,0 +1,43 @@ +category: Histogram +title: "bmc5 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/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: diff --git a/tests/operators/test_colormaps.py b/tests/operators/test_colormaps.py index 7b4ff6ff69..d9a0e902d4 100644 --- a/tests/operators/test_colormaps.py +++ b/tests/operators/test_colormaps.py @@ -463,6 +463,7 @@ def test_colorbar_map_nimrod_wts(cube, tmp_working_dir): cube.rename("Hourly wts accumulation") expected_levels = np.arange(-0.5, 14.5, 1.0) expected_colors = [ + "#dcdcdc", "#d10000", "purple", "#8f00d6", @@ -475,7 +476,6 @@ def test_colorbar_map_nimrod_wts(cube, tmp_working_dir): "#37a648", "#8edc64", "#c5ffc5", - "#dcdcdc", "#ffffff", ] expected_cmap = mpl.colors.ListedColormap(expected_colors)