fix(model_attributes_to_shapefile): use modelgrid kwarg if provided

Author: wpbonelli

autotest/test_shapefile_utils.py

@@ -3,10 +3,12 @@
 """
 
 import numpy as np
+import pytest
 from modflow_devtools.markers import requires_pkg
 
 import flopy
 from flopy.discretization import StructuredGrid, UnstructuredGrid, VertexGrid
+from flopy.discretization.grid import Grid
 from flopy.export.shapefile_utils import model_attributes_to_shapefile, shp2recarray
 from flopy.utils.crs import get_shapefile_crs
 
@@ -49,6 +51,75 @@ def test_model_attributes_to_shapefile(example_data_path, function_tmpdir):
     assert shpfile_path.exists()
 
 
+@requires_pkg("geopandas")
+def test_model_attributes_to_shapefile_modelgrid_kwarg(function_tmpdir):
+    """Repro https://github.com/modflowpy/flopy/issues/2744
+
+    The modelgrid kwarg to model_attributes_to_shapefile should be used as
+    the geometry source if provided, overriding the model's own modelgrid.
+    """
+    import warnings
+
+    nrow, ncol = 3, 4
+    delr = np.ones(ncol) * 10.0
+    delc = np.ones(nrow) * 10.0
+    crs = 26916
+
+    # Model without a DIS package: modelgrid is a bare Grid with no geometry.
+    # Without the fix, this reproduces the reported error:
+    #   TypeError: Grid.to_geodataframe() missing 1 required positional argument: 'features'
+    sim = flopy.mf6.MFSimulation(sim_name="test", sim_ws=str(function_tmpdir))
+    gwf = flopy.mf6.ModflowGwf(sim, modelname="test")
+    mg = StructuredGrid(delr=delr, delc=delc, nlay=1, crs=crs)
+    shpfile = function_tmpdir / "test_no_dis.shp"
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", DeprecationWarning)
+        model_attributes_to_shapefile(shpfile, gwf, modelgrid=mg)
+    assert shpfile.exists()
+
+    # The error in the issue could also have been avoided by retrieving the
+    # modelgrid after attaching the DIS to the model. Before DIS is added,
+    # gwf.modelgrid is the bare base Grid class. After DIS is attached, it
+    # becomes a StructuredGrid whose to_geodataframe() needs no 'features'.
+    sim1 = flopy.mf6.MFSimulation(sim_name="test_b", sim_ws=str(function_tmpdir))
+    gwf1 = flopy.mf6.ModflowGwf(sim1, modelname="test_b")
+    assert isinstance(gwf1.modelgrid, Grid)
+    flopy.mf6.ModflowGwfdis(gwf1, nlay=1, nrow=nrow, ncol=ncol)
+    mg1 = gwf1.modelgrid
+    assert isinstance(mg1, StructuredGrid)
+    shpfile1 = function_tmpdir / "test_dis_first.shp"
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", DeprecationWarning)
+        model_attributes_to_shapefile(
+            shpfile1, gwf1, package_names=["dis"], modelgrid=mg1
+        )
+    assert shpfile1.exists()
+
+    # Without a modelgrid kwarg and no DIS, to_geodataframe raises a clear
+    # AttributeError rather than the cryptic TypeError about 'features'.
+    sim_err = flopy.mf6.MFSimulation(sim_name="test_err", sim_ws=str(function_tmpdir))
+    gwf_err = flopy.mf6.ModflowGwf(sim_err, modelname="test_err")
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", DeprecationWarning)
+        with pytest.raises(AttributeError, match="discretization package"):
+            model_attributes_to_shapefile(function_tmpdir / "test_err.shp", gwf_err)
+
+    # Model with a DIS package but a different modelgrid passed via kwarg:
+    # the kwarg modelgrid's CRS should take precedence.
+    sim2 = flopy.mf6.MFSimulation(sim_name="test2", sim_ws=str(function_tmpdir))
+    gwf2 = flopy.mf6.ModflowGwf(sim2, modelname="test2")
+    flopy.mf6.ModflowGwfdis(gwf2, nlay=1, nrow=nrow, ncol=ncol)
+
+    mg2 = StructuredGrid(delr=delr, delc=delc, nlay=1, crs=crs)
+    shpfile2 = function_tmpdir / "test_with_dis.shp"
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", DeprecationWarning)
+        model_attributes_to_shapefile(
+            shpfile2, gwf2, package_names=["dis"], modelgrid=mg2
+        )
+    assert shpfile2.exists()
+
+
 @requires_pkg("geopandas")
 def test_create_geodataframe(
     minimal_unstructured_grid_info, minimal_vertex_grid_info, function_tmpdir

flopy/export/shapefile_utils.py

@@ -220,10 +220,15 @@ def model_attributes_to_shapefile(
     else:
         package_names = [pak.name[0] for pak in ml.packagelist]
 
-    gdf = ml.to_geodataframe(package_names=package_names, shorten_attr=True)
+    modelgrid = kwargs.pop("modelgrid", None)
+    init_gdf = modelgrid.to_geodataframe() if modelgrid is not None else None
+    gdf = ml.to_geodataframe(
+        gdf=init_gdf, package_names=package_names, shorten_attr=True
+    )
 
     if array_dict:
-        modelgrid = ml.modelgrid
+        if modelgrid is None:
+            modelgrid = ml.modelgrid
         for name, array in array_dict.items():
             if modelgrid.grid_type() == "unstructured":
                 gdf[name] = array.ravel()

flopy/mbase.py

@@ -790,13 +790,20 @@ def to_geodataframe(self, gdf=None, kper=0, package_names=None, shorten_attr=Fal
             gdf : GeoDataFrame
         """
         if gdf is None:
+            from .discretization.grid import Grid
+
             modelgrid = self.modelgrid
-            if modelgrid is not None:
-                gdf = modelgrid.to_geodataframe()
-            else:
+            if modelgrid is None:
                 raise AttributeError(
                     "model does not have a grid instance, please supply a geodataframe"
                 )
+            if type(modelgrid) is Grid:
+                raise AttributeError(
+                    "model does not have a discretization package; cannot build a "
+                    "GeoDataFrame without geometry. Attach a discretization package "
+                    "or supply a gdf."
+                )
+            gdf = modelgrid.to_geodataframe()
 
         if package_names is None:
             package_names = [pak.name[0] for pak in self.packagelist]

flopy/mf6/mfmodel.py

@@ -816,14 +816,21 @@ def to_geodataframe(self, gdf=None, kper=0, package_names=None, shorten_attr=Fal
             gdf : GeoDataFrame
         """
         if gdf is None:
+            from ..discretization.grid import Grid
+
             modelgrid = self.modelgrid
-            if modelgrid is not None:
-                gdf = modelgrid.to_geodataframe()
-            else:
+            if modelgrid is None:
                 raise AttributeError(
                     "model does not have a grid instance, "
                     "please supply a geodataframe"
                 )
+            if type(modelgrid) is Grid:
+                raise AttributeError(
+                    "model does not have a discretization package; cannot build a "
+                    "GeoDataFrame without geometry. Attach a discretization package "
+                    "or supply a gdf."
+                )
+            gdf = modelgrid.to_geodataframe()
 
         if package_names is None:
             package_names = [pak.name[0] for pak in self.packagelist]