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..d3b338d6fb 100755 --- a/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py +++ b/src/CSET/cset_workflow/app/fetch_nimrod/bin/fetch_nimrod.py @@ -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 @@ -94,82 +94,89 @@ def retrieve_nimrod(): * ANALYSIS_LENGTH * ROSE_DATAC """ + # Grab the environment variables required for handling Nimrod data. + def _load_nimrod(path: Path) -> iris.cube.Cube | None: + """Load a Nimrod file, or return None if it's missing.""" + if not path.exists(): + logger.warning("Nimrod file not found: %s", path) + return None + return iris.load_cube(path) + v = _get_needed_environment_variables_nimrod() # Grab the Nimrod handling dictionary. 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"]: if nimrod_field: logger.info("Processing Nimrod field: %s", nimrod_field) - + nimrod_cfg = nimrod_dict[nimrod_field] # Prepare the output directory for the Nimrod field. - nimrod_dir = f"{v['rose_datac']}/data/{nimrod_dict[nimrod_field]['obs_id']}" + nimrod_dir = Path(v["rose_datac"]) / "data" / nimrod_cfg["obs_id"] os.makedirs(nimrod_dir, exist_ok=True) logger.info("Cylc-run Nimrod directory: %s", nimrod_dir) # Prepare the output directory for the Nimrod weights field. if nimrod_dict[nimrod_field]["weights_fname"]: - nimrod_dir_wei = ( - f"{v['rose_datac']}/data/{nimrod_dict[nimrod_field]['wei_id']}" - ) + nimrod_dir_wei = Path(v["rose_datac"]) / "data" / nimrod_cfg["wei_id"] os.makedirs(nimrod_dir_wei, exist_ok=True) logger.info("Cylc-run Nimrod weights directory: %s", nimrod_dir_wei) + # 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 = ( - f"{nimrod_dict[nimrod_field]['basedir']}/" - f"{nimrod_dict[nimrod_field]['obs_dir']}/" - f"{date_use.year}/{date_use.strftime('%Y%m%d%H%M')}" - f"{nimrod_dict[nimrod_field]['obs_fname']}" + + nimrod_obs_path = ( + Path(nimrod_cfg["basedir"]) + / nimrod_cfg["obs_dir"] + / str(date_use.year) + / f"{date_use:%Y%m%d%H%M}{nimrod_cfg['obs_fname']}" ) - try: - nimrod_cube_obs = iris.load_cube(nimrod_obs) - nimrod_obs_exist = "True" - except OSError: - logger.warning("Iris cannot find Nimrod file %s", nimrod_obs) + + nimrod_cube_obs = _load_nimrod(nimrod_obs_path) # Load the Nimrod weights into an Iris cube. - nimrod_weights_exist = "False" - if nimrod_dict[nimrod_field]["weights_fname"] != "": - nimrod_weights = ( - f"{nimrod_dict[nimrod_field]['basedir']}/" - f"{nimrod_dict[nimrod_field]['weights_dir']}/" - f"{date_use.year}/{date_use.strftime('%Y%m%d%H%M')}" - f"{nimrod_dict[nimrod_field]['weights_fname']}" + nimrod_cube_weights = None + if nimrod_cfg["weights_fname"]: + nimrod_weights_path = ( + Path(nimrod_cfg["basedir"]) + / nimrod_cfg["weights_dir"] + / str(date_use.year) + / f"{date_use:%Y%m%d%H%M}{nimrod_cfg['weights_fname']}" ) - try: - nimrod_cube_weights = iris.load_cube(nimrod_weights) - nimrod_weights_exist = "True" - except OSError: - logger.warning( - "Iris cannot find Nimrod weights file %s", nimrod_weights - ) - # Ensure that the weights are unitless. - nimrod_cube_weights.units = "1" + nimrod_cube_weights = _load_nimrod(nimrod_weights_path) + + if nimrod_cube_weights is not None: + # Ensure that the weights are unitless. + nimrod_cube_weights.units = "1" # QC the the Nimrod observations using the weights field. - if nimrod_obs_exist == "True" and nimrod_weights_exist == "True": + if nimrod_cube_obs is not None and nimrod_cube_weights is not None: nimrod_cube_obs, nimrod_cube_weights = apply_radar_weights( nimrod_cube_obs, nimrod_cube_weights ) # Now that the accumulation weights have been applied if required, # write the Nimrod obs to a NetCDF file. - if nimrod_obs_exist == "True": + if nimrod_cube_obs is not None: filename_obs_nc = ( f"{nimrod_dir}/{date_use.strftime('%Y%m%d%H%M')}" f"_{nimrod_dict[nimrod_field]['obs_id']}" @@ -179,7 +186,7 @@ def retrieve_nimrod(): # Write the Nimrod weights to a NetCDF file if the switch # for this is set in v["weights"]. - if v["weights"] == "True" and nimrod_weights_exist == "True": + if v["weights"] == "True" and nimrod_cube_weights is not None: filename_wei_nc = ( f"{nimrod_dir_wei}/{date_use.strftime('%Y%m%d%H%M')}" f"_{nimrod_dict[nimrod_field]['obs_id']}_weights" @@ -194,6 +201,6 @@ def retrieve_nimrod(): ) -# Run the function that fetches the NImrod radar obs. +# Run the function that fetches the Nimrod radar obs. if __name__ == "__main__": retrieve_nimrod() diff --git a/src/CSET/cset_workflow/meta/observations/rose-meta.conf b/src/CSET/cset_workflow/meta/observations/rose-meta.conf index b786240ac2..674ee63c96 100644 --- a/src/CSET/cset_workflow/meta/observations/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/observations/rose-meta.conf @@ -111,6 +111,8 @@ trigger=template variables=NIMROD_COMP_XKM: True; template variables=NIMROD_COMP_2KM: True; template variables=NIMROD_WEIGHTS: True; template variables=NIMROD_COMP_5MIN: True; + template variables=PROCESS_RADAR_2D: True; + template variables=PROCESS_RADAR_HISTOGRAMS: True; [template variables=CARDINGTON_OBS] ns=Observations @@ -223,8 +225,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 +271,28 @@ 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 the Rose tab Diagnostics/Surface (2D) fields. +compulsory=true +type=python_boolean +sort-key=nimrod-radar-9a + ################################################################ # Cardington UK Research Observations [ns=Observations/Cardington] diff --git a/src/CSET/cset_workflow/rose-suite.conf.example b/src/CSET/cset_workflow/rose-suite.conf.example index af76fd0172..a9095bc8e5 100644 --- a/src/CSET/cset_workflow/rose-suite.conf.example +++ b/src/CSET/cset_workflow/rose-suite.conf.example @@ -104,6 +104,8 @@ MODERATE_RAIN_PRESENCE_SPATIAL_PLOT=False !!NIMROD_COMP_2KM=False !!NIMROD_COMP_5MIN=False !!NIMROD_COMP_XKM=False +!!PROCESS_RADAR_2D=False +!!PROCESS_RADAR_HISTOGRAMS=False NIMROD_RADAR_OBS=False !!NIMROD_WEIGHTS=False !!PLEVEL_TRANSECT_AGGREGATION=False,False,False,False diff --git a/src/CSET/loaders/radar.py b/src/CSET/loaders/radar.py index ae74148a0a..2478a9636c 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: list[str]) -> str: + """Select the preferred radar source from a list of sources.""" + # Preference list for Nimrod sources, with most preferred heading the list. + nimrod_preference = ["Nimrod2km", "Nimrodxkm", "Nimrod1km"] + + 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]: @@ -83,8 +96,34 @@ def get_radar_sources(conf) -> list[dict]: return radar_sources +def _model_and_radar_names_ids(models, accum_radars): + """Select the model and radar names and ids.""" + # Model and radar names and ids. + model_names_list = [model["name"] for model in models] + model_ids_list = [model["id"] for model in models] + radar_source = select_radar_source([source["id"] for source in accum_radars]) + radar_obs_ids = [radar_source] + radar_wts_ids = [radar_source + "_weights"] + + # Construct the dictionary to return. + model_and_radar = { + "radar_source": radar_source, + "radar_obs_ids": radar_obs_ids, + "radar_wts_ids": radar_wts_ids, + "model_names_list": model_names_list, + "model_ids_list": model_ids_list, + "combined_names": model_names_list + radar_obs_ids + radar_wts_ids, + "combined_ids": model_ids_list + radar_obs_ids + radar_wts_ids, + } + + return model_and_radar + + def load(conf: Config): """Yield recipes from the given workflow configuration.""" + # Load a list of model detail dictionaries. + models = get_models(conf.asdict()) + # Load the required radar observation sources. radar_sources = get_radar_sources(conf) @@ -95,12 +134,18 @@ def load(conf: Config): if radar["varname"] == "Hourly rain accumulation" ] + # Specify the lateral boundary margin width to use when trimming model fields. + boundary_margin = 16 + # Surface (2D) fields for Nimrod radar rainfall. # + # These loaders produce 2D plots of both the Nimrod + # surface rainfall amount in [mm] and the Nimrod weights field, which gives an indication of data quality. + # # The different sources of Nimrod rainfall accumulation have # different spatial grids. So each source requires its own # recipe to prevent incompatible cubes being created. - if conf.SPATIAL_SURFACE_FIELD: + if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D and (len(radar_sources) > 0): for radar in radar_sources: yield RawRecipe( recipe="generic_surface_spatial_plot_sequence_radar_rainfall.yaml", @@ -118,77 +163,89 @@ def load(conf: Config): aggregation=False, ) - # Histograms for surface (2D) Nimrod radar hourly accumulated rainfall. - # - # To get multiple radar sources plotted on the histogram the - # recipe must be done by passing lists of the radar_ids and - # the radar_names. As this is a multiline plot, all radar sources - # share the same radar variable name. - if conf.HISTOGRAM_SURFACE_FIELD and accum_radars: + # Radar 2D plots using a common domain between model and radar network. + if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D and (len(accum_radars) > 0): + names_ids = _model_and_radar_names_ids(models, accum_radars) yield RawRecipe( - recipe="generic_surface_histogram_series.yaml", - # model_ids -> Becomes $INPUT_PATHS - model_ids=[radar["id"] for radar in accum_radars], + recipe="model_radar_common_domain_2d_radar.yaml", variables={ - "VARNAME": next(radar["varname"] for radar in accum_radars), - "MODEL_NAME": [radar["name"] for radar in accum_radars], - "SEQUENCE": "time" - if conf.HISTOGRAM_SURFACE_FIELD_SEQUENCE - else "realization", + "MODEL_VARNAME": "surface_microphysical_rainfall_rate", + "RADAR_VARNAME": "Hourly rain accumulation", + "MODEL_LABEL": names_ids["model_names_list"][0], + "RADAR_LABEL": names_ids["radar_obs_ids"][0], + "MASK_LABEL": names_ids["radar_wts_ids"][0], + "METHOD": "SEQ", + "BOUNDARY_MARGIN": boundary_margin, "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", }, + model_ids=[ + names_ids["model_ids_list"][0], + names_ids["radar_obs_ids"][0], + names_ids["radar_wts_ids"][0], + ], aggregation=False, ) - # Histograms for surface (2D) Nimrod radar 5 minute rainfall rate. - if conf.HISTOGRAM_SURFACE_FIELD and conf.NIMROD_COMP_5MIN: - yield RawRecipe( - recipe="generic_surface_histogram_series.yaml", - # model_ids -> Becomes $INPUT_PATHS - model_ids="Nimrod5min", - variables={ - "VARNAME": "Rainfall rate Composite", - "MODEL_NAME": "Nimrod_5min", - "SEQUENCE": "time" - if conf.HISTOGRAM_SURFACE_FIELD_SEQUENCE - else "realization", - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, - "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", - }, - aggregation=False, - ) - - # Timeseries plot of Nimrod hourly surface rainfall accumulation. - if conf.TIMESERIES_SURFACE_FIELD and accum_radars: - yield RawRecipe( - recipe="radar_mean_time_series.yaml", - # model_ids -> Becomes $INPUT_PATHS - model_ids=[radar["id"] for radar in accum_radars], - variables={ - "VARNAME": next(radar["varname"] for radar in accum_radars), - "MODEL_NAME": [radar["name"] for radar in accum_radars], - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, - "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", - }, - aggregation=False, - ) + # Model rainfall 2D plots using a common domain between model and radar network. + if conf.NIMROD_RADAR_OBS and conf.PROCESS_RADAR_2D and (len(accum_radars) > 0): + names_ids = _model_and_radar_names_ids(models, accum_radars) + # Loop over the model names. + for idx, model_use in enumerate(names_ids["model_names_list"]): + yield RawRecipe( + recipe="model_radar_common_domain_2d_model.yaml", + variables={ + "MODEL_VARNAME": "surface_microphysical_rainfall_rate", + "RADAR_VARNAME": "Hourly rain accumulation", + "MODEL_LABEL": model_use, + "RADAR_LABEL": names_ids["radar_obs_ids"][0], + "MASK_LABEL": names_ids["radar_wts_ids"][0], + "METHOD": "SEQ", + "BOUNDARY_MARGIN": boundary_margin, + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT + if conf.SELECT_SUBAREA + else None, + "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", + }, + model_ids=[ + names_ids["model_ids_list"][idx], + names_ids["radar_obs_ids"][0], + names_ids["radar_wts_ids"][0], + ], + aggregation=False, + ) - # Timeseries plot of Nimrod 5 minute rainfall rate. - if conf.TIMESERIES_SURFACE_FIELD and conf.NIMROD_COMP_5MIN: - yield RawRecipe( - recipe="radar_mean_time_series.yaml", - # model_ids -> Becomes $INPUT_PATHS - model_ids="Nimrod5min", - variables={ - "VARNAME": "Rainfall rate Composite", - "MODEL_NAME": "Nimrod_5min", - "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, - "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None, - "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", - }, - aggregation=False, - ) + # Histogram sequence rainfall using common domain between + # model and Nimrod radar observations. + if ( + conf.NIMROD_RADAR_OBS + and conf.PROCESS_RADAR_HISTOGRAMS + and (len(accum_radars) > 0) + ): + names_ids = _model_and_radar_names_ids(models, accum_radars) + for sequence, title in ( + ("time", "sequence plots"), + ("realization", "case study"), + ): + yield RawRecipe( + recipe="radar_common_domain_histogram.yaml", + variables={ + "MODEL_VARNAME": "surface_microphysical_rainfall_rate", + "RADAR_VARNAME": "Hourly rain accumulation", + "RADAR_WTS_VARNAME": "Hourly wts accumulation", + "ALL_LABEL": names_ids["combined_names"], + "SEQUENCE": sequence, + "OUTPUTS": "all", + "BOUNDARY_MARGIN": boundary_margin, + "TITLE_STRING": title, + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT + if conf.SELECT_SUBAREA + else None, + "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", + }, + model_ids=names_ids["combined_ids"], + aggregation=False, + ) 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..9e19025615 100644 --- a/src/CSET/operators/_colormaps.py +++ b/src/CSET/operators/_colormaps.py @@ -535,6 +535,10 @@ def custom_colormap_precipitation(cube: iris.cube.Cube, cmap, levels, norm): # Normalize the levels norm = mcolors.BoundaryNorm(levels, cmap.N) logger.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/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..cdfcd4188a --- /dev/null +++ b/src/CSET/operators/radar_filter.py @@ -0,0 +1,456 @@ +# © 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 containing model names and at least one Nimrod hourly + rainfall accumulation source. + Possible radar sources are: + "Nimrod2km", "Nimrod_2km". + "Nimrodxkm", "Nimrod_xkm". + "Nimrod1km", "Nimrod_1km". + + Returns + ------- + list[str] + A list of the Nimrod weights files to use with each of the input + models and radar sources. + + Notes + ----- + At least one of the entries in the input list must be a Nimrod hourly + rainfall accumulation source. + + If just one Nimrod source is specified, then the weights file associated + with this source is used. + + If more than one Nimrod source is in the input list, then each of the these + Nimrod sources is associated with its own weights file e.g. if the input list + contains ["Nimrod1km", "Nimrod2km"] then the weights files for these will + be ["Nimrod1km_weights", "Nimrod2km_weights"]. Any model fields in the input + list will be allocated a weights file according to the order of preference + 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 = f"{model}{wei}" + else: + nimrod_mask = f"{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 radar weights field as a mask. + + Parameters + ---------- + cubes: iris.cube.CubeList + CubeList containing fields to mask and radar weights fields to use as the masks. + model_names: list[str] + A list of the model_names or radar sources to mask. + weights_names: list[str] + A list of radar weights sources to use as masks. There should be an entry + in weights_names to correspond with every entry in model_names. + + Returns + ------- + CubeList: + A CubeList of masked fields. + + Examples + -------- + >>> field_filtered = mask_by_weights(cubelist, model_names) + + """ + # 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. + for model, mask in zip( + iter_maybe(model_names), + iter_maybe(weights_names), + strict=True, + ): + # Grab the field to filter. + model_constraint = iris.AttributeConstraint(model_name=model) + unfiltered_field = cubes.extract_cube(model_constraint) + + # Select the field to use as the mask. + # Nice to do - put in support for a static mask. + mask_constraint = iris.AttributeConstraint(model_name=mask) + mask_field = cubes.extract_cube(mask_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) + + # 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 field 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 filters.generate_mask, you can combine multiple masks in a + recipe using other functions before applying the mask to the data. + + Examples + -------- + >>> radar_domain_only = radar_apply_mask( surface_microphysical_rainfall_rate, Nimrod2km_wts) + """ + # Create an empty cubelist to hold the filtered fields. + masked_fields = iris.cube.CubeList([]) + + # 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() + + # Set the model perimeter to NaN as these gridpoints contain no useful data. + # c.f. boundary_margin in regrid.py + _mask_boundary(masked_field.data, boundary_margin) + + # If the field and mask are on different grids, then regrid the field. + if M[0].shape != masked_field[0].shape: + 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 an attribute to the masked field detailing the mask used. + masked_field.attributes["mask"] = f"mask_of_{F.name()}" + + # Append the masked field to the output list of masked fields. + 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 _mask_boundary(data: np.ndarray, boundary_margin: int) -> None: + """Put a NaN strip along margin of width boundary_margin. + + Parameters + ---------- + data: array + The data array to be masked. + boundary_margin + Number of points to set as NaN at each edge. + """ + if boundary_margin < 0: + raise ValueError(f"boundary_margin must be >= 0, got {boundary_margin}") + if boundary_margin == 0: + return + + m = boundary_margin + data[..., :m, :] = np.nan # first rows + data[..., -m:, :] = np.nan # last rows + data[..., :m] = np.nan # first columns + data[..., -m:] = np.nan # last columns + + +def radar_mask( + model_field: iris.cube.Cube | iris.cube.CubeList, + nimrod_field: iris.cube.Cube | iris.cube.CubeList, + nimrod_mask: iris.cube.Cube | iris.cube.CubeList, + boundary_margin: int = 8, + required_outputs: str = "radar", +) -> iris.cube.Cube | iris.cube.CubeList: + """Apply a mask to given fields using a masked array. + + Parameters + ---------- + model_field: iris.cube.Cube | iris.cube.CubeList + The model field(s) to be masked. + nimrod_field: iris.cube.Cube | iris.cube.CubeList + The Nimrod field(s) to be masked. + nimrod_mask: iris.cube.Cube | iris.cube.CubeList + The Nimrod mask(s) to use. These are normally Nimrod wts fields. + boundary_margin: int, optional + Number of grid points from the domain boundary considered "unreliable". + Defaults to 8. + required_outputs: str, optional + Specifies which outputs are required: + "radar" outputs masked Nimrod rainfall field(s) + "model" outputs masked model field(s) + "all" outputs both masked model and masked Nimrod field(s). + + Returns + ------- + masked_field: iris.cube.Cube | iris.cube.CubeList + A cube or CubeList of the masked field(s). + + Examples + -------- + To mask both model and Nimrod fields using the Nimrod weights in nimrod_wts: + >>> masked_fields = radar_mask( model_fields, nimrod_fields, nimrod_wts, required_outputs="all") + + """ + # Create an empty cubelist to hold the filtered fields. + filtered_fields = iris.cube.CubeList([]) + filtered_radar = iris.cube.CubeList([]) + filtered_model = iris.cube.CubeList([]) + + # Loop over nimrod_field, model_field and nimrod_mask. + for M, F, N in zip( + iter_maybe(nimrod_mask), + iter_maybe(model_field), + iter_maybe(nimrod_field), + strict=True, + ): + # Apply the function radar_apply_mask to generate the re-gridded + # and masked model field. + masked_model_field = radar_apply_mask(F, M, boundary_margin=boundary_margin) + + # Use the masked model field as the mask for the Nimrod field. + # Note: no re-gridding required. + min_timesteps = min(N.shape[0], masked_model_field.shape[0]) + + temp_mask = masked_model_field[0:min_timesteps].copy() + temp_mask.data[~np.isnan(temp_mask.data)] = 1.0 + + masked_nimrod_field = N[0:min_timesteps].copy() + 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 either the masked model or Nimrod fields, or both. + match required_outputs: + case "radar": + filtered_fields.append(filtered_radar.merge_cube()) + case "model": + filtered_fields.append(filtered_model.merge_cube()) + case "all": + filtered_fields.append(filtered_model.merge_cube()) + filtered_fields.append(filtered_radar.merge_cube()) + + # Return either a single cube or a cubelist. + if len(filtered_fields) == 1: + return filtered_fields[0] + else: + 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, + required_outputs: str = "all", +) -> iris.cube.Cube | iris.cube.CubeList: + """Find common domains between a list of models and radar sources. + + Parameters + ---------- + model_field: iris.cube.Cube | iris.cube.CubeList + The model field(s) to be masked. + nimrod_field: iris.cube.Cube | iris.cube.CubeList + The Nimrod field(s) to be masked. + nimrod_mask: iris.cube.Cube | iris.cube.CubeList + The Nimrod mask(s) to use. These are normally Nimrod wts fields. + boundary_margin: int, optional + Number of grid points from the domain boundary considered "unreliable". + Defaults to 8. + required_outputs: str, optional + Specifies which outputs are required: + "radar" outputs masked Nimrod rainfall field(s) + "model" outputs masked model field(s) + "all" outputs both masked model and masked Nimrod field(s). + + Returns + ------- + masked_field: iris.cube.Cube | iris.cube.CubeList + A cube or CubeList of the masked field(s). + + Examples + -------- + To mask both model and Nimrod fields using the Nimrod weights in nimrod_wts: + >>> masked_fields = radar_mask_loop( model_fields, nimrod_fields, nimrod_wts, required_outputs="all") + + """ + # Create an empty cubelist to hold the filtered fields. + filtered_cubes = iris.cube.CubeList([]) + + use_nimrod_field = nimrod_field + use_nimrod_mask = nimrod_mask + + # Loop over the models. + # using required_outputs="all" returns both + # masked model and radar fields + for model in model_field: + filtered_model = radar_mask( + model, + use_nimrod_field, + use_nimrod_mask, + boundary_margin=boundary_margin, + required_outputs=required_outputs, + ) + filtered_cubes.append(filtered_model) + + 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 isinstance(cubes, iris.cube.Cube): + cubes_in = iris.cube.CubeList([cubes]) + else: + cubes_in = cubes + if len(cubes_in) == 1: + return cubes + + # Initialise the list of matched cubes. + new_cubelist = iris.cube.CubeList([]) + + # Use the first cube in the CubeList as the base cube. + base_cube = cubes_in[0] + new_cubelist.append(base_cube) + + # Loop over the cubes matching each to the base cube. + for cube in cubes_in[1:]: + new_cube = cube.copy() + + # 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 + + # Append the matched cube to the output cube list. + new_cubelist.append(new_cube) + + 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 deleted file mode 100644 index 66c8606d8e..0000000000 --- a/src/CSET/recipes/surface_fields/generic_surface_spatial_plot_sequence_radar_rainfall.yaml +++ /dev/null @@ -1,20 +0,0 @@ -category: Surface Spatial Plot -title: $RADAR_NAME $VARNAME $SUBAREA_NAME -description: Extracts and plots radar rainfall data. - -steps: - - operator: read.read_cube - file_paths: $INPUT_PATHS - constraint: - operator: constraints.combine_constraints - varname_constraint: - operator: constraints.generate_var_constraint - varname: $VARNAME - subarea_type: $SUBAREA_TYPE - subarea_extent: $SUBAREA_EXTENT - - - operator: plot.spatial_pcolormesh_plot - sequence_coordinate: time - - - operator: write.write_cube_to_nc - overwrite: True diff --git a/src/CSET/recipes/surface_fields/nimrod/generic_surface_spatial_plot_sequence_radar_rainfall.yaml b/src/CSET/recipes/surface_fields/nimrod/generic_surface_spatial_plot_sequence_radar_rainfall.yaml new file mode 100644 index 0000000000..e3fbb35abc --- /dev/null +++ b/src/CSET/recipes/surface_fields/nimrod/generic_surface_spatial_plot_sequence_radar_rainfall.yaml @@ -0,0 +1,27 @@ +category: Surface Spatial Plot +title: "$RADAR_NAME $VARNAME $SUBAREA_NAME" +description: | + Extracts and plots radar rainfall data. + ATTENTION: the radar data can be 1 hour accumulation amounts in [mm]. Because showers + typically advect, on a 2D plot of accumulated rainfall, showers show + as streaks i.e. lower intensity rainfall over a greater area than a + plot of instantaneous rainfall rate. This results in accumulated rainfall + being damped in intensity compared to rainfall rates so caution is + required when comparing model rainfall rate with observed rainfall accumulation. + +steps: + - operator: read.read_cube + file_paths: $INPUT_PATHS + constraint: + operator: constraints.combine_constraints + varname_constraint: + operator: constraints.generate_var_constraint + varname: $VARNAME + subarea_type: $SUBAREA_TYPE + subarea_extent: $SUBAREA_EXTENT + + - operator: plot.spatial_pcolormesh_plot + sequence_coordinate: time + + - operator: write.write_cube_to_nc + overwrite: True diff --git a/src/CSET/recipes/surface_fields/nimrod/model_radar_common_domain_2d_model.yaml b/src/CSET/recipes/surface_fields/nimrod/model_radar_common_domain_2d_model.yaml new file mode 100644 index 0000000000..92ceb9eefd --- /dev/null +++ b/src/CSET/recipes/surface_fields/nimrod/model_radar_common_domain_2d_model.yaml @@ -0,0 +1,68 @@ +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. + ATTENTION: the radar data are 1 hour accumulations. Because showers + typically advect, on a 2D plot of accumulated rainfall, showers show + as streaks i.e. lower intensity rainfall over a greater area than a + plot of instantaneous rainfall rate. This results in accumulated rainfall + being damped in intensity compared to rainfall rates so caution is + required when comparing model rainfall rate with observed rainfall accumulation. + +steps: + - operator: read.read_cubes + 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: $BOUNDARY_MARGIN + 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/nimrod/model_radar_common_domain_2d_radar.yaml b/src/CSET/recipes/surface_fields/nimrod/model_radar_common_domain_2d_radar.yaml new file mode 100644 index 0000000000..f73b4067d9 --- /dev/null +++ b/src/CSET/recipes/surface_fields/nimrod/model_radar_common_domain_2d_radar.yaml @@ -0,0 +1,68 @@ +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. + ATTENTION: the radar data are 1 hour accumulations. Because showers + typically advect, on a 2D plot of accumulated rainfall, showers show + as streaks i.e. lower intensity rainfall over a greater area than a + plot of instantaneous rainfall rate. This results in accumulated rainfall + being damped in intensity compared to rainfall rates so caution is + required when comparing model rainfall rate with observed rainfall accumulation. + +steps: + - operator: read.read_cubes + 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: $BOUNDARY_MARGIN + 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/nimrod/radar_common_domain_2d_plot.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_2d_plot.yaml new file mode 100644 index 0000000000..4aab5a7768 --- /dev/null +++ b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_2d_plot.yaml @@ -0,0 +1,62 @@ +category: Surface Spatial Plot +title: "bmc2d" +#title: "bmc2c $RADAR_VARNAME\n common domain $MODEL_LABEL and radar" +description: | + Creates a 2-d surface plot of the model field using the + grid and masking appropriate to the Nimrod file $MASK_LABEL. An + 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, $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: $BOUNDARY_MARGIN + required_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/nimrod/radar_common_domain_histogram.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram.yaml new file mode 100644 index 0000000000..b6e4ed506e --- /dev/null +++ b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram.yaml @@ -0,0 +1,68 @@ +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 accumulation. + The model field is regridded onto the radar grid. Both regridded model + and radar grids are then masked so that only data within a domain common + to both are used. + If $SEQUENCE is set to "time" then time sequence plots are produced. + If $SEQUENCE is set to "realization" then a single case-study mean plot is produced. + ATTENTION: the radar data are 1 hour accumulations. Because showers + typically advect, on a 2D plot of accumulated rainfall, showers show + as streaks i.e. lower intensity rainfall over a greater area than a + plot of instantaneous rainfall rate. This results in accumulated rainfall + being damped in intensity compared to rainfall rates so caution is + required when comparing model rainfall rate with observed rainfall accumulation. + +steps: + - operator: read.read_cubes + 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: $BOUNDARY_MARGIN + required_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/nimrod/radar_common_domain_histogram_domain_mean.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_domain_mean.yaml new file mode 100644 index 0000000000..5102312e06 --- /dev/null +++ b/src/CSET/recipes/surface_fields/nimrod/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: $BOUNDARY_MARGIN + required_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/nimrod/radar_common_domain_histogram_sequence.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_common_domain_histogram_sequence.yaml new file mode 100644 index 0000000000..3c9116b9b9 --- /dev/null +++ b/src/CSET/recipes/surface_fields/nimrod/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: $BOUNDARY_MARGIN + required_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/nimrod/radar_dev3.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_dev3.yaml new file mode 100644 index 0000000000..4d62b52abf --- /dev/null +++ b/src/CSET/recipes/surface_fields/nimrod/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/nimrod/radar_masked_mean_time_series.yaml b/src/CSET/recipes/surface_fields/nimrod/radar_masked_mean_time_series.yaml new file mode 100644 index 0000000000..f723f53033 --- /dev/null +++ b/src/CSET/recipes/surface_fields/nimrod/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/nimrod/radar_mean_time_series.yaml similarity index 94% rename from src/CSET/recipes/surface_fields/radar_mean_time_series.yaml rename to src/CSET/recipes/surface_fields/nimrod/radar_mean_time_series.yaml index aaf417f4eb..293dc8c9e8 100644 --- a/src/CSET/recipes/surface_fields/radar_mean_time_series.yaml +++ b/src/CSET/recipes/surface_fields/nimrod/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) diff --git a/tests/operators/test_radar_filter.py b/tests/operators/test_radar_filter.py new file mode 100644 index 0000000000..e6d9c55bd2 --- /dev/null +++ b/tests/operators/test_radar_filter.py @@ -0,0 +1,490 @@ +# © Crown copyright, Met Office (2022-2026) and CSET contributors. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Test filter operators.""" + +import iris +import numpy as np +import pytest + +from CSET.operators import constraints, radar_filter + +constraint_single = constraints.combine_constraints( + constraints.generate_stash_constraint("m01s03i236"), + a=constraints.generate_cell_methods_constraint([]), +) + + +def make_cube(data, name, model_name=None): + """Create a cube with a model_name attribute when requested.""" + cube = iris.cube.Cube(np.asarray(data, dtype=float), long_name=name) + cube.rename(name) + if model_name is not None: + cube.attributes["model_name"] = model_name + return cube + + +class TestMaskList: + """Tests for the Nimrod mask selection helper.""" + + def test_returns_empty_list_when_no_nimrod_source_present(self): + """No Nimrod source should produce no mask names.""" + assert radar_filter.mask_list(["UM_model", "UKV_model"]) == [] + + def test_prefers_highest_resolution_nimrod_source(self): + """The preferred Nimrod weights should match the highest resolution source.""" + result = radar_filter.mask_list(["UM_model", "Nimrod1km", "Nimrod2km"]) + assert result == [ + "Nimrod2km_weights", + "Nimrod1km_weights", + "Nimrod2km_weights", + ] + + +class TestMaskByWeights: + """Tests for applying a radar weight mask to a field.""" + + def test_applies_mask_using_weight_threshold(self): + """Only cells with a valid Nimrod weight should remain in the field.""" + field = make_cube(np.arange(12, dtype=float).reshape(1, 3, 4), "rain_rate") + field.attributes["model_name"] = "UM_model" + + mask_values = np.array( + [ + [ + [0.0, 11.0, 12.0, 4.0], + [13.0, 0.0, 5.0, 14.0], + [10.0, 11.0, 0.0, 21.0], + ] + ], + dtype=float, + ) + weights = make_cube( + mask_values, "Nimrod2km_weights", model_name="Nimrod2km_weights" + ) + + result = radar_filter.mask_by_weights( + iris.cube.CubeList([field, weights]), + ["UM_model"], + ["Nimrod2km_weights"], + ) + + expected = np.full((1, 3, 4), np.nan, dtype=float) + valid = np.array( + [ + [ + [False, True, True, False], + [True, False, False, True], + [False, True, False, True], + ] + ], + dtype=bool, + ) + expected[valid] = field.data[valid] + + np.testing.assert_allclose(result.data, expected, equal_nan=True) + + +class TestRadarApplyMask: + """Tests for the boundary-aware radar mask application.""" + + def test_masks_boundary_points_and_applies_valid_mask(self): + """Boundary cells and mask zeros should both become NaN.""" + field = make_cube(np.arange(120, dtype=float).reshape(1, 10, 12), "rain_rate") + mask = make_cube(np.ones((1, 10, 12), dtype=float), "mask") + + result = radar_filter.radar_apply_mask(field, mask, boundary_margin=2) + + assert np.isnan(result.data[0, 0, 0]) + assert np.isnan(result.data[0, 9, 11]) + assert result.data[0, 2, 2] == 26.0 + assert result.data[0, 5, 5] == 65.0 + + +class TestRadarMask: + """Tests for masking radar and model data together.""" + + def test_returns_masked_radar_cube(self): + """Valid model points should propagate to the radar field output.""" + model_field = make_cube( + np.arange(120, dtype=float).reshape(1, 10, 12), "model_field" + ) + radar_field = make_cube( + np.arange(120, 240, dtype=float).reshape(1, 10, 12), "radar_field" + ) + mask = make_cube(np.ones((1, 10, 12), dtype=float), "nimrod_mask") + + result = radar_filter.radar_mask( + model_field, radar_field, mask, boundary_margin=2, required_outputs="radar" + ) + + assert isinstance(result, iris.cube.Cube) + assert np.isnan(result.data[0, 0, 0]) + assert np.isnan(result.data[0, 9, 11]) + assert result.data[0, 2, 2] == 146.0 + + def test_returns_both_model_and_radar_cubes_for_all_output(self): + """The all-output mode should return both model and radar cubes.""" + model_field = make_cube( + np.arange(120, dtype=float).reshape(1, 10, 12), "model_field" + ) + radar_field = make_cube( + np.arange(120, 240, dtype=float).reshape(1, 10, 12), "radar_field" + ) + mask = make_cube(np.ones((1, 10, 12), dtype=float), "nimrod_mask") + + result = radar_filter.radar_mask( + model_field, radar_field, mask, boundary_margin=2, required_outputs="all" + ) + + assert isinstance(result, iris.cube.CubeList) + assert len(result) == 2 + assert result[0].name() == "model_field" + assert result[1].name() == "radar_field" + + +def test_handles_multiple_field_and_mask_cubes(): + """Test the loop over both the mask and rainfall fields.""" + field1 = make_cube(np.ones((1, 3, 3), dtype=float), "field_a") + field2 = make_cube(np.full((1, 3, 3), 2.0, dtype=float), "field_b") + + mask1 = make_cube(np.ones((1, 3, 3), dtype=float), "mask_a") + mask2 = make_cube(np.ones((1, 3, 3), dtype=float), "mask_b") + + result = radar_filter.radar_apply_mask( + iris.cube.CubeList([field1, field2]), + iris.cube.CubeList([mask1, mask2]), + boundary_margin=0, + ) + + assert isinstance(result, iris.cube.CubeList) + assert len(result) == 2 + assert np.allclose(result[0].data, 1.0) + assert np.allclose(result[1].data, 2.0) + + +def test_returns_masked_model_cube_when_output_is_model(): + """Test the branch that executes if just model output is required.""" + model_field = make_cube( + np.arange(120, dtype=float).reshape(1, 10, 12), "model_field" + ) + radar_field = make_cube( + np.arange(120, 240, dtype=float).reshape(1, 10, 12), "radar_field" + ) + mask = make_cube(np.ones((1, 10, 12), dtype=float), "nimrod_mask") + + result = radar_filter.radar_mask( + model_field, radar_field, mask, boundary_margin=2, required_outputs="model" + ) + + assert isinstance(result, iris.cube.Cube) + assert result.name() == "model_field" + assert np.isnan(result.data[0, 0, 0]) + + +def test_radar_mask_loop_returns_model_and_radar_outputs(): + """Test the loop over all models.""" + model_field = iris.cube.CubeList( + [ + make_cube(np.ones((1, 3, 3), dtype=float), "model_1"), + make_cube(np.full((1, 3, 3), 2.0, dtype=float), "model_2"), + ] + ) + radar_field = make_cube(np.full((1, 3, 3), 9.0, dtype=float), "radar_field") + mask = make_cube(np.ones((1, 3, 3), dtype=float), "nimrod_mask") + + result = radar_filter.radar_mask_loop( + model_field, + radar_field, + mask, + boundary_margin=0, + required_outputs="radar", + ) + + assert isinstance(result, iris.cube.CubeList) + assert len(result) >= 2 + + +def test_boundary_margin_zero_keeps_all_inner_points(): + """Test zero width boundary on model field.""" + field = make_cube(np.arange(9, dtype=float).reshape(1, 3, 3), "rain_rate") + mask = make_cube(np.ones((1, 3, 3), dtype=float), "mask") + + result = radar_filter.radar_apply_mask(field, mask, boundary_margin=0) + + assert np.all(np.isfinite(result.data)) + + +# Session scope fixtures, so the test data only has to be loaded once. +@pytest.fixture(scope="session") +def cube_radar() -> iris.cube.Cube: + """Construct a radar rainfall accumulation cube.""" + shape = (1, 3, 3) + data = np.ones(shape, dtype=np.float32) + + lat_coord = iris.coords.DimCoord( + np.linspace(-90, 90, 3), standard_name="latitude", units="degrees" + ) + lon_coord = iris.coords.DimCoord( + np.linspace(-180, 180, 3), standard_name="longitude", units="degrees" + ) + + time_points = np.array([3], dtype=np.float64) + time_coord = iris.coords.DimCoord( + time_points, + standard_name="time", + units="hours since 2026-09-23 00:00:00", + ) + + fp_coord = iris.coords.DimCoord( + time_points, standard_name="forecast_period", units="hours" + ) + + test_cube = iris.cube.Cube( + data, + long_name="Hourly rain accumulation", + units="mm", + dim_coords_and_dims=[ + (time_coord, 0), + (lat_coord, 1), + (lon_coord, 2), + ], + ) + + test_cube.add_aux_coord(fp_coord, data_dims=0) + + frt_coord = iris.coords.AuxCoord( + [0.0], + standard_name="forecast_reference_time", + units="hours since 2026-09-23 00:00:00", + ) + test_cube.add_aux_coord(frt_coord) + + test_cube.attributes["model_name"] = "RadarA" + + new_data = np.array([[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]]) + test_cube.data = new_data + + return test_cube + + +@pytest.fixture(scope="session") +def cube_radar_wts(cube_radar) -> iris.cube.Cube: + """Construct a radar rainfall accumulation weights cube.""" + test_cube = cube_radar.copy() + new_data = np.array([[[13, 0, 13], [9, 12, 11], [13, 11, 10]]]) + test_cube.data = new_data + test_cube.attributes["model_name"] = "RadarA_wts" + test_cube.long_name = "Hourly wts accumulation" + test_cube.units = "1" + return test_cube + + +@pytest.fixture(scope="session") +def cube_radar_masked_by_wts(cube_radar) -> iris.cube.Cube: + """Construct a radar rainfall accumulation cube that has been masked by the weights.""" + test_cube = cube_radar.copy() + new_data = np.array([[[1.0, np.nan, 3.0], [np.nan, 5.0, 6.0], [7.0, 8.0, np.nan]]]) + test_cube.data = new_data + test_cube.attributes["mask"] = "mask_of_" + test_cube.long_name + return test_cube + + +@pytest.fixture(scope="session") +def cube_model(cube_radar) -> iris.cube.Cube: + """Construct a model rainfall rate cube.""" + test_cube = cube_radar.copy() + new_data = np.array([[[11.0, 12.0, 13.0], [14.0, 15.0, 16.0], [17.0, 18.0, 19.0]]]) + test_cube.data = new_data + test_cube.attributes["model_name"] = "ModelA" + test_cube.long_name = "surface_microphysical_rainfall_rate" + test_cube.units = "mm hr-1" + return test_cube + + +def test_print_cube_radar(cube_radar): + """Print the radar cube.""" + print("----> cube_radar ", cube_radar) + print() + print("----> cube_radar.data ", cube_radar.data) + print() + print( + '----> cube_radar.attributes["model_name"] ', + cube_radar.attributes["model_name"], + ) + + +def test_print_cube_radar_mask(cube_radar_wts): + """Print the radar wts cube.""" + print("----> cube_radar_wts ", cube_radar_wts) + print() + print("----> cube_radar_wts.data ", cube_radar_wts.data) + + +def test_print_cube_radar_masked_by_wts(cube_radar_masked_by_wts): + """Print the radar masked by wts cube.""" + print("----> cube_radar_masked_by_wts ", cube_radar_masked_by_wts) + print() + print("----> cube_radar_masked_by_wts.data ", cube_radar_masked_by_wts.data) + + +def test_mask_list(): + """Test the mask_list function.""" + list_weights = radar_filter.mask_list(["UM_model", "Nimrod1km", "Nimrod2km"]) + expected_list = ["Nimrod2km_weights", "Nimrod1km_weights", "Nimrod2km_weights"] + assert list_weights == expected_list + + +def test_mask_by_weights(cube_radar, cube_radar_wts, cube_radar_masked_by_wts): + """Test the mask_by_weights function.""" + obs_name = cube_radar.attributes["model_name"] + wts_name = cube_radar_wts.attributes["model_name"] + cube_list = iris.cube.CubeList([]) + cube_list.append(cube_radar) + cube_list.append(cube_radar_wts) + masked = radar_filter.mask_by_weights(cube_list, [obs_name], [wts_name]) + print() + print("---->masked ", masked) + print() + print("----> masked.data ", masked.data) + assert np.array_equal(masked.data, cube_radar_masked_by_wts.data, equal_nan=True) + + +def test_radar_apply_mask(cube_model, cube_radar_wts): + """Test the mask_by_weights function.""" + masked = radar_filter.radar_apply_mask( + cube_model, cube_radar_wts, boundary_margin=0 + ) + print() + print("---->masked from test_radar_apply_mask", masked) + print() + print("----> masked.data ", masked.data) + + +def test_match_varname_and_units(cube_model, cube_radar): + """Test the match_varname_and_units function.""" + cube1 = cube_model.copy() + cube2 = cube_model.copy() + cube3 = cube_radar.copy() + + # Construct the CubeList to test. + cube_list = iris.cube.CubeList([]) + cube_list.append(cube1) + cube_list.append(cube2) + cube_list.append(cube3) + + # Push the test CubeList through the function match_varname_and_units. + cubes_matched = radar_filter.match_varname_and_units(cube_list) + + assert cubes_matched[0].var_name == cubes_matched[1].var_name + assert cubes_matched[0].var_name == cubes_matched[2].var_name + + assert cubes_matched[0].long_name == cubes_matched[1].long_name + assert cubes_matched[0].long_name == cubes_matched[2].long_name + + assert cubes_matched[0].units == cubes_matched[1].units + assert cubes_matched[0].units == cubes_matched[2].units + + +def test_match_varname_and_units_single_cube(cube_model): + """Test the match_varname_and_units function.""" + cube1 = cube_model.copy() + + # Push the test Cube through the function match_varname_and_units. + cubes_matched = radar_filter.match_varname_and_units(cube1) + + assert cubes_matched.var_name == cube1.var_name + assert cubes_matched.long_name == cube1.long_name + assert cubes_matched.units == cube1.units + + +def test_match_varname_and_units_single_cube_in_cubelist(cube_model): + """Test the match_varname_and_units function.""" + cube1 = cube_model.copy() + + # Construct the CubeList to test. + cube_list = iris.cube.CubeList([]) + cube_list.append(cube1) + + # Push the test CubeList through the function match_varname_and_units. + cubes_matched = radar_filter.match_varname_and_units(cube_list) + + assert cubes_matched[0].var_name == cube_model.var_name + assert cubes_matched[0].long_name == cube_model.long_name + assert cubes_matched[0].units == cube_model.units + + +def test_mask_boundary_zero_margin_leaves_data_unchanged(): + """Test the mask_boundary function with boundary_margin set to 0.""" + data = np.zeros((2, 5, 6)) + radar_filter._mask_boundary(data, 0) + assert not np.isnan(data).any() + + +def test_mask_boundary_blanks_edges_only(): + """Test the mask_boundary function with boundary_margin set to 1.""" + data = np.zeros((5, 6)) + radar_filter._mask_boundary(data, 1) + expected_nan = np.ones((5, 6), dtype=bool) + expected_nan[1:-1, 1:-1] = False + np.testing.assert_array_equal(np.isnan(data), expected_nan) + + +def test_mask_boundary_is_symmetric(): + """Test the mask_boundary function with boundary_margin set to 3.""" + data = np.zeros((11, 11)) + radar_filter._mask_boundary(data, 3) + nan_rows = np.isnan(data[:, 5]) + np.testing.assert_array_equal(np.flatnonzero(nan_rows), [0, 1, 2, 8, 9, 10]) + + +def test_mask_boundary_applies_to_all_leading_dims(): + """Test the mask_boundary function with boundary_margin set to 1 and with time + realisation coords.""" + data = np.zeros((3, 2, 5, 5)) # e.g. (time, realization, y, x) + radar_filter._mask_boundary(data, 1) + for t in range(3): + for r in range(2): + assert np.isnan(data[t, r, 0, :]).all() + assert not np.isnan(data[t, r, 2, 2]) + + +def test_mask_boundary_negative_margin_raises(): + """Test bad boundary margin raises.""" + with pytest.raises(ValueError): + radar_filter._mask_boundary(np.zeros((4, 4)), -1) + + +# def test_radar_apply_mask(): +# """Test the radar_apply_mask function.""" +# obs = Cube([0,[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]]) +# #wts = Cube([[13, 0, 13], [9, 12, 11], [13, 11, 10]]) +# msk = Cube([0,[[1, np.nan, 1], [np.nan, 1, 1], [1, 1, 1]]]) +# +# var_name_obs = "hourly_rain_accumulation" +# var_name_msk = "hourly_wts_accumulation" +# +# obs_name = "Radar" +# msk_name = "Radar_msk" +# obs.attributes["model_name"] = obs_name +# msk.attributes["model_name"] = msk_name +# +# obs.rename(var_name_obs) +# obs.long_name = var_name_obs +# obs.var_name = var_name_obs +# +# msk.rename(var_name_msk) +# msk.long_name = var_name_msk +# msk.var_name = var_name_msk +# +# masked = radar_filter.radar_apply_mask(obs, msk, boundary_margin = 1) diff --git a/tests/test_fetch_nimrod.py b/tests/test_fetch_nimrod.py index 7c2f2cb52b..d0337a5e8b 100644 --- a/tests/test_fetch_nimrod.py +++ b/tests/test_fetch_nimrod.py @@ -84,11 +84,12 @@ def test_get_needed_environment_variables_nimrod(monkeypatch): def test_retrieve_nimrod_saves_observation_and_weights(tmp_path, monkeypatch): """Test retrieve_nimrod saves both observation and weight files when available.""" # Create a temporary nimrod dictionary file. + basedir = tmp_path / "nimrod" nimrod_dict = { "Nimrod_comp_xkm": { "obs_id": "nimrod_xkm", "wei_id": "nimrod_xkm_weights", - "basedir": "/data/nimrod", + "basedir": str(basedir), "obs_dir": "obs", "obs_fname": "_obs.nc", "weights_dir": "weights", @@ -96,17 +97,22 @@ def test_retrieve_nimrod_saves_observation_and_weights(tmp_path, monkeypatch): "freq": "PT01H", } } + + for subdir, fname in [("obs", "_obs.nc"), ("weights", "_weights.nc")]: + f = basedir / subdir / "2026" / f"202604071300{fname}" + f.parent.mkdir(parents=True) + f.touch() + config_path = tmp_path / "restricted_nimrod_met_office.json" config_path.write_text(json.dumps(nimrod_dict)) - # Set the required environment variables. monkeypatch.setenv("NIMROD_COMP_XKM", "True") monkeypatch.setenv("NIMROD_COMP_1KM", "False") monkeypatch.setenv("NIMROD_COMP_2KM", "False") monkeypatch.setenv("NIMROD_COMP_5MIN", "False") monkeypatch.setenv("NIMROD_WEIGHTS", "True") - monkeypatch.setenv("CYLC_TASK_CYCLE_POINT", "2026-05-07T12:00:00") - monkeypatch.setenv("ANALYSIS_LENGTH", "PT00H") + monkeypatch.setenv("CYLC_TASK_CYCLE_POINT", "2026-04-07T12:00:00") + monkeypatch.setenv("ANALYSIS_LENGTH", "PT01H") rose_datac = tmp_path / "rose_datac" monkeypatch.setenv("ROSE_DATAC", str(rose_datac)) @@ -144,13 +150,13 @@ def save_side_effect(cube, filename): fetch_nimrod.apply_radar_weights.assert_called_once_with(obs_cube, weights_cube) expected_obs_file = str( - Path(rose_datac) / "data" / "nimrod_xkm" / "202605071200_nimrod_xkm.nc" + Path(rose_datac) / "data" / "nimrod_xkm" / "202604071300_nimrod_xkm.nc" ) expected_weight_file = str( Path(rose_datac) / "data" / "nimrod_xkm_weights" - / "202605071200_nimrod_xkm_weights.nc" + / "202604071300_nimrod_xkm_weights.nc" ) saved_paths = [path for _, path in save_calls]