SpatialVx (version 1.0-3)

FeatureAxis: 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 FeatureAxis
plot(x, ..., zoom = FALSE)
# S3 method for FeatureAxis
summary(object, ...)

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.

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.

Author

Eric Gilleland

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.

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)

if (FALSE) {
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)
   }