Package 'SpatialVx'

Spatial Forecast Verification

Description

SpatialVx contains functions to perform many spatial forecast verification methods.

Details

Primary functions include:

0. make.SpatialVx: An object that contains the verification sets and pertinent information.

1. Filter Methods:

1a. Neighborhood Methods:

Neighborhood methods generally apply a convolution kernel smoother to one or both of the fields in the verificaiton set, and then apply the traditional scores. Most of the methods reviewed in Ebert (2008, 2009) are included in this package. The main functions are:

hoods2d, pphindcast2d, kernel2dsmooth, and plot.hoods2d.

1b. Scale Separation Methods:

Scale separation refers to the idea of applying a band-pass filter (and/or doing a multi-resolution analysis, MRA) to the verification set. Typically, skill is assessed on a scale-by-scale basis. However, other techniques are also applied. For example, denoising the field before applying traditional statistics, using the variogram, or applying a statistical test based on the variogram (these last are less similar to the spirit of the scale separation idea, but are at least somewhat related).

There is functionality to do the wavelet methods proposed in Briggs and Levine (1997). In particular, to simply denoise the fields before applying traditional verification statistics, use

wavePurifyVx.

To apply verification statistics to detail fields (in either the wavelet or field space), use:

waverify2d (dyadic fields) or mowaverify2d (non-dyadic or dyadic) fields.

The intensity-scale technique introduced in Casati et al. (2004) and the new developments of the technique proposed in Casati (2010) can be performed with

waveIS.

Although not strictly a “scale separation” method, the structure function (for which the variogram is a special case) is in the same spirit in the sense that it analyzes the field for different separation distances, and these “scales” are separate from each other (i.e., the score does not necessarily improve or decline as the scale increases). This package contains essentially wrapper functions to the vgram.matrix and plot.vgram.matrix functions from the fields package, but there is also a function called

variogram.matrix

that is a modification of vgram.matrix that allows for missing values. The primary function for doing this is called

griddedVgram, which has a plot method function associated with it.

There are also slight modifications of these functions (small modifications of the fields functions) to calculate the structure function of Harris et al. (2001). These functions are called

structurogram (for non-gridded fields) and structurogram.matrix (for gridded fields).

The latter allows for ignoring zero-valued grid points (as detailed in Harris et al., 2001) where the former does not (they must be removed prior to calling the function).

2. Displacement Methods:

In Gilleland et al (2009), this category was broken into two main types as field deformation and features-based. The former lumped together binary image measures/metrics with field deformation techniques because the binary image measures inform about the “similarity” (or dissimilarity) between the spatial extent or pattern of two fields (across the entire field). Here, they are broken down further into those that yield only a single (or small vector of) metric(s) or measure(s) (location measures), and those that have mechanisms for moving grid-point locations to match the fields better spatially (field deformation).

2a. Distance-based and Spatial-Alignment Summary Measures:

Gilleland (2020) introduced a new spatial alignment summary measure that falls between zero and one, with one representing a perfect match and zero a bad match. There is one user-selectable parameter/argument that determines the rate of decrease of the measure towards zero. Another two summary measures also incorporate intensity information. These summaries are available via Gbeta, GbetaIL and G2IL.

In addition to the above new measures, older well-known measures are inlcuded, including: the Hausdorff metric, partial-Hausdorff measure, FQI (Venugopal et al., 2005), Baddeley's delta metric (Baddeley, 1992; Gilleland, 2011; Schwedler et al., 2011), metrV (Zhu et al., 2011), as well as the localization performance measures described in Baddeley, 1992: mean error distance, mean square error distance, and Pratt's Figure of Merit (FOM).

locmeasures2d, metrV, distob, locperf

Image moments can give useful information about location errors, and are used within feature-based methods, particularly MODE, as they give the centroid of an image (or feature), as well as the orientation angle, among other useful properties. See the imomenter function for more details.

2b. Field deformation:

Thanks to Caren Marzban for supplying his optical flow code for this package (it has been modified some). These functions perform the analyses described in Marzban and Sandgathe (2010) and are based on the work of Lucas and Kanade (1981). See the help file for

OF.

Rigid transformations can be estimated using the rigider function. To simply rigidly transform a field (or feature) using specified parameters (x- and y- translations and/or rotations), the rigidTransform function can be used. For these functions, which may result in transformations that do not perfectly fall onto grid points, the function Fint2d can be used to interpolate from nearest grid points. Interpolation options include “round” (simply take the nearest location value), “bilinear” and “bicubic”.

2c. Features-based methods: These methods are also sometimes called object-based methods (the term “features” is used in this package in order to differentiate from an R object), and have many similarities to techniques used in Object-Based Image Analysis (OBIA), a relatively new research area that has emerged primarily as a result of advances in earth observations sensors and GIScience (Blaschke et al., 2008). It is attempted to identify individual features within a field, and subsequently analyze the fields on a feature-by-feature basis. This may involve intensity error information in addition to location-specific error information. Additionally, contingency table verifcation statistics can be found using new definitions for hits, misses and false alarms (correct negatives are more difficult to asses, but can also be done).

Currently, there is functionality for performing the analyses introduced in Davis et al. (2006,2009), including the merge/match algorithm of Gilleland et al (2008), as well as the SAL technique of Wernli et al (2008, 2009). Some functionality for composite analysis (Nachamkin, 2004) is provided by way of placing individual features onto a relative grid so that each shares the same centroid. Shape analysis is partially supported by way of functions to identify boundary points (Micheas et al. 2007; Lack et al. 2010). In particular, see:

Functions to identify features: FeatureFinder

Functions to match/merge features: centmatch, deltamm, minboundmatch

Functions to diagnose features and/or compare matched features:

FeatureAxis, FeatureComps, FeatureMatchAnalyzer, FeatureProps, FeatureTable, interester

See compositer for setting up composited objects, and see hiw (along with distill and summary method functions) for some shape analysis functionality.

The cluster analysis methods of Marzban and Sandgathe (2006; 2008) have been added. The former method was written from scratch by Eric Gilleland

clusterer

and the latter variation was modified from code originally written by Hilary Lyons

CSIsamples.

The Structure, Amplitude and Location (SAL) method can be performed with saller.

2d. Geometrical characterization measures:

Perhaps the measures in this sub-heading are best described as part-and-parcel of 2c. They are certainly useful in that domain, but have been proposed also for entire fields by AghaKouchak et al. (2011); though similar measures have been applied in, e.g., MODE. The measures introduced in AghaKouchak et al. (2011) available here are: connectivity index (Cindex), shape index (Sindex), and area index (Aindex):

Cindex, Sindex, Aindex

3. Statistical inferences for spatial (and/or spatiotemporal) fields:

In addition to the methods categorized in Gilleland et al. (2009), there are also functions for making comparisons between two spatial fields. The field significance approach detailed in Elmore et al. (2006), which requires a temporal dimension as well, involves using a circular block bootstrap (CBB) algorithm (usually for the mean error) at each grid point (or location) individually to determine grid-point significance (null hypothesis that the mean error is zero), and then a semi-parametric Monte Carlo method viz. Livezey and Chen (1983) to determine field significance.

spatbiasFS, LocSig, MCdof

In addition, the spatial prediction comparison test (SPCT) introduced by Hering and Genton (2011) is included via the functions: lossdiff, empiricalVG and flossdiff. Supporting functions for calculating the loss functions include: absolute error (abserrloss), square error (sqerrloss) and correlation skill (corrskill), as well as the distance map loss function (distmaploss) introduced in Gilleland (2013).

4. Other:

The bias corrected TS and ETS (or TS dHdA and ETS dHdA) introduced in Mesinger (2008) are now included within the vxstats function.

The 2-d Gaussian Mixture Model (GMM) approach introduced in Lakshmanan and Kain (2010) can be carried out using the

gmm2d

function (to estimate the GMM) and the associated summary function calculates the parameter comparisons. Also available are plot and predict method functions, but it can be very slow to run. The gmm2d employs an initialization function that takes the K largest object areas (connected components) and uses their centroids as initial estimates for the means, and uses the axes as initial guesses for the standard deviations. However, the user may supply their own initial estimate function.

The S1 score and anomaly correlation (ACC) are available through the functions

S1 and ACC.

See Brown et al. (2012) and Thompson and Carter (1972) for more on these statistics.

Also included is a function to do the geographic box-plot of Willmott et al. (2007). Namely,

GeoBoxPlot.

Datasets:

All of the initial Spatial Forecast Verification Inter-Comparison Project (ICP, https://projects.ral.ucar.edu/icp/) data sets used in the special collection of the Weather and Forecasting journal are included. See the help file for

obs0426,

which gives information on all of these datasets that are included, as well as two examples for plotting them: one that does not preserve projections, but plots the data without modification, and another that preserves the projections, but possibly with some interpolative smoothing.

Ebert (2008) provides a nice review of these methods. Roberts and Lean (2008) describes one of the methods, as well as the primary boxcar kernel smoothing method used throughout this package. Gilleland et al. (2009, 2010) provides an overview of most of the various recently proposed methods, and Ahijevych et al. (2009) describes the data sets included in this package. Some of these have been applied to the ICP test cases in Ebert (2009).

Additionally, one of the NIMROD cases (as provided by the UK Met Office) analyzed in Casati et al (2004) (case 6) is included along with approximate lon/lat locations. See the help file for UKobs6 more information.

A spatio-temporal verification dataset is also included for testing the method of Elmore et al. (2006). See the help file for

GFSNAMfcstEx.

A simulated dataset similar to the one used in Marzban and Sandgathe (2010) is also available and is called

hump.

Author(s)

Eric Gilleland

References

AghaKouchak, A., Nasrollahi, N. Li, J. Imam, B. and Sorooshian, S. (2011) Geometrical characterization of precipitation patterns. J. Hydrometeorology, 12, 274–285, doi:10.1175/2010JHM1298.1.

Ahijevych, D., Gilleland, E., Brown, B. G. and Ebert, E. E. (2009) Application of spatial verification methods to idealized and NWP gridded precipitation forecasts. Wea. Forecasting, 24 (6), 1485–1497.

Baddeley, A. J. (1992) An error metric for binary images. In Robust Computer Vision Algorithms, Forstner, W. and Ruwiedel, S. Eds., Wichmann, 59–78.

Blaschke, T., Lang, S. and Hay, G. (Eds.) (2008) Object-Based Image Analysis. Berlin, Germany: Springer-Verlag, 818 pp.

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

Brown, B. G., Gilleland, E. and Ebert, E. E. (2012) Chapter 6: Forecasts of spatial fields. pp. 95–117, In Forecast Verification: A Practitioner's Guide in Atmospheric Science, 2nd edition. Edts. Jolliffe, I. T. and Stephenson, D. B., Chichester, West Sussex, U.K.: Wiley, 274 pp.

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.

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.

Davis, C. A., Brown, B. G. and Bullock, R. G. (2006) Object-based verification of precipitation forecasts, Part I: Methodology and application to mesoscale rain areas. Mon. Wea. Rev., 134, 1772–1784.

Ebert, E. E. (2008) Fuzzy verification of high resolution gridded forecasts: A review and proposed framework. Meteorol. Appl., 15, 51–64. DOI: 10.1002/met.25

Ebert, E. E. (2009) Neighborhood verification: A strategy for rewarding close forecasts. Wea. Forecasting, 24, 1498–1510, doi:10.1175/2009WAF2222251.1.

Elmore, K. L., Baldwin, M. E. and Schultz, D. M. (2006) Field significance revisited: Spatial bias errors in forecasts as applied to the Eta model. Mon. Wea. Rev., 134, 519–531.

Gilleland, E., 2020. Novel measures for summarizing high-resolution forecast performance. Submitted to Advances in Statistical Climatology, Meteorology and Oceanography on 19 July 2020.

Gilleland, E. (2013) Testing competing precipitation forecasts accurately and efficiently: The spatial prediction comparison test. Mon. Wea. Rev., 141, (1), 340–355.

Gilleland, E. (2011) Spatial forecast verification: Baddeley's delta metric applied to the ICP test cases. Wea. Forecasting, 26, 409–415, doi:10.1175/WAF-D-10-05061.1.

Gilleland, E., Lee, T. C. M., Halley Gotway, J., Bullock, R. G. and Brown, B. G. (2008) Computationally efficient spatial forecast verification using Baddeley's delta image metric. Mon. Wea. Rev., 136, 1747–1757.

Gilleland, E., Ahijevych, D., Brown, B. G., Casati, B. and Ebert, E. E. (2009) Intercomparison of Spatial Forecast Verification Methods. Wea. Forecasting, 24, 1416–1430, doi:10.1175/2009WAF2222269.1.

Gilleland, E., Ahijevych, D. A., Brown, B. G. and Ebert, E. E. (2010) Verifying Forecasts Spatially. Bull. Amer. Meteor. Soc., October, 1365–1373.

Harris, D., Foufoula-Georgiou, E., Droegemeier, K. K. and Levit, J. J. (2001) Multiscale statistical properties of a high-resolution precipitation forecast. J. Hydrometeorol., 2, 406–418.

Hering, A. S. and Genton, M. G. (2011) Comparing spatial predictions. Technometrics 53, (4), 414–425.

Lack, S., Limpert, G. L. and Fox, N. I. (2010) An object-oriented multiscale verification scheme. Wea. Forecasting, 25, 79–92, DOI: 10.1175/2009WAF2222245.1

Lakshmanan, V. and Kain, J. S. (2010) A Gaussian Mixture Model Approach to Forecast Verification. Wea. Forecasting, 25 (3), 908–920.

Livezey, R. E. and Chen, W. Y. (1983) Statistical field significance and its determination by Monte Carlo techniques. Mon. Wea. Rev., 111, 46–59.

Lucas, B D. and Kanade, T. (1981) An iterative image registration technique with an application to stereo vision. Proc. Imaging Understanding Workshop, DARPA, 121–130.

Marzban, C. and Sandgathe, S. (2006) Cluster analysis for verification of precipitation fields. Wea. Forecasting, 21, 824–838.

Marzban, C. and Sandgathe, S. (2008) Cluster Analysis for Object-Oriented Verification of Fields: A Variation. Mon. Wea. Rev., 136, (3), 1013–1025.

Marzban, C. and Sandgathe, S. (2009) Verification with variograms. Wea. Forecasting, 24 (4), 1102–1120, doi: 10.1175/2009WAF2222122.1

Marzban, C. and Sandgathe, S. (2010) Optical flow for verification. Wea. Forecasting, 25, 1479–1494, doi:10.1175/2010WAF2222351.1.

Mesinger, F. (2008) Bias adjusted precipitation threat scores. Adv. Geosci., 16, 137–142.

Micheas, A. C., Fox, N. I., Lack, S. A. and Wikle, C. K. (2007) Cell identification and verification of QPF ensembles using shape analysis techniques. J. of Hydrology, 343, 105–116.

Nachamkin, J. E. (2004) Mesoscale verification using meteorological composites. Mon. Wea. Rev., 132, 941–955.

Roberts, N. M. and Lean, H. W. (2008) Scale-selective verification of rainfall accumulations from high-resolution forecasts of convective events. Mon. Wea. Rev., 136, 78–97. doi:10.1175/2007MWR2123.1.

Schwedler, B. R. J. and Baldwin, M. E. (2011) Diagnosing the sensitivity of binary image measures to bias, location, and event frequency within a forecast verification framework. Wea. Forecasting, 26, 1032–1044, doi:10.1175/WAF-D-11-00032.1.

Thompson, J. C. and Carter, G. M. (1972) On some characteristics of the S1 score. J. Appl. Meteorol., 11, 1384–1385.

Venugopal, V., Basu, S. and Foufoula-Georgiou, E. (2005) A new metric for comparing precipitation patterns with an application to ensemble forecasts. J. Geophys. Res., 110, D08111, doi:10.1029/2004JD005395, 11pp.

Wernli, H., Paulat, M., Hagen, M. and Frei, C. (2008) SAL–A novel quality measure for the verification of quantitative precipitation forecasts. Mon. Wea. Rev., 136, 4470–4487.

Wernli, H., Hofmann, C. and Zimmer, M. (2009) Spatial forecast verification methods intercomparison project: Application of the SAL technique. Wea. Forecasting, 24, 1472–1484, doi:10.1175/2009WAF2222271.1

Willmott, C. J., Robeson, S. M. and Matsuura, K. (2007) Geographic box plots. Physical Geography, 28, 331–344, DOI: 10.2747/0272-3646.28.4.331.

Zhu, M., Lakshmanan, V. Zhang, P. Hong, Y. Cheng, K. and Chen, S. (2011) Spatial verification using a true metric. Atmos. Res., 102, 408–419, doi:10.1016/j.atmosres.2011.09.004.

Examples

## See help files for above named functions and datasets
## for specific examples.

---

Loss functions for the spatial prediction comparison test (SPCT)

Description

Loss functions for applying the spatial prediction comparison test (SPCT) for competing forecasts.

Usage

abserrloss(x, y, ...)
corrskill(x, y, ...)
sqerrloss(x, y, ...)
distmaploss(x, y, threshold = 0, const = Inf, ...)

Arguments

x, y	m by n numeric matrices against which to calculate the loss (or skill) functions.
threshold	numeric giving the threshold over which (and including) binary fields are created from x and y in order to make a distance map.
const	numeric giving the constant beyond which the differences in distance maps between x and y are set to zero. If Inf (default), then no cut-off is taken. The SPCT is probably not powerful for large values of const.
...	Not used by abserrloss or sqerrloss (there for consistency only, and in order to work with lossdiff). For corrskill, these are optional arguments to sd. For distmaploss, these are optional arguments to the distmap function from pacakge spatstat.

Details

These are simple loss functions that can be used in conjunction with lossdiff to carry out the spatial prediction comparison test (SPCT) as introduced in Hering and Genton (2011); see also Gilleland (2013) in particular for details about the distance map loss function.

The distance map loss function does not zero-out well as the other loss functions do. Therefore, zero.out should be FALSE in the call to lossdiff. Further, as pointed out in Gilleland (2013), the distance map loss function can easily be hedged by having a lot of correct negatives. The image warp loss function is probably better for this purpose if, e.g., there are numerous zero-valued grid points in all fields.

Value

numeric m by n matrices containing the value of the loss (or skill) function at each location i of the original set of locations (or grid of points).

Author(s)

Eric Gilleland

References

Gilleland, E. (2013) Testing competing precipitation forecasts accurately and efficiently: The spatial prediction comparison test. Mon. Wea. Rev., 141, (1), 340–355.

Hering, A. S. and Genton, M. G. (2011) Comparing spatial predictions. Technometrics 53, (4), 414–425.

Examples

# See help file for lossdiff for examples.

---

Area Index

Description

Calculate Area index described in AghaKouchak et al. (2011).

Usage

Aindex(x, thresh = NULL, dx = 1, dy = 1, ...)

## Default S3 method:
Aindex(x, thresh = NULL, dx = 1, dy = 1, ...)

## S3 method for class 'SpatialVx'
Aindex(x, thresh = NULL, dx = 1, dy = 1, ...,
    time.point=1, obs = 1, model=1)

Arguments

x	Default: m by n numeric matrix giving the field for which the area index is to be calculated. Aindex.SpatialVx: list object of class “SpatialVx”.
thresh	Values under this threshold are set to zero. If NULL, it will be set to 1e-8 (a very small value).
dx, dy	numeric giving the grid point size in each direction if it is desired to apply such a correction. However, the values are simply canceled out in the index, so these arguments are probably not necessary. If it is desired to only get the area of the non-zero values in the field, or the convex hull, then these make sense.
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.
...	Not used.

Details

The area index introduced in AghaKouchak et al. (2011) is given by

Aindex = A/Aconvex,

where A is the area of the pattern, and Aconvex the area of its convex hull (area.owin from package spatstat is used to calculate this latter area, and the functions as.im and solutionset from spatstat are also used by this function). Values are between 0 and 1. Values closer to unity indicate a more structured pattern, and values closer to zero indicate higher dispersiveness of the pattern, but note that two highly structured patterns far away from each other may also give a low value (see examples below). Because of this property, this measure is perhaps best applied to individual features in a field.

Value

numeric vector (or two-row matrix in the case of Aindex.SpatialVx) with named components (columns):

Aindex	numeric giving the area index.
A, Aconvex	numeric giving the area of th epattern and the convex hull, resp.
dx, dy	the values of dx and dy as input to the function.

Author(s)

Eric Gilleland

References

AghaKouchak, A., Nasrollahi, N., Li, J., Imam, B. and Sorooshian, S. (2011) Geometrical characterization of precipitation patterns. J. Hydrometeorology, 12, 274–285, doi:10.1175/2010JHM1298.1.

See Also

Cindex, Sindex

Examples

# Gemetric shape that is highly structured.
# Re-create Fig. 7a from AghaKouchak et al. (2011).
tmp <- matrix(0, 8, 8)
tmp[3,2:4] <- 1
tmp[5,4:6] <- 1
tmp[7,6:7] <- 1
Aindex(tmp)

---

Bearing from One Spatial Location to Another

Description

Find the bearing from one spatial location to another.

Usage

bearing(point1, point2, deg = TRUE, aty = "compass")

Arguments

point1, point2	two-column numeric matrices giving lon/lat coordinates for the origin point(s) (point1) and the destination point(s) (point2).
deg	logical, should the output be converted from radians to degrees?
aty	character stating either “compass” (default) or “radial”. The former gives the standard compass bearing angle (0 is north, increase clockwise), and the latter is for polar coordinates (0 is East, increase counter-clockwise).

Details

The bearing, beta, of a point B as seen from a point A is given by

beta = atan2(S,T)

where

S = cos(phi_B) * sin(L_A - L_B), and

T = cos(phi_A)*sin(phi_B) - sin(phi_A)*cos(phi_B)*cos(L_A - L_B)

where phi_A (phi_B) is the latitude of point A (B), and L_A (L_B) is the longitude of point A (B).

Note that there is no simple relationship between the bearing of A to B vs. the bearing of B to A. The bearing given here is in the usual R convention for lon/lat information, which gives points east of Greenwich as negative longitude, and south of the equator as negative latitude.

Value

numeric giving the bearing angle.

Author(s)

Eric Gilleland and Randy Bullock, bullock “at” ucar.edu

References

Keay, W. (1995) Land Navigation: Routefinding with Map & Compass, Coventry, UK: Clifford Press Ltd., ISBN 0319008452, 978-0319008454

See Also

atan2, FeatureAxis, fields::rdist.earth

Examples

# Boulder, Colorado and Wallaroo, Australia.
A <- rbind(c(-105.2833, 40.0167), c(137.65, -33.9333))

# Wallaroo, Australia and Boulder, Colorado.
B <- rbind(c(137.65, -33.9333), c(-105.2833, 40.0167))

bearing(A,B)
bearing(A,B,aty="radial")

plot(A, type="n", xlab="", ylab="")
points(A[,1], A[,2], pch="*", col="darkblue")

# Boulder, Colorado to Wallaroo, Australia.
arrows(A[1,1], A[1,2], A[2,1], A[2,2], col="red", lwd=1.5)

---

Create Binary Fields

Description

Convert a spatial field to a binary spatial field via thresholding.

Usage

binarizer(X, Xhat, threshold = NULL,
	  rule = c(">", "<", ">=", "<=", "<>", "><", "=<>",
		   "<>=", "=<>=", "=><", "><=", "=><="),
	  value = c("matrix", "owin"), ...)

Arguments

X, Xhat	matrix or “owin” objects.
threshold	single number, numeric vector of two numbers, or two-by-two matrix; depending on the value of rule. May be missing or null if both X and Xhat are “owin” class objects, in which case binary fields are made wherever these fields are greater than zero, and other arguments are ignored.
rule	character giving the rule for identifying 1-valued grid squares. For example, if rule is the default (“>”), then threshold should either be a single numeric or a vector of length two. If the latter, it specifies a different threshold for X and Xhat, the result is that everywhere X > threshold will have a value of 1 and zero otherwise, etc. The rule “<>” means 1-values whenever X (Xhat) are less than the lower threshold and higher than the higher threshold, etc. For rules requiring two threshold values, if a two-by-two matrix is given, the first column is associated with X and the second with Xhat.
value	character telling whether the returned object be a list with two matrices or a list with two “owin” class objects.
...	Not used.

Details

The binary fields are created by assigning ones according to the rule: 1. ">": if X > threshold, assign 1, zero otherwise. 2. "<": if X < threshold, assign 1, zero otherwise. 3. ">=": if X >= threshold, assign 1, zero otherwise. 4. "<=": if X <= threshold, assign 1, zero otherwise. 5. "<>": if X < threshold[ 1 ] or X > threshold[ 2 ], assign 1, zero otherwise. 6. "><": if threshold[ 1 ] < X < threshold[ 2 ], assign 1, zero otherwise. 7. "=<>": if threshold[ 1 ] <= X or X > threshold[ 2 ], assign 1, zero otherwise. 8. "<>=": if threshold[ 1 ] < X or X >= threshold[ 2 ], assign 1, zero otherwise. 9. "=<>=": if X <= threshold[ 1 ] or X >= threshold[ 2 ], assign 1, zero otherwise. 10. "=><": if threshold[ 1 ] <= X < threshold[ 2 ], assign 1, zero otherwise. 11. "><=": if threshold[ 1 ] < X <= threshold[ 2 ], assign 1, zero otherwise. 12. "=><=": if threshold[ 1 ] <= X <= threshold[ 2 ], assign 1, zero otherwise.

Value

A list object with two components is returned with the first component being the binary version of X and the second that for Xhat. These fields will either be matrices of the same dimension as X and Xhat or they will be owin objects depending on the value argument.

Author(s)

Eric Gilleland

Examples

data( "obs0601" )
data( "wrf4ncar0531" )
bin <- binarizer( X = obs0601, Xhat = wrf4ncar0531, threshold = 2.1 )

image.plot( bin[[ 1 ]] )
image.plot( bin[[ 2 ]] )

bin2 <- binarizer( X = obs0601, Xhat = wrf4ncar0531,
		  threshold = 2.1, value = "owin" )
plot( bin2[[ 1 ]] )
plot( bin2[[ 2 ]] )

---

Calculate dFSS Score Value

Description

Calculates the value for the dFSS binary distance metric. The dFSS uses the Fraction Skill Score to provide a measure of spatial displacement of precipitation in two precipitation fields.

Usage

calculate_dFSS(fbin1, fbin2)

Arguments

fbin1	A numeric matrix representing the first binary field. Only values 0 and 1 are allowed in the matrix.
fbin2	A numeric matrix representing the second binary field. Only values 0 and 1 are allowed in the matrix. The matrix needs to have the same dimensions as fbin1.

Details

The dFSS uses the Fraction Skill Score to provide a measure of spatial displacement of precipitation in two precipitation fields.

The function requires two binary fields as input. A binary field can only have values of 0 or 1 and can be obtained through a thresholding process of the original continuous precipitation field (e.g., by setting all values below a selected precipitation threshold to zero, and all values above the threshold to one).

The dFSS has a requirement that the frequency bias of precipitation needs to be small in order for the metric to work properly (i.e. the number of non-zero grid points has to be similar in both binary fields). The unbiased fields can be obtained from the original continuous precipitation fields via the use of a frequency (percentile) threshold. For example, instead of using a predefined physical threshold (e.g. 1 mm/h), which might produce binary fields with a different number of non-zero points, a frequency threshold (e.g. 5 %) can be used which guarantees that both fields will have the same number of non-zero grid-points and will thus be unbiased (provided that enough grid points in the domain contain non-zero precipitation). Function quantile can be used to determine the value of a physical threshold that corresponds to a prescribed frequency threshold.

If the frequency bias is larger than 1.5 the function will work but produce a warning. If the frequency bias is larger than 2 the function will produce an error. The dFSS value can only be calculated if both fields contain at least one non-zero grid point. For correct interpretation of the results and some other considerations please look at the "Recipe" in the Conclusions section of Skok and Roberts (2018).

The code utilizes the fast method for computing fractions (Faggian et al., 2015) and the Bisection method to arrive more quickly at the correct displacement. Optionally, a significantly faster R code that requires significantly less memory and uses some embedded C++ code is available upon request from the author.

Value

The function returns a single numeric value representing the size of the estimated spatial displacement (expressed as a number of grid points - see the example below).

Author(s)

Gregor Skok ([email protected])

References

Skok, G. and Roberts, N. (2018), Estimating the displacement in precipitation forecasts using the Fractions Skill Score. Q.J.R. Meteorol. Soc. doi:10.1002/qj.3212.

Faggian N., Roux B., Steinle P., Ebert B., 2015: Fast calculation of the Fractions Skill Score, MAUSAM, 66 (3), 457-466.

See Also

hoods2d

Examples

# ---------------------------------------------
# A simple example with two 500 x 500 fields
# ---------------------------------------------

# generate two empty 500 x 500 fields where all values are 0 
fbin1=matrix(0, 500, 500, byrow = FALSE)
fbin2=fbin1

# in the fields define a single 20x20 non-zero region of precipitation
# that is horizontally displaced in the second field by 100 grid points   
fbin1[200:220,200:220]=1
fbin2[200:220,300:320]=1

# calulate dFSS value
dFSS=calculate_dFSS(fbin1, fbin2)

# print dFSS value 
print(dFSS)

# The example should output 97 which means that the spatial displacement
# estimated by dFSS is 97 grid points.

---

Calculate FSS Values

Description

Calculates the value of Fraction Skill Score (FSS) for multiple neighborhood sizes.

Usage

calculate_FSSvector_from_binary_fields(fbin1, fbin2, nvector)

Arguments

fbin1	A numeric matrix representing the first binary field. Only values 0 and 1 are allowed in the matrix.
fbin2	A numeric matrix representing the second binary field. Only values 0 and 1 are allowed in the matrix. The matrix needs to have the same dimensions as fbin1.
nvector	A numeric vector containing neighborhood sizes for which the FSS values are to be calculated. Only positive odd values are allowed in the vector. A square neighborhood shape is assumed and the specified value represents the length of square side.

Details

Fractions Skill Score is a neighborhood-based spatial verification metric frequently used for verifying precipitation (see Roberts and Lean, 2008, for details).

The function requires two binary fields as input. A binary field can only have values of 0 or 1 and can be obtained through a thresholding process of the original continuous precipitation field (e.g., by setting all values below a selected precipitation threshold to zero, and all values above the threshold to one). Either a predefined physical threshold (e.g. 1 mm/h) or a frequency threshold (e.g. 5 %) can be used to produce the binary fields from the original continuous precipitation fields. If a frequency threshold is used the binary fields will be unbiased and the FSS value will asymptote to 1 at large neighborhoods. Function quantile can be used to determine the value of a physical threshold that corresponds to a prescribed frequency threshold.

The code utilizes the fast method for computing fractions (Faggian et al., 2015) that enables fast computation of FSS values at multiple neighborhood sizes. Optionally, a significantly faster R code that requires significantly less memory and uses some embedded C++ code is available upon request from the author.

Value

A numeric vector of the same dimension as nvector that contains the FSS values at corresponding neighborhood sizes.

Author(s)

Gregor Skok ([email protected])

References

Roberts, N.M., Lean, H.W., 2008. Scale-Selective Verification of Rainfall Accumulations from High-Resolution Forecasts of Convective Events. Mon. Wea. Rev. 136, 78-97.

Faggian N., Roux B., Steinle P., Ebert B., 2015: Fast calculation of the Fractions Skill Score, MAUSAM, 66 (3), 457-466.

See Also

hoods2d

Examples

# ---------------------------------------------
# A simple example with two 500 x 500 fields
# ---------------------------------------------

# generate two empty 500 x 500 binary fields where all values are 0 
fbin1=matrix(0, 500, 500, byrow = FALSE)
fbin2=fbin1

# in the fields define a single 20x20 non-zero region of precipitation that
# is horizontally displaced in the second field by 100 grid points 
fbin1[200:220,200:220]=1
fbin2[200:220,300:320]=1

# specify a vector of neighborhood sizes for which the FSS values are to be calculated
nvector = c(1,51,101,201,301,601,901,1501)

# calulate FSS values
FSSvector=calculate_FSSvector_from_binary_fields(fbin1, fbin2, nvector)

# print FSS values 
print(FSSvector)

# The example should output:
#  0.00000000 0.00000000 0.04271484 0.52057596 0.68363656 0.99432823 1.00000000 1.00000000

---

Calculate FSSwind Score Value

Description

Calculates the value for the FSSwind metric that can be used for spatial verification of 2D wind fields.

Usage

calculate_FSSwind(findex1, findex2, nvector)

Arguments

findex1	A numeric matrix representing the first wind class index field. Only integer values larger than 0 are allowed in the matrix. Each integer value corresponds to a certain wind class.
findex2	A numeric matrix representing the second wind class index field. Only integer values larger than 0 are allowed in the matrix. Each integer value corresponds to a certain wind class. The matrix needs to have the same dimensions as findex1.
nvector	A numeric vector containing neighborhood sizes for which the FSSwind score value is to be calculated. Only positive odd values are allowed in the vector.

Details

The FSSwind is based on the idea of the Fractions Skill Score, a neighborhood-based spatial verification metric frequently used for verifying precipitation. The FSSwind avoids some of the problems of traditional non-spatial verification metrics (the "double penalty" problem and the failure to distinguish between a "near miss" and much poorer forecasts) and can distinguish forecasts even when the spatial displacement of wind patterns is large. Moreover, the time-averaged score value in combination with a statistical significance test enables different wind forecasts to be ranked by their performance (see Skok and Hladnik, 2018, for details).

The score can be used to spatially compare two 2D wind vector fields. In order to calculate the FSSwind value, wind classes have first to be defined. The choice of classes will define how the score behaves and influence the results. The definition of classes should reflect what a user wants to verify. The score will evaluate the spatial matching of the areas of the wind classes. The class definitions should cover the whole phase space of possible wind values (i.e. to make sure every wind vector can be assigned a wind class). It does not make sense choosing a class definition with an overly large number of classes or a definition where a single class would totally dominate over all the other classes. Once the class definition is chosen, the wind vector fields need to be converted to wind class index fields with each class assigned an unique integer value. These wind class index fields serve as input to the calculate_FSSwind function. The FSSwind can have values between 0 and 1 with 0 indicating the worst possible forecast and 1 indicating a perfect forecast. For guidance on correct interpretation of the results and other important considerations please refer to Skok and Hladnik (2018).

The code utilizes the fast method for computing fractions (Faggian et al., 2015). Optionally, a significantly faster R code that requires significantly less memory and uses some embedded C++ code is available upon request from the author.

Value

A numeric vector of the same dimension as nvector that contains the FSSwind values at corresponding neighborhood sizes.

Author(s)

Gregor Skok ([email protected])

References

Skok, G. and V. Hladnik, 2018: Verification of Gridded Wind Forecasts in Complex Alpine Terrain: A New Wind Verification Methodology Based on the Neighborhood Approach. Mon. Wea. Rev., 146, 63-75, https://doi.org/10.1175/MWR-D-16-0471.1

Faggian N., Roux B., Steinle P., Ebert B., 2015: Fast calculation of the Fractions Skill Score, MAUSAM, 66 (3), 457-466.

See Also

hoods2d

Examples

# ---------------------------------------------
# A simple example with two 500 x 500 fields
# ---------------------------------------------

# generate two 500 x 500 wind class index fields where all values are 1 (wind class 1)
findex1=matrix(1, 500, 500, byrow = FALSE)
findex2=findex1

# in the fields generate some rectangular areas with other wind classes (classes 2,3 and 4)
findex1[001:220,200:220]=2
findex1[100:220,300:220]=3
findex1[300:500,100:200]=4

findex2[050:220,100:220]=2
findex2[200:320,300:220]=3
findex2[300:500,300:500]=4

# specify a vector of neighborhood sizes for which the FSSwind values are to be calculated
nvector = c(1,51,101,201,301,601,901,1501)

# calulate FSSwind values
FSSwindvector=calculate_FSSwind(findex1, findex2, nvector)

# print FSSwind values 
print(FSSwindvector)

# The example should output:
# 0.6199600 0.6580598 0.7056385 0.8029494 0.8838075 0.9700274 0.9754587 0.9756134

---

Centered on Square Domain Baddeley's Delta Metric

Description

Baddeley's delta metric is sensitive to the position of non-zero grid points within the domain, as well as to the size of the domain. In order to obtain consistent values of the metric across cases, it is recommended to first position the sets to be compared so that they are centered with respect to one another on a square domain; where the square domain is the same for all comparison sets.

Usage

censqdelta(x, y, N, const = Inf, p = 2, ...)

Arguments

x, y	Matrices representing binary images to be compared. If they are not binary, then they will be forced to binary by setting anything above zero to one.
N	The size of the square domain. If missing, it will be the size of the largest side, and if it is even, one will be added to it.
const	single numeric giving the c argument to deltametric, which is the constant value over which the distance map is reduced to this value.
p	single numeric giving the p argument to deltametric, which specifies the type of Lp norm used to calculate the delta maetric.
...	Not used.

Details

Baddeley's delta metric (Baddeley, 1992a,b) is the L_p norm over the absolute difference of distance maps for two binary images, A and B. A concave function (e.g., f(t) = min(t, constant)) may first be applied ot each distance map before taking their absolute differences, which makes the result less sensitive to small changes in one or both images than other similar metrics. The metric is sensitive to size, shape and location differences, which make it very practical for comparing forecasts to observations in terms of the position, area extent, and area shape errors. However, its sensitivity to domain size and position within the domain are undesirable, but are easily fixed by calculating the metric over a consistent, square domain with the combined verification set centerd on that domain. See the example section below to see the issue.

This function essentially takes a window of size N by N and moves so that the centroid of each pair of sets, A and B, is the center of the window before calculating the metric.

Centering and squaring is recommended for carrying out a procedure such as that proposed in Gilleland et al. (2008). Centering on a square domain alleviates the problems discovered by Schwedler and Baldwin (2011) who suggested using a small value of the constant in f(t) = min(t, constant) applied to the distance maps. This solution is not very appealing because of the sensitivity in choice of the constant that generally diminishes as it approaches the domain size (Gilleland, 2011).

After centering the sets on a square domain, the function deltametric from package spatstat is used to calculate the metric.

Value

A single numeric value is returned.

Author(s)

Eric Gilleland

References

Baddeley, A. (1992a) An error metric for binary images. In Robust Computer Vision Algorithms, W. Forstner and S. Ruwiedel, Eds., Wichmann, 59–78.

Baddeley, A. (1992b) Errors in binary images and an Lp version of the Hausdorff metric. Nieuw Arch. Wiskunde, 10, 157–183.

Gilleland, E. (2011) Spatial Forecast Verification: Baddeley's Delta Metric Applied to the ICP Test Cases. Weather Forecast., 26 (3), 409–415.

Gilleland, E. (2017) A new characterization in the spatial verification framework for false alarms, misses, and overall patterns. Weather Forecast., 32 (1), 187–198, DOI: 10.1175/WAF-D-16-0134.1.

Gilleland, E., Lee, T. C. M., Halley Gotway, J., Bullock, R. G. and Brown, B. G. (2008) Computationally efficient spatial forecast verification using Baddeley's delta image metric. Mon. Wea. Rev., 136, 1747–1757.

Schwedler, B. R. J. and Baldwin, M. E. (2011) Diagnosing the sensitivity of binary image measures to bias, location, and event frequency within a forecast verification framework. Weather Forecast., 26, 1032–1044.

See Also

locmeasures2d

Examples

x <- y <- matrix( 0, 100, 200 )
x[ 45, 10 ] <- 1
x <- kernel2dsmooth( x, kernel.type = "disk", r = 4 )

y[ 50, 60 ] <- 1
y <- kernel2dsmooth( y, kernel.type = "disk", r = 10 )

censqdelta( x, y )

## Not run: 
# Example form Gilleland (2017).
#
# I1 = circle with radius = 20 centered at 100, 100
# I2 = circle with radius = 20 centered at 140, 100
# I3 = circle with radius = 20 centered at 180, 100
# I4 = circle with radius = 20 centered at 140, 140

I1 <- I2 <- I3 <- I4 <- matrix( 0, 200, 200 )

I1[ 100, 100 ] <- 1
I1 <- kernel2dsmooth( I1, kernel.type = "disk", r = 20 )
I1[ I1 > 0 ] <- 1
if( any( I1 < 0 ) ) I1[ I1 < 0 ] <- 0

I2[ 140, 100 ] <- 1
I2 <- kernel2dsmooth( I2, kernel.type = "disk", r = 20 )
I2[ I2 > 0 ] <- 1
if( any( I2 < 0 ) ) I2[ I2 < 0 ] <- 0

I3[ 180, 100 ] <- 1
I3 <- kernel2dsmooth( I3, kernel.type = "disk", r = 20 )
I3[ I3 > 0 ] <- 1
if( any( I3 < 0 ) ) I3[ I3 < 0 ] <- 0

I4[ 140, 140 ] <- 1
I4 <- kernel2dsmooth( I4, kernel.type = "disk", r = 20 )
I4[ I4 > 0 ] <- 1
if( any( I4 < 0 ) ) I4[ I4 < 0 ] <- 0

image( I1, col = c("white", "darkblue") )
contour( I2, add = TRUE )
contour( I3, add = TRUE )
contour( I4, add = TRUE )

# Each circle is the same size and shape, and the domain is square.
# I1 and I2, I2 and I3, and I2 and I4 are all the same distance
# away from each other.  I1 and I4 and I3 and I4 are also the same distance
# from each other.  I3 touches the edge of the domain.
# 

# First, calculate the Baddeley delta metric on each 
# comparison.

I1im <- as.im( I1 )
I2im <- as.im( I2 )
I3im <- as.im( I3 )
I4im <- as.im( I4 )

I1im <- solutionset( I1im > 0 )
I2im <- solutionset( I2im > 0 )
I3im <- solutionset( I3im > 0 )
I4im <- solutionset( I4im > 0 )

deltametric( I1im, I2im )
deltametric( I2im, I3im )
deltametric( I2im, I4im )

# Above are all different values.
# Below, they are all 28.84478.
censqdelta( I1, I2 )
censqdelta( I2, I3 )
censqdelta( I2, I4 )

# Similarly for I1 and I4 vs I3 and I4.
deltametric( I1im, I4im )
deltametric( I3im, I4im )

censqdelta( I1, I4 )
censqdelta( I3, I4 )

# To see why this problem exists.
dm1 <- distmap( I1im )
dm1 <- as.matrix( dm1 )
dm2 <- distmap( I2im )
dm2 <- as.matrix( dm2 )

par( mfrow = c( 2, 2 ) )
image.plot( dm1 )
contour( I1, add = TRUE, col = "white" )
image.plot( dm2 )
contour( I2, add = TRUE, col = "white" )

image.plot( abs( dm1 ) - abs( dm2 ) )
contour( I1, add = TRUE, col = "white" )
contour( I2, add = TRUE, col = "white" )


## End(Not run)

---

Centroid Distance Between Two Identified Objects

Description

Find the centroid distance between two identified objects/features.

Usage

centdist(x, y, distfun = "rdist", loc = NULL, ...)

Arguments

x, y	objects of class “owin” (package spatstat) containing binary images of features of interest.
distfun	character string naming a distance function that should take arguments x1 and x2 as 1 by 2 matrices, and return a single numeric value. Default uses the fields function, rdist, where the fields function rdist.earth is an obvious alternative.
loc	two-column matrix giving the location values for which to calculate the centroids. If NULL, indices according to the dimension of the fields are used.
...	optional arguments to distfun.

Details

This is a simple function that calculates the centroid for each of x and y (to get their centroids), and then finds the distance between them according to distfun. The centroids are calculated using FeatureProps.

Value

numeric giving the centroid distance.

Author(s)

Eric Gilleland

Examples

x <- y <- matrix(0, 10, 12)
x[2:3,c(3:6, 8:10)] <- 1
y[c(4:7, 9:10),c(7:9, 11:12)] <- 1

x <- as.im(x)
x <- solutionset(x>0)
y <- as.im(y)
y <- solutionset(y>0)
centdist(x,y)

---

Connectivity Index

Description

Calculate the connectivity index of an image.

Usage

Cindex(x, thresh = NULL, connect.method = "C", ...)

## Default S3 method:
Cindex(x, thresh = NULL, connect.method = "C", ...)

## S3 method for class 'SpatialVx'
Cindex(x, thresh = NULL, connect.method = "C", ...,
    time.point = 1, obs = 1, model = 1)

Arguments

x	Default: m by n numeric matrix giving the field for which the connectivity index is to be calculated. Sindex.SpatialVx: list object of class “SpatialVx”.
thresh	Set values under (strictly less than) this threshold to zero, and calculate the connectivity index for the resulting image. If NULL, no threshold is applied.
connect.method	character string giving the method argument for the connected function of package spatstat. This must be one of “C” or “interpreted”. See the help file for connected for more details.
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.
...	Not used.

Details

The connectivity index is introduced in AghaKouchak et al. (2011), and is designed to automatically determine how connected an image is. It is defined by

Cindex = 1 - (NC - 1)/(sqrt(NP) + NC),

where 0 <= Cindex <= 1 is the connectivity index (values close to zero are less connected, and values close to 1 are more connected), NP is the number of nonzero pixels, and NC is the number of isolated clusters.

The function connected from package spatstat is used to identify the number of isolated clusters.

Value

numeric giving the connectivity index.

Author(s)

Eric Gilleland

References

AghaKouchak, A., Nasrollahi, N., Li, J., Imam, B. and Sorooshian, S. (2011) Geometrical characterization of precipitation patterns. J. Hydrometerology, 12, 274–285, doi:10.1175/2010JHM1298.1.

See Also

Sindex, Aindex

Examples

# Re-create Fig. 7a from AghaKouchak et al. (2011).
tmp <- matrix(0, 8, 8)
tmp[3,2:4] <- 1
tmp[5,4:6] <- 1
tmp[7,6:7] <- 1
Cindex(tmp)

---

Circular Histogram

Description

Make a circle histogram, also known as a wind-rose diagram.

Usage

CircleHistogram(wspd, wdir, numPetals = 12, radians = FALSE,
		COLS = NULL, scale.factor = 3, varwidth = TRUE,
		minW = NULL, maxW = NULL, circFr = 10,
		main = "Wind Rose", cir.ind = 0.05,
		max.perc = NULL, leg = FALSE, units = "units",
		verbose = FALSE, ...)

Arguments

wspd	numeric vector of length n giving the magnitudes for the histograms (i.e., the petal colors). For example, with MODE output, one might use the centroid distances between the paired objects, or the Baddeley Delta metric, or whatever other attribute is of interest. The argument is called wspd because in the original context of windrose graphs, the windspeed was used for the petal colors.
wdir	numeric vector of length n giving the angles (from the north) between two matched objects. In the original context of a windrose diagram, the direction from the north of the wind. See bearing for a way to obtain such vectors.
numPetals	numeric giving the number of petals to use.
radians	logical if TRUE the angles displayed are radians, if FALSE degrees.
COLS	vector defining the colors to be used for the petals. See, for example, topo.colors.
scale.factor	numeric determining the line widths (scaled against the bin sizes), only used if varwidth is TRUE.
varwidth	logical determining whether to vary the widths of the petals or not.
minW, maxW	single numerics giving the minimum and maximum break ranges for the histogram of each petal. If NULL, it will be computed as min( wspd) and max( wspd).
circFr	numeric giving the bin width.
main	character string giving the title to add to the plot.
cir.ind	numeric only used if max.perc is NULL.
max.perc	numeric giving the maximum percentage to show
leg	logical determining whether or not to add a legend to the plot.
units	character string giving the units for use with the legend. Not used if leg is FALSE.
verbose	logical telling whether or not to print information to the screen.
...	optional arguments to the plot or text functions.

Details

The windrose diagram, or circle histogram, is similar to a regular histogram, but adds more information when direction is important. Binned directions are placed in a full circle with frequencies of for each angle shown via the lengths of each pedal. Colors along the pedal show the conditional histogram of another variable along the pedal.

Value

A plot is produced. If assigned to an object, then a list object is returned with the components:

summary	a list object giving the histogram information for each petal.
number.obs	numeric giving the number of observations.
number.calm	numeric giving the number of zero wspd and missing values.

Author(s)

Matt Pocernich

See Also

bearing

Examples

set.seed( 1001 )
wdir <- runif( 1000, 0, 360 )
set.seed( 2002 )
wspd <- rgamma( 1000, 15 )
CircleHistogram( wspd = wspd, wdir = wdir, leg = TRUE )

---

Cluster Analysis Verification

Description

Perform Cluster Analysis (CA) verifcation per Marzban and Sandgathe (2006).

Usage

clusterer(X, Y = NULL, ...)

## Default S3 method:
clusterer(X, Y = NULL, ..., xloc = NULL, xyp = TRUE, threshold = 1e-08, 
    linkage.method = "complete", stand = TRUE, trans = "identity", 
    a = NULL, verbose = FALSE)

## S3 method for class 'SpatialVx'
clusterer(X, Y = NULL, ..., time.point = 1, obs = 1, model = 1, xyp = TRUE, 
    threshold = 1e-08, linkage.method = "complete", stand = TRUE, 
    trans = "identity", verbose = FALSE)

## S3 method for class 'clusterer'
plot(x, ..., mfrow = c(1, 2), col = c("gray", tim.colors(64)), 
    horizontal = FALSE)

## S3 method for class 'summary.clusterer'
plot(x, ...)

## S3 method for class 'clusterer'
print(x, ...)

## S3 method for class 'clusterer'
summary(object, ...)

Arguments

X, Y	clusterer default method, these are m by n matrices giving the verification and forecast fields, resp. “SpatialVx” method function, X is an object of class “SpatialVx” and Y is not used (a warning is given if it is not missing and not NULL).
object, x	list object of class “clusterer” as returned by clusterer (or summary.clusterer in the case of plot.summary.clusterer).
xloc	(optional) numeric mn by 2 matrix giving the gridpoint locations. If NULL, this will be created using 1:m and 1:n.
xyp	logical, should the cluster analysis be performed on the locations and intensities (TRUE) or only the locations (FALSE)?
threshold	numeric of length one or two giving the threshold to apply to each field (>=). If length is two, the first value corresponds to the threshold for the verification field, and the second to the foreast field.
linkage.method	character naming a valid linkage method accepted by hclust.
stand	logical, should the data matrices consisting of xloc and each field first be standardized before performing cluster analysis?
trans	character naming a function to be applied to the field intensities before performing the CA. Only used if xyp is TRUE. Default applies no transformation.
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.
a	(optional) list giving object attributes associated with a “SpatialVx” class object. The clusterer method for “SpatialVx” objects calls the default method function, and uses this argument to pass the attributes through to the final returned object, as well as to grab location information.
mfrow	mfrow parameter (see help file for par). If NULL, then the parameter is not re-set.
col	color vector for image plots of fields after applying the threshold(s).
horizontal	logical, should the image plot color legend be placed horizontally or vertically? Only for image plot sof the fields.
verbose	logical, should progress information be printed to the screen?
...	optional arguments to the hclust function. In the case of the summary method function, z and/or sigma giving a numeric value used to find the cut-off given by median + z*sigma for detemining matched objects (see Marzban and Sandgathe 2006) where defaults of 1 and the standard deviation of minimum inter-cluster distances are used, and silent (logical should information be printed to the screen (FALSE) or not (TRUE); default is to print to the screen. In the case of the plot method functions, these are optional arguments to the summary method function.

Details

This function performs cluster analysis (CA) on positive values from each of two fields in a verification set using the hclust function from package fastcluster. Inter-cluster distances are computed between each cluster of each field at every level of the CA. The function clusterer performs CA on both fields, and finds the inter-cluster distances across fields for every possible combination of objects at each iteration of each CA. The summary method function finishes the analysis by determining hits, misses and false alarms as well as the numbers of clusters. It also computes CSI for each number of cluster combinations. This is the verification approach described in Marzban and Sandgathe (2006).

The plot method function creates a 4 by 2 panel of plots. The top two plots give image plots of the verification and forecast fields with grid points below the threshold(s) showing zero. The next two plots are dendrograms as performed by the plot method function for hclust (dendrogram) objects. The next row gives a histogram of the minimum inter-cluster distances, then box plots showing the hits, misses and false alarms for every possible combination of levels of each CA. Finally, the bottom two plots show, for each combination of CA level (i.e., numbers of clusters), the CSI and average error (inter-cluster distance) for all matched objects. These last three plots are the ones made by the plot method for values returned from the summary method function.

print is currently not very useful here, but it prevents printing a big mess to the screen.

Value

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

linkage.method	character vector of length one or two giving the linkage method as passed into the function. The length is two only if the McQuitty method is chosen in which case this method is used for the CA, but not for the inter-cluster differencs across fields (average is used for that instead).
trans	character naming the transformation function applied to the intensities.
N	numeric giving the size of the fields.
threshold	numeric of length two giving the threshold applied to each field.
NCo, NCf	numeric vectors giving the number of clusters at each iteration of the CA for the verification and forecast fields, resp.
cluster.identifiers	a list with components X and Y giving lists of lists identifying specific CA components at each level of the CA for both fields.
idX, idY	logical vectors describing which grid points were included in the CA for each field (i.e., which grid points were >= threshold and had non-missing values).
cluster.objects	a list with components X and Y giving the objects returned by hclust for each field.
inter.cluster.dist	a list of list objects with NCf by NCo matrix components giving the inter-cluster distances (between verification and forecast fields) for each iteration of CA for each field.
min.intercluster.dists	numeric vector givng the minimum values inter.cluster.dist at each iteration. Used to determine the cut-off for matched objects.

The summary method function returns a list with the same components as above, but also the components:

cutoff	The cut-off value used for determining matches.
csi, AvgErr	NCo by NCf numeric matrix giving the critical success index (CSI) and average intercluster error (distance) based on matched/un-matched objects.
HMF	NCo by NCf by 3 array giving the hits, misses and false alarms based on matched/un-matched objects.

If the argument a is not NULL, then these are returned as attributes of the returned object. In the case of “SpatialVx” objects, the attributes are preserved.

plot and print methods do not return anything.

Warning

Although some effort has been put into making the functions in this package as computationally efficient as possible, there is a lot of bookeeping involved with this approach, and the current functions are probably not as efficient as they could be. In any case, they will likely be slow for large data sets. The function can work quickly on large fields if an adequately high threshold is used (e.g., if threshold=10 is replaced for 16 in the not run example below, the function is VERY slow). Performing the actual cluster analysis on each field is fast because the hclust function from the fastcluster package is used, which works very well. However, bookeeping after the CA is done employs a lot of loops within loops, which possibly can be made more efficient (and maybe someday will be), but for now...

If it is desired to simply look at the CA for the two fields, the function hclust from fastcluster can be used, which essentially replaces the hclust function from the stats package with a faster version, but otherwise operates the same as far as what is returned, etc., and the same method functions can be employed.

Note

Contact Caren Marzban, marzban “at” u.washington.edu, for questions about the method, and Eric Gilleland, ericg “at” ucar.edu, for problems with the code.

Author(s)

Eric Gilleland

References

Marzban, C. and Sandgathe, S. (2006) Cluster analysis for verification of precipitation fields. Wea. Forecasting, 21, 824–838.

See Also

hclust, hclust, as.dendrogram, cutree, make.SpatialVx, CSIsamples

Examples

data( "UKobs6" )
data( "UKfcst6" )
look <- clusterer(X=UKobs6, Y=UKfcst6, threshold=16, trans="log", verbose=TRUE)
plot( look )

## Not run: 
data( "UKloc" )

# Now, do the same thing, but using a "SpatialVx" object.
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" )

look2 <- clusterer(hold, threshold=16, trans="log", verbose=TRUE)
plot( look2 )
# Note that values differ because now we're using the
# actual locations instead of integer indicators of
# positions.

## End(Not run)

---

Combine Features/Matched Objects

Description

Combine two or more features or matched class objects into one object for aggregation purposes.

Usage

combiner(...)

Arguments

...	Two or more objects of class “features” or “matched” (can also be a list of these objects).

Details

Useful for functions such as compositer and/or (coming soon) aggregating results for feature-based methods.

Value

A list object of class “combined” with the same components as the input arguments, but where some components (namely, X, Xhat, X.labeled, Y.labeled) are now arrays containing these values from each combined object. The lists of lists contained in the X.feats and Y.feats include one long list of lists containing all of the individual features from each object.

Author(s)

Eric Gilleland

See Also

Functions that create objects of class “features” and “matched”: FeatureFinder, centmatch, and deltamm.

Examples

# TO DO

---

Create Composite Features

Description

After identifying features in a verification set, re-grid them so that their centroids are all the same, and the new grid is as small as possible to completely contain all of the features in the verification set.

Usage

compositer(x, level = 0, verbose = FALSE, ...)

## S3 method for class 'features'
compositer(x, level = 0, verbose = FALSE, ...)

## S3 method for class 'matched'
compositer(x, level = 0, verbose = FALSE, ...)

## S3 method for class 'combined'
compositer(x, level = 0, verbose = FALSE, ...) 

## S3 method for class 'composited'
plot(x, ..., type = c("all", "X", "Xhat", "X|Xhat", "Xhat|X"),
    dist.crit = 100, FUN = "mean", col = c("gray", tim.colors(64)))

Arguments

x	compositer: an object of class “features”, “matched”, or “combined”. plot: an object of class “composited”.
type	character, stating which composite features should be plotted. Default makes a two by two panel of plots with all of the choices.
dist.crit	maximum value beyond which any minimum centroid distances are considered too far for features to be “present” in the area.
FUN	name of a function to be applied to the composites in order to give the distributional summary.
col	color pallette to be used.
level	numeric used in shrinking the grid to a smaller size.
verbose	logical, should information be printed to the screen?
...	Not used by compositer. plot: Optional arguments to image.plot (only for the scale legend), such as legend.only, legend.lab, legend.mar, horizontal, etc.

Details

This is functionality for performing an analysis similar to the composite verification method of Nachamkin (2004). See also Nachamkin et al. (2005) and Nachamkin (2009). The main difference is that this function centers all features to the same point, then re-sizes the grid to the smallest possible size to contain all features. The "existence" of a feature at the same time point is determined by the centroid distance (because, here, the compositing is done for a large field rather than a small area), but it does not allow for having half of the feature in the domain in order to be considered.

compositer takes an object of class “features” or “matched” and centers all of the identified features onto the same point so that all features have the same centroid. It also then re-grids the composited features so that they are contained on the smallest possible domain that includes all of the features in the verification set.

Generally, because the composite approach is distributional in nature, it makes sense to look at features across multiple time points. The function combiner allows for combining features from more than one object of class “features” or “matched” in order to subsequently run with compositer.

plot takes the composite features and adds them together creating a density of the composite features, then, Depending on the type argument, the verification (type = “X”), model (type = “Xhat”), verification conditioned on the model (type = “X|Xhat”), or the model conditioned on the verification composite features are plotted. In the case of type = “all”, then a panel of four plots are made with all of these choices. In the case of the conditional plots, the sum of composites for one field are masked out so that only the density of the other field is plotted where composited features from the first field exist.

Value

A list object of class “composited” is returned with all of the same components and attributes as the x argument, but with additional components:

distances	List with components X and Xhat giving the minimum centroid distances from each feature in X (Xhat) to a feature in the other field (used for determining the conditional distributions; i.e., a feature is present if its centroid distance is less than some pre-specified amount).
Xcentered, Ycentered	list of “owin” objects containing each feature similar to X.feats and Y.feats, but centered on the same spot and re-gridded

Note

centroids are rounded to the nearest whole number so that interpolation is not necessary. This may introduce a slight bias in results, but it should not be a major issue.

Author(s)

Eric Gilleland

References

Nachamkin, J. E. (2004) Mesoscale verification using meteorological composites. Mon. Wea. Rev., 132, 941–955.

Nachamkin, J. E. (2009) Application of the Composite Method to the Spatial Forecast Verification Methods Intercomparison Dataset. Wea. Forecasting, 24 (5), 1390–1400, DOI: 10.1175/2009WAF2222225.1.

Nachamkin, J. E., Chen, S. and Schmidt, J. S. (2005) Evaluation of heavy precipitation forecasts using composite-based methods: A distributions-oriented approach. Mon. Wea. Rev., 133, 2163–2177.

See Also

Identifying features: FeatureFinder

Examples

x <- y <- matrix(0, 100, 100)
x[2:3,c(3:6, 8:10)] <- 1
y[c(4:7, 9:10),c(7:9, 11:12)] <- 1

x[30:50,45:65] <- 1
y[c(22:24, 99:100),c(50:52, 99:100)] <- 1

hold <- make.SpatialVx( x, y, field.type = "contrived", units = "none",
    data.name = "Example", obs.name = "x", model.name = "y" )

look <- FeatureFinder(hold, smoothpar=0.5)

look2 <- compositer(look)
plot(look2, horizontal = TRUE)

---

Forecast Verification with Cluster Analysis: The Variation

Description

A variation on cluster analysis for forecast verification as proposed by Marzban and Sandgathe (2008).

Usage

CSIsamples(x, ...)

## Default S3 method:
CSIsamples(x, ..., xhat, nbr.csi.samples = 100, threshold = 20, 
    k = 100, width = 25, stand = TRUE, z.mult = 0, hit.threshold = 0.1, 
    max.csi.clust = 100, diss.metric = "euclidean", linkage.method = "average", 
    verbose = FALSE)

## S3 method for class 'SpatialVx'
CSIsamples(x, ..., time.point = 1, obs = 1, model = 1, nbr.csi.samples = 100, 
    threshold = 20, k = 100, width = 25, stand = TRUE, z.mult = 0, 
    hit.threshold = 0.1, max.csi.clust = 100, diss.metric = "euclidean", 
    linkage.method = "average", verbose = FALSE)

## S3 method for class 'CSIsamples'
summary(object, ...)

## S3 method for class 'CSIsamples'
plot(x, ...)

## S3 method for class 'summary.CSIsamples'
plot(x, ...)

## S3 method for class 'CSIsamples'
print(x, ...)

Arguments

x, xhat	default method: matrices giving the verification and forecast fields, resp. “SpatialVx” method: x is an object of class “SpatialVx”. plot, print methods: list object of class “CSIsamples” or “summary.CSIsamples” (in the case of plot).
object	list object of class “CSIsamples”.
nbr.csi.samples	integer giving the number of samples to take at each level of the CA.
threshold	numeric giving a value over which is to be considered an event.
k	numeric giving the value for centers in the call to kmeans.
width	numeric giving the size of the samples for each cluster sample.
stand	logical, should the data first be standardized before applying CA?
z.mult	numeric giving a value by which to multiply the z- component. If zero, then the CA is performed on locations only. Can be used to give more or less weight to the actual values at these locations.
hit.threshold	numeric between zero and one giving the threshold for the proportion of a cluster that is from the verification field vs the forecast field used for determining whether the cluster consitutes a hit (vs false alarm or miss depending).
max.csi.clust	integer giving the maximum number of clusters allowed.
diss.metric	character giving which method to use in the call to dist (which dissimilarity metric should be used?).
linkage.method	character giving the name of a linkage method acceptable to the method argument from the hclust function of package fastcluster.
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.
verbose	logical, should progress information be printed to the screen?
...	Not used by CSIsamples method functions. summary method function: the argument silent may be specified, which is a logical stating whether to print the information to the screen (FALSE) or not (TRUE). If not given, the summary information will be printed to the screen. Not used by the plot method function.

Details

This function carries out the procedure described in Marzban and Sandgathe (2008) for verifying forecasts. Effectively, it combines the verification and forecast fields (keeping track of which values belong to which field) and applies CA to the combined field. Clusters identified with a proportion of values belonging to the verification field within a certain range (defined by the hit.threshold argument) are determined to be hits, misses or false alarms. From this information, the CSI (at each number of clusters; scale) is calculated. A sampling scheme is used to speed up the process.

The plot and summary functions all give the same information, but in different formats: i.e., CSI by number of clusters (scale).

Value

A list is returned by CSIsamples with components:

data.name	character vector giving the names of the verification and forecast fields analyzed, resp.
call	an object of class “call” giving the function call.
results	max.csi.clust by nbr.csi.samples matrix giving the caluclated CSI for each sample and iteration of CA.

The summary method function invisibly returns the same list, but with the additional component:

csi	vector of length max.csi.clust giving the sample average CSI for each iteration of CA.

The plot method functions do not return anything. Plots are created.

Note

Special thanks to Caren Marzban, marzban “at” u.washington.edu, for making the CSIsamples (originally called csi.samples) function available for use with this package.

Author(s)

Hillary Lyons, h.lyons “at” comcast.net, and modified by Eric Gilleland

References

Marzban, C., Sandgathe, S. (2008) Cluster Analysis for Object-Oriented Verification of Fields: A Variation. Mon. Wea. Rev., 136, (3), 1013–1025.

See Also

hclust, hclust, kmeans, clusterer

Examples

## Not run: 
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)
look2 <- sim.rf( obj)

res <- CSIsamples(x=look, xhat=look2, 10, threshold=0, k=100,
                  width=2, z.mult=0, hit.threshold=0.25, max.csi.clust=75)
plot(res)
y <- summary(res)
plot(y)

## End(Not run)
## Not run: 
data( "UKfcst6" )
data( "UKobs6" )
data( "UKloc" )

hold <- make.SpatialVx(UKobs6, UKfcst6, thresholds=0,
    loc=UKloc, map=TRUE, field.type="Rainfall", units="mm/h",
    data.name = "Nimrod", obs.name = "obs 6", model.name = "fcst 6" )

res <- CSIsamples( hold, threshold = 0, k = 200, z.mult = 0.3, hit.threshold = 0.2,
                  max.csi.clust = 150, verbose = TRUE)
plot( res )

summary( res )

y <- summary( res )

plot( y )

## End(Not run)

---

Merge and/or Match Identified Features Within Two Fields

Description

Merge and/or match identified features within two fields using the delta metric method described in Gilleland et al. (2008), or the matching only method of Davis et al. (2006a).

Usage

deltamm(x, p = 2, max.delta = Inf, const = Inf, type = c( "sqcen", "original" ),
    N = NULL, verbose = FALSE, ...)

centmatch(x, criteria = 1, const = 14, distfun = "rdist", areafac = 1,
    verbose = FALSE, ...)

## S3 method for class 'matched'
plot(x, mfrow = c(1, 2), ...)

## S3 method for class 'matched'
print(x, ...)

## S3 method for class 'matched'
summary(object, ...)

Arguments

x	For deltamm and centmatch, an object of class “features” with components X.labeled, Y.labeled (matrices with numbered features), as well as, X.feats and Y.feats, each of which are list objects containing numbered components within which are objects of class “owin” containing logical matrices that define features for the forecast (Y) and verification (X) fields, resp. For example, as returned by the FeatureFinder function. For plot an object of class “matched”. Argument y is used if it is not NULL, otherwise argument x is used (but only one of x or y is used). If x and y are missing, but not object, then object will be used, in which case it must be of class “features”.
object	list object of class “matched”.
p	Baddeley delta metric parameter. A value of 1 gives arithmetic averages, Inf gives the Hausdorff metric and -Inf gives a minimum. The default of 2 is most common.
max.delta	single numeric giving a cut-off value for delta that disallows two features to be merged or matched if the delta between them is larger than this value.
const	deltamm: a constant value over which the shortest distances in a distance map are set. See Gilleland (2011) and Schwedler and Baldwin (2011) for more information about this parameter.
type	character specifying whether Baddeley's delta metric should be calculated after centering object pairs on a new square grid (default) or performed in their original positions on the original grid.
N	If type is “sqcen”, then N is the argument to censqdelta that specifies the common grid size. If not specified the maximum side of the original domain is used (possibly adding one first to make it an odd number). It is possible that values could be pushed off the new grid, and making N larger might alleviate the issue.
centmatch	numeric giving the number of grid squares whereby if the centroid distance (D) is less than this value, a match is declared (only used if criteria is 3.
criteria	1, 2 or 3 telling which criteria for determining a match based on centroid distance, D, to use. The first (1) is a match if D is less than the sum of the sizes of the two features in question (size is the square root of the area of the feature). The second is a match if D is less than the average size of the two features in question. The third is a match if D is less than a constant given by the argument const.
distfun	character string naming a distance function. Default uses rdist from the fields package.
areafac	single numeric used to multiply by grid-space based area in order to at least approximate the correct distance (e.g., using the ICP test cases, 4 would make the areas approximately square km instead of grid points). This should not be used unless distfun is “rdist.earth” in which case it will use the locations given in the call to make.SpatialVx, which are assumed to be lat/lon coordinates.
mfrow	mfrow parameter (see help file for par). If NULL, then the parameter is not re-set.
verbose	logical, should progress information be printed to the screen?
...	For deltamm: additional optional arguments to the distmap function from package spatstat. For centmatch: optional arguments to distfun. For plot.matched, additional arguments to image.plot concerning the color legend bar only.

Details

deltamm:

Gilleland et al. (2008) describe a method for automatically merging, and simultaneously, matching identified features within two fields (a verification set). The method was proposed with the general method for spatial forecast verification introduced by Davis et al. (2006 a,b) in mind. It relies heavily on use of a binary image metric introduced by Baddeley (1992a,b) for comparing binary images; henceforth referred to as the delta metric, or just delta.

The procedure is as follows. Suppose there are m identified forecast features and n identified verification features.

1. Compute delta for each feature identified in the forecast field against each feature identified in the verification field. Store these values in an m by n matrix, Upsilon.

2. For each of the m rows of Upsilon, rank the values of delta to identify the features, j_1, ..., j_n that provide the lowest (best) to highest (worst) value, and do the same for each of the n columns to find the forecast features i1, ...,i_m that yield the lowest to highest values for each verification feature.

3. Create a new m by n matrix, Psi, whose columns contain delta computed between each of the individual features in the forecast and (first column) the corresponding j_1 feature from the verification field, and each successive column, k, has delta between the i-th forecast feature and the union of j_1, j_2, ..., j_k.

4. Create a similar m by n matrix, Ksi, that has delta computed between each individual feature in the verification field and the successively bigger unions i_1, ..., i_l for the l-th column.

5. Let Q=[Upsilon, Psi, Ksi], and merge and match features based on the rankings of delta in Q. That is, find the smallest delta in Q, and determine which mergings (if any) and matchings correspond to this value. Remove the appropriate row(s) and column(s) of Q corresponding to the already determined matchings and/or mergings. Repeat this until all features in at least one field have been exhausted.

The above algorithm suffers from two deficiencies. First, features that are merged in one field cannot be matched to merged features in another field. One possible remedy for this is to run this algorithm twice, though this is not a universally good solution. Second, features can be merged and/or matched to features that are very different from each other. A possible remedy for this is to use the cut-off argument, max.delta, to disallow mergings or matchings between features whose delta value is not <= this cut-off. In practice, these two deficiencies are not likely very problematic.

centmatch:

This function works similarly as deltamm, though it does not merge features. It is based on the method proposed by Davis et al. (2006a). It is possible for more than one object to be matched to the same object in another field. As a result, when plotting, it might appear that features have been merged, but they have not been. For informational purposes, the criteria, appelled criteria.values (as determined by the criteria argument), along with the centroid distance matrix, appelled centroid.distances, are returned.

plot: The plot method function for matched features plots matched features across fields in the same color using rainbow. Unmatched features in either field are all colored gray. Zero values are colored white. The function MergeForce must first be called, however, in order to organize the object into a format that allows the plot method function to determine the correct color coding.

The print method function will tell you which features matched between fields, so one can plot the originally derived features (e.g., from FeatureFinder) to identify matched features.

summary:

The summary method function so far simply reverts the class back to “features” and calls that summary function.

Value

A list object of class “matched” is returned by both centmatch and deltamm containing several components added to the value of x or y passed in, and possibly with attributes inhereted from object.

match.message	A character string stating how features were matched.
match.type	character string naming the matching function used.
matches	two-column matrix with forecast object numbers in the first column and corresponding matched observed features in the second column. If no matches, this will have value integer(0) for each column giving a matrix with dimension 0 by 2.
unmatched	list with components X and Xhat giving the unmatched object numbers, if any, from the observed and forecast fields, resp. If none, the value will be integer(0).
Q	(deltamm only) an array of dimension n by m by 3 giving all of the delta values that were computed in determining the mergings and matchings.
criteria	(centmatch only) 1, 2, or 3 as given by the criteria argument.
criteria.values, centroid.distances	(centmatch only) matrices giving the forecast by observed object criteria and centroid distances.
implicit.merges	(centmatch only) list displaying multiple matches for each field (this could define potential merges). Each component of the list is a unified set of matched features in the form of two-column matrices analogous to the matches component. If there are no implicit mergings or no matched features, this component will be named, but also NULL. Note: such implicit mergings may or may not make physical sense, and are not considered to be merged generally, but will show up as having been merged/clustered when plotted.

If the argument ‘object’ is passed in, then the list object will also contain nearly the same attributes, with the data.name attribute possibly changed to reflect the specific model used. It will also contain a time.point and model attribute.

Author(s)

Eric Gilleland

References

Baddeley, A. (1992a) An error metric for binary images. In Robust Computer Vision Algorithms, W. Forstner and S. Ruwiedel, Eds., Wichmann, 59–78.

Baddeley, A. (1992b) Errors in binary images and an Lp version of the Hausdorff metric. Nieuw Arch. Wiskunde, 10, 157–183.

Davis, C. A., Brown, B. G. and Bullock, R. G. (2006a) Object-based verification of precipitation forecasts, Part I: Methodology and application to mesoscale rain areas. Mon. Wea. Rev., 134, 1772–1784.

Davis, C. A., Brown, B. G. and Bullock, R. G. (2006b) Object-based verification of precipitation forecasts, Part II: Application to convective rain systems. Mon. Wea. Rev., 134, 1785–1795.

Gilleland, E. (2011) Spatial Forecast Verification: Baddeley's Delta Metric Applied to the ICP Test Cases. Wea. Forecasting, 26 (3), 409–415.

Gilleland, E., Lee, T. C. M., Halley Gotway, J., Bullock, R. G. and Brown, B. G. (2008) Computationally efficient spatial forecast verification using Baddeley's delta image metric. Mon. Wea. Rev., 136, 1747–1757.

Schwedler, B. R. J. and Baldwin, M. E. (2011) Diagnosing the sensitivity of binary image measures to bias, location, and event frequency within a forecast verification framework. Wea. Forecasting, 26, 1032–1044.

Examples

## Not run: 
x <- y <- matrix(0, 100, 100)
x[2:3,c(3:6, 8:10)] <- 1
y[c(4:7, 9:10),c(7:9, 11:12)] <- 1

x[30:50,45:65] <- 1
y[c(22:24, 99:100),c(50:52, 99:100)] <- 1

hold <- make.SpatialVx( x, y, field.type = "contrived", units = "none",
    data.name = "Example", obs.name = "x", model.name = "y" )

look <- FeatureFinder( hold, smoothpar = 0.5 ) 

# The next line fails because the centering pushes one object out of the new domain.
# look2 <- deltamm( look )
# Setting N larger fixes the problem.
look2 <- deltamm( look, N = 300 )
look2 <- MergeForce( look2 )

look2

plot( look2 )

FeatureTable(look2)

look3 <- centmatch(look)

FeatureTable(look3)

look3 <- MergeForce( look3 )

plot( look3 )

## End(Not run)

---

Identify Disjoint Sets of Connected Components

Description

Identify disjoint sets of contiguous events in a binary field. In many areas of research, this function finds connected components.

Usage

disjointer(x, method = "C")

Arguments

x	A numeric matrix or other object that as.im from package spatstat works on.
method	Same argument as that in connected from package spatstat.

Details

disjointer essentially follows the help file for connected to produce a list object where each component is an image describing one set of connected components (or blobs). It is essentially a wrapper function to connected. This function is mainly used internally by FeatureFinder and similar, but could be of use outside such functions.

Value

An unnamed list object where each component is an image describing one set of connected components (or blobs).

Author(s)

Eric Gilleland

References

Park, J.-M., Looney, C.G. and Chen, H.-C. (2000) Fast connected component labeling algorithm using a divide and conquer technique. Pages 373-376 in S.Y. Shin (ed) Computers and Their Applications: Proceedings of the ISCA 15th International Conference on Computers and Their Applications, March 29–31, 2000, New Orleans, Louisiana USA. ISCA 2000, ISBN 1-880843-32-3.

Rosenfeld, A. and Pfalz, J.L. (1966) Sequential operations in digital processing. Journal of the Association for Computing Machinery 13 471–494.

Examples

##
## For examples, see FeatureFinder
##

---

Elmore, Baldwin and Schultz Method for Field Significance for Spatial Bias Errors

Description

Apply the method of Elmore, Baldwin and Schultz (2006) for calculating field significance of spatial bias errors.

Usage

EBS(object, model = 1, block.length = NULL, alpha.boot = 0.05,
    field.sig = 0.05, bootR = 1000, ntrials = 1000,
    verbose = FALSE)

## S3 method for class 'EBS'
plot(x, ..., mfrow = c(1, 2), col, horizontal)

Arguments

object	list object of class “SpatialVx”.
x	object of class “EBS” as returned by EBS.
model	number or character describing which model (if more than one in the “SpatialVx” object) to compare.
block.length	numeric giving the block length to be used n the block bootstrap algorithm. If NULL, floor(sqrt(n)) is used.
alpha.boot	numeric between 0 and 1 giving the confidence level desired for the bootstrap algorithm.
field.sig	numeric between 0 and 1 giving the desired field significance level.
bootR	numeric integer giving the number of bootstrap replications to use.
ntrials	numeric integer giving the number of Monte Carol iterations to use.
mfrow	mfrow parameter (see help file for par). If NULL, then the parameter is not re-set.
col, horizontal	optional arguments to image.plot from fields.
verbose	logical, should progress information be printed to the screen?
...	optional arguments to image.plot from fields.

Details

this is a wrapper function for the spatbiasFS function utilizing the “SpatialVx” object class to simplify the arguments.

Value

A list object of class “EBS” with the same attributes as the input object and additional attribute (called “arguments”)that is a named vector giving information provided by the user. Components of the list include:

block.boot.results	object of class “LocSig”.
sig.results	list object containing information about the significance of the results.

Author(s)

Eric Gilleland

References

Elmore, K. L., Baldwin, M. E. and Schultz, D. M. (2006) Field significance revisited: Spatial bias errors in forecasts as applied to the Eta model. Mon. Wea. Rev., 134, 519–531.

See Also

boot::boot, boot:tsboot, spatbiasFS, LocSig, make.SpatialVx

Examples

data( "GFSNAMfcstEx" )
data( "GFSNAMobsEx" )
data( "GFSNAMlocEx" )

id <- GFSNAMlocEx[,"Lon"] >=-95
id <- id & GFSNAMlocEx[,"Lon"] <= -75
id <- id & GFSNAMlocEx[,"Lat"] <= 32

##
## This next step is a bit awkward, but these data
## are not in the format of the SpatialVx class.
## These are being set up with arbitrarily chosen
## dimensions (49 X 48) for the spatial part.  It
## won't matter to the analyses or plots.
##
Vx <- GFSNAMobsEx
Fcst <- GFSNAMfcstEx
Ref <- array(t(Vx), dim=c(49, 48, 361))
Mod <- array(t(Fcst), dim=c(49, 48, 361)) 

hold <- make.SpatialVx(Ref, Mod, loc=GFSNAMlocEx,
    projection=TRUE, map=TRUE, loc.byrow = TRUE, subset=id,
    field.type="Precipitation", units="mm",
    data.name = "GFS/NAM", obs.name = "Reference", model.name = "Model" )

look <- EBS(hold, bootR=500, ntrials=500, verbose=TRUE)
plot( look )

## Not run: 
# Same as above, but now we'll do it for all points.
# A little slower, but not terribly bad.

hold <- make.SpatialVx(Ref, Mod, loc = GFSNAMlocEx,
    projection = TRUE, map = TRUE, loc.byrow = TRUE,
    field.type = "Precipitation", reg.grid = FALSE, units = "mm",
    data.name = "GFS/NAM", obs.name = "Reference", model.name = "Model" )

look <- EBS(hold, bootR=500, ntrials=500, verbose=TRUE)
plot( look )

## End(Not run)

---

Simulated Spatial Verification Set

Description

A simulated spatial verification set for use by various examples for this package.

Usage

data(ExampleSpatialVxSet)

Format

The format is: List of 2 $ vx : num [1:50, 1:50] 0 0 0 0 0 0 0 0 0 0 ... $ fcst: num [1:50, 1:50] 0.0141 0 0 0 0 ...

Details

The data here were generated using the sim.rf function from fields (Furrer et al., 2012):

x <- y <- matrix(0, 10, 12) x[2:3,c(3:6, 8:10)] <- 1 y[c(1:2, 9:10),c(3:6)] <- 1

grid <- list(x=seq(0,5,,50), y=seq(0,5,,50)) obj <- Exp.image.cov(grid=grid, theta=0.5, setup=TRUE) x <- sim.rf(obj) x[x < 0] <- 0 x <- zapsmall(x)

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

References

Reinhard Furrer, Douglas Nychka and Stephen Sain (2012). fields: Tools for spatial data. R package version 6.6.3. http://CRAN.R-project.org/package=fields

Examples

data( "ExampleSpatialVxSet" )
x <- ExampleSpatialVxSet$vx
xhat <- ExampleSpatialVxSet$fcst
par(mfrow=c(1,2))
image.plot(x, col=c("gray",tim.colors(64)))
image.plot(xhat, col=c("gray",tim.colors(64)))

---

Exponential Variogram

Description

Compute the exponential variogram.

Usage

expvg(p, vg, ...)

## S3 method for class 'flossdiff.expvg'
predict(object, newdata, ...)

## S3 method for class 'flossdiff.expvg'
print(x, ...)

Arguments

p	numeric vector of length two. Each component should be positively valued. The first component is the nugget and the second is the range parameter.
vg	A list object with component d giving a numeric vector of distances over which the variogram is to be calculated.
object, x	A list object returned by flossdiff using expvg as the variogram model.
newdata	Numeric giving the distances over which to use the fitted exponential variogram model to make predictions. The default is to go from zero to the maximum lag distance for a given data set, which is not the usual convention for the generic predict, which usually defaults to operate on the lags used in performing the fit.
...	Not used.

Details

A very simple function used mainly internally by flossdiff when fitting the exponential variogram to the empirical one, and by the predict, print and summary method functions for lossdiff objects. For those wishing to use a different variogram model than the exponential, use this function and its method functions as a template. Be sure to create predict and print method functions to operate on objects of class “flossdiff.XXX” where “XXX” is the name of the variogram function you write (so, “expvg” in the current example).

Value

Numeric vector of length equal to that of the d component of vg giving the corresponding exponential variogram values with nugget and range defined by p.

Author(s)

Eric Gilleland

References

Cressie, N. A. (2015) Statistics for Spatial Data. Wiley-Interscience; Revised Edition edition (July 27, 2015), ISBN-10: 1119114616, ISBN-13: 978-1119114611, 928 pp.

See Also

lossdiff, flossdiff

Examples

##
## For examples, see lossdiff and flossdiff
##

---

Exponential Variogram

Description

Calculates the empirical variogram for use with function spct.

Usage

expvgram(p, h, ...)

Arguments

p	numeric vector of length two giving the nugget and range parameter values, resp.
h	numeric vector of separation distances.
...	Not used.

Details

Simple function to work with spct to calculate the exponential variogram for given parameters and separation distances. The exponential variogram employed here is parameterized by

gamma(h) = sigma * ( 1 - exp( - h * theta ) )

where p is c( sigma, theta ).

Value

A numeric vector of variogram values for each separation distance in h.

Author(s)

Eric Gilleland

See Also

spct

Examples

# See help file for spct for examples.

---

Major and Minor Axes of a Feature

Description

Calculate the major and minor axes of a feature and various other properties such as the aspect ratio.

Usage

FeatureAxis(x, fac = 1, flipit = FALSE, twixt = FALSE)

## S3 method for class 'FeatureAxis'
plot(x, ..., zoom = FALSE)

## S3 method for class 'FeatureAxis'
summary(object, ...)

Arguments

x	For FeatureAxis this is an object of class “owin” containing a binary image matrix defining the feature. In the case of plot.FeatureAxis, this is the value returned from FeatureAxis.
object	list object of class “FeatureAxis” as returned by FeatureAxis.
fac	numeric, in determining the lengths of the axes, they are multiplied by a factor of fac (e.g., if the grid points are k by k km each, then one could set this to k so that the resulting lengths are in terms of km rather than grid points.
flipit	logical, should the objects be flipped over x and y? The disjointer function results in images that are flipped, this would flip them back.
twixt	logical, should the major axis angle be forced to be between +/- 90 degrees?
zoom	logical, should the object be plotted on its bounding box (TRUE) or on the original grid (FALSE, default)? Useful if the feature is too small to be seen well on the original gid.
...	For plot.FeatureAxis these are additional arguments to the plot function. Not used by summary.FeatureAxis.

Details

This function attempts to identify the major and minor axes for a pre-defined feature (sometimes referred to as an object). This function relies heavily on the spatstat and smatr packages. First, the convex hull of the feature is determined using the convexhull function from the spatstat package. The major axis is then found using the sma function from package smatr, which is then converted into a psp object (see as.psp from spatstat) from which the axis angle and length are found (using angles.psp and lengths_psp, resp., from spatstat).

The minor axis anlge is easily found after rotating the major axis 90 degrees using rotate.psp from spatstat. The length of the minor axis is more difficult. Here, it is found by rotating the convex hull of the feature by the major axis angle (so that it is upright) using rotate.owin from spatstat, and then computing the bounding box (using boundingbox from spatstat). The differnce is then taken between the range of x- coordinates of the bounding box. This seems to give a reasonable value for the length of the minor axis. A psp object is then created using the mid point of the major axis (which should be close to the centroid of the feature) using as.psp and midpoints.psp from spatstat along with the length and angle already found for the minor axis.

See the help files for the above mentioned functions for references, etc.

Value

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

z	same as the argument x passed in.
MajorAxis, MinorAxis	a psp object with one segment that is the major (minor) axis.
OrientationAngle	list with two components: MajorAxis (the angle in degrees of the major axis wrt the abscissa), MinorAxis (the angle in degrees wrt the abscissa).
aspect.ratio	numeric giving the ratio of the length of the minor axis to that of the major axis (always between 0 and 1).
MidPoint	an object of class “ppp” giving the mid point of the major (minor) axis.
lengths	list object with components: MajorAxis giving the length (possibly multiplied by a factor) of the major axis, and MinorAxis same as MajorAxis but for the minor axis.
sma.fit	The fitted object returned by the sma function. This is useful, e.g., if confidence intervals for the axis are desired. See the sma help file for more details.

No value is returned from the plot or summary method functions.

Author(s)

Eric Gilleland

Examples

data( "ExampleSpatialVxSet" )

x <- ExampleSpatialVxSet$vx

look <- disk2dsmooth(x,5)
u <- quantile(look,0.99)
sIx <- matrix(0, 100, 100)
sIx[ look > u] <- 1
look2 <- disjointer(sIx)[[1]]
look2 <- flipxy(look2)
tmp <- FeatureAxis(look2)
plot(tmp)
summary(tmp)

## Not run: 
data( "pert000" )
data( "pert004" )
data( "ICPg240Locs" )

hold <- make.SpatialVx( pert000, pert004,
    loc = ICPg240Locs, projection = TRUE, map = TRUE,
    loc.byrow = TRUE,
    field.type = "Precipitation", units = "mm/h",
    data.name = "Perturbed ICP Cases", obs.name = "pert000",
    model.name = "pert004" )

look <- FeatureFinder(hold, smoothpar=10.5)
par(mfrow=c(1,2))
plot(look)

par(mfrow=c(2,2))
image.plot(look$X.labeled)
image.plot(look$Y.labeled)

# The next line will likely be very slow.
look2 <- deltamm(x=look, verbose=TRUE)
image.plot(look2$X.labeled)
image.plot(look2$Y.labeled)

look2$mm.new.labels # the first seven features are matched.

ang1 <- FeatureAxis(look2$X.feats[[1]])
ang2 <- FeatureAxis(look2$Y.feats[[1]])
plot(ang1)
plot(ang2)
summary(ang1)
summary(ang2)

ang3 <- FeatureAxis(look2$X.feats[[4]])
ang4 <- FeatureAxis(look2$Y.feats[[4]])
plot(ang3)
plot(ang4)
summary(ang3)
summary(ang4)
   
## End(Not run)

---

Threshold-based Feature Finder

Description

Identify spatial features within a verification set using a threshold-based method.

Usage

FeatureFinder(object, smoothfun = "disk2dsmooth", do.smooth = TRUE,
    smoothpar = 1, smoothfunargs = NULL, thresh = 1e-08, idfun = "disjointer",
    min.size = 1, max.size = Inf, fac = 1, zero.down = FALSE, time.point = 1,
    obs = 1, model = 1, ...)

## S3 method for class 'features'
plot(x, ..., type = c("both", "obs", "model"))

## S3 method for class 'features'
print(x, ...)

## S3 method for class 'features'
summary(object, ...)

## S3 method for class 'summary.features'
plot(x, ...)

Arguments

object	An object of class “SpatialVx”.
x	list object of class “features” as returned by FeatureFinder.
smoothfun	character naming a 2-d smoothing function from package smoothie. Not used if do.smooth is FALSE.
do.smooth	logical, should the field first be smoothed before trying to identify features (resulting field will not be smoothed, this is just for identifying features). Default is to do convolution smoothing using a disc kernel as is recommended by Davis et al (2006a).
smoothpar	numeric of length one or two giving the smoothing parameter for smoothfun. If length is two, the first value is applied to the forecast field and the second to the verification field. The default smooth function (smoothfun argument) is the disk kernel smoother, so this argument gives the radius of the disk. See the help file for hoods2dsmooth from package smoothie for other smoother choices; this argument corresponds to the lambda argument of these functions.
smoothfunargs	list object with named additional arguments to smoothfun.
thresh	numeric vector of length one or two giving the threshold over which (inclusive) features should be identified. If different thresholds are used for the forecast and verification fields, then the first element is the threshold for the forecast, and the second for the verification field.
idfun	character naming the function used to identify (and label) individual features in the thresholded, and possibly smoothed, fields. Must take an argument 'x', the thresholded, and possibly smoothed, field.
min.size	numeric of length one or two giving the minimum number of contiguous grid points exceeding the threshold in order to be included as a feature (can be used to exclude any small features). Default does not exclude any features. If length is two, first value applies to the forecast and second to the verification field.
max.size	numeric of length one or two giving the maximum number of contiguous grid points exceeding the threshold in order to be included as a feature (can be used to exclude large features, if the need be). Default does not exclude any features. If length is two, then the first value applies ot the forecast field, and the second to the verification.
fac	numeric of length one or two giving a factor by which to multiply the R quantile in determining the threshold from the fields. For example, ~ 1/15 is suggested in Wernli et al (2008, 2009). If length is two, then the first value applies to the threshold of the forecast and the second to that of the verification field.
zero.down	logical, should negative values and relatively very small values be set to zero after smoothing the fields? For thresholds larger than such values, this argument is moot. 'zapsmall' is used to set the very small positive values to zero.
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.
type	character string stating which features to plot (observed, forecast or both). If both, a panel of two plots will be made side-by-side.
...	FeatureFinder: additional arguments to idfun. plot: optional arguments to image.plot for the color legend bar only. Not used by the print or summary functions. The 'summary' method function can take the argument: 'silent'-logical, should information be printed to the screen (FALSE) or not (TRUE).

Details

FeatureFinder applies for finding features based on three proposed methods from different papers; and also allows for combinations of the methods. The methods include: the convolution-threshold approach of Davis et al. (2006a,b), which uses a disc kernel convolution smoother to first smooth the fields, then applies a threshold to remove low-intensity areas. Feautres are identified by groups of contiguous “events” (or connected components in the computer vision/image analysis literature) using idfun. Nachamkin (2009) and Lack et al (2010) further require that features have at least min.size connected components in order to be considered a feature (in order to remove very small areas of threshold excesses). Wernli et al. (2009) modify the threshold by a factor (see the fac argument).

In addition to the above options, it is also possible to remove features that are too large, as for some purposes, it is the small-scale features that are of interest, and sometimes the larger features can cause problems when merging and matching features across fields.

Value

FeatureFinder returns a list object of class “features” with comopnents:

data.name	character vector naming the verification and forecast (R object) fields, resp.
X.feats, Y.feats	The identified features for the verification and forecast fields as returned by the idfun function.
X.labeled, Y.labeled	matrices of same dimension as the forecast and verification fields giving the images of the convolved and thresholded verification and forecast fields, but with each individually identified object labeled 1 to the number of objects in each field.
identifier.function, identifier.label	character strings naming the function and giving the long name (for use with plot method function).

An additional attribute, named “call”, is given. This attribute shows the original function call, and is used mainly by the print function..

The plot method functions do not return anything.

The summary method function for objects of class “features” returns a list with components:

X, Y

matrices whose rows are objects and columns are properties: centroidX and centroidY (the x- and y- coordinates for the feature centroids), area (the area of each feature in squared grid points), the orientation angle for the fitted major axis, the aspect ratio, Intensity0.25 and Intensity0.9 (the lower quartile and 0.9 quantile of intensity values for each feature).

Note

This function replaces the now deprecated functions: convthresh, threshsizer and threshfac.

Author(s)

Eric Gilleland

References

Davis, C. A., Brown, B. G. and Bullock, R. G. (2006a) Object-based verification of precipitation forecasts, Part I: Methodology and application to mesoscale rain areas. _Mon. Wea. Rev._, *134*, 1772-1784.

Davis, C. A., Brown, B. G. and Bullock, R. G. (2006b) Object-based verification of precipitation forecasts, Part II: Application to convective rain systems. _Mon. Wea. Rev._, *134*, 1785-1795.

Lack, S. A., Limpert, G. L. and Fox, N. I. (2010) An object-oriented multiscale verification scheme. _Wea. Forecasting_, *25*, 79-92, doi:10.1175/2009WAF2222245.1.

Nachamkin, J. E. (2009) Application of the composite method to the spatial forecast verification methods intercomparison dataset. _Wea. Forecasting_, *24*, 1390-1400, doi:10.1175/2009WAF2222225.1.

Wernli, H., Paulat, M. Hagen, M. and Frei, C. (2008) SAL-A novel quality measure for the verification of quantitative precipitation forecasts. _Mon. Wea. Rev._, *136*, 4470-4487.

Wernli, H., Hofmann, C. and Zimmer, M. (2009) Spatial forecast verification methods intercomparison project: Application of the SAL technique. _Wea. Forecasting_, *24*, 1472-1484, doi:10.1175/2009WAF2222271.1.

Examples

data( "ExampleSpatialVxSet" )

x <- ExampleSpatialVxSet$vx
xhat <- ExampleSpatialVxSet$fcst

hold <- make.SpatialVx( x, xhat, field.type = "simulated",
		       units = "none", data.name = "Example",
		       obs.name = "x", model.name = "xhat" )

look <- FeatureFinder( hold, smoothpar = 0.5, thresh = 1 )

par( mfrow=c(1,2))
image.plot(look$X.labeled)
image.plot(look$Y.labeled)

## Not run: 
x <- y <- matrix(0, 100, 100)
x[2:3,c(3:6, 8:10)] <- 1
y[c(4:7, 9:10),c(7:9, 11:12)] <- 1
     
x[30:50,45:65] <- 1
y[c(22:24, 99:100),c(50:52, 99:100)] <- 1
     
hold <- make.SpatialVx( x, y, field.type = "contrived", units = "none",
         data.name = "Example", obs.name = "x", model.name = "y" )
     
look <- FeatureFinder(hold, smoothpar=0.5) 

par( mfrow=c(1,2))
image.plot(look$X.labeled)
image.plot(look$Y.labeled)
     
look2 <- centmatch(look)
     
FeatureTable(look2)
     
look3 <- deltamm( look, N = 201, verbose = TRUE ) 
FeatureTable( look3 )


# data( "pert000" )
# data( "pert004" )
# data( "ICPg240Locs" )
     
# hold <- make.SpatialVx( pert000, pert004,
#    loc = ICPg240Locs, projection = TRUE, map = TRUE, loc.byrow = TRUE,
#    field.type = "Precipitation", units = "mm/h",
#    data.name = "ICP Perturbed Cases", obs.name = "pert000",
#    model.name = "pert004" )
     
# look <- FeatureFinder(hold, smoothpar=10.5, thresh = 5)
# plot(look)

# look2 <- deltamm( look, N = 701, verbose = TRUE )

# look2 <- MergeForce( look2 )

# plot(look2)

# summary( look2 )

# Now remove smallest features ( those with fewer than 700 grid squares).

# look <- FeatureFinder( hold, smoothpar = 10.5, thresh = 5, min.size = 700 )

# look # Now only two features.

# plot( look )

# Now remove the largest features (those with more than 1000 grid squares). 

# look <- FeatureFinder( hold, smoothpar = 10.5, thresh = 5, max.size = 1000 )

# look

# plot( look )

# Remove any features smaller than 700 and larger than 2000 grid squares).

# look <- FeatureFinder( hold, smoothpar = 10.5, thresh = 5,
 #    min.size = 700, max.size = 2000 )

# look

# plot( look )

# Find features according to Wernli et al. (2008).
# look <- FeatureFinder( hold, thresh = 5, do.smooth = FALSE, fac = 1 / 15 )

# look

# plot( look )

# Now do a mix of the two types of methods.
# look <- FeatureFinder( hold, smoothpar = 10.5, thresh = 5, fac = 1 / 15 )

# look

# plot( look )


## End(Not run)

---

Analyze Features of a Verification Set

Description

Analyze matched features of a verification set.

Usage

FeatureMatchAnalyzer(x, which.comps=c("cent.dist", "angle.diff", "area.ratio", "int.area",
                    "bdelta", "haus", "ph", "med", "msd", "fom", "minsep",
		    "bearing"), sizefac=1, alpha=0.1, k=4, p=2, c=Inf,
		    distfun="distmapfun", ...)

## S3 method for class 'matched.centmatch'
FeatureMatchAnalyzer(x, which.comps=c("cent.dist", "angle.diff",
		    "area.ratio", "int.area", "bdelta", "haus", "ph", "med",
		    "msd", "fom", "minsep", "bearing"), sizefac=1, alpha=0.1, k=4, p=2,
		    c=Inf, distfun="distmapfun", ...)

## S3 method for class 'matched.deltamm'
FeatureMatchAnalyzer(x, which.comps = c("cent.dist", "angle.diff",
		    "area.ratio", "int.area", "bdelta", "haus", "ph", "med", "msd",
		    "fom", "minsep", "bearing"), sizefac = 1, alpha = 0.1, k = 4, p = 2,
		    c = Inf, distfun = "distmapfun", ..., y = NULL, matches = NULL,
		    object = NULL)

## S3 method for class 'FeatureMatchAnalyzer'
summary(object, ...)

## S3 method for class 'FeatureMatchAnalyzer'
plot(x, ..., type = c("all", "ph", "med", "msd", 
    "fom", "minsep", "cent.dist", "angle.diff", "area.ratio",
    "int.area", "bearing", "bdelta", "haus"))

## S3 method for class 'FeatureMatchAnalyzer'
print(x, ...)

FeatureComps(Y, X, which.comps=c("cent.dist", "angle.diff", "area.ratio", "int.area",
    "bdelta", "haus", "ph", "med", "msd", "fom", "minsep", "bearing"),
    sizefac=1, alpha=0.1, k=4, p=2, c=Inf, distfun="distmapfun", deg = TRUE,
    aty = "compass", loc = NULL, ...)

## S3 method for class 'FeatureComps'
distill(x, ...)

Arguments

x, y, matches	x, y and matches are list objects with components as output by deltamm or similar function. Only one is used, and it first checks for matches, then y, and finally x. It expects a component named mm.new.labels that gives the number of matched objects. In the case of the plot and print method functions, x is a list object as returned by FeatureMatchAnalyzer. distill: output from FeatureComps.
X, Y	list object giving a pixel image as output from solutionset from package spatstat for the verification and forecast fields, resp. These arguments are passed directly to the locperf function.
object	list object returned of class “FeatureMatchAnalyzer”, this is the returned value from the self-same function.
which.comps, type	character vector indicating which properties of the features are to be analyzed (which.comps) or plotted (type).
sizefac	single numeric by which area calculations should be multiplied in order to get the desired units. If unity (default) results are in terms of grid squares.
alpha	numeric value for the FOM measure (see the help file for locperf.
k	numeric indicating which quantile to use if the partial Hausdorff measure is to be used.
p	numeric giving the value of the parameter p for the Baddeley metric.
c	numeric giving the cut-off value for the Baddeley metric.
distfun	character naming a distance functions to use in calculating the various binary image measures. Default is Euclidean distance.
deg, aty	optional arguments to the bearing function.
loc	two-column matrix giving location coordinates for centroid distance. If NULL, uses an indices based on the dimension of the field.
...	Additional arguments to deltametric from package spatstat. In case of the summary method function, additional optional arguments may be passed, which include silent (logical, should the information be printed to the screen or not?), interest (numeric vector defining an interest value for calculating total interest for each matched object, if NULL, this is not performed), con (name of function that takes three arguments, the first two are matrices whose rows are objects and columns are matched feature properties, where the former is a matrix of matched feature property values (e.g., angle difference) and the latter is a matrix of interest values determined by the interest argument (whereby each row is identical), the third argument to con must be called which.comps, and it gives the short-form feature property names (i.e., same as which.comps argument); see details section). In the case of the plot method function, these are optional arguments to the function barplot. Not used by distill.

Details

FeatureMatchAnalyzer operates on objects of class “matched”. It is set up to calculate the values discussed in sec. 4 of Davis et al. (2006) for a single verification set (i.e., mean and standard deviation are not computed because it is only a single case). If criteria is 1, then features separated by a distance D < the sum of the sizes of the two features (size of a feature is defined as the square root of its area) are considered a match. If criteria is 2, then a match is made if D < the average of the sizes of the two features. Finally, criteria 3 decides a match as being anything less than a pre-determined constant.

FeatureComps is the primary function called by FeatureMatchAnalyzer, and is designed as a more stand-alone type of function. Several of the measures that can be calculated are simply the binary image measures/metrics available via, e.g., locperf. It calculates comparisons between two matched features (i.e., between the verification and forecast fields).

distill reduces a “FeatureComps” list object to a named numeric vector containing (in this order) the components that exist from "cent.dist", "angle.diff", "area.ratio", "int.area", "bdelta", "haus", "ph", "med", "msd", "fom", and "minsep". This is used, for example, by interester, which is why the order is important.

The summary method function for FeatureMatchAnalyzer allows for passing a function, con, to determine confidence for each interest value. The idea being to set the interest to zero when the particular interest value does not make sense. For example, angle difference makes no sense if both objects are circles. Currently, no functions are included in this package for actually doing this, and so the functionality itself has not been tested.

The print method function for FeatureMatchAnalyzer first converts the object to a simple named matrix, then prints the matrix out. The resulting matrix is returned invisibly.

Value

FeatureMatchAnalyzer returns a list of list objects. The specific components depend on the 'which.comps' argument, and are the same as those returned by FeatureComps. These can be any of the following.

cent.dist	numeric giving the centroid (Euclidean) distance.
angle.diff	numeric giving the orientation (major axis) angle difference.
area.ratio	numeric giving the area ratio, which is always between 0 and 1 because this is defined by Davis et al. (2006) to be the area of the smaller feature divided by that of the larger feature regardless of which field the feature belongs to.
int.area	numeric giving the intersection area of the features.
bdelta	numeric giving Baddeley's delta metric between the two features.
haus, ph, med, msd, fom, minsep	numeric, see locperf for specific information.
bearing	numeric giving the bearing from the forecast object centroid to the observed object centroid.

The summary method for FeatureMatchAnalyzer invisibly returns a matrix with the same information, but where each matched object is a row and each column is the specific statistic. Or, if optional interest argument is passed, a list with components:

print returns a named vector invisibly.

Author(s)

Eric Gilleland

References

Davis, C. A., Brown, B. G. and Bullock, R. G. (2006) Object-based verification of precipitation forecasts, Part I: Methodology and application to mesoscale rain areas. Mon. Wea. Rev., 134, 1772–1784.

Examples

data( "ExampleSpatialVxSet" )

x <- ExampleSpatialVxSet$vx
xhat <- ExampleSpatialVxSet$fcst

hold <- make.SpatialVx( x, xhat, field.type="Example",
    units = "units", data.name = "Example", 
    obs.name = "x", model.name = "xhat" )

look <- FeatureFinder(hold, smoothpar=1.5)
look2 <- centmatch(look)

tmp <- FeatureMatchAnalyzer(look2)
tmp
summary(tmp)
plot(tmp)

---

Single Feature Properties

Description

Calculate properties for an identified feature.

Usage

FeatureProps(x, Im = NULL, which.props = c("centroid", "area", "axis", "intensity"),
    areafac = 1, q = c(0.25, 0.9), loc = NULL, ...)

Arguments

x	object of class “owin” containing a binary image matrix defining the feature.
Im	Matrix giving the original values of the field from which the feature was extracted. Only needed if the feature intensity is desired.
which.props	character vector giving one or more of “centroid”, “area”, “axis” and “intensity”. If “axis” is given, then a call to FeatureAxis is made.
areafac	numeric, in determining the lengths of the axes, they are multiplied by a factor of fac (e.g., if the grid points are k by k km each, then one could set this to k so that the resulting lengths are in terms of km rather than grid points.
q	numeric vector of values between 0 and 1 inclusive giving the quantiles for determining the intensity of the feature.
loc	optional argument giving a two-column matrix of grid locations for finding the centroid. If NULL, indices based on the dimension of x are used.
...	additional arguments to FeatureAxis.

Details

This function takes an owin image and returns several property values for that image, including: centroid, spatial area, major and minor axis angle/length, as well as the overall intensity of the field (cf., Davis et al., 2006a, b).

Value

list object with components depending on the which.props argument. One or more of:

centroid	list with components x and y giving the centroid of the object.
area	numeric giving the area of the feature.
axis	list object of class FeatureAxis as returned by the same-named function.

Author(s)

Eric Gilleland

References

Davis, C. A., Brown, B. G. and Bullock, R. G. (2006a) Object-based verification of precipitation forecasts, Part I: Methodology and application to mesoscale rain areas. Mon. Wea. Rev., 134, 1772–1784.

Davis, C. A., Brown, B. G. and Bullock, R. G. (2006b) Object-based verification of precipitation forecasts, Part II: Application to convective rain systems. Mon. Wea. Rev., 134, 1785–1795.

Examples

data( "ExampleSpatialVxSet" )

x <- ExampleSpatialVxSet$vx

look <- disk2dsmooth(x,5)
u <- quantile(look,0.99)
sIx <- matrix(0, 100, 100)
sIx[ look > u] <- 1
look2 <- disjointer(sIx)[[1]]
look2 <- flipxy(look2)

FeatureProps(look2,
    which.props=c("centroid", "area", "axis"))

---

Feature-Based Contingency Table

Description

Create a feature-based contingency table from a matched object and calculate some summary scores with their standard errors.

Usage

FeatureTable(x, fudge = 1e-08, hits.random = NULL, correct.negatives = NULL, fA = 0.05)

## S3 method for class 'FeatureTable'
ci(x, alpha = 0.05, ...)

## S3 method for class 'FeatureTable'
print(x, ...)

## S3 method for class 'FeatureTable'
summary(object, ...)

Arguments

x	FeatureTable: An object of class “matched” (e.g., from deltamm or centmatch. print: An object of class “FeatureTable”.
object	An object of class “FeatureTable”.
fudge	value added to denominators of scores to ensure no division by zero. Set to zero if this practice is not desired.
hits.random	If a different value for random hits in the caluclation of GSS than provided is desired, it can be given here. Default uses Eq (3) from Davis et al (2009).
correct.negatives	If a different value for correct negatives than provided is desired, it can be given here. Default uses Eq (4) from Davis et al (2009).
fA	numeric between zero and 1 giving the fraction of area occupied for the purpose of matching as in Davis et al (2009).
alpha	numeric between zero and one giving the (1 - alpha) * 100 percent confidence level.
...	Additional arguments to ci.

Details

This function takes an object of class “matched” and calculates a contingency table based on matched and unmatched objects. If no value for correct negatives is given, then it will also determine them based on Eq (3) from Davis et al (2009). The following contingency table scores and their standard errors (based on their usual traditional version) are returned. It should be noted that the standard errors may not be entirely meaningful because they do not capture the uncertainty associated with identiying, merging and matching features within the fields. Neverhteless, they are calculated here for investigative purposes. Note that hits are determined by number of matched objects, which for some matching algorithms can mean that features are matched more than once (e.g., if using centmatch). In essence, this fact may artificially increase thenumber of hits. On the other hand, situations exist where such handling may be more appropriate than not having duplicate matches.

hits are determined by the total number of matched features.

false alarms are the total number of unmatched forecast features.

misses are the total number of unmatched observed features.

correct negatives are less obviously defined. If the user does not supply a value, then these are calculated according to Eq (4) in Davis et al (2009).

GSS: Gilbert skill score (aka Equitable Threat Score) based on Eq (2) of Davis et al (2009).

POD: probability of detecting an event (aka the hit rate).

false alarm rate: (aka probability of false detection) is the ratio of false alarms to the number of false alarms and correct negtives.

FAR: the false alarm ratio is the ratio of false alarms to the total forecast events (in this case, the total number of forecast features in the field).

HSS: Heidke skill score

The print method function simply calls summary, which prints the feature-based contingency table in addition to calling ci. The confidence intervals are based on the normal approximation method using the estimated standard errors, which themselves are suspicious. In any case, the intervals can give a feel for some of the uncertainty associated with the scores, but should not be considered as solid.

Value

A list with inherited attributes from x and components:

estimates	named numeric vector giving the estimated scores.
se	named numeric vector giving the estimated standard errors of the scores.
feature.contingency.table	named numeric vector giving the feature-based contingency table.

Note

Standard error estimates are based on the univariate equivalent formulations, which do not account for uncertainties introduced in the feature identification, merging/clustering and matching. They should not be considered as legitimate, and resulting confidence intervals should be mistrusted.

Author(s)

Eric Gilleland

References

Davis, C. A., Brown, B. G., Bullock, R. G. and Halley Gotway, J. (2009) The Method for Object-based Diagnostic Evaluation (MODE) applied to numerical forecasts from the 2005 NSSL/SPC Spring Program. Wea. Forecsting, 24, 1252–1267, DOI: 10.1175/2009WAF2222241.1.

See Also

To identify features in the fields: FeatureFinder

To match (and merge) features: centmatch, deltamm

Examples

##
## See help file for 'deltamm' for examples.
##

---

2-d Interpolation

Description

Interpolate a function of two variables by rounding (i.e. taking the nearest value), bilinear or bicubic interpolation.

Usage

Fint2d(X, Ws, s, method = c("round", "bilinear", "bicubic"), derivs = FALSE, ...)

Arguments

X	A numeric n by m matrix giving the value of the function at the old coordinates.
Ws	A numeric k by 2 matrix of new grid coordinates where k <= m * n.
s	A numeric k by 2 matrix of old grid coordinates where k <= m * n.
method	character naming one of “round” (default), “bilinear”, or “bicubic” giving the specific interpolation method to use.
derivs	logical, should the gradient interpolatants be returned?
...	Not used.

Details

Method round simply returns the values at each grid point that correspond to the nearest points in the old grid.

Interpolation of a function, say H, is achieved by the following formula (cf. Gilleland et al 2010, sec. 3), where r and s represent the fractional part of their respective coordinate. that is, r = x - g( x ) and s = y - g( y ), where g( x ) is the greatest integer less than x.

sum_k sum_l b_k( r ) * b_l( s ) * H(g( x ) + l, g( y ) + k).

The specific choices for the values of b_l and b_k and their ranges depends on the type of interpolation. For bilinear interpolation, they both range from 0 to 1, and are given by: b_0( x ) = 1 - x and b_1( x ) = x. for bicubic interpolation, they both range from -1 to 2 and are given by:

b_(-1)( t ) = (2 * t^2 - t^3 - t) / 2

b_(0)( t ) = (3 * t^3 - 5 * t^2 + 2) / 2

b_(1)( t ) = (4 * t^2 - 3 * t^3 + t) / 2

b_(2)( t ) = ((t - 1) * t^2) / 2.

Value

If deriv is FALSE, then a matrix is returned whose values correspond to the new coordinates. Otherwise a list is returned with components:

xy	matrix whose values correspond to the new coordinates.
dx, dy	matrices giving the x and y direction gradients of the interpolation.

Author(s)

Eric Gilleland

References

Gilleland and co-authors (2010) Spatial forecast verification: Image warping. NCAR Technical Note, NCAR/TN-482+STR, DOI: 10.5065/D62805JJ.

See Also

rigider, rigidTransform

Examples

# see rigider for an example.

---

Forecast Quality Index

Description

Functions for calculating the Forecast Quality Index (FQI) and its components.

Usage

FQI(object, surr = NULL, k = 4, time.point = 1, obs = 1, model = 1, ...)

UIQI(X, Xhat, ...)

ampstats(X, Xhat, only.nonzero = FALSE)

## S3 method for class 'fqi'
print(x, ...)

## S3 method for class 'fqi'
summary(object, ...)

Arguments

object	list object of class “SpatialVx”. In the case of the summary method, object is the list object returned by FQI.
X, Xhat	numeric matrices giving the fields for the verification set.
x	list object of class “fqi” as returned by FQI.
surr	three-dimesnional array containing surrogate fields for X, e.g. as returned by surrogater2d. If NULL, these will be calculated using surrogater2d.
only.nonzero	logical, should the means and variances of only the non-zero values of the fields be calculated (if so, the covariance is returned as NA)?
k	numeric vector for use with the partial Hausdorff distance. For k that are whole numerics or integers >= 1, then the k-th highest value is returned by locmeasures2d. If 0 <= k < 1, then the corresponding quantile is returned.
time.point	numeric or character indicating which time point from the “SpatialVx” verification set to select for analysis.
obs, model	numeric indicating which forecast model to select for the analysis.
...	In the case of FQI, additional arguments to surrogater2d. Only used if surr is NULL. In the case of UIQI, additional arguments to ampstats. In the case of summary.fqi, these are not used.

Details

The FQI was proposed as a spatial verification metric (a true metric in the mathematical sense) by Venugopal et al. (2005) to combine amplitude and displacement error information in a single summary statistic. It is given by

FQI = (PHD_k(X, Xhat)/mean( PHD_k(X, surr_i); i in 1 to number of surrogates)) / (brightness * distortion)

where the numerator is a normalized partial Hausdorff distance (see help file for locperf), brightness (also called bias) is given by 2*(mu1*mu2)/(mu1^2+mu2^2), where mu1 (mu2) is the mean value of X (Xhat), and the distortion term is given by 2*(sig1*sig2)/(sig1^2+sig2^2), where sig1^2 (sig2^2) is the variance of X (Xhat) values. The denominator is a modified UIQI (Universal Image Quality Index; Wang and Bovik, 2002), which itself is given by

UIQI = cor(X,Xhat)*brightness*distortion.

Note that if only.nonzero is TRUE in the call to UIQI, then the modified UIQI used in the FQI formulation is returned (i.e., without multiplying by the correlation term).

The print method so far just calls the summary method.

Value

FQI returns a list with with the following components:

phd.norm	matrix of normalized partial Hausdorff distances for each value of k (rows) and each threshold (columns).
uiqi.norm	numeric vector of modified UIQI values for each threshold.
fqi	matrix of FQI values for each value of k (rows) and each threshold (columns).

It will also have the same attributes as the “SpatialVx” object with additional attributes defining the arguments specific to parameters used by the function.

UIQI returns a list with components:

data.name	character vector giving the names of the two fields.
cor	single numeric giving the correlation between the two fields.
brightness.bias	single numeric giving the brightness (bias) value.
distortion.variability	single numeric giving the distortion (variability) value.
UIQI	single numeric giving the UIQI (or modified UIQI if only.nonzero is set to TRUE) value.

ampstats returns a list object with components:

mean.fcst, mean.vx	single numerics giving the mean of Xhat and X, resp.
var.fcst, var.vx	single numerics giving the variance of Xhat and X, resp.
cov	single numeric giving the covariance between Xhat and X (if only.nonzero is TRUE, this will be NA).

Author(s)

Eric Gilleland

References

Venugopal, V., Basu, S. and Foufoula-Georgiou, E. (2005) A new metric for comparing precipitation patterns with an application to ensemble forecasts. J. Geophys. Res., 110, D08111, 11 pp., doi:10.1029/2004JD005395.

Wang, Z. and Bovik, A. C. (2002) A universal image quality index. IEEE Signal Process. Lett., 9, 81–84.

See Also

locperf, surrogater2d, locmeasures2d

Examples

data( "ExampleSpatialVxSet" )

 x <- ExampleSpatialVxSet$vx
 xhat <- ExampleSpatialVxSet$fcst

 # Now, find surrogates of the simulated field.
 z <- surrogater2d(x, zero.down=TRUE, n=10)

 u <- list( X = cbind( quantile( c(x), c(0.75, 0.9)) ),
	Xhat = cbind( quantile( c(xhat), c(0.75, 0.9) ) ) )

hold <- make.SpatialVx(x, xhat, thresholds = u,
			field.type = "Example", units = "none",
			data.name = "ExampleSpatialVxSet",
			obs.name = "X", model.name = "Xhat" )

FQI(hold, surr = z, k = c(4, 0.75) )

---

Various Verification Statistics on Possibly Neighborhood-Smoothed Fields.

Description

Functions to calculate various verification statistics on possibly neighborhood smoothed fields. Used by hoods2d, but can be called on their own.

Usage

fss2dfun(sPy, sPx, subset = NULL, verbose = FALSE)

fuzzyjoint2dfun(sPy, sPx, subset = NULL)

MinCvg2dfun(sIy, sIx, subset = NULL)

multicon2dfun(sIy, Ix, subset = NULL)

pragmatic2dfun(sPy, Ix, mIx = NULL, subset = NULL)

upscale2dfun(sYy, sYx, threshold = NULL, which.stats = c("rmse",
                 "bias", "ts", "ets"), rule = ">=", subset = NULL)

Arguments

sPy	n by m matrix giving a smoothed binary forecast field.
sPx	n by m matrix giving a smoothed binary observed field.
sIy	n by m matrix giving a binary forecast field.
sIx	n by m matrix giving a binary observed field (the s indicates that the binary field is obtained from a smoothed field).
Ix	n by m matrix giving a binary observed field.
mIx	(optional) single numeric giving the base rate. If NULL, this will be calculated by the function. Simply a computation saving step if this has already been calculated.
sYy	n by m matrix giving a smoothed forecast field.
sYx	n by m matrix giving a smoothed observed field.
threshold	(optional) numeric vector of length 2 giving the threshold over which to calculate the verification statistics: bias, ts and ets. If NULL, only the rmse will be calculated.
which.stats	character vector naming which statistic(s) should be caluclated for upscale2dfun.
subset	(optional) numeric indicating over which points the summary scores should be calculated. If NULL, all of the points are used.
rule	character string giving the sort of thresholding process desired. See the help file for thresholder for more information.
verbose	logical, should progress information be printed to the screen?

Details

These are modular functions that calculate the neighborhood smoothing method statistics in spatial forecast verification (see, e.g., Ebert, 2008, 2009; Gilleland et al., 2009, 2010; Roberts and Lean,2008). These functions take fields that have already had the neighborhood smoothing applied (e.g., using kernele2d) when appropriate. They are called by hoods2d, so need not be called by the user, but they can be.

Value

In the case of fss2dfun, a single numeric giving the FSS value is returned. In the other cases, list objects are returned with one or more of the following components, depending on the particular function.

fuzzy	fuzzyjoint2dfun returns a list with this list as one component. The list component fuzzy has the components: pod, far and ets.
joint	fuzzyjoint2dfun returns a list with this list as one component. The list component joint has the components: pod, far and ets.
pod	numeric giving the probability of detection, or hit rate.
far	numeric giving the false alarm ratio.
ets	numeric giving the equitable threat score, or Gilbert Skill Score.
f	numeric giving the false alarm rate.
hk	numeric giving the Hanssen-Kuipers statistic.
bs	Brier Score
bss	Brier Skill Score. The pragmatic2dfun returns the bs and bss values. The Brier Skill Score here uses the mean square error between the base rate and the Ix field as the reference forecast.
ts	numeric giving the threat score.
bias	numeric giving the frequency bias.

Author(s)

Eric Gilleland

References

Ebert, E. E. (2008) Fuzzy verification of high resolution gridded forecasts: A review and proposed framework. Meteorol. Appl., 15, 51–64. doi:10.1002/met.25

Ebert, E. E. (2009) Neighborhood verification: A strategy for rewarding close forecasts. Wea. Forecasting, 24, 1498–1510, doi:10.1175/2009WAF2222251.1.

Gilleland, E., Ahijevych, D., Brown, B. G., Casati, B. and Ebert, E. E. (2009) Intercomparison of Spatial Forecast Verification Methods. Wea. Forecasting, 24, 1416–1430, doi:10.1175/2009WAF2222269.1.

Gilleland, E., Ahijevych, D. A., Brown, B. G. and Ebert, E. E. (2010) Verifying Forecasts Spatially. Bull. Amer. Meteor. Soc., October, 1365–1373.

Roberts, N. M. and Lean, H. W. (2008) Scale-selective verification of rainfall accumulations from high-resolution forecasts of convective events. Mon. Wea. Rev., 136, 78–97. doi:10.1175/2007MWR2123.1.

Examples

x <- y <- matrix( 0, 100, 100)
x[ sample(1:100, 10), sample(1:100, 10)] <- 1
y[ sample(1:100, 20), sample(1:100, 20)] <- 1
Px <- kernel2dsmooth( x, kernel.type="boxcar", n=9, xdim=c(100, 100))
Py <- kernel2dsmooth( y, kernel.type="boxcar", n=9, xdim=c(100, 100))
par( mfrow=c(2,2))
image( x, col=c("grey", "darkblue"), main="Simulated Observed Events")
image( y, col=c("grey", "darkblue"), main="Simulated Forecast Events")
image( Px, col=c("grey", tim.colors(256)), main="Forecast Event Frequencies (9 nearest neighbors)")
image( Py, col=c("grey", tim.colors(256)), main="Smoothed Observed Events (9 nearest neighbors)")
fss2dfun( Py, Px)

---

Create Several Graphics for List Objects Returned from hoods2d

Description

Creates several graphics for list objects returned from hoods2d. Mostly quilt and matrix plots for displaying results of smoothing fields over different neighborhood lengths and thresholds.

Usage

fss2dPlot(x, ..., matplotcol = 1:6, mfrow = c(1, 2), add.text = FALSE)

upscale2dPlot(object, args, ..., type = c("all",
                 "gss", "ts", "bias", "rmse"))

Arguments

x	list object with components fss, fss.random and fss.uniform. Effectively, it does the same thing as hoods2dPlot, but adds the fss.random and fss.uniform horizontal lines to the matrix plot.
object	list object with named components: rmse (numeric vector), ets, ts and bias all matrices whose rows represent neighborhood lengths, and whose columns represent thresholds.
args	list object passed to hoods2dPlot, see its help file for more details.
mfrow	mfrow parameter (see help file for par). If NULL, then the parameter is not re-set.
add.text	logical, if TRUE, FSS values will be added to the quilt plot as text (in addition to the color).
type	character string stating which plots to make (default is “all”).
...	Optional arguments to image and image.plot for fss2dPlot, and optional arguments to hoods2dPlot for upscale2dPlot
matplotcol	col argument to function matplot.

Details

makes quilt and matrix plots for output from hoods2d.

Value

No value is returned. A series of plots are created. It may be useful to use this function in conjunction with pdf in order to view all of the plots. See the help file for hoods2dPlot to plot individual results.

Author(s)

Eric Gilleland

Examples

##
## This is effectively an internal function, so the example is commented out
## in order for R's check to run faster.
##
## Not run: 
data( "geom001" )
data( "geom000" )
data( "ICPg240Locs" )

hold <- make.SpatialVx( geom000, geom001, thresholds = c(0.01,50.01),
    loc = ICPg240Locs, map = TRUE, projection = TRUE, loc.byrow = TRUE,
    units = "mm/h", data.name = "Geometric", obs.name = "observation",
    model.name = "case 1" )

look <- hoods2d(hold, levels=c(1, 3, 5, 33, 65),
    verbose=TRUE)
plot( look)

## End(Not run)

---

Spatial-Alignment Summary Measures

Description

Calculates the spatial alignment summary measures from Gilleland (2020)

Usage

Gbeta(X, Xhat, threshold, beta, alpha = 0, rule = ">", ...)

GbetaIL(X, Xhat, threshold, beta, alpha = 0, rule = ">", w = 0.5, ...)

G2IL(X, Xhat, threshold, beta, alpha = 0, rule = ">", ...)

Arguments

X, Xhat	Observed and Forecast fields in the form of a matrix. X and Xhat must have the same dimensions. May be class “owin” objects for Gbeta.
threshold	If X and Xhat are not binary “owin” objects, then the threshold is used to define the binary field. See binarizer for more information about this argument.
beta	single numeric defining the upper limit for y = y1y2 or y = y1y2(1+y3) determining the rate of decrease in the measure. Default is half the domain size squared.
alpha	single numeric defining a lower limit for y = y1y2 or y = y1y2(1+y3) determining what constitutes a perfect match. The default of zero requires the two fields to be identical, and is probably what is wanted in most cases.
rule	See binarizer for more information about this argument.
w	single numeric between zero and one describing how much weight to give to the first term in the definition of GbetaIL. See details section below.
...	Not used.

Details

These summary measures were proposed in Gilleland (2020) and provide an index between zero (bad score) and one (perfect score) describing the closeness in spatial alignment between the two fields. GbetaIL and G2IL also incorporate a distributional summary of the intensity errors. Gbeta applied between two fields A and B is defined as

Gbeta(A,B) = max( 1 - y/beta, 0),

where y = y1 * y2, with y1 a measure of the overlap between A and B (if they overlap completely, y1 = 0); it is the number of points in AB^c and A^cB, where ^c denotes set complement. The term y2 = (MED(A,B) * nB + MED(B,A) * nA), where MED is the mean-error distance (see locperf for more about MED), with nA and nB representing the number of 1-valued grid points in the sets A and B, resp. If alpha != 0, then the term y/beta is replaced with (y - alpha) / (beta - alpha).

GbetaIL is defined to be

GbetaIL(A,B) = w * Gbeta + (1-w) * theta(A,B),

where theta is the maximum of zero and the linear correlation coefficient between the intensity values after having sorted them; this part is carried out via a call to qqplot. If the number of points in the two sets differs (i.e., if nA != nB), then the larger set is linear interpolated to be the same size as the smaller set. If both fields are empty, theta = 1. If field A is empty, then theta = 1 - (nB / N), where N is the size of the domain. Similarly, if field B is empty. The rationale is that if one field is empty and the other has very few nonzero points, then the two fields are more similar than if the other field has many points.

G2betaIL(A,B) = max(1 - (y1 * y2 *(1 + y3 ))/beta, 0),

where y1 and y2 are as above and y3 is the mean-absolute difference between the sorted values from the sets A and B analogous as for GbetaIL. If nA = nB = 0, then y3 = 0. If only one of nA or nB is zero, then the absolute value of the maximum intensity of the other set is used. The rationale being that this value represents the most egregious error so that if it is large, then the difference is penalized more.

See Gilleland (2020) for more details about these measures.

Value

An object of class “Gbeta” giving a single numeric giving the value of the summary (index) measures is returned, but with additional attributes that can be obtained using the attributes function, and are also displayed via the print function. To remove the attributes, the function c can be used. For Gbeta these include:

components	numeric vector giving nA, nB, nAB (the number of points in both sets), nA + nB - 2nAB (i.e., the overlap measure, y1), medAB, medBA, medAB * nnB, medBA * nA, and the two asymmetric versions of this measure (see Gilleland, 2020).
beta	The value of beta used.
alpha	The value of alpha used.
threshold	numeric vector giving the value of the threshold and the rule used.
data.name	character vector giving the object names used for X and Xhat.

The attributes for GbetaIL and G2IL are the same, but the components vector also includes the value of theta. GbetaIL has an additional attribute called weights that gives w and 1 - w.

Author(s)

Eric Gilleland

References

Gilleland, E. (2020) Novel measures for summarizing high-resolution forecast performance. Advances in Statistical Climatology, Meteorology and Oceanography, 7 (1), 13–34, doi: 10.5194/ascmo-7-13-2021.

See Also

TheBigG,qqplot, binarizer, locperf

Examples

data( "obs0601" )
data( "wrf4ncar0531" )
res <- Gbeta( X = obs0601, Xhat = wrf4ncar0531, threshold = 2.1 )
c( res )
res
attributes( res )

GbetaIL( X = obs0601, Xhat = wrf4ncar0531, threshold = 2.1,
	beta = 601 * 501 )

G2IL( X = obs0601, Xhat = wrf4ncar0531, threshold = 2.1,
     beta = 601 * 501 )

---

Geographic Box Plot

Description

Make a geographic box plot as detailed in Willmott et al. (2007).

Usage

GeoBoxPlot(x, areas, ...)

Arguments

x	numeric giving the values to be box-plotted.
areas	numeric of same length as x giving the associated areas for each value.
...	optional arguments to the boxplot function of R. The argument plot is not allowed.

Details

This function makes the geographic box plots described in Willmott et al. (2007) that calculates the five statistics in such a way as to account for the associated areas (e.g., over a grid where each grid box may have differing areas).

Missing values are not handled, and ideally should be handled before calling ths routine.

In future, this function may allow other options for x than currently, but for now, only numeric vectors are allowed.

Value

List with the same components as returned by boxplot.

Author(s)

Eric Gilleland

References

Willmott, C. J., Robeson, S. M. and Matsuura, K. (2007) Geographic box plots. Physical Geography, 28, 331–344, doi:10.2747/0272-3646.28.4.331.

See Also

boxplot

Examples

##
## Reproduce the boxplots of Fig. 1 in Willmott et al. (2007).
##
x <- c(4,9,1,3,10,6,7)
a <- c(rep(1,4),2,1,3)
boxplot( x, at=1, xlim=c(0,3))
GeoBoxPlot(x, a, at=2, add=TRUE)
axis( 1, at=c(1,2), labels=c("Traditional", "Geographic"))

---

Example Verification Set

Description

Example verification set of accumulated precipitation (mm) with 361 time points in addition to 2352 spatial locations on a grid. Taken from a real, but unknown, weather model and observation analysis (one of GFS or NAM). Accumulation is either 3-h or 24-h.

Usage

data(GFSNAMfcstEx)
data(GFSNAMobsEx)
data(GFSNAMlocEx)

Format

The format is: num [1:2352, 1:361] 0 0 0 0 0 ...

The format is: num [1:2352, 1:361] 0 0 0 0 0 ...

The format is: A data frame with 2352 observations on the following 2 variables.

Lat

a numeric vector of latitude coordinates for GFS/NAM example verification set.

Lon

a numeric vector of longitude coordinates for GFS/NAM example verification set.

Details

Example verification set with 2352 spatial locations over the United States, and 361 time points. For both the forecast (GFSNAMfcstEx) and verification (GFSNAMobsEx), these are numeric matrices whose rows represent time, and columns represent space. The associated lon/lat coordinates are provided by GFSNAMlocEx (2352 by 2 data frame with named components giving the lon and lat values).

Note that the available spatial locations are a subset of the original 70 X 100 60-km grid where each time point had no missing observations. This example set is included with the package simply toi demonstrate some functionality that involves both space and time; though this is mostly a spatial-only package.

Examples

data( "GFSNAMfcstEx" )
data( "GFSNAMobsEx" )
data( "GFSNAMlocEx" )

x <- colMeans(GFSNAMfcstEx, na.rm=TRUE)
y <- colMeans(GFSNAMobsEx, na.rm=TRUE)
look <- as.image(x - y, x=GFSNAMlocEx)
image.plot(look)

---

2-d Gaussian Mixture Models Verification

Description

Use 2-d Gaussian Mixture Models (GMM) to assess forecast performance.

Usage

gmm2d(x, ...)

## Default S3 method:
gmm2d(x, ..., xhat, K = 3, gamma = 1, threshold = NULL,
    initFUN = "initGMM", verbose = FALSE)

## S3 method for class 'SpatialVx'
gmm2d(x, ..., time.point = 1, obs = 1, model = 1, K = 3, gamma = 1,
    threshold = NULL, initFUN = "initGMM", verbose = FALSE)

## S3 method for class 'gmm2d'
plot(x, ..., col = c("gray", tim.colors(64)),
    zlim = c(0, 1), horizontal = TRUE)

## S3 method for class 'gmm2d'
predict(object, ..., x)

## S3 method for class 'gmm2d'
print(x, ...)

## S3 method for class 'gmm2d'
summary(object, ...)

Arguments

x, xhat	Default: m by n numeric matrices giving the verification and forecast fields, resp. gmm2d.SpatialVx: object of class “SpatialVx”. plot and print: object returned by gmm2d. predict: k by 2 matrix of lon/lat coordinates on which to predict the model.
object	output from gmm2d.
K	single numeric giving the number of mixture components to use.
gamma	Value of the gamma parameter from Eq (11) of Lakshmanan and Kain (2010). This affects the number of times a location is repeated.
threshold	numeric giving a threshold over which (and including) the GMM is to be fit (zero-valued grid points are not included in the estimation here for speed). If NULL, no thresholding is applied.
initFUN	character naming a function to provide initial estimates for the GMM. Must take an m by n matrix as input, and return a dataframe a component called ind that is a vector indicating the order of the rows for which the first K will be used, a third column giving the x-coordinates of the initial estimate of the mean for the x direction, fourth column giving the initial estimate for the mean of the y-direction, and fifth and sixth columns giving initial estimates for the standard deviations of the x- and y-directions. The default identifies all connected components using the disjointer function, then uses their centroids as the initial estimates of the means, and their axes as initial estimates for the standard deviations. The ind component gives the order of the object areas from largest to smallest so that the K largest objects are used to provide initial estimates. Note that this differs from the initial estimates in Lakshmanan and Kain (2010) where they break the field into different areas first.
verbose	logical, should progress information be printed to the screen?
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.
col, zlim, horizontal	optional arguments to fields function(s) poly.image, image.plot.
...	In the case of gmm2d: optional arguments to initFUN. In the case of plot: not used. In the case of predict: N by 2 matrix of grid point locations on which to predict the probability from the 2-d GMM model. In the case of summary: this can include the arguments: silent, logical stating whether to print summaries to the screen (FALSE) or not (TURE), e1, e2, ..., e5, giving alternative weights in calculating the overall error (Eq 15 in Lakshmanan and Kain, 2010, but see details section below).

Details

These functions carry out the spatial verification approach described in Lakshmanan and Kain (2010), which fits a 2-d Gaussian Mixture Model (GMM) to the locations for each field in the verification set, and makes comparisons using the estimated parameters. In fitting the GMMs, first an initial estimate is provided by using the initFUN argument, which is a function. The default function is relatively fast (it might seem slow, but for what it does, it is very fast!), but is typically the slowest part of the process. Although the EM algorithm is a fairly computationally intensive procedure, acceleration algorithms are employed (via the turboem function of the turboEM package) so that once initial estimates are found, the procedure is very fast.

Because the fit is to the locations only, Lakshmanan and Kain (2010) suggest two ways to incorporate intensity information. The first is to repeat points with higher intensities, and the second is to multiply the results by the total intensities over the fields. The points are repeated M times according to the formula (Eq 11 in Lakshmanan and Kain, 2010):

M = 1 + gamma * round( CFD(I_(xy))/frequency(I_MODE)),

where CFD is the cumulative *frequency* distribution (here estimated from the histogram using the ‘hist’ function), I_(xy) is intensity at grid point (x,y), I_MODE is the mode of intensity values, and gamma is a user-supplied parameter controlling how much to repeat points where higher numbers will result in larger repetitions of high intensity values.

The function gmm2d fits the 2-d GMM to both fields, plot.gmm2d first uses predict.gmm2d to obtain probabilities for each grid point, and then makes a plot similar to those in Lakshmanan and Kain (2010) Figs. 3, 4 and 5, but giving the probabilities instead of the probabilities times A. Note that predict.gmm2d can be very slow to compute so that plot.gmm2d can also be very slow. Less effort was put into speeding these functions up because they are not necessary for obtaining results via the parameters. However, they can give the user an idea of how good the fit is.

The 2-d GMM is given by

G(x,y) = A*sum(lambda*f(x,y))

where lambda and f(x,y) are numeric vectors of length K, lambda components describe the mixing, and f(x,y) is the bivariate normal distribution with mean (mu.x, mu.y) and covariance function. ‘A’ is the total sum of intensities over the field.

Comparisons between forecast and observed fields are carried out finally by the summary method function. In particular, the translation error

e.tr = sqrt((mu.xf - mu.xo)^2 + (mu.yf - mu.yo)^2),

where f means forecast and o verification fields, resp., and mu .x is the mean in the x- direction, and mu.y in the y- direction. The rotation error is given by

e.rot = (180/pi)*acos(theta),

where theta is the dot product between the first eigenvectors of the covariance matrices for the verification and forecast fields. The scaling error is given by

e.sc = Af*lambda.f/Ao*lambda.o,

where lambda is the mixture component and Af/Ao is the forecast/observed total intensity.

The overall error (Eq 15 of Lakshmanana and Kain, 2010) is given by

e.overall = e1 * min(e.tr/e2, 1) + e3*min(e.rot,180 - e.rot)/e4 + e5*(max(e.sc,1/e.sc)-1),

where e1 to e5 can be supplied by the user, but the defaults are those given by Lakshmanan and Kain (2010). Namely, e1 = 0.3, e2 = 100, e3=0.2, e4 = 90, and e5=0.5.

Value

For gmm2d, a list object of class “gmm2d” is returned with components:

fitX, fitY	list objects returned by the turboem function from the turboEM package that describe the EM estimates of the 2-d GMM parameters for the verification and forecast fields, resp.
initX, initY	numeric vectors giving the initial estimates used in the EM algorithm for the verification and forecast fields, resp. The first 2*K values are the initial mean estiamtes for the x- and y- directions, resp. The next 4*K values are the initial estiamtes of the covariances (note that the cross-covariance terms are zero regardless of initialization function employed (maybe this will be improved in the future). The final K values are the initial estimates for lambda.
sX, sY	N by 2 matrix giving the repeated coordinates calculated per M as described in the details section for the verification and forecast fields, resp.
k	single numeric giving the value of K
Ax, Ay	single numerics giving the value of A (the total sum of intensities over the field) for the verifiaction and forecast fields, resp.

For plot.gmm2d no value is returned. A plot is created.

For predict.gmm2d, a list is returned with components:

predX, predY

numeric vectors giving the GMM predicted values for the verification and forecast fields, resp.

For summary.gmm2d, a list is returned invisibly (if silent is FALSE, information is printed to the screen) with components:

meanX, meanY	Estimated mean vectors for each GMM component for the verification and forecast fields, resp.
covX, covY	Estimated covariances for each GMM component for the verification and forecast fields, resp.
lambdasX, lambdasY	Estimated mixture components for each GMM component for the verification and forecast fields, resp.
e.tr, e.rot, e.sc, e.overall	K by K matrices giving the errors between each GMM component in the verification field (rows) to each GMM component in the forecast field (columns). The errors are: translation (e.tr), rotation (e.rot), scaling (e.sc), and overall (e.overall).

Author(s)

Eric Gilleland

References

Lakshmanan, V. and Kain, J. S. (2010) A Gaussian Mixture Model Approach to Forecast Verification. Wea. Forecasting, 25 (3), 908–920.

Examples

## Not run: 
data( "ExampleSpatialVxSet" )

x <- ExampleSpatialVxSet$vx
xhat <- ExampleSpatialVxSet$fcst

u <- min(quantile(c(x[x > 0]), probs = 0.75),
    quantile(c(xhat[xhat > 0]), probs = 0.75))

look <- gmm2d(x, xhat=xhat, threshold=u, verbose=TRUE)
summary(look)
plot(look)

## End(Not run)
## Not run: 
# Alternative method to skin the cat.
hold <- make.SpatialVx( x, xhat, field.type = "MV Gaussian w/ Exp. Cov.",
    units = "units", data.name = "Example", obs.name = "x",
    model.name = "xhat" )

look2 <- gmm2d( hold, threshold = u, verbose = TRUE)
summary(look2)
plot(look2)


## End(Not run)

---

Variograms for a Gridded Verification Set

Description

Find (and plot) variograms for each field in a gridded verification set.

Usage

griddedVgram(object, zero.in = TRUE, zero.out = TRUE, time.point = 1, 
    obs = 1, model = 1, ...)

## S3 method for class 'griddedVgram'
plot(x, ... )

Arguments

object	list object of class “SpatialVx” containing information on the verification set.
zero.in, zero.out	logical, should the variogram be calculated over the entire field (zero.in), and/or over only the non-zero values (zero.out)?
x	list object as returned by griddedVgram.
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.
...	In the case of griddedVgram, these are optional arguments to the vgram.matrix function from package fields. In the case of plot.griddedVgram, these are optional arguments to plot.vgram.matrix, which in turn are optional arguments to image.plot.

Details

Here, the terms semi-variogram and variogram are used interchangeably.

This is a simple wrapper function to vgram.matrix (entire field) from fields and/or variogram.matrix (non-zero grid points only) for finding the variogram between two gridded fields. It calls this function for each of two fields in a verification set. This function allows one to do the diagnostic analysis proposed in Marzban and Sangathe (2009).

Value

A list object containing the entire list passed in by the object argument, and components:

Vx.cgram.matrix, Fcst.vgram.matrix

list objects as returned by vgram.matrix containing the variogram information for each field.

No value is returned by plot.griddedVgram, plots are created showing the empirical variogram (circles), along with directional empirical variograms (dots), and the variogram by direction (image plot).

Author(s)

Eric Gilleland

References

Marzban, C. and Sandgathe, S. (2009) Verification with variograms. Wea. Forecasting, 24 (4), 1102–1120, doi:10.1175/2009WAF2222122.1.

See Also

fields::vgram.matrix, make.SpatialVx

Examples

data( "ExampleSpatialVxSet" )

x <- ExampleSpatialVxSet$vx
xhat <- ExampleSpatialVxSet$fcst

hold <- make.SpatialVx( x, xhat, field.type = "contrived",
    units="none", data.name = "Example", obs.name = "x",
    model.name = "xhat" )

res <- griddedVgram( hold, R = 8 )
plot( res )

---