[R-sig-Geo] Adaptive Savitzky Golay Filter function
Matteo Mattiuzzi
matteo.mattiuzzi at boku.ac.at
Wed Jun 12 12:58:04 CEST 2013
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.
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() 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
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list