MODIS (version 1.2.11)

whittaker.raster: 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",
  ...
)

Value

A Whittaker-smoothed RasterStack.

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').

Author

Matteo Mattiuzzi and Agustin Lobo

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.

References

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

See Also

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

Examples

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