R/mldthermocline.R

In kevans27/aqeco: Basic Functions for Aquatic Ecologists

mldthermocline = function(df, depthname, tempname){

  #Collect depths from the dataframe
  depths <- unlist(df[depthname], use.names = FALSE)
  temps <- unlist(df[tempname], use.names = FALSE)

  #Make a list of the slopes between each consecutive pair of points
  slopes <- c()
  for(i in 1:(length(depths)-1)){
    xs <- c(depths[i], depths[i+1])
    ys <- c(temps[i], temps[i+1])
    slope <- coef(lm(ys~xs))[2]
    slopes <- c(slopes, slope)
  }

  #Absolute value these slopes to allow for
  #inverted profiles (warmer depths and cooler shallows)
  slopes <- abs(unname(slopes))

  #Had to give a 0.5% range in case of slight unequalties
  #Find the indices of the depths at which the slow is at the max value
  inds <- which(slopes < (max(slopes)*1.005) & slopes > (max(slopes)*0.995))

  #Return results including printing a message of the number
  #of thermocline depths and the depth right before the thermocline hits.
  if (length(inds) == 1){
    cat('You have a clear thermocline. I\'ll return a single value for depth.', '\n')
    if (temps[1] < temps[length(temps)]){
      cat('Check your profile. Looks inverted.', '\n')
    }
    return(depths[inds])
  }
  if (length(inds) > 1){
    cat('You have multiple potential thermoclines. I am returning a list of the possible depths.', '\n')
    if (temps[1] < temps[length(temps)]){
      cat('Check your profile. Looks inverted.', '\n')
    }
    return(depths[inds])
  } else{
    cat('Something went wrong. Check your profile.', '\n')
    return()
  }
}