interpolate_raster: Interpolate a Raster* time series to new time stamps...
In antoinestevens/MODISExtra: Download and extract MODIS imagery

Description Usage Arguments Details Value Note Author(s) References See Also Examples

Interpolate Raster* times series to new time stamps

interpolate_raster(x, w = NULL, t = NULL, timeInfo = orgTime(x),
       method = c("spline","linear","constant","nn","whittaker"),
       lambda = 5000, nIter= 3, outlierThreshold = NULL,
       df = 6,
       cl = NULL, filename=rasterTmpFile(),...)

`x`	A `Raster-class` object to interpolate
`w`	Optional, a `Raster-class` object with weight information. See e.g. `makeWeights`.
`t`	Optional, a `Raster-class` object with MODIS 'composite day of the year' data.
`timeInfo`	object returned by `orgTime` defining the input and output time stamps to interpolate to.
`method`	character `vector` of length 1 with one of the following value: 'spline' for smoothing splines ,'linear' for linear approximation ,'constant' for constant approximation, 'nn' for nearest neighbour , 'whittaker' for whittaked smoothing.
`lambda`	Yearly lambda value. Default is 5000. See `whittaker.raster`.
`nIter`	numeric `vector` of length 1 giving the number of iterations for the upper envelope fitting ('whittaker' `method`). Default is 3. See `whittaker.raster`.
`outlierThreshold`	Numeric. If provided, `outlierTreshold` allows to remove outliers before smoothing by fitting a whittaker function and removing values with residuals higher than `outlierTreshold`.
`df`	Yearly degree of freedom for 'spline' `method`. Default is 6. See `smooth.spline.raster`
`cl`	cluster object for parallel processing (do not use it if beginCluster() has been called. Default is `NULL`
`filename`	name of the file where the interpolated raster is saved. Default is created through `rasterTmpFile`
`...`	arguments passed to `writeRaster`

'linear' and 'constant' method use the approx interpolation functions while 'spline' use splinefun.

A Raster-class object with values interpolated to the output time series (as defined in timeInfo)

Most of the code is borrowed from smooth.spline.raster and whittaker.raster. Compared to the MODIS version, code is simplified and made general for reflectance values

Antoine Stevens & Matteo Mattiuzzi

Atzberger, C., and Eilers, P.H.C. (2011). A time series for monitoring vegetation activity and phenology at 10-daily time steps covering large parts of South America. International Journal of Digital Earth 4, 365-386.

gapfill_raster

## Not run: 
# Get and extract NDVI MOD13A2 for 2009
runGdal(product="MOD13A2", begin = "2009.01.01", end = "2009.12.31",
        tileH = 19, tileV = 5 ,
        SDSstring =  "100000000000",job="H19V5",
        gdalPath = "C:/OSGeo4W64/bin")
# process and convert to raster
path <- "~/MODIS_DATA/PROCESSED/H19V5"
r <- convert_modis(path = path,pattern = "16_days_NDVI",
                   begin = "2009.01.01", end = "2009.12.31",extractAll=T)
str(r)
r_NDVI <- r$raster_NDVI
# Smooth by splines
r_smooth <- interpolate_raster(x = r_NDVI, method = "spline")
# plot
x <- r$ti_NDVI$inputLayerDates
y1 <- r_NDVI[1000]
y2 <- r_smooth[1000]
plot(x,y1)
lines(x,y2,col = "red")
# now, using weights and composite day of the year to smooth
# first extract data
runGdal(product="MOD13A2", begin = "2009.01.01", end = "2009.12.31",
        tileH = 19, tileV = 5 ,
        SDSstring =  "001000000010", # extract VI_usefulness and composite_day_of_the_year
        job="H19V5",
        gdalPath = "C:/OSGeo4W64/bin")
cdoy <- convert_modis(path = path ,pattern = "16_days_composite",
                      begin = "2009.01.01", end = "2009.12.31")
QF <- convert_modis(path = path ,pattern = "16_days_VI_Quality",
                   begin = "2009.01.01", end = "2009.12.31", convertDN = FALSE)
# Convert QF to 'usefulness'
usefulness <- convert_qf_modis(r = QF,qf = "VI_usefulness")
# Convert usefulness to weights (see reference section)
w <- (15-usefulness)/15
r_smooth2 <- interpolate_raster(x = r_NDVI, w = w, t = cdoy, method = "spline")
y3 <- r_smooth2[1000]
lines(x,y3,col="green")
# One could aslo interpolate to new time stamps
ti <- orgTime(r_NDVI,nDays = 10,
              begin = "2009.01.01", end = "2009.12.31") # resample every 10 days
r_smooth_resampled <- interpolate_raster(x = r$raster_NDVI, w = w, t = cdoy,
                                         timeInfo = ti, method = "spline")
y4 <- r_smooth_resampled[1000]
points(ti$outputLayerDates,y4,pch=3)

## End(Not run)