Tolerance

Raster tiles from different sources may have tiny floating-point differences in their transform parameters (resolution, origin), causing alignment or concatenation to fail even though the tiles are intended to be compatible.

Rasterix applies a small default relative tolerance (transform_rtol=1e-12) when comparing transform parameters. Use set_options() to adjust this. Both transform_rtol (relative) and transform_atol (absolute) are supported, matching math.isclose() semantics.

import rasterix

rasterix.get_options()

{'transform_rtol': 1e-12, 'transform_atol': 0.0}

Combining

Consider two tiles that should concatenate along X, but have a tiny difference in their X resolution (926.625433054... vs 926.625433055...):

Note

These transforms come from a real dataset MODIS/Terra Vegetation Indices Monthly L3 Global 1km SIN Grid V061!

import pyproj
import numpy as np

transforms = [
    "-8895604.157333 926.6254330549995 0.0 3335851.559 0.0 -926.6254330558334",
    "-7783653.637667 926.6254330558338 0.0 3335851.559 0.0 -926.6254330558334",
]

dsets = [
    xr.Dataset(
        {"temp": (("y", "x"), np.ones((10, 1200)), {"grid_mapping": "spatial_ref"})},
        coords={
            "spatial_ref": (
                (),
                0,
                pyproj.CRS.from_epsg(3857).to_cf() | {"GeoTransform": transform},
            )
        },
    )
    for transform in transforms
]
dsets = list(map(rasterix.assign_index, dsets))

The relative difference here is ~9e-13, which is handled by the default tolerance:

xr.concat(dsets, dim="x")

<xarray.Dataset> Size: 211kB
Dimensions:      (y: 10, x: 2400)
Coordinates:
  * y            (y) float64 80B 3.335e+06 3.334e+06 ... 3.328e+06 3.327e+06
  * x            (x) float64 19kB -8.895e+06 -8.894e+06 ... -6.672e+06
    spatial_ref  int64 8B 0
Data variables:
    temp         (y, x) float64 192kB 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0
Indexes:
  ┌ x        RasterIndex (crs=None)
  └ y

Larger differences

With larger differences (~1e-8 relative), the default tolerance is not enough:

transforms_noisy = [
    "0.0 1.0 0.0 1.0 0.0 -1.0",
    "10.0 1.00000001 0.0 1.0 0.0 -1.0",
]

dsets_noisy = [
    xr.Dataset(
        {"temp": (("y", "x"), np.ones((10, 10)), {"grid_mapping": "spatial_ref"})},
        coords={
            "spatial_ref": (
                (),
                0,
                pyproj.CRS.from_epsg(4326).to_cf() | {"GeoTransform": transform},
            )
        },
    )
    for transform in transforms_noisy
]
dsets_noisy = list(map(rasterix.assign_index, dsets_noisy))

xr.concat(dsets_noisy, dim="x")

ValueError: Transform parameters are not compatible for affine 0: | 1.00, 0.00, 0.00|
| 0.00,-1.00, 1.00|
| 0.00, 0.00, 1.00|, and affine 1 | 1.00, 0.00, 10.00|
| 0.00,-1.00, 1.00|
| 0.00, 0.00, 1.00|

Increase the tolerance using set_options() as a context manager:

with rasterix.set_options(transform_rtol=1e-7):
    result = xr.concat(dsets_noisy, dim="x")
result

<xarray.Dataset> Size: 2kB
Dimensions:      (y: 10, x: 20)
Coordinates:
  * y            (y) float64 80B 0.5 -0.5 -1.5 -2.5 -3.5 ... -5.5 -6.5 -7.5 -8.5
  * x            (x) float64 160B 0.5 1.5 2.5 3.5 4.5 ... 16.5 17.5 18.5 19.5
    spatial_ref  int64 8B 0
Data variables:
    temp         (y, x) float64 2kB 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0
Indexes:
  ┌ x        RasterIndex (crs=None)
  └ y

Alignment

The same tolerance applies to xarray.align():

with rasterix.set_options(transform_rtol=1e-7):
    aligned = xr.align(*dsets_noisy, join="outer")
aligned[0]

<xarray.Dataset> Size: 2kB
Dimensions:      (x: 20, y: 10)
Coordinates:
  * x            (x) float64 160B 0.5 1.5 2.5 3.5 4.5 ... 16.5 17.5 18.5 19.5
  * y            (y) float64 80B 0.5 -0.5 -1.5 -2.5 -3.5 ... -5.5 -6.5 -7.5 -8.5
    spatial_ref  int64 8B 0
Data variables:
    temp         (y, x) float64 2kB 1.0 1.0 1.0 1.0 1.0 ... nan nan nan nan nan
Indexes:
  ┌ x        RasterIndex (crs=None)
  └ y

For join="exact", tolerance helps when two datasets cover the same area but differ only in their resolution. Here both datasets have origin 0.0 and 10 pixels, but slightly different spacing:

transforms_exact = [
    "0.0 1.0 0.0 1.0 0.0 -1.0",
    "0.0 1.00000001 0.0 1.0 0.0 -1.0",
]

dsets_exact = [
    xr.Dataset(
        {"temp": (("y", "x"), np.ones((10, 10)), {"grid_mapping": "spatial_ref"})},
        coords={
            "spatial_ref": (
                (),
                0,
                pyproj.CRS.from_epsg(4326).to_cf() | {"GeoTransform": transform},
            )
        },
    )
    for transform in transforms_exact
]
dsets_exact = list(map(rasterix.assign_index, dsets_exact))

with rasterix.set_options(transform_rtol=1e-7):
    exact = xr.align(*dsets_exact, join="exact")
exact[0]

<xarray.Dataset> Size: 968B
Dimensions:      (y: 10, x: 10)
Coordinates:
  * y            (y) float64 80B 0.5 -0.5 -1.5 -2.5 -3.5 ... -5.5 -6.5 -7.5 -8.5
  * x            (x) float64 80B 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5
    spatial_ref  int64 8B 0
Data variables:
    temp         (y, x) float64 800B 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0
Indexes:
  ┌ x        RasterIndex (crs=None)
  └ y