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()
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
