Periodic axes

RasterIndex accepts x_period (or y_period), in assign_index() or RasterIndex.from_transform() to indicate periodic axes:

ds = rasterix.assign_index(ds, x_period=360)

rasterix does not infer the period from the CRS (see Design Choices).

Show code cell content

Hide code cell content

import numpy as np
import xarray as xr
from affine import Affine
from rasterix import RasterIndex

index = RasterIndex.from_transform(
    Affine(30, 0, -180, 0, -30, 90), width=12, height=2, x_period=360
)
da = xr.DataArray(
    np.broadcast_to(np.arange(12), (2, 12)),
    dims=("y", "x"),
    coords=xr.Coordinates.from_xindex(index),
    name="column",
)
da

<xarray.DataArray 'column' (y: 2, x: 12)> Size: 192B
array([[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11],
       [ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11]])
Coordinates:
  * y        (y) float64 16B 75.0 45.0
  * x        (x) float64 96B -165.0 -135.0 -105.0 -75.0 ... 105.0 135.0 165.0
Indexes:
  ┌ x        RasterIndex (crs=None, x_period=360)
  └ y

xarray.DataArray

'column'

• y: 2
• x: 12

• 0 1 2 3 4 5 6 7 8 9 10 11 0 1 2 3 4 5 6 7 8 9 10 11

  array([[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11],
         [ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11]])

• Coordinates: (2)
  • y
    (y)
    float64
    75.0 45.0

    [2 values with dtype=float64]

  • x
    (x)
    float64
    -165.0 -135.0 ... 135.0 165.0

    [12 values with dtype=float64]

• Indexes: (1)
  • x
    y
    RasterIndex (crs=None, x_period=360)

    RasterIndex(crs=None)
        AxisAffineTransformIndex(AxisAffineTransform(a=30, b=0, c=-165, d=0, e=-30, f=75, axis=X, dim='x', period=360))
        AxisAffineTransformIndex(AxisAffineTransform(a=30, b=0, c=-165, d=0, e=-30, f=75, axis=Y, dim='y'))

Read the period back with RasterIndex.x_period and RasterIndex.y_period:

index.x_period, index.y_period

(360, None)

Requirements

period / |dx| must be an integer number of cells, to within 1e-3 cell by default (set period_atol / period_rtol with rasterix.set_options()). For example, 360 / 0.7 is not:

RasterIndex.from_transform(
    Affine(0.7, 0, -180, 0, -1, 90), width=514, height=2, x_period=360
)

ValueError: period=360 is not an integer number of cells (period / |dx| = 514.2857142857143). If this is round-off in dx, loosen the tolerance with rasterix.set_options(period_atol=...) (now 0.001 cells).

The axis must not have more cells than one period. A grid made with np.linspace(-180, 180, n) has a copy of the first column at the end; drop it first with ds.isel(x=slice(0, -1)):

RasterIndex.from_transform(
    Affine(30, 0, -195, 0, -30, 90), width=13, height=2, x_period=360
)

ValueError: The 'x' axis has 13 cells, but one period has 12. Drop the duplicate seam cells first, e.g. ``.isel(x=slice(0, 12))``.

Behavior

All rules work in cell positions, so they also apply to a decreasing axis (dx < 0). On a decreasing axis, write slices in decreasing label order, as usual in Xarray.

Labels outside the extent wrap

If a label is in the extent i.e. ≤ period, it is used as-is. If not, the label is wrapped by the period. Thus −100° and 260° give the same column:

da.sel(x=-100, method="nearest").values, da.sel(x=260, method="nearest").values

(array([2, 2]), array([2, 2]))

Slices across the period seam

start < stop

When slicing across the period seam (e.g. the anti-meridian), start is wrapped and the width stop - start is kept. The cells come from both ends of the array, and the coordinates of the result are unwrapped and contiguous.

For example, on a global longitude axis, a slice that starts left of the extent, say -180, returns coordianates in that frame. Below, on a 1° grid from −180° to 180°, slice(-190, -170) gives coordinates −189.5 … −169.5.

global_index = RasterIndex.from_transform(
    Affine(1, 0, -180, 0, -1, 90), width=360, height=1, x_period=360
)
global_da = xr.DataArray(
    np.arange(360)[np.newaxis],
    dims=("y", "x"),
    coords=xr.Coordinates.from_xindex(global_index),
)
global_da.sel(x=slice(-190, -170))

<xarray.DataArray (y: 1, x: 21)> Size: 168B
array([[350, 351, 352, 353, 354, 355, 356, 357, 358, 359,   0,   1,   2,
          3,   4,   5,   6,   7,   8,   9,  10]])
Coordinates:
  * y        (y) float64 8B 89.5
  * x        (x) float64 168B -189.5 -188.5 -187.5 ... -171.5 -170.5 -169.5
Indexes:
  ┌ x        RasterIndex (crs=None, x_period=360)
  └ y

xarray.DataArray

• y: 1
• x: 21

• 350 351 352 353 354 355 356 357 358 359 0 1 2 3 4 5 6 7 8 9 10

  array([[350, 351, 352, 353, 354, 355, 356, 357, 358, 359,   0,   1,   2,
            3,   4,   5,   6,   7,   8,   9,  10]])

• Coordinates: (2)
  • y
    (y)
    float64
    89.5

    [1 values with dtype=float64]

  • x
    (x)
    float64
    -189.5 -188.5 ... -170.5 -169.5

    [21 values with dtype=float64]

• Indexes: (1)
  • x
    y
    RasterIndex (crs=None, x_period=360)

    RasterIndex(crs=None)
        AxisAffineTransformIndex(AxisAffineTransform(a=1, b=0, c=-189.5, d=0, e=-1, f=89.5, axis=X, dim='x', period=360))
        AxisAffineTransformIndex(AxisAffineTransform(a=1, b=0, c=-179.5, d=0, e=-1, f=89.5, axis=Y, dim='y'))

stop < start

Behaviour is similar at the other end of the axis:

da.sel(x=slice(150, -150))

<xarray.DataArray 'column' (y: 2, x: 3)> Size: 48B
array([[11,  0,  1],
       [11,  0,  1]])
Coordinates:
  * y        (y) float64 16B 75.0 45.0
  * x        (x) float64 24B 165.0 195.0 225.0
Indexes:
  ┌ x        RasterIndex (crs=None, x_period=360)
  └ y

xarray.DataArray

'column'

• y: 2
• x: 3

• 11 0 1 11 0 1

  array([[11,  0,  1],
         [11,  0,  1]])

• Coordinates: (2)
  • y
    (y)
    float64
    75.0 45.0

    [2 values with dtype=float64]

  • x
    (x)
    float64
    165.0 195.0 225.0

    [3 values with dtype=float64]

• Indexes: (1)
  • x
    y
    RasterIndex (crs=None, x_period=360)

    RasterIndex(crs=None)
        AxisAffineTransformIndex(AxisAffineTransform(a=30, b=0, c=165, d=0, e=-30, f=75, axis=X, dim='x', period=360))
        AxisAffineTransformIndex(AxisAffineTransform(a=30, b=0, c=-165, d=0, e=-30, f=75, axis=Y, dim='y'))

A slice that contains the full extent is not rolled

We never duplicate data of the array by slicing.

da.sel(x=slice(-200, 200)).x.values

array([-165., -135., -105.,  -75.,  -45.,  -15.,   15.,   45.,   75.,
        105.,  135.,  165.])

Period is preserved

da.sel(x=slice(120, 210)) has coordinates 135° … 225° and preserves x_period.

subset = da.sel(x=slice(120, 210))
subset

<xarray.DataArray 'column' (y: 2, x: 4)> Size: 64B
array([[10, 11,  0,  1],
       [10, 11,  0,  1]])
Coordinates:
  * y        (y) float64 16B 75.0 45.0
  * x        (x) float64 32B 135.0 165.0 195.0 225.0
Indexes:
  ┌ x        RasterIndex (crs=None, x_period=360)
  └ y

xarray.DataArray

'column'

• y: 2
• x: 4

• 10 11 0 1 10 11 0 1

  array([[10, 11,  0,  1],
         [10, 11,  0,  1]])

• Coordinates: (2)
  • y
    (y)
    float64
    75.0 45.0

    [2 values with dtype=float64]

  • x
    (x)
    float64
    135.0 165.0 195.0 225.0

    [4 values with dtype=float64]

• Indexes: (1)
  • x
    y
    RasterIndex (crs=None, x_period=360)

    RasterIndex(crs=None)
        AxisAffineTransformIndex(AxisAffineTransform(a=30, b=0, c=135, d=0, e=-30, f=75, axis=X, dim='x', period=360))
        AxisAffineTransformIndex(AxisAffineTransform(a=30, b=0, c=-165, d=0, e=-30, f=75, axis=Y, dim='y'))

So −160° (the same as 200°) still selects column 0:

subset.sel(x=-160, method="nearest")

<xarray.DataArray 'column' (y: 2)> Size: 16B
array([0, 0])
Coordinates:
    y        (y) float64 16B 75.0 45.0
    x        float64 8B 195.0

xarray.DataArray

'column'

• y: 2

• 0 0

  array([0, 0])

• Coordinates: (2)
  • y
    (y)
    float64
    75.0 45.0

    [2 values with dtype=float64]

  • x
    ()
    float64
    195.0

    array(195.)

Alignment and arithmetic need the same frame

Alignment compares labels exactly; it does not wrap them. subset has labels 195 and 225, which are the same cells as −165 and −135 in da. Exact matching would drop these cells (inner join) or duplicate them (outer join), so alignment raises an error. Use the same frame on both sides before you combine.

da + subset

ValueError: Cannot align RasterIndexes on periodic dimension 'x': they share cells at labels that differ by 360 (1 period(s) of 360). Alignment compares labels exactly, so these cells would be dropped or duplicated. Select both objects in the same frame first.

Concatenation across the seam

The end and the start of a periodic axis concatenate to a contiguous result. See Periodic axes for the rule.

xr.concat(
    [da.isel(x=slice(-2, None)), da.isel(x=slice(0, 2))], dim="x"
)

<xarray.DataArray 'column' (y: 2, x: 4)> Size: 64B
array([[10, 11,  0,  1],
       [10, 11,  0,  1]])
Coordinates:
  * y        (y) float64 16B 75.0 45.0
  * x        (x) float64 32B 135.0 165.0 195.0 225.0
Indexes:
  ┌ x        RasterIndex (crs=None, x_period=360)
  └ y

xarray.DataArray

'column'

• y: 2
• x: 4

• 10 11 0 1 10 11 0 1

  array([[10, 11,  0,  1],
         [10, 11,  0,  1]])

• Coordinates: (2)
  • y
    (y)
    float64
    75.0 45.0

    [2 values with dtype=float64]

  • x
    (x)
    float64
    135.0 165.0 195.0 225.0

    [4 values with dtype=float64]

• Indexes: (1)
  • x
    y
    RasterIndex (crs=None, x_period=360)

    RasterIndex(crs=None)
        AxisAffineTransformIndex(AxisAffineTransform(a=30, b=0, c=135, d=0, e=-30, f=75, axis=X, dim='x', period=360))
        AxisAffineTransformIndex(AxisAffineTransform(a=30, b=0, c=135, d=0, e=-30, f=75, axis=Y, dim='y'))

Indexes with different periods do not combine

no_period = da.copy().assign_coords(
    xr.Coordinates.from_xindex(
        RasterIndex.from_transform(
            Affine(30, 0, -180, 0, -30, 90), width=12, height=2
        )
    )
)
xr.align(da, no_period, join="outer")

ValueError: Cannot combine RasterIndexes with different periods: ['(360, None)', '(None, None)']