repDoy: Repair MODIS "composite_day_of_the_year" SDS

In MODIS: Acquisition and Processing of MODIS Products

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)

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'.

Value

A matrix with sequential Julian dates.

Author(s)

Matteo Mattiuzzi

Examples

## Not run: 
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))

## End(Not run)
 

MODIS documentation built on Jan. 6, 2023, 5:10 p.m.