From 4d57f72e05a8793c10d9555a29b4c38e1c1a3aca Mon Sep 17 00:00:00 2001 From: w-bonelli Date: Tue, 10 Dec 2024 15:34:12 -0500 Subject: [PATCH] map view example --- .docs/Notebooks/plot_map_view_example.py | 181 +++++++++++++++++++---- 1 file changed, 149 insertions(+), 32 deletions(-) diff --git a/.docs/Notebooks/plot_map_view_example.py b/.docs/Notebooks/plot_map_view_example.py index 78e3cb961..17362f16c 100644 --- a/.docs/Notebooks/plot_map_view_example.py +++ b/.docs/Notebooks/plot_map_view_example.py @@ -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 @@ -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) @@ -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. @@ -444,7 +518,7 @@ 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", @@ -452,13 +526,13 @@ 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 @@ -483,7 +557,7 @@ 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", @@ -491,11 +565,11 @@ 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 @@ -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() @@ -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: @@ -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] @@ -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: @@ -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") @@ -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)) @@ -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() @@ -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) @@ -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 @@ -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]