Skip to content

Commit

Permalink
map view example
Browse files Browse the repository at this point in the history
  • Loading branch information
@wpbonelli
wpbonelli committed Dec 10, 2024
1 parent de991a4 commit 4d57f72
Showing 1 changed file with 149 additions and 32 deletions.
181 changes: 149 additions & 32 deletions .docs/Notebooks/plot_map_view_example.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,12 +28,14 @@
# +
import os
import sys
from pathlib import Path
from pprint import pformat
from tempfile import TemporaryDirectory

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

import flopy
Expand All @@ -51,21 +53,41 @@
vmf6 = "mf6"
exe_name_mf6 = "mf6"
exe_mp = "mp6"
sim_name = "freyberg"

# Set the paths
loadpth = os.path.join("..", "..", "examples", "data", "freyberg")
tempdir = TemporaryDirectory()
modelpth = tempdir.name
modelpth = Path(tempdir.name)

# + [markdown] pycharm={"name": "#%% md\n"}
# ### Load and Run an Existing MODFLOW-2005 Model
# A model called the "Freyberg Model" is located in the loadpth folder. In the following code block, we load that model, then change into a new workspace (modelpth) where we recreate and run the model. For this to work properly, the MODFLOW-2005 executable (mf2005) must be in the path. We verify that it worked correctly by checking for the presence of freyberg.hds and freyberg.cbc.
# A model called the "Freyberg Model" is located in the modelpth folder. In the following code block, we load that model, then change into a new workspace (modelpth) where we recreate and run the model. For this to work properly, the MODFLOW-2005 executable (mf2005) must be in the path. We verify that it worked correctly by checking for the presence of freyberg.hds and freyberg.cbc.

file_names = {
"freyberg.bas": "63266024019fef07306b8b639c6c67d5e4b22f73e42dcaa9db18b5e0f692c097",
"freyberg.dis": "62d0163bf36c7ee9f7ee3683263e08a0abcdedf267beedce6dd181600380b0a2",
"freyberg.githds": "abe92497b55e6f6c73306e81399209e1cada34cf794a7867d776cfd18303673b",
"freyberg.gitlist": "aef02c664344a288264d5f21e08a748150e43bb721a16b0e3f423e6e3e293056",
"freyberg.lpf": "06500bff979424f58e5e4fbd07a7bdeb0c78f31bd08640196044b6ccefa7a1fe",
"freyberg.nam": "e66321007bb603ef55ed2ba41f4035ba6891da704a4cbd3967f0c66ef1532c8f",
"freyberg.oc": "532905839ccbfce01184980c230b6305812610b537520bf5a4abbcd3bd703ef4",
"freyberg.pcg": "0d1686fac4680219fffdb56909296c5031029974171e25d4304e70fa96ebfc38",
"freyberg.rch": "37a1e113a7ec16b61417d1fa9710dd111a595de738a367bd34fd4a359c480906",
"freyberg.riv": "7492a1d5eb23d6812ec7c8227d0ad4d1e1b35631a765c71182b71e3bd6a6d31d",
"freyberg.wel": "00aa55f59797c02f0be5318a523b36b168fc6651f238f34e8b0938c04292d3e7",
}
for fname, fhash in file_names.items():
pooch.retrieve(
url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/{sim_name}/{fname}",
fname=fname,
path=modelpth,
known_hash=fhash,
)

# +
ml = flopy.modflow.Modflow.load(
"freyberg.nam", model_ws=loadpth, exe_name=exe_name_2005, version=v2005
"freyberg.nam", model_ws=modelpth, exe_name=exe_name_2005, version=v2005
)
ml.change_model_ws(new_pth=modelpth)
ml.write_input()
success, buff = ml.run_model(silent=True, report=True)
assert success, pformat(buff)
Expand Down Expand Up @@ -432,6 +454,58 @@
#
# The shapefile must have intersecting geographic coordinates as the `PlotMapView` object in order for it to overlay correctly on the plot. The `plot_shapefile()` method and function do not use any of the projection information that may be stored with the shapefile. If you reset `xoff`, `yoff`, and `angrot` in the `ml.modelgrid.set_coord_info()` call below, you will see that the grid will no longer overlay correctly with the shapefile.

file_names = {
"bedrock_outcrop_hole.dbf": "c48510bc0b04405e4d3433e6cd892351c8342a7c46215f48332a7e6292249da6",
"bedrock_outcrop_hole.sbn": "48fd1496d84822c9637d7f3065edf4dfa2038406be8fa239cb451b1a3b28127c",
"bedrock_outcrop_hole.sbx": "9a36aee5f3a4bcff0a453ab743a7523ea19acb8841e8273bbda34f27d7237ea5",
"bedrock_outcrop_hole.shp": "25c241ac90dd47be28f761ba60ba94a511744f5219600e35a80a93f19ec99f97",
"bedrock_outcrop_hole.shx": "88b06395fa4c58ea04d300e10e6f6ea81e17fb0baa20d8ac78470d19101430be",
"bedrock_outcrop_hole_rotate14.dbf": "e05bbfc826fc069666a05e949acc833b54de51b14267c9c54b1c129b4a8ab82d",
"bedrock_outcrop_hole_rotate14.sbn": "136d8f86b8a13abc8f0386108228ca398037cf8c28ba6077086fd7e1fd54abf7",
"bedrock_outcrop_hole_rotate14.sbx": "1c2f2f2791db9c752fb1b355f13e46a8740ccd66654ae34d130172a3bdcda805",
"bedrock_outcrop_hole_rotate14.shp": "3e722d8fa9331ab498dbf9544085b30f60d2e38cc82a0955792d11a4e6a4419d",
"bedrock_outcrop_hole_rotate14.shp.xml": "ff6a3e80d10d9e68863ffe224e8130b862c13c2265d3a604342eb20a700d38fd",
"bedrock_outcrop_hole_rotate14.shx": "32a75461fab39b21769c474901254e7cbd24073c53d62b494fd70080cfcd3383",
"cross_section.cpg": "3ad3031f5503a4404af825262ee8232cc04d4ea6683d42c5dd0a2f2a27ac9824",
"cross_section.dbf": "3b050b1d296a7efe1b4f001c78030d5c81f79d3cd101d459e4426944fbd4e8e7",
"cross_section.sbn": "3b6a8f72f78f7b0d12e5823d6e8307040cfd5af88a8fb9427687d027aa805126",
"cross_section.sbx": "72e33139aaa99a8d12922af3774bd6b1a73613fc1bc852d1a1d1426ef48a832a",
"cross_section.shp": "0eb9e37dcbdbb5d932101c4c5bcb971271feb2c1d81d2a5f8dbc0fbf8d799ee5",
"cross_section.shp.xml": "ff99002ecd63a843fe628c107dfb02926b6838132c6f503db38b792644fb368e",
"cross_section.shx": "c6fa1307e1c32c535842796b24b2a0a07865065ace3324b0f6b1b71e9c1a8e1e",
"cross_section_rotate14.cpg": "3ad3031f5503a4404af825262ee8232cc04d4ea6683d42c5dd0a2f2a27ac9824",
"cross_section_rotate14.dbf": "72f8ed25c45a92822fe593862e543ae4167357cbc8fba4f24b889aa2bbf2729a",
"cross_section_rotate14.sbn": "3f7a3b66cf58be8c979353d2c75777303035e19ff58d96a089dde5c95fa8b597",
"cross_section_rotate14.sbx": "7d40bc92b42fde2af01a2805c9205c18c0fe89ae7cf1ba88ac6627b7c6a69b89",
"cross_section_rotate14.shp": "5f0ea7a65b5ddc9a43c874035969e30d58ae578aec9feb6b0e8538b68d5bd0d2",
"cross_section_rotate14.shp.xml": "79e38d9542ce764ace47883c673cf1d9aab16cd7851ae62a8e9bf27ce1091e13",
"cross_section_rotate14.shx": "b750b9d44ef31e0c593e2f78acfc08813667bb73733e6524f1b417e605cae65d",
"model_extent.cpg": "3ad3031f5503a4404af825262ee8232cc04d4ea6683d42c5dd0a2f2a27ac9824",
"model_extent.dbf": "72f8ed25c45a92822fe593862e543ae4167357cbc8fba4f24b889aa2bbf2729a",
"model_extent.sbn": "622376387ac9686e54acc6c57ace348c217d3a82e626274f32911a1d0006a164",
"model_extent.sbx": "2957bc1b5c918e20089fb6f6998d60d4488995d174bac21afa8e3a2af90b3489",
"model_extent.shp": "c72d5a4c703100e98c356c7645ad4b0bcc124c55e0757e55c8cd8663c7bf15c6",
"model_extent.shx": "e8d3b5618f0c248b59284f4f795f5de8207aec5b15ed60ce8da5a021c1043e2f",
"wells_locations.dbf": "965c846ec0b8f0d27570ef0bdaadfbcb6e718ed70ab89c8dda01d3b819e7a7de",
"wells_locations.sbn": "63f8ad670c6ba53ddec13069e42cfd86f27b6d47c5d0b3f2c25dfd6fb6b55825",
"wells_locations.sbx": "8420907d426c44c38315a5bdc0b24fe07a8cd2cc9a7fc60b817500b8cda79a34",
"wells_locations.shp": "ee53a4532b513f5b8bcd37ee3468dc4b2c8f6afab6cfc5110d74362c79e52287",
"wells_locations.shx": "6e816e96ed0726c2acc61392d2a82df5e9265ab5f5b00dd12f765b139840be79",
"wells_locations_rotate14.dbf": "d9b3636b4312c2f76c837e698bcb0d8ef9f4bbaa1765c484787a9f9d7f8bbaae",
"wells_locations_rotate14.sbn": "b436e34b8f145966b18d571b47ebc18e35671ec73fca1abbc737d9e1aa984bfb",
"wells_locations_rotate14.sbx": "24911f8905155882ce76b0330c9ba5ed449ca985d46833ebc45eee11faabbdaf",
"wells_locations_rotate14.shp": "695894af4678358320fb914e872cadb2613fae2e54c2d159f40c02fa558514cf",
"wells_locations_rotate14.shp.xml": "288183eb273c1fc2facb49d51c34bcafb16710189242da48f7717c49412f3e29",
"wells_locations_rotate14.shx": "da3374865cbf864f81dd69192ab616d1093d2159ac3c682fe2bfc4c295a28e42",
}
for fname, fhash in file_names.items():
pooch.retrieve(
url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/{sim_name}/gis/{fname}",
fname=fname,
path=modelpth / "gis",
known_hash=fhash,
)

# + pycharm={"name": "#%%\n"}
# Setup the figure and PlotMapView. Show a very faint map of ibound and
# model grid by specifying a transparency alpha value.
Expand All @@ -444,21 +518,21 @@
mapview = flopy.plot.PlotMapView(model=ml, ax=ax)

# Plot a shapefile of
shp = os.path.join(loadpth, "gis", "bedrock_outcrop_hole")
shp = os.path.join(modelpth, "gis", "bedrock_outcrop_hole")
patch_collection = mapview.plot_shapefile(
shp,
edgecolor="green",
linewidths=2,
alpha=0.5, # facecolor='none',
)
# Plot a shapefile of a cross-section line
shp = os.path.join(loadpth, "gis", "cross_section")
shp = os.path.join(modelpth, "gis", "cross_section")
patch_collection = mapview.plot_shapefile(
shp, radius=0, lw=[3, 1.5], edgecolor=["red", "green"], facecolor="None"
)

# Plot a shapefile of well locations
shp = os.path.join(loadpth, "gis", "wells_locations")
shp = os.path.join(modelpth, "gis", "wells_locations")
patch_collection = mapview.plot_shapefile(shp, radius=100, facecolor="red")

# Plot the grid and boundary conditions over the top
Expand All @@ -483,19 +557,19 @@
mapview = flopy.plot.PlotMapView(model=ml)

# Plot a shapefile of
shp = os.path.join(loadpth, "gis", "bedrock_outcrop_hole_rotate14")
shp = os.path.join(modelpth, "gis", "bedrock_outcrop_hole_rotate14")
patch_collection = mapview.plot_shapefile(
shp,
edgecolor="green",
linewidths=2,
alpha=0.5, # facecolor='none',
)
# Plot a shapefile of a cross-section line
shp = os.path.join(loadpth, "gis", "cross_section_rotate14")
shp = os.path.join(modelpth, "gis", "cross_section_rotate14")
patch_collection = mapview.plot_shapefile(shp, lw=3, edgecolor="red")

# Plot a shapefile of well locations
shp = os.path.join(loadpth, "gis", "wells_locations_rotate14")
shp = os.path.join(modelpth, "gis", "wells_locations_rotate14")
patch_collection = mapview.plot_shapefile(shp, radius=100, facecolor="red")

# Plot the grid and boundary conditions over the top
Expand Down Expand Up @@ -527,16 +601,16 @@

# + pycharm={"name": "#%%\n"}
# lets extract some shapes from our shapefiles
shp = os.path.join(loadpth, "gis", "bedrock_outcrop_hole_rotate14")
shp = os.path.join(modelpth, "gis", "bedrock_outcrop_hole_rotate14")
with shapefile.Reader(shp) as r:
polygon_w_hole = [r.shape(0)]

shp = os.path.join(loadpth, "gis", "cross_section_rotate14")
shp = os.path.join(modelpth, "gis", "cross_section_rotate14")
with shapefile.Reader(shp) as r:
cross_section = r.shapes()

# Plot a shapefile of well locations
shp = os.path.join(loadpth, "gis", "wells_locations_rotate14")
shp = os.path.join(modelpth, "gis", "wells_locations_rotate14")
with shapefile.Reader(shp) as r:
wells = r.shapes()

Expand Down Expand Up @@ -575,14 +649,43 @@

# + pycharm={"name": "#%%\n"}
# load the Freyberg model into mf6-flopy and run the simulation
sim_name = "mfsim.nam"
sim_path = os.path.join("..", "..", "examples", "data", "mf6-freyberg")

sim_name = "mf6-freyberg"
sim_path = modelpth / "mf6"
file_names = {
"bot.asc": "3107f907cb027460fd40ffc16cb797a78babb31988c7da326c9f500fba855b62",
"description.txt": "94093335eec6a24711f86d4d217ccd5a7716dd9e01cb6b732bc7757d41675c09",
"freyberg.cbc": "c8ad843b1da753eb58cf6c462ac782faf0ca433d6dcb067742d8bd698db271e3",
"freyberg.chd": "d8b8ada8d3978daea1758b315be983b5ca892efc7d69bf6b367ceec31e0dd156",
"freyberg.dis": "cac230a207cc8483693f7ba8ae29ce40c049036262eac4cebe17a4e2347a8b30",
"freyberg.dis.grb": "c8c26fb1fa4b210208134b286d895397cf4b3131f66e1d9dda76338502c7e96a",
"freyberg.hds": "926a06411ca658a89db6b5686f51ddeaf5b74ced81239cab1d43710411ba5f5b",
"freyberg.ic": "6efb56ee9cdd704b9a76fb9efd6dae750facc5426b828713f2d2cf8d35194120",
"freyberg.ims": "6dddae087d85417e3cdaa13e7b24165afb7f9575ab68586f3adb6c1b2d023781",
"freyberg.nam": "cee9b7b000fe35d2df26e878d09d465250a39504f87516c897e3fa14dcda081e",
"freyberg.npf": "81104d3546045fff0eddf5059465e560b83b492fa5a5acad1907ce18c2b9c15f",
"freyberg.oc": "c0715acd75eabcc42c8c47260a6c1abd6c784350983f7e2e6009ddde518b80b8",
"freyberg.rch": "a6ec1e0eda14fd2cdf618a5c0243a9caf82686c69242b783410d5abbcf971954",
"freyberg.riv": "a8cafc8c317cbe2acbb43e2f0cfe1188cb2277a7a174aeb6f3e6438013de8088",
"freyberg.sto": "74d748c2f0adfa0a32ee3f2912115c8f35b91011995b70c1ec6ae1c627242c41",
"freyberg.tdis": "9965cbb17caf5b865ea41a4ec04bcb695fe15a38cb539425fdc00abbae385cbe",
"freyberg.wel": "f19847de455598de52c05a4be745698c8cb589e5acfb0db6ab1f06ded5ff9310",
"k11.asc": "b6a8aa46ef17f7f096d338758ef46e32495eb9895b25d687540d676744f02af5",
"mfsim.nam": "6b8d6d7a56c52fb2bff884b3979e3d2201c8348b4bbfd2b6b9752863cbc9975e",
"top.asc": "3ad2b131671b9faca7f74c1dd2b2f41875ab0c15027764021a89f9c95dccaa6a",
}
for fname, fhash in file_names.items():
pooch.retrieve(
url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/{sim_name}/{fname}",
fname=fname,
path=sim_path,
known_hash=fhash,
)

sim = flopy.mf6.MFSimulation.load(
sim_name=sim_name, version=vmf6, exe_name=exe_name_mf6, sim_ws=sim_path
sim_name="mfsim.nam", version=vmf6, exe_name=exe_name_mf6, sim_ws=sim_path
)

newpth = os.path.join(modelpth)
sim.set_sim_path(newpth)
sim.write_simulation()
success, buff = sim.run_simulation()
if not success:
Expand Down Expand Up @@ -662,14 +765,14 @@

# + pycharm={"name": "#%%\n"}
# get the specific discharge from the cell budget file
cbc_file = os.path.join(newpth, "freyberg.cbc")
cbc_file = os.path.join(sim_path, "freyberg.cbc")
cbc = flopy.utils.CellBudgetFile(cbc_file)
spdis = cbc.get_data(text="SPDIS")[0]

qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, ml6)

# get the head from the head file
head_file = os.path.join(newpth, "freyberg.hds")
head_file = os.path.join(sim_path, "freyberg.hds")
head = flopy.utils.HeadFile(head_file)
hdata = head.get_alldata()[0]

Expand Down Expand Up @@ -985,10 +1088,10 @@ def run_vertex_grid_example(ws):
run_vertex_grid_example(modelpth)

# check if model ran properly
modelpth = os.path.join(modelpth, "mp7_ex2", "mf6")
mp7modelpth = os.path.join(modelpth, "mp7_ex2", "mf6")
files = ["mp7p2.hds", "mp7p2.cbb"]
for f in files:
if os.path.isfile(os.path.join(modelpth, f)):
if os.path.isfile(os.path.join(mp7modelpth, f)):
msg = f"Output file located: {f}"
print(msg)
else:
Expand All @@ -1002,7 +1105,7 @@ def run_vertex_grid_example(ws):
sim_name=vertex_sim_name,
version=vmf6,
exe_name=exe_name_mf6,
sim_ws=modelpth,
sim_ws=mp7modelpth,
)
vertex_ml6 = vertex_sim.get_model("mp7p2")

Expand Down Expand Up @@ -1052,7 +1155,7 @@ def run_vertex_grid_example(ws):

# + pycharm={"name": "#%%\n"}
# get the head output for stress period 1 from the modflow6 head file
head = flopy.utils.HeadFile(os.path.join(modelpth, "mp7p2.hds"))
head = flopy.utils.HeadFile(os.path.join(mp7modelpth, "mp7p2.hds"))
hdata = head.get_alldata()[0, :, :, :]

fig = plt.figure(figsize=(12, 12))
Expand Down Expand Up @@ -1093,11 +1196,11 @@ def run_vertex_grid_example(ws):
# + pycharm={"name": "#%%\n"}
# load the MODPATH-7 results
mp_namea = "mp7p2a_mp"
fpth = os.path.join(modelpth, f"{mp_namea}.mppth")
fpth = os.path.join(mp7modelpth, f"{mp_namea}.mppth")
p = flopy.utils.PathlineFile(fpth)
p0 = p.get_alldata()

fpth = os.path.join(modelpth, f"{mp_namea}.timeseries")
fpth = os.path.join(mp7modelpth, f"{mp_namea}.timeseries")
ts = flopy.utils.TimeseriesFile(fpth)
ts0 = ts.get_alldata()

Expand Down Expand Up @@ -1128,7 +1231,7 @@ def run_vertex_grid_example(ws):

# + pycharm={"name": "#%%\n"}
cbb = flopy.utils.CellBudgetFile(
os.path.join(modelpth, "mp7p2.cbb"), precision="double"
os.path.join(mp7modelpth, "mp7p2.cbb"), precision="double"
)
spdis = cbb.get_data(text="SPDIS")[0]
qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, vertex_ml6)
Expand All @@ -1154,8 +1257,22 @@ def run_vertex_grid_example(ws):
# set up the notebook for unstructured grid plotting
from flopy.discretization import UnstructuredGrid

# this is a folder containing some unstructured grids
datapth = os.path.join("..", "..", "examples", "data", "unstructured")
datapth = modelpth / "unstructured"
file_names = {
"TriMesh_local.exp": "0be6a1a1743972ba98c9d9e63ac2e457813c0809bfbda120e09a97b04411a65e",
"TriMesh_usg.exp": "0b450f2b306253a7b2889796e7a4eea52159f509c7b28a1f65929008dd854e08",
"Trimesh_circle.exp": "1efb86bb77060dcec20e752e242076e3bd23046f5e47d20d948bcf4623b3deb7",
"headu.githds": "cbe94655d471470d931923f70c7548b161ea4c5a22333b7fab6e2255450cda89",
"ugrid_iverts.dat": "7e33ec7f7d1fdbeb6cb7bc8dbcdf35f262c82aaa38dc79b4fb3fe7b53f7c7c1b",
"ugrid_verts.dat": "59493b26c8969789bb5a06d999db7a2dac324bffee280925e123007c81e689c7",
}
for fname, fhash in file_names.items():
pooch.retrieve(
url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/unstructured/{fname}",
fname=fname,
path=datapth,
known_hash=fhash,
)


# simple functions to load vertices and incidence lists
Expand Down Expand Up @@ -1259,14 +1376,14 @@ def load_iverts(fname):

# + pycharm={"name": "#%%\n"}
# get the specific discharge from the cell budget file
cbc_file = os.path.join(newpth, "freyberg.cbc")
cbc_file = os.path.join(sim_path, "freyberg.cbc")
cbc = flopy.utils.CellBudgetFile(cbc_file)
spdis = cbc.get_data(text="SPDIS")[0]

qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, ml6)

# get the head from the head file
head_file = os.path.join(newpth, "freyberg.hds")
head_file = os.path.join(sim_path, "freyberg.hds")
head = flopy.utils.HeadFile(head_file)
hdata = head.get_alldata()[0]

Expand Down

0 comments on commit 4d57f72

Please sign in to comment.