From d6da64264a8c24d2d79c33c30818a0293668c18e Mon Sep 17 00:00:00 2001
From: Joshua Larsen <jlarsen@usgs.gov>
Date: Tue, 9 Nov 2021 21:42:40 -0800
Subject: [PATCH 1/2] update(Raster): added min and max resampling for
 hydrologic resampling and streamflow representation

---
 flopy/utils/rasters.py | 16 ++++++++++++----
 1 file changed, 12 insertions(+), 4 deletions(-)

diff --git a/flopy/utils/rasters.py b/flopy/utils/rasters.py
index e5b39662e4..49723f0dc1 100644
--- a/flopy/utils/rasters.py
+++ b/flopy/utils/rasters.py
@@ -399,8 +399,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)):
@@ -455,8 +455,12 @@ def resample_to_grid(
 
                     if method == "median":
                         val = np.nanmedian(rstr_data)
-                    else:
+                    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:
@@ -533,8 +537,12 @@ def __threaded_resampling(
 
         if method == "median":
             val = np.nanmedian(rstr_data)
-        else:
+        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()
 

From a3083efe4ca8233c8887975fb7b3b7a47889d9c9 Mon Sep 17 00:00:00 2001
From: Joshua Larsen <jlarsen@usgs.gov>
Date: Fri, 12 Nov 2021 11:11:53 -0800
Subject: [PATCH 2/2] update(Raster): added trap for nodata within polygons

* added test_raster_sampling_methods to t065_test_gridintersect.py
---
 autotest/t065_test_gridintersect.py | 58 +++++++++++++++++++++++++++++
 flopy/utils/rasters.py              | 41 +++++++++++++-------
 2 files changed, 85 insertions(+), 14 deletions(-)

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 49723f0dc1..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
@@ -453,14 +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)
-                    elif method == "mean":
-                        val = np.nanmean(rstr_data)
-                    elif method == "max":
-                        val = np.nanmax(rstr_data)
+                    if rstr_data.size == 0:
+                        val = self.nodatavals[0]
                     else:
-                        val = np.nanmin(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:
@@ -535,14 +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)
-        elif method == "mean":
-            val = np.nanmean(rstr_data)
-        elif method == "max":
-            val = np.nanmax(rstr_data)
+        if rstr_data.size == 0:
+            val = self.nodatavals[0]
+
         else:
-            val = np.nanmin(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()