R/DAYCENT_utils.R

In bmcnellis/RSFIA: FIA Tree Mortality Analysis

#' @title Utilities for DAYCENT modelling, including data downloading/file formatting
#' @description
#'
#' Multiple utilities for working with DAYCENT.
#'
#' Function parameters, in no particular order:
#'
#' @param sand Fraction of sand
#' @param clay Fraction of clay
#' @param soc Soil organic matter percent
#' @param DF Density factor between 0.9 and 1.3
#' @param gravel Gravel percent by weight
#' @param thetas Saturation water content, w/o compaction
#'
#' Multiple functions came from: https://github.com/ldemaz/rcropmod
#'
#' @author Jeff Kent, Brandon McNellis, Lyndon Estes
#' @examples DAYCENT_utils()
DAYCENT_utils <- function() {
  message('Available functions:')
  cat('daymet_to_daycent', 'ksat', 'field_cap_df', 'theta_sdf',
      'theta_s', 'field_cap', 'ro_df', 'wilt_point', 'bdens', '\n')
}
#' @describeIn DAYCENT_utils Reads .dmw files and converts to .wth
daymet_to_daycent <- function(mydmwfiles) {
  # DayCent files get the same name as the input dmw files, but with swapped extensions
  # following is how I make the list of files:
  # mydmwfiles <- list.files(path = "weather", pattern = "dmw", full.names = T)
  #REQUIRES the lubridate and stringr libraries
  require(lubridate)
  require(stringr)
  for (j in 1:length(mydmwfiles))
  {
    dmw <- read.table(file=mydmwfiles[j], skip = 7, header=T, sep=",") #skips 7-line header in DayMet
    dmw$date <- as.Date(dmw$yday-1, origin=paste(dmw$year,"01","01",sep="-")) #converts year-doy from the DayMet file into a date-class so that other values (day of month, month of year) can be extracted for the .wth file
    wth <- cbind(mday(dmw$date),month(dmw$date),year(dmw$date),yday(dmw$date),dmw$tmax..deg.c.,dmw$tmin..deg.c.,dmw$prcp..mm.day./10) #reorders and converts the DayMet data fields to DayCent format
    wthname <- str_replace(string=mydmwfiles[j], pattern="dmw", replacement="wth") #make the string for naming the .wth output file
    write.table(x=wth,file=wthname,quote=F,row.names=F,col.names=F,sep="\t")  #writes out the file without extra junk
  }
}
#' @describeIn DAYCENT_utils Saturated hydraulic conductivity, including gravel effects.
#' @export
ksat <- function(sand, clay, soc, DF = 1, gravel = 0) {
  # sand/clay - fraction of sand/clay
  # soc - soil organic matter (%)
  # DF - density factor (0.9 - 1.3)
  # gravel Gravel percent by weight
  # Author: https://github.com/ldemaz/rcropmod
  fcdf <- field_cap_df(sand, clay, soc, DF)
  wp <- wilt_point(sand, clay, soc)
  lambda <- (log(fcdf) - log(wp)) / (log(1500) - log(33))  # = 1/Beta
  thetas <- theta_s(sand, clay, soc)  # theta_sdf no density effects
  mdens <- bdens(thetas, DF, gravel = 0)  # BD no gravel to get matric density
  thetasdf <- theta_sdf(sand, clay, soc, DF = DF)  # ThetaSDF w/density effects
  theta_sdf_fcdf <- thetasdf - fcdf
  theta_sdf_fcdf <- ifelse(theta_sdf_fcdf < 0, 0, theta_sdf_fcdf) # FC ! > por.
  kbks <- (1 - gravel) / (1 - gravel * (1 - 1.5 * (mdens / 2.65)))
  ks <- 1930 * (theta_sdf_fcdf)^(3 - lambda) * kbks
  return(ks)
}
#' @describeIn DAYCENT_utils Calculated field capacity accounting for compaction
#' @export
field_cap_df <- function(sand, clay, soc, DF) {
  # Parameters as in ksat
  # Author: https://github.com/ldemaz/rcropmod
  thetas <- theta_sdf(sand, clay, soc, DF = 1)  # Normal theta_s
  thetasdf <- theta_sdf(sand, clay, soc, DF)  # theta_s with compaction
  fcdf <- field_cap(sand, clay, soc) - 0.2 * (thetas - thetasdf)
  return(fcdf)
}
#' @describeIn DAYCENT_utils Calculates saturated water content, accounting for compaction
#' @export
theta_sdf <- function(sand, clay, soc, DF) {
  # Parameters as in ksat
  # Author: https://github.com/ldemaz/rcropmod
  thetas <- theta_s(sand, clay, soc)
  rodf <- ro_df(thetas, DF)
  thetasdf <- 1 - (rodf / 2.65)
  return(thetasdf)
}
#' @describeIn DAYCENT_utils Calculates saturated moisture content
#' @export
theta_s <- function(sand, clay, soc) {
  # Author: https://github.com/ldemaz/rcropmod
  thetas_33t <-  0.278 * sand + 0.034 * clay + 0.022 * soc -
    (0.018 * sand * soc) - (0.027 * clay * soc) - (0.584 * sand * clay) + 0.078
  thetas_33 <- thetas_33t + (0.636 * thetas_33t - 0.107)
  theta33 <- field_cap(sand, clay, soc)
  thetas <- theta33 + thetas_33 - 0.097 * sand + 0.043
  # Returns field capacity at 33 kPa
  return(thetas)
}
#' @describeIn DAYCENT_utils Calculates field capacity
#' @export
field_cap <- function(sand, clay, soc) {
  # Author: https://github.com/ldemaz/rcropmod
  theta33t <- -0.251 * sand + 0.195 * clay + 0.011 * soc +
    (0.006 * sand * soc) - (0.027 * clay * soc) +
    (0.452 * sand * clay) + 0.299
  theta33 <- theta33t + ((1.283 * theta33t^2) - 0.374 * theta33t - 0.015)
  # Returns field capacity at 33 kPa
  return(theta33)
}
#' @describeIn DAYCENT_utils Matric density accounting for compaction
#' @export
ro_df <- function(thetas, DF = 1) {
  # Author: https://github.com/ldemaz/rcropmod
  rodf <- ((1 - thetas) * 2.65) * DF
  # Returns matric density
  return(rodf)
}
#' @describeIn DAYCENT_utils Calculate wilting point
#' @export
wilt_point <- function(sand, clay, soc) {
  # Author: https://github.com/ldemaz/rcropmod
  theta1500t <- -0.024 * sand + 0.487 * clay + 0.006 * soc +
    (0.005 * sand * soc) - (0.013 * clay * soc) + (0.068 * sand * clay) +
    0.031
  theta1500 <- theta1500t + (0.14 * theta1500t - 0.02)
  # Returns wilting point at 1500 kpa
  return(theta1500)
}
#' @describeIn DAYCENT_utils Bulk density accounting for compaction plus gravel
#' @export
bdens <- function(thetas, DF = 1, gravel = 0) {
  rodf <- ro_df(thetas, DF)
  gravel_pctv <- ((rodf / 2.65 ) * gravel) / (1 - gravel * ( 1 - rodf / 2.65))
  ro_b  <- gravel_pctv * 2.65 + (1 - gravel_pctv) * rodf
  return(ro_b)
}