[R-sig-Geo] MODIS EVI time series missing data interpolation (linear or whittaker.raster)

[Matteo Mattiuzzi]

Mohammad,read answers inline

>>> Mohammad Abdel-Razek <abdelrazek at uni-bonn.de> 10/13/13 3:46 PM >>>
Dear All,
I am analyzing MODIS EVI time sereis for 2004-2010. Preprocessing 
included: reprojecting to WGS84 and export EVI to Geotiff (along with 
blue reflectance and view zenith angel). Then I filtered out cloudy  
(blue reflectance more than 10%) and off-nadir (angel > 32.5) pixels. 
Now, I am trying to interpolate the missing values, so I tried the 
following:

- MODIS EVI (MOD13Q1??) eventually has quality layers with a good cloud info! (see the MODIS:::detectBitinfo function:  detectBitInfo("MOD13Q1") # the VI usefulness. and MODIS:::extractBits to access this info.

#first linear interpolation
#convert Tiffs to points (rasterToPoints)
 > r01 <- raster("2004001_evi_filtered.tif")/10000
 >r01_p <- rasterToPoints(r89,fun=function(x){x>-0.9999})
#Then bind EVI columns, along with X,Y coordinates:
 >evi <- rbind ( EVI2004_1_1=r01_p[,3], EVI2004_1_17=r02_p[,3], 
.......,X= r153_p[,1], Y= r153_p[,2] )
#set 0 values to NA
 >evi[evi == 0] <- NA
 >evi=interpNA(evi, method = "linear")

- please be precise, in which package is 'interpNA'

then I transposed the evi to plot it with plotKML


Now I am trying to use alternative interpolation method for better 
results. MODIS package offer whittaker.raster filtering algorithm

 >evi2 = preStack(pattern = "*.tif", path = "path/to/file")
 >timeInfo <- orgTime(evi2,nDays=16,begin="2004001",end="2010241",pillow=40)
 >evi2 <- preStack(files=evi2,timeInfo=timeInfo)
 >beginCluster(type="SOCK",exclude="MODIS") # See: ?beginCluster
 >system.time(whittaker.raster(evi2,timeInfo=timeInfo,lambda=1000))

My questions are:
1. Is the modified whittaker filter functionreally suitable for 
interpolation in my case?
- yes the whittaker implementation in the MODIS package is suitable for all (?) time series where the bias is negative, where you can expect that if data is of bad quality (ie some clouds) the error is negatively biased, or that data with higher values is probably of better quality than lower values. For more details please check the publication in ?whittaker.raster. 

2. How to destack the resulting GeoTiffs (from whittaker.raster)? I want 
to plot the series using plotKML to visualize the results
- yes I know this should be possible, but isn't implemented, but you could do writeRaster(...,bylayer=T) after filtering.   


3. Some researchers suggested wavelet transformation to smooth the data, 
is that comprable to whittaker function?
- Sorry I don't know, but I think recent publications are going more in direction of spline like methods (whittaker, savitzky-golay, smooth spline,...).

Thank you in advance

Mohammad

p.s. I am new to r and would appreciate simple explanations!

Good luck
Matteo



-- 
PhD Candidate
Institute of Crop Science and Resource Protection
Crop Science Research Group
Katzenburgweg 5 - 53115 Bonn - Germany
Tel.: +49 (0) 228 73 2030     Fax: +49 (0) 228 73 2870
abdelrazek at uni-bonn.de        http://www.lap.uni-bonn.de


    [[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