Add grouped_stats for categorical or continuous binning

doc/source/stats.md

@@ -61,7 +61,7 @@ import numpy as np
 
 # Instantiate a raster from a filename on disk
 filename_rast = gu.examples.get_path("exploradores_aster_dem")
-rast = gu.Raster(filename_rast)
+rast = gu.Raster(filename_rast, force_nodata=-9999)
 rast
 ```
 
@@ -93,6 +93,54 @@ inlier_mask = rast > 1500
 rast.get_stats(inlier_mask=inlier_mask)
 ```
 
+## Grouped statistics
+
+GeoUtils provides support for grouped statistics, allowing statistics to be computed independently over subsets of data
+defined by one or more grouping bins. This is particularly useful when analyzing how statistical properties vary
+across classes, bins, or segmentation derived from the data itself.
+
+
+### Example with altitude intervals
+In this example, we will create different altitude classes from a chosen interval [1000, 2000, 3000].
+Once these bins are created, we reapply them and compute the mean, minimum, and maximum values of the same raster
+for each sub-interval. It is also possible to use a reference other than the raster itself for the group_by.
+
+A dictionary containing the masks that have been created during the computation will also be returned by the function.
+Using GeoUtils functions makes it very easy to visualise them.
+
+```{code-cell} ipython3
+group_by = {"raster": rast.data}
+bins = [[1000, 2000, 3000]]
+statistics = ["mean", "min", "max"]
+
+df = rast.grouped_stats(group_by, bins, statistics)
+df
+```
+
+### Example with vector outlines masks
+
+In this example, we will create a mask such as altitude is more than 2000 meters.
+Once these masks are created as a raster, we reapply them and compute the mean, minimum, and maximum values of the same raster
+for masks = True. It is also possible to use a reference other than the raster itself for the group_by.
+
+```{code-cell} ipython3
+import matplotlib.pyplot as plt
+
+rast = gu.Raster(gu.examples.get_path("everest_landsat_b4"))
+vect = gu.Vector(gu.examples.get_path("everest_rgi_outlines"))
+vect_rasterized = vect.create_mask(rast)
+
+vect_rasterized.plot()
+plt.show()
+
+group_by = {"elevation": rast.data}
+bins = [vect_rasterized.data]
+statistics = ["mean", "min", "max"]
+
+df = rast.grouped_stats(group_by, bins, statistics)
+df
+```
+
 ## Subsampling
 
 The {func}`~geoutils.Raster.subsample` method allows to efficiently extract a valid random subsample from a raster or a point cloud. It can conveniently

geoutils/raster/base.py

@@ -20,6 +20,7 @@
 
 from __future__ import annotations
 
+import logging
 import math
 import pathlib
 import struct
@@ -38,6 +39,7 @@
 )
 
 import numpy as np
+import pandas as pd
 import rasterio as rio
 import xarray as xr
 from affine import Affine
@@ -795,6 +797,131 @@ def get_stats(
             else:
                 warnings.warn("Statistic name " + str(stats_name) + " is a not recognized string", category=UserWarning)
 
+    @overload
+    def grouped_stats(
+        self,
+        groupby_arrays: dict[str, RasterType | NDArrayNum],
+        bins: list[Any],
+        statistics: list[str],
+    ) -> pd.DataFrame: ...
+
+    @overload
+    def grouped_stats(
+        self,
+        groupby_arrays: dict[str, RasterType | NDArrayNum],
+        bins: list[Any],
+        statistics: list[str],
+    ) -> pd.DataFrame: ...
+
+    @profiler.profile("geoutils.raster.stats.grouped_stats", memprof=True)  # type: ignore
+    def grouped_stats(
+        self,
+        groupby_arrays: dict[str, NDArrayNum],
+        bins: list[NDArrayNum | list[float] | list[int] | list[bool]],
+        statistics: list[str],
+    ) -> pd.DataFrame:
+        """
+        Compute statistics of the raster values grouped by one or several variables.
+
+        The function groups the values stored in `self.data` according to one or
+        several arrays provided in `groupby_arrays` and computes the requested
+        statistics within bins defined in `bins`. Each groupby variable must have
+        a corresponding bin specification.
+
+        Binning can be defined either by:
+        - a list of bin edges (used with `pandas.cut`), or
+        - a numpy array providing a precomputed bin or mask.
+
+        Notes :
+        - Boolean arrays are interpreted as binary bins (False=0, True=1).
+        - Masked arrays in `self.data` are converted to `NaN` before aggregation.
+        - If a single bin edge is provided (e.g. `[x]`), it is automatically expanded to `[-inf, x, inf]` to create
+        two bins.
+        - All arrays are flattened before grouping and aggregation.
+
+        :param groupby_arrays: Dictionary of arrays used to group the data. Keys correspond to the variable names and
+        values are arrays with the same shape as `self.data`. Each array defines a grouping variable.
+
+        :param bins: List defining the binning strategy for each groupby variable. The list length must match the number
+         of variables in `groupby_arrays`.
+
+            Each element can be:
+            - a list of numeric bin edges (used with `pandas.cut`),
+            - a numpy array providing precomputed bin values,
+            - a boolean array interpreted as binary bins.
+
+        :param statistics: List of aggregation statistics to compute. These must be valid pandas aggregation functions
+        such as `'mean'`, `'median'`,`'min'`, `'max'`, `'std'`, `'count'`, etc.
+        :returns: A DataFrame indexed by the bin combinations of the groupby variables, containing the requested
+        statistics computed on the flattened values of `self.data`
+        """
+
+        if not isinstance(bins, list):
+            raise ValueError("Bins should be provided as a list of arrays or lists, one per groupby variable.")
+
+        if not isinstance(groupby_arrays, dict):
+            raise ValueError(
+                "Groupby arrays should be provided as a dictionary of RasterType, one per groupby variable."
+            )
+        if not isinstance(statistics, list):
+            raise ValueError("Statistics should be provided as a list of strings. Ex: ['mean', 'median']")
+
+        if len(groupby_arrays) != len(bins):
+            raise ValueError("One bins array must be provided per input array.")
+
+        # Flattened data to aggregate
+        if isinstance(self.data, np.ma.MaskedArray):
+            values = self.data.filled(np.nan).ravel()
+        else:
+            values = self.data.ravel()
+
+        # Working on bins when list
+        for idx, one_bin in enumerate(bins):
+            if isinstance(one_bin, list) and len(one_bin) == 1 and isinstance(one_bin[0], (int, float)):
+                logging.info(
+                    f"Only one bin edge provided for groupby_arrays variable {list(groupby_arrays.keys())[idx]}"
+                    " creating two bins with -inf and inf as limits."
+                )
+                bins[idx] = [-np.inf, one_bin[0], np.inf]
+
+        # Transform groupby arrays into a DataFrame for grouping and aggregation
+        df = pd.DataFrame(
+            {**{key_groupby: values_groupby.ravel() for key_groupby, values_groupby in groupby_arrays.items()}}
+        )
+        # Add the values to be aggregated to the DataFrame
+        df["values"] = values
+
+        # Apply binning
+        group_keys = []
+        for key, bin_edges in zip(groupby_arrays.keys(), bins):
+            bin_col = f"bin_{key}"
+            # Test boolean mask case
+            if isinstance(bin_edges, (np.ndarray, np.ma.MaskedArray)):
+                if bin_edges.dtype == bool:
+                    logging.info(
+                        f"Boolean mask provided for groupby variable {key}. "
+                        "It will be treated as binary binning with False(0) and True(1). "
+                        "If this is not intended, provide a list of bin edges instead."
+                    )
+                    df[bin_col] = bin_edges.ravel().astype(int)
+                else:
+                    df[bin_col] = bin_edges.ravel()
+
+            elif isinstance(bin_edges, list):
+                df[bin_col] = pd.cut(df[key], bins=bin_edges, include_lowest=False)
+
+            else:
+                raise ValueError(
+                    f"Bin edges for groupby variable {key} must be either "
+                    "a list of numeric edges or a numpy array (boolean or numeric)."
+                )
+
+            group_keys.append(bin_col)
+        # Perform aggregation
+        result = df.groupby(group_keys, observed=False)["values"].agg(statistics)
+
+        return result
+
     def _raster_equal_allclose(
         self,
         other: RasterType,
@@ -1117,7 +1244,7 @@ def crop(
         current coordinates. To match the extent of another dataset exactly, use reproject().
 
         :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
+            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 inplace: (DEPRECATED. Use rast = rast.crop() instead) Whether to crop in-place or not.
         :returns: A new cropped raster.
@@ -1222,14 +1349,14 @@ def reproject(
         The reprojected raster is written to disk under the path specified in the configuration
 
         :param ref: Reference raster to match resolution, bounds and CRS.
-        :param crs: Destination coordinate reference system as a string or EPSG. If ``ref`` not set,
+        :param crs: Destination coordinate reference system as a string or EPSG. If `ref` not set,
             defaults to this raster's CRS.
         :param res: Destination resolution (pixel size) in units of destination CRS. Single value or (xres, yres).
-            Do not use with ``grid_size``.
-        :param grid_size: Destination grid size as (x, y). Do not use with ``res``.
+            Do not use with `grid_size`.
+        :param grid_size: Destination grid size as (x, y). Do not use with `res`.
         :param bounds: Destination bounds as a Rasterio bounding box, or a dictionary containing left, bottom,
             right, top bounds in the destination CRS.
-        :param nodata: Destination nodata value. If set to ``None``, will use the same as source. If source does
+        :param nodata: Destination nodata value. If set to `None`, will use the same as source. If source does
             not exist, will use GDAL's default.
         :param dtype: Destination data type of array.
         :param resampling: A Rasterio resampling method, can be passed as a string.
@@ -1473,9 +1600,9 @@ def outside_image(self, xi: ArrayLike, yj: ArrayLike, index: bool = True) -> boo
 
         :param xi: Indices (or coordinates) of x direction to check.
         :param yj: Indices (or coordinates) of y direction to check.
-        :param index: Interpret ij as raster indices (default is ``True``). If False, assumes ij is coordinates.
+        :param index: Interpret ij as raster indices (default is `True`). If False, assumes ij is coordinates.
 
-        :returns is_outside: ``True`` if ij is outside the bounds.
+        :returns is_outside: `True` if ij is outside the bounds.
         """
 
         return _outside_bounds(