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.

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.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)