High-resolution gridded forecast verification using discrete wavelet decomposition.
waverify2d(X, ..., Clim = NULL, wavelet.type = "haar", J = NULL)# S3 method for default
waverify2d(X, ..., Y, Clim = NULL, wavelet.type = "haar", J =
NULL, useLL = FALSE, compute.shannon = FALSE,
which.space = "field", verbose = FALSE)
# S3 method for SpatialVx
waverify2d(X, ..., Clim = NULL, wavelet.type = "haar", J = NULL,
useLL = FALSE, compute.shannon = FALSE, which.space = "field",
time.point = 1, obs = 1, model = 1, verbose = FALSE)
mowaverify2d(X, ..., Clim = NULL, wavelet.type = "haar", J = 4)
# S3 method for default
mowaverify2d(X, ..., Clim = NULL, Y, wavelet.type = "haar", J = 4,
useLL = FALSE, compute.shannon = FALSE, which.space = "field", verbose = FALSE)
# S3 method for SpatialVx
mowaverify2d(X, ..., Clim = NULL, wavelet.type = "haar", J = 4,
useLL = FALSE, compute.shannon = FALSE, which.space = "field",
time.point = 1, obs = 1, model = 1, verbose = FALSE)
# S3 method for waverify2d
plot(x, ..., main1 = "X", main2 = "Y", main3 = "Climate",
which.plots = c("all", "dwt2d", "details", "energy", "mse",
"rmse", "acc"), separate = FALSE, col, horizontal = TRUE)
# S3 method for waverify2d
print(x, ...)
A list object of class “waverify2d” with components:
single numeric giving the number of levels.
objects of class “dwt.2d” describing the wavelet decompositions for the verification and forecast fields (and climatology, if applicable), resp. (see the help file for dwt.2d from package waveslim for more about these objects).
numeric matrix giving the Shannon entropy for each field.
numeric matrix giving the energy at each level and field.
numeric vectors of length J giving the MSE/RMSE for each level between the verification and forecast fields.
If a climatology field is supplied, this is a numeric vector giving the ACC for each level.
m by n dyadic matrices (i.e., m = 2^M and n = 2^N, for M, N some integers) giving the verification and forecast fields (and optionally a climatology field), resp. Alternatively, X may be a “SpatialVx” object, in which case, Y is not given and in either case Clim must be provided if it is to be used.
list object of class “waverify2d” as returned by waverify2d.
character naming the type of wavelet to be used. This is given as the wf argument to the dwt.2d function of package waveslim.
(optional) numeric integer giving the pre-determined number of levels to use. If NULL, J is set to be log2(m) = M in waverify2d only.
logical, should the LL submatrix (i.e., the father wavelet or grand mean) be used to find the inverse DWT's for calculating the detail fields?
logical, should the Shannon entropy be calculated for the wavelet decomposition?
character (one of “field” or “wavelet”) naming from which space the detail fields should be used. If “field”, then it is in the original field (or image) space (i.e., the detail reconstruction), and if “wavelet”, it will be done in the wavelet space (i.e., the detail wavelet coefficients).
numeric or character indicating which time point from the “SpatialVx” verification set to select for analysis.
numeric indicating which observation/forecast model to select for the analysis.
optional characters naming each field to be used for the detail field plots and legend labelling on the energy plot.
character vector describing which plots to make. The default is to make all of them. “dwt2d” option uses plot.dwt2d from waveslim. “details” option makes image plots of the detail fields on which the various statistics are calculated. The rest of the options give line plots showing the statistics.
logical, should the plots be on their own devices (TRUE) or should some of them be put onto a single multi-panel device (FALSE, default)?
optional argument specifying the col argument in calls to functions like image. Default is a concatenation of gray with time.colors(64).
logical, should the legend on image plots be horizontal (TRUE, placed at the bottom of the plot) or vertical (FALSE, placed at the right side of the plot)?
logical, should progress information be printed to the screen, including total run time?
optional additonal plot or image.plot parameters. If detail and energy, mse, rmse or acc plots are desired, must be applicable to both types of plots. Not used by print method function.
Eric Gilleland
This is a function to use discrete wavelet decomposition to analyze verification sets along the lines of Briggs and Levine (1997), as well as Casati et al. (2004) and Casati (2009). In the originally proposed formulation of Briggs and Levine (1997), continuous verification statistics (namely, the anomaly correlation coefficient (ACC) and root mean square error (RMSE)) are calculated for detail fields obtained from wavelet decompositions of each of a forecast and verification field (and for ACC a climatology field as well). Casati et al. (2004) introduced an intensity scale approach that applies 2-d DWT to binary (obtained from thresholding) difference fields (Forecast - Verification), and applying a skill score at each level based on the mean square error (MSE). Casati (2009) extended this idea to look at the energy at each level as well.
This function makes use of the dwt.2d and idwt.2d functions from package waveslim, and plot.waverify2d uses the plot.dwt.2d function if dwt2d is selected through the which.plots argument. See the help file for these functions, the references therein and the references herein for more on these approaches.
Generally, it is not necessary to use the father wavelet for the detail fields, but for some purposes, it may be desired.
mowaverify2d is very similar to waverify2d, but it allows fields to be non-dyadic (and may subsequently be slower). It uses the modwt.2d and imodwt.2d functions from the package waveslim. In particular, it performs a maximal overlap discrete wavelet transform on a matrix of arbitrary dimension. See the help file and references therein for modwt.2d for more information, as well as Percival and Guttorp (1994) and Lindsay et al. (1996).
In Briggs and Levine (1997), they state that the calculations can be done in either the data (called field here) space or the wavelet space, and they do their examples in the field space. If the wavelets are orthogonal, then the detail coefficeints (wavelet space), can be analyzed with the assumption that they are independent; whereas in the data space, they typically cannot be assumed to be independent. Therefore, most statistical tests should be performed in the wavelet space to avoid issues arising from spatial dependence.
Briggs, W. M. and Levine, R. A. (1997) Wavelets and field forecast verification. Mon. Wea. Rev., 125, 1329--1341.
Casati, B., Ross, G. and Stephenson, D. B. (2004) A new intensity-scale approach for the verification of spatial precipitation forecasts. Meteorol. Appl. 11, 141--154.
Casati, B. (2010) New Developments of the Intensity-Scale Technique within the Spatial Verification Methods Inter-Comparison Project. Wea. Forecasting 25, (1), 113--143, doi:10.1175/2009WAF2222257.1.
Lindsay, R. W., Percival, D. B. and Rothrock, D. A. (1996) The discrete wavelet transform and the scale analysis of the surface properties of sea ice. IEEE Transactions on Geoscience and Remote Sensing, 34 (3), 771--787.
Percival, D. B. and Guttorp, P. (1994) Long-memory processes, the Allan variance and wavelets. In Wavelets in Geophysics, Foufoula-Georgiou, E. and Kumar, P., Eds., New York: Academic, pp. 325--343.
grid<- list( x= seq( 0,5,,64), y= seq(0,5,,64))
obj<-Exp.image.cov( grid=grid, theta=.5, setup=TRUE)
look<- sim.rf( obj)
look[ look < 0] <- 0
look <- zapsmall( look)
look2 <- sim.rf( obj)
look2[ look2 < 0] <- 0
look2 <- zapsmall( look2)
res <- waverify2d(look, Y=look2)
plot(res)
summary(res)
if (FALSE) {
data( "UKobs6" )
data( "UKfcst6" )
look <- waverify2d(UKobs6, Y=UKfcst6)
plot(look, which.plots="energy")
look2 <- mowaverify2d(UKobs6, UKfcst6, J=8)
plot(look2, which.plots="energy")
plot(look, main1="NIMROD Analysis", main2="NIMROD Forecast")
plot(look2, main1="NIMROD Analysis", main2="NIMROD Forecast")
# Alternative using "SpatialVx" object.
data( "UKloc" )
hold <- make.SpatialVx( UKobs6, UKfcst6, loc = UKloc,
map = TRUE, field.type = "Rainfall", units = "mm/h",
data.name = "Nimrod", obs.name = "Obs 6",
model.name = "Fcst 6" )
look <- waverify2d(hold)
plot(look, which.plots="details")
}
Run the code above in your browser using DataLab