From 67414d2f98fa3e28a9a99f493daa0ddc377f4588 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alice=20de=20Bardonn=C3=A8che-Richard?= Date: Wed, 7 Jan 2026 15:33:25 +0100 Subject: [PATCH 1/4] old vision of stats --- doc/source/stats.md | 86 +++++++- geoutils/stats/grouped_stats.py | 171 ++++++++++++++++ tests/test_stats/test_grouped_stats.py | 267 +++++++++++++++++++++++++ 3 files changed, 523 insertions(+), 1 deletion(-) create mode 100644 geoutils/stats/grouped_stats.py create mode 100644 tests/test_stats/test_grouped_stats.py diff --git a/doc/source/stats.md b/doc/source/stats.md index f640a7881..de3ac32d0 100644 --- a/doc/source/stats.md +++ b/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,90 @@ 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 +from geoutils.stats import grouped_stats +import math +import matplotlib.pyplot as plt + +group_by = {"rast": rast} +bins = {"rast": [400, 1000, 2000, 3000, np.inf]} +to_aggregate = {"rast": rast} +statistics = ["mean", "min", "max"] + +df, masks = grouped_stats.grouped_stats(group_by, bins, to_aggregate, statistics) +df +``` + +```{code-cell} ipython3 +:tags: [hide-input] +:mystnb: +: code_prompt_show: "Show the code for plotting the figure" +: code_prompt_hide: "Hide the code for plotting the figure" + +groups = list(masks["groupby_rast"].keys()) +n = len(groups) + +ncols = 3 +nrows = math.ceil(n / ncols) + +fig, axes = plt.subplots(nrows, ncols, figsize=(5*ncols, 5*nrows)) +axes = axes.flatten() + +for ax, group in zip(axes, groups): + masks["groupby_rast"][group].plot(ax=ax) + ax.set_title(group) + +for ax in axes[n:]: + ax.axis("off") + +plt.tight_layout() +``` + +```{warning} +Bins can be presented in different ways. It is possible to integrate an interval of minimum 2 values, a mask or a +segmentation map in raster format. +``` + + +### Example with altitude 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 +from geoutils.stats import grouped_stats + +group_by = {"rast": rast} +elev_mask = rast > 2000 +bins = {"rast": elev_mask} + +to_aggregate = {"rast": rast} +statistics = ["mean", "min", "max"] + +df, _ = grouped_stats.grouped_stats(group_by, bins, to_aggregate, statistics) +df +``` + +```{code-cell} ipython3 +elev_mask.plot() +plt.show() +``` ## 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 diff --git a/geoutils/stats/grouped_stats.py b/geoutils/stats/grouped_stats.py new file mode 100644 index 000000000..1dcb54862 --- /dev/null +++ b/geoutils/stats/grouped_stats.py @@ -0,0 +1,171 @@ +# Copyright (c) 2025 Centre National d'Etudes Spatiales (CNES). +# Copyright (c) 2025 GeoUtils developers +# +# This file is part of the xDEM project: +# https://github.com/glaciohack/xdem +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +""" +This module allows the user to process grouped statistics easily +""" + +from __future__ import annotations + +import logging +from typing import Any + +import numpy as np +import pandas as pd + +from geoutils._typing import NDArrayNum +from geoutils.raster.raster import Raster, RasterType + + +def from_raster_to_flattened(dict_raster: dict[str, Any]) -> dict[str, Any]: + """ + Transform every raster in flattened data to be processed by pandas + :param: Dictionary containing binning name associated to be grouped by and a raster + """ + dict_arrays = dict_raster.copy() + + for key, raster in dict_raster.items(): + dict_arrays[key] = raster.data.ravel() + + return dict_arrays + + +def is_interval(bin_value: RasterType | list[int] | object) -> bool: + """ + Verifying the interval coherence if binning is a list + :param bin_value: Bin value + """ + if isinstance(bin_value, list): + bins_array = np.asarray(bin_value) + else: + return False + + if bins_array.ndim != 1: + raise ValueError("If bins is an interval, it must be a 1-dimensional array") + elif not np.issubdtype(bins_array.dtype, np.number): + raise TypeError("Bins must be a list of number") + elif bins_array.size < 2: + raise ValueError("Bins must be of size >= 2") + elif not np.all(np.diff(bins_array) > 0): + raise ValueError("Values must be strictly increasing.") + + return True + + +def grouped_stats( + groupby_vars: dict[str, RasterType], + bins: dict[str, RasterType | list[int] | object], + aggregated_vars: dict[str, NDArrayNum], + statistics: list[str], +) -> tuple[Any, dict[str, dict[str, NDArrayNum] | None]]: + """ + Get statistics grouped (=binned) by other variables, whether categorical or continuous. + :param groupby_vars: Dictionary containing Raster to group by. + :param bins: Bins to use. Can be a list of interval, binary mask or segmentation map. + :param aggregated_vars: Values to group, can be a dictionary . + :param statistics: List or dict of statistics to compute, e.g. ["mean", "std"]. + """ + + # Check that groupby_vars and bins as the same values + if set(groupby_vars) != set(bins): + raise ValueError("One bins/mask/segmentation entry per input array required.") + + # Load bins + new_bins = {} + for key, value in bins.items(): + if isinstance(value, str): + new_bins[key] = Raster(value) + logging.info(f"No need to save {value} again") + else: + new_bins[key] = value + + raster_base = next(iter(groupby_vars.values())) + crs = raster_base.crs + shape = raster_base.data.shape + transform = raster_base.transform + nodata = raster_base.nodata + + # Flatten arrays for pandas + groupby_vars_array = from_raster_to_flattened(groupby_vars) + aggregated_vars_arrays = from_raster_to_flattened(aggregated_vars) + + df = pd.DataFrame({**groupby_vars_array, **aggregated_vars_arrays}) + if nodata is not None: + df = df.replace(nodata, np.nan) + + groupby_keys = [] + returned_masks = {} + + for groupby_key in groupby_vars.keys(): + group_col = f"groupby_{groupby_key}" + bins_array = new_bins[groupby_key] + + if isinstance(bins[groupby_key], (str, Raster)): + returned_masks[group_col] = None + else: + returned_masks[group_col] = {} + + # ---------- Interval bins ---------- + if is_interval(bins_array): + cut = pd.cut(df[groupby_key], bins=bins_array) + df[group_col] = cut + + if returned_masks[group_col] is not None: + for interval in cut.cat.categories: + mask_flat = (cut == interval).to_numpy() + mask_to_raster = Raster.from_array(mask_flat.reshape(shape), transform, crs) + + returned_masks[group_col][f"{groupby_key}_{interval.left}_{interval.right}"] = mask_to_raster + + # ---------- Raster (mask or segmentation) ---------- + elif isinstance(bins_array, Raster): + data = np.asarray(bins_array.data.data).flatten() + + if len(data) != len(df): + raise ValueError(f"Mask has invalid length {len(data)}, expected {len(df)}") + + # Binary mask + if bins_array.is_mask: + df[group_col] = pd.Categorical(data, categories=[False, True]) + + if returned_masks[group_col] is not None: + mask_to_raster = Raster.from_array(data.reshape(shape).astype("uint8"), transform, crs) + + returned_masks[group_col][f"{groupby_key}_mask"] = mask_to_raster + + # Segmentation + else: + df[group_col] = pd.Categorical(data) + + if returned_masks[group_col] is not None: + segm_names = np.unique(data[~np.isnan(data)]) + for segm_name in segm_names: + mask_to_raster = Raster.from_array((data == segm_name).reshape(shape), transform, crs) + + returned_masks[group_col][f"{groupby_key}_seg_{int(segm_name)}"] = mask_to_raster + + else: + raise NotImplementedError("This type of bins does not yet work") + + groupby_keys.append(group_col) + + # Compute statistics + result = df.groupby(groupby_keys, observed=True)[list(aggregated_vars.keys())].agg(statistics) + + return result, returned_masks diff --git a/tests/test_stats/test_grouped_stats.py b/tests/test_stats/test_grouped_stats.py new file mode 100644 index 000000000..550e1fd50 --- /dev/null +++ b/tests/test_stats/test_grouped_stats.py @@ -0,0 +1,267 @@ +from typing import Any, Dict, List + +import numpy as np +import pandas as pd +import pytest +from pandas.testing import assert_frame_equal + +from geoutils import Raster +from geoutils.stats import grouped_stats + + +@pytest.fixture +def transform(): # type: ignore + return 30.0, 0.0, 478000.0, 0.0, -30.0, 3108140.0 + + +@pytest.fixture +def slope_raster(transform): # type: ignore + data = np.array([[1, 2, 3, 4, 5], [6, 7, 8, 9, 10], [11, 12, 13, 14, 15], [16, 17, 18, 19, 20]]) + return Raster.from_array(data, transform, 32645) + + +@pytest.fixture +def elev_raster(transform): # type: ignore + data = np.array([[1, 2, 3, 4, 5], [6, 7, 8, 9, 10], [11, 12, 13, 14, 15], [16, 17, 18, 19, 20]]) + return Raster.from_array(data, transform, 32645) + + +@pytest.fixture +def mask_raster(tmp_path, transform): # type: ignore + data_mask = np.array( + [ + [True, False, True, False, False], + [True, False, True, False, False], + [True, True, False, False, False], + [True, True, False, False, False], + ] + ) + mask = Raster.from_array(data_mask, transform, 32645) + return mask + + +@pytest.fixture +def segm_raster(transform): # type: ignore + data_segm = np.array([[1, 1, 2, 3, 1], [2, 2, 3, 1, 1], [1, 2, 1, 1, 3], [2, 2, 2, 3, 3]]) + return Raster.from_array(data_segm, transform, 32645) + + +@pytest.fixture +def statistics(): # type: ignore + return ["mean", "min"] + + +@pytest.fixture +def aggregated_vars(elev_raster): # type: ignore + return {"raster": elev_raster} + + +def test_from_raster_to_flattened(elev_raster: Raster) -> None: + dict_test = {"elev_test": elev_raster} + array_dict_test = grouped_stats.from_raster_to_flattened(dict_test) + assert array_dict_test["elev_test"].shape == (20,) + + +def test_is_interval() -> None: + test_interval = [0, 5, 20, 40, np.inf] + assert grouped_stats.is_interval(test_interval) is True + + +@pytest.mark.parametrize( + "interval, type_error, match", + [ + ( + [[0, 5, 20, 40, np.inf], [0, 5, 20, 40, np.inf]], + ValueError, + "If bins is an interval, it must be a 1-dimensional array", + ), + ([0, "a", 20, 40, np.inf], TypeError, "Bins must be a list of number"), + ([0], ValueError, "Bins must be of size >= 2"), + ([0, 5, -2, 40, np.inf], ValueError, "Values must be strictly increasing."), + ], +) # type: ignore +def test_is_interval_with_error(interval: List[Any], type_error: TypeError, match: str) -> None: + + with pytest.raises(type_error, match=match): + _ = grouped_stats.is_interval(interval) + + +@pytest.mark.parametrize( + "bins_test", + [ + ({"slope1": [0, 1, 2]}), # different size of group_by + ({"slope3": [0, 1, 2], "slope4": [0, 1, 2]}), # different keys + ], +) # type: ignore +def test_grouped_stats_errors( + slope_raster: Raster, aggregated_vars: Dict[str, Any], bins_test: Dict[str, Any], statistics: List[Any] +) -> None: + group_by_test = {"slope1": slope_raster, "slope2": slope_raster} + with pytest.raises(ValueError, match="One bins/mask/segmentation entry per input array required."): + _ = grouped_stats.grouped_stats(group_by_test, bins_test, aggregated_vars, statistics) + + +def test_grouped_stats_interval( + slope_raster: Raster, aggregated_vars: Dict[str, Any], statistics: List[Any], elev_raster: Raster +) -> None: + group_by_test = {"slope1": slope_raster} + expected_df = pd.DataFrame( + { + ("raster", "mean"): [3, 8, 15.5], + ("raster", "min"): [1, 6, 11], + }, + index=pd.CategoricalIndex( + pd.IntervalIndex.from_breaks([0.0, 5.0, 10.0, np.inf], closed="right"), + ordered=True, + name="groupby_slope1", + ), + ) + + crs = elev_raster.crs + transform_test = elev_raster.transform + + expected_mask = { + "groupby_slope1": { + "slope1_0.0_5.0": Raster.from_array( + np.array( + [ + [True, True, True, True, True], + [False, False, False, False, False], + [False, False, False, False, False], + [False, False, False, False, False], + ] + ), + transform=transform_test, + crs=crs, + nodata=None, + ), + "slope1_5.0_10.0": Raster.from_array( + np.array( + [ + [False, False, False, False, False], + [True, True, True, True, True], + [False, False, False, False, False], + [False, False, False, False, False], + ] + ), + transform=transform_test, + crs=crs, + nodata=None, + ), + "slope1_10.0_inf": Raster.from_array( + np.array( + [ + [False, False, False, False, False], + [False, False, False, False, False], + [True, True, True, True, True], + [True, True, True, True, True], + ] + ), + transform=transform_test, + crs=crs, + nodata=None, + ), + } + } + + df_test, masks_test = grouped_stats.grouped_stats( + group_by_test, {"slope1": [0, 5, 10, np.inf]}, aggregated_vars, statistics + ) + assert_frame_equal(df_test, expected_df) + + assert all( + expected_mask["groupby_slope1"][group] == masks_test["groupby_slope1"][group] # type: ignore + for group in expected_mask["groupby_slope1"] + ) + + +def test_grouped_stats_mask( + slope_raster: Raster, aggregated_vars: Dict[str, Any], mask_raster: Raster, statistics: List[Any] # type: ignore +) -> None: + group_by_test = {"slope1": slope_raster} + expected_df = pd.DataFrame( + { + ("raster", "mean"): [11.333333, 9.25], + ("raster", "min"): [2, 1], + }, + index=pd.CategoricalIndex( + [False, True], categories=[False, True], ordered=False, dtype="category", name="groupby_slope1" + ), + ) + + df_test, masks_test = grouped_stats.grouped_stats( + group_by_test, {"slope1": mask_raster}, aggregated_vars, statistics + ) + assert_frame_equal(df_test, expected_df) + + expected_mask = {"groupby_slope1": None} + + assert expected_mask == masks_test + + +def test_grouped_stats_mask_nodata( + aggregated_vars: Dict[str, Any], mask_raster: Raster, statistics: List[Any] # type: ignore +) -> None: + + slope_raster = Raster.from_array( + np.array([[1, 2, 3, 4, 5], [6, 7, 8, 9, 10], [11, 12, 13, 14, 15], [16, 17, 18, 19, -999]]), + transform=mask_raster.transform, + crs=mask_raster.crs, + nodata=-9999, + ) + + group_by_test = {"slope1": slope_raster} + expected_df = pd.DataFrame( + { + ("raster", "mean"): [11.333333, 9.25], + ("raster", "min"): [2, 1], + }, + index=pd.CategoricalIndex( + [False, True], categories=[False, True], ordered=False, dtype="category", name="groupby_slope1" + ), + ) + + df_test, masks_test = grouped_stats.grouped_stats( + group_by_test, {"slope1": mask_raster}, aggregated_vars, statistics + ) + assert_frame_equal(df_test, expected_df) + + +def test_grouped_stats_segm( + slope_raster: Raster, aggregated_vars: Dict[str, Any], segm_raster: Raster, statistics: List[Any] +) -> None: + group_by_test = {"segm": slope_raster} + + expected_df = pd.DataFrame( + { + ("raster", "mean"): [8.125000, 11.285714, 13.200000], + ("raster", "min"): [1, 3, 4], + }, + index=pd.CategoricalIndex( + [1, 2, 3], categories=[1, 2, 3], ordered=False, dtype="category", name="groupby_segm" + ), + ) + df_test, masks_test = grouped_stats.grouped_stats(group_by_test, {"segm": segm_raster}, aggregated_vars, statistics) + assert_frame_equal(df_test, expected_df) + + expected_mask = {"groupby_segm": None} + + assert expected_mask == masks_test + + +def test_intersection_stats( + slope_raster: Raster, aggregated_vars: Dict[str, Any], segm_raster: Raster, statistics: List[Any] +) -> None: + group_by_test = {"slope1": slope_raster, "slope2": slope_raster} + bins_test = {"slope1": [0, 5, 10, np.inf], "slope2": segm_raster} + + df_test, _ = grouped_stats.grouped_stats(group_by_test, bins_test, aggregated_vars, statistics) + + assert df_test.loc[(pd.Interval(10.0, np.inf), 3), ("raster", "mean")] == (15 + 19 + 20) / 3 + assert df_test.loc[(pd.Interval(10.0, np.inf), 3), ("raster", "min")] == 15 + + assert df_test.loc[(pd.Interval(0.0, 5), 1), ("raster", "mean")] == (1 + 2 + 5) / 3 + assert df_test.loc[(pd.Interval(0.0, 5), 1), ("raster", "min")] == 1 + + assert df_test.loc[(pd.Interval(10.0, np.inf), 2), ("raster", "mean")] == (16 + 12 + 17 + 18) / 4 + assert df_test.loc[(pd.Interval(10.0, np.inf), 2), ("raster", "min")] == 12 From c5d2a85f12086359798942bb3abefe90aa829077 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alice=20de=20Bardonn=C3=A8che-Richard?= Date: Tue, 3 Mar 2026 16:41:01 +0100 Subject: [PATCH 2/4] New version --- doc/source/stats.md | 68 ++----- geoutils/raster/base.py | 93 +++++++++ geoutils/stats/grouped_stats.py | 171 ---------------- tests/test_raster/test_base.py | 165 +++++++++++++++ tests/test_stats/test_grouped_stats.py | 267 ------------------------- 5 files changed, 275 insertions(+), 489 deletions(-) delete mode 100644 geoutils/stats/grouped_stats.py delete mode 100644 tests/test_stats/test_grouped_stats.py diff --git a/doc/source/stats.md b/doc/source/stats.md index de3ac32d0..0700088da 100644 --- a/doc/source/stats.md +++ b/doc/source/stats.md @@ -106,77 +106,43 @@ Once these bins are created, we reapply them and compute the mean, minimum, and 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. +Using GeoUtils functions makes it very easy to visualise them. ```{code-cell} ipython3 -from geoutils.stats import grouped_stats -import math -import matplotlib.pyplot as plt - -group_by = {"rast": rast} -bins = {"rast": [400, 1000, 2000, 3000, np.inf]} -to_aggregate = {"rast": rast} +group_by = {"raster": rast.data} +bins = [[1000, 2000, 3000]] statistics = ["mean", "min", "max"] -df, masks = grouped_stats.grouped_stats(group_by, bins, to_aggregate, statistics) +df = rast.grouped_stats(group_by, bins, statistics) df ``` -```{code-cell} ipython3 -:tags: [hide-input] -:mystnb: -: code_prompt_show: "Show the code for plotting the figure" -: code_prompt_hide: "Hide the code for plotting the figure" - -groups = list(masks["groupby_rast"].keys()) -n = len(groups) - -ncols = 3 -nrows = math.ceil(n / ncols) - -fig, axes = plt.subplots(nrows, ncols, figsize=(5*ncols, 5*nrows)) -axes = axes.flatten() - -for ax, group in zip(axes, groups): - masks["groupby_rast"][group].plot(ax=ax) - ax.set_title(group) - -for ax in axes[n:]: - ax.axis("off") - -plt.tight_layout() -``` - -```{warning} -Bins can be presented in different ways. It is possible to integrate an interval of minimum 2 values, a mask or a -segmentation map in raster format. -``` - - -### Example with altitude masks +### 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 -from geoutils.stats import grouped_stats +import matplotlib.pyplot as plt -group_by = {"rast": rast} -elev_mask = rast > 2000 -bins = {"rast": elev_mask} +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) -to_aggregate = {"rast": rast} +vect_rasterized.plot() +plt.show() + +group_by = {"elevation": rast.data} +bins = [vect_rasterized.data] statistics = ["mean", "min", "max"] -df, _ = grouped_stats.grouped_stats(group_by, bins, to_aggregate, statistics) +df = rast.grouped_stats(group_by, bins, statistics) df ``` -```{code-cell} ipython3 -elev_mask.plot() -plt.show() -``` ## 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 diff --git a/geoutils/raster/base.py b/geoutils/raster/base.py index 688f0c3bf..06a6ea3c9 100644 --- a/geoutils/raster/base.py +++ b/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, diff --git a/geoutils/stats/grouped_stats.py b/geoutils/stats/grouped_stats.py deleted file mode 100644 index 1dcb54862..000000000 --- a/geoutils/stats/grouped_stats.py +++ /dev/null @@ -1,171 +0,0 @@ -# Copyright (c) 2025 Centre National d'Etudes Spatiales (CNES). -# Copyright (c) 2025 GeoUtils developers -# -# This file is part of the xDEM project: -# https://github.com/glaciohack/xdem -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -""" -This module allows the user to process grouped statistics easily -""" - -from __future__ import annotations - -import logging -from typing import Any - -import numpy as np -import pandas as pd - -from geoutils._typing import NDArrayNum -from geoutils.raster.raster import Raster, RasterType - - -def from_raster_to_flattened(dict_raster: dict[str, Any]) -> dict[str, Any]: - """ - Transform every raster in flattened data to be processed by pandas - :param: Dictionary containing binning name associated to be grouped by and a raster - """ - dict_arrays = dict_raster.copy() - - for key, raster in dict_raster.items(): - dict_arrays[key] = raster.data.ravel() - - return dict_arrays - - -def is_interval(bin_value: RasterType | list[int] | object) -> bool: - """ - Verifying the interval coherence if binning is a list - :param bin_value: Bin value - """ - if isinstance(bin_value, list): - bins_array = np.asarray(bin_value) - else: - return False - - if bins_array.ndim != 1: - raise ValueError("If bins is an interval, it must be a 1-dimensional array") - elif not np.issubdtype(bins_array.dtype, np.number): - raise TypeError("Bins must be a list of number") - elif bins_array.size < 2: - raise ValueError("Bins must be of size >= 2") - elif not np.all(np.diff(bins_array) > 0): - raise ValueError("Values must be strictly increasing.") - - return True - - -def grouped_stats( - groupby_vars: dict[str, RasterType], - bins: dict[str, RasterType | list[int] | object], - aggregated_vars: dict[str, NDArrayNum], - statistics: list[str], -) -> tuple[Any, dict[str, dict[str, NDArrayNum] | None]]: - """ - Get statistics grouped (=binned) by other variables, whether categorical or continuous. - :param groupby_vars: Dictionary containing Raster to group by. - :param bins: Bins to use. Can be a list of interval, binary mask or segmentation map. - :param aggregated_vars: Values to group, can be a dictionary . - :param statistics: List or dict of statistics to compute, e.g. ["mean", "std"]. - """ - - # Check that groupby_vars and bins as the same values - if set(groupby_vars) != set(bins): - raise ValueError("One bins/mask/segmentation entry per input array required.") - - # Load bins - new_bins = {} - for key, value in bins.items(): - if isinstance(value, str): - new_bins[key] = Raster(value) - logging.info(f"No need to save {value} again") - else: - new_bins[key] = value - - raster_base = next(iter(groupby_vars.values())) - crs = raster_base.crs - shape = raster_base.data.shape - transform = raster_base.transform - nodata = raster_base.nodata - - # Flatten arrays for pandas - groupby_vars_array = from_raster_to_flattened(groupby_vars) - aggregated_vars_arrays = from_raster_to_flattened(aggregated_vars) - - df = pd.DataFrame({**groupby_vars_array, **aggregated_vars_arrays}) - if nodata is not None: - df = df.replace(nodata, np.nan) - - groupby_keys = [] - returned_masks = {} - - for groupby_key in groupby_vars.keys(): - group_col = f"groupby_{groupby_key}" - bins_array = new_bins[groupby_key] - - if isinstance(bins[groupby_key], (str, Raster)): - returned_masks[group_col] = None - else: - returned_masks[group_col] = {} - - # ---------- Interval bins ---------- - if is_interval(bins_array): - cut = pd.cut(df[groupby_key], bins=bins_array) - df[group_col] = cut - - if returned_masks[group_col] is not None: - for interval in cut.cat.categories: - mask_flat = (cut == interval).to_numpy() - mask_to_raster = Raster.from_array(mask_flat.reshape(shape), transform, crs) - - returned_masks[group_col][f"{groupby_key}_{interval.left}_{interval.right}"] = mask_to_raster - - # ---------- Raster (mask or segmentation) ---------- - elif isinstance(bins_array, Raster): - data = np.asarray(bins_array.data.data).flatten() - - if len(data) != len(df): - raise ValueError(f"Mask has invalid length {len(data)}, expected {len(df)}") - - # Binary mask - if bins_array.is_mask: - df[group_col] = pd.Categorical(data, categories=[False, True]) - - if returned_masks[group_col] is not None: - mask_to_raster = Raster.from_array(data.reshape(shape).astype("uint8"), transform, crs) - - returned_masks[group_col][f"{groupby_key}_mask"] = mask_to_raster - - # Segmentation - else: - df[group_col] = pd.Categorical(data) - - if returned_masks[group_col] is not None: - segm_names = np.unique(data[~np.isnan(data)]) - for segm_name in segm_names: - mask_to_raster = Raster.from_array((data == segm_name).reshape(shape), transform, crs) - - returned_masks[group_col][f"{groupby_key}_seg_{int(segm_name)}"] = mask_to_raster - - else: - raise NotImplementedError("This type of bins does not yet work") - - groupby_keys.append(group_col) - - # Compute statistics - result = df.groupby(groupby_keys, observed=True)[list(aggregated_vars.keys())].agg(statistics) - - return result, returned_masks diff --git a/tests/test_raster/test_base.py b/tests/test_raster/test_base.py index 0a98accd3..f4ea6959c 100644 --- a/tests/test_raster/test_base.py +++ b/tests/test_raster/test_base.py @@ -7,6 +7,7 @@ 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 @@ -14,6 +15,7 @@ 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) diff --git a/tests/test_stats/test_grouped_stats.py b/tests/test_stats/test_grouped_stats.py deleted file mode 100644 index 550e1fd50..000000000 --- a/tests/test_stats/test_grouped_stats.py +++ /dev/null @@ -1,267 +0,0 @@ -from typing import Any, Dict, List - -import numpy as np -import pandas as pd -import pytest -from pandas.testing import assert_frame_equal - -from geoutils import Raster -from geoutils.stats import grouped_stats - - -@pytest.fixture -def transform(): # type: ignore - return 30.0, 0.0, 478000.0, 0.0, -30.0, 3108140.0 - - -@pytest.fixture -def slope_raster(transform): # type: ignore - data = np.array([[1, 2, 3, 4, 5], [6, 7, 8, 9, 10], [11, 12, 13, 14, 15], [16, 17, 18, 19, 20]]) - return Raster.from_array(data, transform, 32645) - - -@pytest.fixture -def elev_raster(transform): # type: ignore - data = np.array([[1, 2, 3, 4, 5], [6, 7, 8, 9, 10], [11, 12, 13, 14, 15], [16, 17, 18, 19, 20]]) - return Raster.from_array(data, transform, 32645) - - -@pytest.fixture -def mask_raster(tmp_path, transform): # type: ignore - data_mask = np.array( - [ - [True, False, True, False, False], - [True, False, True, False, False], - [True, True, False, False, False], - [True, True, False, False, False], - ] - ) - mask = Raster.from_array(data_mask, transform, 32645) - return mask - - -@pytest.fixture -def segm_raster(transform): # type: ignore - data_segm = np.array([[1, 1, 2, 3, 1], [2, 2, 3, 1, 1], [1, 2, 1, 1, 3], [2, 2, 2, 3, 3]]) - return Raster.from_array(data_segm, transform, 32645) - - -@pytest.fixture -def statistics(): # type: ignore - return ["mean", "min"] - - -@pytest.fixture -def aggregated_vars(elev_raster): # type: ignore - return {"raster": elev_raster} - - -def test_from_raster_to_flattened(elev_raster: Raster) -> None: - dict_test = {"elev_test": elev_raster} - array_dict_test = grouped_stats.from_raster_to_flattened(dict_test) - assert array_dict_test["elev_test"].shape == (20,) - - -def test_is_interval() -> None: - test_interval = [0, 5, 20, 40, np.inf] - assert grouped_stats.is_interval(test_interval) is True - - -@pytest.mark.parametrize( - "interval, type_error, match", - [ - ( - [[0, 5, 20, 40, np.inf], [0, 5, 20, 40, np.inf]], - ValueError, - "If bins is an interval, it must be a 1-dimensional array", - ), - ([0, "a", 20, 40, np.inf], TypeError, "Bins must be a list of number"), - ([0], ValueError, "Bins must be of size >= 2"), - ([0, 5, -2, 40, np.inf], ValueError, "Values must be strictly increasing."), - ], -) # type: ignore -def test_is_interval_with_error(interval: List[Any], type_error: TypeError, match: str) -> None: - - with pytest.raises(type_error, match=match): - _ = grouped_stats.is_interval(interval) - - -@pytest.mark.parametrize( - "bins_test", - [ - ({"slope1": [0, 1, 2]}), # different size of group_by - ({"slope3": [0, 1, 2], "slope4": [0, 1, 2]}), # different keys - ], -) # type: ignore -def test_grouped_stats_errors( - slope_raster: Raster, aggregated_vars: Dict[str, Any], bins_test: Dict[str, Any], statistics: List[Any] -) -> None: - group_by_test = {"slope1": slope_raster, "slope2": slope_raster} - with pytest.raises(ValueError, match="One bins/mask/segmentation entry per input array required."): - _ = grouped_stats.grouped_stats(group_by_test, bins_test, aggregated_vars, statistics) - - -def test_grouped_stats_interval( - slope_raster: Raster, aggregated_vars: Dict[str, Any], statistics: List[Any], elev_raster: Raster -) -> None: - group_by_test = {"slope1": slope_raster} - expected_df = pd.DataFrame( - { - ("raster", "mean"): [3, 8, 15.5], - ("raster", "min"): [1, 6, 11], - }, - index=pd.CategoricalIndex( - pd.IntervalIndex.from_breaks([0.0, 5.0, 10.0, np.inf], closed="right"), - ordered=True, - name="groupby_slope1", - ), - ) - - crs = elev_raster.crs - transform_test = elev_raster.transform - - expected_mask = { - "groupby_slope1": { - "slope1_0.0_5.0": Raster.from_array( - np.array( - [ - [True, True, True, True, True], - [False, False, False, False, False], - [False, False, False, False, False], - [False, False, False, False, False], - ] - ), - transform=transform_test, - crs=crs, - nodata=None, - ), - "slope1_5.0_10.0": Raster.from_array( - np.array( - [ - [False, False, False, False, False], - [True, True, True, True, True], - [False, False, False, False, False], - [False, False, False, False, False], - ] - ), - transform=transform_test, - crs=crs, - nodata=None, - ), - "slope1_10.0_inf": Raster.from_array( - np.array( - [ - [False, False, False, False, False], - [False, False, False, False, False], - [True, True, True, True, True], - [True, True, True, True, True], - ] - ), - transform=transform_test, - crs=crs, - nodata=None, - ), - } - } - - df_test, masks_test = grouped_stats.grouped_stats( - group_by_test, {"slope1": [0, 5, 10, np.inf]}, aggregated_vars, statistics - ) - assert_frame_equal(df_test, expected_df) - - assert all( - expected_mask["groupby_slope1"][group] == masks_test["groupby_slope1"][group] # type: ignore - for group in expected_mask["groupby_slope1"] - ) - - -def test_grouped_stats_mask( - slope_raster: Raster, aggregated_vars: Dict[str, Any], mask_raster: Raster, statistics: List[Any] # type: ignore -) -> None: - group_by_test = {"slope1": slope_raster} - expected_df = pd.DataFrame( - { - ("raster", "mean"): [11.333333, 9.25], - ("raster", "min"): [2, 1], - }, - index=pd.CategoricalIndex( - [False, True], categories=[False, True], ordered=False, dtype="category", name="groupby_slope1" - ), - ) - - df_test, masks_test = grouped_stats.grouped_stats( - group_by_test, {"slope1": mask_raster}, aggregated_vars, statistics - ) - assert_frame_equal(df_test, expected_df) - - expected_mask = {"groupby_slope1": None} - - assert expected_mask == masks_test - - -def test_grouped_stats_mask_nodata( - aggregated_vars: Dict[str, Any], mask_raster: Raster, statistics: List[Any] # type: ignore -) -> None: - - slope_raster = Raster.from_array( - np.array([[1, 2, 3, 4, 5], [6, 7, 8, 9, 10], [11, 12, 13, 14, 15], [16, 17, 18, 19, -999]]), - transform=mask_raster.transform, - crs=mask_raster.crs, - nodata=-9999, - ) - - group_by_test = {"slope1": slope_raster} - expected_df = pd.DataFrame( - { - ("raster", "mean"): [11.333333, 9.25], - ("raster", "min"): [2, 1], - }, - index=pd.CategoricalIndex( - [False, True], categories=[False, True], ordered=False, dtype="category", name="groupby_slope1" - ), - ) - - df_test, masks_test = grouped_stats.grouped_stats( - group_by_test, {"slope1": mask_raster}, aggregated_vars, statistics - ) - assert_frame_equal(df_test, expected_df) - - -def test_grouped_stats_segm( - slope_raster: Raster, aggregated_vars: Dict[str, Any], segm_raster: Raster, statistics: List[Any] -) -> None: - group_by_test = {"segm": slope_raster} - - expected_df = pd.DataFrame( - { - ("raster", "mean"): [8.125000, 11.285714, 13.200000], - ("raster", "min"): [1, 3, 4], - }, - index=pd.CategoricalIndex( - [1, 2, 3], categories=[1, 2, 3], ordered=False, dtype="category", name="groupby_segm" - ), - ) - df_test, masks_test = grouped_stats.grouped_stats(group_by_test, {"segm": segm_raster}, aggregated_vars, statistics) - assert_frame_equal(df_test, expected_df) - - expected_mask = {"groupby_segm": None} - - assert expected_mask == masks_test - - -def test_intersection_stats( - slope_raster: Raster, aggregated_vars: Dict[str, Any], segm_raster: Raster, statistics: List[Any] -) -> None: - group_by_test = {"slope1": slope_raster, "slope2": slope_raster} - bins_test = {"slope1": [0, 5, 10, np.inf], "slope2": segm_raster} - - df_test, _ = grouped_stats.grouped_stats(group_by_test, bins_test, aggregated_vars, statistics) - - assert df_test.loc[(pd.Interval(10.0, np.inf), 3), ("raster", "mean")] == (15 + 19 + 20) / 3 - assert df_test.loc[(pd.Interval(10.0, np.inf), 3), ("raster", "min")] == 15 - - assert df_test.loc[(pd.Interval(0.0, 5), 1), ("raster", "mean")] == (1 + 2 + 5) / 3 - assert df_test.loc[(pd.Interval(0.0, 5), 1), ("raster", "min")] == 1 - - assert df_test.loc[(pd.Interval(10.0, np.inf), 2), ("raster", "mean")] == (16 + 12 + 17 + 18) / 4 - assert df_test.loc[(pd.Interval(10.0, np.inf), 2), ("raster", "min")] == 12 From 80da2e86ae91a04b24f3fbbb074f268488fdb004 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alice=20de=20Bardonn=C3=A8che-Richard?= Date: Thu, 5 Mar 2026 17:26:19 +0100 Subject: [PATCH 3/4] doc improvment --- doc/source/stats.md | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/doc/source/stats.md b/doc/source/stats.md index 0700088da..7c8a1073b 100644 --- a/doc/source/stats.md +++ b/doc/source/stats.md @@ -101,7 +101,7 @@ 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]. +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. @@ -126,10 +126,8 @@ for masks = True. It is also possible to use a reference other than the raster i ```{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) +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() From df44b998443ae543786fef71ed8948b30efe0dd2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alice=20de=20Bardonn=C3=A8che-Richard?= Date: Fri, 6 Mar 2026 09:22:28 +0100 Subject: [PATCH 4/4] docstring --- geoutils/raster/base.py | 52 ++++++++++++++++++++++++++++++++++------- 1 file changed, 43 insertions(+), 9 deletions(-) diff --git a/geoutils/raster/base.py b/geoutils/raster/base.py index 06a6ea3c9..0d0a1bc9d 100644 --- a/geoutils/raster/base.py +++ b/geoutils/raster/base.py @@ -817,10 +817,44 @@ def grouped_stats( def grouped_stats( self, groupby_arrays: dict[str, NDArrayNum], - bins: list[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.") @@ -1210,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. @@ -1315,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. @@ -1566,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(