update(GeoSpatialCollection): add support for GeoDataFrame objects (#…

…2063)

* update(GeoSpatialCollection): add support for geopandas GeoDataFrame objects

* update test_geospatial_util.py to include GeoDataFrame tests
* add geopandas as optional dependency

* linting

* * Added geo_dataframe property to grid objects
* Added convert_grid method to grid objects

* linting

* add ignore_holes flag to add_polygon

autotest/test_geospatial_util.py

@@ -360,7 +360,7 @@ def test_multilinestring(multilinestring):
         assert gi1 == gi2, "GeoSpatialUtil multilinestring conversion error"
 
 
-@requires_pkg("shapely", "geojson")
+@requires_pkg("shapely", "geojson", "geopandas")
 def test_polygon_collection(polygon, poly_w_hole, multipolygon):
     col = [
         Shape.from_geojson(polygon),
@@ -377,8 +377,9 @@ def test_polygon_collection(polygon, poly_w_hole, multipolygon):
     points = gc1.points
     geojson = gc1.geojson
     fp_geo = gc1.flopy_geometry
+    gdf = gc1.geo_dataframe
 
-    collections = [shp, shply, points, geojson, fp_geo]
+    collections = [shp, shply, points, geojson, fp_geo, gdf]
     for col in collections:
         gc2 = GeoSpatialCollection(col, shapetype)
 
@@ -410,8 +411,9 @@ def test_point_collection(point, multipoint):
     points = gc1.points
     geojson = gc1.geojson
     fp_geo = gc1.flopy_geometry
+    gdf = gc1.geo_dataframe
 
-    collections = [shp, shply, points, geojson, fp_geo]
+    collections = [shp, shply, points, geojson, fp_geo, gdf]
     for col in collections:
         gc2 = GeoSpatialCollection(col, shapetype)
         gi2 = [i.flopy_geometry.__geo_interface__ for i in gc2]
@@ -439,8 +441,9 @@ def test_linestring_collection(linestring, multilinestring):
     points = gc1.points
     geojson = gc1.geojson
     fp_geo = gc1.flopy_geometry
+    gdf = gc1.geo_dataframe
 
-    collections = [shp, shply, points, geojson, fp_geo]
+    collections = [shp, shply, points, geojson, fp_geo, gdf]
     for col in collections:
         gc2 = GeoSpatialCollection(col, shapetype)
         gi2 = [i.flopy_geometry.__geo_interface__ for i in gc2]
@@ -485,8 +488,9 @@ def test_mixed_collection(
     points = gc1.points
     geojson = gc1.geojson
     fp_geo = gc1.flopy_geometry
+    gdf = gc1.geo_dataframe
 
-    collections = [shp, shply, lshply, points, geojson, fp_geo]
+    collections = [shp, shply, lshply, points, geojson, fp_geo, gdf]
     for col in collections:
         gc2 = GeoSpatialCollection(col, shapetype)
 

autotest/test_grid.py

@@ -4,6 +4,7 @@
 from contextlib import nullcontext
 from warnings import warn
 
+import geopandas
 import matplotlib
 import numpy as np
 import pytest
@@ -1233,13 +1234,13 @@ def test_vertex_neighbors(vertex_grid):
 
 def test_unstructured_neighbors(unstructured_grid):
     rook_neighbors = unstructured_grid.neighbors(5)
-    assert np.allclose(rook_neighbors, [0, 10, 1, 2, 3, 7, 12])
+    assert np.allclose(rook_neighbors, [0, 10, 1, 6, 11, 2, 7, 12])
 
     queen_neighbors = unstructured_grid.neighbors(
         5, method="queen", reset=True
     )
     assert np.allclose(
-        queen_neighbors, [0, 10, 1, 6, 11, 2, 3, 4, 7, 8, 12, 13]
+        queen_neighbors, [0, 10, 1, 6, 11, 2, 3, 7, 8, 12, 13]
     )
 
 
@@ -1300,3 +1301,66 @@ def test_get_lni_unstructured(grid):
             )
         )
         assert csum[layer] + i == nn
+
+
+def test_structured_convert(structured_grid):
+    factor = 3
+    new_grid = structured_grid.convert_grid(factor=factor)
+
+    xf = np.sum(new_grid.xvertices) / np.sum(structured_grid.xvertices)
+    yf = np.sum(new_grid.yvertices) / np.sum(structured_grid.yvertices)
+    zf = np.sum(new_grid.zvertices) / np.sum(structured_grid.zvertices)
+    if xf != factor or yf != factor or zf != factor:
+        raise AssertionError(
+            "structured grid conversion is not returning proper vertices"
+        )
+
+
+def test_vertex_convert(vertex_grid):
+    factor = 3
+    new_grid = vertex_grid.convert_grid(factor=factor)
+
+    xf = np.sum(new_grid.xvertices) / np.sum(vertex_grid.xvertices)
+    yf = np.sum(new_grid.yvertices) / np.sum(vertex_grid.yvertices)
+    zf = np.sum(new_grid.zvertices) / np.sum(vertex_grid.zvertices)
+    if xf != factor or yf != factor or zf != factor:
+        raise AssertionError(
+            "structured grid conversion is not returning proper vertices"
+        )
+
+
+def test_unstructured_convert(unstructured_grid):
+    factor = 3
+    new_grid = unstructured_grid.convert_grid(factor=factor)
+
+    xf = np.sum(new_grid.xvertices) / np.sum(unstructured_grid.xvertices)
+    yf = np.sum(new_grid.yvertices) / np.sum(unstructured_grid.yvertices)
+    zf = np.sum(new_grid.zvertices) / np.sum(unstructured_grid.zvertices)
+    if xf != factor or yf != factor or zf != factor:
+        raise AssertionError(
+            "structured grid conversion is not returning proper vertices"
+        )
+
+
+@requires_pkg("geopandas")
+def test_geo_dataframe(structured_grid, vertex_grid, unstructured_grid):
+    grids = (
+        structured_grid,
+        vertex_grid,
+        unstructured_grid
+    )
+
+    for grid in grids:
+        gdf = grid.geo_dataframe
+        if not isinstance(gdf, geopandas.GeoDataFrame):
+            raise TypeError("geo_dataframe not returning GeoDataFrame object")
+
+        geoms = gdf.geometry.values
+        for node, geom in enumerate(geoms):
+            coords = geom.exterior.coords
+            cv = grid.get_cell_vertices(node)
+            for coord in coords:
+                if coord not in cv:
+                    raise AssertionError(
+                        f"Cell vertices incorrect for node={node}"
+                    )

autotest/test_grid_cases.py

@@ -11,6 +11,7 @@ def structured_small(self):
         delc = 1.0 * np.ones(nrow, dtype=float)
         delr = 1.0 * np.ones(ncol, dtype=float)
         top = 10.0 * np.ones((nrow, ncol), dtype=float)
+        idomain = np.ones((nlay, nrow, ncol), dtype=int)
         botm = np.zeros((nlay, nrow, ncol), dtype=float)
         botm[0, :, :] = 5.0
         botm[1, :, :] = 0.0
@@ -23,6 +24,7 @@ def structured_small(self):
             delr=delr,
             top=top,
             botm=botm,
+            idomain=idomain
         )
 
     def structured_cbd_small(self):
@@ -81,6 +83,7 @@ def vertex_small(self):
             [3, 1.5, 1.5, 4, 4, 5, 8, 7],
             [4, 0.5, 0.5, 4, 6, 7, 10, 9],
         ]
+        idomain = np.ones((nlay, ncpl), dtype=int)
         top = np.ones(ncpl, dtype=float) * 10.0
         botm = np.zeros((nlay, ncpl), dtype=float)
         botm[0, :] = 5.0
@@ -93,6 +96,7 @@ def vertex_small(self):
             cell2d=cell2d,
             top=top,
             botm=botm,
+            idomain=idomain
         )
 
     def unstructured_small(self):
@@ -112,21 +116,21 @@ def unstructured_small(self):
             [10, 1.0, 0.0],
         ]
         iverts = [
-            [0, 0, 1, 4, 3],
-            [1, 1, 2, 5, 4],
-            [2, 3, 4, 7, 6],
-            [3, 4, 5, 8, 7],
-            [4, 6, 7, 10, 9],
-            [5, 0, 1, 4, 3],
-            [6, 1, 2, 5, 4],
-            [7, 3, 4, 7, 6],
-            [8, 4, 5, 8, 7],
-            [9, 6, 7, 10, 9],
-            [10, 0, 1, 4, 3],
-            [11, 1, 2, 5, 4],
-            [12, 3, 4, 7, 6],
-            [13, 4, 5, 8, 7],
-            [14, 6, 7, 10, 9],
+            [0, 1, 4, 3],
+            [1, 2, 5, 4],
+            [3, 4, 7, 6],
+            [4, 5, 8, 7],
+            [6, 7, 10, 9],
+            [0, 1, 4, 3],
+            [1, 2, 5, 4],
+            [3, 4, 7, 6],
+            [4, 5, 8, 7],
+            [6, 7, 10, 9],
+            [0, 1, 4, 3],
+            [1, 2, 5, 4],
+            [3, 4, 7, 6],
+            [4, 5, 8, 7],
+            [6, 7, 10, 9],
         ]
         xcenters = [
             0.5,
@@ -142,6 +146,7 @@ def unstructured_small(self):
             1.5,
             0.5,
         ]
+        idomain = np.ones((nlay, 5), dtype=int)
         top = np.ones((nlay, 5), dtype=float)
         top[0, :] = 10.0
         top[1, :] = 5.0
@@ -159,6 +164,7 @@ def unstructured_small(self):
             ncpl=ncpl,
             top=top.flatten(),
             botm=botm.flatten(),
+            idomain=idomain.flatten()
         )
 
     def unstructured_medium(self):

flopy/discretization/grid.py

@@ -618,6 +618,41 @@ def xyzvertices(self):
     def cross_section_vertices(self):
         return self.xyzvertices[0], self.xyzvertices[1]
 
+    def geo_dataframe(self, polys):
+        """
+        Method returns a geopandas GeoDataFrame of the Grid
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        from ..utils.geospatial_utils import GeoSpatialCollection
+
+        gc = GeoSpatialCollection(
+            polys, shapetype=["Polygon" for _ in range(len(polys))]
+        )
+        gdf = gc.geo_dataframe
+        if self.crs is not None:
+            gdf = gdf.set_crs(crs=self.crs)
+
+        return gdf
+
+    def convert_grid(self, factor):
+        """
+        Method to scale the model grid based on user supplied scale factors
+
+        Parameters
+        ----------
+        factor
+
+        Returns
+        -------
+            Grid object
+        """
+        raise NotImplementedError(
+            "convert_grid must be defined in the child class"
+        )
+
     def _set_neighbors(self, reset=False, method="rook"):
         """
         Method to calculate neighbors via shared edges or shared vertices

flopy/discretization/structuredgrid.py

@@ -774,6 +774,44 @@ def map_polygons(self):
 
         return self._polygons
 
+    @property
+    def geo_dataframe(self):
+        """
+        Returns a geopandas GeoDataFrame of the model grid
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        polys = [[list(zip(*i))] for i in zip(*self.cross_section_vertices)]
+        gdf = super().geo_dataframe(polys)
+        return gdf
+
+    def convert_grid(self, factor):
+        """
+        Method to scale the model grid based on user supplied scale factors
+
+        Parameters
+        ----------
+        factor
+
+        Returns
+        -------
+            Grid object
+        """
+        if super().is_complete:
+            return StructuredGrid(
+                delc=self.delc * factor,
+                delr=self.delr * factor,
+                top=self.top * factor,
+                botm=self.botm * factor,
+                idomain=self.idomain,
+            )
+        else:
+            raise AssertionError(
+                "Grid is not complete and cannot be converted"
+            )
+
     ###############
     ### Methods ###
     ###############

flopy/discretization/unstructuredgrid.py

@@ -553,6 +553,49 @@ def map_polygons(self):
 
         return copy.copy(self._polygons)
 
+    @property
+    def geo_dataframe(self):
+        """
+        Returns a geopandas GeoDataFrame of the model grid
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        polys = [[self.get_cell_vertices(nn)] for nn in range(self.nnodes)]
+        gdf = super().geo_dataframe(polys)
+        return gdf
+
+    def convert_grid(self, factor):
+        """
+        Method to scale the model grid based on user supplied scale factors
+
+        Parameters
+        ----------
+        factor
+
+        Returns
+        -------
+            Grid object
+        """
+        if self.is_complete:
+            return UnstructuredGrid(
+                vertices=[
+                    [i[0], i[1] * factor, i[2] * factor]
+                    for i in self._vertices
+                ],
+                iverts=self._iverts,
+                xcenters=self._xc * factor,
+                ycenters=self._yc * factor,
+                top=self.top * factor,
+                botm=self.botm * factor,
+                idomain=self.idomain,
+            )
+        else:
+            raise AssertionError(
+                "Grid is not complete and cannot be converted"
+            )
+
     def intersect(self, x, y, z=None, local=False, forgive=False):
         """
         Get the CELL2D number of a point with coordinates x and y