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)

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

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

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.4 | packaged by conda-forge | (main, Jun 17 2024, 10:23:07) [GCC 12.3.0]
numpy version: 2.0.1
matplotlib version: 3.9.1
flopy version: 3.8.0

write a numpy record array to a shapefile

in this case, we want to visualize output from the checker first make a toy model

temp_dir = TemporaryDirectory()
workspace = temp_dir.name

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

grid = m.modelgrid
grid.set_coord_info(xoff=600000, yoff=5170000, crs="EPSG:26715", angrot=45)

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

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

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"]]]

geoms[0].type

'Polygon'

geoms[0].exterior

((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)))

geoms[0].bounds

(np.float64(600000.0),
 np.float64(5170636.396103068),
 np.float64(600141.4213562373),
 np.float64(5170777.817459305))

geoms[0].plot()  # this feature requires descartes

<Axes: >

write the shapefile

• the projection (.prj) file can be written using an epsg code

• or copied from an existing .prj file

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"),
)

read it back in

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

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

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

ra.geometry[0].plot()

<Axes: >

Other geometry types

Linestring

• create geometry objects for pathlines from a MODPATH simulation

• plot the paths using the built in plotting method

pthfile = PathlineFile("../../examples/data/mp6/EXAMPLE-3.pathline")
pthdata = pthfile._data.view(np.recarray)

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

geoms[0]

<flopy.utils.geometry.LineString at 0x7f08c94f7ec0>

geoms[0].plot()

<Axes: >

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

Points

eptfile = EndpointFile("../../examples/data/mp6/EXAMPLE-3.endpoint")
eptdata = eptfile.get_alldata()

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))]

fig, ax = plt.subplots()
for g in geoms:
    g.plot(ax=ax)
ax.autoscale()
ax.set_aspect(2e-6)

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