CRS handling

RasterIndex composes with xproj to provide CRS handling.

Show code cell source

Hide code cell source

%xmode minimal

import xarray as xr

import rasterix

xr.set_options(display_expand_indexes=True)

Show code cell output

Hide code cell output

Exception reporting mode: Minimal

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

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

To enable CRS-handling, start by assigning a CRS using the .proj.assign_crs function from xproj.

da = xr.open_dataarray(source, engine="rasterio")
da = da.proj.assign_crs(spatial_ref=da.spatial_ref.attrs["crs_wkt"])
da = da.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:
    spatial_ref  CRSIndex (crs=EPSG:3412)
  ┌ x            RasterIndex (crs=EPSG:3412)
  └ y
Attributes:
    AREA_OR_POINT:  Area

xarray.DataArray

'band_data'

• band: 1
• y: 332
• x: 316

• ...

  [104912 values with dtype=float32]

• Coordinates: (4)
  • band
    (band)
    int64
    1

    array([1])

  • y
    (y)
    float64
    4.338e+06 4.312e+06 ... -3.938e+06

    axis :

    Y

    long_name :

    Northing

    standard_name :

    projection_y_coordinate

    units :

    metre

    [332 values with dtype=float64]

  • x
    (x)
    float64
    -3.938e+06 -3.912e+06 ... 3.938e+06

    axis :

    X

    long_name :

    Easting

    standard_name :

    projection_x_coordinate

    units :

    metre

    [316 values with dtype=float64]

  • spatial_ref
    ()
    int64
    ...

    crs_wkt :

    PROJCS["NSIDC Sea Ice Polar Stereographic South",GEOGCS["Hughes 1980",DATUM["Hughes_1980",SPHEROID["Hughes 1980",6378273,298.279411123064,AUTHORITY["EPSG","7058"]],AUTHORITY["EPSG","1359"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","10345"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",-70],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",NORTH],AXIS["Northing",NORTH],AUTHORITY["EPSG","3412"]]

    semi_major_axis :

    6378273.0

    semi_minor_axis :

    6356889.449

    inverse_flattening :

    298.279411123064

    reference_ellipsoid_name :

    Hughes 1980

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    Hughes 1980

    horizontal_datum_name :

    Hughes 1980

    projected_crs_name :

    NSIDC Sea Ice Polar Stereographic South

    grid_mapping_name :

    polar_stereographic

    standard_parallel :

    -70.0

    straight_vertical_longitude_from_pole :

    0.0

    false_easting :

    0.0

    false_northing :

    0.0

    spatial_ref :

    PROJCS["NSIDC Sea Ice Polar Stereographic South",GEOGCS["Hughes 1980",DATUM["Hughes_1980",SPHEROID["Hughes 1980",6378273,298.279411123064,AUTHORITY["EPSG","7058"]],AUTHORITY["EPSG","1359"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","10345"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",-70],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",NORTH],AXIS["Northing",NORTH],AUTHORITY["EPSG","3412"]]

    [1 values with dtype=int64]

• Indexes: (3)
  • band
    PandasIndex

    PandasIndex(Index([1], dtype='int64', name='band'))

  • spatial_ref
    CRSIndex (crs=EPSG:3412)

    CRSIndex
    <Projected CRS: EPSG:3412>
    Name: NSIDC Sea Ice Polar Stereographic South
    Axis Info [cartesian]:
    - [north]: Easting (metre)
    - [north]: Northing (metre)
    Area of Use:
    - undefined
    Coordinate Operation:
    - name: unnamed
    - method: Polar Stereographic (variant B)
    Datum: Hughes 1980
    - Ellipsoid: Hughes 1980
    - Prime Meridian: Greenwich
    

  • x
    y
    RasterIndex (crs=EPSG:3412)

    RasterIndex(crs=EPSG:3412)
        AxisAffineTransformIndex(AxisAffineTransform(a=2.5e+04, b=0, c=-3.938e+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=-2.5e+04, f=4.338e+06, axis=Y, dim='y'))

• Attributes: (1)

  AREA_OR_POINT :

  Area

Note how the x variable has appropriate attributes!

da.x.attrs

{'axis': 'X',
 'long_name': 'Easting',
 'standard_name': 'projection_x_coordinate',
 'units': 'metre'}

Alignment

For demo purposes we’ll create a second copy with a different CRS.

diffcrs = da.copy(deep=True).drop_vars("spatial_ref").drop_indexes(["x", "y"])
diffcrs = diffcrs.proj.assign_crs(spatial_ref="epsg:4326").pipe(rasterix.assign_index)
diffcrs

<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 0
Indexes:
  ┌ x            RasterIndex (crs=EPSG:4326)
  └ y
    spatial_ref  CRSIndex (crs=EPSG:4326)
Attributes:
    AREA_OR_POINT:  Area

xarray.DataArray

'band_data'

• band: 1
• y: 332
• x: 316

• ...

  [104912 values with dtype=float32]

• Coordinates: (4)
  • band
    (band)
    int64
    1

    array([1])

  • y
    (y)
    float64
    4.338e+06 4.312e+06 ... -3.938e+06

    standard_name :

    latitude

    long_name :

    latitude coordinate

    units :

    degrees_north

    axis :

    Y

    [332 values with dtype=float64]

  • x
    (x)
    float64
    -3.938e+06 -3.912e+06 ... 3.938e+06

    standard_name :

    longitude

    long_name :

    longitude coordinate

    units :

    degrees_east

    axis :

    X

    [316 values with dtype=float64]

  • spatial_ref
    ()
    int64
    0

    array(0)

• Indexes: (3)
  • band
    PandasIndex

    PandasIndex(Index([1], dtype='int64', name='band'))

  • x
    y
    RasterIndex (crs=EPSG:4326)

    RasterIndex(crs=EPSG:4326)
        AxisAffineTransformIndex(AxisAffineTransform(a=2.5e+04, b=0, c=-3.938e+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=-2.5e+04, f=4.338e+06, axis=Y, dim='y'))

  • spatial_ref
    CRSIndex (crs=EPSG:4326)

    CRSIndex
    <Geographic 2D CRS: EPSG:4326>
    Name: WGS 84
    Axis Info [ellipsoidal]:
    - Lat[north]: Geodetic latitude (degree)
    - Lon[east]: Geodetic longitude (degree)
    Area of Use:
    - name: World.
    - bounds: (-180.0, -90.0, 180.0, 90.0)
    Datum: World Geodetic System 1984 ensemble
    - Ellipsoid: WGS 84
    - Prime Meridian: Greenwich
    

• Attributes: (1)

  AREA_OR_POINT :

  Area

Again note that the attributes on x have changed appropriately. This is enabled by RasterIndex’s integration with xproj

diffcrs.x.attrs

{'standard_name': 'longitude',
 'long_name': 'longitude coordinate',
 'units': 'degrees_east',
 'axis': 'X'}

Attempting to add the two raises an error (as it should)!

diffcrs + da

ValueError: raster indexes on objects to align do not have the same CRS
first index:
RasterIndex(crs=EPSG:4326)
    AxisAffineTransformIndex(AxisAffineTransform(a=2.5e+04, b=0, c=-3.938e+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=-2.5e+04, f=4.338e+06, axis=Y, dim='y'))

second index:
RasterIndex(crs=EPSG:3412)
    AxisAffineTransformIndex(AxisAffineTransform(a=2.5e+04, b=0, c=-3.938e+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=-2.5e+04, f=4.338e+06, axis=Y, dim='y'))