From de39c35b8a6b7f80aea4383df1176045c398494c Mon Sep 17 00:00:00 2001
From: Romain Hugonnet <romain.hugonnet@gmail.com>
Date: Thu, 22 Aug 2024 17:26:59 -0800
Subject: [PATCH] Add bounds to georeferencing, accept None input for
 interpolate with return_interpolator

---
 geoutils/raster/georeferencing.py |  6 ++++++
 geoutils/raster/interpolate.py    | 20 +++++++++++---------
 geoutils/raster/raster.py         |  4 ++--
 3 files changed, 19 insertions(+), 11 deletions(-)

diff --git a/geoutils/raster/georeferencing.py b/geoutils/raster/georeferencing.py
index c5ce5e33..d64c637d 100644
--- a/geoutils/raster/georeferencing.py
+++ b/geoutils/raster/georeferencing.py
@@ -164,3 +164,9 @@ 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]):
+    """See description of Raster.bounds."""
+
+    return rio.coords.BoundingBox(*rio.transform.array_bounds(height=shape[0], width=shape[1], transform=transform))
\ No newline at end of file
diff --git a/geoutils/raster/interpolate.py b/geoutils/raster/interpolate.py
index 40d9ae70..12e78331 100644
--- a/geoutils/raster/interpolate.py
+++ b/geoutils/raster/interpolate.py
@@ -237,7 +237,7 @@ def _interp_points(
     array: NDArrayNum,
     transform: rio.transform.Affine,
     area_or_point: Literal["Area", "Point"] | None,
-    points: tuple[Number, Number] | tuple[NDArrayNum, NDArrayNum],
+    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,
@@ -250,16 +250,18 @@ def _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
 
-    # Get coordinates
-    x, y = points
-
-    i, j = _xy2ij(x, y, transform=transform, area_or_point=area_or_point, shift_area_or_point=shift_area_or_point)
+    # TODO: Add check about None for "points" depending on "return_interpolator"
+    if not return_interpolator:
+        x, y = points
+        i, j = _xy2ij(x, y, transform=transform, area_or_point=area_or_point,
+                      shift_area_or_point=shift_area_or_point)
 
-    shape: tuple[int, int] = array.shape[0:2]  # type: ignore
-    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)
+        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"
diff --git a/geoutils/raster/raster.py b/geoutils/raster/raster.py
index 7032ab99..04528f90 100644
--- a/geoutils/raster/raster.py
+++ b/geoutils/raster/raster.py
@@ -50,7 +50,7 @@
     reproject_from_latlon,
 )
 from geoutils.raster.array import get_mask_from_array
-from geoutils.raster.georeferencing import _coords, _ij2xy, _outside_image, _res, _xy2ij
+from geoutils.raster.georeferencing import _coords, _ij2xy, _outside_image, _res, _xy2ij, _bounds
 from geoutils.raster.interpolate import _interp_points
 from geoutils.raster.sampling import subsample_array
 from geoutils.vector import Vector
@@ -793,7 +793,7 @@ def res(self) -> tuple[float | int, float | int]:
     @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: