Fail to reproject and reproject_match a dataset with rotation affine.
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