Shape analysis for spatial forecast verification (hiw is OE for shape; yields MnE hue).
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15  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)

x,object 

simplify 

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 1e6 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 blowup 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 
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 usersupplied 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 redetermined 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 nonconvex 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.).
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 twocolumn matrices containing the x and y coordinate centroids for each feature (as determined by centroid.owin). 
intensities 
list with components X and Xhat giving threecolumn 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 
twocolumn matrix whose rows indicate the feature numbers from each field being compared; corresponding to the columns of SS above. 
Eric Gilleland
Lack, S. A., Limpert, G. L. and Fox, N. I. (2010) An objectoriented 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.
To indentify features and create objects of class “features”, see, for example: FeatureFinder
centroid.owin
, as.psp
, psp
, crossing.psp
, lengths.psp
, angles.psp
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44  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")
## Not run:
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 = "ICP Geometric Cases", obs.name = "geom000",
model.name = "geom004" )
look < FeatureFinder(hold, do.smooth = FALSE)
look < hiw(look)
par(mfrow=c(1,2))
plot(look)
plot(look, which = "Xhat")
## End(Not run)

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.