# Rflow package - functions for the BAS MODFLOW package
# read and write BAS package ----------------------------------------------
#' Read basic (BAS) package file
#'
#' Reads information from a MODFLOW-2000 BAS package file. The BAS
#' package contains information about the active region, the head assigned
#' to no-flow cells and starting heads.
#'
#' @param filename
#' character string;
#' file to be read
#' @param dis
#' character string or object of class DIS.MFpackage;
#' a DIS package file (as file path or object), from which model dimensions
#' will be read
#' @param MF96
#' logical [1];
#' whether to expect an MF96 or MODFLOW-SURFACT formatted BAS file
#'
#' @return
#' object of class BAS.MFpackage:\cr
#' \code{$Options}: chr \code{[]}; options such as SHOWPROGRESS, FREE\cr
#' \code{$IBOUND} : num \code{[NCOL, NROW, NLAY]} or \code{[NLAY]}; flags for
#' inactive (0), active (1) and fixed head cells (-1)\cr
#' \code{$HNOFLO} : num \code{[1]}; head assigned to no flow cells\cr
#' \code{$STRT} : num \code{[NCOL, NROW, NLAY]} or \code{[NLAY]}; initial
#' heads\cr
#'
#' @export
#'
#' @examples
#'
read.BAS <- function(filename, dis, MF96 = FALSE){
if(!MF96) dis <- switch(class(dis),
character = read.DIS(dis),
DIS.MFpackage = dis,
stop("invalid dis to read.BAS"))
txt <- readLines(filename)
# remove comment/ title lines
if(MF96) txt <- txt[-(1:2)] else txt <- txt[!is.comment(txt, "#")]
# options
Options <- if(!MF96) scan(text = txt[1L], what = "character", quiet = TRUE)
# extent (only for MF96 - otherwise this information is in the DIS package)
if(MF96){
extent <- as.integer(read.fws(txt[1L], 10L))
names(extent) <- c("NLAY", "NROW", "NCOL", "NPER", "ITMUNI", rep("", length(extent) - 5L))
}else extent <- dis$extent
# unit numbers (only for MF96 - otherwise this information is in the NAM file)
IUNIT <- if(MF96) as.integer(read.fws(txt[2L], 3L))
# IBOUND array
ln <- if(MF96) 4L else 2L
IBOUND <- replicate(extent["NLAY"], {
nIB <- expected.lines.RIARRAY(txt[ln], extent["NCOL"], extent["NROW"], 1L)
lns <- seq(ln, ln + nIB, 1L)
ar <- interpret.RIARRAY(txt[lns], T, extent["NCOL"], 1L)
ln <<- ln + nIB + 1L
ar
}, simplify = FALSE)
lgs <- lengths(IBOUND)
if(all(lgs == 1L)){
IBOUND <- unname(c(IBOUND, recursive = TRUE))
}else{
ar <- array(dim = extent[c("NCOL", "NROW", "NLAY")])
for(l in 1:length(IBOUND)) ar[,, l] <- IBOUND[[l]]
IBOUND <- ar; rm(ar)
}
# no flow cell head label
# - MODFLOW-SURFACT BAS files also have an entry for Datum; read.BAS ignores this
HNOFLO <- as.numeric(substr(txt[ln], 1L, 10L))
# starting head
ln <- ln + 1L
STRT <- replicate(extent["NLAY"], {
nST <- expected.lines.RIARRAY(txt[ln], extent["NCOL"], extent["NROW"], 1L)
lns <- ln:(ln + nST)
ar <- interpret.RIARRAY(txt[lns], T, extent["NCOL"], 1L)
ln <<- ln + nST + 1L
ar
}, simplify = FALSE)
lgs <- lengths(STRT)
if(all(lgs == 1L)){
STRT <- unname(c(STRT, recursive = TRUE))
}else{
ar <- array(dim = extent[c("NCOL", "NROW", "NLAY")])
for(l in 1:length(STRT)) ar[,, l] <- STRT[[l]]
STRT <- ar; rm(ar)
}
structure(mget(if(MF96)
c("extent", "IUNIT", "IBOUND", "HNOFLO", "STRT") else c("Options", "IBOUND", "HNOFLO", "STRT")),
class = "BAS.MFpackage")
}
#' Write MODFLOW BAS package file
#'
#' Writes information from a BAS.MFpackage list object to a MODFLOW-readable
#' BAS package file.
#'
#' @param BAS
#' object of class BAS.MFpackage, as would be returned by
#' \code{\link{read.BAS}}
#' @param filename
#' character string;
#' file to write to, ideally ending in \code{".bas"}
#' @param unit
#' integer \code{[1]};
#' unit number that the BAS file will be assigned to
#' @param title
#' character string;
#' optional title to prefix to the file, for your reference (not used if MF2k =
#' FALSE)
#' @param MF2k
#' logical;
#' should the file be for MODFLOW 2000, or MODFLOW 96 (in which case stress
#' period information is written from sps)
#' @param sps
#' data.frame or object of class DIS.MFpackage;
#' only required if \code{MF2k = FALSE}; a data frame of stress period
#' information with the same columns as described in \code{\link{read.DIS}},
#' except that the \code{TR} column is not required (ignored if provided)
#' @param NCOL,NROW,NLAY
#' integer \code{[1]}; only required if \code{MF2k = FALSE} and
#' \code{class(sps) != "DIS.MFpackage"}
#'
#' @return \code{NULL}
#'
#' @import data.table
#' @export
#'
#' @examples
#'
write.BAS <- function(BAS, filename, unit = 1L, title = "", MF2k = TRUE, sps,
NCOL, NROW, NLAY){
stopifnot(is(BAS, "BAS.MFpackage"))
if(!MF2k){
if(class(sps) == "DIS.MFpackage"){
NCOL <- sps$extent[c("NCOL")]
NROW <- sps$extent[c("NROW")]
NLAY <- sps$extent[c("NLAY")]
}
sps <- switch(class(sps)[1L],
data.frame = sps,
data.table = sps,
DIS.MFpackage = sps$sps,
stop("write.BAS: invalid sps"))
setDF(sps)
sps <- sps[, c("PERLEN", "NSTP", "TSMULT")]
}
con <- file(filename, "wt")
on.exit(close(con))
writeLines(paste("#", {
substr(gsub("\n", "", title), 1L, if(MF2k) 1000L else 78L)
}), con)
writeLines(c(if(MF2k){
"# MODFLOW 2000 basic package file created by write.BAS function, in R"
}else{
"# MODFLOW 1996 basic package file created by write.BAS function, in R"
}), con)
# MF96 item 3
# - NLAY, NROW, NCOL, NPER, ITMUNI
if(!MF2k){
writeLines(paste0(formatC(NLAY, 0L, 10L, "d"),
formatC(NROW, 0L, 10L, "d"),
formatC(NCOL, 0L, 10L, "d"),
formatC(nrow(sps), 0L, 10L, "d"),
formatC(0L, 0L, 10L, "d")), con)
}
# MF2k item 1
# MF96 item 4
# - Options
if(!all(BAS$Options %in% c("XSECTION", "CHTOCH",
"FREE", "SHOWPROGRESS"))){
warning("write.BAS: invalid BAS options are being written to BAS file")
}
writeLines(paste(BAS$Options, collapse = " "), con)
# MF96 item 5
# - IAPART, ISTRT
if(!MF2k){
writeLines(" 1 1", con)
}
# MF2k item 2
# MF96 item 6
# - IBOUND array
writeLines(if(is.vector(BAS$IBOUND)){
RIARRAY.splitlayers(length(BAS$IBOUND), CNSTNT = BAS$IBOUND,
compress = TRUE, FMTIN_type = "I", FMTIN_w = 3L,
flag.no = unit)
}else{
RIARRAY.splitlayers(dim(BAS$IBOUND)[3L], arr = BAS$IBOUND,
compress = TRUE, FMTIN_type = "I", FMTIN_w = 3L,
flag.no = unit)
}, con)
# MF2k item 3
# MF96 item 7
# - HNOFLO
writeLines(formatC(BAS$HNOFLO, 2L, 10L, "e"), con)
# MF2k item 4
# MF96 item 8
# - STRT/ Shead
writeLines(if(is.vector(BAS$STRT)){
RIARRAY.splitlayers(length(BAS$STRT), CNSTNT = BAS$STRT,
compress = TRUE, FMTIN_type = "e", FMTIN_w = 12L,
FMTIN_d = 4L, flag.no = unit)
}else{
RIARRAY.splitlayers(dim(BAS$STRT)[3L], arr = BAS$STRT,
compress = TRUE, FMTIN_type = "e", FMTIN_w = 12L,
FMTIN_d = 4L, flag.no = unit)
}, con)
# MF96 item 9
# - stress period information
# - PERLEN, NSTP, TSMULT
if(!MF2k){
setDT(sps)
on.exit(setDF(sps), add = TRUE)
ff.widths <- c(10L, 10L, 10L)
ff.digit <- c(2L, 0L, 2L)
ff.format <- c("e", "d", "e")
sps[, {
ffs <- matrix("", .N, 4L)
ffs[, 1:3] <- vapply(1:3,
function(col){
formatC(.SD[[col]], width = 10L,
digits = ff.digit[col],
format = ff.format[col])
}, character(.N))
writeLines(apply(ffs, 1L, paste, collapse = ""), con)
}]
}
invisible()
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.