Convert a MODFLOW model file to VTK file

Brief Description of the R code: 
This snippet reads a MODFLOW-2000/2005 DIS file (discretization file) and converts it to a VTK file which can be visualize in 3D in Paraview, VisIt and various VTK viewers.
R code: 
################################################################################
#  This script Reads the discretization file of modflow 2000 / 2005 and outputs
#  a VTK file. Inputs are:
#  1. A MODFLOW 2000/2005 discretization file
#  2. Offset values for the MODFLOW grid coordinates (in the x and y direction)
#     Offsets are simply the coordinates of the lower left corner. Default values
#     are zero.
################################################################################
 
library(tcltk)
 
################################################################################
#  Functions
################################################################################
 
## TclTk input box
inputBox<-function(tmsg="Put your message here")
{
tt<-tktoplevel()
Name <- tclVar("0")
entry.Name <-tkentry(tt,width="20",textvariable=Name)
tkgrid(tklabel(tt,text=tmsg))
tkgrid(entry.Name)
OnOK <- function()
{
	NameVal <- tclvalue(Name)
	tkdestroy(tt)
	}
OK.but <-tkbutton(tt,text="   OK   ",command=OnOK)
tkbind(entry.Name, "<Return>",OnOK)
tkgrid(OK.but)
tkfocus(tt)
tkwait.window(tt)
res<-as.numeric(tclvalue(Name))
return(res)
}
 
 
# 2-dimensional  Array reading function mimicking MODFLOW source code 
u2drel<-function(zz,dummy.line,nrec,n.row){  
  n.col<-nrec/n.row # number of columns in the formatted Fortran record
  flag1<-as.numeric(substr(dummy.line,1,10))
  print("Dunnmy.line")
  print(dummy.line)
  ind1<-regexpr("(",dummy.line,fixed=TRUE)+1
  #print("ind1")
  print(ind1)
  ind2<-regexpr(")",dummy.line,fixed=TRUE)-1
#    flag2<-unlist(strsplit(dummy.line[2],split=")",fixed=TRUE))
#    flag3<-strsplit(unlist(flag2),split="(",fixed=TRUE)
#    flag3<-unlist(flag3)
    fmt<-substr(dummy.line,ind1,ind2)
    print(fmt)
  fmt.col<-rec.extr(fmt)[1]
  nl1<-n.col%/%fmt.col
  if (n.col%%fmt.col > 0) nl1<-nl1+1
  nl2<-nl1*n.row
  if(flag1==0) { 
    dis.array<-vector(length=nrec)
    dis.array[1:nrec]<-as.numeric(substr(dummy.line,11,20))
  } else {
    rec.extr(fmt)[2]->recLength
    #dis.array<-vector()
    dummy.array<-read.fwf(zz,widths=rep(recLength,fmt.col),n=nl2,fill=TRUE)
    dummy.mat<-as.matrix(dummy.array)
    dummy.vec<-c(t(dummy.mat))
    na.ind<-which(is.na(dummy.vec)==FALSE)
    dis.array<-dummy.vec[na.ind]      
  }
  dis.array<-as.numeric(dis.array)   
  return(dis.array)
}   
 
# One dimensional fix format array reading function 
u1drel<-function(zz,dummy.line,nrec){  
  dummy.line<-tolower(dummy.line)
  print(dummy.line)
  if (length(grep("free",dummy.line))>=1)  {
    dis.array<-scan(zz,n=nrec,what=numeric())
  } else { 
     flag1<-as.numeric(substr(dummy.line,1,10))
  print("Dummy.line")
  print(dummy.line)
  ind1<-regexpr("(",dummy.line,fixed=TRUE)+1
  ind2<-regexpr(")",dummy.line,fixed=TRUE)-1
    fmt<-substr(dummy.line,ind1,ind2)
    print(fmt)
    fmt.col<-rec.extr(fmt)[1]
    nl1<-nrec%/%fmt.col
    if (nrec%%fmt.col > 0) nl1<-nl1+1
    if(flag1==0) { 
      dis.array<-vector(length=nrec)
      dis.array[1:nrec]<-as.numeric(substr(dummy.line,11,20))
      } else {
      rec.extr(fmt)[2]->recLength
      dummy.array<-read.fwf(zz,widths=rep(recLength,fmt.col),n=nl1,fill=TRUE)
      dummy.mat<-as.matrix(dummy.array)
      dummy.vec<-c(t(dummy.mat))
      na.ind<-which(is.na(dummy.vec)==FALSE)
      dis.array<-dummy.vec[na.ind]      
    }
  }
  dis.array<-as.numeric(dis.array)   
  return(dis.array)
}
#        *******   Record length extraction from FORTRAN format       ************
rec.extr<-function(fmt) {
  fmt<-tolower(fmt)
  if ((fmt!="free")&&(fmt!="binary")) { 
    fmt<-gsub("e",".",fmt,fixed=TRUE)
    fmt<-gsub("f",".",fmt,fixed=TRUE)
    fmt<-gsub("g",".",fmt,fixed=TRUE)
    fmt<-gsub("i",".",fmt,fixed=TRUE)
    fmt<-unlist(strsplit(fmt,".",fixed=TRUE))
    recLength<-as.numeric(fmt)
  }
  if (fmt[1]=="free") recLength=0
  if (fmt[1]=="binary") recLength=0
  return(recLength)
}
 
 
 
#        *******Dis file reading to get model grid information************
read.grid<-function(dis.f)  {
zz <- file(dis.f, "rt")
amper<-"#"
head.count=0
while (amper=="#") {
     header.line<-readLines(zz,n=1)
     amper<-substr(header.line,1,1)
     head.count<-head.count+1
   }
head.count<-head.count-1
close(zz)
zz <- file(dis.f, "rt")
dim.line<-scan(zz,n=6,sep= "",skip=head.count)
n.col<-dim.line[3]
n.row<-dim.line[2]
n.lay<-dim.line[1]
conf.flag<-scan(zz,n=n.lay,sep= "")
dummy.line1<-readLines(zz,n=1)
print(dummy.line1)
deltax<-u1drel(zz,dummy.line1,n.col)
dummy.line2<-readLines(zz,n=1)
print(dummy.line2)
deltay<-u1drel(zz,unlist(dummy.line2),n.row)
 
grd.info<-list(4+n.lay-1)
 
for (i in 4:(4+n.lay)) {
   dummy.line<-readLines(zz,n=1)
   nrec<-n.row*n.col
   print(dummy.line)
   gen.array<-u2drel(zz,unlist(dummy.line),nrec,n.row)
   mat<-matrix(gen.array,nrow=n.row,byrow=TRUE)
   grd.info[[i]]<-t(mat)[,n.row:1]
}  
image(grd.info[[4]])
close(zz)
 
grd.info[[1]]<-dim.line
grd.info[[2]]<-deltax
grd.info[[3]]<-deltay
return(grd.info)
 
 
}
 
 
################################################################################
#   Main Program
################################################################################
 
 
## Reading Input file
dis.f<-tclvalue(tkgetOpenFile(title="Modflow DIS file",
                filetypes="{{MODFLOW Discretization files} 
                {.dis}} {{All files} *}"))
grd.info<-read.grid(dis.f)
 
# Defining X and Y dimensions
Xseq<-grd.info[[2]]
Yseq<-grd.info[[3]]
 
# Creating grid
Xseq1<-c(0,Xseq)
Yseq1<-c(0,Yseq)
 
Xseq1[1]<-0
for (i in 2:length(Xseq1)) {
   Xseq1[i]<-Xseq1[i-1]+Xseq[i-1]
   }
 
Yseq1[1]<-0  
for (i in 2:length(Yseq1)) {
    Yseq1[i]<-Yseq1[i-1]+Yseq[i-1]
   }
 
#Enter coordinates offset 
xoffset<-inputBox("Input Xoffset")
yoffset<-inputBox("Input Yoffset")
 
## Adding Offset
Xseq1<-Xseq1+xoffset
Yseq1<-Yseq1+yoffset
dim.line<-grd.info[[1]]   
n.col<-dim.line[3]
n.row<-dim.line[2]
n.lay<-dim.line[1]   
init.coord<-expand.grid(Xseq1,Yseq1)
init.coord<-cbind(rep(init.coord[,1],n.lay+1),rep(init.coord[,2],n.lay+1))
 
dummy<-grd.info[[4]]
dummy<-cbind(dummy[,1],dummy)
dummy<-rbind(dummy[1,],dummy)
 
for (i in 5:(n.lay+4)){
   dummy2<-grd.info[[i]]
   dummy2<-cbind(dummy2[,1],dummy2)
   dummy2<-rbind(dummy2[1,],dummy2)
   dummy<-append(dummy,dummy2)
}
final.coord<-cbind(init.coord,dummy)
tclvalue(tkgetSaveFile(title="Output vtk File Name",filetypes="{{VTK files} 
                {.vtk}} {{All files} *}"))->new.vtk
#close(zz)
zz1 <- file(new.vtk, "wt")
cat("# vtk DataFile Version 2.0",file=zz1)
cat("\n",file=zz1)
cat("Regolith data",file=zz1)
cat("\n",file=zz1)
cat("ASCII",file=zz1)
cat("\n",file=zz1)
cat("DATASET STRUCTURED_GRID",file=zz1)
cat("\n",file=zz1)
cat("DIMENSIONS", dim.line[3]+1,dim.line[2]+1,dim.line[1]+1,sep=" ",file=zz1)
cat("\n",file=zz1)
t_pts<-(n.lay+1)*(n.row+1)*(n.col+1)
cat("POINTS",t_pts,"float", sep=" ",file=zz1)
cat("\n",file=zz1)
write.table(final.coord,file=zz1,sep=" ",row.names=FALSE,col.names=FALSE, append=TRUE)
cat("CELL_DATA",n.lay*n.row*n.col,file=zz1,sep=" ")
cat("\n",file=zz1)
cat("SCALARS layerNum int",file=zz1)
cat("\n",file=zz1)
cat("LOOKUP_TABLE default",file=zz1)
cat("\n",file=zz1)
write(rep(1:(n.lay),each=n.row*n.col), file=zz1,ncolumns=1,append=TRUE)
close(zz1)

Created by Pretty R at inside-R.org