Merge spatial_tools.merge_bounding_boxes and projtools.merge_bounds functions

geoutils/projtools.py

@@ -105,20 +105,27 @@ def bounds2poly(
 
 
 def merge_bounds(
-    bounds_list: abc.Iterable[list[float] | rio.io.DatasetReader | gpd.GeoDataFrame],
+    bounds_list: abc.Iterable[
+        list[float] | tuple[float] | rio.coords.BoundingBox | rio.io.DatasetReader | gpd.GeoDataFrame
+    ],
+    resolution: float | None = None,
     merging_algorithm: str = "union",
-) -> tuple[float, ...]:
+    return_rio_bbox: bool = False,
+) -> tuple[float, ...] | rio.coords.BoundingBox:
     """
     Merge a list of bounds into single bounds, using either the union or intersection.
 
-    :param bounds_list: A list of geometries with bounds, i.e. a list of coordinates (xmin, ymin, xmax, ymax), \
-a rasterio/Raster object, a geoPandas/Vector object.
-    :param merging_algorithm: the algorithm to use for merging, either "union" or "intersection"
+    :param bounds_list: List of geometries with bounds, i.e. list of coordinates (xmin, ymin, xmax, ymax),
+        rasterio bounds, a rasterio Dataset (or Raster), a geopandas object (or Vector).
+    :param resolution: (For Rasters) Resolution, to make sure extent is a multiple of it.
+    :param merging_algorithm: Algorithm to use for merging, either "union" or "intersection".
+    :param return_rio_bbox: Whether to return a rio.coords.BoundingBox object instead of a tuple.
 
     :returns: Output bounds (xmin, ymin, xmax, ymax) or empty tuple
     """
     # Check that bounds_list is a list of bounds objects
     assert isinstance(bounds_list, (list, tuple)), "bounds_list must be a list/tuple"
+
     for bounds in bounds_list:
         assert hasattr(bounds, "bounds") or hasattr(bounds, "total_bounds") or isinstance(bounds, (list, tuple)), (
             "bounds_list must be a list of lists/tuples of coordinates or an object with attributes bounds "
@@ -137,8 +144,25 @@ def merge_bounds(
         else:
             raise ValueError("merging_algorithm must be 'union' or 'intersection'")
 
-    new_bounds: tuple[float] = output_poly.bounds
-    return new_bounds
+    new_bounds = output_poly.bounds
+
+    rio_bounds = {"left": new_bounds[0], "bottom": new_bounds[1], "right": new_bounds[2], "top": new_bounds[3]}
+
+    print(rio_bounds)
+
+    # Make sure that extent is a multiple of resolution
+    if resolution is not None:
+        for key1, key2 in zip(("left", "bottom"), ("right", "top")):
+            modulo = (rio_bounds[key2] - rio_bounds[key1]) % resolution
+            rio_bounds[key2] += modulo
+
+    # Format output
+    if return_rio_bbox:
+        final_bounds = rio.coords.BoundingBox(**rio_bounds)
+    else:
+        final_bounds = tuple(rio_bounds.values())
+
+    return final_bounds
 
 
 def align_bounds(

geoutils/spatial_tools.py

@@ -132,7 +132,7 @@ def load_multiple_rasters(
         bounds.append(bound)
 
     # Second get the intersection of all raster bounds
-    intersection = gu.projtools.merge_bounds(bounds, "intersection")
+    intersection = gu.projtools.merge_bounds(bounds, merging_algorithm="intersection")
 
     # Optionally, crop the rasters
     if crop:
@@ -182,22 +182,6 @@ def load_multiple_rasters(
     return output_rst
 
 
-def merge_bounding_boxes(bounds: list[rio.coords.BoundingBox], resolution: float) -> rio.coords.BoundingBox:
-    max_bounds = dict(zip(["left", "right", "top", "bottom"], [np.nan] * 4))
-    for bound in bounds:
-        for key in "right", "top":
-            max_bounds[key] = np.nanmax([max_bounds[key], bound.__getattribute__(key)])
-        for key in "bottom", "left":
-            max_bounds[key] = np.nanmin([max_bounds[key], bound.__getattribute__(key)])
-
-    # Make sure that extent is a multiple of resolution
-    for key1, key2 in zip(("left", "bottom"), ("right", "top")):
-        modulo = (max_bounds[key2] - max_bounds[key1]) % resolution
-        max_bounds[key2] += modulo
-
-    return rio.coords.BoundingBox(**max_bounds)
-
-
 def stack_rasters(
     rasters: list[RasterType],
     reference: int | gu.Raster = 0,
@@ -246,7 +230,9 @@ def stack_rasters(
     if use_ref_bounds:
         dst_bounds = reference_raster.bounds
     else:
-        dst_bounds = merge_bounding_boxes([raster.bounds for raster in rasters], resolution=reference_raster.res[0])
+        dst_bounds = gu.projtools.merge_bounds(
+            [raster.bounds for raster in rasters], resolution=reference_raster.res[0], return_rio_bbox=True
+        )
 
     # Make a data list and add all of the reprojected rasters into it.
     data: list[np.ndarray] = []

tests/test_spatial_tools.py

@@ -198,7 +198,7 @@ def test_stack_rasters(rasters) -> None:  # type: ignore
     assert type(stacked_img) == gu.Raster  # Check output object is always Raster, whatever input was given
     assert np.count_nonzero(np.isnan(stacked_img.data)) == 0  # Check no NaNs introduced
 
-    merged_bounds = gu.spatial_tools.merge_bounding_boxes(
+    merged_bounds = gu.projtools.merge_bounds(
         [rasters.img1.bounds, rasters.img2.bounds], resolution=rasters.img1.res[0]
     )
     assert merged_bounds == stacked_img.bounds