In rasterio.warp.reproject, if the source ndarray (b) is indexing another ndarray (a), the reprojection should operate on b. It appears that is operates on a.
a = np.zeros((5, 5))
b = a[::2, ::2]
reproject(
b,
destination,
src_transform=src_transform,
src_crs=src_crs,
dst_transform=dst_transform,
dst_crs=dst_crs,
resampling=RESAMPLING.nearest)
Linux version 4.10.0-42-generic (buildd@lgw01-amd64-012) (gcc version 6.3.0 20170406 (Ubuntu 6.3.0-12ubuntu2) )
0.36.0 installed with conda install rasterio
So it looks like the ndarray raw data is passed to GDAL here:
https://github.com/mapbox/rasterio/blob/master/rasterio/_shim1.pyx#L109
But np.PyArray_DATA is the pointer to the base array data (a), so it is normal that b's indexing is not taken into account.
A simple fix would be to make a copy of the source ndarray, or at least to warn the user that array views are not taken into account.
I'm preparing a PR.
Thanks for looking into this @davidbrochart ! I think we likely have bugs like this in a couple other places, so a general solution would be very appreciated.
@davidbrochart I took another look at this issue and your PR this morning and realized that I need to ask you if switching to rasterio==1.0a12 addresses the problem. In recent versions we have switched to using GDAL's GDALRasterIOEx() function for its support of array (or view) striding: https://github.com/mapbox/rasterio/blob/master/rasterio/shim_rasterioex.pxi#L42-L43. I have a test of it at https://github.com/mapbox/rasterio/blob/master/tests/test_warp.py#L284.
Would you be willing to create a new conda environment and grab rasterio==1.0a12 from the conda-forge dev channel?
Yes, version 1.0a12 fixes the problem, no copy is needed anymore.
Thanks, this PR can be discarded.