From 4fcecc9bbc9bd7c2377ddf196569623d2ccac903 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Tue, 5 Aug 2025 15:47:52 -0700 Subject: [PATCH 01/13] Checking in progress on geodataframe support throughout code base --- flopy/mf6/data/mfdataarray.py | 48 ++++++++++++++++++++ flopy/mf6/data/mfdatalist.py | 72 ++++++++++++++++++++++++++++++ flopy/pakbase.py | 38 ++++++++++++++++ flopy/utils/util_array.py | 84 +++++++++++++++++++++++++++++++++++ flopy/utils/util_list.py | 53 ++++++++++++++++++++++ 5 files changed, 295 insertions(+) diff --git a/flopy/mf6/data/mfdataarray.py b/flopy/mf6/data/mfdataarray.py index e00c877b1..f4c2b6f09 100644 --- a/flopy/mf6/data/mfdataarray.py +++ b/flopy/mf6/data/mfdataarray.py @@ -323,6 +323,54 @@ def data(self): """Returns array data. Calls get_data with default parameters.""" return self._get_data() + def to_geo_dataframe(self, gdf=None): + """ + Method to add data to a GeoDataFrame for exporting as a geospatial file + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + + Returns + ------- + GeoDataFrame + """ + if self.model is None: + return gdf + else: + modelgrid = self.model.modelgrid + if gdf is None: + if modelgrid is None: + return gdf + gdf = modelgrid.geo_dataframe + + if modelgrid is not None: + if modelgrid.grid_type() != "unstructured": + ncpl = modelgrid.ncpl + else: + ncpl = modelgrid.nnodes + else: + ncpl = len(gdf) + + name = self.name + data = self.array + if data.size == ncpl: + gdf[name] = data.ravel() + + elif data.size % ncpl == 0: + data = data.reshape((-1, ncpl)) + for ix, dat in enumerate(data): + aname = f"{name}_{ix}" + gdf[aname] = dat + else: + raise ValueError( + f"Data size {data.size} not compatible with dataframe length {ncpl}" + ) + + return gdf + def new_simulation(self, sim_data): """Initialize MFArray object for a new simulation diff --git a/flopy/mf6/data/mfdatalist.py b/flopy/mf6/data/mfdatalist.py index 529177c41..22e5490c1 100644 --- a/flopy/mf6/data/mfdatalist.py +++ b/flopy/mf6/data/mfdatalist.py @@ -145,6 +145,78 @@ def to_array(self, kper=0, mask=False): model_grid = self.data_dimensions.get_model_grid() return list_to_array(sarr, model_grid, kper, mask) + def to_geo_dataframe(self, gdf=None, sparse=False): + """ + Method to add data to a GeoDataFrame for exporting as a geospatial file + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + sparse : bool + boolean flag for sparse dataframe construction. Default is False + + Returns + ------- + GeoDataFrame + """ + import pandas as pd + + if self.model is None: + return gdf + else: + modelgrid = self.model.modelgrid + if modelgrid is None: + return gdf + + if gdf is None: + gdf = modelgrid.geo_dataframe + + # name = self.name + recarray = self.array + if "cellid" not in recarray.dtype.names: + return gdf + + cellid = recarray["cellid"] + df = pd.DataFrame.from_records(recarray) + + grid_type = modelgrid.grid_type() + if grid_type != "unstructured": + layer = [cid[0] for cid in cellid] + if grid_type == "structured": + cids = [(0,) + cid[1:] for cid in cellid] + nodes = modelgrid.get_node(cids) + else: + nodes = [cid[1] for cid in cellid] + + df["layer"] = layer + df["node"] = nodes + df = df.drop(columns=["cellid", ]) + for layer in df.layer.unique(): + tmp = df[df.layer == layer] + tmp = tmp.drop(columns=["layer"]) + renames = { + nm: f"{nm}_{layer}" for nm in list(df) if nm not in ("node",) + } + tmp = tmp.set_index("node") + tmp = tmp.rename(columns=renames) + gdf = pd.merge(gdf, tmp, how="left", left_index=True, right_index=True) + + else: + nodes = [cid[0] for cid in cellid] + df["node"] = nodes + df = df.drop(columns=["cellid",]) + df = df.set_index("node") + gdf = pd.merge(gdf, df, how="left", left_index=True, right_index=True) + + if sparse: + cols = [col for col in list(df) if col != "node"] + if cols: + gdf = gdf[~np.isnan(gdf[cols[0]])] + + return gdf + def new_simulation(self, sim_data): """Initialize MFList object for a new simulation. diff --git a/flopy/pakbase.py b/flopy/pakbase.py index 8873254a6..dca6e8e89 100644 --- a/flopy/pakbase.py +++ b/flopy/pakbase.py @@ -673,6 +673,44 @@ def export(self, f, **kwargs): return export.utils.package_export(f, self, **kwargs) + def to_geo_dataframe(self, gdf=None, kper=0): + """ + Method to create a GeoDataFrame from a modflow package + + Parameters + ---------- + gdf : GeoDataFrame + optional geopandas geodataframe object to add data to. Default is None + kper : int + stress period to get transient data from + + Returns + ------- + gdf : GeoDataFrame + """ + from .mbase import BaseModel + if gdf is None: + if isinstance(self.parent, BaseModel): + modelgrid = self.parent.modelgrid + if modelgrid is not None: + gdf = modelgrid.geo_dataframe + else: + raise AttributeError( + "model does not have a grid instance, " + "please supply a geodataframe" + ) + else: + raise AssertionError( + "Package does not have a model instance, " + "please supply a geodataframe" + ) + + for attr, value in self.__dict__.items(): + if callable(getattr(value, "to_geo_dataframe", None)): + gdf = value.to_geo_dataframe(gdf, forgive=True, kper=kper) + + return gdf + def _generate_heading(self): """Generate heading.""" from . import __version__ diff --git a/flopy/utils/util_array.py b/flopy/utils/util_array.py index b83025a90..01f336e39 100644 --- a/flopy/utils/util_array.py +++ b/flopy/utils/util_array.py @@ -633,6 +633,40 @@ def data_type(self): def plottable(self): return True + def to_geo_dataframe(self, gdf=None, forgive=False, **kwargs): + """ + Method to add data to a GeoDataFrame for exporting as a geospatial file + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + forgive : bool + optional flag to continue running and pass data that is not compatible + with the geodataframe shape + + Returns + ------- + GeoDataFrame + """ + if self.model is None: + return gdf + else: + modelgrid = self.model.modelgrid + if modelgrid is None: + return gdf + + if gdf is None: + gdf = modelgrid.geo_dataframe + + names = self.name + for lay, u2d in enumerate(self.util_2ds): + name = f"{names[lay]}_{lay}" + gdf = u2d.to_geo_dataframe(gdf=gdf, name=name, forgive=forgive) + + return gdf + def export(self, f, **kwargs): from .. import export @@ -1875,6 +1909,56 @@ def _decide_how(self): else: self._how = "internal" + def to_geo_dataframe(self, gdf=None, name=None, forgive=False, **kwargs): + """ + Method to add an input array to a geopandas GeoDataFrame + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame object + name : str + optional attribute name, default uses util2d name + forgive : bool + optional flag to continue if data shape not compatible with GeoDataFrame + + Returns + ------- + geopandas GeoDataFrame + """ + if self.model is None: + return gdf + else: + modelgrid = self.model.modelgrid + if gdf is None: + if modelgrid is None: + return gdf + gdf = modelgrid.geo_dataframe + + if modelgrid is not None: + if modelgrid.grid_type != "unstructured": + ncpl = modelgrid.ncpl + else: + ncpl = modelgrid.nnodes + else: + ncpl = len(gdf) + + if name is None: + name = self.name + + data = self.array + + if data.size == ncpl: + gdf[name] = data.ravel() + elif forgive: + return gdf + else: + raise AssertionError( + f"Data size {data.size} not compatible with dataframe length {ncpl}" + ) + + return gdf + def plot( self, title=None, diff --git a/flopy/utils/util_list.py b/flopy/utils/util_list.py index da88d8e1b..edbf71edb 100644 --- a/flopy/utils/util_list.py +++ b/flopy/utils/util_list.py @@ -123,6 +123,59 @@ def get_empty(self, ncell=0): d = create_empty_recarray(ncell, self.dtype, default_value=-1.0e10) return d + def to_geo_dataframe(self, gdf=None, kper=0, sparse=False, **kwargs): + """ + Method to add data to a GeoDataFrame for exporting as a geospatial file + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + kper : int + stress period to export + sparse : bool + boolean flag for sparse dataframe construction. Default is False + + Returns + ------- + GeoDataFrame + """ + if self.model is None: + return gdf + else: + modelgrid = self.model.modelgrid + if modelgrid is None: + return gdf + + if gdf is None: + gdf = modelgrid.geo_dataframe + + data = self.array + + active = [] + for name, array4d in data.items(): + aname = f"{self.name[0].lower()}_{name}" + array = array4d[kper] + if modelgrid.grid_type == "unstructured": + array = array.ravel() + if sparse: + idx = np.where(~np.isnan(array))[0] + active.extend(idx) + gdf[aname] = array + else: + for lay in range(modelgrid.nlay): + arr = array[lay].ravel() + if sparse: + idx = np.where(~np.isnan(arr))[0] + active.extend(idx) + gdf[f"{aname}_{lay}"] = arr.ravel() + + if sparse: + gdf = gdf.iloc[active] + + return gdf + def export(self, f, **kwargs): from .. import export From 1c9e53a0b38038f1226b2791c46b325758d6475a Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Wed, 6 Aug 2025 13:21:00 -0700 Subject: [PATCH 02/13] Updates for geodataframe generation and exporting --- flopy/mbase.py | 34 ++++++++++++++++++ flopy/modflow/mfhfb.py | 75 +++++++++++++++++++++++++++++++++++++++ flopy/modflow/mfsfr2.py | 38 ++++++++++++++++++++ flopy/pakbase.py | 6 ++-- flopy/utils/util_array.py | 58 ++++++++++++++++++++++++++++-- flopy/utils/util_list.py | 13 +++---- 6 files changed, 212 insertions(+), 12 deletions(-) diff --git a/flopy/mbase.py b/flopy/mbase.py index 40cdbf700..15c25b09c 100644 --- a/flopy/mbase.py +++ b/flopy/mbase.py @@ -759,6 +759,40 @@ def _output_msg(self, i, add=True): f"{txt2} the output list." ) + def to_geo_dataframe(self, gdf=None, kper=0): + """ + Method to build a Geodataframe from model inputs. Note: transient data + will only be exported for a sinlge stress period. + + Parameters + ---------- + gdf : GeoDataFrame + optional geopandas geodataframe object to add data to. Default is None + kper : int + stress period to get transient data from + + Returns + ------- + gdf : GeoDataFrame + """ + if gdf is None: + modelgrid = self.modelgrid + if modelgrid is not None: + gdf = modelgrid.geo_dataframe + else: + raise AttributeError( + "model does not have a grid instance, " + "please supply a geodataframe" + ) + + for package in self.packagelist: + if package.package_type in ("hfb6",): + continue + if callable(getattr(package, "to_geo_dataframe", None)): + gdf = package.to_geo_dataframe(gdf, kper=kper, sparse=False) + + return gdf + def add_output_file( self, unit, diff --git a/flopy/modflow/mfhfb.py b/flopy/modflow/mfhfb.py index e018f54d4..ad0696a0f 100644 --- a/flopy/modflow/mfhfb.py +++ b/flopy/modflow/mfhfb.py @@ -179,6 +179,81 @@ def _ncells(self): """ return self.nhfbnp + def _get_hfb_lines(self): + """ + Method to get hfb lines. Lines are ordered by recarray index + + Returns + ------- + list : list of hfb lines + """ + modelgrid = self.parent.modelgrid + hfbs = self.hfb_data + iverts = modelgrid.iverts + verts = modelgrid.verts + nodes = [] + if "k" in self.hfb_data.dtype.names: + for rec in hfbs: + node1 = modelgrid.get_node([(0, rec["irow1"], rec["icol1"])])[0] + node2 = modelgrid.get_node([(0, rec["irow2"], rec["icol2"])])[0] + nodes.append((node1, node2)) + else: + nodes = list(zip(hfbs["node1"], hfbs["node2"])) + + shared_edges = [] + for (node0, node1) in nodes: + iv0 = iverts[node0] + iv1 = iverts[node1] + edges = [] + for ix in range(len(iv0)): + edges.append(tuple(sorted((iv0[ix - 1], iv0[ix])))) + + for ix in range(len(iv1)): + edge = tuple(sorted((iv1[ix - 1], iv1[ix]))) + if edge in edges: + shared_edges.append(edge) + break + + if not shared_edges: + raise AssertionError( + f"No shared cell edges found. Cannot represent HFB for nodes {node0} and {node1}" + ) + + lines = [] + for edge in shared_edges: + lines.append((tuple(verts[edge[0]]), tuple(verts[edge[1]]))) + + return lines + + def to_geo_dataframe(self, **kwargs): + """ + Method to create a LineString based geodataframe of horizontal flow barriers + + Returns + ------- + GeoDataFrame + + """ + import geopandas as gpd + + lines = self._get_hfb_lines() + geo_interface = {"type": "FeatureCollection"} + features = [ + { + "id": f"{ix}", + "geometry": {"coordinates": line, "type": "LineString"}, + "properties": {} + } + for ix, line in enumerate(lines) + ] + geo_interface["features"] = features + gdf = gpd.GeoDataFrame.from_features(geo_interface) + hfbs = self.hfb_data + for name in hfbs.dtype.names: + gdf[name] = hfbs[name] + + return gdf + def write_file(self): """ Write the package file. diff --git a/flopy/modflow/mfsfr2.py b/flopy/modflow/mfsfr2.py index 1abd2ce50..c994fe13a 100644 --- a/flopy/modflow/mfsfr2.py +++ b/flopy/modflow/mfsfr2.py @@ -623,6 +623,44 @@ def paths(self): def df(self): return pd.DataFrame(self.reach_data) + def to_geo_dataframe(self, gdf=None, sparse=True, **kwargs): + """ + Method to export SFR reach data to a GeoDataFrame + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + sparse : bool + boolean flag for sparse dataframe construction. Default is True + """ + modelgrid = self.parent.modelgrid + if modelgrid is None: + return gdf + + if gdf is None: + gdf = modelgrid.geo_dataframe + + df = self.df + if "k" in list(df): + df["node"] = df["node"] - (modelgrid.ncpl * df["k"]) + df = df.drop(columns=["k", "i", "j"]) + + df["node"] += 1 + gdf = gdf.merge(df, how="left", on="node") + + + if sparse: + col_names = [col for col in list(df) if col != "node"] + gdf = gdf.dropna(subset=col_names, how="all") + gdf = gdf.dropna(axis="columns", how="all") + else: + gdf = gdf.drop_duplicates(subset=["node"]) + gdf = gdf.reset_index(drop=True) + + return gdf + def _make_graph(self): # get all segments and their outseg graph = {} diff --git a/flopy/pakbase.py b/flopy/pakbase.py index dca6e8e89..3793f3fdd 100644 --- a/flopy/pakbase.py +++ b/flopy/pakbase.py @@ -673,7 +673,7 @@ def export(self, f, **kwargs): return export.utils.package_export(f, self, **kwargs) - def to_geo_dataframe(self, gdf=None, kper=0): + def to_geo_dataframe(self, gdf=None, kper=0, **kwargs): """ Method to create a GeoDataFrame from a modflow package @@ -707,7 +707,9 @@ def to_geo_dataframe(self, gdf=None, kper=0): for attr, value in self.__dict__.items(): if callable(getattr(value, "to_geo_dataframe", None)): - gdf = value.to_geo_dataframe(gdf, forgive=True, kper=kper) + if isinstance(value, (BaseModel, PackageInterface)): + continue + gdf = value.to_geo_dataframe(gdf, forgive=True, kper=kper, sparse=False) return gdf diff --git a/flopy/utils/util_array.py b/flopy/utils/util_array.py index 01f336e39..498092eb0 100644 --- a/flopy/utils/util_array.py +++ b/flopy/utils/util_array.py @@ -633,7 +633,7 @@ def data_type(self): def plottable(self): return True - def to_geo_dataframe(self, gdf=None, forgive=False, **kwargs): + def to_geo_dataframe(self, gdf=None, forgive=False, name=None, **kwargs): """ Method to add data to a GeoDataFrame for exporting as a geospatial file @@ -645,6 +645,8 @@ def to_geo_dataframe(self, gdf=None, forgive=False, **kwargs): forgive : bool optional flag to continue running and pass data that is not compatible with the geodataframe shape + name : str + optional array name base. If None, method uses the .name attribute Returns ------- @@ -660,7 +662,10 @@ def to_geo_dataframe(self, gdf=None, forgive=False, **kwargs): if gdf is None: gdf = modelgrid.geo_dataframe - names = self.name + if name is not None: + names = [name for _ in range(len(self.util_2ds))] + else: + names = self.name for lay, u2d in enumerate(self.util_2ds): name = f"{names[lay]}_{lay}" gdf = u2d.to_geo_dataframe(gdf=gdf, name=name, forgive=forgive) @@ -1092,6 +1097,32 @@ def data_type(self): def plottable(self): return False + def to_geo_dataframe(self, gdf=None, kper=0, forgive=False, **kwargs): + """ + Method to add data to a GeoDataFrame for exporting as a geospatial file + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + kper : int + stress period to export + forgive : bool + boolean flag for sparse dataframe construction. Default is False + + Returns + ------- + GeoDataFrame + """ + u3d = self.transient_3ds[kper] + # note: may need to provide a pass through name for u3d to avoid s.p. + # number being tacked on. Test this on a model with the LAK package... + name = self.name_base[:-1].lower() + gdf = u3d.to_geo_dataframe(gdf=gdf, forgive=forgive, name=name, **kwargs) + return gdf + + def get_zero_3d(self, kper): name = f"{self.name_base}{kper + 1}(filled zero)" return Util3d( @@ -1432,6 +1463,29 @@ def from_4d(cls, model, pak_name, m4ds): name=name, ) + def to_geo_dataframe(self, gdf=None, kper=0, forgive=False, **kwargs): + """ + Method to add data to a GeoDataFrame for exporting as a geospatial file + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + kper : int + stress period to export + forgive : bool + boolean flag for sparse dataframe construction. Default is False + + Returns + ------- + GeoDataFrame + """ + u2d = self.transient_2ds[kper] + name = self.name_base[:-1] + gdf = u2d.to_geo_dataframe(gdf=gdf, name=name, forgive=forgive, **kwargs) + return gdf + def __setattr__(self, key, value): if hasattr(self, "transient_2ds") and key == "cnstnt": # set cnstnt for each u2d diff --git a/flopy/utils/util_list.py b/flopy/utils/util_list.py index edbf71edb..3b99bc86f 100644 --- a/flopy/utils/util_list.py +++ b/flopy/utils/util_list.py @@ -153,26 +153,23 @@ def to_geo_dataframe(self, gdf=None, kper=0, sparse=False, **kwargs): data = self.array - active = [] + col_names = [] for name, array4d in data.items(): aname = f"{self.name[0].lower()}_{name}" array = array4d[kper] if modelgrid.grid_type == "unstructured": array = array.ravel() - if sparse: - idx = np.where(~np.isnan(array))[0] - active.extend(idx) gdf[aname] = array + col_names.append(aname) else: for lay in range(modelgrid.nlay): arr = array[lay].ravel() - if sparse: - idx = np.where(~np.isnan(arr))[0] - active.extend(idx) gdf[f"{aname}_{lay}"] = arr.ravel() + col_names.append(f"{aname}_{lay}") if sparse: - gdf = gdf.iloc[active] + gdf = gdf.dropna(subset=col_names, how="all") + gdf = gdf.dropna(axis="columns", how="all") return gdf From 65a9d16ab3ce3aef9e8b4b93eb351bf0c8a314e6 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Wed, 6 Aug 2025 14:19:43 -0700 Subject: [PATCH 03/13] Add to_geo_dataframe() method to LayerFile for Binary and Formatted layer file support --- flopy/utils/datafile.py | 58 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 57 insertions(+), 1 deletion(-) diff --git a/flopy/utils/datafile.py b/flopy/utils/datafile.py index 183fb874e..369cc0844 100644 --- a/flopy/utils/datafile.py +++ b/flopy/utils/datafile.py @@ -240,6 +240,59 @@ def __enter__(self): def __exit__(self, *exc): self.close() + def to_geo_dataframe( + self, + gdf=None, + modelgrid=None, + kstpkper=None, + totim=None, + attrib_name=None + ): + """ + Generate a GeoDataFrame with data from a LayerFile instance + + Parameters + ---------- + gdf : GeoDataFrame + optional, existing geodataframe with NCPL geometries + modelgrid : Grid + optional modelgrid instance to generate a GeoDataFrame from + kstpkper : tuple of ints + A tuple containing the time step and stress period (kstp, kper). + These are zero-based kstp and kper values. + totim : float + The simulation time. + attrib_name : str + optional base name of attribute columns. (default is text attribute) + + + Returns + ------- + GeoDataFrame + """ + if gdf is None: + if modelgrid is None: + if self.mg is None: + raise AssertionError( + "A geodataframe or modelgrid instance must be supplied" + ) + modelgrid = self.mg + + gdf = modelgrid.geo_dataframe + + array = np.atleast_3d( + self.get_data(kstpkper=kstpkper, totim=totim).transpose() + ).transpose() + + if attrib_name is None: + attrib_name = self.text.decode() + + for ix, arr in enumerate(array): + name = f"{attrib_name}_{ix}" + gdf[name] = np.ravel(arr) + + return gdf + def to_shapefile( self, filename: Union[str, os.PathLike], @@ -287,7 +340,10 @@ def to_shapefile( >>> times = hdobj.get_times() >>> hdobj.to_shapefile('test_heads_sp6.shp', totim=times[-1]) """ - + warnings.warn( + "to_shapefile() is deprecated and is being replaced by to_geo_dataframe()", + DeprecationWarning + ) plotarray = np.atleast_3d( self.get_data(kstpkper=kstpkper, totim=totim, mflay=mflay).transpose() ).transpose() From 23129fe7963ffd1f10569c56c314a0b78417d5fb Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Wed, 6 Aug 2025 15:07:24 -0700 Subject: [PATCH 04/13] Add geodataframe exporting to CellBudgetFile --- flopy/utils/binaryfile/__init__.py | 51 ++++++++++++++++++++++++++++++ flopy/utils/datafile.py | 20 ++++++------ 2 files changed, 61 insertions(+), 10 deletions(-) diff --git a/flopy/utils/binaryfile/__init__.py b/flopy/utils/binaryfile/__init__.py index 5c9842b74..3000b96b8 100644 --- a/flopy/utils/binaryfile/__init__.py +++ b/flopy/utils/binaryfile/__init__.py @@ -1604,6 +1604,57 @@ def get_data( if t ] + def to_geo_dataframe( + self, + gdf=None, + modelgrid=None, + idx=None, + kstpkper=None, + totim=None, + text=None, + ): + if (idx, kstpkper, totim) == (None, None, None): + raise AssertionError( + "to_geo_dataframe() missing 1 required argument: " + "please provide 'idx', 'kstpkper', or 'totim'" + ) + + + if gdf is None: + if modelgrid is None: + if self.modelgrid is None: + raise AssertionError( + "A geodataframe or modelgrid instance must be supplied" + ) + modelgrid = self.modelgrid + + gdf = modelgrid.geo_dataframe + + col_names = [] + if text is None: + textlist = [i.decode() for i in self.textlist] + for text in textlist: + data = self.get_data( + idx=idx, kstpkper=kstpkper, totim=totim, text=text, full3D=True + ) + + text = text.strip().lower().replace(" ", "_") + for ix, arr in enumerate(data[0]): + name = f"{text}_{ix}" + gdf[name] = arr.ravel() + col_names.append(name) + else: + data = self.get_data( + idx=idx, kstpkper=kstpkper, totim=totim, text=text, full3D=True + ) + text = text.strip().lower().replace(" ", "_") + for ix, arr in enumerate(data[0]): + name = f"{text}_{ix}" + gdf[name] = arr.ravel() + col_names.append(name) + + return gdf + def get_ts(self, idx, text=None, times=None): """ Get a time series from the binary budget file. diff --git a/flopy/utils/datafile.py b/flopy/utils/datafile.py index 369cc0844..e6fd901c4 100644 --- a/flopy/utils/datafile.py +++ b/flopy/utils/datafile.py @@ -197,18 +197,18 @@ def __init__(self, filename: Union[str, os.PathLike], precision, verbose, **kwar self.model = None self.dis = None - self.mg = None + self.modelgrid = None if "model" in kwargs.keys(): self.model = kwargs.pop("model") - self.mg = self.model.modelgrid + self.modelgrid = self.model.modelgrid self.dis = self.model.dis if "dis" in kwargs.keys(): self.dis = kwargs.pop("dis") - self.mg = self.dis.parent.modelgrid + self.modelgrid = self.dis.parent.modelgrid if "tdis" in kwargs.keys(): self.tdis = kwargs.pop("tdis") if "modelgrid" in kwargs.keys(): - self.mg = kwargs.pop("modelgrid") + self.modelgrid = kwargs.pop("modelgrid") if len(kwargs.keys()) > 0: args = ",".join(kwargs.keys()) raise Exception(f"LayerFile error: unrecognized kwargs: {args}") @@ -218,8 +218,8 @@ def __init__(self, filename: Union[str, os.PathLike], precision, verbose, **kwar # now that we read the data and know nrow and ncol, # we can make a generic mg if needed - if self.mg is None: - self.mg = StructuredGrid( + if self.modelgrid is None: + self.modelgrid = StructuredGrid( delc=np.ones((self.nrow,)), delr=np.ones(self.ncol), nlay=self.nlay, @@ -272,11 +272,11 @@ def to_geo_dataframe( """ if gdf is None: if modelgrid is None: - if self.mg is None: + if self.modelgrid is None: raise AssertionError( "A geodataframe or modelgrid instance must be supplied" ) - modelgrid = self.mg + modelgrid = self.modelgrid gdf = modelgrid.geo_dataframe @@ -357,7 +357,7 @@ def to_shapefile( from ..export.shapefile_utils import write_grid_shapefile - write_grid_shapefile(filename, self.mg, attrib_dict, verbose) + write_grid_shapefile(filename, self.modelgrid, attrib_dict, verbose) def plot( self, @@ -469,7 +469,7 @@ def plot( axes=axes, filenames=filenames, mflay=mflay, - modelgrid=self.mg, + modelgrid=self.modelgrid, **kwargs, ) From eec18e1a9231287469a231d17c6717094ba3d262 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Thu, 7 Aug 2025 12:50:12 -0700 Subject: [PATCH 05/13] Begin adding MF6 geodataframe support --- flopy/mf6/data/mfdataarray.py | 25 ++++++--- flopy/mf6/data/mfdatalist.py | 103 ++++++++++++++++++++++------------ flopy/mf6/mfpackage.py | 39 +++++++++++++ 3 files changed, 121 insertions(+), 46 deletions(-) diff --git a/flopy/mf6/data/mfdataarray.py b/flopy/mf6/data/mfdataarray.py index f4c2b6f09..37501a2b6 100644 --- a/flopy/mf6/data/mfdataarray.py +++ b/flopy/mf6/data/mfdataarray.py @@ -323,19 +323,22 @@ def data(self): """Returns array data. Calls get_data with default parameters.""" return self._get_data() - def to_geo_dataframe(self, gdf=None): + def to_geo_dataframe(self, gdf=None, name=None, forgive=False, **kwargs): """ - Method to add data to a GeoDataFrame for exporting as a geospatial file + Method to add an input array to a geopandas GeoDataFrame Parameters ---------- gdf : GeoDataFrame - optional GeoDataFrame instance. If GeoDataFrame is None, one will be - constructed from modelgrid information + optional GeoDataFrame object + name : str + optional attribute name, default uses util2d name + forgive : bool + optional flag to continue if data shape not compatible with GeoDataFrame Returns ------- - GeoDataFrame + geopandas GeoDataFrame """ if self.model is None: return gdf @@ -347,23 +350,27 @@ def to_geo_dataframe(self, gdf=None): gdf = modelgrid.geo_dataframe if modelgrid is not None: - if modelgrid.grid_type() != "unstructured": + if modelgrid.grid_type != "unstructured": ncpl = modelgrid.ncpl else: ncpl = modelgrid.nnodes else: ncpl = len(gdf) - name = self.name + if name is None: + name = self.name + data = self.array if data.size == ncpl: gdf[name] = data.ravel() elif data.size % ncpl == 0: data = data.reshape((-1, ncpl)) - for ix, dat in enumerate(data): + for ix, arr in enumerate(data): aname = f"{name}_{ix}" - gdf[aname] = dat + gdf[aname] = arr + elif forgive: + return gdf else: raise ValueError( f"Data size {data.size} not compatible with dataframe length {ncpl}" diff --git a/flopy/mf6/data/mfdatalist.py b/flopy/mf6/data/mfdatalist.py index 22e5490c1..ae7226e67 100644 --- a/flopy/mf6/data/mfdatalist.py +++ b/flopy/mf6/data/mfdatalist.py @@ -145,7 +145,7 @@ def to_array(self, kper=0, mask=False): model_grid = self.data_dimensions.get_model_grid() return list_to_array(sarr, model_grid, kper, mask) - def to_geo_dataframe(self, gdf=None, sparse=False): + def to_geo_dataframe(self, gdf=None, sparse=False, **kwargs): """ Method to add data to a GeoDataFrame for exporting as a geospatial file @@ -173,47 +173,26 @@ def to_geo_dataframe(self, gdf=None, sparse=False): if gdf is None: gdf = modelgrid.geo_dataframe - # name = self.name - recarray = self.array - if "cellid" not in recarray.dtype.names: + data = self.to_array(mask=True) + if data is None: return gdf - cellid = recarray["cellid"] - df = pd.DataFrame.from_records(recarray) - - grid_type = modelgrid.grid_type() - if grid_type != "unstructured": - layer = [cid[0] for cid in cellid] - if grid_type == "structured": - cids = [(0,) + cid[1:] for cid in cellid] - nodes = modelgrid.get_node(cids) + col_names = [] + for name, array3d in data.items(): + aname = f"{self.path[1].lower()}_{name}" + if modelgrid.grid_type == "unstructured": + array = array3d.ravel() + gdf[aname] = array + col_names.append(aname) else: - nodes = [cid[1] for cid in cellid] - - df["layer"] = layer - df["node"] = nodes - df = df.drop(columns=["cellid", ]) - for layer in df.layer.unique(): - tmp = df[df.layer == layer] - tmp = tmp.drop(columns=["layer"]) - renames = { - nm: f"{nm}_{layer}" for nm in list(df) if nm not in ("node",) - } - tmp = tmp.set_index("node") - tmp = tmp.rename(columns=renames) - gdf = pd.merge(gdf, tmp, how="left", left_index=True, right_index=True) - - else: - nodes = [cid[0] for cid in cellid] - df["node"] = nodes - df = df.drop(columns=["cellid",]) - df = df.set_index("node") - gdf = pd.merge(gdf, df, how="left", left_index=True, right_index=True) + for lay in range(modelgrid.nlay): + arr = array3d[lay].ravel() + gdf[f"{aname}_{lay}"] = arr.ravel() + col_names.append(f"{aname}_{lay}") if sparse: - cols = [col for col in list(df) if col != "node"] - if cols: - gdf = gdf[~np.isnan(gdf[cols[0]])] + gdf = gdf.dropna(subset=col_names, how="all") + gdf = gdf.dropna(axis="columns", how="all") return gdf @@ -1668,6 +1647,56 @@ def to_array(self, kper=0, mask=False): """Returns list data as an array.""" return super().to_array(kper, mask) + def to_geo_dataframe(self, gdf=None, kper=0, sparse=False, **kwargs): + """ + Method to add data to a GeoDataFrame for exporting as a geospatial file + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + kper : int + stress period to export + sparse : bool + boolean flag for sparse dataframe construction. Default is False + + Returns + ------- + GeoDataFrame + """ + if self.model is None: + return gdf + else: + modelgrid = self.model.modelgrid + if modelgrid is None: + return gdf + + if gdf is None: + gdf = modelgrid.geo_dataframe + + data = self.masked_4D_arrays + + col_names = [] + for name, array4d in data.items(): + aname = f"{self.path[1].lower()}_{name}" + array = array4d[kper] + if modelgrid.grid_type == "unstructured": + array = array.ravel() + gdf[aname] = array + col_names.append(aname) + else: + for lay in range(modelgrid.nlay): + arr = array[lay].ravel() + gdf[f"{aname}_{lay}"] = arr.ravel() + col_names.append(f"{aname}_{lay}") + + if sparse: + gdf = gdf.dropna(subset=col_names, how="all") + gdf = gdf.dropna(axis="columns", how="all") + + return gdf + def remove_transient_key(self, transient_key): """Remove transient stress period key. Method is used internally by FloPy and is not intended to the end user. diff --git a/flopy/mf6/mfpackage.py b/flopy/mf6/mfpackage.py index 4ee9d6175..3854c8112 100644 --- a/flopy/mf6/mfpackage.py +++ b/flopy/mf6/mfpackage.py @@ -2086,6 +2086,45 @@ def package_filename_dict(self): ) return self._package_container.package_filename_dict + def to_geo_dataframe(self, gdf=None, kper=0, **kwargs): + """ + Method to create a GeoDataFrame from a modflow package + + Parameters + ---------- + gdf : GeoDataFrame + optional geopandas geodataframe object to add data to. Default is None + kper : int + stress period to get transient data from + + Returns + ------- + gdf : GeoDataFrame + """ + if gdf is None: + if isinstance(self.parent, ModelInterface): + modelgrid = self.parent.modelgrid + if modelgrid is not None: + gdf = modelgrid.geo_dataframe + else: + raise AttributeError( + "model does not have a grid instance, " + "please supply a geodataframe" + ) + else: + raise AssertionError( + "Package does not have a model instance, " + "please supply a geodataframe" + ) + + for attr, value in self.__dict__.items(): + if callable(getattr(value, "to_geo_dataframe", None)): + if isinstance(value, (ModelInterface, PackageInterface)): + continue + gdf = value.to_geo_dataframe(gdf, forgive=True, kper=kper, sparse=False) + + return gdf + def get_package(self, name=None, type_only=False, name_only=False): """ Finds a package by package name, package key, package type, or partial From f117c6410e3b5dd71325507d28c9b32d44e6fd89 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Thu, 7 Aug 2025 14:15:56 -0700 Subject: [PATCH 06/13] Continue linking up code for geodataframe construction from model inputs --- flopy/mf6/data/mfdataarray.py | 47 ++++++++++++++++ flopy/mf6/data/mfdatalist.py | 13 ++--- flopy/mf6/data/mfdataplist.py | 102 +++++++++++++++++++++++++++++++++- flopy/mf6/data/mfdatautil.py | 2 + flopy/mf6/mfmodel.py | 34 ++++++++++++ flopy/mf6/mfpackage.py | 19 +++++-- flopy/pakbase.py | 9 ++- 7 files changed, 209 insertions(+), 17 deletions(-) diff --git a/flopy/mf6/data/mfdataarray.py b/flopy/mf6/data/mfdataarray.py index 37501a2b6..ef3e78508 100644 --- a/flopy/mf6/data/mfdataarray.py +++ b/flopy/mf6/data/mfdataarray.py @@ -361,6 +361,9 @@ def to_geo_dataframe(self, gdf=None, name=None, forgive=False, **kwargs): name = self.name data = self.array + if data is None: + return gdf + if data.size == ncpl: gdf[name] = data.ravel() @@ -1945,6 +1948,50 @@ def _build_period_data( output[sp] = data return output + def to_geo_dataframe(self, gdf=None, kper=0, forgive=False, **kwargs): + """ + + """ + if self.model is None: + return gdf + else: + modelgrid = self.model.modelgrid + if gdf is None: + if modelgrid is None: + return gdf + gdf = modelgrid.geo_dataframe + + if modelgrid is not None: + if modelgrid.grid_type != "unstructured": + ncpl = modelgrid.ncpl + else: + ncpl = modelgrid.nnodes + else: + ncpl = len(gdf) + + if self.array is None: + return gdf + + name = f"{self.path[1]}_{self.name}" + + data = self.get_data(key=kper, apply_mult=True) + if data.size == ncpl: + gdf[name] = data.ravel() + + elif data.size % ncpl == 0: + data = data.reshape((-1, ncpl)) + for ix, arr in enumerate(data): + aname = f"{name}_{ix}" + gdf[aname] = arr + elif forgive: + return gdf + else: + raise ValueError( + f"Data size {data.size} not compatible with dataframe length {ncpl}" + ) + + return gdf + def set_record(self, data_record): """Sets data and metadata at layer `layer` and time `key` to `data_record`. For unlayered data do not pass in `layer`. diff --git a/flopy/mf6/data/mfdatalist.py b/flopy/mf6/data/mfdatalist.py index ae7226e67..3d05e79a2 100644 --- a/flopy/mf6/data/mfdatalist.py +++ b/flopy/mf6/data/mfdatalist.py @@ -161,8 +161,6 @@ def to_geo_dataframe(self, gdf=None, sparse=False, **kwargs): ------- GeoDataFrame """ - import pandas as pd - if self.model is None: return gdf else: @@ -1675,19 +1673,20 @@ def to_geo_dataframe(self, gdf=None, kper=0, sparse=False, **kwargs): if gdf is None: gdf = modelgrid.geo_dataframe - data = self.masked_4D_arrays + data = self.to_array(kper=kper, mask=True) + if data is None: + return gdf col_names = [] - for name, array4d in data.items(): + for name, array3d in data.items(): aname = f"{self.path[1].lower()}_{name}" - array = array4d[kper] if modelgrid.grid_type == "unstructured": - array = array.ravel() + array = array3d.ravel() gdf[aname] = array col_names.append(aname) else: for lay in range(modelgrid.nlay): - arr = array[lay].ravel() + arr = array3d[lay].ravel() gdf[f"{aname}_{lay}"] = arr.ravel() col_names.append(f"{aname}_{lay}") diff --git a/flopy/mf6/data/mfdataplist.py b/flopy/mf6/data/mfdataplist.py index 331177ef0..60c9162cf 100644 --- a/flopy/mf6/data/mfdataplist.py +++ b/flopy/mf6/data/mfdataplist.py @@ -839,6 +839,55 @@ def to_array(self, kper=0, mask=False): model_grid = self.data_dimensions.get_model_grid() return list_to_array(sarr, model_grid, kper, mask) + def to_geo_dataframe(self, gdf=None, sparse=False, **kwargs): + """ + Method to add data to a GeoDataFrame for exporting as a geospatial file + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + sparse : bool + boolean flag for sparse dataframe construction. Default is False + + Returns + ------- + GeoDataFrame + """ + if self.model is None: + return gdf + else: + modelgrid = self.model.modelgrid + if modelgrid is None: + return gdf + + if gdf is None: + gdf = modelgrid.geo_dataframe + + data = self.to_array(mask=True) + if data is None: + return gdf + + col_names = [] + for name, array3d in data.items(): + aname = f"{self.path[1].lower()}_{name}" + if modelgrid.grid_type == "unstructured": + array = array3d.ravel() + gdf[aname] = array + col_names.append(aname) + else: + for lay in range(modelgrid.nlay): + arr = array3d[lay].ravel() + gdf[f"{aname}_{lay}"] = arr.ravel() + col_names.append(f"{aname}_{lay}") + + if sparse: + gdf = gdf.dropna(subset=col_names, how="all") + gdf = gdf.dropna(axis="columns", how="all") + + return gdf + def set_record(self, record, autofill=False, check_data=True): """Sets the contents of the data and metadata to "data_record". Data_record is a dictionary with has the following format: @@ -1935,9 +1984,7 @@ def plot( ) -class MFPandasTransientList( - MFPandasList, mfdata.MFTransient, DataListInterface -): +class MFPandasTransientList(MFPandasList, mfdata.MFTransient, DataListInterface): """ Provides an interface for the user to access and update MODFLOW transient pandas list data. @@ -2019,6 +2066,55 @@ def to_array(self, kper=0, mask=False): """Returns list data as an array.""" return super().to_array(kper, mask) + def to_geo_dataframe(self, gdf=None, kper=0, sparse=False, **kwargs): + """ + Method to add data to a GeoDataFrame for exporting as a geospatial file + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + kper : int + stress period to export + sparse : bool + boolean flag for sparse dataframe construction. Default is False + + Returns + ------- + GeoDataFrame + """ + if self.model is None: + return gdf + else: + modelgrid = self.model.modelgrid + if modelgrid is None: + return gdf + + if gdf is None: + gdf = modelgrid.geo_dataframe + + data = self.to_array(kper=kper, mask=True) + + col_names = [] + for name, array3d in data.items(): + aname = f"{self.path[1].lower()}_{name}" + if modelgrid.grid_type == "unstructured": + array = array3d.ravel() + gdf[aname] = array + col_names.append(aname) + else: + for lay in range(modelgrid.nlay): + arr = array3d[lay].ravel() + gdf[f"{aname}_{lay}"] = arr.ravel() + col_names.append(f"{aname}_{lay}") + + if sparse: + gdf = gdf.dropna(subset=col_names, how="all") + gdf = gdf.dropna(axis="columns", how="all") + + return gdf + def remove_transient_key(self, transient_key): """Remove transient stress period key. Method is used internally by FloPy and is not intended to the end user. diff --git a/flopy/mf6/data/mfdatautil.py b/flopy/mf6/data/mfdatautil.py index 3260ca12c..bcbd600ea 100644 --- a/flopy/mf6/data/mfdatautil.py +++ b/flopy/mf6/data/mfdatautil.py @@ -202,6 +202,8 @@ def list_to_array(sarr, model_grid, kper=0, mask=False): cnt = np.zeros(shape, dtype=np.float64) for sp_rec in sarr: if sp_rec is not None: + if "cellid" not in sp_rec.dtype.names: + return None for rec in sp_rec: arr[rec["cellid"]] += rec[name] cnt[rec["cellid"]] += 1.0 diff --git a/flopy/mf6/mfmodel.py b/flopy/mf6/mfmodel.py index 0abc25cbd..907a33454 100644 --- a/flopy/mf6/mfmodel.py +++ b/flopy/mf6/mfmodel.py @@ -789,6 +789,40 @@ def output(self): except AttributeError: return MF6Output(self, budgetkey=budgetkey) + def to_geo_dataframe(self, gdf=None, kper=0): + """ + Method to build a Geodataframe from model inputs. Note: transient data + will only be exported for a sinlge stress period. + + Parameters + ---------- + gdf : GeoDataFrame + optional geopandas geodataframe object to add data to. Default is None + kper : int + stress period to get transient data from + + Returns + ------- + gdf : GeoDataFrame + """ + if gdf is None: + modelgrid = self.modelgrid + if modelgrid is not None: + gdf = modelgrid.geo_dataframe + else: + raise AttributeError( + "model does not have a grid instance, " + "please supply a geodataframe" + ) + + for package in self.packagelist: + if package.package_type in ("hfb",): + continue + if callable(getattr(package, "to_geo_dataframe", None)): + gdf = package.to_geo_dataframe(gdf, kper=kper, sparse=False) + + return gdf + def export(self, f, **kwargs): """Method to export a model to a shapefile or netcdf file diff --git a/flopy/mf6/mfpackage.py b/flopy/mf6/mfpackage.py index 3854c8112..d7f289143 100644 --- a/flopy/mf6/mfpackage.py +++ b/flopy/mf6/mfpackage.py @@ -2086,7 +2086,7 @@ def package_filename_dict(self): ) return self._package_container.package_filename_dict - def to_geo_dataframe(self, gdf=None, kper=0, **kwargs): + def to_geo_dataframe(self, gdf=None, kper=0, sparse=False, **kwargs): """ Method to create a GeoDataFrame from a modflow package @@ -2117,11 +2117,18 @@ def to_geo_dataframe(self, gdf=None, kper=0, **kwargs): "please supply a geodataframe" ) - for attr, value in self.__dict__.items(): - if callable(getattr(value, "to_geo_dataframe", None)): - if isinstance(value, (ModelInterface, PackageInterface)): - continue - gdf = value.to_geo_dataframe(gdf, forgive=True, kper=kper, sparse=False) + for attr, value in self.__dict__.items(): + if callable(getattr(value, "to_geo_dataframe", None)): + if isinstance(value, (ModelInterface, PackageInterface)): + continue + # do not pass sparse in here, "sparsify" after all data has been + # added to geodataframe + gdf = value.to_geo_dataframe(gdf, forgive=True, kper=kper, sparse=False) + + if sparse: + col_names = [i for i in gdf if i not in ("geometry", "node", "row", "col")] + gdf = gdf.dropna(subset=col_names, how="all") + gdf = gdf.dropna(axis="columns", how="all") return gdf diff --git a/flopy/pakbase.py b/flopy/pakbase.py index 3793f3fdd..b0a55b48b 100644 --- a/flopy/pakbase.py +++ b/flopy/pakbase.py @@ -673,7 +673,7 @@ def export(self, f, **kwargs): return export.utils.package_export(f, self, **kwargs) - def to_geo_dataframe(self, gdf=None, kper=0, **kwargs): + def to_geo_dataframe(self, gdf=None, kper=0, sparse=False, **kwargs): """ Method to create a GeoDataFrame from a modflow package @@ -709,8 +709,15 @@ def to_geo_dataframe(self, gdf=None, kper=0, **kwargs): if callable(getattr(value, "to_geo_dataframe", None)): if isinstance(value, (BaseModel, PackageInterface)): continue + # do not pass sparse in here, make sparse after all data has been + # added to geodataframe gdf = value.to_geo_dataframe(gdf, forgive=True, kper=kper, sparse=False) + if sparse: + col_names = [i for i in gdf if i not in ("geometry", "node", "row", "col")] + gdf = gdf.dropna(subset=col_names, how="all") + gdf = gdf.dropna(axis="columns", how="all") + return gdf def _generate_heading(self): From 40b508c38274382cfda81f9e3f552d11ab2874e9 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Thu, 7 Aug 2025 15:57:59 -0700 Subject: [PATCH 07/13] Add in code for hfb linework exporting --- flopy/mf6/mfpackage.py | 26 +++++++++++++++- flopy/modflow/mfhfb.py | 42 ++++---------------------- flopy/plot/plotutil.py | 67 ++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 97 insertions(+), 38 deletions(-) diff --git a/flopy/mf6/mfpackage.py b/flopy/mf6/mfpackage.py index d7f289143..b2dcacd2c 100644 --- a/flopy/mf6/mfpackage.py +++ b/flopy/mf6/mfpackage.py @@ -2105,7 +2105,31 @@ def to_geo_dataframe(self, gdf=None, kper=0, sparse=False, **kwargs): if isinstance(self.parent, ModelInterface): modelgrid = self.parent.modelgrid if modelgrid is not None: - gdf = modelgrid.geo_dataframe + if self.package_type == "hfb": + import geopandas as gpd + from ..plot.plotutil import hfb_data_to_linework + + recarray = self.stress_period_data.data[kper] + lines = hfb_data_to_linework(recarray, modelgrid) + geo_interface = {"type": "FeatureCollection"} + features = [ + { + "id": f"{ix}", + "geometry": {"coordinates": line, "type": "LineString"}, + "properties": {} + } + for ix, line in enumerate(lines) + ] + geo_interface["features"] = features + gdf = gpd.GeoDataFrame.from_features(geo_interface) + + for name in recarray.dtype.names: + gdf[name] = recarray[name] + + return gdf + + else: + gdf = modelgrid.geo_dataframe else: raise AttributeError( "model does not have a grid instance, " diff --git a/flopy/modflow/mfhfb.py b/flopy/modflow/mfhfb.py index ad0696a0f..fc9dc397d 100644 --- a/flopy/modflow/mfhfb.py +++ b/flopy/modflow/mfhfb.py @@ -187,43 +187,8 @@ def _get_hfb_lines(self): ------- list : list of hfb lines """ - modelgrid = self.parent.modelgrid - hfbs = self.hfb_data - iverts = modelgrid.iverts - verts = modelgrid.verts - nodes = [] - if "k" in self.hfb_data.dtype.names: - for rec in hfbs: - node1 = modelgrid.get_node([(0, rec["irow1"], rec["icol1"])])[0] - node2 = modelgrid.get_node([(0, rec["irow2"], rec["icol2"])])[0] - nodes.append((node1, node2)) - else: - nodes = list(zip(hfbs["node1"], hfbs["node2"])) - - shared_edges = [] - for (node0, node1) in nodes: - iv0 = iverts[node0] - iv1 = iverts[node1] - edges = [] - for ix in range(len(iv0)): - edges.append(tuple(sorted((iv0[ix - 1], iv0[ix])))) - - for ix in range(len(iv1)): - edge = tuple(sorted((iv1[ix - 1], iv1[ix]))) - if edge in edges: - shared_edges.append(edge) - break - - if not shared_edges: - raise AssertionError( - f"No shared cell edges found. Cannot represent HFB for nodes {node0} and {node1}" - ) - - lines = [] - for edge in shared_edges: - lines.append((tuple(verts[edge[0]]), tuple(verts[edge[1]]))) - - return lines + from ..plot.plotutil import hfb_data_to_linework + return hfb_data_to_linework(self.hfb_data, self.parent.modelgrid) def to_geo_dataframe(self, **kwargs): """ @@ -289,6 +254,9 @@ def get_empty(ncells=0, aux_names=None, structured=True): been extended to include aux variables and associated aux names. + Returns + ------- + np.recarray """ dtype = ModflowHfb.get_default_dtype(structured=structured) if aux_names is not None: diff --git a/flopy/plot/plotutil.py b/flopy/plot/plotutil.py index 8c1bfec41..6a2ccb6bd 100644 --- a/flopy/plot/plotutil.py +++ b/flopy/plot/plotutil.py @@ -2954,3 +2954,70 @@ def to_prt_pathlines( return df else: return ret + + +def hfb_data_to_linework(recarray, modelgrid): + """ + Method to get lines representing horizontal flow barriers + + Parameters + ---------- + recarray : np.recarray + recarray of hfb input data + modelgrid : Grid + modelgrid instance + + Returns + ------- + list : list of line segments + """ + iverts = modelgrid.iverts + verts = modelgrid.verts + nodes = [] + if modelgrid.grid_type == "structured": + if "k" in recarray.dtype.names: + for rec in recarray: + node1 = modelgrid.get_node([(0, rec["irow1"], rec["icol1"])])[0] + node2 = modelgrid.get_node([(0, rec["irow2"], rec["icol2"])])[0] + nodes.append((node1, node2)) + else: + for rec in recarray: + node1 = modelgrid.get_node([(0,) + rec["cellid1"][1:]])[0] + node2 = modelgrid.get_node([(0,) + rec["cellid2"][1:]])[0] + nodes.append((node1, node2)) + + elif modelgrid.grid_type == "vertex": + for rec in recarray: + nodes.append((rec["cellid1"][-1], rec["cellid2"][-1])) + + else: + if "node1" in recarray.dtype.names: + nodes = list(zip(recarray["node1"], recarray["node2"])) + else: + for rec in recarray: + nodes.append((rec["cellid1"][0], rec["cellid2"][0])) + + shared_edges = [] + for (node0, node1) in nodes: + iv0 = iverts[node0] + iv1 = iverts[node1] + edges = [] + for ix in range(len(iv0)): + edges.append(tuple(sorted((iv0[ix - 1], iv0[ix])))) + + for ix in range(len(iv1)): + edge = tuple(sorted((iv1[ix - 1], iv1[ix]))) + if edge in edges: + shared_edges.append(edge) + break + + if not shared_edges: + raise AssertionError( + f"No shared cell edges found. Cannot represent HFB for nodes {node0} and {node1}" + ) + + lines = [] + for edge in shared_edges: + lines.append((tuple(verts[edge[0]]), tuple(verts[edge[1]]))) + + return lines \ No newline at end of file From b7860e2f91bf8d368387e35428e7ecdf8382192b Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Thu, 7 Aug 2025 16:55:56 -0700 Subject: [PATCH 08/13] fix testing and lint progress --- autotest/test_binaryfile.py | 2 +- autotest/test_formattedfile.py | 2 +- flopy/mbase.py | 3 +-- flopy/mf6/mfpackage.py | 1 + flopy/modflow/mfhfb.py | 3 ++- flopy/modflow/mfsfr2.py | 1 - flopy/pakbase.py | 1 + flopy/plot/plotutil.py | 7 ++++--- flopy/utils/binaryfile/__init__.py | 1 - flopy/utils/datafile.py | 9 ++------- flopy/utils/util_array.py | 1 - 11 files changed, 13 insertions(+), 18 deletions(-) diff --git a/autotest/test_binaryfile.py b/autotest/test_binaryfile.py index 8f4d265a1..bc9ec5583 100644 --- a/autotest/test_binaryfile.py +++ b/autotest/test_binaryfile.py @@ -275,7 +275,7 @@ def test_load_binary_head_file(example_data_path): def test_plot_binary_head_file(example_data_path): hf = HeadFile(example_data_path / "freyberg" / "freyberg.githds") - hf.mg.set_coord_info(xoff=1000.0, yoff=200.0, angrot=15.0) + hf.modelgrid.set_coord_info(xoff=1000.0, yoff=200.0, angrot=15.0) assert isinstance(hf.plot(), Axes) plt.close() diff --git a/autotest/test_formattedfile.py b/autotest/test_formattedfile.py index 0d31ed308..26ac20f45 100644 --- a/autotest/test_formattedfile.py +++ b/autotest/test_formattedfile.py @@ -69,7 +69,7 @@ def test_headfile_build_index(example_data_path): def test_formattedfile_reference(example_data_path): h = FormattedHeadFile(example_data_path / "mf2005_test" / "test1tr.githds") assert isinstance(h, FormattedHeadFile) - h.mg.set_coord_info(xoff=1000.0, yoff=200.0, angrot=15.0) + h.modelgrid.set_coord_info(xoff=1000.0, yoff=200.0, angrot=15.0) assert isinstance(h.plot(masked_values=[6999.000]), Axes) plt.close() diff --git a/flopy/mbase.py b/flopy/mbase.py index 15c25b09c..58d3145f4 100644 --- a/flopy/mbase.py +++ b/flopy/mbase.py @@ -781,8 +781,7 @@ def to_geo_dataframe(self, gdf=None, kper=0): gdf = modelgrid.geo_dataframe else: raise AttributeError( - "model does not have a grid instance, " - "please supply a geodataframe" + "model does not have a grid instance, please supply a geodataframe" ) for package in self.packagelist: diff --git a/flopy/mf6/mfpackage.py b/flopy/mf6/mfpackage.py index b2dcacd2c..c98d687d9 100644 --- a/flopy/mf6/mfpackage.py +++ b/flopy/mf6/mfpackage.py @@ -2107,6 +2107,7 @@ def to_geo_dataframe(self, gdf=None, kper=0, sparse=False, **kwargs): if modelgrid is not None: if self.package_type == "hfb": import geopandas as gpd + from ..plot.plotutil import hfb_data_to_linework recarray = self.stress_period_data.data[kper] diff --git a/flopy/modflow/mfhfb.py b/flopy/modflow/mfhfb.py index fc9dc397d..d695dbedd 100644 --- a/flopy/modflow/mfhfb.py +++ b/flopy/modflow/mfhfb.py @@ -188,6 +188,7 @@ def _get_hfb_lines(self): list : list of hfb lines """ from ..plot.plotutil import hfb_data_to_linework + return hfb_data_to_linework(self.hfb_data, self.parent.modelgrid) def to_geo_dataframe(self, **kwargs): @@ -207,7 +208,7 @@ def to_geo_dataframe(self, **kwargs): { "id": f"{ix}", "geometry": {"coordinates": line, "type": "LineString"}, - "properties": {} + "properties": {}, } for ix, line in enumerate(lines) ] diff --git a/flopy/modflow/mfsfr2.py b/flopy/modflow/mfsfr2.py index c994fe13a..27e791b5e 100644 --- a/flopy/modflow/mfsfr2.py +++ b/flopy/modflow/mfsfr2.py @@ -650,7 +650,6 @@ def to_geo_dataframe(self, gdf=None, sparse=True, **kwargs): df["node"] += 1 gdf = gdf.merge(df, how="left", on="node") - if sparse: col_names = [col for col in list(df) if col != "node"] gdf = gdf.dropna(subset=col_names, how="all") diff --git a/flopy/pakbase.py b/flopy/pakbase.py index b0a55b48b..8718f188d 100644 --- a/flopy/pakbase.py +++ b/flopy/pakbase.py @@ -689,6 +689,7 @@ def to_geo_dataframe(self, gdf=None, kper=0, sparse=False, **kwargs): gdf : GeoDataFrame """ from .mbase import BaseModel + if gdf is None: if isinstance(self.parent, BaseModel): modelgrid = self.parent.modelgrid diff --git a/flopy/plot/plotutil.py b/flopy/plot/plotutil.py index 6a2ccb6bd..17393f6f7 100644 --- a/flopy/plot/plotutil.py +++ b/flopy/plot/plotutil.py @@ -2998,7 +2998,7 @@ def hfb_data_to_linework(recarray, modelgrid): nodes.append((rec["cellid1"][0], rec["cellid2"][0])) shared_edges = [] - for (node0, node1) in nodes: + for node0, node1 in nodes: iv0 = iverts[node0] iv1 = iverts[node1] edges = [] @@ -3013,11 +3013,12 @@ def hfb_data_to_linework(recarray, modelgrid): if not shared_edges: raise AssertionError( - f"No shared cell edges found. Cannot represent HFB for nodes {node0} and {node1}" + f"No shared cell edges found. Cannot represent HFB " + f"for nodes {node0} and {node1}" ) lines = [] for edge in shared_edges: lines.append((tuple(verts[edge[0]]), tuple(verts[edge[1]]))) - return lines \ No newline at end of file + return lines diff --git a/flopy/utils/binaryfile/__init__.py b/flopy/utils/binaryfile/__init__.py index 3000b96b8..50ea1a278 100644 --- a/flopy/utils/binaryfile/__init__.py +++ b/flopy/utils/binaryfile/__init__.py @@ -1619,7 +1619,6 @@ def to_geo_dataframe( "please provide 'idx', 'kstpkper', or 'totim'" ) - if gdf is None: if modelgrid is None: if self.modelgrid is None: diff --git a/flopy/utils/datafile.py b/flopy/utils/datafile.py index e6fd901c4..78da71c2f 100644 --- a/flopy/utils/datafile.py +++ b/flopy/utils/datafile.py @@ -241,12 +241,7 @@ def __exit__(self, *exc): self.close() def to_geo_dataframe( - self, - gdf=None, - modelgrid=None, - kstpkper=None, - totim=None, - attrib_name=None + self, gdf=None, modelgrid=None, kstpkper=None, totim=None, attrib_name=None ): """ Generate a GeoDataFrame with data from a LayerFile instance @@ -342,7 +337,7 @@ def to_shapefile( """ warnings.warn( "to_shapefile() is deprecated and is being replaced by to_geo_dataframe()", - DeprecationWarning + DeprecationWarning, ) plotarray = np.atleast_3d( self.get_data(kstpkper=kstpkper, totim=totim, mflay=mflay).transpose() diff --git a/flopy/utils/util_array.py b/flopy/utils/util_array.py index 498092eb0..2e4175209 100644 --- a/flopy/utils/util_array.py +++ b/flopy/utils/util_array.py @@ -1122,7 +1122,6 @@ def to_geo_dataframe(self, gdf=None, kper=0, forgive=False, **kwargs): gdf = u3d.to_geo_dataframe(gdf=gdf, forgive=forgive, name=name, **kwargs) return gdf - def get_zero_3d(self, kper): name = f"{self.name_base}{kper + 1}(filled zero)" return Util3d( From 1d799b99eae80a9a02eedd3d9a63b6ed1eb2398f Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Thu, 7 Aug 2025 16:58:10 -0700 Subject: [PATCH 09/13] spelling bee updates --- flopy/mbase.py | 2 +- flopy/mf6/mfmodel.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/flopy/mbase.py b/flopy/mbase.py index 58d3145f4..ce7589f03 100644 --- a/flopy/mbase.py +++ b/flopy/mbase.py @@ -762,7 +762,7 @@ def _output_msg(self, i, add=True): def to_geo_dataframe(self, gdf=None, kper=0): """ Method to build a Geodataframe from model inputs. Note: transient data - will only be exported for a sinlge stress period. + will only be exported for a single stress period. Parameters ---------- diff --git a/flopy/mf6/mfmodel.py b/flopy/mf6/mfmodel.py index 907a33454..0cb7e321a 100644 --- a/flopy/mf6/mfmodel.py +++ b/flopy/mf6/mfmodel.py @@ -792,7 +792,7 @@ def output(self): def to_geo_dataframe(self, gdf=None, kper=0): """ Method to build a Geodataframe from model inputs. Note: transient data - will only be exported for a sinlge stress period. + will only be exported for a single stress period. Parameters ---------- From f3089e6feba6376a90e6abca716957cbd1e32d3d Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Tue, 12 Aug 2025 16:24:48 -0700 Subject: [PATCH 10/13] Geodataframe support for modpath output --- flopy/mf6/mfpackage.py | 4 +- flopy/modflow/mfhfb.py | 3 +- flopy/utils/modpathfile.py | 102 +++++++++----- flopy/utils/particletrackfile.py | 224 +++++++++++++++++-------------- 4 files changed, 197 insertions(+), 136 deletions(-) diff --git a/flopy/mf6/mfpackage.py b/flopy/mf6/mfpackage.py index c98d687d9..d2b0af8cd 100644 --- a/flopy/mf6/mfpackage.py +++ b/flopy/mf6/mfpackage.py @@ -12,6 +12,7 @@ from ..pakbase import PackageInterface from ..utils import datautil from ..utils.check import mf6check +from ..utils.utl_import import import_optional_dependency from ..version import __version__ from .coordinates import modeldimensions from .data import ( @@ -2106,8 +2107,7 @@ def to_geo_dataframe(self, gdf=None, kper=0, sparse=False, **kwargs): modelgrid = self.parent.modelgrid if modelgrid is not None: if self.package_type == "hfb": - import geopandas as gpd - + gpd = import_optional_dependency("geopandas") from ..plot.plotutil import hfb_data_to_linework recarray = self.stress_period_data.data[kper] diff --git a/flopy/modflow/mfhfb.py b/flopy/modflow/mfhfb.py index d695dbedd..81df8a55c 100644 --- a/flopy/modflow/mfhfb.py +++ b/flopy/modflow/mfhfb.py @@ -14,6 +14,7 @@ from ..pakbase import Package from ..utils.flopy_io import line_parse from ..utils.recarray_utils import create_empty_recarray +from ..utils.utl_import import import_optional_dependency from .mfparbc import ModflowParBc as mfparbc @@ -200,7 +201,7 @@ def to_geo_dataframe(self, **kwargs): GeoDataFrame """ - import geopandas as gpd + gpd = import_optional_dependency("geopandas") lines = self._get_hfb_lines() geo_interface = {"type": "FeatureCollection"} diff --git a/flopy/utils/modpathfile.py b/flopy/utils/modpathfile.py index 6fbb930d2..7e10521af 100644 --- a/flopy/utils/modpathfile.py +++ b/flopy/utils/modpathfile.py @@ -12,6 +12,7 @@ from flopy.utils.particletrackfile import ParticleTrackFile from ..utils.flopy_io import loadtxt +from ..utils.utl_import import import_optional_dependency class ModpathFile(ParticleTrackFile): @@ -671,6 +672,61 @@ def get_destination_endpoint_data(self, dest_cells, source=False): inds = np.isin(raslice, dest_cells) return data[inds].copy().view(np.recarray) + def to_geo_dataframe( + self, + modelgrid, + data=None, + direction="ending", + ): + """ + Create a geodataframe of particle starting / ending locations. + + Parameters + ---------- + modelgrid : flopy.discretization.grid instance + Used to scale and rotate Global x,y,z values in MODPATH Endpoint + file. + data : np.recarray + Record array of same form as that returned by EndpointFile.get_alldata. + (if none, EndpointFile.get_alldata() is exported). + direction : str + String defining if 'starting' or 'ending' particle locations should be + considered. (default is 'ending') + """ + from ..utils import geometry + gpd = import_optional_dependency("geopandas") + shapely_geo = import_optional_dependency("shapely.geometry") + if data is None: + data = self.get_alldata() + + if direction.lower() == "ending": + xcol, ycol, zcol = "x", "y", "z" + elif direction.lower() == "starting": + xcol, ycol, zcol = "x0", "y0", "z0" + else: + raise Exception( + 'flopy.map.plot_endpoint direction must be "ending" or "starting".' + ) + x, y = geometry.transform( + data[xcol], + data[ycol], + xoff=modelgrid.xoffset, + yoff=modelgrid.yoffset, + angrot_radians=modelgrid.angrot_radians, + ) + z = data[zcol] + + geoms = [shapely_geo.Point(p) for p in zip(z, y, z)] + gdf = gpd.GeoDataFrame(data, geometry=geoms, crs=modelgrid.crs) + + # adjust to 1 based node numbers + for col in list(gdf): + if col in self.kijnames: + gdf[col] += 1 + + return gdf + + def write_shapefile( self, data=None, @@ -715,44 +771,20 @@ def write_shapefile( - ``epsg`` (int): use ``crs`` instead. """ - from ..discretization import StructuredGrid - from ..export.shapefile_utils import recarray2shp - from ..utils import geometry - from ..utils.geometry import Point - + import warnings + warnings.warn( + "write_shapefile is Deprecated, please use to_geo_dataframe() in the future" + ) epd = (data if data is not None else endpoint_data).copy() - if epd is None: - epd = self.get_alldata() + gdf = self.to_geo_dataframe(modelgrid=mg, data=epd, direction=direction) - if direction.lower() == "ending": - xcol, ycol, zcol = "x", "y", "z" - elif direction.lower() == "starting": - xcol, ycol, zcol = "x0", "y0", "z0" - else: - raise Exception( - 'flopy.map.plot_endpoint direction must be "ending" or "starting".' - ) - if mg is None: - raise ValueError("A modelgrid object was not provided.") - - if isinstance(mg, StructuredGrid): - x, y = geometry.transform( - epd[xcol], - epd[ycol], - xoff=mg.xoffset, - yoff=mg.yoffset, - angrot_radians=mg.angrot_radians, - ) - else: - x, y = mg.get_coords(epd[xcol], epd[ycol]) - z = epd[zcol] + if crs is not None: + if gdf.crs is None: + gdf = gdf.set_crs(crs) + else: + gdf = gdf.to_crs(crs) - geoms = [Point(x[i], y[i], z[i]) for i in range(len(epd))] - # convert back to one-based - for n in self.kijnames: - if n in epd.dtype.names: - epd[n] += 1 - recarray2shp(epd, geoms, shpname=shpname, crs=crs, **kwargs) + gdf.to_file(shpname) class TimeseriesFile(ModpathFile): diff --git a/flopy/utils/particletrackfile.py b/flopy/utils/particletrackfile.py index c5df347c0..42777bf3b 100644 --- a/flopy/utils/particletrackfile.py +++ b/flopy/utils/particletrackfile.py @@ -9,6 +9,7 @@ import numpy as np from numpy.lib.recfunctions import stack_arrays +from .utl_import import import_optional_dependency MIN_PARTICLE_TRACK_DTYPE = np.dtype( [ @@ -183,19 +184,20 @@ def intersect(self, cells, to_recarray) -> np.recarray: """Find intersection of pathlines with cells.""" pass - def write_shapefile( + def to_geo_dataframe( self, + modelgrid, data=None, one_per_particle=True, direction="ending", - shpname="endpoints.shp", - mg=None, - crs=None, - **kwargs, ): """ - Write particle track data to a shapefile. + Create a geodataframe of particle tracks. + Parameters + ---------- + modelgrid : flopy.discretization.Grid instance + Used to scale and rotate Global x, y, z values. data : np.recarray Record array of same form as that returned by get_alldata(). (if none, get_alldata() is exported). @@ -209,120 +211,146 @@ def write_shapefile( String defining if starting or ending particle locations should be included in shapefile attribute information. Only used if one_per_particle=False. (default is 'ending') - shpname : str - File path for shapefile - mg : flopy.discretization.grid instance - Used to scale and rotate Global x,y,z values. - crs : pyproj.CRS, int, str, optional - Coordinate reference system (CRS) for the model grid - (must be projected; geographic CRS are not supported). - The value can be anything accepted by - :meth:`pyproj.CRS.from_user_input() `, - such as an authority string (eg "EPSG:26916") or a WKT string. - kwargs : keyword arguments to flopy.export.shapefile_utils.recarray2shp - - .. deprecated:: 3.5 - The following keyword options will be removed for FloPy 3.6: - - - ``epsg`` (int): use ``crs`` instead. + Returns + ------- + GeoDataFrame """ - from ..export.shapefile_utils import recarray2shp from . import geometry - from .geometry import LineString + shapely_geo = import_optional_dependency("shapely.geometry") + gpd = import_optional_dependency("geopandas") - series = data - if series is None: - series = self._data.view(np.recarray) + if data is None: + data = self._data.view(np.recarray) else: # convert pathline list to a single recarray - if isinstance(series, list): - s = series[0] + if isinstance(data, list): + s = data[0] print(s.dtype) - for n in range(1, len(series)): - s = stack_arrays((s, series[n])) - series = s.view(np.recarray) - - series = series.copy() - series.sort(order=["particleid", "time"]) + for n in range(1, len(data)): + s = stack_arrays((s, data[n])) + data = s.view(np.recarray) - if mg is None: - raise ValueError("A modelgrid object was not provided.") + data = data.copy() + data.sort(order=["particleid", "time"]) - particles = np.unique(series.particleid) + particles = np.unique(data.particleid) geoms = [] - # create dtype with select attributes in pth - names = series.dtype.names - dtype = [] - atts = ["particleid", "particlegroup", "time", "k", "i", "j", "node"] - for att in atts: - if att in names: - t = np.int32 - if att == "time": - t = np.float32 - dtype.append((att, t)) - dtype = np.dtype(dtype) - - # reset names to the selected names in the created dtype - names = dtype.names - - # 1 geometry for each path + # create a dict of attrs? + headings = ["particleid", "particlegroup", "time", "k", "i", "j", "node"] + attrs = [] + for h in headings: + if h in data.dtype.names: + attrs.append(h) + if one_per_particle: - loc_inds = 0 + dfdata = {a: [] for a in attrs} if direction == "ending": - loc_inds = -1 - - sdata = [] - for pid in particles: - ra = series[series.particleid == pid] + idx = -1 + else: + idx = 0 + + for p in particles: + ra = data[data.particleid == p] + for k in dfdata.keys(): + if k == "time": + dfdata[k].append(np.max(ra[k])) + else: + dfdata[k].append(ra[k][idx]) x, y = geometry.transform( - ra.x, ra.y, mg.xoffset, mg.yoffset, mg.angrot_radians + ra.x, ra.y, modelgrid.xoffset, modelgrid.yoffset, modelgrid.angrot_radians ) z = ra.z - geoms.append(LineString(list(zip(x, y, z)))) - - t = [pid] - if "particlegroup" in names: - t.append(ra.particlegroup[0]) - t.append(ra.time.max()) - if "k" in names: - t.append(ra.k[loc_inds]) - if "node" in names: - t.append(ra.node[loc_inds]) - else: - if "i" in names: - t.append(ra.i[loc_inds]) - if "j" in names: - t.append(ra.j[loc_inds]) - sdata.append(tuple(t)) - - sdata = np.array(sdata, dtype=dtype).view(np.recarray) - - # geometry for each row in PathLine file + + line = list(zip(x, y, z)) + geoms.append(shapely_geo.LineString(line)) + else: - dtype = series.dtype - sdata = [] + dfdata = {a: [] for a in attrs} for pid in particles: - ra = series[series.particleid == pid] - if mg is not None: - x, y = geometry.transform( - ra.x, ra.y, mg.xoffset, mg.yoffset, mg.angrot_radians - ) - else: - x, y = geometry.transform(ra.x, ra.y, 0, 0, 0) + ra = data[data.particleid == pid] + x, y = geometry.transform( + ra.x, ra.y, modelgrid.xoffset, modelgrid.yoffset, modelgrid.angrot_radians + ) z = ra.z geoms += [ - LineString([(x[i - 1], y[i - 1], z[i - 1]), (x[i], y[i], z[i])]) + shapely_geo.LineString([(x[i - 1], y[i - 1], z[i - 1]), (x[i], y[i], z[i])]) for i in np.arange(1, (len(ra))) ] - sdata += ra[1:].tolist() - sdata = np.array(sdata, dtype=dtype).view(np.recarray) + for k in dfdata.keys(): + dfdata[k].extend(ra[k][1:]) + + # now create a geodataframe + gdf = gpd.GeoDataFrame(dfdata, geometry=geoms, crs=modelgrid.crs) + + # adjust to 1 based node numbers + for col in list(gdf): + if col in self.kijnames: + gdf[col] += 1 + + return gdf - # convert back to one-based - for n in set(self.kijnames).intersection(set(sdata.dtype.names)): - sdata[n] += 1 + def write_shapefile( + self, + data=None, + one_per_particle=True, + direction="ending", + shpname="endpoints.shp", + mg=None, + crs=None, + **kwargs, + ): + """ + Write particle track data to a shapefile. + + data : np.recarray + Record array of same form as that returned by + get_alldata(). (if none, get_alldata() is exported). + one_per_particle : boolean (default True) + True writes a single LineString with a single set of attribute + data for each particle. False writes a record/geometry for each + pathline segment (each row in the PathLine file). This option can + be used to visualize attribute information (time, model layer, + etc.) across a pathline in a GIS. + direction : str + String defining if starting or ending particle locations should be + included in shapefile attribute information. Only used if + one_per_particle=False. (default is 'ending') + shpname : str + File path for shapefile + mg : flopy.discretization.grid instance + Used to scale and rotate Global x,y,z values. + crs : pyproj.CRS, int, str, optional + Coordinate reference system (CRS) for the model grid + (must be projected; geographic CRS are not supported). + The value can be anything accepted by + :meth:`pyproj.CRS.from_user_input() `, + such as an authority string (eg "EPSG:26916") or a WKT string. + kwargs : keyword arguments to flopy.export.shapefile_utils.recarray2shp + + .. deprecated:: 3.5 + The following keyword options will be removed for FloPy 3.6: + + - ``epsg`` (int): use ``crs`` instead. + + """ + import warnings + warnings.warn( + "write_shapefile will be Deprecated, please use to_geo_dataframe()", + DeprecationWarning + ) + gdf = self.to_geo_dataframe( + modelgrid=mg, + data=data, + one_per_particle=one_per_particle, + direction=direction + ) + if crs is not None: + if gdf.crs is None: + gdf = gdf.set_crs(crs) + else: + gdf = gdf.to_crs(crs) - # write the final recarray to a shapefile - recarray2shp(sdata, geoms, shpname=shpname, crs=crs, **kwargs) + gdf.to_file(shpname) \ No newline at end of file From 34424af2971ef55a00a34fb17249da98090c1689 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Tue, 12 Aug 2025 16:27:56 -0700 Subject: [PATCH 11/13] woof! --- flopy/utils/modpathfile.py | 3 ++- flopy/utils/particletrackfile.py | 25 +++++++++++++++++++------ 2 files changed, 21 insertions(+), 7 deletions(-) diff --git a/flopy/utils/modpathfile.py b/flopy/utils/modpathfile.py index 7e10521af..29f8d60c2 100644 --- a/flopy/utils/modpathfile.py +++ b/flopy/utils/modpathfile.py @@ -694,6 +694,7 @@ def to_geo_dataframe( considered. (default is 'ending') """ from ..utils import geometry + gpd = import_optional_dependency("geopandas") shapely_geo = import_optional_dependency("shapely.geometry") if data is None: @@ -726,7 +727,6 @@ def to_geo_dataframe( return gdf - def write_shapefile( self, data=None, @@ -772,6 +772,7 @@ def write_shapefile( """ import warnings + warnings.warn( "write_shapefile is Deprecated, please use to_geo_dataframe() in the future" ) diff --git a/flopy/utils/particletrackfile.py b/flopy/utils/particletrackfile.py index 42777bf3b..0e3a4ead4 100644 --- a/flopy/utils/particletrackfile.py +++ b/flopy/utils/particletrackfile.py @@ -9,6 +9,7 @@ import numpy as np from numpy.lib.recfunctions import stack_arrays + from .utl_import import import_optional_dependency MIN_PARTICLE_TRACK_DTYPE = np.dtype( @@ -217,6 +218,7 @@ def to_geo_dataframe( GeoDataFrame """ from . import geometry + shapely_geo = import_optional_dependency("shapely.geometry") gpd = import_optional_dependency("geopandas") @@ -260,7 +262,11 @@ def to_geo_dataframe( dfdata[k].append(ra[k][idx]) x, y = geometry.transform( - ra.x, ra.y, modelgrid.xoffset, modelgrid.yoffset, modelgrid.angrot_radians + ra.x, + ra.y, + modelgrid.xoffset, + modelgrid.yoffset, + modelgrid.angrot_radians, ) z = ra.z @@ -272,11 +278,17 @@ def to_geo_dataframe( for pid in particles: ra = data[data.particleid == pid] x, y = geometry.transform( - ra.x, ra.y, modelgrid.xoffset, modelgrid.yoffset, modelgrid.angrot_radians + ra.x, + ra.y, + modelgrid.xoffset, + modelgrid.yoffset, + modelgrid.angrot_radians, ) z = ra.z geoms += [ - shapely_geo.LineString([(x[i - 1], y[i - 1], z[i - 1]), (x[i], y[i], z[i])]) + shapely_geo.LineString( + [(x[i - 1], y[i - 1], z[i - 1]), (x[i], y[i], z[i])] + ) for i in np.arange(1, (len(ra))) ] for k in dfdata.keys(): @@ -337,15 +349,16 @@ def write_shapefile( """ import warnings + warnings.warn( "write_shapefile will be Deprecated, please use to_geo_dataframe()", - DeprecationWarning + DeprecationWarning, ) gdf = self.to_geo_dataframe( modelgrid=mg, data=data, one_per_particle=one_per_particle, - direction=direction + direction=direction, ) if crs is not None: if gdf.crs is None: @@ -353,4 +366,4 @@ def write_shapefile( else: gdf = gdf.to_crs(crs) - gdf.to_file(shpname) \ No newline at end of file + gdf.to_file(shpname) From d8311e235e3cd815c26c38c02f8952e364e32188 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Wed, 13 Aug 2025 14:47:56 -0700 Subject: [PATCH 12/13] updates for existing modpath export testing --- autotest/test_modpathfile.py | 26 ++++++++++---------------- flopy/utils/modpathfile.py | 2 +- flopy/utils/particletrackfile.py | 2 +- 3 files changed, 12 insertions(+), 18 deletions(-) diff --git a/autotest/test_modpathfile.py b/autotest/test_modpathfile.py index 2bedeeb3f..a891c6b8d 100644 --- a/autotest/test_modpathfile.py +++ b/autotest/test_modpathfile.py @@ -306,11 +306,11 @@ def test_get_destination_endpoint_data( ) -@pytest.mark.parametrize("longfieldname", [True, False]) + @requires_exe("mf6", "mp7") -@requires_pkg("pyshp", "shapely", name_map={"pyshp": "shapefile"}) -def test_write_shapefile(function_tmpdir, mp7_small, longfieldname): - from shapefile import Reader +@requires_pkg("geopandas", "shapely") +def test_write_shapefile(function_tmpdir, mp7_small): + import geopandas as gpd # setup and run model, then copy outputs to function_tmpdir sim, forward_model_name, _, _, _ = mp7_small @@ -330,13 +330,7 @@ def test_write_shapefile(function_tmpdir, mp7_small, longfieldname): # define shapefile path shp_file = ws / "pathlines.shp" - # add a column to the pathline recarray - fieldname = "newfield" + ("longname" if longfieldname else "") - fieldval = "x" - pathlines = [ - rfn.append_fields(pl, fieldname, list(repeat(fieldval, len(pl))), dtypes="|S1") - for pl in pathlines - ] + pline_names = [name[:10] for name in pathlines[0].dtype.names] # write the pathline recarray to shapefile pathline_file.write_shapefile( @@ -350,8 +344,8 @@ def test_write_shapefile(function_tmpdir, mp7_small, longfieldname): assert shp_file.is_file() # load shapefile - with Reader(shp_file) as reader: - fieldnames = [f[0] for f in reader.fields[1:]] - fieldname = "newfiname_" if longfieldname else fieldname - assert fieldname in fieldnames - assert all(r[fieldname] == fieldval for r in reader.iterRecords()) + gdf = gpd.read_file(shp_file) + fieldnames = list(gdf) + for fname in pline_names: + assert fname in fieldnames + diff --git a/flopy/utils/modpathfile.py b/flopy/utils/modpathfile.py index 29f8d60c2..938206de8 100644 --- a/flopy/utils/modpathfile.py +++ b/flopy/utils/modpathfile.py @@ -717,7 +717,7 @@ def to_geo_dataframe( ) z = data[zcol] - geoms = [shapely_geo.Point(p) for p in zip(z, y, z)] + geoms = [shapely_geo.Point(p) for p in zip(x, y, z)] gdf = gpd.GeoDataFrame(data, geometry=geoms, crs=modelgrid.crs) # adjust to 1 based node numbers diff --git a/flopy/utils/particletrackfile.py b/flopy/utils/particletrackfile.py index 0e3a4ead4..483ce1570 100644 --- a/flopy/utils/particletrackfile.py +++ b/flopy/utils/particletrackfile.py @@ -274,7 +274,7 @@ def to_geo_dataframe( geoms.append(shapely_geo.LineString(line)) else: - dfdata = {a: [] for a in attrs} + dfdata = {a: [] for a in data.dtype.names} for pid in particles: ra = data[data.particleid == pid] x, y = geometry.transform( From 862ec36e8bcabad6463ddc48d7f9ec79b5760ba5 Mon Sep 17 00:00:00 2001 From: jlarsen-usgs Date: Wed, 13 Aug 2025 14:49:38 -0700 Subject: [PATCH 13/13] woof --- autotest/test_modpathfile.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/autotest/test_modpathfile.py b/autotest/test_modpathfile.py index a891c6b8d..9f667fac8 100644 --- a/autotest/test_modpathfile.py +++ b/autotest/test_modpathfile.py @@ -306,7 +306,6 @@ def test_get_destination_endpoint_data( ) - @requires_exe("mf6", "mp7") @requires_pkg("geopandas", "shapely") def test_write_shapefile(function_tmpdir, mp7_small): @@ -348,4 +347,3 @@ def test_write_shapefile(function_tmpdir, mp7_small): fieldnames = list(gdf) for fname in pline_names: assert fname in fieldnames -