Reproducible Demo: Writing a broken TIF in rasterio

[8ern4ard]

This issue was first reported with rasterio but the error is just pasted through from gdal, so I closed it there.

I stumbled over this multiple times in various configurations, e.g. the window size, value range and band number seem to be important. As you can see I nest the code in an endless while loop since the exception doesn't occur every time for me.

Notably, the error occurs while writing, not reading. In my real-world example it actually throws an error while reading.

The problem itself disappears when i) not packbits is used, ii) band count of the target is different (e.g. 3 bands), iii) the window size is different, iii) the value range is different (e.g. random * 100) or iv) sometimes just randomly.

I'm pretty sure I saw this problem also with other compression styles, but this is a first example with packbits.

I run Big Sur
rasterio 1.2.6 py39h906574e_2 conda-forge
gdal 3.3.1 py39h0530131_1 conda-forge
numpy 1.21.1 py39h7eed0ac_0 conda-forge

from math import ceil

import numpy as np
import rasterio
from rasterio.crs import CRS
from rasterio.transform import Affine
from rasterio.windows import Window, crop

# SETUP
profile = {
    'driver': 'GTiff',
    'dtype': 'uint16',
    'nodata': 0.0,
    'width': 1857,
    'height': 1221,
    'count': 32,
    'crs': CRS.from_epsg(32633),
    'transform': Affine(10.0, 0.0, 557950.0, 0.0, -10.0, 5417590.0),
    'blockxsize': 512,
    'blockysize': 128,
    'tiled': True,
    'compress': 'packbits',
    'interleave': 'pixel'}

win_sz = 100
windows = [crop(Window(win_sz * w, win_sz * h, win_sz, win_sz), profile['height'], profile['width'])
           for w in range(ceil(profile['width'] / win_sz))
           for h in range(ceil(profile['height'] / win_sz))]

while True:
    
    # INIT
    target = '/Users/bug.tif'
    with rasterio.open(target, 'w', **profile) as trg:
        pass

    # LOOP
    for window in windows:
        with rasterio.open(target, 'r+') as trg:
            data = np.random.rand(profile['count'], window.height, window.width) * 10000
            trg.write(data.astype(profile['dtype']), window=window)

    with rasterio.open(target, 'r') as trg:
        for w in windows:
            test = trg.read(window=w)

ERROR 1: PackBitsDecode:Not enough data for scanline 256
ERROR 1: TIFFReadEncodedTile() failed.
...
Traceback (most recent call last):
File "rasterio/_io.pyx", line 1382, in rasterio._io.DatasetWriterBase.write
File "rasterio/shim_rasterioex.pxi", line 142, in rasterio._shim.io_multi_band
File "rasterio/_err.pyx", line 190, in rasterio._err.exc_wrap_int
rasterio._err.CPLE_AppDefinedError: /Users/bug.tif, band 1: IReadBlock failed at X offset 2, Y offset 2: TIFFReadEncodedTile() failed.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "/Users/Library/Application Support/JetBrains/PyCharmCE2021.2/scratches/bug.py", line 41, in
trg.write(data.astype(profile['dtype']), window=window)
File "rasterio/_io.pyx", line 1384, in rasterio._io.DatasetWriterBase.write
rasterio.errors.RasterioIOError: Read or write failed. /Users/bug.tif, band 1: IReadBlock failed at X offset 2, Y offset 2: TIFFReadEncodedTile() failed.

Process finished with exit code 1

[rouault]

Minimum deterministic reproducer with

from math import ceil

import numpy as np
import rasterio
from rasterio.crs import CRS
from rasterio.transform import Affine
from rasterio.windows import Window, crop

from rasterio._env import set_gdal_config
set_gdal_config('GDAL_CACHEMAX', 1000)


# SETUP
profile = {
    'driver': 'GTiff',
    'dtype': 'uint16',
    'nodata': 0.0,
    'width': 512,
    'height': 768,
    'count': 12,
    'blockxsize': 128,
    'blockysize': 128,
    'tiled': True,
    'compress': 'packbits',
    'interleave': 'pixel'}

win_sz = 100
windows = [crop(Window(win_sz * w, win_sz * h, win_sz, win_sz), profile['height'], profile['width'])
           for w in range(ceil(profile['width'] / win_sz))
           for h in range(ceil(profile['height'] / win_sz))]
np.random.seed(seed=0)

# INIT
target = '/tmp/bug.tif'
with rasterio.open(target, 'w', **profile) as trg:
    pass

# LOOP
for window in windows:
    with rasterio.open(target, 'r+') as trg:
        print(window)
        data = np.random.rand(profile['count'], window.height, window.width) * 10000
        trg.write(data.astype(profile['dtype']), window=window)

output:

Window(col_off=300, row_off=600, width=100, height=100)
Traceback (most recent call last):
  File "rasterio/_io.pyx", line 1606, in rasterio._io.DatasetWriterBase.write
    io_multi_band(self._hds, 1, xoff, yoff, width, height, arr, indexes_arr)
  File "rasterio/_io.pyx", line 174, in rasterio._io.io_multi_band
    return exc_wrap_int(retval)
  File "rasterio/_err.pyx", line 187, in rasterio._err.exc_wrap_int
    raise exc
rasterio._err.CPLE_AppDefinedError: /tmp/bug.tif, band 1: IReadBlock failed at X offset 3, Y offset 4: TIFFReadEncodedTile() failed.

[rouault]

GDAL only somewhat minimal reproducer:

from osgeo import gdal
import array

ds = gdal.GetDriverByName('GTiff').Create('bug.tif', 144*2, 128, 1, options = ['TILED=YES', 'COMPRESS=PACKBITS', 'BLOCKXSIZE=144', 'BLOCKYSIZE=128'])
x = ((144*128)//2) - 645
ds.GetRasterBand(1).WriteRaster(0, 0, 144, 128, b'\x00' * x  + array.array('B', [i % 255 for i in range(144*128-x)]))
ds.GetRasterBand(1).WriteRaster(144, 0, 144, 128,b'\x00' * (x + 8) + array.array('B', [i % 255 for i in range(144*128-(x+8))]))
ds = None
ds = gdal.Open('bug.tif', gdal.GA_Update)
ds.GetRasterBand(1).ReadRaster(144, 0, 144, 128)
ds.GetRasterBand(1).WriteRaster(0, 0, 144, 128, array.array('B', [i % 255 for i in range(144*128)]))
ds.FlushCache()
assert ds.GetRasterBand(1).ReadRaster(144, 0, 144, 128) is not None