check if the Y bounds are inverted and flip the image on the Y axis

[vincentsarago]

possible solution for developmentseed/titiler-cmr#34

import earthaccess
import matplotlib.pyplot as plt
import rio_tiler

earthaccess.login()

results = earthaccess.search_data(
    concept_id="C2036881735-POCLOUD", count=1, temporal=("2024-09-20", "2024-09-20")
)
fp = earthaccess.download(results, "earthaccess_data")[0]
ds = xarray.open_dataset(fp, decode_coords="all")

da = ds["analysed_sst"]
if not da.rio.crs:
    da = da.rio.write_crs("epsg:4326")

with XarrayReader(da) as src:
    img = src.part([-90.0,22.5,-67.5,45.0])

print(img.crs)
plt.imshow(img.data_as_image())

cc @maxrjones @hrodmn

[vincentsarago]

We check if the unordered bounds are in form of xmin, ymax, xmax, ymin and if the transform's y resolution is positive.

Those 2 clues tells us the image if flipped

I don't know enough about NetCDF/rioxarray to tell if this will be enough

[vincentsarago]

Just a quick note: This PR seems to fix our original issue but I'm really not sure about the logic. I'm going to pause on this because I want to see first if it's not something in rioxarray

[vincentsarago]

changed my mind! I think it's not in a better shape!

from why I understand, in some case, rioxarray consider somehow the origin of the image to be the bottom left, while in rasterio/gdal world it's the top left. Once we apply reprojection, GDAL/Rasterio will deal with it (using the transform) and flip the image.

When we do not have reprojection, the image might be upside down so we check the unordered_bounds and the Y resolution to see if we need to flip the image vertically.

[hrodmn]

Thanks for figuring out a quick solution @vincentsarago! I left a few small suggestions but nothing blocking.

[maxrjones]

I checked what happens if we write this DataArray to a GeoTiff and try with the rasterio reader:

import numpy as np
import xarray as xr
import rioxarray # noqa
import rio_tiler
from datetime import datetime
import matplotlib.pyplot as plt
from odc.geo.geom import Geometry
from shapely.geometry import box

arr = np.arange(0, 33 * 35).reshape(1, 33, 35)
da = xr.DataArray(
    arr,
    dims=("time", "y", "x"),
    coords={
        "x": np.arange(-170, 180, 10),
        "y": np.arange(-80, 85, 5),
        "time": [datetime(2022, 1, 1)],
    },
)
da.attrs.update({"valid_min": arr.min(), "valid_max": arr.max()})
da.rio.write_crs("epsg:4326", inplace=True)

da.rio.to_raster("test.tif")

# Set up GeoJSON for reader
bbox = box(-99.524, 19.043, -78.647, 31.126)
gbox = Geometry(bbox, "EPSG:4326")
shape = gbox.geojson()

r = rio_tiler.io.rasterio.Reader(f"test.tif")
im = r.feature(shape)
plt.imshow(im.data.squeeze())
plt.colorbar()
plt.show()

result:

WindowError: Bounds and transform are inconsistent

This is where the error is raised in rasterio if you want to use the same logic for checking for inverted coordinates here:
https://github.com/rasterio/rasterio/blob/74ccaf126d08fc6eca3eacd7cb20ac8bb155ee3b/rasterio/windows.py#L319-L324

However, as I mentioned in slack I am concerned about the slippery slope of handling these cases within rio-tiler. It is extremely common to expect users to transform their coordinates before handing the data to your library. For example, rioxarray and pyresample expects the user to set the x and y dimension/coordinate names, crs, and the proper axis order. Rio-tiler could similarly expect the bounds and transform to be consistent. Obviously your call, but I would lean towards raising an error in this case and fixing it at the titiler or custom reader level rather than here for consistency with the rasterio reader.

If you do want to flip the data in rio-tiler, here is the most common method for making sure your 'y' coordinates are in the expected order (which can vary depending on how you're trying to use the data) using xarray rather than numpy, which can help readability:
https://github.com/carbonplan/cmip6-downscaling/blob/2ed6c5e2d92a4307f7161d9c2431b01d3c10798b/cmip6_downscaling/data/cmip.py#L57-L58

[vincentsarago]

I checked what happens if we write this DataArray to a GeoTiff and try with the rasterio reader:

This tells me that rioxarray does something weird TBH. I feel that rioxarray expect origin to be the bottom right, while it's common to assume the origin is top left 🤷

If we don't fix this here, we would add a IF test to check if the data is in correct order 🤷 so I feel that doing this is OK (I may regret it)