diff --git a/autotest/t065_test_gridintersect.py b/autotest/t065_test_gridintersect.py
index 63d2f65966..1ec78ed685 100644
--- a/autotest/t065_test_gridintersect.py
+++ b/autotest/t065_test_gridintersect.py
@@ -1361,5 +1361,63 @@ def test_rasters():
     del rio
 
 
+# %% test raster sampling methods
+
+
+def test_raster_sampling_methods():
+    from flopy.utils import Raster
+    import os
+    import flopy as fp
+
+    ws = os.path.join("..", "examples", "data", "options")
+    raster_name = "dem.img"
+
+    try:
+        rio = Raster.load(os.path.join(ws, "dem", raster_name))
+    except:
+        return
+
+    ml = fp.modflow.Modflow.load(
+        "sagehen.nam", version="mfnwt", model_ws=os.path.join(ws, "sagehen")
+    )
+    xoff = 214110
+    yoff = 4366620
+    ml.modelgrid.set_coord_info(xoff, yoff)
+
+    x0, x1, y0, y1 = rio.bounds
+
+    x0 += 3000
+    y0 += 3000
+    x1 -= 3000
+    y1 -= 3000
+    shape = np.array([(x0, y0), (x0, y1), (x1, y1), (x1, y0), (x0, y0)])
+
+    rio.crop(shape)
+
+    methods = {
+        "min": 2088.52343,
+        "max": 2103.54882,
+        "mean": 2097.05053,
+        "median": 2097.36254,
+        "nearest": 2097.81079,
+        "linear": 2097.81079,
+        "cubic": 2097.81079
+    }
+
+    for method, value in methods.items():
+        data = rio.resample_to_grid(
+            ml.modelgrid,
+            band=rio.bands[0],
+            method=method
+        )
+
+        print(data[30, 37])
+        if np.abs(data[30, 37] - value) > 1e-05:
+            raise AssertionError(
+                f"{method} resampling returning incorrect values"
+            )
+
+
 if __name__ == "__main__":
     test_rasters()
+    test_raster_sampling_methods()
\ No newline at end of file
diff --git a/flopy/utils/rasters.py b/flopy/utils/rasters.py
index e5b39662e4..3aa6a64559 100644
--- a/flopy/utils/rasters.py
+++ b/flopy/utils/rasters.py
@@ -351,6 +351,11 @@ def resample_to_grid(
             ``mean`` for mean sampling
 
             ``median`` for median sampling
+
+            ``min`` for minimum sampling
+
+            ``max`` for maximum sampling
+
         multithread : bool
             boolean flag indicating if multithreading should be used with
             the ``mean`` and ``median`` sampling methods
@@ -399,8 +404,8 @@ def resample_to_grid(
                 method=method,
             )
 
-        elif method in ("median", "mean"):
-            # these methods are slow and could use a speed u
+        elif method in ("median", "mean", "min", "max"):
+            # these methods are slow and could use a speed up
             ncpl = modelgrid.ncpl
             data_shape = modelgrid.xcellcenters.shape
             if isinstance(ncpl, (list, np.ndarray)):
@@ -453,10 +458,17 @@ def resample_to_grid(
                     msk = np.in1d(rstr_data, self.nodatavals)
                     rstr_data[msk] = np.nan
 
-                    if method == "median":
-                        val = np.nanmedian(rstr_data)
+                    if rstr_data.size == 0:
+                        val = self.nodatavals[0]
                     else:
-                        val = np.nanmean(rstr_data)
+                        if method == "median":
+                            val = np.nanmedian(rstr_data)
+                        elif method == "mean":
+                            val = np.nanmean(rstr_data)
+                        elif method == "max":
+                            val = np.nanmax(rstr_data)
+                        else:
+                            val = np.nanmin(rstr_data)
 
                     data[node] = val
         else:
@@ -531,10 +543,19 @@ def __threaded_resampling(
         msk = np.in1d(rstr_data, self.nodatavals)
         rstr_data[msk] = np.nan
 
-        if method == "median":
-            val = np.nanmedian(rstr_data)
+        if rstr_data.size == 0:
+            val = self.nodatavals[0]
+
         else:
-            val = np.nanmean(rstr_data)
+            if method == "median":
+                val = np.nanmedian(rstr_data)
+            elif method == "mean":
+                val = np.nanmean(rstr_data)
+            elif method == "max":
+                val = np.nanmax(rstr_data)
+            else:
+                val = np.nanmin(rstr_data)
+
         q.put((node, val))
         container.release()