Shape analysis for spatial forecast verification (hiw is OE for shape; yields MnE hue).
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 hiw
distill(x, ...)
# S3 method for hiw
plot(x, ..., which = c("X", "Xhat"), ftr.num = 1, zoom = TRUE,
seg.col = "darkblue")
# S3 method for hiw
print(x, ...)
# S3 method for hiw
summary(object, ..., silent = FALSE)
A list object of class “hiw” is returned with components the same as in the original “features” class object, as well as:
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.
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).
list with components X and Xhat giving three-column matrices that contain the mean, min and max intensities for each feature.
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:
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.
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.
two-column matrix whose rows indicate the feature numbers from each field being compared; corresponding to the columns of SS above.
hiw: object of class “features”.
distill, plot, summary: object of class “hiw”.
dmin argument in call to simplify.owin from spatstat. If 0 (default), then no call is made to simplify.owin.
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.
logical, should progress information be printed to the screen?
character string naming whether to plot a feature from the obsevation field (default) or the forecast field.
integer stating which feature number to plot.
logical, should the feature be plotted within its original domain, or a blow-up of the feature (default)?
color for the line segments.
logical, should the summary information be printed to the screen?
Not used by hiw, distill, summary. Optional arguments to plot.
Eric Gilleland
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.).
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.
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 = "ICP Geometric Cases", obs.name = "geom000",
model.name = "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")
Run the code above in your browser using DataLab