diff --git a/src/rasterix/raster_index.py b/src/rasterix/raster_index.py
index 3d263ed..31b2e73 100644
--- a/src/rasterix/raster_index.py
+++ b/src/rasterix/raster_index.py
@@ -895,7 +895,7 @@ def concat(
 
     def _new_with_bbox(self, bbox: BoundingBox) -> RasterIndex:
         affine = self.transform()
-        new_affine, Nx, Ny = bbox_to_affine(bbox, rx=affine.a, ry=affine.e)
+        new_affine, Nx, Ny = bbox_to_affine(bbox, affine)
         # TODO: set xdim, ydim explicitly
         new_index = self.from_transform(new_affine, width=Nx, height=Ny)
         assert new_index.bbox == bbox
@@ -984,18 +984,28 @@ def get_indexer(off, our_off, start, stop, spacing, tol, size) -> np.ndarray:
     return idxr
 
 
-def bbox_to_affine(bbox: BoundingBox, rx, ry) -> tuple[Affine, int, int]:
+def bbox_to_affine(bbox: BoundingBox, affine: Affine) -> tuple[Affine, int, int]:
     # 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
 
-    offx, nx = snap_grid(bbox.left, bbox.right, rx, 0, tol=tol)
-    offy, ny = snap_grid(bbox.bottom, bbox.top, ry, 0, tol=tol)
+    rx = affine.a
+    ry = affine.e
 
-    affine = Affine.translation(offx, offy) * Affine.scale(rx, ry)
+    # Calculate the pixel fraction offset from the affine's grid alignment.
+    # This ensures the new grid aligns with the original grid rather than snapping
+    # to a global grid at x=0, y=0.
+    # See https://github.com/xarray-contrib/rasterix/issues/67
+    off_pix_x = (affine.c / abs(rx)) % 1.0
+    off_pix_y = (affine.f / abs(ry)) % 1.0
 
-    return affine, nx, ny
+    offx, nx = snap_grid(bbox.left, bbox.right, rx, off_pix_x, tol=tol)
+    offy, ny = snap_grid(bbox.bottom, bbox.top, ry, off_pix_y, tol=tol)
+
+    new_affine = Affine.translation(offx, offy) * Affine.scale(rx, ry)
+
+    return new_affine, nx, ny
 
 
 def as_compatible_bboxes(*indexes: RasterIndex, concat_dim: Hashable | None) -> tuple[BoundingBox, ...]:
diff --git a/tests/test_raster_index.py b/tests/test_raster_index.py
index 03e8c05..ab48aa4 100644
--- a/tests/test_raster_index.py
+++ b/tests/test_raster_index.py
@@ -408,6 +408,39 @@ def test_align():
         aligned = xr.align(*datasets, join="exact")
 
 
+def test_join_partially_overlapping_indexes():
+    """Test joining two partially overlapping indexes with non-grid-aligned coordinates.
+
+    Regression test for https://github.com/xarray-contrib/rasterix/issues/67
+
+    This tests the case where the transform coordinates are not aligned to a global
+    grid at x=0, y=0 (e.g., 399955.0 / 10.0 = 39995.5, which has a 0.5 pixel offset).
+    """
+    transform1 = Affine.from_gdal(399955.0, 10.0, 0.0, 5500025.0, 0.0, -10.0)
+    transform2 = Affine.from_gdal(399955.0, 10.0, 0.0, 5497625.0, 0.0, -10.0)
+
+    index1 = RasterIndex.from_transform(transform1, width=256, height=256)
+    index2 = RasterIndex.from_transform(transform2, width=256, height=256)
+
+    # Inner join should produce a 256x16 pixel index (the overlap region)
+    result = index1.join(index2, how="inner")
+
+    assert result.xy_shape == (256, 16)
+    assert result.bbox.left == 399955.0
+    assert result.bbox.right == 402515.0
+    assert result.bbox.bottom == 5497465.0
+    assert result.bbox.top == 5497625.0
+
+    # Outer join should produce a 256x496 pixel index (combined coverage)
+    result_outer = index1.join(index2, how="outer")
+
+    assert result_outer.xy_shape == (256, 496)  # 256 + 256 - 16 = 496
+    assert result_outer.bbox.left == 399955.0
+    assert result_outer.bbox.right == 402515.0
+    assert result_outer.bbox.bottom == 5495065.0
+    assert result_outer.bbox.top == 5500025.0
+
+
 def test_repr_inline() -> None:
     index1 = RasterIndex.from_transform(Affine.identity(), width=12, height=10)
     ds1 = xr.Dataset(coords=xr.Coordinates.from_xindex(index1))