diff --git a/doc/source/stats.md b/doc/source/stats.md index f640a788..7c8a1073 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,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 diff --git a/geoutils/raster/base.py b/geoutils/raster/base.py index 688f0c3b..0d0a1bc9 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,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( diff --git a/tests/test_raster/test_base.py b/tests/test_raster/test_base.py index 0a98accd..f4ea6959 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)