[R-sig-Geo] Adaptive Savitzky Golay Filter function

[Matteo Mattiuzzi]

Dear Nicolas,


The MODIS package on R-forge is intended to do exactly this. There are
implemented a non iterative smooth.spline filtering and the iterative
Whittaker filter (Atzberger and Eilers 2012).
MODIS offers also a lot of data acquisition and pre processing
capabilities but if you intend to use those capabilities you have to do
a little bit of set up first. 


Install the package and run MODISoptions(), this function is intended to
guide you through. Be aware that on Linux it is much simpler, but give
it a try (and give me some feedback pls).


install.packages("MODIS", repos="http://R-Forge.R-project.org")


########
# on Linux avoid to run the next part as a non typical user ie as root,
because some directory is created and you will not access it as normal
after.
########


library(MODIS)
MODISoptions()


# You are on a Windows machine you need to install a HDF supporting GDAL
 or MRT (run MODIS:::checkTools() for help). If MODIS continues to
bother you with thinks redo the steps.
# The entire processing should look approximately like that:


getProduct() # to list available once


runGdal(product="M.D13Q1",begin="2005001",end="2006001", job="testJob")
# this downloads, converts (HDF to GDAL formats), crops, warps.
see?runGdal for details. You can limit the data sung the SDSstring
argument.
# if you have MRT
runMrt(product="M.D13Q1",begin="2005001",end="2006001", job="testJob") #
this downloads, converts (HDF to GDAL formats), crops, warps. see?runMrt
for details, but I suggest to use gdal!


setwd(paste0(options()$MODIS_outDirPath,"/testJob"))
dir() # your data


 # look for data
vi     <- preStack(pattern="*_NDVI.tif$")
qa   <- preStack(pattern="*_VI_Quality.tif$")
time <- preStack(pattern="*_composite_day_of_the_year.tif$")
    
crono <- orgTime(vi) # here you decide the time steps for the output


# order data by time, and use only what is defined by orgTime    
vi      <- stack(preStack(files=vi,timeInfo=crono))
time <- stack(preStack(files=time,timeInfo=crono))
qa    <- stack(preStack(files=qa,timeInfo=crono))


# the weighting is generated internally in whittaker.raster but you
could also do it separately for more control see ?makeWeights
# qa <-  makeWeights(qa, bitShift=2, bitMask=15, threshold=6) # for
bitShift and bitMask info on other products run detectBitInfo("MOD13Q1")



# filtering 
whittaker.raster(vi=vi, wt=qa, inT=time, timeInfo=crono, lambda=1000,
overwrite=TRUE) # 


Good luck,
Matteo



>>> Nicolas Bories <nicolas.bories at orange.fr> 06/07/13 4:29 PM >>>
Dear all,

I'm trying to use the Adaptive Savitzky Golay Filterfor smoothing and
gap filling a time series of remote sensing observations  (MODIS).
For this exercise, I need to use a vector weighting strategy for a
smoothing that takes into account the pixelreliability. The TIMESAT code
(L. Eklundh, P. Jönsson,2011) seems to be ideal for this. But, it is
limiting my specific needs because it is available as an executable and
I cannot customize it according to my needs. So, I am wondering if there
is any other alternative function/tools to do this task.

Thanks in advance for any suggestions !

Nicolas.





    [[alternative HTML version deleted]]