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

Matteo Mattiuzzi matteo.mattiuzzi at boku.ac.at
Wed Jun 12 12:58:04 CEST 2013


using the meuse dataset you had a shift of 100km? or was it using your
'zone'? Is the problem in your 'zone' as I have no shift using the
Internally the specified extent is reprojected to LatLon using
rgdal:::spTransform() so maybe you have to check the CRS of your zone?
If you want to follow deeper in this argument you have to reveal more
about your 'zone'  
zoneLatLon <- spTransform(zone,CRS("+proj=longlat +datum=WGS84
+ellps=WGS84 +towgs84=0,0,0"))


# do these two shp files match visualized in a GIS? (with enabled 'on
the fly' reprojection!)

# Filtering
there are two helper function intended to organize the input data for
whittaker.raster() and smooth.spline.raster().
preStack and orgTime

the first ('preStack') has two rules:
1. to look for data in 'path' with a given pattern
2. to sort files by date and remove superfluous data (this info comes
from orgTime)

the secound ('orgTime') checks the available input data (from the first
preStack run), and generates a output structure based on your needs
(mainly these args: 'nDays','begin','end' and 'pillow') and the
available input data.

Example (in ?whittaker.raster):
path <- 'dir to your data'
# with the first run of preStack all files in 'path' with the ending
*_NDVI.tif$ are listened
vi <- preStack(path=path, pattern="*_NDVI.tif$")

# with orgTime you say: retrieve and sort ascending by time, all
filenames (vi), that are needed to process the period begin to end,
adding (_if available_) a pillow of data on the begin and on the end you
your time serie to ensure a high quality filtering also on the temporal
timeInfo <- orgTime(vi,nDays=16,begin="2005001",end="2005365",pillow=40)

# now just get the data in vi, selected and sorted by timeInfo     
vi <- preStack(files=vi,timeInfo=timeInfo)

# now the simples run of whittaker.raster
# 'vi' is a vector of filenames, but you can also stack it:
# vi <- stack(vi)
whittaker.raster(vi=vi,timeInfo=timeInfo) # the result is written to a
file in the current dir as specified by the default argument in
whittaker.raster outDirPath="./"

# if you what 
use also the weighting info from VI_Quality and the precise day of the
year for the "X-Axis" [t], you have to organize also the following files
qa   <- preStack(path=path, pattern="*_VI_Quality.tif$")
time <- preStack(path=path, pattern="*_composite_day_of_the_year.tif$")

qa <- preStack(files=qa,timeInfo=timeInfo)

qa <- stack(qa)

time <- preStack(files=time,timeInfo=timeInfo)

time <- stack(time)

whittaker.raster(vi=vi, wt=qa, inT=time, timeInfo=timeInfo) # this
overwrites the result of the first filtering!
# in whittaker.raster and smooth.spline.raster there are functions that
are serializing the 'time' info and also extracting the bit-coded info
in qa, if you want to generate your own weighting rules based on qa you
need to run makeWeights() yourself and use it as input:

qa2 <- makeWeights(qa) # see ?makeWeights for the how to (important is:
threshold=XX, and decodeOnly=FALSE) see in the example what 'res' is
whittaker.raster(vi=vi, wt=qa2, inT=time, timeInfo=timeInfo) 

# NOTE you can also speedup the filtering using 'beginCluster()'
good luck, 

>>> Nicolas Bories <nicolas.bories at orange.fr> 06/11/13 6:51 PM >>>

Dear Matteo,

Thanks a lot for your help !
I installed MODIS 0.9-3.
First of all, I followed your example with meuse dataset.
The projection was kept but there is a shift biger than 100 km with the 
real position.
In a second time I tried your solution to avoid any resampling.
It's worked perfectly : no shift, no resampling !
Thank you !

Now, I have to understand how to get a long time serie smoothed and gap 
filled for each pixel.
It's tricky for me, because I don't understand clearly how to use 
whittaker.rasteRe: [R-sig-Geo] Adaptive Savitzky Golay Filter function
From:	 Matteo Mattiuzzi
To:	 Nicolas Bories; r-sig-geo at r-project.org
Date:	 Tuesday - June 11, 2013 4:44 PM
Subject: 	Re: [R-sig-Geo] Adaptive Savitzky Golay Filter function

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!


     Reply  Reply All  Forward   Move   Mark Unread    Delete   Print
Mail Properties
From:	Nicolas Bories <nicolas.bories at orange.fr>	Tuesday - June
11, 2013 1:21 AM
To:	r-sig-geo at r-project.org
Subject:	Re: [R-sig-Geo] Adaptive Savitzky Golay Filter function
Part.001 (1 KB)	   View
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() f> 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, converts (HDF to GDAL formats), crops, warps.
> 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
> # 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
> (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