MODIS (version 1.2.11)

repDoy: Repair MODIS "composite_day_of_the_year" SDS

Description

Currently works only for MODIS 16 days composites! In MODIS composites, the Julian dates inside the 'composite_day_of_the_year' SDS are referring always to the year they are effectively in. The problem is that the layer/SDS name from the last files from Terra and Aqua within a year can include dates from the following year and so starting again with 1. The problem occurs if you want to sort values of a time series by date (e.g. for precise time series functions). This function generates a sequential vector beginning always with the earliest SDS/layer date and ending with the total sum of days of the time series length.

Usage

repDoy(pixX, layerDate = NULL, bias = 0)

Value

A matrix with sequential Julian dates.

Arguments

pixX

matrix of values, usually derived from raster::as.matrix().

layerDate

If NULL (default), try to autodetect layer dates. If you want to be sure, use the result from extractDate() or orgTime().

bias

integer. Bias applied to all values in 'pixX'.

Author

Matteo Mattiuzzi

Examples

if (FALSE) {
tfs <- runGdal(product="M.D13A2", begin="2010350", end="2011016"
               , extent="Luxembourg", job="deleteme", SDSstring="100000000010")

ndviFiles <- grep("NDVI.tif$", unlist(tfs, use.names = FALSE), value = TRUE)
ndviFiles <- preStack(files = ndviFiles, timeInfo = orgTime(ndviFiles))
ndvi <- stack(ndviFiles)

doyFiles <- grep("composite_day_of_the_year.tif$"
                 , unlist(tfs, use.names = FALSE), value = TRUE)
doyFiles <- preStack(files = doyFiles, timeInfo = orgTime(doyFiles))
doy <- stack(doyFiles)

layerDates <- extractDate(doyFiles)

pixX <- 169

y <- ndvi[pixX]
print(x1 <- doy[pixX])
print(x2 <- repDoy(x1,layerDates))

# the plotting example is not really good. 
# To create a figurative example it would be necessary to dolwnload to much data! 
plot("",xlim=c(1,max(x1,x2)),ylim=c(0,2000),xlab="time",ylab="NDVI*10000")
lines(y=y,x=x1,col="red",lwd=3)
lines(y=y,x=x2,col="green",lwd=2)

# repDoy function is thought to be embedded in something like that:
tr <- blockSize(ndvi)

doyOk <- brick(doy)
doyOk <- writeStart(doyOk, filename='test.tif',  overwrite=TRUE)

for (i in 1:tr$n)
{
  pixX  <- getValues(doy,tr$row[i],tr$nrows[i])
  ok    <- repDoy(pixX,layerDates)
  doyOk <- writeValues(x=doyOk,v=ok,start=tr$row[i])
}
doyOk <- writeStop(doyOk)

unlink(filename(doyOk))
}