whittaker.raster: Filter Vegetation Index with Modified Whittaker Approach

In MODIS: Acquisition and Processing of MODIS Products

Filter Vegetation Index with Modified Whittaker Approach

Description

Use a modified Whittaker filter function (see References) from package ptw to filter a vegetation index (VI) time series of satellite data.

Usage

whittaker.raster(
  vi,
  w = NULL,
  t = NULL,
  timeInfo = orgTime(vi),
  lambda = 5000,
  nIter = 3,
  outputAs = "single",
  collapse = FALSE,
  prefixSuffix = c("MCD", "ndvi"),
  outDirPath = ".",
  outlierThreshold = NULL,
  mergeDoyFun = "max",
  ...
)

Arguments

vi	Raster* or character file names, sorted VI. Use preStack() functionality to ensure the right input.
w	Raster* or character file names. In case of MODIS composite, the sorted 'VI_Quality' layers.
t	Raster* or character file names. In case of MODIS composite, the sorted 'composite_day_of_the_year' layers. If missing, the date is determined using 'timeInfo'.
timeInfo	Output from orgTime().
lambda	character or integer. Yearly lambda value passed to ptw::whit2(). If set as character (i.e., lambda = "600"), it is not adapted to the time series length, but used as a fixed value (see Details). High values = stiff/rigid spline.
nIter	integer. Number of iterations for the upper envelope fitting.
outputAs	character, organization of output files. "single" (default) means each date one RasterLayer; "yearly" a RasterBrick for each year, and "one" one RasterBrick for the entire time series.
collapse	logical. Collapse input data of multiple years into one single year before filtering.
prefixSuffix	character, file naming. Names are dot-separated: paste0(prefixSuffix[1], "YYYDDD", lambda, prefixSuffix[2], ".defaultFileExtension").
outDirPath	character, output path. Defaults to the current working directory.
outlierThreshold	numeric in the same unit as 'vi', used for outlier removal (see Details).
mergeDoyFun	Especially when using collapse = TRUE, multiple measurements for one day can be present. Here you can choose how those values are merged to one single value: "max" uses the highest value, "mean" or "weighted.mean" use mean() or stats::weighted.mean().
...	Arguments passed to raster::writeRaster() (except for 'filename').

Details

The argument 'lambda' is passed to MODIS:::miwhitatzb1. You can set it as yearly 'lambda', which means that it doesn't matter how long the input time series is because 'lambda' is adapted to it with: lambda * ('length of input time series in days' / 365). The input length can differ from the output because of the 'pillow' argument in orgTime(). But it can also be set as character (i.e., lambda = "1000"). In this case, the adaption to the time series length is not performed.

Value

A Whittaker-smoothed RasterStack.

Note

Currently tested on MODIS and Landsat data. Using M*D13, it is also possible to use the 'composite_day_of_the_year' and the 'VI_Quality' layers.

Author(s)

Matteo Mattiuzzi and Agustin Lobo

References

Modified Whittaker smoother, according to Atzberger & Eilers 2011 International Journal of Digital Earth 4(5):365-386, doi: 10.1080/17538947.2010.505664. Implementation in R: Agustin Lobo 2012

See Also

smooth.spline.raster(), raster::raster().

Examples

## Not run: 
# The following function will download bit more than 1 year of MOD13A1 (~180mB) and therefore
# take while to execute! Data will be downloaded to options("MODIS_localArcPath") and processed 
# to 'paste0(options("MODIS_outDirPath"),"fullCapa")'
# You need to extract: 'vegetation index', 'VI_Quality layer', and 'composite day of the year',
# this is expressed by the argument 'SDSstring'
runGdal(product="MOD13A2",begin="2004340",extent="ireland",end="2006020", job="fullCapa",
SDSstring="101000000010")
path <- paste0(options("MODIS_outDirPath"),"fullCapa")

# the only obligatory dataset is the vegetatino index 
# get the 'vi' data in the source directory: 
vi <- preStack(path=path, pattern="*_NDVI.tif$")

# "orgTime" detects timing information of the input data and generates based on the arguments
# the output date information. 
# For spline functions (in general) the beginning and the end of the time series
# is always problematic. So there is the argument "pillow" (default 75 days) that adds
# (if available) some more layers on the two endings.
timeInfo <- orgTime(vi,nDays=16,begin="2005001",end="2005365",pillow=40)

# now re-run "preStack" with two differences, 'files' (output of the first 'preStack' call)
# and the 'timeInfo'
# Here only the data needed for the filtering is extracted:
vi <- preStack(files=vi,timeInfo=timeInfo)

whittaker.raster(vi,timeInfo=timeInfo,lambda=5000)

# if the files are M*D13 you can use also Quality layers and the composite day of the year:
wt <- preStack(path=path, pattern="*_VI_Quality.tif$", timeInfo=timeInfo)
# can also be already stacked:
inT <- preStack(path=path, pattern="*_composite_day_of_the_year.tif$", timeInfo=timeInfo)

whittaker.raster(vi=vi, wt=wt, inT=inT, timeInfo=timeInfo, lambda=5000, overwrite=TRUE)

## End(Not run)

MODIS documentation built on Jan. 6, 2023, 5:10 p.m.