pyvista-modflow

Read MODFLOW-USG unstructured grids into PyVista as exact, flat-layered volumes.

Built and maintained by CoDimensional PBC.

_{FloPy's Freyberg MODFLOW-USG model, layers pulled apart, 60× vertical exaggeration (scripts/render_readme.py).}

Install

pip install pyvista-modflow

It depends only on NumPy and PyVista (0.49 or later, with VTK 9.4 or later), and it registers itself as the PyVista reader for .dis and .disu files.

Quickstart

import pyvista as pv

grid = pv.read('model.disu')  # finds model.gsf or model.shp beside it
grid = pv.read('model.disu', footprints='cells.shp')  # or name the outlines
grid.plot(show_edges=True)  # colored by layer

pyvista_modflow.read(...) takes the same arguments. The grid carries two cell arrays: layer, and node, the 1-based MODFLOW node number. Use node to attach heads or any other per-node result:

import numpy as np

heads = np.loadtxt('heads.txt')  # one value per MODFLOW node, from any post-processor
grid.cell_data['head'] = heads[grid.cell_data['node'] - 1]

Command line

Installing the package also lets the pyvista command read DISU files. The footprint file must sit beside the DISU with the same name (model.gsf or model.shp).

pyvista convert model.disu model.vtu  # one model
pyvista convert "*.disu" .vtu         # every model in the folder
pyvista plot model.disu

Recipes

Aquifers are thin, so exaggerate z before looking at them:

grid.scale([1, 1, 20]).plot(scalars='head', show_edges=True)

One layer, a cross-section, or a cutaway of whole cells:

layer2 = grid.threshold([2, 2], scalars='layer')
section = grid.slice(normal='y')
cutaway = grid.clip(normal='x', crinkle=True)

Save it for ParaView or any other VTK reader:

grid.save('model.vtu')

Background

A MODFLOW-USG DISU file holds each cell's top, bottom and area, but no x or y coordinates. The cell outlines come from a sidecar file, and the usual conversions build the cells from that sidecar's vertices. Those vertices are shared between neighboring cells and carry a single z, so every cell is bent to follow an interpolated surface instead of its own flat top and bottom.

This reader takes each cell's outline from the sidecar and its elevations from the DISU, so every cell is a flat-topped prism. Neighboring columns break at different elevations, so the reader splits every side wall wherever a neighbor's top or bottom meets it, and inserts the extra corner a refined (quadtree or nested) neighbor puts part way along a coarse cell's edge. A wall shared by two cells is the same face in both, so the surface PyVista extracts has no interior walls. Tops and bottoms are not split the same way: where a layer is refined differently from the one above it, the coarse cell's top and the fine cells' bottoms touch without sharing faces.

Inputs

Each outline's area is checked against the DISU AREA, which catches outlines matched to the wrong nodes; area_rtol=None skips it. A node whose TOP is not above its BOT has no volume and is skipped with a warning, and the node cell array maps every cell back to its MODFLOW node.

MODFLOW 6 (DISV, DISU, binary grid files) and structured DIS files are not read yet.

Development

just sync   # uv sync with the dev group
just test   # pytest
just lint   # pre-commit on every file
just figure # re-render the figure above

The test models in tests/data/ come from FloPy (public domain, CC0); see tests/data/README.md.