fix(model_attributes_to_shapefile): pass modelgrid

Author: wpbonelli

autotest/test_shapefile_utils.py

@@ -7,6 +7,7 @@
 
 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 +50,66 @@ 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()
+
+    # 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()