This page was generated from shapefile_export_example.py. It's also available as a notebook.

Shapefile export demo

The goal of this notebook is to demonstrate ways to export model information to shapefiles. This example will cover: * basic exporting of information for a model, individual package, or dataset * custom exporting of combined data from different packages * general exporting and importing of geographic data from other sources

import os

import sys
from tempfile import TemporaryDirectory

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

import flopy

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

# temporary directory
temp_dir = TemporaryDirectory()
outdir = os.path.join(temp_dir.name, "shapefile_export")

# load an existing model
model_ws = "../../examples/data/freyberg"
m = flopy.modflow.Modflow.load(
    "freyberg.nam",
    model_ws=model_ws,
    verbose=False,
    check=False,
    exe_name="mfnwt",
)

m.get_package_list()

['DIS', 'BAS6', 'LPF', 'WEL', 'RIV', 'RCH', 'OC', 'PCG']

set the model coordinate information

the coordinate information where the grid is located in a projected coordinate system (e.g. UTM)

grid = m.modelgrid
grid.set_coord_info(xoff=273170, yoff=5088657, crs=26916)

grid.extent

(np.float64(273170.0),
 np.float64(278170.0),
 np.float64(5088657.0),
 np.float64(5098657.0))

Declarative export using attached .export() methods

Export the whole model to a single shapefile

fname = f"{outdir}/model.shp"
m.export(fname)

/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for delc array, DIS package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for delr array, DIS package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for laycbd array, DIS package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for nstp array, DIS package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for perlen array, DIS package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for steady array, DIS package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for tsmult array, DIS package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for chani array, LPF package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for layavg array, LPF package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for laytyp array, LPF package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for layvka array, LPF package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for laywet array, LPF package
  warn(

ax = plt.subplot(1, 1, 1, aspect="equal")
extents = grid.extent
pc = flopy.plot.plot_shapefile(fname, ax=ax, edgecolor="k", facecolor="none")
ax.set_xlim(extents[0], extents[1])
ax.set_ylim(extents[2], extents[3])
ax.set_title(fname)

Text(0.5, 1.0, '/tmp/tmpqrba6fea/shapefile_export/model.shp')

fname = f"{outdir}/wel.shp"
m.wel.export(fname)

Export a package to a shapefile

Export a FloPy list or array object

m.lpf.hk

<flopy.utils.util_array.Util3d at 0x7f3af2b51d00>

fname = f"{outdir}/hk.shp"
m.lpf.hk.export(f"{outdir}/hk.shp")

ax = plt.subplot(1, 1, 1, aspect="equal")
extents = grid.extent
a = m.lpf.hk.array.ravel()
pc = flopy.plot.plot_shapefile(fname, ax=ax, a=a)
ax.set_xlim(extents[0], extents[1])
ax.set_ylim(extents[2], extents[3])
ax.set_title(fname)

Text(0.5, 1.0, '/tmp/tmpqrba6fea/shapefile_export/hk.shp')

m.riv.stress_period_data

<flopy.utils.util_list.MfList at 0x7f3af5f6c920>

m.riv.stress_period_data.export(f"{outdir}/riv_spd.shp")

MfList.export() exports the whole grid by default, regardless of the locations of the boundary cells

sparse=True only exports the boundary cells in the MfList

m.riv.stress_period_data.export(f"{outdir}/riv_spd.shp", sparse=True)

m.wel.stress_period_data.export(f"{outdir}/wel_spd.shp", sparse=True)

Ad-hoc exporting using recarray2shp

• The main idea is to create a recarray with all of the attribute information, and a list of geometry features (one feature per row in the recarray)

• each geometry feature is an instance of the Point, LineString or Polygon classes in flopy.utils.geometry. The shapefile format requires all the features to be of the same type.

• We will use pandas dataframes for these examples because they are easy to work with, and then convert them to recarrays prior to exporting.

from flopy.export.shapefile_utils import recarray2shp

combining data from different packages

write a shapefile of RIV and WEL package cells

wellspd = pd.DataFrame(m.wel.stress_period_data[0])
rivspd = pd.DataFrame(m.riv.stress_period_data[0])
spd = pd.concat([wellspd, rivspd])
spd.head()

	k	i	j	flux	iface	stage	cond	rbot
0	0	8	15	-0.00820	0.0	NaN	NaN	NaN
1	0	10	12	-0.00410	0.0	NaN	NaN	NaN
2	0	19	13	-0.00390	0.0	NaN	NaN	NaN
3	0	25	9	-0.00083	0.0	NaN	NaN	NaN
4	0	28	5	-0.00072	0.0	NaN	NaN	NaN

from flopy.utils.geometry import Polygon

vertices = []
for row, col in zip(spd.i, spd.j):
    vertices.append(grid.get_cell_vertices(row, col))
polygons = [Polygon(vrt) for vrt in vertices]
polygons

[<flopy.utils.geometry.Polygon at 0x7f3aea604ef0>,
 <flopy.utils.geometry.Polygon at 0x7f3aea604470>,
 <flopy.utils.geometry.Polygon at 0x7f3aea604e90>,
 <flopy.utils.geometry.Polygon at 0x7f3aea605520>,
 <flopy.utils.geometry.Polygon at 0x7f3aea605130>,
 <flopy.utils.geometry.Polygon at 0x7f3aea605a00>,
 <flopy.utils.geometry.Polygon at 0x7f3aea6059a0>,
 <flopy.utils.geometry.Polygon at 0x7f3aea6058b0>,
 <flopy.utils.geometry.Polygon at 0x7f3aea6055b0>,
 <flopy.utils.geometry.Polygon at 0x7f3aea605820>,
 <flopy.utils.geometry.Polygon at 0x7f3aea605880>,
 <flopy.utils.geometry.Polygon at 0x7f3aea604110>,
 <flopy.utils.geometry.Polygon at 0x7f3aea604830>,
 <flopy.utils.geometry.Polygon at 0x7f3aea605b50>,
 <flopy.utils.geometry.Polygon at 0x7f3aea605be0>,
 <flopy.utils.geometry.Polygon at 0x7f3aea605c70>,
 <flopy.utils.geometry.Polygon at 0x7f3aea605d00>,
 <flopy.utils.geometry.Polygon at 0x7f3aea605d90>,
 <flopy.utils.geometry.Polygon at 0x7f3aea605e20>,
 <flopy.utils.geometry.Polygon at 0x7f3aea605eb0>,
 <flopy.utils.geometry.Polygon at 0x7f3aea605f40>,
 <flopy.utils.geometry.Polygon at 0x7f3aea605fd0>,
 <flopy.utils.geometry.Polygon at 0x7f3aea606060>,
 <flopy.utils.geometry.Polygon at 0x7f3aea6060f0>,
 <flopy.utils.geometry.Polygon at 0x7f3aea606180>,
 <flopy.utils.geometry.Polygon at 0x7f3aea606210>,
 <flopy.utils.geometry.Polygon at 0x7f3aea6062a0>,
 <flopy.utils.geometry.Polygon at 0x7f3aea606330>,
 <flopy.utils.geometry.Polygon at 0x7f3aea6063c0>,
 <flopy.utils.geometry.Polygon at 0x7f3aea606450>,
 <flopy.utils.geometry.Polygon at 0x7f3aea6064e0>,
 <flopy.utils.geometry.Polygon at 0x7f3aea606570>,
 <flopy.utils.geometry.Polygon at 0x7f3aea606600>,
 <flopy.utils.geometry.Polygon at 0x7f3aea606690>,
 <flopy.utils.geometry.Polygon at 0x7f3aea606720>,
 <flopy.utils.geometry.Polygon at 0x7f3aea6067b0>,
 <flopy.utils.geometry.Polygon at 0x7f3aea606840>,
 <flopy.utils.geometry.Polygon at 0x7f3aea6068d0>,
 <flopy.utils.geometry.Polygon at 0x7f3aea606960>,
 <flopy.utils.geometry.Polygon at 0x7f3aea6069f0>,
 <flopy.utils.geometry.Polygon at 0x7f3aea606a80>,
 <flopy.utils.geometry.Polygon at 0x7f3aea606b10>,
 <flopy.utils.geometry.Polygon at 0x7f3aea606ba0>,
 <flopy.utils.geometry.Polygon at 0x7f3aea606c30>,
 <flopy.utils.geometry.Polygon at 0x7f3aea606cc0>,
 <flopy.utils.geometry.Polygon at 0x7f3aea606d50>]

fname = f"{outdir}/bcs.shp"
recarray2shp(spd.to_records(), geoms=polygons, shpname=fname, crs=grid.epsg)

ax = plt.subplot(1, 1, 1, aspect="equal")
extents = grid.extent
pc = flopy.plot.plot_shapefile(fname, ax=ax)
ax.set_xlim(extents[0], extents[1])
ax.set_ylim(extents[2], extents[3])
ax.set_title(fname)

Text(0.5, 1.0, '/tmp/tmpqrba6fea/shapefile_export/bcs.shp')

exporting other data

Suppose we have some well data with actual locations that we want to export to a shapefile

welldata = pd.DataFrame(
    {
        "wellID": np.arange(0, 10),
        "q": np.random.randn(10) * 100 - 1000,
        "x_utm": np.random.rand(10) * 5000 + grid.xoffset,
        "y_utm": grid.yoffset + np.random.rand(10) * 10000,
    }
)
welldata.head()

	wellID	q	x_utm	y_utm
0	0	-941.734584	276971.747573	5.095699e+06
1	1	-959.901006	276618.502072	5.096199e+06
2	2	-900.913792	275399.858091	5.093298e+06
3	3	-1029.627011	277402.031861	5.096737e+06
4	4	-958.939576	276254.012449	5.094008e+06

from flopy.utils.geometry import Point

geoms = [Point(x, y) for x, y in zip(welldata.x_utm, welldata.y_utm)]

fname = f"{outdir}/wel_data.shp"
recarray2shp(welldata.to_records(), geoms=geoms, shpname=fname, crs=grid.epsg)

ax = plt.subplot(1, 1, 1, aspect="equal")
extents = grid.extent
pc = flopy.plot.plot_shapefile(fname, ax=ax, radius=100)
ax.set_xlim(extents[0], extents[1])
ax.set_ylim(extents[2], extents[3])
ax.set_title(fname)

Text(0.5, 1.0, '/tmp/tmpqrba6fea/shapefile_export/wel_data.shp')

Adding attribute data to an existing shapefile

Suppose we have a GIS coverage representing the river in the riv package

from flopy.utils.geometry import LineString

### make up a linestring shapefile of the river reaches
i, j = m.riv.stress_period_data[0].i, m.riv.stress_period_data[0].j
x0 = grid.xyzcellcenters[0][i[0], j[0]]
x1 = grid.xyzcellcenters[0][i[-1], j[-1]]
y0 = grid.xyzcellcenters[1][i[0], j[0]]
y1 = grid.xyzcellcenters[1][i[-1], j[-1]]
x = np.linspace(x0, x1, m.nrow + 1)
y = np.linspace(y0, y1, m.nrow + 1)
l0 = zip(list(zip(x[:-1], y[:-1])), list(zip(x[1:], y[1:])))
lines = [LineString(l) for l in l0]

rivdata = pd.DataFrame(m.riv.stress_period_data[0])
rivdata["reach"] = np.arange(len(lines))
lines_shapefile = f"{outdir}/riv_reaches.shp"
recarray2shp(
    rivdata.to_records(index=False),
    geoms=lines,
    shpname=lines_shapefile,
    crs=grid.epsg,
)

ax = plt.subplot(1, 1, 1, aspect="equal")
extents = grid.extent
pc = flopy.plot.plot_shapefile(lines_shapefile, ax=ax, radius=25)
ax.set_xlim(extents[0], extents[1])
ax.set_ylim(extents[2], extents[3])
ax.set_title(lines_shapefile)

Text(0.5, 1.0, '/tmp/tmpqrba6fea/shapefile_export/riv_reaches.shp')

shp2recarray reads a shapefile into a numpy record array, which can easily be converted to a DataFrame

from flopy.export.shapefile_utils import shp2recarray

linesdata = shp2recarray(lines_shapefile)
linesdata = pd.DataFrame(linesdata)
linesdata.head()

	k	i	j	stage	cond	rbot	iface	reach	geometry
0	0	0	14	20.100000	0.05	20.00	0.0	0	<flopy.utils.geometry.LineString object at 0x7...
1	0	1	14	19.870001	0.05	19.75	0.0	1	<flopy.utils.geometry.LineString object at 0x7...
2	0	2	14	19.650000	0.05	19.50	0.0	2	<flopy.utils.geometry.LineString object at 0x7...
3	0	3	14	19.420000	0.05	19.25	0.0	3	<flopy.utils.geometry.LineString object at 0x7...
4	0	4	14	19.190001	0.05	19.00	0.0	4	<flopy.utils.geometry.LineString object at 0x7...

# make up some fluxes between the river and aquifer at each reach
q = np.random.randn(len(linesdata)) + 1
q

array([ 0.13783958, -0.01509178,  1.5108926 ,  0.97161615,  2.95655034,
        3.61131162,  3.13825095,  1.16022603,  2.29182296,  1.8203779 ,
        1.80463882,  0.49242326,  0.68449548,  0.066337  ,  0.9028356 ,
        1.83232057,  1.70588022,  0.98326336,  2.34521421,  1.88817166,
       -0.12852692,  2.96876627,  2.11589205,  1.04840417, -0.05795394,
        2.13247322,  2.181918  ,  1.99001218,  2.59941009, -0.08517741,
        0.25763747,  0.47859263,  0.65536752,  1.00432287,  0.69606944,
        1.12176559,  2.07923593,  3.45060047, -0.81698554,  1.23316581])

linesdata["qreach"] = q
linesdata["qstream"] = np.cumsum(q)

recarray2shp(
    linesdata.drop("geometry", axis=1).to_records(),
    geoms=linesdata.geometry.values,
    shpname=lines_shapefile,
    crs=grid.epsg,
)

ax = plt.subplot(1, 1, 1, aspect="equal")
extents = grid.extent
pc = flopy.plot.plot_shapefile(lines_shapefile, ax=ax, radius=25)
ax.set_xlim(extents[0], extents[1])
ax.set_ylim(extents[2], extents[3])
ax.set_title(lines_shapefile)

Text(0.5, 1.0, '/tmp/tmpqrba6fea/shapefile_export/riv_reaches.shp')

Overriding the model’s modelgrid with a user supplied modelgrid

In some cases it may be necessary to override the model’s modelgrid instance with a seperate modelgrid. An example of this is if the model discretization is in feet and the user would like it projected in meters. Exporting can be accomplished by supplying a modelgrid as a kwarg in any of the export() methods within flopy. Below is an example:

mg0 = m.modelgrid

# build a new modelgrid instance with discretization in meters
modelgrid = flopy.discretization.StructuredGrid(
    delc=mg0.delc * 0.3048,
    delr=mg0.delr * 0.3048,
    top=mg0.top,
    botm=mg0.botm,
    idomain=mg0.idomain,
    xoff=mg0.xoffset * 0.3048,
    yoff=mg0.yoffset * 0.3048,
)

# exporting an entire model
m.export(f"{outdir}/freyberg.shp", modelgrid=modelgrid)

No CRS information for writing a .prj file.
Supply an valid coordinate system reference to the attached modelgrid object or .export() method.
No CRS information for writing a .prj file.
Supply an valid coordinate system reference to the attached modelgrid object or .export() method.

/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for delc array, DIS package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for delr array, DIS package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for laycbd array, DIS package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for nstp array, DIS package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for perlen array, DIS package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for steady array, DIS package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for tsmult array, DIS package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for chani array, LPF package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for layavg array, LPF package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for laytyp array, LPF package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for layvka array, LPF package
  warn(
/home/runner/micromamba/envs/flopy/lib/python3.12/site-packages/flopy/export/shapefile_utils.py:319: UserWarning: Failed to get data for laywet array, LPF package
  warn(

And for a specific parameter the method is the same

fname = f"{outdir}/hk.shp"
m.lpf.hk.export(fname, modelgrid=modelgrid)

No CRS information for writing a .prj file.
Supply an valid coordinate system reference to the attached modelgrid object or .export() method.

ax = plt.subplot(1, 1, 1, aspect="equal")
extents = modelgrid.extent
a = m.lpf.hk.array.ravel()
pc = flopy.plot.plot_shapefile(fname, ax=ax, a=a)
ax.set_xlim(extents[0], extents[1])
ax.set_ylim(extents[2], extents[3])
ax.set_title(fname)

Text(0.5, 1.0, '/tmp/tmpqrba6fea/shapefile_export/hk.shp')

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