R/package-huf.R
In rogiersbart/RMODFLOW: Pre- and post-processing of MODFLOW files
Defines functions
rmf_write_lvda
rmf_read_lvda
rmf_create_lvda
rmf_write_kdep
rmf_read_kdep
rmf_create_kdep
rmf_write_huf
rmf_read_huf
rmf_create_huf
Documented in
rmf_create_huf
rmf_create_kdep
rmf_create_lvda
rmf_read_huf
rmf_read_kdep
rmf_read_lvda
rmf_write_huf
rmf_write_kdep
rmf_write_lvda
#' Create an \code{RMODFLOW} huf object
#'
#' \code{rmf_create_huf} creates an \code{RMODFLOW} huf object.
#'
#' @param dis RMODFLOW dis object
#' @param ihufcb flag and unit number for writing cell-by-cell flow terms; defaults to 0
#' @param hdry head assigned to cells that are converted to dry cells; defaults to -888
#' @param nhuf number of hydrogeological units (hgu's)
#' @param iohufheads flag and unit number indicating whether interpolated heads should be written. Defaults to 0 (interpolated heads are not written). See details.
#' @param iohufflows flag and unit number indicating whether interpolated cell-by-cell flows should be written. Defaults to 0 (interpolated flows are not written). See details.
#' @param lthuf vector length \code{dis$nlay}, specifying layer type; defaults to all confined (0) except the first layer (1)
#' @param laywt vector length \code{dis$nlay},, indicating if wetting is active; defaults to 0 for each layer
#' @param wetfct is a factor that is included in the calculation of the head that is initially established at a cell when it is converted from dry to wet; defaults to 0.1
#' @param iwetit is the iteration interval for attempting to wet cells; defaults to 1
#' @param ihdwet is a flag that determines which equation is used to define the initial head at cells that become wet; defaults to 0
#' @param wetdry 3d array with a wetting threshold and flag indicating which neighboring cells can cause a cell to become wet; defaults to NULL. If not read for a specific layer, set all values in that layer to NA.
#' @param hgunam character vector of length nhuf specifying the names of the hgu's
#' @param top 3d array with the top elevations of each hgu (hgu is the third dimension)
#' @param thck 3d array with the thicknesses of each hgu (hgu is the third dimension)
#' @param hguhani vector of flags indicating whether HANI parameters (0) or hguhani (> 0) is used to specify horizontal anisotropies for each hgu; defaults to 1 for each hgu
#' @param hguvani vector of flags indicating whether vk parameters will be specified (0) or HGUVANI is used for vertical anisotropy (> 0). VANI parameters may overwrite the the latter value; defaults to 0 for each hgu
#' @param parameters list of \code{rmf_parameter} as created by \code{\link{rmf_create_parameter}}. See details; defaults to NULL
#' @param hk (named) vector with hydraulic conductivity along rows; one value per hgu specified in the order of \code{hgunam}; defaults to 0.0001. See details.
#' @param hani (named) vector with the ratio of hydraulic conductivity along columns to that along rows; one value per hgu specified in the order of \code{hgunam}; defaults to 1. See details.
#' @param vka (named) vector with vertical hydraulic conductivity or the ratio of horizontal to vertical conductivity depending on the value of hguvani; one value per hgu specified in the order of \code{hgunam}; defaults to NULL. See details.
#' @param ss (named) vector with specific storage; only required when there are transient stress periods; one value per hgu specified in the order of \code{hgunam}; defaults to 1e-5. See details.
#' @param sy (named) vector with specific yield; only required when there are transient stress periods; one value per hgu specified in the order of \code{hgunam}; defaults to 0.15. See details.
#' @param sytp numeric value defining sytp; only required when there are transient stress periods and all layers are confined; defaults to NULL. See details.
#' @details If iohufheads (iohufflows) > 0, heads (cell-by-cell flows) will be interpolated to each hgu and saved, resulting in a single value per hgu at every row-column combination.
#' Flow variables are any of \code{HK, HANI, VK, VANI, SS, SY and SYTP}. Variables can be specified either through the use of parameters or direct vector input through hk, hani, vk, ss, sy and sytp arguments.
#' Specifying directly with parameters is useful if the parametrization is complex. If only a single value per hgu per partyp is needed, the direct vector input can be used.
#' If any parameters are specified, all flow variables have to be specified by parameters. If hk, hani, vk, ss, sy or sytp are not named, the hgu number will be added as a suffix to the type (e.g. 'HK_1').
#' Flow variable SYTP can not be specified per hgu, only for all hgu's. Using direct vector input, VANI can not be specified, only VK. VANI can be set through hguvani.
#'
#' \code{\link{rmf_convert_huf_to_grid}} can be used to convert parameters defined on the HUF grid to the numerical grid
#'
#' @return Object of class huf
#' @export
#' @seealso \code{\link{rmf_create_parameter}}, \code{\link{rmf_convert_huf_to_grid}}, \code{\link{rmf_read_huf}}, \code{\link{rmf_write_huf}} and \url{http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?huf.htm}
rmf_create_huf <- function(dis,
ihufcb = 0,
hdry = -888,
nhuf,
iohufheads = 0,
iohufflows = 0,
lthuf = rmfi_ifelse0(dis$nlay == 1, 1, c(1,rep(0, dis$nlay - 1))),
laywt = rep(0, dis$nlay),
wetfct = 0.1,
iwetit = 1,
ihdwet = 0,
wetdry = NULL,
hgunam,
top,
thck,
hguhani = rep(1.0, nhuf),
hguvani = rep(0, nhuf),
parameters = NULL,
hk = rep(0.0001, nhuf),
hani = rep(1, nhuf),
vk = hk,
ss = rep(1e-5, nhuf),
sy = rep(0.15, nhuf),
sytp = NULL) {
huf <- list()
# data set 0
# to provide comments, use ?comment on the resulting huf object
# data set 1
huf$ihufcb <- ihufcb
huf$hdry <- hdry
huf$nhuf <- nhuf
huf$nphuf <- 0
huf$iohufheads <- iohufheads
huf$iohufflows <- iohufflows
# data set 2
huf$lthuf <- rmfi_ifelse0(length(lthuf) == 1, rep(lthuf, dis$nlay), lthuf)
# data set 3
huf$laywt <- rmfi_ifelse0(length(laywt) == 1, rep(laywt, dis$nlay), laywt)
# data set 4
if(any(huf$laywet != 0)) {
huf$wetfct <- wetfct
huf$iwetit <- iwetit
huf$ihdwet <- ihdwet
# data set 5
huf$wetdry <- wetdry
}
# data set 6
huf$hgunam <- hgunam
# data set 7
huf$top <- top
# data set 8
huf$thck <- thck
# data set 9
huf$hguhani <- rmfi_ifelse0(length(hguhani) == 1, rep(hguhani, huf$nhuf), hguhani)
huf$hguvani <- rmfi_ifelse0(length(hguvani) == 1, rep(hguvani, huf$nhuf), hguvani)
# data set 10-11
if(!is.null(parameters)) {
huf$nphuf <- length(parameters)
# error check
if(any(vapply(parameters, function(i) is.null(attr(i, 'partyp')) || is.null(attr(i, 'hgunam')) || is.null(attr(i, 'parnam')) || is.null(attr(i, 'parval')), TRUE))) {
stop('Please make sure all parameters have a parnam, parval, partyp and hgunam attribute', call. = FALSE)
}
huf$parameters <- list()
huf$parameter_values <- NULL
for(i in 1:huf$nphuf) {
attrb <- attributes(parameters[[i]])
parnam <- attrb$parnam
huf$parameter_values[parnam] <- attrb$parval
huf$parameters[[parnam]] <- parameters[[i]]
}
} else {
# simple parameters
create_huf_parameter <- function(value, index, partyp, hgunam){
if(is.null(names(value))) names(value) <- paste(partyp, index, sep = '_')
rmf_create_parameter(dis = dis, parnam = names(value), parval = value, partyp = partyp, hgunam = hgunam)
}
# hk
if(length(hk) == 1) hk <- rep(hk, huf$nhuf)
hk_p <- lapply(1:huf$nhuf, function(i) create_huf_parameter(value = hk[i], index = i, partyp = 'HK', hgunam = huf$hgunam[i]))
huf$parameters <- hk_p
# hani
if(any(huf$hguhani == 0)) {
if(length(hani) == 1) hani <- rep(hani, huf$nhuf)
hani_p <- lapply(1:huf$nhuf, function(i) rmfi_ifelse0(huf$hguhani[i] == 0, create_huf_parameter(value = hani[i], index = i, partyp = 'HANI', hgunam = huf$hgunam[i]), NULL))
huf$parameters <- c(huf$parameters, hani_p[huf$hguhani == 0])
}
# vk
if(any(huf$hguvani == 0)) {
if(length(vk) == 1) vk <- rep(vk, huf$nhuf)
vk_p <- lapply(1:huf$nhuf, function(i) rmfi_ifelse0(huf$hguvani[i] == 0, create_huf_parameter(value = vk[i], index = i, partyp = 'VK', hgunam = huf$hgunam[i]), NULL))
huf$parameters <- c(huf$parameters, vk_p[huf$hguvani == 0])
}
# ss
if('TR' %in% dis$sstr) {
if(length(ss) == 1) ss <- rep(ss, huf$nhuf)
ss_p <- lapply(1:huf$nhuf, function(i) create_huf_parameter(value = ss[i], index = i, partyp = 'SS', hgunam = huf$hgunam[i]))
huf$parameters <- c(huf$parameters, ss_p)
}
# sy
if('TR' %in% dis$sstr && any(huf$lthuf != 0)) {
if(length(sy) == 1) sy <- rep(sy, huf$nhuf)
sy_p <- lapply(1:huf$nhuf, function(i) rmfi_ifelse0(huf$lthuf[i] != 0, create_huf_parameter(value = sy[i], index = i, partyp = 'SY', hgunam = huf$hgunam[i]), NULL))
huf$parameters <- c(huf$parameters, sy_p[huf$lthuf != 0])
}
# sytp
if('TR' %in% dis$sstr && all(huf$lthuf == 0) && !is.null(sytp)) {
if(is.null(names(sytp))) names(sytp) <- 'SYTP'
sytp_p <- create_huf_parameter(value = sytp, index = 0, partyp = 'SYTP', hgunam = 'SYTP')
huf$parameters <- c(huf$parameters, sytp_p)
}
huf$nphuf <- length(huf$parameters)
huf$parameter_values <- vapply(huf$parameters, function(i) attr(i, 'parval'), 1.0)
names(huf$parameter_values) <- vapply(huf$parameters, function(i) attr(i, 'parnam'), "name")
}
# data set 12
# not supported yet
class(huf) <- c('huf','rmf_package')
return(huf)
}
#' Read a MODFLOW hydrogeologic unit flow file
#'
#' \code{read_huf} reads in a MODFLOW hydrogeologic unit flow file and returns it as an \code{\link{RMODFLOW}} huf object.
#'
#' @param file filename; typically '*.huf'
#' @param dis discretization file object; defaults to that with the same filename but with extension '.dis'
#' @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.
#'
#' @details \code{\link{rmf_convert_huf_to_grid}} can be used to convert parameters defined on the HUF grid to the numerical grid
#' @return object of class huf
#' @export
#' @seealso \code{\link{rmf_convert_huf_to_grid}}, \code{\link{rmf_create_huf}}, \code{\link{rmf_write_huf}} and \url{http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?huf.htm}
rmf_read_huf <- function(file = {cat('Please select huf file ...\n'); file.choose()},
dis = {cat('Please select corresponding dis file ...\n'); rmf_read_dis(file.choose())},
mlt = NULL,
zon = NULL,
...) {
huf_lines <- readr::read_lines(file, lazy = FALSE)
huf <- list()
# data set 0
data_set_0 <- rmfi_parse_comments(huf_lines)
comment(huf) <- data_set_0$comments
huf_lines <- data_set_0$remaining_lines
rm(data_set_0)
# data set 1
data_set_1 <- rmfi_parse_variables(huf_lines)
huf$ihufcb <- as.numeric(data_set_1$variables[1])
huf$hdry <- as.numeric(data_set_1$variables[2])
huf$nhuf <- as.numeric(data_set_1$variables[3])
huf$nphuf <- as.numeric(data_set_1$variables[4])
huf$iohufheads <- ifelse(is.na(data_set_1$variables[5]) || any(is.na(suppressWarnings(as.numeric(data_set_1$variables[5])))) , 0, as.numeric(data_set_1$variables[5]))
huf$iohufflows <- ifelse(is.na(data_set_1$variables[6]) || any(is.na(suppressWarnings(as.numeric(data_set_1$variables[5:6])))), 0, as.numeric(data_set_1$variables[6]))
huf_lines <- data_set_1$remaining_lines
rm(data_set_1)
# data set 2
data_set_2 <- rmfi_parse_variables(huf_lines, nlay = dis$nlay)
huf$lthuf <- as.numeric(data_set_2$variables[1:dis$nlay])
huf_lines <- data_set_2$remaining_lines
rm(data_set_2)
# data set 3
data_set_3 <- rmfi_parse_variables(huf_lines, nlay = dis$nlay)
huf$laywt <- as.numeric(data_set_3$variables[1:dis$nlay])
huf_lines <- data_set_3$remaining_lines
rm(data_set_3)
# data set 4
if(any(huf$laywt > 0)) {
data_set_4 <- rmfi_parse_variables(huf_lines)
huf$wetfct <- as.numeric(data_set_4$variables[1])
huf$iwetit <- as.numeric(data_set_4$variables[2])
huf$ihdwet <- as.numeric(data_set_4$variables[3])
huf_lines <- data_set_4$remaining_lines
rm(data_set_4)
# data set 5
data_set_5 <- rmfi_parse_array(huf_lines,dis$nrow,dis$ncol,sum(which(huf$laywt!=0)), ndim = 3, file = file, ...)
huf$wetdry <- rmf_create_array(dim = c(dis$nrow, dis$ncol, dis$nlay))
huf$wetdry[,,which(huf$laywt != 0)] <- data_set_5$array
huf_lines <- data_set_5$remaining_lines
rm(data_set_5)
}
# data set 6-8
huf$hgunam <- vector(mode='character',length=huf$nhuf)
huf$top <- rmf_create_array(dim=c(dis$nrow, dis$ncol, huf$nhuf))
huf$thck <- rmf_create_array(dim=c(dis$nrow, dis$ncol, huf$nhuf))
for(i in 1:huf$nhuf) {
data_set_6 <- rmfi_parse_variables(huf_lines, character = TRUE)
huf$hgunam[i] <- as.character(data_set_6$variables[1])
huf_lines <- data_set_6$remaining_lines
rm(data_set_6)
data_set_7 <- rmfi_parse_array(huf_lines,dis$nrow,dis$ncol, 1, ndim = 2, file = file, ...)
huf$top[,,i] <- data_set_7$array
huf_lines <- data_set_7$remaining_lines
rm(data_set_7)
data_set_8 <- rmfi_parse_array(huf_lines,dis$nrow,dis$ncol, 1, ndim = 2, file = file, ...)
huf$thck[,,i] <- data_set_8$array
huf_lines <- data_set_8$remaining_lines
rm(data_set_8)
}
# data set 9
huf$hguhani <- vector(mode='numeric',length=huf$nhuf)
huf$hguvani <- vector(mode='numeric',length=huf$nhuf)
data_set_9 <- rmfi_parse_variables(huf_lines, character = TRUE)
if(toupper(data_set_9$variables[1]) == 'ALL') {
huf$hguhani <- rep(as.numeric(data_set_9$variables[2]),huf$nhuf)
huf$hguvani <- rep(as.numeric(data_set_9$variables[3]),huf$nhuf)
huf_lines <- data_set_9$remaining_lines
} else {
k <- which(huf$hgunam == data_set_9$variables[1])
huf$hguhani[k] <- as.numeric(data_set_9$variables[2])
huf$hguvani[k] <- as.numeric(data_set_9$variables[3])
huf_lines <- data_set_9$remaining_lines
for(i in 2:huf$nhuf) {
data_set_9 <- rmfi_parse_variables(huf_lines, character = TRUE)
k <- which(huf$hgunam == data_set_9$variables[1])
huf$hguhani[k] <- as.numeric(data_set_9$variables[2])
huf$hguvani[k] <- as.numeric(data_set_9$variables[3])
huf_lines <- data_set_9$remaining_lines
}
}
rm(data_set_9)
# data set 10-11
types <- NULL
if(huf$nphuf > 0) {
huf$parameters <- list()
huf$parameter_values <- NULL
for(i in 1:huf$nphuf) {
data_set_10 <- rmfi_parse_variables(huf_lines, character = TRUE)
parnam <- data_set_10$variables[1]
partyp <- data_set_10$variables[2]
parval <- as.numeric(data_set_10$variables[3])
nclu <- as.numeric(data_set_10$variables[4])
huf_lines <- data_set_10$remaining_lines
rm(data_set_10)
ds11 <- list(layer = NULL, mltarr = NULL, zonarr = NULL, iz = list())
for(j in 1:nclu) {
data_set_11 <- rmfi_parse_variables(huf_lines, character = TRUE)
ds11$hgunam[j] <- data_set_11$variables[1]
ds11$mltarr[j] <- data_set_11$variables[2]
ds11$zonarr[j] <- data_set_11$variables[3]
# zero or character entry terminates IZ
if(ds11$zonarr[j] == 'ALL') {
ds11$iz[[j]] <- NULL
} else {
iz <- suppressWarnings(as.numeric(data_set_11$variables[4:length(data_set_11$variables)]))
ds11$iz[[j]] <- iz[1:min(length(iz), which(is.na(iz))[1] - 1, which(iz == 0)[1] - 1, na.rm = TRUE)]
}
huf_lines <- data_set_11$remaining_lines
if(toupper(data_set_11$variables[2]) != 'NONE') {
if(is.null(mlt)) stop('Please provide a mlt object', call. = FALSE)
}
if(toupper(data_set_11$variables[3]) != 'ALL') {
if(is.null(zon)) stop('Please provide a zon object', call. = FALSE)
}
rm(data_set_11)
}
huf$parameter_values[parnam] <- parval
huf$parameters[[parnam]] <- rmf_create_parameter(dis = dis, parnam = parnam, partyp = partyp, parval = parval, hgunam = ds11$hgunam, mltnam = ds11$mltarr, zonnam = ds11$zonarr, iz = ds11$iz, mlt = mlt, zon = zon)
# if(is.null(huf[[tolower(partyp)]])) huf[[tolower(partyp)]] <- rmf_create_array(NA, dim = c(dis$nrow, dis$ncol, dis$nlay))
# huf[[tolower(partyp)]][,,unique(ds11$layer)] <- c(huf$parameters[[parnam]])
#
# if(!(toupper(partyp) %in% types)) types <- append(types, toupper(partyp))
}
}
# data set 12
# These are print options, not implemented yet...
class(huf) <- c('huf','rmf_package')
return(huf)
}
#' Write a MODFLOW hydrogeologic unit flow file
#'
#' @param huf an \code{\link{RMODFLOW}} huf object
#' @param file filename to write to; typically '*.huf'
#' @param iprn format code for printing arrays in the listing file; defaults to -1 (no printing)
#' @param ... arguments passed to \code{rmfi_write_array}. Can be ignored when arrays are INTERNAL or CONSTANT.
#' @return \code{NULL}
#' @export
#' @seealso \code{\link{rmf_convert_huf_to_grid}}, \code{\link{rmf_create_huf}}, \code{\link{rmf_read_huf}} and \url{http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?huf.htm}
rmf_write_huf <- function(huf,
file = {cat('Please select huf file to overwrite or provide new filename ...\n'); file.choose()},
iprn=-1,
...) {
# data set 0
v <- packageDescription("RMODFLOW")$Version
cat(paste('# MODFLOW Hydrogeologic Unit Flow Package created by RMODFLOW, version',v,'\n'), file = file)
cat(paste('#', comment(huf)), sep='\n', file=file, append=TRUE)
# data set 1
rmfi_write_variables(as.integer(huf$ihufcb),huf$hdry,as.integer(huf$nhuf),as.integer(huf$nphuf),as.integer(huf$iohufheads),as.integer(huf$iohufflows), file = file)
# data set 2
rmfi_write_variables(huf$lthuf, file=file, integer = TRUE)
# data set 3
rmfi_write_variables(huf$laywt, file=file, integer = TRUE)
# data set 4
if(any(huf$laywt > 0)) {
rmfi_write_variables(huf$wetfct, as.integer(huf$iwetit), as.integer(huf$ihdwet), file = file)
# data set 5
if(dim(huf$wetdry)[3]>0) {
rmfi_write_array(huf$wetdry, file = file, iprn = iprn, ...)
}
}
# data set 6-8
for(i in 1:huf$nhuf) {
rmfi_write_variables(huf$hgunam[i], file=file)
rmfi_write_array(huf$top[,,i], file = file, iprn = iprn, ...)
rmfi_write_array(huf$thck[,,i], file = file, iprn = iprn, ...)
}
# data set 9
for(i in 1:huf$nhuf) {
rmfi_write_variables(huf$hgunam[i],huf$hguhani[i],huf$hguvani[i], file=file)
}
# data set 10-11
for(i in 1:huf$nphuf) {
attrb <- attributes(huf$parameters[[i]])
rmfi_write_variables(attrb$parnam, attrb$partyp, attrb$parval, as.integer(length(attrb$mlt)), file=file)
for(j in 1:length(attrb$mlt)) {
rmfi_write_variables(attrb$hgunam[j], attrb$mlt[j], attrb$zon[j], rmfi_ifelse0(attrb$zon[j] != 'ALL', as.integer(attrb$iz[[j]]), ''), file=file)
}
}
# data set 12
# Print options, not implemented
}
#' Create an \code{RMODFLOW} kdep object
#'
#' \code{rmf_create_kdep} creates an \code{RMODFLOW} kdep object.
#'
#' @param parameters either a single \code{rmf_parameter} or a list of \code{rmf_parameters} specifying the depth-dependency coefficients. See details.
#' @param rs optional 2d array specifying the reference elevation surface
#' @param dis \code{RMODFLOW} dis object. Only used if \code{rs} is supplied.
#'
#' @return a \code{RMODFLOW} kdep object
#' @export
#' @details All parameters should have a hgunam attribute and their partyp attribute set to 'KDEP'.
#' The KDEP package can only be used in conjunction with the HUF package.
#' Note that the parameters are defined on the HUF grid, not the numerical grid
#' \code{\link{rmf_convert_huf_to_grid}} can be used to convert parameters defined on the HUF grid to the numerical grid
#' @seealso \code{link{rmf_convert_huf_to_grid}}, \code{\link{rmf_read_kdep}}, \code{\link{rmf_write_kdep}} and \url{http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?kdep.htm}
rmf_create_kdep <- function(parameters,
rs = NULL,
dis) {
if(length(parameters) == 1 && inherits(parameters[[1]], 'list') && !(inherits(parameters[[1]], 'rmf_parameter'))) parameters <- parameters[[1]]
kdep <- list()
# data set 0
# to provide comments, use ?comment on the resulting kdep object
# data set 1
kdep$npkdep <- length(parameters)
kdep$ifkdep <- ifelse(is.null(rs), 0, 1)
# data set 2
if(!is.null(rs)) kdep$rs <- rmf_create_array(rs, dim = c(dis$nrow, dis$ncol))
# data set 3 & 4
# error check
if(any(vapply(parameters, function(i) is.null(attr(i, 'partyp')) || is.null(attr(i, 'hgunam')) || is.null(attr(i, 'parnam')) || is.null(attr(i, 'parval')), TRUE))) {
stop('Please make sure all parameters have a parnam, parval, partyp and hgunam attribute', call. = FALSE)
}
kdep$parameters <- list()
kdep$parameter_values <- NULL
for(i in 1:kdep$npkdep) {
attrb <- attributes(parameters[[i]])
parnam <- attrb$parnam
kdep$parameter_values[parnam] <- attrb$parval
kdep$parameters[[parnam]] <- parameters[[i]]
}
class(kdep) <- c('kdep', 'rmf_package')
return(kdep)
}
#' Read a MODFLOW hydraulic conductivity depth-dependence capability file
#'
#' \code{rmf_read_kdep} reads in a MODFLOW Hydraulic-Conductivity Depth-Dependence Capability file and returns it as an \code{RMODFLOW} kdep object.
#'
#' @param filename filename; typically *.kdep
#' @param dis \code{RMODFLOW} dis object
#' @param ... arguments passed to \code{rmfi_parse_array}. Can be ignored when input arrays are free-format and INTERNAL or CONSTANT.
#' @return object of class kdep
#' @details Note that the parameters are defined on the HUF grid, not the numerical grid
#' \code{\link{rmf_convert_huf_to_grid}} can be used to convert parameters defined on the huf grid to the numerical grid
#' @export
#' @seealso \code{\link{rmf_convert_huf_to_grid}}, \code{\link{rmf_create_kdep}}, \code{\link{rmf_write_kdep}} and \url{http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?kdep.htm}
rmf_read_kdep <- function(file = {cat('Please select kdep file ...\n'); file.choose()},
dis = {cat('Please select corresponding dis file ...\n'); rmf_read_dis(file.choose())},
...) {
kdep_lines <- readr::read_lines(file, lazy = FALSE)
kdep <- list()
# data set 0
data_set_0 <- rmfi_parse_comments(kdep_lines)
comments <- data_set_0$comments
kdep_lines <- data_set_0$remaining_lines
rm(data_set_0)
# data set 1
data_set_1 <- rmfi_parse_variables(kdep_lines)
kdep$npkdep <- as.numeric(data_set_1$variables[1])
kdep$ifkdep <- as.numeric(data_set_1$variables[2])
kdep_lines <- data_set_1$remaining_lines
rm(data_set_1)
# data set 2
if(kdep$ifkdep > 0) {
data_set2 <- rmfi_parse_array(kdep_lines,dis$nrow,dis$ncol,1, ndim = 2, file = file, ...)
kdep_lines <- data_set2$remaining_lines
kdep$rs <- data_set2$array
rm(data_set2)
}
# data set 3-4
for(i in 1:kdep$npkdep) {
data_set_3 <- rmfi_parse_variables(kdep_lines)
parnam <- data_set_3$variables[1]
partyp <- data_set_3$variables[2]
parval <- as.numeric(data_set_3$variables[3])
nclu <- as.numeric(data_set_3$variables[4])
kdep_lines <- data_set_3$remaining_lines
rm(data_set_3)
ds4 <- list(layer = NULL, mltarr = NULL, zonarr = NULL, iz = list())
for(j in 1:nclu) {
data_set_4 <- rmfi_parse_variables(kdep_lines)
ds4$hgunam[j] <- data_set_4$variables[1]
ds4$mltarr[j] <- data_set_4$variables[2]
ds4$zonarr[j] <- data_set_4$variables[3]
# zero or character entry terminates IZ
if(toupper(ds4$zonarr[j]) == 'ALL') {
ds4$iz[[j]] <- NULL
} else {
iz <- suppressWarnings(as.numeric(data_set_4$variables[4:length(data_set_4$variables)]))
ds4$iz[[j]] <- iz[1:min(length(iz), which(is.na(iz))[1] - 1, which(iz == 0)[1] - 1, na.rm = TRUE)]
}
kdep_lines <- data_set_4$remaining_lines
if(toupper(data_set_4$variables[2]) != 'NONE') {
if(is.null(mlt)) stop('Please provide a mlt object', call. = FALSE)
}
if(toupper(data_set_4$variables[3]) != 'ALL') {
if(is.null(zon)) stop('Please provide a zon object', call. = FALSE)
}
rm(data_set_4)
}
kdep$parameter_values[parnam] <- parval
kdep$parameters[[parnam]] <- rmf_create_parameter(dis = dis, parnam = parnam, partyp = partyp, parval = parval, hgunam = ds4$hgunam, mltnam = ds4$mltarr, zonnam = ds4$zonarr, iz = ds4$iz, mlt = mlt, zon = zon)
# if(is.null(kdep[[tolower(partyp)]])) kdep[[tolower(partyp)]] <- rmf_create_array(NA, dim = c(dis$nrow, dis$ncol, dis$nlay))
# kdep[[tolower(partyp)]][,,unique(ds4$layer)] <- c(kdep$parameters[[parnam]])
#
# if(!(toupper(partyp) %in% types)) types <- append(types, toupper(partyp))
}
comment(kdep) <- comments
class(kdep) <- c('kdep','rmf_package')
return(kdep)
}
#' Write a MODFLOW hydraulic conductivity depth-dependence capability file
#'
#' @param kdep an \code{RMODFLOW} kdep object
#' @param file filename to write to; typically '*.kdep'
#' @param iprn format code for printing arrays in the listing file; defaults to -1 (no printing)
#' @param ... arguments passed to \code{rmfi_write_array}. Can be ignored when arrays are INTERNAL or CONSTANT.
#'
#' @return \code{NULL}
#' @export
#'
#' @seealso \code{\link{rmf_create_kdep}}, \code{\link{rmf_read_kdep}} and \url{http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?kdep.htm}
rmf_write_kdep <- function(kdep,
file = {cat('Please select kdep file to overwrite or provide new filename ...\n'); file.choose()},
iprn = -1,
...) {
# data set 0
v <- packageDescription("RMODFLOW")$Version
cat(paste('# MODFLOW Hydraulic Conductivity Depth-Dependence Capability Package created by RMODFLOW, version',v,'\n'), file = file)
cat(paste('#', comment(kdep)), sep='\n', file=file, append=TRUE)
# data set 1
rmfi_write_variables(kdep$npkdep, kdep$ifkdep, file = file, integer = TRUE)
# data set 2
if(kdep$ifkdep > 0) rmfi_write_array(kdep$rs, file = file, iprn = iprn, ...)
# data set 3-4
for(i in 1:kdep$npkdep) {
attrb <- attributes(kdep$parameters[[i]])
rmfi_write_variables(attrb$parnam, attrb$partyp, attrb$parval, as.integer(length(attrb$mlt)), file=file)
for(j in 1:length(attrb$mlt)) {
rmfi_write_variables(attrb$hgunam[j], attrb$mlt[j], attrb$zon[j], rmfi_ifelse0(attrb$zon[j] != 'ALL', as.integer(attrb$iz[[j]]), ''), file=file)
}
}
}
#' Create an \code{RMODFLOW} lvda object
#'
#' \code{rmf_create_lvda} creates an \code{RMODFLOW} lvda object.
#'
#' @param parameters either a single \code{rmf_parameter} or a list of \code{rmf_parameters} specifying the angle between the grid axis and the principal direction of horizontal hydraulic conductivity. See details.
#' @param dis \code{RMODFLOW} dis object
#'
#' @return a \code{RMODFLOW} lvda object
#' @export
#' @details All parameters should have a layer attribute and their partyp attribute set to 'LVDA'.
#' The LVDA package can only be used in conjunction with the HUF package.
#' Note that the parameters are defined on the numerical grid, not the HUF grid
#' @seealso \code{\link{rmf_read_lvda}}, \code{\link{rmf_write_lvda}} and \url{http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?lvda.htm}
rmf_create_lvda <- function(parameters, dis) {
if(length(parameters) == 1 && inherits(parameters[[1]], 'list') && !(inherits(parameters[[1]], 'rmf_parameter'))) parameters <- parameters[[1]]
lvda <- list()
# data set 0
# to provide comments, use ?comment on the resulting lvda object
# data set 1
lvda$nplvda <- length(parameters)
# data set 2 & 3
# error check
if(any(vapply(parameters, function(i) is.null(attr(i, 'partyp')) || is.null(attr(i, 'layer')) || is.null(attr(i, 'parnam')) || is.null(attr(i, 'parval')), TRUE))) {
stop('Please make sure all parameters have a parnam, parval, partyp and layer attribute', call. = FALSE)
}
lvda$parameters <- list()
lvda$parameter_values <- NULL
lvda$lvda <- rmf_create_array(0, dim = c(dis$nrow, dis$ncol, dis$nlay))
for(i in 1:lvda$nplvda) {
attrb <- attributes(parameters[[i]])
parnam <- attrb$parnam
lvda$parameter_values[parnam] <- attrb$parval
lvda$parameters[[parnam]] <- parameters[[i]]
lvda$lvda[,,unique(attrb$layer)] <- lvda$lvda[,,unique(attrb$layer)] + c(lvda$parameters[[parnam]])
}
class(lvda) <- c('lvda', 'rmf_package')
return(lvda)
}
#' Read a MODFLOW model-layer variable-direction horizontal anisotropy capability file
#'
#' \code{rmf_read_lvda} reads in a MODFLOW Model-Layer Variable-Direction Horizontal Anisotropy Capability file and returns it as an \code{RMODFLOW} lvda object.
#'
#' @param filename filename; typically *.lvda
#' @param dis \code{RMODFLOW} dis object
#' @param ... ignored
#' @details Note that the parameters are defined on the numerical grid, not the HUF grid
#' @return object of class lvda
#' @export
#' @seealso \code{\link{rmf_create_lvda}}, \code{\link{rmf_write_lvda}} and \url{http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?lvda.htm}
rmf_read_lvda <- function(file,
dis = {cat('Please select corresponding dis file ...\n'); rmf_read_dis(file.choose())},
...) {
lvda_lines <- readr::read_lines(file, lazy = FALSE)
lvda <- list()
# data set 0
data_set_0 <- rmfi_parse_comments(lvda_lines)
comments <- data_set_0$comments
lvda_lines <- data_set_0$remaining_lines
rm(data_set_0)
# data set 1
data_set_1 <- rmfi_parse_variables(lvda_lines)
lvda$nplvda <- as.numeric(data_set_1$variables[1])
lvda_lines <- data_set_1$remaining_lines
rm(data_set_1)
lvda$lvda <- rmf_create_array(0, dim = c(dis$nrow, dis$ncol, dis$nlay))
# data set 2-3
for(i in 1:lvda$nplvda) {
data_set_2 <- rmfi_parse_variables(lvda_lines)
parnam <- data_set_2$variables[1]
partyp <- data_set_2$variables[2]
parval <- as.numeric(data_set_2$variables[3])
nclu <- as.numeric(data_set_2$variables[4])
lvda_lines <- data_set_2$remaining_lines
rm(data_set_2)
ds3 <- list(layer = NULL, mltarr = NULL, zonarr = NULL, iz = list())
for(j in 1:nclu) {
data_set_3 <- rmfi_parse_variables(lvda_lines)
ds3$layer[j] <- data_set_3$variables[1]
ds3$mltarr[j] <- data_set_3$variables[2]
ds3$zonarr[j] <- data_set_3$variables[3]
# zero or character entry terminates IZ
if(toupper(ds3$zonarr[j]) == 'ALL') {
ds3$iz[[j]] <- NULL
} else {
iz <- suppressWarnings(as.numeric(data_set_3$variables[4:length(data_set_3$variables)]))
ds3$iz[[j]] <- iz[1:min(length(iz), which(is.na(iz))[1] - 1, which(iz == 0)[1] - 1, na.rm = TRUE)]
}
lvda_lines <- data_set_3$remaining_lines
if(toupper(data_set_3$variables[2]) != 'NONE') {
if(is.null(mlt)) stop('Please provide a mlt object', call. = FALSE)
}
if(toupper(data_set_3$variables[3]) != 'ALL') {
if(is.null(zon)) stop('Please provide a zon object', call. = FALSE)
}
rm(data_set_3)
}
lvda$parameter_values[parnam] <- parval
lvda$parameters[[parnam]] <- rmf_create_parameter(dis = dis, parnam = parnam, partyp = partyp, parval = parval, layer = ds3$layer, mltnam = ds3$mltarr, zonnam = ds3$zonarr, iz = ds3$iz, mlt = mlt, zon = zon)
# if(is.null(lvda[[tolower(partyp)]])) lvda[[tolower(partyp)]] <- rmf_create_array(NA, dim = c(dis$nrow, dis$ncol, dis$nlay))
# lvda[[tolower(partyp)]][,,unique(ds3$layer)] <- c(lvda$parameters[[parnam]])
#
# if(!(toupper(partyp) %in% types)) types <- append(types, toupper(partyp))
lvda$lvda[,,unique(ds3$layer)] <- lvda$lvda[,,unique(ds3$layer)] + c(lvda$parameters[[parnam]])
}
comment(lvda) <- comments
class(lvda) <- c('lvda','rmf_package')
return(lvda)
}
#' Write a MODFLOW model-layer variable-direction horizontal anisotropy capability file
#'
#' @param lvda an \code{RMODFLOW} lvda object
#' @param file filename to write to; typically '*.lvda'
#' @param ... ignored
#'
#' @return \code{NULL}
#' @export
#'
#' @seealso \code{\link{rmf_create_lvda}}, \code{\link{rmf_read_lvda}} and \url{http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?lvda.htm}
rmf_write_lvda <- function(lvda,
file = {cat('Please select lvda file to overwrite or provide new filename ...\n'); file.choose()},
...) {
# data set 0
v <- packageDescription("RMODFLOW")$Version
cat(paste('# MODFLOW Model-Layer Variable-Direction Horizontal Anisotropy Capability Package created by RMODFLOW, version',v,'\n'), file = file)
cat(paste('#', comment(lvda)), sep='\n', file=file, append=TRUE)
# data set 1
rmfi_write_variables(lvda$nplvda, file = file, integer = TRUE)
# data set 3-4
for(i in 1:lvda$nplvda) {
attrb <- attributes(lvda$parameters[[i]])
rmfi_write_variables(attrb$parnam, attrb$partyp, attrb$parval, as.integer(length(attrb$mlt)), file=file)
for(j in 1:length(attrb$mlt)) {
rmfi_write_variables(as.integer(attrb$layer[j]), attrb$mlt[j], attrb$zon[j], rmfi_ifelse0(attrb$zon[j] != 'ALL', as.integer(attrb$iz[[j]]), ''), file=file)
}
}
}
rogiersbart/RMODFLOW documentation built on Jan. 14, 2023, 4:21 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 rmf_write_lvda rmf_read_lvda rmf_create_lvda rmf_write_kdep rmf_read_kdep rmf_create_kdep rmf_write_huf rmf_read_huf rmf_create_huf
Documented in rmf_create_huf rmf_create_kdep rmf_create_lvda rmf_read_huf rmf_read_kdep rmf_read_lvda rmf_write_huf rmf_write_kdep rmf_write_lvda
#' Create an \code{RMODFLOW} huf object
#'
#' \code{rmf_create_huf} creates an \code{RMODFLOW} huf object.
#'
#' @param dis RMODFLOW dis object
#' @param ihufcb flag and unit number for writing cell-by-cell flow terms; defaults to 0
#' @param hdry head assigned to cells that are converted to dry cells; defaults to -888
#' @param nhuf number of hydrogeological units (hgu's)
#' @param iohufheads flag and unit number indicating whether interpolated heads should be written. Defaults to 0 (interpolated heads are not written). See details.
#' @param iohufflows flag and unit number indicating whether interpolated cell-by-cell flows should be written. Defaults to 0 (interpolated flows are not written). See details.
#' @param lthuf vector length \code{dis$nlay}, specifying layer type; defaults to all confined (0) except the first layer (1)
#' @param laywt vector length \code{dis$nlay},, indicating if wetting is active; defaults to 0 for each layer
#' @param wetfct is a factor that is included in the calculation of the head that is initially established at a cell when it is converted from dry to wet; defaults to 0.1
#' @param iwetit is the iteration interval for attempting to wet cells; defaults to 1
#' @param ihdwet is a flag that determines which equation is used to define the initial head at cells that become wet; defaults to 0
#' @param wetdry 3d array with a wetting threshold and flag indicating which neighboring cells can cause a cell to become wet; defaults to NULL. If not read for a specific layer, set all values in that layer to NA.
#' @param hgunam character vector of length nhuf specifying the names of the hgu's
#' @param top 3d array with the top elevations of each hgu (hgu is the third dimension)
#' @param thck 3d array with the thicknesses of each hgu (hgu is the third dimension)
#' @param hguhani vector of flags indicating whether HANI parameters (0) or hguhani (> 0) is used to specify horizontal anisotropies for each hgu; defaults to 1 for each hgu
#' @param hguvani vector of flags indicating whether vk parameters will be specified (0) or HGUVANI is used for vertical anisotropy (> 0). VANI parameters may overwrite the the latter value; defaults to 0 for each hgu
#' @param parameters list of \code{rmf_parameter} as created by \code{\link{rmf_create_parameter}}. See details; defaults to NULL
#' @param hk (named) vector with hydraulic conductivity along rows; one value per hgu specified in the order of \code{hgunam}; defaults to 0.0001. See details.
#' @param hani (named) vector with the ratio of hydraulic conductivity along columns to that along rows; one value per hgu specified in the order of \code{hgunam}; defaults to 1. See details.
#' @param vka (named) vector with vertical hydraulic conductivity or the ratio of horizontal to vertical conductivity depending on the value of hguvani; one value per hgu specified in the order of \code{hgunam}; defaults to NULL. See details.
#' @param ss (named) vector with specific storage; only required when there are transient stress periods; one value per hgu specified in the order of \code{hgunam}; defaults to 1e-5. See details.
#' @param sy (named) vector with specific yield; only required when there are transient stress periods; one value per hgu specified in the order of \code{hgunam}; defaults to 0.15. See details.
#' @param sytp numeric value defining sytp; only required when there are transient stress periods and all layers are confined; defaults to NULL. See details.
#' @details If iohufheads (iohufflows) > 0, heads (cell-by-cell flows) will be interpolated to each hgu and saved, resulting in a single value per hgu at every row-column combination.
#' Flow variables are any of \code{HK, HANI, VK, VANI, SS, SY and SYTP}. Variables can be specified either through the use of parameters or direct vector input through hk, hani, vk, ss, sy and sytp arguments.
#' Specifying directly with parameters is useful if the parametrization is complex. If only a single value per hgu per partyp is needed, the direct vector input can be used.
#' If any parameters are specified, all flow variables have to be specified by parameters. If hk, hani, vk, ss, sy or sytp are not named, the hgu number will be added as a suffix to the type (e.g. 'HK_1').
#' Flow variable SYTP can not be specified per hgu, only for all hgu's. Using direct vector input, VANI can not be specified, only VK. VANI can be set through hguvani.
#'
#' \code{\link{rmf_convert_huf_to_grid}} can be used to convert parameters defined on the HUF grid to the numerical grid
#'
#' @return Object of class huf
#' @export
#' @seealso \code{\link{rmf_create_parameter}}, \code{\link{rmf_convert_huf_to_grid}}, \code{\link{rmf_read_huf}}, \code{\link{rmf_write_huf}} and \url{http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?huf.htm}
rmf_create_huf <- function(dis,
ihufcb = 0,
hdry = -888,
nhuf,
iohufheads = 0,
iohufflows = 0,
lthuf = rmfi_ifelse0(dis$nlay == 1, 1, c(1,rep(0, dis$nlay - 1))),
laywt = rep(0, dis$nlay),
wetfct = 0.1,
iwetit = 1,
ihdwet = 0,
wetdry = NULL,
hgunam,
top,
thck,
hguhani = rep(1.0, nhuf),
hguvani = rep(0, nhuf),
parameters = NULL,
hk = rep(0.0001, nhuf),
hani = rep(1, nhuf),
vk = hk,
ss = rep(1e-5, nhuf),
sy = rep(0.15, nhuf),
sytp = NULL) {
huf <- list()
# data set 0
# to provide comments, use ?comment on the resulting huf object
# data set 1
huf$ihufcb <- ihufcb
huf$hdry <- hdry
huf$nhuf <- nhuf
huf$nphuf <- 0
huf$iohufheads <- iohufheads
huf$iohufflows <- iohufflows
# data set 2
huf$lthuf <- rmfi_ifelse0(length(lthuf) == 1, rep(lthuf, dis$nlay), lthuf)
# data set 3
huf$laywt <- rmfi_ifelse0(length(laywt) == 1, rep(laywt, dis$nlay), laywt)
# data set 4
if(any(huf$laywet != 0)) {
huf$wetfct <- wetfct
huf$iwetit <- iwetit
huf$ihdwet <- ihdwet
# data set 5
huf$wetdry <- wetdry
}
# data set 6
huf$hgunam <- hgunam
# data set 7
huf$top <- top
# data set 8
huf$thck <- thck
# data set 9
huf$hguhani <- rmfi_ifelse0(length(hguhani) == 1, rep(hguhani, huf$nhuf), hguhani)
huf$hguvani <- rmfi_ifelse0(length(hguvani) == 1, rep(hguvani, huf$nhuf), hguvani)
# data set 10-11
if(!is.null(parameters)) {
huf$nphuf <- length(parameters)
# error check
if(any(vapply(parameters, function(i) is.null(attr(i, 'partyp')) || is.null(attr(i, 'hgunam')) || is.null(attr(i, 'parnam')) || is.null(attr(i, 'parval')), TRUE))) {
stop('Please make sure all parameters have a parnam, parval, partyp and hgunam attribute', call. = FALSE)
}
huf$parameters <- list()
huf$parameter_values <- NULL
for(i in 1:huf$nphuf) {
attrb <- attributes(parameters[[i]])
parnam <- attrb$parnam
huf$parameter_values[parnam] <- attrb$parval
huf$parameters[[parnam]] <- parameters[[i]]
}
} else {
# simple parameters
create_huf_parameter <- function(value, index, partyp, hgunam){
if(is.null(names(value))) names(value) <- paste(partyp, index, sep = '_')
rmf_create_parameter(dis = dis, parnam = names(value), parval = value, partyp = partyp, hgunam = hgunam)
}
# hk
if(length(hk) == 1) hk <- rep(hk, huf$nhuf)
hk_p <- lapply(1:huf$nhuf, function(i) create_huf_parameter(value = hk[i], index = i, partyp = 'HK', hgunam = huf$hgunam[i]))
huf$parameters <- hk_p
# hani
if(any(huf$hguhani == 0)) {
if(length(hani) == 1) hani <- rep(hani, huf$nhuf)
hani_p <- lapply(1:huf$nhuf, function(i) rmfi_ifelse0(huf$hguhani[i] == 0, create_huf_parameter(value = hani[i], index = i, partyp = 'HANI', hgunam = huf$hgunam[i]), NULL))
huf$parameters <- c(huf$parameters, hani_p[huf$hguhani == 0])
}
# vk
if(any(huf$hguvani == 0)) {
if(length(vk) == 1) vk <- rep(vk, huf$nhuf)
vk_p <- lapply(1:huf$nhuf, function(i) rmfi_ifelse0(huf$hguvani[i] == 0, create_huf_parameter(value = vk[i], index = i, partyp = 'VK', hgunam = huf$hgunam[i]), NULL))
huf$parameters <- c(huf$parameters, vk_p[huf$hguvani == 0])
}
# ss
if('TR' %in% dis$sstr) {
if(length(ss) == 1) ss <- rep(ss, huf$nhuf)
ss_p <- lapply(1:huf$nhuf, function(i) create_huf_parameter(value = ss[i], index = i, partyp = 'SS', hgunam = huf$hgunam[i]))
huf$parameters <- c(huf$parameters, ss_p)
}
# sy
if('TR' %in% dis$sstr && any(huf$lthuf != 0)) {
if(length(sy) == 1) sy <- rep(sy, huf$nhuf)
sy_p <- lapply(1:huf$nhuf, function(i) rmfi_ifelse0(huf$lthuf[i] != 0, create_huf_parameter(value = sy[i], index = i, partyp = 'SY', hgunam = huf$hgunam[i]), NULL))
huf$parameters <- c(huf$parameters, sy_p[huf$lthuf != 0])
}
# sytp
if('TR' %in% dis$sstr && all(huf$lthuf == 0) && !is.null(sytp)) {
if(is.null(names(sytp))) names(sytp) <- 'SYTP'
sytp_p <- create_huf_parameter(value = sytp, index = 0, partyp = 'SYTP', hgunam = 'SYTP')
huf$parameters <- c(huf$parameters, sytp_p)
}
huf$nphuf <- length(huf$parameters)
huf$parameter_values <- vapply(huf$parameters, function(i) attr(i, 'parval'), 1.0)
names(huf$parameter_values) <- vapply(huf$parameters, function(i) attr(i, 'parnam'), "name")
}
# data set 12
# not supported yet
class(huf) <- c('huf','rmf_package')
return(huf)
}
#' Read a MODFLOW hydrogeologic unit flow file
#'
#' \code{read_huf} reads in a MODFLOW hydrogeologic unit flow file and returns it as an \code{\link{RMODFLOW}} huf object.
#'
#' @param file filename; typically '*.huf'
#' @param dis discretization file object; defaults to that with the same filename but with extension '.dis'
#' @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.
#'
#' @details \code{\link{rmf_convert_huf_to_grid}} can be used to convert parameters defined on the HUF grid to the numerical grid
#' @return object of class huf
#' @export
#' @seealso \code{\link{rmf_convert_huf_to_grid}}, \code{\link{rmf_create_huf}}, \code{\link{rmf_write_huf}} and \url{http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?huf.htm}
rmf_read_huf <- function(file = {cat('Please select huf file ...\n'); file.choose()},
dis = {cat('Please select corresponding dis file ...\n'); rmf_read_dis(file.choose())},
mlt = NULL,
zon = NULL,
...) {
huf_lines <- readr::read_lines(file, lazy = FALSE)
huf <- list()
# data set 0
data_set_0 <- rmfi_parse_comments(huf_lines)
comment(huf) <- data_set_0$comments
huf_lines <- data_set_0$remaining_lines
rm(data_set_0)
# data set 1
data_set_1 <- rmfi_parse_variables(huf_lines)
huf$ihufcb <- as.numeric(data_set_1$variables[1])
huf$hdry <- as.numeric(data_set_1$variables[2])
huf$nhuf <- as.numeric(data_set_1$variables[3])
huf$nphuf <- as.numeric(data_set_1$variables[4])
huf$iohufheads <- ifelse(is.na(data_set_1$variables[5]) || any(is.na(suppressWarnings(as.numeric(data_set_1$variables[5])))) , 0, as.numeric(data_set_1$variables[5]))
huf$iohufflows <- ifelse(is.na(data_set_1$variables[6]) || any(is.na(suppressWarnings(as.numeric(data_set_1$variables[5:6])))), 0, as.numeric(data_set_1$variables[6]))
huf_lines <- data_set_1$remaining_lines
rm(data_set_1)
# data set 2
data_set_2 <- rmfi_parse_variables(huf_lines, nlay = dis$nlay)
huf$lthuf <- as.numeric(data_set_2$variables[1:dis$nlay])
huf_lines <- data_set_2$remaining_lines
rm(data_set_2)
# data set 3
data_set_3 <- rmfi_parse_variables(huf_lines, nlay = dis$nlay)
huf$laywt <- as.numeric(data_set_3$variables[1:dis$nlay])
huf_lines <- data_set_3$remaining_lines
rm(data_set_3)
# data set 4
if(any(huf$laywt > 0)) {
data_set_4 <- rmfi_parse_variables(huf_lines)
huf$wetfct <- as.numeric(data_set_4$variables[1])
huf$iwetit <- as.numeric(data_set_4$variables[2])
huf$ihdwet <- as.numeric(data_set_4$variables[3])
huf_lines <- data_set_4$remaining_lines
rm(data_set_4)
# data set 5
data_set_5 <- rmfi_parse_array(huf_lines,dis$nrow,dis$ncol,sum(which(huf$laywt!=0)), ndim = 3, file = file, ...)
huf$wetdry <- rmf_create_array(dim = c(dis$nrow, dis$ncol, dis$nlay))
huf$wetdry[,,which(huf$laywt != 0)] <- data_set_5$array
huf_lines <- data_set_5$remaining_lines
rm(data_set_5)
}
# data set 6-8
huf$hgunam <- vector(mode='character',length=huf$nhuf)
huf$top <- rmf_create_array(dim=c(dis$nrow, dis$ncol, huf$nhuf))
huf$thck <- rmf_create_array(dim=c(dis$nrow, dis$ncol, huf$nhuf))
for(i in 1:huf$nhuf) {
data_set_6 <- rmfi_parse_variables(huf_lines, character = TRUE)
huf$hgunam[i] <- as.character(data_set_6$variables[1])
huf_lines <- data_set_6$remaining_lines
rm(data_set_6)
data_set_7 <- rmfi_parse_array(huf_lines,dis$nrow,dis$ncol, 1, ndim = 2, file = file, ...)
huf$top[,,i] <- data_set_7$array
huf_lines <- data_set_7$remaining_lines
rm(data_set_7)
data_set_8 <- rmfi_parse_array(huf_lines,dis$nrow,dis$ncol, 1, ndim = 2, file = file, ...)
huf$thck[,,i] <- data_set_8$array
huf_lines <- data_set_8$remaining_lines
rm(data_set_8)
}
# data set 9
huf$hguhani <- vector(mode='numeric',length=huf$nhuf)
huf$hguvani <- vector(mode='numeric',length=huf$nhuf)
data_set_9 <- rmfi_parse_variables(huf_lines, character = TRUE)
if(toupper(data_set_9$variables[1]) == 'ALL') {
huf$hguhani <- rep(as.numeric(data_set_9$variables[2]),huf$nhuf)
huf$hguvani <- rep(as.numeric(data_set_9$variables[3]),huf$nhuf)
huf_lines <- data_set_9$remaining_lines
} else {
k <- which(huf$hgunam == data_set_9$variables[1])
huf$hguhani[k] <- as.numeric(data_set_9$variables[2])
huf$hguvani[k] <- as.numeric(data_set_9$variables[3])
huf_lines <- data_set_9$remaining_lines
for(i in 2:huf$nhuf) {
data_set_9 <- rmfi_parse_variables(huf_lines, character = TRUE)
k <- which(huf$hgunam == data_set_9$variables[1])
huf$hguhani[k] <- as.numeric(data_set_9$variables[2])
huf$hguvani[k] <- as.numeric(data_set_9$variables[3])
huf_lines <- data_set_9$remaining_lines
}
}
rm(data_set_9)
# data set 10-11
types <- NULL
if(huf$nphuf > 0) {
huf$parameters <- list()
huf$parameter_values <- NULL
for(i in 1:huf$nphuf) {
data_set_10 <- rmfi_parse_variables(huf_lines, character = TRUE)
parnam <- data_set_10$variables[1]
partyp <- data_set_10$variables[2]
parval <- as.numeric(data_set_10$variables[3])
nclu <- as.numeric(data_set_10$variables[4])
huf_lines <- data_set_10$remaining_lines
rm(data_set_10)
ds11 <- list(layer = NULL, mltarr = NULL, zonarr = NULL, iz = list())
for(j in 1:nclu) {
data_set_11 <- rmfi_parse_variables(huf_lines, character = TRUE)
ds11$hgunam[j] <- data_set_11$variables[1]
ds11$mltarr[j] <- data_set_11$variables[2]
ds11$zonarr[j] <- data_set_11$variables[3]
# zero or character entry terminates IZ
if(ds11$zonarr[j] == 'ALL') {
ds11$iz[[j]] <- NULL
} else {
iz <- suppressWarnings(as.numeric(data_set_11$variables[4:length(data_set_11$variables)]))
ds11$iz[[j]] <- iz[1:min(length(iz), which(is.na(iz))[1] - 1, which(iz == 0)[1] - 1, na.rm = TRUE)]
}
huf_lines <- data_set_11$remaining_lines
if(toupper(data_set_11$variables[2]) != 'NONE') {
if(is.null(mlt)) stop('Please provide a mlt object', call. = FALSE)
}
if(toupper(data_set_11$variables[3]) != 'ALL') {
if(is.null(zon)) stop('Please provide a zon object', call. = FALSE)
}
rm(data_set_11)
}
huf$parameter_values[parnam] <- parval
huf$parameters[[parnam]] <- rmf_create_parameter(dis = dis, parnam = parnam, partyp = partyp, parval = parval, hgunam = ds11$hgunam, mltnam = ds11$mltarr, zonnam = ds11$zonarr, iz = ds11$iz, mlt = mlt, zon = zon)
# if(is.null(huf[[tolower(partyp)]])) huf[[tolower(partyp)]] <- rmf_create_array(NA, dim = c(dis$nrow, dis$ncol, dis$nlay))
# huf[[tolower(partyp)]][,,unique(ds11$layer)] <- c(huf$parameters[[parnam]])
#
# if(!(toupper(partyp) %in% types)) types <- append(types, toupper(partyp))
}
}
# data set 12
# These are print options, not implemented yet...
class(huf) <- c('huf','rmf_package')
return(huf)
}
#' Write a MODFLOW hydrogeologic unit flow file
#'
#' @param huf an \code{\link{RMODFLOW}} huf object
#' @param file filename to write to; typically '*.huf'
#' @param iprn format code for printing arrays in the listing file; defaults to -1 (no printing)
#' @param ... arguments passed to \code{rmfi_write_array}. Can be ignored when arrays are INTERNAL or CONSTANT.
#' @return \code{NULL}
#' @export
#' @seealso \code{\link{rmf_convert_huf_to_grid}}, \code{\link{rmf_create_huf}}, \code{\link{rmf_read_huf}} and \url{http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?huf.htm}
rmf_write_huf <- function(huf,
file = {cat('Please select huf file to overwrite or provide new filename ...\n'); file.choose()},
iprn=-1,
...) {
# data set 0
v <- packageDescription("RMODFLOW")$Version
cat(paste('# MODFLOW Hydrogeologic Unit Flow Package created by RMODFLOW, version',v,'\n'), file = file)
cat(paste('#', comment(huf)), sep='\n', file=file, append=TRUE)
# data set 1
rmfi_write_variables(as.integer(huf$ihufcb),huf$hdry,as.integer(huf$nhuf),as.integer(huf$nphuf),as.integer(huf$iohufheads),as.integer(huf$iohufflows), file = file)
# data set 2
rmfi_write_variables(huf$lthuf, file=file, integer = TRUE)
# data set 3
rmfi_write_variables(huf$laywt, file=file, integer = TRUE)
# data set 4
if(any(huf$laywt > 0)) {
rmfi_write_variables(huf$wetfct, as.integer(huf$iwetit), as.integer(huf$ihdwet), file = file)
# data set 5
if(dim(huf$wetdry)[3]>0) {
rmfi_write_array(huf$wetdry, file = file, iprn = iprn, ...)
}
}
# data set 6-8
for(i in 1:huf$nhuf) {
rmfi_write_variables(huf$hgunam[i], file=file)
rmfi_write_array(huf$top[,,i], file = file, iprn = iprn, ...)
rmfi_write_array(huf$thck[,,i], file = file, iprn = iprn, ...)
}
# data set 9
for(i in 1:huf$nhuf) {
rmfi_write_variables(huf$hgunam[i],huf$hguhani[i],huf$hguvani[i], file=file)
}
# data set 10-11
for(i in 1:huf$nphuf) {
attrb <- attributes(huf$parameters[[i]])
rmfi_write_variables(attrb$parnam, attrb$partyp, attrb$parval, as.integer(length(attrb$mlt)), file=file)
for(j in 1:length(attrb$mlt)) {
rmfi_write_variables(attrb$hgunam[j], attrb$mlt[j], attrb$zon[j], rmfi_ifelse0(attrb$zon[j] != 'ALL', as.integer(attrb$iz[[j]]), ''), file=file)
}
}
# data set 12
# Print options, not implemented
}
#' Create an \code{RMODFLOW} kdep object
#'
#' \code{rmf_create_kdep} creates an \code{RMODFLOW} kdep object.
#'
#' @param parameters either a single \code{rmf_parameter} or a list of \code{rmf_parameters} specifying the depth-dependency coefficients. See details.
#' @param rs optional 2d array specifying the reference elevation surface
#' @param dis \code{RMODFLOW} dis object. Only used if \code{rs} is supplied.
#'
#' @return a \code{RMODFLOW} kdep object
#' @export
#' @details All parameters should have a hgunam attribute and their partyp attribute set to 'KDEP'.
#' The KDEP package can only be used in conjunction with the HUF package.
#' Note that the parameters are defined on the HUF grid, not the numerical grid
#' \code{\link{rmf_convert_huf_to_grid}} can be used to convert parameters defined on the HUF grid to the numerical grid
#' @seealso \code{link{rmf_convert_huf_to_grid}}, \code{\link{rmf_read_kdep}}, \code{\link{rmf_write_kdep}} and \url{http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?kdep.htm}
rmf_create_kdep <- function(parameters,
rs = NULL,
dis) {
if(length(parameters) == 1 && inherits(parameters[[1]], 'list') && !(inherits(parameters[[1]], 'rmf_parameter'))) parameters <- parameters[[1]]
kdep <- list()
# data set 0
# to provide comments, use ?comment on the resulting kdep object
# data set 1
kdep$npkdep <- length(parameters)
kdep$ifkdep <- ifelse(is.null(rs), 0, 1)
# data set 2
if(!is.null(rs)) kdep$rs <- rmf_create_array(rs, dim = c(dis$nrow, dis$ncol))
# data set 3 & 4
# error check
if(any(vapply(parameters, function(i) is.null(attr(i, 'partyp')) || is.null(attr(i, 'hgunam')) || is.null(attr(i, 'parnam')) || is.null(attr(i, 'parval')), TRUE))) {
stop('Please make sure all parameters have a parnam, parval, partyp and hgunam attribute', call. = FALSE)
}
kdep$parameters <- list()
kdep$parameter_values <- NULL
for(i in 1:kdep$npkdep) {
attrb <- attributes(parameters[[i]])
parnam <- attrb$parnam
kdep$parameter_values[parnam] <- attrb$parval
kdep$parameters[[parnam]] <- parameters[[i]]
}
class(kdep) <- c('kdep', 'rmf_package')
return(kdep)
}
#' Read a MODFLOW hydraulic conductivity depth-dependence capability file
#'
#' \code{rmf_read_kdep} reads in a MODFLOW Hydraulic-Conductivity Depth-Dependence Capability file and returns it as an \code{RMODFLOW} kdep object.
#'
#' @param filename filename; typically *.kdep
#' @param dis \code{RMODFLOW} dis object
#' @param ... arguments passed to \code{rmfi_parse_array}. Can be ignored when input arrays are free-format and INTERNAL or CONSTANT.
#' @return object of class kdep
#' @details Note that the parameters are defined on the HUF grid, not the numerical grid
#' \code{\link{rmf_convert_huf_to_grid}} can be used to convert parameters defined on the huf grid to the numerical grid
#' @export
#' @seealso \code{\link{rmf_convert_huf_to_grid}}, \code{\link{rmf_create_kdep}}, \code{\link{rmf_write_kdep}} and \url{http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?kdep.htm}
rmf_read_kdep <- function(file = {cat('Please select kdep file ...\n'); file.choose()},
dis = {cat('Please select corresponding dis file ...\n'); rmf_read_dis(file.choose())},
...) {
kdep_lines <- readr::read_lines(file, lazy = FALSE)
kdep <- list()
# data set 0
data_set_0 <- rmfi_parse_comments(kdep_lines)
comments <- data_set_0$comments
kdep_lines <- data_set_0$remaining_lines
rm(data_set_0)
# data set 1
data_set_1 <- rmfi_parse_variables(kdep_lines)
kdep$npkdep <- as.numeric(data_set_1$variables[1])
kdep$ifkdep <- as.numeric(data_set_1$variables[2])
kdep_lines <- data_set_1$remaining_lines
rm(data_set_1)
# data set 2
if(kdep$ifkdep > 0) {
data_set2 <- rmfi_parse_array(kdep_lines,dis$nrow,dis$ncol,1, ndim = 2, file = file, ...)
kdep_lines <- data_set2$remaining_lines
kdep$rs <- data_set2$array
rm(data_set2)
}
# data set 3-4
for(i in 1:kdep$npkdep) {
data_set_3 <- rmfi_parse_variables(kdep_lines)
parnam <- data_set_3$variables[1]
partyp <- data_set_3$variables[2]
parval <- as.numeric(data_set_3$variables[3])
nclu <- as.numeric(data_set_3$variables[4])
kdep_lines <- data_set_3$remaining_lines
rm(data_set_3)
ds4 <- list(layer = NULL, mltarr = NULL, zonarr = NULL, iz = list())
for(j in 1:nclu) {
data_set_4 <- rmfi_parse_variables(kdep_lines)
ds4$hgunam[j] <- data_set_4$variables[1]
ds4$mltarr[j] <- data_set_4$variables[2]
ds4$zonarr[j] <- data_set_4$variables[3]
# zero or character entry terminates IZ
if(toupper(ds4$zonarr[j]) == 'ALL') {
ds4$iz[[j]] <- NULL
} else {
iz <- suppressWarnings(as.numeric(data_set_4$variables[4:length(data_set_4$variables)]))
ds4$iz[[j]] <- iz[1:min(length(iz), which(is.na(iz))[1] - 1, which(iz == 0)[1] - 1, na.rm = TRUE)]
}
kdep_lines <- data_set_4$remaining_lines
if(toupper(data_set_4$variables[2]) != 'NONE') {
if(is.null(mlt)) stop('Please provide a mlt object', call. = FALSE)
}
if(toupper(data_set_4$variables[3]) != 'ALL') {
if(is.null(zon)) stop('Please provide a zon object', call. = FALSE)
}
rm(data_set_4)
}
kdep$parameter_values[parnam] <- parval
kdep$parameters[[parnam]] <- rmf_create_parameter(dis = dis, parnam = parnam, partyp = partyp, parval = parval, hgunam = ds4$hgunam, mltnam = ds4$mltarr, zonnam = ds4$zonarr, iz = ds4$iz, mlt = mlt, zon = zon)
# if(is.null(kdep[[tolower(partyp)]])) kdep[[tolower(partyp)]] <- rmf_create_array(NA, dim = c(dis$nrow, dis$ncol, dis$nlay))
# kdep[[tolower(partyp)]][,,unique(ds4$layer)] <- c(kdep$parameters[[parnam]])
#
# if(!(toupper(partyp) %in% types)) types <- append(types, toupper(partyp))
}
comment(kdep) <- comments
class(kdep) <- c('kdep','rmf_package')
return(kdep)
}
#' Write a MODFLOW hydraulic conductivity depth-dependence capability file
#'
#' @param kdep an \code{RMODFLOW} kdep object
#' @param file filename to write to; typically '*.kdep'
#' @param iprn format code for printing arrays in the listing file; defaults to -1 (no printing)
#' @param ... arguments passed to \code{rmfi_write_array}. Can be ignored when arrays are INTERNAL or CONSTANT.
#'
#' @return \code{NULL}
#' @export
#'
#' @seealso \code{\link{rmf_create_kdep}}, \code{\link{rmf_read_kdep}} and \url{http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?kdep.htm}
rmf_write_kdep <- function(kdep,
file = {cat('Please select kdep file to overwrite or provide new filename ...\n'); file.choose()},
iprn = -1,
...) {
# data set 0
v <- packageDescription("RMODFLOW")$Version
cat(paste('# MODFLOW Hydraulic Conductivity Depth-Dependence Capability Package created by RMODFLOW, version',v,'\n'), file = file)
cat(paste('#', comment(kdep)), sep='\n', file=file, append=TRUE)
# data set 1
rmfi_write_variables(kdep$npkdep, kdep$ifkdep, file = file, integer = TRUE)
# data set 2
if(kdep$ifkdep > 0) rmfi_write_array(kdep$rs, file = file, iprn = iprn, ...)
# data set 3-4
for(i in 1:kdep$npkdep) {
attrb <- attributes(kdep$parameters[[i]])
rmfi_write_variables(attrb$parnam, attrb$partyp, attrb$parval, as.integer(length(attrb$mlt)), file=file)
for(j in 1:length(attrb$mlt)) {
rmfi_write_variables(attrb$hgunam[j], attrb$mlt[j], attrb$zon[j], rmfi_ifelse0(attrb$zon[j] != 'ALL', as.integer(attrb$iz[[j]]), ''), file=file)
}
}
}
#' Create an \code{RMODFLOW} lvda object
#'
#' \code{rmf_create_lvda} creates an \code{RMODFLOW} lvda object.
#'
#' @param parameters either a single \code{rmf_parameter} or a list of \code{rmf_parameters} specifying the angle between the grid axis and the principal direction of horizontal hydraulic conductivity. See details.
#' @param dis \code{RMODFLOW} dis object
#'
#' @return a \code{RMODFLOW} lvda object
#' @export
#' @details All parameters should have a layer attribute and their partyp attribute set to 'LVDA'.
#' The LVDA package can only be used in conjunction with the HUF package.
#' Note that the parameters are defined on the numerical grid, not the HUF grid
#' @seealso \code{\link{rmf_read_lvda}}, \code{\link{rmf_write_lvda}} and \url{http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?lvda.htm}
rmf_create_lvda <- function(parameters, dis) {
if(length(parameters) == 1 && inherits(parameters[[1]], 'list') && !(inherits(parameters[[1]], 'rmf_parameter'))) parameters <- parameters[[1]]
lvda <- list()
# data set 0
# to provide comments, use ?comment on the resulting lvda object
# data set 1
lvda$nplvda <- length(parameters)
# data set 2 & 3
# error check
if(any(vapply(parameters, function(i) is.null(attr(i, 'partyp')) || is.null(attr(i, 'layer')) || is.null(attr(i, 'parnam')) || is.null(attr(i, 'parval')), TRUE))) {
stop('Please make sure all parameters have a parnam, parval, partyp and layer attribute', call. = FALSE)
}
lvda$parameters <- list()
lvda$parameter_values <- NULL
lvda$lvda <- rmf_create_array(0, dim = c(dis$nrow, dis$ncol, dis$nlay))
for(i in 1:lvda$nplvda) {
attrb <- attributes(parameters[[i]])
parnam <- attrb$parnam
lvda$parameter_values[parnam] <- attrb$parval
lvda$parameters[[parnam]] <- parameters[[i]]
lvda$lvda[,,unique(attrb$layer)] <- lvda$lvda[,,unique(attrb$layer)] + c(lvda$parameters[[parnam]])
}
class(lvda) <- c('lvda', 'rmf_package')
return(lvda)
}
#' Read a MODFLOW model-layer variable-direction horizontal anisotropy capability file
#'
#' \code{rmf_read_lvda} reads in a MODFLOW Model-Layer Variable-Direction Horizontal Anisotropy Capability file and returns it as an \code{RMODFLOW} lvda object.
#'
#' @param filename filename; typically *.lvda
#' @param dis \code{RMODFLOW} dis object
#' @param ... ignored
#' @details Note that the parameters are defined on the numerical grid, not the HUF grid
#' @return object of class lvda
#' @export
#' @seealso \code{\link{rmf_create_lvda}}, \code{\link{rmf_write_lvda}} and \url{http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?lvda.htm}
rmf_read_lvda <- function(file,
dis = {cat('Please select corresponding dis file ...\n'); rmf_read_dis(file.choose())},
...) {
lvda_lines <- readr::read_lines(file, lazy = FALSE)
lvda <- list()
# data set 0
data_set_0 <- rmfi_parse_comments(lvda_lines)
comments <- data_set_0$comments
lvda_lines <- data_set_0$remaining_lines
rm(data_set_0)
# data set 1
data_set_1 <- rmfi_parse_variables(lvda_lines)
lvda$nplvda <- as.numeric(data_set_1$variables[1])
lvda_lines <- data_set_1$remaining_lines
rm(data_set_1)
lvda$lvda <- rmf_create_array(0, dim = c(dis$nrow, dis$ncol, dis$nlay))
# data set 2-3
for(i in 1:lvda$nplvda) {
data_set_2 <- rmfi_parse_variables(lvda_lines)
parnam <- data_set_2$variables[1]
partyp <- data_set_2$variables[2]
parval <- as.numeric(data_set_2$variables[3])
nclu <- as.numeric(data_set_2$variables[4])
lvda_lines <- data_set_2$remaining_lines
rm(data_set_2)
ds3 <- list(layer = NULL, mltarr = NULL, zonarr = NULL, iz = list())
for(j in 1:nclu) {
data_set_3 <- rmfi_parse_variables(lvda_lines)
ds3$layer[j] <- data_set_3$variables[1]
ds3$mltarr[j] <- data_set_3$variables[2]
ds3$zonarr[j] <- data_set_3$variables[3]
# zero or character entry terminates IZ
if(toupper(ds3$zonarr[j]) == 'ALL') {
ds3$iz[[j]] <- NULL
} else {
iz <- suppressWarnings(as.numeric(data_set_3$variables[4:length(data_set_3$variables)]))
ds3$iz[[j]] <- iz[1:min(length(iz), which(is.na(iz))[1] - 1, which(iz == 0)[1] - 1, na.rm = TRUE)]
}
lvda_lines <- data_set_3$remaining_lines
if(toupper(data_set_3$variables[2]) != 'NONE') {
if(is.null(mlt)) stop('Please provide a mlt object', call. = FALSE)
}
if(toupper(data_set_3$variables[3]) != 'ALL') {
if(is.null(zon)) stop('Please provide a zon object', call. = FALSE)
}
rm(data_set_3)
}
lvda$parameter_values[parnam] <- parval
lvda$parameters[[parnam]] <- rmf_create_parameter(dis = dis, parnam = parnam, partyp = partyp, parval = parval, layer = ds3$layer, mltnam = ds3$mltarr, zonnam = ds3$zonarr, iz = ds3$iz, mlt = mlt, zon = zon)
# if(is.null(lvda[[tolower(partyp)]])) lvda[[tolower(partyp)]] <- rmf_create_array(NA, dim = c(dis$nrow, dis$ncol, dis$nlay))
# lvda[[tolower(partyp)]][,,unique(ds3$layer)] <- c(lvda$parameters[[parnam]])
#
# if(!(toupper(partyp) %in% types)) types <- append(types, toupper(partyp))
lvda$lvda[,,unique(ds3$layer)] <- lvda$lvda[,,unique(ds3$layer)] + c(lvda$parameters[[parnam]])
}
comment(lvda) <- comments
class(lvda) <- c('lvda','rmf_package')
return(lvda)
}
#' Write a MODFLOW model-layer variable-direction horizontal anisotropy capability file
#'
#' @param lvda an \code{RMODFLOW} lvda object
#' @param file filename to write to; typically '*.lvda'
#' @param ... ignored
#'
#' @return \code{NULL}
#' @export
#'
#' @seealso \code{\link{rmf_create_lvda}}, \code{\link{rmf_read_lvda}} and \url{http://water.usgs.gov/nrp/gwsoftware/modflow2000/MFDOC/index.html?lvda.htm}
rmf_write_lvda <- function(lvda,
file = {cat('Please select lvda file to overwrite or provide new filename ...\n'); file.choose()},
...) {
# data set 0
v <- packageDescription("RMODFLOW")$Version
cat(paste('# MODFLOW Model-Layer Variable-Direction Horizontal Anisotropy Capability Package created by RMODFLOW, version',v,'\n'), file = file)
cat(paste('#', comment(lvda)), sep='\n', file=file, append=TRUE)
# data set 1
rmfi_write_variables(lvda$nplvda, file = file, integer = TRUE)
# data set 3-4
for(i in 1:lvda$nplvda) {
attrb <- attributes(lvda$parameters[[i]])
rmfi_write_variables(attrb$parnam, attrb$partyp, attrb$parval, as.integer(length(attrb$mlt)), file=file)
for(j in 1:length(attrb$mlt)) {
rmfi_write_variables(as.integer(attrb$layer[j]), attrb$mlt[j], attrb$zon[j], rmfi_ifelse0(attrb$zon[j] != 'ALL', as.integer(attrb$iz[[j]]), ''), file=file)
}
}
}
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.