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

[Matteo Mattiuzzi]

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


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




>>> 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", 
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.

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


    [[alternative HTML version deleted]]