R and MODFLOW
[This article was first published on Bart Rogiers - Sreigor, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Here are some functions for reading and writing MODFLOW files from R. I hope to update this in the future!
################################################################################
### read.modflow.pval ##########################################################
################################################################################
read.modflow.pval <- function(filename, read.all=F)
{
pval.lines <- scan(filename, what=character(), sep='\n')
pval <- NULL
# Data set 0
pval.lines <- remove.comments.from.lines(pval.lines)
# Data set 1
ifelse(read.all, pval$NP <- length(pval.lines)-1, pval$NP <- as.numeric(pval.lines[1]))
pval.lines <- pval.lines[-1]
# Data set 2
for(i in 1:pval$NP)
{
#print(strsplit(pval.lines[1],' '))
pval$PARNAM[i] <- as.character(strsplit(pval.lines[1],' ')[[1]][1])
pval$Parval[i] <- as.numeric(remove.empty.strings(strsplit(pval.lines[1],' ')[[1]])[2])
pval.lines <- pval.lines[-1]
}
return(pval)
}
################################################################################
### read.modflow.dis ###########################################################
################################################################################
read.modflow.dis <- function(filename)
{
dis.lines <- scan(filename, what=character(), sep='\n')
# Data set 0
dis.lines <- remove.comments.from.lines(dis.lines)
# Data set 1
dis.dataset1 <- strsplit(dis.lines[1],' ')
dis <- NULL
dis$NLAY <- as.numeric(dis.dataset1[[1]][1])
dis$NROW <- as.numeric(dis.dataset1[[1]][2])
dis$NCOL <- as.numeric(dis.dataset1[[1]][3])
dis$NPER <- as.numeric(dis.dataset1[[1]][4])
dis$ITMUNI <- as.numeric(dis.dataset1[[1]][5])
dis$LENUNI <- as.numeric(dis.dataset1[[1]][6])
dis.lines <- dis.lines[-1]
# Data set 2
dis$LAYCBD <- as.numeric(strsplit(dis.lines[1],' ')[[1]])
dis.lines <- dis.lines[-1]
# Data set 3
dis.lines <- dis.lines[-1]
dis$DELR <- as.numeric(strsplit(dis.lines[1],' ')[[1]])
dis.lines <- dis.lines[-1]
# Data set 4
dis.lines <- dis.lines[-1]
dis$DELC <- as.numeric(strsplit(dis.lines[1],' ')[[1]])
dis.lines <- dis.lines[-1]
# Data set 5
dis.lines <- dis.lines[-1]
dis$TOP <- matrix(nrow=dis$NROW, ncol=dis$NCOL)
for(i in 1:dis$NROW)
{
dis$TOP[i,] <- as.numeric(strsplit(dis.lines[1],' ')[[1]])
dis.lines <- dis.lines[-1]
}
# Data set 6
dis$BOTM <- list(rep(matrix(nrow=180, ncol=250),dis$NLAY+sum(dis$LAYCBD)))
for(i in 1:(dis$NLAY+sum(dis$LAYCBD)))
{
dis.lines <- dis.lines[-1]
botmatrix <- matrix(nrow=dis$NROW, ncol=dis$NCOL)
for(j in 1:dis$NROW)
{
botmatrix[j,] <- as.numeric(strsplit(dis.lines[1],' ')[[1]])
dis.lines <- dis.lines[-1]
}
dis$BOTM[[i]] <- botmatrix
}
# Data set 7
dis.dataset7 <- strsplit(dis.lines[1],' ')[[1]]
dis$PERLEN <- as.numeric(dis.dataset7[1])
dis$NSTP <- as.numeric(dis.dataset7[2])
dis$TSMULT <- as.numeric(dis.dataset7[3])
dis$SSTR <- as.character(dis.dataset7[4])
return(dis)
}
################################################################################
### read.modflow.mlt ###########################################################
################################################################################
read.modflow.mlt <- function(filename, dis)
{
mlt <- NULL
mlt.lines <- scan(filename, what=character(), sep='\n')
# Data set 0
mlt.lines <- remove.comments.from.lines(mlt.lines)
# Data set 1
mlt$NML <- as.numeric(mlt.lines[1])
mlt.lines <- mlt.lines[-1]
# Data set 2 + 3
mlt$RMLT <- list()
for(i in 1:mlt$NML)
{
mlt$MLTNAM[i] <- as.character(strsplit(mlt.lines[1],' ')[1])
mlt.lines <- mlt.lines[-1]
if(strsplit(mlt.lines[1],' ')[[1]][1]=='CONSTANT') {mlt$RMLT[[i]] <- as.numeric(strsplit(mlt.lines[1],' ')[[1]][2]);mlt.lines <- mlt.lines[-1]}
else if(strsplit(mlt.lines[1],' ')[[1]][1]=='INTERNAL')
{
mlt.lines <- mlt.lines[-1]
mlt$RMLT[[i]] <- matrix(nrow=dis$NROW, ncol=dis$NCOL)
for(j in 1:dis$NROW)
{
mlt$RMLT[[i]][j,] <- as.numeric(strsplit(mlt.lines[1],' ')[[1]])
mlt.lines <- mlt.lines[-1]
}
}
}
return(mlt)
}
################################################################################
### write.modflow.mlt ##########################################################
################################################################################
write.modflow.mlt <- function(mlt, filename, info='No further information provided')
{
# Data set 0
cat('# MODFLOW Multiplier File created in R\n', file=filename)
cat(paste('#', info, '\n'), file=filename, append=TRUE)
# Data set 1
cat(paste(mlt$NML, '\n', sep=''), file=filename, append=TRUE)
# Data set 2 + 3
for(i in 1:mlt$NML)
{
cat(paste(mlt$MLTNAM[i], '\n', sep=''), file=filename, append=TRUE)
if(length(mlt$RMLT[[i]])==1)
{
cat(paste('CONSTANT ', mlt$RMLT[[i]], '\n', sep=''), file=filename, append=TRUE)
}
else
{
cat(paste('INTERNAL 1.0 (free) 0', '\n', sep=''), file=filename, append=TRUE)
write.table(mlt$RMLT[[i]], file=filename, append=TRUE, sep=' ', col.names=FALSE, row.names=FALSE)
}
}
}
################################################################################
### read.gms.2dgrid ############################################################
################################################################################
read.gms.grid2d <- function(filename)
{
grid2d <- NULL
grid2d.lines <- scan(filename, what=character(), sep='\n')
#2dgrid.lines <- remove.comments.from.lines(mlt.lines)
grid2d.lines <- grid2d.lines[-1]
grid2d$objtype <- as.character(strsplit(grid2d.lines[1],'\"')[[1]][2])
grid2d.lines <- grid2d.lines[-1]
grid2d.lines <- grid2d.lines[-1]
grid2d$nd <- as.numeric(strsplit(grid2d.lines[1],' ')[[1]][3])
grid2d.lines <- grid2d.lines[-1]
grid2d$nc <- as.numeric(strsplit(grid2d.lines[1],' ')[[1]][3])
grid2d.lines <- grid2d.lines[-1]
grid2d$n <- as.character(strsplit(grid2d.lines[1],'\"')[[1]][2])
grid2d.lines <- grid2d.lines[-1]
grid2d.lines <- grid2d.lines[-1]
grid2d$nddata <- as.numeric(grid2d.lines[1:grid2d$nd])
grid2d$ncdata <- as.numeric(grid2d.lines[(1+grid2d$nd):(grid2d$nd + grid2d$nc)])
return(grid2d)
}
################################################################################
### Utilities ##################################################################
################################################################################
# Remove beginning comment lines.
remove.comments.from.lines <- function(lines)
{
i <- 0
while(i==0) ifelse(substr(lines[1], 1,1)=='#', lines <- lines[-1], i<-1)
#if(i==1) cat('Comments removed\n')
return(lines)
}
# Return column information from dis file
node.info <- function(rownumber, colnumber, dis)
{
cat('Column width = ',dis$DELR[colnumber], '\n')
cat('Row width = ', dis$DELC[rownumber], '\n')
cat('Vertical boundaries:\n')
# layers: top bottom thickness
cat('\t\t Top \t\t Bottom \t Thickness\n', sep='')
cat('Layer 1:\t', dis$TOP[rownumber, colnumber], '\t', dis$BOTM[[1]][rownumber,colnumber],'\t', dis$TOP[rownumber, colnumber]-dis$BOTM[[1]][rownumber,colnumber],'\n', sep='')
for(i in 2:dis$NLAY)
{
cat('Layer ',i,':\t', dis$BOTM[[i-1]][rownumber, colnumber], '\t', dis$BOTM[[i]][rownumber,colnumber],'\t', dis$BOTM[[i-1]][rownumber, colnumber]-dis$BOTM[[i]][rownumber,colnumber],'\n', sep='')
}
}
# Remove empty strings from string array
remove.empty.strings <- function(stringArray)
{
newStringArray <- NULL
for(i in 1:length(stringArray))
{
#print(stringArray[i])
if(stringArray[i] != '') {newStringArray <- c(newStringArray, stringArray[i])}
}
return(newStringArray)
}
To leave a comment for the author, please follow the link and comment on their blog: Bart Rogiers - Sreigor.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.