#' Define stem-cyclic phases
#'
#' @description The function identifies and assigns each timestamp to one of the three distinct phases of contraction, expansion and stem-radius increment (Deslauriers et al. 2011) for dendrometer series from a \code{data.frame} with gap-free dendrometer data.
#'
#' @usage phase_def(dm.gpf, resolution = dendro.resolution(dm.gpf),
#' shapeSensitivity = 0.6, minmaxDist = 0.2, minmaxSD = 2,
#' radialIncrease = "max")
#'
#' @param dm.gpf a \code{data.frame} with either gap-free or gap-filled dendrometer series as produced by \code{\link{fill_gaps}}.
#' @param resolution a \code{numeric} specifying the resolution of the dendrometer data in seconds. Defaults to the resolution of \code{dm.gpf} as calculated using \code{\link{dendro.resolution}}.
#' @param shapeSensitivity a \code{numeric} specifying a time window, defined as proportion of a single day. Within this time window possible extrema points (i.e. minimum and maximum) in dendrometer measurements are searched for. Defaults to 0.6 (60\% of a day). See details for further explanation.
#' @param minmaxDist a \code{numeric} specifying the minimum temporal distance between consecutive minimum and maximum points (i.e. in the x direction). Defaults to 0.2 (20\% of a day).
#' @param minmaxSD a \code{numeric} specifying the minimum difference between consecutive minimum and maximum points expressed as a number of standard deviations (i.e. in the y direction). Defaults to 2.
#' @param radialIncrease a \code{character} string of \code{"max", "min", "mid"}, specifying when the stem-radius increment phase should start, with \code{"max"} as the most, and \code{"min"} as the least conservative approach; \code{"mid"} is in between. See details for further explanation.
#'
#' @details The function defines the stem-cyclic phases of contraction, expansion, and stem-radius increment, as described in Deslauriers et al. (2011). The function is a more robust version of the original SAS routine, as its architecture allows to handle noisy data as well.
#'
#' First, the function searches for minimum and maximum points within a daily time window as specified by \code{shapeSensitivity}. Then, the original dendrometer series are offset by \code{(1 - shapeSensitivity) / 2} in both directions to assure whether the identified extrema are indeed the extrema of cyclic phases. A comparison between the original and offset series allows to select all appropriate minimum and maximum values.
#'
#' The arguments \code{minmaxDist} and \code{minmaxSD} specify the temporal distance and the minimum difference between consecutive minimum and maximum points (i.e. in x and y direction), respectively. The argument \code{radialIncrease} determines from which moment on data points should be assigned to the stem-radius increment phase: when points are continuously above the previous maximum (\code{"max"}), when a single data point is above the previous maximum (\code{"min"}), or right in between \code{"min"} and \code{"max"} (\code{"mid"}).
#'
#' @return
#' The function returns a \code{data.frame} with numbers indicating the different stem-cyclic phases: (1) contraction, (2) expansion, (3) stem-radius increment for each timestamp.
#'
#' @author Marko Smiljanic
#'
#' @references Deslauriers, A., Rossi, S., Turcotte, A., Morin, H. and Krause, C. (2011) A three-step procedure in SAS to analyze the time series from automatic dendrometers. \emph{Dendrochronologia} 29: 151-161.
#'
#' @examples
#' data(dmCD)
#' dm.phase <- phase_def(dmCD)
#'
#' @export phase_def
#'
phase_def <- function(dm.gpf, resolution = dendro.resolution(dm.gpf), shapeSensitivity = 0.6, minmaxDist = 0.2, minmaxSD = 2, radialIncrease = "max")
{
### Recursive call if there is more than one dendrometer series ###
if(!inherits(dm.gpf, "numeric")) {
dm.gpf <- as.data.frame(dm.gpf)
output <- as.data.frame(matrix(nrow = nrow(dm.gpf), ncol = ncol(dm.gpf)))
if(is.dendro(dm.gpf)) {
for(d in 1:ncol(dm.gpf)) {
den <- dm.gpf[[d]]
res <- phase_def(den,resolution,shapeSensitivity, minmaxDist, minmaxSD, radialIncrease)
rowRange <- (nrow(output)-length(res)+1):nrow(output)
output[rowRange, d] <- res
}
names(output) <- names(dm.gpf)
rownames(output) <- rownames(dm.gpf)
}
else {
stop("there is a problem with the input data")
}
return(output)
}
if(any(is.na(dm.gpf))) {
stop("NA values detected in dendrometer series. Run fill_gaps first")
}
### Detect local minima and maxima ###
local_extrema <- function(dendro = parent.frame()$dm.gpf) {
minima <- vector()
maxima <- vector()
i <- 2
while(i < length(dendro)) {
dat <- dendro[i]
datAfter <- dendro[i+1]
datBefore <- dendro[i-1]
if(dat == datAfter) {
dat <- dendro[i] <- dendro[i] + (dat-datBefore)/2
}
if(dat == datBefore) {
iStart <- i
while(dendro[i] == dat) {
i <- i+1
}
pos <- ((iStart-1):i)
vals <- (pos-1)*(dendro[i]-dat)/(i-iStart+1)+dat
}
if(!is.na(datAfter) && !is.na(datBefore)) {
if(dat > datAfter && dat > datBefore) {
maxima <- c(maxima, i)
}
else if(dat < datAfter && dat < datBefore) {
minima <- c(minima, i)
}
}
else {
stop("NA values detected in dendrometer series")
}
i <- i+1
}
result <- list()
result$minima <- minima
result$maxima <- maxima
return(result)
}
### select minima/maxima points for every day depending on the resolution ###
extrema <- function(dendro = parent.frame()$dm.gpf, minima = parent.frame()$minima, maxima = parent.frame()$maxima, shapeSensitivity = parent.frame()$shapeSensitivity) {
dayLen <- 86400/resolution
dayNumb <- round(length(dendro)/dayLen)
dayShape <- vector()
minimaNew <- vector()
maximaNew <- vector()
for(day in 1:dayNumb) {
dayStart <- (day-1)*dayLen+1
dayEnd <- min(dayStart+dayLen, length(dendro))
dendroDay <- dendro[dayStart:dayEnd]
minDendroDay <- min(dendroDay, na.rm=TRUE)
minDendroDay <- which(dendroDay == minDendroDay)
minDendroDay <- minDendroDay[[1]]
maxDendroDay <- max(dendroDay, na.rm=TRUE)
maxDendroDay <- which(dendroDay == maxDendroDay)
maxDendroDay <- maxDendroDay[[1]]
lowSens <- dayLen*shapeSensitivity
highSens <- (1-shapeSensitivity)*dayLen
if((minDendroDay > lowSens) & (minDendroDay < highSens)) {
dayShp <- 0
}
else if((maxDendroDay > lowSens) & (maxDendroDay < highSens)) {
dayShp <- 1
}
else {
dayShp <- NA
}
dayShape <- c(dayShape, dayShp)
if(!is.na(dayShp)) {
dendroMinimaID <- minDendroDay+dayStart-1
dendroMaximaID <- maxDendroDay+dayStart-1
dendroMinimaID <- dendroMinimaID[[1]]
dendroMaximaID <- dendroMaximaID[[1]]
minimaNew <- c(minimaNew, minima[which(minima==dendroMinimaID)])
maximaNew <- c(maximaNew, maxima[which(maxima==dendroMaximaID)])
}
}
result <- list(minimaNew,maximaNew,dayShape)
names(result) <- c("minimaNew", "maximaNew", "dayShape")
return(result)
}
shapeSensitivity = (1-shapeSensitivity)/2
locExt <- local_extrema()
minima <- locExt$minima
maxima <- locExt$maxima
extVals <- extrema()
# offset original series depending on the shapeSensitivity variable and redo the extrema detection
# TODO: Offset might be too big currently. Needs further testing
dayLen <- 86400/resolution
locExtOffset <- local_extrema(dm.gpf[floor(dayLen*shapeSensitivity):length(dm.gpf)])
extValsOffset <- extrema(dm.gpf[floor(dayLen*shapeSensitivity):length(dm.gpf)], locExtOffset$minima, locExtOffset$maxima)
### Depending on the number of indeterminate day curve shapes select either original or offset series
origShape <- extVals$dayShape
offsetShape <- extValsOffset$dayShape
origShapeNAs <- sum(is.na(origShape))
offsetShapeNAs <- sum(is.na(offsetShape))
lenConc <- length(origShape) - length(which(origShape == 0))
lenConv <- length(origShape) - length(which(origShape == 1))
lenConcOff <- length(offsetShape) - length(which(offsetShape == 0))
lenConvOff <- length(offsetShape) - length(which(offsetShape == 1))
diffCC <- abs(lenConc - lenConv)
diffCCOff <- abs(lenConcOff - lenConvOff)
dm.orig <- dm.gpf
if(offsetShapeNAs < origShapeNAs) {
dm.gpf <- dm.gpf[floor(dayLen*shapeSensitivity):length(dm.gpf)]
minima <- locExtOffset$minima
maxima <- locExtOffset$maxima
extVals <- extValsOffset
}
else if(offsetShapeNAs == origShapeNAs) {
if(diffCCOff > diffCC) {
dm.gpf <- dm.gpf[floor(dayLen*shapeSensitivity):length(dm.gpf)]
minima <- locExtOffset$minima
maxima <- locExtOffset$maxima
extVals <- extValsOffset
}
}
if(length(dm.orig) > length(dm.gpf)) {
reps <- length(dm.orig)-length(dm.gpf)
dm.gpf <- c(rep(NA, reps), dm.gpf)
extVals$minimaNew <- extVals$minimaNew+reps
extVals$maximaNew <- extVals$maximaNew+reps
minima <- minima+reps
maxima <- maxima+reps
}
minimaNew <- extVals$minimaNew
maximaNew <- extVals$maximaNew
dayShape <- extVals$dayShape
# par(ask=FALSE) #DEBUG
# plot(dm.gpf, type="l") #DEBUG
### Depending on the ratio of convex to concave day curve shapes maxima or minima will hold true values.
### Fill the other vector from its local minima values
lenConc <- length(dayShape) - length(which(dayShape == 0))
lenConv <- length(dayShape) - length(which(dayShape == 1))
ratioCC <- lenConc/lenConv
if(ratioCC <= 1) {
maximaNew <- vector()
for(mini in 1:(length(minimaNew)-1)) {
minimaStart <- minimaNew[mini]
minimaEnd <- minimaNew[mini+1]
maximas <- maxima[which(maxima > minimaStart & maxima < minimaEnd)]
dendroMax <- dm.gpf[maximas]
maximaNew <- c(maximaNew, maximas[which(dendroMax == max(dendroMax))][1])
}
}
else if(ratioCC > 1) {
minimaNew <- vector()
for(maxi in 1:(length(maximaNew)-1)) {
maximaStart <- maximaNew[maxi]
maximaEnd <- maximaNew[maxi+1]
minimas <- minima[which(minima > maximaStart & minima < maximaEnd)]
dendroMin <- dm.gpf[minimas]
minimaNew <- c(minimaNew, minimas[which(dendroMin == min(dendroMin))][1])
}
}
### if the minima maxima pair is too close to each other delete both
for(i in 1:length(minimaNew)) {
mini <- minimaNew[i]
potMaxiID <- which(abs(maximaNew-mini) < dayLen*minmaxDist)
if(length(potMaxiID) != 0 ) {
potMaxi <- maximaNew[potMaxiID]
if(mini != length(dm.gpf) && potMaxi != length(dm.gpf) && mini != 1 && potMaxi != 1) {
minimaNew[i] <- maximaNew[potMaxiID] <- NA
minimaNew <- na.omit(minimaNew)
maximaNew <- na.omit(maximaNew)
}
}
}
### if the difference between the minima and maxima in minima maxima pair is too small delete both
extremaPoints <- as.data.frame(rbind(cbind(maximaNew,0),cbind(minimaNew, 1)))
names(extremaPoints) <- c("Position", "type")
extremaPoints <- extremaPoints[order(extremaPoints$Position),]
extremaValues <- as.data.frame(cbind(dm.gpf[extremaPoints$Position], extremaPoints$type))
extremaDiff <- abs(diff(extremaValues[,1]))
extremaSD <- (extremaDiff-mean(extremaDiff))/sd(extremaDiff)
extremaPot <- which(extremaSD < (-minmaxSD))
for(ext in extremaPot) {
extPot <- extremaPoints[ext:(ext+1),]
if(diff(extPot$type) == 1) {
for(p in extPot$Position){
maximaNew[which(maximaNew == p)] <- NA
minimaNew[which(minimaNew == p)] <- NA
}
}
}
minimaNew <- na.omit(minimaNew)
maximaNew <- na.omit(maximaNew)
#points(y=dm.gpf[minimaNew],x=minimaNew, col=2) #DEBUG
#points(x=maximaNew, y=dm.gpf[maximaNew], col=3) #DEBUG
### Assign the growth phases ###
minima <- minimaNew
maxima <- maximaNew
minimaLen <- length(minima)
maximaLen <- length(maxima)
minimaFirst <- minima[1] < maxima[1]
maximaLast <- minima[length(minima)] < maxima[length(maxima)]
phase <- vector()
# Depending on the order which comes first in dendrometer series (minima or maxima)
# this if statement goes through extrema points and assigns the phases to vectors
# currently in the if statement one block of code is copied and minima replaced with maxima
# TODO: Avoid this copy with minima or maxima being a reference. If that does not work use a function
if(minimaFirst) {
for(i in 1:length(minima)) {
if(!is.na(maxima[i])) {
phase[(minima[i]+1):maxima[i]] <- 2
if(i > 1) {
ref <- dm.gpf[maxima[i-1]]
dotsPhase3 <- which(dm.gpf[minima[i]:maxima[i]] > ref)
phaseStart <- NULL
if(length(dotsPhase3) > 0 && radialIncrease == "min") {
phaseStart <- min(dotsPhase3)
}
else if(length(dotsPhase3) > 0 && radialIncrease == "max") {
if(length(dotsPhase3) == 1) {
phaseStart <- dotsPhase3
}
else {
while(length(unique(diff(dotsPhase3))) != 1) {
dotsPhase3 <- dotsPhase3[2:length(dotsPhase3)]
}
if(unique(diff(dotsPhase3)) == 1) {
phaseStart <- min(dotsPhase3)
}
else {
phaseStart <- max(dotsPhase3)
}
}
}
else if(length(dotsPhase3) > 0 && radialIncrease == "mid") {
phaseMin <- min(dotsPhase3)
if(length(dotsPhase3) == 1) {
phaseMax <- dotsPhase3
}
else {
while(length(unique(diff(dotsPhase3))) != 1) {
dotsPhase3 <- dotsPhase3[2:length(dotsPhase3)]
}
if(unique(diff(dotsPhase3)) == 1) {
phaseMax <- min(dotsPhase3)
}
else {
phaseMax <- max(dotsPhase3)
}
}
phaseStart <- round(mean(c(phaseMin, phaseMax)))
}
if(!is.null(phaseStart)) {
phaseStart <- phaseStart + minima[i] - 1
phase[phaseStart:maxima[i]] <- 3
}
}
}
if(!is.na(minima[i+1]) && !is.na(maxima[i])) {
phase[(maxima[i]+1):minima[i+1]] <- 1
}
}
}
else {
if(maxima[1] < 1) {
phase[1:maxima[1]] <- NA
}
for(i in 1:length(maxima)) {
if(!is.na(minima[i])) {
phase[(maxima[i]+1):minima[i]] <- 1
}
if(!is.na(maxima[i+1])) {
phase[(minima[i]+1):maxima[i+1]] <- 2
ref <- dm.gpf[maxima[i]]
dotsPhase3 <- which(dm.gpf[minima[i]:maxima[i+1]] > ref)
phaseStart <- NULL
if(length(dotsPhase3) > 0 && radialIncrease == "min") {
phaseStart <- min(dotsPhase3)
}
else if(length(dotsPhase3) > 0 && radialIncrease == "max") {
if(length(dotsPhase3) == 1) {
phaseStart <- dotsPhase3
}
else {
while(length(unique(diff(dotsPhase3))) != 1) {
dotsPhase3 <- dotsPhase3[2:length(dotsPhase3)]
}
if(unique(diff(dotsPhase3)) == 1) {
phaseStart <- min(dotsPhase3)
}
else {
phaseStart <- max(dotsPhase3)
}
}
}
else if(length(dotsPhase3) > 0 && radialIncrease == "mid") {
phaseMin <- min(dotsPhase3)
if(length(dotsPhase3) == 1) {
phaseMax <- dotsPhase3
}
else {
while(length(unique(diff(dotsPhase3))) != 1) {
dotsPhase3 <- dotsPhase3[2:length(dotsPhase3)]
}
if(unique(diff(dotsPhase3)) == 1) {
phaseMax <- min(dotsPhase3)
}
else {
phaseMax <- max(dotsPhase3)
}
}
phaseStart <- round(mean(c(phaseMin, phaseMax)))
}
if(!is.null(phaseStart)) {
phaseStart <- phaseStart + minima[i] - 1
phase[phaseStart:maxima[i+1]] <- 3
}
}
}
}
if(length(phase) != length(dm.gpf)) {
phase[(length(phase)+1):length(dm.gpf)] <- NA
}
#DEBUG: points(x=1:length(phase), y=dendro[1:length(phase)], pch=19, col=phase, cex=1)
return(phase)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.