From d7edbc61febda718254a33bbfc14e39ed8059173 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Tue, 12 May 2026 13:46:38 -0700 Subject: [PATCH 01/13] update(intersect): update days calculation to fix dropped precision * replace timedelta.days with timedelta.total_seconds() / 86400 --- flopy/discretization/modeltime.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/flopy/discretization/modeltime.py b/flopy/discretization/modeltime.py index f072d11070..2e92453031 100644 --- a/flopy/discretization/modeltime.py +++ b/flopy/discretization/modeltime.py @@ -618,15 +618,14 @@ def intersect( "objects, set time units or use totim for intersection" ) - elif self.time_units == "days": - totim = timedelta.days - - elif self.time_units in {"hours", "minutes", "seconds"}: + elif self.time_units in {"days", "hours", "minutes", "seconds"}: totim = timedelta.total_seconds() if self.time_units == "minutes": totim /= 60 elif self.time_units == "hours": totim /= 3600 + elif self.time_units == "days": + totim /= 86400. else: # years condition @@ -647,7 +646,7 @@ def intersect( ) timedelta = datetime_obj - dt_iyear - days = timedelta.days + days = timedelta.total_seconds() / 86400. yr_frac = days / ndays totim += yr_frac From 8d2811b25f59946e4a659f9581b05eff4f1054d0 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Tue, 12 May 2026 14:36:36 -0700 Subject: [PATCH 02/13] linting --- flopy/discretization/modeltime.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/flopy/discretization/modeltime.py b/flopy/discretization/modeltime.py index 2e92453031..2056c0110f 100644 --- a/flopy/discretization/modeltime.py +++ b/flopy/discretization/modeltime.py @@ -625,7 +625,7 @@ def intersect( elif self.time_units == "hours": totim /= 3600 elif self.time_units == "days": - totim /= 86400. + totim /= 86400.0 else: # years condition @@ -646,7 +646,7 @@ def intersect( ) timedelta = datetime_obj - dt_iyear - days = timedelta.total_seconds() / 86400. + days = timedelta.total_seconds() / 86400.0 yr_frac = days / ndays totim += yr_frac From 523991459a32ebc7d4dbd26208e4ff731588cc89 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Tue, 12 May 2026 15:57:06 -0700 Subject: [PATCH 03/13] update test_modeltime.py, remove note in hfb_util.py --- autotest/test_modeltime.py | 4 ++-- flopy/utils/hfb_util.py | 1 - 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/autotest/test_modeltime.py b/autotest/test_modeltime.py index 966d21d6be..fbdb54e335 100644 --- a/autotest/test_modeltime.py +++ b/autotest/test_modeltime.py @@ -270,7 +270,7 @@ def test_mf2005_modeltime(): raise AssertionError("start_datetime improperly stored") result = modeltime.intersect("3/06/2024 23:59:59") - if result != (2, 0): + if result != (2, 1): raise AssertionError("ModelTime intersect not working correctly") @@ -317,7 +317,7 @@ def test_mf6_modeltime(): raise AssertionError("start_datetime improperly stored") result = modeltime.intersect("3/06/2024 23:59:59") - if result != (2, 0): + if result != (2, 1): raise AssertionError("ModelTime intersect not working correctly") diff --git a/flopy/utils/hfb_util.py b/flopy/utils/hfb_util.py index 97d6c83983..036d33c87d 100644 --- a/flopy/utils/hfb_util.py +++ b/flopy/utils/hfb_util.py @@ -232,7 +232,6 @@ def make_hfb_array(modelgrid, geom): for record in result: node = record.cellid ixshp = record.ixshapes - # todo: numpy array this and Transpose coords = np.array(ixshp.coords.xy).T x0, y0 = coords[0, 0], coords[0, 1] x1, y1 = coords[-1, 0], coords[-1, 1] From f8c61255924e3d06e8e4de723b0811d78eece77f Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Wed, 13 May 2026 15:13:26 -0700 Subject: [PATCH 04/13] update(ModflowWel): add basic support for tabfile based well files --- flopy/modflow/mfwel.py | 46 ++++++++++++++++++++++++++++++------------ flopy/pakbase.py | 13 ++++++++++++ 2 files changed, 46 insertions(+), 13 deletions(-) diff --git a/flopy/modflow/mfwel.py b/flopy/modflow/mfwel.py index 86c815c373..d0a68ced0a 100644 --- a/flopy/modflow/mfwel.py +++ b/flopy/modflow/mfwel.py @@ -295,24 +295,44 @@ def add_record(self, kper, index, values): raise Exception(f"mfwel error adding record to list: {e!s}") @staticmethod - def get_default_dtype(structured=True): - if structured: - dtype = np.dtype( - [ - ("k", int), - ("i", int), - ("j", int), - ("flux", np.float32), - ] - ) + def get_default_dtype(structured=True, tabfiles=False): + if not tabfiles: + if structured: + dtype = np.dtype( + [ + ("k", int), + ("i", int), + ("j", int), + ("flux", np.float32), + ] + ) + else: + dtype = np.dtype([("node", int), ("flux", np.float32)]) else: - dtype = np.dtype([("node", int), ("flux", np.float32)]) + if structured: + dtype = np.dtype( + [ + ("tabunit", int), + ("tabval", int), + ("k", int), + ("i", int), + ("j", int) + ] + ) + else: + dtype = np.dtype( + [ + ("tabunit", int), + ("tabval", int), + ("node", int) + ] + ) return dtype @staticmethod - def get_empty(ncells=0, aux_names=None, structured=True): + def get_empty(ncells=0, aux_names=None, structured=True, tabfiles=False): # get an empty recarray that corresponds to dtype - dtype = ModflowWel.get_default_dtype(structured=structured) + dtype = ModflowWel.get_default_dtype(structured=structured, tabfiles=tabfiles) if aux_names is not None: dtype = Package.add_to_dtype(dtype, aux_names, np.float32) return create_empty_recarray(ncells, dtype, default_value=-1.0e10) diff --git a/flopy/pakbase.py b/flopy/pakbase.py index b9bb0d9fa1..7cd85358e9 100644 --- a/flopy/pakbase.py +++ b/flopy/pakbase.py @@ -961,6 +961,8 @@ def load( nwt_options = OptionBlock.load_options(f, pak_type) line = f.readline() + nwt_tabfiles = False + # check for parameters nppak = 0 if "parameter" in line.lower(): @@ -1062,6 +1064,11 @@ def load( options = nwt_options else: f.seek(ipos) + + if options: + if options.tabfiles is not None: + nwt_tabfiles = True + elif "flopy.modflow.mfchd.modflowchd".lower() in pak_type_str: partype = ["shead", "ehead"] @@ -1085,6 +1092,12 @@ def load( bnd_output_cln = None stress_period_data_cln = {} current_cln = None + + if nwt_tabfiles: + # pass tabfile flag using the existing usg_args dict, change nper to 1 + nper = 1 + usg_args["tabfiles"] = True + for iper in range(nper): if model.verbose: msg = f" loading {pak_type} for kper {iper + 1:5d}" From 2a9dd95f612c408924fa9392fc862acd13455666 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Wed, 13 May 2026 15:17:15 -0700 Subject: [PATCH 05/13] Linting --- flopy/modflow/mfwel.py | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/flopy/modflow/mfwel.py b/flopy/modflow/mfwel.py index d0a68ced0a..be48526291 100644 --- a/flopy/modflow/mfwel.py +++ b/flopy/modflow/mfwel.py @@ -316,17 +316,11 @@ def get_default_dtype(structured=True, tabfiles=False): ("tabval", int), ("k", int), ("i", int), - ("j", int) + ("j", int), ] ) else: - dtype = np.dtype( - [ - ("tabunit", int), - ("tabval", int), - ("node", int) - ] - ) + dtype = np.dtype([("tabunit", int), ("tabval", int), ("node", int)]) return dtype @staticmethod From 9f1cde243c2a767477b509467f84bd8112f6be22 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Wed, 13 May 2026 15:25:10 -0700 Subject: [PATCH 06/13] remove "flux" sfac_column for tabfile wel records --- flopy/pakbase.py | 1 + 1 file changed, 1 insertion(+) diff --git a/flopy/pakbase.py b/flopy/pakbase.py index 7cd85358e9..29c81ae518 100644 --- a/flopy/pakbase.py +++ b/flopy/pakbase.py @@ -1096,6 +1096,7 @@ def load( if nwt_tabfiles: # pass tabfile flag using the existing usg_args dict, change nper to 1 nper = 1 + sfac_columns = [] usg_args["tabfiles"] = True for iper in range(nper): From 8f3f63586c4b266db4527eb7b617d960301ea322 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Wed, 13 May 2026 15:31:35 -0700 Subject: [PATCH 07/13] update tabfile option check --- flopy/pakbase.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/flopy/pakbase.py b/flopy/pakbase.py index 29c81ae518..d34a6cfe0a 100644 --- a/flopy/pakbase.py +++ b/flopy/pakbase.py @@ -1066,7 +1066,7 @@ def load( f.seek(ipos) if options: - if options.tabfiles is not None: + if options.tabfiles: nwt_tabfiles = True elif "flopy.modflow.mfchd.modflowchd".lower() in pak_type_str: From 1367c9d8b59ffdd0c57466ce5c5522fbe4b5d1b4 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Wed, 13 May 2026 15:52:36 -0700 Subject: [PATCH 08/13] update aux variable check for tabfile type list records --- flopy/modflow/mfwel.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/flopy/modflow/mfwel.py b/flopy/modflow/mfwel.py index be48526291..be63248196 100644 --- a/flopy/modflow/mfwel.py +++ b/flopy/modflow/mfwel.py @@ -163,6 +163,7 @@ def __init__( self.url = "wel.html" self.np = 0 + tabfiles = False if options is None: options = [] self.specify = False @@ -178,6 +179,9 @@ def __init__( self.phiramp = self.options.phiramp self.iunitramp = self.options.iunitramp # this is to grab the aux variables... + if self.options.tabfiles: + tabfiles = True + options = [] else: @@ -196,7 +200,9 @@ def __init__( self.dtype = self.get_default_dtype(structured=self.parent.structured) # determine if any aux variables in dtype - dt = self.get_default_dtype(structured=self.parent.structured) + dt = self.get_default_dtype( + structured=self.parent.structured, tabfiles=tabfiles + ) if len(self.dtype.names) > len(dt.names): for name in self.dtype.names[len(dt.names) :]: ladd = True From 11ec71b1dcba13d20572c4e9cae0836e31198994 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Thu, 14 May 2026 08:59:05 -0700 Subject: [PATCH 09/13] update option block check for tabfiles --- flopy/pakbase.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/flopy/pakbase.py b/flopy/pakbase.py index d34a6cfe0a..64c744e944 100644 --- a/flopy/pakbase.py +++ b/flopy/pakbase.py @@ -1065,7 +1065,7 @@ def load( else: f.seek(ipos) - if options: + if isinstance(options, OptionBlock): if options.tabfiles: nwt_tabfiles = True From fa10b00bccf34882882d14db337e65da214d0ec3 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Thu, 16 Jul 2026 09:35:40 -0700 Subject: [PATCH 10/13] Updates for multiband rasters and unsigned integer arrays --- flopy/utils/rasters.py | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/flopy/utils/rasters.py b/flopy/utils/rasters.py index 1c8989d895..02600712b2 100644 --- a/flopy/utils/rasters.py +++ b/flopy/utils/rasters.py @@ -81,7 +81,10 @@ def __init__( elif array.dtype in Raster.FLOAT64: dtype = "float64" elif array.dtype in Raster.INT8: - dtype = "int8" + if np.max(array) > 127: + dtype = "uint8" + else: + dtype = "int8" elif array.dtype in Raster.INT16: dtype = "int16" elif array.dtype in Raster.INT32: @@ -996,6 +999,8 @@ def plot(self, ax=None, contour=False, **kwargs): import_optional_dependency("rasterio") from rasterio.plot import show + band = kwargs.pop("band", None) + if self._dataset is not None: ax = show(self._dataset, ax=ax, contour=contour, **kwargs) @@ -1008,11 +1013,15 @@ def plot(self, ax=None, contour=False, **kwargs): if d1 is None: raise AssertionError("No plottable arrays found") - data = np.zeros((d0, d1, d2), dtype=float) - i = 0 - for _, arr in sorted(self.__arr_dict.items()): - data[i, :, :] = arr - i += 1 + if band is not None: + data = np.zeros((1, d1, d2), dtype=float) + data[0] = self.__arr_dict[band] + else: + data = np.zeros((d0, d1, d2), dtype=float) + i = 0 + for _, arr in sorted(self.__arr_dict.items()): + data[i, :, :] = arr + i += 1 data = np.ma.masked_where(data == self.nodatavals, data) ax = show( From b62d42c60adbe0ccff91a88fb927d7a1f69aab5e Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Thu, 16 Jul 2026 14:08:01 -0700 Subject: [PATCH 11/13] update(PlotCrossSection): add additional flexibility to cross sectional plotting * new feature: view, users can optionally provide `view="x"` or `view="y"` to force the view angle of the cross section * add in initial and end line points to the intersection routine which can include an additional cell on each side of cross section * provide uint support for Raster class --- flopy/plot/crosssection.py | 11 ++++++++++- flopy/plot/plotutil.py | 10 ++++++++++ flopy/utils/rasters.py | 4 ++-- 3 files changed, 22 insertions(+), 3 deletions(-) diff --git a/flopy/plot/crosssection.py b/flopy/plot/crosssection.py index e50ad5ed87..6d8c820465 100644 --- a/flopy/plot/crosssection.py +++ b/flopy/plot/crosssection.py @@ -51,6 +51,13 @@ class PlotCrossSection: minimum width of a grid cell polygon to be plotted. Cells with a cross-sectional width less than min_segment_length will be ignored and not included in the plot. Default is 1e-02. + view : str + view can be used to force the view of the cross section when a line is provided. + "auto" is default and mimics the long term behavior of flopy, which decides + on the view by taking the maximum of the x and y direction of the cross + sectional line. "x" forces the view to be plotted from the x direction + (bottom of the unrotated grid). and "y" forces the view to be plotted from the + y-direction "left" or "right" side of the unrotated grid. """ def __init__( @@ -62,7 +69,9 @@ def __init__( extent=None, geographic_coords=False, min_segment_length=1e-02, + view="auto" ): + view = view.lower() self.ax = ax self.geographic_coords = geographic_coords self.model = model @@ -159,7 +168,7 @@ def __init__( yp.append(v2) xp, yp = self.mg.get_local_coords(xp, yp) - if np.max(xp) - np.min(xp) > np.max(yp) - np.min(yp): + if (np.max(xp) - np.min(xp) > np.max(yp) - np.min(yp)) or view not in ("auto", "y"): # this is x-projection and we should buffer x by small amount idx0 = np.argmax(xp) idx1 = np.argmin(xp) diff --git a/flopy/plot/plotutil.py b/flopy/plot/plotutil.py index f79b133212..e87eda7ad2 100644 --- a/flopy/plot/plotutil.py +++ b/flopy/plot/plotutil.py @@ -1685,6 +1685,15 @@ def line_intersect_grid(ptsin, xgrid, ygrid): if t: vdict[cell] = t + # check if line vertex 0 or -1 should be included in the crosssection + ptschk = [ptsin[0], ptsin[-1]] + for cell, vrts in vdict.items(): + if len(vrts) == 1: + for pt in ptschk: + if np.min(xgrid[cell]) < pt[0] < np.max(xgrid[cell]): + if np.min(ygrid[cell]) < pt[1] < np.max(ygrid[cell]): + vdict[cell] = [vrts[0], tuple(pt)] + return vdict @staticmethod @@ -1713,6 +1722,7 @@ def filter_line_segments(vdict, threshold=1e-2): for node in nodes: points = vdict[node] if len(points) < 2: + # todo: need to pass in cross section points.... dist = 0 else: pt0 = points[0] diff --git a/flopy/utils/rasters.py b/flopy/utils/rasters.py index c568866ab7..10489e3084 100644 --- a/flopy/utils/rasters.py +++ b/flopy/utils/rasters.py @@ -998,8 +998,8 @@ def plot(self, ax=None, contour=False, **kwargs): """ import_optional_dependency("rasterio") from rasterio.plot import show - - band = kwargs.pop("band", None) + + band = kwargs.pop("band", None) # rasterio >= 1.5 defaults to adjusting data to [0, 1], which # can clash with user-supplied vmin/vmax. don't adjust unless # the caller explicitly asks for it From 6da1ee2040d528224a50238018ff84a723e4d8fd Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Thu, 16 Jul 2026 14:10:26 -0700 Subject: [PATCH 12/13] linting --- flopy/plot/crosssection.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/flopy/plot/crosssection.py b/flopy/plot/crosssection.py index 6d8c820465..5ae25e5182 100644 --- a/flopy/plot/crosssection.py +++ b/flopy/plot/crosssection.py @@ -69,7 +69,7 @@ def __init__( extent=None, geographic_coords=False, min_segment_length=1e-02, - view="auto" + view="auto", ): view = view.lower() self.ax = ax @@ -168,7 +168,10 @@ def __init__( yp.append(v2) xp, yp = self.mg.get_local_coords(xp, yp) - if (np.max(xp) - np.min(xp) > np.max(yp) - np.min(yp)) or view not in ("auto", "y"): + if (np.max(xp) - np.min(xp) > np.max(yp) - np.min(yp)) or view not in ( + "auto", + "y", + ): # this is x-projection and we should buffer x by small amount idx0 = np.argmax(xp) idx1 = np.argmin(xp) From 07b83ddf36a36c149937c577527249df311ace82 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Fri, 17 Jul 2026 07:49:38 -0700 Subject: [PATCH 13/13] remove todo note --- flopy/plot/plotutil.py | 1 - 1 file changed, 1 deletion(-) diff --git a/flopy/plot/plotutil.py b/flopy/plot/plotutil.py index e87eda7ad2..23f6d00f3c 100644 --- a/flopy/plot/plotutil.py +++ b/flopy/plot/plotutil.py @@ -1722,7 +1722,6 @@ def filter_line_segments(vdict, threshold=1e-2): for node in nodes: points = vdict[node] if len(points) < 2: - # todo: need to pass in cross section points.... dist = 0 else: pt0 = points[0]