Geospatial Referentiation of a MODFLOW Model with Flopy - Tutorial

Groundwater models are geospatial referenced (unless you are in a laboratory) since we represent the actual and future conditions of a certain porous / fractured media, however the actual model matrix resolution is spatially independent since it deals with a hydrogeological conceptualized array of columns, rows and layers. The nexus in between the matrix and the groundwater flow system has been a topic in the model development, even later versions of Model Muse ask for the model system of reference (as EPSG or Proj4 code), however the user has to keep in mind the water heads and where those water heads are located.  The USGS has developed a Python package called Flopy for the model construction, simulation and output representation; this package has interesting features for the interaction with the input and output data and for the georeferentiation of model data. This tutorial shows the complete procedure to geospatial reference a MODFLOW NWT model with some lines of Python code in a Jupyter Notebook and the representation of geospatial referenced model discretization data in QGIS 3.

Tutorial

Input Data

You can download the input data from this link.

Code

This is the Python code:

# coding: utf-8

# In[1]:

get_ipython().magic('matplotlib inline')
import time
import os
import sys
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import flopy

print(sys.version)
print('numpy version: {}'.format(np.__version__))
print('matplotlib version: {}'.format(mpl.__version__))
print('pandas version: {}'.format(pd.__version__))
print('flopy version: {}'.format(flopy.__version__))


# In[2]:

modflownwt = '../bin/MODFLOW-NWT_64.exe'


# In[3]:

model_ws = '../model'
ml = flopy.modflow.Modflow.load("Model1.nam", model_ws=model_ws, verbose=False,
                               check=False, exe_name=modflownwt)


# In[4]:

ml.dis.sr


# In[5]:

import re
dislines = open('../model/Model1.dis').readlines()
#extracting features from the upper left corner
upperleft = re.sub('\(|\)',',',dislines[1]).split(',')[1:3]
modelxul = float(upperleft[0])
modelyul = float(upperleft[1])

#extracting gridrotation
gridangle = float(re.sub('\\n',' ',dislines[5]).split(' ')[6])

print(modelxul,modelyul,gridangle)


# In[6]:

from flopy.utils.reference import SpatialReference
ml.sr = SpatialReference(delr=ml.dis.delr, delc=ml.dis.delc, 
                         xul=modelxul, yul=modelyul, rotation=gridangle, epsg=32717)


# In[7]:

ml.dis.sr


# In[8]:

ml.dis.plot(mflay=1)


# In[9]:

#Remember to have pyshp installed
#in mac: pip install pyshp
#in windows on anaconda prompt; python -m pip install pyshp
ml.dis.export(os.path.join('../shp',"Model1_dis.shp"))


# In[10]:

ml.upw.export(os.path.join('../shp',"Model1_upw.shp"))