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

[Matteo Mattiuzzi]

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