fix(plot_pathline): update modpath intersection routines (#1112)

* fix(plot_pathline): update modpath intersection routines for multilayered grids

* added test_pathline_plot_xc to t049_test.py

* updates to crosssection.py

* updates to t049_test.py

* update(test_pathline_plot_xc): update exe_name for unix compatibility

autotest/t049_test.py

@@ -247,6 +247,52 @@ def test_pathline_plot():
     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,
+                                    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=mpth_exe,
+                                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")
@@ -444,6 +490,7 @@ def eval_timeseries(file):
 if __name__ == "__main__":
     test_modpath()
     test_pathline_plot()
+    test_pathline_plot_xc()
     test_mp5_load()
     test_mp5_timeseries_load()
     test_mp6_timeseries_load()

flopy/plot/crosssection.py

@@ -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,
         )
 

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: