Type:	Package
Title:	Shape Selection for Landsat Time Series of Forest Dynamics
Version:	1.7
Date:	2023-08-19
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.
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
Depends:	coneproj (≥ 1.6), raster (≥ 2.3-40), R (≥ 3.0.2)
NeedsCompilation:	no
Suggests:	stats, MASS, graphics, grDevices, utils
Author:	Xiyue Liao [aut, cre], Mary Meyer [aut], Elizabeth Freeman [aut], Gretchen Moisen [aut]
Maintainer:	Xiyue Liao <xliao@sdsu.edu>
Repository:	CRAN
Packaged:	2023-08-19 16:05:16 UTC; xliao
Date/Publication:	2023-08-19 16:42:35 UTC

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 <xiyue@rams.colostate.edu>

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")

Format

A 21 by 7 matrix.

See Also

shape, ShapeSelectForest-package

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)

Arguments

Details

Creates one jpeg for each year in years.

Value

Returns nothing

Author(s)

Liz Freeman

See Also

f2a.raster

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)

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)

Arguments

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

See Also

flat2annual

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)

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)

Arguments

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:

Value

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

Author(s)

Liz Freeman

See Also

flat2annual

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)

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)

Arguments

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)

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)

Arguments

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)

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)

Arguments

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.

See Also

shape

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)

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)

Arguments

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

See Also

shape

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)

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)

Arguments

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

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.

See Also

plotshape, edf0s

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)

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)

Arguments

Value

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.

See Also

shape

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)

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")

Format

A 26 by 36 matrix.

Source

US Forest Service

See Also

shape, plotshape