MODIS EVI time series missing data interpolation (linear or whittaker.raster)
Mohammad:you still leave open a lot of important details especially on how you processed the data. The function whittaker.raster() has 3 important data inputs: vi: vegetation index layer w: VI qualitiy layer, derivate of it or any custom indicator t: composite day of the year layer if you provide 't' the function 'repDoy' generates (internally) a vector where each pixel value is placed on the exact DOY of it acquisition (if not the date in the filename is used), the output is always equally spaced depending on your 'timeInfo' argument. So the linear interpolation should not be needed, and I would also say a less good solution. I don't know what parameters you have used: lambda? niter? outlier? pillow? I have no idea what the problem could be that causes the difference. Maybe post a self containing example (restricted area!) (but I currently have very little time to share). Matteo
Mohammad Abdel-Razek <abdelrazek at uni-bonn.de> 11/05/13 3:58 PM >>>
Matteo, Thanks for your reply 1. Yes it is MOD13Q1 series, 2. interpNA function I used is from timeSeries package 3. I re-did the filtering this time with MODIS pixel quality as you suggested, which is much easier than my method :) I used the whittaker on both filtered and raw evi time series, still it dosn't give me a clear indication of land use change (compared to published papers for the same point coordinates) What seems to be the problem is the interpolation /before /whittaker. In this research (http://www.dsr.inpe.br/laf/series/artigos/Freitas_Ramon_M_et_al_2011_JCIS_Virtual_laboratory_of_remote_sensing_time_series.pdf) the authors did the following: "The filtered data were then linear interpolated based on the date of the pixel of the image composition to provide equally spaced time series." They did it with Matlab, is there a way to do it in R? Thanks again Mohammad
On 21.10.2013 12:49, Matteo Mattiuzzi wrote:
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