diff --git a/.docs/Notebooks/modelgrid_examples.py b/.docs/Notebooks/modelgrid_examples.py
index 92e4d2b1b9..810290c53f 100644
--- a/.docs/Notebooks/modelgrid_examples.py
+++ b/.docs/Notebooks/modelgrid_examples.py
@@ -937,20 +937,19 @@ def load_iverts(fname):
 plt.colorbar(pc)
 # -
 
-# #### `write_shapefile()`
+# #### `to_geodataframe()`
+#
+# Method to create a geopandas GeoDataFrame of the grid with just the cellid attributes.
 #
-# Method to write a shapefile of the grid with just the cellid attributes. Input parameters include:
 #
-#    - `filename` : shapefile name
-#    - `crs` : either an epsg integer code of model coordinate system, a proj4 string or pyproj CRS instance
-#    - `prjfile` : path to ".prj" projection file that describes the model coordinate system
 
 # +
-# write a shapefile
+# Create a GeoDataFrame and write it to file
 shp_name = os.path.join(gridgen_ws, "freyberg-6_grid.shp")
-epsg = 32614
+modelgrid.crs = 32614
 
-ml.modelgrid.write_shapefile(shp_name, crs=epsg)
+gdf = ml.modelgrid.to_geodataframe()
+gdf.to_file(shp_name)
 # -
 
 try:
diff --git a/.docs/Notebooks/shapefile_export_example.py b/.docs/Notebooks/shapefile_export_example.py
index c542a71a84..b73aec1652 100644
--- a/.docs/Notebooks/shapefile_export_example.py
+++ b/.docs/Notebooks/shapefile_export_example.py
@@ -6,15 +6,25 @@
 #       extension: .py
 #       format_name: light
 #       format_version: '1.5'
-#       jupytext_version: 1.14.5
+#       jupytext_version: 1.17.2
 #   kernelspec:
 #     display_name: Python 3 (ipykernel)
 #     language: python
 #     name: python3
+#   language_info:
+#     codemirror_mode:
+#       name: ipython
+#       version: 3
+#     file_extension: .py
+#     mimetype: text/x-python
+#     name: python
+#     nbconvert_exporter: python
+#     pygments_lexer: ipython3
+#     version: 3.13.5
 #   metadata:
-#     section: export
 #     authors:
-#       - name: Andy Leaf
+#     - name: Andy Leaf
+#     section: export
 # ---
 
 # # Shapefile export demo
@@ -31,6 +41,7 @@
 from pathlib import Path
 from tempfile import TemporaryDirectory
 
+import geopandas as gpd
 import git
 import matplotlib as mpl
 import matplotlib.pyplot as plt
@@ -48,7 +59,8 @@
 # +
 # temporary directory
 temp_dir = TemporaryDirectory()
-outdir = os.path.join(temp_dir.name, "shapefile_export")
+outdir = Path(temp_dir.name) / "shapefile_export"
+outdir.mkdir(exist_ok=True)
 
 # Check if we are in the repository and define the data path.
 
@@ -103,11 +115,12 @@
 
 grid.extent
 
-# ## Declarative export using attached `.export()` methods
+# ## Declarative export using `.to_geodataframe()` method
 # #### Export the whole model to a single shapefile
 
 fname = f"{outdir}/model.shp"
-m.export(fname)
+gdf = m.to_geodataframe(truncate_attrs=True)
+gdf.to_file(fname)
 
 ax = plt.subplot(1, 1, 1, aspect="equal")
 extents = grid.extent
@@ -117,7 +130,8 @@
 ax.set_title(fname)
 
 fname = f"{outdir}/wel.shp"
-m.wel.export(fname)
+gdf = m.wel.to_geodataframe()
+gdf.to_file(fname)
 
 # ### Export a package to a shapefile
 
@@ -126,7 +140,8 @@
 m.lpf.hk
 
 fname = f"{outdir}/hk.shp"
-m.lpf.hk.export(f"{outdir}/hk.shp")
+gdf = m.lpf.hk.to_geodataframe()
+gdf.to_file(fname)
 
 ax = plt.subplot(1, 1, 1, aspect="equal")
 extents = grid.extent
@@ -138,47 +153,31 @@
 
 m.riv.stress_period_data
 
-m.riv.stress_period_data.export(f"{outdir}/riv_spd.shp")
+gdf = m.riv.stress_period_data.to_geodataframe()
+gdf.to_file(f"{outdir}/riv_spd.shp")
 
-# ### MfList.export() exports the whole grid by default, regardless of the locations of the boundary cells
+# ### MfList.to_geodataframe() exports the whole grid by default, regardless of the locations of the boundary cells
 # `sparse=True` only exports the boundary cells in the MfList
 
-m.riv.stress_period_data.export(f"{outdir}/riv_spd.shp", sparse=True)
+gdf = m.riv.stress_period_data.to_geodataframe(sparse=True)
+gdf.to_file(f"{outdir}/riv_spd.shp")
 
-m.wel.stress_period_data.export(f"{outdir}/wel_spd.shp", sparse=True)
+gdf = m.wel.stress_period_data.to_geodataframe(sparse=True)
+gdf.to_file(f"{outdir}/wel_spd.shp")
 
-# ## Ad-hoc exporting using `recarray2shp`
-# * The main idea is to create a recarray with all of the attribute information, and a list of geometry features (one feature per row in the recarray)
-# * each geometry feature is an instance of the `Point`, `LineString` or `Polygon` classes in `flopy.utils.geometry`. The shapefile format requires all the features to be of the same type.
-# * We will use pandas dataframes for these examples because they are easy to work with, and then convert them to recarrays prior to exporting.
-#
-
-from flopy.export.shapefile_utils import recarray2shp
+# ## Ad-hoc exporting using `to_geodataframe()`
 
 # ### combining data from different packages
 # write a shapefile of RIV and WEL package cells
 
-wellspd = pd.DataFrame(m.wel.stress_period_data[0])
-rivspd = pd.DataFrame(m.riv.stress_period_data[0])
-spd = pd.concat([wellspd, rivspd])
-spd.head()
-
-# ##### Create a list of Polygon features from the cell vertices stored in the modelgrid object
-
-# +
-from flopy.utils.geometry import Polygon
-
-vertices = []
-for row, col in zip(spd.i, spd.j):
-    vertices.append(grid.get_cell_vertices(row, col))
-polygons = [Polygon(vrt) for vrt in vertices]
-polygons
-# -
+gdf = m.wel.stress_period_data.to_geodataframe(truncate_attrs=True)
+gdf = m.riv.stress_period_data.to_geodataframe(gpd=gpd, truncate_attrs=True)
+gdf.head()
 
 # ##### write the shapefile
 
 fname = f"{outdir}/bcs.shp"
-recarray2shp(spd.to_records(), geoms=polygons, shpname=fname, crs=grid.epsg)
+gdf.to_file(fname)
 
 ax = plt.subplot(1, 1, 1, aspect="equal")
 extents = grid.extent
@@ -204,11 +203,17 @@
 
 # +
 from flopy.utils.geometry import Point
+from flopy.utils.geospatial_utils import GeoSpatialCollection
 
 geoms = [Point(x, y) for x, y in zip(welldata.x_utm, welldata.y_utm)]
+geoms = GeoSpatialCollection(geoms, "Point")
+gdf = gpd.GeoDataFrame.from_features(geoms)
+
+for col in list(welldata):
+    gdf[col] = welldata[col]
+
+gdf.to_file(f"{outdir}/wel_data.shp")
 
-fname = f"{outdir}/wel_data.shp"
-recarray2shp(welldata.to_records(), geoms=geoms, shpname=fname, crs=grid.epsg)
 # -
 
 ax = plt.subplot(1, 1, 1, aspect="equal")
@@ -235,15 +240,13 @@
 l0 = zip(list(zip(x[:-1], y[:-1])), list(zip(x[1:], y[1:])))
 lines = [LineString(l) for l in l0]
 
-rivdata = pd.DataFrame(m.riv.stress_period_data[0])
-rivdata["reach"] = np.arange(len(lines))
+gdf = gpd.GeoDataFrame(data=m.riv.stress_period_data[0], geometry=lines)
+gdf["reach"] = np.arange(len(lines))
+gdf = gdf.set_crs(epsg=grid.epsg)
+
 lines_shapefile = f"{outdir}/riv_reaches.shp"
-recarray2shp(
-    rivdata.to_records(index=False),
-    geoms=lines,
-    shpname=lines_shapefile,
-    crs=grid.epsg,
-)
+gdf.to_file(lines_shapefile)
+
 # -
 
 ax = plt.subplot(1, 1, 1, aspect="equal")
@@ -253,13 +256,10 @@
 ax.set_ylim(extents[2], extents[3])
 ax.set_title(lines_shapefile)
 
-# #### read in the GIS coverage using `shp2recarray`
-# `shp2recarray` reads a shapefile into a numpy record array, which can easily be converted to a DataFrame
+# #### read in the GIS coverage using geopandas
+# `read_file()` reads a geospatial file into a GeoDataFrame
 
-from flopy.export.shapefile_utils import shp2recarray
-
-linesdata = shp2recarray(lines_shapefile)
-linesdata = pd.DataFrame(linesdata)
+linesdata = gpd.read_file(lines_shapefile)
 linesdata.head()
 
 # ##### Suppose we have some flow information that we read in from the cell budget file
@@ -272,13 +272,7 @@
 
 linesdata["qreach"] = q
 linesdata["qstream"] = np.cumsum(q)
-
-recarray2shp(
-    linesdata.drop("geometry", axis=1).to_records(),
-    geoms=linesdata.geometry.values,
-    shpname=lines_shapefile,
-    crs=grid.epsg,
-)
+linesdata.to_file(lines_shapefile)
 
 ax = plt.subplot(1, 1, 1, aspect="equal")
 extents = grid.extent
@@ -289,7 +283,7 @@
 
 # ## Overriding the model's modelgrid with a user supplied modelgrid
 #
-# In some cases it may be necessary to override the model's modelgrid instance with a seperate modelgrid. An example of this is if the model discretization is in feet and the user would like it projected in meters. Exporting can be accomplished by supplying a modelgrid as a `kwarg` in any of the `export()` methods within flopy. Below is an example:
+# In some cases it may be necessary to override the model's modelgrid instance with a seperate modelgrid. An example of this is if the model discretization is in feet and the user would like it projected in meters. Exporting can be accomplished by supplying a GeoDataFrame in the `to_geodataframe()` methods within flopy. Below is an example:
 
 # +
 mg0 = m.modelgrid
@@ -306,13 +300,18 @@
 )
 
 # exporting an entire model
-m.export(f"{outdir}/freyberg.shp", modelgrid=modelgrid)
+gdf = modelgrid.to_geodataframe()
+
+gdf = m.to_geodataframe(gdf=gdf)
+gdf.to_file(f"{outdir}/freyberg.shp")
 # -
 
 # And for a specific parameter the method is the same
 
 fname = f"{outdir}/hk.shp"
-m.lpf.hk.export(fname, modelgrid=modelgrid)
+gdf = modelgrid.to_geodataframe()
+gdfhk = m.lpf.hk.to_geodataframe(gdf=gdf)
+gdfhk.to_file(fname)
 
 ax = plt.subplot(1, 1, 1, aspect="equal")
 extents = modelgrid.extent
diff --git a/.docs/Notebooks/shapefile_feature_examples.py b/.docs/Notebooks/shapefile_feature_examples.py
index 84a9050cc6..4af64a3702 100644
--- a/.docs/Notebooks/shapefile_feature_examples.py
+++ b/.docs/Notebooks/shapefile_feature_examples.py
@@ -6,23 +6,31 @@
 #       extension: .py
 #       format_name: light
 #       format_version: '1.5'
-#       jupytext_version: 1.14.5
+#       jupytext_version: 1.17.2
 #   kernelspec:
 #     display_name: Python 3 (ipykernel)
 #     language: python
 #     name: python3
+#   language_info:
+#     codemirror_mode:
+#       name: ipython
+#       version: 3
+#     file_extension: .py
+#     mimetype: text/x-python
+#     name: python
+#     nbconvert_exporter: python
+#     pygments_lexer: ipython3
+#     version: 3.13.5
 #   metadata:
-#     section: flopy
 #     authors:
-#       - name: Andy Leaf
+#     - name: Andy Leaf
+#     section: flopy
 # ---
 
 # # Working with shapefiles
 #
 # This notebook shows some lower-level functionality in `flopy` for working with shapefiles
 # including:
-# * `recarray2shp` convience function for writing a numpy record array to a shapefile
-# * `shp2recarray` convience function for quickly reading a shapefile into a numpy recarray
 # * `utils.geometry` classes for writing shapefiles of model input/output. For example, quickly writing a shapefile of model cells with errors identified by the checker
 # * examples of how the `Point` and `LineString` classes can be used to quickly plot pathlines and endpoints from MODPATH (these are also used by the `PathlineFile` and `EndpointFile` classes to write shapefiles of this output)
 
@@ -34,6 +42,7 @@
 from pathlib import Path
 from tempfile import TemporaryDirectory
 
+import geopandas as gpd
 import git
 import matplotlib as mpl
 import matplotlib.pyplot as plt
@@ -41,7 +50,6 @@
 import pooch
 
 import flopy
-from flopy.export.shapefile_utils import recarray2shp, shp2recarray
 from flopy.utils import geometry
 from flopy.utils.geometry import LineString, Point, Polygon
 from flopy.utils.modpathfile import EndpointFile, PathlineFile
@@ -104,32 +112,24 @@
 
 geoms[0].plot()  # this feature requires descartes
 
-# ### write the shapefile
-# * the projection (.prj) file can be written using an epsg code
-# * or copied from an existing .prj file
+# ### create a GeoDataFrame and write it to shapefile
+# * the projection (.prj) file can be written using an epsg code or crs
 
 # +
 from pathlib import Path
 
-recarray2shp(chk.summary_array, geoms, os.path.join(workspace, "test.shp"), crs=26715)
-shape_path = os.path.join(workspace, "test.prj")
+gdf = gpd.GeoDataFrame(data=chk.summary_array, geometry=geoms)
+gdf.to_file(os.path.join(workspace, "test.shp"))
 
-shutil.copy(shape_path, os.path.join(workspace, "26715.prj"))
-recarray2shp(
-    chk.summary_array,
-    geoms,
-    os.path.join(workspace, "test.shp"),
-    prjfile=os.path.join(workspace, "26715.prj"),
-)
 # -
 
 # ### read it back in
-# * flopy geometry objects representing the shapes are stored in the 'geometry' field
+# * We can use geopandas to read it back into memory
 
-ra = shp2recarray(os.path.join(workspace, "test.shp"))
-ra
+gdf = gpd.read_file(os.path.join(workspace, "test.shp"))
+gdf.head()
 
-ra.geometry[0].plot()
+gdf.geometry.plot()
 # -
 
 # ## Other geometry types
diff --git a/autotest/test_binaryfile.py b/autotest/test_binaryfile.py
index e6da83e5a6..59558d98a2 100644
--- a/autotest/test_binaryfile.py
+++ b/autotest/test_binaryfile.py
@@ -275,7 +275,7 @@ def test_load_binary_head_file(example_data_path):
 
 def test_plot_binary_head_file(example_data_path):
     hf = HeadFile(example_data_path / "freyberg" / "freyberg.githds")
-    hf.mg.set_coord_info(xoff=1000.0, yoff=200.0, angrot=15.0)
+    hf.modelgrid.set_coord_info(xoff=1000.0, yoff=200.0, angrot=15.0)
 
     assert isinstance(hf.plot(), Axes)
     plt.close()
diff --git a/autotest/test_copy.py b/autotest/test_copy.py
index e466d5007f..3a99ff16f3 100644
--- a/autotest/test_copy.py
+++ b/autotest/test_copy.py
@@ -39,7 +39,11 @@ def model_is_copy(m1, m2):
             return False
     for k, v in m1.__dict__.items():
         v2 = m2.__dict__[k]
-        if v2 is v and type(v) not in [bool, str, type(None), float, int]:
+        # Allow identity sharing for immutable types including NumPy scalars
+        is_immutable = type(v) in [bool, str, type(None), float, int] or isinstance(
+            v, np.generic
+        )
+        if v2 is v and not is_immutable:
             # some mf6 objects aren't copied with deepcopy
             if isinstance(v, MFSimulationData):
                 continue
@@ -78,7 +82,15 @@ def package_is_copy(pk1, pk2):
     """
     for k, v in pk1.__dict__.items():
         v2 = pk2.__dict__[k]
-        if v2 is v and type(v) not in [bool, str, type(None), float, int, tuple]:
+        is_immutable = type(v) in [
+            bool,
+            str,
+            type(None),
+            float,
+            int,
+            tuple,
+        ] or isinstance(v, np.generic)
+        if v2 is v and not is_immutable:
             # Deep copy doesn't work for ModflowUtltas
             if not inspect.isclass(v):
                 return False
@@ -139,11 +151,10 @@ def package_is_copy(pk1, pk2):
                     if not isinstance(a1, np.ndarray):
                         if a1 != a2:
                             return False
-                    # TODO: this may return False if there are nans
-                    elif not np.allclose(v.array, v2.array):
+                    elif not np.allclose(v.array, v2.array, equal_nan=True):
                         return False
         elif isinstance(v, np.ndarray):
-            if not np.allclose(v, v2):
+            if not np.allclose(v, v2, equal_nan=True):
                 return False
         elif v != v2:
             return False
diff --git a/autotest/test_export.py b/autotest/test_export.py
index f6c0506332..525f94c0ed 100644
--- a/autotest/test_export.py
+++ b/autotest/test_export.py
@@ -239,32 +239,6 @@ def test_freyberg_export(function_tmpdir, example_data_path):
     assert m.drn.stress_period_data.mg.yoffset == m.modelgrid.yoffset
     assert m.drn.stress_period_data.mg.angrot == m.modelgrid.angrot
 
-    # get wkt text from pyproj
-    wkt = m.modelgrid.crs.to_wkt()
-
-    # if wkt text was fetched from pyproj
-    if wkt is not None:
-        # test default package export
-        shape = function_tmpdir / f"{name}_dis.shp"
-        m.dis.export(shape)
-        for suffix in [".dbf", ".prj", ".shp", ".shx"]:
-            part = shape.with_suffix(suffix)
-            assert part.exists()
-            if suffix == ".prj":
-                assert part.read_text() == wkt
-            part.unlink()
-
-        # test default package export to higher level dir ?
-
-        # test sparse package export
-        shape = function_tmpdir / f"{name}_drn_sparse.shp"
-        m.drn.stress_period_data.export(shape, sparse=True)
-        for suffix in [".dbf", ".prj", ".shp", ".shx"]:
-            part = shape.with_suffix(suffix)
-            assert part.exists()
-            if suffix == ".prj":
-                assert part.read_text() == wkt
-
 
 @requires_pkg("pyshp", name_map={"pyshp": "shapefile"})
 @pytest.mark.parametrize("missing_arrays", [True, False])
@@ -319,7 +293,7 @@ def test_export_output(crs, function_tmpdir, example_data_path):
     assert read_crs == get_authority_crs(4326)
 
 
-@requires_pkg("pyshp", name_map={"pyshp": "shapefile"})
+@requires_pkg("pyshp", "geopandas", name_map={"pyshp": "shapefile"})
 def test_write_gridlines_shapefile(function_tmpdir):
     import shapefile
 
@@ -334,7 +308,10 @@ def test_write_gridlines_shapefile(function_tmpdir):
         crs=26715,
     )
     outshp = function_tmpdir / "gridlines.shp"
-    write_gridlines_shapefile(outshp, sg)
+    gdf = sg.grid_line_geodataframe()
+    gdf.to_file(outshp)
+
+    # write_gridlines_shapefile(outshp, sg)
 
     for suffix in [".dbf", ".shp", ".shx"]:
         assert outshp.with_suffix(suffix).exists()
@@ -373,7 +350,9 @@ def test_export_shapefile_polygon_closed(function_tmpdir):
     shp.close()
 
 
-@requires_pkg("rasterio", "pyshp", "scipy", name_map={"pyshp": "shapefile"})
+@requires_pkg(
+    "rasterio", "pyshp", "scipy", "geopandas", name_map={"pyshp": "shapefile"}
+)
 def test_export_array(function_tmpdir, example_data_path):
     import rasterio
     from scipy.ndimage import rotate
@@ -396,7 +375,8 @@ def test_export_array(function_tmpdir, example_data_path):
     )
     arr = np.loadtxt(function_tmpdir / "fb.asc", skiprows=6)
 
-    m.modelgrid.write_shapefile(function_tmpdir / "grid.shp")
+    gdf = m.modelgrid.to_geodataframe()
+    gdf.to_file(function_tmpdir / "grid.shp")
 
     # check bounds
     with open(function_tmpdir / "fb.asc") as src:
@@ -507,7 +487,7 @@ def test_shapefile(function_tmpdir, namfile):
     )
 
 
-@requires_pkg("pyshp", name_map={"pyshp": "shapefile"})
+@requires_pkg("pyshp", "geopandas", name_map={"pyshp": "shapefile"})
 @pytest.mark.slow
 @pytest.mark.parametrize("namfile", namfiles())
 def test_shapefile_export_modelgrid_override(function_tmpdir, namfile):
@@ -820,7 +800,7 @@ def test_export_contours(function_tmpdir, example_data_path):
 
 
 @pytest.mark.mf6
-@requires_pkg("pyshp", "shapely", name_map={"pyshp": "shapefile"})
+@requires_pkg("pyshp", "shapely", "geopandas", name_map={"pyshp": "shapefile"})
 def test_export_mf6_shp(function_tmpdir):
     from shapefile import Reader
 
@@ -909,12 +889,15 @@ def test_export_mf6_shp(function_tmpdir):
         )
         pass
 
-    m.export(function_tmpdir / "mfnwt.shp")
-    gwf.export(function_tmpdir / "mf6.shp")
+    gdf = m.to_geodataframe()
+    gdf.to_file(function_tmpdir / "mfnwt.shp")
 
-    # check that the shapefiles are the same
-    ra = shp2recarray(function_tmpdir / "mfnwt.shp")
-    ra6 = shp2recarray(function_tmpdir / "mf6.shp")
+    gdf2 = gwf.to_geodataframe()
+    gdf2.to_file(function_tmpdir / "mf6.shp")
+
+    # check that the geodataframes are the same
+    ra = gdf.to_records()
+    ra6 = gdf2.to_records()
 
     # check first and last exported cells
     assert ra.geometry[0] == ra6.geometry[0]
@@ -922,11 +905,12 @@ def test_export_mf6_shp(function_tmpdir):
     # fields
     different_fields = list(set(ra.dtype.names).difference(ra6.dtype.names))
     different_fields = [
-        f for f in different_fields if "thick" not in f and "rech" not in f
+        f
+        for f in different_fields
+        if "thick" not in f and "rech" not in f and "model" not in f
     ]
     assert len(different_fields) == 0
-    for lay in np.arange(m.nlay) + 1:
-        assert np.sum(np.abs(ra[f"rech_{lay}"] - ra6[f"rechar{lay}"])) < 1e-6
+
     common_fields = set(ra.dtype.names).intersection(ra6.dtype.names)
     common_fields.remove("geometry")
     # array values
@@ -938,15 +922,13 @@ def test_export_mf6_shp(function_tmpdir):
                 assert np.abs(it - it6) < 1e-6
 
     # Compare exported riv shapefiles
-    riv.export(function_tmpdir / "riv.shp")
-    riv6.export(function_tmpdir / "riv6.shp")
-    with (
-        Reader(function_tmpdir / "riv.shp") as riv_shp,
-        Reader(function_tmpdir / "riv6.shp") as riv6_shp,
-    ):
-        assert list(riv_shp.shapeRecord(-1).record) == list(
-            riv6_shp.shapeRecord(-1).record
-        )
+    rdf = riv.to_geodataframe()
+    rdf.to_file(function_tmpdir / "riv.shp")
+
+    rdf6 = riv6.to_geodataframe()
+    rdf6.to_file(function_tmpdir / "riv6.shp")
+
+    assert rdf6.geometry.values[-1] == rdf.geometry.values[-1]
 
     # Check wel export with timeseries
     wel_spd_0 = flopy.mf6.ModflowGwfwel.stress_period_data.empty(
@@ -958,10 +940,11 @@ def test_export_mf6_shp(function_tmpdir):
         maxbound=1,
         stress_period_data={0: wel_spd_0[0]},
     )
-    wel.export(function_tmpdir / "wel_test.shp")
+    gdf = wel.to_geodataframe()
+    gdf.to_file(function_tmpdir / "wel_test.shp")
 
 
-@requires_pkg("pyshp", name_map={"pyshp": "shapefile"})
+@requires_pkg("geopandas")
 @pytest.mark.slow
 def test_export_huge_shapefile(function_tmpdir):
     nlay = 2
@@ -987,8 +970,8 @@ def test_export_huge_shapefile(function_tmpdir):
         top=top,
         botm=botm,
     )
-
-    m.export(function_tmpdir / "huge.shp")
+    gdf = m.to_geodataframe()
+    gdf.to_file(function_tmpdir / "huge.shp")
 
 
 @requires_pkg("netCDF4", "pyproj")
@@ -2053,7 +2036,7 @@ def test_vtk_export_disu2_grid(function_tmpdir, example_data_path):
 
 @pytest.mark.mf6
 @requires_exe("mf6", "gridgen")
-@requires_pkg("vtk", "pyshp", "shapely", name_map={"pyshp": "shapefile"})
+@requires_pkg("vtk", "pyshp", "shapely", "geopandas", name_map={"pyshp": "shapefile"})
 def test_vtk_export_disu_model(function_tmpdir):
     from vtkmodules.util.numpy_support import vtk_to_numpy
 
@@ -2177,7 +2160,7 @@ def test_mf6_chd_shapefile_export_structured(function_tmpdir, use_pandas, sparse
             assert set(chd_vals) == {1.0, 2.0}
 
 
-@requires_pkg("pyshp", name_map={"pyshp": "shapefile"})
+@requires_pkg("pyshp", "geopandas", name_map={"pyshp": "shapefile"})
 @pytest.mark.parametrize("use_pandas", [True])  # TODO: test non-pandas
 @pytest.mark.parametrize("sparse", [True])  # TODO: test non-sparse
 def test_mf6_chd_shapefile_export_unstructured(function_tmpdir, use_pandas, sparse):
@@ -2264,7 +2247,7 @@ def disv_sim(name, tmpdir):
     return sim
 
 
-@requires_pkg("pyshp", name_map={"pyshp": "shapefile"})
+@requires_pkg("geopandas")
 @pytest.mark.parametrize("use_pandas", [True])  # TODO: test non-pandas
 @pytest.mark.parametrize("sparse", [True])  # TODO: test non-sparse
 def test_mf6_chd_shapefile_export_vertex(function_tmpdir, use_pandas, sparse):
@@ -2286,39 +2269,37 @@ def test_mf6_chd_shapefile_export_vertex(function_tmpdir, use_pandas, sparse):
     chd = ModflowGwfchd(gwf, stress_period_data=chd_cells)
 
     shpfile = function_tmpdir / f"chd_disv_{use_pandas}_{sparse}.shp"
-    gwf.chd.stress_period_data.export(shpfile, sparse=sparse)
+    gdf = gwf.chd.stress_period_data.to_geodataframe(sparse=sparse)
+    gdf.to_file(shpfile)
 
     # Check that shapefile and sidecar files exist
     for ext in [".shp", ".shx", ".dbf"]:
         assert shpfile.with_suffix(ext).exists(), f"{shpfile.with_suffix(ext)} missing"
 
-    # Read shapefile and check records
-    with shapefile.Reader(str(shpfile)) as sf:
-        records = list(sf.records())
-        fields = [f[0] for f in sf.fields[1:]]  # skip DeletionFlag
+    if sparse:
+        # Only CHD cells should be exported
+        assert len(gdf) == len(chd_cells)
+    else:
+        # All grid cells should be exported
+        assert len(gdf) == gwf.modelgrid.nnodes
 
-        if sparse:
-            # Only CHD cells should be exported
-            assert len(records) == len(chd_cells)
-        else:
-            # All grid cells should be exported
-            assert len(records) == gwf.modelgrid.nnodes
+    # Check attribute values for CHD cells
+    cnt = 0
+    chd_vals = []
+    for _, row in gdf.iterrows():
+        col = f"chd_head_{cnt}_0"
+        if not np.isnan(row[col]):
+            chd_vals.append(row[col])
+            cnt += 1
 
-        # Check attribute values for CHD cells
-        chd_vals = [
-            rec[fields.index(next(f for f in fields if "head" in f))] for rec in records
-        ]
-        if sparse:
-            # Should match the CHD values
-            assert set(chd_vals) == {1.0, 2.0, 3.0}
+    if sparse:
+        # Should match the CHD values
+        assert set(chd_vals) == {1.0, 2.0, 3.0}
 
 
-@requires_pkg("pyshp", name_map={"pyshp": "shapefile"})
+@requires_pkg("geopandas")
 def test_issue_2628(function_tmpdir):
-    import shapefile
-
     from flopy.discretization import StructuredGrid
-    from flopy.export.shapefile_utils import write_grid_shapefile
 
     nrow, ncol = 3, 3
     idomain = np.ones((nrow, ncol), dtype=int)
@@ -2339,23 +2320,20 @@ def test_issue_2628(function_tmpdir):
     props = {"idomain": idomain, "sy": sy, "thickness": thickness, "k": k}
     prop_names = list(props.keys())
 
+    gdf = grid.to_geodataframe()
+    for prop, array in props.items():
+        gdf[prop] = array.ravel()
+
     outshp = function_tmpdir / "test_dict_handling.shp"
-    write_grid_shapefile(outshp, grid, props)
+    gdf.to_file(outshp)
 
     for suffix in [".dbf", ".shp", ".shx"]:
         assert outshp.with_suffix(suffix).exists()
 
-    with shapefile.Reader(str(outshp)) as sf:
-        records = list(sf.records())
-        fields = [f[0] for f in sf.fields[1:]]  # skip DeletionFlag
-
-        assert len(records) == nrow * ncol
-        assert [f for f in fields if f in prop_names] == prop_names
+    assert len(gdf) == nrow * ncol
+    assert [f for f in list(gdf) if f in prop_names] == prop_names
 
-        for record in records:
-            row = record[fields.index("row")] - 1
-            col = record[fields.index("column")] - 1
-            assert record[fields.index("idomain")] == idomain[row, col]
-            assert np.isclose(record[fields.index("sy")], sy[row, col])
-            assert np.isclose(record[fields.index("thickness")], thickness[row, col])
-            assert np.isclose(record[fields.index("k")], k[row, col])
+    assert np.allclose(gdf.idomain.values, idomain.ravel())
+    assert np.allclose(gdf.sy.values, sy.ravel())
+    assert np.allclose(gdf.thickness.values, thickness.ravel())
+    assert np.allclose(gdf.k.values, k.ravel())
diff --git a/autotest/test_formattedfile.py b/autotest/test_formattedfile.py
index 0d31ed308c..26ac20f451 100644
--- a/autotest/test_formattedfile.py
+++ b/autotest/test_formattedfile.py
@@ -69,7 +69,7 @@ def test_headfile_build_index(example_data_path):
 def test_formattedfile_reference(example_data_path):
     h = FormattedHeadFile(example_data_path / "mf2005_test" / "test1tr.githds")
     assert isinstance(h, FormattedHeadFile)
-    h.mg.set_coord_info(xoff=1000.0, yoff=200.0, angrot=15.0)
+    h.modelgrid.set_coord_info(xoff=1000.0, yoff=200.0, angrot=15.0)
 
     assert isinstance(h.plot(masked_values=[6999.000]), Axes)
     plt.close()
diff --git a/autotest/test_geospatial_util.py b/autotest/test_geospatial_util.py
index decbcf3706..135d8b6957 100644
--- a/autotest/test_geospatial_util.py
+++ b/autotest/test_geospatial_util.py
@@ -374,7 +374,7 @@ def test_polygon_collection(polygon, poly_w_hole, multipolygon):
     points = gc1.points
     geojson = gc1.geojson
     fp_geo = gc1.flopy_geometry
-    gdf = gc1.geo_dataframe
+    gdf = gc1.geodataframe
 
     collections = [shp, shply, points, geojson, fp_geo, gdf]
     for col in collections:
@@ -408,7 +408,7 @@ def test_point_collection(point, multipoint):
     points = gc1.points
     geojson = gc1.geojson
     fp_geo = gc1.flopy_geometry
-    gdf = gc1.geo_dataframe
+    gdf = gc1.geodataframe
 
     collections = [shp, shply, points, geojson, fp_geo, gdf]
     for col in collections:
@@ -436,7 +436,7 @@ def test_linestring_collection(linestring, multilinestring):
     points = gc1.points
     geojson = gc1.geojson
     fp_geo = gc1.flopy_geometry
-    gdf = gc1.geo_dataframe
+    gdf = gc1.geodataframe
 
     collections = [shp, shply, points, geojson, fp_geo, gdf]
     for col in collections:
@@ -481,7 +481,7 @@ def test_mixed_collection(
     points = gc1.points
     geojson = gc1.geojson
     fp_geo = gc1.flopy_geometry
-    gdf = gc1.geo_dataframe
+    gdf = gc1.geodataframe
 
     collections = [shp, shply, lshply, points, geojson, fp_geo, gdf]
     for col in collections:
@@ -525,7 +525,7 @@ def test_geopandas_dtypes(
     col = Collection(col)
 
     gc1 = GeoSpatialCollection(col)
-    gdf = gc1.geo_dataframe
+    gdf = gc1.geodataframe
 
     collections = [gdf, gdf.geometry, gdf.geometry.values]
     for col in collections:
diff --git a/autotest/test_grid.py b/autotest/test_grid.py
index c390f9b151..82d32ba069 100644
--- a/autotest/test_grid.py
+++ b/autotest/test_grid.py
@@ -1587,14 +1587,14 @@ def test_unstructured_convert(unstructured_grid):
 
 
 @requires_pkg("geopandas")
-def test_geo_dataframe(structured_grid, vertex_grid, unstructured_grid):
+def test_geodataframe(structured_grid, vertex_grid, unstructured_grid):
     geopandas = import_optional_dependency("geopandas")
     grids = (structured_grid, vertex_grid, unstructured_grid)
 
     for grid in grids:
-        gdf = grid.geo_dataframe
+        gdf = grid.to_geodataframe()
         if not isinstance(gdf, geopandas.GeoDataFrame):
-            raise TypeError("geo_dataframe not returning GeoDataFrame object")
+            raise TypeError("geodataframe not returning GeoDataFrame object")
 
         geoms = gdf.geometry.values
         for node, geom in enumerate(geoms):
diff --git a/autotest/test_gridgen.py b/autotest/test_gridgen.py
index 65a80b3412..d369a21cd9 100644
--- a/autotest/test_gridgen.py
+++ b/autotest/test_gridgen.py
@@ -134,7 +134,7 @@ def test_add_refinement_feature(function_tmpdir, grid_type):
 
 @pytest.mark.slow
 @requires_exe("mf6", "gridgen")
-@requires_pkg("shapely")
+@requires_pkg("shapely", "geopandas")
 def test_mf6disv(function_tmpdir):
     from shapely.geometry import Polygon
 
@@ -202,9 +202,12 @@ def test_mf6disv(function_tmpdir):
 
     # write grid and model shapefiles
     fname = function_tmpdir / "grid.shp"
-    gwf.modelgrid.write_shapefile(fname)
+    gdf = gwf.modelgrid.to_geodataframe()
+    gdf.to_file(fname)
+
     fname = function_tmpdir / "model.shp"
-    gwf.export(fname)
+    gdf = gwf.to_geodataframe()
+    gdf.to_file(fname)
 
     sim.run_simulation(silent=True)
     head = gwf.output.head().get_data()
@@ -348,11 +351,16 @@ def test_mf6disu(sim_disu_diff_layers):
 
     # write grid and model shapefiles
     fname = ws / "grid.shp"
-    gwf.modelgrid.write_shapefile(fname)
+    gdf = gwf.modelgrid.to_geodataframe()
+    gdf.to_file(fname)
+
     fname = ws / "model.shp"
-    gwf.export(fname)
+    gdf = gwf.to_geodataframe()
+    gdf.to_file(fname)
+
     fname = ws / "chd.shp"
-    gwf.chd.export(fname)
+    gdf = gwf.chd.to_geodataframe()
+    gdf.to_file(fname)
 
     sim.run_simulation(silent=True)
     head = gwf.output.head().get_data()
@@ -561,10 +569,17 @@ def test_mfusg(function_tmpdir):
     ), msg
 
     # test disu, bas6, lpf shapefile export for mfusg unstructured models
-    m.disu.export(function_tmpdir / f"{name}_disu.shp")
-    m.bas6.export(function_tmpdir / f"{name}_bas6.shp")
-    m.lpf.export(function_tmpdir / f"{name}_lpf.shp")
-    m.export(function_tmpdir / f"{name}.shp")
+    gdf = m.disu.to_geodataframe()
+    gdf.to_file(function_tmpdir / f"{name}_disu.shp")
+
+    gdf = m.bas6.to_geodataframe()
+    gdf.to_file(function_tmpdir / f"{name}_bas6.shp")
+
+    gdf = m.lpf.to_geodataframe()
+    gdf.to_file(function_tmpdir / f"{name}_lpf.shp")
+
+    gdf = m.to_geodataframe()
+    gdf.to_file(function_tmpdir / f"{name}.shp")
 
 
 @pytest.mark.slow
diff --git a/autotest/test_modpathfile.py b/autotest/test_modpathfile.py
index 652ad4f8b5..883d11f1b7 100644
--- a/autotest/test_modpathfile.py
+++ b/autotest/test_modpathfile.py
@@ -308,9 +308,9 @@ def test_get_destination_endpoint_data(
 
 @pytest.mark.parametrize("longfieldname", [True, False])
 @requires_exe("mf6", "mp7")
-@requires_pkg("pyshp", "shapely", name_map={"pyshp": "shapefile"})
+@requires_pkg("geopandas", "shapely")
 def test_write_shapefile(function_tmpdir, mp7_small, longfieldname):
-    from shapefile import Reader
+    import geopandas as gpd
 
     # setup and run model, then copy outputs to function_tmpdir
     sim, forward_model_name, _, _, _ = mp7_small
@@ -330,13 +330,7 @@ def test_write_shapefile(function_tmpdir, mp7_small, longfieldname):
     # define shapefile path
     shp_file = ws / "pathlines.shp"
 
-    # add a column to the pathline recarray
-    fieldname = "newfield" + ("longname" if longfieldname else "")
-    fieldval = "x"
-    pathlines = [
-        rfn.append_fields(pl, fieldname, list(repeat(fieldval, len(pl))), dtypes="|U1")
-        for pl in pathlines
-    ]
+    pline_names = [name[:10] for name in pathlines[0].dtype.names]
 
     # write the pathline recarray to shapefile
     pathline_file.write_shapefile(
@@ -350,8 +344,7 @@ def test_write_shapefile(function_tmpdir, mp7_small, longfieldname):
     assert shp_file.is_file()
 
     # load shapefile
-    with Reader(shp_file) as reader:
-        fieldnames = [f[0] for f in reader.fields[1:]]
-        fieldname = "newfiname_" if longfieldname else fieldname
-        assert fieldname in fieldnames
-        assert all(r[fieldname] == fieldval for r in reader.iterRecords())
+    gdf = gpd.read_file(shp_file)
+    fieldnames = list(gdf)
+    for fname in pline_names:
+        assert fname in fieldnames
diff --git a/autotest/test_mp6.py b/autotest/test_mp6.py
index d2d99e0b0a..10aec0a69f 100644
--- a/autotest/test_mp6.py
+++ b/autotest/test_mp6.py
@@ -123,8 +123,10 @@ def test_mpsim(function_tmpdir, mp6_test_path):
         assert stllines[6].strip().split()[-1] == "p2"
 
 
-@requires_pkg("pyshp", "shapely", name_map={"pyshp": "shapefile"})
+@requires_pkg("geopandas", "shapely")
 def test_get_destination_data(function_tmpdir, mp6_test_path):
+    import geopandas as gpd
+
     copy_modpath_files(mp6_test_path, function_tmpdir, "EXAMPLE.")
     copy_modpath_files(mp6_test_path, function_tmpdir, "EXAMPLE-3.")
 
@@ -169,74 +171,59 @@ def test_get_destination_data(function_tmpdir, mp6_test_path):
     )
     assert np.all(np.isin(starting_locs, pathline_locs))
 
-    # test writing a shapefile of endpoints
-    epd.write_shapefile(
-        well_epd,
-        direction="starting",
-        shpname=function_tmpdir / "starting_locs.shp",
-        mg=m.modelgrid,
+    gdf = epd.to_geodataframe(
+        modelgrid=m.modelgrid, data=well_epd, direction="starting"
     )
+    gdf.to_file(function_tmpdir / "starting_locs.shp")
 
     # test writing shapefile of pathlines
-    fpth = function_tmpdir / "pathlines_1per.shp"
-    pthld.write_shapefile(
-        well_pthld,
-        one_per_particle=True,
-        direction="starting",
-        mg=m.modelgrid,
-        shpname=fpth,
-    )
-    fpth = function_tmpdir / "pathlines_1per_end.shp"
-    pthld.write_shapefile(
-        well_pthld,
-        one_per_particle=True,
-        direction="ending",
-        mg=m.modelgrid,
-        shpname=fpth,
+    gdf = pthld.to_geodataframe(
+        m.modelgrid, data=well_pthld, one_per_particle=True, direction="starting"
     )
-    # test writing shapefile of pathlines
-    fpth = function_tmpdir / "pathlines_1per2.shp"
-    pthld.write_shapefile(
-        well_pthld,
-        one_per_particle=True,
-        direction="starting",
-        mg=mg,
-        shpname=fpth,
+    gdf.to_file(function_tmpdir / "pathlines_1per.shp")
+
+    gdf = pthld.to_geodataframe(
+        m.modelgrid, data=well_pthld, one_per_particle=True, direction="ending"
     )
+    gdf.to_file(function_tmpdir / "pathlines_1per_end.shp")
+
     # test writing shapefile of pathlines
-    fpth = function_tmpdir / "pathlines_1per2_ll.shp"
-    pthld.write_shapefile(
-        well_pthld,
-        one_per_particle=True,
-        direction="starting",
-        mg=mg,
-        shpname=fpth,
+    gdf = pthld.to_geodataframe(
+        mg, data=well_pthld, one_per_particle=True, direction="starting"
     )
-    fpth = function_tmpdir / "pathlines.shp"
-    pthld.write_shapefile(
-        well_pthld, one_per_particle=False, mg=m.modelgrid, shpname=fpth
+    gdf.to_file(function_tmpdir / "pathlines_1per2.shp")
+
+    # test writing shapefile of pathlines
+    gdf = pthld.to_geodataframe(
+        mg, data=well_pthld, one_per_particle=True, direction="starting"
     )
+    gdf.to_file(function_tmpdir / "pathlines_1per2_ll.shp")
+
+    # test shapefile construction with one_per_particle=False
+    gdf = pthld.to_geodataframe(m.modelgrid, data=well_pthld, one_per_particle=False)
+    gdf.to_file(function_tmpdir / "pathlines.shp")
 
     # test that endpoints were rotated and written correctly
-    ra = shp2recarray(function_tmpdir / "starting_locs.shp")
-    p3 = ra.geometry[ra.particleid == 4][0]
+    ra = gpd.read_file(function_tmpdir / "starting_locs.shp")
+    p3 = ra[ra.particleid == 4]
+    p3 = p3.geometry.values[0]
     xorig, yorig = m.modelgrid.get_coords(well_epd.x0[0], well_epd.y0[0])
     assert p3.x - xorig + p3.y - yorig < 1e-4
     xorig, yorig = mg1.xcellcenters[3, 4], mg1.ycellcenters[3, 4]
     assert np.abs(p3.x - xorig + p3.y - yorig) < 1e-4  # this also checks for 1-based
 
     # test that particle attribute information is consistent with pathline file
-    ra = shp2recarray(function_tmpdir / "pathlines.shp")
-    inds = (ra.particleid == 8) & (ra.i == 12) & (ra.j == 12)
-    assert ra.time[inds][0] - 20181.7 < 0.1
-    assert ra.xloc[inds][0] - 0.933 < 0.01
+    ra = gpd.read_file(function_tmpdir / "pathlines.shp")
+    tst = ra[(ra.particleid == 8) & (ra.i == 12) & (ra.j == 12)]
+    assert tst.time.values[0] - 20181.7 < 0.1
+    assert tst.xloc.values[0] - 0.933 < 0.01
 
     # test that k, i, j are correct for single geometry pathlines, forwards
     # and backwards
-    ra = shp2recarray(function_tmpdir / "pathlines_1per.shp")
-    assert ra.i[0] == 4, ra.j[0] == 5
-    ra = shp2recarray(function_tmpdir / "pathlines_1per_end.shp")
-    assert ra.i[0] == 13, ra.j[0] == 13
+    ra = gpd.read_file(function_tmpdir / "pathlines_1per.shp")
+    assert ra.i.values[0] == 4, ra.j.values[0] == 5
+    ra = gpd.read_file(function_tmpdir / "pathlines_1per_end.shp")
+    assert ra.i.values[0] == 13, ra.j.values[0] == 13
 
     # test use of arbitrary spatial reference and offset
     mg1.set_coord_info(
@@ -245,21 +232,31 @@ def test_get_destination_data(function_tmpdir, mp6_test_path):
         angrot=mg.angrot,
         crs=mg.epsg,
     )
-    ra = shp2recarray(function_tmpdir / "pathlines_1per2.shp")
-    p3_2 = ra.geometry[ra.particleid == 4][0]
-    test1 = mg1.xcellcenters[3, 4]
-    test2 = mg1.ycellcenters[3, 4]
+    ra = gpd.read_file(function_tmpdir / "pathlines_1per2.shp")
+    p3_2 = ra[ra.particleid == 4]
+    p3_2 = p3_2.geometry.values[0]
     assert (
-        np.abs(p3_2.x[0] - mg1.xcellcenters[3, 4] + p3_2.y[0] - mg1.ycellcenters[3, 4])
+        np.abs(
+            p3_2.xy[0][0]
+            - mg1.xcellcenters[3, 4]
+            + p3_2.xy[1][0]
+            - mg1.ycellcenters[3, 4]
+        )
         < 1e-4
     )
 
     # arbitrary spatial reference with ll specified instead of ul
-    ra = shp2recarray(function_tmpdir / "pathlines_1per2_ll.shp")
-    p3_2 = ra.geometry[ra.particleid == 4][0]
+    ra = gpd.read_file(function_tmpdir / "pathlines_1per2_ll.shp")
+    p3_2 = ra[ra.particleid == 4]
+    p3_2 = p3_2.geometry.values[0]
     mg.set_coord_info(xoff=mg.xoffset, yoff=mg.yoffset, angrot=30.0)
     assert (
-        np.abs(p3_2.x[0] - mg.xcellcenters[3, 4] + p3_2.y[0] - mg.ycellcenters[3, 4])
+        np.abs(
+            p3_2.xy[0][0]
+            - mg.xcellcenters[3, 4]
+            + p3_2.xy[1][0]
+            - mg.ycellcenters[3, 4]
+        )
         < 1e-4
     )
 
@@ -277,10 +274,14 @@ def test_get_destination_data(function_tmpdir, mp6_test_path):
     )
 
     fpth = function_tmpdir / "dis2.shp"
-    m.dis.export(fpth)
+    gdf = m.dis.to_geodataframe()
+    gdf.to_file(fpth)
+
     pthobj = PathlineFile(function_tmpdir / "EXAMPLE-3.pathline")
     fpth = function_tmpdir / "pathlines_1per3.shp"
-    pthobj.write_shapefile(shpname=fpth, direction="ending", mg=mg4)
+
+    gdf = pthobj.to_geodataframe(mg4, direction="ending")
+    gdf.to_file(fpth)
 
 
 def test_loadtxt(function_tmpdir, mp6_test_path):
diff --git a/autotest/test_sfr.py b/autotest/test_sfr.py
index 64ed8bca0f..999342f618 100644
--- a/autotest/test_sfr.py
+++ b/autotest/test_sfr.py
@@ -359,8 +359,10 @@ def test_const(sfr_data):
     assert True
 
 
-@requires_pkg("pyshp", "shapely", name_map={"pyshp": "shapefile"})
+@requires_pkg("geopandas")
 def test_export(function_tmpdir, sfr_data):
+    import geopandas as gpd
+
     m = Modflow()
     dis = ModflowDis(m, 1, 10, 10, lenuni=2, itmuni=4)
 
@@ -373,28 +375,13 @@ def test_export(function_tmpdir, sfr_data):
     sfr.export_outlets(function_tmpdir / "outlets.shp")
     sfr.export_transient_variable(function_tmpdir / "inlets.shp", "flow")
 
-    from flopy.export.shapefile_utils import shp2recarray
-
-    ra = shp2recarray(function_tmpdir / "inlets.shp")
-    assert ra.flow0[0] == 1e4
-    ra = shp2recarray(function_tmpdir / "outlets.shp")
-    assert ra.iseg[0] + ra.ireach[0] == 5
-    ra = shp2recarray(function_tmpdir / "linkages.shp")
-    crds = np.array(list(ra.geometry[2].coords))
+    gdf = gpd.read_file(function_tmpdir / "inlets.shp")
+    assert gdf.flow0.values[0] == 1e4
+    gdf = gpd.read_file(function_tmpdir / "outlets.shp")
+    assert gdf.iseg.values[0] + gdf.ireach.values[0] == 5
+    gdf = gpd.read_file(function_tmpdir / "linkages.shp")
+    crds = np.array(list(gdf.geometry.values[2].coords))
     assert np.array_equal(crds, np.array([[2.5, 4.5], [3.5, 5.5]]))
-    ra = shp2recarray(function_tmpdir / "sfr.shp")
-    assert ra.iseg.sum() == sfr.reach_data.iseg.sum()
-    assert ra.ireach.sum() == sfr.reach_data.ireach.sum()
-    y = np.concatenate([np.array(g.exterior)[:, 1] for g in ra.geometry])
-    x = np.concatenate([np.array(g.exterior)[:, 0] for g in ra.geometry])
-
-    assert (x.min(), x.max(), y.min(), y.max()) == m.modelgrid.extent
-    assert ra[(ra.iseg == 2) & (ra.ireach == 1)]["geometry"][0].bounds == (
-        6.0,
-        2.0,
-        7.0,
-        3.0,
-    )
 
 
 def test_example(mf2005_model_path):
diff --git a/autotest/test_shapefile_utils.py b/autotest/test_shapefile_utils.py
index a81eec6f2e..a96440abc5 100644
--- a/autotest/test_shapefile_utils.py
+++ b/autotest/test_shapefile_utils.py
@@ -14,7 +14,7 @@
 from .test_grid import minimal_unstructured_grid_info, minimal_vertex_grid_info
 
 
-@requires_pkg("pyshp", "shapely", name_map={"pyshp": "shapefile"})
+@requires_pkg("geopandas")
 def test_model_attributes_to_shapefile(example_data_path, function_tmpdir):
     # freyberg mf2005 model
     name = "freyberg"
@@ -22,8 +22,9 @@ def test_model_attributes_to_shapefile(example_data_path, function_tmpdir):
     ws = example_data_path / name
     m = flopy.modflow.Modflow.load(namfile, model_ws=ws, check=False, verbose=False)
     shpfile_path = function_tmpdir / f"{name}.shp"
-    pakg_names = ["DIS", "BAS6", "LPF", "WEL", "RIV", "RCH", "OC", "PCG"]
-    model_attributes_to_shapefile(shpfile_path, m, pakg_names)
+    pkg_names = ["DIS", "BAS6", "LPF", "WEL", "RIV", "RCH", "OC", "PCG"]
+    gdf = m.to_geodataframe(package_names=pkg_names)
+    gdf.to_file(shpfile_path)
     assert shpfile_path.exists()
 
     # freyberg mf6 model
@@ -31,8 +32,9 @@ def test_model_attributes_to_shapefile(example_data_path, function_tmpdir):
     sim = flopy.mf6.MFSimulation.load(sim_name=name, sim_ws=example_data_path / name)
     m = sim.get_model()
     shpfile_path = function_tmpdir / f"{name}.shp"
-    pakg_names = ["dis", "bas6", "npf", "wel", "riv", "rch", "oc", "pcg"]
-    model_attributes_to_shapefile(shpfile_path, m, pakg_names)
+    pkg_names = ["dis", "bas6", "npf", "wel", "riv", "rch", "oc", "pcg"]
+    gdf = m.to_geodataframe(package_names=pkg_names)
+    gdf.to_file(shpfile_path)
     assert shpfile_path.exists()
 
     # model with a DISU grid with no angldegx arrays
@@ -41,43 +43,48 @@ def test_model_attributes_to_shapefile(example_data_path, function_tmpdir):
     sim = disu_sim(name, function_tmpdir, missing_arrays=True)
     m = sim.get_model(name)
     shpfile_path = function_tmpdir / f"{name}.shp"
-    pakg_names = ["dis"]
-    model_attributes_to_shapefile(shpfile_path, m, pakg_names)
+    pkg_names = ["dis"]
+    gdf = m.to_geodataframe(package_names=pkg_names)
+    gdf.to_file(shpfile_path)
     assert shpfile_path.exists()
 
 
-@requires_pkg("pyproj", "pyshp", "shapely", name_map={"pyshp": "shapefile"})
-def test_write_grid_shapefile(
+@requires_pkg("geopandas")
+def test_create_geodataframe(
     minimal_unstructured_grid_info, minimal_vertex_grid_info, function_tmpdir
 ):
-    import pyproj
-
     d = minimal_unstructured_grid_info
     delr = np.ones(10)
     delc = np.ones(10)
     crs = 26916
     shapefilename = function_tmpdir / "grid.shp"
     sg = StructuredGrid(delr=delr, delc=delc, nlay=1, crs=crs)
-    sg.write_shapefile(shapefilename)
-    data = shp2recarray(shapefilename)
+    gdf = sg.to_geodataframe()
+    gdf.to_file(shapefilename)
+
+    data = gdf.to_records()
     # check that row and column appear as integers in recarray
     assert np.issubdtype(data.dtype["row"], np.integer)
-    assert np.issubdtype(data.dtype["column"], np.integer)
+    assert np.issubdtype(data.dtype["col"], np.integer)
     assert len(data) == sg.nnodes
     written_crs = get_shapefile_crs(shapefilename)
     assert written_crs.to_epsg() == crs
 
     usg = UnstructuredGrid(**d, crs=crs)
-    usg.write_shapefile(shapefilename)
-    data = shp2recarray(shapefilename)
+    gdf = usg.to_geodataframe()
+    gdf.to_file(shapefilename)
+
+    data = gdf.to_records()
     assert len(data) == usg.nnodes
     written_crs = get_shapefile_crs(shapefilename)
     assert written_crs.to_epsg() == crs
 
     d = minimal_vertex_grid_info
     vg = VertexGrid(**d, crs=crs)
-    vg.write_shapefile(shapefilename)
-    data = shp2recarray(shapefilename)
+    gdf = vg.to_geodataframe()
+    gdf.to_file(shapefilename)
+
+    data = gdf.to_records()
     assert len(data) == vg.nnodes
     written_crs = get_shapefile_crs(shapefilename)
     assert written_crs.to_epsg() == crs
diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py
index 1d59801b05..0e14b6f330 100644
--- a/flopy/discretization/grid.py
+++ b/flopy/discretization/grid.py
@@ -593,7 +593,7 @@ def xyzvertices(self):
     def cross_section_vertices(self):
         return self.xyzvertices[0], self.xyzvertices[1]
 
-    def geo_dataframe(self, features, featuretype="Polygon"):
+    def to_geodataframe(self, features, featuretype="Polygon"):
         """
         Method returns a geopandas GeoDataFrame of the Grid
 
@@ -606,26 +606,13 @@ def geo_dataframe(self, features, featuretype="Polygon"):
         gc = GeoSpatialCollection(
             features, shapetype=[featuretype for _ in range(len(features))]
         )
-        gdf = gc.geo_dataframe
+        gdf = gc.geodataframe
         gdf["node"] = gdf.index + 1
         if self.crs is not None:
             gdf = gdf.set_crs(crs=self.crs)
 
         return gdf
 
-    @property
-    def grid_line_geo_dataframe(self):
-        """
-        Method to get a GeoDataFrame of grid lines
-
-        Returns
-        -------
-            GeoDataFrame
-        """
-        gdf = self.geo_dataframe(self.grid_lines, featuretype="LineString")
-        gdf = gdf.rename(columns={"node": "number"})
-        return gdf
-
     def convert_grid(self, factor):
         """
         Method to scale the model grid based on user supplied scale factors
diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py
index 0d2df58fdf..e27a678b1f 100644
--- a/flopy/discretization/structuredgrid.py
+++ b/flopy/discretization/structuredgrid.py
@@ -759,8 +759,7 @@ def map_polygons(self):
 
         return self._polygons
 
-    @property
-    def geo_dataframe(self):
+    def to_geodataframe(self):
         """
         Returns a geopandas GeoDataFrame of the model grid
 
@@ -769,11 +768,40 @@ def geo_dataframe(self):
             GeoDataFrame
         """
         polys = [[list(zip(*i))] for i in zip(*self.cross_section_vertices)]
-        gdf = super().geo_dataframe(polys)
+        gdf = super().to_geodataframe(polys)
         gdf["row"] = sorted(list(range(1, self.nrow + 1)) * self.ncol)
         gdf["col"] = list(range(1, self.ncol + 1)) * self.nrow
         return gdf
 
+    def grid_line_geodataframe(self):
+        """
+        Method to get a GeoDataFrame of grid lines
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        gdf = super().to_geodataframe(self.grid_lines, featuretype="LineString")
+        gdf = gdf.rename(columns={"node": "number"})
+        return gdf
+
+    @property
+    def geo_dataframe(self):
+        """
+        Returns a geopandas GeoDataFrame of the model grid
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        import warnings
+
+        warnings.warn(
+            "geo_dataframe has been deprecated, use to_geodataframe() instead",
+            DeprecationWarning,
+        )
+        return self.to_geodataframe()
+
     def convert_grid(self, factor):
         """
         Method to scale the model grid based on user supplied scale factors
diff --git a/flopy/discretization/unstructuredgrid.py b/flopy/discretization/unstructuredgrid.py
index f34764dd6f..508a4d2a31 100644
--- a/flopy/discretization/unstructuredgrid.py
+++ b/flopy/discretization/unstructuredgrid.py
@@ -580,6 +580,23 @@ def map_polygons(self):
 
         return copy.copy(self._polygons)
 
+    def to_geodataframe(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().to_geodataframe(polys)
+        if self.nlay > 1:
+            lays = []
+            for ix, ncpl in enumerate(self.ncpl):
+                lays.extend([ix + 1] * ncpl)
+            gdf["layer"] = lays
+        return gdf
+
     @property
     def geo_dataframe(self):
         """
@@ -589,8 +606,24 @@ def geo_dataframe(self):
         -------
             GeoDataFrame
         """
-        polys = [[self.get_cell_vertices(nn)] for nn in range(self.nnodes)]
-        gdf = super().geo_dataframe(polys)
+        import warnings
+
+        warnings.warn(
+            "geo_dataframe has been deprecated, use to_geodataframe() instead",
+            DeprecationWarning,
+        )
+        return self.to_geodataframe()
+
+    def grid_line_geodataframe(self):
+        """
+        Method to get a GeoDataFrame of grid lines
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        gdf = super().to_geodataframe(self.grid_lines, featuretype="LineString")
+        gdf = gdf.rename(columns={"node": "number"})
         return gdf
 
     def neighbors(self, node=None, **kwargs):
diff --git a/flopy/discretization/vertexgrid.py b/flopy/discretization/vertexgrid.py
index dc5be87109..136216283d 100644
--- a/flopy/discretization/vertexgrid.py
+++ b/flopy/discretization/vertexgrid.py
@@ -300,8 +300,7 @@ def map_polygons(self):
 
         return copy.copy(self._polygons)
 
-    @property
-    def geo_dataframe(self):
+    def to_geodataframe(self):
         """
         Returns a geopandas GeoDataFrame of the model grid
 
@@ -313,9 +312,38 @@ def geo_dataframe(self):
         featuretype = "Polygon"
         if self._cell1d is not None:
             featuretype = "multilinestring"
-        gdf = super().geo_dataframe(cells, featuretype)
+        gdf = super().to_geodataframe(cells, featuretype)
         return gdf
 
+    def grid_line_geodataframe(self):
+        """
+        Method to get a GeoDataFrame of grid lines
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        gdf = super().to_geodataframe(self.grid_lines, featuretype="LineString")
+        gdf = gdf.rename(columns={"node": "number"})
+        return gdf
+
+    @property
+    def geo_dataframe(self):
+        """
+        Returns a geopandas GeoDataFrame of the model grid
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        import warnings
+
+        warnings.warn(
+            "geo_dataframe has been deprecated, use to_geodataframe() instead",
+            DeprecationWarning,
+        )
+        return self.to_geodataframe()
+
     def convert_grid(self, factor):
         """
         Method to scale the model grid based on user supplied scale factors
diff --git a/flopy/export/shapefile_utils.py b/flopy/export/shapefile_utils.py
index 3c0e188279..d1b419cd5b 100644
--- a/flopy/export/shapefile_utils.py
+++ b/flopy/export/shapefile_utils.py
@@ -17,12 +17,11 @@
 import numpy as np
 
 from ..datbase import DataInterface, DataType
-from ..discretization.grid import Grid
 from ..utils import Util3d, flopy_io, import_optional_dependency
 from ..utils.crs import get_crs
 
 
-def write_gridlines_shapefile(filename: Union[str, PathLike], mg):
+def write_gridlines_shapefile(filename: Union[str, PathLike], modelgrid):
     """
     Write a polyline shapefile of the grid lines - a lightweight alternative
     to polygons.
@@ -31,7 +30,7 @@ def write_gridlines_shapefile(filename: Union[str, PathLike], mg):
     ----------
     filename : str or PathLike
         path of the shapefile to write
-    mg : flopy.discretization.grid.Grid object
+    modelgrid : flopy.discretization.grid.Grid object
         flopy model grid
 
     Returns
@@ -39,29 +38,27 @@ def write_gridlines_shapefile(filename: Union[str, PathLike], mg):
     None
 
     """
-    shapefile = import_optional_dependency("shapefile")
-    if not isinstance(mg, Grid):
+    warnings.warn(
+        "the write_gridlines_shapefile() utility will be deprecated and is "
+        "replaced by modelgrid.to_geodataframe()",
+        DeprecationWarning,
+    )
+
+    from ..discretization.grid import Grid
+
+    if not isinstance(modelgrid, Grid):
         raise ValueError(
-            f"'mg' must be a flopy Grid subclass instance; found '{type(mg)}'"
+            f"'modelgrid' must be a flopy Grid subclass instance; "
+            f"found '{type(modelgrid)}'"
         )
-    wr = shapefile.Writer(str(filename), shapeType=shapefile.POLYLINE)
-    wr.field("number", "N", 18, 0)
-    grid_lines = mg.grid_lines
-    for i, line in enumerate(grid_lines):
-        wr.line([line])
-        wr.record(i)
 
-    wr.close()
-    try:
-        write_prj(filename, modelgrid=mg)
-    except ImportError:
-        pass
-    return
+    gdf = modelgrid.grid_line_geodataframe()
+    gdf.write_file(filename)
 
 
 def write_grid_shapefile(
-    path: Union[str, PathLike],
-    mg,
+    filename: Union[str, PathLike],
+    modelgrid,
     array_dict,
     nan_val=np.nan,
     crs=None,
@@ -74,9 +71,9 @@ def write_grid_shapefile(
 
     Parameters
     ----------
-    path : str or PathLike
+    filename : str or PathLike
         shapefile file path
-    mg : flopy.discretization.grid.Grid object
+    modelgrid : flopy.discretization.grid.Grid object
         flopy model grid
     array_dict : dict
         dictionary of model input arrays
@@ -105,127 +102,54 @@ def write_grid_shapefile(
     None
 
     """
-    shapefile = import_optional_dependency("shapefile")
-    w = shapefile.Writer(str(path), shapeType=shapefile.POLYGON)
-    w.autoBalance = 1
+    warnings.warn(
+        "the write_grid_shapefile() utility will be deprecated and is "
+        "replaced by modelgrid.to_geodataframe()",
+        DeprecationWarning,
+    )
+
+    from ..discretization.grid import Grid
 
-    if not isinstance(mg, Grid):
+    if not isinstance(modelgrid, Grid):
         raise ValueError(
-            f"'mg' must be a flopy Grid subclass instance; found '{type(mg)}'"
+            f"'modelgrid' must be a flopy Grid subclass instance; "
+            f"found '{type(modelgrid)}'"
         )
-    elif mg.grid_type == "structured":
-        verts = [
-            mg.get_cell_vertices(i, j) for i in range(mg.nrow) for j in range(mg.ncol)
-        ]
-    elif mg.grid_type == "vertex":
-        verts = [mg.get_cell_vertices(cellid) for cellid in range(mg.ncpl)]
-    elif mg.grid_type == "unstructured":
-        verts = [mg.get_cell_vertices(cellid) for cellid in range(mg.nnodes)]
-    else:
-        raise NotImplementedError(f"Grid type {mg.grid_type} not supported.")
-
-    # set up the attribute fields and arrays of attributes
-    if mg.grid_type == "structured":
-        names = ["node", "row", "column"] + list(array_dict.keys())
-        dtypes = [
-            ("node", np.dtype("int")),
-            ("row", np.dtype("int")),
-            ("column", np.dtype("int")),
-        ] + [
-            (enforce_10ch_limit([name])[0], array_dict[name].dtype)
-            for name in names[3:]
-        ]
-        node = list(range(1, mg.ncol * mg.nrow + 1))
-        col = list(range(1, mg.ncol + 1)) * mg.nrow
-        row = sorted(list(range(1, mg.nrow + 1)) * mg.ncol)
-        at = np.vstack(
-            [node, row, col] + [array_dict[name].ravel() for name in names[3:]]
-        ).transpose()
-
-        names = enforce_10ch_limit(names)
-
-    elif mg.grid_type == "vertex":
-        names = ["node"] + list(array_dict.keys())
-        dtypes = [("node", np.dtype("int"))] + [
-            (enforce_10ch_limit([name])[0], array_dict[name].dtype)
-            for name in names[1:]
-        ]
-        node = list(range(1, mg.ncpl + 1))
-        at = np.vstack(
-            [node] + [array_dict[name].ravel() for name in names[1:]]
-        ).transpose()
-
-        names = enforce_10ch_limit(names)
-
-    elif mg.grid_type == "unstructured":
-        if mg.nlay is None:
-            names = ["node"] + list(array_dict.keys())
-            dtypes = [("node", np.dtype("int"))] + [
-                (enforce_10ch_limit([name])[0], array_dict[name].dtype)
-                for name in names[1:]
-            ]
-            node = list(range(1, mg.nnodes + 1))
-            at = np.vstack(
-                [node] + [array_dict[name].ravel() for name in names[1:]]
-            ).transpose()
+    gdf = modelgrid.to_geodataframe()
+    names = list(array_dict.keys())
+    at = np.vstack([array_dict[name].ravel() for name in names])
+    names = enforce_10ch_limit(names)
+
+    for ix, name in enumerate(names):
+        gdf[name] = at[ix]
+
+    gdf = gdf.fillna(nan_val)
+
+    if "epsg" in kwargs:
+        epsg = kwargs.pop("epsg")
+        crs = f"EPSG:{epsg}"
+
+    if crs is not None:
+        if gdf.crs is None:
+            gdf = gdf.set_crs(crs)
         else:
-            names = ["node", "layer"] + list(array_dict.keys())
-            dtypes = [("node", np.dtype("int")), ("layer", np.dtype("int"))] + [
-                (enforce_10ch_limit([name])[0], array_dict[name].dtype)
-                for name in names[2:]
-            ]
-            node = list(range(1, mg.nnodes + 1))
-            layer = np.zeros(mg.nnodes)
-            for ilay in range(mg.nlay):
-                istart, istop = mg.get_layer_node_range(ilay)
-                layer[istart:istop] = ilay + 1
-            at = np.vstack(
-                [node] + [layer] + [array_dict[name].ravel() for name in names[2:]]
-            ).transpose()
-
-        names = enforce_10ch_limit(names)
-
-    # flag nan values and explicitly set the dtypes
-    if at.dtype in [float, np.float32, np.float64]:
-        at[np.isnan(at)] = nan_val
-    at = np.array([tuple(i) for i in at], dtype=dtypes)
-
-    # write field information
-    fieldinfo = {name: get_pyshp_field_info(dtype.name) for name, dtype in dtypes}
-    for n in names:
-        w.field(n, *fieldinfo[n])
-
-    for i, r in enumerate(at):
-        # check if polygon is closed, if not close polygon for QGIS
-        if verts[i][-1] != verts[i][0]:
-            verts[i] = verts[i] + [verts[i][0]]
-        w.poly([verts[i]])
-        w.record(*r)
-
-    # close
-    w.close()
+            gdf = gdf.to_crs(crs)
+
+    gdf.to_file(filename)
+
     if verbose:
-        print(f"wrote {flopy_io.relpath_safe(path)}")
+        print(f"wrote {flopy_io.relpath_safe(os.getcwd(), filename)}")
 
-    # handle deprecated projection kwargs; warnings are raised in crs.py
-    write_prj_args = {}
-    if "epsg" in kwargs:
-        write_prj_args["epsg"] = kwargs.pop("epsg")
-    if "prj" in kwargs:
-        write_prj_args["prj"] = kwargs.pop("prj")
-    if kwargs:
-        raise TypeError(f"unhandled keywords: {kwargs}")
-    # write the projection file
-    try:
-        write_prj(path, mg, crs=crs, prjfile=prjfile, **write_prj_args)
-    except ImportError:
-        if verbose:
-            print("projection file not written")
-    return
+    if "prj" in kwargs or "prjfile" in kwargs or "wkt_string" in kwargs:
+        try:
+            write_prj(filename, modelgrid, crs=crs, prjfile=prjfile, **kwargs)
+        except ImportError:
+            if verbose:
+                print("projection file not written")
 
 
 def model_attributes_to_shapefile(
-    path: Union[str, PathLike],
+    filename: Union[str, PathLike],
     ml,
     package_names=None,
     array_dict=None,
@@ -239,7 +163,7 @@ def model_attributes_to_shapefile(
 
     Parameters
     ----------
-    path : str or PathLike
+    filename : str or PathLike
         path to write the shapefile to
     ml : flopy.mbase
         model instance
@@ -276,6 +200,11 @@ def model_attributes_to_shapefile(
     >>> flopy.utils.model_attributes_to_shapefile('model.shp', m)
 
     """
+    warnings.warn(
+        "model_attributes_to_shapefile is deprecated, please use the built in "
+        "to_geodataframe() method on the model object",
+        DeprecationWarning,
+    )
 
     if array_dict is None:
         array_dict = {}
@@ -286,143 +215,40 @@ def model_attributes_to_shapefile(
     else:
         package_names = [pak.name[0] for pak in ml.packagelist]
 
-    if "modelgrid" in kwargs:
-        grid = kwargs.pop("modelgrid")
-    else:
-        grid = ml.modelgrid
-
-    horz_shape = grid.get_plottable_layer_shape()
-    for pname in package_names:
-        pak = ml.get_package(pname)
-        attrs = dir(pak)
-        if pak is not None:
-            if "start_datetime" in attrs:
-                attrs.remove("start_datetime")
-            for attr in attrs:
-                a = getattr(pak, attr)
-                if a is None or not hasattr(a, "data_type") or a.name == "thickness":
-                    continue
-                if a.data_type == DataType.array2d:
-                    if a.array is None or a.array.shape != horz_shape:
-                        warn(
-                            "Failed to get data for "
-                            f"{a.name} array, {pak.name[0]} package"
-                        )
-                        continue
-                    name = shape_attr_name(a.name, keep_layer=True)
-                    array_dict[name] = a.array
-                elif a.data_type == DataType.array3d:
-                    # Not sure how best to check if an object has array data
-                    if a.array is None:
-                        warn(
-                            "Failed to get data for "
-                            f"{a.name} array, {pak.name[0]} package"
-                        )
-                        continue
-                    if isinstance(a.name, list) and a.name[0] == "thickness":
-                        continue
-                    if a.array.shape == horz_shape:
-                        if hasattr(a, "shape"):
-                            if a.shape[1] is None:  # usg unstructured Util3d
-                                # return a flattened array,
-                                # with a.name[0] (a per-layer list)
-                                array_dict[a.name[0]] = a.array
-                            else:
-                                array_dict[a.name] = a.array
-                        else:
-                            array_dict[a.name] = a.array
-                    else:
-                        # array is not the same shape as the layer shape
-                        for ilay in range(a.array.shape[0]):
-                            try:
-                                arr = a.array[ilay]
-                            except:
-                                arr = a[ilay]
-
-                            if isinstance(a, Util3d):
-                                aname = shape_attr_name(a[ilay].name)
-                            else:
-                                aname = a.name
-
-                            if arr.shape == (1,) + horz_shape:
-                                # fix for mf6 case
-                                arr = arr[0]
-                            assert arr.shape == horz_shape
-                            name = f"{aname}_{ilay + 1}"
-                            array_dict[name] = arr
-                elif a.data_type == DataType.transient2d:
-                    # Not sure how best to check if an object has array data
-                    try:
-                        assert a.array is not None
-                    except:
-                        warn(
-                            "Failed to get data for "
-                            f"{a.name} array, {pak.name[0]} package"
-                        )
-                        continue
-                    for kper in range(a.array.shape[0]):
-                        name = f"{shape_attr_name(a.name)}{kper + 1}"
-                        arr = a.array[kper][0]
-                        assert arr.shape == horz_shape
-                        array_dict[name] = arr
-                elif a.data_type == DataType.transientlist:
-                    # Skip empty transientlist
-                    if not a.data:
-                        continue
-
-                    # Use first recarray kper to check transientlist
-                    for kper in a.data.keys():
-                        if isinstance(a.data[kper], np.recarray):
-                            break
-                    # Skip transientlist if all elements are of object type
-                    if all(
-                        dtype == np.object_
-                        for dtype, _ in a.data[kper].dtype.fields.values()
-                    ):
-                        continue
-
-                    for name, array in a.masked_4D_arrays_itr():
-                        n = shape_attr_name(name, length=4)
-                        for kper in range(array.shape[0]):
-                            # guard clause for disu case
-                            # array is (kper, node) only
-                            if len(array.shape) == 2:
-                                aname = f"{n}{kper + 1}"
-                                arr = array[kper]
-                                assert arr.shape == horz_shape
-                                if np.all(np.isnan(arr)):
-                                    continue
-                                array_dict[aname] = arr
-                                continue
-                            # non-disu case
-                            for k in range(array.shape[1]):
-                                aname = f"{n}{k + 1}{kper + 1}"
-                                arr = array[kper][k]
-                                assert arr.shape == horz_shape
-                                if np.all(np.isnan(arr)):
-                                    continue
-                                array_dict[aname] = arr
-                elif isinstance(a, list):
-                    for v in a:
-                        if (
-                            isinstance(a, DataInterface)
-                            and v.data_type == DataType.array3d
-                        ):
-                            for ilay in range(a.model.modelgrid.nlay):
-                                u2d = a[ilay]
-                                name = f"{shape_attr_name(u2d.name)}_{ilay + 1}"
-                                arr = u2d.array
-                                assert arr.shape == horz_shape
-                                array_dict[name] = arr
-
-    # write data arrays to a shapefile
-    write_grid_shapefile(path, grid, array_dict)
+    gdf = ml.to_geodataframe(package_names=package_names, truncate_attrs=True)
+
+    if array_dict:
+        modelgrid = ml.modelgrid
+        for name, array in array_dict.items():
+            if modelgrid.grid_type() == "unstructured":
+                gdf[name] = array.ravel()
+            else:
+                if array.size == modelgrid.ncpl:
+                    gdf[name] = array.ravel()
+                elif array.size % modelgrid.ncpl == 0:
+                    array = array.reshape((-1, modelgrid.ncpl))
+                    for ix, lay in enumerate(array):
+                        gdf[f"{name}_{ix}"] = lay
+                else:
+                    raise ValueError(
+                        f"{name} array size is not a multiple of ncpl {modelgrid.ncpl}"
+                    )
+
     crs = kwargs.get("crs", None)
+    if crs is not None:
+        if gdf.crs is None:
+            gdf = gdf.set_crs(crs)
+        else:
+            gdf = gdf.to_crs(crs)
+
+    gdf.to_file(filename)
+
     prjfile = kwargs.get("prjfile", None)
-    try:
-        write_prj(path, grid, crs=crs, prjfile=prjfile)
-    except ImportError:
-        pass
+    if prjfile is not None:
+        try:
+            write_prj(filename, ml.modelgrid, crs=crs, prjfile=prjfile)
+        except ImportError:
+            pass
 
 
 def shape_attr_name(name, length=6, keep_layer=False):
@@ -506,42 +332,12 @@ def truncate(s):
     return names
 
 
-def get_pyshp_field_info(dtypename):
-    """Get pyshp dtype information for a given numpy dtype."""
-    fields = {
-        "int": ("N", 18, 0),
-        "<i": ("N", 18, 0),
-        "float": ("F", 20, 12),
-        "<f": ("F", 20, 12),
-        "bool": ("L", 1),
-        "b1": ("L", 1),
-        "str": ("C", 50),
-        "object": ("C", 50),
-    }
-    k = [k for k in fields.keys() if k in dtypename.lower()]
-    if len(k) == 1:
-        return fields[k[0]]
-    else:
-        return fields["str"]
-
-
-def get_pyshp_field_dtypes(code):
-    """Returns a numpy dtype for a pyshp field type."""
-    dtypes = {
-        "N": int,
-        "F": float,
-        "L": bool,
-        "C": object,
-    }
-    return dtypes.get(code, object)
-
-
-def shp2recarray(shpname: Union[str, PathLike]):
+def shp2recarray(filename: Union[str, PathLike]):
     """Read a shapefile into a numpy recarray.
 
     Parameters
     ----------
-    shpname : str or PathLike
+    filename : str or PathLike
         ESRI Shapefile path
 
     Returns
@@ -549,26 +345,39 @@ def shp2recarray(shpname: Union[str, PathLike]):
     np.recarray
 
     """
+    warnings.warn(
+        "shp2recarray will be deprecated, shapefiles can be read in using "
+        "geopandas standard methods. e.g., gpd.read_file(filename)",
+        DeprecationWarning,
+    )
+
     from ..utils.geospatial_utils import GeoSpatialCollection
 
-    sf = import_optional_dependency("shapefile")
+    gpd = import_optional_dependency("geopandas")
+    gdf = gpd.read_file(filename)
+    recarray = gdf.to_records()
+    return recarray
 
-    sfobj = sf.Reader(str(shpname))
+    """
+    sfobj = sf.Reader(str(filename))
     dtype = [(str(f[0]), get_pyshp_field_dtypes(f[1])) for f in sfobj.fields[1:]]
 
     geoms = GeoSpatialCollection(sfobj).flopy_geometry
     records = [tuple(r) + (geoms[i],) for i, r in enumerate(sfobj.iterRecords())]
     dtype += [("geometry", object)]
 
+
+
     recarray = np.array(records, dtype=dtype).view(np.recarray)
     return recarray
+    """
 
 
 def recarray2shp(
     recarray,
     geoms,
-    shpname: Union[str, PathLike] = "recarray.shp",
-    mg=None,
+    filename: Union[str, PathLike] = "recarray.shp",
+    modelgrid=None,
     crs=None,
     prjfile: Union[str, PathLike, None] = None,
     verbose=False,
@@ -590,9 +399,9 @@ def recarray2shp(
             list of shapefile.Shape objects, or geojson geometry collection
         The number of geometries in geoms must equal the number of records in
         recarray.
-    shpname : str or PathLike, default "recarray.shp"
+    filename : str or PathLike, default "recarray.shp"
         Path for the output shapefile
-    mg : flopy.discretization.Grid object
+    modelgrid : flopy.discretization.Grid object
         flopy model grid
     crs : pyproj.CRS, int, str, optional if `prjfile` is specified
         Coordinate reference system (CRS) for the model grid
@@ -624,77 +433,43 @@ def recarray2shp(
     if len(recarray) == 0:
         raise Exception("Recarray is empty")
 
-    geomtype = None
+    # set up for geopandas
+    gdf = GeoSpatialCollection(geoms).geodataframe
+    names = enforce_10ch_limit(recarray.dtype.names)
 
-    geoms = GeoSpatialCollection(geoms).flopy_geometry
+    for ix, name in enumerate(names):
+        ra_name = recarray.dtype.names[ix]
+        gdf[name] = recarray[ra_name]
 
-    for g in geoms:
-        try:
-            geomtype = g.shapeType
-        except AttributeError:
-            continue
+    if "epsg" in kwargs:
+        epsg = kwargs.pop("epsg")
+        crs = f"EPSG:{epsg}"
 
-    # set up for pyshp 2
-    shapefile = import_optional_dependency("shapefile")
-    w = shapefile.Writer(str(shpname), shapeType=geomtype)
-    w.autoBalance = 1
+    if crs is not None:
+        gdf = gdf.set_crs(crs)
+
+    gdf.to_file(filename)
 
-    # set up the attribute fields
-    names = enforce_10ch_limit(recarray.dtype.names)
-    for i, npdtype in enumerate(recarray.dtype.descr):
-        key = names[i]
-        if not isinstance(key, str):
-            key = str(key)
-        w.field(key, *get_pyshp_field_info(npdtype[1]))
-
-    # write the geometry and attributes for each record
-    ralist = recarray.tolist()
-    if geomtype == shapefile.POLYGON:
-        for i, r in enumerate(ralist):
-            w.poly(geoms[i].pyshp_parts)
-            w.record(*r)
-    elif geomtype == shapefile.POLYLINE:
-        for i, r in enumerate(ralist):
-            w.line(geoms[i].pyshp_parts)
-            w.record(*r)
-    elif geomtype == shapefile.POINT:
-        # pyshp version 2.x w.point() method can only take x and y
-        # code will need to be refactored in order to write POINTZ
-        # shapes with the z attribute.
-        for i, r in enumerate(ralist):
-            w.point(*geoms[i].pyshp_parts[:2])
-            w.record(*r)
-
-    w.close()
     if verbose:
-        print(f"wrote {flopy_io.relpath_safe(os.getcwd(), shpname)}")
+        print(f"wrote {flopy_io.relpath_safe(os.getcwd(), filename)}")
 
-    # handle deprecated projection kwargs; warnings are raised in crs.py
-    write_prj_args = {}
-    if "epsg" in kwargs:
-        write_prj_args["epsg"] = kwargs.pop("epsg")
-    if "prj" in kwargs:
-        write_prj_args["prj"] = kwargs.pop("prj")
-    if kwargs:
-        raise TypeError(f"unhandled keywords: {kwargs}")
-    # write the projection file
-    try:
-        write_prj(shpname, mg, crs=crs, prjfile=prjfile, **write_prj_args)
-    except ImportError:
-        if verbose:
-            print("projection file not written")
-    return
+    if "prj" in kwargs or "prjfile" in kwargs or "wkt_string" in kwargs:
+        try:
+            write_prj(filename, modelgrid, crs=crs, prjfile=prjfile, **kwargs)
+        except ImportError:
+            if verbose:
+                print("projection file not written")
 
 
 def write_prj(
-    shpname,
+    filename,
     modelgrid=None,
     crs=None,
     prjfile=None,
     **kwargs,
 ):
     # projection file name
-    output_projection_file = Path(shpname).with_suffix(".prj")
+    output_projection_file = Path(filename).with_suffix(".prj")
 
     # handle deprecated projection kwargs; warnings are raised in crs.py
     get_crs_args = {}
diff --git a/flopy/export/utils.py b/flopy/export/utils.py
index 6a99e550fe..f62d9d45fe 100644
--- a/flopy/export/utils.py
+++ b/flopy/export/utils.py
@@ -554,7 +554,15 @@ def output_helper(
                             attrib_dict[name] = plotarray[per][k]
 
         if attrib_dict:
-            shapefile_utils.write_grid_shapefile(f, ml.modelgrid, attrib_dict)
+            gdf = ml.modelgrid.to_geodataframe()
+            names = list(attrib_dict.keys())
+            at = np.vstack([attrib_dict[name].ravel() for name in names])
+            names = shapefile_utils.enforce_10ch_limit(names)
+
+            for ix, name in enumerate(names):
+                gdf[name] = at[ix]
+
+            gdf.to_file(f)
 
     else:
         msg = f"unrecognized export argument:{f}"
@@ -603,9 +611,13 @@ def model_export(f: Union[str, PathLike, NetCdf, dict], ml, fmt=None, **kwargs):
         f = NetCdf(f, ml, **kwargs)
 
     if isinstance(f, (str, PathLike)) and Path(f).suffix.lower() == ".shp":
-        shapefile_utils.model_attributes_to_shapefile(
-            f, ml, package_names=package_names, **kwargs
-        )
+        modelgrid = kwargs.pop("modelgrid", None)
+        if modelgrid is None:
+            gdf = None
+        else:
+            gdf = modelgrid.to_geodataframe()
+        gdf = ml.to_geodataframe(gdf=gdf, package_names=package_names, **kwargs)
+        gdf.to_file(f)
 
     elif isinstance(f, NetCdf):
         for pak in ml.packagelist:
@@ -692,9 +704,26 @@ def package_export(
         f = NetCdf(f, pak.parent, **kwargs)
 
     if isinstance(f, (str, PathLike)) and Path(f).suffix.lower() == ".shp":
-        shapefile_utils.model_attributes_to_shapefile(
-            f, pak.parent, package_names=pak.name, verbose=verbose, **kwargs
-        )
+        from .shapefile_utils import write_prj
+
+        gdf = pak.to_geodataframe()
+
+        crs = kwargs.get("crs", None)
+        if crs is not None:
+            if gdf.crs is None:
+                gdf = gdf.set_crs(crs)
+            else:
+                gdf = gdf.to_crs(crs)
+
+        gdf.to_file(f)
+
+        prjfile = kwargs.get("prjfile", None)
+
+        if prjfile is not None:
+            try:
+                write_prj(f, pak.parent.modelgrid, crs=crs, prjfile=prjfile)
+            except ImportError:
+                pass
 
     elif isinstance(f, NetCdf) or isinstance(f, dict):
         for a in pak.data_list:
@@ -869,87 +898,15 @@ def mflist_export(f: Union[str, PathLike, NetCdf], mfl, **kwargs):
         kper = kwargs.get("kper", 0)
         squeeze = kwargs.get("squeeze", True)
 
-        if modelgrid is None:
-            raise Exception("MfList.to_shapefile: modelgrid is not set")
-        elif modelgrid.grid_type == "USG-Unstructured":
-            raise Exception(
-                "Flopy does not support exporting to shapefile "
-                "from a MODFLOW-USG unstructured grid."
-            )
-
-        if kper is None:
-            keys = mfl.data.keys()
-            keys.sort()
-        else:
-            keys = [kper]
-        if not sparse:
-            array_dict = {}
-            for kk in keys:
-                arrays = mfl.to_array(kk)
-                for name, array in arrays.items():
-                    for k in range(array.shape[0]):
-                        n = shapefile_utils.shape_attr_name(name, length=4)
-                        aname = f"{n}{k + 1}{int(kk) + 1}"
-                        array_dict[aname] = array[k]
-            shapefile_utils.write_grid_shapefile(f, modelgrid, array_dict)
+        if modelgrid is not None:
+            gdf = modelgrid.to_geodataframe()
         else:
-            from ..export.shapefile_utils import recarray2shp
-            from ..utils.geometry import Polygon
+            gdf = None
 
-            df = (
-                mfl.get_dataframe(squeeze=squeeze)
-                if isinstance(mfl, MfList)
-                else mfl.get_dataframe()
-            )
-            if isinstance(df, dict):
-                df = df[kper]
-            if "kper" in kwargs or df is None:
-                ra = mfl[kper]
-                verts = np.array(modelgrid.get_cell_vertices(ra.i, ra.j))
-            elif df is not None:
-                if modelgrid.grid_type == "unstructured":
-                    node_index_name = (
-                        "cellid_node" if "cellid_node" in df.columns else "node"
-                    )
-                    verts = np.array(
-                        [
-                            modelgrid.get_cell_vertices(node)
-                            for node in df[node_index_name].values
-                        ]
-                    )
-                elif modelgrid.grid_type == "vertex":
-                    layer_index_name = (
-                        "cellid_layer" if "cellid_layer" in df.columns else "layer"
-                    )
-                    cell_index_name = (
-                        "cellid_cell" if "cellid_cell" in df.columns else "cell"
-                    )
-                    verts = np.array(
-                        [
-                            modelgrid.get_cell_vertices(layer * modelgrid.ncpl + cell)
-                            for layer, cell in zip(
-                                df[layer_index_name].values, df[cell_index_name].values
-                            )
-                        ]
-                    )
-                else:
-                    row_index_name = "i" if "i" in df.columns else "cellid_row"
-                    col_index_name = "j" if "j" in df.columns else "cellid_column"
-                    verts = np.array(
-                        [
-                            modelgrid.get_cell_vertices(
-                                i, df[col_index_name].values[ix]
-                            )
-                            for ix, i in enumerate(df[row_index_name].values)
-                        ]
-                    )
-                ra = df.to_records(index=False)
-            crs = kwargs.get("crs", None)
-            prjfile = kwargs.get("prjfile", None)
-            polys = np.array([Polygon(v) for v in verts])
-            recarray2shp(
-                ra, geoms=polys, shpname=f, mg=modelgrid, crs=crs, prjfile=prjfile
-            )
+        gdf = mfl.to_geodataframe(
+            gdf=gdf, kper=kper, sparse=sparse, truncate_attrs=True
+        )
+        gdf.to_file(f)
 
     elif isinstance(f, NetCdf) or isinstance(f, dict):
         base_name = mfl.package.name[0].lower()
@@ -1185,24 +1142,13 @@ def array3d_export(f: Union[str, PathLike], u3d, fmt=None, **kwargs):
         f = NetCdf(f, u3d.model, **kwargs)
 
     if isinstance(f, (str, PathLike)) and Path(f).suffix.lower() == ".shp":
-        array_dict = {}
-        array_shape = u3d.array.shape
-
-        if len(array_shape) == 1:
-            name = shapefile_utils.shape_attr_name(u3d.name)
-            array_dict[name] = u3d.array
+        if modelgrid is not None:
+            gdf = modelgrid.to_geodataframe()
         else:
-            for ilay in range(array_shape[0]):
-                u2d = u3d[ilay]
-                if isinstance(u2d, np.ndarray):
-                    dname = u3d.name
-                    array = u2d
-                else:
-                    dname = u2d.name
-                    array = u2d.array
-                name = f"{shapefile_utils.shape_attr_name(dname)}_{ilay + 1}"
-                array_dict[name] = array
-        shapefile_utils.write_grid_shapefile(f, modelgrid, array_dict)
+            gdf = None
+
+        gdf = u3d.to_geodataframe(gdf=gdf, truncate_attrs=True)
+        gdf.to_file(f)
 
     elif isinstance(f, NetCdf) or isinstance(f, dict):
         var_name = u3d.name
@@ -1613,6 +1559,7 @@ def export_contours(
     contours,
     fieldname="level",
     verbose=False,
+    crs=None,
     **kwargs,
 ):
     """
@@ -1628,7 +1575,13 @@ def export_contours(
         gis attribute table field name
     verbose : bool, optional, default False
         whether to show verbose output
-    **kwargs : key-word arguments to flopy.export.shapefile_utils.recarray2shp
+    crs : pyproj.CRS, int, str, optional if `prjfile` is specified
+        Coordinate reference system (CRS) for the model grid
+        (must be projected; geographic CRS are not supported).
+        The value can be anything accepted by
+        :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+        such as an authority string (eg "EPSG:26916") or a WKT string.
+    **kwargs :
 
     Returns
     -------
@@ -1640,7 +1593,10 @@ def export_contours(
     from matplotlib.path import Path
 
     from ..utils.geometry import LineString
-    from .shapefile_utils import recarray2shp
+    from ..utils.geospatial_utils import GeoSpatialCollection
+    from ..utils.utl_import import import_optional_dependency
+
+    gpd = import_optional_dependency("geopandas")
 
     if not isinstance(contours, list):
         contours = [contours]
@@ -1702,12 +1658,18 @@ def export_contours(
     if verbose:
         print(f"Writing {len(level)} contour lines")
 
-    ra = np.array(level, dtype=[(fieldname, float)]).view(np.recarray)
-
-    recarray2shp(ra, geoms, filename, **kwargs)
+    geoms = GeoSpatialCollection(geoms, "LineString")
+    gdf = gpd.GeoDataFrame.from_features(geoms)
+    gdf[fieldname] = level
+    if crs is not None:
+        gdf = gdf.set_crs(crs)
+    gdf.to_file(filename)
+    return gdf
 
 
-def export_contourf(filename, contours, fieldname="level", verbose=False, **kwargs):
+def export_contourf(
+    filename, contours, fieldname="level", verbose=False, crs=None, **kwargs
+):
     """
     Write matplotlib filled contours to shapefile.
 
@@ -1723,11 +1685,17 @@ def export_contourf(filename, contours, fieldname="level", verbose=False, **kwar
         the range represented by the polygon.  Default is 'level'.
     verbose : bool, optional, default False
         whether to show verbose output
+    crs : pyproj.CRS, int, str, optional if `prjfile` is specified
+        Coordinate reference system (CRS) for the model grid
+        (must be projected; geographic CRS are not supported).
+        The value can be anything accepted by
+        :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+        such as an authority string (eg "EPSG:26916") or a WKT string.
     **kwargs : keyword arguments to flopy.export.shapefile_utils.recarray2shp
 
     Returns
     -------
-    None
+    gdf : GeoDataFrame of matplotlib contours
 
     Examples
     --------
@@ -1744,7 +1712,10 @@ def export_contourf(filename, contours, fieldname="level", verbose=False, **kwar
     from matplotlib.path import Path
 
     from ..utils.geometry import Polygon, is_clockwise
-    from .shapefile_utils import recarray2shp
+    from ..utils.geospatial_utils import GeoSpatialCollection
+    from ..utils.utl_import import import_optional_dependency
+
+    gpd = import_optional_dependency("geopandas")
 
     if not isinstance(contours, list):
         contours = [contours]
@@ -1829,9 +1800,13 @@ def export_contourf(filename, contours, fieldname="level", verbose=False, **kwar
     if verbose:
         print(f"Writing {len(level)} polygons")
 
-    ra = np.array(level, dtype=[(fieldname, float)]).view(np.recarray)
-
-    recarray2shp(ra, geoms, filename, **kwargs)
+    geoms = GeoSpatialCollection(geoms, "Polygon")
+    gdf = gpd.GeoDataFrame.from_features(geoms)
+    gdf[fieldname] = level
+    if crs is not None:
+        gdf = gdf.set_crs(crs)
+    gdf.to_file(filename)
+    return gdf
 
 
 def export_array_contours(
@@ -1842,6 +1817,7 @@ def export_array_contours(
     interval=None,
     levels=None,
     maxlevels=1000,
+    crs=None,
     **kwargs,
 ):
     """
@@ -1863,7 +1839,13 @@ def export_array_contours(
         list of contour levels
     maxlevels : int
         maximum number of contour levels
-    **kwargs : keyword arguments to flopy.export.shapefile_utils.recarray2shp
+    crs : pyproj.CRS, int, str, optional if `prjfile` is specified
+        Coordinate reference system (CRS) for the model grid
+        (must be projected; geographic CRS are not supported).
+        The value can be anything accepted by
+        :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+        such as an authority string (eg "EPSG:26916") or a WKT string.
+    **kwargs :
 
     """
     import matplotlib.pyplot as plt
@@ -1880,7 +1862,7 @@ def export_array_contours(
     ctr = contour_array(modelgrid, ax, a, layer, levels=levels)
 
     kwargs["mg"] = modelgrid
-    export_contours(filename, ctr, fieldname, **kwargs)
+    export_contours(filename, ctr, fieldname, crs=crs, **kwargs)
     plt.close()
 
 
diff --git a/flopy/mbase.py b/flopy/mbase.py
index ef5a26b4be..60b6566bb2 100644
--- a/flopy/mbase.py
+++ b/flopy/mbase.py
@@ -760,6 +760,56 @@ def _output_msg(self, i, add=True):
             f"{txt2} the output list."
         )
 
+    def to_geodataframe(
+        self, gdf=None, kper=0, package_names=None, truncate_attrs=False
+    ):
+        """
+        Method to build a Geodataframe from model inputs. Note: transient data
+        will only be exported for a single stress period.
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional geopandas geodataframe object to add data to. Default is None
+        kper : int
+            stress period to get transient data from
+        package_names : list
+            optional list of package names
+        truncate_attrs : bool
+            method to truncate attribute names for shapefile attribute name length
+            restrictions
+
+        Returns
+        -------
+            gdf : GeoDataFrame
+        """
+        if gdf is None:
+            modelgrid = self.modelgrid
+            if modelgrid is not None:
+                gdf = modelgrid.to_geodataframe()
+            else:
+                raise AttributeError(
+                    "model does not have a grid instance, please supply a geodataframe"
+                )
+
+        if package_names is None:
+            package_names = [pak.name[0] for pak in self.packagelist]
+        else:
+            pass
+
+        for package in self.packagelist:
+            if package.package_type in ("hfb6",):
+                continue
+            if package.name[0] not in package_names:
+                # todo: this needs testing
+                continue
+            if callable(getattr(package, "to_geodataframe", None)):
+                gdf = package.to_geodataframe(
+                    gdf, kper=kper, sparse=False, truncate_attrs=truncate_attrs
+                )
+
+        return gdf
+
     def add_output_file(
         self,
         unit,
diff --git a/flopy/mf6/data/mfdataarray.py b/flopy/mf6/data/mfdataarray.py
index e00c877b1a..dc9b70ea3f 100644
--- a/flopy/mf6/data/mfdataarray.py
+++ b/flopy/mf6/data/mfdataarray.py
@@ -323,6 +323,71 @@ def data(self):
         """Returns array data.  Calls get_data with default parameters."""
         return self._get_data()
 
+    def to_geodataframe(self, gdf=None, name=None, forgive=False, truncate_attrs=False, **kwargs):
+        """
+        Method to add an input array to a geopandas GeoDataFrame
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame object
+        name : str
+            optional attribute name, default uses util2d name
+        forgive : bool
+            optional flag to continue if data shape not compatible with GeoDataFrame
+        truncate_attrs : bool
+            method to truncate attribute names for shapefile restrictions
+
+        Returns
+        -------
+            geopandas GeoDataFrame
+        """
+        from ...export.shapefile_utils import shape_attr_name
+
+        if self.model is None:
+            return gdf
+        else:
+            modelgrid = self.model.modelgrid
+            if gdf is None:
+                if modelgrid is None:
+                    return gdf
+                gdf = modelgrid.to_geodataframe()
+
+            if modelgrid is not None:
+                if modelgrid.grid_type != "unstructured":
+                    ncpl = modelgrid.ncpl
+                else:
+                    ncpl = modelgrid.nnodes
+            else:
+                ncpl = len(gdf)
+
+            if name is None:
+                name = self.name
+
+            data = self.array
+            if data is None:
+                return gdf
+
+            if truncate_attrs:
+                name = shape_attr_name(name=name)
+
+            if data.size == ncpl:
+                gdf[name] = data.ravel()
+
+            elif data.size % ncpl == 0:
+                data = data.reshape((-1, ncpl))
+                for ix, arr in enumerate(data):
+                    aname = f"{name}_{ix}"
+                    gdf[aname] = arr
+            elif forgive:
+                return gdf
+            else:
+                raise ValueError(
+                    f"Data size {data.size} not compatible with dataframe length {ncpl}"
+                )
+
+            return gdf
+
     def new_simulation(self, sim_data):
         """Initialize MFArray object for a new simulation
 
@@ -1890,6 +1955,76 @@ def _build_period_data(
                     output[sp] = data
         return output
 
+    def to_geodataframe(self, gdf=None, kper=0, forgive=False, truncate_attrs=False, **kwargs):
+        """
+        Method to add an input array to a geopandas GeoDataFrame
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame object
+        name : str
+            optional attribute name, default uses util2d name
+        kper : int
+            stress period number
+        forgive : bool
+            optional flag to continue if data shape not compatible with GeoDataFrame
+        truncate_attrs : bool
+            method to truncate attribute names for shapefile restrictions
+
+        Returns
+        -------
+            geopandas GeoDataFrame
+        """
+        from ...export.shapefile_utils import shape_attr_name
+
+        if self.model is None:
+            return gdf
+        else:
+            modelgrid = self.model.modelgrid
+            if gdf is None:
+                if modelgrid is None:
+                    return gdf
+                gdf = modelgrid.to_geodataframe()
+
+            if modelgrid is not None:
+                if modelgrid.grid_type != "unstructured":
+                    ncpl = modelgrid.ncpl
+                else:
+                    ncpl = modelgrid.nnodes
+            else:
+                ncpl = len(gdf)
+
+            if self.array is None:
+                return gdf
+
+            if truncate_attrs:
+                name = shape_attr_name(self.name, length=4)
+            else:
+                name = f"{self.path[1]}_{self.name}"
+
+            data = self.get_data(key=kper, apply_mult=True)
+            if data.size == ncpl:
+                name = f"{name}_{kper}"
+                gdf[name] = data.ravel()
+
+            elif data.size % ncpl == 0:
+                data = data.reshape((-1, ncpl))
+                for ix, arr in enumerate(data):
+                    if truncate_attrs:
+                        aname = f"{name}{ix}{kper}"
+                    else:
+                        aname = f"{name}_{ix}_{kper}"
+                    gdf[aname] = arr
+            elif forgive:
+                return gdf
+            else:
+                raise ValueError(
+                    f"Data size {data.size} not compatible with dataframe length {ncpl}"
+                )
+
+            return gdf
+
     def set_record(self, data_record):
         """Sets data and metadata at layer `layer` and time `key` to
         `data_record`.  For unlayered data do not pass in `layer`.
diff --git a/flopy/mf6/data/mfdatalist.py b/flopy/mf6/data/mfdatalist.py
index 529177c41c..d089caa747 100644
--- a/flopy/mf6/data/mfdatalist.py
+++ b/flopy/mf6/data/mfdatalist.py
@@ -145,6 +145,63 @@ def to_array(self, kper=0, mask=False):
         model_grid = self.data_dimensions.get_model_grid()
         return list_to_array(sarr, model_grid, kper, mask)
 
+    def to_geodataframe(self, gdf=None, sparse=False, truncate_attrs=False, **kwargs):
+        """
+        Method to add data to a GeoDataFrame for exporting as a geospatial file
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame instance. If GeoDataFrame is None, one will be
+            constructed from modelgrid information
+        sparse : bool
+            boolean flag for sparse dataframe construction. Default is False
+        truncate_attrs : bool
+            method to truncate attribute names for shapefile restrictions
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        from ...export.shapefile_utils import shape_attr_name
+
+        if self.model is None:
+            return gdf
+        else:
+            modelgrid = self.model.modelgrid
+            if modelgrid is None:
+                return gdf
+
+            if gdf is None:
+                gdf = modelgrid.to_geodataframe()
+
+            data = self.to_array(mask=True)
+            if data is None:
+                return gdf
+
+            col_names = []
+            for name, array3d in data.items():
+                if truncate_attrs:
+                    aname = shape_attr_name(name)
+                else:
+                    aname = f"{self.path[1].lower()}_{name}"
+
+                if modelgrid.grid_type == "unstructured":
+                    array = array3d.ravel()
+                    gdf[aname] = array
+                    col_names.append(aname)
+                else:
+                    for lay in range(modelgrid.nlay):
+                        arr = array3d[lay].ravel()
+                        gdf[f"{aname}_{lay}"] = arr.ravel()
+                        col_names.append(f"{aname}_{lay}")
+
+            if sparse:
+                gdf = gdf.dropna(subset=col_names, how="all")
+                gdf = gdf.dropna(axis="columns", how="all")
+
+            return gdf
+
     def new_simulation(self, sim_data):
         """Initialize MFList object for a new simulation.
 
@@ -1596,6 +1653,69 @@ def to_array(self, kper=0, mask=False):
         """Returns list data as an array."""
         return super().to_array(kper, mask)
 
+    def to_geodataframe(self, gdf=None, kper=0, sparse=False, truncate_attrs=False, **kwargs):
+        """
+        Method to add data to a GeoDataFrame for exporting as a geospatial file
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame instance. If GeoDataFrame is None, one will be
+            constructed from modelgrid information
+        kper : int
+            stress period to export
+        sparse : bool
+            boolean flag for sparse dataframe construction. Default is False
+        truncate_attrs : bool
+            method to truncate attribute names for shapefile restrictions
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        from ...export.shapefile_utils import shape_attr_name
+
+        if self.model is None:
+            return gdf
+        else:
+            modelgrid = self.model.modelgrid
+            if modelgrid is None:
+                return gdf
+
+            if gdf is None:
+                gdf = modelgrid.to_geodataframe()
+
+            data = self.to_array(kper=kper, mask=True)
+            if data is None:
+                return gdf
+
+            col_names = []
+            for name, array3d in data.items():
+                if truncate_attrs:
+                    name = shape_attr_name(name, length=4)
+                else:
+                    name = f"{self.path[1].lower()}_{name}"
+
+                if modelgrid.grid_type == "unstructured":
+                    array = array3d.ravel()
+                    gdf[name] = array
+                    col_names.append(name)
+                else:
+                    for lay in range(modelgrid.nlay):
+                        arr = array3d[lay].ravel()
+                        if truncate_attrs:
+                            aname = f"{name}{lay}{kper}"
+                        else:
+                            aname = f"{name}_{lay}_{kper}"
+                        gdf[aname] = arr.ravel()
+                        col_names.append(aname)
+
+            if sparse:
+                gdf = gdf.dropna(subset=col_names, how="all")
+                gdf = gdf.dropna(axis="columns", how="all")
+
+            return gdf
+
     def remove_transient_key(self, transient_key):
         """Remove transient stress period key.  Method is used
         internally by FloPy and is not intended to the end user.
diff --git a/flopy/mf6/data/mfdataplist.py b/flopy/mf6/data/mfdataplist.py
index a7aef4625e..bb6c1e2af0 100644
--- a/flopy/mf6/data/mfdataplist.py
+++ b/flopy/mf6/data/mfdataplist.py
@@ -842,6 +842,63 @@ def to_array(self, kper=0, mask=False):
         model_grid = self.data_dimensions.get_model_grid()
         return list_to_array(sarr, model_grid, kper, mask)
 
+    def to_geodataframe(self, gdf=None, sparse=False, truncate_attrs=False, **kwargs):
+        """
+        Method to add data to a GeoDataFrame for exporting as a geospatial file
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame instance. If GeoDataFrame is None, one will be
+            constructed from modelgrid information
+        sparse : bool
+            boolean flag for sparse dataframe construction. Default is False
+        truncate_attrs : bool
+            method to truncate attribute names for shapefile restrictions
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        from ...export.shapefile_utils import shape_attr_name
+
+        if self.model is None:
+            return gdf
+        else:
+            modelgrid = self.model.modelgrid
+            if modelgrid is None:
+                return gdf
+
+            if gdf is None:
+                gdf = modelgrid.to_geodataframe()
+
+            data = self.to_array(mask=True)
+            if data is None:
+                return gdf
+
+            col_names = []
+            for name, array3d in data.items():
+                if truncate_attrs:
+                    aname = shape_attr_name(name)
+                else:
+                    aname = f"{self.path[1].lower()}_{name}"
+
+                if modelgrid.grid_type == "unstructured":
+                    array = array3d.ravel()
+                    gdf[aname] = array
+                    col_names.append(aname)
+                else:
+                    for lay in range(modelgrid.nlay):
+                        arr = array3d[lay].ravel()
+                        gdf[f"{aname}_{lay}"] = arr.ravel()
+                        col_names.append(f"{aname}_{lay}")
+
+            if sparse:
+                gdf = gdf.dropna(subset=col_names, how="all")
+                gdf = gdf.dropna(axis="columns", how="all")
+
+            return gdf
+
     def set_record(self, record, autofill=False, check_data=True):
         """Sets the contents of the data and metadata to "data_record".
         Data_record is a dictionary with has the following format:
@@ -2019,6 +2076,68 @@ def __init__(
         self.repeating = True
         self.empty_keys = {}
 
+    def to_geodataframe(self, gdf=None, kper=0, sparse=False, truncate_attrs=False, **kwargs):
+        """
+        Method to add data to a GeoDataFrame for exporting as a geospatial file
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame instance. If GeoDataFrame is None, one will be
+            constructed from modelgrid information
+        kper : int
+            stress period to export
+        sparse : bool
+            boolean flag for sparse dataframe construction. Default is False
+        truncate_attrs : bool
+            method to truncate attribute names for shapefile restrictions
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        from ...export.shapefile_utils import shape_attr_name
+
+        if self.model is None:
+            return gdf
+        else:
+            modelgrid = self.model.modelgrid
+            if modelgrid is None:
+                return gdf
+
+            if gdf is None:
+                gdf = modelgrid.to_geodataframe()
+
+            data = self.to_array(kper=kper, mask=True)
+
+            col_names = []
+            for name, array3d in data.items():
+                if truncate_attrs:
+                    name = shape_attr_name(name, length=4)
+                else:
+                    name = f"{self.path[1].lower()}_{name}"
+                if modelgrid.grid_type == "unstructured":
+                    array = array3d.ravel()
+                    aname = f"{name}_{kper}"
+                    gdf[aname] = array
+                    col_names.append(aname)
+                else:
+                    for lay in range(modelgrid.nlay):
+                        arr = array3d[lay].ravel()
+                        if truncate_attrs:
+                            aname = f"{name}{lay}{kper}"
+                        else:
+                            aname = f"{name}_{lay}_{kper}"
+                        gdf[aname] = arr.ravel()
+                        col_names.append(aname)
+
+            if sparse:
+                gdf = gdf.dropna(subset=col_names, how="all")
+                gdf = gdf.dropna(axis="columns", how="all")
+
+            return gdf
+
+
     @property
     def data_type(self):
         return DataType.transientlist
diff --git a/flopy/mf6/data/mfdatautil.py b/flopy/mf6/data/mfdatautil.py
index 9b0a885119..17ee03ab8c 100644
--- a/flopy/mf6/data/mfdatautil.py
+++ b/flopy/mf6/data/mfdatautil.py
@@ -202,6 +202,8 @@ def list_to_array(sarr, model_grid, kper=0, mask=False):
         cnt = np.zeros(shape, dtype=np.float64)
         for sp_rec in sarr:
             if sp_rec is not None:
+                if "cellid" not in sp_rec.dtype.names:
+                    return None
                 for rec in sp_rec:
                     arr[rec["cellid"]] += rec[name]
                     cnt[rec["cellid"]] += 1.0
diff --git a/flopy/mf6/mfmodel.py b/flopy/mf6/mfmodel.py
index d17695d209..cb387582de 100644
--- a/flopy/mf6/mfmodel.py
+++ b/flopy/mf6/mfmodel.py
@@ -793,6 +793,54 @@ def output(self):
         except AttributeError:
             return MF6Output(self, budgetkey=budgetkey)
 
+    def to_geodataframe(self, gdf=None, kper=0, package_names=None, truncate_attrs=False):
+        """
+        Method to build a Geodataframe from model inputs. Note: transient data
+        will only be exported for a single stress period.
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional geopandas geodataframe object to add data to. Default is None
+        kper : int
+            stress period to get transient data from
+        package_names : list
+            optional list of package names
+        truncate_attrs : bool
+            method to truncate attribute names for shapefile attribute name length
+            restrictions
+
+        Returns
+        -------
+            gdf : GeoDataFrame
+        """
+        if gdf is None:
+            modelgrid = self.modelgrid
+            if modelgrid is not None:
+                gdf = modelgrid.to_geodataframe()
+            else:
+                raise AttributeError(
+                    "model does not have a grid instance, "
+                    "please supply a geodataframe"
+                )
+
+        if package_names is None:
+            package_names = [pak.name[0] for pak in self.packagelist]
+        else:
+            pass
+
+        for package in self.packagelist:
+            if package.package_type in ("hfb",):
+                continue
+            if package.name[0] not in package_names and package.package_type not in package_names:
+                continue
+            if callable(getattr(package, "to_geodataframe", None)):
+                gdf = package.to_geodataframe(
+                    gdf, kper=kper, sparse=False, truncate_attrs=truncate_attrs
+                )
+
+        return gdf
+
     def export(self, f, **kwargs):
         """Method to export a model to a shapefile or netcdf file
 
diff --git a/flopy/mf6/mfpackage.py b/flopy/mf6/mfpackage.py
index 735ca3e416..eaa3abf08c 100644
--- a/flopy/mf6/mfpackage.py
+++ b/flopy/mf6/mfpackage.py
@@ -13,6 +13,7 @@
 from ..pakbase import PackageInterface
 from ..utils import datautil
 from ..utils.check import mf6check
+from ..utils.utl_import import import_optional_dependency
 from ..version import __version__
 from .coordinates import modeldimensions
 from .data import (
@@ -2102,6 +2103,84 @@ def package_filename_dict(self):
         )
         return self._package_container.package_filename_dict
 
+    def to_geodataframe(self, gdf=None, kper=0, sparse=False, truncate_attrs=False, **kwargs):
+        """
+        Method to create a GeoDataFrame from a modflow package
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional geopandas geodataframe object to add data to. Default is None
+        kper : int
+            stress period to get transient data from
+        sparse : bool
+            optional parameter to create a sparse geodataframe
+        truncate_attrs : bool
+            method to truncate attribute names for shapefile attribute name length
+            restrictions
+
+        Returns
+        -------
+            gdf : GeoDataFrame
+        """
+        if gdf is None:
+            if isinstance(self.parent, ModelInterface):
+                modelgrid = self.parent.modelgrid
+                if modelgrid is not None:
+                    if self.package_type == "hfb":
+                        gpd = import_optional_dependency("geopandas")
+                        from ..plot.plotutil import hfb_data_to_linework
+
+                        recarray = self.stress_period_data.data[kper]
+                        lines = hfb_data_to_linework(recarray, modelgrid)
+                        geo_interface = {"type": "FeatureCollection"}
+                        features = [
+                            {
+                                "id": f"{ix}",
+                                "geometry": {"coordinates": line, "type": "LineString"},
+                                "properties": {}
+                            }
+                            for ix, line in enumerate(lines)
+                        ]
+                        geo_interface["features"] = features
+                        gdf = gpd.GeoDataFrame.from_features(geo_interface)
+
+                        for name in recarray.dtype.names:
+                            gdf[name] = recarray[name]
+
+                        return gdf
+
+                    else:
+                        gdf = modelgrid.to_geodataframe()
+                else:
+                    raise AttributeError(
+                        "model does not have a grid instance, "
+                        "please supply a geodataframe"
+                    )
+            else:
+                raise AssertionError(
+                    "Package does not have a model instance, "
+                    "please supply a geodataframe"
+                )
+
+        for attr, value in self.__dict__.items():
+            if callable(getattr(value, "to_geodataframe", None)):
+                if isinstance(value, (ModelInterface, PackageInterface)):
+                    continue
+                # do not pass sparse in here, "sparsify" after all data has been
+                #  added to geodataframe
+                gdf = value.to_geodataframe(
+                    gdf, forgive=True, kper=kper, sparse=False, truncate_attrs=truncate_attrs
+                )
+
+        if sparse:
+            col_names = [i for i in gdf if i not in ("geometry", "node", "row", "col")]
+            gdf = gdf.dropna(subset=col_names, how="all")
+            gdf = gdf.dropna(axis="columns", how="all")
+
+        return gdf
+
+
     def get_package(self, name=None, type_only=False, name_only=False):
         """
         Finds a package by package name, package key, package type, or partial
diff --git a/flopy/modflow/mfhfb.py b/flopy/modflow/mfhfb.py
index 468cf1a421..70cb7cae4a 100644
--- a/flopy/modflow/mfhfb.py
+++ b/flopy/modflow/mfhfb.py
@@ -16,6 +16,7 @@
 from ..pakbase import Package
 from ..utils.flopy_io import line_parse
 from ..utils.recarray_utils import create_empty_recarray
+from ..utils.utl_import import import_optional_dependency
 from .mfparbc import ModflowParBc as mfparbc
 
 
@@ -181,6 +182,47 @@ def _ncells(self):
         """
         return self.nhfbnp
 
+    def _get_hfb_lines(self):
+        """
+        Method to get hfb lines. Lines are ordered by recarray index
+
+        Returns
+        -------
+            list : list of hfb lines
+        """
+        from ..plot.plotutil import hfb_data_to_linework
+
+        return hfb_data_to_linework(self.hfb_data, self.parent.modelgrid)
+
+    def to_geodataframe(self, **kwargs):
+        """
+        Method to create a LineString based geodataframe of horizontal flow barriers
+
+        Returns
+        -------
+            GeoDataFrame
+
+        """
+        gpd = import_optional_dependency("geopandas")
+
+        lines = self._get_hfb_lines()
+        geo_interface = {"type": "FeatureCollection"}
+        features = [
+            {
+                "id": f"{ix}",
+                "geometry": {"coordinates": line, "type": "LineString"},
+                "properties": {},
+            }
+            for ix, line in enumerate(lines)
+        ]
+        geo_interface["features"] = features
+        gdf = gpd.GeoDataFrame.from_features(geo_interface)
+        hfbs = self.hfb_data
+        for name in hfbs.dtype.names:
+            gdf[name] = hfbs[name]
+
+        return gdf
+
     def write_file(self):
         """
         Write the package file.
diff --git a/flopy/modflow/mfsfr2.py b/flopy/modflow/mfsfr2.py
index 306d77a661..57df8872d5 100644
--- a/flopy/modflow/mfsfr2.py
+++ b/flopy/modflow/mfsfr2.py
@@ -622,6 +622,43 @@ def paths(self):
     def df(self):
         return pd.DataFrame(self.reach_data)
 
+    def to_geodataframe(self, gdf=None, sparse=True, **kwargs):
+        """
+        Method to export SFR reach data to a GeoDataFrame
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame instance. If GeoDataFrame is None, one will be
+            constructed from modelgrid information
+        sparse : bool
+            boolean flag for sparse dataframe construction. Default is True
+        """
+        modelgrid = self.parent.modelgrid
+        if modelgrid is None:
+            return gdf
+
+        if gdf is None:
+            gdf = modelgrid.to_geodataframe()
+
+        df = self.df
+        if "k" in list(df):
+            df["node"] = df["node"] - (modelgrid.ncpl * df["k"])
+            df = df.drop(columns=["k", "i", "j"])
+
+        df["node"] += 1
+        gdf = gdf.merge(df, how="left", on="node")
+
+        if sparse:
+            col_names = [col for col in list(df) if col != "node"]
+            gdf = gdf.dropna(subset=col_names, how="all")
+            gdf = gdf.dropna(axis="columns", how="all")
+        else:
+            gdf = gdf.drop_duplicates(subset=["node"])
+            gdf = gdf.reset_index(drop=True)
+
+        return gdf
+
     def _make_graph(self):
         # get all segments and their outseg
         graph = {}
@@ -1895,8 +1932,10 @@ def export_linkages(self, f, **kwargs):
         reaches can be used to filter for the longest connections in a GIS.
 
         """
-        from ..export.shapefile_utils import recarray2shp
         from ..utils.geometry import LineString
+        from ..utils.utl_import import import_optional_dependency
+
+        gpd = import_optional_dependency("geopandas")
 
         rd = self.reach_data.copy()
         m = self.parent
@@ -1922,11 +1961,16 @@ def export_linkages(self, f, **kwargs):
             lengths.append(np.sqrt((x1 - x0) ** 2 + (y1 - y0) ** 2))
         lengths = np.array(lengths)
 
-        # append connection lengths for filtering in GIS
-        rd = recfunctions.append_fields(
-            rd, names=["length"], data=[lengths], usemask=False, asrecarray=True
-        )
-        recarray2shp(rd, geoms, f, **kwargs)
+        gdf = gpd.GeoDataFrame(data=rd, geometry=geoms)
+        gdf["length"] = lengths
+        crs = kwargs.pop("crs", None)
+        if crs is None:
+            try:
+                crs = self.parent.modelgrid.crs
+            except AttributeError:
+                pass
+        gdf = gdf.set_crs(crs)
+        gdf.to_file(f)
 
     def export_outlets(self, f, **kwargs):
         """
@@ -1934,8 +1978,10 @@ def export_outlets(self, f, **kwargs):
         the model (outset=0).
 
         """
-        from ..export.shapefile_utils import recarray2shp
         from ..utils.geometry import Point
+        from ..utils.utl_import import import_optional_dependency
+
+        gpd = import_optional_dependency("geopandas")
 
         rd = self.reach_data
         if np.min(rd.outreach) == np.max(rd.outreach):
@@ -1949,7 +1995,15 @@ def export_outlets(self, f, **kwargs):
         x0 = mg.xcellcenters[rd.i, rd.j]
         y0 = mg.ycellcenters[rd.i, rd.j]
         geoms = [Point(x, y) for x, y in zip(x0, y0)]
-        recarray2shp(rd, geoms, f, **kwargs)
+        gdf = gpd.GeoDataFrame(data=rd, geometry=geoms)
+        crs = kwargs.pop("crs", None)
+        if crs is None:
+            try:
+                crs = self.parent.modelgrid.crs
+            except AttributeError:
+                pass
+        gdf = gdf.set_crs(crs)
+        gdf.to_file(f)
 
     def export_transient_variable(self, f, varname, **kwargs):
         """
@@ -1968,8 +2022,10 @@ def export_transient_variable(self, f, varname, **kwargs):
             Variable in SFR Package dataset 6a (see SFR package documentation)
 
         """
-        from ..export.shapefile_utils import recarray2shp
         from ..utils.geometry import Point
+        from ..utils.utl_import import import_optional_dependency
+
+        gpd = import_optional_dependency("geopandas")
 
         rd = self.reach_data
         if np.min(rd.outreach) == np.max(rd.outreach):
@@ -1982,7 +2038,15 @@ def export_transient_variable(self, f, varname, **kwargs):
         x0 = mg.xcellcenters[ra.i, ra.j]
         y0 = mg.ycellcenters[ra.i, ra.j]
         geoms = [Point(x, y) for x, y in zip(x0, y0)]
-        recarray2shp(ra, geoms, f, **kwargs)
+        gdf = gpd.GeoDataFrame(data=ra, geometry=geoms)
+        crs = kwargs.pop("crs", None)
+        if crs is None:
+            try:
+                crs = self.parent.modelgrid.crs
+            except AttributeError:
+                pass
+        gdf = gdf.set_crs(crs)
+        gdf.to_file(f)
 
     @staticmethod
     def _ftype():
diff --git a/flopy/pakbase.py b/flopy/pakbase.py
index 71a7f38e69..2dd06af72e 100644
--- a/flopy/pakbase.py
+++ b/flopy/pakbase.py
@@ -674,6 +674,67 @@ def export(self, f, **kwargs):
 
         return export.utils.package_export(f, self, **kwargs)
 
+    def to_geodataframe(
+        self, gdf=None, kper=0, sparse=False, truncate_attrs=False, **kwargs
+    ):
+        """
+        Method to create a GeoDataFrame from a modflow package
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional geopandas geodataframe object to add data to. Default is None
+        kper : int
+            stress period to get transient data from
+        sparse : bool
+            optional parameter to create a sparse geodataframe
+        truncate_attrs : bool
+            method to truncate attribute names for shapefile attribute name length
+            restrictions
+
+        Returns
+        -------
+            gdf : GeoDataFrame
+        """
+        from .mbase import BaseModel
+
+        if gdf is None:
+            if isinstance(self.parent, BaseModel):
+                modelgrid = self.parent.modelgrid
+                if modelgrid is not None:
+                    gdf = modelgrid.to_geodataframe()
+                else:
+                    raise AttributeError(
+                        "model does not have a grid instance, "
+                        "please supply a geodataframe"
+                    )
+            else:
+                raise AssertionError(
+                    "Package does not have a model instance, "
+                    "please supply a geodataframe"
+                )
+
+        for attr, value in self.__dict__.items():
+            if callable(getattr(value, "to_geodataframe", None)):
+                if isinstance(value, (BaseModel, PackageInterface)):
+                    continue
+                # do not pass sparse in here, make sparse after all data has been
+                #  added to geodataframe
+                gdf = value.to_geodataframe(
+                    gdf,
+                    forgive=True,
+                    kper=kper,
+                    sparse=False,
+                    truncate_attrs=truncate_attrs,
+                )
+
+        if sparse:
+            col_names = [i for i in gdf if i not in ("geometry", "node", "row", "col")]
+            gdf = gdf.dropna(subset=col_names, how="all")
+            gdf = gdf.dropna(axis="columns", how="all")
+
+        return gdf
+
     def _generate_heading(self):
         """Generate heading."""
         from . import __version__
diff --git a/flopy/plot/plotutil.py b/flopy/plot/plotutil.py
index e8085cc8f2..d6ae28f694 100644
--- a/flopy/plot/plotutil.py
+++ b/flopy/plot/plotutil.py
@@ -2947,3 +2947,71 @@ def to_prt_pathlines(
         return df
     else:
         return ret
+
+
+def hfb_data_to_linework(recarray, modelgrid):
+    """
+    Method to get lines representing horizontal flow barriers
+
+    Parameters
+    ----------
+    recarray : np.recarray
+        recarray of hfb input data
+    modelgrid : Grid
+        modelgrid instance
+
+    Returns
+    -------
+        list : list of line segments
+    """
+    iverts = modelgrid.iverts
+    verts = modelgrid.verts
+    nodes = []
+    if modelgrid.grid_type == "structured":
+        if "k" in recarray.dtype.names:
+            for rec in recarray:
+                node1 = modelgrid.get_node([(0, rec["irow1"], rec["icol1"])])[0]
+                node2 = modelgrid.get_node([(0, rec["irow2"], rec["icol2"])])[0]
+                nodes.append((node1, node2))
+        else:
+            for rec in recarray:
+                node1 = modelgrid.get_node([(0,) + rec["cellid1"][1:]])[0]
+                node2 = modelgrid.get_node([(0,) + rec["cellid2"][1:]])[0]
+                nodes.append((node1, node2))
+
+    elif modelgrid.grid_type == "vertex":
+        for rec in recarray:
+            nodes.append((rec["cellid1"][-1], rec["cellid2"][-1]))
+
+    else:
+        if "node1" in recarray.dtype.names:
+            nodes = list(zip(recarray["node1"], recarray["node2"]))
+        else:
+            for rec in recarray:
+                nodes.append((rec["cellid1"][0], rec["cellid2"][0]))
+
+    shared_edges = []
+    for node0, node1 in nodes:
+        iv0 = iverts[node0]
+        iv1 = iverts[node1]
+        edges = []
+        for ix in range(len(iv0)):
+            edges.append(tuple(sorted((iv0[ix - 1], iv0[ix]))))
+
+        for ix in range(len(iv1)):
+            edge = tuple(sorted((iv1[ix - 1], iv1[ix])))
+            if edge in edges:
+                shared_edges.append(edge)
+                break
+
+        if not shared_edges:
+            raise AssertionError(
+                f"No shared cell edges found. Cannot represent HFB "
+                f"for nodes {node0} and {node1}"
+            )
+
+    lines = []
+    for edge in shared_edges:
+        lines.append((tuple(verts[edge[0]]), tuple(verts[edge[1]])))
+
+    return lines
diff --git a/flopy/utils/binaryfile/__init__.py b/flopy/utils/binaryfile/__init__.py
index 7d7ab8e8a2..7b4875318b 100644
--- a/flopy/utils/binaryfile/__init__.py
+++ b/flopy/utils/binaryfile/__init__.py
@@ -456,7 +456,7 @@ def get_ts(self, idx):
         result = self._init_result(nstation)
 
         # Determine grid type
-        grid_type = "structured" if self.mg is None else self.mg.grid_type
+        grid_type = "structured" if self.modelgrid is None else self.modelgrid.grid_type
 
         istat = 1
         for item in kijlist:
@@ -1641,6 +1641,56 @@ def get_data(
             if t
         ]
 
+    def to_geodataframe(
+        self,
+        gdf=None,
+        modelgrid=None,
+        idx=None,
+        kstpkper=None,
+        totim=None,
+        text=None,
+    ):
+        if (idx, kstpkper, totim) == (None, None, None):
+            raise AssertionError(
+                "to_geodataframe() missing 1 required argument: "
+                "please provide 'idx', 'kstpkper', or 'totim'"
+            )
+
+        if gdf is None:
+            if modelgrid is None:
+                if self.modelgrid is None:
+                    raise AssertionError(
+                        "A geodataframe or modelgrid instance must be supplied"
+                    )
+                modelgrid = self.modelgrid
+
+            gdf = modelgrid.to_geodataframe()
+
+        col_names = []
+        if text is None:
+            textlist = [i.decode() for i in self.textlist]
+            for text in textlist:
+                data = self.get_data(
+                    idx=idx, kstpkper=kstpkper, totim=totim, text=text, full3D=True
+                )
+
+                text = text.strip().lower().replace(" ", "_")
+                for ix, arr in enumerate(data[0]):
+                    name = f"{text}_{ix}"
+                    gdf[name] = arr.ravel()
+                    col_names.append(name)
+        else:
+            data = self.get_data(
+                idx=idx, kstpkper=kstpkper, totim=totim, text=text, full3D=True
+            )
+            text = text.strip().lower().replace(" ", "_")
+            for ix, arr in enumerate(data[0]):
+                name = f"{text}_{ix}"
+                gdf[name] = arr.ravel()
+                col_names.append(name)
+
+        return gdf
+
     def get_ts(self, idx, text=None, times=None, variable="q"):
         """
         Get a time series from the binary budget file.
diff --git a/flopy/utils/datafile.py b/flopy/utils/datafile.py
index ce803bdd4d..cf4960fb15 100644
--- a/flopy/utils/datafile.py
+++ b/flopy/utils/datafile.py
@@ -197,18 +197,18 @@ def __init__(self, filename: Union[str, PathLike], precision, verbose, **kwargs)
 
         self.model = None
         self.dis = None
-        self.mg = None
+        self.modelgrid = None
         if "model" in kwargs.keys():
             self.model = kwargs.pop("model")
-            self.mg = self.model.modelgrid
+            self.modelgrid = self.model.modelgrid
             self.dis = self.model.dis
         if "dis" in kwargs.keys():
             self.dis = kwargs.pop("dis")
-            self.mg = self.dis.parent.modelgrid
+            self.modelgrid = self.dis.parent.modelgrid
         if "tdis" in kwargs.keys():
             self.tdis = kwargs.pop("tdis")
         if "modelgrid" in kwargs.keys():
-            self.mg = kwargs.pop("modelgrid")
+            self.modelgrid = kwargs.pop("modelgrid")
         if len(kwargs.keys()) > 0:
             args = ",".join(kwargs.keys())
             raise ValueError(f"LayerFile error: unrecognized kwargs: {args}")
@@ -217,9 +217,9 @@ def __init__(self, filename: Union[str, PathLike], precision, verbose, **kwargs)
         self._build_index()
 
         # now that we read the data and know nrow and ncol,
-        # we can make a generic mg if needed
-        if self.mg is None:
-            self.mg = StructuredGrid(
+        # we can make a generic modelgrid if needed
+        if self.modelgrid is None:
+            self.modelgrid = StructuredGrid(
                 delc=np.ones((self.nrow,)),
                 delr=np.ones(self.ncol),
                 nlay=self.nlay,
@@ -240,6 +240,54 @@ def __enter__(self):
     def __exit__(self, *exc):
         self.close()
 
+    def to_geodataframe(
+        self, gdf=None, modelgrid=None, kstpkper=None, totim=None, attrib_name=None
+    ):
+        """
+        Generate a GeoDataFrame with data from a LayerFile instance
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional, existing geodataframe with NCPL geometries
+        modelgrid : Grid
+            optional modelgrid instance to generate a GeoDataFrame from
+        kstpkper : tuple of ints
+            A tuple containing the time step and stress period (kstp, kper).
+            These are zero-based kstp and kper values.
+        totim : float
+            The simulation time.
+        attrib_name : str
+            optional base name of attribute columns. (default is text attribute)
+
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        if gdf is None:
+            if modelgrid is None:
+                if self.modelgrid is None:
+                    raise AssertionError(
+                        "A geodataframe or modelgrid instance must be supplied"
+                    )
+                modelgrid = self.modelgrid
+
+            gdf = modelgrid.to_geodataframe()
+
+        array = np.atleast_3d(
+            self.get_data(kstpkper=kstpkper, totim=totim).transpose()
+        ).transpose()
+
+        if attrib_name is None:
+            attrib_name = self.text.decode()
+
+        for ix, arr in enumerate(array):
+            name = f"{attrib_name}_{ix}"
+            gdf[name] = np.ravel(arr)
+
+        return gdf
+
     def to_shapefile(
         self,
         filename: Union[str, PathLike],
@@ -287,7 +335,10 @@ def to_shapefile(
         >>> times = hdobj.get_times()
         >>> hdobj.to_shapefile('test_heads_sp6.shp', totim=times[-1])
         """
-
+        warnings.warn(
+            "to_shapefile() is deprecated and is being replaced by to_geodataframe()",
+            DeprecationWarning,
+        )
         plotarray = np.atleast_3d(
             self.get_data(kstpkper=kstpkper, totim=totim, mflay=mflay).transpose()
         ).transpose()
@@ -301,7 +352,7 @@ def to_shapefile(
 
         from ..export.shapefile_utils import write_grid_shapefile
 
-        write_grid_shapefile(filename, self.mg, attrib_dict, verbose)
+        write_grid_shapefile(filename, self.modelgrid, attrib_dict, verbose)
 
     def plot(
         self,
@@ -413,7 +464,7 @@ def plot(
             axes=axes,
             filenames=filenames,
             mflay=mflay,
-            modelgrid=self.mg,
+            modelgrid=self.modelgrid,
             **kwargs,
         )
 
@@ -631,7 +682,7 @@ def _build_kijlist(self, idx):
         - DISU: list of integers (node)
         """
         # Determine grid type
-        grid_type = "structured" if self.mg is None else self.mg.grid_type
+        grid_type = "structured" if self.modelgrid is None else self.modelgrid.grid_type
 
         # Normalize idx to a list
         if isinstance(idx, int):
@@ -681,9 +732,9 @@ def _build_kijlist(self, idx):
                     )
                 if k < 0 or k >= self.nlay:
                     raise ValueError(f"Layer index {k} out of range [0, {self.nlay})")
-                if cell < 0 or cell >= self.mg.ncpl:
+                if cell < 0 or cell >= self.modelgrid.ncpl:
                     raise ValueError(
-                        f"Cell index {cell} out of range [0, {self.mg.ncpl})"
+                        f"Cell index {cell} out of range [0, {self.modelgrid.ncpl})"
                     )
                 # Store as 2-tuple for DISV
                 kijlist.append((k, cell))
@@ -708,9 +759,9 @@ def _build_kijlist(self, idx):
                         f"DISU unstructured grid requires integer node index "
                         f"or 3-tuple (dummy, dummy, node), got: {item}"
                     )
-                if node < 0 or node >= self.mg.nnodes:
+                if node < 0 or node >= self.modelgrid.nnodes:
                     raise ValueError(
-                        f"Node index {node} out of range [0, {self.mg.nnodes})"
+                        f"Node index {node} out of range [0, {self.modelgrid.nnodes})"
                     )
                 # Store as integer for DISU
                 kijlist.append(node)
diff --git a/flopy/utils/geospatial_utils.py b/flopy/utils/geospatial_utils.py
index a60a16aca5..36cfde7b1e 100644
--- a/flopy/utils/geospatial_utils.py
+++ b/flopy/utils/geospatial_utils.py
@@ -453,7 +453,7 @@ def shapely(self):
         return self._shapely
 
     @property
-    def geo_dataframe(self):
+    def geodataframe(self):
         """
         Property that returns a geopandas DataFrame
 
diff --git a/flopy/utils/gridgen.py b/flopy/utils/gridgen.py
index efe66c82c8..1857e8457f 100644
--- a/flopy/utils/gridgen.py
+++ b/flopy/utils/gridgen.py
@@ -465,6 +465,8 @@ def build(self, verbose=False):
         None
 
         """
+        gpd = import_optional_dependency("geopandas")
+
         fname = os.path.join(self.model_ws, "_gridgen_build.dfn")
         f = open(fname, "w")
 
@@ -507,8 +509,9 @@ def build(self, verbose=False):
         )
 
         # Create a recarray of the grid polygon shapefile
-        shapename = os.path.join(self.model_ws, "qtgrid")
-        self.qtra = shp2recarray(shapename)
+        shapename = os.path.join(self.model_ws, "qtgrid.shp")
+        gdf = gpd.read_file(shapename)
+        self.qtra = gdf.to_records()
 
     def get_vertices(self, nodenumber):
         """
diff --git a/flopy/utils/modpathfile.py b/flopy/utils/modpathfile.py
index a5bc077886..5a34175d5f 100644
--- a/flopy/utils/modpathfile.py
+++ b/flopy/utils/modpathfile.py
@@ -12,6 +12,7 @@
 from flopy.utils.particletrackfile import ParticleTrackFile
 
 from ..utils.flopy_io import loadtxt
+from ..utils.utl_import import import_optional_dependency
 
 
 class ModpathFile(ParticleTrackFile):
@@ -671,6 +672,61 @@ def get_destination_endpoint_data(self, dest_cells, source=False):
         inds = np.isin(raslice, dest_cells)
         return data[inds].copy().view(np.recarray)
 
+    def to_geodataframe(
+        self,
+        modelgrid,
+        data=None,
+        direction="ending",
+    ):
+        """
+        Create a geodataframe of particle starting / ending locations.
+
+        Parameters
+        ----------
+        modelgrid : flopy.discretization.grid instance
+            Used to scale and rotate Global x,y,z values in MODPATH Endpoint
+            file.
+        data : np.recarray
+            Record array of same form as that returned by EndpointFile.get_alldata.
+            (if none, EndpointFile.get_alldata() is exported).
+        direction : str
+            String defining if 'starting' or 'ending' particle locations should be
+            considered. (default is 'ending')
+        """
+        from ..utils import geometry
+
+        gpd = import_optional_dependency("geopandas")
+        shapely_geo = import_optional_dependency("shapely.geometry")
+        if data is None:
+            data = self.get_alldata()
+
+        if direction.lower() == "ending":
+            xcol, ycol, zcol = "x", "y", "z"
+        elif direction.lower() == "starting":
+            xcol, ycol, zcol = "x0", "y0", "z0"
+        else:
+            raise Exception(
+                'flopy.map.plot_endpoint direction must be "ending" or "starting".'
+            )
+        x, y = geometry.transform(
+            data[xcol],
+            data[ycol],
+            xoff=modelgrid.xoffset,
+            yoff=modelgrid.yoffset,
+            angrot_radians=modelgrid.angrot_radians,
+        )
+        z = data[zcol]
+
+        geoms = [shapely_geo.Point(p) for p in zip(x, y, z)]
+        gdf = gpd.GeoDataFrame(data, geometry=geoms, crs=modelgrid.crs)
+
+        # adjust to 1 based node numbers
+        for col in list(gdf):
+            if col in self.kijnames:
+                gdf[col] += 1
+
+        return gdf
+
     def write_shapefile(
         self,
         data=None,
@@ -715,44 +771,21 @@ def write_shapefile(
                 - ``epsg`` (int): use ``crs`` instead.
 
         """
-        from ..discretization import StructuredGrid
-        from ..export.shapefile_utils import recarray2shp
-        from ..utils import geometry
-        from ..utils.geometry import Point
+        import warnings
 
+        warnings.warn(
+            "write_shapefile is Deprecated, please use to_geodataframe() in the future"
+        )
         epd = (data if data is not None else endpoint_data).copy()
-        if epd is None:
-            epd = self.get_alldata()
+        gdf = self.to_geodataframe(modelgrid=mg, data=epd, direction=direction)
 
-        if direction.lower() == "ending":
-            xcol, ycol, zcol = "x", "y", "z"
-        elif direction.lower() == "starting":
-            xcol, ycol, zcol = "x0", "y0", "z0"
-        else:
-            raise Exception(
-                'flopy.map.plot_endpoint direction must be "ending" or "starting".'
-            )
-        if mg is None:
-            raise ValueError("A modelgrid object was not provided.")
-
-        if isinstance(mg, StructuredGrid):
-            x, y = geometry.transform(
-                epd[xcol],
-                epd[ycol],
-                xoff=mg.xoffset,
-                yoff=mg.yoffset,
-                angrot_radians=mg.angrot_radians,
-            )
-        else:
-            x, y = mg.get_coords(epd[xcol], epd[ycol])
-        z = epd[zcol]
+        if crs is not None:
+            if gdf.crs is None:
+                gdf = gdf.set_crs(crs)
+            else:
+                gdf = gdf.to_crs(crs)
 
-        geoms = [Point(x[i], y[i], z[i]) for i in range(len(epd))]
-        # convert back to one-based
-        for n in self.kijnames:
-            if n in epd.dtype.names:
-                epd[n] += 1
-        recarray2shp(epd, geoms, shpname=shpname, crs=crs, **kwargs)
+        gdf.to_file(shpname)
 
 
 class TimeseriesFile(ModpathFile):
diff --git a/flopy/utils/particletrackfile.py b/flopy/utils/particletrackfile.py
index 00e9a7a141..6a7c838993 100644
--- a/flopy/utils/particletrackfile.py
+++ b/flopy/utils/particletrackfile.py
@@ -10,6 +10,8 @@
 import numpy as np
 from numpy.lib.recfunctions import stack_arrays
 
+from .utl_import import import_optional_dependency
+
 MIN_PARTICLE_TRACK_DTYPE = np.dtype(
     [
         ("x", np.float32),
@@ -183,6 +185,125 @@ def intersect(self, cells, to_recarray) -> np.recarray:
         """Find intersection of pathlines with cells."""
         pass
 
+    def to_geodataframe(
+        self,
+        modelgrid,
+        data=None,
+        one_per_particle=True,
+        direction="ending",
+    ):
+        """
+        Create a geodataframe of particle tracks.
+
+        Parameters
+        ----------
+        modelgrid : flopy.discretization.Grid instance
+            Used to scale and rotate Global x, y, z values.
+        data : np.recarray
+            Record array of same form as that returned by
+            get_alldata(). (if none, get_alldata() is exported).
+        one_per_particle : boolean (default True)
+            True writes a single LineString with a single set of attribute
+            data for each particle. False writes a record/geometry for each
+            pathline segment (each row in the PathLine file). This option can
+            be used to visualize attribute information (time, model layer,
+            etc.) across a pathline in a GIS.
+        direction : str
+            String defining if starting or ending particle locations should be
+            included in shapefile attribute information. Only used if
+            one_per_particle=False. (default is 'ending')
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        from . import geometry
+
+        shapely_geo = import_optional_dependency("shapely.geometry")
+        gpd = import_optional_dependency("geopandas")
+
+        if data is None:
+            data = self._data.view(np.recarray)
+        else:
+            # convert pathline list to a single recarray
+            if isinstance(data, list):
+                s = data[0]
+                print(s.dtype)
+                for n in range(1, len(data)):
+                    s = stack_arrays((s, data[n]))
+                data = s.view(np.recarray)
+
+        data = data.copy()
+        data.sort(order=["particleid", "time"])
+
+        particles = np.unique(data.particleid)
+        geoms = []
+
+        # create a dict of attrs?
+        headings = ["particleid", "particlegroup", "time", "k", "i", "j", "node"]
+        attrs = []
+        for h in headings:
+            if h in data.dtype.names:
+                attrs.append(h)
+
+        if one_per_particle:
+            dfdata = {a: [] for a in attrs}
+            if direction == "ending":
+                idx = -1
+            else:
+                idx = 0
+
+            for p in particles:
+                ra = data[data.particleid == p]
+                for k, _ in dfdata.items():
+                    if k == "time":
+                        dfdata[k].append(np.max(ra[k]))
+                    else:
+                        dfdata[k].append(ra[k][idx])
+
+                x, y = geometry.transform(
+                    ra.x,
+                    ra.y,
+                    modelgrid.xoffset,
+                    modelgrid.yoffset,
+                    modelgrid.angrot_radians,
+                )
+                z = ra.z
+
+                line = list(zip(x, y, z))
+                geoms.append(shapely_geo.LineString(line))
+
+        else:
+            dfdata = {a: [] for a in data.dtype.names}
+            for pid in particles:
+                ra = data[data.particleid == pid]
+                x, y = geometry.transform(
+                    ra.x,
+                    ra.y,
+                    modelgrid.xoffset,
+                    modelgrid.yoffset,
+                    modelgrid.angrot_radians,
+                )
+                z = ra.z
+                geoms += [
+                    shapely_geo.LineString(
+                        [(x[i - 1], y[i - 1], z[i - 1]), (x[i], y[i], z[i])]
+                    )
+                    for i in np.arange(1, (len(ra)))
+                ]
+                for k in dfdata.keys():
+                    dfdata[k].extend(ra[k][1:])
+
+        # now create a geodataframe
+        gdf = gpd.GeoDataFrame(dfdata, geometry=geoms, crs=modelgrid.crs)
+
+        # adjust to 1 based node numbers
+        for col in list(gdf):
+            if col in self.kijnames:
+                gdf[col] += 1
+
+        return gdf
+
     def write_shapefile(
         self,
         data=None,
@@ -227,102 +348,22 @@ def write_shapefile(
                - ``epsg`` (int): use ``crs`` instead.
 
         """
-        from ..export.shapefile_utils import recarray2shp
-        from . import geometry
-        from .geometry import LineString
-
-        series = data
-        if series is None:
-            series = self._data.view(np.recarray)
-        else:
-            # convert pathline list to a single recarray
-            if isinstance(series, list):
-                s = series[0]
-                print(s.dtype)
-                for n in range(1, len(series)):
-                    s = stack_arrays((s, series[n]))
-                series = s.view(np.recarray)
-
-        series = series.copy()
-        series.sort(order=["particleid", "time"])
+        import warnings
 
-        if mg is None:
-            raise ValueError("A modelgrid object was not provided.")
-
-        particles = np.unique(series.particleid)
-        geoms = []
-
-        # create dtype with select attributes in pth
-        names = series.dtype.names
-        dtype = []
-        atts = ["particleid", "particlegroup", "time", "k", "i", "j", "node"]
-        for att in atts:
-            if att in names:
-                t = np.int32
-                if att == "time":
-                    t = np.float32
-                dtype.append((att, t))
-        dtype = np.dtype(dtype)
-
-        # reset names to the selected names in the created dtype
-        names = dtype.names
-
-        # 1 geometry for each path
-        if one_per_particle:
-            loc_inds = 0
-            if direction == "ending":
-                loc_inds = -1
-
-            sdata = []
-            for pid in particles:
-                ra = series[series.particleid == pid]
-
-                x, y = geometry.transform(
-                    ra.x, ra.y, mg.xoffset, mg.yoffset, mg.angrot_radians
-                )
-                z = ra.z
-                geoms.append(LineString(list(zip(x, y, z))))
-
-                t = [pid]
-                if "particlegroup" in names:
-                    t.append(ra.particlegroup[0])
-                t.append(ra.time.max())
-                if "k" in names:
-                    t.append(ra.k[loc_inds])
-                if "node" in names:
-                    t.append(ra.node[loc_inds])
-                else:
-                    if "i" in names:
-                        t.append(ra.i[loc_inds])
-                    if "j" in names:
-                        t.append(ra.j[loc_inds])
-                sdata.append(tuple(t))
-
-            sdata = np.array(sdata, dtype=dtype).view(np.recarray)
-
-        # geometry for each row in PathLine file
-        else:
-            dtype = series.dtype
-            sdata = []
-            for pid in particles:
-                ra = series[series.particleid == pid]
-                if mg is not None:
-                    x, y = geometry.transform(
-                        ra.x, ra.y, mg.xoffset, mg.yoffset, mg.angrot_radians
-                    )
-                else:
-                    x, y = geometry.transform(ra.x, ra.y, 0, 0, 0)
-                z = ra.z
-                geoms += [
-                    LineString([(x[i - 1], y[i - 1], z[i - 1]), (x[i], y[i], z[i])])
-                    for i in np.arange(1, (len(ra)))
-                ]
-                sdata += ra[1:].tolist()
-            sdata = np.array(sdata, dtype=dtype).view(np.recarray)
-
-        # convert back to one-based
-        for n in set(self.kijnames).intersection(set(sdata.dtype.names)):
-            sdata[n] += 1
+        warnings.warn(
+            "write_shapefile will be Deprecated, please use to_geo_dataframe()",
+            DeprecationWarning,
+        )
+        gdf = self.to_geodataframe(
+            modelgrid=mg,
+            data=data,
+            one_per_particle=one_per_particle,
+            direction=direction,
+        )
+        if crs is not None:
+            if gdf.crs is None:
+                gdf = gdf.set_crs(crs)
+            else:
+                gdf = gdf.to_crs(crs)
 
-        # write the final recarray to a shapefile
-        recarray2shp(sdata, geoms, shpname=shpname, crs=crs, **kwargs)
+        gdf.to_file(shpname)
diff --git a/flopy/utils/util_array.py b/flopy/utils/util_array.py
index 3823260615..afef0d8567 100644
--- a/flopy/utils/util_array.py
+++ b/flopy/utils/util_array.py
@@ -634,6 +634,55 @@ def data_type(self):
     def plottable(self):
         return True
 
+    def to_geodataframe(
+        self, gdf=None, forgive=False, name=None, truncate_attrs=False, **kwargs
+    ):
+        """
+        Method to add data to a GeoDataFrame for exporting as a geospatial file
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame instance. If GeoDataFrame is None, one will be
+            constructed from modelgrid information
+        forgive : bool
+            optional flag to continue running and pass data that is not compatible
+            with the geodataframe shape
+        name : str
+            optional array name base. If None, method uses the .name attribute
+        truncate_attrs : bool
+            method to truncate attribute names for shapefile restrictions
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        if self.model is None:
+            return gdf
+        else:
+            modelgrid = self.model.modelgrid
+            if modelgrid is None:
+                return gdf
+
+            if gdf is None:
+                gdf = modelgrid.to_geodataframe()
+
+            if name is not None:
+                names = [name for _ in range(len(self.util_2ds))]
+            else:
+                names = self.name
+            for lay, u2d in enumerate(self.util_2ds):
+                name = f"{names[lay]}_{lay}"
+                gdf = u2d.to_geodataframe(
+                    gdf=gdf,
+                    name=name,
+                    forgive=forgive,
+                    truncate_attrs=truncate_attrs,
+                    **kwargs,
+                )
+
+            return gdf
+
     def export(self, f, **kwargs):
         from .. import export
 
@@ -1059,6 +1108,37 @@ def data_type(self):
     def plottable(self):
         return False
 
+    def to_geodataframe(
+        self, gdf=None, kper=0, forgive=False, truncate_attrs=False, **kwargs
+    ):
+        """
+        Method to add data to a GeoDataFrame for exporting as a geospatial file
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame instance. If GeoDataFrame is None, one will be
+            constructed from modelgrid information
+        kper : int
+            stress period to export
+        forgive : bool
+            boolean flag for sparse dataframe construction. Default is False
+        truncate_attrs : bool
+            method to truncate attribute names for shapefile restrictions
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        u3d = self.transient_3ds[kper]
+        # note: may need to provide a pass through name for u3d to avoid s.p.
+        # number being tacked on. Test this on a model with the LAK package...
+        name = self.name_base[:-1].lower()
+        gdf = u3d.to_geodataframe(
+            gdf=gdf, forgive=forgive, name=name, truncate_attrs=truncate_attrs, **kwargs
+        )
+        return gdf
+
     def get_zero_3d(self, kper):
         name = f"{self.name_base}{kper + 1}(filled zero)"
         return Util3d(
@@ -1399,6 +1479,35 @@ def from_4d(cls, model, pak_name, m4ds):
             name=name,
         )
 
+    def to_geodataframe(
+        self, gdf=None, kper=0, forgive=False, truncate_attrs=False, **kwargs
+    ):
+        """
+        Method to add data to a GeoDataFrame for exporting as a geospatial file
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame instance. If GeoDataFrame is None, one will be
+            constructed from modelgrid information
+        kper : int
+            stress period to export
+        forgive : bool
+            boolean flag for sparse dataframe construction. Default is False
+        truncate_attrs : bool
+            method to truncate attribute names for shapefile restrictions
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        u2d = self.transient_2ds[kper]
+        name = self.name_base[:-1]
+        gdf = u2d.to_geodataframe(
+            gdf=gdf, name=name, forgive=forgive, truncate_attrs=truncate_attrs, **kwargs
+        )
+        return gdf
+
     def __setattr__(self, key, value):
         if hasattr(self, "transient_2ds") and key == "cnstnt":
             # set cnstnt for each u2d
@@ -1876,6 +1985,65 @@ def _decide_how(self):
         else:
             self._how = "internal"
 
+    def to_geodataframe(
+        self, gdf=None, name=None, forgive=False, truncate_attrs=False, **kwargs
+    ):
+        """
+        Method to add an input array to a geopandas GeoDataFrame
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame object
+        name : str
+            optional attribute name, default uses util2d name
+        forgive : bool
+            optional flag to continue if data shape not compatible with GeoDataFrame
+        truncate_attrs : bool
+            method to truncate attribute names for shapefile restrictions
+
+        Returns
+        -------
+            geopandas GeoDataFrame
+        """
+        from ..export.shapefile_utils import shape_attr_name
+
+        if self.model is None:
+            return gdf
+        else:
+            modelgrid = self.model.modelgrid
+            if gdf is None:
+                if modelgrid is None:
+                    return gdf
+                gdf = modelgrid.to_geodataframe()
+
+            if modelgrid is not None:
+                if modelgrid.grid_type != "unstructured":
+                    ncpl = modelgrid.ncpl
+                else:
+                    ncpl = modelgrid.nnodes
+            else:
+                ncpl = len(gdf)
+
+            if name is None:
+                name = self.name
+
+            if truncate_attrs:
+                name = shape_attr_name(name, keep_layer=True)
+
+            data = self.array
+
+            if data.size == ncpl:
+                gdf[name] = data.ravel()
+            elif forgive:
+                return gdf
+            else:
+                raise AssertionError(
+                    f"Data size {data.size} not compatible with dataframe length {ncpl}"
+                )
+
+            return gdf
+
     def plot(
         self,
         title=None,
diff --git a/flopy/utils/util_list.py b/flopy/utils/util_list.py
index 3ddbb550fa..fb89228d16 100644
--- a/flopy/utils/util_list.py
+++ b/flopy/utils/util_list.py
@@ -124,6 +124,76 @@ def get_empty(self, ncell=0):
         d = create_empty_recarray(ncell, self.dtype, default_value=-1.0e10)
         return d
 
+    def to_geodataframe(
+        self, gdf=None, kper=0, sparse=False, truncate_attrs=False, **kwargs
+    ):
+        """
+        Method to add data to a GeoDataFrame for exporting as a geospatial file
+
+        Parameters
+        ----------
+        gdf : GeoDataFrame
+            optional GeoDataFrame instance. If GeoDataFrame is None, one will be
+            constructed from modelgrid information
+        kper : int
+            stress period to export
+        sparse : bool
+            boolean flag for sparse dataframe construction. Default is False
+        truncate_attrs : bool
+            method to truncate attribute names for shapefile restrictions
+
+        Returns
+        -------
+            GeoDataFrame
+        """
+        from ..export.shapefile_utils import shape_attr_name
+
+        if self.model is None:
+            return gdf
+        else:
+            modelgrid = self.model.modelgrid
+            if modelgrid is None:
+                return gdf
+
+            if gdf is None:
+                gdf = modelgrid.to_geodataframe()
+
+            data = self.array
+
+            col_names = []
+            for name, array4d in data.items():
+                if truncate_attrs:
+                    name = shape_attr_name(name, length=4)
+                else:
+                    name = f"{self.name[0].lower()}_{name}"
+                array = array4d[kper]
+                if modelgrid.grid_type == "unstructured":
+                    array = array.ravel()
+                    aname = name
+                    if truncate_attrs:
+                        aname = f"{name}{kper}"
+                    gdf[aname] = array
+                    col_names.append(aname)
+                else:
+                    if modelgrid.nlay is None:
+                        nlay = 1
+                    else:
+                        nlay = modelgrid.nlay
+                    for lay in range(nlay):
+                        arr = array[lay].ravel()
+                        if truncate_attrs:
+                            aname = f"{name}{lay}{kper}"
+                        else:
+                            aname = f"{name}_{lay}_{kper}"
+                        gdf[aname] = arr.ravel()
+                        col_names.append(aname)
+
+            if sparse:
+                gdf = gdf.dropna(subset=col_names, how="all")
+                gdf = gdf.dropna(axis="columns", how="all")
+
+            return gdf
+
     def export(self, f, **kwargs):
         from .. import export