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,56 @@ 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 [400, 1000, 2000, 3000, >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
+
+filename_rast = gu.examples.get_path("everest_landsat_b4")
+filename_vect = gu.examples.get_path("everest_rgi_outlines")
+rast = gu.Raster(filename_rast)
+vect = gu.Vector(filename_vect)
+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,97 @@ 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],
+        statistics: list[str],
+    ) -> pd.DataFrame:
+        """ """
+
+        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,

tests/test_raster/test_base.py

@@ -7,13 +7,15 @@
 import geopandas as gpd
 import numpy as np
 import pandas as pd
+import pandas.testing as pdt
 import pytest
 import rasterio as rio
 import xarray as xr
 from packaging.version import Version
 from pandas.testing import assert_frame_equal
 from pyproj import CRS
 
+from geoutils import examples  # flake8 error
 from geoutils import Raster, Vector, open_raster
 from geoutils.raster import MultiprocConfig
 from geoutils.raster.base import RasterBase
@@ -54,6 +56,11 @@ def assert_output_equal(output1: Any, output2: Any, use_allclose: bool = False,
         df1 = pd.DataFrame(index=[0], data=output1)
         df2 = pd.DataFrame(index=[0], data=output2)
         assert_frame_equal(df1, df2, check_dtype=False)
+
+    # For pandas DataFrame
+    elif isinstance(output1, pd.DataFrame) and isinstance(output2, pd.DataFrame):
+        pdt.assert_frame_equal(output1, output2)
+
     # For any other object type
     else:
         assert output1 == output2
@@ -219,6 +226,14 @@ def test_properties__equality_and_loading(self, path_index: int, prop: str, lazy
         ("subsample", {"subsample": 1000, "random_state": 42}),
         ("filter", {"method": "median", "size": 7}),
         ("get_stats", {}),
+        (
+            "grouped_stats",
+            {
+                "groupby_arrays": {"self"},
+                "bins": [[0, 5, 20]],
+                "statistics": ["mean"],
+            },
+        ),
         # 2.2. In-place methods
         ("load", {}),
     ]
@@ -267,6 +282,8 @@ def test_methods__equality_and_loading(
             args.update({"other": ds.copy(deep=False)})
         elif method == "copy" and "new_array" in args:
             args.update({"new_array": np.ones(ds.shape)})
+        elif method == "grouped_stats":
+            args.update({"groupby_arrays": {"raster_data": raster.data}})
 
         # Apply method for each class
         output_raster = getattr(raster, method)(**args)
@@ -504,3 +521,151 @@ def test_methods__match_raster(self) -> None:
 
         # TODO: Finalize after consistent input check function #850
         assert True
+
+
+class TestGroupedStatsRaster:
+    """Test the grouped_stats method of Raster"""
+
+    data = np.array([[1, 2, 3, 4, 5], [6, 7, 8, 9, 10], [11, 12, 13, 14, 15], [16, 17, 18, 19, 20]])
+    transform = 30.0, 0.0, 478000.0, 0.0, -30.0, 3108140.0
+    test_raster = Raster.from_array(data, transform, 32645)
+
+    def test_wrong_inputs(self) -> None:
+        """Test that wrong inputs raise the expected errors."""
+
+        raster = Raster(examples.get_path_test("exploradores_aster_dem"))
+
+        with pytest.raises(
+            ValueError, match="Bins should be provided as a list of arrays or lists, one per groupby variable."
+        ):
+            raster.grouped_stats(groupby_arrays={"slope": np.array([1, 2, 3])}, bins="not_a_list", statistics=["mean"])
+
+        with pytest.raises(
+            ValueError,
+            match="Groupby arrays should be provided as a dictionary of RasterType, one per groupby variable.",
+        ):
+            raster.grouped_stats(groupby_arrays=["slope", "aspect"], bins=[np.array([0, 1, 2])], statistics=["mean"])
+
+        with pytest.raises(ValueError, match="Statistics should be provided as a list of strings."):
+            raster.grouped_stats(
+                groupby_arrays={"slope": np.array([1, 2, 3])},
+                bins=[np.array([0, 1, 2])],
+                statistics="mean",
+            )
+
+        with pytest.raises(ValueError, match="One bins array must be provided per input array."):
+            raster.grouped_stats(
+                groupby_arrays={"slope": np.array([1, 2, 3]), "aspect": np.array([4, 5, 6])},
+                bins=[np.array([0, 1, 2])],
+                statistics=["mean"],
+            )
+
+    def test_interval_bins(self) -> None:
+        """Test that the method works with interval bins."""
+
+        bins = [[0, 10, 20]]
+        arrays = {"test_interval": self.test_raster.data}
+
+        # Ground truth
+        expected = pd.DataFrame(
+            {
+                "mean": [5.5, 15.5],
+                "std": [3.02765, 3.02765],
+            },
+            index=pd.CategoricalIndex(
+                pd.IntervalIndex.from_tuples([(0, 10), (10, 20)], closed="right"),
+                ordered=True,
+                name="bin_test_interval",
+            ),
+        )
+
+        grouped_stats_test = self.test_raster.grouped_stats(
+            groupby_arrays=arrays, bins=bins, statistics=["mean", "std"]
+        )
+
+        pdt.assert_frame_equal(expected, grouped_stats_test)
+
+    def test_boolean_mask(self) -> None:
+        """Test grouped stats with a boolean mask as bins."""
+
+        arrays = {"test_bool_mask": self.test_raster.data}
+
+        # Create mask for odd values
+        bins = [self.test_raster.data % 2 == 1]
+
+        # Ground truth
+        expected = pd.DataFrame(
+            {
+                "mean": [11.0, 10.0],
+                "std": [6.055301, 6.055301],
+            },
+            index=pd.Index([0, 1], name="bin_test_bool_mask"),
+        )
+
+        grouped_stats_test = self.test_raster.grouped_stats(
+            groupby_arrays=arrays, bins=bins, statistics=["mean", "std"]
+        )
+
+        pdt.assert_frame_equal(expected, grouped_stats_test)
+
+    def test_multimodal_mask(self) -> None:
+        """Test grouped stats with a boolean mask as bins."""
+
+        arrays = {"test_multimodal_mask": self.test_raster.data}
+
+        # Create boolean mask for odd and prime values
+        bins = [np.array([[1, 0, 2, 0, 2], [0, 2, 0, 1, 0], [2, 0, 2, 0, 1], [0, 2, 0, 2, 0]])]
+
+        # Ground truth
+        expected = pd.DataFrame(
+            {
+                "mean": [11.0, 8.333333, 10.714286],
+                "std": [6.055301, 7.023769, 6.047432],
+            },
+            index=pd.Index([0, 1, 2], name="bin_test_multimodal_mask"),
+        )
+
+        grouped_stats_test = self.test_raster.grouped_stats(
+            groupby_arrays=arrays, bins=bins, statistics=["mean", "std"]
+        )
+
+        pdt.assert_frame_equal(expected, grouped_stats_test)
+
+    def test_multi_bins(self) -> None:
+        """"""
+
+        arrays = {"test_multimodal_mask": self.test_raster.data, "test_elev": self.test_raster.data}
+
+        # Create boolean mask for odd and prime values
+        bins = [
+            np.array([[1, 0, 2, 0, 2], [0, 2, 0, 1, 0], [2, 0, 2, 0, 1], [0, 2, 0, 2, 0]]),
+            self.test_raster.data > 10,
+        ]
+
+        # Ground truth
+        index = pd.MultiIndex.from_tuples(
+            [
+                (0, 0),
+                (0, 1),
+                (1, 0),
+                (1, 1),
+                (2, 0),
+                (2, 1),
+            ],
+            names=["bin_test_multimodal_mask", "bin_test_elev"],
+        )
+
+        # Ground truth
+        expected = pd.DataFrame(
+            {
+                "mean": [6.0, 16.0, 5.0, 15.0, 5.0, 15.0],
+                "std": [3.162278, 3.162278, 5.656854, np.nan, 2.0, 3.651484],
+            },
+            index=index,
+        )
+
+        grouped_stats_test = self.test_raster.grouped_stats(
+            groupby_arrays=arrays, bins=bins, statistics=["mean", "std"]
+        )
+
+        pdt.assert_frame_equal(expected, grouped_stats_test)