update(GeoSpatialCollection, Gridgen): Added support for lists of sha…

…pely objects (#1565)

* update(GeoSpatialCollection, Gridgen): Added lists of shapely objects to input type 

* updated Gridgen to accept StructuredGrid objects. Begin deprecation of DIS object support in Gridgen
* updated flopy3_gridgen.ipynb

* Fix typo in Gridgen, update autotests

* update to test_geospatial_util.py

autotest/test_export.py

@@ -1936,7 +1936,7 @@ def test_vtk_export_disu_model(tmpdir):
         botm=botm,
     )
 
-    g = Gridgen(dis5, model_ws=str(tmpdir))
+    g = Gridgen(ml5.modelgrid, model_ws=str(tmpdir))
 
     rf0shp = os.path.join(str(tmpdir), "rf0")
     xmin = 7 * delr

autotest/test_geospatial_util.py

@@ -469,11 +469,12 @@ def test_mixed_collection(
     shapetype = gc1.shapetype
     shp = gc1.shape
     shply = gc1.shapely
+    lshply = list(gc1.shapely)
     points = gc1.points
     geojson = gc1.geojson
     fp_geo = gc1.flopy_geometry
 
-    collections = [shp, shply, points, geojson, fp_geo]
+    collections = [shp, shply, lshply, points, geojson, fp_geo]
     for col in collections:
         gc2 = GeoSpatialCollection(col, shapetype)
 
@@ -489,3 +490,4 @@ def test_mixed_collection(
                 is_equal = gi2 == gi1[ix]
 
             assert is_equal, "GeoSpatialCollection conversion error"
+

autotest/test_gridgen.py

@@ -46,7 +46,7 @@ def test_mf6disv(tmpdir):
     )
 
     # Create and build the gridgen model with a refined area in the middle
-    g = Gridgen(dis, model_ws=str(tmpdir))
+    g = Gridgen(gwf.modelgrid, model_ws=str(tmpdir))
     polys = [Polygon([(4, 4), (6, 4), (6, 6), (4, 6)])]
     g.add_refinement_features(polys, "polygon", 3, range(nlay))
     g.build()
@@ -182,7 +182,7 @@ def test_mf6disu(tmpdir):
     )
 
     # Create and build the gridgen model with a refined area in the middle
-    g = Gridgen(dis, model_ws=str(tmpdir))
+    g = Gridgen(gwf.modelgrid, model_ws=str(tmpdir))
     polys = [Polygon([(4, 4), (6, 4), (6, 6), (4, 6)])]
     g.add_refinement_features(polys, "polygon", 3, layers=[0])
     g.build()
@@ -346,7 +346,7 @@ def test_mfusg(tmpdir):
     )
 
     # Create and build the gridgen model with a refined area in the middle
-    g = Gridgen(dis, model_ws=str(tmpdir))
+    g = Gridgen(m.modelgrid, model_ws=str(tmpdir))
     polys = [Polygon([(4, 4), (6, 4), (6, 6), (4, 6)])]
     g.add_refinement_features(polys, "polygon", 3, layers=[0])
     g.build()
@@ -514,10 +514,10 @@ def test_gridgen(tmpdir):
 
     gridgen = Path(which("gridgen")).name
     ws = str(tmpdir)
-    g = Gridgen(dis5, model_ws=ws, exe_name=gridgen)
-    g6 = Gridgen(dis6, model_ws=ws, exe_name=gridgen)
+    g = Gridgen(ms.modelgrid, model_ws=ws, exe_name=gridgen)
+    g6 = Gridgen(gwf.modelgrid, model_ws=ws, exe_name=gridgen)
     gu = Gridgen(
-        dis_usg,
+        ms_u.modelgrid,
         model_ws=ws,
         exe_name=gridgen,
         vertical_pass_through=True,

autotest/test_headufile.py

@@ -45,7 +45,7 @@ def mfusg_model(module_tmpdir):
     )
 
     # Create and build the gridgen model with a refined area in the middle
-    g = Gridgen(dis, model_ws=str(module_tmpdir))
+    g = Gridgen(m.modelgrid, model_ws=str(module_tmpdir))
 
     polys = [Polygon([(4, 4), (6, 4), (6, 6), (4, 6)])]
     g.add_refinement_features(polys, "polygon", 3, layers=[0])

flopy/utils/geospatial_utils.py

@@ -260,6 +260,10 @@ def __init__(self, obj, shapetype=None):
         self.__shapefile = import_optional_dependency(
             "shapefile", errors="silent"
         )
+        shapely_geo = import_optional_dependency(
+            "shapely.geometry", errors="silent"
+        )
+
         self.__obj = obj
         self.__collection = []
         self._geojson = None
@@ -283,6 +287,21 @@ def __init__(self, obj, shapetype=None):
                     for shape in obj:
                         self.__collection.append(GeoSpatialUtil(shape))
 
+            if shapely_geo is not None:
+                if isinstance(
+                    obj[0],
+                    (
+                        shapely_geo.Point,
+                        shapely_geo.MultiPoint,
+                        shapely_geo.Polygon,
+                        shapely_geo.MultiPolygon,
+                        shapely_geo.LineString,
+                        shapely_geo.MultiLineString,
+                    ),
+                ):
+                    for shape in obj:
+                        self.__collection.append(GeoSpatialUtil(shape))
+
             if not self.__collection:
                 if shapetype is None:
                     err = "a list of shapetypes must be provided"
@@ -324,7 +343,6 @@ def __init__(self, obj, shapetype=None):
                 for geom in obj.geometries:
                     self.__collection.append(GeoSpatialUtil(geom))
 
-        shapely_geo = import_optional_dependency("shapely.geometry")
         if shapely_geo is not None:
             if isinstance(
                 obj,

flopy/utils/gridgen.py

@@ -1,12 +1,15 @@
 import os
 import subprocess
+import warnings
 
 import numpy as np
 
+from ..discretization import StructuredGrid
 from ..export.shapefile_utils import shp2recarray
 from ..mbase import which
 from ..mf6.modflow import ModflowGwfdis
 from ..mfusg.mfusgdisu import MfUsgDisU
+from ..modflow import ModflowDis
 from ..utils import import_optional_dependency
 from .util_array import Util2d  # read1d,
 
@@ -172,8 +175,10 @@ class Gridgen:
 
     Parameters
     ----------
-    dis : flopy.modflow.ModflowDis
-        Flopy discretization object
+    modelgrid : flopy.discretization.StructuredGrid
+        Flopy StructuredGrid object. Note this also accepts ModflowDis and
+        ModflowGwfdis objects, however it is deprecated and support will be
+        removed in version 3.3.7
     model_ws : str
         workspace location for creating gridgen files (default is '.')
     exe_name : str
@@ -204,28 +209,41 @@ class Gridgen:
 
     def __init__(
         self,
-        dis,
+        modelgrid,
         model_ws=".",
         exe_name="gridgen",
         surface_interpolation="replicate",
         vertical_pass_through=False,
         **kwargs,
     ):
-        self.dis = dis
-        if isinstance(dis, ModflowGwfdis):
-            self.nlay = self.dis.nlay.get_data()
-            self.nrow = self.dis.nrow.get_data()
-            self.ncol = self.dis.ncol.get_data()
-            self.modelgrid = self.dis.parent.modelgrid
+        if isinstance(modelgrid, StructuredGrid):
+            if modelgrid.top is None or modelgrid.botm is None:
+                raise AssertionError(
+                    "A complete modelgrid must be supplied to use Gridgen"
+                )
+
+            self.modelgrid = modelgrid
+
+        elif isinstance(modelgrid, (ModflowGwfdis, ModflowDis)):
+            warnings.warn(
+                "Supplying a dis object is deprecated, and support will be "
+                "removed in version 3.3.7. Please supply StructuredGrid."
+            )
+            # this is actually a DIS file
+            self.modelgrid = modelgrid.parent.modelgrid
+
         else:
-            self.nlay = self.dis.nlay
-            self.nrow = self.dis.nrow
-            self.ncol = self.dis.ncol
-            self.modelgrid = self.dis.parent.modelgrid
+            raise TypeError(
+                "A StructuredGrid object must be supplied to Gridgen"
+            )
+
+        self.nlay = self.modelgrid.nlay
+        self.nrow = self.modelgrid.nrow
+        self.ncol = self.modelgrid.ncol
 
         self.nodes = 0
         self.nja = 0
-        self.nodelay = np.zeros((self.nlay), dtype=int)
+        self.nodelay = np.zeros((self.nlay,), dtype=int)
         self._vertdict = {}
         self.model_ws = model_ws
         exe_name = which(exe_name)
@@ -1674,6 +1692,8 @@ def get_gridprops_unstructuredgrid(self):
         gridprops["ycenters"] = ycenters
         gridprops["top"] = top
         gridprops["botm"] = bot
+        gridprops["iac"] = self.get_iac()
+        gridprops["ja"] = self.get_ja()
 
         return gridprops
 
@@ -1783,7 +1803,7 @@ def _mfgrid_block(self):
         s += f"  NCOL = {self.ncol}\n"
 
         # delr
-        delr = self.dis.delr.array
+        delr = self.modelgrid.delr
         if delr.min() == delr.max():
             s += f"  DELR = CONSTANT {delr.min()}\n"
         else:
@@ -1792,7 +1812,7 @@ def _mfgrid_block(self):
             np.savetxt(fname, np.atleast_2d(delr))
 
         # delc
-        delc = self.dis.delc.array
+        delc = self.modelgrid.delc
         if delc.min() == delc.max():
             s += f"  DELC = CONSTANT {delc.min()}\n"
         else:
@@ -1801,7 +1821,7 @@ def _mfgrid_block(self):
             np.savetxt(fname, np.atleast_2d(delc))
 
         # top
-        top = self.dis.top.array
+        top = self.modelgrid.top
         if top.min() == top.max():
             s += f"  TOP = CONSTANT {top.min()}\n"
         else:
@@ -1810,12 +1830,9 @@ def _mfgrid_block(self):
             np.savetxt(fname, top)
 
         # bot
-        botm = self.dis.botm.array
+        botm = self.modelgrid.botm
         for k in range(self.nlay):
-            if isinstance(self.dis, ModflowGwfdis):
-                bot = botm[k]
-            else:
-                bot = botm[k]
+            bot = botm[k]
             if bot.min() == bot.max():
                 s += f"  BOTTOM LAYER {k + 1} = CONSTANT {bot.min()}\n"
             else: