Dear Matteo,
You're right, I did a mistake. Now I have no shift in position and all
is ok.
I followed all your explanations and it's worked perfectly.
Thanks a lot again,
Nicolas.
On 12/06/2013 12:58, Matteo Mattiuzzi wrote:
> 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.bor...@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@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.bor...@orange.fr> Tuesday - June
> 11, 2013 1:21 AM
> To: r-sig-geo@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.bor...@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@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>
> [[alternative HTML version deleted]]
>
>
>
>
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
