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:

recarray2shp convience function for writing a numpy record array to a shapefile
shp2recarray convience function for quickly reading a shapefile into a numpy recarray
utils.geometry classes for writing shapefiles of model input/output. For example, quickly writing a shapefile of model cells with errors identified by the checker
examples of how the Point and LineString classes can be used to quickly plot pathlines and endpoints from MODPATH (these are also used by the PathlineFile and EndpointFile classes to write shapefiles of this output)

[1]:

import os
import shutil
import sys
import warnings
from pathlib import Path
from tempfile import TemporaryDirectory

import git
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pooch

import flopy
from flopy.export.shapefile_utils import recarray2shp, shp2recarray
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.12.10 | packaged by conda-forge | (main, Apr 10 2025, 22:21:13) [GCC 13.3.0]
numpy version: 2.2.5
matplotlib version: 3.10.3
flopy version: 3.9.3

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: >

../_images/Notebooks_shapefile_feature_examples_12_1.png

write the shapefile

the projection (.prj) file can be written using an epsg code
or copied from an existing .prj file

[10]:

from pathlib import Path

recarray2shp(chk.summary_array, geoms, os.path.join(workspace, "test.shp"), crs=26715)
shape_path = os.path.join(workspace, "test.prj")

shutil.copy(shape_path, os.path.join(workspace, "26715.prj"))
recarray2shp(
    chk.summary_array,
    geoms,
    os.path.join(workspace, "test.shp"),
    prjfile=os.path.join(workspace, "26715.prj"),
)

/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/geopandas/_compat.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

read it back in

flopy geometry objects representing the shapes are stored in the ‘geometry’ field

[11]:

ra = shp2recarray(os.path.join(workspace, "test.shp"))
ra

[11]:

rec.array([('Error', 'DIS', 1, 5, 5, -2.5, 'zero or negative thickness', <flopy.utils.geometry.Polygon object at 0x7f99c221b920>),
           ('Error', 'DIS', 1, 6, 6, -2.5, 'zero or negative thickness', <flopy.utils.geometry.Polygon object at 0x7f99c221b260>)],
          dtype=[('type', 'O'), ('package', 'O'), ('k', '<i8'), ('i', '<i8'), ('j', '<i8'), ('value', '<f8'), ('desc', 'O'), ('geometry', 'O')])

[12]:

ra.geometry[0].plot()
# -

[12]:

<Axes: >

../_images/Notebooks_shapefile_feature_examples_17_1.png

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 0x7f99c20f0170>

[16]:

geoms[0].plot()

[16]:

<Axes: >

../_images/Notebooks_shapefile_feature_examples_22_1.png

[17]:

fig, ax = plt.subplots()
for g in geoms:
    g.plot(ax=ax)
ax.autoscale()
ax.set_aspect(1)

../_images/Notebooks_shapefile_feature_examples_23_0.png

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)

../_images/Notebooks_shapefile_feature_examples_27_0.png

[21]:

try:
    # ignore PermissionError on Windows
    temp_dir.cleanup()
except:
    pass