SpatialVx (version 1.0-3)

variographier: Variography Score

Description

Calculate the variography score between two spatial fields based on the fitted exponential variogram.

Usage

variographier(x, init, zero.out = FALSE, ...)
# S3 method for default
variographier( x, init, zero.out = FALSE, ..., y )
# S3 method for SpatialVx
variographier( x, init, zero.out = FALSE, ...,
    obs = 1, model = 1, time.point = 1 )

Value

A list object of class “variographied” is returned with components:

obs.vg, mod.vg

Empirical variogram objects as returned by either vgram.matrix or variogram.matrix

obs.parvg, mod.parvg

objects returned by nlminb containing the fitted exponential variogram model parameters and some information about the optimization.

variography

single numeric giving the variography measure.

Arguments

x, y

matrices giving the fields on which to calculate the variography or a “SpatialVx” class object (x only).

init

list with components px and py that give initial values for parameter estimates (sill + nugget and range). If missing, default will attempt to find reasonable starting values.

zero.out

logical should the variogram be calculated over all grid points or just ones where one or both fields are non-zero? See variogram.matrix.

time.point

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

obs, model

numeric indicating which observation/forecast model to select for the analysis.

...

optional arguments to vgram.matrix or variogram.matrix (if zero.out is TRUE). Can also have optional arguments to nlminb (but not lower or upper).

Author

Eric Gilleland

Details

The variography score calculated here is that from Ekstrom (2016). So far, only the exponential variogram is allowed.

Note that in the fitting, the model g(h) = c * ( 1 - exp( -a * h ) ) is used, but the variography is calculated for theta = 3 / a. Therefore, the values in the par component of the returned fitted variograms correspond to a, while the variography score corresponds to theta. The score is given by:

v = 1 / sqrt( c_0^2 + c_m^2 + ( theta_0 - theta_m )^2 )

where c_0 and c_m are the sill + nugget terms for the observation and model, resp., and similarly for theta_0 and theta_m.

The parameters are *not* currently normalized, here, to give equal weight between sill + nugget and range. If several fields are analyzed (e.g., an ensemble), then the fitted parameters could be gathered, and one could use that information to calculate the score based on a normalized version.

References

Ekstrom, M. (2016) Metrics to identify meaningful downscaling skill in WRF simulations of intense rainfall events. Environmental Modelling and Software, 79, 267--284, DOI: 10.1016/j.envsoft.2016.01.012.

Examples

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

hold <- make.SpatialVx( UKobs6, UKfcst6, thresholds = c(0.01, 20.01),
    loc = UKloc, field.type = "Precipitation", units = "mm/h",
    data.name = "Nimrod", obs.name = "Observations 6", model.name = "Forecast 6",
    map = TRUE)

look <- variographier( hold )
look
plot( look )