Rasterio: Array-to-array resampling fails in Python 2 - long long warning?

[j08lue]

I am trying to use rasterio.warp.reproject to resample data from a ndarray to a ndarray like this

import numpy as np
import rasterio
import rasterio.warp
import affine

def test_resample():

    arr = (np.random.random((20, 30)) * 100).astype('int')
    fill_value = 42
    newarr = np.full((200, 300), fill_value=fill_value, dtype='int')

    src_crs = {'init': 'EPSG:32632'}
    src_transform = affine.Affine(600.0, 0.0, 399960.0, 0.0, -600.0, 6100020.0)
    dst_transform = affine.Affine(60.0, 0.0, 399960.0, 0.0, -60.0, 6100020.0)

    rasterio.warp.reproject(arr, newarr,
        src_transform=src_transform,
        dst_transform=dst_transform,
        src_crs=src_crs,
        dst_crs=src_crs,
        resample=rasterio.warp.Resampling.nearest)

    assert not (newarr == fill_value).all(), 'target array contains interpolated data'
````

It works like a charm in Python 3, but in Python 2, it seems like no data is written to the `destination` array, as shown in [this notebook](https://gist.github.com/j08lue/f577d557dfd63d4af5db562f5abcfee4).

Also only in Python 2, I get the warning

Exception ValueError: "Buffer dtype mismatch, expected 'long' but got 'long long'" in 'rasterio._io.io_auto' ignored
```

My assumption is that the same issue that causes the warning also means that the reference to the destination array is lost and the data is written to some copy of it instead.

I know Python 3 merged long and int, so that is why the warning disappears. But is that also the reason why the reproject function works?

By the way, I get the warning / issue for any dtype I use, so it cannot be related to the source / destination arrays.

[j08lue]

NB: Resampling also fails (in Python 2) if I use a MemFile as destination like this:

from rasterio.io import MemoryFile

profile = dict(
    driver='GTiff', 
    count=1, 
    crs=src_crs, 
    transform=dst_transform,
    width=newarr.shape[0],
    height=newarr.shape[1],
    dtype=newarr.dtype)

with MemoryFile() as memfile:
    with memfile.open(**profile) as dst:
        rasterio.warp.reproject(
            source=arr, 
            destination=rasterio.band(dst, 1),
            src_transform=src_transform,
            src_crs=src_crs,
            resample=rasterio.warp.Resampling.nearest)
        memfile.seek(0)
    with memfile.open() as out:
        newarr = out.read(1)

[citterio]

Also reproduced by following the docs examples.

[sgillies]

Thanks for the update @citterio. I have no time to fix this bug until we get to our first 1.0 beta, I'm busy making sure all the breaking changes are done and fully explained. I can review and help with a bug fix in the meantime. Thanks for your patience!

[sgillies]

@j08lue @citterio would there be any problems on Python 2.7 if you used 'uint8' or 'int32' as a dtype instead of 'int'? I think sticking to the GDAL supported data types is probably the work around here. I'm going to move this to post-1.0.

[j08lue]

would there be any problems on Python 2.7 if you used 'uint8' or 'int32' as a dtype instead of 'int'?

Yes, same thing: I get a warnings about long long and no data is written to the destination array. No matter if I use int, int16, int32, or uint8 for input and output array data type.

[j08lue]

If my data is int64, I get

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-37-884d0ca767db> in <module>()
      4     src_crs=src_crs,
      5     dst_crs=src_crs,
----> 6     resample=rasterio.warp.Resampling.nearest)

C:\Users\josl\AppData\Local\Continuum\Miniconda3\envs\py2\lib\site-packages\rasterio\env.pyc in wrapper(*args, **kwds)
    346         else:
    347             with Env.from_defaults():
--> 348                 return f(*args, **kwds)
    349     return wrapper
    350 

C:\Users\josl\AppData\Local\Continuum\Miniconda3\envs\py2\lib\site-packages\rasterio\env.pyc in wrapper(*args, **kwds)
    520                                 param, full_kwds[param], inequality, version, reason))
    521 
--> 522             return f(*args, **kwds)
    523 
    524         return wrapper

C:\Users\josl\AppData\Local\Continuum\Miniconda3\envs\py2\lib\site-packages\rasterio\warp.pyc in reproject(source, destination, src_transform, gcps, src_crs, src_nodata, dst_transform, dst_crs, dst_nodata, resampling, init_dest_nodata, **kwargs)
    278     _reproject(source, destination, src_transform, gcps, src_crs, src_nodata,
    279                dst_transform, dst_crs, dst_nodata, False, resampling,
--> 280                init_dest_nodata, **kwargs)
    281 
    282 

rasterio/_warp.pyx in rasterio._warp._reproject()

KeyError: 'int64'

[sgillies]

Well, the "KeyError: 'int64'" error is not something we can do anything about except for raise a better exception because GDAL doesn't have 64-bit int types.

>>> arr = (np.random.random((20, 30)) * 100).astype('int')
>>> arr.dtype
dtype('int64')

For what it's worth, I don't see the buffer type mismatch warning, so it might be specific to the Anaconda build.

[j08lue]

I don't see the buffer type mismatch warning

The origin of the warning may actually be a Windows thing (where sizeof(long) == sizeof(int)). Sorry for not mentioning that. Similar to https://github.com/numba/numba/issues/2578.

But the original issue here also occurs on Unix, does it not?

[sgillies]

@j08lue yes. This issue, GDAL's lack of support for 64-bit integer data, affects every operating system and version of GDAL.

[alessioarena]

I also can confirm this error. I'm using rasterio 1.0.11, gdal 2.2.3 with ECW support and python 2.7.

When running the reproject function as array-to-array (but to resample only), and I'm using the logging.DEBUG level I get:

Entering env context: <rasterio.env.Env object at 0x0000000002E70A20>
Starting outermost env
No GDAL environment exists
New GDAL environment <rasterio._env.GDALEnv object at 0x000000000D6EBB70> created
GDAL_DATA found in environment: 'C:\\Temp\\planned-burns-toolkit\\external_libs\\Lib\\site-packages\\osgeo\\gdal-data'.
Started GDALEnv <rasterio._env.GDALEnv object at 0x000000000D6EBB70>.
Got a copy of environment <rasterio._env.GDALEnv object at 0x000000000D6EBB70> options
Entered env context: <rasterio.env.Env object at 0x0000000002E70A20>
PROJ.4 to be imported: '+ellps=GRS80 +no_defs +proj=utm +south +units=m +zone=54'
Set CRS on temp dataset: PROJCS["UTM Zone 54, Southern Hemisphere",GEOGCS["GRS 1980(IUGG, 1980)",DATUM["unknown",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",141],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["Meter",1]]
PROJ.4 to be imported: '+ellps=GRS80 +no_defs +proj=utm +south +units=m +zone=54'
Set CRS on temp dataset: PROJCS["UTM Zone 54, Southern Hemisphere",GEOGCS["GRS 1980(IUGG, 1980)",DATUM["unknown",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",141],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["Meter",1]]
Created temp destination dataset.
Created transformer and options.
Setting NUM_THREADS option: 1
Configured to warp src band 1 to destination band 1
Configured to warp src band 2 to destination band 2
Configured to warp src band 3 to destination band 3
Set transformer options
Chunk and warp window: 0, 0, 800, 800.
Exiting env context: <rasterio.env.Env object at 0x0000000002E70A20>
Cleared existing <rasterio._env.GDALEnv object at 0x000000000D6EBB70> options
Stopped GDALEnv <rasterio._env.GDALEnv object at 0x000000000D6EBB70>.
Exiting outermost env
Exited env context: <rasterio.env.Env object at 0x0000000002E70A20>

In my case the code runs till the end but produces an empty array having correct size, crs and dtype.
However, if I force a RuntimeError before completion I get:

ValueError: Buffer dtype mismatch, expected 'long' but got 'long long'
Exception ValueError: "Buffer dtype mismatch, expected 'long' but got 'long long'" in 'rasterio._io.io_auto' ignored
ValueError: Buffer dtype mismatch, expected 'long' but got 'long long'
Exception ValueError: "Buffer dtype mismatch, expected 'long' but got 'long long'" in 'rasterio._io.io_auto' ignored
 ValueError: Buffer dtype mismatch, expected 'long' but got 'long long'
Exception ValueError: "Buffer dtype mismatch, expected 'long' but got 'long long'" in 'rasterio._io.io_auto' ignored
Traceback (most recent call last):
File "preprocessing_helper.py", line 185, in <module>
execute(args.input_data, outpath=args.outpath, target_res=args.res, resampling_method=args.res_m, ndvi=(args.ndvi_input if args.ndvi else False))
File "preprocessing_helper.py", line 153, in executedata.append(read_img(i, target_res, resampling_method))
File "preprocessing_helper.py", line 31, in read_img
data = resample(img, data, target_res, resampling_method)
File "preprocessing_helper.py", line 98, in resample
raise RuntimeError
RuntimeError

[sgillies]

@alessioarena thanks for the report!

I'm digging into this today. I've started by enabling propagation of exceptions raised in key Cython functions and am now able to reproduce the buffer mismatch. I think it comes down to passing wrongly typed arrays of band indexes and should be a relatively simple fix for 1.0.12.

[sgillies]

We don't have builds for the platforms where this issue bites, but I'm certain #1573 will help if not outright fix the issue.

[alessioarena]

I tested version 1.0.12 today and the bug seems fixed. Thanks a lot for that

[sgillies]

Thanks for the confirmation @alessioarena !