Unable to reproject epsg:32645 to epsg:4326

[zxdawn]

I have one tif file whose projection is epsg:32645 and it has transform:

import rioxarray
import xarray

test = xr.open_dataset('/Users/xinz/Downloads/test.tif')

print(test.rio.crs)
print(test.rio.transform())

Output:

EPSG:32645
| 321540.00, 30.00, 0.00|
| 3084990.00, 0.00,-30.00|
| 0.00, 0.00, 1.00|

I have tried to reproject it into epsg:4326:

test.rio.reproject('EPSG:4326')

Got this error:

---------------------------------------------------------------------------
CPLE_AppDefinedError                      Traceback (most recent call last)
Cell In[14], line 3
      1 test = xr.open_dataset('[/Users/xinz/Downloads/test.tif](https://file+.vscode-resource.vscode-cdn.net/Users/xinz/Downloads/test.tif)')
----> 3 test.rio.reproject('EPSG:4326')

File [~/miniconda3/envs/enpt/lib/python3.10/site-packages/rioxarray/raster_dataset.py:121](https://file+.vscode-resource.vscode-cdn.net/Users/xinz/Documents/github/enmap/notebooks/~/miniconda3/envs/enpt/lib/python3.10/site-packages/rioxarray/raster_dataset.py:121), in RasterDataset.reproject(self, dst_crs, resolution, shape, transform, resampling, nodata, **kwargs)
    118 try:
    119     x_dim, y_dim = _get_spatial_dims(self._obj, var)
    120     resampled_dataset[var] = (
--> 121         self._obj[var]
    122         .rio.set_spatial_dims(x_dim=x_dim, y_dim=y_dim, inplace=True)
    123         .rio.reproject(
    124             dst_crs,
    125             resolution=resolution,
    126             shape=shape,
    127             transform=transform,
    128             resampling=resampling,
    129             nodata=nodata,
    130             **kwargs,
    131         )
    132     )
    133 except MissingSpatialDimensionError:
    134     if len(self._obj[var].dims) >= 2 and not get_option(
    135         SKIP_MISSING_SPATIAL_DIMS
...
File rasterio/_warp.pyx:690, in rasterio._warp._calculate_default_transform()

File rasterio/_err.pyx:221, in rasterio._err.exc_wrap_pointer()

CPLE_AppDefinedError: Cannot invert geotransform

Any idea how to reproject it correctly?

[Kirill888]

@zxdawn looks like messed up transform to me. Probably stored in row major order inside TIFF instead of GDAL order.

[zxdawn]

Thanks, but I got the same error after updating the transform to row order:

import rioxarray
import xarray as xr
import rasterio

test = xr.open_dataset('/Users/xinz/Downloads/test.tif')

row_order = test.rio.transform().to_shapely()

print('Wrong transform', test.rio.transform())
print('Row order: ', row_order)

row_transform = rasterio.Affine(row_order[0], row_order[1], row_order[2], row_order[3], row_order[4], row_order[5])
test.rio.write_transform(row_transform, inplace=True)

print('New transform: ', test.rio.transform())

test.rio.reproject('EPSG:4326')

Wrong transform | 321540.00, 30.00, 0.00|
| 3084990.00, 0.00,-30.00|
| 0.00, 0.00, 1.00|
Row order:  (321540.0, 30.0, 3084990.0, 0.0, 0.0, -30.0)
New transform:  | 321540.00, 30.00, 3084990.00|
| 0.00, 0.00,-30.00|
| 0.00, 0.00, 1.00|

[Kirill888]

row major would be like this:

| 30.0,  0.00, 321540.00|
| 0.00, -30.0, 3084990.00|
| 0.00,  0.00, 1.00|

sx, 0, tx, 0, sy, ty, 0,0,1

[zxdawn]

Thanks a lot, @Kirill888. I finally found out that's caused by the wrong definition of transform. I should use Affine.from_gdal() instead of Affine() when writing the transform.

# write crs and geotransform
dataset.rio.write_transform(Affine.from_gdal(*geotransform), inplace=True)
dataset.rio.write_crs(CRS.from_wkt(spatial_ref), inplace=True)