Indexing

We will demonstrate the raster index with a simple dataset with rectilinear coordinates and no rotation or skew. Both x and y coordinates are 1-dimensional.

Hide code cell source

 
%xmode minimal

import numpy as np
import xarray as xr

import rasterix

np.set_printoptions(threshold=10, edgeitems=2)
xr.set_options(display_expand_indexes=True)

Hide code cell output

 
Exception reporting mode: Minimal

<xarray.core.options.set_options at 0x7fcb24015400>

source = "../data/S_20240101_concentration_v4.0.tif"

da = xr.open_dataarray(source, engine="rasterio").pipe(rasterix.assign_index)
da

<xarray.DataArray 'band_data' (band: 1, y: 332, x: 316)> Size: 420kB
[104912 values with dtype=float32]
Coordinates:
  * band         (band) int64 8B 1
  * y            (y) float64 3kB 4.338e+06 4.312e+06 ... -3.912e+06 -3.938e+06
  * x            (x) float64 3kB -3.938e+06 -3.912e+06 ... 3.912e+06 3.938e+06
    spatial_ref  int64 8B ...
Indexes:
  ┌ x        RasterIndex (crs=None)
  └ y
Attributes:
    AREA_OR_POINT:  Area

Positional Indexing (isel)

See also

Now we will demonstrate indexing with a RasterIndex. See the Xarray docs and tutorial to understand the concepts underlying this section.

Slicing both x and y

Slicing preserves the laziness of coordinates since the new dataset can be represented by a new affine transform.

da_sliced = da.isel(x=slice(1, 4), y=slice(None, None, 2))
da_sliced

<xarray.DataArray 'band_data' (band: 1, y: 166, x: 3)> Size: 2kB
[498 values with dtype=float32]
Coordinates:
  * band         (band) int64 8B 1
  * y            (y) float64 1kB 4.338e+06 4.288e+06 ... -3.862e+06 -3.912e+06
  * x            (x) float64 24B -3.912e+06 -3.888e+06 -3.862e+06
    spatial_ref  int64 8B ...
Indexes:
  ┌ x        RasterIndex (crs=None)
  └ y
Attributes:
    AREA_OR_POINT:  Area

Compare the transforms before/after slicing.

(da.xindexes["x"].transform(), da_sliced.xindexes["x"].transform())

(Affine(25000.0, 0.0, -3950000.0,
        0.0, -25000.0, 4350000.0),
 Affine(25000.0, 0.0, -3925000.0,
        0.0, -50000.0, 4362500.0))

print(da_sliced.xindexes["x"])
print(da_sliced.xindexes["y"])

RasterIndex(crs=None)
    AxisAffineTransformIndex(AxisAffineTransform(a=2.5e+04, b=0, c=-3.912e+06, d=0, e=-2.5e+04, f=4.338e+06, axis=X, dim='x'))
    AxisAffineTransformIndex(AxisAffineTransform(a=2.5e+04, b=0, c=-3.938e+06, d=0, e=-5e+04, f=4.338e+06, axis=Y, dim='y'))
RasterIndex(crs=None)
    AxisAffineTransformIndex(AxisAffineTransform(a=2.5e+04, b=0, c=-3.912e+06, d=0, e=-2.5e+04, f=4.338e+06, axis=X, dim='x'))
    AxisAffineTransformIndex(AxisAffineTransform(a=2.5e+04, b=0, c=-3.938e+06, d=0, e=-5e+04, f=4.338e+06, axis=Y, dim='y'))

Outer indexing with array-like indexers

da_outer = da.isel(x=[0, 2, 4], y=[0, 0, 1])
da_outer

<xarray.DataArray 'band_data' (band: 1, y: 3, x: 3)> Size: 36B
[9 values with dtype=float32]
Coordinates:
  * band         (band) int64 8B 1
    y            (y) float64 24B 4.338e+06 4.338e+06 4.312e+06
    x            (x) float64 24B -3.938e+06 -3.888e+06 -3.838e+06
    spatial_ref  int64 8B ...
Attributes:
    AREA_OR_POINT:  Area

We cannot compute a new affine transform given arbitrary array positions, so the resulting dataset has no indexes associated with x and y.

da_outer.xindexes

Indexes:
    band     PandasIndex

Basic indexing with scalars

da_scalar = da.isel(x=0, y=1)
da_scalar

<xarray.DataArray 'band_data' (band: 1)> Size: 4B
[1 values with dtype=float32]
Coordinates:
  * band         (band) int64 8B 1
    spatial_ref  int64 8B ...
    x            float64 8B -3.938e+06
    y            float64 8B 4.312e+06
Attributes:
    AREA_OR_POINT:  Area
da_xscalar = da.isel(x=0)
da_xscalar

<xarray.DataArray 'band_data' (band: 1, y: 332)> Size: 1kB
[332 values with dtype=float32]
Coordinates:
  * band         (band) int64 8B 1
    y            (y) float64 3kB 4.338e+06 4.312e+06 ... -3.912e+06 -3.938e+06
    spatial_ref  int64 8B ...
    x            float64 8B -3.938e+06
Attributes:
    AREA_OR_POINT:  Area
# da_xscalar.xindexes["y"]  # should return an index

Vectorized (fancy) indexing

See also

See the Xarray tutorial for more on this topic.

Indexing the spatial coordinates with Xarray Variable or DataArray objects returns an unindexed object. The result cannot be represented by a RasterIndex in general.

da_points = da.isel(x=xr.DataArray([0, 1], dims="z"), y=xr.DataArray([1, 1], dims="z"))
da_points

<xarray.DataArray 'band_data' (band: 1, z: 2)> Size: 8B
[2 values with dtype=float32]
Coordinates:
  * band         (band) int64 8B 1
    x            (z) float64 16B -3.938e+06 -3.912e+06
    y            (z) float64 16B 4.312e+06 4.312e+06
    spatial_ref  int64 8B ...
Dimensions without coordinates: z
Attributes:
    AREA_OR_POINT:  Area
da_points2d = da.isel(
    x=xr.Variable(("u", "v"), [[0, 1], [2, 3]]),
    y=xr.Variable(("u", "v"), [[1, 1], [2, 2]]),
)
da_points2d

<xarray.DataArray 'band_data' (band: 1, u: 2, v: 2)> Size: 16B
[4 values with dtype=float32]
Coordinates:
  * band         (band) int64 8B 1
    x            (u, v) float64 32B -3.938e+06 -3.912e+06 -3.888e+06 -3.862e+06
    y            (u, v) float64 32B 4.312e+06 4.312e+06 4.288e+06 4.288e+06
    spatial_ref  int64 8B ...
Dimensions without coordinates: u, v
Attributes:
    AREA_OR_POINT:  Area

Label-based Indexing (sel)

Label-based indexing also preserves the RasterIndex where possible, like positional indexing.

da.sel(x=slice(-2e6, 2e6), y=slice(4e6, -2e6))

<xarray.DataArray 'band_data' (band: 1, y: 241, x: 161)> Size: 155kB
[38801 values with dtype=float32]
Coordinates:
  * band         (band) int64 8B 1
  * y            (y) float64 2kB 3.988e+06 3.962e+06 ... -1.988e+06 -2.012e+06
  * x            (x) float64 1kB -1.988e+06 -1.962e+06 ... 1.988e+06 2.012e+06
    spatial_ref  int64 8B ...
Indexes:
  ┌ x        RasterIndex (crs=None)
  └ y
Attributes:
    AREA_OR_POINT:  Area