diff --git a/autotest/t074_test_geospatial_util.py b/autotest/t074_test_geospatial_util.py
index 90492d110f..0da761bcfe 100644
--- a/autotest/t074_test_geospatial_util.py
+++ b/autotest/t074_test_geospatial_util.py
@@ -3,12 +3,12 @@
     "coordinates": (
         (
             (-121.389308, 38.560816),
-            (-121.363391, 38.568835),
-            (-121.358641, 38.565972),
-            (-121.359327, 38.562767),
-            (-121.369932, 38.560575),
-            (-121.370609, 38.557232),
             (-121.385435, 38.555018),
+            (-121.370609, 38.557232),
+            (-121.369932, 38.560575),
+            (-121.359327, 38.562767),
+            (-121.358641, 38.565972),
+            (-121.363391, 38.568835),
             (-121.389308, 38.560816),
         ),
     ),
@@ -19,19 +19,19 @@
     "coordinates": (
         (
             (-121.383097, 38.565764),
-            (-121.342866, 38.579086),
-            (-121.342739, 38.578995),
-            (-121.323309, 38.578953),
-            (-121.358295, 38.561163),
-            (-121.379047, 38.559053),
             (-121.382318, 38.562934),
+            (-121.379047, 38.559053),
+            (-121.358295, 38.561163),
+            (-121.323309, 38.578953),
+            (-121.342739, 38.578995),
+            (-121.342866, 38.579086),
             (-121.383097, 38.565764),
         ),
         (
             (-121.367281, 38.567214),
-            (-121.362633, 38.562622),
-            (-121.345857, 38.570301),
             (-121.352168, 38.572258),
+            (-121.345857, 38.570301),
+            (-121.362633, 38.562622),
             (-121.367281, 38.567214),
         ),
     ),
@@ -43,17 +43,17 @@
         [
             (
                 (-121.433775, 38.544254),
-                (-121.424263, 38.547474),
                 (-121.422917, 38.540376),
+                (-121.424263, 38.547474),
                 (-121.433775, 38.544254),
             )
         ],
         [
             (
                 (-121.456113, 38.552220),
-                (-121.440092, 38.548303),
-                (-121.440053, 38.537820),
                 (-121.459991, 38.541350),
+                (-121.440053, 38.537820),
+                (-121.440092, 38.548303),
                 (-121.456113, 38.552220),
             )
         ],
@@ -113,16 +113,16 @@
 def test_import_geospatial_utils():
 
     from flopy.utils.geospatial_utils import (
-        GeoSpatialCollection,
         GeoSpatialUtil,
+        GeoSpatialCollection,
     )
 
     return
 
 
 def test_polygon():
-    from flopy.utils.geometry import Polygon, Shape
     from flopy.utils.geospatial_utils import GeoSpatialUtil
+    from flopy.utils.geometry import Shape, Polygon
 
     poly = Shape.from_geojson(polygon)
     gi1 = poly.__geo_interface__
@@ -149,12 +149,16 @@ def test_polygon():
         is_equal = gi1 == gi2
 
         if not is_equal:
-            raise AssertionError("GeoSpatialUtil polygon conversion error")
+            # pyshp < 2.2.0 sorts coordinates in opposite direction
+            gi2['coordinates'] = (gi2['coordinates'][0][::-1],)
+            is_equal = gi1 == gi2
+            if not is_equal:
+                raise AssertionError("GeoSpatialUtil polygon conversion error")
 
 
 def test_polygon_with_hole():
-    from flopy.utils.geometry import Polygon, Shape
     from flopy.utils.geospatial_utils import GeoSpatialUtil
+    from flopy.utils.geometry import Shape, Polygon
 
     poly = Shape.from_geojson(poly_w_hole)
     gi1 = poly.__geo_interface__
@@ -181,12 +185,17 @@ def test_polygon_with_hole():
         is_equal = gi1 == gi2
 
         if not is_equal:
-            raise AssertionError("GeoSpatialUtil polygon conversion error")
+            # pyshp < 2.2.0 sorts coordinates in opposite direction
+            t = reversed(t)
+            gi2 = t.__geo_interface__
+            is_equal = gi1 == gi2
+            if not is_equal:
+                raise AssertionError("GeoSpatialUtil polygon conversion error")
 
 
 def test_multipolygon():
-    from flopy.utils.geometry import MultiPolygon, Shape
     from flopy.utils.geospatial_utils import GeoSpatialUtil
+    from flopy.utils.geometry import Shape, MultiPolygon
 
     poly = Shape.from_geojson(multipolygon)
     gi1 = poly.__geo_interface__
@@ -213,14 +222,19 @@ def test_multipolygon():
         is_equal = gi1 == gi2
 
         if not is_equal:
-            raise AssertionError(
-                "GeoSpatialUtil multipolygon conversion error"
-            )
+            # pyshp < 2.2.0 sorts coordinates in opposite direction
+            t = reversed(t)
+            gi2 = t.__geo_interface__
+            is_equal = gi1 == gi2
+            if not is_equal:
+                raise AssertionError(
+                    "GeoSpatialUtil multipolygon conversion error"
+                )
 
 
 def test_point():
-    from flopy.utils.geometry import Point, Shape
     from flopy.utils.geospatial_utils import GeoSpatialUtil
+    from flopy.utils.geometry import Shape, Point
 
     pt = Shape.from_geojson(point)
     gi1 = pt.__geo_interface__
@@ -251,8 +265,8 @@ def test_point():
 
 
 def test_multipoint():
-    from flopy.utils.geometry import MultiPoint, Shape
     from flopy.utils.geospatial_utils import GeoSpatialUtil
+    from flopy.utils.geometry import Shape, MultiPoint
 
     mpt = Shape.from_geojson(multipoint)
     gi1 = mpt.__geo_interface__
@@ -283,8 +297,8 @@ def test_multipoint():
 
 
 def test_linestring():
-    from flopy.utils.geometry import LineString, Shape
     from flopy.utils.geospatial_utils import GeoSpatialUtil
+    from flopy.utils.geometry import Shape, LineString
 
     lstr = Shape.from_geojson(linestring)
     gi1 = lstr.__geo_interface__
@@ -315,8 +329,8 @@ def test_linestring():
 
 
 def test_multilinestring():
-    from flopy.utils.geometry import MultiLineString, Shape
     from flopy.utils.geospatial_utils import GeoSpatialUtil
+    from flopy.utils.geometry import Shape, MultiLineString
 
     mlstr = Shape.from_geojson(multilinestring)
     gi1 = mlstr.__geo_interface__
@@ -349,8 +363,8 @@ def test_multilinestring():
 
 
 def test_polygon_collection():
-    from flopy.utils.geometry import Collection, Shape
     from flopy.utils.geospatial_utils import GeoSpatialCollection
+    from flopy.utils.geometry import Shape, Collection
 
     col = [
         Shape.from_geojson(polygon),
@@ -375,20 +389,26 @@ def test_polygon_collection():
             continue
 
         gc2 = GeoSpatialCollection(col, shapetype)
-        gi2 = [i.flopy_geometry.__geo_interface__ for i in gc2]
 
-        for ix, gi in enumerate(gi2):
-            is_equal = gi == gi1[ix]
+        for ix, gi in enumerate(gc2):
+            t = gi.flopy_geometry
+            gi2 = t.__geo_interface__
+            is_equal = gi2 == gi1[ix]
 
             if not is_equal:
-                raise AssertionError(
-                    "GeoSpatialCollection Polygon conversion error"
-                )
+                # pyshp < 2.2.0 sorts coordinates in opposite direction
+                t = reversed(t)
+                gi2 = t.__geo_interface__
+                is_equal = gi2 == gi1[ix]
+                if not is_equal:
+                    raise AssertionError(
+                        "GeoSpatialCollection Polygon conversion error"
+                    )
 
 
 def test_point_collection():
-    from flopy.utils.geometry import Collection, Shape
     from flopy.utils.geospatial_utils import GeoSpatialCollection
+    from flopy.utils.geometry import Shape, Collection
 
     col = [Shape.from_geojson(point), Shape.from_geojson(multipoint)]
 
@@ -421,8 +441,8 @@ def test_point_collection():
 
 
 def test_linestring_collection():
-    from flopy.utils.geometry import Collection, Shape
     from flopy.utils.geospatial_utils import GeoSpatialCollection
+    from flopy.utils.geometry import Shape, Collection
 
     col = [Shape.from_geojson(linestring), Shape.from_geojson(multilinestring)]
 
@@ -455,8 +475,8 @@ def test_linestring_collection():
 
 
 def test_mixed_collection():
-    from flopy.utils.geometry import Collection, Shape
     from flopy.utils.geospatial_utils import GeoSpatialCollection
+    from flopy.utils.geometry import Shape, Collection
 
     col = [
         Shape.from_geojson(polygon),
@@ -485,13 +505,22 @@ def test_mixed_collection():
             continue
 
         gc2 = GeoSpatialCollection(col, shapetype)
-        gi2 = [i.flopy_geometry.__geo_interface__ for i in gc2]
 
-        for ix, gi in enumerate(gi2):
-            is_equal = gi == gi1[ix]
+        for ix, gi in enumerate(gc2):
+            t = gi.flopy_geometry
+            gi2 = t.__geo_interface__
+
+            is_equal = gi2 == gi1[ix]
 
             if not is_equal:
-                raise AssertionError("GeoSpatialCollection conversion error")
+                t = reversed(t)
+                gi2 = t.__geo_interface__
+                is_equal = gi2 == gi1[ix]
+
+                if not is_equal:
+                    raise AssertionError(
+                        "GeoSpatialCollection conversion error"
+                    )
 
 
 if __name__ == "__main__":
diff --git a/flopy/modflow/mfag.py b/flopy/modflow/mfag.py
index 5eae021880..786e2ec37c 100644
--- a/flopy/modflow/mfag.py
+++ b/flopy/modflow/mfag.py
@@ -360,7 +360,7 @@ def write_file(self, check=False):
 
             # check if item 12 exists and write item 12 - 14
             if self.segment_list is not None:
-                foo.write("# segment list for irriagation diversions\n")
+                foo.write("# segment list for irrigation diversions\n")
                 foo.write("SEGMENT LIST\n")
                 for iseg in self.segment_list:
                     foo.write(f"{iseg}\n")
diff --git a/flopy/utils/geometry.py b/flopy/utils/geometry.py
index e095a51897..490fbd82d6 100644
--- a/flopy/utils/geometry.py
+++ b/flopy/utils/geometry.py
@@ -199,6 +199,11 @@ def bounds(self):
 
         return xmin, ymin, xmax, ymax
 
+    def __reversed__(self):
+        for shp in self:
+            reversed(shp)
+        return self
+
     def plot(self, ax=None, **kwargs):
         """
         Plotting method for collection
@@ -409,6 +414,19 @@ def pyshp_parts(self):
     def patch(self):
         return self.get_patch()
 
+    def __reversed__(self):
+        # method to reverse the sorting on polygon points, patch for
+        # differences in pyshp 2.2.0 vs pyshp < 2.2.0
+        if self.exterior:
+            self.exterior = self.exterior[::-1]
+        if self.interiors:
+            interiors = []
+            for i in self.interiors:
+                interiors.append(i[::-1])
+            self.interiors = interiors
+
+        return self
+
     def get_patch(self, **kwargs):
         descartes = import_optional_dependency("descartes")
         from descartes import PolygonPatch
@@ -515,6 +533,10 @@ def bounds(self):
     def pyshp_parts(self):
         return [self.coords]
 
+    def __reversed__(self):
+        self.coords = self.coords[::-1]
+        return self
+
     def plot(self, ax=None, **kwargs):
         import matplotlib.pyplot as plt
 
@@ -611,6 +633,9 @@ def bounds(self):
     def pyshp_parts(self):
         return self.coords
 
+    def __reversed__(self):
+        return self
+
     def plot(self, ax=None, **kwargs):
         import matplotlib.pyplot as plt
 
diff --git a/flopy/utils/rasters.py b/flopy/utils/rasters.py
index f67b072f00..df56397d45 100644
--- a/flopy/utils/rasters.py
+++ b/flopy/utils/rasters.py
@@ -265,7 +265,7 @@ def sample_point(self, *point, band=1):
 
         return value
 
-    def sample_polygon(self, polygon, band, invert=False):
+    def sample_polygon(self, polygon, band, invert=False, **kwargs):
         """
         Method to get an unordered list of raster values that are located
         within a arbitrary polygon
@@ -309,7 +309,7 @@ def sample_polygon(self, polygon, band, invert=False):
                     arr_dict[b] = t
 
         else:
-            mask = self._intersection(polygon, invert)
+            mask = self._intersection(polygon, invert, **kwargs)
 
             arr_dict = {}
             for b, arr in self.__arr_dict.items():
@@ -455,7 +455,9 @@ def resample_to_grid(
             else:
                 for node in range(ncpl):
                     verts = modelgrid.get_cell_vertices(node)
-                    rstr_data = self.sample_polygon(verts, band).astype(float)
+                    rstr_data = self.sample_polygon(
+                        verts, band, convert=False
+                    ).astype(float)
                     msk = np.in1d(rstr_data, self.nodatavals)
                     rstr_data[msk] = np.nan
 
@@ -540,7 +542,9 @@ def __threaded_resampling(
         """
         container.acquire()
         verts = modelgrid.get_cell_vertices(node)
-        rstr_data = self.sample_polygon(verts, band).astype(float)
+        rstr_data = self.sample_polygon(verts, band, convert=False).astype(
+            float
+        )
         msk = np.in1d(rstr_data, self.nodatavals)
         rstr_data[msk] = np.nan
 
@@ -715,7 +719,7 @@ def _sample_rio_dataset(self, polygon, invert):
 
         return arr_dict, rstr_crp_meta
 
-    def _intersection(self, polygon, invert):
+    def _intersection(self, polygon, invert, **kwargs):
         """
         Internal method to create an intersection mask, used for cropping
         arrays and sampling arrays.
@@ -740,14 +744,18 @@ def _intersection(self, polygon, invert):
             mask : np.ndarray (dtype = bool)
 
         """
-        from .geospatial_utils import GeoSpatialUtil
+        # the convert kwarg is to speed up the resample_to_grid method
+        #  which already provides the proper datatype to _intersect()
+        convert = kwargs.pop("convert", True)
+        if convert:
+            from .geospatial_utils import GeoSpatialUtil
 
-        if isinstance(polygon, (list, tuple, np.ndarray)):
-            polygon = [polygon]
+            if isinstance(polygon, (list, tuple, np.ndarray)):
+                polygon = [polygon]
 
-        geom = GeoSpatialUtil(polygon, shapetype="Polygon")
+            geom = GeoSpatialUtil(polygon, shapetype="Polygon")
 
-        polygon = geom.points[0]
+            polygon = geom.points[0]
 
         # step 2: create a grid of centoids
         xc = self.xcenters