Alignment

See also

See the Xarray tutorial on Alignment.

Simple forms of alignment are supported.

Consider this non-overlapping pair of tiles with two transform:

transforms = [
    "-50.0 5 0.0 0.0 0.0 -0.25",
    "-40.0 5 0.0 0.0 0.0 -0.25",
]

Hide code cell source

 
import pyproj
import numpy as np
import xarray as xr
import rasterix

dsets = [
    (i+1) * xr.Dataset(
        {"foo": (("y", "x"), np.ones((4, 2)), {"grid_mapping": "spatial_ref"})},
        coords={
            "spatial_ref": (
                (),
                0,
                pyproj.CRS.from_epsg(4326).to_cf() | {"GeoTransform": transform},
            )
        },
    )
    for i, transform in enumerate(transforms)
]
datasets = list(map(rasterix.assign_index, dsets))
datasets[0]

<xarray.Dataset> Size: 120B
Dimensions:      (y: 4, x: 2)
Coordinates:
  * y            (y) float64 32B -0.125 -0.375 -0.625 -0.875
  * x            (x) float64 16B -47.5 -42.5
    spatial_ref  int64 8B 0
Data variables:
    foo          (y, x) float64 64B 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
Indexes:
  ┌ x        RasterIndex (crs=None)
  └ y
datasets[1]

<xarray.Dataset> Size: 120B
Dimensions:      (y: 4, x: 2)
Coordinates:
  * y            (y) float64 32B -0.125 -0.375 -0.625 -0.875
  * x            (x) float64 16B -37.5 -32.5
    spatial_ref  int64 8B 0
Data variables:
    foo          (y, x) float64 64B 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
Indexes:
  ┌ x        RasterIndex (crs=None)
  └ y

join="outer"

For outer joins, we reindex so that both datasets are expanded to the union of the input bounding boxes:

outer = xr.align(*datasets, join="outer")
outer[0]

<xarray.Dataset> Size: 200B
Dimensions:      (x: 4, y: 4)
Coordinates:
  * x            (x) float64 32B -47.5 -42.5 -37.5 -32.5
  * y            (y) float64 32B -0.125 -0.375 -0.625 -0.875
    spatial_ref  int64 8B 0
Data variables:
    foo          (y, x) float64 128B 1.0 1.0 nan nan 1.0 ... nan 1.0 1.0 nan nan
Indexes:
  ┌ x        RasterIndex (crs=None)
  └ y
outer[1]

<xarray.Dataset> Size: 200B
Dimensions:      (x: 4, y: 4)
Coordinates:
  * x            (x) float64 32B -47.5 -42.5 -37.5 -32.5
  * y            (y) float64 32B -0.125 -0.375 -0.625 -0.875
    spatial_ref  int64 8B 0
Data variables:
    foo          (y, x) float64 128B nan nan 2.0 2.0 nan ... 2.0 nan nan 2.0 2.0
Indexes:
  ┌ x        RasterIndex (crs=None)
  └ y

join="inner"

For outer joins, we reindex so that both datasets are restricted to the intersection of the input bounding boxes:

There is no overlap in x for this set of non-overlapping tiles

inner = xr.align(*datasets, join="inner")
inner[0]

<xarray.Dataset> Size: 40B
Dimensions:      (x: 0, y: 4)
Coordinates:
  * x            (x) float64 0B 
  * y            (y) float64 32B -0.125 -0.375 -0.625 -0.875
    spatial_ref  int64 8B 0
Data variables:
    foo          (y, x) float64 0B 
Indexes:
  ┌ x        RasterIndex (crs=None)
  └ y

join="exact"

For exact joins, we compare transforms. Since these are not identical, we see an AlignmentError.

xr.align(*datasets, join="exact")

AlignmentError: cannot align objects with join='exact' where index/labels/sizes are not equal along these coordinates (dimensions): 'x' ('x',), 'y' ('y',)