Coordinate Reference System Management

xarray “… is particularly tailored to working with netCDF files, which were the source of xarray’s data model…” (http://xarray.pydata.org).

For netCDF files, the GIS community uses CF conventions (http://cfconventions.org/).

Additionally, GDAL also supports these attributes:

  • spatial_ref (Well Known Text)

  • GeoTransform (GeoTransform array)

References:

Operations on xarray objects can cause data loss. Due to this, rioxarray writes and expects the spatial reference information to exist in the coordinates.

Accessing the CRS object

If you have opened a dataset and the Coordinate Reference System (CRS) can be determined, you can access it via the rio.crs accessor.

Search order for the CRS (DataArray and Dataset):

  1. Look in attributes (attrs) of your data array for the grid_mapping coordinate name. Inside the grid_mapping coordinate first look for spatial_ref then crs_wkt and lastly the CF grid mapping attributes. This is in line with the Climate and Forecast (CF) conventions for storing the CRS as well as GDAL netCDF conventions.

  2. Look in the crs attribute and load in the CRS from there. This is for backwards compatibility with xarray.open_rasterio, which is deprecated since version 0.20.0. We recommend using rioxarray.open_rasterio instead.

The value for the crs is anything accepted by rasterio.crs.CRS.from_user_input()

Search order for the CRS for Dataset:

If the CRS is not found using the search methods above, it also searches the data_vars and uses the first valid CRS found.

decode_coords=”all”

If you use one of xarray’s open methods such as xarray.open_dataset to load netCDF files with the default engine, it is recommended to use decode_coords="all". This will load the grid mapping variable into coordinates for compatibility with rioxarray.

API Documentation

[1]:

import rioxarray  # activate the rio accessor
import xarray
from affine import Affine

[2]:

rds = xarray.open_dataset("../../test/test_data/input/PLANET_SCOPE_3D.nc", decode_coords="all")

[3]:

rds.green.attrs

[3]:

{'units': 'DN', 'nodata': 0.0}

[4]:

rds.green.spatial_ref

[4]:

<xarray.DataArray 'spatial_ref' ()>
array(0)
Coordinates:
    spatial_ref  int64 0
Attributes:
    spatial_ref:  PROJCS["WGS 84 / UTM zone 22S",GEOGCS["WGS 84",DATUM["WGS_1...
[5]:

rds.green.rio.crs

[5]:

CRS.from_epsg(32722)

Setting the CRS

Use the rio.write_crs method to set the CRS on your xarray.Dataset or xarray.DataArray. This modifies the xarray.Dataset or xarray.DataArray and sets the CRS in a CF compliant manner.

Note: It is recommended to use rio.write_crs() if you want the CRS to persist on the Dataset/DataArray and to write the CRS CF compliant metadata. Calling only rio.set_crs() CRS storage method is lossy and will not modify the Dataset/DataArray metadata.

[6]:

xda = xarray.DataArray(1)
xda.rio.write_crs(4326, inplace=True)
xda.spatial_ref

[6]:

<xarray.DataArray 'spatial_ref' ()>
array(0)
Coordinates:
    spatial_ref  int64 0
Attributes:
    crs_wkt:                      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["...
    semi_major_axis:              6378137.0
    semi_minor_axis:              6356752.314245179
    inverse_flattening:           298.257223563
    reference_ellipsoid_name:     WGS 84
    longitude_of_prime_meridian:  0.0
    prime_meridian_name:          Greenwich
    geographic_crs_name:          WGS 84
    grid_mapping_name:            latitude_longitude
    spatial_ref:                  GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["...
[7]:

xda.rio.crs

[7]:

CRS.from_epsg(4326)

Spatial dimensions

Only 1-dimensional X and Y dimensions are supported.

The expected X/Y dimension names searched for in the coords are:

  • x | y

  • longitude | latitude

  • Coordinates (coords) with the CF attributes in attrs:

    • axis: X | Y

    • standard_name: longitude | latitude or projection_x_coordinate | projection_y_coordinate

Option 1: Write the CF attributes for non-standard dimension names

If you don’t want to rename your dimensions/coordinates, you can write the CF attributes so the coordinates can be found.

[ ]:

rds.rio.write_crs(
    4326
    inplace=True,
).rio.set_spatial_dims(
    x_dim="lon",
    y_dim="lat"
    inplace=True,
).rio.write_coordinate_system(inplace=True)

Option 2: Rename your coordinates

xarray.Dataset.rename

[ ]:

rds = rds.rename(lon=longitude, lat=latitude)

Setting the transform of the dataset

The transform can be calculated from the coordinates of your data. This method is useful if your netCDF file does not have coordinates present. Use the rio.write_transform method to set the transform on your xarray.Dataset or xarray.DataArray.

[8]:

transform = Affine(3.0, 0.0, 466266.0, 0.0, -3.0, 8084700.0)
xda.rio.write_transform(transform, inplace=True)
xda.spatial_ref.GeoTransform

[8]:

'466266.0 3.0 0.0 8084700.0 0.0 -3.0'

[9]:

xda.rio.transform()

[9]:

Affine(3.0, 0.0, 466266.0,
       0.0, -3.0, 8084700.0)