Example of interaction with Flopy and Modflow
Flopy provides a python interface for the USGS Modflow model. In this example, we illustrate the potential of Qgridder as a Qgis based preprocessor for Flopy. We consider the propagation of the tidal signal in a theoretical island. The simulation is transient, sea level is imposed at the coast of the island.
You can download the dataset with python scripts and follow the instructions below.
- Qgis 1.8 with Numpy and Matplotlib python packages
- Qgridder
- Flopy
- The example dataset.
There is no need to install an independent python distribution, Qgis contains its own.
In Qgis, load "island.shp" and start the Qgridder plugin. From the Qgridder "Build New Grid" tab, select the layer "island" and click "Update Extents from layer". Set grid resolution to 50 m. Set the output path and click "Build grid". The grid is ready.
Go to the Qgridder "Pre-processing" tab. Check that "Model type" is set to "Modflow" at the top of the window. Click "Proceed to numbering". This will add fields ROW and COL to the grid attribute table.
- IBOUND of active cells must be set to 1. In Qgis, open the menu Vector / Research Tools / Select by Location. Select features in the grid layer that intersect features in the island layer. Open the grid layer attribute table, start editing, and click the calculator icon. Check "Create new field", set "Name" to IBOUND, "Type" to Whole number (integer), and "Width" to 10. Type 1 in the "expression" text box below. Click OK in the calculator, remain in the attribute table.
- IBOUND of inactive cells must be set to 0. From the grid attribute table, click "invert selection". Open the calculator, check "Update existing field" and set IBOUND of selected features to 0.
Subsequent pre-rocessing steps must be completed from the Qgis python console. Available from the menu Plugin / Python Console. The whole python script is available here, but the principal Qgridder python commands are detailed hereafter.
Before running the commands detailed hereafter, you must set a variable pointing to the grid layer. First select the grid layer in the layer list panel, and type :
gridLayer = qgis.utils.iface.activeLayer()
nrow, ncol = qgridder_utils.get_rgrid_nrow_ncol(gridLayer)
delr, delc = qgridder_utils.get_rgrid_delr_delc(gridLayer)
Following a flopy standard, delr and delc will be scalars for regular grid, list otherwise.
This can be used to set any modflow distributed variable (IBOUND, TRANS, ...)
ibound = qgridder_utils.get_param(gridLayer, output_type = 'array', fieldName = 'IBOUND')
This can be used to set boundary conditions. Note that you must first select the cells of interest in the grid layer.
bcells = qgridder_utils.get_param(gridLayer, output_type = 'list', layer = 1)
This can be used for post-processing purposes.
cx = qgridder_utils.get_param(gridLayer, output_type = 'array', fieldName = 'CX')
cy = qgridder_utils.get_param(gridLayer, output_type = 'array', fieldName = 'CY')
The result may be visualized with an animated gif. The python script is avaiable here.
