Rasterio: [Bug] int > None not supported

[dazza-codes]

Objective

Easily convert any numpy.ndarray into a default rasterio.DataSet using a MemoryFile (or something).

Backtrace

  File "/home/joe/miniconda3/envs/project/lib/python3.7/site-packages/rasterio/io.py", line 136, in open
    nodata=nodata, sharing=sharing, **kwargs)
  File "rasterio/_io.pyx", line 1091, in rasterio._io.DatasetWriterBase.__init__
TypeError: '>' not supported between instances of 'int' and 'NoneType'

Reproduce

# script.py

import numpy as np
from rasterio.io import MemoryFile
# from rasterio.profiles import DefaultGTiffProfile

def create_raster(tif_band: np.ndarray, dtype: str = "int32"):

    profile = {
        'interleave': 'band',
        'tiled': True,
        'blockxsize': 256,
        'blockysize': 256,
        'compress': 'lzw',
        'nodata': 0,
        'dtype': dtype
    }

    with MemoryFile() as mem_file:
        with mem_file.open(driver='GTiff', count=1, **profile) as ds:
            # ds.write(tif_band.tobytes())  # same error
            ds.write(tif_band)
        with mem_file.open() as ds:
            return ds

if __name__ == '__main__':

    dtype = "int32"
    shape = (10, 10)
    x = np.linspace(-5, 5, shape[0])
    y = np.linspace(-5, 5, shape[1])
    xy = np.abs(np.sinc(np.outer(x, y))) * 10000
    xy = xy.astype(dtype)

    ds = create_raster(xy, dtype)

Run it

$ python3 script.py 
Traceback (most recent call last):
  File "script.py", line 34, in <module>
    ds = create_raster(xy, dtype)
  File "script.py", line 19, in create_raster
    with mem_file.open(driver='GTiff', count=1, **profile) as ds:
  File "/home/joe/miniconda3/envs/project/lib/python3.7/site-packages/rasterio/env.py", line 394, in wrapper
    return f(*args, **kwds)
  File "/home/joe/miniconda3/envs/project/lib/python3.7/site-packages/rasterio/io.py", line 136, in open
    nodata=nodata, sharing=sharing, **kwargs)
  File "rasterio/_io.pyx", line 1091, in rasterio._io.DatasetWriterBase.__init__
TypeError: '>' not supported between instances of 'int' and 'NoneType'

versions

$ python --version
Python 3.7.3
$ rio --version
1.0.25

references

• https://rasterio.readthedocs.io/en/latest/topics/memory-files.html
• https://rasterio.readthedocs.io/en/latest/topics/profiles.html
• https://stackoverflow.com/questions/43708541/python-typeerror-typeerror-not-supported-between-instances-of-nonetype

[vincentsarago]

@darrenleeweber The error is because you are trying to create a tiled dataset with blocksize == 256 while your data is 10x10.
ref: https://github.com/cogeotiff/rio-cogeo/issues/80#issuecomment-509359041

if you remove

        'tiled': True,
        'blockxsize': 256,
        'blockysize': 256,

it should work.

[vincentsarago]

I think the real issue is why rasterio behve differently than GDAL here

in
https://github.com/mapbox/rasterio/blob/096cb4878fd3490bad52e17835285de8a2aeeebf/rasterio/_io.pyx#L1090

we force the dataset to be bigger than the blocksize while in GDAL this is somehow supported

cc @sgillies

[sgillies]

@vincentsarago you are not quite correct.

@darrenleeweber in your case MemoryFile.open must be called with height and width parameters, otherwise a dataset can't be creaed. That's where the problem begins. Rasterio does fail mysteriously, that is something we can fix.

[dazza-codes]

working script

import numpy as np
import rasterio
from rasterio.io import MemoryFile
from rasterio.transform import Affine

# from rasterio.profiles import DefaultGTiffProfile


def create_raster(
    tif_band: np.ndarray, crs: str = "+proj=latlong", affine_xfm: Affine = None
) -> rasterio.DatasetReader:

    if affine_xfm is None:
        affine_xfm = Affine.identity()

    profile = {
        "driver": "GTiff",
        "height": tif_band.shape[0],
        "width": tif_band.shape[1],
        "count": 1,
        "dtype": tif_band.dtype,
        "crs": crs,
        "transform": affine_xfm,
    }

    with MemoryFile() as mem_file:
        with rasterio.open(mem_file, "w", **profile) as dataset:
            dataset.write(tif_band, 1)

        src = mem_file.open()
    return src


if __name__ == "__main__":

    shape = (10, 10)
    x = np.linspace(-5, 5, shape[0])
    y = np.linspace(-5, 5, shape[1])
    xy = np.abs(np.sinc(np.outer(x, y))) * 10000

    ds = create_raster(xy)

    import ipdb
    ipdb.set_trace()

run it

In : ds
<open DatasetReader name='/vsimem/6f1ee4d8-7bfd-4e06-8aa5-a3beeb590cf8.' mode='r'>
In : pp(ds.profile)
{'count': 1,
 'crs': CRS.from_wkt('GEOGCS["GRS 1980(IUGG, 1980)",DATUM["unknown",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]'),
 'driver': 'GTiff',
 'dtype': 'float64',
 'height': 10,
 'interleave': 'band',
 'nodata': None,
 'tiled': False,
 'transform': Affine(1.0, 0.0, 0.0,
       0.0, 1.0, 0.0),
 'width': 10}
In : xy == ds.read(1)
array([[ True,  True,  True,  True,  True,  True,  True,  True,  True,
         True],
       [ True,  True,  True,  True,  True,  True,  True,  True,  True,
         True],
       [ True,  True,  True,  True,  True,  True,  True,  True,  True,
         True],
       [ True,  True,  True,  True,  True,  True,  True,  True,  True,
         True],
       [ True,  True,  True,  True,  True,  True,  True,  True,  True,
         True],
       [ True,  True,  True,  True,  True,  True,  True,  True,  True,
         True],
       [ True,  True,  True,  True,  True,  True,  True,  True,  True,
         True],
       [ True,  True,  True,  True,  True,  True,  True,  True,  True,
         True],
       [ True,  True,  True,  True,  True,  True,  True,  True,  True,
         True],
       [ True,  True,  True,  True,  True,  True,  True,  True,  True,
         True]])

[sgillies]

This will be fixed in 1.0.26, unscheduled, but maybe early next week.