Rasterio: Boundless reads now throw NotGeoreferencedWarning and have degraded performance

[dionhaefner]

Consider the following example of reading a birds-eye view of a raster dataset into a fixed tile size:

import time
import os
import tempfile

import affine
import numpy as np
import rasterio
import rasterio.windows
from rasterio.enums import Resampling
import matplotlib.pyplot as plt

raster_data = np.arange(-128 * 256, 128 * 256, dtype='int16').reshape(256, 256)

profile = {
    'driver': 'GTiff',
    'dtype': 'int16',
    'nodata': 10000,
    'width': raster_data.shape[1],
    'height': raster_data.shape[0],
    'count': 1,
    'crs': {'init': 'epsg:32637'},
    'transform': affine.Affine(
        2.0, 0.0, 694920.0,
        0.0, -2.0, 2055666.0
    )
}

raster_file = tempfile.NamedTemporaryFile(delete=False)
raster_file.close()

with rasterio.open(raster_file.name, 'w', **profile) as dst:
    dst.write(raster_data, 1)


# start testing
with rasterio.open(raster_file.name) as src:
    w, s, e, n = src.bounds
    window = rasterio.windows.from_bounds(w - 2000, s - 2000, e + 2000, n + 2000,
                                          transform=src.transform)

    for resampling in ('nearest', 'bilinear', 'average', 'cubic'):
        resampling_enum = getattr(Resampling, resampling)
        start = time.clock()
        out = src.read(1, window=window, boundless=True, resampling=resampling_enum,
                       fill_value=10000, out_shape=(256, 256))
        end = time.clock()
        print(f'Boundless read took {end-start:.3f}s ({resampling})')

plt.figure()
plt.imshow(raster_data)

plt.figure()
plt.imshow(out)
plt.show()

os.remove(raster_file.name)

With Rasterio 1.0.6, this works nicely and prints

$ python rio_bug.py
Boundless read took 0.021s (nearest)
Boundless read took 0.032s (bilinear)
Boundless read took 0.028s (average)
Boundless read took 0.024s (cubic)

However, from Rasterio 1.0.7 I get

$ python rio_bug.py
ERROR 4: No such file or directory
rio_bug.py:45: NotGeoreferencedWarning: Dataset has no geotransform set. The identity matrix may be returned.
  fill_value=10000, out_shape=(256, 256))
Boundless read took 0.015s (nearest)
ERROR 4: No such file or directory
Boundless read took 0.171s (bilinear)
ERROR 4: No such file or directory
Boundless read took 0.166s (average)
ERROR 4: No such file or directory
Boundless read took 0.178s (cubic)

Note the extra warnings and ~5x worse performance for any resampling methods that is not nearest.

Platform

Python 3.6, Windows 64 bit, Rasterio installed from conda-forge

[jpprins1]

Might be related to: https://github.com/mapbox/rasterio/issues/1478 ?
I'm only using nearest.

[dionhaefner]

Yes, I'm pretty sure that is related. It's nice and fast if I remove the fill_value argument, so I assume the fix for #1478 made something worse when using non-default resampling.

[sgillies]

@dionhaefner can you try upgrading to 1.0.8?

[dionhaefner]

I was on 1.0.8 when I discovered the bug. Same effect:

ERROR 4: No such file or directory
rio_bug.py:45: NotGeoreferencedWarning: Dataset has no geotransform set. The identity matrix may be returned.
  fill_value=10000, out_shape=(256, 256))
Boundless read took 0.023s (nearest)
ERROR 4: No such file or directory
Boundless read took 0.205s (bilinear)
ERROR 4: No such file or directory
Boundless read took 0.153s (average)
ERROR 4: No such file or directory
Boundless read took 0.214s (cubic)

[sgillies]

Ah, ok, I read 1.0.7 above. I'll look into this. If you've added overviews to your dataset, it's possible that they are being missed by the VRT-using code.

[dionhaefner]

If you've added overviews to your dataset

Of course, as a law-abiding citizen I only use cloud-optimized GeoTIFF for my geospatial needs 💪

Edit: Actually, in the test case I didn't use overviews (see above).

[sgillies]

@dionhaefner thanks again for the script. It works perfectly, by which I mean it exposes the problem.

In version 1.0.7 we added an in-memory file to represent the fill value, the background: https://github.com/mapbox/rasterio/compare/1.0.6...1.0.7. The warning arises from the lack of a geo transform, which I am fixing.

The poor performance I am not sure about, but it's related to this new file. I don't have a bead on it quite yet.

The ERROR 4 comes from the fact that most of your code executes outside of any GDAL error handler. Inside of a with rasterio.Env(): block, GDAL errors are dispatched to Python's logging facility instead of directly to stderr.

[sgillies]

I see more like a 9x slow down, except for nearest which is faster in 1.0.10.

1.0.10

$ python issue1499.py 
Boundless read took 0.002s (nearest)
Boundless read took 0.091s (bilinear)
Boundless read took 0.083s (average)
Boundless read took 0.095s (cubic)

1.0.6

$ python issue1499.py 
Boundless read took 0.016s (nearest)
Boundless read took 0.010s (bilinear)
Boundless read took 0.009s (average)
Boundless read took 0.008s (cubic)

[sgillies]

I'm a little stumped by this one. I think our current implementation is much more correct than the old one but clearly there's a slowdown for resampling reads when we add the background layer to the VRT. I'm going to have to let the performance issue slide to 1.0.12. The needless warning will be fixed in 1.0.11.

[dionhaefner]

Thanks for the update. I just used a VRT for now (VRT behavior seems a bit more predictable), but it would be nice to have for cases where you don't want to change the CRS.

[sgillies]

@dionhaefner can you help me understand how you get a fill value with your VRT? I ask because Rasterio uses a VRT to implement boundless reads and we're taking the approach that Frank Warmerdam suggests at https://trac.osgeo.org/gdal/ticket/3327#comment:5. It turns out that this approach is slow when you resample because time is wasted resampling the background layer.

[sgillies]

I've taken this issue upstream.

[dionhaefner]

We currently use add_alpha=True, and then read that alpha band:

https://github.com/DHI-GRAS/terracotta/blob/master/terracotta/drivers/raster_base.py#L437-L532

It's not a pretty solution and not as fast as it could be, but at least it seems somewhat stable (consistent performance and results).

[sgillies]

Got a killer optimization tip from Even that I'm testing in PR #1571. It seems to completely reverse the regression you reported @dionhaefner .