SpatialVx (version 0.3)

hiw: Spatial Forecast Verification Shape Analysis

Description

Shape analysis for spatial forecast verification (hiw is OE for shape; yields MnE hue).

Usage

hiw(x, simplify = 0, A = pi * c(0, 1/16, 1/8, 1/6, 1/4, 1/2, 9/16, 5/8, 2/3, 3/4),
    verbose = FALSE, ...)
## S3 method for class 'hiw':
distill(x, ...)
## S3 method for class 'hiw':
plot(x, ..., which = c("X", "Xhat"), ftr.num = 1, zoom = TRUE, 
    seg.col = "darkblue")
## S3 method for class 'hiw':
print(x, ...)
## S3 method for class 'hiw':
summary(object, ..., silent = FALSE)

Arguments

x,object

hiw: object of class features.

distill, plot, summary: object of class hiw.

simplify

dmin argument in call to simplify.owin from spatstat. If 0 (default), then no call is made to simplify.owin.

A

numeric vector of angles for which to apply shape analysis. Note that this vector will be rounded to 6 digits. If values are less than that, might be prudent to add 1e-6 to them.

verbose

logical, should progress information be printed to the screen?

which

character string naming whether to plot a feature from the obsevation field (default) or the forecast field.

ftr.num

integer stating which feature number to plot.

zoom

logical, should the feature be plotted within its original domain, or a blow-up of the feature (default)?

seg.col

color for the line segments.

silent

logical, should the summary information be printed to the screen?

...

Not used by hiw, distill, summary. Optional arguments to plot.

Value

• A list object of class hiw is returned with components the same as in the original features class object, as well as:
• radial.segments
  a list with components X and Xhat each giving lists of the psp class (i.e., line segment) object for each feature containing the radial segments from the feature centroids to the boundaries.
• centers
  list with components X and Xhat giving two-column matrices containing the x- and y- coordinate centroids for each feature (as determined by centroid.owin).
• intensities
  list with components X and Xhat giving three-column matrices that contain the mean, min and max intensities for each feature.
• angles,lengths
  list with components X and Xhat each giving lists containing the lengths of the line segments and their respective angles. Missing values are possible here.
• distill returns an array whose dimensions are the number of landmarks (i.e., boundary points) by two by the number of observed and forecast features. An attribute called field.identifier is also given that is a character vector containing repeated X and Xhat values identiying which of the third dimension are associated with the observed field (X) and those identified with the forecast field (Xhat). Note that missing values may be present, which may need to be dealt with (by the user) before using functions from the shapes package.

  summary invisibly returns a list object with components:

• X,Xhat
  matrices whose rows represent features and whose columns give their centroids (note that x refers to the columns and y to the rows), as well as the average, min and max intensities.
• SS
  matrix with four rows and columns equal to the number of possible combinations of feature matchings between fields. Gives the sum of square translation/location error (i.e., squared centroid distance), as well as the average, min and max squared differences between each combination of features.
• ind
  two-column matrix whose rows indicate the feature numbers from each field being compared; corresponding to the columns of SS above.

Details

This function is an attempt to approximate the technique described first in Micheas et al. (2007) and as modified in Lack et al. (2010). It will only find the centroids, rays extending from them to the boundaries, and the boundary points. Use distill to convert this output into an object readable by, for example, procGPA from package shapes.

First, identified features (which may be identified by any feature identification function that yields an object of class features) are taken, the centroid is found (the centroid is found via centroid.owin so that the x- and y- coordinates are fliped from what you might expect) and very long line segments are found radiating out in both directions from the center. They are then clipped by where they cross the boundaries of the features.

The spatstat package is used heavily by this function. In particular, the function as.polygonal is applied to the owin objects (possible after first calling simplify.owin). Line segments are created using the feature centroids, as found by centroid.owin, and the user-supplied angles, along with a very long length (equal to the domain size). Boundary crossings are found using crossing.psp, and new line segment patterns are created using the centroids and boundary crossing information (extra points along line segments are subsequently removed through a painstaking process, and as.psp is called again, and any missing line segments are subsequently accounted for, for later calculations). Additionally, lengths of line segments are found via lengths.psp. Angles must also be re-determined and corresponded to the originally passed angles. Therefore, it is necessary to round the angles to 6 digits, or equal angles may not be considered equal, which will cause problems.

The hiw function merely does the above step, as well as finds the lengths of the resulting line segments. For non-convex objects, the longest line segment is returned, and if the boundary crossings do not lie on opposite sides of the centroid, then the negative of the shortest segment is returned for that particular value. Also returned are the mean, min and maximum intensities for each feature, as well as the final angles returned. It is possible to have missing values for some of these components.

The summary function computes SSloc, SSavg, SSmin, and SSmax between each pair of features between fields. distill may be used to create an object that can be further analyzed (for shape) using the shapes package.

While any feature identification function may be used, it is recommended to throw out small sized features as the results may be misleading (e.g., comparisons between features consisting of single points, etc.).

References

Lack, S. A., Limpert, G. L. and Fox, N. I. (2010) An object-oriented multiscale verification scheme. Wea. Forecasting, 25, 79--92.

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. Hydrology, 343, 105--116.

See Also

To indentify features and create objects of class features, see, for example: FeatureFinder

centroid.owin, as.psp, psp, crossing.psp, lengths.psp, angles.psp

Examples

data(geom000)
data(geom001)
data(geom004)
data(ICPg240Locs)

hold <- make.SpatialVx(geom000, geom001, thresholds=c(0.01,50.01),
    projection=TRUE, map=TRUE, loc=ICPg240Locs, loc.byrow = TRUE,
    field.type="Geometric Objects Pretending to be Precipitation",
    units="mm/h", data.name=c("ICP Geometric Cases", "geom000", "geom001"))

look <- FeatureFinder(hold, do.smooth = FALSE, thresh = 2, min.size = 200)

look <- hiw(look)

distill.hiw(look)

# Actually, you just need to type:
# distill(look)

summary(look)

# Note: procGPA will not allow missing values.

par(mfrow=c(1,2))
plot(look)
plot(look, which = "Xhat")

hold <- make.SpatialVx(geom000, geom004, thresholds=c(0.01,50.01),
    projection=TRUE, map=TRUE, loc=ICPg240Locs, loc.byrow = TRUE,
    field.type="Geometric Objects Pretending to be Precipitation",
    units="mm/h", data.name=c("ICP Geometric Cases", "geom000", "geom004"))

look <- FeatureFinder(hold, do.smooth = FALSE)

look <- hiw(look)

par(mfrow=c(1,2))
plot(look)
plot(look, which = "Xhat")