# Author: Matteo Mattiuzzi, matteo.mattiuzzi@boku.ac.at
# Date : August 2012
# Licence GPL v3
# maybe to add: derivate=1
# vx=val[index,u];wx=wtu[index,u];tx=inTu[index,u]
# tmtr <- cbind(vx,wx,tx)
unifyDoubleWM <- function(vx,wx,tx)
{
tx <- as.numeric(tx)
vx <- as.numeric(vx)
wx <- as.numeric(wx)
double <- tx[duplicated(tx)]
if(length(double)>0)
{
double <- unique(double)
for(i in seq_along(double))
{
inx <- which(tx==double[i])
vx[inx[1]] <- weighted.mean(vx[inx],w=wx[inx])
wx[inx[1]] <- max(wx[inx])
vx <- vx[-inx[-1]]
wx <- wx[-inx[-1]]
tx <- tx[-inx[-1]]
}
}
list(vx=vx,wx=wx,tx=tx)
}
unifyDoubleMX <- function(vx,wx,tx)
{
tx <- as.numeric(tx)
vx <- as.numeric(vx)
wx <- as.numeric(wx)
double <- tx[duplicated(tx)]
if(length(double)>0)
{
double <- unique(double)
for(i in seq_along(double))
{
inx <- which(tx==double[i])
mx <- which.max(wx[inx])
vx <- vx[-inx[-mx]]
wx <- wx[-inx[-mx]]
tx <- tx[-inx[-mx]]
}
}
list(vx=vx,wx=wx,tx=tx)
}
doCollapse <- function(tx,pillow)
{
ord <- order(tx)
txS <- tx[ord]
t0 <- 365 - pillow
tS <- ord[txS >= t0]
tE <- ord[txS <= pillow]
s0 <- txS[txS >= t0] - t0
s1 <- txS + pillow
s2 <- txS[txS <= pillow] + 365 + pillow
list(order=c(tS,ord,tE),sequence=c(s0,s1,s2)+1)
}
whittaker.raster <- function(vi, w=NULL, t=NULL, timeInfo = orgTime(vi), lambda = 5000, nIter= 3, outputAs=c("single","yearly","one"), collapse=FALSE, prefixSuffix=c("MCD","ndvi"), outDirPath=getwd(), outlierThreshold=NULL, mergeDoyFun="max", ...)
{
if(!require(ptw))
{
stop("In order to use the the whittaker filter please install 'ptw': install.package('ptw')")
}
# debug
# w=wt; t=inT; groupYears=TRUE; lambda = 5000; nIter= 3; outDirPath = "./"; collapse=FALSE; opt <- MODIS:::combineOptions();removeOutlier=FALSE;mergeDoyFun="max";opts$outDirPath <- MODIS:::setPath(outDirPath)
# w=NULL; t=NULL; groupYears=TRUE; lambda = 5000; nIter= 5; outDirPath = "./SUB/";collapse=TRUE; opts <- MODIS:::combineOptions();removeOutlier=FALSE;mergeDoyFun="max"
# opt <- list(bitShift=2,bitMask=15,threshold=6)
# opt <- list(bitShift=2,bitMask=15)
# opt <- list()
opts <- combineOptions(...)
outDirPath <- setPath(outDirPath)
bitShift <- opts$bitShift
bitMask <- opts$bitMask
threshold <- opts$threshold
dataFormat <- opts$dataFormat
rasterOut <- toupper(writeFormats())
if(toupper(dataFormat) %in% rasterOut[,"name"])
{
dataFormat <- raster:::.defaultExtension(dataFormat)
} else
{
stop("Argument dataFormat='",dataFormat,"' is unknown/not supported. Please run 'writeFormats()' (column 'name') so list available dataFormat's")
}
if (is.null(opts$minDat))
{
minDat <- 3 # 3 is very small!
} else
{
minDat <- opts$minDat
}
if(collapse)
{
if(timeInfo$call$nDays=="asIn")
{
stop("'orgTime' argument nDays='asIn' is not allowed when using collapse=TRUE")
}
fitt <- seq(as.numeric(format(min(timeInfo$outputLayerDates),"%j")),as.numeric(format(max(timeInfo$outputLayerDates),"%j")),by=timeInfo$call$nDays) + timeInfo$call$pillow
} else
{
fitt <- timeInfo$outSeq
}
inlam <- lambda
if (is.character(lambda))
{
cat("Using fixed 'lambda':",lambda,"\n")
nameL <- "fL"
} else
{
if (collapse)
{
lambda <- lambda*((365 + 2*timeInfo$call$pillow)/365)
cat("Yearly 'lambda' is:",inlam,"\nNow changed with lambda*((365+2*pillow)/365) to:",lambda,"\n")
} else
{
lambda <- lambda*((max(timeInfo$inSeq) - min(timeInfo$inSeq) -1)/365)
cat("Yearly 'lambda' is:",inlam,"\nNow changed with lambda*('length of input data period in days'/365) to:",lambda,"\n")
}
nameL <- "yL"
}
inlam <- round(inlam) #from here on used only in outputfilename
lambda <- as.numeric(lambda)
if(!inherits(vi,"Raster"))
{
vi <- stack(vi,quick=TRUE)
}
if(!inherits(w,"Raster") & !is.null(w))
{
w <- stack(w,quick=TRUE)
}
if(!inherits(t,"Raster") & !is.null(t))
{
t <- stack(t,quick=TRUE)
}
if(is.null(opts$datatype))
{
opts$datatype <- dataType(vi[[1]])
}
if(opts$datatype == "FLT")
{
doround <- FALSE
} else
{
doround <- TRUE
}
if (is.null(opts$overwrite))
{
opts$overwrite <- raster:::.overwrite()
}
outputAs <- tolower(outputAs)
if (collapse) # use only DOY obmitt year (all years into one)
{
d <- sprintf("%03d",seq(as.numeric(format(min(timeInfo$outputLayerDates),"%j")),as.numeric(format(max(timeInfo$outputLayerDates),"%j")),by=timeInfo$call$nDays))
if(outputAs=="single")
{
oname <- paste0(outDirPath,prefixSuffix[1],".",d,".",prefixSuffix[2],dataFormat)
b <- vector(mode="list",length=length(oname))
for(a in seq_along(oname))
{
b[[a]] <- raster(vi)
b[[a]] <- writeStart(b[[a]], filename=oname[a], datatype=opts$datatype, overwrite=opts$overwrite)
}
} else
{
b <- list()
b[[1]] <- writeStart(brick(raster(vi),nl=length(d)), filename=oname, datatype=opts$datatype, overwrite=opts$overwrite)
names(b[[1]]) <- paste0("doy",d)
}
} else if (outputAs=="yearly")
{
b <- vector(mode="list",length=length(unique(format(timeInfo$outputLayerDates,"%Y"))))
for (a in seq_along(unique(format(timeInfo$outputLayerDates,"%Y"))))
{
y <- unique(format(timeInfo$outputLayerDates,"%Y"))[a]
oname <- paste0(outDirPath,prefixSuffix[1],".year",y,".",nameL,inlam,".",prefixSuffix[2],dataFormat)
names <- timeInfo$outputLayerDates[format(timeInfo$outputLayerDates,"%Y")==y]
b[[a]] <- brick(raster(vi),nl=as.integer(sum(format(timeInfo$outputLayerDates,"%Y")==y)), values=FALSE)
b[[a]] <- writeStart(b[[a]], filename=oname, datatype=opts$datatype, overwrite=opts$overwrite)
names(b[[a]]) <- names
}
} else if(outputAs=="one")
{
y <- unique(format(timeInfo$outputLayerDates,"%Y"))
oname <- paste0(outDirPath,prefixSuffix[1],"_from",paste0(y,collapse="to"),".",nameL,inlam,".",prefixSuffix[2],dataFormat)
b <- list()
b[[1]] <- brick(raster(vi),nl=length(fitt), values=FALSE)
b[[1]] <- writeStart(b[[1]], filename=oname, datatype=opts$datatype, overwrite=opts$overwrite)
names(b[[1]]) <- timeInfo$outputLayerDates
} else if (outputAs=="single")
{
d <- sort(format(timeInfo$outputLayerDates,"%Y%j"))
oname <- paste0(outDirPath,prefixSuffix[1],".",d,".",nameL,inlam,".",prefixSuffix[2],dataFormat)
b <- vector(mode="list",length=length(oname))
for(a in seq_along(oname))
{
b[[a]] <- raster(vi)
b[[a]] <- writeStart(b[[a]], filename=oname[a], datatype=opts$datatype, overwrite=opts$overwrite)
}
}
tr <- blockSize(vi)
cluster <- raster:::.doCluster()
if (cluster)
{
# beginCluster()
cl <- getCluster()
on.exit(endCluster())
nodes <- getOption("rasterClusterCores")
# MODIS fails to load if not done like that ...
clF <- function(i){require(MODIS)}
for (i in 1:nodes)
{
sendCall(cl[[i]], clF, i, tag=i)
recvOneData(cl)
}
# better to be save than sorry:
clusterEvalQ(cl,require(bitops))
clusterEvalQ(cl,require(rgdal))
clusterEvalQ(cl,require(ptw))
tr <- blockSizeCluster(vi)
}
cat("Data is in, start processing!\n")
if(tolower(mergeDoyFun)=="max")
{
mergeFun <- unifyDoubleMX
} else if (tolower(mergeDoyFun)=="weighted.mean" | tolower(mergeDoyFun)== "mean")
{
mergeFun <- unifyDoubleWM
}
clFun <- function(l)
{
val <- getValues(vi, row=tr$row[l], nrows=tr$nrows[l])
val <- t(val)
mtrdim <- dim(val)
set0 <- matrix(FALSE,nrow=mtrdim[1], ncol=mtrdim[2])
set0[is.na(val)] <- TRUE
if (!is.null(w))
{
wtu <- getValues(w, row=tr$row[l], nrows=tr$nrows[l])
# is it not a weight info [0-1]?
if(max(wtu,na.rm=TRUE) > 1)
{
if(is.null(bitShift) | is.null(bitMask))
{
# try to detect VI usefulness layer
bits <- detectBitInfo(vi,"VI usefulness",warn=FALSE)
bitShift <- bits$bitShift
bitMask <- bits$bitMask
}
if(is.null(bitShift) | is.null(bitMask))
{
stop("Could not extract 'bits' for weighting from this product. Use '?makeWeights' function to generate weights manually!")
}
wtu <- makeWeights(wtu, bitShift = bitShift, bitMask = bitMask, threshold = threshold, decodeOnly = FALSE)
}
wtu <- t(wtu)
set0[wtu==0] <- TRUE
} else
{
# if no weighting info is available then weight = 1
wtu <- matrix(1,nrow=mtrdim[1],ncol=mtrdim[2])
}
if (inherits(t,"Raster"))
{
inTu <- getValues(t, row=tr$row[l], nrows=tr$nrows[l])
inTu <- t(inTu)
set0[is.na(inTu)] <- TRUE
set0[ inTu <= 0 ] <- TRUE
t0 <- min(timeInfo$inDoys[1]) - 1
if(!collapse)
{
inTu <- t(repDoy(t(inTu),layerDate=timeInfo,bias=-t0))
}
inTu[set0] <- 0
} else
{
if (collapse)
{
inTu <- matrix(timeInfo$inDoys,nrow=length(timeInfo$inDoys),ncol=mtrdim[2])
} else
{
inTu <- matrix(timeInfo$inSeq,nrow=length(timeInfo$inSeq),ncol=mtrdim[2])
}
}
# the entire info to use or not a pix is in "wtu"
wtu[set0] <- 0
val[set0] <- 0
out <- matrix(NA, nrow=length(fitt), ncol=mtrdim[2])
if(!is.null(outlierThreshold))
{
kickOutlier <- function(vals,weights,lambda,threshold)
{
fTS <- whit2(vals,w=weights,lambda=lambda)
weights[weights==1][abs(vals[weights==1]-fTS[weights==1]) > threshold] <- 0
return(weights)
}
} else
{
# if is.null(outlierThreshold) generate a fake function to avoid a per pixel "if"
kickOutlier <- function(vals,weights,lambda,threshold)
{
return(weights)
}
}
if (collapse)
{
vec0 <- rep(0,365 + (2*timeInfo$call$pillow) + 30) # add a save length of data (because layer doy + effectice composite doy)
} else
{
vec0 <- rep(0,max(timeInfo$inSeq,timeInfo$outSeq) - min(timeInfo$inSeq,timeInfo$outSeq) - 1 + 30)
}
# minimum "minDat" input values for filtering
Cvec <- (colSums(wtu > 0) >= minDat)
Cvec <- (1:mtrdim[2])[Cvec]
ind <- inTu > 0
win <- options("warn")
options(warn=-1)
for (u in Cvec)
{
index <- ind[,u]
use <- mergeFun(vx=val[index,u],wx=wtu[index,u],tx=inTu[index,u])
valVec <- wtVec <- vec0
if(!collapse)
{
valVec[use$tx] <- use$vx
wtVec[use$tx] <- use$wx
} else
{
newOrder <- doCollapse(tx=use$tx,pillow=timeInfo$call$pillow)
valVec[newOrder$sequence] <- use$vx[newOrder$order]
wtVec[newOrder$sequence] <- use$wx[newOrder$order]
}
wtVec <- kickOutlier(vals=valVec,weights=wtVec,lambda=lambda,threshold=outlierThreshold)
#plot(valVec,ylim=c(-1000,9000))
for(i in 1:nIter)
{
fTS <- whit2(valVec,w=wtVec,lambda=lambda)
valVec[valVec < fTS] <- fTS[valVec < fTS]
}
out[,u] <- fTS[fitt]
#lines(fTS,col=2)
}
options(warn=win$warn)
out[,colSums(abs(out))==0] <- NA
return(t(out))
}
if (!cluster)
{
for (i in seq_along(tr$row))
{
res <- clFun(i)
if(doround)
{
res <- round(res)
}
b <- writeValuesMODIS(b,res,tr$row[i],timeInfo,collapse,outputAs)
}
} else
{
for (i in 1:nodes)
{
sendCall(cl[[i]], fun = clFun, args = i, tag=i)
}
for (i in 1:tr$n)
{
d <- recvOneData(cl)
if (!d$value$success)
{
stop("cluster error")
}
ni <- nodes + i
if (ni <= tr$n)
{
sendCall(cl[[d$node]], fun = clFun, args = ni, tag=ni)
}
if(doround)
{
d$value$value <- round(d$value$value)
}
#####
b <- writeValuesMODIS(b,res,tr$row[d$value$tag],timeInfo,collapse,outputAs)
}
}
writeStopMODIS(b,timeInfo,outputAs,collapse)
return(b)
}
#####################################################
writeValuesMODIS <- function(b,val,row,timeInfo,collapse,outputAs)
{
if(collapse)
{
d <- seq_along(sprintf("%03d",seq(as.numeric(format(min(timeInfo$outputLayerDates),"%j")),as.numeric(format(max(timeInfo$outputLayerDates),"%j")),by=timeInfo$call$nDays)))
} else
{
d <- seq_along(b)
}
if(outputAs=="single")
{
for (a in seq_along(d))
{
b[[a]] <- writeValues(b[[a]], val[,a], row)
}
} else if(outputAs=="one")
{
b[[1]] <- writeValues(b[[1]], val, row)
} else
{
for (a in seq_along(d))
{
y <- unique(format(timeInfo$outputLayerDates,"%Y"))[a]
b[[a]] <- writeValues(b[[a]], val[,format(timeInfo$outputLayerDates,"%Y")==y], row)
}
}
return(b)
}
writeStopMODIS <- function(b,timeInfo,outputAs,collapse)
{
for (a in seq_along(b))
{
b[[a]] <- writeStop(b[[a]])
nam <- filename(b[[a]])
extension(nam) <- ""
if (collapse & outputAs!="single")
{
write.table(x=unique(format(timeInfo$outputLayerDates,"%j")),
file=nam, col.names=FALSE, row.names=FALSE)
} else if(outputAs=="one")
{
write.table(x=timeInfo$outputLayerDates, file=nam,
col.names=FALSE, row.names=FALSE)
} else if (outputAs=="yearly")
{
y <- unique(format(timeInfo$outputLayerDates,"%Y"))
for (ax in seq_along(y))
{
ind <- format(timeInfo$outputLayerDates,"%Y")==y[ax]
write.table(x=timeInfo$outputLayerDates[ind], file=nam,
col.names=FALSE, row.names=FALSE)
}
}
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.