[R-sig-Geo] Adaptive Savitzky Golay Filter function
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
# do these two shp files match visualized in a GIS? (with enabled 'on
the fly' reprojection!)
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
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
# 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()'
>>> Nicolas Bories <nicolas.bories at orange.fr> 06/11/13 6:51 PM >>>
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
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
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"
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
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
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
The MODIS package seems to be very useful and very powerful !
Nonetheless I encounter some problems with the runGdal() function and
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 :
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
+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
> # 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
> # this downloads, converts (HDF to GDAL formats), crops, warps.
> see?runGdal for details. You can limit the data sung the SDSstring
> # 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!
> 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,
>>>> 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 !
> [[alternative HTML version deleted]]
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
[[alternative HTML version deleted]]
More information about the R-sig-Geo