From be229b68f3411b27c202c9f98887cfb6d1fd485d Mon Sep 17 00:00:00 2001 From: Joshua Larsen Date: Wed, 12 May 2021 12:59:35 -0700 Subject: [PATCH 1/4] fix(plot_pathline): update modpath intersection routines for multilayered grids * added test_pathline_plot_xc to t049_test.py --- autotest/t049_test.py | 322 ++++++++++++++++++------------------- flopy/plot/crosssection.py | 14 +- flopy/plot/plotutil.py | 80 ++++++--- 3 files changed, 226 insertions(+), 190 deletions(-) diff --git a/autotest/t049_test.py b/autotest/t049_test.py index 23b662a829..867e6894a5 100644 --- a/autotest/t049_test.py +++ b/autotest/t049_test.py @@ -8,20 +8,20 @@ try: import pymake except ImportError: - print("could not import pymake") + print('could not import pymake') pymake = False -cpth = os.path.join("temp", "t049") +cpth = os.path.join('temp', 't049') # delete the directory if it exists if os.path.isdir(cpth): shutil.rmtree(cpth) # make the directory os.makedirs(cpth) -mf2005_exe = "mf2005" +mf2005_exe = 'mf2005' v = flopy.which(mf2005_exe) -mpth_exe = "mp6" +mpth_exe = 'mp6' v2 = flopy.which(mpth_exe) rung = True @@ -30,8 +30,8 @@ def test_modpath(): - pth = os.path.join("..", "examples", "data", "freyberg") - mfnam = "freyberg.nam" + pth = os.path.join('..', 'examples', 'data', 'freyberg') + mfnam = 'freyberg.nam' if pymake: run = rung @@ -41,9 +41,8 @@ def test_modpath(): run = False lpth = pth - m = flopy.modflow.Modflow.load( - mfnam, model_ws=lpth, verbose=True, exe_name=mf2005_exe - ) + m = flopy.modflow.Modflow.load(mfnam, model_ws=lpth, verbose=True, + exe_name=mf2005_exe) assert m.load_fail is False if run: @@ -51,24 +50,18 @@ def test_modpath(): success, buff = m.run_model(silent=False) except: success = False - assert success, "modflow model run did not terminate successfully" + assert success, 'modflow model run did not terminate successfully' # create the forward modpath file - mpnam = "freybergmp" - mp = flopy.modpath.Modpath6( - mpnam, exe_name=mpth_exe, modflowmodel=m, model_ws=lpth - ) - mpbas = flopy.modpath.Modpath6Bas( - mp, - hnoflo=m.bas6.hnoflo, - hdry=m.lpf.hdry, - ibound=m.bas6.ibound.array, - prsity=0.2, - prsityCB=0.2, - ) - sim = mp.create_mpsim( - trackdir="forward", simtype="endpoint", packages="RCH" - ) + mpnam = 'freybergmp' + mp = flopy.modpath.Modpath6(mpnam, exe_name=mpth_exe, modflowmodel=m, + model_ws=lpth) + mpbas = flopy.modpath.Modpath6Bas(mp, hnoflo=m.bas6.hnoflo, + hdry=m.lpf.hdry, + ibound=m.bas6.ibound.array, prsity=0.2, + prsityCB=0.2) + sim = mp.create_mpsim(trackdir='forward', simtype='endpoint', + packages='RCH') # write forward particle track files mp.write_input() @@ -78,25 +71,18 @@ def test_modpath(): success, buff = mp.run_model(silent=False) except: success = False - assert success, ( - "forward modpath model run " + "did not terminate successfully" - ) - - mpnam = "freybergmpp" - mpp = flopy.modpath.Modpath6( - mpnam, exe_name=mpth_exe, modflowmodel=m, model_ws=lpth - ) - mpbas = flopy.modpath.Modpath6Bas( - mpp, - hnoflo=m.bas6.hnoflo, - hdry=m.lpf.hdry, - ibound=m.bas6.ibound.array, - prsity=0.2, - prsityCB=0.2, - ) - sim = mpp.create_mpsim( - trackdir="backward", simtype="pathline", packages="WEL" - ) + assert success, 'forward modpath model run ' + \ + 'did not terminate successfully' + + mpnam = 'freybergmpp' + mpp = flopy.modpath.Modpath6(mpnam, exe_name=mpth_exe, + modflowmodel=m, model_ws=lpth) + mpbas = flopy.modpath.Modpath6Bas(mpp, hnoflo=m.bas6.hnoflo, + hdry=m.lpf.hdry, + ibound=m.bas6.ibound.array, prsity=0.2, + prsityCB=0.2) + sim = mpp.create_mpsim(trackdir='backward', simtype='pathline', + packages='WEL') # write backward particle track files mpp.write_input() @@ -106,44 +92,41 @@ def test_modpath(): success, buff = mpp.run_model(silent=False) except: success = False - assert success, ( - "backward modpath model run " + "did not terminate successfully" - ) + assert success, 'backward modpath model run ' + \ + 'did not terminate successfully' # load modpath output files if run: endfile = os.path.join(lpth, mp.sim.endpoint_file) pthfile = os.path.join(lpth, mpp.sim.pathline_file) else: - endfile = os.path.join( - "..", "examples", "data", "mp6_examples", "freybergmp.gitmpend" - ) - pthfile = os.path.join( - "..", "examples", "data", "mp6_examples", "freybergmpp.gitmppth" - ) + endfile = os.path.join('..', 'examples', 'data', 'mp6_examples', + 'freybergmp.gitmpend') + pthfile = os.path.join('..', 'examples', 'data', 'mp6_examples', + 'freybergmpp.gitmppth') # load the endpoint data try: endobj = flopy.utils.EndpointFile(endfile) except: - assert False, "could not load endpoint file" + assert False, 'could not load endpoint file' ept = endobj.get_alldata() - assert ept.shape == (695,), "shape of endpoint file is not (695,)" + assert ept.shape == (695,), 'shape of endpoint file is not (695,)' # load the pathline data try: pthobj = flopy.utils.PathlineFile(pthfile) except: - assert False, "could not load pathline file" + assert False, 'could not load pathline file' plines = pthobj.get_alldata() - assert len(plines) == 576, "there are not 576 particle pathlines in file" + assert len(plines) == 576, 'there are not 576 particle pathlines in file' return def test_pathline_plot(): - pth = os.path.join("..", "examples", "data", "freyberg") - mfnam = "freyberg.nam" + pth = os.path.join('..', 'examples', 'data', 'freyberg') + mfnam = 'freyberg.nam' run = rung try: @@ -154,39 +137,35 @@ def test_pathline_plot(): nampath = os.path.join(lpth, mfnam) assert os.path.exists(nampath), "namefile {} doesn't exist.".format( - nampath - ) + nampath) # load the modflow files for model map - m = flopy.modflow.Modflow.load( - mfnam, model_ws=lpth, verbose=True, forgive=False, exe_name=mf2005_exe - ) + m = flopy.modflow.Modflow.load(mfnam, model_ws=lpth, verbose=True, + forgive=False, + exe_name=mf2005_exe) # load modpath output files if run: - pthfile = os.path.join(lpth, "freybergmpp.mppth") + pthfile = os.path.join(lpth, 'freybergmpp.mppth') else: - pthfile = os.path.join( - "..", "examples", "data", "mp6_examples", "freybergmpp.gitmppth" - ) + pthfile = os.path.join('..', 'examples', 'data', 'mp6_examples', + 'freybergmpp.gitmppth') # load the pathline data try: pthobj = flopy.utils.PathlineFile(pthfile) except: - assert False, "could not load pathline file" + assert False, 'could not load pathline file' # determine version ver = pthobj.version - assert ver == 6, "{} is not a MODPATH version 6 pathline file".format( - pthfile - ) + assert ver == 6, '{} is not a MODPATH version 6 pathline file'.format(pthfile) # get all pathline data plines = pthobj.get_alldata() mm = flopy.plot.PlotMapView(model=m) try: - mm.plot_pathline(plines, colors="blue", layer="all") + mm.plot_pathline(plines, colors='blue', layer='all') except: assert False, 'could not plot pathline with layer="all"' @@ -195,83 +174,127 @@ def test_pathline_plot(): mm.plot_grid() mm.plot_ibound() except: - assert False, "could not plot grid and ibound" + assert False, 'could not plot grid and ibound' try: - fpth = os.path.join(lpth, "pathline.png") + fpth = os.path.join(lpth, 'pathline.png') plt.savefig(fpth) plt.close() except: - assert False, "could not save plot as {}".format(fpth) + assert False, 'could not save plot as {}'.format(fpth) mm = flopy.plot.PlotMapView(model=m) try: - mm.plot_pathline(plines, colors="green", layer=0) + mm.plot_pathline(plines, colors='green', layer=0) except: - assert False, "could not plot pathline with layer=0" + assert False, 'could not plot pathline with layer=0' # plot the grid and ibound array try: mm.plot_grid() mm.plot_ibound() except: - assert False, "could not plot grid and ibound" + assert False, 'could not plot grid and ibound' try: - fpth = os.path.join(lpth, "pathline2.png") + fpth = os.path.join(lpth, 'pathline2.png') plt.savefig(fpth) plt.close() except: - assert False, "could not save plot as {}".format(fpth) + assert False, 'could not save plot as {}'.format(fpth) mm = flopy.plot.PlotMapView(model=m) try: - mm.plot_pathline(plines, colors="red") + mm.plot_pathline(plines, colors='red') except: - assert False, "could not plot pathline" + assert False, 'could not plot pathline' # plot the grid and ibound array try: mm.plot_grid() mm.plot_ibound() except: - assert False, "could not plot grid and ibound" + assert False, 'could not plot grid and ibound' try: - fpth = os.path.join(lpth, "pathline3.png") + fpth = os.path.join(lpth, 'pathline3.png') plt.savefig(fpth) plt.close() except: - assert False, "could not save plot as {}".format(fpth) + assert False, 'could not save plot as {}'.format(fpth) return +def test_pathline_plot_xc(): + from matplotlib.collections import LineCollection + + # test with multi-layer example + model_ws = os.path.join("..", 'examples', 'data', 'mp6') + + ml = flopy.modflow.Modflow.load("EXAMPLE.nam", model_ws=model_ws) + ml.change_model_ws(os.path.join('.', 'temp')) + ml.write_input() + ml.run_model() + + mp = flopy.modpath.Modpath6(modelname='ex6', + exe_name='mp6', + modflowmodel=ml, + model_ws=os.path.join('.', 'temp'), + dis_file=ml.name + '.DIS', + head_file=ml.name + '.hed', + budget_file=ml.name + '.bud') + + mpb = flopy.modpath.Modpath6Bas(mp, hdry=ml.lpf.hdry, laytyp=ml.lpf.laytyp, + ibound=1, prsity=0.1) + + sim = mp.create_mpsim(trackdir='forward', simtype='pathline', + packages='RCH', start_time=(2, 0, 1.)) + mp.write_input() + + mp.run_model(silent=False) + + pthobj = flopy.utils.PathlineFile(os.path.join('temp', 'ex6.mppth')) + well_pathlines = pthobj.get_destination_pathline_data( + dest_cells=[(4, 12, 12)]) + + mx = flopy.plot.PlotCrossSection(model=ml, line={'row': 4}) + mx.plot_bc("WEL", kper=2, color='blue') + pth = mx.plot_pathline(well_pathlines, method='cell', colors='red') + + if not isinstance(pth, LineCollection): + raise AssertionError() + + if len(pth._paths) != 6: + raise AssertionError() + + plt.close() + + def test_mp5_load(): # load the base freyberg model - pth = os.path.join("..", "examples", "data", "freyberg") + pth = os.path.join('..', 'examples', 'data', 'freyberg') # load the modflow files for model map - m = flopy.modflow.Modflow.load( - "freyberg.nam", model_ws=pth, check=False, verbose=True, forgive=False - ) + m = flopy.modflow.Modflow.load('freyberg.nam', model_ws=pth, check=False, + verbose=True, forgive=False) # load the pathline data - fpth = os.path.join("..", "examples", "data", "mp5", "m.ptl") + fpth = os.path.join('..', 'examples', 'data', 'mp5', 'm.ptl') try: pthobj = flopy.utils.PathlineFile(fpth) except: - assert False, "could not load pathline file" + assert False, 'could not load pathline file' # load endpoint data - fpth = os.path.join("..", "examples", "data", "mp5", "m.ept") + fpth = os.path.join('..', 'examples', 'data', 'mp5', 'm.ept') try: endobj = flopy.utils.EndpointFile(fpth, verbose=True) except: - assert False, "could not load endpoint file" + assert False, 'could not load endpoint file' # determine version ver = pthobj.version - assert ver == 5, "{} is not a MODPATH version 5 pathline file".format(fpth) + assert ver == 5, '{} is not a MODPATH version 5 pathline file'.format(fpth) # read all of the pathline and endpoint data plines = pthobj.get_alldata() @@ -279,9 +302,9 @@ def test_mp5_load(): # determine the number of particles in the pathline file nptl = pthobj.nid.shape[0] - assert nptl == 64, "number of MODPATH 5 particles does not equal 64" + assert nptl == 64, 'number of MODPATH 5 particles does not equal 64' - hsv = plt.get_cmap("hsv") + hsv = plt.get_cmap('hsv') colors = hsv(np.linspace(0, 1.0, nptl)) # plot the pathlines one pathline at a time @@ -290,44 +313,37 @@ def test_mp5_load(): p = pthobj.get_data(partid=n) e = endobj.get_data(partid=n) try: - mm.plot_pathline(p, colors=colors[n], layer="all") + mm.plot_pathline(p, colors=colors[n], layer='all') except: - assert False, ( - "could not plot pathline {} ".format(n + 1) - + 'with layer="all"' - ) + assert False, 'could not plot pathline {} '.format(n + 1) + \ + 'with layer="all"' try: mm.plot_endpoint(e) except: - assert False, ( - "could not plot endpoint {} ".format(n + 1) - + 'with layer="all"' - ) + assert False, 'could not plot endpoint {} '.format(n + 1) + \ + 'with layer="all"' # plot the grid and ibound array try: mm.plot_grid(lw=0.5) mm.plot_ibound() except: - assert False, "could not plot grid and ibound" + assert False, 'could not plot grid and ibound' try: - fpth = os.path.join(cpth, "mp5.pathline.png") + fpth = os.path.join(cpth, 'mp5.pathline.png') plt.savefig(fpth, dpi=300) plt.close() except: - assert False, "could not save plot as {}".format(fpth) + assert False, 'could not save plot as {}'.format(fpth) return def test_mp5_timeseries_load(): - pth = os.path.join("..", "examples", "data", "mp5") - files = [ - os.path.join(pth, name) - for name in sorted(os.listdir(pth)) - if ".timeseries" in name - ] + pth = os.path.join('..', 'examples', 'data', 'mp5') + files = [os.path.join(pth, name) for name in sorted(os.listdir(pth)) + if '.timeseries' in name] for file in files: print(file) eval_timeseries(file) @@ -335,12 +351,9 @@ def test_mp5_timeseries_load(): def test_mp6_timeseries_load(): - pth = os.path.join("..", "examples", "data", "mp6") - files = [ - os.path.join(pth, name) - for name in sorted(os.listdir(pth)) - if ".timeseries" in name - ] + pth = os.path.join('..', 'examples', 'data', 'mp6') + files = [os.path.join(pth, name) for name in sorted(os.listdir(pth)) + if '.timeseries' in name] for file in files: print(file) eval_timeseries(file) @@ -349,10 +362,8 @@ def test_mp6_timeseries_load(): def eval_timeseries(file): ts = flopy.utils.TimeseriesFile(file) - msg = ( - "{} ".format(os.path.basename(file)) - + "is not an instance of flopy.utils.TimeseriesFile" - ) + msg = '{} '.format(os.path.basename(file)) + \ + 'is not an instance of flopy.utils.TimeseriesFile' assert isinstance(ts, flopy.utils.TimeseriesFile), msg # get the all of the data @@ -360,9 +371,8 @@ def eval_timeseries(file): tsd = ts.get_alldata() except: pass - msg = "could not load data using get_alldata() from " + "{}.".format( - os.path.basename(file) - ) + msg = 'could not load data using get_alldata() from ' + \ + '{}.'.format(os.path.basename(file)) assert len(tsd) > 0, msg # get the data for the last particleid @@ -371,47 +381,37 @@ def eval_timeseries(file): partid = ts.get_maxid() except: pass - msg = ( - "could not get maximum particleid using get_maxid() from " - + "{}.".format(os.path.basename(file)) - ) + msg = 'could not get maximum particleid using get_maxid() from ' + \ + '{}.'.format(os.path.basename(file)) assert partid is not None, msg try: tsd = ts.get_data(partid=partid) except: pass - msg = ( - "could not load data for particleid {} ".format(partid) - + "using get_data() from " - + "{}. ".format(os.path.basename(file)) - + "Maximum partid = " - + "{}.".format(ts.get_maxid()) - ) + msg = 'could not load data for particleid {} '.format(partid) + \ + 'using get_data() from ' + \ + '{}. '.format(os.path.basename(file)) + 'Maximum partid = ' + \ + '{}.'.format(ts.get_maxid()) assert tsd.shape[0] > 0, msg timemax = None try: - timemax = ts.get_maxtime() / 2.0 + timemax = ts.get_maxtime() / 2. except: pass - msg = ( - "could not get maximum time using get_maxtime() from " - + "{}.".format(os.path.basename(file)) - ) + msg = 'could not get maximum time using get_maxtime() from ' + \ + '{}.'.format(os.path.basename(file)) assert timemax is not None, msg try: tsd = ts.get_alldata(totim=timemax) except: pass - msg = ( - "could not load data for totim>={} ".format(timemax) - + "using get_alldata() from " - + "{}. ".format(os.path.basename(file)) - + "Maximum totim = " - + "{}.".format(ts.get_maxtime()) - ) + msg = 'could not load data for totim>={} '.format(timemax) + \ + 'using get_alldata() from ' + \ + '{}. '.format(os.path.basename(file)) + 'Maximum totim = ' + \ + '{}.'.format(ts.get_maxtime()) assert len(tsd) > 0, msg timemax = None @@ -419,31 +419,27 @@ def eval_timeseries(file): timemax = ts.get_maxtime() except: pass - msg = ( - "could not get maximum time using get_maxtime() from " - + "{}.".format(os.path.basename(file)) - ) + msg = 'could not get maximum time using get_maxtime() from ' + \ + '{}.'.format(os.path.basename(file)) assert timemax is not None, msg try: tsd = ts.get_alldata(totim=timemax, ge=False) except: pass - msg = ( - "could not load data for totim<={} ".format(timemax) - + "using get_alldata() from " - + "{}. ".format(os.path.basename(file)) - + "Maximum totim = " - + "{}.".format(ts.get_maxtime()) - ) + msg = 'could not load data for totim<={} '.format(timemax) + \ + 'using get_alldata() from ' + \ + '{}. '.format(os.path.basename(file)) + 'Maximum totim = ' + \ + '{}.'.format(ts.get_maxtime()) assert len(tsd) > 0, msg return -if __name__ == "__main__": +if __name__ == '__main__': test_modpath() test_pathline_plot() + test_pathline_plot_xc() test_mp5_load() test_mp5_timeseries_load() test_mp6_timeseries_load() diff --git a/flopy/plot/crosssection.py b/flopy/plot/crosssection.py index 13ee99a207..a6b99c3873 100644 --- a/flopy/plot/crosssection.py +++ b/flopy/plot/crosssection.py @@ -631,7 +631,7 @@ def plot_inactive(self, ibound=None, color_noflow="black", **kwargs): else: ibound = self.mg.idomain - plotarray = np.zeros(ibound.shape, dtype=int) + plotarray = np.zeros(ibound.shape, dtype=np.int) idx1 = ibound == 0 plotarray[idx1] = 1 plotarray = np.ma.masked_equal(plotarray, 0) @@ -685,7 +685,7 @@ def plot_ibound( ibound = self.mg.idomain - plotarray = np.zeros(ibound.shape, dtype=int) + plotarray = np.zeros(ibound.shape, dtype=np.int) idx1 = ibound == 0 idx2 = ibound < 0 plotarray[idx1] = 1 @@ -823,11 +823,11 @@ def plot_bc( idx = mflist["node"] if len(self.mg.shape) != 3: - plotarray = np.zeros((self._nlay, self._ncpl), dtype=int) + plotarray = np.zeros((self._nlay, self._ncpl), dtype=np.int) plotarray[tuple(idx)] = 1 else: plotarray = np.zeros( - (self.mg.nlay, self.mg.nrow, self.mg.ncol), dtype=int + (self.mg.nlay, self.mg.nrow, self.mg.ncol), dtype=np.int ) plotarray[idx[0], idx[1], idx[2]] = 1 @@ -1155,7 +1155,7 @@ def plot_discharge( # thickness by setting laytyp to zeros if head is None or laytyp is None: head = np.zeros(botm.shape, np.float32) - laytyp = np.zeros((nlay,), dtype=int) + laytyp = np.zeros((nlay,), dtype=np.int) head[0, :, :] = top if nlay > 1: head[1:, :, :] = botm[:-1, :, :] @@ -1170,7 +1170,7 @@ def plot_discharge( ) if qz is None: - qz = np.zeros(qx.shape, dtype=float) + qz = np.zeros(qx.shape, dtype=np.float) qx = qx.ravel() qy = qy.ravel() @@ -1258,6 +1258,7 @@ def plot_pathline( self.xvertices, self.yvertices, self.direction, + self._ncpl, method=method, ) plines = plotutil.reproject_modpath_to_crosssection( @@ -1266,6 +1267,7 @@ def plot_pathline( self.xypts, self.direction, self.mg, + self._ncpl, self.geographic_coords, ) diff --git a/flopy/plot/plotutil.py b/flopy/plot/plotutil.py index 781f3f5ad1..ed15c3525f 100644 --- a/flopy/plot/plotutil.py +++ b/flopy/plot/plotutil.py @@ -2712,19 +2712,35 @@ def intersect_modpath_with_crosssection( xvertices, yvertices, projection, + ncpl, method="cell", starting=False, ): """ + Method to intersect modpath output with a cross-section - :param recarrays: - :param projpts: - :param xvertices: - :param yvertices: - :param projection: - :param method: - :param starting: - :return: + Parameters + ---------- + recarrays : list + list of numpy recarrays + projpts : dict + dict of crossectional cell vertices + xvertices : np.array + array of modelgrid xvertices + yvertices : np.array + array of modelgrid yvertices + projection : str + projection direction (x or y) + ncpl : int + number of cells per layer (cross sectional version) + method : str + intersection method ('cell' or 'all') + starting : bool + modpath starting location flag + + Returns + ------- + dict : dictionary of intersecting recarrays """ from ..utils.geometry import point_in_polygon @@ -2754,12 +2770,15 @@ def intersect_modpath_with_crosssection( nppts = {} for cell, verts in projpts.items(): + tcell = cell + while tcell >= ncpl: + tcell -= ncpl zmin = np.min(np.array(verts)[:, 1]) zmax = np.max(np.array(verts)[:, 1]) - nmin = np.min(v_norm[cell]) - nmax = np.max(v_norm[cell]) - omin = np.min(v_opp[cell]) - omax = np.max(v_opp[cell]) + nmin = np.min(v_norm[tcell]) + nmax = np.max(v_norm[tcell]) + omin = np.min(v_opp[tcell]) + omax = np.max(v_opp[tcell]) oppts[cell] = np.array( [ [omin, zmax], @@ -2818,19 +2837,35 @@ def reproject_modpath_to_crosssection( xypts, projection, modelgrid, + ncpl, geographic_coords, starting=False, ): """ + Method to reproject modpath points onto cross sectional line - :param idict: - :param projpts: - :param xypts: - :param projection: - :param modelgrid: - :param geographic_coords: - :param starting: - :return: + Parameters + ---------- + idict : dict + dictionary of intersecting points + projpts : dict + dictionary of cross sectional cells + xypts : dict + dictionary of cross sectional line + projection : str + projection direction (x or y) + modelgrid : Grid object + flopy modelgrid object + ncpl : int + number of cells per layer (cross sectional version) + geographic_coords : bool + flag for plotting in geographic coordinates + starting : bool + flag for modpath position + + Returns + ------- + dictionary of projected modpath lines or points """ from ..utils import geometry @@ -2845,7 +2880,10 @@ def reproject_modpath_to_crosssection( ptdict = {} if not geographic_coords: for cell, recarrays in idict.items(): - line = xypts[cell] + tcell = cell + while tcell >= ncpl: + tcell -= ncpl + line = xypts[tcell] if projection == "x": d0 = np.min([i[0] for i in projpts[cell]]) else: From 3d9eb23102dab9b2d5ff9cecbb89037c6374b963 Mon Sep 17 00:00:00 2001 From: Joshua Larsen Date: Wed, 12 May 2021 14:10:03 -0700 Subject: [PATCH 2/4] updates to crosssection.py --- flopy/plot/crosssection.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/flopy/plot/crosssection.py b/flopy/plot/crosssection.py index a6b99c3873..06d8a90999 100644 --- a/flopy/plot/crosssection.py +++ b/flopy/plot/crosssection.py @@ -631,7 +631,7 @@ def plot_inactive(self, ibound=None, color_noflow="black", **kwargs): else: ibound = self.mg.idomain - plotarray = np.zeros(ibound.shape, dtype=np.int) + plotarray = np.zeros(ibound.shape, dtype=int) idx1 = ibound == 0 plotarray[idx1] = 1 plotarray = np.ma.masked_equal(plotarray, 0) @@ -685,7 +685,7 @@ def plot_ibound( ibound = self.mg.idomain - plotarray = np.zeros(ibound.shape, dtype=np.int) + plotarray = np.zeros(ibound.shape, dtype=int) idx1 = ibound == 0 idx2 = ibound < 0 plotarray[idx1] = 1 @@ -823,11 +823,11 @@ def plot_bc( idx = mflist["node"] if len(self.mg.shape) != 3: - plotarray = np.zeros((self._nlay, self._ncpl), dtype=np.int) + plotarray = np.zeros((self._nlay, self._ncpl), dtype=int) plotarray[tuple(idx)] = 1 else: plotarray = np.zeros( - (self.mg.nlay, self.mg.nrow, self.mg.ncol), dtype=np.int + (self.mg.nlay, self.mg.nrow, self.mg.ncol), dtype=int ) plotarray[idx[0], idx[1], idx[2]] = 1 @@ -1155,7 +1155,7 @@ def plot_discharge( # thickness by setting laytyp to zeros if head is None or laytyp is None: head = np.zeros(botm.shape, np.float32) - laytyp = np.zeros((nlay,), dtype=np.int) + laytyp = np.zeros((nlay,), dtype=int) head[0, :, :] = top if nlay > 1: head[1:, :, :] = botm[:-1, :, :] @@ -1170,7 +1170,7 @@ def plot_discharge( ) if qz is None: - qz = np.zeros(qx.shape, dtype=np.float) + qz = np.zeros(qx.shape, dtype=float) qx = qx.ravel() qy = qy.ravel() From 7b6b28216c40cde8b9a420e335eb0972aa027151 Mon Sep 17 00:00:00 2001 From: Joshua Larsen Date: Wed, 12 May 2021 14:10:30 -0700 Subject: [PATCH 3/4] updates to t049_test.py --- autotest/t049_test.py | 276 +++++++++++++++++++++++++----------------- 1 file changed, 163 insertions(+), 113 deletions(-) diff --git a/autotest/t049_test.py b/autotest/t049_test.py index 867e6894a5..d4794ef6ed 100644 --- a/autotest/t049_test.py +++ b/autotest/t049_test.py @@ -8,20 +8,20 @@ try: import pymake except ImportError: - print('could not import pymake') + print("could not import pymake") pymake = False -cpth = os.path.join('temp', 't049') +cpth = os.path.join("temp", "t049") # delete the directory if it exists if os.path.isdir(cpth): shutil.rmtree(cpth) # make the directory os.makedirs(cpth) -mf2005_exe = 'mf2005' +mf2005_exe = "mf2005" v = flopy.which(mf2005_exe) -mpth_exe = 'mp6' +mpth_exe = "mp6" v2 = flopy.which(mpth_exe) rung = True @@ -30,8 +30,8 @@ def test_modpath(): - pth = os.path.join('..', 'examples', 'data', 'freyberg') - mfnam = 'freyberg.nam' + pth = os.path.join("..", "examples", "data", "freyberg") + mfnam = "freyberg.nam" if pymake: run = rung @@ -41,8 +41,9 @@ def test_modpath(): run = False lpth = pth - m = flopy.modflow.Modflow.load(mfnam, model_ws=lpth, verbose=True, - exe_name=mf2005_exe) + m = flopy.modflow.Modflow.load( + mfnam, model_ws=lpth, verbose=True, exe_name=mf2005_exe + ) assert m.load_fail is False if run: @@ -50,18 +51,24 @@ def test_modpath(): success, buff = m.run_model(silent=False) except: success = False - assert success, 'modflow model run did not terminate successfully' + assert success, "modflow model run did not terminate successfully" # create the forward modpath file - mpnam = 'freybergmp' - mp = flopy.modpath.Modpath6(mpnam, exe_name=mpth_exe, modflowmodel=m, - model_ws=lpth) - mpbas = flopy.modpath.Modpath6Bas(mp, hnoflo=m.bas6.hnoflo, - hdry=m.lpf.hdry, - ibound=m.bas6.ibound.array, prsity=0.2, - prsityCB=0.2) - sim = mp.create_mpsim(trackdir='forward', simtype='endpoint', - packages='RCH') + mpnam = "freybergmp" + mp = flopy.modpath.Modpath6( + mpnam, exe_name=mpth_exe, modflowmodel=m, model_ws=lpth + ) + mpbas = flopy.modpath.Modpath6Bas( + mp, + hnoflo=m.bas6.hnoflo, + hdry=m.lpf.hdry, + ibound=m.bas6.ibound.array, + prsity=0.2, + prsityCB=0.2, + ) + sim = mp.create_mpsim( + trackdir="forward", simtype="endpoint", packages="RCH" + ) # write forward particle track files mp.write_input() @@ -71,18 +78,25 @@ def test_modpath(): success, buff = mp.run_model(silent=False) except: success = False - assert success, 'forward modpath model run ' + \ - 'did not terminate successfully' - - mpnam = 'freybergmpp' - mpp = flopy.modpath.Modpath6(mpnam, exe_name=mpth_exe, - modflowmodel=m, model_ws=lpth) - mpbas = flopy.modpath.Modpath6Bas(mpp, hnoflo=m.bas6.hnoflo, - hdry=m.lpf.hdry, - ibound=m.bas6.ibound.array, prsity=0.2, - prsityCB=0.2) - sim = mpp.create_mpsim(trackdir='backward', simtype='pathline', - packages='WEL') + assert success, ( + "forward modpath model run " + "did not terminate successfully" + ) + + mpnam = "freybergmpp" + mpp = flopy.modpath.Modpath6( + mpnam, exe_name=mpth_exe, modflowmodel=m, model_ws=lpth + ) + mpbas = flopy.modpath.Modpath6Bas( + mpp, + hnoflo=m.bas6.hnoflo, + hdry=m.lpf.hdry, + ibound=m.bas6.ibound.array, + prsity=0.2, + prsityCB=0.2, + ) + sim = mpp.create_mpsim( + trackdir="backward", simtype="pathline", packages="WEL" + ) # write backward particle track files mpp.write_input() @@ -92,41 +106,44 @@ def test_modpath(): success, buff = mpp.run_model(silent=False) except: success = False - assert success, 'backward modpath model run ' + \ - 'did not terminate successfully' + assert success, ( + "backward modpath model run " + "did not terminate successfully" + ) # load modpath output files if run: endfile = os.path.join(lpth, mp.sim.endpoint_file) pthfile = os.path.join(lpth, mpp.sim.pathline_file) else: - endfile = os.path.join('..', 'examples', 'data', 'mp6_examples', - 'freybergmp.gitmpend') - pthfile = os.path.join('..', 'examples', 'data', 'mp6_examples', - 'freybergmpp.gitmppth') + endfile = os.path.join( + "..", "examples", "data", "mp6_examples", "freybergmp.gitmpend" + ) + pthfile = os.path.join( + "..", "examples", "data", "mp6_examples", "freybergmpp.gitmppth" + ) # load the endpoint data try: endobj = flopy.utils.EndpointFile(endfile) except: - assert False, 'could not load endpoint file' + assert False, "could not load endpoint file" ept = endobj.get_alldata() - assert ept.shape == (695,), 'shape of endpoint file is not (695,)' + assert ept.shape == (695,), "shape of endpoint file is not (695,)" # load the pathline data try: pthobj = flopy.utils.PathlineFile(pthfile) except: - assert False, 'could not load pathline file' + assert False, "could not load pathline file" plines = pthobj.get_alldata() - assert len(plines) == 576, 'there are not 576 particle pathlines in file' + assert len(plines) == 576, "there are not 576 particle pathlines in file" return def test_pathline_plot(): - pth = os.path.join('..', 'examples', 'data', 'freyberg') - mfnam = 'freyberg.nam' + pth = os.path.join("..", "examples", "data", "freyberg") + mfnam = "freyberg.nam" run = rung try: @@ -137,35 +154,39 @@ def test_pathline_plot(): nampath = os.path.join(lpth, mfnam) assert os.path.exists(nampath), "namefile {} doesn't exist.".format( - nampath) + nampath + ) # load the modflow files for model map - m = flopy.modflow.Modflow.load(mfnam, model_ws=lpth, verbose=True, - forgive=False, - exe_name=mf2005_exe) + m = flopy.modflow.Modflow.load( + mfnam, model_ws=lpth, verbose=True, forgive=False, exe_name=mf2005_exe + ) # load modpath output files if run: - pthfile = os.path.join(lpth, 'freybergmpp.mppth') + pthfile = os.path.join(lpth, "freybergmpp.mppth") else: - pthfile = os.path.join('..', 'examples', 'data', 'mp6_examples', - 'freybergmpp.gitmppth') + pthfile = os.path.join( + "..", "examples", "data", "mp6_examples", "freybergmpp.gitmppth" + ) # load the pathline data try: pthobj = flopy.utils.PathlineFile(pthfile) except: - assert False, 'could not load pathline file' + assert False, "could not load pathline file" # determine version ver = pthobj.version - assert ver == 6, '{} is not a MODPATH version 6 pathline file'.format(pthfile) + assert ver == 6, "{} is not a MODPATH version 6 pathline file".format( + pthfile + ) # get all pathline data plines = pthobj.get_alldata() mm = flopy.plot.PlotMapView(model=m) try: - mm.plot_pathline(plines, colors='blue', layer='all') + mm.plot_pathline(plines, colors="blue", layer="all") except: assert False, 'could not plot pathline with layer="all"' @@ -174,54 +195,54 @@ def test_pathline_plot(): mm.plot_grid() mm.plot_ibound() except: - assert False, 'could not plot grid and ibound' + assert False, "could not plot grid and ibound" try: - fpth = os.path.join(lpth, 'pathline.png') + fpth = os.path.join(lpth, "pathline.png") plt.savefig(fpth) plt.close() except: - assert False, 'could not save plot as {}'.format(fpth) + assert False, "could not save plot as {}".format(fpth) mm = flopy.plot.PlotMapView(model=m) try: - mm.plot_pathline(plines, colors='green', layer=0) + mm.plot_pathline(plines, colors="green", layer=0) except: - assert False, 'could not plot pathline with layer=0' + assert False, "could not plot pathline with layer=0" # plot the grid and ibound array try: mm.plot_grid() mm.plot_ibound() except: - assert False, 'could not plot grid and ibound' + assert False, "could not plot grid and ibound" try: - fpth = os.path.join(lpth, 'pathline2.png') + fpth = os.path.join(lpth, "pathline2.png") plt.savefig(fpth) plt.close() except: - assert False, 'could not save plot as {}'.format(fpth) + assert False, "could not save plot as {}".format(fpth) mm = flopy.plot.PlotMapView(model=m) try: - mm.plot_pathline(plines, colors='red') + mm.plot_pathline(plines, colors="red") except: - assert False, 'could not plot pathline' + assert False, "could not plot pathline" # plot the grid and ibound array try: mm.plot_grid() mm.plot_ibound() except: - assert False, 'could not plot grid and ibound' + assert False, "could not plot grid and ibound" try: - fpth = os.path.join(lpth, 'pathline3.png') + fpth = os.path.join(lpth, "pathline3.png") plt.savefig(fpth) plt.close() except: - assert False, 'could not save plot as {}'.format(fpth) + assert False, "could not save plot as {}".format(fpth) return @@ -273,28 +294,29 @@ def test_pathline_plot_xc(): def test_mp5_load(): # load the base freyberg model - pth = os.path.join('..', 'examples', 'data', 'freyberg') + pth = os.path.join("..", "examples", "data", "freyberg") # load the modflow files for model map - m = flopy.modflow.Modflow.load('freyberg.nam', model_ws=pth, check=False, - verbose=True, forgive=False) + m = flopy.modflow.Modflow.load( + "freyberg.nam", model_ws=pth, check=False, verbose=True, forgive=False + ) # load the pathline data - fpth = os.path.join('..', 'examples', 'data', 'mp5', 'm.ptl') + fpth = os.path.join("..", "examples", "data", "mp5", "m.ptl") try: pthobj = flopy.utils.PathlineFile(fpth) except: - assert False, 'could not load pathline file' + assert False, "could not load pathline file" # load endpoint data - fpth = os.path.join('..', 'examples', 'data', 'mp5', 'm.ept') + fpth = os.path.join("..", "examples", "data", "mp5", "m.ept") try: endobj = flopy.utils.EndpointFile(fpth, verbose=True) except: - assert False, 'could not load endpoint file' + assert False, "could not load endpoint file" # determine version ver = pthobj.version - assert ver == 5, '{} is not a MODPATH version 5 pathline file'.format(fpth) + assert ver == 5, "{} is not a MODPATH version 5 pathline file".format(fpth) # read all of the pathline and endpoint data plines = pthobj.get_alldata() @@ -302,9 +324,9 @@ def test_mp5_load(): # determine the number of particles in the pathline file nptl = pthobj.nid.shape[0] - assert nptl == 64, 'number of MODPATH 5 particles does not equal 64' + assert nptl == 64, "number of MODPATH 5 particles does not equal 64" - hsv = plt.get_cmap('hsv') + hsv = plt.get_cmap("hsv") colors = hsv(np.linspace(0, 1.0, nptl)) # plot the pathlines one pathline at a time @@ -313,37 +335,44 @@ def test_mp5_load(): p = pthobj.get_data(partid=n) e = endobj.get_data(partid=n) try: - mm.plot_pathline(p, colors=colors[n], layer='all') + mm.plot_pathline(p, colors=colors[n], layer="all") except: - assert False, 'could not plot pathline {} '.format(n + 1) + \ - 'with layer="all"' + assert False, ( + "could not plot pathline {} ".format(n + 1) + + 'with layer="all"' + ) try: mm.plot_endpoint(e) except: - assert False, 'could not plot endpoint {} '.format(n + 1) + \ - 'with layer="all"' + assert False, ( + "could not plot endpoint {} ".format(n + 1) + + 'with layer="all"' + ) # plot the grid and ibound array try: mm.plot_grid(lw=0.5) mm.plot_ibound() except: - assert False, 'could not plot grid and ibound' + assert False, "could not plot grid and ibound" try: - fpth = os.path.join(cpth, 'mp5.pathline.png') + fpth = os.path.join(cpth, "mp5.pathline.png") plt.savefig(fpth, dpi=300) plt.close() except: - assert False, 'could not save plot as {}'.format(fpth) + assert False, "could not save plot as {}".format(fpth) return def test_mp5_timeseries_load(): - pth = os.path.join('..', 'examples', 'data', 'mp5') - files = [os.path.join(pth, name) for name in sorted(os.listdir(pth)) - if '.timeseries' in name] + pth = os.path.join("..", "examples", "data", "mp5") + files = [ + os.path.join(pth, name) + for name in sorted(os.listdir(pth)) + if ".timeseries" in name + ] for file in files: print(file) eval_timeseries(file) @@ -351,9 +380,12 @@ def test_mp5_timeseries_load(): def test_mp6_timeseries_load(): - pth = os.path.join('..', 'examples', 'data', 'mp6') - files = [os.path.join(pth, name) for name in sorted(os.listdir(pth)) - if '.timeseries' in name] + pth = os.path.join("..", "examples", "data", "mp6") + files = [ + os.path.join(pth, name) + for name in sorted(os.listdir(pth)) + if ".timeseries" in name + ] for file in files: print(file) eval_timeseries(file) @@ -362,8 +394,10 @@ def test_mp6_timeseries_load(): def eval_timeseries(file): ts = flopy.utils.TimeseriesFile(file) - msg = '{} '.format(os.path.basename(file)) + \ - 'is not an instance of flopy.utils.TimeseriesFile' + msg = ( + "{} ".format(os.path.basename(file)) + + "is not an instance of flopy.utils.TimeseriesFile" + ) assert isinstance(ts, flopy.utils.TimeseriesFile), msg # get the all of the data @@ -371,8 +405,9 @@ def eval_timeseries(file): tsd = ts.get_alldata() except: pass - msg = 'could not load data using get_alldata() from ' + \ - '{}.'.format(os.path.basename(file)) + msg = "could not load data using get_alldata() from " + "{}.".format( + os.path.basename(file) + ) assert len(tsd) > 0, msg # get the data for the last particleid @@ -381,37 +416,47 @@ def eval_timeseries(file): partid = ts.get_maxid() except: pass - msg = 'could not get maximum particleid using get_maxid() from ' + \ - '{}.'.format(os.path.basename(file)) + msg = ( + "could not get maximum particleid using get_maxid() from " + + "{}.".format(os.path.basename(file)) + ) assert partid is not None, msg try: tsd = ts.get_data(partid=partid) except: pass - msg = 'could not load data for particleid {} '.format(partid) + \ - 'using get_data() from ' + \ - '{}. '.format(os.path.basename(file)) + 'Maximum partid = ' + \ - '{}.'.format(ts.get_maxid()) + msg = ( + "could not load data for particleid {} ".format(partid) + + "using get_data() from " + + "{}. ".format(os.path.basename(file)) + + "Maximum partid = " + + "{}.".format(ts.get_maxid()) + ) assert tsd.shape[0] > 0, msg timemax = None try: - timemax = ts.get_maxtime() / 2. + timemax = ts.get_maxtime() / 2.0 except: pass - msg = 'could not get maximum time using get_maxtime() from ' + \ - '{}.'.format(os.path.basename(file)) + msg = ( + "could not get maximum time using get_maxtime() from " + + "{}.".format(os.path.basename(file)) + ) assert timemax is not None, msg try: tsd = ts.get_alldata(totim=timemax) except: pass - msg = 'could not load data for totim>={} '.format(timemax) + \ - 'using get_alldata() from ' + \ - '{}. '.format(os.path.basename(file)) + 'Maximum totim = ' + \ - '{}.'.format(ts.get_maxtime()) + msg = ( + "could not load data for totim>={} ".format(timemax) + + "using get_alldata() from " + + "{}. ".format(os.path.basename(file)) + + "Maximum totim = " + + "{}.".format(ts.get_maxtime()) + ) assert len(tsd) > 0, msg timemax = None @@ -419,24 +464,29 @@ def eval_timeseries(file): timemax = ts.get_maxtime() except: pass - msg = 'could not get maximum time using get_maxtime() from ' + \ - '{}.'.format(os.path.basename(file)) + msg = ( + "could not get maximum time using get_maxtime() from " + + "{}.".format(os.path.basename(file)) + ) assert timemax is not None, msg try: tsd = ts.get_alldata(totim=timemax, ge=False) except: pass - msg = 'could not load data for totim<={} '.format(timemax) + \ - 'using get_alldata() from ' + \ - '{}. '.format(os.path.basename(file)) + 'Maximum totim = ' + \ - '{}.'.format(ts.get_maxtime()) + msg = ( + "could not load data for totim<={} ".format(timemax) + + "using get_alldata() from " + + "{}. ".format(os.path.basename(file)) + + "Maximum totim = " + + "{}.".format(ts.get_maxtime()) + ) assert len(tsd) > 0, msg return -if __name__ == '__main__': +if __name__ == "__main__": test_modpath() test_pathline_plot() test_pathline_plot_xc() From b8dbd5bd56e0d1dc9705084df9661f063d91a548 Mon Sep 17 00:00:00 2001 From: Joshua Larsen Date: Wed, 12 May 2021 15:25:45 -0700 Subject: [PATCH 4/4] update(test_pathline_plot_xc): update exe_name for unix compatibility --- autotest/t049_test.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/autotest/t049_test.py b/autotest/t049_test.py index d4794ef6ed..e4597d9980 100644 --- a/autotest/t049_test.py +++ b/autotest/t049_test.py @@ -253,13 +253,14 @@ def test_pathline_plot_xc(): # test with multi-layer example model_ws = os.path.join("..", 'examples', 'data', 'mp6') - ml = flopy.modflow.Modflow.load("EXAMPLE.nam", model_ws=model_ws) + ml = flopy.modflow.Modflow.load("EXAMPLE.nam", model_ws=model_ws, + exe_name=mf2005_exe) ml.change_model_ws(os.path.join('.', 'temp')) ml.write_input() ml.run_model() mp = flopy.modpath.Modpath6(modelname='ex6', - exe_name='mp6', + exe_name=mpth_exe, modflowmodel=ml, model_ws=os.path.join('.', 'temp'), dis_file=ml.name + '.DIS',