From f320528f54b4192144021e359efc00b3d10a165d Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Wed, 16 Apr 2025 21:28:23 -0600
Subject: [PATCH 01/25] WIP align implementation

---
 src/rasterix/raster_index.py | 76 ++++++++++++++++++++++++++++++++++++
 1 file changed, 76 insertions(+)

diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index 25bfe3c..263608c 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -32,6 +32,50 @@ def assign_index(obj: T_Xarray, *, x_dim: str | None = None, y_dim: str | None =
     return obj.assign_coords(coords)
 
 
+def _assert_transforms_are_compatible(A1: Affine, A2: Affine) -> None:
+    assert A1.a == A2.a
+    assert A1.b == A2.b
+    assert A1.d == A2.d
+    assert A1.e == A2.e
+
+
+def compute_bounding_box_affine(
+    A1: Affine, shape1: tuple[int, int], A2: Affine, shape2: tuple[int, int], *, how: str
+):
+    _assert_transforms_are_compatible(A1, A2)
+
+    # Apply both affine transformations to the corners
+    transformed1 = np.stack([A1 * point for point in [(0, 0), shape1]])
+    transformed2 = np.stack([A2 * point for point in [(0, 0), shape2]])
+
+    # Concatenate transformed points
+    all_points = np.stack([transformed1, transformed2])
+
+    # Compute bounding box min/max
+    if how == "outer":
+        min_x, min_y = all_points.min(axis=(0, 1)).tolist()
+        max_x, max_y = all_points.max(axis=(0, 1)).tolist()
+
+        assert max_x > min_x
+        assert max_y > min_y
+
+        # FIXME: floating point error here
+        Nx = None if A1.a == 0 else int(np.abs((max_x - min_x) / A1.a))
+        Ny = None if A1.e == 0 else int(np.abs((max_y - min_y) / A1.e))
+    else:
+        raise NotImplementedError
+
+    Anew = Affine(
+        A1.a,
+        A1.b,
+        min_x if A1.a > 0 else max_x,
+        A1.d,
+        A1.e,
+        min_y if A1.e > 0 else max_y,  # TODO: handle orientation here
+    )
+    return Anew, Nx, Ny
+
+
 class AffineTransform(CoordinateTransform):
     """Affine 2D transform wrapper."""
 
@@ -304,6 +348,7 @@ class RasterIndex(Index):
     """
 
     _wrapped_indexes: dict[WrappedIndexCoords, WrappedIndex]
+    _shape: dict[str, int]
 
     def __init__(self, indexes: Mapping[WrappedIndexCoords, WrappedIndex]):
         idx_keys = list(indexes)
@@ -321,6 +366,13 @@ def __init__(self, indexes: Mapping[WrappedIndexCoords, WrappedIndex]):
         assert axis_dependent ^ axis_independent
 
         self._wrapped_indexes = dict(indexes)
+        if axis_independent:
+            self._shape = {
+                "x": self._wrapped_indexes["x"].axis_transform.size,
+                "y": self._wrapped_indexes["y"].axis_transform.size,
+            }
+        else:
+            self._shape = next(iter(self._wrapped_indexes.values())).dim_size
 
     @classmethod
     def from_transform(
@@ -436,3 +488,27 @@ def transform(self) -> Affine:
             index = next(iter(self._wrapped_indexes.values()))
             aff = index.affine
         return aff * Affine.translation(-0.5, -0.5)
+
+    def join(self, other, how) -> "RasterIndex":
+        if not isinstance(other, RasterIndex):
+            raise ValueError(
+                f"Alignment is only supported between RasterIndexes. Received RasterIndex and {type(other)!r} instead"
+            )
+
+        if set(self._wrapped_indexes.keys()) != set(other._wrapped_indexes.keys()):
+            # TODO: better error message
+            raise ValueError(
+                "Alignment is only supported between RasterIndexes, when both contain compatible transforms."
+            )
+
+        # move affines back to top-left corner
+        trans = Affine.translation(-0.5, -0.5)
+        ours = self.transform() * trans
+        theirs = other.transform() * trans
+
+        our_shape = tuple(self._shape[k] for k in ("x", "y"))
+        their_shape = tuple(self._shape[k] for k in ("x", "y"))
+        new_affine, Nx, Ny = compute_bounding_box_affine(ours, our_shape, theirs, their_shape, how=how)
+
+        # TODO: set xdim, ydim explicitly
+        return self.from_transform(new_affine, width=Nx, height=Ny)

From 3665bc7fe11c20b15b74f12cbd65a3ccdfc8a88d Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Thu, 17 Apr 2025 09:49:14 -0600
Subject: [PATCH 02/25] Try using BoundingBox

---
 src/rasterix/raster_index.py | 77 ++++++++++++++++++++++++++++++++++--
 1 file changed, 74 insertions(+), 3 deletions(-)

diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index 263608c..75511b8 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -3,11 +3,16 @@
 import textwrap
 from collections.abc import Hashable, Mapping
 from typing import Any, TypeVar
+from typing import Any, Self
 
 import numpy as np
 import pandas as pd
 from affine import Affine
 from xarray import Coordinates, DataArray, Dataset, Index, Variable
+from xarray import DataArray, Index, Variable
+from odc.geo import BoundingBox
+from odc.geo.geom import bbox_intersection, bbox_union
+from xarray import DataArray, Index, Variable
 from xarray.core.coordinate_transform import CoordinateTransform
 
 # TODO: import from public API once it is available
@@ -489,6 +494,14 @@ def transform(self) -> Affine:
             aff = index.affine
         return aff * Affine.translation(-0.5, -0.5)
 
+    @property
+    def bbox(self) -> BoundingBox:
+        return BoundingBox.from_transform(
+            shape=tuple(self._shape[k] for k in ("x", "y")),
+            transform=self.combined_affine_transform() * Affine.translation(-0.5, -0.5),
+            crs=None,
+        )
+
     def join(self, other, how) -> "RasterIndex":
         if not isinstance(other, RasterIndex):
             raise ValueError(
@@ -505,10 +518,68 @@ def join(self, other, how) -> "RasterIndex":
         trans = Affine.translation(-0.5, -0.5)
         ours = self.transform() * trans
         theirs = other.transform() * trans
+        ours, theirs = as_compatible_bboxes(self, other)
+
+        if how == "outer":
+            new_bbox = bbox_union([ours, theirs])
+            print(new_bbox, ours, theirs)
+        elif how == "inner":
+            new_bbox = bbox_intersection([ours, theirs])
 
-        our_shape = tuple(self._shape[k] for k in ("x", "y"))
-        their_shape = tuple(self._shape[k] for k in ("x", "y"))
-        new_affine, Nx, Ny = compute_bounding_box_affine(ours, our_shape, theirs, their_shape, how=how)
+        affine = self.combined_affine_transform()
+        new_affine, Nx, Ny = bbox_to_affine(new_bbox, rx=affine.a, ry=affine.e)
 
         # TODO: set xdim, ydim explicitly
         return self.from_transform(new_affine, width=Nx, height=Ny)
+
+    def reindex_like(self, other: Self, method=None, tolerance=None) -> dict[Hashable, Any]:
+        affine = self.combined_affine_transform()
+        ours, theirs = as_compatible_bboxes(self, other)
+        inter = bbox_intersection([ours, theirs])
+        dx = affine.a
+        dy = affine.e
+
+        # Fraction of a pixel that can be ignored, defaults to 1/100. Bounding box of the output
+        # geobox is allowed to be smaller than supplied bounding box by that amount.
+        # FIXME: translate user-provided `tolerance` to `tol`
+        tol: float = 0.01
+
+        indexers = {}
+        indexers["x"] = slice_bounds(theirs.left, inter.left, inter.right, spacing=dx, tol=tol)
+        indexers["y"] = slice_bounds(theirs.top, inter.bottom, inter.top, spacing=dy, tol=tol)
+        return indexers
+
+
+def slice_bounds(off, start, stop, spacing, tol):
+    from math import ceil
+
+    from odc.geo.math import maybe_int
+
+    istart = ceil(maybe_int((start - off) / spacing, tol))
+    istop = ceil(maybe_int((stop - off) / spacing, tol))
+
+    return slice(istart, istop)
+
+
+def bbox_to_affine(bbox: BoundingBox, rx, ry) -> Affine:
+    # Fraction of a pixel that can be ignored, defaults to 1/100. Bounding box of the output
+    # geobox is allowed to be smaller than supplied bounding box by that amount.
+    # FIXME: translate user-provided `tolerance` to `tol`
+    tol: float = 0.01
+
+    from odc.geo.math import snap_grid
+
+    offx, nx = snap_grid(bbox.left, bbox.right, rx, 0, tol=tol)
+    offy, ny = snap_grid(bbox.bottom, bbox.top, ry, 0, tol=tol)
+
+    affine = Affine.translation(offx, offy) * Affine.scale(rx, ry)
+
+    return affine, nx, ny
+
+
+def as_compatible_bboxes(r1: RasterIndex, r2: RasterIndex) -> tuple[BoundingBox, BoundingBox]:
+    r1_transform = r1.combined_affine_transform()
+    r2_transform = r2.combined_affine_transform()
+    _assert_transforms_are_compatible(r1_transform, r2_transform)
+
+    return r1.bbox, r2.bbox

From 56c6308ac6a95be0390afe670a641585573bc40b Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Thu, 17 Apr 2025 11:32:16 -0600
Subject: [PATCH 03/25] bugfix

---
 src/rasterix/raster_index.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index 75511b8..163e869 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -497,7 +497,7 @@ def transform(self) -> Affine:
     @property
     def bbox(self) -> BoundingBox:
         return BoundingBox.from_transform(
-            shape=tuple(self._shape[k] for k in ("x", "y")),
+            shape=tuple(self._shape[k] for k in ("y", "x")),
             transform=self.combined_affine_transform() * Affine.translation(-0.5, -0.5),
             crs=None,
         )
@@ -530,7 +530,9 @@ def join(self, other, how) -> "RasterIndex":
         new_affine, Nx, Ny = bbox_to_affine(new_bbox, rx=affine.a, ry=affine.e)
 
         # TODO: set xdim, ydim explicitly
-        return self.from_transform(new_affine, width=Nx, height=Ny)
+        new_index = self.from_transform(new_affine, width=Nx, height=Ny)
+        assert new_index.bbox == new_bbox
+        return new_index
 
     def reindex_like(self, other: Self, method=None, tolerance=None) -> dict[Hashable, Any]:
         affine = self.combined_affine_transform()

From bf847183cbc66fcfc9cc4bb2ae19098c35840683 Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Thu, 17 Apr 2025 12:14:19 -0600
Subject: [PATCH 04/25] Some more progress

---
 src/rasterix/raster_index.py | 39 +++++++++++++++++++++++++++++++-----
 1 file changed, 34 insertions(+), 5 deletions(-)

diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index 163e869..709a801 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -3,6 +3,8 @@
 import textwrap
 from collections.abc import Hashable, Mapping
 from typing import Any, TypeVar
+from collections.abc import Hashable, Mapping
+from collections.abc import Hashable, Iterable, Mapping, Sequence
 from typing import Any, Self
 
 import numpy as np
@@ -502,6 +504,17 @@ def bbox(self) -> BoundingBox:
             crs=None,
         )
 
+    @classmethod
+    def concat(
+        cls,
+        indexes: Sequence[Self],
+        dim: Hashable,
+        positions: Iterable[Iterable[int]] | None = None,
+    ) -> Self:
+        if len(indexes) == 1:
+            return next(iter(indexes))
+        raise NotImplementedError
+
     def join(self, other, how) -> "RasterIndex":
         if not isinstance(other, RasterIndex):
             raise ValueError(
@@ -547,20 +560,36 @@ def reindex_like(self, other: Self, method=None, tolerance=None) -> dict[Hashabl
         tol: float = 0.01
 
         indexers = {}
-        indexers["x"] = slice_bounds(theirs.left, inter.left, inter.right, spacing=dx, tol=tol)
-        indexers["y"] = slice_bounds(theirs.top, inter.bottom, inter.top, spacing=dy, tol=tol)
+        indexers["x"] = get_indexer(
+            theirs.left, ours.left, inter.left, inter.right, spacing=dx, tol=tol, size=other._shape["x"]
+        )
+        indexers["y"] = get_indexer(
+            theirs.top, ours.top, inter.top, inter.bottom, spacing=dy, tol=tol, size=other._shape["y"]
+        )
+        import ipdb
+
+        ipdb.set_trace()
         return indexers
 
 
-def slice_bounds(off, start, stop, spacing, tol):
+def get_indexer(off, our_off, start, stop, spacing, tol, size) -> np.ndarray:
     from math import ceil
 
     from odc.geo.math import maybe_int
 
     istart = ceil(maybe_int((start - off) / spacing, tol))
-    istop = ceil(maybe_int((stop - off) / spacing, tol))
 
-    return slice(istart, istop)
+    ours_istart = ceil(maybe_int((start - our_off) / spacing, tol))
+    ours_istop = ceil(maybe_int((stop - our_off) / spacing, tol))
+
+    idxr = np.concatenate(
+        [
+            np.full((istart,), fill_value=-1),
+            np.arange(ours_istart, ours_istop),
+            np.full((size - istart - (ours_istop - ours_istart),), fill_value=-1),
+        ]
+    )
+    return idxr
 
 
 def bbox_to_affine(bbox: BoundingBox, rx, ry) -> Affine:

From a10bfe7e2390f202c821accc888fa8e406994230 Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Fri, 16 May 2025 09:23:25 -0600
Subject: [PATCH 05/25] cleanup

---
 src/rasterix/raster_index.py | 23 ++++++++++-------------
 1 file changed, 10 insertions(+), 13 deletions(-)

diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index 709a801..8a2e569 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -1,20 +1,15 @@
 from __future__ import annotations
 
 import textwrap
-from collections.abc import Hashable, Mapping
-from typing import Any, TypeVar
-from collections.abc import Hashable, Mapping
 from collections.abc import Hashable, Iterable, Mapping, Sequence
-from typing import Any, Self
+from typing import Any, Self, TypeVar
 
 import numpy as np
 import pandas as pd
 from affine import Affine
-from xarray import Coordinates, DataArray, Dataset, Index, Variable
-from xarray import DataArray, Index, Variable
 from odc.geo import BoundingBox
 from odc.geo.geom import bbox_intersection, bbox_union
-from xarray import DataArray, Index, Variable
+from xarray import Coordinates, DataArray, Dataset, Index, Variable
 from xarray.core.coordinate_transform import CoordinateTransform
 
 # TODO: import from public API once it is available
@@ -500,7 +495,7 @@ def transform(self) -> Affine:
     def bbox(self) -> BoundingBox:
         return BoundingBox.from_transform(
             shape=tuple(self._shape[k] for k in ("y", "x")),
-            transform=self.combined_affine_transform() * Affine.translation(-0.5, -0.5),
+            transform=self.transform() * Affine.translation(-0.5, -0.5),
             crs=None,
         )
 
@@ -515,7 +510,7 @@ def concat(
             return next(iter(indexes))
         raise NotImplementedError
 
-    def join(self, other, how) -> "RasterIndex":
+    def join(self, other, how) -> RasterIndex:
         if not isinstance(other, RasterIndex):
             raise ValueError(
                 f"Alignment is only supported between RasterIndexes. Received RasterIndex and {type(other)!r} instead"
@@ -538,8 +533,10 @@ def join(self, other, how) -> "RasterIndex":
             print(new_bbox, ours, theirs)
         elif how == "inner":
             new_bbox = bbox_intersection([ours, theirs])
+        else:
+            raise NotImplementedError
 
-        affine = self.combined_affine_transform()
+        affine = self.transform()
         new_affine, Nx, Ny = bbox_to_affine(new_bbox, rx=affine.a, ry=affine.e)
 
         # TODO: set xdim, ydim explicitly
@@ -548,7 +545,7 @@ def join(self, other, how) -> "RasterIndex":
         return new_index
 
     def reindex_like(self, other: Self, method=None, tolerance=None) -> dict[Hashable, Any]:
-        affine = self.combined_affine_transform()
+        affine = self.transform()
         ours, theirs = as_compatible_bboxes(self, other)
         inter = bbox_intersection([ours, theirs])
         dx = affine.a
@@ -609,8 +606,8 @@ def bbox_to_affine(bbox: BoundingBox, rx, ry) -> Affine:
 
 
 def as_compatible_bboxes(r1: RasterIndex, r2: RasterIndex) -> tuple[BoundingBox, BoundingBox]:
-    r1_transform = r1.combined_affine_transform()
-    r2_transform = r2.combined_affine_transform()
+    r1_transform = r1.transform()
+    r2_transform = r2.transform()
     _assert_transforms_are_compatible(r1_transform, r2_transform)
 
     return r1.bbox, r2.bbox

From 0ee5616d59cde1e307ffe8ff7af2c8864461eb7a Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Fri, 16 May 2025 09:23:51 -0600
Subject: [PATCH 06/25] remove old bbox code

---
 src/rasterix/raster_index.py | 37 ------------------------------------
 1 file changed, 37 deletions(-)

diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index 8a2e569..c10af3d 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -41,43 +41,6 @@ def _assert_transforms_are_compatible(A1: Affine, A2: Affine) -> None:
     assert A1.e == A2.e
 
 
-def compute_bounding_box_affine(
-    A1: Affine, shape1: tuple[int, int], A2: Affine, shape2: tuple[int, int], *, how: str
-):
-    _assert_transforms_are_compatible(A1, A2)
-
-    # Apply both affine transformations to the corners
-    transformed1 = np.stack([A1 * point for point in [(0, 0), shape1]])
-    transformed2 = np.stack([A2 * point for point in [(0, 0), shape2]])
-
-    # Concatenate transformed points
-    all_points = np.stack([transformed1, transformed2])
-
-    # Compute bounding box min/max
-    if how == "outer":
-        min_x, min_y = all_points.min(axis=(0, 1)).tolist()
-        max_x, max_y = all_points.max(axis=(0, 1)).tolist()
-
-        assert max_x > min_x
-        assert max_y > min_y
-
-        # FIXME: floating point error here
-        Nx = None if A1.a == 0 else int(np.abs((max_x - min_x) / A1.a))
-        Ny = None if A1.e == 0 else int(np.abs((max_y - min_y) / A1.e))
-    else:
-        raise NotImplementedError
-
-    Anew = Affine(
-        A1.a,
-        A1.b,
-        min_x if A1.a > 0 else max_x,
-        A1.d,
-        A1.e,
-        min_y if A1.e > 0 else max_y,  # TODO: handle orientation here
-    )
-    return Anew, Nx, Ny
-
-
 class AffineTransform(CoordinateTransform):
     """Affine 2D transform wrapper."""
 

From 5909fa8ce58df4e0eae00ec4826d06bf47736590 Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Fri, 16 May 2025 09:56:19 -0600
Subject: [PATCH 07/25] Adapt to new equals

---
 src/rasterix/raster_index.py | 18 +++++++++++-------
 tests/test_raster_index.py   | 19 +++++++++++++++++++
 2 files changed, 30 insertions(+), 7 deletions(-)

diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index c10af3d..a0bc4c4 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -85,7 +85,9 @@ def reverse(self, coord_labels):
 
         return results
 
-    def equals(self, other):
+    def equals(self, other, *, exclude):
+        if exclude is not None:
+            raise NotImplementedError
         if not isinstance(other, AffineTransform):
             return False
         return self.affine == other.affine and self.dim_size == other.dim_size
@@ -138,9 +140,11 @@ def reverse(self, coord_labels: dict[Hashable, Any]) -> dict[str, Any]:
 
         return {self.dim: positions}
 
-    def equals(self, other):
+    def equals(self, other, *, exclude):
         if not isinstance(other, AxisAffineTransform):
             return False
+        if exclude is not None:
+            raise NotImplementedError
 
         # only compare the affine parameters of the relevant axis
         if self.is_xaxis:
@@ -413,7 +417,9 @@ def sel(self, labels: dict[Any, Any], method=None, tolerance=None) -> IndexSelRe
 
         return merge_sel_results(results)
 
-    def equals(self, other: Index) -> bool:
+    def equals(self, other: Index, *, exclude=None) -> bool:
+        if exclude is None:
+            exclude = {}
         if not isinstance(other, RasterIndex):
             return False
         if set(self._wrapped_indexes) != set(other._wrapped_indexes):
@@ -422,6 +428,7 @@ def equals(self, other: Index) -> bool:
         return all(
             index.equals(other._wrapped_indexes[k])  # type: ignore[arg-type]
             for k, index in self._wrapped_indexes.items()
+            if k not in exclude
         )
 
     def to_pandas_index(self) -> pd.Index:
@@ -473,7 +480,7 @@ def concat(
             return next(iter(indexes))
         raise NotImplementedError
 
-    def join(self, other, how) -> RasterIndex:
+    def join(self, other: RasterIndex, how) -> RasterIndex:
         if not isinstance(other, RasterIndex):
             raise ValueError(
                 f"Alignment is only supported between RasterIndexes. Received RasterIndex and {type(other)!r} instead"
@@ -526,9 +533,6 @@ def reindex_like(self, other: Self, method=None, tolerance=None) -> dict[Hashabl
         indexers["y"] = get_indexer(
             theirs.top, ours.top, inter.top, inter.bottom, spacing=dy, tol=tol, size=other._shape["y"]
         )
-        import ipdb
-
-        ipdb.set_trace()
         return indexers
 
 
diff --git a/tests/test_raster_index.py b/tests/test_raster_index.py
index 045ae3a..663f1b7 100644
--- a/tests/test_raster_index.py
+++ b/tests/test_raster_index.py
@@ -1,4 +1,5 @@
 import numpy as np
+import pyproj
 import rioxarray  # noqa
 import xarray as xr
 from affine import Affine
@@ -31,3 +32,21 @@ def test_sel_slice():
 
     assert actual_transform == actual.rio.transform()
     assert actual_transform == (transform * Affine.translation(0, 3))
+
+
+def test_combine_nested():
+    transforms = [
+        "-50.0 0.5 0.0 0.0 0.0 -0.25",
+        "-40.0 0.5 0.0 0.0 0.0 -0.25",
+    ]
+    crs_attrs = pyproj.CRS.from_epsg(4326).to_cf()
+
+    datasets = [
+        xr.Dataset(
+            {"foo": (("y", "x"), np.ones((4, 2)), {"grid_mapping": "spatial_ref"})},
+            coords={"spatial_ref": ((), 0, crs_attrs | {"GeoTransform": transform})},
+        )
+        for transform in transforms
+    ]
+    datasets = list(map(assign_index, datasets))
+    xr.combine_nested(datasets, concat_dim="x")

From fbe3b371cee2a5ed51dceaf83eb0f99c3789550e Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Fri, 16 May 2025 10:13:32 -0600
Subject: [PATCH 08/25] Vendor odc's BoundingBox

---
 src/rasterix/odc_compat.py   | 241 +++++++++++++++++++++++++++++++++++
 src/rasterix/raster_index.py |   4 +-
 2 files changed, 242 insertions(+), 3 deletions(-)
 create mode 100644 src/rasterix/odc_compat.py

diff --git a/src/rasterix/odc_compat.py b/src/rasterix/odc_compat.py
new file mode 100644
index 0000000..c67ecaa
--- /dev/null
+++ b/src/rasterix/odc_compat.py
@@ -0,0 +1,241 @@
+import itertools
+import math
+from collections.abc import Iterable, Sequence
+from typing import Any
+
+from affine import Affine
+
+
+class BoundingBox(Sequence[float]):
+    """Bounding box, defining extent in cartesian coordinates."""
+
+    __slots__ = "_box"
+
+    def __init__(self, left: float, bottom: float, right: float, top: float):
+        self._box = (left, bottom, right, top)
+
+    @property
+    def left(self):
+        return self._box[0]
+
+    @property
+    def bottom(self):
+        return self._box[1]
+
+    @property
+    def right(self):
+        return self._box[2]
+
+    @property
+    def top(self):
+        return self._box[3]
+
+    @property
+    def bbox(self) -> tuple[float, float, float, float]:
+        return self._box
+
+    def __iter__(self):
+        return self._box.__iter__()
+
+    def __eq__(self, other: Any) -> bool:
+        if isinstance(other, BoundingBox):
+            return self._box == other._box
+        return self._box == other
+
+    def __hash__(self) -> int:
+        return hash(self._box)
+
+    def __len__(self) -> int:
+        return 4
+
+    def __getitem__(self, idx):
+        return self._box[idx]
+
+    def __repr__(self) -> str:
+        return f"BoundingBox(left={self.left}, bottom={self.bottom}, right={self.right}, top={self.top})"
+
+    def __str__(self) -> str:
+        return self.__repr__()
+
+    def __and__(self, other: "BoundingBox") -> "BoundingBox":
+        return bbox_intersection([self, other])
+
+    def __or__(self, other: "BoundingBox") -> "BoundingBox":
+        return bbox_union([self, other])
+
+    def buffered(self, xbuff: float, ybuff: float | None = None) -> "BoundingBox":
+        """
+        Return a new BoundingBox, buffered in the x and y dimensions.
+
+        :param xbuff: X dimension buffering amount
+        :param ybuff: Y dimension buffering amount
+        :return: new BoundingBox
+        """
+        if ybuff is None:
+            ybuff = xbuff
+
+        return BoundingBox(
+            left=self.left - xbuff,
+            bottom=self.bottom - ybuff,
+            right=self.right + xbuff,
+            top=self.top + ybuff,
+        )
+
+    @property
+    def span_x(self) -> float:
+        """Span of the bounding box along x axis."""
+        return self.right - self.left
+
+    @property
+    def span_y(self) -> float:
+        """Span of the bounding box along y axis."""
+        return self.top - self.bottom
+
+    @property
+    def aspect(self) -> float:
+        """Aspect ratio."""
+        return self.span_x / self.span_y
+
+    @property
+    def width(self) -> int:
+        """``int(span_x)``"""
+        return int(self.right - self.left)
+
+    @property
+    def height(self) -> int:
+        """``int(span_y)``"""
+        return int(self.top - self.bottom)
+
+    @property
+    def shape(self) -> tuple[int, int]:
+        """``(int(span_y), int(span_x))``."""
+        return (self.height, self.width)
+
+    @property
+    def range_x(self) -> tuple[float, float]:
+        """``left, right``"""
+        return (self.left, self.right)
+
+    @property
+    def range_y(self) -> tuple[float, float]:
+        """``bottom, top``"""
+        return (self.bottom, self.top)
+
+    @property
+    def points(self) -> list[tuple[float, float]]:
+        """Extract four corners of the bounding box."""
+        x0, y0, x1, y1 = self._box
+        return list(itertools.product((x0, x1), (y0, y1)))
+
+    def transform(self, transform: Affine) -> "BoundingBox":
+        """
+        Map bounding box through a linear transform.
+
+        Apply linear transform on four points of the bounding box and compute bounding box of these
+        four points.
+        """
+        pts = [transform * pt for pt in self.points]
+        xx = [x for x, _ in pts]
+        yy = [y for _, y in pts]
+        return BoundingBox(min(xx), min(yy), max(xx), max(yy))
+
+    def map_bounds(self) -> tuple[tuple[float, float], tuple[float, float]]:
+        """
+        Convert to bounds in folium/ipyleaflet style.
+
+        Returns SW, and NE corners in lat/lon order.
+        ``((lat_w, lon_s), (lat_e, lon_n))``.
+        """
+        x0, y0, x1, y1 = self._box
+        return (y0, x0), (y1, x1)
+
+    @staticmethod
+    def from_xy(x: tuple[float, float], y: tuple[float, float]) -> "BoundingBox":
+        """
+        Construct :py:class:`BoundingBox` from x and y ranges.
+
+        :param x: (left, right)
+        :param y: (bottom, top)
+        """
+        x1, x2 = sorted(x)
+        y1, y2 = sorted(y)
+        return BoundingBox(x1, y1, x2, y2)
+
+    @staticmethod
+    def from_points(p1: tuple[float, float], p2: tuple[float, float]) -> "BoundingBox":
+        """
+        Construct :py:class:`BoundingBox` from two points.
+
+        :param p1: (x, y)
+        :param p2: (x, y)
+        """
+        return BoundingBox.from_xy((p1[0], p2[0]), (p1[1], p2[1]))
+
+    @staticmethod
+    def from_transform(shape: tuple[int, int], transform: Affine) -> "BoundingBox":
+        """
+        Construct :py:class:`BoundingBox` from image shape and transform.
+
+        :param shape: image dimensions
+        :param transform: Affine mapping from pixel to world
+        """
+        ny, nx = shape
+        pts = [(0, 0), (nx, 0), (nx, ny), (0, ny)]
+        transform.itransform(pts)
+        xx, yy = list(zip(*pts))
+        return BoundingBox.from_xy((min(xx), max(xx)), (min(yy), max(yy)))
+
+    def round(self) -> "BoundingBox":
+        """
+        Expand bounding box to nearest integer on all sides.
+        """
+        x0, y0, x1, y1 = self._box
+        return BoundingBox(math.floor(x0), math.floor(y0), math.ceil(x1), math.ceil(y1))
+
+
+def bbox_union(bbs: Iterable[BoundingBox]) -> BoundingBox:
+    """
+    Compute union of bounding boxes.
+
+    Given a stream of bounding boxes compute enclosing :py:class:`~odc.geo.geom.BoundingBox`.
+    """
+    # pylint: disable=invalid-name
+    try:
+        bb, *bbs = bbs
+    except ValueError:
+        raise ValueError("Union of empty stream is undefined") from None
+
+    L, B, R, T = bb
+
+    for bb in bbs:
+        l, b, r, t = bb  # noqa: E741
+        L = min(l, L)
+        B = min(b, B)
+        R = max(r, R)
+        T = max(t, T)
+
+    return BoundingBox(L, B, R, T)
+
+
+def bbox_intersection(bbs: Iterable[BoundingBox]) -> BoundingBox:
+    """
+    Compute intersection of bounding boxes.
+
+    Given a stream of bounding boxes compute the overlap :py:class:`~odc.geo.geom.BoundingBox`.
+    """
+    # pylint: disable=invalid-name
+    try:
+        bb, *bbs = bbs
+    except ValueError:
+        raise ValueError("Intersection of empty stream is undefined") from None
+
+    L, B, R, T = bb
+
+    for bb in bbs:
+        l, b, r, t = bb  # noqa: E741
+        L = max(l, L)
+        B = max(b, B)
+        R = min(r, R)
+        T = min(t, T)
+
+    return BoundingBox(L, B, R, T)
diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index a0bc4c4..0d0805f 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -7,8 +7,6 @@
 import numpy as np
 import pandas as pd
 from affine import Affine
-from odc.geo import BoundingBox
-from odc.geo.geom import bbox_intersection, bbox_union
 from xarray import Coordinates, DataArray, Dataset, Index, Variable
 from xarray.core.coordinate_transform import CoordinateTransform
 
@@ -16,6 +14,7 @@
 from xarray.core.indexes import CoordinateTransformIndex, PandasIndex
 from xarray.core.indexing import IndexSelResult, merge_sel_results
 
+from rasterix.odc_compat import BoundingBox, bbox_intersection, bbox_union
 from rasterix.rioxarray_compat import guess_dims
 
 T_Xarray = TypeVar("T_Xarray", "DataArray", "Dataset")
@@ -466,7 +465,6 @@ def bbox(self) -> BoundingBox:
         return BoundingBox.from_transform(
             shape=tuple(self._shape[k] for k in ("y", "x")),
             transform=self.transform() * Affine.translation(-0.5, -0.5),
-            crs=None,
         )
 
     @classmethod

From 5cfd0afca453170f2212c31c22b32565be9c6667 Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Fri, 16 May 2025 10:22:50 -0600
Subject: [PATCH 09/25] Vendor odc math funcs

---
 src/rasterix/odc_compat.py   | 89 ++++++++++++++++++++++++++++++++++++
 src/rasterix/raster_index.py |  4 +-
 2 files changed, 90 insertions(+), 3 deletions(-)

diff --git a/src/rasterix/odc_compat.py b/src/rasterix/odc_compat.py
index c67ecaa..a6b5ce7 100644
--- a/src/rasterix/odc_compat.py
+++ b/src/rasterix/odc_compat.py
@@ -6,6 +6,95 @@
 from affine import Affine
 
 
+def _snap_edge_pos(x0: float, x1: float, res: float, tol: float) -> tuple[float, int]:
+    assert res > 0
+    assert x1 >= x0
+    _x0 = math.floor(maybe_int(x0 / res, tol))
+    _x1 = math.ceil(maybe_int(x1 / res, tol))
+    nx = max(1, _x1 - _x0)
+    return _x0 * res, nx
+
+
+def _snap_edge(x0: float, x1: float, res: float, tol: float) -> tuple[float, int]:
+    assert x1 >= x0
+    if res > 0:
+        return _snap_edge_pos(x0, x1, res, tol)
+    _tx, nx = _snap_edge_pos(x0, x1, -res, tol)
+    tx = _tx + nx * (-res)
+    return tx, nx
+
+
+def snap_grid(
+    x0: float, x1: float, res: float, off_pix: float | None = 0, tol: float = 1e-6
+) -> tuple[float, int]:
+    """
+    Compute grid snapping for single axis.
+
+    :param x0: In point ``x0 <= x1``
+    :param x1: Out point ``x0 <= x1``
+    :param res: Pixel size and direction (can be negative)
+    :param off_pix:
+       Pixel fraction to align to ``x=0``.
+       0 - edge aligned
+       0.5 - center aligned
+       None - don't snap
+
+    :return: ``tx, nx`` that defines 1-d grid, such that ``x0`` and ``x1`` are within edge pixels.
+    """
+    assert (off_pix is None) or (0 <= off_pix < 1)
+    if off_pix is None:
+        if res > 0:
+            nx = math.ceil(maybe_int((x1 - x0) / res, tol))
+            return x0, max(1, nx)
+        nx = math.ceil(maybe_int((x1 - x0) / (-res), tol))
+        return x1, max(nx, 1)
+
+    off = off_pix * abs(res)
+    _tx, nx = _snap_edge(x0 - off, x1 - off, res, tol)
+    return _tx + off, nx
+
+
+def split_float(x: float) -> tuple[float, float]:
+    """
+    Split float number into whole and fractional parts.
+
+    Adding the two numbers back together should result in the original value.
+    Fractional part is always in the ``(-0.5, +0.5)`` interval, and whole part
+    is equivalent to ``round(x)``.
+
+    :param x: floating point number
+    :return: ``whole, fraction``
+    """
+    if not math.isfinite(x):
+        return (x, 0)
+
+    x_part = math.fmod(x, 1.0)
+    x_whole = x - x_part
+    if x_part > 0.5:
+        x_part -= 1
+        x_whole += 1
+    elif x_part < -0.5:
+        x_part += 1
+        x_whole -= 1
+    return (x_whole, x_part)
+
+
+def maybe_int(x: float, tol: float) -> int | float:
+    """
+    Turn almost ints to actual ints.
+
+    pass through other values unmodified.
+    """
+    if not math.isfinite(x):
+        return x
+
+    x_whole, x_part = split_float(x)
+
+    if abs(x_part) < tol:  # almost int
+        return int(x_whole)
+    return x
+
+
 class BoundingBox(Sequence[float]):
     """Bounding box, defining extent in cartesian coordinates."""
 
diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index 0d0805f..ded01b8 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -14,7 +14,7 @@
 from xarray.core.indexes import CoordinateTransformIndex, PandasIndex
 from xarray.core.indexing import IndexSelResult, merge_sel_results
 
-from rasterix.odc_compat import BoundingBox, bbox_intersection, bbox_union
+from rasterix.odc_compat import BoundingBox, bbox_intersection, bbox_union, snap_grid
 from rasterix.rioxarray_compat import guess_dims
 
 T_Xarray = TypeVar("T_Xarray", "DataArray", "Dataset")
@@ -560,8 +560,6 @@ def bbox_to_affine(bbox: BoundingBox, rx, ry) -> Affine:
     # FIXME: translate user-provided `tolerance` to `tol`
     tol: float = 0.01
 
-    from odc.geo.math import snap_grid
-
     offx, nx = snap_grid(bbox.left, bbox.right, rx, 0, tol=tol)
     offy, ny = snap_grid(bbox.bottom, bbox.top, ry, 0, tol=tol)
 

From c5e7e73105541451c4d6f2d3192205985a76d384 Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Fri, 16 May 2025 14:49:11 -0600
Subject: [PATCH 10/25] Add new dim concat test

---
 src/rasterix/raster_index.py |  2 ++
 tests/test_raster_index.py   | 21 ++++++++++++++++++++-
 2 files changed, 22 insertions(+), 1 deletion(-)

diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index ded01b8..64b3d54 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -34,6 +34,8 @@ def assign_index(obj: T_Xarray, *, x_dim: str | None = None, y_dim: str | None =
 
 
 def _assert_transforms_are_compatible(A1: Affine, A2: Affine) -> None:
+    # TODO: offsets could be multiples
+    #       offset should be compatible too
     assert A1.a == A2.a
     assert A1.b == A2.b
     assert A1.d == A2.d
diff --git a/tests/test_raster_index.py b/tests/test_raster_index.py
index 663f1b7..3a54d3f 100644
--- a/tests/test_raster_index.py
+++ b/tests/test_raster_index.py
@@ -49,4 +49,23 @@ def test_combine_nested():
         for transform in transforms
     ]
     datasets = list(map(assign_index, datasets))
-    xr.combine_nested(datasets, concat_dim="x")
+    # xr.combine_nested(datasets, concat_dim="x")
+    xr.concat(datasets, concat_dim="x")
+
+
+def test_concat_new_dim():
+    transforms = [
+        "-50.0 0.5 0.0 0.0 0.0 -0.25",
+        "-50.0 0.5 0.0 0.0 0.0 -0.25",
+    ]
+    crs_attrs = pyproj.CRS.from_epsg(4326).to_cf()
+
+    datasets = [
+        xr.Dataset(
+            {"foo": (("y", "x"), np.ones((4, 2)), {"grid_mapping": "spatial_ref"})},
+            coords={"spatial_ref": ((), 0, crs_attrs | {"GeoTransform": transform})},
+        )
+        for transform in transforms
+    ]
+    datasets = list(map(assign_index, datasets))
+    xr.concat(datasets, dim="time", join="exact")

From 9fddc8488b1f47a9312bb6248a0a0737c86d5319 Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Fri, 16 May 2025 15:01:37 -0600
Subject: [PATCH 11/25] more vendoring

---
 src/rasterix/raster_index.py | 13 +++++--------
 1 file changed, 5 insertions(+), 8 deletions(-)

diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index 64b3d54..5202833 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -1,5 +1,6 @@
 from __future__ import annotations
 
+import math
 import textwrap
 from collections.abc import Hashable, Iterable, Mapping, Sequence
 from typing import Any, Self, TypeVar
@@ -14,7 +15,7 @@
 from xarray.core.indexes import CoordinateTransformIndex, PandasIndex
 from xarray.core.indexing import IndexSelResult, merge_sel_results
 
-from rasterix.odc_compat import BoundingBox, bbox_intersection, bbox_union, snap_grid
+from rasterix.odc_compat import BoundingBox, bbox_intersection, bbox_union, maybe_int, snap_grid
 from rasterix.rioxarray_compat import guess_dims
 
 T_Xarray = TypeVar("T_Xarray", "DataArray", "Dataset")
@@ -537,14 +538,10 @@ def reindex_like(self, other: Self, method=None, tolerance=None) -> dict[Hashabl
 
 
 def get_indexer(off, our_off, start, stop, spacing, tol, size) -> np.ndarray:
-    from math import ceil
+    istart = math.ceil(maybe_int((start - off) / spacing, tol))
 
-    from odc.geo.math import maybe_int
-
-    istart = ceil(maybe_int((start - off) / spacing, tol))
-
-    ours_istart = ceil(maybe_int((start - our_off) / spacing, tol))
-    ours_istop = ceil(maybe_int((stop - our_off) / spacing, tol))
+    ours_istart = math.ceil(maybe_int((start - our_off) / spacing, tol))
+    ours_istop = math.ceil(maybe_int((stop - our_off) / spacing, tol))
 
     idxr = np.concatenate(
         [

From a3691c68598e54685b3c67c24bd80211f5283291 Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Wed, 28 May 2025 22:12:59 -0600
Subject: [PATCH 12/25] fix transforms

---
 src/rasterix/raster_index.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index 5202833..9ea767e 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -454,6 +454,9 @@ def __repr__(self):
 
     def transform(self) -> Affine:
         """Returns Affine transform for top-left corners."""
+        return self.center_transform() * Affine.translation(-0.5, -0.5)
+
+    def center_transform(self) -> Affine:
         if len(self._wrapped_indexes) > 1:
             x = self._wrapped_indexes["x"].axis_transform.affine
             y = self._wrapped_indexes["y"].axis_transform.affine
@@ -461,13 +464,13 @@ def transform(self) -> Affine:
         else:
             index = next(iter(self._wrapped_indexes.values()))
             aff = index.affine
-        return aff * Affine.translation(-0.5, -0.5)
+        return aff
 
     @property
     def bbox(self) -> BoundingBox:
         return BoundingBox.from_transform(
             shape=tuple(self._shape[k] for k in ("y", "x")),
-            transform=self.transform() * Affine.translation(-0.5, -0.5),
+            transform=self.transform(),
         )
 
     @classmethod
@@ -493,15 +496,12 @@ def join(self, other: RasterIndex, how) -> RasterIndex:
                 "Alignment is only supported between RasterIndexes, when both contain compatible transforms."
             )
 
-        # move affines back to top-left corner
-        trans = Affine.translation(-0.5, -0.5)
-        ours = self.transform() * trans
-        theirs = other.transform() * trans
+        ours = self.transform()
+        theirs = other.transform()
         ours, theirs = as_compatible_bboxes(self, other)
 
         if how == "outer":
             new_bbox = bbox_union([ours, theirs])
-            print(new_bbox, ours, theirs)
         elif how == "inner":
             new_bbox = bbox_intersection([ours, theirs])
         else:

From a7a9069243e103da60e59cee360064bf30f27f7c Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Wed, 28 May 2025 22:27:10 -0600
Subject: [PATCH 13/25] Fix test

---
 src/rasterix/raster_index.py | 50 ++++++++++++++++++++----------------
 tests/test_raster_index.py   | 15 ++++++++---
 2 files changed, 39 insertions(+), 26 deletions(-)

diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index 9ea767e..36d3c5f 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -34,13 +34,15 @@ def assign_index(obj: T_Xarray, *, x_dim: str | None = None, y_dim: str | None =
     return obj.assign_coords(coords)
 
 
-def _assert_transforms_are_compatible(A1: Affine, A2: Affine) -> None:
+def _assert_transforms_are_compatible(*affines) -> None:
     # TODO: offsets could be multiples
     #       offset should be compatible too
-    assert A1.a == A2.a
-    assert A1.b == A2.b
-    assert A1.d == A2.d
-    assert A1.e == A2.e
+    A1 = affines[0]
+    for A2 in affines[1:]:
+        assert A1.a == A2.a
+        assert A1.b == A2.b
+        assert A1.d == A2.d
+        assert A1.e == A2.e
 
 
 class AffineTransform(CoordinateTransform):
@@ -482,7 +484,22 @@ def concat(
     ) -> Self:
         if len(indexes) == 1:
             return next(iter(indexes))
-        raise NotImplementedError
+
+        if positions is not None:
+            raise NotImplementedError
+
+        # Note: I am assuming that xarray has calling `align(..., exclude="x" [or 'y'])` already
+        # and checked for equality along "y" [or 'x']
+        new_bbox = bbox_union(as_compatible_bboxes(*indexes))
+        return indexes[0]._new_with_bbox(new_bbox)
+
+    def _new_with_bbox(self, bbox: BoundingBox) -> Self:
+        affine = self.transform()
+        new_affine, Nx, Ny = bbox_to_affine(bbox, rx=affine.a, ry=affine.e)
+        # TODO: set xdim, ydim explicitly
+        new_index = self.from_transform(new_affine, width=Nx, height=Ny)
+        assert new_index.bbox == bbox
+        return new_index
 
     def join(self, other: RasterIndex, how) -> RasterIndex:
         if not isinstance(other, RasterIndex):
@@ -496,10 +513,7 @@ def join(self, other: RasterIndex, how) -> RasterIndex:
                 "Alignment is only supported between RasterIndexes, when both contain compatible transforms."
             )
 
-        ours = self.transform()
-        theirs = other.transform()
         ours, theirs = as_compatible_bboxes(self, other)
-
         if how == "outer":
             new_bbox = bbox_union([ours, theirs])
         elif how == "inner":
@@ -507,13 +521,7 @@ def join(self, other: RasterIndex, how) -> RasterIndex:
         else:
             raise NotImplementedError
 
-        affine = self.transform()
-        new_affine, Nx, Ny = bbox_to_affine(new_bbox, rx=affine.a, ry=affine.e)
-
-        # TODO: set xdim, ydim explicitly
-        new_index = self.from_transform(new_affine, width=Nx, height=Ny)
-        assert new_index.bbox == new_bbox
-        return new_index
+        return self._new_with_bbox(new_bbox)
 
     def reindex_like(self, other: Self, method=None, tolerance=None) -> dict[Hashable, Any]:
         affine = self.transform()
@@ -567,9 +575,7 @@ def bbox_to_affine(bbox: BoundingBox, rx, ry) -> Affine:
     return affine, nx, ny
 
 
-def as_compatible_bboxes(r1: RasterIndex, r2: RasterIndex) -> tuple[BoundingBox, BoundingBox]:
-    r1_transform = r1.transform()
-    r2_transform = r2.transform()
-    _assert_transforms_are_compatible(r1_transform, r2_transform)
-
-    return r1.bbox, r2.bbox
+def as_compatible_bboxes(*indexes: RasterIndex) -> tuple[BoundingBox, ...]:
+    transforms = tuple(i.transform() for i in indexes)
+    _assert_transforms_are_compatible(*transforms)
+    return tuple(i.bbox for i in indexes)
diff --git a/tests/test_raster_index.py b/tests/test_raster_index.py
index 3a54d3f..44572bb 100644
--- a/tests/test_raster_index.py
+++ b/tests/test_raster_index.py
@@ -3,6 +3,7 @@
 import rioxarray  # noqa
 import xarray as xr
 from affine import Affine
+from xarray.testing import assert_identical
 
 from rasterix import RasterIndex, assign_index
 
@@ -36,8 +37,8 @@ def test_sel_slice():
 
 def test_combine_nested():
     transforms = [
-        "-50.0 0.5 0.0 0.0 0.0 -0.25",
-        "-40.0 0.5 0.0 0.0 0.0 -0.25",
+        "-50.0 5 0.0 0.0 0.0 -0.25",
+        "-40.0 5 0.0 0.0 0.0 -0.25",
     ]
     crs_attrs = pyproj.CRS.from_epsg(4326).to_cf()
 
@@ -49,8 +50,14 @@ def test_combine_nested():
         for transform in transforms
     ]
     datasets = list(map(assign_index, datasets))
-    # xr.combine_nested(datasets, concat_dim="x")
-    xr.concat(datasets, concat_dim="x")
+
+    expected = xr.Dataset(
+        {"foo": (("y", "x"), np.ones((4, 2 * len(datasets))), {"grid_mapping": "spatial_ref"})},
+        coords={"spatial_ref": ((), 0, crs_attrs | {"GeoTransform": transforms[0]})},
+    ).pipe(assign_index)
+
+    assert_identical(xr.combine_nested(datasets, concat_dim="x", combine_attrs="override"), expected)
+    assert_identical(xr.concat(datasets, dim="x"), expected)
 
 
 def test_concat_new_dim():

From 919e3535fa71b064a2501f27318aeb1f3c52338e Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Wed, 28 May 2025 22:31:05 -0600
Subject: [PATCH 14/25] dump 3.10

---
 .github/workflows/release.yml | 2 +-
 .github/workflows/test.yml    | 2 +-
 pyproject.toml                | 5 ++---
 3 files changed, 4 insertions(+), 5 deletions(-)

diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
index c847e5e..e84faf9 100644
--- a/.github/workflows/release.yml
+++ b/.github/workflows/release.yml
@@ -18,7 +18,7 @@ jobs:
       - uses: actions/setup-python@v5.2.0
         name: Install Python
         with:
-          python-version: "3.10"
+          python-version: "3.12"
 
       - name: Install PyBuild
         run: |
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index f546fbf..a69f7d0 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -20,7 +20,7 @@ jobs:
 
     strategy:
       matrix:
-        python-version: ["3.10", "3.13"]
+        python-version: ["3.11", "3.13"]
         os: ["ubuntu-latest"]
     runs-on: ${{ matrix.os }}
 
diff --git a/pyproject.toml b/pyproject.toml
index 7ef820d..9eca752 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -14,14 +14,13 @@ name = "rasterix"
 description = "Raster extensions for Xarray"
 license = "Apache-2.0"
 readme = "README.md"
-requires-python = ">=3.10"
+requires-python = ">=3.11"
 keywords = ["xarray"]
 classifiers = [
     "Development Status :: 4 - Beta",
     "Natural Language :: English",
     "Operating System :: OS Independent",
     "Programming Language :: Python",
-    "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
     "Programming Language :: Python :: 3.13",
@@ -75,7 +74,7 @@ dependencies = [
 features = ["test"]
 
 [[tool.hatch.envs.test.matrix]]
-python = ["3.10", "3.13"]
+python = ["3.11", "3.13"]
 
 [tool.hatch.envs.test.scripts]
 run-coverage = "pytest --cov-config=pyproject.toml --cov=pkg --cov-report xml --cov=src --junitxml=junit.xml -o junit_family=legacy"

From 7fad664343db4361fc3f0ad4941d92eac0cee02f Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Tue, 3 Jun 2025 09:02:23 -0600
Subject: [PATCH 15/25] cleanup

---
 src/rasterix/odc_compat.py   | 63 ------------------------------------
 src/rasterix/raster_index.py |  4 +--
 2 files changed, 2 insertions(+), 65 deletions(-)

diff --git a/src/rasterix/odc_compat.py b/src/rasterix/odc_compat.py
index a6b5ce7..c192970 100644
--- a/src/rasterix/odc_compat.py
+++ b/src/rasterix/odc_compat.py
@@ -152,64 +152,11 @@ def __and__(self, other: "BoundingBox") -> "BoundingBox":
     def __or__(self, other: "BoundingBox") -> "BoundingBox":
         return bbox_union([self, other])
 
-    def buffered(self, xbuff: float, ybuff: float | None = None) -> "BoundingBox":
-        """
-        Return a new BoundingBox, buffered in the x and y dimensions.
-
-        :param xbuff: X dimension buffering amount
-        :param ybuff: Y dimension buffering amount
-        :return: new BoundingBox
-        """
-        if ybuff is None:
-            ybuff = xbuff
-
-        return BoundingBox(
-            left=self.left - xbuff,
-            bottom=self.bottom - ybuff,
-            right=self.right + xbuff,
-            top=self.top + ybuff,
-        )
-
-    @property
-    def span_x(self) -> float:
-        """Span of the bounding box along x axis."""
-        return self.right - self.left
-
-    @property
-    def span_y(self) -> float:
-        """Span of the bounding box along y axis."""
-        return self.top - self.bottom
-
-    @property
-    def aspect(self) -> float:
-        """Aspect ratio."""
-        return self.span_x / self.span_y
-
-    @property
-    def width(self) -> int:
-        """``int(span_x)``"""
-        return int(self.right - self.left)
-
-    @property
-    def height(self) -> int:
-        """``int(span_y)``"""
-        return int(self.top - self.bottom)
-
     @property
     def shape(self) -> tuple[int, int]:
         """``(int(span_y), int(span_x))``."""
         return (self.height, self.width)
 
-    @property
-    def range_x(self) -> tuple[float, float]:
-        """``left, right``"""
-        return (self.left, self.right)
-
-    @property
-    def range_y(self) -> tuple[float, float]:
-        """``bottom, top``"""
-        return (self.bottom, self.top)
-
     @property
     def points(self) -> list[tuple[float, float]]:
         """Extract four corners of the bounding box."""
@@ -228,16 +175,6 @@ def transform(self, transform: Affine) -> "BoundingBox":
         yy = [y for _, y in pts]
         return BoundingBox(min(xx), min(yy), max(xx), max(yy))
 
-    def map_bounds(self) -> tuple[tuple[float, float], tuple[float, float]]:
-        """
-        Convert to bounds in folium/ipyleaflet style.
-
-        Returns SW, and NE corners in lat/lon order.
-        ``((lat_w, lon_s), (lat_e, lon_n))``.
-        """
-        x0, y0, x1, y1 = self._box
-        return (y0, x0), (y1, x1)
-
     @staticmethod
     def from_xy(x: tuple[float, float], y: tuple[float, float]) -> "BoundingBox":
         """
diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index 36d3c5f..5084985 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -515,9 +515,9 @@ def join(self, other: RasterIndex, how) -> RasterIndex:
 
         ours, theirs = as_compatible_bboxes(self, other)
         if how == "outer":
-            new_bbox = bbox_union([ours, theirs])
+            new_bbox = ours | theirs
         elif how == "inner":
-            new_bbox = bbox_intersection([ours, theirs])
+            new_bbox = ours & theirs
         else:
             raise NotImplementedError
 

From 5ab308a1265587bddaf1ce5d4fc68e974becfe71 Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Thu, 12 Jun 2025 08:44:58 -0600
Subject: [PATCH 16/25] PR comment

---
 src/rasterix/raster_index.py | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index 5084985..f3767d5 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -502,11 +502,6 @@ def _new_with_bbox(self, bbox: BoundingBox) -> Self:
         return new_index
 
     def join(self, other: RasterIndex, how) -> RasterIndex:
-        if not isinstance(other, RasterIndex):
-            raise ValueError(
-                f"Alignment is only supported between RasterIndexes. Received RasterIndex and {type(other)!r} instead"
-            )
-
         if set(self._wrapped_indexes.keys()) != set(other._wrapped_indexes.keys()):
             # TODO: better error message
             raise ValueError(

From aa58457cf042e887a670931f06305dc62b70a9cc Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Thu, 12 Jun 2025 08:46:25 -0600
Subject: [PATCH 17/25] Add odc license

---
 src/rasterix/odc_compat.py | 204 +++++++++++++++++++++++++++++++++++++
 1 file changed, 204 insertions(+)

diff --git a/src/rasterix/odc_compat.py b/src/rasterix/odc_compat.py
index c192970..60cdf8a 100644
--- a/src/rasterix/odc_compat.py
+++ b/src/rasterix/odc_compat.py
@@ -1,3 +1,207 @@
+# The code below is reproduced from the odc-geo project under terms of their Apache-2.0 license (reproduced below)
+#
+#                                 Apache License
+#                           Version 2.0, January 2004
+#                        http://www.apache.org/licenses/
+#
+#   TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION
+#
+#   1. Definitions.
+#
+#      "License" shall mean the terms and conditions for use, reproduction,
+#      and distribution as defined by Sections 1 through 9 of this document.
+#
+#      "Licensor" shall mean the copyright owner or entity authorized by
+#      the copyright owner that is granting the License.
+#
+#      "Legal Entity" shall mean the union of the acting entity and all
+#      other entities that control, are controlled by, or are under common
+#      control with that entity. For the purposes of this definition,
+#      "control" means (i) the power, direct or indirect, to cause the
+#      direction or management of such entity, whether by contract or
+#      otherwise, or (ii) ownership of fifty percent (50%) or more of the
+#      outstanding shares, or (iii) beneficial ownership of such entity.
+#
+#      "You" (or "Your") shall mean an individual or Legal Entity
+#      exercising permissions granted by this License.
+#
+#      "Source" form shall mean the preferred form for making modifications,
+#      including but not limited to software source code, documentation
+#      source, and configuration files.
+#
+#      "Object" form shall mean any form resulting from mechanical
+#      transformation or translation of a Source form, including but
+#      not limited to compiled object code, generated documentation,
+#      and conversions to other media types.
+#
+#      "Work" shall mean the work of authorship, whether in Source or
+#      Object form, made available under the License, as indicated by a
+#      copyright notice that is included in or attached to the work
+#      (an example is provided in the Appendix below).
+#
+#      "Derivative Works" shall mean any work, whether in Source or Object
+#      form, that is based on (or derived from) the Work and for which the
+#      editorial revisions, annotations, elaborations, or other modifications
+#      represent, as a whole, an original work of authorship. For the purposes
+#      of this License, Derivative Works shall not include works that remain
+#      separable from, or merely link (or bind by name) to the interfaces of,
+#      the Work and Derivative Works thereof.
+#
+#      "Contribution" shall mean any work of authorship, including
+#      the original version of the Work and any modifications or additions
+#      to that Work or Derivative Works thereof, that is intentionally
+#      submitted to Licensor for inclusion in the Work by the copyright owner
+#      or by an individual or Legal Entity authorized to submit on behalf of
+#      the copyright owner. For the purposes of this definition, "submitted"
+#      means any form of electronic, verbal, or written communication sent
+#      to the Licensor or its representatives, including but not limited to
+#      communication on electronic mailing lists, source code control systems,
+#      and issue tracking systems that are managed by, or on behalf of, the
+#      Licensor for the purpose of discussing and improving the Work, but
+#      excluding communication that is conspicuously marked or otherwise
+#      designated in writing by the copyright owner as "Not a Contribution."
+#
+#      "Contributor" shall mean Licensor and any individual or Legal Entity
+#      on behalf of whom a Contribution has been received by Licensor and
+#      subsequently incorporated within the Work.
+#
+#   2. Grant of Copyright License. Subject to the terms and conditions of
+#      this License, each Contributor hereby grants to You a perpetual,
+#      worldwide, non-exclusive, no-charge, royalty-free, irrevocable
+#      copyright license to reproduce, prepare Derivative Works of,
+#      publicly display, publicly perform, sublicense, and distribute the
+#      Work and such Derivative Works in Source or Object form.
+#
+#   3. Grant of Patent License. Subject to the terms and conditions of
+#      this License, each Contributor hereby grants to You a perpetual,
+#      worldwide, non-exclusive, no-charge, royalty-free, irrevocable
+#      (except as stated in this section) patent license to make, have made,
+#      use, offer to sell, sell, import, and otherwise transfer the Work,
+#      where such license applies only to those patent claims licensable
+#      by such Contributor that are necessarily infringed by their
+#      Contribution(s) alone or by combination of their Contribution(s)
+#      with the Work to which such Contribution(s) was submitted. If You
+#      institute patent litigation against any entity (including a
+#      cross-claim or counterclaim in a lawsuit) alleging that the Work
+#      or a Contribution incorporated within the Work constitutes direct
+#      or contributory patent infringement, then any patent licenses
+#      granted to You under this License for that Work shall terminate
+#      as of the date such litigation is filed.
+#
+#   4. Redistribution. You may reproduce and distribute copies of the
+#      Work or Derivative Works thereof in any medium, with or without
+#      modifications, and in Source or Object form, provided that You
+#      meet the following conditions:
+#
+#      (a) You must give any other recipients of the Work or
+#          Derivative Works a copy of this License; and
+#
+#      (b) You must cause any modified files to carry prominent notices
+#          stating that You changed the files; and
+#
+#      (c) You must retain, in the Source form of any Derivative Works
+#          that You distribute, all copyright, patent, trademark, and
+#          attribution notices from the Source form of the Work,
+#          excluding those notices that do not pertain to any part of
+#          the Derivative Works; and
+#
+#      (d) If the Work includes a "NOTICE" text file as part of its
+#          distribution, then any Derivative Works that You distribute must
+#          include a readable copy of the attribution notices contained
+#          within such NOTICE file, excluding those notices that do not
+#          pertain to any part of the Derivative Works, in at least one
+#          of the following places: within a NOTICE text file distributed
+#          as part of the Derivative Works; within the Source form or
+#          documentation, if provided along with the Derivative Works; or,
+#          within a display generated by the Derivative Works, if and
+#          wherever such third-party notices normally appear. The contents
+#          of the NOTICE file are for informational purposes only and
+#          do not modify the License. You may add Your own attribution
+#          notices within Derivative Works that You distribute, alongside
+#          or as an addendum to the NOTICE text from the Work, provided
+#          that such additional attribution notices cannot be construed
+#          as modifying the License.
+#
+#      You may add Your own copyright statement to Your modifications and
+#      may provide additional or different license terms and conditions
+#      for use, reproduction, or distribution of Your modifications, or
+#      for any such Derivative Works as a whole, provided Your use,
+#      reproduction, and distribution of the Work otherwise complies with
+#      the conditions stated in this License.
+#
+#   5. Submission of Contributions. Unless You explicitly state otherwise,
+#      any Contribution intentionally submitted for inclusion in the Work
+#      by You to the Licensor shall be under the terms and conditions of
+#      this License, without any additional terms or conditions.
+#      Notwithstanding the above, nothing herein shall supersede or modify
+#      the terms of any separate license agreement you may have executed
+#      with Licensor regarding such Contributions.
+#
+#   6. Trademarks. This License does not grant permission to use the trade
+#      names, trademarks, service marks, or product names of the Licensor,
+#      except as required for reasonable and customary use in describing the
+#      origin of the Work and reproducing the content of the NOTICE file.
+#
+#   7. Disclaimer of Warranty. Unless required by applicable law or
+#      agreed to in writing, Licensor provides the Work (and each
+#      Contributor provides its Contributions) on an "AS IS" BASIS,
+#      WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
+#      implied, including, without limitation, any warranties or conditions
+#      of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A
+#      PARTICULAR PURPOSE. You are solely responsible for determining the
+#      appropriateness of using or redistributing the Work and assume any
+#      risks associated with Your exercise of permissions under this License.
+#
+#   8. Limitation of Liability. In no event and under no legal theory,
+#      whether in tort (including negligence), contract, or otherwise,
+#      unless required by applicable law (such as deliberate and grossly
+#      negligent acts) or agreed to in writing, shall any Contributor be
+#      liable to You for damages, including any direct, indirect, special,
+#      incidental, or consequential damages of any character arising as a
+#      result of this License or out of the use or inability to use the
+#      Work (including but not limited to damages for loss of goodwill,
+#      work stoppage, computer failure or malfunction, or any and all
+#      other commercial damages or losses), even if such Contributor
+#      has been advised of the possibility of such damages.
+#
+#   9. Accepting Warranty or Additional Liability. While redistributing
+#      the Work or Derivative Works thereof, You may choose to offer,
+#      and charge a fee for, acceptance of support, warranty, indemnity,
+#      or other liability obligations and/or rights consistent with this
+#      License. However, in accepting such obligations, You may act only
+#      on Your own behalf and on Your sole responsibility, not on behalf
+#      of any other Contributor, and only if You agree to indemnify,
+#      defend, and hold each Contributor harmless for any liability
+#      incurred by, or claims asserted against, such Contributor by reason
+#      of your accepting any such warranty or additional liability.
+#
+#   END OF TERMS AND CONDITIONS
+#
+#   APPENDIX: How to apply the Apache License to your work.
+#
+#      To apply the Apache License to your work, attach the following
+#      boilerplate notice, with the fields enclosed by brackets "[]"
+#      replaced with your own identifying information. (Don't include
+#      the brackets!)  The text should be enclosed in the appropriate
+#      comment syntax for the file format. We also recommend that a
+#      file or class name and description of purpose be included on the
+#      same "printed page" as the copyright notice for easier
+#      identification within third-party archives.
+#
+#   Copyright [yyyy] [name of copyright owner]
+#
+#   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.
+
 import itertools
 import math
 from collections.abc import Iterable, Sequence

From 8649541cfe0e1cfe8ddbac7c9176ab022d245e90 Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Thu, 12 Jun 2025 09:15:37 -0600
Subject: [PATCH 18/25] Add more tests

---
 tests/test_raster_index.py | 141 ++++++++++++++++++++++++++++++++++---
 1 file changed, 133 insertions(+), 8 deletions(-)

diff --git a/tests/test_raster_index.py b/tests/test_raster_index.py
index 44572bb..93e28c4 100644
--- a/tests/test_raster_index.py
+++ b/tests/test_raster_index.py
@@ -1,5 +1,6 @@
 import numpy as np
 import pyproj
+import pytest
 import rioxarray  # noqa
 import xarray as xr
 from affine import Affine
@@ -35,11 +36,39 @@ def test_sel_slice():
     assert actual_transform == (transform * Affine.translation(0, 3))
 
 
-def test_combine_nested():
-    transforms = [
-        "-50.0 5 0.0 0.0 0.0 -0.25",
-        "-40.0 5 0.0 0.0 0.0 -0.25",
-    ]
+@pytest.mark.parametrize(
+    "transforms, concat_dim",
+    [
+        (
+            [
+                "-50.0 5 0.0 0.0 0.0 -0.25",
+                "-40.0 5 0.0 0.0 0.0 -0.25",
+                "-30.0 5 0.0 0.0 0.0 -0.25",
+            ],
+            "x",
+        ),
+        (
+            [
+                # decreasing Δy
+                "-40.0 5 0.0 2.0 0.0 -0.5",
+                "-40.0 5 0.0 0.0 0.0 -0.5",
+                "-40.0 5 0.0 -2.0 0.0 -0.5",
+            ],
+            "y",
+        ),
+        (
+            [
+                # increasing Δy
+                "-40.0 5 0.0 -2.0 0.0 0.5",
+                "-40.0 5 0.0 0.0 0.0 0.5",
+                "-40.0 5 0.0 2.0 0.0 0.5",
+            ],
+            "y",
+        ),
+    ],
+)
+def test_concat_and_combine_nested_1D(transforms, concat_dim):
+    """Models two side-by-side tiles in"""
     crs_attrs = pyproj.CRS.from_epsg(4326).to_cf()
 
     datasets = [
@@ -51,16 +80,112 @@ def test_combine_nested():
     ]
     datasets = list(map(assign_index, datasets))
 
+    if concat_dim == "x":
+        new_data = np.ones((4, 2 * len(transforms)))
+    else:
+        new_data = np.ones((4 * len(transforms), 2))
     expected = xr.Dataset(
-        {"foo": (("y", "x"), np.ones((4, 2 * len(datasets))), {"grid_mapping": "spatial_ref"})},
+        {"foo": (("y", "x"), new_data, {"grid_mapping": "spatial_ref"})},
         coords={"spatial_ref": ((), 0, crs_attrs | {"GeoTransform": transforms[0]})},
     ).pipe(assign_index)
 
-    assert_identical(xr.combine_nested(datasets, concat_dim="x", combine_attrs="override"), expected)
-    assert_identical(xr.concat(datasets, dim="x"), expected)
+    for actual in [
+        xr.combine_nested(datasets, concat_dim=concat_dim, combine_attrs="override"),
+        xr.concat(datasets, dim=concat_dim),
+    ]:
+        assert_identical(actual, expected)
+        assert_identical(actual, expected)
+        concat_coord = xr.concat([ds[concat_dim] for ds in datasets], dim=concat_dim)
+        assert_identical(actual[concat_dim], concat_coord)
+
+
+@pytest.mark.parametrize(
+    "transforms, concat_dim",
+    [
+        (
+            [
+                # out-of-order for Y
+                "-40.0 5 0.0 -2.0 0.0 -0.5",
+                "-40.0 5 0.0 0.0 0.0 -0.5",
+                "-40.0 5 0.0 2.0 0.0 -0.5",
+            ],
+            "y",
+        ),
+        (
+            [
+                # incompatible, different origins
+                "-50.0 5 0.0 -2.0 0.0 -0.5",
+                "-40.0 5 0.0 0.0 0.0 -0.5",
+            ],
+            "x",
+        ),
+        (
+            [
+                # incompatible, different Δx
+                "-50.0 2 0.0 0.0 0.0 -0.5",
+                "-40.0 5 0.0 0.0 0.0 -0.5",
+            ],
+            "x",
+        ),
+        (
+            [
+                # incompatible, different Δy
+                "-50.0 2 0.0 0.0 0.0 -0.5",
+                "-40.0 5 0.0 0.0 0.0 -0.25",
+            ],
+            "x",
+        ),
+        (
+            [
+                # exact same transform, makes no sense to concat
+                "-50.0 5 0.0 0.0 0.0 -0.5",
+                "-50.0 5 0.0 0.0 0.0 -0.5",
+            ],
+            "x",
+        ),
+    ],
+)
+def test_concat_errors(transforms, concat_dim):
+    """Models two side-by-side tiles in"""
+    crs_attrs = pyproj.CRS.from_epsg(4326).to_cf()
+
+    datasets = [
+        xr.Dataset(
+            {"foo": (("y", "x"), np.ones((4, 2)), {"grid_mapping": "spatial_ref"})},
+            coords={"spatial_ref": ((), 0, crs_attrs | {"GeoTransform": transform})},
+        )
+        for transform in transforms
+    ]
+    datasets = list(map(assign_index, datasets))
+
+    with pytest.raises(ValueError):
+        xr.combine_nested(datasets, concat_dim=concat_dim, combine_attrs="override")
+    with pytest.raises(ValueError):
+        xr.concat(datasets, dim=concat_dim, combine_attrs="override")
+
+
+def test_concat_different_shape_compatible_transform_error():
+    crs_attrs = pyproj.CRS.from_epsg(4326).to_cf()
+    concat_dim = "x"
+
+    ds1 = xr.Dataset(
+        {"foo": (("y", "x"), np.ones((4, 3)), {"grid_mapping": "spatial_ref"})},
+        coords={"spatial_ref": ((), 0, crs_attrs | {"GeoTransform": "-50.0 5 0.0 0.0 0.0 -0.5"})},
+    )
+    ds2 = xr.Dataset(
+        {"foo": (("y", "x"), np.ones((4, 2)), {"grid_mapping": "spatial_ref"})},
+        coords={"spatial_ref": ((), 0, crs_attrs | {"GeoTransform": "-40.0 5 0.0 0.0 0.0 -0.5"})},
+    )
+
+    datasets = list(map(assign_index, [ds1, ds2]))
+    with pytest.raises(ValueError):
+        xr.combine_nested(datasets, concat_dim=concat_dim, combine_attrs="override")
+    with pytest.raises(ValueError):
+        xr.concat(datasets, dim=concat_dim, combine_attrs="override")
 
 
 def test_concat_new_dim():
+    """models concat along `time` for two tiles with same transform."""
     transforms = [
         "-50.0 0.5 0.0 0.0 0.0 -0.25",
         "-50.0 0.5 0.0 0.0 0.0 -0.25",

From 06af30a6f0a5216074f2fdccab85b713cc441caa Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Thu, 12 Jun 2025 17:12:02 -0600
Subject: [PATCH 19/25] Better errors

---
 src/rasterix/raster_index.py | 30 +++++++++++++++++++++++++-----
 tests/test_raster_index.py   |  3 +--
 2 files changed, 26 insertions(+), 7 deletions(-)

diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index f3767d5..5623485 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -38,11 +38,11 @@ def _assert_transforms_are_compatible(*affines) -> None:
     # TODO: offsets could be multiples
     #       offset should be compatible too
     A1 = affines[0]
-    for A2 in affines[1:]:
-        assert A1.a == A2.a
-        assert A1.b == A2.b
-        assert A1.d == A2.d
-        assert A1.e == A2.e
+    for index, A2 in enumerate(affines[1:]):
+        if A1.a != A2.a or A1.b != A2.b or A1.d != A2.d or A1.e != A2.e:
+            raise ValueError(
+                f"Transform parameters are not compatible for affine 0: {A1}, and affine {index + 1} {A2}"
+            )
 
 
 class AffineTransform(CoordinateTransform):
@@ -573,4 +573,24 @@ def bbox_to_affine(bbox: BoundingBox, rx, ry) -> Affine:
 def as_compatible_bboxes(*indexes: RasterIndex) -> tuple[BoundingBox, ...]:
     transforms = tuple(i.transform() for i in indexes)
     _assert_transforms_are_compatible(*transforms)
+
+    expected_off_x = tuple(t.c + i._shape["x"] * t.a for i, t in zip(indexes, transforms))
+    expected_off_y = tuple(t.f + i._shape["y"] * t.d for i, t in zip(indexes, transforms))
+
+    off_x = tuple(t.c for t in transforms)
+    off_y = tuple(t.c for t in transforms)
+
+    if all(o == off_x[0] for o in off_x[1:]) and all(o == off_y[0] for o in off_y[1:]):
+        raise ValueError("Attempting to concatenate arrays with same transform along X or Y.")
+
+    if any(a != b for a, b in zip(off_x, expected_off_x)):
+        raise ValueError(
+            f"X offsets are incompatible. Provided offsets {off_x}, expected offsets: {expected_off_x}"
+        )
+
+    if any(a != b for a, b in zip(off_y, expected_off_y)):
+        raise ValueError(
+            f"Y offsets are incompatible. Provided offsets {off_y}, expected offsets: {expected_off_y}"
+        )
+
     return tuple(i.bbox for i in indexes)
diff --git a/tests/test_raster_index.py b/tests/test_raster_index.py
index 93e28c4..90f2acc 100644
--- a/tests/test_raster_index.py
+++ b/tests/test_raster_index.py
@@ -68,7 +68,7 @@ def test_sel_slice():
     ],
 )
 def test_concat_and_combine_nested_1D(transforms, concat_dim):
-    """Models two side-by-side tiles in"""
+    """Models two side-by-side tiles"""
     crs_attrs = pyproj.CRS.from_epsg(4326).to_cf()
 
     datasets = [
@@ -146,7 +146,6 @@ def test_concat_and_combine_nested_1D(transforms, concat_dim):
     ],
 )
 def test_concat_errors(transforms, concat_dim):
-    """Models two side-by-side tiles in"""
     crs_attrs = pyproj.CRS.from_epsg(4326).to_cf()
 
     datasets = [

From 2320fed018212155c77b37da5ab34c3a52ba8b7f Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Thu, 12 Jun 2025 17:16:24 -0600
Subject: [PATCH 20/25] Add 2D combine nested test

---
 src/rasterix/raster_index.py | 45 ++++++++++++++++++++++--------------
 tests/test_raster_index.py   | 26 +++++++++++++++++++++
 2 files changed, 54 insertions(+), 17 deletions(-)

diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index 5623485..2777a69 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -490,7 +490,7 @@ def concat(
 
         # Note: I am assuming that xarray has calling `align(..., exclude="x" [or 'y'])` already
         # and checked for equality along "y" [or 'x']
-        new_bbox = bbox_union(as_compatible_bboxes(*indexes))
+        new_bbox = bbox_union(as_compatible_bboxes(*indexes, concat_dim=dim))
         return indexes[0]._new_with_bbox(new_bbox)
 
     def _new_with_bbox(self, bbox: BoundingBox) -> Self:
@@ -508,13 +508,13 @@ def join(self, other: RasterIndex, how) -> RasterIndex:
                 "Alignment is only supported between RasterIndexes, when both contain compatible transforms."
             )
 
-        ours, theirs = as_compatible_bboxes(self, other)
+        ours, theirs = as_compatible_bboxes(self, other, concat_dim=None)
         if how == "outer":
             new_bbox = ours | theirs
         elif how == "inner":
             new_bbox = ours & theirs
         else:
-            raise NotImplementedError
+            raise NotImplementedError(f"{how=!r} not implemented yet for RasterIndex.")
 
         return self._new_with_bbox(new_bbox)
 
@@ -570,27 +570,38 @@ def bbox_to_affine(bbox: BoundingBox, rx, ry) -> Affine:
     return affine, nx, ny
 
 
-def as_compatible_bboxes(*indexes: RasterIndex) -> tuple[BoundingBox, ...]:
+def as_compatible_bboxes(*indexes: RasterIndex, concat_dim: Hashable | None) -> tuple[BoundingBox, ...]:
     transforms = tuple(i.transform() for i in indexes)
     _assert_transforms_are_compatible(*transforms)
 
-    expected_off_x = tuple(t.c + i._shape["x"] * t.a for i, t in zip(indexes, transforms))
-    expected_off_y = tuple(t.f + i._shape["y"] * t.d for i, t in zip(indexes, transforms))
+    expected_off_x = (transforms[0].c,) + tuple(
+        t.c + i._shape["x"] * t.a for i, t in zip(indexes[:-1], transforms[:-1])
+    )
+    expected_off_y = (transforms[0].f,) + tuple(
+        t.f + i._shape["y"] * t.e for i, t in zip(indexes[:-1], transforms[:-1])
+    )
 
     off_x = tuple(t.c for t in transforms)
-    off_y = tuple(t.c for t in transforms)
+    off_y = tuple(t.f for t in transforms)
 
-    if all(o == off_x[0] for o in off_x[1:]) and all(o == off_y[0] for o in off_y[1:]):
-        raise ValueError("Attempting to concatenate arrays with same transform along X or Y.")
+    if concat_dim is not None:
+        if all(o == off_x[0] for o in off_x[1:]) and all(o == off_y[0] for o in off_y[1:]):
+            raise ValueError("Attempting to concatenate arrays with same transform along X or Y.")
 
-    if any(a != b for a, b in zip(off_x, expected_off_x)):
-        raise ValueError(
-            f"X offsets are incompatible. Provided offsets {off_x}, expected offsets: {expected_off_x}"
-        )
+    if concat_dim == "x":
+        if any(off_y[0] != o for o in off_y[1:]):
+            raise ValueError("offsets must be identical in X when concatenating along Y")
+        if any(a != b for a, b in zip(off_x, expected_off_x)):
+            raise ValueError(
+                f"X offsets are incompatible. Provided offsets {off_x}, expected offsets: {expected_off_x}"
+            )
+    elif concat_dim == "y":
+        if any(off_x[0] != o for o in off_x[1:]):
+            raise ValueError("offsets must be identical in X when concatenating along Y")
 
-    if any(a != b for a, b in zip(off_y, expected_off_y)):
-        raise ValueError(
-            f"Y offsets are incompatible. Provided offsets {off_y}, expected offsets: {expected_off_y}"
-        )
+        if any(a != b for a, b in zip(off_y, expected_off_y)):
+            raise ValueError(
+                f"Y offsets are incompatible. Provided offsets {off_y}, expected offsets: {expected_off_y}"
+            )
 
     return tuple(i.bbox for i in indexes)
diff --git a/tests/test_raster_index.py b/tests/test_raster_index.py
index 90f2acc..d3cb1a0 100644
--- a/tests/test_raster_index.py
+++ b/tests/test_raster_index.py
@@ -200,3 +200,29 @@ def test_concat_new_dim():
     ]
     datasets = list(map(assign_index, datasets))
     xr.concat(datasets, dim="time", join="exact")
+
+
+def test_combine_nested_2d():
+    """models 2d tiling"""
+    transforms = [
+        # row 1
+        "-50.0 5 0.0 0.0 0.0 -0.25",
+        "-40.0 5 0.0 0.0 0.0 -0.25",
+        "-30.0 5 0.0 0.0 0.0 -0.25",
+        # row 2
+        "-50.0 5 0.0 -1 0.0 -0.25",
+        "-40.0 5 0.0 -1 0.0 -0.25",
+        "-30.0 5 0.0 -1 0.0 -0.25",
+    ]
+    crs_attrs = pyproj.CRS.from_epsg(4326).to_cf()
+
+    datasets = [
+        xr.Dataset(
+            {"foo": (("y", "x"), np.ones((4, 2)), {"grid_mapping": "spatial_ref"})},
+            coords={"spatial_ref": ((), 0, crs_attrs | {"GeoTransform": transform})},
+        )
+        for transform in transforms
+    ]
+    datasets = list(map(assign_index, datasets))
+    datasets = [datasets[:3], datasets[3:]]
+    xr.combine_nested(datasets, concat_dim=["y", "x"])

From b2cbf896748a4493df3e6effac46c5fade704fb9 Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Thu, 12 Jun 2025 22:40:22 -0600
Subject: [PATCH 21/25] cleanup

---
 src/rasterix/raster_index.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index 2777a69..d25a902 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -35,8 +35,6 @@ def assign_index(obj: T_Xarray, *, x_dim: str | None = None, y_dim: str | None =
 
 
 def _assert_transforms_are_compatible(*affines) -> None:
-    # TODO: offsets could be multiples
-    #       offset should be compatible too
     A1 = affines[0]
     for index, A2 in enumerate(affines[1:]):
         if A1.a != A2.a or A1.b != A2.b or A1.d != A2.d or A1.e != A2.e:

From e1d3623cad815b92b7db834096fd9d452c3a1aaf Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Thu, 12 Jun 2025 22:54:41 -0600
Subject: [PATCH 22/25] Cleanup

---
 tests/test_raster_index.py | 55 ++++++++++----------------------------
 1 file changed, 14 insertions(+), 41 deletions(-)

diff --git a/tests/test_raster_index.py b/tests/test_raster_index.py
index d3cb1a0..06a6006 100644
--- a/tests/test_raster_index.py
+++ b/tests/test_raster_index.py
@@ -8,6 +8,15 @@
 
 from rasterix import RasterIndex, assign_index
 
+CRS_ATTRS = pyproj.CRS.from_epsg(4326).to_cf()
+
+
+def dataset_from_transform(transform: str) -> xr.Dataset:
+    return xr.Dataset(
+        {"foo": (("y", "x"), np.ones((4, 2)), {"grid_mapping": "spatial_ref"})},
+        coords={"spatial_ref": ((), 0, CRS_ATTRS | {"GeoTransform": transform})},
+    ).pipe(assign_index)
+
 
 def test_rectilinear():
     source = "/vsicurl/https://noaadata.apps.nsidc.org/NOAA/G02135/south/daily/geotiff/2024/01_Jan/S_20240101_concentration_v3.0.tif"
@@ -69,16 +78,7 @@ def test_sel_slice():
 )
 def test_concat_and_combine_nested_1D(transforms, concat_dim):
     """Models two side-by-side tiles"""
-    crs_attrs = pyproj.CRS.from_epsg(4326).to_cf()
-
-    datasets = [
-        xr.Dataset(
-            {"foo": (("y", "x"), np.ones((4, 2)), {"grid_mapping": "spatial_ref"})},
-            coords={"spatial_ref": ((), 0, crs_attrs | {"GeoTransform": transform})},
-        )
-        for transform in transforms
-    ]
-    datasets = list(map(assign_index, datasets))
+    datasets = list(map(dataset_from_transform, transforms))
 
     if concat_dim == "x":
         new_data = np.ones((4, 2 * len(transforms)))
@@ -86,7 +86,7 @@ def test_concat_and_combine_nested_1D(transforms, concat_dim):
         new_data = np.ones((4 * len(transforms), 2))
     expected = xr.Dataset(
         {"foo": (("y", "x"), new_data, {"grid_mapping": "spatial_ref"})},
-        coords={"spatial_ref": ((), 0, crs_attrs | {"GeoTransform": transforms[0]})},
+        coords={"spatial_ref": ((), 0, CRS_ATTRS | {"GeoTransform": transforms[0]})},
     ).pipe(assign_index)
 
     for actual in [
@@ -146,17 +146,7 @@ def test_concat_and_combine_nested_1D(transforms, concat_dim):
     ],
 )
 def test_concat_errors(transforms, concat_dim):
-    crs_attrs = pyproj.CRS.from_epsg(4326).to_cf()
-
-    datasets = [
-        xr.Dataset(
-            {"foo": (("y", "x"), np.ones((4, 2)), {"grid_mapping": "spatial_ref"})},
-            coords={"spatial_ref": ((), 0, crs_attrs | {"GeoTransform": transform})},
-        )
-        for transform in transforms
-    ]
-    datasets = list(map(assign_index, datasets))
-
+    datasets = list(map(dataset_from_transform, transforms))
     with pytest.raises(ValueError):
         xr.combine_nested(datasets, concat_dim=concat_dim, combine_attrs="override")
     with pytest.raises(ValueError):
@@ -189,16 +179,7 @@ def test_concat_new_dim():
         "-50.0 0.5 0.0 0.0 0.0 -0.25",
         "-50.0 0.5 0.0 0.0 0.0 -0.25",
     ]
-    crs_attrs = pyproj.CRS.from_epsg(4326).to_cf()
-
-    datasets = [
-        xr.Dataset(
-            {"foo": (("y", "x"), np.ones((4, 2)), {"grid_mapping": "spatial_ref"})},
-            coords={"spatial_ref": ((), 0, crs_attrs | {"GeoTransform": transform})},
-        )
-        for transform in transforms
-    ]
-    datasets = list(map(assign_index, datasets))
+    datasets = list(map(dataset_from_transform, transforms))
     xr.concat(datasets, dim="time", join="exact")
 
 
@@ -214,15 +195,7 @@ def test_combine_nested_2d():
         "-40.0 5 0.0 -1 0.0 -0.25",
         "-30.0 5 0.0 -1 0.0 -0.25",
     ]
-    crs_attrs = pyproj.CRS.from_epsg(4326).to_cf()
 
-    datasets = [
-        xr.Dataset(
-            {"foo": (("y", "x"), np.ones((4, 2)), {"grid_mapping": "spatial_ref"})},
-            coords={"spatial_ref": ((), 0, crs_attrs | {"GeoTransform": transform})},
-        )
-        for transform in transforms
-    ]
-    datasets = list(map(assign_index, datasets))
+    datasets = list(map(dataset_from_transform, transforms))
     datasets = [datasets[:3], datasets[3:]]
     xr.combine_nested(datasets, concat_dim=["y", "x"])

From 1fdb66e9d627c3aade42c8123f5e84458b830a8d Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Thu, 12 Jun 2025 23:11:07 -0600
Subject: [PATCH 23/25] Add align test

---
 src/rasterix/odc_compat.py   |  2 ++
 src/rasterix/raster_index.py |  2 +-
 tests/test_raster_index.py   | 19 +++++++++++++++++++
 3 files changed, 22 insertions(+), 1 deletion(-)

diff --git a/src/rasterix/odc_compat.py b/src/rasterix/odc_compat.py
index 60cdf8a..53c8552 100644
--- a/src/rasterix/odc_compat.py
+++ b/src/rasterix/odc_compat.py
@@ -255,6 +255,8 @@ def snap_grid(
 
     off = off_pix * abs(res)
     _tx, nx = _snap_edge(x0 - off, x1 - off, res, tol)
+    if x1 == x0:
+        nx = 0
     return _tx + off, nx
 
 
diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index d25a902..e89bde7 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -518,7 +518,7 @@ def join(self, other: RasterIndex, how) -> RasterIndex:
 
     def reindex_like(self, other: Self, method=None, tolerance=None) -> dict[Hashable, Any]:
         affine = self.transform()
-        ours, theirs = as_compatible_bboxes(self, other)
+        ours, theirs = as_compatible_bboxes(self, other, concat_dim=None)
         inter = bbox_intersection([ours, theirs])
         dx = affine.a
         dy = affine.e
diff --git a/tests/test_raster_index.py b/tests/test_raster_index.py
index 06a6006..e522974 100644
--- a/tests/test_raster_index.py
+++ b/tests/test_raster_index.py
@@ -199,3 +199,22 @@ def test_combine_nested_2d():
     datasets = list(map(dataset_from_transform, transforms))
     datasets = [datasets[:3], datasets[3:]]
     xr.combine_nested(datasets, concat_dim=["y", "x"])
+
+
+def test_align():
+    transforms = [
+        "-50.0 5 0.0 0.0 0.0 -0.25",
+        "-40.0 5 0.0 0.0 0.0 -0.25",
+    ]
+    expected_affine = Affine(5, 0, -50, 0, -0.25, 0)
+    datasets = list(map(dataset_from_transform, transforms))
+
+    aligned = xr.align(*datasets, join="outer")
+    assert all(a.sizes == {"x": 4, "y": 4} for a in aligned)
+    assert all(a.xindexes["x"].transform() == expected_affine for a in aligned)
+
+    aligned = xr.align(*datasets, join="inner")
+    assert all(a.sizes["x"] == 0 for a in aligned)
+
+    with pytest.raises(xr.AlignmentError):
+        aligned = xr.align(*datasets, join="exact")

From b468b659cd171e3b6181e16f4cca0124e8db1fe8 Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Fri, 13 Jun 2025 08:43:47 -0600
Subject: [PATCH 24/25] Switch to AssertionError

---
 tests/test_raster_index.py | 27 ++++++++++++++++++---------
 1 file changed, 18 insertions(+), 9 deletions(-)

diff --git a/tests/test_raster_index.py b/tests/test_raster_index.py
index e522974..8f2c398 100644
--- a/tests/test_raster_index.py
+++ b/tests/test_raster_index.py
@@ -111,14 +111,6 @@ def test_concat_and_combine_nested_1D(transforms, concat_dim):
             ],
             "y",
         ),
-        (
-            [
-                # incompatible, different origins
-                "-50.0 5 0.0 -2.0 0.0 -0.5",
-                "-40.0 5 0.0 0.0 0.0 -0.5",
-            ],
-            "x",
-        ),
         (
             [
                 # incompatible, different Δx
@@ -137,7 +129,7 @@ def test_concat_and_combine_nested_1D(transforms, concat_dim):
         ),
         (
             [
-                # exact same transform, makes no sense to concat
+                # exact same transform, makes no sense to concat in X or Y
                 "-50.0 5 0.0 0.0 0.0 -0.5",
                 "-50.0 5 0.0 0.0 0.0 -0.5",
             ],
@@ -153,6 +145,23 @@ def test_concat_errors(transforms, concat_dim):
         xr.concat(datasets, dim=concat_dim, combine_attrs="override")
 
 
+def test_concat_error_alignment():
+    transforms, concat_dim = (
+        [
+            # incompatible, different origins
+            "-50.0 5 0.0 -2.0 0.0 -0.5",
+            "-40.0 5 0.0 0.0 0.0 -0.5",
+        ],
+        "x",
+    )
+
+    datasets = list(map(dataset_from_transform, transforms))
+    with pytest.raises(AssertionError):
+        xr.combine_nested(datasets, concat_dim=concat_dim, combine_attrs="override")
+    with pytest.raises(AssertionError):
+        xr.concat(datasets, dim=concat_dim, combine_attrs="override")
+
+
 def test_concat_different_shape_compatible_transform_error():
     crs_attrs = pyproj.CRS.from_epsg(4326).to_cf()
     concat_dim = "x"

From 8b9b155c48a701e801dafcd8232658bf63d08576 Mon Sep 17 00:00:00 2001
From: Deepak Cherian <deepak@cherian.net>
Date: Fri, 13 Jun 2025 08:45:13 -0600
Subject: [PATCH 25/25] Add combine_by_coords skipped test

---
 tests/test_raster_index.py | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/tests/test_raster_index.py b/tests/test_raster_index.py
index 8f2c398..1bfe331 100644
--- a/tests/test_raster_index.py
+++ b/tests/test_raster_index.py
@@ -210,6 +210,24 @@ def test_combine_nested_2d():
     xr.combine_nested(datasets, concat_dim=["y", "x"])
 
 
+@pytest.mark.skip(reason="xarray converts to PandasIndex")
+def test_combine_by_coords():
+    """models 2d tiling"""
+    transforms = [
+        # row 1
+        "-50.0 5 0.0 0.0 0.0 -0.25",
+        "-40.0 5 0.0 0.0 0.0 -0.25",
+        "-30.0 5 0.0 0.0 0.0 -0.25",
+        # row 2
+        "-50.0 5 0.0 -1 0.0 -0.25",
+        "-40.0 5 0.0 -1 0.0 -0.25",
+        "-30.0 5 0.0 -1 0.0 -0.25",
+    ]
+
+    datasets = list(map(dataset_from_transform, transforms))
+    xr.combine_by_coords(datasets)
+
+
 def test_align():
     transforms = [
         "-50.0 5 0.0 0.0 0.0 -0.25",