This page was generated from
shapefile_feature_examples.py.
It's also available as a notebook.
Working with shapefiles
This notebook shows some lower-level functionality in flopy for working with shapefiles including:
utils.geometryclasses for writing shapefiles of model input/output. For example, quickly writing a shapefile of model cells with errors identified by the checkerexamples of how the
PointandLineStringclasses can be used to quickly plot pathlines and endpoints from MODPATH (these are also used by thePathlineFileandEndpointFileclasses to write shapefiles of this output)
[1]:
import os
import shutil
import sys
import warnings
from pathlib import Path
from tempfile import TemporaryDirectory
import geopandas as gpd
import git
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pooch
import flopy
from flopy.utils import geometry
from flopy.utils.geometry import LineString, Point, Polygon
from flopy.utils.modpathfile import EndpointFile, PathlineFile
warnings.simplefilter("ignore", UserWarning)
print(sys.version)
print(f"numpy version: {np.__version__}")
print(f"matplotlib version: {mpl.__version__}")
print(f"flopy version: {flopy.__version__}")
3.10.18 | packaged by conda-forge | (main, Jun 4 2025, 14:45:41) [GCC 13.3.0]
numpy version: 2.2.6
matplotlib version: 3.10.6
flopy version: 3.11.0.dev0
write a numpy record array to a shapefile
in this case, we want to visualize output from the checker first make a toy model
[2]:
temp_dir = TemporaryDirectory()
workspace = temp_dir.name
# Check if we are in the repository and define the data path.
try:
root = Path(git.Repo(".", search_parent_directories=True).working_dir)
except:
root = None
data_path = root / "examples" / "data" if root else Path.cwd()
m = flopy.modflow.Modflow("toy_model", model_ws=workspace)
botm = np.zeros((2, 10, 10))
botm[0, :, :] = 1.5
botm[1, 5, 5] = 4 # negative layer thickness!
botm[1, 6, 6] = 4
dis = flopy.modflow.ModflowDis(
nrow=10, ncol=10, nlay=2, delr=100, delc=100, top=3, botm=botm, model=m
)
set coordinate information
[3]:
grid = m.modelgrid
grid.set_coord_info(xoff=600000, yoff=5170000, crs="EPSG:26715", angrot=45)
[4]:
chk = dis.check()
chk.summary_array
DIS PACKAGE DATA VALIDATION:
2 Errors:
2 instances of zero or negative thickness
Checks that passed:
thin cells (less than checker threshold of 1.0)
nan values in top array
nan values in bottom array
[4]:
rec.array([(np.str_('Error'), np.str_('DIS'), 1, 5, 5, -2.5, np.str_('zero or negative thickness')),
(np.str_('Error'), np.str_('DIS'), 1, 6, 6, -2.5, np.str_('zero or negative thickness'))],
dtype=[('type', 'O'), ('package', 'O'), ('k', '<i8'), ('i', '<i8'), ('j', '<i8'), ('value', '<f8'), ('desc', 'O')])
make geometry objects for the cells with errors
geometry objects allow the shapefile writer to be simpler and agnostic about the kind of geometry
[5]:
get_vertices = (
m.modelgrid.get_cell_vertices
) # function to get the referenced vertices for a model cell
geoms = [Polygon(get_vertices(i, j)) for i, j in chk.summary_array[["i", "j"]]]
[6]:
geoms[0].type
[6]:
'Polygon'
[7]:
geoms[0].exterior
[7]:
((np.float64(600000.0), np.float64(5170707.106781187)),
(np.float64(600070.7106781186), np.float64(5170777.817459305)),
(np.float64(600141.4213562373), np.float64(5170707.106781187)),
(np.float64(600070.7106781186), np.float64(5170636.396103068)))
[8]:
geoms[0].bounds
[8]:
(np.float64(600000.0),
np.float64(5170636.396103068),
np.float64(600141.4213562373),
np.float64(5170777.817459305))
[9]:
geoms[0].plot() # this feature requires descartes
[9]:
<Axes: >
create a GeoDataFrame and write it to shapefile
the projection (.prj) file can be written using an epsg code or crs
[10]:
from pathlib import Path
gdf = gpd.GeoDataFrame(data=chk.summary_array, geometry=geoms)
gdf.to_file(os.path.join(workspace, "test.shp"))
/home/runner/work/flopy/flopy/modflow6/.pixi/envs/rtd/lib/python3.10/site-packages/pyogrio/__init__.py:7: DeprecationWarning: The 'shapely.geos' module is deprecated, and will be removed in a future version. All attributes of 'shapely.geos' are available directly from the top-level 'shapely' namespace (since shapely 2.0.0).
import shapely.geos # noqa: F401
read it back in
We can use geopandas to read it back into memory
[11]:
gdf = gpd.read_file(os.path.join(workspace, "test.shp"))
gdf.head()
[11]:
| type | package | k | i | j | value | desc | geometry | |
|---|---|---|---|---|---|---|---|---|
| 0 | Error | DIS | 1 | 5 | 5 | -2.5 | zero or negative thickness | POLYGON ((600000 5170707.107, 600070.711 51707... |
| 1 | Error | DIS | 1 | 6 | 6 | -2.5 | zero or negative thickness | POLYGON ((600141.421 5170707.107, 600212.132 5... |
[12]:
gdf.geometry.plot()
# -
[12]:
<Axes: >
Other geometry types
Linestring
create geometry objects for pathlines from a MODPATH simulation
plot the paths using the built in plotting method
[13]:
fname = "EXAMPLE-3.pathline"
pthfile = pooch.retrieve(
url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/mp6/{fname}",
fname=fname,
path=data_path / "mp6",
known_hash=None,
)
pthfile = PathlineFile(pthfile)
pthdata = pthfile._data.view(np.recarray)
[14]:
length_mult = 1.0 # multiplier to convert coordinates from model to real world
rot = 0 # grid rotation
particles = np.unique(pthdata.particleid)
geoms = []
for pid in particles:
ra = pthdata[pthdata.particleid == pid]
x, y = geometry.rotate(
ra.x * length_mult, ra.y * length_mult, grid.xoffset, grid.yoffset, rot
)
z = ra.z
geoms.append(LineString(list(zip(x, y, z))))
[15]:
geoms[0]
[15]:
<flopy.utils.geometry.LineString at 0x7f7cf944c910>
[16]:
geoms[0].plot()
[16]:
<Axes: >
[17]:
fig, ax = plt.subplots()
for g in geoms:
g.plot(ax=ax)
ax.autoscale()
ax.set_aspect(1)
Points
[18]:
fname = "EXAMPLE-3.endpoint"
eptfile = pooch.retrieve(
url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/mp6/{fname}",
fname=fname,
path=data_path / "mp6",
known_hash=None,
)
eptfile = EndpointFile(eptfile)
eptdata = eptfile.get_alldata()
[19]:
x, y = geometry.rotate(
eptdata["x0"] * length_mult,
eptdata["y0"] * length_mult,
grid.xoffset,
grid.yoffset,
rot,
)
z = eptdata["z0"]
geoms = [Point(x[i], y[i], z[i]) for i in range(len(eptdata))]
[20]:
fig, ax = plt.subplots()
for g in geoms:
g.plot(ax=ax)
ax.autoscale()
ax.set_aspect(2e-6)
[21]:
try:
# ignore PermissionError on Windows
temp_dir.cleanup()
except:
pass