useqmap: Quantile Mapping Using the 'qmap' Package

In jonasbhend/biascorrection: Bias Correct Ensemble Forecasts

Quantile Mapping Using the qmap Package

Description

Computes bias correction with quantile mapping from the qmap package

Usage

useqmap(fcst, obs, fcst.out = fcst, nn = nrow(fcst), exact = FALSE, ...)

Arguments

fcst	n x m x k array of n lead times, m forecasts, of k ensemble members
obs	n x m matrix of veryfing observations
fcst.out	array of forecast values to which bias correction should be applied (defaults to fcst)
nn	number of lead times to be included around the lead time to be calibrated (i.e. ceiling((nn - 1)/2) before and floor((nn - 1)/2) after). The interval to estimate quantile corrections is held constant, that is, the first nn lead times are used to estimate the quantile correction for lead time one. The default is to lump all lead times together nn = nrow(fcst).
exact	logical, should surrounding lead-times be used for quantiles of forecast as well?
...	additional arguments passed to fitQmap and doQmap (see help(fitQmap) and help(doQmap))

Details

This function provides a wrapper to use the quantile mapping functionality from the qmap package. All lead times are lumped together, a potential seasonal cycle (lead-time dependency of quantile correction) is not explicitly taken into account. Similarly, for the bias corrected forecasts, all ensemble members are lumped together.

The quantile mapping is lead time dependent, parameter nn is used to select the number of lead times on either side of the lead time that is to be corrected to be included in the quantile estimation. For the beginning and end of the series, the lead-time interval is kept constant, so that to estimate the quantile correction for the first lead time, the first nn lead times are used. If exact = FALSE, the lead time dependent quantiles for the forecast are directly estimated from single lead times without the surrounding 2*nn lead times. Note, however, that unlike for qqmap, this does not speed up computation.

Examples

## initialise forcast observation pairs
nens <- 51
signal <- outer(sin(seq(0,4,length=215)), sort(rnorm(30, sd=0.2)), '+')
fcst <- array(rgamma(length(signal)*nens, shape=3), c(dim(signal), nens)) +
  c(signal)
obs <- array(rnorm(length(signal), mean=2), dim(signal)) + 
  signal
fcst.debias <- biascorrection:::useqmap(fcst[,1:20,], 
  obs[,1:20], fcst.out=fcst[,21:30,], method='RQUANT', qstep=0.05, wet.day=FALSE)
oprob <- (seq(obs[,21:30]) - 1/3) / (length(obs[,21:30]) + 1/3)
plot(quantile(obs[,21:30], type=8, oprob), 
  quantile(fcst[,21:30,], type=8, oprob),
  type='l', lwd=2, xlab='Observed quantiles',
  ylab='Forecast quantiles',
  main='Out-of-sample validation for qqmap')
abline(c(0,1), lwd=2, lty=2)
lines(quantile(obs[,21:30], type=8, oprob),
  quantile(fcst.debias, type=8, oprob), lwd=2, col=2)
minprob <- 0.05
abline(v=quantile(obs[,1:20], type=8, prob=c(minprob, 1-minprob)), 
  lwd=2, lty=3)
text(quantile(obs[,1:20], type=8, 1-minprob) + 0.1, par('usr')[3] + 0.5, 
  'Extrapolated\ncorrection', adj=c(0,0), cex=0.67)
text(quantile(obs[,1:20], type=8, 1-minprob) - 0.1, par('usr')[3] + 0.5, 
  'Explicit quantile\ncorrection', adj=c(1,0), cex=0.67)
legend('topleft', c('No bias correction', 'useqmap (RQUANT)'), 
  lwd=2, col=1:2, inset=0.05)

jonasbhend/biascorrection documentation built on Nov. 11, 2023, 1:16 a.m.