R/AASHGS_Execute.R

In SGronsdahl/AASHGS: At-a-station hydraulic geometry simulator

#' A function to execute the ASHGS model
#'
#' This function executes the at-a-station hydraulic geometry simulator, which is a method of determining weighted usable area.
#' For more information about the theoretical component of this model see:
#' McParland et al. (2016). At-a-station hydrualic geometry simulator. River Res. Applic. 32: 399 - 410.
#' This function requires that the 'birk' package is loaded in your directory.
#'
#'
#' @param S channel gradient (m/m)
#' @param wb reach averaged bankfull width (m).
#' @param db reach averaged bankfull depth (m).
#' @param db_max Defaults to NULL. reach averaged maximum bankfull depth (m).  Specifying db_max is preferred.
#' @param max_Q maximum discharge (m3/s) to simualted WUA for.  Defaults to simulated bankfull dsicharge.
#' @param d_curve A dataframe of depth suitability curves with four columns as follows:  column 1: "species"  - This identified species (eg "species"); column 2: "life_stage" - This identifies life stage (eg "fry"); column 3: "depth" (m); column 4: "suit" (0 - 1)
#' @param v_curve A dataframe of velocity suitability curves with four columns as follows:  column 1: "species"  - This identified species (eg "species"); column 2: "life_stage" - This identifies life stage (eg "fry"); column 3: "velocity" (m/s); column 4: "suit" (0 - 1)
#' @param s_curve Defaults to NULL> A dataframe of substrate suitability curves with five columns as follows:  column 1: "species"  - This identified species (eg "species"); column 2: "life_stage" - This identifies life stage (eg "fry"); column 3: "lower" (mm) - This identifies the lower limit of each class; column 4: "upper" (mm) - This identifies the upper limit of each class; column 5: "suit" (0 - 1)
#' @param D84 grain size (mm)
#' @param gsd Defaults to NULL. A dataframe of the grain size distribution with one column: individual grain sizes (mm). Intended for use with output of the Find_D84 function.
#' @export
#' @return .csv of hydraulic data
#' @return .csv of WUA suitability data 'suitability.csv'
#' @return .pdf figure of simulated cross section
#' @return .pdf figure of the WUA - Q curve
#' AASHGS_Execute()

AASHGS_Execute = function(S, wb, db, db_max, max_Q, d_curve, v_curve, s_curve = NULL,
                          D84, gsd = NULL) {

  ###########################################################

  # load libraries
  library(birk)
  library(dplyr)

  ##### Define Find_U Function #####

  findU = function(Wb, S, D84, depths) {

    deltaX = 0.0001
    Xgrid = Wb * seq(0, 1, deltaX)

    wet.vert = depths[depths >= 0]
    Wi = length(wet.vert) * deltaX * Wb
    Ai = sum(wet.vert * deltaX * Wb)
    di = Ai / Wi
    Pi = sum((diff(wet.vert)^2 + (max(Xgrid) * deltaX) ^2) ^ (1/2))
    Ri = Ai/Pi

    # Ferguson's continuously varying power law
    D.84 = D84 / 1000
    g = 9.81
    a1 = 6.5
    a2 = 2.5
    Res = a1 * a2 * (Ri / D.84) /
      (a1^2 + a2^2 * (Ri / D.84) ^ (5/3)) ^ (1/2)
    Ui = Res * sqrt(g * Ri * S) # Velocity (m/s)

    # Formatting the outputs in a dataframe
    df = data.frame(Ai, Wi, di, Ui)

    return(df)
  }

  #############################################################
  ##### Simulate Hydraulics #####

  # define grid
  deltaX = 0.0001 # Resolution of grid upon which to calculate Q
  # (as a proportion of wb)
  deltaY = 0.001  # increment by which to change depths when estimating HG

  # Ferguson model's shape factor b
  if (is.na(db_max) == TRUE){
    b = (wb / db) / 100
  } else {
    b = 1 - (db / db_max) # preferred method
  }

  # Generate the cross section
  dmax = (1 + b) / (1 - b) * db
  X = c(0, b * wb, 0.99 * wb, wb)
  Y = 5 * db- c(0, db, dmax, 0)

  # Interpolate the distribution onto an XS raster
  Xgrid = wb * seq(0, 1, deltaX)
  Ygrid = matrix(unlist(approx(X, Y, Xgrid)), ncol = 10001, byrow = TRUE)[2,]

  # Specify the values of the water surface elevation for which to calculate Wi
  Zw = 5 * db - dmax + seq(0.1 * dmax, dmax, deltaY * dmax)

  # For loop to calculate the width and discharge for each chosen water level
  output = data.frame()

  results = list()
  for (j in 1:length(Zw)) {
    #j = 3
    depths = Zw[j] - Ygrid # Calculate the depths, for each vertical
    results = findU(wb, S, D84, depths)
    results = c(Q = results[1, 4] * results[1, 1], results)
    output = rbind(output, results)
  }
  write.csv(output, "simulated_hydraulics.csv", row.names = FALSE)

  ###### Prepare graph of cross section #####
  # select y-values based on db_max input
  if(is.na(db_max) == TRUE){
    plot_y = c(0, (db * -1), (db / b * - 1), 0)
  } else {
    plot_y = c(0, (db * -1), (db_max * - 1), 0)
  }

  # set up x-values to plot
  plot_x = c(0, (b * wb), (0.99 * wb), wb )

  # output a plot of the cross section
  pdf("cross_section.pdf", width = 6, height = 4)
  par(mar = c(4.5, 4.5, 1, 1))
  plot(plot_x, plot_y, type = "l",
       xlab = "Width (m)",
       ylab = "Depth (m)",
       ylim = c((min(plot_y) * 1.2), 0)
  )
  abline(a = (db * -1), 0, lty = 2, col = "grey")
  legend("bottomleft", col = c("black", "grey"), bty = "n",
         lty = c(1, 2), legend = c("Channel cross section",
                                   "Average depth"))
  dev.off()



  ########### Distribution Settings ###########

  # set up as per original code put together by McParland et al. (2016)

  bins = seq(0, 3, 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

  unv = 1 # mean of normal velocity distribution
  sdnv = 0.52 # standard deviation of normal velocity distribution
  nv = dnorm(bins, unv, sdnv) # normal velocity denisty distribution

  ulnv = 0 # mean of lognormal velocity distribution
  sdlnv = 1.19 # standard deviation of lognormal velocity distribution
  lnv = dlnorm(bins, ulnv, sdlnv) # lognormal velocity density distribution

  # determine substrate suitability criteria
  ###########

  # if-else criteria based on sub_suit definition at start
  if(missing(s_curve) | missing(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])
  }

  # read in files for use in the for-loop to simulate WUA
  simulated = read.csv("simulated_hydraulics.csv")
  simulated = subset(simulated, Q <= max_Q)

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

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

    # Calculate the Froude number and mixing factor for the given flow
    Fr = simulated$Ui[i] / sqrt(9.81 * simulated$di[i])
    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)
    # absolute depths
    abs.depths = simulated$di[i] * bins

    # relative velocity distribution
    vel.dist = ((1 - Smix) * nd + Smix * lnd)
    # absolute velocities
    abs.velocities = simulated$Ui[i] * bins

    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 = 61
      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
    #if(sub_suit == "yes"){
    WUA.out[i, 1] = simulated$Q[i]
    WUA.out[i, 2] = d.suit
    WUA.out[i, 3] = v.suit
    WUA.out[i, 4] = sub.suit
    WUA.out[i, 5] = simulated$Wi[i]
    WUA.out[ , 6] = WUA.out[, 2] * WUA.out[, 3] * WUA.out[, 4] * WUA.out[, 5]
    colnames(WUA.out) = c("Q", "d.suit", "v.suit", "s.suit", "w", "WUA")
  }
  write.csv(WUA.out, "wua_q.csv", row.names = FALSE)

  WUA.out$MAD = WUA.out[ , 1] / max(WUA.out[ , 1]) * 100
  pdf("WUA_Q.pdf", width = 7, height = 5)
  par(mar = c(4.5, 4.5, 1, 1))
  plot(WUA.out$MAD, WUA.out$WUA, type = "l",
       xlab = "Discharge (% of MAD)",
       ylab = expression("WUA ("*m^{2}*"/m)"),
       xaxs = "i", yaxs = "i",
       ylim = c(0, (1.1 * max(WUA.out$WUA)))
  )
  grid()

  dev.off()

  ## End of Function ##
}