Adaptive Savitzky Golay Filter function

Nicolas,


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
'meuse'. 
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'  
ie:
bbox(zone)
proj4string(zone)
zoneLatLon <- spTransform(zone,CRS("+proj=longlat +datum=WGS84
+ellps=WGS84 +towgs84=0,0,0"))

bbox(zoneLatLon)
proj4string(zoneLatLon)


# do these two shp files match visualized in a GIS? (with enabled 'on
the fly' reprojection!)
require(raster)
shapefile("z1",zone)
shapefile("z2",zoneLatLon)


#####
# 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$")
vi


# 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
borders. 
timeInfo <- orgTime(vi,nDays=16,begin="2005001",end="2005365",pillow=40)
timeInfo


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

# 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
plotting
whittaker.raster(vi=vi, wt=qa2, inT=time, timeInfo=timeInfo) 


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

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
BC:	
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:
http://ivfl-arc.boku.ac.at/owncloud/public.php?service=shorty_relay&id=QR5FZSxe0g)


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

require(rgdal)
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",
extent=zone)

####################
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
runGdal(product="MOD13Q1",begin="2002001",end="2002001",
job="testJob",tileH=tileH,tileV=tileV)

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


Matteo



     Reply  Reply All  Forward   Move   Mark Unread    Delete   Print
View   
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
 Attachments:	
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", 
job="testJob",extent=zone)
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 
4925025.8738480


2) runGdal(product="MOD13Q1",begin="2002001",end="2002001", 
job="testJob",extent=zone,pixelSize="asIn",resamplingType="near",outProj="+proj=sinu

+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
+datum=WGS84"


This produces TIFF files reprojected in LongLat


Can anyone please help me?


Many thanks in advance,


Nicolas.

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() 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",

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]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo