Support deterministic forecast source when clustering

improver/clustering/realization_clustering.py

@@ -535,6 +535,27 @@ def _convert_hours_to_seconds(hours: list[int]) -> list[int]:
         """
         return [h * 3600 for h in hours]
 
+    @staticmethod
+    def _ensure_realization_coord(cube: Cube) -> Cube:
+        """Add a scalar realization dimension coordinate if absent.
+
+        Deterministic input cubes carry no realization coordinate. This method
+        adds a realization DimCoord with value 0 as the leading dimension so
+        that downstream matching code can treat deterministic and ensemble inputs
+        uniformly.
+
+        Args:
+            cube: The input cube, which may or may not have a realization coordinate.
+
+        Returns:
+            The cube with a realization dimension coordinate as the leading axis.
+            If the cube already has a realization coordinate, it is returned unchanged.
+        """
+        if not cube.coords("realization"):
+            cube.add_aux_coord(DimCoord(0, standard_name="realization", units="1"))
+            cube = new_axis(cube, "realization")
+        return cube
+
     def cluster_primary_input(
         self, primary_cube: Cube, target_grid_cube: Cube | None
     ) -> tuple[Cube, Cube]:
@@ -666,6 +687,9 @@ def _categorise_secondary_inputs(
             UserWarning: If a secondary input has an inconsistent realization count
                 compared to its earliest valid forecast period; all forecast periods
                 from the first mismatch onwards are dropped for that input.
+            UserWarning: If a hierarchy secondary input has no cubes matching the
+                model_id_attr value in the specified forecast period range; that
+                input is ignored.
         """
         full_realization_inputs = []
         partial_realization_inputs = []
@@ -683,6 +707,12 @@ def _categorise_secondary_inputs(
             fp_constr = iris.Constraint(forecast_period=fp_seconds_range)
             model_cubes = cubes.extract(model_id_constr & fp_constr)
             if not model_cubes:
+                warnings.warn(
+                    f"Secondary input '{candidate_name}' has no cubes matching "
+                    f"{self.model_id_attr}='{candidate_name}' in the forecast period "
+                    f"range {fp_range}. This input will be ignored.",
+                    UserWarning,
+                )
                 continue  # No cubes found in this range for this model
 
             # Get forecast period and cube pairs for this model.
@@ -721,12 +751,21 @@ def _categorise_secondary_inputs(
             # forecast period and truncate this secondary input if the count changes
             # at later lead times. This avoids a source pulsing in/out and keeps
             # all periods for this source mergeable for consistent multi-period
-            # matching.
+            # matching. If no realization coordinate exists, treat it as a single
+            # realization (deterministic).
             first_fp, first_cube = valid_fp_cube_pairs[0]
-            n_realizations = len(first_cube.coord("realization").points)
+            n_realizations = (
+                len(first_cube.coord("realization").points)
+                if first_cube.coords("realization")
+                else 1
+            )
             forecast_periods_in_range = []
             for idx, (fp, cube) in enumerate(valid_fp_cube_pairs):
-                n_realizations_at_fp = len(cube.coord("realization").points)
+                n_realizations_at_fp = (
+                    len(cube.coord("realization").points)
+                    if cube.coords("realization")
+                    else 1
+                )
                 if n_realizations_at_fp != n_realizations:
                     dropped_fps = [
                         future_fp for future_fp, _ in valid_fp_cube_pairs[idx:]
@@ -1083,6 +1122,7 @@ def _extract_merge_and_match(
         )
         fp_constr = iris.Constraint(forecast_period=fps)
         candidate_cube = MergeCubes()(cubes.extract(model_id_constr & fp_constr))
+        candidate_cube = self._ensure_realization_coord(candidate_cube)
         enforce_coordinate_ordering(candidate_cube, ["realization"])
         if ensure_fp_dim:
             candidate_cube = self._ensure_forecast_period_is_dimension(candidate_cube)
@@ -1261,6 +1301,7 @@ def _process_partial_realization_inputs(
             for fp in forecast_periods:
                 fp_constr = iris.Constraint(forecast_period=fp)
                 candidate_cube = cubes.extract_cube(model_id_constr & fp_constr)
+                candidate_cube = self._ensure_realization_coord(candidate_cube)
 
                 # Index the candidate cube using the realization indices determined
                 # from the combined match across all forecast periods.
@@ -1347,7 +1388,7 @@ def process(self, cubes: CubeList) -> Cube:
             ValueError: If no primary cube is found with the specified
                 model_id_attr.
 
-        Warnings:
+        Warns:
             UserWarning: If primary cubes have different realization numbering schemes
                 when renumber_primary_realizations=False, which may cause merge
                 failures.
@@ -1356,6 +1397,8 @@ def process(self, cubes: CubeList) -> Cube:
                 be returned.
             UserWarning: If secondary inputs have forecast periods not present in the
                 primary input, which will be ignored.
+            UserWarning: If input cubes have model_id_attr values not referenced in
+                the hierarchy; those cubes will be ignored.
         """
         if self.cycletime is not None:
             for cube in cubes:
@@ -1396,6 +1439,23 @@ def process(self, cubes: CubeList) -> Cube:
                 f"{self.hierarchy['primary_input']}"
             )
 
+        # Warn about cubes with model_id_attr values not referenced in the hierarchy.
+        hierarchy_names = {self.hierarchy["primary_input"]} | set(
+            self.hierarchy["secondary_inputs"].keys()
+        )
+        input_model_ids = {
+            cube.attributes[self.model_id_attr]
+            for cube in cubes
+            if self.model_id_attr in cube.attributes
+        }
+        unreferenced = input_model_ids - hierarchy_names
+        if unreferenced:
+            warnings.warn(
+                f"Input cubes have {self.model_id_attr} values not referenced in the "
+                f"hierarchy: {sorted(unreferenced)}. These cubes will be ignored.",
+                UserWarning,
+            )
+
         target_grid_cube = None
         if self.regrid_for_clustering:
             try:
@@ -1457,7 +1517,6 @@ def process(self, cubes: CubeList) -> Cube:
             cluster_sources[cluster_idx][primary_name] = list(
                 clustered_primary_cube.coord("forecast_period").points
             )
-
         # Create a mapping to track which realizations from secondary inputs correspond
         # to which clusters.
         secondary_input_realizations_to_clusters = {}
@@ -1491,7 +1550,6 @@ def process(self, cubes: CubeList) -> Cube:
         result_cube = MergeCubes()(
             CubeList([iris.util.squeeze(c) for c in matched_cubes])
         )
-
         # Use json.dumps to store dictionary as attribute.
         result_cube.attributes["primary_input_realizations_to_clusters"] = json.dumps(
             primary_input_realizations_to_clusters

improver_tests/clustering/test_realization_clustering.py

@@ -3052,6 +3052,220 @@ def test_clusterandmatch_secondary_input_missing_primary_forecast_period(
             )
 
 
+def test_clusterandmatch_secondary_no_matching_cubes_warns():
+    """Test that a warning is issued when a hierarchy secondary input has no cubes
+    matching the model_id_attr value in the specified forecast period range.
+
+    This verifies that when the hierarchy references a secondary model name that does
+    not match any cube's model_id attribute, the plugin warns and skips that input
+    rather than silently ignoring it.
+    """
+    pytest.importorskip("kmedoids")
+    pytest.importorskip("esmf_regrid")
+
+    cubes = CubeList()
+    spatial_shape = (3, 3)
+
+    cubes.extend(
+        _create_4d_realization_cube(
+            n_realizations=3,
+            forecast_periods=[0, 6],
+            y_dim=spatial_shape[0],
+            x_dim=spatial_shape[1],
+            base_value=100.0,
+            model_id="primary_model",
+            merge=False,
+        )
+    )
+    cubes.append(_create_target_grid_cube(spatial_shape=spatial_shape))
+
+    hierarchy = {
+        "primary_input": "primary_model",
+        "secondary_inputs": {"nonexistent_model": [0, 6]},
+    }
+
+    plugin = RealizationClusterAndMatch(
+        hierarchy=hierarchy,
+        model_id_attr="model_id",
+        clustering_method="KMedoids",
+        target_grid_name="target_grid",
+        n_clusters=2,
+        random_state=42,
+    )
+
+    with pytest.warns(
+        UserWarning,
+        match=(
+            r"Secondary input 'nonexistent_model' has no cubes matching "
+            r"model_id='nonexistent_model' in the forecast period range \[0, 6\]\."
+        ),
+    ):
+        plugin.process(cubes)
+
+
+def test_clusterandmatch_unreferenced_model_id_warns():
+    """Test that a warning is issued when input cubes contain model_id_attr values
+    not referenced in the hierarchy.
+
+    This verifies that when extra cubes are supplied with a model_id attribute value
+    that does not appear in the hierarchy (neither as primary_input nor as a secondary
+    input key), the plugin warns that those cubes will be ignored.
+    """
+    pytest.importorskip("kmedoids")
+    pytest.importorskip("esmf_regrid")
+
+    cubes = CubeList()
+    spatial_shape = (3, 3)
+
+    cubes.extend(
+        _create_4d_realization_cube(
+            n_realizations=3,
+            forecast_periods=[0, 6],
+            y_dim=spatial_shape[0],
+            x_dim=spatial_shape[1],
+            base_value=100.0,
+            model_id="primary_model",
+            merge=False,
+        )
+    )
+    # Extra cube whose model_id is not in the hierarchy.
+    cubes.extend(
+        _create_4d_realization_cube(
+            n_realizations=3,
+            forecast_periods=[0, 6],
+            y_dim=spatial_shape[0],
+            x_dim=spatial_shape[1],
+            base_value=200.0,
+            model_id="unlisted_model",
+            merge=False,
+        )
+    )
+    cubes.append(_create_target_grid_cube(spatial_shape=spatial_shape))
+
+    hierarchy = {
+        "primary_input": "primary_model",
+        "secondary_inputs": {},
+    }
+
+    plugin = RealizationClusterAndMatch(
+        hierarchy=hierarchy,
+        model_id_attr="model_id",
+        clustering_method="KMedoids",
+        target_grid_name="target_grid",
+        n_clusters=2,
+        random_state=42,
+    )
+
+    with pytest.warns(
+        UserWarning,
+        match=(
+            r"Input cubes have model_id values not referenced in the hierarchy: "
+            r"\['unlisted_model'\]\. These cubes will be ignored\."
+        ),
+    ):
+        result = plugin.process(cubes)
+
+    assert all(
+        "unlisted_model" not in str(value)
+        for value in result.attributes.values()
+    )
+
+
+def test_clusterandmatch_deterministic_secondary_input():
+    """Test that a deterministic (no-realization) secondary input is supported.
+
+    The primary input provides the baseline clustered forecast for all forecast periods.
+    When a secondary input cube has no realization coordinate it should be treated
+    as having a single realization and follow the partial-realization path, replacing
+    the best-matching cluster at the relevant forecast periods with its data.
+    """
+    pytest.importorskip("kmedoids")
+    pytest.importorskip("esmf_regrid")
+
+    cubes = CubeList()
+    spatial_shape = (3, 3)
+
+    # Primary ensemble input — 4 realizations, base value 100.
+    cubes.extend(
+        _create_4d_realization_cube(
+            n_realizations=4,
+            forecast_periods=[0, 6],
+            y_dim=spatial_shape[0],
+            x_dim=spatial_shape[1],
+            base_value=100.0,
+            model_id="primary_model",
+            merge=False,
+        )
+    )
+
+    # Deterministic secondary input — 2D cubes with no realization coordinate,
+    # base value 500 (far from primary so it will match the nearest cluster).
+    for fp_hours in [0, 6]:
+        det_data = np.full(spatial_shape, 500.0 + fp_hours, dtype=np.float32)
+        det_cube = set_up_variable_cube(
+            det_data,
+            name="air_temperature",
+            units="K",
+            spatial_grid="equalarea",
+        )
+        det_cube.attributes["model_id"] = "det_model"
+        # Set time coordinates to match the primary cubes.
+        det_cube.coord("forecast_period").points = [fp_hours * 3600]
+        det_cube.coord("time").points = [
+            det_cube.coord("forecast_reference_time").points[0] + fp_hours * 3600
+        ]
+        assert not det_cube.coords("realization"), (
+            "Deterministic cube should have no realization coordinate"
+        )
+        cubes.append(det_cube)
+
+    cubes.append(_create_target_grid_cube(spatial_shape=spatial_shape))
+
+    hierarchy = {
+        "primary_input": "primary_model",
+        "secondary_inputs": {"det_model": [0, 6]},
+    }
+
+    plugin = RealizationClusterAndMatch(
+        hierarchy=hierarchy,
+        model_id_attr="model_id",
+        clustering_method="KMedoids",
+        target_grid_name="target_grid",
+        n_clusters=2,
+        random_state=42,
+    )
+
+    result = plugin.process(cubes)
+
+    # Result should have realization as a dimension coordinate.
+    assert result.coords("realization", dim_coords=True)
+    assert result.coord("realization").points.size == 2
+
+    # Both forecast periods should be present.
+    np.testing.assert_array_equal(
+        result.coord("forecast_period").points, [0, 6 * 3600]
+    )
+
+    # The deterministic value (500 / 506) is far from the primary (100 / 106), so
+    # exactly one cluster should carry the deterministic data at each lead time.
+    for fp_hours in [0, 6]:
+        fp_data = result.extract(
+            iris.Constraint(forecast_period=fp_hours * 3600)
+        ).data
+        expected_det = 500.0 + fp_hours
+        expected_primary = 100.0 + fp_hours
+        det_clusters = np.isclose(fp_data, expected_det, atol=5.0).any(axis=(-1, -2))
+        primary_clusters = np.isclose(fp_data, expected_primary, atol=5.0).any(
+            axis=(-1, -2)
+        )
+        assert det_clusters.sum() == 1, (
+            f"fp={fp_hours}h: exactly one cluster should carry deterministic data"
+        )
+        assert primary_clusters.sum() == 1, (
+            f"fp={fp_hours}h: exactly one cluster should carry primary data"
+        )
+
+
 def test_select_realizations_for_kmedoid_clusters_too_many_clusters():
     """Test that ValueError is raised if number of clusters > number of realizations."""
     # Create a cube with 2 realizations