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(f"numpy version: {np.__version__}")
print(f"matplotlib version: {mpl.__version__}")
print(f"flopy version: {flopy.__version__}")
3.12.2 | packaged by conda-forge | (main, Feb 16 2024, 20:50:58) [GCC 12.3.0]
numpy version: 1.26.4
matplotlib version: 3.8.4
flopy version: 3.7.0.dev0
[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, crs=26916)
[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 = 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(
[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/tmp337w87bg/shapefile_export/model.shp')
[9]:
fname = f"{outdir}/wel.shp"
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 0x7f06041cb830>
[11]:
fname = f"{outdir}/hk.shp"
m.lpf.hk.export(f"{outdir}/hk.shp")
[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/tmp337w87bg/shapefile_export/hk.shp')
[13]:
m.riv.stress_period_data
[13]:
<flopy.utils.util_list.MfList at 0x7f0606bcf230>
[14]:
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
[15]:
m.riv.stress_period_data.export(f"{outdir}/riv_spd.shp", sparse=True)
[16]:
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
orPolygon
classes inflopy.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 0x7f06041398b0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1c710>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1c6b0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1c0e0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1c680>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1c5f0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1c170>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1c1a0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1c200>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1c3b0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1c4d0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1c7d0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1c830>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1c8c0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1c950>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1c9e0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1ca70>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1cb00>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1cb90>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1cc20>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1ccb0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1cd40>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1cdd0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1ce60>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1cef0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1cf80>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1d010>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1d0a0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1d130>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1d1c0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1d250>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1d2e0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1d370>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1d400>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1d490>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1d520>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1d5b0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1d640>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1d6d0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1d760>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1d7f0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1d880>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1d910>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1d9a0>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1da30>,
<flopy.utils.geometry.Polygon at 0x7f05fbc1dac0>]
[20]:
fname = f"{outdir}/bcs.shp"
recarray2shp(spd.to_records(), geoms=polygons, shpname=fname, crs=grid.epsg)
[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/tmp337w87bg/shapefile_export/bcs.shp')
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 | -1000.818465 | 277493.287550 | 5.094084e+06 |
1 | 1 | -1016.089049 | 278077.582741 | 5.093017e+06 |
2 | 2 | -1067.506321 | 275639.299205 | 5.090580e+06 |
3 | 3 | -929.310395 | 276093.262016 | 5.098619e+06 |
4 | 4 | -1079.372817 | 277544.916201 | 5.095658e+06 |
[23]:
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)
[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/tmp337w87bg/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
[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 = f"{outdir}/riv_reaches.shp"
recarray2shp(
rivdata.to_records(index=False),
geoms=lines,
shpname=lines_shapefile,
crs=grid.epsg,
)
[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/tmp337w87bg/shapefile_export/riv_reaches.shp')
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([ 1.2142898 , 2.1521107 , 1.31569332, 1.29411128, 1.68710723,
1.29470387, -0.47810135, 1.66789849, 2.16907214, 1.67434153,
1.33297642, 0.93729245, 0.03731037, 1.1365753 , 0.53530465,
1.32906057, 1.0978084 , 1.49603452, -0.20803492, -0.12902618,
1.18766483, 1.10388585, 1.0483132 , 0.6448886 , 1.33253496,
0.697038 , 2.59569236, 1.7574997 , 1.44184797, 2.1769406 ,
0.77563506, 2.64575455, 0.97887366, 0.39593215, 0.75126892,
1.18819337, -0.54491212, 0.63702057, -0.74875019, 1.37972087])
[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,
crs=grid.epsg,
)
[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/tmp337w87bg/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:
[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(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
[34]:
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.
[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/tmp337w87bg/shapefile_export/hk.shp')
[36]:
try:
# ignore PermissionError on Windows
temp_dir.cleanup()
except:
pass