improver.calibration.rainforest_calibration module

RainForests calibration Plugins.

RainForests calibration

RainForests calibration is a situation dependent non-parametric method for the calibration of rainfall based on the ECPoint method of Hewson and Pillosu Hewson & Pillosu, 2021.

Sub-grid variability as a means of calibration

ECPoint is based on the principle of using sub-grid variability as a means to calibrate grid-scale rainfall forecasts. Sub-grid variability in this context is the connection between the point observations one would measure within a grid-box when presented with a given areal average grid-scale forecast value. This relationship is encapsulated in the mapping that takes a grid-scale forecast value and maps it onto a distribution of typical point observation values, or equivalently a set of forecast error values, and is central to this calibration method.

Naturally the relationship between this distribution and the forecast value is contingent on the underlying rainfall processes at play. For instance, the distribution associated with a rainband will be characteristically different to that associated with post-frontal showers and different still to that associated with deep-tropical convection. To this end, a suitable set of meteorological variables can be used to distinguish different rainfall regimes and identify the associated distribution that describes the underlying sub-grid variability.

Equipped with the relevant distribution that describes the sub-grid variability for a given forecast, the grid-scale input forecast value can be mapped onto a series of realisable point-scale values, yielding a form of conditional bias-correction. This process is done on a per realization basis, using the mapping function most consistent with the input ensemble member forecast. In this way we calibrate to point scale and correct model bias within each ensemble member independently, with each grid-scale ensemble member producing a series of realisable point-scale forecast values.

The output produced from each ensemble member can be considered a pseudo-ensemble which represents a conditional probability distribution that describes the likelihood of observing a given outcome when the realised atmospheric state is consistent with that represented in the input ensemble member forecast.

Combining these ensemble member pseudo-ensembles, we can produce an intra-model “super-ensemble”. Usage of the term super-ensemble here is distinct from the more common usage which refers to the inter-model super-ensemble formed by combining multiple NWP ensemble models. Herein the usage of super-ensemble refers to the intra-model super-ensemble. Rather than containing a range of atmospheric states to which a single outcome is associated, the super-ensemble describes a range of possible outcomes given sourced from each atmospheric state.

This super-ensemble is represented by a series of realization values, and the associated probability distribution can be sourced using the same approach one would apply to the input forecast ensemble.

Hewson & Pillosu propose a framework (ECPoint) to determine and apply these distributions across a variety of weather types through the use of a manually constructed decision tree model taking an appropriate set of feature parameters (meteorological variables) as inputs. Details of this method can be found in Hewson & Pillosu, 2021.

One advantage of the ECPoint approach is that the calibration is inherently non-local. As the calibration is done by identifying distinct weather types, the model bias and scale difference should be independent of any given location and time as the underlying physical process should be identical. Thus a grid-point can be calibrated using data from any location, provided the underlying weather type is the consistent. This enables effective calibration to be applied to areas that are typically lacking sufficient cases to calibrate against.

The RainForests method

RainForests is an adaptation of the ECPoint method that seeks to use machine learning based tree models to replace the manually constructed tree model of ECPoint. Here we use gradient-boosted decision tree (GBDT) ensemble models.

The underlying principle of the calibration methodology is essentially the same, namely using a set of feature variables to map onto error distributions which can be applied to the forecast values to produce calibrated forecast values. However, the way in which these distributions are constructed differs somewhat within the RainForests framework.

Our approach is to use a series of GBDT models, taking the feature variables as inputs, to produce exceedance probability values for representative error thresholds. Collectively these represent the error distributions as a cumulative distribution function.

Each CDF is then mapped onto a series of equispaced percentile values to yield a series of representative error values which can be applied to the grid-scale forecast value to produce the series of possible point-scale (calibrated) values.

Calibration of the ensemble forecast proceeds by determining the underlying distribution on a per realization basis and applying this to each forecast value to produce a series of distinct calibrated forecast values for each realization. Collectively these values form the so-called super-ensemble, which we subsequently sample to produce the calibrated forecast ensemble output.

Another point of difference between RainForests and ECPoint is in the choice of error value underlying the error distribution. For RainForests we have chosen to use additive error in place of a multiplicative error (forecast error ratio in ECPoint) to allow us to calibrate input forecast values of zero-rainfall.

GBDT vs manually constructed DT

The choice of using GBDT models in place of the manually constructed DT of ECPoint comes with some advantages, but at the expense of some trade-offs:

Advantages:

• GBDT is an ensemble of many trees rather than a single tree. This means outputs are near-continuous relative to the inputs.

• Trees are built algorithmically, not manually, with each branch split automatically chosen to be optimal relative to some loss function. In principle this gives better accuracy, and makes it easier to retrain on new data.

trade-offs:

• By using an ensemble of trees, the intuitive connection between weather type and feature variables becomes obscured.

• Using a series of decision tree ensembles in place of a single decision tree increases the computational demand significantly.

• Some initial effort is required to select a good set of model hyper-parameters that neither under- or over-fit. This process is more challenging and less transparent than ECPoint, however is required only once rather than each time the decision tree(s) are constructed.

Implementation details

Model training

The model training process is relatively simple and involves collating a series of forecast-observation pairs with the associated feature variables into a single pandas dataframe.

Using this dataframe, a series of lightGBM binary-classification models are trained against the truth values for exceedances of the forecast error relative to a series of error-threshold values, resulting in one model per error threshold.

As the calibration method works on a per realization basis, the errors used in training must be representative of the underlying systematic model biases for a single ensemble member rather than those of the ensemble as a whole. This requires the use of a single representative ensemble forecast member in training. It also requires the underlying weather type associated with the forecast and observation in the forecast-observation pair be as consistent as possible.

For the forecast values used in training, we use the control at the shortest available lead-times that distinctly represent of the underlying weather types. Specifically, we take the earliest 24-hour period for daily accumulations, and all lead-times over the earliest 24-hour period for sub-daily accumulations to capture the full cycle of diurnal variations.

Currently model training is done offline, using a minimum 12-month period to capture the full seasonal cycle.

Forecast calibration

Forecast calibration uses the trained GBDT models, along with the forecast cube and associated feature cubes. The tree-models are passed in via a model-config json which identifies the appropriate tree-model file for each error-threshold.

Forecast calibration proceeds via a 2-step process:

1. Evaluate the error CDF defined over the series of error-thresholds used in model training. Each exceedance probability is evaluated using the corresponding tree-model, and the feature variables as inputs.

2. Interpolate the CDF to extract a series of percentile values for the error distributions. The error percentiles are then added to each associated ensemble realization from the forecast variable to produce a series of realisable forecast values.

Collectively these series form the calibrated super-ensemble which is obtained by collapsing the two realization dimensions into one. This is then sampled to provide the calibrated ensemble forecast.

Deterministic forecasts can also be calibrated using the same approach to produce a calibrated pseudo-ensemble; in this case inputs are treated as an ensemble of size 1.

A typical usage example: we typically use around 25 error threshold values to construct the CDF for the distribution of forecast errors. For each error threshold we have an associated GBDT model which is used to evaluate the exceedance probabilities that describe the CDF. So starting with an input ensemble forecast consisting of 50 realizations, we evaluate 25 threshold probability values for each realization to construct a forecast error CDF for each realization (50 distributions in total, each containing 25 threshold values).

For each distribution, we then interpolate between the threshold probabilities to extract 20 evenly-spaced percentiles. These are then applied to each forecast realization to produce 20 calibrated forecast realizations, resulting in 50 * 20 (1000) forecast values which collective for the calibrated “super-ensemble”. Finally, we sample the super-ensemble by taking 100 equispaced percentile values to be representative realizations for the calibrated forecast ensemble.

This final step is not required, but ensures efficient processing in downstream CLIs.

The number of error-percentiles used in the interim step, and the number of output ensemble realizations are taken as input parameters.

class ApplyRainForestsCalibration(model_config_dict: dict, threads: int = 1)

Bases: PostProcessingPlugin

Generic class to calibrate input forecast via RainForests.

The choice of tree-model library is determined from package availability, and whether all required models files are available. Treelite is the preferred option, defaulting to lightGBM if requirements are missing.

_abc_impl = <_abc_data object>

process()

Subclasses should override this function.

Return type

None

class ApplyRainForestsCalibrationLightGBM(model_config_dict, threads=1)

Bases: ApplyRainForestsCalibration

Class to calibrate input forecast given via RainForests approach using light-GBM tree models

__init__(model_config_dict, threads=1)

Initialise the tree model variables used in the application of RainForests Calibration. LightGBM Boosters are used for tree model predictors.

Parameters

• model_config_dict (dict) – Dictionary containing Rainforests model configuration variables.

• threads (int) – Number of threads to use during prediction with tree-model objects.

Dictionary is of format:

{
    "-50.0" : {
        "lightgbm_model" : "<path_to_lightgbm_model_object>"
    },
    "-25.0" : {
        "lightgbm_model" : "<path_to_lightgbm_model_object>"
    },
    ...,
    "50.0" : {
        "lightgbm_model" : "<path_to_lightgbm_model_object>"
    }
}

The keys specify the error threshold value, while the associated values are the path to the corresponding tree-model objects for that threshold.

_abc_impl = <_abc_data object>

_align_feature_variables(feature_cubes, forecast_cube)

Ensure that feature cubes have consistent dimension coordinates. If realization dimension present in any cube, all cubes lacking this dimension will have realization dimension added and broadcast along this new dimension.

This situation occurs when derived fields (such as accumulated solar radiation) are used as predictors. As these fields do not contain a realization dimension, they must be broadcast to match the NWP fields that do contain realization, so that all features have consistent shape.

In the case of deterministic models (those without a realization dimension), a realization dimension is added to allow consistent behaviour between ensemble and deterministic models.

Parameters

• feature_cubes (CubeList) – Cubelist containing feature variables to align.

• forecast_cube (Cube) – Cube containing the forecast variable to align.

Return type

Tuple[CubeList, Cube]

Returns

• feature_cubes with realization coordinate added to each cube if absent

• forecast_cube with realization coordinate added if absent

Raises

ValueError – if feature/forecast variables have inconsistent dimension coordinates (excluding realization dimension), or if feature/forecast variables have different length realization coordinate over cubes containing a realization coordinate.

_apply_error_to_forecast(forecast_cube, error_percentiles_cube)

Apply the error distributions (as error percentiles) to the forecast cube. The result is a series (sub-ensemble) of values for each forecast realization.

Note

Within the RainForests approach we work with an additive error correction as opposed to a multiplicative correction used in ECPoint. The advantage of using an additive error is that we are also able to calibrate zero-values in the input forecast.

Warning

After applying the error distributions to the forecast cube, values outside the expected bounds of the forecast parameter can arise. These values occur when when the input forecast value is between error thresholds and there exists a lower bound on the observable value (eg. 0 in the case of rainfall).

In this situation, error thresholds below the residual value (min(obs) - fcst) must have a probability of exceedance of 1, whereas as error thresholds above this value can take on any value between [0, 1]. In the subsequent step where error percentile values are extracted, the linear interpolation in mapping from probabilities to percentiles can give percentile values that lie below the residual value; when these are applied to the forecast value, they result in forecast values outside the expected bounds of the forecast parameter in the resultant sub-ensemble.

To address this, we remap all values outside of the expected bounds to nearest bound (eg. negative values are mapped to 0 in the case of rainfall).

Parameters

• forecast_cube (Cube) – Cube containing the forecast to be calibrated.

• error_percentiles_cube (Cube) – Cube containing percentile values for the error distributions.

Return type

Cube

Returns

Cube containing the forecast sub-ensembles.

_calculate_error_probabilities(forecast_cube, feature_cubes)

Evaluate the error exceedence probabilities for forecast_cube using the tree_models, with the associated feature_cubes taken as inputs to the tree_model predictors.

Parameters

• forecast_cube (Cube) – Cube containing the variable to be calibrated.

• feature_cubes (CubeList) – Cubelist containing the independent feature variables for prediction.

Return type

Cube

Returns

A cube containing error exceedence probabilities.

Raises

ValueError – If an unsupported model object is passed. Expects lightgbm Booster, or treelite_runtime Predictor (if treelite dependency is available).

_check_num_features(features)

Check that the correct number of features has been passed into the model.

Parameters

features (CubeList) – Cubelist containing feature variables.

Return type

None

_combine_subensembles(forecast_subensembles, output_realizations_count)

Combine the forecast sub-ensembles into a single ensemble. This is done by first stacking the sub-ensembles into a single super-ensemble and then resampling the super-ensemble to produce a subset of output realizations.

Parameters

• forecast_subensembles (Cube) – Cube containing a series of forecast sub-ensembles.

• output_realizations_count (Optional[int]) – The number of ensemble realizations that will be extracted from the super-ensemble. If realizations_count is None, all realizations will be returned.

Return type

Cube

Returns

Cube containing single realization dimension.

_evaluate_probabilities(forecast_data, input_data, forecast_variable, forecast_variable_unit, output_data)

Evaluate probability that error in forecast exceeds thresholds, setting the result to 1 when forecast + threshold is less than or equal to the lower bound of forecast_variable, as defined in constants.BOUNDS_FOR_ECDF`.

Parameters

• forecast_data (ndarray) – 1-d containing data for the variable to be calibrated.

• input_data (ndarray) – 2-d array of data for the feature variables of the model

• forecast_variable (str) – name of forecast variable

• forecast_variable_unit (str) – unit of forecast variable

• output_data (ndarray) – array to populate with output; will be modified in place

Return type

None

_extract_error_percentiles(error_probability_cube, error_percentiles_count)

Extract error percentile values from the error exceedence probabilities.

Parameters

• error_probability_cube – A cube containing error exceedence probabilities.

• error_percentiles_count – The number of error percentiles to extract. The resulting percentiles will be evenly spaced over the interval (0, 100).

Returns

Cube containing percentile values for the error distributions.

_make_decreasing(probability_data)

Enforce monotonicity on the error CDF data, where threshold dimension is assumed to be the leading dimension.

This is achieved by identifying the minimum value progressively along the leading dimension by comparing to all preceding probability values along this dimension. The same is done for maximum values, comparing to all succeeding values along the leading dimension. Averaging these resulting arrays results in an array decreasing monotonically in the threshold dimension.

Parameters

probability_data (ndarray) – The error probability data as exceedence probabilities.

Return type

ndarray

Returns

The error probability data, enforced to be monotonically decreasing along the leading dimension.

_prepare_error_probability_cube(forecast_cube)

Initialise a cube with the same dimensions as the input forecast_cube, with an additional threshold dimension added as the leading dimension.

Parameters

forecast_cube – Cube containing the forecast to be calibrated.

Returns

An empty probability cube.

_prepare_features_array(feature_cubes)

Convert gridded feature cubes into a numpy array, with feature variables sorted alphabetically.

Note: It is expected that feature_cubes has been aligned using _align_feature_variables prior to calling this function.

Parameters

feature_cubes (CubeList) – Cubelist containing the independent feature variables for prediction.

Return type

ndarray

Returns

Array containing flattened feature variables,

Raises

ValueError – If flattened cubes have differing length.

_stack_subensembles(forecast_subensembles)

Stacking the realization and percentile dimensions in forecast_subensemble into a single realization dimension. Realization and percentile are assumed to be the first and second dimensions respectively.

Parameters

input_cube – Cube containing the forecast_subensembles.

Return type

Cube

Returns

Cube containing single realization dimension in place of the realization and percentile dimensions in forecast_subensemble.

Raises

ValueError – if realization and percentile are not the first and second dimensions.

process(forecast_cube, feature_cubes, error_percentiles_count=19, output_realizations_count=100)

Apply rainforests calibration to forecast cube.

Ensemble forecasts must be in realization representation. Deterministic forecasts can be processed to produce a pseudo-ensemble; a realization dimension will be added to deterministic forecast cubes if one is not already present.

The calibration is done in a situation dependent fashion using a series of decision-tree models to construct representative error distributions which are then used to map each input ensemble member onto a series of realisable values.

These error distributions are formed in a two-step process:

1. Evaluate error CDF defined over the specified error_thresholds. Each exceedence probability is evaluated using the corresponding decision-tree model.

2. Interpolate the CDF to extract a series of percentiles for the error distribution. The error percentiles are then applied to each associated ensemble realization to produce a series of realisable values; collectively these series form a calibrated super-ensemble, which is then sub-sampled to provide the calibrated forecast.

Parameters

• forecast_cube (Cube) – Cube containing the forecast to be calibrated; must be as realizations.

• feature_cubes (CubeList) – Cubelist containing the feature variables (physical parameters) used as inputs to the tree-models for the generation of the associated error distributions. Feature cubes are expected to have the same dimensions as forecast_cube, with the exception of the realization dimension. Where the feature_cube contains a realization dimension this is expected to be consistent, otherwise the cube will be broadcast along the realization dimension.

• error_percentiles_count (int) – The number of error percentiles to extract from the associated error CDFs evaluated via the tree-models. These error percentiles are applied to each ensemble realization to produce a series of values, which collectively form the super-ensemble. The resulting super-ensemble will be of size = forecast.realization.size * error_percentiles_count.

• output_realizations_count (int) – The number of ensemble realizations that will be extracted from the super-ensemble. If realizations_count is None, all realizations will be returned.

Return type

Cube

Returns

The calibrated forecast cube.

Raises

RuntimeError – If the number of tree-models is inconsistent with the number of error thresholds.

class ApplyRainForestsCalibrationTreelite(model_config_dict, threads=1)

Bases: ApplyRainForestsCalibrationLightGBM

Class to calibrate input forecast given via RainForests approach using treelite compiled tree models

__init__(model_config_dict, threads=1)

Initialise the tree model variables used in the application of RainForests Calibration. Treelite Predictors are used for tree model predictors.

Parameters

• model_config_dict (dict) – Dictionary containing Rainforests model configuration variables.

• threads (int) – Number of threads to use during prediction with tree-model objects.

Dictionary is of format:

{
    "-50.0" : {
        "treelite_model" : "<path_to_treelite_model_object>"
    },
    "-25.0" : {
        "treelite_model" : "<path_to_treelite_model_object>"
    },
    ...,
    "50.0" : {
        "treelite_model" : "<path_to_treelite_model_object>"
    }
}

The keys specify the error threshold value, while the associated values are the path to the corresponding tree-model objects for that threshold.

_abc_impl = <_abc_data object>

_check_num_features(features)

Check that the correct number of features has been passed into the model. :type features: CubeList :param features: Cubelist containing feature variables.

Return type

None