R/calculateRelativeBedStability.R
In jasonelaw/nrsa: Calculate metrics for NRSA protocol data.
Defines functions
calculateRefinedCriticalDiameter
calculateCrudeCriticalDiameter
calculateCriticalDiameterKaufmann
calculateRelativeBedStability
fillMissingLWD
calculateS_RP100
forceSlopeGt0
fixBoatableDepthUnits
Documented in
calculateCriticalDiameterKaufmann
calculateCrudeCriticalDiameter
calculateRefinedCriticalDiameter
calculateRelativeBedStability
calculateS_RP100
fillMissingLWD
fixBoatableDepthUnits
forceSlopeGt0
#' Calculate several versions of critical diameter: ltest, ldmb_bw4, ldmb_bw5, s_ldmb_bw5 and ldcbf_g08
#'
#' These functions can be used to calculate several versions of the critical particle
#' diameter:
#' \itemize{
#' \item a crude estimate via \code{calculateCrudeCriticalDiameter}
#' \item several more refined ones via \code{calculateRefinedCriticalDiameter}
#' \item a research based one due to Kauffman 2007 via \code{calculateCriticalDiameterKaufmann}
#' }
#'
#' \code{calculateRefinedCriticalDiameter} is used to calculate the metrics
#' ldmb_bw4 (with old.version = T), ldmb_bw5, and s_ldmb_bw5 (when the rp100
#' argument is as returned by \link{calculateS_RP100}).
#' @param xdepth average depth in centimeters
#' @param xslope average slope in percent
#' @param rp100 average residual depth in centimeters
#' @param xbkf_h average bankfull height in meters
#' @param v1w_msq volume lwd in cubic meters / 100 square meeters
#' @param old.version Calculate old version of logged critical diameter that were based on the
#' wood density when it accidentally used the wetted widths instead of the
#' bankfull widths.
#' @param xbkf_w average bankfull width in meters; must be supplied when old.version = T.
#' @param xwidth average width in meters; must be supplied when old.version = T.
#' @param lsub_dmm mean log10(substrate diameter) in millimeters
#' @return the log10(critical diameter)
#' @references
#' Kaufmann, P.R. et al., A roughness-corrected index of relative bed stability
#' for regional stream surveys. Geomorphology (2007)
#' @export
#' @examples
#' library(reshape2)
#'x <- structure(list(UID = c("1", "10", "11", "2", "3", "4", "5", "6", "7", "8", "9"),
#' lsub_dmm = c(1, -1.31451, 1.18991, 1, 1, 1, 0.5408, 3.204, 0.75694, 0.37664, 0.31911),
#' lsub2dmm = c(2, -1.31451, NA, 2, 2, 2, 0.5408, 3.19253, 0.7244, 0.37664, 0.29885),
#' rp100 = c(1, 10.3489, 112.73, 1, 1, 1, 19.0993, 28.3768, 6.4191, 27.2566, 22.2393),
#' sddepth = c(3.95494014, 8.8661, 1.69386, 3.95494014, 3.95494014,
#' 3.95494014, 27.7557, 44.2636, 10.4878, 22.6112, 28.7622),
#' v1w_msq = c(1, NA, 0.005818, NA, NA, 0.3, 4e-05, 0.01461, 0.37977, 0.00015, NA),
#' xbkf_h = c(1, 0.33636, 1.47273, 1, 1, 1, 0.75545, 0.61818, 0.3, 0.4, 0),
#' xbkf_w = c(1, 3.718, 255.273, 1, 1, 1, 21.818, 6.673, 2.927, 3.4818, 8.745),
#' xdepth = c(1, 42.67, 3.29555, 1, 1, 1, 51.16, 71.19, 24.007, 31.033, 66.51),
#' xfc_lwd = c(1, 0.00455, 0.00455, 1, 0, NA, 0, 0.01364, 0.11818, 0, 0),
#' xslope = c(1, 0.775, 0.036, 1, 1, 1, 1.22, 9.75, 8.05, 0, 1.57),
#' xwidth = c(1, 3.085, 163.364, 1, 1, 1, 10.0105, 5.255, 2.445, 2.385, 9.1474),
#' PROTOCOL = c("WADEABLE", "WADEABLE", "BOATABLE", "WADEABLE",
#' "WADEABLE", "WADEABLE", "WADEABLE",
#' "WADEABLE", "WADEABLE", "WADEABLE", "WADEABLE")),
#' .Names = c("UID", "lsub_dmm", "lsub2dmm", "rp100", "sddepth", "v1w_msq", "xbkf_h",
#' "xbkf_w", "xdepth", "xfc_lwd", "xslope", "xwidth", "PROTOCOL"),
#' row.names = c(NA, -11L),
#' class = "data.frame")
#' x$xslope <- forceSlopeGt0(x$xslope)
#' x$xdepth <- fixBoatableDepthUnits(x$xdepth, x$PROTOCOL == 'WADEABLE')
#' x$sddepth <- fixBoatableDepthUnits(x$sddepth, x$PROTOCOL == 'WADEABLE')
#' x$s_rp100 <- calculateS_RP100(x$xslope, x$sddepth)
#' x$v1w_msq <- fillMissingLWD(x$v1w_msq, x$xfc_lwd)
#' x$ltest <- with(x, calculateCrudeCriticalDiameter(xdepth, xslope))
#' x$ldmb_bw5 <- with(x, calculateRefinedCriticalDiameter(xdepth, xslope, rp100, xbkf_h, v1w_msq))
#' x$s_ldmb_bw5 <- with(x, calculateRefinedCriticalDiameter(xdepth, xslope, s_rp100, xbkf_h, v1w_msq))
#' x$ldmb_bw4 <- with(x, calculateRefinedCriticalDiameter(xdepth, xslope, rp100, xbkf_h, v1w_msq, T, xbkf_w, xwidth))
#' x <- cbind(x, with(x, calculateCriticalDiameterKaufmann(xdepth, xslope, rp100, xbkf_h, v1w_msq, lsub_dmm)))
#' x <- cbind(x, with(x, calculateRelativeBedStability(lsub_dmm, lsub2dmm, ltest, ldmb_bw4, ldmb_bw5, s_ldmb_bw5, ldcbf_g08)))
#' metsExpected <-
#' structure(list(UID = c("1", "1", "1", "1", "1", "1", "1", "1",
#' "1", "1", "1", "1", "1", "1", "2", "2", "2", "2", "2", "2", "2",
#' "2", "2", "2", "2", "2", "2", "2", "3", "3", "3", "3", "3", "3",
#' "3", "3", "3", "3", "3", "3", "3", "3", "4", "4", "4", "4", "4",
#' "4", "4", "4", "4", "4", "4", "4", "4", "4", "5", "5", "5", "5",
#' "5", "5", "5", "5", "5", "5", "5", "5", "5", "5", "6", "6", "6",
#' "6", "6", "6", "6", "6", "6", "6", "6", "6", "6", "6", "7", "7",
#' "7", "7", "7", "7", "7", "7", "7", "7", "7", "7", "7", "7", "8",
#' "8", "8", "8", "8", "8", "8", "8", "8", "8", "8", "8", "8", "8",
#' "9", "9", "9", "9", "9", "9", "9", "9", "9", "9", "9", "9", "9",
#' "9", "10", "10", "10", "10", "10", "10", "10", "10", "10", "10",
#' "10", "10", "10", "10", "11", "11", "11", "11", "11", "11", "11",
#' "11", "11", "11", "11", "11", "11", "11"),
#' METRIC = c("ltest",
#' "lrbs_tst", "ldmb_bw5", "s_ldmb_bw5", "ldmb_bw4", "lrbs_bw4",
#' "lrbs_bw5", "s_lrbs_bw5", "lrbs_bw6", "s_lrbs_bw6", "Dcbf_g08",
#' "ldcbf_g08", "lrbs_g08", "s_rp100", "ltest", "lrbs_tst", "ldmb_bw5",
#' "s_ldmb_bw5", "ldmb_bw4", "lrbs_bw4", "lrbs_bw5", "s_lrbs_bw5",
#' "lrbs_bw6", "s_lrbs_bw6", "Dcbf_g08", "ldcbf_g08", "lrbs_g08",
#' "s_rp100", "ltest", "lrbs_tst", "ldmb_bw5", "s_ldmb_bw5", "ldmb_bw4",
#' "lrbs_bw4", "lrbs_bw5", "s_lrbs_bw5", "lrbs_bw6", "s_lrbs_bw6",
#' "Dcbf_g08", "ldcbf_g08", "lrbs_g08", "s_rp100", "ltest", "lrbs_tst",
#' "ldmb_bw5", "s_ldmb_bw5", "ldmb_bw4", "lrbs_bw4", "lrbs_bw5",
#' "s_lrbs_bw5", "lrbs_bw6", "s_lrbs_bw6", "Dcbf_g08", "ldcbf_g08",
#' "lrbs_g08", "s_rp100", "ltest", "lrbs_tst", "ldmb_bw5", "s_ldmb_bw5",
#' "ldmb_bw4", "lrbs_bw4", "lrbs_bw5", "s_lrbs_bw5", "lrbs_bw6",
#' "s_lrbs_bw6", "Dcbf_g08", "ldcbf_g08", "lrbs_g08", "s_rp100",
#' "ltest", "lrbs_tst", "ldmb_bw5", "s_ldmb_bw5", "ldmb_bw4", "lrbs_bw4",
#' "lrbs_bw5", "s_lrbs_bw5", "lrbs_bw6", "s_lrbs_bw6", "Dcbf_g08",
#' "ldcbf_g08", "lrbs_g08", "s_rp100", "ltest", "lrbs_tst", "ldmb_bw5",
#' "s_ldmb_bw5", "ldmb_bw4", "lrbs_bw4", "lrbs_bw5", "s_lrbs_bw5",
#' "lrbs_bw6", "s_lrbs_bw6", "Dcbf_g08", "ldcbf_g08", "lrbs_g08",
#' "s_rp100", "ltest", "lrbs_tst", "ldmb_bw5", "s_ldmb_bw5", "ldmb_bw4",
#' "lrbs_bw4", "lrbs_bw5", "s_lrbs_bw5", "lrbs_bw6", "s_lrbs_bw6",
#' "Dcbf_g08", "ldcbf_g08", "lrbs_g08", "s_rp100", "ltest", "lrbs_tst",
#' "ldmb_bw5", "s_ldmb_bw5", "ldmb_bw4", "lrbs_bw4", "lrbs_bw5",
#' "s_lrbs_bw5", "lrbs_bw6", "s_lrbs_bw6", "Dcbf_g08", "ldcbf_g08",
#' "lrbs_g08", "s_rp100", "ltest", "lrbs_tst", "ldmb_bw5", "s_ldmb_bw5",
#' "ldmb_bw4", "lrbs_bw4", "lrbs_bw5", "s_lrbs_bw5", "lrbs_bw6",
#' "s_lrbs_bw6", "Dcbf_g08", "ldcbf_g08", "lrbs_g08", "s_rp100",
#' "ltest", "lrbs_tst", "ldmb_bw5", "s_ldmb_bw5", "ldmb_bw4", "lrbs_bw4",
#' "lrbs_bw5", "s_lrbs_bw5", "lrbs_bw6", "s_lrbs_bw6", "Dcbf_g08",
#' "ldcbf_g08", "lrbs_g08", "s_rp100"),
#' RESULT = c(-0.164309429,
#' 1.164309429, 0.835690571, 0.831325766, -0.164309429, 1.164309429,
#' 0.164309429, 0.168674234, 1.164309429, 1.168674234, 113.9277151,
#' 2.056629387, -1.056629387, 2, -0.164309429, 1.164309429, NA,
#' NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 2, -0.164309429, 1.164309429,
#' 1.835690571, 1.831325766, 1.835690571, -0.835690571, -0.835690571,
#' -0.831325766, 0.164309429, 0.168674234, NA, NA, NA, 2, -0.164309429,
#' 1.164309429, 1.437750563, 1.426755179, 0.437750563, 0.562249437,
#' -0.437750563, -0.426755179, 0.562249437, 0.573244821, 146.4599268,
#' 2.165718813, -1.165718813, 2, 1.630980938, -1.090180938, 1.95385339,
#' 1.941619453, 1.953783007, -1.412983007, -1.41305339, -1.400819453,
#' -1.41305339, -1.400819453, 105.9303106, 2.025020246, -1.484220246,
#' 22.08840009, 2.67711418, 0.52688582, 2.832055426, 2.842904539,
#' 2.815915646, 0.388084354, 0.371944574, 0.361095461, 0.360474574,
#' 0.349625461, 1588.155449, 3.200893009, 0.003106991, 25.80377156,
#' 2.121824344, -1.364884344, 1.418982992, 1.436906832, 0.418982992,
#' 0.337957008, -0.662042992, -0.679966832, -0.694582992, -0.712506832,
#' 148.8617733, 2.172783188, -1.415843188, 4.41366392, -0.672485667,
#' 1.049125667, -0.523367109, -0.768424143, -0.523802116, 0.900442116,
#' 0.900007109, 1.145064143, 0.900007109, 1.145064143, 0.32636017,
#' -0.486302849, 0.862942849, 46.120606, 1.854477172, -1.535367172,
#' 1.677706613, 1.680788736, 1.677706613, -1.358596613, -1.358596613,
#' -1.361678736, -1.378856613, -1.381938736, NA, NA, NA, 21.92365464,
#' 1.355114917, -2.669624917, NA, NA, NA, NA, NA, NA, NA, NA, NA,
#' NA, NA, 5.80084933, 0.909920977, 0.279989023, 0.951821206, -0.080588008,
#' 0.949639991, 0.240270009, 0.238088794, 1.2704980082, NA, NA,
#' 17.97017283, 1.254552254, -0.064642254, 441.9807675)),
#' .Names = c("UID", "METRIC", "RESULT"), row.names = c(NA, 154L), class = "data.frame")
#' mets <- melt(subset(x, select = -PROTOCOL), id.var = 'UID', variable.name = 'METRIC', value.name = 'RESULT')
#' res <- merge(metsExpected, mets, by = c('UID', 'METRIC'), all.x = T)
#' res$diff <- abs(res$RESULT.x - res$RESULT.y)
#' res
calculateRefinedCriticalDiameter <- function(xdepth, xslope, rp100, xbkf_h,
v1w_msq, old.version = F, xbkf_w = NULL, xwidth = NULL){
if (old.version){
stopifnot(!is.null(xbkf_w), !is.null(xwidth))
Rw <- v1w_msq * 1000 * xbkf_w / xwidth
is.na(Rw) <- xwidth <= 0
} else {
Rw <- v1w_msq * 1000
}
Rbf <- 0.5 * ((xdepth - rp100) * 10 + xbkf_h * 1000)
Rbf <- pmax(0.1 * Rbf, Rbf - Rw)
critdia <- 13.7 * Rbf * xslope / 100
is.na(critdia) <- critdia <= 0
return(log10(critdia))
}
#' @rdname calculateRefinedCriticalDiameter
#' @export
calculateCrudeCriticalDiameter <- function(xdepth, xslope){
stopifnot(xslope <= 100, xslope >= 0)
critdia <- 13.7 * (0.5 * xdepth * 10) * (xslope / 100)
is.na(critdia) <- critdia <= 0
return(log10(critdia))
}
#' @rdname calculateRefinedCriticalDiameter
#' @export
calculateCriticalDiameterKaufmann <- function(xdepth, xslope, rp100, xbkf_h, v1w_msq, lsub_dmm){
rho <- 998
rhoSed <- 2650
g <- 9.807
Dbf_th <- xbkf_h + (xdepth / 100)
Rb3 <- 0.65 * Dbf_th
rp <- rp100 / 100
s <- xslope / 100
viscosity <- 1.02e-6
sub_dm <- 10 ^ (lsub_dmm - 3)# is (10^lsub_dmm) / 1000 because 10^lsub_dmm * 10^-3 = 10^(lsub_dmm - 3)
# Total hydraulic resistance, eqn 13b commented out. Original parameterization
# of eqn 13a is used.
#Ct_rpwd <- 1.21 * (rp^1.08) * ((rp + v1w_msq)^0.638) * (Dbf_th^-3.32)
reg.mat <- log10(cbind(10, rp, rp + v1w_msq, Dbf_th))
reg.pars <- c(0.0835, 1.08, 0.638, -3.32)
Ct_rpwd <- 10 ^ drop(reg.mat %*% reg.pars)
# Hydraulic resistance due to particles
Cp3_mill <- pmax(0.002, (1/8) * (2.03 * log10(12.2 * Rb3 / sub_dm))^(-2))
# Hydraulic resistance due to particles, eqn 10b, with restriction as described for eqn 16
Rrpw3 <- Rb3 * (pmin(1, Cp3_mill / Ct_rpwd) ^ (1 / 3))
# Reynolds number at bankfull, eqn 14
ReyP3 <- sqrt(g * Rb3 * s) * sub_dm / viscosity
# Shields parameter, eqn 15a, 15b. Note 15a uses an abbreviated value of
# the exponent in the case of small Reynolds numbers
Shld_Px3 <- ifelse(ReyP3 < 26,
0.04 * ReyP3^(-0.24),
0.5 * ((0.22 * (ReyP3^(-0.6))) + 0.06 * (10^(-7.7 * (ReyP3^(-0.6))))))
is.na(Shld_Px3) <- ReyP3 <= 0
# Bed surface particle critical diameter (mm) Dcbf* from eqn 16.
Dcbf_g08 <- 1000 * (rho * g * Rrpw3 * s) / (Shld_Px3 * (rhoSed - rho) * g)
ldcbf_g08 <- log10(Dcbf_g08)
return(cbind(Dcbf_g08, ldcbf_g08))
}
#' Calculate relative bed stability metrics.
#'
#' Calculates several version of relative bed stability metrics using two versions
#' of substrate diamter and several critical diameter versions.
#' @param lsub_dmm mean log10(substrate diameter)
#' @param lsub2dmm mean log2(substrate diameter)
#' @param ltest critical diameter (crude)
#' @param ldmb_bw4 critical diameter (refined) old version
#' @param ldmb_bw5 critical diameter (refined)
#' @param s_ldmb_bw5 critical diameter (refined) with s_rp100 rather than rp100
#' @param ldcbf_g08 critial diameter from Kauffman 2007
#' @return a matrix of relative bed stability metrics.
#' @export
calculateRelativeBedStability <- function(lsub_dmm, lsub2dmm, ltest, ldmb_bw4, ldmb_bw5, s_ldmb_bw5, ldcbf_g08){
ld <- cbind(ltest, ldmb_bw4, ldmb_bw5, s_ldmb_bw5, ldcbf_g08)
ans1 <- lsub_dmm - ld
ans2 <- lsub2dmm - ld[, 3:4, drop = F]
ans <- cbind(ans1, ans2)
colnames(ans) <- c('lrbs_tst', 'lrbs_bw4', 'lrbs_bw5', 's_lrbs_bw5', 'lrbs_g08', 'lrbs_bw6', 's_lrbs_bw6')
progressReport('Finished calculating relative bed stability metrics.')
return(ans)
}
#' Fill in missing LWD volume.
#'
#' Fill in missing LWD volume as 0 when no LWD fishcover (i.e., xfc_lwd == 0) otherwise we really can't guess.
#' @param v1w_msq a vector of v1w_msq metrics for each site; m^3 of lwd / m^2 channel
#' @param xfc_lwd a vector of xfc_lwd metrics for each site
#' @return a modified vector of v1w_msq metrics
#' @export
#' @examples
#' fillMissingLWD(NA, 1)
#' fillMissingLWD(NA, 0)
#' fillMissingLWD(3, 1)
fillMissingLWD <- function(v1w_msq, xfc_lwd){
nas <- is.na(v1w_msq)
no.lwd <- xfc_lwd == 0
v1w_msq[nas & no.lwd] <- 0
return(v1w_msq)
}
#' Not even sure what this modified metric is...
#'
#' Wish I knew.
#' @param xslope mean site slope
#' @param sddepth standard deviation of thalweg depth
#' @return the s_rp100 metric
#' @export
calculateS_RP100 <- function(xslope, sddepth){
k1 <- -0.44767
k2 <- 1.25381
k3 <- -0.20675
s_rp100 <- 10 ^ (k1 + k2 * log10(sddepth) + k3 * log10(xslope))
return(s_rp100)
}
#' Make slope slightly positive
#'
#' Make slopes slightly positive if they are less than 0.01.
#' @param slope a vector of slopes
#' @return a vector of slopes with slopes < .01, replaced by 0.01
#' @export
forceSlopeGt0 <- function(slope){
pmax(0.01, slope)
}
#' Convert boatable depth units from m to cm
#'
#' Convert boatable depth units from m to cm.
#' @param x a depth based metric measured in m
#' @param x logical, \code{TRUE} if the site was wadeable
#' @return a vector of depth metrics in cm
#' @export
fixBoatableDepthUnits <- function(x, is.wadeable){
i <- !is.wadeable
x[i] <- x[i] * 100
return(x)
}
jasonelaw/nrsa documentation built on Nov. 8, 2019, 11:34 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions calculateRefinedCriticalDiameter calculateCrudeCriticalDiameter calculateCriticalDiameterKaufmann calculateRelativeBedStability fillMissingLWD calculateS_RP100 forceSlopeGt0 fixBoatableDepthUnits
Documented in calculateCriticalDiameterKaufmann calculateCrudeCriticalDiameter calculateRefinedCriticalDiameter calculateRelativeBedStability calculateS_RP100 fillMissingLWD fixBoatableDepthUnits forceSlopeGt0
#' Calculate several versions of critical diameter: ltest, ldmb_bw4, ldmb_bw5, s_ldmb_bw5 and ldcbf_g08
#'
#' These functions can be used to calculate several versions of the critical particle
#' diameter:
#' \itemize{
#' \item a crude estimate via \code{calculateCrudeCriticalDiameter}
#' \item several more refined ones via \code{calculateRefinedCriticalDiameter}
#' \item a research based one due to Kauffman 2007 via \code{calculateCriticalDiameterKaufmann}
#' }
#'
#' \code{calculateRefinedCriticalDiameter} is used to calculate the metrics
#' ldmb_bw4 (with old.version = T), ldmb_bw5, and s_ldmb_bw5 (when the rp100
#' argument is as returned by \link{calculateS_RP100}).
#' @param xdepth average depth in centimeters
#' @param xslope average slope in percent
#' @param rp100 average residual depth in centimeters
#' @param xbkf_h average bankfull height in meters
#' @param v1w_msq volume lwd in cubic meters / 100 square meeters
#' @param old.version Calculate old version of logged critical diameter that were based on the
#' wood density when it accidentally used the wetted widths instead of the
#' bankfull widths.
#' @param xbkf_w average bankfull width in meters; must be supplied when old.version = T.
#' @param xwidth average width in meters; must be supplied when old.version = T.
#' @param lsub_dmm mean log10(substrate diameter) in millimeters
#' @return the log10(critical diameter)
#' @references
#' Kaufmann, P.R. et al., A roughness-corrected index of relative bed stability
#' for regional stream surveys. Geomorphology (2007)
#' @export
#' @examples
#' library(reshape2)
#'x <- structure(list(UID = c("1", "10", "11", "2", "3", "4", "5", "6", "7", "8", "9"),
#' lsub_dmm = c(1, -1.31451, 1.18991, 1, 1, 1, 0.5408, 3.204, 0.75694, 0.37664, 0.31911),
#' lsub2dmm = c(2, -1.31451, NA, 2, 2, 2, 0.5408, 3.19253, 0.7244, 0.37664, 0.29885),
#' rp100 = c(1, 10.3489, 112.73, 1, 1, 1, 19.0993, 28.3768, 6.4191, 27.2566, 22.2393),
#' sddepth = c(3.95494014, 8.8661, 1.69386, 3.95494014, 3.95494014,
#' 3.95494014, 27.7557, 44.2636, 10.4878, 22.6112, 28.7622),
#' v1w_msq = c(1, NA, 0.005818, NA, NA, 0.3, 4e-05, 0.01461, 0.37977, 0.00015, NA),
#' xbkf_h = c(1, 0.33636, 1.47273, 1, 1, 1, 0.75545, 0.61818, 0.3, 0.4, 0),
#' xbkf_w = c(1, 3.718, 255.273, 1, 1, 1, 21.818, 6.673, 2.927, 3.4818, 8.745),
#' xdepth = c(1, 42.67, 3.29555, 1, 1, 1, 51.16, 71.19, 24.007, 31.033, 66.51),
#' xfc_lwd = c(1, 0.00455, 0.00455, 1, 0, NA, 0, 0.01364, 0.11818, 0, 0),
#' xslope = c(1, 0.775, 0.036, 1, 1, 1, 1.22, 9.75, 8.05, 0, 1.57),
#' xwidth = c(1, 3.085, 163.364, 1, 1, 1, 10.0105, 5.255, 2.445, 2.385, 9.1474),
#' PROTOCOL = c("WADEABLE", "WADEABLE", "BOATABLE", "WADEABLE",
#' "WADEABLE", "WADEABLE", "WADEABLE",
#' "WADEABLE", "WADEABLE", "WADEABLE", "WADEABLE")),
#' .Names = c("UID", "lsub_dmm", "lsub2dmm", "rp100", "sddepth", "v1w_msq", "xbkf_h",
#' "xbkf_w", "xdepth", "xfc_lwd", "xslope", "xwidth", "PROTOCOL"),
#' row.names = c(NA, -11L),
#' class = "data.frame")
#' x$xslope <- forceSlopeGt0(x$xslope)
#' x$xdepth <- fixBoatableDepthUnits(x$xdepth, x$PROTOCOL == 'WADEABLE')
#' x$sddepth <- fixBoatableDepthUnits(x$sddepth, x$PROTOCOL == 'WADEABLE')
#' x$s_rp100 <- calculateS_RP100(x$xslope, x$sddepth)
#' x$v1w_msq <- fillMissingLWD(x$v1w_msq, x$xfc_lwd)
#' x$ltest <- with(x, calculateCrudeCriticalDiameter(xdepth, xslope))
#' x$ldmb_bw5 <- with(x, calculateRefinedCriticalDiameter(xdepth, xslope, rp100, xbkf_h, v1w_msq))
#' x$s_ldmb_bw5 <- with(x, calculateRefinedCriticalDiameter(xdepth, xslope, s_rp100, xbkf_h, v1w_msq))
#' x$ldmb_bw4 <- with(x, calculateRefinedCriticalDiameter(xdepth, xslope, rp100, xbkf_h, v1w_msq, T, xbkf_w, xwidth))
#' x <- cbind(x, with(x, calculateCriticalDiameterKaufmann(xdepth, xslope, rp100, xbkf_h, v1w_msq, lsub_dmm)))
#' x <- cbind(x, with(x, calculateRelativeBedStability(lsub_dmm, lsub2dmm, ltest, ldmb_bw4, ldmb_bw5, s_ldmb_bw5, ldcbf_g08)))
#' metsExpected <-
#' structure(list(UID = c("1", "1", "1", "1", "1", "1", "1", "1",
#' "1", "1", "1", "1", "1", "1", "2", "2", "2", "2", "2", "2", "2",
#' "2", "2", "2", "2", "2", "2", "2", "3", "3", "3", "3", "3", "3",
#' "3", "3", "3", "3", "3", "3", "3", "3", "4", "4", "4", "4", "4",
#' "4", "4", "4", "4", "4", "4", "4", "4", "4", "5", "5", "5", "5",
#' "5", "5", "5", "5", "5", "5", "5", "5", "5", "5", "6", "6", "6",
#' "6", "6", "6", "6", "6", "6", "6", "6", "6", "6", "6", "7", "7",
#' "7", "7", "7", "7", "7", "7", "7", "7", "7", "7", "7", "7", "8",
#' "8", "8", "8", "8", "8", "8", "8", "8", "8", "8", "8", "8", "8",
#' "9", "9", "9", "9", "9", "9", "9", "9", "9", "9", "9", "9", "9",
#' "9", "10", "10", "10", "10", "10", "10", "10", "10", "10", "10",
#' "10", "10", "10", "10", "11", "11", "11", "11", "11", "11", "11",
#' "11", "11", "11", "11", "11", "11", "11"),
#' METRIC = c("ltest",
#' "lrbs_tst", "ldmb_bw5", "s_ldmb_bw5", "ldmb_bw4", "lrbs_bw4",
#' "lrbs_bw5", "s_lrbs_bw5", "lrbs_bw6", "s_lrbs_bw6", "Dcbf_g08",
#' "ldcbf_g08", "lrbs_g08", "s_rp100", "ltest", "lrbs_tst", "ldmb_bw5",
#' "s_ldmb_bw5", "ldmb_bw4", "lrbs_bw4", "lrbs_bw5", "s_lrbs_bw5",
#' "lrbs_bw6", "s_lrbs_bw6", "Dcbf_g08", "ldcbf_g08", "lrbs_g08",
#' "s_rp100", "ltest", "lrbs_tst", "ldmb_bw5", "s_ldmb_bw5", "ldmb_bw4",
#' "lrbs_bw4", "lrbs_bw5", "s_lrbs_bw5", "lrbs_bw6", "s_lrbs_bw6",
#' "Dcbf_g08", "ldcbf_g08", "lrbs_g08", "s_rp100", "ltest", "lrbs_tst",
#' "ldmb_bw5", "s_ldmb_bw5", "ldmb_bw4", "lrbs_bw4", "lrbs_bw5",
#' "s_lrbs_bw5", "lrbs_bw6", "s_lrbs_bw6", "Dcbf_g08", "ldcbf_g08",
#' "lrbs_g08", "s_rp100", "ltest", "lrbs_tst", "ldmb_bw5", "s_ldmb_bw5",
#' "ldmb_bw4", "lrbs_bw4", "lrbs_bw5", "s_lrbs_bw5", "lrbs_bw6",
#' "s_lrbs_bw6", "Dcbf_g08", "ldcbf_g08", "lrbs_g08", "s_rp100",
#' "ltest", "lrbs_tst", "ldmb_bw5", "s_ldmb_bw5", "ldmb_bw4", "lrbs_bw4",
#' "lrbs_bw5", "s_lrbs_bw5", "lrbs_bw6", "s_lrbs_bw6", "Dcbf_g08",
#' "ldcbf_g08", "lrbs_g08", "s_rp100", "ltest", "lrbs_tst", "ldmb_bw5",
#' "s_ldmb_bw5", "ldmb_bw4", "lrbs_bw4", "lrbs_bw5", "s_lrbs_bw5",
#' "lrbs_bw6", "s_lrbs_bw6", "Dcbf_g08", "ldcbf_g08", "lrbs_g08",
#' "s_rp100", "ltest", "lrbs_tst", "ldmb_bw5", "s_ldmb_bw5", "ldmb_bw4",
#' "lrbs_bw4", "lrbs_bw5", "s_lrbs_bw5", "lrbs_bw6", "s_lrbs_bw6",
#' "Dcbf_g08", "ldcbf_g08", "lrbs_g08", "s_rp100", "ltest", "lrbs_tst",
#' "ldmb_bw5", "s_ldmb_bw5", "ldmb_bw4", "lrbs_bw4", "lrbs_bw5",
#' "s_lrbs_bw5", "lrbs_bw6", "s_lrbs_bw6", "Dcbf_g08", "ldcbf_g08",
#' "lrbs_g08", "s_rp100", "ltest", "lrbs_tst", "ldmb_bw5", "s_ldmb_bw5",
#' "ldmb_bw4", "lrbs_bw4", "lrbs_bw5", "s_lrbs_bw5", "lrbs_bw6",
#' "s_lrbs_bw6", "Dcbf_g08", "ldcbf_g08", "lrbs_g08", "s_rp100",
#' "ltest", "lrbs_tst", "ldmb_bw5", "s_ldmb_bw5", "ldmb_bw4", "lrbs_bw4",
#' "lrbs_bw5", "s_lrbs_bw5", "lrbs_bw6", "s_lrbs_bw6", "Dcbf_g08",
#' "ldcbf_g08", "lrbs_g08", "s_rp100"),
#' RESULT = c(-0.164309429,
#' 1.164309429, 0.835690571, 0.831325766, -0.164309429, 1.164309429,
#' 0.164309429, 0.168674234, 1.164309429, 1.168674234, 113.9277151,
#' 2.056629387, -1.056629387, 2, -0.164309429, 1.164309429, NA,
#' NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 2, -0.164309429, 1.164309429,
#' 1.835690571, 1.831325766, 1.835690571, -0.835690571, -0.835690571,
#' -0.831325766, 0.164309429, 0.168674234, NA, NA, NA, 2, -0.164309429,
#' 1.164309429, 1.437750563, 1.426755179, 0.437750563, 0.562249437,
#' -0.437750563, -0.426755179, 0.562249437, 0.573244821, 146.4599268,
#' 2.165718813, -1.165718813, 2, 1.630980938, -1.090180938, 1.95385339,
#' 1.941619453, 1.953783007, -1.412983007, -1.41305339, -1.400819453,
#' -1.41305339, -1.400819453, 105.9303106, 2.025020246, -1.484220246,
#' 22.08840009, 2.67711418, 0.52688582, 2.832055426, 2.842904539,
#' 2.815915646, 0.388084354, 0.371944574, 0.361095461, 0.360474574,
#' 0.349625461, 1588.155449, 3.200893009, 0.003106991, 25.80377156,
#' 2.121824344, -1.364884344, 1.418982992, 1.436906832, 0.418982992,
#' 0.337957008, -0.662042992, -0.679966832, -0.694582992, -0.712506832,
#' 148.8617733, 2.172783188, -1.415843188, 4.41366392, -0.672485667,
#' 1.049125667, -0.523367109, -0.768424143, -0.523802116, 0.900442116,
#' 0.900007109, 1.145064143, 0.900007109, 1.145064143, 0.32636017,
#' -0.486302849, 0.862942849, 46.120606, 1.854477172, -1.535367172,
#' 1.677706613, 1.680788736, 1.677706613, -1.358596613, -1.358596613,
#' -1.361678736, -1.378856613, -1.381938736, NA, NA, NA, 21.92365464,
#' 1.355114917, -2.669624917, NA, NA, NA, NA, NA, NA, NA, NA, NA,
#' NA, NA, 5.80084933, 0.909920977, 0.279989023, 0.951821206, -0.080588008,
#' 0.949639991, 0.240270009, 0.238088794, 1.2704980082, NA, NA,
#' 17.97017283, 1.254552254, -0.064642254, 441.9807675)),
#' .Names = c("UID", "METRIC", "RESULT"), row.names = c(NA, 154L), class = "data.frame")
#' mets <- melt(subset(x, select = -PROTOCOL), id.var = 'UID', variable.name = 'METRIC', value.name = 'RESULT')
#' res <- merge(metsExpected, mets, by = c('UID', 'METRIC'), all.x = T)
#' res$diff <- abs(res$RESULT.x - res$RESULT.y)
#' res
calculateRefinedCriticalDiameter <- function(xdepth, xslope, rp100, xbkf_h,
v1w_msq, old.version = F, xbkf_w = NULL, xwidth = NULL){
if (old.version){
stopifnot(!is.null(xbkf_w), !is.null(xwidth))
Rw <- v1w_msq * 1000 * xbkf_w / xwidth
is.na(Rw) <- xwidth <= 0
} else {
Rw <- v1w_msq * 1000
}
Rbf <- 0.5 * ((xdepth - rp100) * 10 + xbkf_h * 1000)
Rbf <- pmax(0.1 * Rbf, Rbf - Rw)
critdia <- 13.7 * Rbf * xslope / 100
is.na(critdia) <- critdia <= 0
return(log10(critdia))
}
#' @rdname calculateRefinedCriticalDiameter
#' @export
calculateCrudeCriticalDiameter <- function(xdepth, xslope){
stopifnot(xslope <= 100, xslope >= 0)
critdia <- 13.7 * (0.5 * xdepth * 10) * (xslope / 100)
is.na(critdia) <- critdia <= 0
return(log10(critdia))
}
#' @rdname calculateRefinedCriticalDiameter
#' @export
calculateCriticalDiameterKaufmann <- function(xdepth, xslope, rp100, xbkf_h, v1w_msq, lsub_dmm){
rho <- 998
rhoSed <- 2650
g <- 9.807
Dbf_th <- xbkf_h + (xdepth / 100)
Rb3 <- 0.65 * Dbf_th
rp <- rp100 / 100
s <- xslope / 100
viscosity <- 1.02e-6
sub_dm <- 10 ^ (lsub_dmm - 3)# is (10^lsub_dmm) / 1000 because 10^lsub_dmm * 10^-3 = 10^(lsub_dmm - 3)
# Total hydraulic resistance, eqn 13b commented out. Original parameterization
# of eqn 13a is used.
#Ct_rpwd <- 1.21 * (rp^1.08) * ((rp + v1w_msq)^0.638) * (Dbf_th^-3.32)
reg.mat <- log10(cbind(10, rp, rp + v1w_msq, Dbf_th))
reg.pars <- c(0.0835, 1.08, 0.638, -3.32)
Ct_rpwd <- 10 ^ drop(reg.mat %*% reg.pars)
# Hydraulic resistance due to particles
Cp3_mill <- pmax(0.002, (1/8) * (2.03 * log10(12.2 * Rb3 / sub_dm))^(-2))
# Hydraulic resistance due to particles, eqn 10b, with restriction as described for eqn 16
Rrpw3 <- Rb3 * (pmin(1, Cp3_mill / Ct_rpwd) ^ (1 / 3))
# Reynolds number at bankfull, eqn 14
ReyP3 <- sqrt(g * Rb3 * s) * sub_dm / viscosity
# Shields parameter, eqn 15a, 15b. Note 15a uses an abbreviated value of
# the exponent in the case of small Reynolds numbers
Shld_Px3 <- ifelse(ReyP3 < 26,
0.04 * ReyP3^(-0.24),
0.5 * ((0.22 * (ReyP3^(-0.6))) + 0.06 * (10^(-7.7 * (ReyP3^(-0.6))))))
is.na(Shld_Px3) <- ReyP3 <= 0
# Bed surface particle critical diameter (mm) Dcbf* from eqn 16.
Dcbf_g08 <- 1000 * (rho * g * Rrpw3 * s) / (Shld_Px3 * (rhoSed - rho) * g)
ldcbf_g08 <- log10(Dcbf_g08)
return(cbind(Dcbf_g08, ldcbf_g08))
}
#' Calculate relative bed stability metrics.
#'
#' Calculates several version of relative bed stability metrics using two versions
#' of substrate diamter and several critical diameter versions.
#' @param lsub_dmm mean log10(substrate diameter)
#' @param lsub2dmm mean log2(substrate diameter)
#' @param ltest critical diameter (crude)
#' @param ldmb_bw4 critical diameter (refined) old version
#' @param ldmb_bw5 critical diameter (refined)
#' @param s_ldmb_bw5 critical diameter (refined) with s_rp100 rather than rp100
#' @param ldcbf_g08 critial diameter from Kauffman 2007
#' @return a matrix of relative bed stability metrics.
#' @export
calculateRelativeBedStability <- function(lsub_dmm, lsub2dmm, ltest, ldmb_bw4, ldmb_bw5, s_ldmb_bw5, ldcbf_g08){
ld <- cbind(ltest, ldmb_bw4, ldmb_bw5, s_ldmb_bw5, ldcbf_g08)
ans1 <- lsub_dmm - ld
ans2 <- lsub2dmm - ld[, 3:4, drop = F]
ans <- cbind(ans1, ans2)
colnames(ans) <- c('lrbs_tst', 'lrbs_bw4', 'lrbs_bw5', 's_lrbs_bw5', 'lrbs_g08', 'lrbs_bw6', 's_lrbs_bw6')
progressReport('Finished calculating relative bed stability metrics.')
return(ans)
}
#' Fill in missing LWD volume.
#'
#' Fill in missing LWD volume as 0 when no LWD fishcover (i.e., xfc_lwd == 0) otherwise we really can't guess.
#' @param v1w_msq a vector of v1w_msq metrics for each site; m^3 of lwd / m^2 channel
#' @param xfc_lwd a vector of xfc_lwd metrics for each site
#' @return a modified vector of v1w_msq metrics
#' @export
#' @examples
#' fillMissingLWD(NA, 1)
#' fillMissingLWD(NA, 0)
#' fillMissingLWD(3, 1)
fillMissingLWD <- function(v1w_msq, xfc_lwd){
nas <- is.na(v1w_msq)
no.lwd <- xfc_lwd == 0
v1w_msq[nas & no.lwd] <- 0
return(v1w_msq)
}
#' Not even sure what this modified metric is...
#'
#' Wish I knew.
#' @param xslope mean site slope
#' @param sddepth standard deviation of thalweg depth
#' @return the s_rp100 metric
#' @export
calculateS_RP100 <- function(xslope, sddepth){
k1 <- -0.44767
k2 <- 1.25381
k3 <- -0.20675
s_rp100 <- 10 ^ (k1 + k2 * log10(sddepth) + k3 * log10(xslope))
return(s_rp100)
}
#' Make slope slightly positive
#'
#' Make slopes slightly positive if they are less than 0.01.
#' @param slope a vector of slopes
#' @return a vector of slopes with slopes < .01, replaced by 0.01
#' @export
forceSlopeGt0 <- function(slope){
pmax(0.01, slope)
}
#' Convert boatable depth units from m to cm
#'
#' Convert boatable depth units from m to cm.
#' @param x a depth based metric measured in m
#' @param x logical, \code{TRUE} if the site was wadeable
#' @return a vector of depth metrics in cm
#' @export
fixBoatableDepthUnits <- function(x, is.wadeable){
i <- !is.wadeable
x[i] <- x[i] * 100
return(x)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.