How to import a Water Table from MODFLOW in QGIS with Python - Tutorial

Advances in groundwater modeling with MODFLOW allow us to have higher refinements on the representation of the water heads and water table as well as more capabilities in the representation of physical process related to groundwater flow. On a regional scale, we can deal with models of more than 500K elements and most times we need to represent this data on a GIS software for further study or the creation of figures for the end users, stakeholders and reports. By the use of Python scripts we can speed up the process of model output representation on a GIS software as QGIS.

Python scripts can be a little bit long and very declarative, but the process time is much smaller than the traditional clicking process on the GUI interface. The purpose is to store these scripts and use then every time one have to process the MODFLOW output data. 

 

Tutorial

 

Code

This is the complete code on the Python console in QGIS.

#####
#Part1
#####

#import the required libraries
import matplotlib.pyplot as plt
import os, re
import numpy as np
from osgeo import gdal, gdal_array, osr

#change directory to the model files path
os.chdir('C:\Users\Saul\Documents\Ih_MODFLOWResultsinQGIS\Model')

#open the *.DIS file and get the model boundaries and discretization
dis= open('Modelo1.dis').readlines()

lowerleft, upperright = re.sub('\(|\)' , ',' , dis[2]), re.sub('\(|\)' , ',' , dis[3])
lowerleft, upperright = lowerleft.split(','), upperright.split(',')

#border coordinates
xmin, ymin = float(lowerleft[1]), float(lowerleft[2])
xmax, ymax = float(upperright[1]), float(upperright[2])

#number of layers, rows and columns
nlay, nrows, ncols = int(dis[6].split()[0]), int(dis[6].split()[1]), int(dis[6].split()[2])

#grid resolution
xres, yres = (xmax-xmin)/ncols, (ymax-ymin)/nrows



#####
#Part2
#####

#open the *.FHD file to read the heads, beware that the script runs only on static models
hds = open('Modelo1.fhd').readlines()

#lines that separate heads from a layer
breakers=[x for x in range(len(hds)) if hds[x][:4]=='    ']

#counter plus a empty numpy array for model results
i=0
headmatrix = np.zeros((nlay,nrows,ncols))

#loop to read all the heads in one layer, reshape them to the grid dimensions and add to the numpy array
for salto in breakers:
    rowinicio = salto + 1
    rowfinal = salto + 2545 #here we insert the len of lines per layer
    matrixlist = []
    
    for linea in range(rowinicio,rowfinal,1):
        listaitem = [float(item) for item in hds[linea].split()]
        for item in listaitem: matrixlist.append(item)
 
    matrixarray = np.asarray(matrixlist)
    matrixarray = matrixarray.reshape(nrows,ncols)
    
    headmatrix[i,:,:]=matrixarray
    i+=1

#empty numpy array for the water table
wt = np.zeros((nrows,ncols))

#obtain the first positive or real head from the head array
for row in range(nrows):
    for col in range(ncols):
        a = headmatrix[:,row,col]
        if a[a>-1000] != []:                    #just in case there are some inactive zones
            wt[row,col] = a[a>-1000][0]
        else: wt[row,col] = None



#####
#Part3
#####

#definition of the raster name
outputraster = "waterlevel.tif"

#delete output raster if exits
try:
    os.remove(outputraster)
except OSError:
    pass

#creation of a raster and parameters for the array geotransformation
geotransform = [xmin,xres,0,ymax,0,-yres]
raster = gdal.GetDriverByName('GTiff').Create(outputraster,ncols,nrows,1,gdal.GDT_Float32)
raster.SetGeoTransform(geotransform)

#assign of a spatial reference by EPSG code
srs=osr.SpatialReference()
srs.ImportFromEPSG(32717)

#write results on the created raster
raster.SetProjection(srs.ExportToWkt())
raster.GetRasterBand(1).WriteArray(wt)

#add the raster to the canvas
iface.addRasterLayer(outputraster,"WaterTable","gdal")

 

Input data

You can download the required data here.