feat(CellBudgetFile): add get_pathlines method to CellBudgetFile

autotest/test_binaryfile.py

@@ -1,4 +1,5 @@
 import os
+from typing import List
 
 import numpy as np
 import pytest
@@ -545,3 +546,14 @@ def test_cellbudgetfile_readrecord_waux(example_data_path):
                 record
             )
     v.close()
+
+
+def test_cellbudgetfile_get_pathlines(example_data_path):
+    prt_model_path = example_data_path / "prt_test"
+    cbb_fname = prt_model_path / "mp7-p04_prt.cbb"
+    cbb = CellBudgetFile(cbb_fname)
+    assert isinstance(cbb, CellBudgetFile)
+
+    pls = cbb.get_pathlines()
+    assert isinstance(pls, list)
+    assert all(isinstance(pl, np.recarray) for pl in pls)

flopy/utils/binaryfile.py

@@ -11,7 +11,7 @@
 import os
 import warnings
 from pathlib import Path
-from typing import Union
+from typing import List, Union
 
 import numpy as np
 
@@ -1365,7 +1365,7 @@ def get_data(
         text=None,
         paknam=None,
         full3D=False,
-    ):
+    ) -> Union[List, np.ndarray]:
         """
         Get data from the binary budget file.
 
@@ -1963,6 +1963,56 @@ def get_residual(self, totim, scaled=False):
 
         return residual
 
+    def get_pathlines(self, prpnam=None):
+        dtype = {
+            "names": ["x", "y", "z", "time", "k", "particleid"],
+            "formats": ["<f4", "<f4", "<f4", "<f4", "<i4", "<i4"],
+        }
+        times = self.get_times()
+
+        # Extract pathline data
+        pdict = {}
+        # For each time, read and parse particle data
+        # and add it to a dictionary keyed on a particle identifier consisting
+        # of release point number and release time
+        for totim in times:
+            data = self.get_data(text="DATA-PRTCL", paknam=prpnam, totim=totim)
+            for ploc in data[0]:
+                # particle data read from the binary output file are in terms
+                # of one-based indexing; will convert to zero-based indexing
+                irpt, node, xp, yp, zp, trelease, ttrack = [
+                    ploc[i] for i in (0, 1, 3, 4, 5, 6, 7)
+                ]
+                irpt -= 1
+                node -= 1
+
+                # todo: how to retrieve ncpl here unless it's passed as argument?
+                #  do we need to include k (layer number) in the pathpoint tuple?
+                # kp = node // ncpl
+
+                # 2nd to last element was kp in Alden's original implementation
+                # last element: issue(mf6): particle ids not assigned yet in mf6
+                pathpoint = (xp, yp, zp, ttrack, 0, 0)
+
+                key = (irpt, trelease)
+                # Note that repeat locations written for terminated particles
+                # are not filtered out here, though they could be
+                if key not in pdict:
+                    pdict[key] = [pathpoint]
+                else:
+                    pdict[key].append(pathpoint)
+
+        # Reformat into mf6pathlines array
+        mf6pathlines = []
+        for key in pdict:
+            path = np.core.records.fromarrays(
+                np.array(pdict[key]).T, dtype=dtype
+            )
+            mf6pathlines.append(path)
+
+        # Return mf6pathlines
+        return mf6pathlines
+
     def close(self):
         """
         Close the file handle