R/Habitat.R

In SGronsdahl/GIFT: Geomorphic instream flow tool for conducting hydraulic-habitat modelling.

### Geomorphic Instream Flow Tool ###

#' A function to execute the habitat model
#'
#' This function executes the hydraulic geometry simulator to evaluate reach-averaged depths and velocities
#' generated at flows less than bankfull conditions.
#' For more information about this model see: McParland et al. (2016) and Gronsdahl et al. (XXXX)

#' @param hydraulics dataframe output generated by the AvgHydraulics function
#' @param d_curve dataframe of depth suitabilities with columns: 'depth' (m) and 'suit' (dimensionless suitability from 0-1)
#' @param v_curve dataframe of velocity suitabilities with columns: 'velocity' (m) and 'suit' (dimensionless suitability from 0-1)
#' @param s_curve dataframe of substrate suitabilities with columns: 'lower' (mm) and 'upper' (mm) of grain size classes and 'suit' (dimensionless suitability from 0-1)
#' @param gsd vector of grain size distribution (mm)
#' @param wua_output defaults to TRUE. Expression to specify a .jpeg output of streamflow-WUA relationship
#' @export
#' @return data frame of streamflow, suitabilities, channel width, and WUA
#' Habitat()

Habitat = function(hydraulics,
                          d_curve, v_curve, s_curve = NULL,
                          gsd = NULL, wua_output = TRUE) {

  # load libraries
  library(birk)
  library(dplyr)
  library(zoo)

  ##################################################################
  #input validation
  if(!is.data.frame(hydraulics)) {stop("The 'hydraulics' input parameter should be a data frame")}
  if(nrow(hydraulics) < 3) {warning("The 'hydraulics' input parameter is a data frame with less than 3 rows")}
  if(!"Q" %in% colnames(hydraulics)) {stop("The 'hydraulics' input parameter should be a data frame with a column named 'Q'")}
  if(!"Ai" %in% colnames(hydraulics)) {stop("The 'hydraulics' input parameter should be a data frame with a column named 'Ai'")}
  if(!"Wi" %in% colnames(hydraulics)) {stop("The 'hydraulics' input parameter should be a data frame with a column named 'Wi'")}
  if(!"di" %in% colnames(hydraulics)) {stop("The 'hydraulics' input parameter should be a data frame with a column named 'di'")}
  if(!"Ui" %in% colnames(hydraulics)) {stop("The 'hydraulics' input parameter should be a data frame with a column named 'Ui'")}
  if(!is.numeric(hydraulics$Q) | !is.numeric(hydraulics$Ai) | !is.numeric(hydraulics$Wi) |
     !is.numeric(hydraulics$di) | !is.numeric(hydraulics$Ui)) {stop('A field in the hydraulics data frame is not numeric')}

  if(!is.data.frame(d_curve)) {stop("The 'd_curve' input parameter should be a data frame")}
  if(!"depth" %in% colnames(d_curve)) {stop("The 'd_curve' input parameter should be a data frame with a column named 'depth'")}
  if(!"suit" %in% colnames(d_curve)) {stop("The 'd_curve' input parameter should be a data frame with a column named 'suit")}
  if(!is.numeric(d_curve$depth)) {stop("The depths specified in d_curve are not numeric")}
  if(!is.numeric(d_curve$suit)) {stop("The suitabilities specified in d_curve are not numeric")}
  if(min(d_curve$depth) < 0) {stop("The 'd_curve' data frame specifies depths < 0")}
  if(min(d_curve$suit) < 0) {stop("The 'd_curve' data frame specifies suitabilities < 0")}
  if(max(d_curve$suit) > 1) {stop("The 'd_curve' data frame specifies suitabilities > 1")}

  if(!is.data.frame(v_curve)) {stop("The 'v_curve' input parameter should be a data frame")}
  if(!"velocity" %in% colnames(v_curve)) {stop("The 'v_curve' input parameter should be a data frame with a column named 'velocity'")}
  if(!"suit" %in% colnames(v_curve)) {stop("The 'v_curve' input parameter should be a data frame with a column named 'suit")}
  if(!is.numeric(v_curve$velocity)) {stop("The velocities specified in v_curve are not numeric")}
  if(!is.numeric(v_curve$suit)) {stop("The suitabilities specified in v_curve are not numeric")}
  if(min(v_curve$velocity) < 0) {stop("The 'v_curve' data frame specifies velocities < 0")}
  if(min(v_curve$suit) < 0) {stop("The 'v_curve' data frame specifies suitabilities < 0")}
  if(max(v_curve$suit) > 1) {stop("The 'v_curve' data frame specifies suitabilities > 1")}

  if(is.null(s_curve) == TRUE) {
  } else {
    if(!is.data.frame(s_curve)) {stop("If specified, the 's_curve' parameter should be a data frame.")}
    if(!"upper" %in% colnames(s_curve)) {stop("The 's_curve' input parameter should be a data frame with a column named 'upper'")}
    if(!"lower" %in% colnames(s_curve)) {stop("The 's_curve' input parameter should be a data frame with a column named 'lower'")}
    if(!"suit" %in% colnames(s_curve)) {stop("The 's_curve' input parameter should be a data frame with a column named 'suit")}
    if(!is.numeric(s_curve$upper) | !is.numeric(s_curve$lower) | !is.numeric(s_curve$suit)) {stop("One of 'upper', 'lower', or 'suit' specified in 's_curve' is not numeric")}
    if(min(s_curve$lower) < 0 | min(s_curve$upper) < 0) {stop("The 's_curve' data frame specifies grain sizes < 0")}
    if(min(s_curve$suit) < 0) {stop("The 's_curve' data frame specifies suitabilities < 0")}
    if(max(s_curve$suit) > 1) {stop("The 's_curve' data frame specifies suitabilities > 1")}
    if(max(s_curve$upper) < 1 | max(s_curve$lower) < 1) {warning("Check the 's_curve' 'upper' and 'lower' grain sizes are in mm.")}
  }

  if(is.null(gsd) == TRUE) {
  } else {
      if(!is.vector(gsd)){stop("If specified, 'gsd' should be a vector")}
      if(!is.numeric(gsd)){stop("If specified, 'gsd' should be a numeric")}
      if(min(gsd) < 0 ){stop(" 'gsd' includes grain sizes < 0 mm ")}
      if(max(gsd) < 2 ){warning("Ensure 'gsd' are specified in mm ")}
      if(length(gsd) < 10 ){warning(" 'gsd' includes < 10 observations ")}
    }
  ################################################################################

  # input hydraulics
  mod_hyd = hydraulics

  # Substrate suitability criteria
  #########################################################
  # if-else function to evaluate a reach-averaged substrate suitability
  # value based on a grain-size distribution specified in the model.
  # if no suitability curve or GSD are provided the substrate suitability
  # defaults to 1.0

  if(is.null(s_curve) | is.null(gsd)) {
    sub.suit = 1
  } else {
    n.class = nrow(s_curve)

    sub.out = data.frame()
    for( i in 1 : n.class) {
      #i = 1
      class.obvs = length(which(gsd >= s_curve$lower[i] &
                                  gsd < s_curve$upper[i]))
      to.add = c(i, class.obvs, s_curve$suit[i])
      sub.out = rbind(sub.out, to.add)
    }

    sub.suit = sum((sub.out[ , 2] * sub.out[ , 3])) /
      sum(sub.out[ , 2])
  }


  ###########################################
  # Estimate statistical distributions

  # set up as per original code put together by McParland et al. (2016)
  bins = seq(0, 6, 0.05) # bins that will be used for depth and velocity distributions

  und = 1 # mean of normal depth distribution
  sdnd = 0.52 # standard deviation of normal depth distribution
  nd = dnorm(bins, und, sdnd) # normal depth density distribution

  ulnd = 0 # mean of lognormal depth distribution
  sdlnd = 1.09 # standard deviation of lognormal depth distribution
  lnd = dlnorm(bins, ulnd, sdlnd) # lognormal depth density distribution

  ########################################################
  # convert depth, velocity, and substrate suitabilities into a WUA value

  # generate objects for use in the loop
  n.mod = nrow(mod_hyd)
  WUA.out = data.frame()

  for (i in 1:n.mod){
    #i = 9

    ### Depths ###
    # Calculate the Froude number and mixing factor for the given flow
    Fr = mod_hyd$Ui[i] / sqrt(9.81 * mod_hyd$di[i])

    # Calculate mixing parameter (6a of Schweizer)
    Smix = (exp(-4.72 - 2.84 * (log(Fr)))) / (1 + (exp(-4.72 - 2.84 * (log(Fr)))))

    # relative depth distribution
    depth.dist = ((1 - Smix) * nd + Smix * lnd)

    # overwirte any negative values to zero
    depth.dist = ifelse(depth.dist < 0, 0, depth.dist)

    # absolute depths
    abs.depths = mod_hyd$di[i] * bins

    ### Velocities ###
    # calculate scale parameter for velocity
    s = -0.150 - (0.252 * log(Fr))

    # calculate velocity distributions (Saraeva and Hardy)
    vel.dist = s*((3.33 * exp(-bins / 0.693)) +
                    (0.117 * exp( - ((bins - 8) / 1.73)^2))) +
      ((1 - s) * (0.653 * exp(- ((bins - 1) / 0.664)^2 )))

    # overwirte any negative values to zero
    vel.dist = ifelse(vel.dist < 0, 0, vel.dist)

    # absolute velocities
    abs.velocities = mod_hyd$Ui[i] * bins

    ######################################
    # apply suitability curves to outputs
    n.bins <- length(abs.depths)

    # find the indexed depth/velocity suitability that nearest matches the
    # simulated bins
    sum.d.suit = 0
    sum.v.suit = 0

    for (j in 1 : n.bins){

      # identify index of nearest suitability value
      j.depth <- which.closest(d_curve$depth, abs.depths[j])
      j.vel <- which.closest(v_curve$velocity, abs.velocities[j])

      sum.d.suit = sum.d.suit + depth.dist[j]*d_curve$suit[j.depth]
      sum.v.suit = sum.v.suit + vel.dist[j]*v_curve$suit[j.vel]
    }

    d.suit = sum.d.suit/sum(depth.dist)
    v.suit = sum.v.suit/sum(vel.dist)

    # calculate the composite suitability and output WUA using if-else
    # statement to determine if substrate suitability is applied

    WUA.out[i, 1] <- mod_hyd$Q[i]
    WUA.out[i, 2] <- d.suit
    WUA.out[i, 3] <- v.suit
    WUA.out[i, 4] <- sub.suit
    WUA.out[i, 5] <- mod_hyd$Wi[i]
    WUA.out[i , 6] <- WUA.out[i, 2] * WUA.out[i, 3] * WUA.out[i, 4] * WUA.out[i, 5]
  }

  # name columns
  colnames(WUA.out) = c("Q", "d.suit", "v.suit", "s.suit", "w", "WUA")

  # output WUA curve and figure
  if(wua_output == TRUE){
    ####################################################
    # output WUA figure

    # apply smoothing filter
    plot_y = rollmean(WUA.out$WUA, k = 5, na.pad = TRUE)

    jpeg("WUA_Q.jpeg", width = 7, height = 5, units = "in", res = 300)
    par(mar = c(4.5, 4.5, 1, 1))
    plot(WUA.out$Q, plot_y, type = "l", lwd = 2,
         xlab = expression("Discharge ("*m^3*s^{-1}*")"),
         ylab = expression("WUA ("*m^{2}*m^{-1}*")"),
         xaxs = "i", yaxs = "i",
         ylim = c(0, (1.1 * max(plot_y, na.rm = T))))
    grid()
    box()
    dev.off()

  } else {
  }

  return(WUA.out)
}