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

Matteo Mattiuzzi matteo.mattiuzzi at boku.ac.at
Tue Jun 11 16:44:38 CEST 2013

Dear Nicolas,

Fixed in MODIS 0.9-3. Thanks! (If you like to build it yourself you can
get the source from:

example with the meuse dataset:
coordinates(meuse) <- ~x+y
proj4string(meuse) <- CRS("+init=epsg:28992")

zone <- spTransform(meuse,CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0
+a=6371007.181 +b=6371007.181 +units=m +no_defs"))

# make sure your extent has an assigned projection string, if not do:

# proj4string(zone) <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
+b=6371007.181 +units=m +no_defs"

runGdal(product="MOD13Q1",begin="2002001",end="2002001", job="testJob",

Be aware, that if you use an 'extent' argument a certain re-sampling, in
form of a pixel shift, is nearly always performed, as a specified
'extent' may not matches the original MODIS grid! If you want absolutely
avoid any resampling you have to do something like that (processing the
entire TILE(s)):

tileH <- getTile(zone)$tileH
tileV <- getTile(zone)$tileV

and then you can use raster:::crop() to reduce the extents if you like!


>>> Nicolas Bories  06/11/13 1:21 AM >>>
Dear all,

The MODIS package seems to be very useful and very powerful !
Nonetheless I encounter some problems with the runGdal() function and 
extent argument.

I need to extract a part of "MOD13Q1" product between 2002 and 2005.
I also need to keep data in the same projection and the same resolution 
than hdf files, without resampling.
I tried 2 cases :

1) runGdal(product="MOD13Q1",begin="2002001",end="2006001", 
zone is a SpatialPointsDataFrame with same projection than hdf files

This returns :
Error in ReplProj4string(obj, CRS(value)) :
   Geographical CRS given to non-conformant data:  -94112.4108437 

2) runGdal(product="MOD13Q1",begin="2002001",end="2002001", 

+lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
zone is a SpatialPointsDataFrame  with CRS : "+proj=longlat

This produces TIFF files reprojected in LongLat

Can anyone please help me?

Many thanks in advance,


On 07/06/2013 17:24, Matteo Mattiuzzi wrote:
> 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
> a little bit of set up first.
> Install the package and run MODISoptions(), this function is intended
> 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
> 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
>   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",
> # 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, conv> 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
> # filtering
> whittaker.raster(vi=vi, wt=qa, inT=time, timeInfo=crono, lambda=1000,
> overwrite=TRUE) #
> Good luck,
> Matteo
>>>> Nicolas Bories  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
> (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
> I cannot customize it according to my needs. So, I am wondering if
> is any other alternative function/tools to do this task.
> Thanks in advance for any suggestions !
> Nicolas.
>      [[alternative HTML version deleted]]
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

    [[alternative HTML version deleted]]

More information about the R-sig-Geo mailing list