Note

This page was generated from Notebooks/shapefile_export_example. Download as a Jupyter notebook (.ipynb) or Python script (.py) or launch interactively with Binder Binder badge.

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

[1]:

import os

[2]:

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("numpy version: {}".format(np.__version__))
print("matplotlib version: {}".format(mpl.__version__))
print("flopy version: {}".format(flopy.__version__))

3.8.18 (default, Aug 28 2023, 08:27:22)
[GCC 11.4.0]
numpy version: 1.24.4
matplotlib version: 3.7.3
flopy version: 3.4.3

[3]:

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

[4]:

m.get_package_list()

[4]:

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

[5]:

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

/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/utils/crs.py:119: PendingDeprecationWarning: the epsg argument will be deprecated and will be removed in version 3.4. Use crs instead.
  warnings.warn(

[6]:

grid.extent

[6]:

(273170.0, 278170.0, 5088657.0, 5098657.0)

Declarative export using attached .export() methods

Export the whole model to a single shapefile

[7]:

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

/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for delc array, DIS package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for delr array, DIS package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for laycbd array, DIS package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for nstp array, DIS package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for perlen array, DIS package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for steady array, DIS package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for tsmult array, DIS package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for chani array, LPF package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for layavg array, LPF package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for laytyp array, LPF package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for layvka array, LPF package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for laywet array, LPF package
  warn(

[8]:

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)

[8]:

Text(0.5, 1.0, '/tmp/tmptjfsm5qe/shapefile_export/model.shp')
../_images/Notebooks_shapefile_export_example_10_1.png 
[9]:

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

Export a package to a shapefile

Export a FloPy list or array object

[10]:

m.lpf.hk

[10]:

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

[11]:

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

[12]:

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)

[12]:

Text(0.5, 1.0, '/tmp/tmptjfsm5qe/shapefile_export/hk.shp')
../_images/Notebooks_shapefile_export_example_16_1.png 
[13]:

m.riv.stress_period_data

[13]:

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

[14]:

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

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

[15]:

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

wrote ../../../../home/runner/work/flopy/flopy/.docs/Notebooks

[16]:

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

wrote ../../../../home/runner/work/flopy/flopy/.docs/Notebooks

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.

[17]:

from flopy.export.shapefile_utils import recarray2shp

combining data from different packages

write a shapefile of RIV and WEL package cells

[18]:

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

[18]:
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 
[19]:

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

[19]:

[<flopy.utils.geometry.Polygon at 0x7ff0c54bde80>,
 <flopy.utils.geometry.Polygon at 0x7ff0c54ad730>,
 <flopy.utils.geometry.Polygon at 0x7ff0c54ade50>,
 <flopy.utils.geometry.Polygon at 0x7ff0c54ad160>,
 <flopy.utils.geometry.Polygon at 0x7ff0c54ad940>,
 <flopy.utils.geometry.Polygon at 0x7ff0c54ad3d0>,
 <flopy.utils.geometry.Polygon at 0x7ff0c54adbb0>,
 <flopy.utils.geometry.Polygon at 0x7ff0c54adee0>,
 <flopy.utils.geometry.Polygon at 0x7ff0c54adc40>,
 <flopy.utils.geometry.Polygon at 0x7ff0c54ad6a0>,
 <flopy.utils.geometry.Polygon at 0x7ff0c54ada00>,
 <flopy.utils.geometry.Polygon at 0x7ff0c54ad070>,
 <flopy.utils.geometry.Polygon at 0x7ff0c54ada90>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299070>,
 <flopy.utils.geometry.Polygon at 0x7ff0c32990d0>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299220>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299340>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299310>,
 <flopy.utils.geometry.Polygon at 0x7ff0c32994f0>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299580>,
 <flopy.utils.geometry.Polygon at 0x7ff0c32996a0>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299730>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299700>,
 <flopy.utils.geometry.Polygon at 0x7ff0c32998e0>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299820>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299a00>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299a60>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299b20>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299b80>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299c40>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299ca0>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299d60>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299df0>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299dc0>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299f10>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299fa0>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299160>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299040>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299670>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299790>,
 <flopy.utils.geometry.Polygon at 0x7ff0c3299940>,
 <flopy.utils.geometry.Polygon at 0x7ff0c32b72b0>,
 <flopy.utils.geometry.Polygon at 0x7ff0c32b7340>,
 <flopy.utils.geometry.Polygon at 0x7ff0c32b7430>,
 <flopy.utils.geometry.Polygon at 0x7ff0c32b74c0>,
 <flopy.utils.geometry.Polygon at 0x7ff0c32b7550>]

[20]:

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

wrote ../../../../home/runner/work/flopy/flopy/.docs/Notebooks

/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/utils/crs.py:119: PendingDeprecationWarning: the epsg argument will be deprecated and will be removed in version 3.4. Use crs instead.
  warnings.warn(

[21]:

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)

[21]:

Text(0.5, 1.0, '/tmp/tmptjfsm5qe/shapefile_export/bcs.shp')
../_images/Notebooks_shapefile_export_example_30_1.png

exporting other data

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

[22]:

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

[22]:
wellID q x_utm y_utm
0 0 -1092.359884 276660.473089 5.095095e+06
1 1 -974.519847 276306.978258 5.093062e+06
2 2 -889.823255 277969.021802 5.096678e+06
3 3 -1062.366654 277663.656335 5.090236e+06
4 4 -968.146430 277554.297277 5.090216e+06 
[23]:

from flopy.utils.geometry import Point

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

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

wrote ../../../../home/runner/work/flopy/flopy/.docs/Notebooks

/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/utils/crs.py:119: PendingDeprecationWarning: the epsg argument will be deprecated and will be removed in version 3.4. Use crs instead.
  warnings.warn(

[24]:

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)

[24]:

Text(0.5, 1.0, '/tmp/tmptjfsm5qe/shapefile_export/wel_data.shp')
../_images/Notebooks_shapefile_export_example_35_1.png

Adding attribute data to an existing shapefile

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

[25]:

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 = "{}/riv_reaches.shp".format(outdir)
recarray2shp(
    rivdata.to_records(index=False),
    geoms=lines,
    shpname=lines_shapefile,
    epsg=grid.epsg,
)

wrote ../../../../home/runner/work/flopy/flopy/.docs/Notebooks

/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/utils/crs.py:119: PendingDeprecationWarning: the epsg argument will be deprecated and will be removed in version 3.4. Use crs instead.
  warnings.warn(

[26]:

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)

[26]:

Text(0.5, 1.0, '/tmp/tmptjfsm5qe/shapefile_export/riv_reaches.shp')
../_images/Notebooks_shapefile_export_example_38_1.png

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

[27]:

from flopy.export.shapefile_utils import shp2recarray

[28]:

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

[28]:
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... 
[29]:

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

[29]:

array([-0.18403005,  2.35045997, -0.21018577,  2.56692767,  2.3528112 ,
        1.05254341,  2.98830938,  1.34309542,  0.99476571,  1.11428521,
        0.94711756,  1.24274964,  0.98880093, -1.35575662, -0.08662407,
        1.58616024, -0.03223463,  1.41355697, -0.52118678,  1.81407137,
        2.24146616,  1.20358632,  0.60193288,  1.76019454, -0.8627548 ,
        0.18114315, -0.49923485,  0.15026523, -0.19594543,  0.87682404,
       -0.33079539, -0.31925346,  1.1363268 , -0.23418252,  1.03831157,
        2.35007046,  0.90669645,  1.42527292, -0.14586778,  1.52432627])

[30]:

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

[31]:

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

wrote ../../../../home/runner/work/flopy/flopy/.docs/Notebooks

/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/utils/crs.py:119: PendingDeprecationWarning: the epsg argument will be deprecated and will be removed in version 3.4. Use crs instead.
  warnings.warn(

[32]:

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)

[32]:

Text(0.5, 1.0, '/tmp/tmptjfsm5qe/shapefile_export/riv_reaches.shp')
../_images/Notebooks_shapefile_export_example_47_1.png

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:

[33]:

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("{}/freyberg.shp".format(outdir), 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.

/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for delc array, DIS package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for delr array, DIS package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for laycbd array, DIS package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for nstp array, DIS package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for perlen array, DIS package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for steady array, DIS package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for tsmult array, DIS package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for chani array, LPF package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for layavg array, LPF package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for laytyp array, LPF package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for layvka array, LPF package
  warn(
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/flopy/export/shapefile_utils.py:307: UserWarning: Failed to get data for laywet array, LPF package
  warn(

And for a specific parameter the method is the same

[34]:

fname = "{}/hk.shp".format(outdir)
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.

[35]:

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)

[35]:

Text(0.5, 1.0, '/tmp/tmptjfsm5qe/shapefile_export/hk.shp')
../_images/Notebooks_shapefile_export_example_52_1.png 
[36]:

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