Latest edits.

Author: gavinevans

improver/calibration/load_and_apply_quantile_regression_random_forest.py

@@ -97,12 +97,6 @@ def _get_inputs(
             except ValueError:
                 qrf_model = joblib.load(file_path)
 
-        if not qrf_model:
-            msg = (
-                "No QRF model found in the provided file paths. "
-                "A trained QRF model must be provided."
-            )
-            raise ValueError(msg)
         if not cube_inputs:
             msg = (
                 "No features found in the provided file paths. "
@@ -113,18 +107,34 @@ def _get_inputs(
         # Extract all additional cubes which are associated with a feature in the
         # feature_config.
         forecast_constraint = iris.Constraint(name=self.target_cube_name)
-        forecast_cube = cube_inputs.extract_cube(forecast_constraint)
+        forecast_cube = cube_inputs.extract(forecast_constraint)
+
+        if forecast_cube:
+            (forecast_cube,) = forecast_cube
+        else:
+            msg = (
+                "No target forecast provided. An input file representing the target "
+                "must be provided, even if the target will not be used as a feature. "
+                f"The target is '{self.target_cube_name}'."
+            )
+            raise ValueError(msg)
+
+        if not qrf_model:
+            return None, forecast_cube, None
+
+        if len(cube_inputs) != len(self.feature_config.keys()):
+            msg = (
+                "The number of cubes loaded does not match the number of features "
+                "expected. The number of cubes loaded was: "
+                f"{len(cube_inputs)}. The number of features expected was: "
+                f"{len(self.feature_config.keys())}."
+            )
+            raise ValueError(msg)
 
         # If target diagnostic not a feature in the training then remove.
         if self.target_cube_name not in self.feature_config.keys():
             cube_inputs.remove(forecast_cube)
 
-        if len(cube_inputs) + 1 != len(file_paths):
-            raise ValueError(
-                "Unable to identify the correct number of inputs. "
-                f"The number of file paths provided was: {len(file_paths)}. "
-                f"The number of items loaded into cubes was: {len(cube_inputs)}."
-            )
         return cube_inputs, forecast_cube, qrf_model
 
     @staticmethod
@@ -232,6 +242,8 @@ def process(
                 The calibrated forecast cube.
         """
         cube_inputs, forecast_cube, qrf_model = self._get_inputs(file_paths)
+        if not qrf_model:
+            return forecast_cube
         if forecast_cube.coords("percentile"):
             percentiles = self._compute_percentiles(forecast_cube, "percentile")
             cube_inputs = self._percentiles_to_realizations(cube_inputs)

improver/calibration/load_and_train_quantile_regression_random_forest.py

@@ -10,6 +10,7 @@
 import iris
 import numpy as np
 import pandas as pd
+import pyarrow as pa
 import pyarrow.parquet as pq
 
 from improver import PostProcessingPlugin
@@ -67,7 +68,7 @@ def _split_cubes_and_parquet_files(self, file_paths):
             try:
                 cube = load_cube(str(file_path))
                 cube_inputs.append(cube)
-            except (OSError, IsADirectoryError):
+            except IsADirectoryError:
                 # For loop here because the read_schema must read a .parquet file rather than a directory.
                 for file in Path(file_path).glob("**/*.parquet"):
                     try:
@@ -78,6 +79,19 @@ def _split_cubes_and_parquet_files(self, file_paths):
                     if forecast_table_path and truth_table_path:
                         break
 
+        if len(self.feature_config.keys()) not in [
+            len(cube_inputs),
+            len(cube_inputs) + 1,
+        ]:
+            msg = (
+                "The number of cubes loaded does not match the number of features "
+                "expected. These can mismatch if the some features are coming from the "
+                "historic forecast. The number of cubes loaded was: "
+                f"{len(cube_inputs)}. The number of features expected was: "
+                f"{len(self.feature_config.keys())}."
+            )
+            raise ValueError(msg)
+
         return forecast_table_path, truth_table_path, cube_inputs
 
     def _read_parquet_files(
@@ -108,12 +122,31 @@ def _read_parquet_files(
                 ("experiment", "==", self.experiment),
             ]
         ]
+        for file in Path(forecast_table_path).glob("**/*.parquet"):
+            if pq.read_schema(file).get_all_field_indices("percentile"):
+                altered_schema = FORECAST_SCHEMA
+            elif pq.read_schema(file).get_all_field_indices("realization"):
+                altered_schema = FORECAST_SCHEMA.remove(
+                    FORECAST_SCHEMA.get_field_index("percentile")
+                )
+                altered_schema = altered_schema.append(
+                    pa.field("realization", pa.int64())
+                )
+            else:
+                msg = (
+                    "The forecast parquet file is expected to contain either a "
+                    "'percentile' or 'realization' field. Neither was found."
+                )
+                raise ValueError(msg)
+            break
+
         forecast_df = pd.read_parquet(
             forecast_table_path,
             filters=filters,
-            schema=FORECAST_SCHEMA,
+            schema=altered_schema,
             engine="pyarrow",
         )
+
         # Convert df columns from ms to pandas timestamp object to work with existing code
         for column in ["time", "forecast_reference_time", "blend_time"]:
             forecast_df[column] = pd.to_datetime(
@@ -294,12 +327,6 @@ def process(
             forecast_df, truth_df, forecast_periods
         )
 
-        # Check the number of features provided as separate files, plus 2
-        # (the forecasts and truths for training) are equal to the total number of
-        # file paths provided.
-        if len(cube_inputs) + 2 != len(file_paths):
-            raise ValueError("Unable to identify the correct number of inputs")
-
         # If target_forecast is also a dynamic feature in the feature config then
         # add it to cube_inputs
         for feature_name in self.feature_config.keys():

improver/calibration/quantile_regression_random_forest.py

@@ -15,6 +15,7 @@
 
 from improver import BasePlugin, PostProcessingPlugin
 from improver.constants import DAYS_IN_YEAR, HOURS_IN_DAY
+from improver.utilities.cube_manipulation import enforce_coordinate_ordering
 
 
 def _remove_item_from_list(alist: list, items: list):
@@ -364,8 +365,13 @@ def _organise_truth_data(forecast_cube: Cube, truth_cube: Cube) -> list[np.ndarr
                 )
         else:
             # Forecast reference time and forecast period are different dimensions.
-            for frt in list(frt_coord.cells()):
-                for fp in fp_coord.points:
+            enforce_coordinate_ordering(
+                forecast_cube,
+                ["forecast_period", "forecast_reference_time"],
+                anchor_start=True,
+            )
+            for fp in fp_coord.points:
+                for frt in list(frt_coord.cells()):
                     time_datetimes.append(
                         frt.point._to_real_datetime()
                         + pd.to_timedelta(

improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_apply_quantile_regression_random_forest.py

@@ -31,7 +31,7 @@
         (100, 2, 55, 5, None, 0, {}, False, [0.5, 0.9], [[4.1, 5.1], [4.2, 5.1]]),  # noqa Multiple quantiles
         (1, 1, 55, 5, None, 0, {}, False, [0.5], [6.2, 6.2]),  # noqa Fewer estimators and reduced depth
         (1, 1, 73, 5, None, 0, {}, False, [0.5], [4.2, 6.2]),  # Different random state
-        (2, 2, 55, 5, "log", 10, {}, False, [0.5], [4.1, 5.1]),  # Log transformation
+        (2, 2, 55, 5, "log", 10, {}, False, [0.5], [5.11, 5.64]),  # Log transformation
         (
             2,
             2,
@@ -42,7 +42,7 @@
             {},
             False,
             [0.5],
-            [4.1, 5.1],
+            [5.11, 5.64],
         ),  # Log10 transformation
         (
             2,
@@ -54,7 +54,7 @@
             {},
             False,
             [0.5],
-            [4.1, 5.1],
+            [5.11, 5.64],
         ),  # Square root transformation
         (
             2,
@@ -66,7 +66,7 @@
             {},
             False,
             [0.5],
-            [4.1, 5.1],
+            [5.13, 5.64],
         ),  # Cube root transformation
         (2, 2, 55, 5, None, 0, {"max_samples_leaf": 0.5}, False, [0.5], [5.15, 6.2]),  # noqa # Different criterion
         (2, 5, 55, 5, None, 0, {}, True, [0.5], [5.15, 5.65]),  # Include static data
@@ -155,7 +155,16 @@ def test_load_and_apply_qrf(
     assert result.units == "m s-1"
 
 
-@pytest.mark.parametrize("exception", ["no_model_output", "no_features"])
+@pytest.mark.parametrize(
+    "exception",
+    [
+        "no_model_output",
+        "no_features",
+        "missing_target_feature",
+        "missing_static_feature",
+        "missing_dynamic_feature",
+    ],
+)
 def test_exceptions(
     tmp_path,
     exception,
@@ -207,6 +216,10 @@ def test_exceptions(
         forecast_cube, day_of_training_period, "forecast_reference_time"
     )
 
+    ancil_cube = _create_ancil_file()
+    ancil_filepath = tmp_path / "ancil.nc"
+    save_netcdf(ancil_cube, ancil_filepath)
+
     features_dir = tmp_path / "features"
     features_dir.mkdir(parents=True)
     forecast_filepath = str(features_dir / "forecast.nc")
@@ -219,12 +232,40 @@ def test_exceptions(
         pre_transform_addition=pre_transform_addition,
     )
 
-    if exception == "no_qrf_model":
+    if exception == "no_model_output":
         file_paths = [forecast_filepath]
-        with pytest.raises(ValueError, match="No QRF model found"):
-            plugin.process(file_paths=file_paths)
-
-    if exception == "no_features":
+        result = plugin.process(file_paths=file_paths)
+        assert isinstance(result, Cube)
+        # assert result == forecast_cube
+        assert result.name() == "wind_speed_at_10m"
+        assert result.units == "m s-1"
+        assert result.data.shape == forecast_cube.data.shape
+        assert np.allclose(result.data, forecast_cube.data)
+    elif exception == "no_features":
         file_paths = [model_output]
         with pytest.raises(ValueError, match="No features found"):
             plugin.process(file_paths=file_paths)
+    elif exception == "missing_target_feature":
+        file_paths = [model_output, ancil_filepath]
+        with pytest.raises(ValueError, match="No target forecast provided."):
+            plugin.process(file_paths=file_paths)
+    elif exception == "missing_static_feature":
+        feature_config = {
+            "wind_speed_at_10m": ["mean", "std"],
+            "distance_to_water": ["static"],
+        }
+        plugin.feature_config = feature_config
+        file_paths = [model_output, forecast_filepath]
+        with pytest.raises(ValueError, match="The number of cubes loaded."):
+            plugin.process(file_paths=file_paths)
+    elif exception == "missing_dynamic_feature":
+        feature_config = {
+            "wind_speed_at_10m": ["mean", "std"],
+            "air_temperature": ["mean", "std"],
+        }
+        plugin.feature_config = feature_config
+        file_paths = [model_output, forecast_filepath]
+        with pytest.raises(ValueError, match="The number of cubes loaded."):
+            plugin.process(file_paths=file_paths)
+    else:
+        raise ValueError(f"Unknown exception type: {exception}")