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. 





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


#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

#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


#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
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)

#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


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

#delete output raster if exits
except OSError:

#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)

#assign of a spatial reference by EPSG code

#write results on the created raster

#add the raster to the canvas


Input data

You can download the required data here.

Saul Montoya

Saul Montoya es Ingeniero Civil graduado de la Pontificia Universidad Católica del Perú en Lima con estudios de postgrado en Manejo e Ingeniería de Recursos Hídricos (Programa WAREM) de la Universidad de Stuttgart con mención en Ingeniería de Aguas Subterráneas y Hidroinformática.

Smiley face

Subscribe to our free e-newsletter for tutorials, articles, webminars, courses and more.