diff --git a/improver/calibration/__init__.py b/improver/calibration/__init__.py
index 6b288ea29a..723395a9ae 100644
--- a/improver/calibration/__init__.py
+++ b/improver/calibration/__init__.py
@@ -24,7 +24,7 @@
         ("altitude", pa.float32()),
         ("blend_time", pa.timestamp("s", "utc")),
         ("forecast_period", pa.int64()),
-        ("forecast_reference_time", pa.int64()),
+        ("forecast_reference_time", pa.timestamp("s", "utc")),
         ("latitude", pa.float32()),
         ("longitude", pa.float32()),
         ("time", pa.timestamp("s", "utc")),
@@ -44,7 +44,7 @@
         ("latitude", pa.float32()),
         ("longitude", pa.float32()),
         ("altitude", pa.float32()),
-        ("time", pa.int64()),
+        ("time", pa.timestamp("s", "utc")),
         ("wmo_id", pa.string()),
         ("ob_value", pa.float32()),
     ]
diff --git a/improver/calibration/dataframe_utilities.py b/improver/calibration/dataframe_utilities.py
index bc103eb983..857442390b 100644
--- a/improver/calibration/dataframe_utilities.py
+++ b/improver/calibration/dataframe_utilities.py
@@ -76,8 +76,7 @@ def _dataframe_column_check(df: DataFrame, compulsory_columns: Sequence) -> None
     if not set(compulsory_columns).issubset(df.columns):
         diff = set(compulsory_columns).difference(df.columns)
         msg = (
-            "The following compulsory column(s) are missing from the "
-            f"DataFrame: {diff}"
+            f"The following compulsory column(s) are missing from the DataFrame: {diff}"
         )
         raise ValueError(msg)
 
@@ -183,11 +182,7 @@ def _fill_missing_entries(df, combi_cols, static_cols, site_id_col):
     df = df.set_index(combi_cols).reindex(new_index).reset_index(level=combi_cols)
 
     # Fill the NaNs within the static columns for each wmo_id.
-    filled_df = (
-        df.groupby(site_id_col)[combi_cols + static_cols]
-        .fillna(method="ffill")
-        .fillna(method="bfill")
-    )
+    filled_df = df.groupby(site_id_col)[combi_cols + static_cols].ffill().bfill()
     df = df.drop(columns=static_cols)
     df = df.merge(filled_df, on=combi_cols)
 
@@ -361,8 +356,7 @@ def get_forecast_representation(df: DataFrame) -> str:
         )
     if len(representations) == 0:
         raise ValueError(
-            f"None of the columns {REPRESENTATION_COLUMNS} "
-            "exist in the input dataset"
+            f"None of the columns {REPRESENTATION_COLUMNS} exist in the input dataset"
         )
     return representations.pop()
 
diff --git a/improver/calibration/load_and_apply_quantile_regression_random_forest.py b/improver/calibration/load_and_apply_quantile_regression_random_forest.py
index c28d9896e6..5647f64d3e 100644
--- a/improver/calibration/load_and_apply_quantile_regression_random_forest.py
+++ b/improver/calibration/load_and_apply_quantile_regression_random_forest.py
@@ -12,7 +12,9 @@
 import iris
 import joblib
 import numpy as np
+import pandas as pd
 from iris.cube import Cube, CubeList
+from iris.pandas import as_data_frame
 from quantile_forest import RandomForestQuantileRegressor
 
 from improver import PostProcessingPlugin
@@ -22,11 +24,12 @@
 )
 from improver.ensemble_copula_coupling.ensemble_copula_coupling import (
     RebadgePercentilesAsRealizations,
-    RebadgeRealizationsAsPercentiles,
 )
 from improver.ensemble_copula_coupling.utilities import choose_set_of_percentiles
 from improver.utilities.cube_checker import assert_spatial_coords_match
 
+iris.FUTURE.pandas_ndim = True
+
 
 class LoadAndApplyQRF(PostProcessingPlugin):
     """Load and apply the trained Quantile Regression Random Forest (QRF) model."""
@@ -34,7 +37,7 @@ class LoadAndApplyQRF(PostProcessingPlugin):
     def __init__(
         self,
         feature_config: dict[str, list[str]],
-        target_cube_name: str,
+        target_cf_name: str,
         transformation: Optional[str] = None,
         pre_transform_addition: Optional[float] = None,
     ):
@@ -52,16 +55,16 @@ def __init__(
                 water cube, then the value should state "static". This will ensure
                 the cube's data is used as the feature.
                 The config will have the structure:
-                "DYNAMIC_VARIABLE_NAME": ["FEATURE1", "FEATURE2"] e.g:
+                "DYNAMIC_VARIABLE_CF_NAME": ["FEATURE1", "FEATURE2"] e.g:
                 {
                 "air_temperature": ["mean", "std", "altitude"],
                 "visibility_at_screen_level": ["mean", "std"]
                 "distance_to_water": ["static"],
                 }
-            target_cube_name (str):
-                A string containing the cube name of the forecast to be
-                calibrated. This will be used to separate it from the rest of the
-                dynamic predictors, if present.
+            target_cf_name (str):
+                A string containing the CF name of diagnostic to be calibrated. This
+                will be used to separate it from the rest of the dynamic predictors,
+                if present.
             transformation (str):
                 Transformation to be applied to the data before fitting.
             pre_transform_addition (float):
@@ -69,7 +72,7 @@ def __init__(
 
         """
         self.feature_config = feature_config
-        self.target_cube_name = target_cube_name
+        self.target_cf_name = target_cf_name
         self.transformation = transformation
         self.pre_transform_addition = pre_transform_addition
 
@@ -110,7 +113,7 @@ 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_constraint = iris.Constraint(name=self.target_cf_name)
         forecast_cube = cube_inputs.extract(forecast_constraint)
 
         if forecast_cube:
@@ -119,7 +122,7 @@ def _get_inputs(
             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}'."
+                f"The target is '{self.target_cf_name}'."
             )
             raise ValueError(msg)
 
@@ -141,7 +144,7 @@ def _get_inputs(
             return forecast_cube
 
         # If target diagnostic not a feature in the training then remove.
-        if self.target_cube_name not in self.feature_config.keys():
+        if self.target_cf_name not in self.feature_config.keys():
             cube_inputs.remove(forecast_cube)
 
         return cube_inputs, forecast_cube, qrf_model
@@ -189,42 +192,40 @@ def _percentiles_to_realizations(cube_inputs: CubeList) -> CubeList:
         return realization_cube_inputs
 
     @staticmethod
-    def _organise_cubes(
-        cube_inputs: CubeList, forecast_cube: Cube
-    ) -> tuple[CubeList, Cube]:
-        """Promote the forecast period and forecast reference time coordinates to be
-        dimension coordinates, if present, on the feature cubes and the template
-        forecast cube.
+    def _cube_to_dataframe(cube_inputs: CubeList) -> pd.DataFrame:
+        """Convert cube inputs to a pandas DataFrame.
 
         Args:
-            cube_inputs: CubeList of feature cubes, which may include the forecast to be
-            forecast_cube: Forecast cube for use as a template.
-
+            cube_inputs: List of cubes containing the features and the forecast to be
+                calibrated.
         Returns:
-            Feature cubes and template cube with forecast period and
-            forecast reference time promoted to dimension coordinates.
+            DataFrame containing the data from the cubes, with auxiliary coordinates
+            included as columns.
         """
-        # Ensure that forecast_period is a dimension on all cubes.
-        fp_dim_cube_inputs = iris.cube.CubeList([])
-        for feature_cube in cube_inputs:
-            if feature_cube.coords("forecast_period", dim_coords=False):
-                feature_cube = iris.util.new_axis(feature_cube, "forecast_period")
-            if feature_cube.coords("forecast_reference_time", dim_coords=False):
-                feature_cube = iris.util.new_axis(
-                    feature_cube, "forecast_reference_time"
-                )
-            fp_dim_cube_inputs.append(feature_cube)
-        cube_inputs = fp_dim_cube_inputs
-
-        # Ensure the forecast cube has the same dimensions as the features
-        template_forecast_cube = iris.util.new_axis(forecast_cube, "forecast_period")
-        template_forecast_cube = iris.util.new_axis(
-            template_forecast_cube, "forecast_reference_time"
-        )
+        # Convert the first cube to a DataFrame.
+        df = as_data_frame(cube_inputs[0], add_aux_coords=True).reset_index()
+
+        # Iteratively convert remaining cubes to DataFrame and merge.
+        for cube in cube_inputs[1:]:
+            temporary_df = as_data_frame(cube, add_aux_coords=True).reset_index()
+            possible_columns = [
+                "wmo_id",
+                "time",
+                "forecast_reference_time",
+                "forecast_period",
+            ]
+            merge_columns = [
+                col for col in possible_columns if col in temporary_df.columns
+            ]
+            df = df.merge(
+                temporary_df[merge_columns + [cube.name()]],
+                on=merge_columns,
+                how="left",
+            )
 
-        # Check that the grids are the same for all dynamic predictors and the forecast
-        assert_spatial_coords_match(cube_inputs)
-        return cube_inputs, template_forecast_cube
+        for column in ["forecast_reference_time", "time"]:
+            df[column] = df[column].apply(lambda x: x._to_real_datetime())
+        return df
 
     def process(
         self,
@@ -254,24 +255,21 @@ def process(
         template_forecast_cube = forecast_cube.copy()
         if forecast_cube.coords("percentile"):
             percentiles = self._compute_percentiles(forecast_cube.copy(), "percentile")
-            cube_inputs = self._percentiles_to_realizations(cube_inputs.copy())
-            template_forecast_cube = RebadgePercentilesAsRealizations()(
-                template_forecast_cube
-            )
         elif forecast_cube.coords("realization"):
             percentiles = self._compute_percentiles(forecast_cube.copy(), "realization")
 
-        cube_inputs, template_forecast_cube = self._organise_cubes(
-            cube_inputs, template_forecast_cube
-        )
+        assert_spatial_coords_match(cube_inputs)
+        df = self._cube_to_dataframe(cube_inputs)
 
-        result = ApplyQuantileRegressionRandomForests(
+        calibrated_forecast = ApplyQuantileRegressionRandomForests(
+            target_name=self.target_cf_name,
             feature_config=self.feature_config,
             quantiles=percentiles,
             transformation=self.transformation,
             pre_transform_addition=self.pre_transform_addition,
-        )(qrf_model, cube_inputs, template_forecast_cube)
+        )(qrf_model, df)
+        calibrated_forecast_cube = template_forecast_cube.copy(
+            data=np.broadcast_to(calibrated_forecast.T, template_forecast_cube.shape)
+        )
 
-        if forecast_cube.coords("percentile"):
-            result = RebadgeRealizationsAsPercentiles()(result)
-        return result
+        return calibrated_forecast_cube
diff --git a/improver/calibration/load_and_train_quantile_regression_random_forest.py b/improver/calibration/load_and_train_quantile_regression_random_forest.py
index 8dd13c4390..3498eaab8e 100644
--- a/improver/calibration/load_and_train_quantile_regression_random_forest.py
+++ b/improver/calibration/load_and_train_quantile_regression_random_forest.py
@@ -14,17 +14,17 @@
 import pandas as pd
 import pyarrow as pa
 import pyarrow.parquet as pq
+from iris.pandas import as_data_frame
 
 from improver import PostProcessingPlugin
 from improver.calibration import FORECAST_SCHEMA, TRUTH_SCHEMA
-from improver.calibration.dataframe_utilities import (
-    forecast_and_truth_dataframes_to_cubes,
-)
 from improver.calibration.quantile_regression_random_forest import (
     TrainQuantileRegressionRandomForests,
 )
 from improver.utilities.load import load_cube
 
+iris.FUTURE.pandas_ndim = True
+
 
 class LoadAndTrainQRF(PostProcessingPlugin):
     """Plugin to load and train a Quantile Regression Random Forest (QRF) model."""
@@ -33,26 +33,30 @@ def __init__(
         self,
         feature_config: dict[str, list[str]],
         target_diagnostic_name: str,
+        target_cf_name: str,
         forecast_periods: str,
         cycletime: str,
         training_length: int,
-        experiment: str = None,
+        experiment: Optional[str] = None,
         n_estimators: int = 100,
-        max_depth: int = None,
-        random_state: int = None,
-        transformation: str = None,
+        max_depth: Optional[int] = None,
+        max_samples: Optional[float] = None,
+        random_state: Optional[int] = None,
+        transformation: Optional[str] = None,
         pre_transform_addition: float = 0,
         compression: int = 5,
     ):
         """Initialise the LoadAndTrainQRF plugin."""
         self.feature_config = feature_config
         self.target_diagnostic_name = target_diagnostic_name
+        self.target_cf_name = target_cf_name
         self.forecast_periods = forecast_periods
         self.cycletime = cycletime
         self.training_length = training_length
         self.experiment = experiment
         self.n_estimators = n_estimators
         self.max_depth = max_depth
+        self.max_samples = max_samples
         self.random_state = random_state
         self.transformation = transformation
         self.pre_transform_addition = pre_transform_addition
@@ -83,24 +87,24 @@ def _split_cubes_and_parquet_files(
 
         # file extraction loop:
         for file_path in file_paths:
-            try:
-                cube = load_cube(str(file_path))
-                cube_inputs.append(cube)
-            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:
-                        pq.read_schema(file).field("forecast_period")
-                        forecast_table_path = file_path
-                    except KeyError:
-                        truth_table_path = file_path
-                    if forecast_table_path and truth_table_path:
-                        break
-            except OSError:
+            if not Path(file_path).exists():
                 # This will occur when the filepath does not exist. In this case,
                 # return None.
                 return None, None, None
+            elif Path(file_path).is_dir():
+                try:
+                    example_file_path = next(Path(file_path).glob("**/*.parquet"))
+                except StopIteration:
+                    # If no parquet files are found, return None.
+                    return None, None, None
+                try:
+                    pq.read_schema(example_file_path).field("forecast_period")
+                    forecast_table_path = file_path
+                except KeyError:
+                    truth_table_path = file_path
+            else:
+                cube = load_cube(str(file_path))
+                cube_inputs.append(cube)
 
         if len(self.feature_config.keys()) not in [
             len(cube_inputs),
@@ -142,7 +146,6 @@ def _read_parquet_files(
                 fields.
         """
         cycletimes = []
-
         for forecast_period in forecast_periods:
             # Load forecasts from parquet file filtering by diagnostic and blend_time.
             forecast_period_td = pd.Timedelta(int(forecast_period), unit="seconds")
@@ -166,23 +169,20 @@ def _read_parquet_files(
             ]
         ]
 
-        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
+        example_file_path = next(Path(forecast_table_path).glob("**/*.parquet"))
+        if pq.read_schema(example_file_path).get_all_field_indices("percentile"):
+            altered_schema = FORECAST_SCHEMA
+        elif pq.read_schema(example_file_path).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)
 
         forecast_df = pd.read_parquet(
             forecast_table_path,
@@ -191,16 +191,18 @@ def _read_parquet_files(
             engine="pyarrow",
         )
 
-        # Convert df columns from ms to pandas timestamp object to work with existing
-        # code
+        forecast_df = forecast_df[
+            forecast_df["forecast_period"].isin(np.array(forecast_periods) * 1e9)
+        ]
+        # Convert df columns from ns to pandas timestamp object.
         for column in ["time", "forecast_reference_time", "blend_time"]:
             forecast_df[column] = pd.to_datetime(
                 forecast_df[column], unit="ns", utc=True
             )
-        forecast_df["forecast_period"] = pd.to_timedelta(
-            forecast_df["forecast_period"], unit="ns"
-        )
-        forecast_df["period"] = pd.to_timedelta(forecast_df["period"], unit="ns")
+        for column in ["forecast_period", "period"]:
+            forecast_df[column] = pd.to_timedelta(forecast_df[column], unit="ns")
+
+        forecast_df = forecast_df.rename(columns={"forecast": self.target_cf_name})
 
         # Load truths from parquet file filtering by diagnostic.
         filters = [[("diagnostic", "==", self.target_diagnostic_name)]]
@@ -218,111 +220,76 @@ def _read_parquet_files(
             raise IOError(msg)
         return forecast_df, truth_df
 
-    def _dataframe_to_cubes(
-        self,
-        forecast_df: pd.DataFrame,
-        truth_df: pd.DataFrame,
-        forecast_periods: list[int],
-    ) -> tuple[iris.cube.Cube, iris.cube.Cube]:
-        """Convert the forecast and truth dataframes to cubes at each forecast period
-        required.
+    @staticmethod
+    def _check_matching_times(
+        forecast_df: pd.DataFrame, truth_df: pd.DataFrame
+    ) -> list[pd.Timestamp]:
+        """Find the intersecting times available within the forecast and truth
+        DataFrames.
 
         Args:
             forecast_df: DataFrame containing the forecast data.
             truth_df: DataFrame containing the truth data.
-            forecast_periods: List of forecast periods in seconds.
-
         Returns:
-            Tuple containing:
-                - Cube containing the forecast data.
-                - Cube containing the truth data.
-
-        Raises:
-            ValueError: The forecast has failed to concatenate into a single cube.
+            List of intersecting times as pandas Timestamp objects.
         """
-        forecast_cubes = iris.cube.CubeList([])
-        truth_cubes = iris.cube.CubeList([])
-
-        for forecast_period in forecast_periods:
-            forecast_cube, truth_cube = forecast_and_truth_dataframes_to_cubes(
-                forecast_df,
-                truth_df,
-                self.cycletime,
-                forecast_period,
-                self.training_length,
-                experiment=self.experiment,
+        # Calling unique() on the time column is quicker than relying upon set() to
+        # find the unique times.
+        return list(
+            set(forecast_df["time"].unique()).intersection(
+                set(truth_df["time"].unique())
             )
+        )
 
-            if forecast_cube is None or truth_cube is None:
-                continue
-
-            if not forecast_cube.coords("realization", dim_coords=True):
-                forecast_cube = iris.util.new_axis(forecast_cube, "realization")
-
-            for forecast_slice in forecast_cube.slices_over("forecast_reference_time"):
-                forecast_cubes.append(forecast_slice)
-
-            for truth_slice in truth_cube.slices_over("time"):
-                truth_slice = iris.util.new_axis(truth_slice, "time")
-                
-                # Multiple forecasts can match to the same observation. This check
-                # ensures that we do not add the same truth slice multiple times.
-                if truth_slice not in truth_cubes:
-                    truth_cubes.append(truth_slice)
-
-        if not forecast_cubes or not truth_cubes:
-            return None, None
-
-        truth_cube = truth_cubes.concatenate_cube()
-        forecast_cube = forecast_cubes.merge()
+    def _add_features_to_df(
+        self, forecast_df: pd.DataFrame, cube_inputs: iris.cube.CubeList
+    ) -> pd.DataFrame:
+        """Add features to the forecast DataFrame based on the feature configuration.
 
-        # concatenate_cube() can fail for the forecast_cube, even though calling
-        # concatenate() results in a single cube. This check ensures the concatenation
-        # was successful.
-        if len(forecast_cube) == 1:
-            forecast_cube = forecast_cube[0]
-        else:
-            msg = "Concatenating the forecast has failed to create a single cube."
-            raise ValueError(msg)
-        
-        # Promote the forecast_reference_time coord to a dimension coordinate if
-        # the forecast_period is one.
-        if forecast_cube.coord_dims("forecast_period") and not forecast_cube.coord_dims("forecast_reference_time"):
-            forecast_cube = iris.util.new_axis(forecast_cube, "forecast_reference_time")
+        Args:
+            forecast_df: DataFrame containing the forecast data.
+            cube_inputs: List of cubes containing additional features.
 
-        return forecast_cube, truth_cube
+        Returns:
+            DataFrame with additional features added.
+        """
+        for feature_name, feature_list in self.feature_config.items():
+            for feature in feature_list:
+                if feature == "static":
+                    # Use the cube's data directly as a feature.
+                    constr = iris.Constraint(name=feature_name)
+                    feature_cube = cube_inputs.extract_cube(constr)
+                    feature_df = as_data_frame(feature_cube, add_aux_coords=True)
+                    forecast_df = forecast_df.merge(
+                        feature_df[["wmo_id", feature_name]], on=["wmo_id"], how="left"
+                    )
+        return forecast_df
 
     @staticmethod
     def filter_bad_sites(
-        forecast_cube: iris.cube.Cube,
-        truth_cube: iris.cube.Cube,
-        cube_inputs: iris.cube.CubeList,
-    ) -> tuple[iris.cube.Cube, iris.cube.Cube, iris.cube.CubeList]:
+        forecast_df: pd.DataFrame,
+        truth_df: pd.DataFrame,
+    ) -> tuple[pd.DataFrame, pd.DataFrame]:
         """Remove sites that have NaNs in the data.
 
         Args:
-            forecast_cube: Cube containing the forecast data.
-            truth_cube: Cube containing the truth data.
-            cube_inputs: List of additional feature cubes.
+            feature_df: DataFrame containing the forecast data with features.
+            truth_df: DataFrame containing the truth data.
 
         Returns:
             Tuple containing:
-                - Cube containing the forecast data with bad sites removed.
-                - Cube containing the truth data with bad sites removed.
-                - List of additional feature cubes with bad sites removed.
+                - DataFrame containing the forecast data with bad sites removed.
+                - DataFrame containing the truth data with bad sites removed.
         """
-        nan_mask = np.any(np.isnan(truth_cube.data), axis=truth_cube.coord_dims("time"))
-        all_site_ids = truth_cube.coord("wmo_id").points
-        bad_site_ids = all_site_ids[nan_mask]
-        constr = iris.Constraint(wmo_id=lambda cell: cell not in bad_site_ids.tolist())
-        truth_cube = truth_cube.extract(constr)
-        forecast_cube = forecast_cube.extract(constr)
-        feature_cube_inputs = iris.cube.CubeList([])
-        for cube in cube_inputs:
-            cube = cube.extract(constr)
-            feature_cube_inputs.append(cube)
-
-        return forecast_cube, truth_cube, feature_cube_inputs
+        truth_df.dropna(
+            subset=["latitude", "longitude", "altitude", "ob_value"], inplace=True
+        )
+
+        wmo_ids = set(forecast_df["wmo_id"]).intersection(set(truth_df["wmo_id"]))
+
+        forecast_df = forecast_df[forecast_df["wmo_id"].isin(wmo_ids)]
+        truth_df = truth_df[truth_df["wmo_id"].isin(wmo_ids)]
+        return forecast_df, truth_df
 
     def process(
         self,
@@ -344,55 +311,6 @@ def process(
                 - The path to a Parquet file containing the forecasts to be used
                 for calibration.
                 - Optionally, paths to NetCDF files containing additional predictors.
-            feature_config (dict):
-                Feature configuration defining the features to be used for quantile
-                regression. The configuration is a dictionary of strings, where the
-                keys are the names of the input cube(s) supplied, and the values are
-                a list. This list can contain both computed features, such as the mean
-                or standard deviation (std), or static features, such as the altitude.
-                The computed features will be computed using the cube defined in the
-                dictionary key. If the key is the feature itself e.g. a distance to
-                water cube, then the value should state "static". This will ensure
-                the cube's data is used as the feature. The config will have the
-                structure:
-                    "DYNAMIC_VARIABLE_NAME": ["FEATURE1", "FEATURE2"] e.g:
-                    {
-                    "air_temperature": ["mean", "std", "altitude"],
-                    "visibility_at_screen_level": ["mean", "std"]
-                    "distance_to_water": ["static"],
-                    }
-            target_diagnostic_name (str):
-                A string containing the diagnostic name of the forecast to be
-                calibrated. This will be used to filter the target forecast and truth
-                dataframes.
-            forecast_period (int):
-                Range of forecast periods to be calibrated in hours in the form:
-                "start:end:interval" e.g. "6:18:6" or a single forecast period e.g. "6".
-                The end value is exclusive, so "6:18:6" will calibrate the 6 and 12
-                hours.
-            cycletime (str):
-                Cycletime of the forecast to be calibrated in a format similar to
-                20170109T0000Z. This is used to filter the correct blendtimes from
-                the dataframe on load.
-            training_length (int):
-                The length of the training period in days.
-            experiment (str):
-                The name of the experiment (step) that calibration is applied to.
-            n_estimators (int):
-                Number of trees in the forest.
-            max_depth (int):
-                Maximum depth of the tree.
-            random_state (int):
-                Random seed for reproducibility.
-            transformation (str):
-                Transformation to be applied to the data before fitting.
-                Supported transformations are "log", "log10", "sqrt", and "cbrt".
-            pre_transform_addition (float):
-                Value to be added before transformation. This is useful for ensuring
-                that the data is positive before applying a transformation such as log.
-            compression (int):
-                Compression level for saving the model. Please see the joblib
-                documentation for more information on compression levels.
             model_output (str):
                 Full path including model file name that will store the pickled model.
 
@@ -421,30 +339,22 @@ def process(
         forecast_df, truth_df = self._read_parquet_files(
             forecast_table_path, truth_table_path, forecast_periods
         )
-
-        forecast_cube, truth_cube = self._dataframe_to_cubes(
-            forecast_df, truth_df, forecast_periods
-        )
-        if forecast_cube is None or truth_cube is None:
+        intersecting_times = self._check_matching_times(forecast_df, truth_df)
+        if len(intersecting_times) == 0:
             return None
 
-        # 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():
-            if feature_name == forecast_cube[0].name():
-                cube_inputs.append(forecast_cube)
-
-        forecast_cube, truth_cube, feature_cube_inputs = self.filter_bad_sites(
-            forecast_cube, truth_cube, cube_inputs
-        )
+        forecast_df = self._add_features_to_df(forecast_df, cube_inputs)
+        forecast_df, truth_df = self.filter_bad_sites(forecast_df, truth_df)
 
         TrainQuantileRegressionRandomForests(
+            target_name=self.target_cf_name,
             feature_config=self.feature_config,
             n_estimators=self.n_estimators,
             max_depth=self.max_depth,
+            max_samples=self.max_samples,
             random_state=self.random_state,
             transformation=self.transformation,
             pre_transform_addition=self.pre_transform_addition,
             compression=self.compression,
             model_output=model_output,
-        )(forecast_cube, truth_cube, feature_cube_inputs)
+        )(forecast_df, truth_df)
diff --git a/improver/calibration/quantile_regression_random_forest.py b/improver/calibration/quantile_regression_random_forest.py
index d24542b6fd..cf43788b27 100644
--- a/improver/calibration/quantile_regression_random_forest.py
+++ b/improver/calibration/quantile_regression_random_forest.py
@@ -6,284 +6,161 @@
 
 from typing import Optional
 
-import iris
 import joblib
 import numpy as np
 import pandas as pd
-from iris.cube import Cube, CubeList
 from quantile_forest import RandomForestQuantileRegressor
 
 from improver import BasePlugin, PostProcessingPlugin
 from improver.constants import DAYS_IN_YEAR, HOURS_IN_DAY
-from improver.utilities.cube_manipulation import enforce_coordinate_ordering
-from improver.ensemble_copula_coupling.ensemble_copula_coupling import (
-            RebadgeRealizationsAsPercentiles,
-        )
-
-
-def _expand_dims(
-    template_cube: Cube, expansion_array: np.ndarray, features: list[str]
-) -> np.ndarray:
-    """Expand dimensions of the expansion_array to match the dimensions of the template
-    cube. For example, the dimensions of the expansion_array may only span a portion
-    of the dimensions present in the template cube. If the expansion_array has
-    singleton dimensions that are not present in the template cube, these singleton
-    dimension will be squeezed away.
-
-    Args:
-        template_cube: Cube that acts as a template for the output. The dimension
-            coordinates and coordinates associated with the dimension coordinates on
-            this cube will be used to expand the expansion_array.
-        expansion_array: Array that will be expanded to match the dimensions of the
-            template_cube.
-        features: List of feature names that are present in the template_cube. These
-            may correspond to dimension coordinates or coordinates associated with
-            dimension coordinates or scalar coordinates.
-
-    Returns:
-        Array with dimensions expanded to match the template_cube.
-    """
-    # The aim here is to find the names of the dimension coordinates but check
-    # whether the name of any coordinate associated with a dimension coordinate matches
-    # the features provided. If there is a match, then use the feature name, rather than
-    # the dimension coordinate name.
-    dim_coord_names = [c.name() for c in template_cube.coords(dim_coords=True)]
-    refined_dim_coord_names = []
-    for coord in template_cube.coords():
-        # Ignore scalar coordinates.
-        if len(template_cube.coord_dims(coord)) > 0:
-            associated_coords = [
-                c.name()
-                for c in template_cube.coords(
-                    dimensions=template_cube.coord_dims(coord)
-                )
-            ]
-            feature_associated_coords = list(
-                set(associated_coords).intersection(features)
-            )
-            if feature_associated_coords:
-                refined_dim_coord_names.append(feature_associated_coords[0])
-            elif coord.name() in dim_coord_names:
-                refined_dim_coord_names.append(coord.name())
-    refined_dim_coord_names = list(set(refined_dim_coord_names))
-
-    expansion_dim_names = list(set(refined_dim_coord_names) - set(features))
-
-    # If the expansion_dims and refined_dim_coord_names match, this implies that
-    # the features provided are not dimension coordinates nor associated with any
-    # dimension coordinate. In this case, any singleton dimensions on the
-    # expansion_array should be removed, as it isn't relevant for defining the shape
-    # of the template cube.
-    if expansion_dim_names == refined_dim_coord_names:
-        expansion_array = np.squeeze(expansion_array)
-
-    dims = []
-    for dim_name in expansion_dim_names:
-        dims.extend(template_cube.coord_dims(dim_name))
-
-    return np.expand_dims(expansion_array, dims)
 
 
 def prep_feature(
-    template_cube: Cube,
-    feature_cube: Cube,
-    feature: str,
-) -> np.ndarray:
-    """Prepare the feature values for the quantile regression random forest model.
-    Each expected feature is handled separately. The key steps are:
-        - Collapse the template cube over the realization dimension to get a shape
-        template.
-        - Extract the feature values from the feature cube provided.
-        - Handle different types of features separately (mean, std, spatial coordinates,
-        forecast period, model weights, day of year, hour of day, and static features)
-        to ensure they are broadcasted correctly to match the template cube shape.
-        Time coordinates are the most complex to handle as they can be scalar, span a
-        single dimension or span multiple dimensions.
-        - Flatten the feature values to create a 1D array.
-        - Set the output dtype based on the feature type.
-        - Return the flattened array of feature values.
+    df: pd.DataFrame,
+    variable_name: str,
+    feature_name: str,
+) -> pd.DataFrame:
+    """Prepare features that require computation from the input DataFrame. Options
+    available are mean and standard deviation of the input feature, the
+    day of year, sine of day of year, cosine of day of year, hour of day,
+    sine of hour of day and cosine of hour of day. When computing the mean or standard
+    deviation, these will be computed over either the percentile or realization column,
+    depending upon which is available.
 
     Args:
-        template_cube (cube):
-            The forecast cube that acts only as a template.
-        feature_cube (cube):
-            Feature cube that is either static or dynamic.
-        feature:
-            The feature to be extracted from the associated feature_cube.
-            If "static" then the cube itself acts as the feature.
+        df: Input DataFrame.
+        variable_name: Name of the variable to be used for the computation.
+        feature_name: Feature to be computed. Options are "mean", "std", "day_of_year",
+            "day_of_year_sin", "day_of_year_cos", "hour_of_day",
+            "hour_of_day_sin" and "hour_of_day_cos".
     Returns:
-        feature_values (numpy.ndarray):
-            Flattened array of feature values.
+        df: DataFrame with the computed feature added.
     """
-    # Collapse the template cube to get a shape template.
-    collapsed_cube = template_cube.collapsed(["realization"], iris.analysis.MEAN)
-
-    # Handle each expected feature type separately.
-    if "mean" == feature:
-        feature_values = feature_cube.collapsed(
-            ["realization"], iris.analysis.MEAN
-        ).data.flatten()
-    elif "std" == feature:
-        feature_values = feature_cube.collapsed(
-            ["realization"], iris.analysis.STD_DEV
-        ).data.flatten()
-    elif feature in ["latitude", "longitude", "altitude"]:
-        coord_multidim = _expand_dims(
-            collapsed_cube, feature_cube.coord(feature).points, ["spot_index"]
+    if feature_name in ["mean", "std"]:
+        representation_name = [
+            n for n in ["percentile", "realization"] if n in df.columns
+        ][0]
+        groupby_cols = ["forecast_reference_time", "forecast_period", "wmo_id"]
+        subset_cols = [*groupby_cols] + [
+            representation_name,
+            variable_name,
+        ]
+        # For a subset of the input DataFrame compute the mean or standard deviation
+        # over the representation column, grouped by the groupby columns.
+        if feature_name == "mean":
+            subset_df = df[subset_cols].groupby(groupby_cols).mean()
+        elif feature_name == "std":
+            subset_df = df[subset_cols].groupby(groupby_cols).std()
+
+        subset_df = subset_df.reset_index()
+        # Rename the column to distinguish the computed feature from the original.
+        subset_df.rename(
+            columns={variable_name: f"{variable_name}_{feature_name}"}, inplace=True
         )
-        feature_values = np.broadcast_to(coord_multidim, collapsed_cube.shape).flatten()
-    elif feature == "forecast_period":
-        if len(feature_cube.coord_dims("forecast_period")) == 1:
-            coord_multidim = _expand_dims(
-                collapsed_cube,
-                feature_cube.coord("forecast_period").points,
-                ["forecast_period"],
-            )
-        else:
-            coord_multidim = feature_cube.coord("forecast_period").points
-
-        feature_values = np.broadcast_to(coord_multidim, collapsed_cube.shape).flatten()
-    elif feature == "model_weights":
-        if len(feature_cube.coord_dims("forecast_period")) == 1:
-            coord_multidim = _expand_dims(
-                collapsed_cube,
-                feature_cube.coord("model_weights").points,
-                ["forecast_period"],
-            )
+        # Merge the computed feature back into the original DataFrame.
+        df = df.merge(
+            subset_df[groupby_cols + [f"{variable_name}_{feature_name}"]],
+            on=groupby_cols,
+            how="left",
+        )
+
+    elif feature_name in ["day_of_year", "day_of_year_sin", "day_of_year_cos"]:
+        # For a large DataFrame, the strftime("%j") computation can take a noticeable
+        # amount of time, so this computation is done once for each unique time
+        # and then merged back into the DataFrame.
+        doy_df = pd.DataFrame({"time": df["time"].unique()})
+        doy_df["day_of_year"] = np.array(doy_df["time"].dt.strftime("%j"), np.int32)
+
+        if feature_name == "day_of_year_sin":
+            doy_df[feature_name] = np.sin(
+                2 * np.pi * doy_df["day_of_year"].values / (DAYS_IN_YEAR + 1)
+            ).astype(np.float32)
+        elif feature_name == "day_of_year_cos":
+            doy_df[feature_name] = np.cos(
+                2 * np.pi * doy_df["day_of_year"].values / (DAYS_IN_YEAR + 1)
+            ).astype(np.float32)
+        df = df.merge(doy_df[["time", feature_name]], on="time", how="left")
+    elif feature_name in ["hour_of_day", "hour_of_day_sin", "hour_of_day_cos"]:
+        # For hour_of_day, unlike day_of_year, the hour attribute doesn't require
+        # computation, therefore there is no benefit to creating the separate DataFrame
+        # and merging it back into the DataFrame.
+        if df["time"].nunique() == 1:
+            hour_of_day = np.int32(df["time"].iloc[0].hour)
         else:
-            coord_multidim = feature_cube.coord("model_weights").points
-
-        feature_values = np.broadcast_to(coord_multidim, collapsed_cube.shape).flatten()
-    elif feature in ["day_of_year", "day_of_year_sin", "day_of_year_cos"]:
-        frt_dims = feature_cube.coord_dims("forecast_reference_time")
-        fp_dims = feature_cube.coord_dims("forecast_period")
-        frt_coord = feature_cube.coord("forecast_reference_time")
-        fp_coord = feature_cube.coord("forecast_period")
-
-        if len(frt_dims) == 0 and len(fp_dims) == 0:
-            time_point = frt_coord.cell(0).point._to_real_datetime() + pd.to_timedelta(
-                fp_coord.points[0], unit=str(fp_coord.units)
-            )
-            day_of_year = np.array([np.int32(time_point.strftime("%j"))])
-        elif frt_dims == fp_dims:
-            # Forecast reference time and forecast period share a dimension coordinate.
-            # Forecast reference time and forecast period must be mixed together
-            # along one dimension.
-            day_of_year = []
-            for frt_cell, fp_point in zip(frt_coord.cells(), fp_coord.points):
-                time_point = frt_cell.point._to_real_datetime() + pd.to_timedelta(
-                    fp_point, unit=str(fp_coord.units)
-                )
-                day_of_year.append(np.int32(time_point.strftime("%j")))
-            day_of_year = np.array(day_of_year)
-            # frt_dims and fp_dims are the same, so we can choose one of them to pop.
-            day_of_year = _expand_dims(
-                collapsed_cube, day_of_year, ["forecast_reference_time"]
+            hour_of_day = np.array(df["time"].dt.hour, dtype=np.int32)
+        if feature_name == "hour_of_day":
+            feature_values = hour_of_day
+        elif feature_name == "hour_of_day_sin":
+            feature_values = np.sin(2 * np.pi * hour_of_day / HOURS_IN_DAY).astype(
+                np.float32
             )
-        else:
-            # Forecast reference time and forecast period are different dimensions.
-            day_of_year = np.zeros((len(frt_coord.points), len(fp_coord.points)))
-            for i, frt_cell in enumerate(frt_coord.cells()):
-                for j, fp_point in enumerate(fp_coord.points):
-                    time_point = frt_cell.point._to_real_datetime() + pd.to_timedelta(
-                        fp_point, unit=str(fp_coord.units)
-                    )
-                    day_of_year[i, j] = time_point.strftime("%j")
-
-            day_of_year = _expand_dims(
-                collapsed_cube,
-                day_of_year,
-                ["forecast_reference_time", "forecast_period"],
+        elif feature_name == "hour_of_day_cos":
+            feature_values = np.cos(2 * np.pi * hour_of_day / HOURS_IN_DAY).astype(
+                np.float32
             )
+        df[feature_name] = feature_values
+    return df
 
-        feature_values = np.broadcast_to(day_of_year, collapsed_cube.shape).flatten()
-        if feature == "day_of_year_sin":
-            feature_values = np.sin(2 * np.pi * feature_values / (DAYS_IN_YEAR + 1))
-        elif feature == "day_of_year_cos":
-            feature_values = np.cos(2 * np.pi * feature_values / (DAYS_IN_YEAR + 1))
-    elif feature in ["hour_of_day", "hour_of_day_sin", "hour_of_day_cos"]:
-        frt_dims = feature_cube.coord_dims("forecast_reference_time")
-        fp_dims = feature_cube.coord_dims("forecast_period")
-        frt_coord = feature_cube.coord("forecast_reference_time")
-        fp_coord = feature_cube.coord("forecast_period")
-
-        if len(frt_dims) == 0 and len(fp_dims) == 0:
-            time_point = frt_coord.cell(0).point._to_real_datetime() + pd.to_timedelta(
-                fp_coord.points[0], unit=str(fp_coord.units)
-            )
-            hour_of_day = np.array([np.int32(time_point.hour)])
-        elif frt_dims == fp_dims:
-            # Forecast reference time and forecast period share a dimension coordinate.
-            # Forecast reference time and forecast period must be mixed together
-            # along one dimension.
-            hour_of_day = []
-            for frt_cell, fp_point in zip(frt_coord.cells(), fp_coord.points):
-                time_point = frt_cell.point._to_real_datetime() + pd.to_timedelta(
-                    fp_point, unit=str(fp_coord.units)
-                )
-                hour_of_day.append(np.int32(time_point.hour))
-            hour_of_day = np.array(hour_of_day)
 
-            # frt_dims and fp_dims are the same, so we can choose one of them to pop.
-            hour_of_day = _expand_dims(
-                collapsed_cube, hour_of_day, ["forecast_reference_time"]
-            )
-        else:
-            # Forecast reference time and forecast period are different dimensions.
-            hour_of_day = np.zeros((len(frt_coord.points), len(fp_coord.points)))
-            for i, frt_cell in enumerate(frt_coord.cells()):
-                for j, fp_point in enumerate(fp_coord.points):
-                    time_point = frt_cell.point._to_real_datetime() + pd.to_timedelta(
-                        fp_point, unit=str(fp_coord.units)
-                    )
-                    hour_of_day[i, j] = time_point.hour
-
-            hour_of_day = _expand_dims(
-                collapsed_cube,
-                hour_of_day,
-                ["forecast_reference_time", "forecast_period"],
-            )
+def sanitise_forecast_dataframe(
+    df: pd.DataFrame, feature_config: dict[str, list[str]]
+) -> pd.DataFrame:
+    """Sanitise the forecast DataFrame by removing columns that are no longer
+    required. Following the computation of e.g. the mean or standard deviation,
+    the original feature can be removed. The column over which the mean or
+    standard deviation has been computed (e.g. the percentile or realization column)
+    is also removed.
 
-        feature_values = np.broadcast_to(hour_of_day, collapsed_cube.shape).flatten()
-        if feature == "hour_of_day_sin":
-            feature_values = np.sin(2 * np.pi * feature_values / HOURS_IN_DAY)
-        elif feature == "hour_of_day_cos":
-            feature_values = np.cos(2 * np.pi * feature_values / HOURS_IN_DAY)
-    elif feature == "day_of_training_period":
-        if len(feature_cube.coord_dims("day_of_training_period")) == 1:
-            coord_multidim = _expand_dims(
-                collapsed_cube,
-                feature_cube.coord("day_of_training_period").points,
-                ["day_of_training_period"],
-            )
-        else:
-            coord_multidim = feature_cube.coord("day_of_training_period").points
-
-        feature_values = np.broadcast_to(coord_multidim, collapsed_cube.shape).flatten()
-    elif feature == "static":
-        coord_multidim = _expand_dims(collapsed_cube, feature_cube.data, ["spot_index"])
-        feature_values = np.broadcast_to(coord_multidim, collapsed_cube.shape).flatten()
-    else:
-        # If the feature is not one of the expected types, raise an error.
-        msg = f"Feature {feature} is not supported."
-        raise ValueError(msg)
+    Args:
+        df: Input DataFrame, potentially including some computed features.
+        feature_config: Feature configuration defining the features to be used for QRF.
+    """
+    representation_name = [n for n in ["percentile", "realization"] if n in df.columns][
+        0
+    ]
+    collapsed_features = []
+    for key, values in feature_config.items():
+        collapsed_features.extend([key for v in values if v in ["mean", "std"]])
+    collapsed_features = list(set(collapsed_features))
+    # Subset the dataframe by the first value of the representation column
+    # and drop the representation column and any features where the original variable
+    # is no longer required. This reduces the size of the DataFrame e.g. if there are
+    # 3 percentiles initially, the subsetted dataframe will be 1/3 of the size.
+    df = df[df[representation_name] == df[representation_name].iloc[0]]
+    df = df.drop(columns=[representation_name, *collapsed_features])
+    return df
+
+
+def get_required_column_names(
+    df: pd.DataFrame, feature_config: dict[str, list[str]]
+) -> list[str]:
+    """Process the feature_config to return the expected column names that will be
+    used as features with the QRF.
 
-    # Ensure the feature values are returned with the correct dtype.
-    if feature in ["mean", "std", "static"]:
-        feature_values = feature_values.astype(feature_cube.dtype)
-    elif feature in [
-        "day_of_year",
-        "day_of_year_sin",
-        "day_of_year_cos",
-        "hour_of_day",
-        "hour_of_day_sin",
-        "hour_of_day_cos",
-    ]:
-        feature_values = feature_values.astype(np.float32)
-    else:
-        feature_values = feature_values.astype(feature_cube.coord(feature).dtype)
+    Args:
+        df: Input DataFrame.
+        feature_config: Feature configuration defining the features to be used for QRF.
+    Returns:
+        List of expected column names that will be used as features with the QRF.
+    Raises:
+        ValueError: If a feature expected in the feature_config is not present in
+        the DataFrame.
+    """
+    feature_column_names = []
+    for variable_name in feature_config.keys():
+        for feature in feature_config[variable_name]:
+            if feature in ["mean", "std"]:
+                feature_column_names.append(f"{variable_name}_{feature}")
+            elif feature in ["static"]:
+                feature_column_names.append(variable_name)
+            else:
+                feature_column_names.append(feature)
+
+    if len(list(set(feature_column_names) - set(df.columns))) > 0:
+        msg = f"Feature '{feature}' is not supported."
+        raise ValueError(msg)
 
-    return feature_values
+    return feature_column_names
 
 
 def _check_valid_transformation(transformation: str):
@@ -306,9 +183,11 @@ class TrainQuantileRegressionRandomForests(BasePlugin):
 
     def __init__(
         self,
+        target_name: str,
         feature_config: dict[str, list[str]],
         n_estimators: int,
         max_depth: Optional[int] = None,
+        max_samples: Optional[float] = None,
         random_state: Optional[int] = None,
         transformation: Optional[str] = None,
         pre_transform_addition: np.float32 = 0,
@@ -317,28 +196,35 @@ def __init__(
         **kwargs,
     ) -> None:
         """Initialise the plugin.
+
         Args:
+            target_name (str):
+                Name of the target variable to be calibrated e.g. 'air_temperature'.
             feature_config (dict):
                 Feature configuration defining the features to be used for quantile
                 regression. The configuration is a dictionary of strings, where the
-                keys are the names of the input cube(s) supplied, and the values are
-                a list. This list can contain both computed features, such as the
-                mean or standard deviation (std), or static features, such as the
-                altitude. The computed features will be computed using the cube defined
-                in the dictionary key. If the key is the feature itself e.g. a distance
-                to water cube, then the value should state "static". This will ensure
-                the cube's data is used as the feature. The config will have the
-                structure:
-                    "DYNAMIC_VARIABLE_NAME": ["FEATURE1", "FEATURE2"] e.g:
-                    {
-                    "air_temperature": ["mean", "std", "altitude"],
-                    "visibility_at_screen_level": ["mean", "std"]
-                    "distance_to_water": ["static"],
-                    }
+                keys are the names of the columns within the dataframe. Some
+                features may be used as initially provided within the dataframe,
+                whilst others may be computed from the data e.g. mean, std.
+                If the key is the feature itself e.g. distance to water, then the value
+                should state "static". In this case, the name of feature e.g.
+                'distance_to_water' is expected to be a column name in the input
+                dataframe. The config will have the structure:
+                "DYNAMIC_VARIABLE_CF_NAME": ["FEATURE1", "FEATURE2"] e.g.
+                {
+                "air_temperature": ["mean", "std", "altitude"],
+                "visibility_at_screen_level": ["mean", "std"]
+                "distance_to_water": ["static"],
+                }
             n_estimators (int):
                 Number of trees in the forest.
             max_depth (int):
                 Maximum depth of the tree.
+            max_samples (float):
+                If an int, then it is the number of samples to draw to train
+                each tree. If a float, then it is the fraction of samples to draw
+                to train each tree. If None, then each tree contains the same
+                total number of samples as originally provided.
             random_state (int):
                 Random seed for reproducibility.
             transformation (str):
@@ -351,11 +237,13 @@ def __init__(
                 Full path including model file name that will store the pickled model.
             kwargs:
                 Additional keyword arguments for the quantile regression model.
-        """
 
+        """
+        self.target_name = target_name
         self.feature_config = feature_config
         self.n_estimators = n_estimators
         self.max_depth = max_depth
+        self.max_samples = max_samples
         self.random_state = random_state
         self.transformation = transformation
         _check_valid_transformation(self.transformation)
@@ -363,6 +251,7 @@ def __init__(
         self.compression = compression
         self.output = model_output
         self.kwargs = kwargs
+        self.expected_coordinate_order = ["forecast_reference_time", "forecast_period"]
 
     def fit_qrf(
         self, forecast_features: np.ndarray, target: np.ndarray
@@ -381,96 +270,26 @@ def fit_qrf(
         qrf_model = RandomForestQuantileRegressor(
             n_estimators=self.n_estimators,
             max_depth=self.max_depth,
+            max_samples=self.max_samples,
             random_state=self.random_state,
             **self.kwargs,
         )
         qrf_model.fit(forecast_features, target)
         return qrf_model
 
-    @staticmethod
-    def _organise_truth_data(forecast_cube: Cube, truth_cube: Cube) -> list[np.ndarray]:
-        """Organise the truth data, so that the validity time matches the validity
-        time of the forecast. This might mean that the truth data is repeated, if, for
-        example, the forecast has multiple forecast reference times and multiple
-        forecast periods that have the same validity time."""
-        frt_dims = forecast_cube.coord_dims("forecast_reference_time")
-        fp_dims = forecast_cube.coord_dims("forecast_period")
-
-        frt_coord = forecast_cube.coord("forecast_reference_time")
-        fp_coord = forecast_cube.coord("forecast_period")
-
-        time_datetimes = []
-
-        # Handle the case where both forecast reference time and forecast period are
-        # non-dimensional coordinates (i.e., the cube has no time dimensions).
-        if frt_dims is None and fp_dims is None:
-            time_datetimes.append(
-                frt_coord.cell(0).point._to_real_datetime()
-                + pd.to_timedelta(fp_coord.points, unit=str(fp_coord.units))
-            )
-        elif frt_dims == fp_dims:
-            # Handle the case where forecast reference time and forecast period share
-            # the same dimension. This means both are mixed together along one
-            # dimension.
-            for frt, fp in zip(
-                list(frt_coord.cells()),
-                fp_coord.points,
-            ):
-                time_datetimes.append(
-                    frt.point._to_real_datetime()
-                    + pd.to_timedelta(
-                        fp, unit=str(forecast_cube.coord("forecast_period").units)
-                    )
-                )
-        else:
-            # Handle the case where forecast reference time and forecast period are on
-            # different dimensions. This requires iterating over all combinations of
-            # forecast period and forecast reference time.
-            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(
-                            fp, unit=str(forecast_cube.coord("forecast_period").units)
-                        )
-                    )
-
-        truth_data_list = []
-        for time_datetime in time_datetimes:
-            constr = iris.Constraint(time=lambda cell: cell.point == time_datetime)
-            truth_data_list.append(truth_cube.extract(constr).data)
-        return truth_data_list
-
     def process(
         self,
-        forecast_cube: Cube,
-        truth_cube: Cube,
-        feature_cubes: Optional[CubeList] = None,
+        forecast_df: pd.DataFrame,
+        truth_df: pd.DataFrame,
     ) -> None:
         """Train a quantile regression random forests model.
+
         Args:
-            forecast_cube:
-                Cube containing the realization forecasts. If the cube provided
-                contains multiple forecast periods, then the cube is expected to have
-                forecast period, forecast reference time, realization
-                and spot_index as the dimensions. This cube is only used as a template.
-                If the forecast cube is a feature cube, then it should also be provided
-                within the feature_cubes list.
-            truth_cube:
+            forecast_df:
+                DataFrame containing the forecast information and features.
+            truth_df:
                 Cube containing the truths. The truths should have the same validity
-                times as the forecast. If the same validity time occurs multiple times
-                within the forecast cube (i.e. due to e.g. a forecast reference time
-                of 20170102T0000Z and a lead time of T+36 having the same validity time
-                as a forecast reference time of 20170103T0000Z and a lead time of T+6),
-                then the truth data will be repeated.
-            feature_cubes:
-                List of additional feature cubes. The name of the cube should match a
-                key in the feature_config dictionary.
+                times as the forecasts.
 
         References:
             Johnson. (2024). quantile-forest: A Python Package for Quantile
@@ -487,55 +306,37 @@ def process(
             operational ensemble post-processing in France using machine learning,
             Nonlin. Processes Geophys., 27, 329–347,
             https://doi.org/10.5194/npg-27-329-2020, 2020.
-        """
 
+        """
         if self.transformation:
-            forecast_cube.data = getattr(np, self.transformation)(
-                forecast_cube.data + self.pre_transform_addition
-            )
-            truth_cube.data = getattr(np, self.transformation)(
-                truth_cube.data + self.pre_transform_addition
+            forecast_df[self.target_name] = getattr(np, self.transformation)(
+                forecast_df[self.target_name] + self.pre_transform_addition
             )
-
-        # Ensure the forecast cube has the correct dimension ordering for prep_feature.
-        
-        coord_dims = [forecast_cube.coord_dims(c) for c in ["forecast_reference_time", "forecast_period"]]
-        coord_names = []
-        for dim in coord_dims:
-            if len(dim) == 0:
-                continue
-            if len(forecast_cube.coords(dimensions=dim, dim_coords=True)) == 0:
-                continue
-            coord_names.append(
-                forecast_cube.coord(dimensions=dim, dim_coords=True).name()
-            )
-
-        # If the temporal coord is 2D then this ensures sometimes that the coordinates are ordered correctly.
-        # However, if the temporal coord is 1D then iris currently ensures the dimension is always the first dimension.
-        enforce_coordinate_ordering(
-                forecast_cube,
-                [coord_names],
-                anchor_start=True,
+            truth_df["ob_value"] = getattr(np, self.transformation)(
+                truth_df["ob_value"] + self.pre_transform_addition
             )
 
-        feature_values = []
-
-        for feature_name in self.feature_config.keys():
-            feature_cube = feature_cubes.extract(iris.Constraint(feature_name))
-            if not feature_cube:
+        for variable_name in self.feature_config.keys():
+            if variable_name not in forecast_df.columns:
                 msg = (
-                    f"Feature cube for {feature_name} not found in the provided "
-                    "feature cubes."
+                    f"Feature '{variable_name}' is not present in the "
+                    "forecast DataFrame."
                 )
                 raise ValueError(msg)
-            for feature in self.feature_config[feature_name]:
-                feature_values.append(
-                    prep_feature(forecast_cube, feature_cube[0], feature)
-                )
-        feature_values = np.array(feature_values).T
+            for feature_name in self.feature_config[variable_name]:
+                forecast_df = prep_feature(forecast_df, variable_name, feature_name)
+
+        forecast_df = sanitise_forecast_dataframe(forecast_df, self.feature_config)
+        feature_column_names = get_required_column_names(
+            forecast_df, self.feature_config
+        )
+        merge_columns = ["wmo_id", "time"]
+        combined_df = forecast_df.merge(
+            truth_df[merge_columns + ["ob_value"]], on=merge_columns, how="inner"
+        )
+        feature_values = np.array(combined_df[feature_column_names])
+        target_values = combined_df["ob_value"].values
 
-        truth_data_list = self._organise_truth_data(forecast_cube, truth_cube)
-        target_values = np.array(truth_data_list).flatten()
         # Fit the quantile regression model
         qrf_model = self.fit_qrf(feature_values, target_values)
 
@@ -547,31 +348,35 @@ class ApplyQuantileRegressionRandomForests(PostProcessingPlugin):
 
     def __init__(
         self,
+        target_name: str,
         feature_config: dict[str, list[str]],
         quantiles: list[np.float32],
         transformation: str = None,
         pre_transform_addition: np.float32 = 0,
     ) -> None:
         """Initialise the plugin.
+
         Args:
+            target_name (str):
+                Name of the target variable to be calibrated.
             feature_config (dict):
-                Feature configuration defining the features to be used for quantile regression.
-                The configuration is a dictionary of strings, where the keys are the names of
-                the input cube(s) supplied, and the values are a list. This list can contain both
-                computed features, such as the mean or standard deviation (std), or static
-                features, such as the altitude. The computed features will be computed using
-                the cube defined in the dictionary key. If the key is the feature itself e.g.
-                a distance to water cube, then the value should state "static". This will ensure
-                the cube's data is used as the feature.
-                The config will have the structure:
-                "DYNAMIC_VARIABLE_NAME": ["FEATURE1", "FEATURE2"] e.g:
+                Feature configuration defining the features to be used for quantile
+                regression. The configuration is a dictionary of strings, where the
+                keys are the names of the columns within the dataframe. Some
+                features may be used as initially provided within the dataframe,
+                whilst others may be computed from the data e.g. mean, std.
+                If the key is the feature itself e.g. distance to water, then the value
+                should state "static". In this case, the name of feature e.g.
+                'distance_to_water' is expected to be a column name in the input
+                dataframe. The config will have the structure:
+                "DYNAMIC_VARIABLE_CF_NAME": ["FEATURE1", "FEATURE2"] e.g.
                 {
                 "air_temperature": ["mean", "std", "altitude"],
                 "visibility_at_screen_level": ["mean", "std"]
                 "distance_to_water": ["static"],
                 }
             quantiles (float):
-                Quantiles used for prediction (values ranging from 0 to 1)
+                Quantiles used for prediction (values ranging from 0 to 1).
             transformation (str):
                 Transformation to be applied to the data before fitting.
             pre_transform_addition (float):
@@ -579,89 +384,78 @@ def __init__(
 
         Raises:
             ValueError: If the transformation is not one of the supported types.
-        """
 
+        """
+        self.target_name = target_name
         self.feature_config = feature_config
         self.quantiles = quantiles
         self.transformation = transformation
         _check_valid_transformation(self.transformation)
         self.pre_transform_addition = pre_transform_addition
 
-    def _reverse_transformation(self, forecast_cube: Cube):
+    def _reverse_transformation(self, forecast: np.ndarray) -> np.ndarray:
         """Reverse the transformation applied to the data prior to fitting the QRF.
-        The forecast cube provided is modified in place.
 
         Args:
-            forecast_cube: Forecast to be calibrated.
+            forecast: Calibrated forecast.
+        Returns:
+            forecast: Forecast with the transformation reversed.
         """
         if self.transformation:
             if self.transformation == "log":
-                forecast_cube.data = (
-                    np.exp(forecast_cube.data) - self.pre_transform_addition
-                )
+                forecast = np.exp(forecast) - self.pre_transform_addition
             elif self.transformation == "log10":
-                forecast_cube.data = (
-                    10 ** (forecast_cube.data) - self.pre_transform_addition
-                )
+                forecast = 10 ** (forecast) - self.pre_transform_addition
             elif self.transformation == "sqrt":
-                forecast_cube.data = forecast_cube.data**2 - self.pre_transform_addition
+                forecast = forecast**2 - self.pre_transform_addition
             elif self.transformation == "cbrt":
-                forecast_cube.data = forecast_cube.data**3 - self.pre_transform_addition
+                forecast = forecast**3 - self.pre_transform_addition
+        return forecast
 
     def process(
         self,
         qrf_model: RandomForestQuantileRegressor,
-        feature_cubes: CubeList,
-        template_forecast_cube: Cube,
-    ) -> Cube:
+        forecast_df: pd.DataFrame,
+    ) -> np.ndarray:
         """Apply a quantile regression random forests model.
 
         Args:
             qrf_model: A trained QRF model.
-            feature_cubes: CubeList of features. This should include the forecast to be
-                calibrated, if that has been used as a feature in the training, and
-                any features that can be provided as cubes.
-            template_forecast_cube: Template forecast cube that provides the required
-                metadata and shape for the output cube. The data from this cube will not
-                be used.
+            forecast_df: DataFrame containing the forecast information and features.
 
         Returns:
-            Calibrated forecast cube with the same metadata as the template
-            forecast cube.
+            Calibrated forecast as a numpy array.
 
         """
         feature_values = []
 
-        for feature_name in self.feature_config.keys():
-            feature_cube = feature_cubes.extract_cube(iris.Constraint(feature_name))
-
+        for variable_name in self.feature_config.keys():
             # Transform the feature cube data if a transformation is specified.
             if (
                 self.transformation
-                and set(["mean", "std"]).intersection(self.feature_config[feature_name])
-                and feature_cube.name() == template_forecast_cube.name()
+                and set(["mean", "std"]).intersection(
+                    self.feature_config[variable_name]
+                )
+                and self.target_name in forecast_df.columns
             ):
-                feature_cube.data = getattr(np, self.transformation)(
-                    feature_cube.data + self.pre_transform_addition
+                forecast_df[self.target_name] = getattr(np, self.transformation)(
+                    forecast_df[self.target_name] + self.pre_transform_addition
                 )
 
-            for feature in self.feature_config[feature_name]:
-                feature_values.append(
-                    prep_feature(template_forecast_cube, feature_cube, feature)
-                )
+            for feature_name in self.feature_config[variable_name]:
+                forecast_df = prep_feature(forecast_df, variable_name, feature_name)
+
+        forecast_df = sanitise_forecast_dataframe(forecast_df, self.feature_config)
+        feature_column_names = get_required_column_names(
+            forecast_df, self.feature_config
+        )
+
+        feature_values = np.array(forecast_df[feature_column_names])
 
-        feature_values = np.array(feature_values).T
         calibrated_forecast = qrf_model.predict(
             feature_values, quantiles=self.quantiles
         )
         calibrated_forecast = np.float32(calibrated_forecast)
 
-        calibrated_forecast_cube = template_forecast_cube.copy(
-            data=np.broadcast_to(calibrated_forecast.T, template_forecast_cube.shape)
-        )
-        self._reverse_transformation(calibrated_forecast_cube)
-        
-        # Rebadge the percentiles as realizations.
-        calibrated_forecast_cube = RebadgeRealizationsAsPercentiles()(calibrated_forecast_cube)
-
-        return calibrated_forecast_cube
+        calibrated_forecast = self._reverse_transformation(calibrated_forecast)
+        return calibrated_forecast
diff --git a/improver/cli/apply_quantile_regression_random_forest.py b/improver/cli/apply_quantile_regression_random_forest.py
index 92ad7a5fb6..c7a6fec8d4 100644
--- a/improver/cli/apply_quantile_regression_random_forest.py
+++ b/improver/cli/apply_quantile_regression_random_forest.py
@@ -13,7 +13,7 @@
 def process(
     *file_paths: cli.inputpath,
     feature_config: cli.inputjson,
-    target_cube_name: str,
+    target_cf_name: str,
     transformation: str = None,
     pre_transform_addition: float = 0,
 ):
@@ -32,23 +32,23 @@ def process(
             - The path to a NetCDF file containing the forecast to be calibrated.
             - Optionally, paths to NetCDF files containing additional predictors.
         feature_config (dict):
-            Feature configuration defining the features to be used for quantile regression.
-                The configuration is a dictionary of strings, where the keys are the names of
-                the input cube(s) supplied, and the values are a list. This list can contain both
-                computed features, such as the mean or standard deviation (std), or static
-                features, such as the altitude. The computed features will be computed using
-                the cube defined in the dictionary key. If the key is the feature itself e.g.
-                a distance to water cube, then the value should state "static". This will ensure
-                the cube's data is used as the feature.
-                The config will have the structure:
-                "DYNAMIC_VARIABLE_NAME": ["FEATURE1", "FEATURE2"] e.g:
-                {
-                "air_temperature": ["mean", "std", "altitude"],
-                "visibility_at_screen_level": ["mean", "std"]
-                "distance_to_water": ["static"],
-                }
-        target_cube_name (str):
-            A string containing the cube name of the forecast to be
+            Feature configuration defining the features to be used for quantile
+            regression. The configuration is a dictionary of strings, where the keys
+            are the names of the input cube(s) supplied, and the values are a list.
+            This list can contain both computed features, such as the mean or
+            standard deviation (std), or static features, such as the altitude. The
+            computed features will be computed using the cube defined in the
+            dictionary key. If the key is the feature itself e.g. a distance to water
+            cube, then the value should state "static". This will ensure the cube's
+            data is used as the feature. The config will have the structure:
+            "DYNAMIC_VARIABLE_CF_NAME": ["FEATURE1", "FEATURE2"] e.g.
+            {
+            "air_temperature": ["mean", "std", "altitude"],
+            "visibility_at_screen_level": ["mean", "std"]
+            "distance_to_water": ["static"],
+            }
+        target_cf_name (str):
+            A string containing the CF name of the forecast to be
             calibrated. This will be used to separate it from the rest of the
             feature cubes, if present.
         transformation (str):
@@ -66,7 +66,7 @@ def process(
 
     result = LoadAndApplyQRF(
         feature_config=feature_config,
-        target_cube_name=target_cube_name,
+        target_cf_name=target_cf_name,
         transformation=transformation,
         pre_transform_addition=pre_transform_addition,
     )(
diff --git a/improver/cli/train_quantile_regression_random_forest.py b/improver/cli/train_quantile_regression_random_forest.py
index 27ea3cf586..2c10a44421 100644
--- a/improver/cli/train_quantile_regression_random_forest.py
+++ b/improver/cli/train_quantile_regression_random_forest.py
@@ -13,12 +13,14 @@ def process(
     *file_paths: cli.inputpath,
     feature_config: cli.inputjson,
     target_diagnostic_name: str,
+    target_cf_name: str,
     forecast_periods: str,
     cycletime: str,
     training_length: int,
     experiment: str = None,
     n_estimators: int = 100,
     max_depth: int = None,
+    max_samples: float = None,
     random_state: int = None,
     transformation: str = None,
     pre_transform_addition: float = 0,
@@ -52,7 +54,7 @@ def process(
             If the key is the feature itself e.g. a distance to water cube, then the
             value should state "static". This will ensure the cube's data is used as
             the feature. The config will have the structure:
-            "DYNAMIC_VARIABLE_NAME": ["FEATURE1", "FEATURE2"] e.g:
+            "DYNAMIC_VARIABLE_CF_NAME": ["FEATURE1", "FEATURE2"] e.g:
             {
             "air_temperature": ["mean", "std", "altitude"],
             "visibility_at_screen_level": ["mean", "std"]
@@ -62,6 +64,8 @@ def process(
             A string containing the diagnostic name of the forecast to be
             calibrated. This will be used to filter the target forecast and truth
             dataframes.
+        target_cf_name (str):
+            A string containing the CF name of the forecast to be calibrated.
         forecast_periods (str):
             Range of forecast periods to be calibrated in hours in the form:
             "start:end:interval" e.g. "6:18:6" or a single forecast period e.g. "6".
@@ -77,6 +81,11 @@ def process(
             Number of trees in the forest.
         max_depth (int):
             Maximum depth of the tree.
+        max_samples (float):
+            If an int, then it is the number of samples to draw to train
+            each tree. If a float, then it is the fraction of samples to draw
+            to train each tree. If None, then each tree contains the same
+            total number of samples as originally provided.
         random_state (int):
             Random seed for reproducibility.
         transformation (str):
@@ -100,11 +109,13 @@ def process(
         experiment=experiment,
         feature_config=feature_config,
         target_diagnostic_name=target_diagnostic_name,
+        target_cf_name=target_cf_name,
         forecast_periods=forecast_periods,
         cycletime=cycletime,
         training_length=training_length,
         n_estimators=n_estimators,
         max_depth=max_depth,
+        max_samples=max_samples,
         random_state=random_state,
         transformation=transformation,
         pre_transform_addition=pre_transform_addition,
diff --git a/improver/utilities/compare.py b/improver/utilities/compare.py
index 75b943ae0c..598349c85e 100644
--- a/improver/utilities/compare.py
+++ b/improver/utilities/compare.py
@@ -406,10 +406,20 @@ def raise_reporter(message):
         reporter(str(exc))
         return
 
-    assert output.n_features_in_ == kgo.n_features_in_
-    assert output.n_outputs_ == kgo.n_outputs_
-    assert output.max_depth == kgo.max_depth
-    assert output.n_estimators == kgo.n_estimators
-    assert output.random_state == kgo.random_state
-    for output_estimator, kgo_estimator in zip(output.estimators_, kgo.estimators_):
-        assert (output_estimator.tree_.value == kgo_estimator.tree_.value).all()
+    difference_found = False
+    try:
+        np.testing.assert_equal(output.n_features_in_, kgo.n_features_in_)
+        np.testing.assert_equal(output.n_outputs_, kgo.n_outputs_)
+        np.testing.assert_equal(output.max_depth, kgo.max_depth)
+        np.testing.assert_equal(output.n_estimators, kgo.n_estimators)
+        np.testing.assert_equal(output.random_state, kgo.random_state)
+        for output_estimator, kgo_estimator in zip(output.estimators_, kgo.estimators_):
+            np.testing.assert_allclose(
+                output_estimator.tree_.value, kgo_estimator.tree_.value
+            )
+    except (AssertionError, ValueError):
+        difference_found = True
+    # call the reporter function outside the except block to avoid nested
+    # exceptions if the reporter function is raising an exception
+    if difference_found:
+        reporter("different pickled forest")
diff --git a/improver_tests/acceptance/SHA256SUMS b/improver_tests/acceptance/SHA256SUMS
index ed087c4f5d..4757fb1a02 100644
--- a/improver_tests/acceptance/SHA256SUMS
+++ b/improver_tests/acceptance/SHA256SUMS
@@ -87,9 +87,9 @@ ed8ab78a9f55b54bf0a49f191d3eb33daae30e0d05b9051275a70c0e697aac71  ./apply-night-
 d2f7d8389b33cde359dd253def2aaf31afbc27557f389bf0f744a28b24e145dd  ./apply-quantile-regression-random-forest/config.json
 8db1b35bde734c16340a6e42454b9ac146ea32f7c012c17c5e64815069af41bf  ./apply-quantile-regression-random-forest/input_forecast.nc
 793f17191d2421ead8d546319c7da34143ef01e2db88ca6cc5faf46e92b4b1de  ./apply-quantile-regression-random-forest/with_transformation_input.pickle
-e5221acb444360ee802dc48e4da1ab3a874a4179ceff1a611138eb2af3f6591c  ./apply-quantile-regression-random-forest/with_transformation_kgo.nc
+d8acadfe782e2c2562dfcc31ada8ad381b320fbdf3854d6b361c38308e6b8226  ./apply-quantile-regression-random-forest/with_transformation_kgo.nc
 31eab8a83a75eb8cee4c0140749a4d4e9d9b98292c28b3c53ee27756c1f3a265  ./apply-quantile-regression-random-forest/without_transformation_input.pickle
-12ae824d2c23362eae37f41a3951c6d3fc48150e45a5998a9eafe467e34f0c1d  ./apply-quantile-regression-random-forest/without_transformation_kgo.nc
+65397673caf0f1205d18a85c08ce81bd3543a71c4101f40f4869a93ddd0cf0f8  ./apply-quantile-regression-random-forest/without_transformation_kgo.nc
 74a40c508b68294652f5dde52d8fe4b93ea002bd4b72ca8d6cc11a83309f28f1  ./apply-rainforests-calibration/basic/kgo.nc
 b4678f178e0e1df0c242c61b269b71dd9a3d7aca64663e2e81e5bfeb71fe97da  ./apply-rainforests-calibration/features/20200802T0000Z-PT0000H00M-clearsky_solar_radiation-PT24H.nc
 dce2d4a28cec494c94b154dc3a9203b46b7cef2ad54d04049663166f4b42ecf5  ./apply-rainforests-calibration/features/20200802T0000Z-PT0024H00M-cape-PT24H.nc
@@ -1003,8 +1003,8 @@ e3d2315b24bf769bcfb137033950d7100c8795878635e0ec9491413a1349904a  ./train-quanti
 1aec562fc0c39c7695521482b889d335a3becb10a76b7eae965447cbd9faf051  ./train-quantile-regression-random-forest/spot_observation_tables/20250804/diagnostic=temperature_at_screen_level/0600Z_0.parquet
 e81c04f9f2d8dc14b8d8cd9bc9e8827ebd0da5ba016db48cf60ab738d327f968  ./train-quantile-regression-random-forest/spot_observation_tables/20250804/diagnostic=temperature_at_screen_level/1200Z_0.parquet
 99e6656982672ea2d4cb9fa1dfb5731408dd41b2c9052b7e0df8d8004864241e  ./train-quantile-regression-random-forest/spot_observation_tables/20250804/diagnostic=temperature_at_screen_level/1800Z_0.parquet
-793f17191d2421ead8d546319c7da34143ef01e2db88ca6cc5faf46e92b4b1de  ./train-quantile-regression-random-forest/with_transformation_kgo.pickle
-31eab8a83a75eb8cee4c0140749a4d4e9d9b98292c28b3c53ee27756c1f3a265  ./train-quantile-regression-random-forest/without_transformation_kgo.pickle
+d634d639dda9a4a2789eed69bfdec59954ce63f51734a71f1d9a20c3958d2786  ./train-quantile-regression-random-forest/with_transformation_kgo.pickle
+e3f9fcd7544028436e34415db38d3c5841a6e22abbc3fc635e6ff87236764fcc  ./train-quantile-regression-random-forest/without_transformation_kgo.pickle
 e2cfb5e19ef5ebcfd73dc1c8504b66e9a55ad76b7b18cb79318659224079f234  ./uv-index/basic/20181210T0600Z-PT0000H00M-radiation_flux_in_uv_downward_at_surface.nc
 e48c3b07bd14214a1a10c0bbf1e0acdf54cfedf93b2c6d766f594ce83a45c2b8  ./uv-index/basic/kgo.nc
 2226a5e95eb29664e21f8c0658101a53d65945a1450a60a19955563a2693d22d  ./vertical-updraught/cape.nc
diff --git a/improver_tests/acceptance/acceptance.py b/improver_tests/acceptance/acceptance.py
index 0cf644c71f..3616f225be 100644
--- a/improver_tests/acceptance/acceptance.py
+++ b/improver_tests/acceptance/acceptance.py
@@ -320,6 +320,8 @@ def compare(
         rtol (float): Relative tolerance
         exclude_vars (Iterable[str]): Variables to exclude from comparison
         exclude_attributes (Iterable[str]): Attributes to exclude from comparison
+        file_type (str): Name of file type to compare, either "netCDF" or
+            "pickled_forest".
 
     Returns:
         None
diff --git a/improver_tests/acceptance/test_apply_quantile_regression_random_forest.py b/improver_tests/acceptance/test_apply_quantile_regression_random_forest.py
index 68e8d9cb40..5cd27f5f17 100644
--- a/improver_tests/acceptance/test_apply_quantile_regression_random_forest.py
+++ b/improver_tests/acceptance/test_apply_quantile_regression_random_forest.py
@@ -30,7 +30,7 @@ def test_basic(tmp_path, transformation):
         qrf_path,
         "--feature-config",
         config_path,
-        "--target-cube-name",
+        "--target-cf-name",
         "air_temperature",
         "--output",
         output_path,
diff --git a/improver_tests/acceptance/test_train_quantile_regression_random_forest.py b/improver_tests/acceptance/test_train_quantile_regression_random_forest.py
index 009852f4d3..64feccd97e 100644
--- a/improver_tests/acceptance/test_train_quantile_regression_random_forest.py
+++ b/improver_tests/acceptance/test_train_quantile_regression_random_forest.py
@@ -43,6 +43,8 @@ def test_basic(
         config_path,
         "--target-diagnostic-name",
         "temperature_at_screen_level",
+        "--target-cf-name",
+        "air_temperature",
         "--forecast-periods",
         "6:18:6",
         "--cycletime",
@@ -59,7 +61,7 @@ def test_basic(
         "42",
         "--compression",
         "5",
-        "--model-output",
+        "--output",
         output_path,
     ]
     if transformation == "with_transformation":
diff --git a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_apply_quantile_regression_random_forest.py b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_apply_quantile_regression_random_forest.py
index 906cd41332..6d237f560b 100644
--- a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_apply_quantile_regression_random_forest.py
+++ b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_apply_quantile_regression_random_forest.py
@@ -6,6 +6,7 @@
 
 import numpy as np
 import pytest
+from iris.coords import AuxCoord
 from iris.cube import Cube
 
 from improver.calibration.load_and_apply_quantile_regression_random_forest import (
@@ -16,27 +17,45 @@
 )
 from improver.utilities.save import save_netcdf
 from improver_tests.calibration.quantile_regression_random_forests_calibration.test_quantile_regression_random_forest import (
-    _add_day_of_training_period,
     _create_ancil_file,
     _create_forecasts,
     _run_train_qrf,
 )
 
 
+def _add_day_of_training_period_to_cube(cube, day_of_training_period, secondary_coord):
+    """Add day of training period coordinate to the cube.
+    Args:
+        cube: Cube to which the day of training period coordinate will be added.
+        day_of_training_period: Day of training period to be added.
+        secondary_coord: Coordinate to associate the day of training period with.
+    Returns:
+        Cube with the day of training period coordinate added.
+    """
+    dims = cube.coord_dims(secondary_coord)
+    day_of_training_period_coord = AuxCoord(
+        np.array(day_of_training_period, dtype=np.int32),
+        long_name="day_of_training_period",
+        units="1",
+    )
+    cube.add_aux_coord(day_of_training_period_coord, data_dims=dims)
+    return cube
+
+
 @pytest.mark.parametrize("percentile_input", [True, False])
 @pytest.mark.parametrize(
     "n_estimators,max_depth,random_state,compression,transformation,pre_transform_addition,extra_kwargs,include_static,quantiles,expected",
     [
-        (2, 2, 55, 5, None, 0, {}, False, [0.5], [5.15, 5.65]),  # noqa Basic test case
-        (100, 2, 55, 5, None, 0, {}, False, [1 / 3, 2 / 3], [[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
+        (2, 2, 55, 5, None, 0, {}, False, [0.5], [4.1, 5.65]),  # noqa Basic test case
+        (100, 2, 55, 5, None, 0, {}, False, [1 / 3, 2 / 3], [[4.1, 5.1], [5.1, 5.1]]),  # noqa Multiple quantiles
+        (1, 1, 55, 5, None, 0, {}, False, [0.5], [4.1, 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], [5.11, 5.64]),  # Log transformation
-        (2, 2, 55, 5, "log10", 10, {}, False, [0.5], [5.11, 5.64]),  # noqa Log10 transformation
-        (2, 2, 55, 5, "sqrt", 10, {}, False, [0.5], [5.11, 5.64]),  # noqa Square root transformation
-        (2, 2, 55, 5, "cbrt", 10, {}, False, [0.5], [5.13, 5.64]),  # noqa 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]),  # noqa Include an additional static feature
+        (2, 2, 55, 5, "log", 10, {}, False, [0.5], [4.1, 5.64]),  # Log transformation
+        (2, 2, 55, 5, "log10", 10, {}, False, [0.5], [4.1, 5.64]),  # noqa Log10 transformation
+        (2, 2, 55, 5, "sqrt", 10, {}, False, [0.5], [4.1, 5.64]),  # noqa Square root transformation
+        (2, 2, 55, 5, "cbrt", 10, {}, False, [0.5], [4.1, 5.64]),  # noqa Cube root transformation
+        (2, 2, 55, 5, None, 0, {"max_samples_leaf": 0.5}, False, [0.5], [4.1, 6.2]),  # noqa # Different criterion
+        (2, 5, 55, 5, None, 0, {}, True, [0.5], [4.1, 5.65]),  # noqa Include an additional static feature
     ],
 )
 def test_load_and_apply_qrf(
@@ -75,20 +94,20 @@ def test_load_and_apply_qrf(
             "20170101T1200Z",
             "20170102T1200Z",
         ],
-        day_of_training_period=[0, 1],
         realization_data=[2, 6, 10],
-        truth_data=[[4.2, 6.2], [4.1, 5.1]],
+        truth_data=[4.2, 6.2, 4.1, 5.1],
     )
 
     frt = "20170103T0000Z"
     vt = "20170103T1200Z"
-    day_of_training_period = 2
     data = np.arange(6, (len(quantiles) * 6) + 1, 6)
+    day_of_training_period = 2
+    forecast_cube = _create_forecasts(frt, vt, data, return_cube=True)
 
-    forecast_cube = _create_forecasts(frt, vt, data)
-    forecast_cube = _add_day_of_training_period(
+    forecast_cube = _add_day_of_training_period_to_cube(
         forecast_cube, day_of_training_period, "forecast_reference_time"
     )
+
     if percentile_input:
         forecast_cube = RebadgeRealizationsAsPercentiles()(forecast_cube)
 
@@ -100,7 +119,7 @@ def test_load_and_apply_qrf(
     file_paths = [model_output, forecast_filepath]
 
     if include_static:
-        ancil_cube = _create_ancil_file()
+        ancil_cube = _create_ancil_file(return_cube=True)
         ancil_filepath = features_dir / "ancil.nc"
         save_netcdf(ancil_cube, ancil_filepath)
         file_paths.append(str(ancil_filepath))
@@ -114,7 +133,7 @@ def test_load_and_apply_qrf(
     result = plugin.process(file_paths=file_paths)
     assert isinstance(result, Cube)
 
-    assert result.data.shape == (1, 1, len(quantiles), 2)
+    assert result.data.shape == (len(quantiles), 2)
     assert np.allclose(result.data, expected, rtol=1e-2)
 
     # Check that the metadata is as expected
@@ -178,22 +197,20 @@ def test_exceptions(
             "20170101T1200Z",
             "20170102T1200Z",
         ],
-        day_of_training_period=[0, 1],
         realization_data=[2, 6, 10],
-        truth_data=[[4.2, 6.2], [4.1, 5.1]],
+        truth_data=[4.2, 6.2, 4.1, 5.1],
     )
 
     frt = "20170103T0000Z"
     vt = "20170103T1200Z"
-    day_of_training_period = 2
     data = np.arange(6, (len(quantiles) * 6) + 1, 6)
-
-    forecast_cube = _create_forecasts(frt, vt, data)
-    forecast_cube = _add_day_of_training_period(
+    day_of_training_period = 2
+    forecast_cube = _create_forecasts(frt, vt, data, return_cube=True)
+    forecast_cube = _add_day_of_training_period_to_cube(
         forecast_cube, day_of_training_period, "forecast_reference_time"
     )
 
-    ancil_cube = _create_ancil_file()
+    ancil_cube = _create_ancil_file(return_cube=True)
     ancil_filepath = tmp_path / "ancil.nc"
     save_netcdf(ancil_cube, ancil_filepath)
 
diff --git a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_train_quantile_regression_random_forest.py b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_train_quantile_regression_random_forest.py
index 51925d984a..0891ec6f82 100644
--- a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_train_quantile_regression_random_forest.py
+++ b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_load_and_train_quantile_regression_random_forest.py
@@ -60,9 +60,6 @@ def _create_multi_site_forecast_parquet_file(tmp_path, representation="percentil
     wind_speed_df = pd.DataFrame(wind_speed_dict)
     joined_df = pd.concat([data_df, wind_speed_df], ignore_index=True)
 
-    joined_df["forecast_period"] = joined_df["forecast_period"].astype(
-        "timedelta64[ms]"
-    )
     output_dir = tmp_path / "forecast_parquet_files"
     output_dir.mkdir(parents=True, exist_ok=True)
     output_path = str(output_dir / "forecast.parquet")
@@ -101,9 +98,6 @@ def _create_multi_percentile_forecast_parquet_file(tmp_path, representation=None
     data_df = pd.DataFrame(data_dict)
     wind_speed_df = pd.DataFrame(wind_speed_dict)
     joined_df = pd.concat([data_df, wind_speed_df], ignore_index=True)
-    joined_df["forecast_period"] = joined_df["forecast_period"].astype(
-        "timedelta64[ms]"
-    )
 
     output_dir = tmp_path / "forecast_parquet_files"
     output_dir.mkdir(parents=True, exist_ok=True)
@@ -152,9 +146,6 @@ def _create_multi_forecast_period_forecast_parquet_file(tmp_path, representation
     data_df = pd.DataFrame(data_dict)
     wind_speed_df = pd.DataFrame(wind_speed_dict)
     joined_df = pd.concat([data_df, wind_speed_df], ignore_index=True)
-    joined_df["forecast_period"] = joined_df["forecast_period"].astype(
-        "timedelta64[ms]"
-    )
 
     output_dir = tmp_path / "forecast_parquet_files"
     output_dir.mkdir(parents=True, exist_ok=True)
@@ -353,7 +344,7 @@ def _create_ancil_file(tmp_path, wmo_ids):
             False,
             False,
             "realization",  # Provide realization input
-            5.6,
+            5.62,
         ),
         (
             _create_multi_forecast_period_forecast_parquet_file,
@@ -391,7 +382,7 @@ def test_load_and_train_qrf(
     model_output = str(model_output_dir / "qrf_model.pkl")
 
     if include_static:
-        ancil_path = _create_ancil_file(tmp_path, list(set(wmo_ids)))
+        ancil_path = _create_ancil_file(tmp_path, sorted(list(set(wmo_ids))))
         file_paths.append(ancil_path)
         feature_config["distance_to_water"] = ["static"]
 
@@ -403,6 +394,7 @@ def test_load_and_train_qrf(
         experiment="latestblend",
         feature_config=feature_config,
         target_diagnostic_name="temperature_at_screen_level",
+        target_cf_name="air_temperature",
         forecast_periods=forecast_periods,
         cycletime="20170103T0000Z",
         training_length=2,
@@ -423,7 +415,7 @@ def test_load_and_train_qrf(
     if remove_target:
         current_forecast = []
     else:
-        current_forecast = [279, 3, 55]
+        current_forecast = [5.64, 3, 55]
 
     if include_static:
         current_forecast.append(2.5)
@@ -459,6 +451,7 @@ def test_load_and_train_qrf_no_paths(tmp_path, make_files):
         experiment="latestblend",
         feature_config=feature_config,
         target_diagnostic_name="temperature_at_screen_level",
+        target_cf_name="air_temperature",
         forecast_periods="6:12:6",
         cycletime="20170102T0000Z",
         training_length=2,
@@ -503,6 +496,7 @@ def test_load_and_train_qrf_mismatches(tmp_path, cycletime, forecast_periods):
         experiment="latestblend",
         feature_config=feature_config,
         target_diagnostic_name="temperature_at_screen_level",
+        target_cf_name="air_temperature",
         forecast_periods=forecast_periods,
         cycletime=cycletime,
         training_length=2,
@@ -586,6 +580,7 @@ def test_exceptions(
         experiment="latestblend",
         feature_config=feature_config,
         target_diagnostic_name="temperature_at_screen_level",
+        target_cf_name="air_temperature",
         forecast_periods=forecast_periods,
         cycletime="20170103T0000Z",
         training_length=2,
diff --git a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py
index b7ff4e6825..d8fc7af795 100644
--- a/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py
+++ b/improver_tests/calibration/quantile_regression_random_forests_calibration/test_quantile_regression_random_forest.py
@@ -7,33 +7,40 @@
 from datetime import datetime as dt
 
 import iris
+import itertools
 import joblib
 import numpy as np
 import pandas as pd
 import pytest
-from iris.coords import AuxCoord
-from iris.cube import Cube, CubeList
+from iris.cube import Cube
+from iris.pandas import as_data_frame
+from pandas.testing import assert_frame_equal
 
 from improver.calibration.quantile_regression_random_forest import (
     ApplyQuantileRegressionRandomForests,
     TrainQuantileRegressionRandomForests,
+    _check_valid_transformation,
+    get_required_column_names,
     prep_feature,
+    sanitise_forecast_dataframe,
 )
 from improver.metadata.constants.time_types import DT_FORMAT
 from improver.synthetic_data.set_up_test_cubes import set_up_spot_variable_cube
-from improver.utilities.temporal import datetime_to_iris_time
 
 ALTITUDE = [10, 20]
 LATITUDE = [50, 60]
 LONGITUDE = [0, 10]
 WMO_ID = ["00001", "00002"]
 
+iris.FUTURE.pandas_ndim = True
+
 
 def _create_forecasts(
     forecast_reference_time: str,
     validity_time: str,
     data: list[int],
-) -> Cube:
+    return_cube: bool = False,
+) -> Cube | pd.DataFrame:
     """Create site forecast cube with realizations.
 
     Args:
@@ -45,8 +52,7 @@ def _create_forecasts(
     Returns:
         Forecast cube containing three percentiles and two sites.
     """
-
-    data = np.array(data, dtype=np.float32).repeat(2).reshape(len(data), 2)
+    data = np.array([data, data + 2], dtype=np.float32).T
     cube = set_up_spot_variable_cube(
         data,
         realizations=range(len(data)),
@@ -59,36 +65,48 @@ def _create_forecasts(
         time=dt.strptime(validity_time, DT_FORMAT),
         frt=dt.strptime(forecast_reference_time, DT_FORMAT),
     )
-    cube.remove_coord("time")
-    return cube
+    if return_cube:
+        return cube
+    df = as_data_frame(
+        cube,
+        add_aux_coords=[
+            "altitude",
+            "latitude",
+            "longitude",
+            "wmo_id",
+            "forecast_period",
+            "forecast_reference_time",
+            "time",
+        ],
+    ).reset_index()
+    df["time"] = df["time"].apply(lambda x: x._to_real_datetime())
+    df["forecast_reference_time"] = df["forecast_reference_time"].apply(
+        lambda x: x._to_real_datetime()
+    )
+    return df
 
 
-def _add_day_of_training_period(cube, day_of_training_period, secondary_coord):
-    """Add day of training period coordinate to the cube.
+def _add_day_of_training_period(df):
+    """Add day of training period coordinate to the dataframe.
 
     Args:
-        cube: Cube to which the day of training period coordinate will be added.
+        df: DataFrame to which the day of training period coordinate will be added.
         day_of_training_period: Day of training period to be added.
-        secondary_coord: Coordinate to associate the day of training period with.
 
     Returns:
         Cube with the day of training period coordinate added.
     """
-    dims = cube.coord_dims(secondary_coord)
-    day_of_training_period_coord = AuxCoord(
-        np.array(day_of_training_period, dtype=np.int32),
-        long_name="day_of_training_period",
-        units="1",
+    df["day_of_training_period"] = df["time"].dt.dayofyear - np.min(
+        df["time"].dt.dayofyear
     )
-    cube.add_aux_coord(day_of_training_period_coord, data_dims=dims)
-    return cube
+    return df
 
 
-def _create_ancil_file():
+def _create_ancil_file(return_cube=False):
     """Create an ancillary file for testing.
 
     Returns:
-        An ancillary cube with a single value.
+        An ancillary DataFrame without temporal columns.
     """
     data = np.array([2, 3], dtype=np.float32)
     template_cube = set_up_spot_variable_cube(
@@ -103,7 +121,17 @@ def _create_ancil_file():
     cube = template_cube.copy()
     for coord in ["time", "forecast_reference_time", "forecast_period"]:
         cube.remove_coord(coord)
-    return cube
+    if return_cube:
+        return cube
+    return as_data_frame(
+        cube,
+        add_aux_coords=[
+            "altitude",
+            "latitude",
+            "longitude",
+            "wmo_id",
+        ],
+    ).reset_index()
 
 
 def _run_train_qrf(
@@ -129,65 +157,49 @@ def _run_train_qrf(
         "20170102T1200Z",
         "20170102T1800Z",
     ],
-    day_of_training_period=[0, 1],
     realization_data=[2, 6, 10],
-    truth_data=[[4.2, 3.8], [5.8, 6], [7, 7.3], [9.1, 9.5]],
+    truth_data=[4.2, 3.8, 5.8, 6, 7, 7.3, 9.1, 9.5],
 ):
     realization_data = np.array(realization_data, dtype=np.float32)
-    forecast_cubes = CubeList()
+    forecast_dfs = []
     for index, (frt, vt) in enumerate(zip(forecast_reference_times, validity_times)):
-        forecast_cube = _create_forecasts(frt, vt, realization_data + index)
-        forecast_cubes.append(forecast_cube)
-    forecast_cube = forecast_cubes.merge_cube()
-
-    forecast_cube = _add_day_of_training_period(
-        forecast_cube, day_of_training_period, "forecast_reference_time"
+        forecast_df = _create_forecasts(frt, vt, realization_data + index)
+        forecast_dfs.append(forecast_df)
+    forecast_df = pd.concat(forecast_dfs)
+    forecast_df = _add_day_of_training_period(forecast_df)
+
+    truth_df = forecast_df.copy()
+    truth_df.drop(
+        columns=[
+            "forecast_period",
+            "forecast_reference_time",
+            "day_of_training_period",
+            "realization",
+            "wind_speed_at_10m",
+        ],
+        inplace=True,
     )
+    truth_df.drop_duplicates(inplace=True)
+    truth_df["ob_value"] = np.array(truth_data, dtype=np.float32)
 
-    truth_cube = next(forecast_cube.slices_over("realization")).copy()
-    for coord in [
-        "day_of_training_period",
-        "realization",
-    ]:
-        truth_cube.remove_coord(coord)
-    truth_slices = CubeList()
-    for truth_slice in truth_cube.slices_over(
-        ["forecast_period", "forecast_reference_time"]
-    ):
-        frt_coord = truth_slice.coord("forecast_reference_time")
-        fp_coord = truth_slice.coord("forecast_period")
-        time_point = datetime_to_iris_time(
-            frt_coord.cell(0).point._to_real_datetime()
-            + pd.to_timedelta(fp_coord.points[0], unit=str(fp_coord.units))
-        )
-        time_coord = iris.coords.AuxCoord(
-            time_point, standard_name="time", units=frt_coord.units
-        )
-        truth_slice.add_aux_coord(time_coord)
-        truth_slice.remove_coord("forecast_period")
-        truth_slice.remove_coord("forecast_reference_time")
-        truth_slices.append(truth_slice)
-
-    truth_cube = truth_slices.merge_cube()
-    truth_cube.data = np.array(truth_data).astype(np.float32)
-
-    feature_cubes = CubeList([forecast_cube])
     if include_static:
-        ancil_cube = _create_ancil_file()
-        feature_cubes.append(ancil_cube)
-        feature_config[ancil_cube.name()] = ["static"]
+        ancil_df = _create_ancil_file()
+        forecast_df = forecast_df.merge(
+            ancil_df[["wmo_id", "distance_to_water"]], on=["wmo_id"], how="left"
+        )
+        feature_config["distance_to_water"] = ["static"]
     if "air_temperature" in feature_config.keys():
-        dynamic_cube = forecast_cube.copy(data=forecast_cube.data + 10)
-        dynamic_cube.rename("air_temperature")
-        dynamic_cube.units = "Celsius"
-        feature_cubes.append(dynamic_cube)
+        forecast_df["air_temperature"] = np.array(
+            forecast_df["wind_speed_at_10m"] + 10, dtype=np.float32
+        )
 
     model_output_dir = tmp_path / "train_qrf"
     model_output_dir.mkdir(parents=True)
     model_output = str(model_output_dir / "qrf_model.pkl")
 
     plugin = TrainQuantileRegressionRandomForests(
-        feature_config,
+        target_name="wind_speed_at_10m",
+        feature_config=feature_config,
         n_estimators=n_estimators,
         max_depth=max_depth,
         random_state=random_state,
@@ -197,172 +209,129 @@ def _run_train_qrf(
         model_output=model_output,
         **extra_kwargs,
     )
-    plugin.process(forecast_cube, truth_cube, feature_cubes=feature_cubes)
+    plugin.process(forecast_df, truth_df)
     return model_output
 
 
 @pytest.mark.parametrize(
     "feature,expected,expected_dtype",
     [
-        ("mean", np.array([6, 6], dtype=np.float32), np.float32),
-        ("std", np.array([4, 4], dtype=np.float32), np.float32),
-        ("latitude", np.array([50, 60], dtype=np.float32), np.float32),
-        ("longitude", np.array([0, 10], dtype=np.float32), np.float32),
-        ("altitude", np.array([10, 20], dtype=np.float32), np.float32),
+        ("mean", np.tile(np.array([6, 8], dtype=np.float32), 3), np.float32),
+        ("std", np.repeat(4, 6).astype(np.float32), np.float32),
+        ("latitude", np.tile(np.array([50, 60], dtype=np.float32), 3), np.float32),
+        ("longitude", np.tile(np.array([0, 10], dtype=np.float32), 3), np.float32),
+        ("altitude", np.tile(np.array([10, 20], dtype=np.float32), 3), np.float32),
         (
             "day_of_year",
-            np.array([1, 1], dtype=np.float32),
-            np.float32,
+            np.tile(np.array([1, 1], dtype=np.int32), 3),
+            np.int32,
         ),
         (
             "day_of_year_sin",
-            np.array([0.01716633, 0.01716633], dtype=np.float32),
+            np.tile(np.array([0.01716633, 0.01716633], dtype=np.float32), 3),
             np.float32,
         ),
         (
             "day_of_year_cos",
-            np.array([0.99985266, 0.99985266], dtype=np.float32),
+            np.tile(np.array([0.99985266, 0.99985266], dtype=np.float32), 3),
             np.float32,
         ),
-        ("hour_of_day", np.array([12, 12], dtype=np.float32), np.float32),
-        ("hour_of_day_sin", np.array([0, 0], dtype=np.float32), np.float32),
-        ("hour_of_day_cos", np.array([-1, -1], dtype=np.float32), np.float32),
-        ("forecast_period", np.array([43200, 43200], dtype=np.int32), np.int32),
-        ("day_of_training_period", np.array([0, 0], dtype=np.int32), np.int32),
-        ("static", np.array([2, 3], dtype=np.float32), np.float32),
+        ("hour_of_day", np.repeat(12, 6).astype(np.int32), np.int32),
+        ("hour_of_day_sin", np.repeat(0, 6).astype(np.float32), np.float32),
+        ("hour_of_day_cos", np.repeat(-1, 6).astype(np.float32), np.float32),
+        ("forecast_period", np.repeat(43200, 6).astype(np.int32), np.int32),
+        ("day_of_training_period", np.repeat(0, 6).astype(np.int32), np.int32),
+        ("static", np.tile(np.array([2, 3], dtype=np.float32), 3), np.float32),
     ],
 )
 def test_prep_feature_single_time(feature, expected, expected_dtype):
     """Test the prep_feature function for a single time."""
+    feature_name = "wind_speed_at_10m"
+
     forecast_reference_time = "20170101T0000Z"
     validity_time = "20170101T1200Z"
-    data = [2, 6, 10]
-    day_of_training_period = [0]
-    forecast_cube = _create_forecasts(forecast_reference_time, validity_time, data)
-    forecast_cube = _add_day_of_training_period(
-        forecast_cube, day_of_training_period, "forecast_reference_time"
-    )
+    data = np.array([2, 6, 10])
+    forecast_df = _create_forecasts(forecast_reference_time, validity_time, data)
+    forecast_df = _add_day_of_training_period(forecast_df)
 
     if feature == "static":
-        feature_cube = _create_ancil_file()
+        feature_df = _create_ancil_file()
+        forecast_df = forecast_df.merge(
+            feature_df[["wmo_id", "distance_to_water"]], on=["wmo_id"], how="left"
+        )
+    else:
+        forecast_df = forecast_df.copy()
+
+    result = prep_feature(forecast_df, feature_name, feature)
+
+    if feature in ["mean", "std"]:
+        assert result.shape == (6, 12)
+        feature_name_modified = f"{feature_name}_{feature}"
+    elif feature in [
+        "day_of_year",
+        "day_of_year_sin",
+        "day_of_year_cos",
+        "hour_of_day",
+        "hour_of_day_sin",
+        "hour_of_day_cos",
+    ]:
+        assert result.shape == (6, 12)
+        feature_name_modified = feature
+    elif feature in ["static"]:
+        assert result.shape == (6, 12)
+        feature_name_modified = "distance_to_water"
     else:
-        feature_cube = forecast_cube.copy()
+        assert result.shape == (6, 11)
+        feature_name_modified = feature
 
-    result = prep_feature(forecast_cube, feature_cube, feature)
-    assert result.shape == (2,)
-    assert result.dtype == expected_dtype
-    np.testing.assert_allclose(result, expected, atol=1e-6)
+    assert result[feature_name_modified].dtype == expected_dtype
+    np.testing.assert_allclose(result[feature_name_modified], expected, atol=1e-6)
 
 
 @pytest.mark.parametrize(
     "feature,expected,expected_dtype",
     [
-        ("mean", np.repeat(6, 8).astype(np.float32), np.float32),
-        ("std", np.repeat(4, 8).astype(np.float32), np.float32),
-        ("latitude", np.tile([50, 60], 4).astype(np.float32), np.float32),
-        ("longitude", np.tile([0, 10], 4).astype(np.float32), np.float32),
-        ("altitude", np.tile([10, 20], 4).astype(np.float32), np.float32),
+        ("mean", np.tile([6, 8], 18).astype(np.float32), np.float32),
+        ("std", np.repeat(4, 36).astype(np.float32), np.float32),
+        ("latitude", np.tile([50, 60], 18).astype(np.float32), np.float32),
+        ("longitude", np.tile([0, 10], 18).astype(np.float32), np.float32),
+        ("altitude", np.tile([10, 20], 18).astype(np.float32), np.float32),
         (
             "day_of_year",
-            np.repeat([1, 2], 4).astype(np.float32),
-            np.float32,
+            np.repeat([1, 2, 3], 12).astype(np.int32),
+            np.int32,
         ),
         (
             "day_of_year_sin",
-            np.repeat([0.017166, 0.034328], 4).astype(np.float32),
+            np.repeat([0.017166, 0.034328, 0.051479], 12).astype(np.float32),
             np.float32,
         ),
         (
             "day_of_year_cos",
-            np.repeat([0.999853, 0.999411], 4).astype(np.float32),
+            np.repeat([0.999853, 0.999411, 0.998674], 12).astype(np.float32),
             np.float32,
         ),
-        ("hour_of_day", np.repeat([12, 12], 4).astype(np.float32), np.float32),
-        ("hour_of_day_sin", np.repeat([0, 0], 4).astype(np.float32), np.float32),
-        ("hour_of_day_cos", np.repeat([-1, -1], 4).astype(np.float32), np.float32),
+        ("hour_of_day", np.tile(np.repeat([6, 12], 6), 3).astype(np.int32), np.int32),
         (
-            "forecast_period",
-            np.tile([43200, 43200, 21600, 21600], 2).astype(np.int32),
-            np.int32,
-        ),
-        ("day_of_training_period", np.repeat([0, 1], 4).astype(np.int32), np.int32),
-        ("static", np.tile([2, 3], 4).astype(np.float32), np.float32),
-    ],
-)
-def test_prep_feature_1d_time_dimension(feature, expected, expected_dtype):
-    """Test the prep_feature function for multiple times."""
-    forecast_reference_times = [
-        "20170101T0000Z",
-        "20170101T0600Z",
-        "20170102T0000Z",
-        "20170102T0600Z",
-    ]
-    validity_times = [
-        "20170101T1200Z",
-        "20170101T1200Z",
-        "20170102T1200Z",
-        "20170102T1200Z",
-    ]
-
-    data = [2, 6, 10]
-
-    forecast_cubes = CubeList()
-    for frt, vt in zip(forecast_reference_times, validity_times):
-        forecast_cubes.append(_create_forecasts(frt, vt, data))
-    forecast_cube = forecast_cubes.merge_cube()
-    day_of_training_period = [0, 0, 1, 1]
-    forecast_cube = _add_day_of_training_period(
-        forecast_cube, day_of_training_period, "forecast_reference_time"
-    )
-
-    if feature == "static":
-        feature_cube = _create_ancil_file()
-    else:
-        feature_cube = forecast_cube.copy()
-
-    result = prep_feature(forecast_cube, feature_cube, feature)
-
-    assert result.shape == (8,)
-    assert result.dtype == expected_dtype
-    np.testing.assert_allclose(result, expected, atol=1e-6)
-
-
-@pytest.mark.parametrize(
-    "feature,expected,expected_dtype",
-    [
-        ("mean", np.repeat(6, 12).astype(np.float32), np.float32),
-        ("std", np.repeat(4, 12).astype(np.float32), np.float32),
-        ("latitude", np.tile([50, 60], 6).astype(np.float32), np.float32),
-        ("longitude", np.tile([0, 10], 6).astype(np.float32), np.float32),
-        ("altitude", np.tile([10, 20], 6).astype(np.float32), np.float32),
-        (
-            "day_of_year",
-            np.repeat([1, 2, 3], 4).astype(np.float32),
-            np.float32,
-        ),
-        (
-            "day_of_year_sin",
-            np.repeat([0.017166, 0.034328, 0.051479], 4).astype(np.float32),
+            "hour_of_day_sin",
+            np.tile(np.repeat([1, 0], 6), 3).astype(np.float32),
             np.float32,
         ),
         (
-            "day_of_year_cos",
-            np.repeat([0.999853, 0.999411, 0.998674], 4).astype(np.float32),
+            "hour_of_day_cos",
+            np.tile(np.repeat([0, -1], 6), 3).astype(np.float32),
             np.float32,
         ),
-        ("hour_of_day", np.tile([6, 6, 12, 12], 3).astype(np.float32), np.float32),
-        ("hour_of_day_sin", np.tile([1, 1, 0, 0], 3).astype(np.float32), np.float32),
-        ("hour_of_day_cos", np.tile([0, 0, -1, -1], 3).astype(np.float32), np.float32),
         (
             "forecast_period",
-            np.tile([21600, 21600, 43200, 43200], 3).astype(np.int32),
+            np.tile(np.repeat([21600, 43200], 6), 3).astype(np.int32),
             np.int32,
         ),
-        ("day_of_training_period", np.repeat([0, 1, 2], 4).astype(np.int32), np.int32),
-        ("static", np.tile([2, 3], 6).astype(np.float32), np.float32),
+        ("day_of_training_period", np.repeat([0, 1, 2], 12).astype(np.int32), np.int32),
+        ("static", np.tile([2, 3], 18).astype(np.float32), np.float32),
     ],
 )
-def test_prep_feature_2d_time_dimension(feature, expected, expected_dtype):
+def test_prep_feature_more_times(feature, expected, expected_dtype):
     """Test the prep_feature function for multiple times."""
     forecast_reference_times = [
         "20170101T0000Z",
@@ -381,43 +350,170 @@ def test_prep_feature_2d_time_dimension(feature, expected, expected_dtype):
         "20170103T1200Z",
     ]
 
-    data = [2, 6, 10]
-    forecast_cubes = CubeList()
+    data = np.array([2, 6, 10])
+
+    feature_name = "wind_speed_at_10m"
+
+    data = np.array([2, 6, 10])
+
+    forecast_dfs = []
     for frt, vt in zip(forecast_reference_times, validity_times):
-        forecast_cubes.append(_create_forecasts(frt, vt, data))
-    forecast_cube = forecast_cubes.merge_cube()
-    # Ensure coordinate order is forecast_reference_time, forecast_period,
-    # realization, spot_index
-    forecast_cube.transpose([1, 0, 2, 3])
-    day_of_training_period = [0, 1, 2]
-    forecast_cube = _add_day_of_training_period(
-        forecast_cube, day_of_training_period, "forecast_reference_time"
-    )
+        forecast_dfs.append(_create_forecasts(frt, vt, data))
+    forecast_df = pd.concat(forecast_dfs)
+    forecast_df = _add_day_of_training_period(forecast_df)
 
     if feature == "static":
-        feature_cube = _create_ancil_file()
+        feature_df = _create_ancil_file()
+        forecast_df = forecast_df.merge(
+            feature_df[["wmo_id", "distance_to_water"]], on=["wmo_id"], how="left"
+        )
+
+    else:
+        forecast_df = forecast_df.copy()
+
+    result = prep_feature(forecast_df, feature_name, feature)
+
+    if feature in ["mean", "std"]:
+        assert result.shape == (36, 12)
+        feature_name_modified = f"{feature_name}_{feature}"
+    elif feature in [
+        "day_of_year",
+        "day_of_year_sin",
+        "day_of_year_cos",
+        "hour_of_day",
+        "hour_of_day_sin",
+        "hour_of_day_cos",
+    ]:
+        assert result.shape == (36, 12)
+        feature_name_modified = feature
+    elif feature in ["static"]:
+        assert result.shape == (36, 12)
+        feature_name_modified = "distance_to_water"
     else:
-        feature_cube = forecast_cube.copy()
+        assert result.shape == (36, 11)
+        feature_name_modified = feature
 
-    result = prep_feature(forecast_cube, feature_cube, feature)
+    assert result[feature_name_modified].dtype == expected_dtype
+    np.testing.assert_allclose(result[feature_name_modified], expected, atol=1e-6)
+
+
+@pytest.mark.parametrize(
+    "feature_config",
+    [
+        {"wind_speed_at_10m": ["mean", "std", "latitude", "longitude"]},
+        {"wind_speed_at_10m": ["latitude", "longitude"]},
+        {
+            "wind_speed_at_10m": ["latitude", "longitude"],
+            "distance_to_water": ["static"],
+        },
+        {"wind_speed_at_10m": ["latitude", "longitude", "height"]},
+    ],
+)
+def test_sanitise_forecast_dataframe(feature_config):
+    """Test sanitise_forecast_dataframe function."""
+    data_dict = {
+        "wmo_id": np.tile(WMO_ID, 3),
+        "latitude": np.tile(LATITUDE, 3),
+        "longitude": np.tile(LONGITUDE, 3),
+        "altitude": np.tile(ALTITUDE, 3),
+        "wind_speed_at_10m_mean": np.repeat(5, 6),
+        "wind_speed_at_10m_std": np.repeat(1, 6),
+        "wind_speed_at_10m": np.tile([4, 6], 3),
+        "realization": np.tile([1, 2], 3),
+        "distance_to_water": np.tile([2.0, 3.0], 3),
+    }
+    df = pd.DataFrame(data_dict)
+
+    expected = df.copy()
+    if (
+        "mean" in feature_config["wind_speed_at_10m"]
+        or "std" in feature_config["wind_speed_at_10m"]
+    ):
+        expected.drop(
+            columns=[
+                "wind_speed_at_10m",
+            ],
+            inplace=True,
+        )
+    expected = expected[expected["realization"] == 1]
+    expected = expected.drop(columns=["realization"])
 
-    assert result.shape == (12,)
-    assert result.dtype == expected_dtype
-    np.testing.assert_allclose(result, expected, atol=1e-6)
+    result = sanitise_forecast_dataframe(df, feature_config)
+    assert_frame_equal(result, expected)
+
+
+@pytest.mark.parametrize(
+    "feature_config",
+    [
+        {"wind_speed_at_10m": ["mean", "std", "latitude", "longitude"]},
+        {"wind_speed_at_10m": ["latitude", "longitude"]},
+        {
+            "wind_speed_at_10m": ["latitude", "longitude"],
+            "distance_to_water": ["static"],
+        },
+        {"wind_speed_at_10m": ["latitude", "longitude", "height"]},
+    ],
+)
+def test_get_required_column_names(feature_config):
+    """Test the get_required_column_names function."""
+    data_dict = {
+        "wmo_id": np.tile(WMO_ID, 3),
+        "latitude": np.tile(LATITUDE, 3),
+        "longitude": np.tile(LONGITUDE, 3),
+        "altitude": np.tile(ALTITUDE, 3),
+        "wind_speed_at_10m_mean": np.repeat(5, 6),
+        "wind_speed_at_10m_std": np.repeat(1, 6),
+        "wind_speed_at_10m": np.tile([4, 6], 3),
+        "realization": np.tile([1, 2], 3),
+        "distance_to_water": np.tile([2.0, 3.0], 3),
+    }
+    df = pd.DataFrame(data_dict)
+
+    # expected = feature_config["wind_speed_at_10m"].copy()
+    expected = list(itertools.chain.from_iterable(feature_config.values()))
+    if "mean" in expected:
+        expected = ["wind_speed_at_10m_mean" if e == "mean" else e for e in expected]
+    if "std" in expected:
+        expected = ["wind_speed_at_10m_std" if e == "std" else e for e in expected]
+    if "static" in expected:
+        expected = ["distance_to_water" if e == "static" else e for e in expected]
+
+    if "height" in expected:
+        with pytest.raises(ValueError, match="Feature 'height' is not supported."):
+            get_required_column_names(df, feature_config)
+    else:
+        result = get_required_column_names(df, feature_config)
+        assert result == expected
+
+
+@pytest.mark.parametrize(
+    "transformation", ["log", "log10", "sqrt", "cbrt", None, "yeojohnson"]
+)
+def test_check_valid_transformation(transformation):
+    """Test the _check_valid_transformation function."""
+
+    if transformation == "yeojohnson":
+        with pytest.raises(
+            ValueError, match="Currently the only supported transformations"
+        ):
+            _check_valid_transformation(transformation)
+    else:
+        result = _check_valid_transformation(transformation)
+        assert result is None
 
 
 @pytest.mark.parametrize(
     "n_estimators,max_depth,random_state,compression,transformation,pre_transform_addition,extra_kwargs,include_static,expected",
     [
-        (2, 2, 55, 5, None, 0, {}, False, 4.15),  # Basic test case
-        (2, 2, 54, 5, None, 0, {}, False, 4.2),  # Different random state
-        (1, 1, 55, 5, None, 0, {}, False, 4.2),  # Fewer estimators and reduced depth
+        (2, 2, 55, 5, None, 0, {}, False, 4.2),  # Basic test case
+        (2, 2, 54, 5, None, 0, {}, False, 4.15),  # Different random state
+        (1, 1, 54, 5, None, 0, {}, False, 4.2),  # Fewer estimators and reduced depth
         (2, 2, 55, 5, "log", 10, {}, False, 2.65),  # Log transformation
         (2, 2, 55, 5, "log10", 10, {}, False, 1.15),  # Log10 transformation
         (2, 2, 55, 5, "sqrt", 10, {}, False, 3.76),  # Square root transformation
         (2, 2, 55, 5, "cbrt", 10, {}, False, 2.42),  # Cube root transformation
-        (2, 2, 55, 5, None, 0, {"criterion": "absolute_error"}, False, 4.15),  # noqa # Different criterion
-        (2, 5, 55, 5, None, 0, {}, True, 4.15),  # Include static data
+        (2, 2, 55, 5, None, 0, {"criterion": "absolute_error"}, False, 4.2),  # noqa # Different criterion
+        (2, 5, 55, 5, None, 0, {}, True, 4.2),  # Include static data
     ],
 )
 def test_train_qrf_single_lead_times(
@@ -456,9 +552,8 @@ def test_train_qrf_single_lead_times(
             "20170101T1200Z",
             "20170102T1200Z",
         ],
-        day_of_training_period=[0, 1],
         realization_data=[2, 6, 10],
-        truth_data=[[4.2, 4.2], [4.1, 4.1]],
+        truth_data=[4.2, 4.1, 4.2, 4.1],
     )
     qrf_model = joblib.load(model_output)
 
@@ -466,7 +561,7 @@ def test_train_qrf_single_lead_times(
     assert qrf_model.max_depth == max_depth
     assert qrf_model.random_state == random_state
 
-    current_forecast = [5, 3, 55, 5]
+    current_forecast = [7.5, 3, 55, 5]
     if include_static:
         current_forecast.append(2.5)
     result = qrf_model.predict(
@@ -478,15 +573,15 @@ def test_train_qrf_single_lead_times(
 @pytest.mark.parametrize(
     "n_estimators,max_depth,random_state,compression,transformation,pre_transform_addition,extra_kwargs,include_static,expected",
     [
-        (2, 2, 55, 5, None, 0, {}, False, 4.0),  # Basic test case
+        (2, 2, 55, 5, None, 0, {}, False, 5.8),  # Basic test case
         (1, 1, 55, 5, None, 0, {}, False, 3.8),  # Fewer estimators and reduced depth
-        (1, 1, 73, 5, None, 0, {}, False, 7),  # Different random state
-        (2, 2, 55, 5, "log", 10, {}, False, 2.743),  # Log transformation
-        (2, 2, 55, 5, "log10", 10, {}, False, 1.191),  # Log10 transformation
-        (2, 2, 55, 5, "sqrt", 10, {}, False, 3.946),  # Square root transformation
-        (2, 2, 55, 5, "cbrt", 10, {}, False, 2.496),  # Cube root transformation
-        (2, 2, 55, 5, None, 0, {"criterion": "absolute_error"}, False, 4),  # noqa # Different criterion
-        (2, 5, 55, 5, None, 0, {}, True, 4),  # Include static data
+        (1, 1, 73, 5, None, 0, {}, False, 5.8),  # Different random state
+        (2, 2, 55, 5, "log", 10, {}, False, 2.752),  # Log transformation
+        (2, 2, 55, 5, "log10", 10, {}, False, 1.195),  # Log10 transformation
+        (2, 2, 55, 5, "sqrt", 10, {}, False, 3.964),  # Square root transformation
+        (2, 2, 55, 5, "cbrt", 10, {}, False, 2.504),  # Cube root transformation
+        (2, 2, 55, 5, None, 0, {"criterion": "absolute_error"}, False, 5.8),  # noqa # Different criterion
+        (1, 1, 73, 5, None, 0, {}, True, 3.8),  # Include static data
     ],
 )
 def test_train_qrf_multiple_lead_times(
@@ -524,9 +619,9 @@ def test_train_qrf_multiple_lead_times(
     assert qrf_model.max_depth == max_depth
     assert qrf_model.random_state == random_state
 
-    current_forecast = [5, 3, 55, 5]
+    current_forecast = [7.5, 3, 55, 5]
     if include_static:
-        current_forecast.append(2.5)
+        current_forecast.append(3.2)
     result = qrf_model.predict(
         np.expand_dims(np.array(current_forecast), 0), quantiles=[0.5]
     )
@@ -536,20 +631,20 @@ def test_train_qrf_multiple_lead_times(
 @pytest.mark.parametrize(
     "feature_config,data,include_static,expected",
     [
-        ({"wind_speed_at_10m": ["mean"]}, [5], False, [4]),  # One feature
-        ({"wind_speed_at_10m": ["latitude"]}, [55], False, [6.65]),  # noqa Without the target
+        ({"wind_speed_at_10m": ["mean"]}, [5], False, [5]),  # One feature
+        ({"wind_speed_at_10m": ["latitude"]}, [61], False, [7.75]),  # noqa Without the target
         ({"wind_speed_at_10m": ["mean"]}, [5], True, [4]),  # With static data
         (
             {"wind_speed_at_10m": ["mean"], "air_temperature": ["mean"]},
             [5],
             False,
-            [4],
+            [5],
         ),  # Multiple dynamic features
         (
             {"wind_speed_at_10m": ["mean"], "pressure_at_mean_sea_level": ["mean"]},
             [5],
             False,
-            "Feature cube for pressure_at_mean_sea_level",
+            "Feature 'pressure_at_mean_sea_level' is not present",
         ),  # Multiple dynamic features
     ],
 )
@@ -606,7 +701,7 @@ def test_alternative_feature_configs(
 
     current_forecast = [*data]
     if include_static:
-        current_forecast.append(2.5)
+        current_forecast.append(3.2)
     if "air_temperature" in feature_config:
         current_forecast.append(15.0)
 
@@ -619,31 +714,31 @@ def test_alternative_feature_configs(
 @pytest.mark.parametrize(
     "quantiles,transformation,pre_transform_addition,include_static,expected",
     [
-        ([0.5], None, 0, False, [4, 4]),  # One quantile
-        # ([0.1, 0.5, 0.9], None, 0, False, [[9.1, 9.14], [9.1, 9.3], [9.1, 9.46]]),  # noqa Multiple quantiles
-        # ([0.1, 0.5, 0.9], "log", 10, False, [[4.46, 4.46], [5.54, 5.54], [6.7, 6.7]]),  # noqa Log transformation
-        # (
-        #     [0.1, 0.5, 0.9],
-        #     "log10",
-        #     10,
-        #     False,
-        #     [[4.46, 4.46], [5.54, 5.54], [6.7, 6.7]],
-        # ),  # Log10 transformation
-        # (
-        #     [0.1, 0.5, 0.9],
-        #     "sqrt",
-        #     10,
-        #     False,
-        #     [[4.47, 4.47], [5.57, 5.57], [6.71, 6.71]],
-        # ),  # Square root transformation
-        # (
-        #     [0.1, 0.5, 0.9],
-        #     "cbrt",
-        #     10,
-        #     False,
-        #     [[4.47, 4.47], [5.56, 5.56], [6.7, 6.7]],
-        # ),  # Cube root transformation
-        # ([0.1, 0.5, 0.9], None, 0, True, [[9.1, 9.14], [9.1, 9.3], [9.1, 9.46]]),  # noqa Include static data
+        ([0.5], None, 0, False, [5, 4.9]),  # One quantile
+        ([0.1, 0.5, 0.9], None, 0, False, [[7.25, 8.25, 9.25], [6.35, 7.75, 9.14]]),  # noqa Multiple quantiles
+        ([0.1, 0.5, 0.9], "log", 10, False, [[4.48, 5.67, 6.96], [4.48, 5.67, 6.96]]),  # noqa Log transformation
+        (
+            [0.1, 0.5, 0.9],
+            "log10",
+            10,
+            False,
+            [[4.48, 5.67, 6.96], [4.48, 5.67, 6.96]],
+        ),  # Log10 transformation
+        (
+            [0.1, 0.5, 0.9],
+            "sqrt",
+            10,
+            False,
+            [[4.5, 5.71, 6.98], [4.5, 5.71, 6.98]],
+        ),  # Square root transformation
+        (
+            [0.1, 0.5, 0.9],
+            "cbrt",
+            10,
+            False,
+            [[4.49, 5.7, 6.97], [4.49, 5.7, 6.97]],
+        ),  # Cube root transformation
+        ([0.1, 0.5, 0.9], None, 0, True, [[7.25, 8.25, 9.25], [6.35, 7.75, 9.15]]),  # noqa Include static data
     ],
 )
 def test_apply_qrf(
@@ -678,30 +773,31 @@ def test_apply_qrf(
 
     frt = "20170103T0000Z"
     vt = "20170103T1200Z"
-    day_of_training_period = 2
     data = np.arange(6, (len(quantiles) * 6) + 1, 6)
 
-    forecast_cube = _create_forecasts(frt, vt, data)
-    forecast_cube = _add_day_of_training_period(
-        forecast_cube, day_of_training_period, "forecast_reference_time"
-    )
+    forecast_df = _create_forecasts(frt, vt, data)
+    forecast_df = _add_day_of_training_period(forecast_df)
 
-    feature_cubes = CubeList([forecast_cube])
     if include_static:
-        ancil_cube = _create_ancil_file()
-        feature_cubes.append(ancil_cube)
+        ancil_df = _create_ancil_file()
+        forecast_df = forecast_df.merge(
+            ancil_df[["wmo_id", "distance_to_water"]], on=["wmo_id"], how="left"
+        )
 
     plugin = ApplyQuantileRegressionRandomForests(
-        feature_config, quantiles, transformation, pre_transform_addition
+        "air_temperature",
+        feature_config,
+        quantiles,
+        transformation,
+        pre_transform_addition,
     )
-    result = plugin.process(qrf_model, feature_cubes, forecast_cube)
+    result = plugin.process(qrf_model, forecast_df)
 
-    assert isinstance(result, Cube)
-    assert result.shape == (len(quantiles), 2)
+    assert isinstance(result, np.ndarray)
+    if len(quantiles) == 3:
+        assert result.shape == (2, 3)
+    else:
+        assert result.shape == (2,)
     assert result.dtype == np.float32
-    expected_cube = forecast_cube.copy(
-        data=np.full((len(quantiles), 2), expected, dtype=np.float32)
-    )
-    assert result.name() == expected_cube.name()
-    assert result.units == expected_cube.units
-    np.testing.assert_almost_equal(result.data, expected_cube.data, decimal=2)
+    # expected_data = np.full((len(quantiles), 2), expected, dtype=np.float32)
+    np.testing.assert_almost_equal(result, expected, decimal=2)