Package 'ShapeSelectForest' reference manual

Title:	Shape Selection for Landsat Time Series of Forest Dynamics
Description:	Landsat satellites collect important data about global forest conditions. Documentation about Landsat's role in forest disturbance estimation is available at the site <https://landsat.gsfc.nasa.gov/>. By constrained quadratic B-splines, this package delivers an optimal shape-restricted trajectory to a time series of Landsat imagery for the purpose of modeling annual forest disturbance dynamics to behave in an ecologically sensible manner assuming one of seven possible "shapes", namely, flat, decreasing, one-jump (decreasing, jump up, decreasing), inverted vee (increasing then decreasing), vee (decreasing then increasing), linear increasing, and double-jump (decreasing, jump up, decreasing, jump up, decreasing). The main routine selects the best shape according to the minimum Bayes information criterion (BIC) or the cone information criterion (CIC), which is defined as the log of the estimated predictive squared error. The package also provides parameters summarizing the temporal pattern including year(s) of inflection, magnitude of change, pre- and post-inflection rates of growth or recovery. In addition, it contains routines for converting a flat map of disturbance agents to time-series disturbance maps and a graphical routine displaying the fitted trajectory of Landsat imagery.
Authors:	Xiyue Liao [aut, cre] , Mary Meyer [aut], Elizabeth Freeman [aut], Gretchen Moisen [aut]
Maintainer:	Xiyue Liao <[email protected]>
License:	GPL (>= 2)
Version:	1.7
Built:	2024-12-11 06:47:13 UTC
Source:	CRAN

Shape Selection for Landsat Time Series of Forest Dynamics

Description

Given a scatterplot of $(x_i, y_i)$ , $i = 1,\ldots,n$ , where $\bold{x}$ could be a vector of years and $\bold{y}$ could be a vector of Landsat signals, constrained least-squares spline fits are obtained for the following shapes:

1. flat
2. decreasing
3. one-jump, i.e., decreasing, jump up, decreasing
4. inverted vee (increasing then decreasing)
5. vee (decreasing then increasing)
6. linear increasing
7. double-jump, i.e., decreasing, jump up, decreasing, jump up, decreasing.

The shape with the smallest information criterion may be considered a "best" fit. This shape-selection problem was motivated by a need to identify types of disturbances to areas of forest, given Landsat signals over a number of years. The satellite signal is constant or slowly decreasing for a healthy forest, with a jump upward in the signal caused by mass destruction of trees.

The main routine to select the shape for a scatterplot is "shape". See shape for more details.

Author(s)

Mary C. Meyer, Xiyue Liao, Elizabeth Freeman, Gretchen G. Moisen

Maintainer: Xiyue Liao <[email protected]>

References

Meyer, M. C. and Woodroofe M (2000) On the Degrees of Freedom in Shape-Restricted Regression. The Annals of Statistics 28, 1083–1104.

Meyer, M. C. (2013a) Semi-parametric additive constrained regression. Journal of Nonparametric Statistics 25(3), 715.

Meyer, M. C. (2013b) A simple new algorithm for quadratic programming with applications in statistics. Communications in Statistics 42(5), 1126–1139.

Liao, X. and M. C. Meyer (2014) coneproj: An R package for the primal or dual cone projections with routines for constrained regression. Journal of Statistical Software 61(12), 1–22.

A 21 by 7 Matrix Storing Edf0 Vectors

Description

The object "edf0s" is a $21$ by $7$ matrix. Each row is an edf0 vector of $7$ elements corresponding to the $7$ shapes in this package. Such a vector will be used in the main routine "shape" to select the best shape for a scatterplot of Landsat signals. Each edf0 vector is simulated through a subroutine called "getedf0", using a total of $1000$ simulations with the random seed being set to $123$ . Each row is an edf0 vector for an equally spaced $\bold{x}$ vector of $n$ elements (e.g., years). From the first row to the last row, the edf0 vector is for a predictor vector $\bold{x}$ of length $n$ which is an integer ranging from $20$ to $40$ . The matrix is built for the convenience of users when they call the routine "shape".

If the $\bold{x}$ vector is equally spaced and its number of elements $n$ is between $20$ and $40$ , then a corresponding edf0 vector will be extracted directly from this matrix and no simulation will be done, which saves a lot of time; otherwise, the subroutine "getedf0" will be called inside the routine "shape" to get an edf0 vector for $\bold{x}$ . The timing depends on the number of elements in $\bold{x}$ and the shape options allowed by the user. For example, when $\bold{x}$ is an equally spaced vector of $26$ elements, the timing is about $167$ seconds if the user allows a double-jump shape and the timing is about $12$ seconds if the user doesn't allow a double-jump shape. Also, when $\bold{x}$ is not an equally spaced vector, no matter how many elements it has, "getedf0" should be called.

Usage

data("edf0s")data("edf0s")

Format

A 21 by 7 matrix.

f2a.map.jpeg

Description

Creates a time series of jpegs. One jpeg is created for each year of the map output of f2a.raster.

Usage

f2a.map.jpeg(years, folder, OUTPUT.fn, height = 10, 
width = 10 * (dim(mapgrid.dist)[2] / dim(mapgrid.dist)[1]), 
units = "in", res = 400)
f2a.map.jpeg(years, folder, OUTPUT.fn, height = 10, 
width = 10 * (dim(mapgrid.dist)[2] / dim(mapgrid.dist)[1]), 
units = "in", res = 400)

Arguments

`years`	Vector of the years included in the time series data.
`folder`	Folder (full path) containing the input rasters. This is also the folder where the output will be written.
`OUTPUT.fn`	Filename of output from `f2a.raster`.
`height`	The height of the device.
`width`	The width of the device. The default is the width that will give accurate ratio of height to width for a raster.
`units`	The units in which `height` and `width` are given, which can be `px` (pixels, the default), `in` (inches), `cm` or `mm`.
`res`	The nominal resolution in ppi which will be recorded in the bitmap file, if a positive integer. Also used for units other than the default. If not specified, taken as 300 ppi to set the size of text and line widths.

Details

Creates one jpeg for each year in years.

Value

Returns nothing

Author(s)

Liz Freeman

Examples

## Not run: 
	# define years
	years <- c(1984, 1986:2010)
	
	# define a folder for all outputs
	folder.in <- paste(system.file(package = "ShapeSelectForest"), 
	"extdata", "helpexamples", sep = "/")
	folder.out <- getwd()

	# define filenames
	flat.pred.fn <- "MINI_FLATPRED.img"
	b5.fn <- "MINI_B5.img"
	fi.fn <- "MINI_FI.img"
	nbr.fn <- "MINI_NBR.img"
	ndvi.fn <- "MINI_NDVI.img"
	INPUT.bands <- c(b5.fn, fi.fn, nbr.fn, ndvi.fn)

	# call f2a.raster
	ans1 <- f2a.raster(years = years, folder.in = folder.in, folder.out = folder.out, 
	OUTPUT.fn <- "f2a_example.img", flat.pred.fn = flat.pred.fn, INPUT.bands = INPUT.bands)

	# create jpegs 
	ans2 <- f2a.map.jpeg(years = years, folder = folder.out, OUTPUT.fn = "f2a_example.img")

## End(Not run)
## Not run: 
	# define years
	years <- c(1984, 1986:2010)
	
	# define a folder for all outputs
	folder.in <- paste(system.file(package = "ShapeSelectForest"), 
	"extdata", "helpexamples", sep = "/")
	folder.out <- getwd()

	# define filenames
	flat.pred.fn <- "MINI_FLATPRED.img"
	b5.fn <- "MINI_B5.img"
	fi.fn <- "MINI_FI.img"
	nbr.fn <- "MINI_NBR.img"
	ndvi.fn <- "MINI_NDVI.img"
	INPUT.bands <- c(b5.fn, fi.fn, nbr.fn, ndvi.fn)

	# call f2a.raster
	ans1 <- f2a.raster(years = years, folder.in = folder.in, folder.out = folder.out, 
	OUTPUT.fn <- "f2a_example.img", flat.pred.fn = flat.pred.fn, INPUT.bands = INPUT.bands)

	# create jpegs 
	ans2 <- f2a.map.jpeg(years = years, folder = folder.out, OUTPUT.fn = "f2a_example.img")

## End(Not run)

Raster based flat-to-annual function.

Description

Applies flat2annual to each pixel in a raster to produce time series maps of disturbance.

Usage

f2a.raster(years, folder.in, folder.out, OUTPUT.fn, flat.pred.fn, INPUT.bands,
layer.shape = 1, layer.dyr = 2, layer.dur = 5)
f2a.raster(years, folder.in, folder.out, OUTPUT.fn, flat.pred.fn, INPUT.bands,
layer.shape = 1, layer.dyr = 2, layer.dur = 5)

Arguments

years

Vector of the years included in the time series data.

folder.in

Folder (full path) containing the input rasters.

folder.out

Folder where the output will be written.

OUTPUT.fn

Filename for output. The extension of this filename will specify the output file type. For image files, OUTPUT.fn should end in .img.

flat.pred.fn

Filename of a single layer .img file of the disturbance agents. The agents are specified by the integers from 0 to 6.

			`0`	Unclassified
			`1`	Conversion
			`2`	Fire
			`3`	Harvest
			`4`	Stable
			`5`	Stress
			`6`	Recovery

INPUT.bands

Filenames of the multi-layer .img files of the shapes to be used.

layer.shape

Number giving the layer of the shape files containing the shape data.

layer.dyr

Number giving the layer of the shape files containing the disturbance year data.

layer.dur

Number giving the layer of the shape files containing the disturbance duration data.

Details

The function writes a multi-layer raster with one layer for each year in years given the predicted agent for each pixel at each year.

The layers for shape, dyr and dur need to be the same in all files named in INPUT.bands. The default is layer.shape = 1, layer.dyr = 2, and layer.dur = 5.

Value

The function does not return a value. Instead, a multi-band .img map file is created.

Author(s)

Liz Freeman

Examples

## Not run: 
	# define years
	years <- c(1984, 1986:2010)

	# define a folder for all output
	folder.in <- paste(system.file(package = "ShapeSelectForest"), 
	"extdata", "helpexamples", sep = "/")
	folder.out <- getwd()	

	# define filenames
	flat.pred.fn <- "MINI_FLATPRED.img"
	b5.fn <- "MINI_B5.img"
	fi.fn <- "MINI_FI.img"
	nbr.fn <- "MINI_NBR.img"
	ndvi.fn <- "MINI_NDVI.img"
	INPUT.bands <- c(b5.fn, fi.fn, nbr.fn, ndvi.fn)

	# call f2a.raster
	ans <- f2a.raster(years = years, folder.in = folder.in, folder.out = folder.out, 
	OUTPUT.fn = "f2a_example.img", flat.pred.fn = flat.pred.fn, INPUT.bands = INPUT.bands)

## End(Not run)
## Not run: 
	# define years
	years <- c(1984, 1986:2010)

	# define a folder for all output
	folder.in <- paste(system.file(package = "ShapeSelectForest"), 
	"extdata", "helpexamples", sep = "/")
	folder.out <- getwd()	

	# define filenames
	flat.pred.fn <- "MINI_FLATPRED.img"
	b5.fn <- "MINI_B5.img"
	fi.fn <- "MINI_FI.img"
	nbr.fn <- "MINI_NBR.img"
	ndvi.fn <- "MINI_NDVI.img"
	INPUT.bands <- c(b5.fn, fi.fn, nbr.fn, ndvi.fn)

	# call f2a.raster
	ans <- f2a.raster(years = years, folder.in = folder.in, folder.out = folder.out, 
	OUTPUT.fn = "f2a_example.img", flat.pred.fn = flat.pred.fn, INPUT.bands = INPUT.bands)

## End(Not run)

Raster based flat-to-parameter function.

Description

Applies flat2parameter to each pixel in a raster to produce maps of disturbance parameters.

Usage

f2p.raster(years, folder.in, folder.out, OUTPUT.fn, flat.pred.fn, INPUT.bands,
layer.shape = 1, layer.dyr = 2, layer.dur = 5, layer.mag = 3)
f2p.raster(years, folder.in, folder.out, OUTPUT.fn, flat.pred.fn, INPUT.bands,
layer.shape = 1, layer.dyr = 2, layer.dur = 5, layer.mag = 3)

Arguments

years

Vector of the years included in the time series data.

folder.in

Folder (full path) containing the input rasters.

folder.out

Folder where the output will be written.

OUTPUT.fn

Filename for output. The extension of this filename will specify the output file type. For image files, OUTPUT.fn should end in .img

flat.pred.fn

Filename of a single layer .img file of the disturbance agents. The agents are specified by the integers from 0 to 6.

			`0`	Unclassified
			`1`	Conversion
			`2`	Fire
			`3`	Harvest
			`4`	Stable
			`5`	Stress
			`6`	Recovery

INPUT.bands

Filenames of the multi-layer .img files of the shapes to be used.

layer.shape

Number giving the layer of the shape files containing the shape data.

layer.dyr

Number giving the layer of the shape files containing the disturbance year data.

layer.dur

Number giving the layer of the shape files containing the disturbance duration data.

layer.mag

Number giving the layer of the shape files containing the disturbance magnitude data.

Details

The function writes a seven-layer raster with layers for disturbance agent, median disturbance year, median disturbance duration and magnitude of all image files given in INPUT.bands, in the order the filenames are listed in INPUT.bands.

The layers for shape, dyr and dur need to be the same in all files named in INPUT.bands. The default is layer.shape = 1, layer.dyr = 2, layer.dur = 5, and layer.mag = 3.

If, for example, INPUT.bands = c(b5.fn, fi.fn, nbr.fn, and ndvi.fn), then the layers of the output file are:

			`1`	Disturbance Agent
			`2`	Disturbance Year
			`3`	Disturbance Duration
			`4`	Magnitude B5
			`5`	Magnitude FI
			`6`	Magnitude NBR
			`7`	Magnitude NDVI

Value

The function does not return a value. Instead, a multi-band .img map file is created.

Author(s)

Liz Freeman

Examples

## Not run: 
	# define years
	years <- c(1984, 1986:2010)

	# define a folder for all output
	folder.in <- paste(system.file(package = "ShapeSelectForest"), 
	"extdata", "helpexamples", sep = "/")
	folder.out <- getwd()	

	# define filenames
	flat.pred.fn <- "MINI_FLATPRED.img"
	b5.fn <- "MINI_B5.img"
	fi.fn <- "MINI_FI.img"
	nbr.fn <- "MINI_NBR.img"
	ndvi.fn <- "MINI_NDVI.img"
	INPUT.bands <- c(b5.fn, fi.fn, nbr.fn, ndvi.fn)

	# call f2p.raster
	ans <- f2p.raster(years = years, folder.in = folder.in, folder.out = folder.out, 
	OUTPUT.fn = "f2p_example.img", flat.pred.fn = flat.pred.fn, INPUT.bands = INPUT.bands)

## End(Not run)
## Not run: 
	# define years
	years <- c(1984, 1986:2010)

	# define a folder for all output
	folder.in <- paste(system.file(package = "ShapeSelectForest"), 
	"extdata", "helpexamples", sep = "/")
	folder.out <- getwd()	

	# define filenames
	flat.pred.fn <- "MINI_FLATPRED.img"
	b5.fn <- "MINI_B5.img"
	fi.fn <- "MINI_FI.img"
	nbr.fn <- "MINI_NBR.img"
	ndvi.fn <- "MINI_NDVI.img"
	INPUT.bands <- c(b5.fn, fi.fn, nbr.fn, ndvi.fn)

	# call f2p.raster
	ans <- f2p.raster(years = years, folder.in = folder.in, folder.out = folder.out, 
	OUTPUT.fn = "f2p_example.img", flat.pred.fn = flat.pred.fn, INPUT.bands = INPUT.bands)

## End(Not run)

flat2annual

Description

Identifies the disturbance year of a single pixel or a plot location.

Usage

flat2annual(years, all.shapes, all.durs, all.dyrs, mtbs, flat.pred)
flat2annual(years, all.shapes, all.durs, all.dyrs, mtbs, flat.pred)

Arguments

`years`	Vector of the years included in the time series data.
`all.shapes`	Vector (length 4) of the shapes of the four remote sensing bands. Shape values range from 1 to 7.
`all.durs`	Vector (length 4) of the duration of each shape.
`all.dyrs`	Vector (length 4) of the year of each shape.
`mtbs`	MTBS
`flat.pred`	Predicted disturbance Agent. Agent values range from 0 to 6.

Details

flat2annual can be used on either a single pixel or on a single data point in a data frame.

Value

Returns a vector of the same length as years with a predicted agent for each year in years.

Author(s)

Liz Freeman

Examples

	# define years
	years <- c(2001:2010)

	# define parameters
	all.shapes <- c(1, 4, 5, 3)
	all.dyrs <- c(2001, 0, 2004, 2004)
	all.durs <- c(1, 0, 3, 5)
	flat.pred <- 5

	# call flat2annual
	ans <- flat2annual(years = years, all.shapes = all.shapes, all.durs = all.durs,
	all.dyrs = all.dyrs, mtbs = mtbs, flat.pred = flat.pred)
# define years
	years <- c(2001:2010)

	# define parameters
	all.shapes <- c(1, 4, 5, 3)
	all.dyrs <- c(2001, 0, 2004, 2004)
	all.durs <- c(1, 0, 3, 5)
	flat.pred <- 5

	# call flat2annual
	ans <- flat2annual(years = years, all.shapes = all.shapes, all.durs = all.durs,
	all.dyrs = all.dyrs, mtbs = mtbs, flat.pred = flat.pred)

flat2parameter

Description

Based on the disturbance type and the shapes of the four remote sensing bands, a vector of averaged parameters is produced. Which bands are included in the average depends on both the predicted disturbance type and the shapes of the four bands.

Usage

flat2parameter(years, all.shapes, all.durs, all.dyrs, all.mags, mtbs, flat.pred)
flat2parameter(years, all.shapes, all.durs, all.dyrs, all.mags, mtbs, flat.pred)

Arguments

`years`	Vector of the years included in the time series data.
`all.shapes`	Vector (length 4) of the shapes of the four remote sensing bands. Shape values range from 1 to 7.
`all.durs`	Vector (length 4) of the duration of each shape.
`all.dyrs`	Vector (length 4) of the year of each shape.
`all.mags`	Vector (length 4) of the magnitude of each shape.
`mtbs`	MTBS.
`flat.pred`	Predicted disturbance Agent. Agent values range from 0 to 6.

Details

flat2parameter can be used on either a single pixel or on a single data point in a data frame.

Value

Returns a vector of length 4 containing the average disturbance year, the duration, the magnitude, and the predicted type.

Author(s)

Liz Freeman

Examples

	# define years
	years <- c(2001:2010)

	# define parameters
	all.shapes <- c(1, 4, 5, 3)
	all.dyrs <- c(2001, 0, 2004, 2004)
	all.durs <- c(1, 0, 3, 5)
	all.mags <- c(100, 0, 1000, 1500)
	flat.pred <- 5

	# call flat2parameter
	ans <- flat2parameter(years = years, all.shapes = all.shapes, all.durs = all.durs,
	all.dyrs = all.dyrs, all.mags = all.mags, mtbs = mtbs, flat.pred = flat.pred)
# define years
	years <- c(2001:2010)

	# define parameters
	all.shapes <- c(1, 4, 5, 3)
	all.dyrs <- c(2001, 0, 2004, 2004)
	all.durs <- c(1, 0, 3, 5)
	all.mags <- c(100, 0, 1000, 1500)
	flat.pred <- 5

	# call flat2parameter
	ans <- flat2parameter(years = years, all.shapes = all.shapes, all.durs = all.durs,
	all.dyrs = all.dyrs, all.mags = all.mags, mtbs = mtbs, flat.pred = flat.pred)

Get the Edf0 Vector for the Shape Routine

Description

An edf0 vector is the estimated "null expected degrees of freedom" for shapes allowed by the user. It is an input of the main routine "shape" and it is used to select the best shape for a scatterplot. See Meyer (2013a) and Meyer (2013b) for further details.

Usage

getedf0(x, flat = TRUE, dec = TRUE, jp = TRUE, invee = TRUE, 
vee = TRUE, inc = TRUE, db = TRUE, nsim = 1e+3, random = FALSE, msg = FALSE)getedf0(x, flat = TRUE, dec = TRUE, jp = TRUE, invee = TRUE, 
vee = TRUE, inc = TRUE, db = TRUE, nsim = 1e+3, random = FALSE, msg = FALSE)

Arguments

`x`	A $n$ by $1$ predictor vector, for example, years.
`flat`	A logical flag. If it is TRUE, there is a flat shape choice; otherwise, there is no such a shape option.
`dec`	A logical flag. If it is TRUE, there is a decreasing shape choice; otherwise, there is no such a shape option.
`jp`	A logical flag. If it is TRUE, there is a one-jump shape choice; otherwise, there is no such a shape option.
`invee`	A logical flag. If it is TRUE, there is an inverted-vee shape choice; otherwise, there is no such a shape option.
`vee`	A logical flag. If it is TRUE, there is a vee shape choice; otherwise, there is no such a shape option.
`inc`	A logical flag. If it is TRUE, there is an increasing shape choice; otherwise, there is no such a shape option.
`db`	A logical flag. If it is TRUE, there is a double-jump option; otherwise, there is no such a shape option.
`nsim`	Number of simulations used to get the edf0 vector. The default is nsim = 1e+3.
`random`	A parameter used by the maintainer to test if each shape option can be both included and excluded.
`msg`	A logical flag. If msg is TRUE, then a warning message will be printed when there is a non-convergence problem; otherwise no warning message will be printed. The default is msg = FALSE

Details

Because the calculations for the edf0 vector for a given set of $\bold{x}$ values (e.g., years) is time-consuming, this is accomplished in the subroutine "getedf0", and the edf0 vector is an input to the main routine "shape". In this way the edf0 values can be determined for one set of years and used for many scatterplots.

Value

The edf0 values for all shape options allowed by the user.

Author(s)

Mary C. Meyer and Xiyue Liao

References

Meyer, M. C. (2013a) Semi-parametric additive constrained regression. Journal of Nonparametric Statistics 25(3), 715

Meyer, M. C. (2013b) A simple new algorithm for quadratic programming with applications in statistics. Communications in Statistics 42(5), 1126–1139.

Examples

## Not run: 
	# define the predictor vector: the year 1985 to the year 2010
	x <- 1985:2010

	# call the getedf0 routine without a double-jump option
	edf0 <- getedf0(x, db = FALSE)

## End(Not run)
## Not run: 
	# define the predictor vector: the year 1985 to the year 2010
	x <- 1985:2010

	# call the getedf0 routine without a double-jump option
	edf0 <- getedf0(x, db = FALSE)

## End(Not run)

The Plot Routine for an Object of the Shape Routine

Description

This routine can plot the best shape selected by the shape routine for each scatterplot of Landsat signals. It can also plot the "BIC" or "CIC" values against shapes for each scatterplot, which is a way to verify that the best shape selected has the smallest "BIC" or "CIC" value.

Usage

plotshape(object, ids = 1, color = "mediumorchid4", lty = 1, lwd = 1, 
cex = .83, cex.main = .93, form = TRUE, icpic = FALSE, both = TRUE, 
tt = NULL, transpose = FALSE, plot = graphics::plot)plotshape(object, ids = 1, color = "mediumorchid4", lty = 1, lwd = 1, 
cex = .83, cex.main = .93, form = TRUE, icpic = FALSE, both = TRUE, 
tt = NULL, transpose = FALSE, plot = graphics::plot)

Arguments

`object`	An object of the shape routine.
`ids`	An integer vector representing a subset of the columns of ymat in the shape routine. Each column of ymat is a time series of Landsat imagery. Suppose that the dimension of ymat is $n$ by $k$ , then each element of "ids" should be an integer between $1$ and $k$ . The default is ids = 1.
`color`	The col argument inherited from the plot routine. The default is color = "mediumorchid4".
`lty`	The lty argument inherited from the lines routine. The default is lty = 1.
`lwd`	The lwd argument inherited from the lines routine. The default is lwd = 1.
`cex`	The cex argument inherited from the plot routine. The default is cex = .83.
`cex.main`	The cex.main argument inherited from the par routine. The default is cex.main = .93.
`form`	A logical flag. If it is TRUE, the user will let the plotshape routine to decide the layout of pictures corresponding to the elements in the "ids" vector; otherwise, the user needs to define the layout of pictures before they call the plotshape routine using par(mfrow = ...) or par(mfcol = ...). The default is form = TRUE.
`icpic`	A logical flag. Given an "ids" vector, if it is TRUE, "BIC" or "CIC" values will be plotted against all shapes allowed by the user and the fitted trajectory will not be plotted; otherwise, the fitted trajectory will also be plotted. The default is icpic = FALSE.
`both`	A logical flag. If it is TRUE, then for each element of the "ids" vector, both the fitted trajectory and the "BIC" or "CIC" plot will be made; otherwise, the "BIC" or "CIC" values will not be plotted. The default is both = TRUE.
`tt`	A vector of titles. The user can define its element as a name of the location for a scatterplot, like "Pixel 420". The default is tt = NULL.
`transpose`	A logical flag which can be used only when form is TRUE. If form is TRUE, then the user can transpose the layout the plotshape routine uses to arrange the pictures. For example, the user wants to make 2 pictures, and the default layout of this routine is 2 rows and 1 column. By setting transpose equal to TRUE, the layout will be 1 row and 2 columns. The default is transpose = TRUE.
`plot`	The genetic plot routine in graphics.

Value

A plot showing the fitted trajectory of the best shape, or showing the "BIC" or "CIC" values against shapes of an object of the shape routine.

Author(s)

Xiyue Liao

Examples

## Not run: 
	# import the matrix of Landsat signals 
	data("ymat")	

	# define the predictor vector: the year 1985 to the year 2010
	x <- 1985:2010

	# make a fit by the shape routine using "CIC" 
	# and not allow a double jump shape.
	ans <- shape(x, ymat, "CIC", db = FALSE)

	# make a plot for the 1st column of ymat 
	plotshape(ans, ids = 1, both = TRUE, form = TRUE, tt = "Pixel 420")

	# transpose the layout
	plotshape(ans, ids = 1, both = TRUE, form = TRUE, tt = "Pixel 420", transpose = TRUE)

	# make a plot for each of the first 6 columns of ymat 
	# showing the best shape 
	# and "CIC" values against the 7 shapes for each plot.
	par(mfrow = c(3, 2))
	plotshape(ans, ids = 1:6)

	# make a plot for each of the first 6 columns of ymat
	# showing both the best shape
	# and "CIC" values against the 7 shapes for each plot. 		
	# Let the routine make the layout.  
	plotshape(ans, ids = 1:6, form = TRUE, col = 2)

	# plot the ic values only
	plotshape(ans, ids = 1:6, form = TRUE, col = 5, icpic = TRUE)

	# make a title vector
	tts <- paste('Pixel', 1:36, sep = " ")

	# make all plots for the 36 scatterplots with the title vector 
	plotshape(ans, ids = 1:15, both = TRUE, form = TRUE, tt = tts[1:15], cex = .5)
	plotshape(ans, ids = 16:30, both = TRUE, form = TRUE, tt = tts[16:30], lty = 2, cex = .3)
	plotshape(ans, ids = 31:36, both = TRUE, form = TRUE, tt = tts[31:36], lty = 2, cex = .1)

## End(Not run)
## Not run: 
	# import the matrix of Landsat signals 
	data("ymat")	

	# define the predictor vector: the year 1985 to the year 2010
	x <- 1985:2010

	# make a fit by the shape routine using "CIC" 
	# and not allow a double jump shape.
	ans <- shape(x, ymat, "CIC", db = FALSE)

	# make a plot for the 1st column of ymat 
	plotshape(ans, ids = 1, both = TRUE, form = TRUE, tt = "Pixel 420")

	# transpose the layout
	plotshape(ans, ids = 1, both = TRUE, form = TRUE, tt = "Pixel 420", transpose = TRUE)

	# make a plot for each of the first 6 columns of ymat 
	# showing the best shape 
	# and "CIC" values against the 7 shapes for each plot.
	par(mfrow = c(3, 2))
	plotshape(ans, ids = 1:6)

	# make a plot for each of the first 6 columns of ymat
	# showing both the best shape
	# and "CIC" values against the 7 shapes for each plot. 		
	# Let the routine make the layout.  
	plotshape(ans, ids = 1:6, form = TRUE, col = 2)

	# plot the ic values only
	plotshape(ans, ids = 1:6, form = TRUE, col = 5, icpic = TRUE)

	# make a title vector
	tts <- paste('Pixel', 1:36, sep = " ")

	# make all plots for the 36 scatterplots with the title vector 
	plotshape(ans, ids = 1:15, both = TRUE, form = TRUE, tt = tts[1:15], cex = .5)
	plotshape(ans, ids = 16:30, both = TRUE, form = TRUE, tt = tts[16:30], lty = 2, cex = .3)
	plotshape(ans, ids = 31:36, both = TRUE, form = TRUE, tt = tts[31:36], lty = 2, cex = .1)

## End(Not run)

Shape Selection

Description

Given a predictor vector $\bold{x}$ , e.g., years, and a matrix $\bold{ymat}$ whose columns are response vectors, e.g., Landsat signals. The shape routine will select a shape that is the best fit for each response vector according to the Bayes information criterion (BIC) or the cone information criterion (CIC).

Usage

shape(x, ymat, infocrit = "CIC", flat = TRUE, dec = TRUE, jp = TRUE, 
invee = TRUE, vee = TRUE, inc = TRUE, db = TRUE, nsim = 1e+3, 
edf0 = NULL, get.edf0 = FALSE, random = FALSE, msg = FALSE)shape(x, ymat, infocrit = "CIC", flat = TRUE, dec = TRUE, jp = TRUE, 
invee = TRUE, vee = TRUE, inc = TRUE, db = TRUE, nsim = 1e+3, 
edf0 = NULL, get.edf0 = FALSE, random = FALSE, msg = FALSE)

Arguments

`x`	A $n$ by $1$ predictor vector, for example, years.
`ymat`	A $n$ by $N$ matrix whose columns are response vectors corresponding to x, for example, Landsat signals.
`infocrit`	The criterion used to select the best shape for a scatterplot. It can either be the Bayes information criterion (BIC) or the cone information criterion (CIC).
`flat`	A logical flag. If it is TRUE, there is a flat shape choice; otherwise, there is no such a shape option.
`dec`	A logical flag. If it is TRUE, there is a decreasing shape choice; otherwise, there is no such a shape option.
`jp`	A logical flag. If it is TRUE, there is a one-jump shape choice; otherwise, there is no such a shape option.
`invee`	A logical flag. If it is TRUE, there is an inverted-vee shape choice; otherwise, there is no such a shape option.
`vee`	A logical flag. If it is TRUE, there is a vee shape choice; otherwise, there is no such a shape option.
`inc`	A logical flag. If it is TRUE, there is an increasing shape choice; otherwise, there is no such a shape option.
`db`	A logical flag. If it is TRUE, there is a double-jump shape choice; otherwise, there is no such a shape option. The routine is usually slower when there is a double-jump shape choice than it is when there is no such a choice.
`nsim`	Number of simulations used to get the edf0 vector. The default is nsim = 1e+3. See references in this section for more details about edf0.
`edf0`	The edf0 given by the user. When $\bold{x}$ is an equally spaced vector whose number of elements is between $20$ and $40$ . The user doesn't need to provide an edf0 vector; otherwise, the user has to set get.edf0 to be TRUE such that the shape routine will simulate an edf0 vector, or the user can choose to simulate an edf0 vector by the getedf0 routine and provide the edf0 vector to the shape routine with this argument. The default is edf0 = NULL.
`get.edf0`	A logical flag. When $\bold{x}$ is not an equally spaced vector whose number of elements is between $20$ and $40$ . The user has to set get.edf0 to be TRUE such that the shape routine will simulate an edf0 vector, or the user can choose to simulate an edf0 vector by the “getedf0” routine and provide the edf0 vector to the shape routine with the edf0 argument. The default is get.edf0 = FALSE.
`random`	A parameter used by the maintainer to test if each shape option can be both included and excluded.
`msg`	A logical flag. If msg is TRUE, then a warning message will be printed when there is a non-convergence problem; otherwise no warning message will be printed. The default is msg = FALSE

Details

Given a scatterplot of $(x_i, y_i)$ , $i=1,\ldots,n$ , where $\bold{x}$ could be a vector of years and $\bold{y}$ could be a vector of Landsat signals, constrained least-squares spline fits are obtained for the following shapes:

1. flat
2. decreasing
3. one-jump, i.e., decreasing, jump up, decreasing
4. inverted vee (increasing then decreasing)
5. vee (decreasing then increasing)
6. linear increasing
7. double-jump, i.e., decreasing, jump up, decreasing, jump up, decreasing.

The "shape" routine chooses one of the shapes allowed by the user based on the minimum Bayes information criterion (BIC) or the cone information criterion (CIC). It also returns the information criterion (IC) values for shapes allowed by the user. Fitting method is constrained quadratic B-splines, number of knots depends on number of observations. The cone projection algorithm used in this routine is implemented by the R package coneproj.

See references cited in this section and the official manual (https://cran.r-project.org/package=coneproj) for the R package coneproj for more details.

Value

`shape`	A $N$ by $1$ vector. The $i$ th element is the best shape for each of the $i$ th scatterplot.
`ic`	A $k$ by $N$ matrix where the $i$ th column is the vector of "BIC" or "CIC" values used to choose the best shape for the $i$ th scatterplot. $k$ is the number of shapes allowed by the user.
`thetab`	A $n$ by $N$ matrix where the $i$ th column is the vector of predicted values for the chosen shape for the $i$ th scatterplot.
`x`	The argument x.
`ymat`	The argument ymat.
`infocrit`	The argument infocrit.
`k`	The number of knots used.
`bs`	A list of coefficient vectors. Each vector is the vector of coefficients for regression basis functions for each scatterplot.
`ijps`	A list storing the position of the first jump for scatterplots whose best shape is one-jump or double-jump. It also stores the position of the knot from where $\bold{f}$ starts increasing (decreasing) for scatterplots whose best shape is vee (inverted vee).
`jjps`	A list storing the position of the second jump for scatterplots whose best shape is double-jump.
`m_is`	A vector storing the centering values for the first ramp edge for scatterplots whose best shape is one-jump or double-jump.
`m_js`	A vector storing the centering values for the second ramp edge for scatterplots whose best shape is double-jump.
`tm`	Total cpu running time.

Author(s)

Mary C. Meyer and Xiyue Liao

References

Meyer, M. C. (2013a) Semi-parametric additive constrained regression. Journal of Nonparametric Statistics 25(3), 715.

Meyer, M. C. (2013b) A simple new algorithm for quadratic programming with applications in statistics. Communications in Statistics 42(5), 1126–1139.

Liao, X. and M. C. Meyer (2014) coneproj: An R package for the primal or dual cone projections with routines for constrained regression. Journal of Statistical Software 61(12), 1–22.

Examples

	# import the matrix of Landsat signals 
	data("ymat")

	# define the predictor vector: the year 1985 to the year 2010	
	x <- 1985:2010
## Not run: 
	# Example 1:	
	# call the shape routine allowing a double jump shape using "BIC"
	ans <- shape(x, ymat, "BIC")
	plotshape(ans, ids = 1:6, both = TRUE, form = TRUE)

## End(Not run)
## Not run: 
	# Example 2:
	# call the shape routine not allowing a double jump shape using "CIC"
	ans <- shape(x, ymat, "CIC", db = FALSE)
	plotshape(ans, ids = 1:6, both = TRUE, form = TRUE)

## End(Not run)
# import the matrix of Landsat signals 
	data("ymat")

	# define the predictor vector: the year 1985 to the year 2010	
	x <- 1985:2010
## Not run: 
	# Example 1:	
	# call the shape routine allowing a double jump shape using "BIC"
	ans <- shape(x, ymat, "BIC")
	plotshape(ans, ids = 1:6, both = TRUE, form = TRUE)

## End(Not run)
## Not run: 
	# Example 2:
	# call the shape routine not allowing a double jump shape using "CIC"
	ans <- shape(x, ymat, "CIC", db = FALSE)
	plotshape(ans, ids = 1:6, both = TRUE, form = TRUE)

## End(Not run)

Shape Parameters

Description

Given the output from the shape function (including the chosen shape, chosen information criteria value ic, vector of fitted values thetab, and corresponding $\bold{x}$ , e.g., years), this routine calculates a set of parameters that describe the behavior of the fitted trajectory.

Usage

shapeparams(shapenum, ic, thetab, x)shapeparams(shapenum, ic, thetab, x)

Arguments

`shapenum`	A number with the index $1$ to $7$ .
`ic`	A $k$ by $N$ matrix where the $i$ th column is the vector of "BIC" or "CIC" values used to choose the best shape for the $i$ th scatterplot. $k$ is the number of shapes allowed by the user.
`thetab`	A $n$ by $N$ matrix where the $i$ th column is the vector of predicted values for the chosen shape for the $i$ th scatterplot.
`x`	A $n$ by $1$ predictor vector, e.g., years.

Value

`shapenum`	the shapenum argument
`pre.rate`	annual rate of decline prior to the primary change point
`pre.rate2`	annual rate of decline prior to the secondary change point
`dist.yr`	year of the primary change points
`dist2.yr`	year of the secondary change points
`dist.mag`	difference in predicted values before and after primary change events
`dist2.mag`	difference in predicted values before and after secondary change events
`dist.mag2`	difference in predicted values before and after primary change points scaled by starting value
`dist2.mag2`	difference in predicted values before and after secondary change points scaled by starting value
`dist.dur`	duration of the change event before resuming a downward turn
`dist2.dur`	duration of the change event before resuming a downward turn
`post.rate`	annual rate of decline after the end of the primary change event
`post2.rate`	annual rate of decline after the end of the secondary change event
`my.ic`	information criteria value for the chosen shape

Author(s)

Gretchen G. Moisen

References

Moisen, G.G., M. Meyer, T.A. Schroeder, C. Toney, X. Liao, E.A. Freeman, K. Schleeweis. Shape-selection in Landsat time series: A tool for monitoring forest dynamics (In Review). Global Change Biology.

Examples

## Not run: 
	# import the matrix of Landsat signals 
	data("ymat")
	
	# define the predictor vector: the year 1985 to the year 2010
	x <- 1985:2010
 	
	# call the shape routine allowing a double-jump shape using "CIC"
	ans <- shape(x, ymat, "CIC")

	# Example 1: parameters for a flat shape
	flat_id <- which(ans$shape == 1)
	i <- flat_id[1]
	ans_flat <- shapeparams(ans$shape[i], ans$ic[, i], ans$thetab[, i], x)	

	# Example 2: parameters for a one-jump shape
	jp_id <- which(ans$shape == 3)
	i <- jp_id[1]
	ans_jp <- shapeparams(ans$shape[i], ans$ic[, i], ans$thetab[, i], x)	

	# Example 3: parameters for a double-jump shape
	db_id <- which(ans$shape == 7)
	i <- db_id[1]
	ans_db <- shapeparams(ans$shape[i], ans$ic[, i], ans$thetab[, i], x)		

## End(Not run)
## Not run: 
	# import the matrix of Landsat signals 
	data("ymat")
	
	# define the predictor vector: the year 1985 to the year 2010
	x <- 1985:2010
 	
	# call the shape routine allowing a double-jump shape using "CIC"
	ans <- shape(x, ymat, "CIC")

	# Example 1: parameters for a flat shape
	flat_id <- which(ans$shape == 1)
	i <- flat_id[1]
	ans_flat <- shapeparams(ans$shape[i], ans$ic[, i], ans$thetab[, i], x)	

	# Example 2: parameters for a one-jump shape
	jp_id <- which(ans$shape == 3)
	i <- jp_id[1]
	ans_jp <- shapeparams(ans$shape[i], ans$ic[, i], ans$thetab[, i], x)	

	# Example 3: parameters for a double-jump shape
	db_id <- which(ans$shape == 7)
	i <- db_id[1]
	ans_db <- shapeparams(ans$shape[i], ans$ic[, i], ans$thetab[, i], x)		

## End(Not run)

Response Variable Matrix

Description

This is a $26$ by $36$ matrix. Each column is a trajectory of Landsat signals corresponding to $26$ consecutive years ranging from $1985$ to $2010$ . It will be used in some examples of this package. They are trajectories from $36$ pixels in South Carolina.

Usage

data("ymat")data("ymat")

Format

A $26$ by $36$ matrix.

Source

US Forest Service

Package 'ShapeSelectForest'

Help Index

Shape Selection for Landsat Time Series of Forest Dynamics

Description

Author(s)

References

A 21 by 7 Matrix Storing Edf0 Vectors

Description

Usage

Format

See Also

f2a.map.jpeg

Description

Usage

Arguments

Details

Value

Author(s)

See Also

Examples

Raster based flat-to-annual function.

Description

Usage

Arguments

Details

Value

Author(s)

See Also

Examples

Raster based flat-to-parameter function.

Description

Usage

Arguments

Details

Value

Author(s)

See Also

Examples

flat2annual

Description

Usage

Arguments

Details

Value

Author(s)

Examples

flat2parameter

Description

Usage

Arguments

Details

Value

Author(s)

Examples

Get the Edf0 Vector for the Shape Routine

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Examples

The Plot Routine for an Object of the Shape Routine

Description

Usage

Arguments

Value

Author(s)

See Also

Examples

Shape Selection

Description

Usage

Arguments

Details

Value

Author(s)

References