Dear Matteo,
I just found out something interesting. As I mentioned in my previous email, I
tried the code you proposed and could see the diffference from the old code,
basically a smoother plot. However, the plot didn't change corresponding to
differenct quality thresohold (1, 2, 3, 4 and 5).
Then I realized that the smoother time series was due to larger Lambda. It was
misspled in the code, hence MODIS package took a default value of 5000 rather
than the 500 I used in the code.
It might be due to the presence of another threshold (outlier removal) used in
MODIS package 0.10-18, hence the package ignored the value. In your newst
package (0.10-23) outlier threshold variable name is changed. But the new
package doesn't allow for using cluster (it throughs an error), so I couldn't
try.
Hope the remark could be helpful in any future improvments :D
Best,Amit
Amit Boshale <amit.bosh...@yahoo.com> schrieb am 13:49 Mittwoch, 20.Mai
2015:
Dear Matteo,
Many thanks for the help!
I just tried it on a small raster stack and it worked. The resulting time
series looks better. I will keep trying with other threshold values(2,3,etc.)to
compare the resulting graphs for improved visual interpretation.
Best,Amit
Matteo Mattiuzzi <matteo.mattiu...@boku.ac.at> schrieb am 21:53 Dienstag,
19.Mai 2015:
Dear Amit,
Thats correct, if you want to use only pixel of 0000 quality you have to set 0
as threshold (everything > than threshold is set to 0 weight). But for what I
remember, and this might explains your result with all 0s, is that 0000 does
not exist in MODIS data despite the fact that it is described so on the
website. You can use the decodeOnly=TRUE option to see the exact values of your
quality layers, but I guess you will not find any values lower than 1.
If you run whittaker.raster, you can feed also argument of the makeWeights
function. If you use a cluster this step is parallelized automatically.
I gave a short look to makeWeights now, probably I can do some little speed
improvements.
It should be possible to do:
library("MODIS")
path=("E:/NDVI/Indonesia")
#stack ndvi images, following MODIS package example
ndvi = preStack(pattern = "*_NDVI.tif", path = path)
timeInfo <- orgTime(ndvi,nDays=16,begin="2000049",end="2015049",pillow=40)
#Restack evi
ndvi <- preStack(files=ndvi,timeInfo=timeInfo)
# You don't need to stack here, it is done with stack(...quick=TRUE) within
whittaker.raster.
wt <- preStack(path=path, pattern="*_VI_Quality.tif$", timeInfo=timeInfo)
inT <- preStack(path=path, pattern="*_composite_day_of_the_year.tif",
timeInfo=timeInfo)
nodes <- 4
library(snow)
beginCluster(nodes)
system.time(whittaker.raster(vi=ndvi, wt = wt, inT = inT, timeInfo = timeInfo,
lamba = 500, bitShift = 2, bitMask = 15, threshold = 1))
endCluster()
Matteo
>>> Amit Boshale <amit.bosh...@yahoo.com> 05/19/15 4:54 PM >>>
Hello!
I use MODIS package to smooth MODIS NDVI raster time series.
Just wanted to try how would it look like if I used the best qulity pixels
only.makeWeights function took about two days and resulted in raster stack of
zeros. According to MODIS, the best qulity is 0, hence I used 0 as the
threshold, is that correct? Is there a multicore version of makeWeights
function?
Best,Amit
library("MODIS")path=("E:/NDVI/Indonesia")
#stack ndvi images, following MODIS package example
ndvi = preStack(pattern = "*_NDVI.tif", path = path)
timeInfo <- orgTime(ndvi,nDays=16,begin="2000049",end="2015049",pillow=40)
#Restack evi
ndvi <- preStack(files=ndvi,timeInfo=timeInfo)
#stack Vegtation index quality layers
VIqual <- stack(preStack(path=path, pattern="*_VI_Quality.tif",
timeInfo=timeInfo))
bitShift = 2
bitMask = 15
threshold=0 # based on MODIS documentation, quality 0 is the best
wt <- makeWeights(VIqual,bitShift,bitMask,threshold,decodeOnly=FALSE)
#stack composite day of the year layers
inT <- stack(preStack(path=path, pattern="*_composite_day_of_the_year.tif",
timeInfo=timeInfo))
nodes <- 4
library(snow)
beginCluster(nodes)
system.time(whittaker.raster(vi=ndvi, wt=wt, inT=inT, timeInfo=timeInfo,
lambda=500))
endCluster()
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
