Handling a ROI when loading raster data (#642)

doc/Makefile

@@ -15,9 +15,9 @@ help:
 
 clean:
 	echo "Removing build files..."
-	if [[ -d "$(BUILDDIR)" ]]; then rm -r "$(BUILDDIR)"; fi
-	if [[ -d "$(SOURCEDIR)/auto_examples" ]]; then rm -r "$(SOURCEDIR)/auto_examples"; fi
-	if [[ -d "$(SOURCEDIR)/gen_modules" ]]; then rm -r "$(SOURCEDIR)/gen_modules"; fi
+	if [ -d "$(BUILDDIR)" ]; then rm -r "$(BUILDDIR)"; fi
+	if [ -d "$(SOURCEDIR)/auto_examples" ]; then rm -r "$(SOURCEDIR)/auto_examples"; fi
+	if [ -d "$(SOURCEDIR)/gen_modules" ]; then rm -r "$(SOURCEDIR)/gen_modules"; fi
 
 # Catch-all target: route all unknown targets to Sphinx using the new
 # "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).

doc/source/api.md

@@ -90,6 +90,7 @@ documentation.
     :toctree: gen_modules/
 
     Raster.crop
+    Raster.icrop
     Raster.reproject
     Raster.polygonize
     Raster.proximity

doc/source/index.md

@@ -135,7 +135,6 @@ config
 :maxdepth: 2
 
 background
-authors
 license
 ```
 

doc/source/raster_class.md

@@ -244,6 +244,8 @@ Resampling methods are listed in **[the dedicated section of Rasterio's API](htt
 ## Crop
 
 Cropping a {class}`~geoutils.Raster` is done through the {func}`~geoutils.Raster.crop` function, which enforces new {attr}`~geoutils.Raster.bounds`.
+Additionally, you can use the {func}`~geoutils.Raster.icrop` method to crop the raster using pixel coordinates instead of geographic bounds.
+Both cropping methods can be used before loading the raster's data into memory. This optimization can prevent loading unnecessary parts of the data, which is particularly useful when working with large rasters.
 
 ```{important}
 As with all geospatial handling methods, the {func}`~geoutils.Raster.crop` function can be passed only a {class}`~geoutils.Raster` or {class}`~geoutils.Vector`
@@ -254,14 +256,24 @@ See {ref}`core-match-ref` for more details.
 
 The {func}`~geoutils.Raster.crop` function can also be passed a {class}`list` or {class}`tuple` of bounds (`xmin`, `ymin`, `xmax`, `ymax`). By default,
 {func}`~geoutils.Raster.crop` returns a new Raster.
+The {func}`~geoutils.Raster.icrop` function accepts only a bounding box in pixel coordinates (xmin, ymin, xmax, ymax) and crop the raster accordingly.
+By default, {func}`~geoutils.Raster.crop` and {func}`~geoutils.Raster.icrop` return a new Raster unless the inplace parameter is set to True, in which case the cropping operation is performed directly on the original raster object.
 For more details, see the {ref}`specific section and function descriptions in the API<api-geo-handle>`.
 
+### Example for `crop`
 ```{code-cell} ipython3
 # Crop raster to smaller bounds
-rast_crop = rast.crop(crop_geom=(0.3, 0.3, 1, 1))
+rast_crop = rast.crop(bbox=(0.3, 0.3, 1, 1))
 print(rast_crop.bounds)
 ```
 
+### Example for `icrop`
+```{code-cell} ipython3
+# Crop raster using pixel coordinates
+rast_icrop = rast.icrop(bbox=(2, 2, 6, 6))
+print(rast_icrop.bounds)
+```
+
 ## Polygonize
 
 Polygonizing a {class}`~geoutils.Raster` is done through the {func}`~geoutils.Raster.polygonize` function, which converts target pixels into a multi-polygon

geoutils/raster/geotransformations.py

@@ -443,7 +443,7 @@ def _reproject(
 
 def _crop(
     source_raster: gu.Raster,
-    crop_geom: gu.Raster | gu.Vector | list[float] | tuple[float, ...],
+    bbox: gu.Raster | gu.Vector | list[float] | tuple[float, ...],
     mode: Literal["match_pixel"] | Literal["match_extent"] = "match_pixel",
 ) -> tuple[MArrayNum, affine.Affine]:
     """Crop raster. See details in Raster.crop()."""
@@ -453,14 +453,14 @@ def _crop(
         "match_pixel",
     ], "mode must be one of 'match_pixel', 'match_extent'"
 
-    if isinstance(crop_geom, (gu.Raster, gu.Vector)):
+    if isinstance(bbox, (gu.Raster, gu.Vector)):
         # For another Vector or Raster, we reproject the bounding box in the same CRS as self
-        xmin, ymin, xmax, ymax = crop_geom.get_bounds_projected(out_crs=source_raster.crs)
-        if isinstance(crop_geom, gu.Raster):
+        xmin, ymin, xmax, ymax = bbox.get_bounds_projected(out_crs=source_raster.crs)
+        if isinstance(bbox, gu.Raster):
             # Raise a warning if the reference is a raster that has a different pixel interpretation
-            _cast_pixel_interpretation(source_raster.area_or_point, crop_geom.area_or_point)
-    elif isinstance(crop_geom, (list, tuple)):
-        xmin, ymin, xmax, ymax = crop_geom
+            _cast_pixel_interpretation(source_raster.area_or_point, bbox.area_or_point)
+    elif isinstance(bbox, (list, tuple)):
+        xmin, ymin, xmax, ymax = bbox
     else:
         raise ValueError("cropGeom must be a Raster, Vector, or list of coordinates.")
 
@@ -479,11 +479,8 @@ def _crop(
         if source_raster.is_loaded:
             # In case data is loaded on disk, can extract directly from np array
             (rowmin, rowmax), (colmin, colmax) = final_window.toranges()
+            crop_img = source_raster.data[..., rowmin:rowmax, colmin:colmax]
 
-            if source_raster.count == 1:
-                crop_img = source_raster.data[rowmin:rowmax, colmin:colmax]
-            else:
-                crop_img = source_raster.data[:, rowmin:rowmax, colmin:colmax]
         else:
 
             assert source_raster._disk_shape is not None  # This should not be the case, sanity check to make mypy happy

geoutils/raster/raster.py

@@ -1,4 +1,5 @@
 # Copyright (c) 2025 GeoUtils developers
+# Copyright (c) 2025 Centre National d'Etudes Spatiales (CNES)
 #
 # This file is part of the GeoUtils project:
 # https://github.com/glaciohack/geoutils
@@ -218,7 +219,7 @@ def _load_rio(
     * resampling : to set the resampling algorithm
     """
     # If out_shape is passed, no need to account for transform and shape
-    if kwargs["out_shape"] is not None:
+    if kwargs.get("out_shape") is not None:
         window = None
         # If multi-band raster, the out_shape needs to contain the count
         if out_count is not None and out_count > 1:
@@ -233,8 +234,6 @@ def _load_rio(
             window = rio.windows.Window(col_off, row_off, *shape[::-1])
         elif sum(param is None for param in [shape, transform]) == 1:
             raise ValueError("If 'shape' or 'transform' is provided, BOTH must be given.")
-        else:
-            window = None
 
     if indexes is None:
         if only_mask:
@@ -441,7 +440,7 @@ def __init__(
         if isinstance(filename_or_dataset, Raster):
             for key in filename_or_dataset.__dict__:
                 setattr(self, key, filename_or_dataset.__dict__[key])
-            return
+
         # Image is a file on disk.
         elif isinstance(filename_or_dataset, (str, pathlib.Path, rio.io.DatasetReader, rio.io.MemoryFile)):
             # ExitStack is used instead of "with rio.open(filename_or_dataset) as ds:".
@@ -492,6 +491,7 @@ def __init__(
             # Downsampled image size
             if not isinstance(downsample, (int, float)):
                 raise TypeError("downsample must be of type int or float.")
+
             if downsample == 1:
                 out_shape = (self.height, self.width)
             else:
@@ -2387,7 +2387,7 @@ def __array_function__(
     @overload
     def crop(
         self: RasterType,
-        crop_geom: RasterType | Vector | list[float] | tuple[float, ...],
+        bbox: RasterType | Vector | list[float] | tuple[float, ...],
         mode: Literal["match_pixel"] | Literal["match_extent"] = "match_pixel",
         *,
         inplace: Literal[False] = False,
@@ -2396,7 +2396,7 @@ def crop(
     @overload
     def crop(
         self: RasterType,
-        crop_geom: RasterType | Vector | list[float] | tuple[float, ...],
+        bbox: RasterType | Vector | list[float] | tuple[float, ...],
         mode: Literal["match_pixel"] | Literal["match_extent"] = "match_pixel",
         *,
         inplace: Literal[True],
@@ -2405,15 +2405,15 @@ def crop(
     @overload
     def crop(
         self: RasterType,
-        crop_geom: RasterType | Vector | list[float] | tuple[float, ...],
+        bbox: RasterType | Vector | list[float] | tuple[float, ...],
         mode: Literal["match_pixel"] | Literal["match_extent"] = "match_pixel",
         *,
         inplace: bool = False,
     ) -> RasterType | None: ...
 
     def crop(
         self: RasterType,
-        crop_geom: RasterType | Vector | list[float] | tuple[float, ...],
+        bbox: RasterType | Vector | list[float] | tuple[float, ...],
         mode: Literal["match_pixel"] | Literal["match_extent"] = "match_pixel",
         *,
         inplace: bool = False,
@@ -2425,8 +2425,8 @@ def crop(
 
         Reprojection is done on the fly if georeferenced objects have different projections.
 
-        :param crop_geom: Geometry to crop raster to. Can use either a raster or vector as match-reference, or a list of
-            coordinates. If ``crop_geom`` is a raster or vector, will crop to the bounds. If ``crop_geom`` is a
+        :param bbox: Geometry to crop raster to. Can use either a raster or vector as match-reference, or a list of
+            coordinates. If ``bbox`` is a raster or vector, will crop to the bounds. If ``bbox`` is a
             list of coordinates, the order is assumed to be [xmin, ymin, xmax, ymax].
         :param mode: Whether to match within pixels or exact extent. ``'match_pixel'`` will preserve the original pixel
             resolution, cropping to the extent that most closely aligns with the current coordinates. ``'match_extent'``
@@ -2436,7 +2436,7 @@ def crop(
         :returns: A new raster (or None if inplace).
         """
 
-        crop_img, tfm = _crop(source_raster=self, crop_geom=crop_geom, mode=mode)
+        crop_img, tfm = _crop(source_raster=self, bbox=bbox, mode=mode)
 
         if inplace:
             self._data = crop_img
@@ -2446,6 +2446,34 @@ def crop(
             newraster = self.from_array(crop_img, tfm, self.crs, self.nodata, self.area_or_point)
             return newraster
 
+    def icrop(
+        self: RasterType,
+        bbox: list[int] | tuple[int, ...],
+        *,
+        inplace: bool = False,
+    ) -> RasterType | None:
+        """
+        Crop raster based on pixel indices (bbox), converting them into georeferenced coordinates.
+
+        :param bbox: Pixel-based bounding box as (xmin, ymin, xmax, ymax) in pixel coordinates.
+        :param inplace: If True, modify the raster in place. Otherwise, return a new cropped raster.
+
+        :returns: Cropped raster or None (if inplace=True).
+        """
+        # Ensure bbox contains four elements (xmin, ymin, xmax, ymax)
+        if len(bbox) != 4:
+            raise ValueError("bbox must be a list or tuple of four integers: (xmin, ymin, xmax, ymax).")
+
+        # Convert pixel coordinates to georeferenced coordinates using the transform
+        bottom_left = self.transform * (bbox[0], self.height - bbox[1])
+        top_right = self.transform * (bbox[2], self.height - bbox[3])
+
+        # Create new georeferenced crop geometry (xmin, ymin, xmax, ymax)
+        new_bbox = (bottom_left[0], bottom_left[1], top_right[0], top_right[1])
+
+        # Call the existing crop() method with the new georeferenced crop geometry
+        return self.crop(bbox=new_bbox, inplace=inplace)
+
     @overload
     def reproject(
         self: RasterType,
@@ -3952,7 +3980,7 @@ def reproject(
     @overload
     def crop(
         self: Mask,
-        crop_geom: Mask | Vector | list[float] | tuple[float, ...],
+        bbox: Mask | Vector | list[float] | tuple[float, ...],
         mode: Literal["match_pixel"] | Literal["match_extent"] = "match_pixel",
         *,
         inplace: Literal[False] = False,
@@ -3961,7 +3989,7 @@ def crop(
     @overload
     def crop(
         self: Mask,
-        crop_geom: Mask | Vector | list[float] | tuple[float, ...],
+        bbox: Mask | Vector | list[float] | tuple[float, ...],
         mode: Literal["match_pixel"] | Literal["match_extent"] = "match_pixel",
         *,
         inplace: Literal[True],
@@ -3970,15 +3998,15 @@ def crop(
     @overload
     def crop(
         self: Mask,
-        crop_geom: Mask | Vector | list[float] | tuple[float, ...],
+        bbox: Mask | Vector | list[float] | tuple[float, ...],
         mode: Literal["match_pixel"] | Literal["match_extent"] = "match_pixel",
         *,
         inplace: bool = False,
     ) -> Mask | None: ...
 
     def crop(
         self: Mask,
-        crop_geom: Mask | Vector | list[float] | tuple[float, ...],
+        bbox: Mask | Vector | list[float] | tuple[float, ...],
         mode: Literal["match_pixel"] | Literal["match_extent"] = "match_pixel",
         *,
         inplace: bool = False,
@@ -3988,16 +4016,16 @@ def crop(
             raise ValueError(NotImplementedError)
             # self._data = self.data.astype("float32")
             # if inplace:
-            #     super().crop(crop_geom=crop_geom, mode=mode, inplace=inplace)
+            #     super().crop(bbox=bbox, mode=mode, inplace=inplace)
             #     self._data = self.data.astype(bool)
             #     return None
             # else:
-            #     output = super().crop(crop_geom=crop_geom, mode=mode, inplace=inplace)
+            #     output = super().crop(bbox=bbox, mode=mode, inplace=inplace)
             #     output._data = output.data.astype(bool)
             #     return output
         # Otherwise, run a classic crop
         else:
-            return super().crop(crop_geom=crop_geom, mode=mode, inplace=inplace)
+            return super().crop(bbox=bbox, mode=mode, inplace=inplace)
 
     def polygonize(
         self,