[R-sig-Geo] Adaptive Savitzky Golay Filter function
matteo.mattiuzzi at boku.ac.at
Tue Jun 11 16:44:38 CEST 2013
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!
>>> Nicolas Bories 06/11/13 1:21 AM >>>
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() 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
> # 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, 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,
>>>> 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 !
> [[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