2113 nimrod regrid

src/CSET/loaders/radar.py

@@ -85,6 +85,9 @@ def get_radar_sources(conf) -> list[dict]:
 
 def load(conf: Config):
     """Yield recipes from the given workflow configuration."""
+    # Load a list of model detail dictionaries.
+    # models = get_models(conf.asdict())
+
     # Load the required radar observation sources.
     radar_sources = get_radar_sources(conf)
 
@@ -95,6 +98,98 @@ def load(conf: Config):
         if radar["varname"] == "Hourly rain accumulation"
     ]
 
+    # Form the list of accumulated hourly weights for Nimrod radar sources.
+    wts_radars = [
+        radar
+        for radar in radar_sources
+        if radar["varname"] == "Hourly wts accumulation"
+    ]
+
+    #    # Radar masking based on sea mask.
+    #    if conf.SPATIAL_SURFACE_FIELD:
+    #        for field in conf.SURFACE_FIELDS:
+    #            yield RawRecipe(
+    #                recipe="sea_mask_for_surface_domain_mean_time_series.yaml",
+    #                variables={
+    #                    "VARNAME": field,
+    #                    "MODEL_NAME": [model["name"] for model in models],
+    #                    "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None,
+    #                    "SUBAREA_EXTENT": conf.SUBAREA_EXTENT
+    #                    if conf.SELECT_SUBAREA
+    #                    else None,
+    #                    "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "",
+    #                },
+    #                model_ids=[model["id"] for model in models],
+    #                aggregation=False,
+    #            )
+
+    # Radar masking of radar obs based on sea mask.
+    if conf.SPATIAL_SURFACE_FIELD:
+        field = "Hourly rain accumulation"
+        yield RawRecipe(
+            recipe="radar_mask_model.yaml",
+            variables={
+                "VARNAME": field,
+                "MODEL_LABEL": "Nimrod2km",
+                "MASK_LABEL": "Nimrod2km",
+                "METHOD": "SEQ",
+                "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None,
+                "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None,
+                "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "",
+            },
+            model_ids=["Nimrod2km", "Nimrod2km_weights"],
+            aggregation=False,
+        )
+
+    # Radar masking of model rainfall based on sea mask.
+    if conf.SPATIAL_SURFACE_FIELD:
+        field = "surface_microphysical_rainfall_rate"
+        yield RawRecipe(
+            recipe="radar_mask_model.yaml",
+            variables={
+                "VARNAME": field,
+                "MODEL_LABEL": "ModelA",
+                "MASK_LABEL": "Nimrod2km",
+                "METHOD": "SEQ",
+                "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None,
+                "SUBAREA_EXTENT": conf.SUBAREA_EXTENT if conf.SELECT_SUBAREA else None,
+                "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "",
+            },
+            model_ids=["1", "Nimrod2km_weights"],
+            aggregation=False,
+        )
+
+    # Surface (2D) fields for model rainfall masked by Nimrod radar.
+    #
+    # The different sources of Nimrod rainfall accumulation have
+    # different spatial grids. So each source requires its own
+    # recipe to prevent incompatible cubes being created.
+    #    if conf.SPATIAL_SURFACE_FIELD:
+    #        radar_source = ["Nimrod_2km"]
+    #        for radar in radar_source:
+    #            model_labels = [model["id"] for model in models]
+    #            radar_label = ["Nimrod2km_weights"]
+    #            combined_ids = [model_labels[0]] + radar_label
+    #            print("Combined ids is: ", combined_ids)
+    #            yield RawRecipe(
+    #                recipe="radar_plot_sequence_rainfall.yaml",
+    ##                model_ids=radar["id"],  # -> Becomes $INPUT_PATHS
+    #                model_ids=combined_ids,
+    #                variables={
+    ##                    "VARNAME": radar["varname"],
+    ##                    "RADAR_NAME": radar["name"],
+    #                    "RADAR_NAME": "Nimrod_2km_weights",
+    ##                    "MODEL_NAME": [model["name"] for model in models],
+    #                    "MODEL_NAME": "ModelA",
+    #                    "METHOD": "SEQ",
+    #                    "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None,
+    #                    "SUBAREA_EXTENT": conf.SUBAREA_EXTENT
+    #                    if conf.SELECT_SUBAREA
+    #                    else None,
+    #                },
+    #                aggregation=False,
+    #            )
+
     # Surface (2D) fields for Nimrod radar rainfall.
     #
     # The different sources of Nimrod rainfall accumulation have
@@ -134,6 +229,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,
             },
@@ -158,6 +287,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:
         yield RawRecipe(

src/CSET/operators/__init__.py

@@ -41,6 +41,7 @@
     plot,
     precipitation,
     pressure,
+    radar_filter,
     read,
     regrid,
     temperature,
@@ -68,6 +69,7 @@
     "plot",
     "precipitation",
     "pressure",
+    "radar_filter",
     "read",
     "regrid",
     "temperature",

src/CSET/operators/plot.py

@@ -800,7 +800,7 @@ def _plot_and_save_spatial_plot(
     # Add watermark with min/max/mean. Currently not user togglable.
     # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
     axes.annotate(
-        f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}",
+        f"Min: {np.nanmin(cube.data):.3g} Max: {np.nanmax(cube.data):.3g} Mean: {np.nanmean(cube.data):.3g}",
         xy=(0.025, yinfopad),
         xycoords="axes fraction",
         xytext=(-5, 5),
@@ -2270,7 +2270,6 @@ def _custom_colourmap_precipitation(cube: iris.cube.Cube, cmap, levels, norm):
     )
 
     if is_rainfall_var:
-        print("varnames_lower ", varnames_lower)
         levels = [0, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256]
         colors = [
             "w",
@@ -2293,6 +2292,9 @@ def _custom_colourmap_precipitation(cube: iris.cube.Cube, cmap, levels, norm):
         norm = mcolors.BoundaryNorm(levels, cmap.N)
         logging.info("Using custom rainfall colourmap.")
 
+        # Set any Nan values to be plotted a light grey.
+        cmap.set_bad("#dcdcdc")
+
     return cmap, levels, norm
 
 
@@ -2362,6 +2364,8 @@ def _custom_colourmap_nimrod_weights(cube: iris.cube.Cube, cmap, levels, norm):
         # levels = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]
         norm = mcolors.BoundaryNorm(levels, cmap.N)
         colours = [
+            "#dcdcdc",
+            #            "darkgray",
             "#d10000",
             "purple",
             "#8f00d6",
@@ -2374,7 +2378,6 @@ def _custom_colourmap_nimrod_weights(cube: iris.cube.Cube, cmap, levels, norm):
             "#37a648",
             "#8edc64",
             "#c5ffc5",
-            "#dcdcdc",
             "#ffffff",
         ]
         # Create a custom colormap