Fail to reproject and reproject_match a dataset with rotation affine.

Floabsolut commented 4 months ago

Dear all, I hope this description will be clear, and that I am not raising a pointless issue. I am available for any questions, Sincerely

Code Sample, a copy-pastable example

let's say you define a funtion to create a dataset with geographical information like the one bellow: create_gis_dataset()

import numpy as np
import xarray as xr
from affine import Affine
from rasterio.crs import CRS
from rioxarray.rioxarray import _generate_spatial_coords

def create_gis_dataset(
    image: np.ndarray,
    dst_affine: Affine,
    crs: CRS = 4326,
) -> xr.Dataset:
    """Create a rioxarray dataset object with coordinates vectors in good agreement with the input affine.

    Args:
    ----
        image (np.ndarray): An image in 2 or 3 dimensions, in uint16 dtpye
        dst_affine (Affine): The transform of the dataset, from which will be calculated the coordinates vectors
        crs (CRS, optional):  The coordinate reference system of the dataset, from which will be calculated the coordinates vectors.
        Defaults to 4326.

    Returns:
    -------
        xr.Dataset: with a proper coordinate synced with the dataset.rio.transform() and
        dataset.spatial_ref.GeoTransform

    """
    dims = ['band', 'y', 'x']
    xcs = (
        xr.DataArray(
            image,
            dims=dims,
            coords=_generate_spatial_coords(
                affine=dst_affine,
                width=image.shape[2],
                height=image.shape[1],
            ),
        )
        .astype('uint16')
        .rio.write_nodata(0, inplace=True)
        .rio.write_crs(crs, inplace=True)
        .rio.write_transform(dst_affine, inplace=True)
        .rio.write_coordinate_system(inplace=True)
    )
    output: xr.Dataset = xcs.to_dataset(name='image', promote_attrs=True)
    return output

Now you create a dataset with this method, using an affine that does not satisfied the condition: if affine.is_rectilinear and not _affine_has_rotation(affine)

for example

image = np.random.randint(0, 255, size=(1, 100, 100), dtype=np.uint8)
affine = Affine(1, 0.2, 0, 0, 1, 0)

dataset = create_gis_dataset(image, affine)

Then you want to reproject this dataset to the affine Affine(1, 0, 0, 0, 1, 0)

squared_dataset = dataset.copy(deep=True)
squared_dataset = squared_dataset.rio.write_transform(transform=Affine(1, 0, 0, 0, 1, 0), inplace=True)

squared_dataset = squared_dataset.rio.reproject_match(dataset)

This reproject_match will fail with this error:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/opt/homebrew/lib/python3.11/site-packages/rioxarray/raster_dataset.py", line 121, in reproject
    self._obj[var]
  File "/opt/homebrew/lib/python3.11/site-packages/rioxarray/raster_array.py", line 482, in reproject
    coords=_make_coords(
           ^^^^^^^^^^^^^
  File "/opt/homebrew/lib/python3.11/site-packages/rioxarray/rioxarray.py", line 166, in _make_coords
    coords = _get_nonspatial_coords(src_data_array)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/homebrew/lib/python3.11/site-packages/rioxarray/rioxarray.py", line 144, in _get_nonspatial_coords
    coords[coord] = xarray.IndexVariable(
                    ^^^^^^^^^^^^^^^^^^^^^
  File "/opt/homebrew/lib/python3.11/site-packages/xarray/core/variable.py", line 2637, in __init__
    raise ValueError(f"{type(self).__name__} objects must be 1-dimensional")
ValueError: IndexVariable objects must be 1-dimensional

Problem description

My understanding of why it fails:

_generate_spatial_coords will output a coords with a value which is not 1 dimension.

Then, trying to reproject will fail with the error in this line here cause IndexVariable is expecting a 1-Dimensional object. My understanding so far is that rioxarray might not handle well "not squared" affines?

Environment Information

python -c "import rioxarray; rioxarray.show_versions()"

rioxarray version 0.15.1
rasterio version 1.3.9
GDAL version 3.6.4
Python version 3.11.0 (main, Mar 1 2023, 18:26:19) [GCC 11.2.0]
Operation System Information Linux-5.15.133.1-microsoft-standard-WSL2-x86_64-with-glibc2.35

Installation method

pip install rioxarray

snowman2 commented 2 months ago

Thank you for this issue - it has all of the elements needed to re-produce the issue. Yes, rotated raster support could definitely be improved.

Floabsolut commented 2 months ago

Dear @snowman2, thanks for the answer and validation of the issue. Let me know if I could help in any way.

snowman2 commented 2 months ago

See #772

Let me know if I could help in any way.

Thank you. This issue was very helpful. One step closer to better support rotated rasters.

corteva / rioxarray