cogeotiff / rio-tiler

User friendly Rasterio plugin to read raster datasets.
https://cogeotiff.github.io/rio-tiler/
BSD 3-Clause "New" or "Revised" License
505 stars 106 forks source link

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

Closed j08lue closed 1 week ago

j08lue commented 1 month ago

I created an rio_tiler.io.xarray.XarrayReader for a dataset that has longitude range 200 to 330 degrees East (which is not uncommon in NetCDF datasets) and tried to get a global tile at zoom level 0.

That would fail, saying the tile was out of bounds.

Once I subtracted 360 to get the longitudes into the -180, 180 range, it worked.

Should we handle datasets with wrapping longitude axes better in XarrayReader? Or document the tight criteria that files need to fulfill to work?

Here is a minimal example: https://gist.github.com/j08lue/6ec9ecc2a853e499e4fca98d8aa5fb78

The gist:

image

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

fails with

---------------------------------------------------------------------------
TileOutsideBounds                         Traceback (most recent call last)
Cell In[6], line 2
      1 with XarrayReader(da) as reader:
----> 2     img = reader.tile(0, 0, 0)

File ~/Code/titiler-xarray/.venv/lib/python3.12/site-packages/rio_tiler/io/xarray.py:235, in XarrayReader.tile(self, tile_x, tile_y, tile_z, tilesize, resampling_method, reproject_method, auto_expand, nodata)
    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
    238 if nodata is not None:

TileOutsideBounds: Tile 0/0/0 is outside bounds
j08lue commented 1 month ago

I see we have some code in TiTiler-Xarray that makes sure datasets conform to the format we expect, e.g.

    if crs == "epsg:4326" and (da.x > 180).any():
        # Adjust the longitude coordinates to the -180 to 180 range
        da = da.assign_coords(x=(da.x + 180) % 360 - 180)

        # Sort the dataset by the updated longitude coordinates
        da = da.sortby(da.x)

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

I guess we should at least add some checks to the XarrayReader, though, so it fails better?

vincentsarago commented 3 weeks ago

that seems reasonable