MODIS (version 1.0.0)

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 serie 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 filenames, sorted VI. Use preStack functionality to ensure the right input.

w

Raster* or character filenames. In case of MODIS composite, the sorted 'VI_Quality' layers.

t

Raster* or character filenames. 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 whit2. If set as character (i.e., lambda = "600"), it is not adapted to the time serie length, but used as a fixed value (see Details). High values = stiff/rigid spline.

nIter

integer. Number of iteration for the upper envelope fitting.

outputAs

character, organisation 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 composed dot-separated: paste0(prefixSuffix[1], "YYYDDD", lambda, refixSuffix[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 the mean if no weighting "w" is available, and weighted.mean if it is.

...

Arguments passed to writeRaster (except for filename).

Value

The filtered data and a text file with the dates of the output layers.

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 serie is because lambda is adapted to it with: lambda*('length of input timeserie 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. Implementation in R: Agustin Lobo 2012

See Also

smooth.spline.raster, raster.

Examples

# The full capacity of the following functions is currently avaliable only with M*D13 data.
# !! The function is very new, double check the result !!

# You need to extract: 'vegetation index', 'VI_Quality layer', and 'composite day of the year'.
# runGdal(product="MOD13A2",begin="2004340",extent="sicily",end="2006070",
# job="fullCapa",SDSstring="101000000010")
# You can download this dataset from (2.6 MB): 
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)