From f1885830ecbeda82b3481c029c7ba9aac1ae65cc Mon Sep 17 00:00:00 2001
From: Romain Hugonnet <romain.hugonnet@gmail.com>
Date: Mon, 26 Aug 2024 20:13:49 -0700
Subject: [PATCH] Fix resampling not being passed to `stack_rasters` and add
 tests (#601)

---
 geoutils/raster/multiraster.py        |  1 +
 tests/test_raster/test_multiraster.py | 12 +++++++++++-
 2 files changed, 12 insertions(+), 1 deletion(-)

diff --git a/geoutils/raster/multiraster.py b/geoutils/raster/multiraster.py
index 5d7f6b2c..6d1c424b 100644
--- a/geoutils/raster/multiraster.py
+++ b/geoutils/raster/multiraster.py
@@ -181,6 +181,7 @@ def stack_rasters(
             crs=reference_raster.crs,
             dtype=reference_raster.data.dtype,
             nodata=reference_raster.nodata,
+            resampling=resampling_method,
             silent=True,
         )
         reprojected_raster.set_nodata(nodata)
diff --git a/tests/test_raster/test_multiraster.py b/tests/test_raster/test_multiraster.py
index 78fc47b7..2f9e9949 100644
--- a/tests/test_raster/test_multiraster.py
+++ b/tests/test_raster/test_multiraster.py
@@ -52,7 +52,7 @@ def __init__(
             inplace=True,
         )
         if different_crs:
-            self.img2 = self.img2.reproject(crs=different_crs)
+            self.img2 = self.img2.reproject(crs=different_crs, resampling="nearest")
 
         # To check that use_ref_bounds work - create a img that do not cover the whole extent
         self.img3 = img.copy()
@@ -160,6 +160,16 @@ def test_stack_rasters(self, rasters) -> None:  # type: ignore
         stacked_img2 = gu.raster.stack_rasters([rasters.img1, rasters.img3], reference=rasters.img, use_ref_bounds=True)
         assert stacked_img2.bounds == rasters.img.bounds
 
+        # This case should preserve unique data values through "nearest" resampling
+        rasters.img1[:] = 5
+        rasters.img1[0:5, 0:5] = 1
+        rasters.img2 = rasters.img1.translate(0.5, 0.5, distance_unit="pixel")
+        stacked_img = gu.raster.stack_rasters([rasters.img1, rasters.img2], resampling_method="nearest")
+        assert np.array_equal(np.unique(stacked_img.data.compressed()), np.array([1, 5]))
+        # But not this case with a shifted raster resampled with "bilinear"
+        stacked_img = gu.raster.stack_rasters([rasters.img1, rasters.img2], resampling_method="bilinear")
+        assert not np.array_equal(np.unique(stacked_img.data.compressed()), np.array([1, 5]))
+
     @pytest.mark.parametrize(
         "rasters",
         [