diff --git a/autotest/test_geospatial_util.py b/autotest/test_geospatial_util.py
index 3ac9f42a9a..2d0508c1ed 100644
--- a/autotest/test_geospatial_util.py
+++ b/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)
 
diff --git a/autotest/test_grid.py b/autotest/test_grid.py
index 214018777e..523441514f 100644
--- a/autotest/test_grid.py
+++ b/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}"
+                    )
diff --git a/autotest/test_grid_cases.py b/autotest/test_grid_cases.py
index 2721754726..1db81a37e4 100644
--- a/autotest/test_grid_cases.py
+++ b/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):
diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py
index 0fe69caf55..6d845e8cd6 100644
--- a/flopy/discretization/grid.py
+++ b/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
diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py
index c850fa303d..4c943f979f 100644
--- a/flopy/discretization/structuredgrid.py
+++ b/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 ###
     ###############
diff --git a/flopy/discretization/unstructuredgrid.py b/flopy/discretization/unstructuredgrid.py
index 2ebef4f9ff..27d4db69f5 100644
--- a/flopy/discretization/unstructuredgrid.py
+++ b/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
diff --git a/flopy/discretization/vertexgrid.py b/flopy/discretization/vertexgrid.py
index f22aec2960..44641ccb91 100644
--- a/flopy/discretization/vertexgrid.py
+++ b/flopy/discretization/vertexgrid.py
@@ -289,6 +289,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.ncpl)]
+        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 VertexGrid(
+                vertices=[
+                    [i[0], i[1] * factor, i[2] * factor]
+                    for i in self._vertices
+                ],
+                cell2d=[
+                    [i[0], i[1] * factor, i[2] * factor] + i[3:]
+                    for i in self._cell2d
+                ],
+                top=self.top * factor,
+                botm=self.botm * factor,
+            )
+        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
diff --git a/flopy/utils/geospatial_utils.py b/flopy/utils/geospatial_utils.py
index 4cd5620a8b..7e9457f71d 100644
--- a/flopy/utils/geospatial_utils.py
+++ b/flopy/utils/geospatial_utils.py
@@ -263,6 +263,8 @@ def __init__(self, obj, shapetype=None):
         self.__shapefile = import_optional_dependency(
             "shapefile", errors="silent"
         )
+        gpd = import_optional_dependency("geopandas", errors="silent")
+
         shapely_geo = import_optional_dependency(
             "shapely.geometry", errors="silent"
         )
@@ -275,6 +277,7 @@ def __init__(self, obj, shapetype=None):
         self._flopy_geometry = None
         self._points = None
         self.__shapetype = None
+        self.__attributes = None
 
         if isinstance(obj, Collection):
             for shape in obj:
@@ -361,6 +364,16 @@ def __init__(self, obj, shapetype=None):
                 for geom in obj.geoms:
                     self.__collection.append(GeoSpatialUtil(geom))
 
+        if gpd is not None:
+            if isinstance(obj, gpd.GeoDataFrame):
+                self.__attributes = {}
+                for geom in obj.geometry.values:
+                    self.__collection.append(GeoSpatialUtil(geom))
+
+                for k in list(obj):
+                    if k != "geometry":
+                        self.__attributes[k] = obj[k].values
+
         if not self.__collection:
             raise AssertionError(
                 f"Reader is not installed for collection type: {type(obj)}"
@@ -419,6 +432,23 @@ def shapely(self):
         )
         return self._shapely
 
+    @property
+    def geo_dataframe(self):
+        """
+        Property that returns a geopandas DataFrame
+
+        Returns
+        -------
+            geopandas.GeoDataFrame
+        """
+        gpd = import_optional_dependency("geopandas")
+        data = {"geometry": self.shapely.geoms}
+        if self.__attributes is not None:
+            for k, v in self.__attributes.items():
+                data[k] = v
+        gdf = gpd.GeoDataFrame(data)
+        return gdf
+
     @property
     def geojson(self):
         """
diff --git a/flopy/utils/triangle.py b/flopy/utils/triangle.py
index 9d7c5bb7b4..a0d86d18f0 100644
--- a/flopy/utils/triangle.py
+++ b/flopy/utils/triangle.py
@@ -58,7 +58,7 @@ def __init__(
         self.additional_args = additional_args
         self._initialize_vars()
 
-    def add_polygon(self, polygon):
+    def add_polygon(self, polygon, ignore_holes=False):
         """
         Add a polygon
 
@@ -72,6 +72,9 @@ def add_polygon(self, polygon):
             shapely Polygon object
             shapefile Polygon shape
             flopy.utils.geometry.Polygon object
+        ignore_holes : bool
+            method to ignore holes in polygon and only use the exterior
+            coordinates
 
         Returns
         -------
@@ -86,9 +89,10 @@ def add_polygon(self, polygon):
         if polygon[0][0] == polygon[0][-1]:
             polygon[0] = polygon[0][:-1]
         self._polygons.append(polygon[0])
-        if len(polygon) > 1:
-            for hole in polygon[1:]:
-                self.add_hole(hole)
+        if not ignore_holes:
+            if len(polygon) > 1:
+                for hole in polygon[1:]:
+                    self.add_hole(hole)
 
     def add_hole(self, hole):
         """
diff --git a/pyproject.toml b/pyproject.toml
index a7304055fb..0bceeb856b 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -65,6 +65,7 @@ optional = [
     "descartes",
     "fiona",
     "geojson",
+    "geopandas",
     "imageio",
     "netcdf4",
     "pymetis ; platform_system != 'Windows'",