R/Habitat.R

In SGronsdahl/Geomorphic-Approach:

### Geomorphic Instream Flow Tool ###

#' A function to execute the ASHGS 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 max_Q maximum discharge (m3/s) to simualted WUA for.  Defaults to 1 m3/s.
#' @param d_curve dataframe of depth suitabilities with two columns: 'depth' (m) and 'suit' (dimensionless suitability from 0-1)
#' @param v_curve dataframe of velocity suitabilities with two columns: 'velocity' (m) and 'suit' (dimensionless suitability from 0-1)
#' @param s_curve dataframe of substrate suitabilities with three 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, max_Q = 10, 
                          d_curve, v_curve, s_curve = NULL,
                          gsd = NULL, wua_output = TRUE) {
  
  # load libraries
  library(birk)
  library(dplyr)
  library(zoo)
  
  # input hydraulics
  mod_hyd = hydraulics %>% filter(Q <= max_Q) 
  
  # 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
    
    ######################################
    # attach suitability data to outputs
    d.bin.suits = data.frame()
    v.bin.suits = data.frame()
    n.bins = length(abs.depths)
    
    # find the indexed depth/velocity suitability that nearest matches the
    # simulated bins
    for (j in 1 : n.bins){
      
      # j = 31
      
      #depth_eval = approx(d_curve$depth, d_curve$suit, xout = abs.depths[j])
      #vel_eval = approx(v_curve$velocity, v_curve$suit, xout = abs.velocities[j])
      
      j.depth = which.closest(d_curve$depth, abs.depths[j])
      j.vel = which.closest(v_curve$velocity, abs.velocities[j])
      
      # output suitability values into data frames
      d.bin.suits[j, 1] = depth.dist[j]
      d.bin.suits[j, 2] = abs.depths[j]
      d.bin.suits[j, 3] = d_curve$suit[j.depth]
      v.bin.suits[j, 1] = vel.dist[j]
      v.bin.suits[j, 2] = abs.velocities[j]
      v.bin.suits[j, 3] = v_curve$suit[j.vel]
    }
    
    colnames(d.bin.suits) = c("density", "depth", "suitability")
    colnames(v.bin.suits) = c("density", "velocity", "suitability")
    
    # aggregate the depth and velocity suitabilities for the simulated flow
    d.suit = sum(d.bin.suits$density * d.bin.suits$suitability) /
      sum(d.bin.suits$density)
    v.suit = sum(v.bin.suits$density * v.bin.suits$suitability) /
      sum(v.bin.suits$density)
    
    # 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)
}