2.1 Loading and harmonizing data

First, we load the required dataset resulting from the study in Mirones et al.(2022) located at Zenodo open repository [https://zenodo.org/record/7014397/files/South_Pacific_precipitation.zip].

###download of the dataset from Zenodo repository
temp <- tempfile(fileext = ".zip")
download.file("https://zenodo.org/record/7014397/files/South_Pacific_precipitation.zip",temp)

The downloaded dataset comprises a collection of stations located in the South Pacific region, obtained from the PACRAIN (Pacific Rainfall Database) [http://pacrain.ou.edu/]. Additionally, precipitation series from the TRMM (Tropical Rainfall Measuring Mission) satellite [https://trmm.gsfc.nasa.gov/] and the ERA5 reanalysis [https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5] have been extracted for the following stations: Kolopelu (Wallis and Futuna), Alofi (Niue), Rarotonga (Cook Islands), Raoul Island (New Zealand), Port Vila (Vanuatu), Aoloau, and Nuu’uli (American Samoa). Daily variables are included for the stations with the following IDs: NZ75400, NZ82400, NZ84317, NZ99701, SP00646, US14000, and US14690.

• TRMM raw satellite precipitation

• ERA5 raw reanalysis precipitation

• PACRAIN rain gauge precipitation

• Weather type associated

• TRMM calibrated (scaling)

• TRMM calibrated (empirical quantile mapping)

• TRMM calibrated (scaling conditioned)

• TRMM calibrated (empirical quantile mapping conditioned)

where the last four are a consequence of the application of 4 different calibration methods.

The required variables are rain gauge precipitation (pr), satellite precipitation (pp_trmm) and the WTs (wt). Thus, we set the corresponding variable names in the argument var when calling to loadStationData from loadeR:

###extracting the variables from the dataset
obs <- loadStationData(dataset = temp, var = "pr")
trmm <- loadStationData(dataset = temp, var = "pp_trmm")
wt <- loadStationData(dataset = temp, var = "wt")

To select the station for analysis, you can use the station.id argument in the subsetGrid function from the transformeR package. The station IDs can be found in $Metadata. We replicate the methodology outlined in the article for the Kolopelu station, as the code for implementing the adaptive methodology is the same for all stations. The only difference lies in the choice of the station.id argument in the subsetGrid function.

###selecting the station to be analyzed
obs <- subsetGrid(obs, station.id = obs$Metadata$station_id[1])
trmm <- subsetGrid(trmm, station.id = obs$Metadata$station_id)
wt <- subsetGrid(wt, station.id = obs$Metadata$station_id)

Lastly, we perform an intersection of the rain gauge and TRMM data to account for potential differences in their time coverage. To achieve this, we utilize the intersectGrid function from the transformeR package. The which.return argument specifies which grid should be returned, encompassing the overlapping time period. In this case, the observation series is returned as the target period for calibration is based on the TRMM data.

obs <- intersectGrid(obs, trmm, which.return = 1)

2.2 Adaptive calibration methodology and evaluation

The computation of the adaptive calibration for the Kolopelu TRMM series is carried out using a function called adaptiveCalibration, which encapsulates the adaptive calibration methodology presented in the article.

For a proper evaluation and comparison of different calibration methods (scaling, eQM, pQM, gpQM-95, and gpQM-75), a consistent scoring methodology is employed, as outlined by Kotlarski et al. (2019).

To ensure code reproducibility, the ./utils directory has been provided, which contains the auxiliary files adaptiveCalibration.R, RFScore.R, and MaxReturnValue.R necessary for implementing the adaptive calibration methodology.

In this methodology, the calibration methods are applied individually for each weather type (WT), and the final calibrated series is determined by selecting the best calibration method for each WT. This process results in a single time series that covers the entire period.

To determine the best calibration method for each WT, a score calculation based on Kotlarski et al. (2019) is utilized. The score computation follows these steps:

1. Climate indices from the VALUE database (http://www.value-cost.eu/validationportal/app/#!indices) are calculated using the climate4R.value package, which serves as a validation reference.

2. The absolute difference between the observed series and the climate indices is computed.

3. The difference in absolute value for each index is normalized across all methods. The resulting score for each method is the average of the normalized values of the climate indices.

This process of the score calculation is encapsulated in a atomic function called RFScore. In RFScore is possible to assign arbitrary weights through weights argument to the indices involved in the RF score, giving more importance to specific precipitation characteristics, such as the representation of extremes. The weighting aims to guide the calibration towards better adjustment in the upper tail of the distribution, thereby achieving improved correction for extreme precipitation events beyond a certain threshold. This allows obtaining a score that focuses more on specific indices, such as extreme precipitation (e.g., R20p or P98WetAmount).

The cross-validation scheme and calibration methods are implemented within the calibration and empirical statistical downscaling tools provided by the downscaleRpackage (Bedia et al., 2020) in the climate4R framework.

The validation of the results is based on the VALUE climatic indices. However, the adaptiveCalibration function allows the inclusion of other user-defined indices using the custom_function argument, such as the MaxReturnValue index.

Finally, the adaptive calibration is performed using the scaling, eQM, pQM, and gpQM (with 0.95 and 0.75 thresholds) techniques. These methods are also computed independently to evaluate whether the adaptive calibration results in improved performance compared to the standard calibration methods.

A few auxiliary functions have been prepared for paper calculations. These are next sourced from the utils directory:

###load of auxiliary functions
source("./utils/adaptiveCalibration.R")
source("./utils/RFscore.R")
source("./utils/MaxReturnValue.R")

###adaptive calibration computations
cal <- adaptiveCalibration(obs = obs, sim = trmm, 
                             clustering = wt$Data,
                             methods = c("scaling", "eqm", "pqm", 
                                         "gpqm", "gpqm"),
                             index = c("Skewness","Mean","SDII",
                                       "R10","R10p","R20","R20p",
                                       "P98Wet","P98WetAmount"),
                             custom_function = list(MaxReturnValue),
                             scaling.type = "multiplicative", 
                             window = NULL, 
                             theta = c(.95, .75))

###adaptive calibration including weigths
w <- c(.05, .05, .05, .05, .05, .15, .15, .15, .2, .1)
cal_w <- adaptiveCalibration(obs = obs, sim = trmm, 
                             clustering = wt$Data,
                             methods = c("scaling", "eqm", "pqm",
                                         "gpqm", "gpqm"),
                             index = c("Skewness","Mean","SDII",
                                       "R10","R10p",
                                       "R20","R20p","P98Wet",
                                       "P98WetAmount"),
                             weights = w,
                             custom_function = list(MaxReturnValue),
                             scaling.type = "multiplicative", 
                             window = NULL, 
                             theta = c(.95, .75))

2.3 Computation of standard calibration techniques

After completing the calculation of the series calibrated using the adaptive methodology, it is crucial to compute a series of standard calibrations to evaluate their performance and compare them with the adaptive methodology. For this purpose, the downscaleR package includes the biasCorrection function, which provides specific arguments depending on the calibration method selected.

1. The multiplicative argument determines the type of scaling method employed. It can be set to additive for unbounded variables like temperature or multiplicative for (lower)bounded variables like precipitation. If the bias correction method is not set to scaling, this argument is disregarded.

2. The fitdistr.args argument allows us to specify additional arguments passed to the fitdistr function, such as densfun, start, and others. This argument is only used when applying the pqm method to account for the different theoretical distribution function fitting.

3. The theta argument defines the upper threshold (and lower threshold for the left tail, if necessary) above which values are fitted to a Generalized Pareto Distribution (GPD). Values below this threshold are fitted to a gamma distribution. By default, theta corresponds to the 95th percentile (and 5th percentile for the left tail). This argument is exclusive to the gpqmmethod.

Hence, we compute the standard calibration for the scaling, eQM, pQM, and gpQM methods. In the case of the gpQM method, two series are generated: gpQM-95 and gpQM-75. The gpQM-95 series fits the data to a GPD above the 95th percentile, while the gpQM-75 series fits the data above the 75th percentile.

obs <- setGridDates.asPOSIXlt(obs,
                              tz = "UTC")
trmm <- setGridDates.asPOSIXlt(trmm,
                               tz = "UTC")

scaling <- biasCorrection(y = obs, x = trmm,
                          precipitation = TRUE,
                          method = "scaling",
                          scaling.type = "multiplicative",
                          window = NULL, 
                          cross.val = "loo")

eqm <- biasCorrection(y = obs, x = trmm,
                      precipitation = TRUE,
                      method = "eqm",
                      cross.val = "loo")

pqm <- biasCorrection(y = obs, x = trmm,
                      precipitation = TRUE,
                      method = "pqm",
                      cross.val = "loo")

gpqm95 <- biasCorrection(y = obs, x = trmm,
                         precipitation = TRUE,
                         method = "gpqm",
                         theta = .95,
                         cross.val = "loo")

gpqm75 <- biasCorrection(y = obs, x = trmm, 
                         precipitation = TRUE,
                         method = "gpqm",
                         theta = .75,
                         cross.val = "loo")