opendatacube / odc-geo

GeoBox and geometry utilities extracted from datacube-core
https://odc-geo.readthedocs.io/en/latest/
Apache License 2.0
80 stars 12 forks source link

`xr_reproject`: CPLE_AppDefinedError: Unable to compute a BILINEAR based transformation #159

Closed maawoo closed 3 months ago

maawoo commented 4 months ago

Hi @Kirill888,

I'm running into the following issue when trying to downsample the resolution of a DataArray by reprojecting in the same CRS:

test_40 = xr_reproject(test, how=test.odc.crs, resolution=40, method='bilinear')
CPLE_AppDefinedError: Unable to compute a BILINEAR based transformation between pixel/line and georeferenced coordinates for MEM:::DATAPOINTER=13593378816,PIXELS=4680,LINES=11820,BANDS=1,DATATYPE=Float32.
Full traceback:
```python --------------------------------------------------------------------------- CPLE_AppDefinedError Traceback (most recent call last) Cell In[64], line 1 ----> 1 test_40.compute() File ~/Documents/pypypy/woody/.pixi/envs/default/lib/python3.12/site-packages/xarray/core/dataarray.py:1178, in DataArray.compute(self, **kwargs) 1153 """Manually trigger loading of this array's data from disk or a 1154 remote source into memory and return a new array. 1155 (...) 1175 dask.compute 1176 """ 1177 new = self.copy(deep=False) -> 1178 return new.load(**kwargs) File ~/Documents/pypypy/woody/.pixi/envs/default/lib/python3.12/site-packages/xarray/core/dataarray.py:1146, in DataArray.load(self, **kwargs) 1126 def load(self, **kwargs) -> Self: 1127 """Manually trigger loading of this array's data from disk or a 1128 remote source into memory and return this array. 1129 (...) 1144 dask.compute 1145 """ -> 1146 ds = self._to_temp_dataset().load(**kwargs) 1147 new = self._from_temp_dataset(ds) 1148 self._variable = new._variable File ~/Documents/pypypy/woody/.pixi/envs/default/lib/python3.12/site-packages/xarray/core/dataset.py:862, in Dataset.load(self, **kwargs) 859 chunkmanager = get_chunked_array_type(*lazy_data.values()) 861 # evaluate all the chunked arrays simultaneously --> 862 evaluated_data: tuple[np.ndarray[Any, Any], ...] = chunkmanager.compute( 863 *lazy_data.values(), **kwargs 864 ) 866 for k, data in zip(lazy_data, evaluated_data): 867 self.variables[k].data = data File ~/Documents/pypypy/woody/.pixi/envs/default/lib/python3.12/site-packages/xarray/namedarray/daskmanager.py:86, in DaskManager.compute(self, *data, **kwargs) 81 def compute( 82 self, *data: Any, **kwargs: Any 83 ) -> tuple[np.ndarray[Any, _DType_co], ...]: 84 from dask.array import compute ---> 86 return compute(*data, **kwargs) File ~/Documents/pypypy/woody/.pixi/envs/default/lib/python3.12/site-packages/dask/base.py:661, in compute(traverse, optimize_graph, scheduler, get, *args, **kwargs) 658 postcomputes.append(x.__dask_postcompute__()) 660 with shorten_traceback(): --> 661 results = schedule(dsk, keys, **kwargs) 663 return repack([f(r, *a) for r, (f, a) in zip(results, postcomputes)]) File ~/Documents/pypypy/woody/.pixi/envs/default/lib/python3.12/site-packages/odc/geo/_dask.py:59, in _do_chunked_reproject() 56 src = ba.extract(src_nodata, dtype=dtype, casting=casting, roi=src_roi) 57 dst_roi = ba.with_yx(src_roi, np.s_[:, :]) ---> 59 _ = _rio_reproject( 60 src, 61 dst[dst_roi], 62 src_gbox, 63 dst_gbox, 64 resampling=resampling, 65 src_nodata=src_nodata, 66 dst_nodata=dst_nodata, 67 **kwargs, 68 ) 70 return dst File ~/Documents/pypypy/woody/.pixi/envs/default/lib/python3.12/site-packages/odc/geo/warp.py:215, in _rio_reproject() 212 src, src_is_bool = _alias_or_convert(src) 213 _dst, _ = _alias_or_convert(dst) --> 215 rasterio.warp.reproject( 216 src, 217 _dst, 218 src_transform=src_transform, 219 gcps=gcps, 220 src_crs=str(s_gbox.crs), 221 dst_transform=d_gbox.transform, 222 dst_crs=str(d_gbox.crs), 223 resampling=resampling, 224 src_nodata=src_nodata, 225 dst_nodata=dst_nodata, 226 **kwargs, 227 ) 229 if dst is not _dst: 230 # int8 workaround copy pixels back to int8 231 if src_is_bool: 232 # undo [0, 1] to [0, 255] stretching of the src File ~/Documents/pypypy/woody/.pixi/envs/default/lib/python3.12/site-packages/rasterio/env.py:401, in wrapper() 399 else: 400 with Env.from_defaults(): --> 401 return f(*args, **kwds) File ~/Documents/pypypy/woody/.pixi/envs/default/lib/python3.12/site-packages/rasterio/warp.py:344, in reproject() 340 destination = np.empty((int(dst_count), int(dst_height), int(dst_width)), 341 dtype=source.dtype) 343 # Call the function in our extension module. --> 344 _reproject( 345 source, destination, src_transform=src_transform, gcps=gcps, rpcs=rpcs, 346 src_crs=src_crs, src_nodata=src_nodata, dst_transform=dst_transform, 347 dst_crs=dst_crs, dst_nodata=dst_nodata, dst_alpha=dst_alpha, 348 src_alpha=src_alpha, resampling=resampling, 349 init_dest_nodata=init_dest_nodata, num_threads=num_threads, 350 warp_mem_limit=warp_mem_limit, **kwargs) 352 return destination, dst_transform File rasterio/_warp.pyx:520, in rasterio._warp._reproject() File rasterio/_warp.pyx:492, in rasterio._warp._reproject() File rasterio/_err.pyx:221, in rasterio._err.exc_wrap_pointer() CPLE_AppDefinedError: Unable to compute a BILINEAR based transformation between pixel/line and georeferenced coordinates for MEM:::DATAPOINTER=13593378816,PIXELS=4680,LINES=11820,BANDS=1,DATATYPE=Float32. ```

Something seems to work (until it doesn't...) as the reprojected geobox looks fine:

image

When computing the result into memory it fails as mentioned above. Same problem occurs if the input DataArray test is computed before reprojecting, so it's not only a problem with dask-backed DataArrays/Datasets.

No problem with rioxarray:

image

I found this related discussion in rasterio, which might help: https://github.com/rasterio/rasterio/discussions/2822 Maybe a problem related to caching?

Kirill888 commented 4 months ago

@maawoo thanks for the report, but that's not enough info to debug.

Error from GDAL looks like it didn't get src transform information.

maawoo commented 3 months ago

Sorry for taking so long to respond and not providing enough info in the first place! I just got back to checking this again and realized that it was just a mistake on my side. Not sure why I was using method= instead of resampling=. With the latter it works without a problem!