Make interp_points robust to non-equal grids and impractical defaul…

…t SciPy arguments (#484)

Signed-off-by: dependabot[bot] <support@github.com>
Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com>

doc/source/raster_class.md

@@ -317,7 +317,7 @@ rast_reproj.value_at_coords(x=0.5, y=0.5, window=3, reducer_function=np.ma.media
 
 ```{code-cell} ipython3
 # Interpolate coordinate value with quintic algorithm
-rast_reproj.interp_points([(0.5, 0.5)], mode="quintic")
+rast_reproj.interp_points([(0.5, 0.5)], method="quintic")
 ```
 
 ```{note}

examples/analysis/point_extraction/interpolation.py

@@ -50,7 +50,7 @@
 # %%
 # We can interpolate again by shifting according to our interpretation, and changing the resampling algorithm (default to "linear").
 
-vals_shifted = rast.interp_points(points=list(zip(x_coords, y_coords)), shift_area_or_point=True, mode="quintic")
+vals_shifted = rast.interp_points(points=list(zip(x_coords, y_coords)), shift_area_or_point=True, method="quintic")
 np.nanmean(vals - vals_shifted)
 
 # %%

geoutils/raster/raster.py

@@ -27,6 +27,7 @@
 from rasterio.enums import Resampling
 from rasterio.features import shapes
 from rasterio.plot import show as rshow
+from scipy.interpolate import interpn
 from scipy.ndimage import distance_transform_edt, map_coordinates
 
 import geoutils.vector as gv
@@ -3230,43 +3231,40 @@ def outside_image(self, xi: ArrayLike, yj: ArrayLike, index: bool = True) -> boo
     def interp_points(
         self,
         points: tuple[list[float], list[float]],
-        input_latlon: bool = False,
-        mode: str = "linear",
+        method: Literal["nearest", "linear", "cubic", "quintic"] = "linear",
         band: int = 1,
+        input_latlon: bool = False,
         shift_area_or_point: bool = False,
+        force_scipy_function: Literal["map_coordinates", "interpn"] | None = None,
         **kwargs: Any,
     ) -> NDArrayNum:
         """
          Interpolate raster values at a set of points.
 
+         Uses scipy.ndimage.map_coordinates if the Raster is on an equal grid, otherwise uses scipy.interpn
+         on a regular grid.
+
          Optionally, user can enforce the interpretation of pixel coordinates in self.tags['AREA_OR_POINT']
          to ensure that the interpolation of points is done at the right location. See parameter description
          of shift_area_or_point for more details.
 
         :param points: Point(s) at which to interpolate raster value. If points fall outside of image, value
             returned is nan. Shape should be (N,2).
+        :param method: Interpolation method, one of 'nearest', 'linear', 'cubic', or 'quintic'. For more information,
+            see scipy.ndimage.map_coordinates and scipy.interpolate.interpn. Default is linear.
+        :param band: Band to use (from 1 to self.count).
         :param input_latlon: Whether the input is in latlon, unregarding of Raster CRS
-        :param mode: One of 'linear', 'cubic', or 'quintic'. Determines what type of spline is used to
-            interpolate the raster value at each point. For more information, see scipy.interpolate.interp2d.
-            Default is linear.
-        :param band: The band to use (from 1 to self.count).
         :param shift_area_or_point: Shifts index to center pixel coordinates if GDAL's AREA_OR_POINT
             attribute (in self.tags) is "Point", keeps the corner pixel coordinate for "Area".
+        :param force_scipy_function: Force to use either map_coordinates or interpn. Mainly for testing purposes.
 
         :returns rpoints: Array of raster value(s) for the given points.
         """
-        assert mode in [
-            "mean",
-            "linear",
-            "cubic",
-            "quintic",
-            "nearest",
-        ], "mode must be mean, linear, cubic, quintic or nearest."
 
         # Get coordinates
         x, y = list(zip(*points))
 
-        # If those are in latlon, convert to Raster crs
+        # If those are in latlon, convert to Raster CRS
         if input_latlon:
             init_crs = pyproj.CRS(4326)
             dest_crs = pyproj.CRS(self.crs)
@@ -3277,10 +3275,42 @@ def interp_points(
 
         ind_invalid = np.vectorize(lambda k1, k2: self.outside_image(k1, k2, index=True))(j, i)
 
-        if self.count == 1:
-            rpoints = map_coordinates(self.data.astype(np.float32), [i, j], **kwargs)
+        array = self.get_nanarray()
+        if self.count != 1:
+            array = array[band - 1, :, :]
+
+        # If the raster is on an equal grid, use scipy.ndimage.map_coordinates
+        force_map_coords = force_scipy_function is not None and force_scipy_function == "map_coordinates"
+        force_interpn = force_scipy_function is not None and force_scipy_function == "interpn"
+
+        if (self.res[0] == self.res[1] or force_map_coords) and not force_interpn:
+
+            # Convert method name into order
+            method_to_order = {"nearest": 0, "linear": 1, "cubic": 3, "quintic": 5}
+            order = method_to_order[method]
+
+            # Remove default spline pre-filtering that is activated by default
+            if "prefilter" not in kwargs.keys():
+                kwargs.update({"prefilter": False})
+            # Change default constant value to NaN for interpolation outside the image bounds
+            if "cval" not in kwargs.keys():
+                kwargs.update({"cval": np.nan})
+
+            rpoints = map_coordinates(array, [i, j], order=order, **kwargs)
+
+        # Otherwise, use scipy.interpolate.interpn
         else:
-            rpoints = map_coordinates(self.data[band - 1, :, :].astype(np.float32), [i, j], **kwargs)
+
+            xycoords = self.coords(offset="corner", grid=False)
+
+            # Let interpolation outside the bounds not raise any error by default
+            if "bounds_error" not in kwargs.keys():
+                kwargs.update({"bounds_error": False})
+            # Return NaN outside image bounds
+            if "fill_value" not in kwargs.keys():
+                kwargs.update({"fill_value": np.nan})
+
+            rpoints = interpn(xycoords, self.get_nanarray(), np.array([i, j]).T, method=method, **kwargs)
 
         rpoints = np.array(rpoints, dtype=np.float32)
         rpoints[np.array(ind_invalid)] = np.nan

tests/test_raster.py

@@ -1766,8 +1766,8 @@ def test_ij2xy_xy2ij(self, example: str) -> None:
         # r.ds.index(x, y, op=np.float32)
         # Out[34]: (75.0, 302.0)
 
-    def test_xy2ij_and_interp(self) -> None:
-        """Test xy2ij with shift_area_or_point argument, and related interp function"""
+    def test_xy2ij(self) -> None:
+        """Test xy2ij with shift_area_or_point argument, and compare to interp_points function for consistency."""
 
         # First, we try on a Raster with a Point interpretation in its "AREA_OR_POINT" metadata: values interpolated
         # at the center of pixel
@@ -1834,7 +1834,7 @@ def test_xy2ij_and_interp(self) -> None:
             list_z_ind.append(z_ind)
 
         # First order interpolation
-        rpts = r.interp_points(pts, order=1)
+        rpts = r.interp_points(pts, method="linear")
         # The values interpolated should be equal
         assert np.array_equal(np.array(list_z_ind, dtype=np.float32), rpts, equal_nan=True)
 
@@ -1866,23 +1866,179 @@ def test_xy2ij_and_interp(self) -> None:
             z_ind = img[int(i[k]), int(j[k])]
             list_z_ind.append(z_ind)
 
-        rpts = r.interp_points(pts, order=1)
+        rpts = r.interp_points(pts, method="linear")
 
         assert np.array_equal(np.array(list_z_ind, dtype=np.float32), rpts, equal_nan=True)
 
         # Test for an invidiual point (shape can be tricky at 1 dimension)
         x = 493120.0
         y = 3101000.0
         i, j = r.xy2ij(x, y)
-        val = r.interp_points([(x, y)], order=1)[0]
+        val = r.interp_points([(x, y)], method="linear")[0]
         val_img = img[int(i[0]), int(j[0])]
         assert val_img == val
 
         # Finally, check that interp convert to latlon
         lat, lon = gu.projtools.reproject_to_latlon([x, y], in_crs=r.crs)
-        val_latlon = r.interp_points([(lat, lon)], order=1, input_latlon=True)[0]
+        val_latlon = r.interp_points([(lat, lon)], method="linear", input_latlon=True)[0]
         assert val == pytest.approx(val_latlon, abs=0.0001)
 
+    @pytest.mark.parametrize("tag_aop", [None, "Area", "Point"])  # type: ignore
+    @pytest.mark.parametrize("shift_aop", [True, False])  # type: ignore
+    def test_interp_points__synthetic(self, tag_aop: str | None, shift_aop: bool) -> None:
+        """Test interp_points function with synthetic data."""
+
+        # We flip the array up/down to facilitate index comparison of Y axis
+        arr = np.flipud(np.array([1, 2, 3, 4, 5, 6, 7, 8, 9]).reshape((3, 3)))
+        transform = rio.transform.from_bounds(0, 0, 3, 3, 3, 3)
+        raster = gu.Raster.from_array(data=arr, transform=transform, crs=None, nodata=-9999)
+
+        # Define the AREA_OR_POINT attribute
+        raster.tags = {"AREA_OR_POINT": tag_aop}
+
+        # Check interpolation falls right on values for points (1, 1), (1, 2) etc...
+        index_x = [0, 1, 2, 0, 1, 2, 0, 1, 2]
+        index_y = [0, 0, 0, 1, 1, 1, 2, 2, 2]
+
+        # The actual X/Y coords will be offset by one because Y axis is inverted and pixel coords is upper-left corner
+
+        # For the single case of "Point" and shift=True, shift the point coordinates to match the calculations below
+        # All other cases should pass without this shift (for tag=None, "Area" is the default)
+        if shift_aop and tag_aop == "Point":
+            reindex_x = [ix - 0.5 for ix in index_x]
+            reindex_y = [iy - 0.5 for iy in index_y]
+            points_x, points_y = raster.ij2xy(i=reindex_x, j=reindex_y)
+        else:
+            points_x, points_y = raster.ij2xy(i=index_x, j=index_y)
+
+        points = np.array((points_x, points_y)).T
+
+        # The following 4 methods should yield the same result because:
+        # Nearest = Linear interpolation at the location of a data point
+        # Regular grid = Equal grid interpolation at the location of a data point
+
+        # Check warning is raised if AREA_OR_POINT is None
+        if tag_aop is None and shift_aop:
+            with pytest.warns(UserWarning, match="Attribute AREA_OR_POINT undefined in self.tags*"):
+                raster_points = raster.interp_points(points, method="nearest", shift_area_or_point=shift_aop)
+        else:
+            raster_points = raster.interp_points(points, method="nearest", shift_area_or_point=shift_aop)
+
+        # Warnings should be raised when AREA_OR_POINT is None everywhere, so we ignore them from now on
+        if tag_aop is None and shift_aop:
+            warnings.filterwarnings(
+                "ignore", category=UserWarning, message="Attribute AREA_OR_POINT undefined in self.tags*"
+            )
+
+        raster_points_lin = raster.interp_points(points, method="linear", shift_area_or_point=shift_aop)
+        raster_points_interpn = raster.interp_points(
+            points, method="nearest", force_scipy_function="interpn", shift_area_or_point=shift_aop
+        )
+        raster_points_interpn_lin = raster.interp_points(
+            points, method="linear", force_scipy_function="interpn", shift_area_or_point=shift_aop
+        )
+
+        assert np.array_equal(raster_points, raster_points_lin)
+        assert np.array_equal(raster_points, raster_points_interpn)
+        assert np.array_equal(raster_points, raster_points_interpn_lin)
+
+        for i in range(3):
+            for j in range(3):
+                ind = 3 * i + j
+                assert raster_points[ind] == arr[index_x[ind], index_y[ind]]
+
+        # Check bilinear interpolation values inside the grid (same here, offset by 1 between X and Y)
+        index_x_in = [0.5, 0.5, 1.5, 1.5]
+        index_y_in = [0.5, 1.5, 0.5, 1.5]
+
+        # Again, exception for "Point"
+        if shift_aop and tag_aop == "Point":
+            reindex_x_in = [ix - 0.5 for ix in index_x_in]
+            reindex_y_in = [iy - 0.5 for iy in index_y_in]
+            points_x_in, points_y_in = raster.ij2xy(i=reindex_x_in, j=reindex_y_in)
+        else:
+            points_x_in, points_y_in = raster.ij2xy(i=index_x_in, j=index_y_in)
+
+        points_in = np.array((points_x_in, points_y_in)).T
+
+        # Here again compare methods
+        raster_points_in = raster.interp_points(points_in, method="linear", shift_area_or_point=shift_aop)
+        raster_points_in_interpn = raster.interp_points(
+            points_in, method="linear", force_scipy_function="interpn", shift_area_or_point=shift_aop
+        )
+
+        assert np.array_equal(raster_points_in, raster_points_in_interpn)
+
+        for i in range(len(points_in)):
+
+            xlow = int(index_x_in[i] - 0.5)
+            xupp = int(index_x_in[i] + 0.5)
+            ylow = int(index_y_in[i] - 0.5)
+            yupp = int(index_y_in[i] + 0.5)
+
+            # Check the bilinear interpolation matches the mean value of those 4 points (equivalent as its the middle)
+            assert raster_points_in[i] == np.mean([arr[xlow, ylow], arr[xupp, ylow], arr[xupp, yupp], arr[xlow, yupp]])
+
+        # Check bilinear extrapolation for points at 1 spacing outside from the input grid
+        points_out = (
+            [(-1, i) for i in np.arange(1, 4)]
+            + [(i, -1) for i in np.arange(1, 4)]
+            + [(4, i) for i in np.arange(1, 4)]
+            + [(i, 4) for i in np.arange(4, 1)]
+        )
+        raster_points_out = raster.interp_points(points_out)
+        assert all(~np.isfinite(raster_points_out))
+
+        # To use cubic or quintic, we need a larger grid (minimum 6x6, but let's aim bigger with 50x50)
+        arr = np.flipud(np.arange(1, 2501).reshape((50, 50)))
+        transform = rio.transform.from_bounds(0, 0, 50, 50, 50, 50)
+        raster = gu.Raster.from_array(data=arr, transform=transform, crs=None, nodata=-9999)
+        raster.tags = {"AREA_OR_POINT": tag_aop}
+
+        # For this, get random points
+        np.random.seed(42)
+        index_x_in_rand = np.random.randint(low=8, high=42, size=(10,)) + np.random.normal(scale=0.3)
+        index_y_in_rand = np.random.randint(low=8, high=42, size=(10,)) + np.random.normal(scale=0.3)
+        points_x_rand, points_y_rand = raster.ij2xy(i=index_x_in_rand, j=index_y_in_rand)
+        points_in_rand = np.array((points_x_rand, points_y_rand)).T
+
+        for method in ["nearest", "linear", "cubic", "quintic"]:
+            raster_points_mapcoords = raster.interp_points(
+                points_in_rand, method=method, force_scipy_function="map_coordinates", shift_area_or_point=shift_aop
+            )
+            raster_points_interpn = raster.interp_points(
+                points_in_rand, method=method, force_scipy_function="interpn", shift_area_or_point=shift_aop
+            )
+
+            assert np.array_equal(raster_points_mapcoords, raster_points_interpn)
+
+        # Check that, outside the edge, the interpolation fails and returns a NaN
+        np.random.seed(42)
+        index_x_edge_rand = [-0.5, -0.5, -0.5, 25, 25, 49.5, 49.5, 49.5]
+        index_y_edge_rand = [-0.5, 25, 49.5, -0.5, 49.5, -0.5, 25, 49.5]
+
+        # Again, exception for "Point"
+        if shift_aop and tag_aop == "Point":
+            reindex_x_edge_rand = [ix - 0.5 for ix in index_x_edge_rand]
+            reindex_y_edge_rand = [iy - 0.5 for iy in index_y_edge_rand]
+            points_x_rand, points_y_rand = raster.ij2xy(i=reindex_x_edge_rand, j=reindex_y_edge_rand)
+        else:
+            points_x_rand, points_y_rand = raster.ij2xy(i=index_x_edge_rand, j=index_y_edge_rand)
+
+        points_edge_rand = np.array((points_x_rand, points_y_rand)).T
+
+        # Nearest doesn't apply, just linear and above
+        for method in ["cubic", "quintic"]:
+            raster_points_mapcoords_edge = raster.interp_points(
+                points_edge_rand, method=method, force_scipy_function="map_coordinates", shift_area_or_point=shift_aop
+            )
+            raster_points_interpn_edge = raster.interp_points(
+                points_edge_rand, method=method, force_scipy_function="interpn", shift_area_or_point=shift_aop
+            )
+
+            assert all(~np.isfinite(raster_points_mapcoords_edge))
+            assert all(~np.isfinite(raster_points_interpn_edge))
+
     def test_value_at_coords(self) -> None:
         """
         Test that value at coords works as intended