fitqmapptf: Quantile mapping using parametric transformations
In qmap: Statistical Transformations for Post-Processing Climate Model Output

Description Usage Arguments Details Value Author(s) References See Also Examples

fitQmapPTF fits a parametric transformations to the quantile-quantile relation of observed and modelled values. doQmapPTF uses the transformation to adjust the distribution of the modelled data to match the distribution of the observations.

fitQmapPTF(obs, mod, ...)
## Default S3 method:
fitQmapPTF(obs, mod, transfun=c("power","linear","expasympt",
"scale","power.x0","expasympt.x0"), wet.day=TRUE,
cost=c("RSS","MAE"), qstep=0.001,opar,...)
## S3 method for class 'matrix'
fitQmapPTF(obs, mod, ...)
## S3 method for class 'data.frame'
fitQmapPTF(obs, mod, ...)

doQmapPTF(x,fobj,...)
## Default S3 method:
doQmapPTF(x,fobj,...)
## S3 method for class 'matrix'
doQmapPTF(x,fobj,...)
## S3 method for class 'data.frame'
doQmapPTF(x,fobj,...)

`obs`	`numeric` vector, column `matrix` or `data.frame` with observed time series.
`mod`	`numeric` vector, column `matrix` or `data.frame` with modelled time series, corresponding to `obs`.
`transfun`	either a character string specifying a predefined function used for the transformation (see Details) or a function with `x` as first argument e.g. `function(x,a,b){a*x^b}`
`wet.day`	`logical` indicating whether to perform wet day correction or not. OR a `numeric` threshold below which all values are set to zero. See Details.
`cost`	Criterion for optimisation. "RSS" minimises the residual sum of squares and produces a least square fit. "MAE" minimises the mean absolute error, which is less sensitive to outliers.
`qstep`	`NULL` or a numeric value between 0 and 1. See Details.
`opar`	a named list with arguments passed to `optim`. Note that `method` is chosen automatically. If `transfun` is a character string default values for `par` are available (but can be overwritten). See examples.
`x`	`numeric` vector or a column `matrix` of modelled time series
`fobj`	output from `fitQmapDIST`
`...`	Further arguments passed to methods

Before further computations the empirical cumulative distribution functions (CDF) of the observed (obs) and modelled (mod) are estimated. If !is.null(qstep) than mod and obs are aggregated to quantiles before model identification as: quantile(x,probs=seq(0,1,by=qstep). If !is.null(qstep) than mod and obs are sorted to produce an estimate of the empirical CDF. In case of different length of mod and obs than quantile(x,probs=seq(0,1,len=n)] is used, where n <- min(length(obs),length(mod)). NOTE that large values of qstep effectively reduce the sample-size and can be used to speedup computations - but may render estimates less reliable.

wet.day is intended for the use for precipitation data. Wet day correction attempts to equalise the fraction of days with precipitation between the observed and the modelled data. If wet.day=TRUE the empirical probability of nonzero observations is found (obs>=0) and the corresponding modelled value is selected as a threshold. All modelled values below this threshold are set to zero. If wet.day is numeric the same procedure is performed after setting all obs<wet.day to zero. The transformations are then only fitted to the portion of the distributions corresponding to observed wet days. See Piani et. al (2010) for further explanations.

Transformations (transfun):

NOTE: If wet day correction is performed (see wet.day), the transformations are only fitted to the portion of the empirical CDF with nonzero observations.

A series of predefined transformations are available and can be accessed by setting transfun to one of the following options (P_o refers to observed and P_m to modelled CDFs):

"power":

P_o=b*P_m^c

"linear":

P_o=a+b*P_m

"expasympt" (exponential tendency to an asymptote):

P_o=(a+b*P_m)*(1-exp(-P_m/tau))

"scale":

P_o=b*P_m

"power.x0":

P_o=b*(P_m-x0)^c

"expasympt.x0" (exponential tendency to an asymptote):

P_o=(a+b*P_m)*(1-exp(-(P_m-x0)/tau))

fitQmapPTF returns an object of class fitQmapPTF containing following elements:

`tfun`	The function used to transform the distribution of the modelled values to match the distribution of the observations.
`par`	A matrix. The (named) columns correspond to the parameters of the transfer function. The rows correspond to pairs of time series in `obs` and `mod`.
`wet.day`	`logical`, indicating whether to perform wet day correction or not. OR a `numeric` threshold below which all values are set to zero.

doQmapPTF returns a numeric vector, matrix or data.frame depending on the format of x.

Lukas Gudmundsson

The implementation is closely related to the methods published in:

Piani, C.; Weedon, G.; Best, M.; Gomes, S.; Viterbo, P.; Hagemann, S. & Haerter, J. Statistical bias correction of global simulated daily precipitation and temperature for the application of hydrological models. Journal of Hydrology, 2010, 395, 199 - 215, doi:10.1016/j.jhydrol.2010.10.024.

Dosio, A. & Paruolo, P. Bias correction of the ENSEMBLES high-resolution climate change projections for use by impact models: Evaluation on the present climate. J. Geophys. Res., AGU, 2011, 116, D16106, doi:10.1029/2011JD015934.

For a general assessment of the methods see:

Gudmundsson, L.; Bremnes, J. B.; Haugen, J. E. & Engen-Skaugen, T. Technical Note: Downscaling RCM precipitation to the station scale using statistical transformations - a comparison of methods. Hydrology and Earth System Sciences, 2012, 16, 3383-3390, doi:10.5194/hess-16-3383-2012.

fitQmap, optim

data(obsprecip)
data(modprecip)

## data.frame example
qm.fit <- fitQmapPTF(obsprecip,modprecip,
                     transfun="power.x0",
                     cost="RSS",wet.day=TRUE,
                     qstep=0.001)
qm <- doQmapPTF(modprecip,qm.fit)

## application to "single time series"
qm.b.fit <- fitQmapPTF(obsprecip[,1],modprecip[,1],
                     transfun="expasympt.x0",
                     cost="RSS",wet.day=0.1,
                     qstep=0.001)
qm.b <- doQmapPTF(modprecip[,1],qm.b.fit)
qm.c.fit <- fitQmapPTF(obsprecip[,1],modprecip[,1],
                     transfun="expasympt",
                     cost="RSS",wet.day=TRUE,
                     qstep=0.001)
qm.c <- doQmapPTF(modprecip[,1],qm.c.fit)

## user defined transfer function
## and usage of the 'opar' argument
## (same as transfun="power")
myff <- function(x,a,b) a*x^b

qm3.fit <- fitQmapPTF(obsprecip[,1],modprecip[,1],
                     transfun=myff,
                     opar=list(par=c(a=1,b=1)),
                     cost="RSS",wet.day=TRUE,
                     qstep=0.001)
qm3 <- doQmapPTF(modprecip[,1],qm3.fit)


sqrtquant <- function(x,qstep=0.01){
  qq <- quantile(x,prob=seq(0,1,by=qstep))
  sqrt(qq)
}
plot(sqrtquant(modprecip[,1]),
     sqrtquant(obsprecip[,1]))
lines(sqrtquant(modprecip[,1]),
      sqrtquant(qm[,1]),col="red")
lines(sqrtquant(modprecip[,1]),
      sqrtquant(qm.b),col="blue")
lines(sqrtquant(modprecip[,1]),
      sqrtquant(qm.c),col="green")
lines(sqrtquant(modprecip[,1]),
      sqrtquant(qm3),col="orange")
legend("topleft",
       legend=c("power.x0","expasympt.x0",
         "expasympt","myff"),
       col=c("red","blue","green","orange"),lty=1)