R/package-rch.R

In rogiersbart/RMODFLOW: Pre- and post-processing of MODFLOW files

#' Create an \code{RMODFLOW} rch object.
#' 
#' \code{rmf_create_rch} creates an \code{RMODFLOW} rch object
#' 
#' @param ... \code{rmf_2d_arrays} (possibly of class \code{rmf_parameter}) or a single \code{list} with \code{rmf_2d_arrays} (possibly of class \code{rmf_parameter}) elements; defines the recharge values. See details.
#' @param dis \code{RMODFLOW} dis object
#' @param nrchop recharge option code; defaults to 3 (recharge is applied to the highest active cell in each vertical column)
#' @param irchcb flag and unit number for writing cell-by-cell flow terms; defaults to 0 
#' @param irch a single \code{rmf_2d_array} or a list of \code{rmf_2d_arrays} specifying the layer numbers defining in which layer recharge is applied. The \code{'kper'} attribute of the arrays define the stress period in which the array is active, see details. Only used when \code{nrchop = 2}. Defaults to NULL
#' @param irchpf numeric of length 1 or length \code{dis$nper}; optional format code for printing the \code{RECH} variable it has been defined by parameters; defaults to -1 (no printing) for all stress periods
#' @details the \code{rmf_2d_arrays} should have \code{kper} attributes specifying the stress period in which they are active. This is also true for the irch arrays. There can be only one non-parameter array active per stress periods. Multiple parameters are however allowed per stress period.
#' @return \code{RMODFLOW} rch object
#' @export
#' @seealso \code{\link{rmf_read_rch}}, \code{\link{rmf_write_rch}}, \url{https://water.usgs.gov/ogw/modflow/MODFLOW-2005-Guide/index.html?rch.htm}

rmf_create_rch <- function(...,
                           dis,
                           nrchop = 3,
                           irchcb = 0,
                           irch = NULL,
                           irchpf = -1
) {
  
  arg <- rmfi_create_bc_array(arg = list(...), dis = dis)
  
  # create rch object
  obj <- arg[c("np", "instances")]
  obj$nrchop <- nrchop
  obj$irchcb <- irchcb
  obj$recharge <- arg$data
  if(arg$np > 0) obj$parameter_values <- arg$parameter_values
  obj$kper <- arg$kper
  
  # irch
  if(nrchop == 2) {
    if(is.null(irch)) stop('Please supply a irch argument when nrchop = 2', call. = FALSE)
    if(!inherits(irch, 'list')) irch <- list(irch)
    irch <- lapply(irch, function(i) {r <- apply(i, MARGIN = 1:length(dim(i)), function(x) as.integer(x)); attributes(r) <- attributes(i); r})
    obj$irch <- irch
    names(obj$irch) <- paste('irch', 1:length(irch), sep = '_')
    obj$kper$irch <- NA_character_
    for(i in 1:length(irch)) {
      obj$kper$irch[c(1:dis$nper) %in% attr(irch[[i]],'kper')] <- names(obj$irch)[i]
    }

    # check if multiple irch arrays are active
    irch_err <- unlist(lapply(irch, function(i) attr(i, 'kper')))
    if(any(duplicated(irch_err))) stop(paste('There can be only 1 active irch array per stress period. Stress period(s)', sort(unique(irch_err[duplicated(irch_err)])), 'have multiple active arrays.'), call. = FALSE)
  } 
  obj$irchpf <- irchpf
  
  class(obj) <- c('rch', 'rmf_package')
  return(obj)
  
}

#' Read a MODFLOW recharge file
#' 
#' \code{rmf_read_rch} reads in a MODFLOW recharge file and returns it as an \code{RMODFLOW} rch object.
#'
#' @param file filename; typically '*.rch'
#' @param dis an \code{RMODFLOW} dis object
#' @param mlt a \code{RMODFLOW} mlt object. Only needed when reading parameter arrays defined by multiplier arrays
#' @param zon a \code{RMODFLOW} zon object. Only needed when reading parameter arrays defined by zone arrays
#' @param ... arguments passed to \code{rmfi_parse_array}. Can be ignored when input arrays are free-format and INTERNAL or CONSTANT.
#' @return \code{RMODFLOW} rch object
#' @export
#' @seealso \code{\link{rmf_write_rch}}, \code{\link{rmf_create_rch}}, \url{https://water.usgs.gov/ogw/modflow/MODFLOW-2005-Guide/index.html?rch.htm}

rmf_read_rch <-  function(file = {cat('Please select rch file ...\n'); file.choose()},
                          dis = {cat('Please select corresponding dis file ...\n'); rmf_read_dis(file.choose())},
                          mlt = NULL,
                          zon = NULL,
                          ... ){
  
  lines <- readr::read_lines(file, lazy = FALSE)
  rmf_arrays <- list()
  
  # data set 0
  data_set_0 <-  rmfi_parse_comments(lines)
  comments <-  data_set_0$comments
  lines <-  data_set_0$remaining_lines
  rm(data_set_0)
  
  # data set 1
  data_set_1 <- rmfi_parse_variables(lines, character = TRUE)
  
  if('PARAMETER' %in% toupper(data_set_1$variables)) {
    np <-  as.numeric(data_set_1$variables[2])
    lines <-  data_set_1$remaining_lines
  }  else {
    np <- 0
  }
  rm(data_set_1)
  
  # data set 2
  data_set_2 <-  rmfi_parse_variables(lines, n=2, ...)
  nrchop <- as.numeric(data_set_2$variables[1])
  irchcb <- as.numeric(data_set_2$variables[2])
  lines <-  data_set_2$remaining_lines
  rm(data_set_2)
  irch <- rmfi_ifelse0(nrchop == 2, list(), NULL)
  
  # parameters: data set 3 & 4
  if(np > 0) {
    data_set_3 <- rmfi_parse_array_parameters(lines, dis = dis, np = np, mlt = mlt, zon = zon)
    rmf_arrays <- data_set_3$parameters
    lines <- data_set_3$remaining_lines
    rm(data_set_3)
  }
  
  # function for setting kper attribute of parameter the same as previous kper
  previous_kper <- function(k, kper) {
    set_previous_kper <- function(kk, kper) {
      if(!is.null(attr(kk, 'kper')) && kper-1 %in% attr(kk, 'kper')) {
        attr(kk, 'kper') <- c(attr(kk, 'kper'), kper)
      }
      return(kk)
    }
    
    if(is.list(k) && !is.null(attr(k[[1]], 'instnam'))) {
      k <- lapply(k, set_previous_kper, kper = kper)
    } else {
      k <- set_previous_kper(k, kper)
    }

    return(k)
  }
  
  for(i in 1:dis$nper){
    # data set 5
    data_set_5 <-  rmfi_parse_variables(lines, n=2, ...)
    inrech <- as.numeric(data_set_5$variables[1])
    if(nrchop == 2) inirch <- as.numeric(data_set_5$variables[2])
    lines <- data_set_5$remaining_lines
    rm(data_set_5)
    
    # data set 6-7
    irchpf <- NULL
    
    if(np == 0) {
      
      if(inrech >= 0) {
        data_set_6 <- rmfi_parse_array(lines, dis$nrow, dis$ncol, 1, ndim = 2, file = file, ...)
        rmf_arrays[[length(rmf_arrays) + 1]] <- structure(data_set_6$array, kper = i)
        lines <- data_set_6$remaining_lines
        rm(data_set_6)
      } else if(inrech < 0 && i > 1) {
        attr(rmf_arrays[[length(rmf_arrays)]], 'kper') <- c(attr(rmf_arrays[[length(rmf_arrays)]], 'kper'), i)
      }
      
    } else {
      # parameters
      if(inrech > 0) {
        for(j in 1:inrech){
          # data set 7
          data_set_7 <-  rmfi_parse_variables(lines, character = TRUE)
          p_name <-  as.character(data_set_7$variables[1])
          if(is.list(rmf_arrays[[p_name]]) && !is.null(attr(rmf_arrays[[p_name]][[1]], 'instnam'))) {
            i_name <- data_set_7$variables[2]
            if(length(data_set_7$variables) > 2 && !is.na(suppressWarnings(as.numeric(data_set_7$variables[3])))) {
              irchpf[i] <- as.numeric(data_set_7$variables[3])
            }
            attr(rmf_arrays[[p_name]][[i_name]], 'kper') <- c(attr(rmf_arrays[[p_name]][[i_name]], 'kper'), i)
          } else {
            i_name <- NULL
            if(length(data_set_7$variables) > 1 && !is.na(suppressWarnings(as.numeric(data_set_7$variables[2])))) {
              irchpf[i] <- as.numeric(data_set_7$variables[2])
            }
            attr(rmf_arrays[[p_name]], 'kper') <- c(attr(rmf_arrays[[p_name]], 'kper'), i)
          }
          
          # rmf_arrays <- lapply(rmf_arrays, set_kper, p_name = p_name, i_name = i_name, kper = i)
          
          lines <- data_set_7$remaining_lines
          rm(data_set_7)
          
        }
        
      } else if(inrech < 0 && i > 1) {
        rmf_arrays <- lapply(rmf_arrays, previous_kper, kper = i)
      }
      
    }
    
    # data set 8
    if(nrchop == 2) {
      if(inirch >= 0) {
        data_set_8 <- rmfi_parse_array(lines, dis$nrow, dis$ncol, 1, ndim = 2, file = file, integer = TRUE, ...)
        irch[[length(irch) + 1]] <- rmf_create_array(structure(apply(data_set_8$array, 1:length(dim(data_set_8$array)), function(i) as.integer(i)), kper = i))
        lines <- data_set_8$remaining_lines
        rm(data_set_8)
      } else if(inirch < 0 && i > 1) {
        attr(irch[[length(irch)]], 'kper') <- c(attr(irch[[length(irch)]], 'kper'), i)
      }
    }
  }
  list_arrays <- function(i) {
    if(is.list(i) && !is.null(attr(i[[1]], 'instnam'))) {
      return(i)
    } else {
      return(list(i))
    }
  }
  rmf_arrays <- lapply(rmf_arrays, list_arrays)
  rmf_arrays <- do.call(c, rmf_arrays)
  rmf_arrays <- lapply(rmf_arrays, function(i) rmfi_ifelse0(is.null(attr(i, 'kper')), structure(i, kper = 0), i))
  
  rch <- rmf_create_rch(rmf_arrays, dis = dis, nrchop = nrchop, irchcb = irchcb, irch = irch, irchpf = rmfi_ifelse0(is.null(irchpf), -1, irchpf))
  comment(rch) <- comments
  return(rch)
}

#' Write a MODFLOW recharge file
#'
#' \code{rmf_write_rch} writes a MODFLOW recharge file based on an \code{RMODFLOW} rch object
#' 
#' @param rch an \code{RMODFLOW} rch object
#' @param dis an \code{RMODFLOW} dis object
#' @param file filename to write to; typically '*.rch'
#' @param iprn format code for printing arrays in the listing file; defaults to -1 (no printing)
#' @param ... arguments passed to \code{rmfi_write_variables} and \code{rmfi_write_array}
#' 
#' @return \code{NULL}
#' @export
#' @seealso \code{\link{rmf_read_rch}}, \code{\link{rmf_create_rch}}, \url{https://water.usgs.gov/ogw/modflow/MODFLOW-2005-Guide/index.html?rch.htm}

rmf_write_rch <-  function(rch,
                           dis = {cat('Please select corresponding dis file ...\n'); rmf_read_dis(file.choose())},
                           file={cat('Please choose rch file to overwrite or provide new filename ...\n'); file.choose()},
                           iprn = -1,
                           ...){
  
  # data set 0
  v <- packageDescription("RMODFLOW")$Version
  cat(paste(paste('# MODFLOW Recharge Package created by RMODFLOW, version'),v,'\n'), file = file)
  cat(paste('#', comment(rch)), sep='\n', file=file, append=TRUE)
  
  # data set 1
  if(rch$np > 0) rmfi_write_variables('PARAMETER', as.integer(rch$np), file=file)
  
  # data set 2
  rmfi_write_variables(rch$nrchop, rch$irchcb, file=file, integer = TRUE, ...)
  
  # parameters
  partyp <- 'RCH'
  if(rch$np > 0) {
    parm_names <- names(rch$parameter_values)
    tv_parm <- rep(FALSE, rch$np)
    if(!is.null(rch$instances)) tv_parm <- rch$instances > 0
    rmfi_write_array_parameters(obj = rch, arrays = rch$recharge, file = file, partyp = 'RCH', ...)
  }
  
  # stress periods
  for (i in 1:dis$nper){
    
    check_prev <- function(kper, i) {
      drop_names <- which(names(kper) %in% c('kper', 'irch'))
      df <- kper[c(i-1,i), -drop_names, drop = FALSE]
      identical(c(df[2,]), c(df[1,]))
    }
    
    # data set 5
    # inrech
    # drop_id <- which(colnames(rch$kper) %in% c('kper', 'irch'))
    # names_act <- colnames(rch$kper)[which(rch$kper[i,which(!is.na(rch$kper[i,]))] != FALSE)[-drop_id]]
    names_act <- colnames(rch$kper)[which(rch$kper[i,which(!is.na(rch$kper[i,]))] != FALSE)]
    names_act <- names_act[which(!(names_act %in% c('kper', 'irch')))]
    if(i > 1 && check_prev(rch$kper, i)) {
      inrech <- -1
    } else {
      inrech <- length(names_act)
    }
    # inirch
    inirch <- 0
    if(rch$nrchop == 2) {
      irch_act <-  rch$kper$irch[i]
      if(!is.na(irch_act)) {
        if(i > 1 && identical(irch_act, rch$kper$irch[i-1])) {
          inirch <- -1
        } else {
          inirch <- length(irch_act)
        }
      } 
    }
    
    if(rch$np > 0) {
      parm_names_active <- parm_names[parm_names %in% names_act]
      np <- length(parm_names_active)
    } else {
      np <- 0
    }
    
    rmfi_write_variables(inrech, ifelse(rch$nrchop == 2, inirch, ''), file=file, integer = TRUE, ...)
    
    # data set 6
    if(np == 0 && inrech >= 0) rmfi_write_array(rch$recharge[[names_act]], file = file, iprn = iprn, ...)
    
    # data set 7
    if(np > 0){
      for(j in 1:np){
        rmfi_write_variables(parm_names_active[j], ifelse(tv_parm[j], rch$kper[i,parm_names_active[j]], ''), as.integer(ifelse(length(rch$irchpf) == 1, rch$irchpf, rch$irchpf[j])), file=file)
      }
    }
    
    # data set 8
    if(rch$nrchop == 2 && inirch >= 0) {
      rmfi_write_array(rch$irch[[irch_act]], file = file, iprn = iprn, ...)
    }
  }
  
}