diff --git a/geoutils/raster/georeferencing.py b/geoutils/raster/georeferencing.py
new file mode 100644
index 00000000..da02d927
--- /dev/null
+++ b/geoutils/raster/georeferencing.py
@@ -0,0 +1,172 @@
+from __future__ import annotations
+
+from typing import Iterable, Literal
+
+import numpy as np
+import rasterio as rio
+
+from geoutils._config import config
+from geoutils._typing import ArrayLike, NDArrayNum
+
+
+def _ij2xy(
+    i: ArrayLike,
+    j: ArrayLike,
+    transform: rio.transform.Affine,
+    area_or_point: Literal["Area", "Point"] | None,
+    shift_area_or_point: bool | None = None,
+    force_offset: str | None = None,
+) -> tuple[NDArrayNum, NDArrayNum]:
+    """See description of Raster.ij2xy."""
+
+    # If undefined, default to the global system config
+    if shift_area_or_point is None:
+        shift_area_or_point = config["shift_area_or_point"]
+
+    # Shift by half a pixel back for "Point" interpretation
+    if shift_area_or_point and force_offset is None:
+        if area_or_point is not None and area_or_point == "Point":
+            i = np.asarray(i) - 0.5
+            j = np.asarray(j) - 0.5
+
+    # Default offset is upper-left for raster coordinates
+    if force_offset is None:
+        force_offset = "ul"
+
+    x, y = rio.transform.xy(transform, i, j, offset=force_offset)
+
+    return x, y
+
+
+def _xy2ij(
+    x: ArrayLike,
+    y: ArrayLike,
+    transform: rio.transform.Affine,
+    area_or_point: Literal["Area", "Point"] | None,
+    op: type = np.float32,
+    precision: float | None = None,
+    shift_area_or_point: bool | None = None,
+) -> tuple[NDArrayNum, NDArrayNum]:
+    """See description of Raster.xy2ij."""
+
+    # If undefined, default to the global system config
+    if shift_area_or_point is None:
+        shift_area_or_point = config["shift_area_or_point"]
+
+    # Input checks
+    if op not in [np.float32, np.float64, float]:
+        raise UserWarning(
+            "Operator is not of type float: rio.Dataset.index might "
+            "return unreliable indexes due to rounding issues."
+        )
+
+    i, j = rio.transform.rowcol(transform, x, y, op=op, precision=precision)
+
+    # Necessary because rio.Dataset.index does not return abc.Iterable for a single point
+    if not isinstance(i, Iterable):
+        i, j = (
+            np.asarray(
+                [
+                    i,
+                ]
+            ),
+            np.asarray(
+                [
+                    j,
+                ]
+            ),
+        )
+    else:
+        i, j = (np.asarray(i), np.asarray(j))
+
+    # AREA_OR_POINT GDAL attribute, i.e. does the value refer to the upper left corner "Area" or
+    # the center of pixel "Point". This normally has no influence on georeferencing, it's only
+    # about the interpretation of the raster values, and thus can affect sub-pixel interpolation,
+    # for more details see: https://gdal.org/user/raster_data_model.html#metadata
+
+    # If the user wants to shift according to the interpretation
+    if shift_area_or_point:
+
+        # Shift by half a pixel if the AREA_OR_POINT attribute is "Point", otherwise leave as is
+        if area_or_point is not None and area_or_point == "Point":
+            if not isinstance(i.flat[0], (np.floating, float)):
+                raise ValueError("Operator must return np.floating values to perform pixel interpretation shifting.")
+
+            i += 0.5
+            j += 0.5
+
+    # Convert output indexes to integer if they are all whole numbers
+    if np.all(np.mod(i, 1) == 0) and np.all(np.mod(j, 1) == 0):
+        i = i.astype(int)
+        j = j.astype(int)
+
+    return i, j
+
+
+def _coords(
+    transform: rio.transform.Affine,
+    shape: tuple[int, int],
+    area_or_point: Literal["Area", "Point"] | None,
+    grid: bool = True,
+    shift_area_or_point: bool | None = None,
+    force_offset: str | None = None,
+) -> tuple[NDArrayNum, NDArrayNum]:
+    """See description of Raster.coords."""
+
+    # The coordinates are extracted from indexes 0 to shape
+    _, yy = _ij2xy(
+        i=np.arange(shape[0] - 1, -1, -1),
+        j=0,
+        transform=transform,
+        area_or_point=area_or_point,
+        shift_area_or_point=shift_area_or_point,
+        force_offset=force_offset,
+    )
+    xx, _ = _ij2xy(
+        i=0,
+        j=np.arange(shape[1]),
+        transform=transform,
+        area_or_point=area_or_point,
+        shift_area_or_point=shift_area_or_point,
+        force_offset=force_offset,
+    )
+
+    # If grid is True, return coordinate grids
+    if grid:
+        meshgrid = tuple(np.meshgrid(xx, np.flip(yy)))
+        return meshgrid  # type: ignore
+    else:
+        return np.asarray(xx), np.asarray(yy)
+
+
+def _outside_image(
+    xi: ArrayLike,
+    yj: ArrayLike,
+    transform: rio.transform.Affine,
+    shape: tuple[int, int],
+    area_or_point: Literal["Area", "Point"] | None,
+    index: bool = True,
+) -> bool:
+    """See description of Raster.outside_image."""
+
+    if not index:
+        xi, xj = _xy2ij(xi, yj, transform=transform, area_or_point=area_or_point)
+
+    if np.any(np.array((xi, yj)) < 0):
+        return True
+    elif np.asanyarray(xi) > shape[1] or np.asanyarray(yj) > shape[0]:
+        return True
+    else:
+        return False
+
+
+def _res(transform: rio.transform.Affine) -> tuple[float, float]:
+    """See description of Raster.res"""
+
+    return transform[0], abs(transform[4])
+
+
+def _bounds(transform: rio.transform.Affine, shape: tuple[int, int]) -> rio.coords.BoundingBox:
+    """See description of Raster.bounds."""
+
+    return rio.coords.BoundingBox(*rio.transform.array_bounds(height=shape[0], width=shape[1], transform=transform))
diff --git a/geoutils/raster/interpolate.py b/geoutils/raster/interpolate.py
new file mode 100644
index 00000000..153d9955
--- /dev/null
+++ b/geoutils/raster/interpolate.py
@@ -0,0 +1,329 @@
+from __future__ import annotations
+
+from typing import Any, Callable, Literal, overload
+
+import numpy as np
+import rasterio as rio
+from scipy.interpolate import RectBivariateSpline, RegularGridInterpolator
+from scipy.ndimage import binary_dilation, distance_transform_edt, map_coordinates
+
+from geoutils._typing import NDArrayNum, Number
+from geoutils.raster.georeferencing import _coords, _outside_image, _res, _xy2ij
+
+method_to_order = {"nearest": 0, "linear": 1, "cubic": 3, "quintic": 5, "slinear": 1, "pchip": 3, "splinef2d": 3}
+
+
+def _get_dist_nodata_spread(order: int, dist_nodata_spread: Literal["half_order_up", "half_order_down"] | int) -> int:
+    """
+    Derive distance of nodata spreading based on interpolation order.
+
+    :param order: Interpolation order.
+    :param dist_nodata_spread: Spreading distance of nodata, either half-order rounded up (default), rounded down, or
+        fixed integer.
+    """
+
+    if dist_nodata_spread == "half_order_up":
+        dist_nodata_spread = int(np.ceil(order / 2))
+    elif dist_nodata_spread == "half_order_down":
+        dist_nodata_spread = int(np.floor(order / 2))
+
+    return dist_nodata_spread
+
+
+def _interpn_interpolator(
+    points: tuple[NDArrayNum, NDArrayNum],
+    values: NDArrayNum,
+    fill_value: Number = np.nan,
+    bounds_error: bool = False,
+    dist_nodata_spread: Literal["half_order_up", "half_order_down"] | int = "half_order_up",
+    method: Literal["nearest", "linear", "cubic", "quintic", "slinear", "pchip", "splinef2d"] = "linear",
+) -> Callable[[tuple[NDArrayNum, NDArrayNum]], NDArrayNum]:
+    """
+    Create SciPy interpolator with nodata spreading. Default is spreading at distance of half the method order
+    rounded up (i.e., linear spreads 1 nodata in each direction, cubic spreads 2, quintic 3).
+
+    Gives the exact same result as scipy.interpolate.interpn, and allows interpolator to be re-used if required (
+    for speed).
+    In practice, returns either a NaN-modified RegularGridInterpolator or a NaN-modified RectBivariateSpline object,
+    both expecting a tuple of X/Y coordinates to be evaluated.
+
+    For input arguments, see scipy.interpolate.RegularGridInterpolator.
+    For additional argument "dist_nodata_spread", see description of Raster.interp_points.
+
+    Adapted from:
+    https://github.com/scipy/scipy/blob/44e4ebaac992fde33f04638b99629d23973cb9b2/scipy/interpolate/_rgi.py#L743.
+    """
+
+    # Derive distance to spread nodata to depending on method order
+    order = method_to_order[method]
+    d = _get_dist_nodata_spread(order=order, dist_nodata_spread=dist_nodata_spread)
+
+    # We compute the nodata mask and dilate it to the distance to spread nodatas
+    mask_nan = ~np.isfinite(values)
+    if d != 0:
+        new_mask = binary_dilation(mask_nan, iterations=d).astype("uint8")
+    # Zero iterations has a different behaviour in binary_dilation than doing nothing, we want the original array
+    else:
+        new_mask = mask_nan.astype("uint8")
+
+    # We create an interpolator for the nodata mask using nearest
+    interp_mask = RegularGridInterpolator(points, new_mask, method="nearest", bounds_error=bounds_error, fill_value=1)
+
+    # Most methods (cubic, quintic, etc) do not support NaNs and require an array full of valid values
+    # We replace thus replace all NaN values by nearest neighbours to give surrounding values of the same order of
+    # magnitude and minimize interpolation errors near NaNs (errors of 10e-2/e-5 relative to the values)
+    # Elegant solution from: https://stackoverflow.com/questions/5551286/filling-gaps-in-a-numpy-array for a fast
+    # nearest neighbour fill
+    indices = distance_transform_edt(mask_nan, return_distances=False, return_indices=True)
+    values = values[tuple(indices)]
+
+    # For the RegularGridInterpolator
+    if method in RegularGridInterpolator._ALL_METHODS:
+
+        # We create the classic interpolator
+        interp = RegularGridInterpolator(
+            points, values, method=method, bounds_error=bounds_error, fill_value=fill_value
+        )
+
+        # We create a new interpolator callable that propagates nodata as defined above
+        def regulargrid_interpolator_with_nan(xi: tuple[NDArrayNum, NDArrayNum]) -> NDArrayNum:
+
+            # Get results
+            results = interp(xi)
+            # Get invalids
+            invalids = interp_mask(xi)
+            results[invalids.astype(bool)] = np.nan
+
+            return results
+
+        return regulargrid_interpolator_with_nan
+
+    # For the RectBivariateSpline
+    else:
+
+        # The coordinates must be in ascending order, which requires flipping the array too (more costly)
+        interp = RectBivariateSpline(np.flip(points[0]), points[1], np.flip(values[:], axis=0))
+
+        # We create a new interpolator callable that propagates nodata as defined above, and supports fill_value
+        def rectbivariate_interpolator_with_fillvalue(xi: tuple[NDArrayNum, NDArrayNum]) -> NDArrayNum:
+
+            # Get invalids
+            invalids = interp_mask(xi)
+
+            # RectBivariateSpline doesn't support fill_value, so we need to wrap here to add them
+            xi_arr = np.array(xi).T
+            xi_shape = xi_arr.shape
+            xi_arr = xi_arr.reshape(-1, xi_arr.shape[-1])
+            idx_valid = np.all(
+                (
+                    points[0][-1] <= xi_arr[:, 0],
+                    xi_arr[:, 0] <= points[0][0],
+                    points[1][0] <= xi_arr[:, 1],
+                    xi_arr[:, 1] <= points[1][-1],
+                ),
+                axis=0,
+            )
+            # Make a copy of values for RectBivariateSpline
+            result = np.empty_like(xi_arr[:, 0])
+            result[idx_valid] = interp.ev(xi_arr[idx_valid, 0], xi_arr[idx_valid, 1])
+            result[np.logical_not(idx_valid)] = fill_value
+
+            # Add back NaNs from dilated mask
+            results = np.atleast_1d(result.reshape(xi_shape[:-1]))
+            results[invalids.astype(bool)] = np.nan
+
+            return results
+
+        return rectbivariate_interpolator_with_fillvalue
+
+
+def _map_coordinates_nodata_propag(
+    values: NDArrayNum,
+    indices: tuple[NDArrayNum, NDArrayNum],
+    order: int,
+    dist_nodata_spread: Literal["half_order_up", "half_order_down"] | int = "half_order_up",
+    **kwargs: Any,
+) -> NDArrayNum:
+    """
+    Perform map_coordinates with nodata spreading. Default is spreading at distance of half the method order rounded
+    up (i.e., linear spreads 1 nodata in each direction, cubic spreads 2, quintic 3).
+
+    For map_coordinates, only nearest and linear are used.
+
+    For input arguments, see scipy.ndimage.map_coordinates.
+    For additional argument "dist_nodata_spread", see description of Raster.interp_points.
+    """
+
+    # Derive distance of nodata spreading
+    d = _get_dist_nodata_spread(order=order, dist_nodata_spread=dist_nodata_spread)
+
+    # We compute the mask and dilate it to the distance to spread nodatas
+    mask_nan = ~np.isfinite(values)
+    if d != 0:
+        new_mask = binary_dilation(mask_nan, iterations=d).astype("uint8")
+    # Zero iterations has a different behaviour in binary_dilation than doing nothing, here we want the original array
+    else:
+        new_mask = mask_nan.astype("uint8")
+
+    # We replace all NaN values by nearest neighbours to minimize interpolation errors near NaNs
+    # Elegant solution from: https://stackoverflow.com/questions/5551286/filling-gaps-in-a-numpy-array
+    ind = distance_transform_edt(mask_nan, return_distances=False, return_indices=True)
+    values = values[tuple(ind)]
+
+    # We interpolate the dilated array at the coordinates with nearest, and transform it back to a boolean to mask NaNs
+    rmask = map_coordinates(new_mask, indices, order=0, cval=1, prefilter=False)
+
+    # Interpolate at indices
+    rpoints = map_coordinates(values, indices, order=order, **kwargs)
+
+    # Set to NaNs based on spreading distance
+    rpoints[rmask.astype(bool)] = np.nan
+
+    return rpoints
+
+
+@overload
+def _interp_points(
+    array: NDArrayNum,
+    transform: rio.transform.Affine,
+    area_or_point: Literal["Area", "Point"] | None,
+    points: tuple[Number, Number] | tuple[NDArrayNum, NDArrayNum],
+    method: Literal["nearest", "linear", "cubic", "quintic", "slinear", "pchip", "splinef2d"] = "linear",
+    dist_nodata_spread: Literal["half_order_up", "half_order_down"] | int = "half_order_up",
+    shift_area_or_point: bool | None = None,
+    force_scipy_function: Literal["map_coordinates", "interpn"] | None = None,
+    *,
+    return_interpolator: Literal[False] = False,
+    **kwargs: Any,
+) -> NDArrayNum:
+    ...
+
+
+@overload
+def _interp_points(
+    array: NDArrayNum,
+    transform: rio.transform.Affine,
+    area_or_point: Literal["Area", "Point"] | None,
+    points: tuple[Number, Number] | tuple[NDArrayNum, NDArrayNum],
+    method: Literal["nearest", "linear", "cubic", "quintic", "slinear", "pchip", "splinef2d"] = "linear",
+    dist_nodata_spread: Literal["half_order_up", "half_order_down"] | int = "half_order_up",
+    shift_area_or_point: bool | None = None,
+    force_scipy_function: Literal["map_coordinates", "interpn"] | None = None,
+    *,
+    return_interpolator: Literal[True],
+    **kwargs: Any,
+) -> Callable[[tuple[NDArrayNum, NDArrayNum]], NDArrayNum]:
+    ...
+
+
+@overload
+def _interp_points(
+    array: NDArrayNum,
+    transform: rio.transform.Affine,
+    area_or_point: Literal["Area", "Point"] | None,
+    points: tuple[Number, Number] | tuple[NDArrayNum, NDArrayNum],
+    method: Literal["nearest", "linear", "cubic", "quintic", "slinear", "pchip", "splinef2d"] = "linear",
+    dist_nodata_spread: Literal["half_order_up", "half_order_down"] | int = "half_order_up",
+    shift_area_or_point: bool | None = None,
+    force_scipy_function: Literal["map_coordinates", "interpn"] | None = None,
+    *,
+    return_interpolator: bool = False,
+    **kwargs: Any,
+) -> NDArrayNum | Callable[[tuple[NDArrayNum, NDArrayNum]], NDArrayNum]:
+    ...
+
+
+def _interp_points(
+    array: NDArrayNum,
+    transform: rio.transform.Affine,
+    area_or_point: Literal["Area", "Point"] | None,
+    points: tuple[Number, Number] | tuple[NDArrayNum, NDArrayNum] | None,
+    method: Literal["nearest", "linear", "cubic", "quintic", "slinear", "pchip", "splinef2d"] = "linear",
+    dist_nodata_spread: Literal["half_order_up", "half_order_down"] | int = "half_order_up",
+    shift_area_or_point: bool | None = None,
+    force_scipy_function: Literal["map_coordinates", "interpn"] | None = None,
+    return_interpolator: bool = False,
+    **kwargs: Any,
+) -> NDArrayNum | Callable[[tuple[NDArrayNum, NDArrayNum]], NDArrayNum]:
+    """See description of Raster.interp_points."""
+
+    # If array is not a floating dtype (to support NaNs), convert dtype
+    if not np.issubdtype(array.dtype, np.floating):
+        array = array.astype(np.float32)
+    shape: tuple[int, int] = array.shape[0:2]  # type: ignore
+
+    # TODO: Add check about None for "points" depending on "return_interpolator"
+    if not return_interpolator:
+        if points is None:
+            raise ValueError("Input 'points' cannot be None if 'return_interpolator' is False.")
+        x, y = points
+        i, j = _xy2ij(x, y, transform=transform, area_or_point=area_or_point, shift_area_or_point=shift_area_or_point)
+
+        ind_invalid = np.vectorize(
+            lambda k1, k2: _outside_image(
+                k1, k2, transform=transform, area_or_point=area_or_point, shape=shape, index=True
+            )
+        )(j, i)
+
+    # 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"
+
+    # Map method name to spline order in map_coordinates, and use only is method compatible
+    method_to_order_mapcoords = {"nearest": 0, "linear": 1}
+    mapcoords_supported = method in method_to_order_mapcoords.keys()
+
+    res = _res(transform)
+    if (res[0] == res[1] or force_map_coords) and not force_interpn and mapcoords_supported and not return_interpolator:
+
+        # Convert method name into order
+        order = method_to_order_mapcoords[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})
+
+        # Use map coordinates with nodata propagation
+        rpoints = _map_coordinates_nodata_propag(
+            values=array, indices=(i, j), order=order, dist_nodata_spread=dist_nodata_spread, **kwargs
+        )
+
+    # Otherwise, use scipy.interpolate.interpn
+    else:
+        # Get lower-left corner coordinates
+        xycoords = _coords(
+            transform=transform,
+            shape=shape,
+            area_or_point=area_or_point,
+            grid=False,
+            shift_area_or_point=shift_area_or_point,
+        )
+
+        # 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})
+
+        # Using direct coordinates, Y is the first axis, and we need to flip it
+        interpolator = _interpn_interpolator(
+            points=(np.flip(xycoords[1], axis=0), xycoords[0]),
+            values=array,
+            method=method,
+            dist_nodata_spread=dist_nodata_spread,
+            bounds_error=kwargs["bounds_error"],
+            fill_value=kwargs["fill_value"],
+        )
+        if return_interpolator:
+            return interpolator
+        else:
+            rpoints = interpolator((y, x))  # type: ignore
+
+    rpoints = np.array(np.atleast_1d(rpoints), dtype=np.float32)
+    rpoints[np.array(ind_invalid)] = np.nan
+
+    return rpoints
diff --git a/geoutils/raster/raster.py b/geoutils/raster/raster.py
index 696e8187..2d615dab 100644
--- a/geoutils/raster/raster.py
+++ b/geoutils/raster/raster.py
@@ -17,7 +17,6 @@
 import matplotlib
 import matplotlib.pyplot as plt
 import numpy as np
-import pyproj
 import rasterio as rio
 import rasterio.warp
 import rasterio.windows
@@ -30,8 +29,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
+from scipy.ndimage import distance_transform_edt
 
 import geoutils.vector as gv
 from geoutils._config import config
@@ -52,6 +50,15 @@
     reproject_from_latlon,
 )
 from geoutils.raster.array import get_mask_from_array
+from geoutils.raster.georeferencing import (
+    _bounds,
+    _coords,
+    _ij2xy,
+    _outside_image,
+    _res,
+    _xy2ij,
+)
+from geoutils.raster.interpolate import _interp_points
 from geoutils.raster.sampling import subsample_array
 from geoutils.vector import Vector
 
@@ -788,12 +795,12 @@ def shape(self) -> tuple[int, int]:
     @property
     def res(self) -> tuple[float | int, float | int]:
         """Resolution (X, Y) of the raster in georeferenced units."""
-        return self.transform[0], abs(self.transform[4])
+        return _res(self.transform)
 
     @property
     def bounds(self) -> rio.coords.BoundingBox:
         """Bounding coordinates of the raster."""
-        return rio.coords.BoundingBox(*rio.transform.array_bounds(self.height, self.width, self.transform))
+        return _bounds(transform=self.transform, shape=self.shape)
 
     @property
     def footprint(self) -> Vector:
@@ -3515,60 +3522,15 @@ def xy2ij(
         :returns i, j: Indices of (x,y) in the image.
         """
 
-        # If undefined, default to the global system config
-        if shift_area_or_point is None:
-            shift_area_or_point = config["shift_area_or_point"]
-
-        # Input checks
-        if op not in [np.float32, np.float64, float]:
-            raise UserWarning(
-                "Operator is not of type float: rio.Dataset.index might "
-                "return unreliable indexes due to rounding issues."
-            )
-
-        i, j = rio.transform.rowcol(self.transform, x, y, op=op, precision=precision)
-
-        # Necessary because rio.Dataset.index does not return abc.Iterable for a single point
-        if not isinstance(i, abc.Iterable):
-            i, j = (
-                np.asarray(
-                    [
-                        i,
-                    ]
-                ),
-                np.asarray(
-                    [
-                        j,
-                    ]
-                ),
-            )
-        else:
-            i, j = (np.asarray(i), np.asarray(j))
-
-        # AREA_OR_POINT GDAL attribute, i.e. does the value refer to the upper left corner "Area" or
-        # the center of pixel "Point". This normally has no influence on georeferencing, it's only
-        # about the interpretation of the raster values, and thus can affect sub-pixel interpolation,
-        # for more details see: https://gdal.org/user/raster_data_model.html#metadata
-
-        # If the user wants to shift according to the interpretation
-        if shift_area_or_point:
-
-            # Shift by half a pixel if the AREA_OR_POINT attribute is "Point", otherwise leave as is
-            if self.area_or_point is not None and self.area_or_point == "Point":
-                if not isinstance(i.flat[0], (np.floating, float)):
-                    raise ValueError(
-                        "Operator must return np.floating values to perform pixel interpretation shifting."
-                    )
-
-                i += 0.5
-                j += 0.5
-
-        # Convert output indexes to integer if they are all whole numbers
-        if np.all(np.mod(i, 1) == 0) and np.all(np.mod(j, 1) == 0):
-            i = i.astype(int)
-            j = j.astype(int)
-
-        return i, j
+        return _xy2ij(
+            x=x,
+            y=y,
+            transform=self.transform,
+            area_or_point=self.area_or_point,
+            op=op,
+            precision=precision,
+            shift_area_or_point=shift_area_or_point,
+        )
 
     def ij2xy(
         self, i: ArrayLike, j: ArrayLike, shift_area_or_point: bool | None = None, force_offset: str | None = None
@@ -3593,23 +3555,14 @@ def ij2xy(
         :returns x, y: x,y coordinates of i,j in reference system.
         """
 
-        # If undefined, default to the global system config
-        if shift_area_or_point is None:
-            shift_area_or_point = config["shift_area_or_point"]
-
-        # Shift by half a pixel back for "Point" interpretation
-        if shift_area_or_point and force_offset is None:
-            if self.area_or_point is not None and self.area_or_point == "Point":
-                i = np.asarray(i) - 0.5
-                j = np.asarray(j) - 0.5
-
-        # Default offset is upper-left for raster coordinates
-        if force_offset is None:
-            force_offset = "ul"
-
-        x, y = rio.transform.xy(self.transform, i, j, offset=force_offset)
-
-        return x, y
+        return _ij2xy(
+            i=i,
+            j=j,
+            transform=self.transform,
+            area_or_point=self.area_or_point,
+            shift_area_or_point=shift_area_or_point,
+            force_offset=force_offset,
+        )
 
     def coords(
         self, grid: bool = True, shift_area_or_point: bool | None = None, force_offset: str | None = None
@@ -3627,23 +3580,14 @@ def coords(
         :returns x,y: Arrays of the (x,y) coordinates.
         """
 
-        # The coordinates are extracted from indexes 0 to shape
-        _, yy = self.ij2xy(
-            i=np.arange(self.height - 1, -1, -1),
-            j=0,
+        return _coords(
+            transform=self.transform,
+            shape=self.shape,
+            area_or_point=self.area_or_point,
+            grid=grid,
             shift_area_or_point=shift_area_or_point,
             force_offset=force_offset,
         )
-        xx, _ = self.ij2xy(
-            i=0, j=np.arange(self.width), shift_area_or_point=shift_area_or_point, force_offset=force_offset
-        )
-
-        # If grid is True, return coordinate grids
-        if grid:
-            meshgrid = tuple(np.meshgrid(xx, np.flip(yy)))
-            return meshgrid  # type: ignore
-        else:
-            return np.asarray(xx), np.asarray(yy)
 
     def outside_image(self, xi: ArrayLike, yj: ArrayLike, index: bool = True) -> bool:
         """
@@ -3655,20 +3599,16 @@ def outside_image(self, xi: ArrayLike, yj: ArrayLike, index: bool = True) -> boo
 
         :returns is_outside: ``True`` if ij is outside the image.
         """
-        if not index:
-            xi, xj = self.xy2ij(xi, yj)
 
-        if np.any(np.array((xi, yj)) < 0):
-            return True
-        elif np.asanyarray(xi) > self.width or np.asanyarray(yj) > self.height:
-            return True
-        else:
-            return False
+        return _outside_image(
+            xi=xi, yj=yj, transform=self.transform, shape=self.shape, area_or_point=self.area_or_point, index=index
+        )
 
     def interp_points(
         self,
         points: tuple[Number, Number] | tuple[NDArrayNum, NDArrayNum],
-        method: Literal["nearest", "linear", "cubic", "quintic"] = "linear",
+        method: Literal["nearest", "linear", "cubic", "quintic", "slinear", "pchip", "splinef2d"] = "linear",
+        dist_nodata_spread: Literal["half_order_up", "half_order_down"] | int = "half_order_up",
         band: int = 1,
         input_latlon: bool = False,
         shift_area_or_point: bool | None = None,
@@ -3678,8 +3618,8 @@ def interp_points(
         """
          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.
+         Uses scipy.ndimage.map_coordinates if the Raster is on an equal grid using "nearest" or "linear" (for speed),
+         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
@@ -3687,8 +3627,12 @@ def interp_points(
 
         :param points: Point(s) at which to interpolate raster value (tuple of X/Y array-likes). If points fall
             outside of image, value returned is nan.
-        :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 method: Interpolation method, one of 'nearest', 'linear', 'cubic', 'quintic', 'slinear', 'pchip' or
+            'splinef2d'. For more information, see scipy.ndimage.map_coordinates and scipy.interpolate.interpn.
+            Default is linear.
+        :param dist_nodata_spread: Distance of nodata spreading during interpolation, either half-interpolation order
+            rounded up (default; equivalent to 0 for nearest, 1 for linear methods, 2 for cubic methods and 3 for
+            quintic method), or rounded down, or a fixed integer.
         :param band: Band to use (from 1 to self.count).
         :param input_latlon: Whether the input is in latlon, unregarding of Raster CRS.
         :param shift_area_or_point: Whether to shift with pixel interpretation, which shifts to center of pixel
@@ -3699,64 +3643,26 @@ def interp_points(
         :returns rpoints: Array of raster value(s) for the given points.
         """
 
-        # Get coordinates
-        x, y = points
-
-        # If those are in latlon, convert to Raster CRS
-        if input_latlon:
-            init_crs = pyproj.CRS(4326)
-            dest_crs = pyproj.CRS(self.crs)
-            transformer = pyproj.Transformer.from_crs(init_crs, dest_crs)
-            x, y = transformer.transform(x, y)
-
-        i, j = self.xy2ij(x, y, shift_area_or_point=shift_area_or_point)
-
-        ind_invalid = np.vectorize(lambda k1, k2: self.outside_image(k1, k2, index=True))(j, i)
-
+        # Extract array supporting NaNs
         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:
-            # Get lower-left corner coordinates
-            xycoords = self.coords(grid=False, shift_area_or_point=shift_area_or_point)
-
-            # 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})
-
-            # Using direct coordinates, Y is the first axis, and we need to flip it
-            rpoints = interpn(
-                (np.flip(xycoords[1], axis=0), xycoords[0]), self.get_nanarray(), (y, x), method=method, **kwargs
-            )
-
-        rpoints = np.array(rpoints, dtype=np.float32)
-        rpoints[np.array(ind_invalid)] = np.nan
+        # If those are in latlon, convert to Raster CRS
+        if input_latlon:
+            points = reproject_from_latlon(points, out_crs=self.crs)  # type: ignore
 
-        return rpoints
+        return _interp_points(
+            array,
+            transform=self.transform,
+            area_or_point=self.area_or_point,
+            points=points,
+            method=method,
+            shift_area_or_point=shift_area_or_point,
+            dist_nodata_spread=dist_nodata_spread,
+            force_scipy_function=force_scipy_function,
+            **kwargs,
+        )
 
     def split_bands(self: RasterType, copy: bool = False, bands: list[int] | int | None = None) -> list[Raster]:
         """
diff --git a/geoutils/vector.py b/geoutils/vector.py
index 1e270be4..f4c61ab8 100644
--- a/geoutils/vector.py
+++ b/geoutils/vector.py
@@ -1138,7 +1138,7 @@ def reproject(
                     try:
                         ds_ref = Vector(ref)
                     except pyogrio.errors.DataSourceError:
-                        raise ValueError("Could not open raster or vector with rasterio or fiona.")
+                        raise ValueError("Could not open raster or vector with rasterio or pyogrio.")
             else:
                 raise TypeError("Type of ref must be string path to file, Raster or Vector.")
 
diff --git a/tests/test_raster/test_georeferencing.py b/tests/test_raster/test_georeferencing.py
new file mode 100644
index 00000000..fe2fcba9
--- /dev/null
+++ b/tests/test_raster/test_georeferencing.py
@@ -0,0 +1,211 @@
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+import geoutils as gu
+from geoutils import examples
+
+
+class TestGeoreferencing:
+
+    landsat_b4_path = examples.get_path("everest_landsat_b4")
+    aster_dem_path = examples.get_path("exploradores_aster_dem")
+    landsat_rgb_path = examples.get_path("everest_landsat_rgb")
+    landsat_b4_crop_path = examples.get_path("everest_landsat_b4_cropped")
+
+    @pytest.mark.parametrize("example", [landsat_b4_path, aster_dem_path, landsat_rgb_path])  # type: ignore
+    def test_ij2xy(self, example: str) -> None:
+        """Test the outputs of ij2xy and that the function is reversible with xy2ij."""
+
+        # Open raster
+        rst = gu.Raster(example)
+        xmin, ymin, xmax, ymax = rst.bounds
+
+        # Check ij2xy manually for the four corners and center
+        # With "force_offset", no considerations of pixel interpolation
+
+        # With offset="center", should be pixel center
+        xmin_center = xmin + rst.res[0] / 2
+        ymin_center = ymin + rst.res[1] / 2
+        xmax_center = xmax - rst.res[0] / 2
+        ymax_center = ymax - rst.res[1] / 2
+        assert rst.ij2xy([0], [0], force_offset="center") == ([xmin_center], [ymax_center])
+        assert rst.ij2xy([rst.shape[0] - 1], [0], force_offset="center") == ([xmin_center], [ymin_center])
+        assert rst.ij2xy([0], [rst.shape[1] - 1], force_offset="center") == ([xmax_center], [ymax_center])
+        assert rst.ij2xy([rst.shape[0] - 1], [rst.shape[1] - 1], force_offset="center") == (
+            [xmax_center],
+            [ymin_center],
+        )
+
+        # Same checks for offset="ll", "ul", "ur", "lr
+        lims_ll = [xmin, ymin, xmax - rst.res[0], ymax - rst.res[1]]
+        lims_ul = [xmin, ymin + rst.res[1], xmax - rst.res[0], ymax]
+        lims_lr = [xmin + rst.res[0], ymin, xmax, ymax - rst.res[1]]
+        lims_ur = [xmin + rst.res[0], ymin + rst.res[1], xmax, ymax]
+        offsets = ["ll", "ul", "lr", "ur"]
+        list_lims = [lims_ll, lims_ul, lims_lr, lims_ur]
+        for i in range(len(list_lims)):
+            offset = offsets[i]
+            lim = list_lims[i]
+            assert rst.ij2xy([0], [0], force_offset=offset) == ([lim[0]], [lim[3]])
+            assert rst.ij2xy([rst.shape[0] - 1], [0], force_offset=offset) == ([lim[0]], [lim[1]])
+            assert rst.ij2xy([0], [rst.shape[1] - 1], force_offset=offset) == ([lim[2]], [lim[3]])
+            assert rst.ij2xy([rst.shape[0] - 1], [rst.shape[1] - 1], force_offset=offset) == ([lim[2]], [lim[1]])
+
+        # Check that default coordinate is upper-left
+
+        # With shift from pixel interpretation, coordinates will be half a pixel back for "Point"
+        if rst.area_or_point is not None and rst.area_or_point == "Point" and gu.config["shift_area_or_point"]:
+            # Shift is backward in X, forward in Y
+            lims_ul[0] = lims_ul[0] - 0.5 * rst.res[0]
+            lims_ul[1] = lims_ul[1] + 0.5 * rst.res[1]
+            lims_ul[2] = lims_ul[2] - 0.5 * rst.res[0]
+            lims_ul[3] = lims_ul[3] + 0.5 * rst.res[1]
+
+        assert rst.ij2xy([0], [0]) == ([lims_ul[0]], [lims_ul[3]])
+        assert rst.ij2xy([rst.shape[0] - 1], [0]) == ([lims_ul[0]], [lims_ul[1]])
+        assert rst.ij2xy([0], [rst.shape[1] - 1]) == ([lims_ul[2]], [lims_ul[3]])
+        assert rst.ij2xy([rst.shape[0] - 1], [rst.shape[1] - 1]) == ([lims_ul[2]], [lims_ul[1]])
+
+    @pytest.mark.parametrize("example", [landsat_b4_path, aster_dem_path, landsat_rgb_path])  # type: ignore
+    def test_xy2ij_ij2xy_reversible(self, example: str):
+
+        # Open raster
+        rst = gu.Raster(example)
+        xmin, ymin, xmax, ymax = rst.bounds
+
+        # We generate random points within the boundaries of the image
+        rng = np.random.default_rng(42)
+        xrand = rng.integers(low=0, high=rst.width, size=(10,)) * list(rst.transform)[0] + xmin
+        yrand = ymax + rng.integers(low=0, high=rst.height, size=(10,)) * list(rst.transform)[4]
+
+        # Test reversibility for any point or area interpretation
+        i, j = rst.xy2ij(xrand, yrand)
+        xnew, ynew = rst.ij2xy(i, j)
+        assert all(xnew == xrand)
+        assert all(ynew == yrand)
+
+        # TODO: clarify this weird behaviour of rasterio.index with floats?
+        # r.ds.index(x, y)
+        # Out[33]: (75, 301)
+        # r.ds.index(x, y, op=np.float32)
+        # Out[34]: (75.0, 302.0)
+
+    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
+        r = gu.Raster(self.landsat_b4_path)
+        assert r.area_or_point == "Point"
+        xmin, ymin, xmax, ymax = r.bounds
+
+        # We generate random points within the boundaries of the image
+        rng = np.random.default_rng(42)
+        xrand = rng.integers(low=0, high=r.width, size=(10,)) * list(r.transform)[0] + xmin
+        yrand = ymax + rng.integers(low=0, high=r.height, size=(10,)) * list(r.transform)[4]
+
+        # Get decimal indexes based on "Point", should refer to the corner still (shift False by default)
+        i, j = r.xy2ij(xrand, yrand, shift_area_or_point=False)
+        assert np.all(i % 1 == 0)
+        assert np.all(j % 1 == 0)
+
+        # Those should all be .5 because values refer to the center and are shifted
+        i, j = r.xy2ij(xrand, yrand, shift_area_or_point=True)
+        assert np.all(i % 1 == 0.5)
+        assert np.all(j % 1 == 0.5)
+
+        # Force "Area", should refer to corner
+        r.set_area_or_point("Area", shift_area_or_point=False)
+        i, j = r.xy2ij(xrand, yrand, shift_area_or_point=True)
+        assert np.all(i % 1 == 0)
+        assert np.all(j % 1 == 0)
+
+        # Now, we calculate the mean of values in each 2x2 slices of the data, and compare with interpolation at order 1
+        list_z_ind = []
+        img = r.data
+        for k in range(len(xrand)):
+            # 2x2 slices
+            z_ind = np.mean(
+                img[
+                    slice(int(np.floor(i[k])), int(np.ceil(i[k])) + 1),
+                    slice(int(np.floor(j[k])), int(np.ceil(j[k])) + 1),
+                ]
+            )
+            list_z_ind.append(z_ind)
+
+        # First order interpolation
+        rpts = r.interp_points((xrand, yrand), method="linear")
+        # The values interpolated should be equal
+        assert np.array_equal(np.array(list_z_ind, dtype=np.float32), rpts, equal_nan=True)
+
+        # Test there is no failure with random coordinates (edge effects, etc)
+        xrand = rng.uniform(low=xmin, high=xmax, size=(1000,))
+        yrand = rng.uniform(low=ymin, high=ymax, size=(1000,))
+        r.interp_points((xrand, yrand))
+
+        # Second, test after a crop: the Raster now has an Area interpretation, those should fall right on the integer
+        # pixel indexes
+        r2 = gu.Raster(self.landsat_b4_crop_path)
+        r.crop(r2)
+        assert r.area_or_point == "Area"
+
+        xmin, ymin, xmax, ymax = r.bounds
+
+        # We can test with several method for the exact indexes: interp, and simple read should
+        # give back the same values that fall right on the coordinates
+        xrand = rng.integers(low=0, high=r.width, size=(10,)) * list(r.transform)[0] + xmin
+        yrand = ymax + rng.integers(low=0, high=r.height, size=(10,)) * list(r.transform)[4]
+        # By default, i and j are returned as integers
+        i, j = r.xy2ij(xrand, yrand, op=np.float32)
+        list_z_ind = []
+        img = r.data
+        for k in range(len(xrand)):
+            # We directly sample the values
+            z_ind = img[int(i[k]), int(j[k])]
+            list_z_ind.append(z_ind)
+
+        rpts = r.interp_points((xrand, yrand), method="linear")
+
+        assert np.array_equal(np.array(list_z_ind, dtype=np.float32), rpts, equal_nan=True)
+
+    @pytest.mark.parametrize("example", [landsat_b4_path, aster_dem_path, landsat_rgb_path])  # type: ignore
+    def test_coords(self, example: str) -> None:
+
+        img = gu.Raster(self.landsat_b4_path)
+
+        # With lower left argument
+        xx0, yy0 = img.coords(grid=False, force_offset="ll")
+
+        assert len(xx0) == img.width
+        assert len(yy0) == img.height
+
+        assert xx0[0] == pytest.approx(img.bounds.left)
+        assert xx0[-1] == pytest.approx(img.bounds.right - img.res[0])
+        if img.res[1] > 0:
+            assert yy0[0] == pytest.approx(img.bounds.bottom)
+            assert yy0[-1] == pytest.approx(img.bounds.top - img.res[1])
+        else:
+            # Currently not covered by test image
+            assert yy0[0] == pytest.approx(img.bounds.top)
+            assert yy0[-1] == pytest.approx(img.bounds.bottom + img.res[1])
+
+        # With center argument
+        xx, yy = img.coords(grid=False, force_offset="center")
+        hx = img.res[0] / 2
+        hy = img.res[1] / 2
+        assert xx[0] == pytest.approx(img.bounds.left + hx)
+        assert xx[-1] == pytest.approx(img.bounds.right - hx)
+        if img.res[1] > 0:
+            assert yy[0] == pytest.approx(img.bounds.bottom + hy)
+            assert yy[-1] == pytest.approx(img.bounds.top - hy)
+        else:
+            # Currently not covered by test image
+            assert yy[-1] == pytest.approx(img.bounds.top + hy)
+            assert yy[0] == pytest.approx(img.bounds.bottom - hy)
+
+        # With grid argument (default argument, repeated here for clarity)
+        xxgrid, yygrid = img.coords(grid=True, force_offset="ll")
+        assert np.array_equal(xxgrid, np.repeat(xx0[np.newaxis, :], img.height, axis=0))
+        assert np.array_equal(yygrid, np.flipud(np.repeat(yy0[:, np.newaxis], img.width, axis=1)))
diff --git a/tests/test_raster/test_interpolate.py b/tests/test_raster/test_interpolate.py
new file mode 100644
index 00000000..d2c743d7
--- /dev/null
+++ b/tests/test_raster/test_interpolate.py
@@ -0,0 +1,601 @@
+from __future__ import annotations
+
+import re
+from typing import Literal
+
+import numpy as np
+import pytest
+import rasterio as rio
+from scipy.interpolate import interpn
+from scipy.ndimage import binary_dilation
+
+import geoutils as gu
+from geoutils import examples
+from geoutils.projtools import reproject_to_latlon
+from geoutils.raster.interpolate import (
+    _get_dist_nodata_spread,
+    _interp_points,
+    _interpn_interpolator,
+    method_to_order,
+)
+
+
+class TestInterpolate:
+
+    landsat_b4_path = examples.get_path("everest_landsat_b4")
+    aster_dem_path = examples.get_path("exploradores_aster_dem")
+    landsat_b4_crop_path = examples.get_path("everest_landsat_b4_cropped")
+    landsat_rgb_path = examples.get_path("everest_landsat_rgb")
+
+    def test_dist_nodata_spread(self) -> None:
+        """Test distance of nodata spreading computation based on interpolation order."""
+
+        assert _get_dist_nodata_spread(0, "half_order_up") == 0
+        assert _get_dist_nodata_spread(0, "half_order_down") == 0
+        assert _get_dist_nodata_spread(0, 5) == 5
+
+        assert _get_dist_nodata_spread(1, "half_order_up") == 1
+        assert _get_dist_nodata_spread(1, "half_order_down") == 0
+        assert _get_dist_nodata_spread(1, 5) == 5
+
+        assert _get_dist_nodata_spread(3, "half_order_up") == 2
+        assert _get_dist_nodata_spread(3, "half_order_down") == 1
+        assert _get_dist_nodata_spread(3, 5) == 5
+
+        assert _get_dist_nodata_spread(5, "half_order_up") == 3
+        assert _get_dist_nodata_spread(5, "half_order_down") == 2
+        assert _get_dist_nodata_spread(5, 5) == 5
+
+    @pytest.mark.parametrize(
+        "method", ["nearest", "linear", "cubic", "quintic", "slinear", "pchip", "splinef2d"]
+    )  # type: ignore
+    def test_interpn_interpolator_accuracy(
+        self, method: Literal["nearest", "linear", "cubic", "quintic", "slinear", "pchip", "splinef2d"]
+    ):
+        """Test that _interpn_interpolator (used by interp_points) returns exactly the same result as scipy.interpn."""
+
+        # Create synthetic 2D array with non-aligned coordinates in X/Y, and X in descending order to mirror a raster's
+        # coordinates
+        shape = (50, 20)
+        coords = (np.linspace(10, 0, shape[0]), np.linspace(20, 30, shape[1]))
+        values = np.random.default_rng(42).normal(size=shape)
+
+        # Get 10 random points in the array boundaries
+        i = np.random.default_rng(42).uniform(coords[0][-1], coords[0][0], size=10)
+        j = np.random.default_rng(42).uniform(coords[1][0], coords[1][-1], size=10)
+
+        # Compare interpn and interpolator
+
+        # Method splinef2d is expecting strictly ascending coordinates (while other methods support desc or asc)
+        if method != "splinef2d":
+            vals = interpn(points=coords, values=values, xi=(i, j), method=method)
+        else:
+            vals = interpn(
+                points=(np.flip(coords[0]), coords[1]), values=np.flip(values[:], axis=0), xi=(i, j), method=method
+            )
+        # With the interpolator (coordinates are re-ordered automatically, as it happens often for rasters)
+        interpolator = _interpn_interpolator(points=coords, values=values, method=method)
+        vals2 = interpolator((i, j))
+
+        assert np.array_equal(vals, vals2, equal_nan=True)
+
+    @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 select known points and compare to the expected interpolation results across all methods, and all pixel
+        interpretations (area_or_point, with or without shift).
+
+        The synthetic data is a 3x3 array of values, which we interpolate at different synthetic points:
+        1/ Points falling right on grid coordinates are compared to their value in the grid,
+        2/ Points falling halfway between grid coordinates are compared to the mean of the surrounding values in the
+        grid (as it should equal the linear interpolation on an equal grid),
+        3/ Random points in the grid, compared between methods (forcing to use either scipy.ndimage.map_coordinates or
+        scipy.interpolate.interpn under-the-hood, to ensure results are consistent).
+
+        These tests also check the behaviour when returning interpolated for points outside of valid values.
+        """
+
+        # 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 without re-transforming
+        raster.set_area_or_point(tag_aop, shift_area_or_point=False)
+
+        # 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
+        points_x, points_y = raster.ij2xy(i=index_x, j=index_y, shift_area_or_point=shift_aop)
+
+        # 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
+
+        raster_points = raster.interp_points((points_x, points_y), method="nearest", shift_area_or_point=shift_aop)
+        raster_points_lin = raster.interp_points((points_x, points_y), method="linear", shift_area_or_point=shift_aop)
+        raster_points_interpn = raster.interp_points(
+            (points_x, points_y), method="nearest", force_scipy_function="interpn", shift_area_or_point=shift_aop
+        )
+        raster_points_interpn_lin = raster.interp_points(
+            (points_x, points_y), 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]
+
+        points_x_in, points_y_in = raster.ij2xy(i=index_x_in, j=index_y_in, shift_area_or_point=shift_aop)
+
+        # Here again compare methods
+        raster_points_in = raster.interp_points(
+            (points_x_in, points_y_in), method="linear", shift_area_or_point=shift_aop
+        )
+        raster_points_in_interpn = raster.interp_points(
+            (points_x_in, points_y_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_x_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)]
+        )
+        points_out_xy = list(zip(*points_out))
+        raster_points_out = raster.interp_points(points_out_xy)
+        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.set_area_or_point(tag_aop, shift_area_or_point=False)
+
+        # For this, get random points
+        rng = np.random.default_rng(42)
+        index_x_in_rand = rng.integers(low=8, high=42, size=(10,)) + rng.normal(scale=0.3)
+        index_y_in_rand = rng.integers(low=8, high=42, size=(10,)) + rng.normal(scale=0.3)
+        points_x_rand, points_y_rand = raster.ij2xy(i=index_x_in_rand, j=index_y_in_rand, shift_area_or_point=shift_aop)
+
+        for method in ["nearest", "linear"]:
+            raster_points_mapcoords = raster.interp_points(
+                (points_x_rand, points_y_rand),
+                method=method,
+                force_scipy_function="map_coordinates",
+                shift_area_or_point=shift_aop,
+            )
+            raster_points_interpn = raster.interp_points(
+                (points_x_rand, points_y_rand),
+                method=method,
+                force_scipy_function="interpn",
+                shift_area_or_point=shift_aop,
+            )
+
+            # Not exactly equal in floating point precision since changes in Scipy 1.13.0,
+            # see https://github.com/GlacioHack/geoutils/issues/533
+            assert np.allclose(raster_points_mapcoords, raster_points_interpn)
+
+        # Check that, outside the edge, the interpolation fails and returns a NaN
+        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]
+
+        points_x_rand, points_y_rand = raster.ij2xy(
+            i=index_x_edge_rand, j=index_y_edge_rand, shift_area_or_point=shift_aop
+        )
+
+        # Test across all methods
+        for method in ["nearest", "linear", "cubic", "quintic"]:
+            raster_points_mapcoords_edge = raster.interp_points(
+                (points_x_rand, points_y_rand),
+                method=method,
+                force_scipy_function="map_coordinates",
+                shift_area_or_point=shift_aop,
+            )
+            raster_points_interpn_edge = raster.interp_points(
+                (points_x_rand, points_y_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))
+
+    @pytest.mark.parametrize("example", [landsat_b4_path, aster_dem_path])  # type: ignore
+    @pytest.mark.parametrize(
+        "method", ["nearest", "linear", "cubic", "quintic", "slinear", "pchip", "splinef2d"]
+    )  # type: ignore
+    def test_interp_points__real(
+        self, example: str, method: Literal["nearest", "linear", "cubic", "quintic", "slinear", "pchip", "splinef2d"]
+    ) -> None:
+        """
+        Test interp_points for real data, checking in particular the propagation of nodata.
+
+        For a random point (dimension 0) and a group of random points (dimension 1) in a real raster, we check the
+        consistency of the output forcing to use either scipy.ndimage.map_coordinates or scipy.interpolate.interpn
+        under-the-hood, or returning a regular-grid interpolator.
+        """
+
+        # Check the accuracy of the interpolation at an exact point, and between methods
+
+        # Open and crop for speed
+        r = gu.Raster(example)
+        r = r.crop((r.bounds.left, r.bounds.bottom, r.bounds.left + r.res[0] * 50, r.bounds.bottom + r.res[1] * 50))
+        r.set_area_or_point("Area", shift_area_or_point=False)
+
+        # 1/ Test for an individual point (shape can be tricky in 1 dimension)
+        itest = 10
+        jtest = 10
+        x, y = r.ij2xy(itest, jtest)
+        val = r.interp_points((x, y), method=method, force_scipy_function="map_coordinates")[0]
+        val_img = r.data[itest, jtest]
+        # For a point exactly at a grid coordinate, only nearest and linear will match
+        # (cubic modifies values at a grid coordinate)
+        if method in ["nearest", "linear"]:
+            assert val_img == val
+
+        # Check the result is exactly the same for both methods
+        val2 = r.interp_points((x, y), method=method, force_scipy_function="interpn")[0]
+        assert val2 == pytest.approx(val)
+
+        # 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), method=method, input_latlon=True)[0]
+        assert val == pytest.approx(val_latlon, abs=0.0001)
+
+        # 2/ Test for multiple points
+        i = np.random.default_rng(42).integers(1, 49, size=10)
+        j = np.random.default_rng(42).integers(1, 49, size=10)
+        x, y = r.ij2xy(i, j)
+        vals = r.interp_points((x, y), method=method, force_scipy_function="map_coordinates")
+        vals2 = r.interp_points((x, y), method=method, force_scipy_function="interpn")
+
+        assert np.array_equal(vals, vals2, equal_nan=True)
+
+        # 3/ Test return_interpolator is consistent with above
+        interp = _interp_points(
+            r.get_nanarray(),
+            transform=r.transform,
+            area_or_point=r.area_or_point,
+            points=(x, y),
+            method=method,
+            return_interpolator=True,
+        )
+        vals3 = interp((y, x))
+        vals3 = np.array(np.atleast_1d(vals3), dtype=np.float32)
+
+        assert np.array_equal(vals2, vals3, equal_nan=True)
+
+    @pytest.mark.parametrize("example", [landsat_b4_path, aster_dem_path])  # type: ignore
+    @pytest.mark.parametrize(
+        "method", ["nearest", "linear", "cubic", "quintic", "slinear", "pchip", "splinef2d"]
+    )  # type: ignore
+    @pytest.mark.parametrize("dist", ["half_order_up", "half_order_down", 0, 1, 5])  # type: ignore
+    def test_interp_point__nodata_propag(
+        self,
+        example: str,
+        method: Literal["nearest", "linear", "cubic", "quintic", "slinear", "pchip", "splinef2d"],
+        dist: Literal["half_order_up", "half_order_down"] | int,
+    ) -> None:
+        """
+        Tests for nodata propagation of interp_points.
+
+        We create artificial NaNs at certain pixels of real rasters (one with integer type, one floating type), then
+        verify that interpolated values propagate these NaNs at the right distances for all methods, and across all
+        nodata propagation distances.
+
+        We also assess the precision of the interpolation under the default nodata propagation "half_order_up", by
+        comparing interpolated values near nodata edges with and without the gap-filling of a constant placeholder
+        value (= ensures the gap-filling of the array required for most methods that do not support NaN has little
+        effect of interpolation).
+        """
+        # Open and crop for speed
+        r = gu.Raster(example)
+        r = r.crop((r.bounds.left, r.bounds.bottom, r.bounds.left + r.res[0] * 100, r.bounds.bottom + r.res[1] * 100))
+
+        # 1/ Check the propagation of NaNs
+
+        # 1.1/ Manual check for a specific point
+        # Convert raster to float
+        r = r.astype(np.float32)
+
+        # Create a NaN at a given pixel (we know the landsat example has no NaNs to begin with)
+        i0, j0 = (10, 10)
+        r[i0, j0] = np.nan
+
+        # Create a big NaN area in the bottom right (for more complex NaN propagation below)
+        r[80:90, 80:90] = np.nan
+
+        # All surrounding pixels with distance half the method order rounded up should be NaNs
+        order = method_to_order[method]
+        d = _get_dist_nodata_spread(order=order, dist_nodata_spread=dist)
+
+        # Get indices of raster pixels within the right distance from NaNs
+        indices_nan = [
+            (i0 + i, j0 + j) for i in np.arange(-d, d + 1) for j in np.arange(-d, d + 1) if (np.abs(i) + np.abs(j)) <= d
+        ]
+        i, j = list(zip(*indices_nan))
+        x, y = r.ij2xy(i, j)
+        vals = r.interp_points((x, y), method=method, force_scipy_function="map_coordinates", dist_nodata_spread=dist)
+        vals2 = r.interp_points((x, y), method=method, force_scipy_function="interpn", dist_nodata_spread=dist)
+
+        assert all(np.isnan(np.atleast_1d(vals))) and all(np.isnan(np.atleast_1d(vals2)))
+
+        # Same check forrandom coordinates within half a pixel of the above coordinates (falling exactly on grid points)
+        xoffset = np.random.default_rng(42).uniform(low=-0.5, high=0.5, size=len(x))
+        yoffset = np.random.default_rng(42).uniform(low=-0.5, high=0.5, size=len(x))
+
+        vals = r.interp_points(
+            (x + xoffset, y + yoffset), method=method, force_scipy_function="map_coordinates", dist_nodata_spread=dist
+        )
+        vals2 = r.interp_points(
+            (x + xoffset, y + yoffset), method=method, force_scipy_function="interpn", dist_nodata_spread=dist
+        )
+
+        assert all(np.isnan(np.atleast_1d(vals))) and all(np.isnan(np.atleast_1d(vals2)))
+
+        # 1.2/ Check for all NaNs in the raster
+
+        # We create the mask of dilated NaNs
+        mask_nan = ~np.isfinite(r.get_nanarray())
+        if d != 0:
+            mask_nan_dilated = binary_dilation(mask_nan, iterations=d).astype("uint8")
+        # (Zero iteration triggers a different behaviour than just "doing nothing" in binary_dilation, we override here)
+        else:
+            mask_nan_dilated = mask_nan.astype("uint8")
+        # Get indices of the related pixels, convert to coordinates
+        i, j = np.where(mask_nan_dilated)
+        x, y = r.ij2xy(i, j)
+        # And interpolate at those coordinates
+        vals = r.interp_points((x, y), method=method, force_scipy_function="map_coordinates", dist_nodata_spread=dist)
+        vals2 = r.interp_points((x, y), method=method, force_scipy_function="interpn", dist_nodata_spread=dist)
+
+        assert all(np.isnan(np.atleast_1d(vals))) and all(np.isnan(np.atleast_1d(vals2)))
+
+        # 2/ Check that interpolated values at the edge of NaNs are valid + have small errors due to filling NaNs
+        # with the nearest value during interpolation (thanks to the spreading of the nodata mask)
+
+        # We compare values interpolated right at the edge of valid values near a NaN between
+        # a/ Implementation of interp_points (that replaces NaNs by the nearest neighbour during interpolation)
+        # b/ Raster filled with placeholder value then running interp_points
+
+        # 2.1/ Manual check for a specific point
+
+        # We get the indexes of valid pixels just at the edge of NaNs
+        indices_edge = [
+            (i0 + i, j0 + j)
+            for i in np.arange(-d - 1, d + 2)
+            for j in np.arange(-d - 1, d + 2)
+            if (np.abs(i) + np.abs(j)) == d + 1
+        ]
+        i, j = list(zip(*indices_edge))
+        x, y = r.ij2xy(i, j)
+        # And get their interpolated value
+        vals = r.interp_points((x, y), method=method, force_scipy_function="map_coordinates", dist_nodata_spread=dist)
+        vals2 = r.interp_points((x, y), method=method, force_scipy_function="interpn", dist_nodata_spread=dist)
+
+        # Then we fill the NaNs in the raster with a placeholder value of the raster mean
+        r_arr = r.get_nanarray()
+        r_arr[~np.isfinite(r_arr)] = np.nanmean(r_arr)
+        r2 = r.copy(new_array=r_arr)
+
+        # All raster values should be valid now
+        assert np.all(np.isfinite(r_arr))
+
+        # Only check the accuracy with the default NaN spreading (half-order rounded up), otherwise anything can happen
+        if dist == "half_order_up":
+
+            # Get the interpolated values
+            vals_near = r2.interp_points(
+                (x, y), method=method, force_scipy_function="map_coordinates", dist_nodata_spread=dist
+            )
+            vals2_near = r2.interp_points(
+                (x, y), method=method, force_scipy_function="interpn", dist_nodata_spread=dist
+            )
+
+            # Both sets of values should be valid + within a relative tolerance of 0.01%
+            assert np.allclose(vals, vals_near, equal_nan=False, rtol=10e-4)
+            assert np.allclose(vals2, vals2_near, equal_nan=False, rtol=10e-4)
+
+            # Same check for values within one pixel of the exact coordinates
+            xoffset = np.random.default_rng(42).uniform(low=-0.5, high=0.5, size=len(x))
+            yoffset = np.random.default_rng(42).uniform(low=-0.5, high=0.5, size=len(x))
+
+            vals = r.interp_points(
+                (x + xoffset, y + yoffset),
+                method=method,
+                force_scipy_function="map_coordinates",
+                dist_nodata_spread=dist,
+            )
+            vals2 = r.interp_points(
+                (x + xoffset, y + yoffset), method=method, force_scipy_function="interpn", dist_nodata_spread=dist
+            )
+            vals_near = r2.interp_points(
+                (x + xoffset, y + yoffset),
+                method=method,
+                force_scipy_function="map_coordinates",
+                dist_nodata_spread=dist,
+            )
+            vals2_near = r2.interp_points(
+                (x + xoffset, y + yoffset), method=method, force_scipy_function="interpn", dist_nodata_spread=dist
+            )
+
+            # Both sets of values should be exactly the same, without any NaNs
+            assert np.allclose(vals, vals_near, equal_nan=False, rtol=10e-4)
+            assert np.allclose(vals2, vals2_near, equal_nan=False, rtol=10e-4)
+
+            # 2.2/ Repeat the same for all edges of NaNs in the raster
+            mask_dilated_plus_one = binary_dilation(mask_nan_dilated, iterations=1).astype(bool)
+            mask_edge_dilated = np.logical_and(mask_dilated_plus_one, ~mask_nan_dilated.astype(bool))
+
+            # Get indices of the related pixels, convert to coordinates
+            i, j = np.where(mask_edge_dilated)
+            x, y = r.ij2xy(i, j)
+            # And interpolate at those coordinates
+            vals = r.interp_points(
+                (x, y), method=method, force_scipy_function="map_coordinates", dist_nodata_spread=dist
+            )
+            vals2 = r.interp_points((x, y), method=method, force_scipy_function="interpn", dist_nodata_spread=dist)
+            vals_near = r2.interp_points(
+                (x, y), method=method, force_scipy_function="map_coordinates", dist_nodata_spread=dist
+            )
+            vals2_near = r2.interp_points(
+                (x, y), method=method, force_scipy_function="interpn", dist_nodata_spread=dist
+            )
+
+            # Both sets of values should be exactly the same, without any NaNs
+            assert np.allclose(vals, vals_near, equal_nan=False, rtol=10e-4)
+            assert np.allclose(vals2, vals2_near, equal_nan=False, rtol=10e-4)
+
+    def test_reduce_points(self) -> None:
+        """
+        Test reduce points.
+        """
+
+        # -- Tests 1: Check based on indexed values --
+
+        # Open raster
+        r = gu.Raster(self.landsat_b4_crop_path)
+
+        # A pixel center where all neighbouring coordinates are different:
+        # array([[[237, 194, 239],
+        #          [250, 173, 164],
+        #          [255, 192, 128]]]
+        itest0 = 120
+        jtest0 = 451
+        # This is the center of the pixel
+        xtest0 = 496975.0
+        ytest0 = 3099095.0
+
+        # Verify coordinates match indexes
+        x_out, y_out = r.ij2xy(itest0, jtest0, force_offset="center")
+        assert x_out == xtest0
+        assert y_out == ytest0
+
+        # Check that the value at this coordinate is the same as when indexing
+        z_val = r.reduce_points((xtest0, ytest0))
+        z = r.data.data[itest0, jtest0]
+        assert z == z_val
+
+        # Check that the value is the same the other 4 corners of the pixel
+        assert z == r.reduce_points((xtest0 + 0.49 * r.res[0], ytest0 - 0.49 * r.res[1]))
+        assert z == r.reduce_points((xtest0 - 0.49 * r.res[0], ytest0 + 0.49 * r.res[1]))
+        assert z == r.reduce_points((xtest0 - 0.49 * r.res[0], ytest0 - 0.49 * r.res[1]))
+        assert z == r.reduce_points((xtest0 + 0.49 * r.res[0], ytest0 + 0.49 * r.res[1]))
+
+        # -- Tests 2: check arguments work as intended --
+
+        # 1/ Lat-lon argument check by getting the coordinates of our last test point
+        lat, lon = reproject_to_latlon(points=[xtest0, ytest0], in_crs=r.crs)
+        z_val_2 = r.reduce_points((lon, lat), input_latlon=True)
+        assert z_val == z_val_2
+
+        # 2/ Band argument
+        # Get the band indexes for the multi-band Raster
+        r_multi = gu.Raster(self.landsat_rgb_path)
+        itest, jtest = r_multi.xy2ij(xtest0, ytest0)
+        itest = int(itest[0])
+        jtest = int(jtest[0])
+        # Extract the values
+        z_band1 = r_multi.reduce_points((xtest0, ytest0), band=1)
+        z_band2 = r_multi.reduce_points((xtest0, ytest0), band=2)
+        z_band3 = r_multi.reduce_points((xtest0, ytest0), band=3)
+        # Compare to the Raster array slice
+        assert list(r_multi.data[:, itest, jtest]) == [z_band1, z_band2, z_band3]
+
+        # 3/ Masked argument
+        r_multi.data[:, itest, jtest] = np.ma.masked
+        z_not_ma = r_multi.reduce_points((xtest0, ytest0), band=1)
+        assert not np.ma.is_masked(z_not_ma)
+        z_ma = r_multi.reduce_points((xtest0, ytest0), band=1, masked=True)
+        assert np.ma.is_masked(z_ma)
+
+        # 4/ Window argument
+        val_window, z_window = r_multi.reduce_points(
+            (xtest0, ytest0), band=1, window=3, masked=True, return_window=True
+        )
+        assert (
+            val_window
+            == np.ma.mean(r_multi.data[0, itest - 1 : itest + 2, jtest - 1 : jtest + 2])
+            == np.ma.mean(z_window)
+        )
+        assert np.array_equal(z_window, r_multi.data[0, itest - 1 : itest + 2, jtest - 1 : jtest + 2])
+
+        # 5/ Reducer function argument
+        val_window2 = r_multi.reduce_points(
+            (xtest0, ytest0), band=1, window=3, masked=True, reducer_function=np.ma.median
+        )
+        assert val_window2 == np.ma.median(r_multi.data[0, itest - 1 : itest + 2, jtest - 1 : jtest + 2])
+
+        # -- Tests 3: check that errors are raised when supposed for non-boolean arguments --
+
+        # Verify that passing a window that is not a whole number fails
+        with pytest.raises(ValueError, match=re.escape("Window must be a whole number.")):
+            r.reduce_points((xtest0, ytest0), window=3.5)  # type: ignore
+        # Same for an odd number
+        with pytest.raises(ValueError, match=re.escape("Window must be an odd number.")):
+            r.reduce_points((xtest0, ytest0), window=4)
+        # But a window that is a whole number as a float works
+        r.reduce_points((xtest0, ytest0), window=3.0)  # type: ignore
+
+        # -- Tests 4: check that passing an array-like object works
+
+        # For simple coordinates
+        x_coords = [xtest0, xtest0 + 100]
+        y_coords = [ytest0, ytest0 - 100]
+        vals = r_multi.reduce_points((x_coords, y_coords))
+        val0, win0 = r_multi.reduce_points((x_coords[0], y_coords[0]), return_window=True)
+        val1, win1 = r_multi.reduce_points((x_coords[1], y_coords[1]), return_window=True)
+
+        assert len(vals) == len(x_coords)
+        assert np.array_equal(vals[0], val0, equal_nan=True)
+        assert np.array_equal(vals[1], val1, equal_nan=True)
+
+        # With a return window argument
+        vals, windows = r_multi.reduce_points((x_coords, y_coords), return_window=True)
+        assert len(windows) == len(x_coords)
+        assert np.array_equal(windows[0], win0, equal_nan=True)
+        assert np.array_equal(windows[1], win1, equal_nan=True)
+
+        # -- Tests 5 -- Check image corners and latlon argument
+
+        # Lower right pixel
+        x, y = [r.bounds.right - r.res[0] / 2, r.bounds.bottom + r.res[1] / 2]
+        lat, lon = reproject_to_latlon([x, y], r.crs)
+        assert r.reduce_points((x, y)) == r.reduce_points((lon, lat), input_latlon=True) == r.data[-1, -1]
+
+        # One pixel above
+        x, y = [r.bounds.right - r.res[0] / 2, r.bounds.bottom + 3 * r.res[1] / 2]
+        lat, lon = reproject_to_latlon([x, y], r.crs)
+        assert r.reduce_points((x, y)) == r.reduce_points((lon, lat), input_latlon=True) == r.data[-2, -1]
+
+        # One pixel left
+        x, y = [r.bounds.right - 3 * r.res[0] / 2, r.bounds.bottom + r.res[1] / 2]
+        lat, lon = reproject_to_latlon([x, y], r.crs)
+        assert r.reduce_points((x, y)) == r.reduce_points((lon, lat), input_latlon=True) == r.data[-1, -2]
diff --git a/tests/test_raster/test_raster.py b/tests/test_raster/test_raster.py
index da6334b2..ed16c0ae 100644
--- a/tests/test_raster/test_raster.py
+++ b/tests/test_raster/test_raster.py
@@ -20,11 +20,9 @@
 from pylint.reporters.text import TextReporter
 
 import geoutils as gu
-import geoutils.projtools as pt
 from geoutils import examples
 from geoutils._typing import MArrayNum, NDArrayNum
 from geoutils.misc import resampling_method_from_str
-from geoutils.projtools import reproject_to_latlon
 from geoutils.raster.raster import _default_nodata, _default_rio_attrs
 
 DO_PLOT = False
@@ -1981,447 +1979,6 @@ def test_intersection(self, example: list[str]) -> None:
             intersection = r.intersection(r_nonoverlap)
             assert intersection == (0.0, 0.0, 0.0, 0.0)
 
-    @pytest.mark.parametrize("example", [landsat_b4_path, aster_dem_path, landsat_rgb_path])  # type: ignore
-    def test_ij2xy(self, example: str) -> None:
-        """Test the outputs of ij2xy and that the function is reversible with xy2ij."""
-
-        # Open raster
-        rst = gu.Raster(example)
-        xmin, ymin, xmax, ymax = rst.bounds
-
-        # Check ij2xy manually for the four corners and center
-        # With "force_offset", no considerations of pixel interpolation
-
-        # With offset="center", should be pixel center
-        xmin_center = xmin + rst.res[0] / 2
-        ymin_center = ymin + rst.res[1] / 2
-        xmax_center = xmax - rst.res[0] / 2
-        ymax_center = ymax - rst.res[1] / 2
-        assert rst.ij2xy([0], [0], force_offset="center") == ([xmin_center], [ymax_center])
-        assert rst.ij2xy([rst.shape[0] - 1], [0], force_offset="center") == ([xmin_center], [ymin_center])
-        assert rst.ij2xy([0], [rst.shape[1] - 1], force_offset="center") == ([xmax_center], [ymax_center])
-        assert rst.ij2xy([rst.shape[0] - 1], [rst.shape[1] - 1], force_offset="center") == (
-            [xmax_center],
-            [ymin_center],
-        )
-
-        # Same checks for offset="ll", "ul", "ur", "lr
-        lims_ll = [xmin, ymin, xmax - rst.res[0], ymax - rst.res[1]]
-        lims_ul = [xmin, ymin + rst.res[1], xmax - rst.res[0], ymax]
-        lims_lr = [xmin + rst.res[0], ymin, xmax, ymax - rst.res[1]]
-        lims_ur = [xmin + rst.res[0], ymin + rst.res[1], xmax, ymax]
-        offsets = ["ll", "ul", "lr", "ur"]
-        list_lims = [lims_ll, lims_ul, lims_lr, lims_ur]
-        for i in range(len(list_lims)):
-            offset = offsets[i]
-            lim = list_lims[i]
-            assert rst.ij2xy([0], [0], force_offset=offset) == ([lim[0]], [lim[3]])
-            assert rst.ij2xy([rst.shape[0] - 1], [0], force_offset=offset) == ([lim[0]], [lim[1]])
-            assert rst.ij2xy([0], [rst.shape[1] - 1], force_offset=offset) == ([lim[2]], [lim[3]])
-            assert rst.ij2xy([rst.shape[0] - 1], [rst.shape[1] - 1], force_offset=offset) == ([lim[2]], [lim[1]])
-
-        # Check that default coordinate is upper-left
-
-        # With shift from pixel interpretation, coordinates will be half a pixel back for "Point"
-        if rst.area_or_point is not None and rst.area_or_point == "Point" and gu.config["shift_area_or_point"]:
-            # Shift is backward in X, forward in Y
-            lims_ul[0] = lims_ul[0] - 0.5 * rst.res[0]
-            lims_ul[1] = lims_ul[1] + 0.5 * rst.res[1]
-            lims_ul[2] = lims_ul[2] - 0.5 * rst.res[0]
-            lims_ul[3] = lims_ul[3] + 0.5 * rst.res[1]
-
-        assert rst.ij2xy([0], [0]) == ([lims_ul[0]], [lims_ul[3]])
-        assert rst.ij2xy([rst.shape[0] - 1], [0]) == ([lims_ul[0]], [lims_ul[1]])
-        assert rst.ij2xy([0], [rst.shape[1] - 1]) == ([lims_ul[2]], [lims_ul[3]])
-        assert rst.ij2xy([rst.shape[0] - 1], [rst.shape[1] - 1]) == ([lims_ul[2]], [lims_ul[1]])
-
-    @pytest.mark.parametrize("example", [landsat_b4_path, aster_dem_path, landsat_rgb_path])  # type: ignore
-    def test_xy2ij_ij2xy_reversible(self, example: str):
-
-        # Open raster
-        rst = gu.Raster(example)
-        xmin, ymin, xmax, ymax = rst.bounds
-
-        # We generate random points within the boundaries of the image
-        rng = np.random.default_rng(42)
-        xrand = rng.integers(low=0, high=rst.width, size=(10,)) * list(rst.transform)[0] + xmin
-        yrand = ymax + rng.integers(low=0, high=rst.height, size=(10,)) * list(rst.transform)[4]
-
-        # Test reversibility for any point or area interpretation
-        i, j = rst.xy2ij(xrand, yrand)
-        xnew, ynew = rst.ij2xy(i, j)
-        assert all(xnew == xrand)
-        assert all(ynew == yrand)
-
-        # TODO: clarify this weird behaviour of rasterio.index with floats?
-        # r.ds.index(x, y)
-        # Out[33]: (75, 301)
-        # r.ds.index(x, y, op=np.float32)
-        # Out[34]: (75.0, 302.0)
-
-    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
-        r = gu.Raster(self.landsat_b4_path)
-        assert r.area_or_point == "Point"
-        xmin, ymin, xmax, ymax = r.bounds
-
-        # We generate random points within the boundaries of the image
-        rng = np.random.default_rng(42)
-        xrand = rng.integers(low=0, high=r.width, size=(10,)) * list(r.transform)[0] + xmin
-        yrand = ymax + rng.integers(low=0, high=r.height, size=(10,)) * list(r.transform)[4]
-
-        # Get decimal indexes based on "Point", should refer to the corner still (shift False by default)
-        i, j = r.xy2ij(xrand, yrand, shift_area_or_point=False)
-        assert np.all(i % 1 == 0)
-        assert np.all(j % 1 == 0)
-
-        # Those should all be .5 because values refer to the center and are shifted
-        i, j = r.xy2ij(xrand, yrand, shift_area_or_point=True)
-        assert np.all(i % 1 == 0.5)
-        assert np.all(j % 1 == 0.5)
-
-        # Force "Area", should refer to corner
-        r.set_area_or_point("Area", shift_area_or_point=False)
-        i, j = r.xy2ij(xrand, yrand, shift_area_or_point=True)
-        assert np.all(i % 1 == 0)
-        assert np.all(j % 1 == 0)
-
-        # Now, we calculate the mean of values in each 2x2 slices of the data, and compare with interpolation at order 1
-        list_z_ind = []
-        img = r.data
-        for k in range(len(xrand)):
-            # 2x2 slices
-            z_ind = np.mean(
-                img[
-                    slice(int(np.floor(i[k])), int(np.ceil(i[k])) + 1),
-                    slice(int(np.floor(j[k])), int(np.ceil(j[k])) + 1),
-                ]
-            )
-            list_z_ind.append(z_ind)
-
-        # First order interpolation
-        rpts = r.interp_points((xrand, yrand), method="linear")
-        # The values interpolated should be equal
-        assert np.array_equal(np.array(list_z_ind, dtype=np.float32), rpts, equal_nan=True)
-
-        # Test there is no failure with random coordinates (edge effects, etc)
-        xrand = rng.uniform(low=xmin, high=xmax, size=(1000,))
-        yrand = rng.uniform(low=ymin, high=ymax, size=(1000,))
-        r.interp_points((xrand, yrand))
-
-        # Second, test after a crop: the Raster now has an Area interpretation, those should fall right on the integer
-        # pixel indexes
-        r2 = gu.Raster(self.landsat_b4_crop_path)
-        r.crop(r2)
-        assert r.area_or_point == "Area"
-
-        xmin, ymin, xmax, ymax = r.bounds
-
-        # We can test with several method for the exact indexes: interp, and simple read should
-        # give back the same values that fall right on the coordinates
-        xrand = rng.integers(low=0, high=r.width, size=(10,)) * list(r.transform)[0] + xmin
-        yrand = ymax + rng.integers(low=0, high=r.height, size=(10,)) * list(r.transform)[4]
-        # By default, i and j are returned as integers
-        i, j = r.xy2ij(xrand, yrand, op=np.float32)
-        list_z_ind = []
-        img = r.data
-        for k in range(len(xrand)):
-            # We directly sample the values
-            z_ind = img[int(i[k]), int(j[k])]
-            list_z_ind.append(z_ind)
-
-        rpts = r.interp_points((xrand, yrand), method="linear")
-
-        assert np.array_equal(np.array(list_z_ind, dtype=np.float32), rpts, equal_nan=True)
-
-    @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 without re-transforming
-        raster.set_area_or_point(tag_aop, shift_area_or_point=False)
-
-        # 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
-        points_x, points_y = raster.ij2xy(i=index_x, j=index_y, shift_area_or_point=shift_aop)
-
-        # 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
-
-        raster_points = raster.interp_points((points_x, points_y), method="nearest", shift_area_or_point=shift_aop)
-        raster_points_lin = raster.interp_points((points_x, points_y), method="linear", shift_area_or_point=shift_aop)
-        raster_points_interpn = raster.interp_points(
-            (points_x, points_y), method="nearest", force_scipy_function="interpn", shift_area_or_point=shift_aop
-        )
-        raster_points_interpn_lin = raster.interp_points(
-            (points_x, points_y), 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]
-
-        points_x_in, points_y_in = raster.ij2xy(i=index_x_in, j=index_y_in, shift_area_or_point=shift_aop)
-
-        # Here again compare methods
-        raster_points_in = raster.interp_points(
-            (points_x_in, points_y_in), method="linear", shift_area_or_point=shift_aop
-        )
-        raster_points_in_interpn = raster.interp_points(
-            (points_x_in, points_y_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_x_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)]
-        )
-        points_out_xy = list(zip(*points_out))
-        raster_points_out = raster.interp_points(points_out_xy)
-        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.set_area_or_point(tag_aop, shift_area_or_point=False)
-
-        # For this, get random points
-        rng = np.random.default_rng(42)
-        index_x_in_rand = rng.integers(low=8, high=42, size=(10,)) + rng.normal(scale=0.3)
-        index_y_in_rand = rng.integers(low=8, high=42, size=(10,)) + rng.normal(scale=0.3)
-        points_x_rand, points_y_rand = raster.ij2xy(i=index_x_in_rand, j=index_y_in_rand, shift_area_or_point=shift_aop)
-
-        for method in ["nearest", "linear", "cubic", "quintic"]:
-            raster_points_mapcoords = raster.interp_points(
-                (points_x_rand, points_y_rand),
-                method=method,
-                force_scipy_function="map_coordinates",
-                shift_area_or_point=shift_aop,
-            )
-            raster_points_interpn = raster.interp_points(
-                (points_x_rand, points_y_rand),
-                method=method,
-                force_scipy_function="interpn",
-                shift_area_or_point=shift_aop,
-            )
-
-            # Not exactly equal in floating point precision since changes in Scipy 1.13.0,
-            # see https://github.com/GlacioHack/geoutils/issues/533
-            assert np.allclose(raster_points_mapcoords, raster_points_interpn, atol=0.01)
-
-        # Check that, outside the edge, the interpolation fails and returns a NaN
-        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]
-
-        points_x_rand, points_y_rand = raster.ij2xy(
-            i=index_x_edge_rand, j=index_y_edge_rand, shift_area_or_point=shift_aop
-        )
-
-        # Nearest doesn't apply, just linear and above
-        for method in ["cubic", "quintic"]:
-            raster_points_mapcoords_edge = raster.interp_points(
-                (points_x_rand, points_y_rand),
-                method=method,
-                force_scipy_function="map_coordinates",
-                shift_area_or_point=shift_aop,
-            )
-            raster_points_interpn_edge = raster.interp_points(
-                (points_x_rand, points_y_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_interp_points__real(self) -> None:
-        """Test interp_points for real data."""
-
-        r = gu.Raster(self.landsat_b4_path)
-        r.set_area_or_point("Area", shift_area_or_point=False)
-
-        # 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), method="linear")[0]
-        val_img = r.data[int(i[0]), int(j[0])]
-        assert val_img == val
-
-        # Check the result is exactly the same for both methods
-        val2 = r.interp_points((x, y), method="linear", force_scipy_function="interpn")[0]
-        assert val2 == pytest.approx(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), method="linear", input_latlon=True)[0]
-        assert val == pytest.approx(val_latlon, abs=0.0001)
-
-    def test_reduce_points(self) -> None:
-        """
-        Test reduce points.
-        """
-
-        # -- Tests 1: Check based on indexed values --
-
-        # Open raster
-        r = gu.Raster(self.landsat_b4_crop_path)
-
-        # A pixel center where all neighbouring coordinates are different:
-        # array([[[237, 194, 239],
-        #          [250, 173, 164],
-        #          [255, 192, 128]]]
-        itest0 = 120
-        jtest0 = 451
-        # This is the center of the pixel
-        xtest0 = 496975.0
-        ytest0 = 3099095.0
-
-        # Verify coordinates match indexes
-        x_out, y_out = r.ij2xy(itest0, jtest0, force_offset="center")
-        assert x_out == xtest0
-        assert y_out == ytest0
-
-        # Check that the value at this coordinate is the same as when indexing
-        z_val = r.reduce_points((xtest0, ytest0))
-        z = r.data.data[itest0, jtest0]
-        assert z == z_val
-
-        # Check that the value is the same the other 4 corners of the pixel
-        assert z == r.reduce_points((xtest0 + 0.49 * r.res[0], ytest0 - 0.49 * r.res[1]))
-        assert z == r.reduce_points((xtest0 - 0.49 * r.res[0], ytest0 + 0.49 * r.res[1]))
-        assert z == r.reduce_points((xtest0 - 0.49 * r.res[0], ytest0 - 0.49 * r.res[1]))
-        assert z == r.reduce_points((xtest0 + 0.49 * r.res[0], ytest0 + 0.49 * r.res[1]))
-
-        # -- Tests 2: check arguments work as intended --
-
-        # 1/ Lat-lon argument check by getting the coordinates of our last test point
-        lat, lon = reproject_to_latlon(points=[xtest0, ytest0], in_crs=r.crs)
-        z_val_2 = r.reduce_points((lon, lat), input_latlon=True)
-        assert z_val == z_val_2
-
-        # 2/ Band argument
-        # Get the band indexes for the multi-band Raster
-        r_multi = gu.Raster(self.landsat_rgb_path)
-        itest, jtest = r_multi.xy2ij(xtest0, ytest0)
-        itest = int(itest[0])
-        jtest = int(jtest[0])
-        # Extract the values
-        z_band1 = r_multi.reduce_points((xtest0, ytest0), band=1)
-        z_band2 = r_multi.reduce_points((xtest0, ytest0), band=2)
-        z_band3 = r_multi.reduce_points((xtest0, ytest0), band=3)
-        # Compare to the Raster array slice
-        assert list(r_multi.data[:, itest, jtest]) == [z_band1, z_band2, z_band3]
-
-        # 3/ Masked argument
-        r_multi.data[:, itest, jtest] = np.ma.masked
-        z_not_ma = r_multi.reduce_points((xtest0, ytest0), band=1)
-        assert not np.ma.is_masked(z_not_ma)
-        z_ma = r_multi.reduce_points((xtest0, ytest0), band=1, masked=True)
-        assert np.ma.is_masked(z_ma)
-
-        # 4/ Window argument
-        val_window, z_window = r_multi.reduce_points(
-            (xtest0, ytest0), band=1, window=3, masked=True, return_window=True
-        )
-        assert (
-            val_window
-            == np.ma.mean(r_multi.data[0, itest - 1 : itest + 2, jtest - 1 : jtest + 2])
-            == np.ma.mean(z_window)
-        )
-        assert np.array_equal(z_window, r_multi.data[0, itest - 1 : itest + 2, jtest - 1 : jtest + 2])
-
-        # 5/ Reducer function argument
-        val_window2 = r_multi.reduce_points(
-            (xtest0, ytest0), band=1, window=3, masked=True, reducer_function=np.ma.median
-        )
-        assert val_window2 == np.ma.median(r_multi.data[0, itest - 1 : itest + 2, jtest - 1 : jtest + 2])
-
-        # -- Tests 3: check that errors are raised when supposed for non-boolean arguments --
-
-        # Verify that passing a window that is not a whole number fails
-        with pytest.raises(ValueError, match=re.escape("Window must be a whole number.")):
-            r.reduce_points((xtest0, ytest0), window=3.5)  # type: ignore
-        # Same for an odd number
-        with pytest.raises(ValueError, match=re.escape("Window must be an odd number.")):
-            r.reduce_points((xtest0, ytest0), window=4)
-        # But a window that is a whole number as a float works
-        r.reduce_points((xtest0, ytest0), window=3.0)  # type: ignore
-
-        # -- Tests 4: check that passing an array-like object works
-
-        # For simple coordinates
-        x_coords = [xtest0, xtest0 + 100]
-        y_coords = [ytest0, ytest0 - 100]
-        vals = r_multi.reduce_points((x_coords, y_coords))
-        val0, win0 = r_multi.reduce_points((x_coords[0], y_coords[0]), return_window=True)
-        val1, win1 = r_multi.reduce_points((x_coords[1], y_coords[1]), return_window=True)
-
-        assert len(vals) == len(x_coords)
-        assert np.array_equal(vals[0], val0, equal_nan=True)
-        assert np.array_equal(vals[1], val1, equal_nan=True)
-
-        # With a return window argument
-        vals, windows = r_multi.reduce_points((x_coords, y_coords), return_window=True)
-        assert len(windows) == len(x_coords)
-        assert np.array_equal(windows[0], win0, equal_nan=True)
-        assert np.array_equal(windows[1], win1, equal_nan=True)
-
-        # -- Tests 5 -- Check image corners and latlon argument
-
-        # Lower right pixel
-        x, y = [r.bounds.right - r.res[0] / 2, r.bounds.bottom + r.res[1] / 2]
-        lat, lon = pt.reproject_to_latlon([x, y], r.crs)
-        assert r.reduce_points((x, y)) == r.reduce_points((lon, lat), input_latlon=True) == r.data[-1, -1]
-
-        # One pixel above
-        x, y = [r.bounds.right - r.res[0] / 2, r.bounds.bottom + 3 * r.res[1] / 2]
-        lat, lon = pt.reproject_to_latlon([x, y], r.crs)
-        assert r.reduce_points((x, y)) == r.reduce_points((lon, lat), input_latlon=True) == r.data[-2, -1]
-
-        # One pixel left
-        x, y = [r.bounds.right - 3 * r.res[0] / 2, r.bounds.bottom + r.res[1] / 2]
-        lat, lon = pt.reproject_to_latlon([x, y], r.crs)
-        assert r.reduce_points((x, y)) == r.reduce_points((lon, lat), input_latlon=True) == r.data[-1, -2]
-
     @pytest.mark.parametrize("example", [landsat_b4_path, aster_dem_path])  # type: ignore
     def test_set_nodata(self, example: str) -> None:
         """
@@ -2890,46 +2447,6 @@ def test_save(self, example: str) -> None:
         except (NotADirectoryError, PermissionError):
             pass
 
-    @pytest.mark.parametrize("example", [landsat_b4_path, aster_dem_path, landsat_rgb_path])  # type: ignore
-    def test_coords(self, example: str) -> None:
-
-        img = gu.Raster(self.landsat_b4_path)
-
-        # With lower left argument
-        xx0, yy0 = img.coords(grid=False, force_offset="ll")
-
-        assert len(xx0) == img.width
-        assert len(yy0) == img.height
-
-        assert xx0[0] == pytest.approx(img.bounds.left)
-        assert xx0[-1] == pytest.approx(img.bounds.right - img.res[0])
-        if img.res[1] > 0:
-            assert yy0[0] == pytest.approx(img.bounds.bottom)
-            assert yy0[-1] == pytest.approx(img.bounds.top - img.res[1])
-        else:
-            # Currently not covered by test image
-            assert yy0[0] == pytest.approx(img.bounds.top)
-            assert yy0[-1] == pytest.approx(img.bounds.bottom + img.res[1])
-
-        # With center argument
-        xx, yy = img.coords(grid=False, force_offset="center")
-        hx = img.res[0] / 2
-        hy = img.res[1] / 2
-        assert xx[0] == pytest.approx(img.bounds.left + hx)
-        assert xx[-1] == pytest.approx(img.bounds.right - hx)
-        if img.res[1] > 0:
-            assert yy[0] == pytest.approx(img.bounds.bottom + hy)
-            assert yy[-1] == pytest.approx(img.bounds.top - hy)
-        else:
-            # Currently not covered by test image
-            assert yy[-1] == pytest.approx(img.bounds.top + hy)
-            assert yy[0] == pytest.approx(img.bounds.bottom - hy)
-
-        # With grid argument (default argument, repeated here for clarity)
-        xxgrid, yygrid = img.coords(grid=True, force_offset="ll")
-        assert np.array_equal(xxgrid, np.repeat(xx0[np.newaxis, :], img.height, axis=0))
-        assert np.array_equal(yygrid, np.flipud(np.repeat(yy0[:, np.newaxis], img.width, axis=1)))
-
     @pytest.mark.parametrize("example", [landsat_b4_path, aster_dem_path, landsat_rgb_path])  # type: ignore
     def test_from_array(self, example: str) -> None:
 
diff --git a/tests/test_vector.py b/tests/test_vector.py
index cabae434..08d848d3 100644
--- a/tests/test_vector.py
+++ b/tests/test_vector.py
@@ -166,7 +166,7 @@ def test_reproject(self) -> None:
         with pytest.raises(ValueError, match=re.escape("Reference raster or vector path does not exist.")):
             v0.reproject(ref="tmp.lol")
         # If it exists but cannot be opened by rasterio or fiona
-        with pytest.raises(ValueError, match=re.escape("Could not open raster or vector with rasterio or fiona.")):
+        with pytest.raises(ValueError, match=re.escape("Could not open raster or vector with rasterio or pyogrio.")):
             v0.reproject(ref="geoutils/examples.py")
         # If input of wrong type
         with pytest.raises(TypeError, match=re.escape("Type of ref must be string path to file, Raster or Vector.")):