SpatialVx (version 0.1-7)

wavePurifyVx: Apply Traditional Forecast Verification After Wavelet Denoising

Description

Apply traditional forecast verification after wavelet denoising ala Briggs and Levine (1997).

Usage

wavePurifyVx(x, y, object = NULL, climate = NULL, which.stats =
                 c("bias", "ts", "ets", "pod", "far", "f", "hk",
                 "mse"), thresholds = NULL, return.fields = FALSE,
                 time.point = 1, model = 1, verbose = FALSE, ...)
## S3 method for class 'wavePurifyVx':
plot(x, ..., set.pw = FALSE, col, zlim, horizontal = TRUE, 
    loc.byrow = TRUE, type = c("stats", "fields"))
## S3 method for class 'wavePurifyVx':
summary(object, ...)

Arguments

x,y

For wavePurifyVx, m by n matrices giving the verification and forecast fields, resp. If either is not supplied, then object must be supplied, and the fields given by object are used instead (even if one of x

object

For wavePurifyVx, a list object of class SpatialVx. Not used if both x and y are given. For summary.wavePurifyVx, list object as returned by wavePurifyVx.

climate

m by n matrix defining a climatology field. If not NULL, then the anamoly correlation coefficient will be applied to the wavelet denoised fields.

which.stats

character describing which traditional verification statistics to calculate on the wavelet denoised fields. This is the argument passed to the argument of the same name in vxstats.

thresholds

If object is used, then this argument is ignored. Otherwise, it is used to define events for all of the verification statistics except MSE. However, if supplied or other statistics are to be computed, then MSE will be calculated for the fie

return.fields

logical, should the denoised fields be returned (e.g., for subsequent plotting)?

time.point

numeric or character indicating which time point from the SpatialVx verification set to select for analysis.

model

numeric indicating which forecast model to select for the analysis.

verbose

logical, should progress information (including total run time) be printed to the screen?

set.pw

logical, should a panel of plots be determined and set by the function (generally ill advised).

col, zlim, horizontal

optional arguments to image, and/or fields functions poly.image and image.plot

loc.byrow

logical, only used if field is a projection, this determines how the locations should be put into matrices.

type

character string stating whether to plot the resulting statistics or the original fields along with their de-noised counter parts.

...

For wavePurifyVx, optional additional arguments to denoise.dwt.2d (or denoise.modwt.2d) from package waveslim. Note that if the argument J is not passed, then it will be determined as J=log2(

Value

• A list object of class wavePurifyVx is returned with possible components (depending on what is supplied in the arguments, etc.):
• X2,Y2
  m by n matrices of the denoised verification and forecast fields, resp. (only if return.fields is TRUE).
• thresholds
  q by 2 matrix of thresholds applied to the forecast (first column) and verification (second column) fields, resp. If climate is not NULL, then the same thresholds for the forecast field are applied to the climatology.
• qs
  If object and thresholds are NULL, and statistics other than MSE or ACC are desired, then this will be created along with the thresholds, and is just a character version of the trhesholds.
• args
  list object containing all the optional arguments passed into ..., and the value of J used (e.g., even if not passed into ...).
• bias,ts,ets,pod,far,f,hk,mse,acc
  numeric vectors of length q (i.e., the number of thresholds) giving the associated verification statistics.

Details

If the fields are dyadic, then the denoise.dwt.2d function from package waveslim is applied to each field before calculating the chosen verification statistics. Otherwise denoise.modwt.2d from the same package is used. The result is that high-frequency fluctuations in the two fields are removed before calculating verification statistics so that the resulting statistics are less susceptible to small-scale errors (see Briggs and Levine, 1997). See Percival and Guttorp (1994) and Lindsay et al. (1996) for more on this type of wavelet analysis including maximal overlap DWT.

References

Briggs, W. M. and Levine, R. A. (1997) Wavelets and field forecast verification. Mon. Wea. Rev., 125, 1329--1341.

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.

See Also

denoise.dwt.2d, denoise.modwt.2d, waverify2d, mowaverify2d, waveIS, vxstats

Examples

Run this code

grid<- list( x= seq( 0,5,,100), y= seq(0,5,,100))
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)

look3 <- sim.rf( obj) 
look3[ look3 < 0] <- 0 
look3 <- zapsmall( look3)

hold <- make.SpatialVx(look2, look2, thresholds=c(0.1,1),
            field.type="random", units="units")
plot(hold, set.pw=TRUE)
res <- wavePurifyVx(object=hold, climate=look3, return.fields=TRUE)
plot(res, type="fields")
plot(res, type="stats")
summary(res)

data(UKobs6)
data(UKfcst6)
data(UKloc)

hold <- surrogater2d(UKobs6, n=1, maxiter=50, verbose=TRUE)
hold <- matrix(hold, 256, 256)

UKobj <- make.SpatialVx(UKobs6, UKfcst6, thresholds=c(0.1,2,5,10), loc=UKloc,
    map=TRUE, field.type="Rainfall", units="mm/h",
    data.name=c("Nimrod", "obs 6", "fcst 6"))

plot(UKobj, set.pw=TRUE) 

look <- wavePurifyVx(object=UKobj, climate=hold, return.fields=TRUE, verbose=TRUE)
plot(look, type="fields", set.pw=TRUE)
plot(look, type="stats", set.pw=TRUE)
summary(look)

data(pert004)
data(pert000)
hold <- make.SpatialVx( pert000, pert004, thresholds=c(1,10,50), loc=ICPg240Locs,
    projection=TRUE, map=TRUE, field.type="Precipitation", units="mm/h",
    data.name=c("Perturbed ICP Cases", "pert000", "pert004"))
plot(hold, set.pw=TRUE)
look <- wavePurifyVx(object=hold, return.fields=TRUE, verbose=TRUE)
plot(look, type="fields", set.pw=TRUE)
plot(look, type="stats", set.pw=TRUE)
summary(look)

Run the code above in your browser using DataLab