get_cover_difference: Calculate quantiles for change in vegetation cover between...

In mbedward/CERMBlidar: A collection of functions to process aerial LIDAR data

Calculate quantiles for change in vegetation cover between two times

Description

This function uses the same approach for Bayesian estimation of vegetation cover as described for function get_cover_quantiles, but applied to two rasters representing LiDAR point counts for the same area at two different times. Argument rcounts0 is the raster for the earlier time t0, while argument rcounts1 is the raster for the later time t1. Each input raster should have two or more bands, with the first band representing ground point counts and subsequent bands representing point counts for vegetation strata arranged in increasing height order (i.e. ground is band 1). For a given pixel and stratum, the distribution of credible cover change values from time t0 to time t1 is found via a simulation approach (see Details). This is a very time consuming process, so rather than computing change quantiles for all pixels there is an option to restrict the calculations to a set of sample locations which can either be specified as either a two-column matrix or data frame of X-Y coordinate values, or an sf spatial data frame of point features.

Usage

get_cover_difference(
  rcounts0,
  rcounts1,
  probs = c(0.05, 0.5, 0.95),
  samplelocs = NULL,
  backgroundNA = TRUE,
  beta_prior = c(1, 1),
  nsim = 1e+05,
  seed = NULL
)

Arguments

rcounts0	Point counts raster for the first (earlier) time. Normally, this will be a terra package SpatRaster object generated by function get_stratum_counts(), but it may also be a RasterStack or RasterBrick (raster package) object. Important: the bands must be arranged in increasing height order, i.e. ground layer, then understorey layer(s), then overstorey layer(s).
rcounts1	Point counts raster for the second (later) time. This must have the same dimensions (rows, columns and number of bands) as rcounts0. Generally it will also have the same spatial bounds, but this is not checked to allow the function to be used for other types of comparisons (e.g. space for time substitution).
probs	A numeric vector of one or more probability values specifying the quantiles of vegetation cover change to calculate. The default c(0.05, 0.5, 0.95) returns quantiles for the median and the central 90% interval of cover change.
samplelocs	Spatial locations of raster pixels to sample. Can be a two-column matrix or data frame of X-Y coordinates, or an sf data frame of point features. Note that the default value (NULL) will result in change estimates being calculated for all pixels. Only do this if you have a lot of time on your hands! If providing sample locations as a two-column matrix or data frame, the two point count rasters must be in the same projection otherwise an error is raised. The projection is then assumed to apply to the sample coordinates.
backgroundNA	Controls what pixel value should be returned for a stratum when there are no LiDAR points for the stratum , any lower strata or ground level in one or both point counts rasters. If set to TRUE (the default), such cases will have missing (NA) pixel values for all quantiles for the stratum. If FALSE, quantiles for change values will be based on the prior beta distribution.
beta_prior	(Expert use only!) Either a numeric vector of two values defining the prior beta distribution to use for all strata, or a named list of two-element vectors where names match raster layer names for vegetation strata (case-insensitive). The default is c(1,1) for a beta(1,1) distribution which is equivalent to a uniform distribution on the (0,1) interval for all strata. This is probably what you want unless (a) you really know what you are doing and (b) there are particular reasons, i.e. previous studies or expert knowledge, to expect cover values for one or more strata to follow alternative distributions. For the list option, the default c(1,1) prior parameters will be used any strata not matched to list element names. The function beta_prior_helper can be used to explore candidate distributions given an expected cover value. See the examples and the Details section for further explanation.
nsim	Number of iterations for simulating the difference between the two beta distributions representing possible cover values at times t0 and t1. Default is 1e5.
seed	An optional integer value to use as an initial seed for random number generation when simulating the differences between beta distributions representing possible cover values for a pixel and stratum at the two times. This is useful when exactly reproducible results are preferred. If NULL (default), results will vary slightly between runs with the same input data.

Details

As described for function get_cover_quantiles, the estimation of stratum vegetation cover is based on treating the discrete LiDAR returns as a binomial sampling process, where the probability of a return from vegetation in the stratum being considered corresponds to vegetation cover. Under this model, the distribution of possible cover values that could have generated the observed point counts for a pixel will follow a beta distribution. Change in vegetation cover can therefore be estimated by taking the distribution of differences between the beta distribution for time 2 minus that for time 1. This function uses a simulation approach to approximate the distribution of differences.

Value

If estimating change for all pixels (samplelocs=NULL), the function returns a SpatRast (terra package) raster object with a layer for each requested quantile of vegetation cover change within each stratum. For example, calculating the default median plus 90% bounds (0.05 and 0.95 quantiles) of cover change for 5 vegetation strata would return a raster object with 15 layers. Layers are arranged by stratum, then by quantile. Layer names have the form stratum_dq_prob, e.g. 'TallShrub_dq_0.5' for the median (50%) quantile of cover change. If estimating change for a subset of pixels specified via the samplelocs argument, an sf spatial data frame is returned with a row for each sample location and columns for stratum x quantile change estimates. Column names have the same form as that described for raster band names.

See Also

get_cover_quantiles

Examples

## Not run: 
# Given two point count rasters, r0 (earlier time t0) and r1 (later time t1),
# calculate quantiles for difference in stratum cover at t1 compared to t0.
#
# Randomly select 100 pixels and create a set of point features
library(sf)
library(terra)

icells <- sample(1:ncells(r0), size = 100)
xy <- as.data.frame( xyFromCell(icells) )
xy <- st_as_sf(xy, coords = 1:2)

# Set the coordinate reference system for sample points. This is optional
# if r0 and r1 are in the same projection.
st_crs(xy) <- st_crs(r0)

# Quantiles of cover change (time t1 relative to t0) corresponding to
# median, 50\
# column for each stratum x quantile.
#
probs <- c(0.05, 0.25, 0.5, 0.75, 0.95)
dat_change <- get_cover_difference(r0, r1, probs, sample = xy, seed = 42)


# If the point count rasters are not too large, or you are happy to wait,
# you can calculate change estimates for all pixels by leaving the
# sample argument at its default of NULL. In this case, a terra::SpatRaster
# object is returned with a band for each stratum x quantile.
#
r_change <- get_cover_difference(r0, r1, probs, seed = 42)

## End(Not run)

mbedward/CERMBlidar documentation built on April 10, 2024, 2:05 p.m.