[R-sig-Geo] Processing best quality MODIS pixels only (MODIS package)

[Amit Boshale]

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.boshale at 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.mattiuzzi at 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.boshale at 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]]