Handle Xarray dataset without CRS

[j08lue]

A pretty common case is that a DataArray passed to XarrayReader does not have a CRS defined (via rioxarray).

E.g. with the xarray tutorial data:

#!pip install rio_tiler xarray rioxarray pooch netCDF4

import xarray
from rio_tiler.io.xarray import XarrayReader

ds = xarray.tutorial.open_dataset("air_temperature")
da = ds["air"]
print(da.rio.crs)  # -> None

with XarrayReader(da) as reader:
    img = reader.tile(0, 0, 0)

CRSError: CRS is invalid: None

---------------------------------------------------------------------------
CRSError                                  Traceback (most recent call last)
File ~/Code/titiler-xarray/.venv/lib/python3.12/site-packages/rio_tiler/io/base.py:106, in SpatialMixin.tile_exists(self, tile_x, tile_y, tile_z)
    105 try:
--> 106     tile_bounds = transform_bounds(
    107         self.tms.rasterio_crs,
    108         self.crs,
    109         *tile_bounds,
    110         densify_pts=21,
    111     )
    112 except:  # noqa
    113     # HACK: gdal will first throw an error for invalid transformation
    114     # but if retried it will then pass.
    115     # Note: It might return `+/-inf` values

File ~/Code/titiler-xarray/.venv/lib/python3.12/site-packages/rasterio/warp.py:148, in transform_bounds(src_crs, dst_crs, left, bottom, right, top, densify_pts)
    147 src_crs = CRS.from_user_input(src_crs)
--> 148 dst_crs = CRS.from_user_input(dst_crs)
    149 return _transform_bounds(
    150     src_crs,
    151     dst_crs,
   (...)
    156     densify_pts,
    157 )

File rasterio/crs.pyx:783, in rasterio.crs.CRS.from_user_input()

CRSError: CRS is invalid: None

During handling of the above exception, another exception occurred:

CRSError                                  Traceback (most recent call last)
Cell In[8], line 11
      8 print(da.rio.crs)  # -> None
     10 with XarrayReader(da) as reader:
---> 11     img = reader.tile(0, 0, 0)

File ~/Code/titiler-xarray/.venv/lib/python3.12/site-packages/rio_tiler/io/xarray.py:234, in XarrayReader.tile(self, tile_x, tile_y, tile_z, tilesize, resampling_method, reproject_method, auto_expand, nodata)
    228     warnings.warn(
    229         "`resampling_method` is deprecated and will be removed in rio-tiler 7.0, please use `reproject_method`",
    230         DeprecationWarning,
    231     )
    232     reproject_method = resampling_method
--> 234 if not self.tile_exists(tile_x, tile_y, tile_z):
    235     raise TileOutsideBounds(f"Tile {tile_z}/{tile_x}/{tile_y} is outside bounds")
    237 ds = self.input

File ~/Code/titiler-xarray/.venv/lib/python3.12/site-packages/rio_tiler/io/base.py:116, in SpatialMixin.tile_exists(self, tile_x, tile_y, tile_z)
    106         tile_bounds = transform_bounds(
    107             self.tms.rasterio_crs,
    108             self.crs,
    109             *tile_bounds,
    110             densify_pts=21,
    111         )
    112     except:  # noqa
    113         # HACK: gdal will first throw an error for invalid transformation
    114         # but if retried it will then pass.
    115         # Note: It might return `+/-inf` values
--> 116         tile_bounds = transform_bounds(
    117             self.tms.rasterio_crs,
    118             self.crs,
    119             *tile_bounds,
    120             densify_pts=21,
    121         )
    123 # If tile_bounds has non-finite value in the dataset CRS we return True
    124 if not all(numpy.isfinite(tile_bounds)):

File ~/Code/titiler-xarray/.venv/lib/python3.12/site-packages/rasterio/warp.py:148, in transform_bounds(src_crs, dst_crs, left, bottom, right, top, densify_pts)
    121 """Transform bounds from src_crs to dst_crs.
    122 
    123 Optionally densifying the edges (to account for nonlinear transformations
   (...)
    145     Outermost coordinates in target coordinate reference system.
    146 """
    147 src_crs = CRS.from_user_input(src_crs)
--> 148 dst_crs = CRS.from_user_input(dst_crs)
    149 return _transform_bounds(
    150     src_crs,
    151     dst_crs,
   (...)
    156     densify_pts,
    157 )

File rasterio/crs.pyx:783, in rasterio.crs.CRS.from_user_input()

CRSError: CRS is invalid: None

In this case, da.rio.crs is None and XarrayReader.crs inherits that in this line

rio-tiler/rio_tiler/io/xarray.py

Line 74 in 8573dbe

self.crs = self.input.rio.crs

What is meant to happen instead? Should we default to WGS84? Or fail earlier, asking the user to make sure that da.rio.crs is defined?

We have this geographic_crs attribute, but it is not used anywhere, afaics.

https://github.com/cogeotiff/rio-tiler/blob/8573dbe7b3ea7325e7db084f29bc47f9797d0bbe/rio_tiler/io/xarray.py#L61C5-L61C19

[j08lue]

Same as

• XarrayReader considers tiles out of bounds when longitudes > 180? #718

Should we add a check like this to XarrayReader?

    # Make sure we have a valid CRS
    crs = da.rio.crs or "epsg:4326"
    da.rio.write_crs(crs, inplace=True)

https://github.com/developmentseed/titiler-xarray/blob/2be19eb9905d94cd6ac21f73947ad9802f9e92cb/titiler/xarray/reader.py#L170-L172

[vincentsarago]

What is meant to happen instead? Should we default to WGS84? Or fail earlier, asking the user to make sure that da.rio.crs is defined?

I think we should assert the crs is set. I would love to avoid setting WGS84 crs by default because we might get weird result if people use something else but forgot to add it correctly in the file.

Maybe we could have something based on env variable, like RIO_TILER_XARRAY_USE_WGS84=TRUE/FALSE