Rasterio: WarpedVRT resampling behavior

[ungarj]

WarpedVRT shows unexpected behavior when using other resampling methods than Resampling.nearest.

Expected behavior and actual behavior.

Resampling.bilinear returns the same output than Resampling.nearest. However, when resampling with GDAL, the output image is smoothed.

GDAL bilinear:

Resampling.nearest:

Resampling.bilinear:

Steps to reproduce the problem.

Here is a zip file containing the example data (a CleanTOPO extract) and the script that reproduces all outputs: warpedvrt_resampling.zip

this is the code snipped that reads and resamples the window:

    with rasterio.open(input_file, "r") as src:
        with WarpedVRT(src, dst_crs=tile.crs, resampling=Resampling.bilinear) as vrt:
            bilinear = vrt.read(
                window=vrt.window(*tile.bounds()),
                out_shape=tile.shape(),
                indexes=1
            )

Not sure whether I used all of the parameters correctly. I used the example from the docs.

Operating system

Ubuntu 16.04

Rasterio version and provenance

1.0a11 installed from PyPI using pip 9.0.1, using Python 2.7

[vincentsarago]

👋 @ungarj can you try to add the resampling method to the vrt.read as in https://github.com/mapbox/rio-tiler/blob/f1c3d5f86650158707fad5843d736d3e9edb082f/rio_tiler/utils.py#L134

[ungarj]

That was the trick, thanks! I must have been quite blind yesterday to not try this out.

It was a bit confusing though because the WarpedVRT object also has a resampling attribute.

[ungarj]

One sidenote: adding boundless=True again leads to the nearest resampled output. But I guess, the keyword is deprecated anyway?

[ungarj]

@vincentsarago resampling seems to fall back to nearest when using boundless keyword. Should I open a fresh issue or append info to this issue and reopen?

[vincentsarago]

@ungarj please do

cc @sgillies

[ungarj]

Ok, I did a little bit of digging:

Shouldn't here be the resampling kwarg in the read() function?
https://github.com/mapbox/rasterio/blob/44b9372c13d9696a6b355da6013e9d0a762ff47c/rasterio/_io.pyx#L362

On my local build adding it seems to work. I'm not sure though whether it's the cleanest fix because I'm not familiar enough with rasterio's internals.

Here is what I suspect: without boundless the internal _read() function is invoked with resampling and everything behaves as expected. Having boundless however triggers WarpedVRT._read() but without resampling.

It seems it is intended that the WarpedVRT.resampling attribute the VRT is initialized with should be passed on to its own internal read() function but it isn't.

Not sure how to solve that in a clean way. One option would be the solution mentioned above but this means the resampling attribute would become obsolete because it shall be passed on to read() anyway. The second option would be to overwrite the from DatasetReaderBase inherited read() and use the own resampling attirbute as default. The latter option seems to be better but I'm not sure how to properly implement this.

In any case the documentation should reflect the correct usage whether to use resampling when creating WarpedVRT and if using it with read() which resampling method is actually applied.

I can prepare a PR with the proposed changes if you want.

[ungarj]

@vincentsarago done: https://github.com/mapbox/rasterio/pull/1238

[sgillies]

@ungarj I added a long comment at #1238 and will follow up here, too.

The case of a boundless read from a WarpedVRT hasn't been well tested. We should try to avoid it and instead create instances of WarpedVRT with large enough dimensions to cover all the anticipated reads.

Rasterio's read() gets pixels from the GDAL block cache and can resample those. The resampling attribute of the WarpedVRT affects the resampling used when warping blocks of imagery that go into the cache. For better or worse, these are two different and independent operations. I will work on improving the documentation on this topic.

[ungarj]

@sgillies thanks!

As mentioned in #1238, I'll try to avoid boundless read by initializing the WarpedVRT with large enough dimensions (width, height) so a boundless read won't be required.

[vincentsarago]

👋 here is another demonstration of boundless=True + resampling unexpected behavior

import rasterio

src = rasterio.open('s3://remotepixel/data/image/LC82330572016015LGN00.tif')
data_expected = src.read(1, window=((0,100), (0,100)), resampling=Resampling.bilinear)

# Reading the data with boundless=True, resampling=Resampling.bilinear
data_boundless = src.read(1, window=((-100,100), (-100,100)), boundless=True, resampling=Resampling.bilinear)

# Reading the data withing the image bounds 
data_inBounds = src.read(1, window=((0,100), (0,100)), boundless=True,  resampling=Resampling.bilinear)

data_boundless != data_expected
data_inBounds != data_expected
data_inBounds == data_boundless

https://gist.github.com/vincentsarago/a58c088c7f7de6fd0806a2ba0106e147

cc @sgillies

[sgillies]

This one has been closed in 1.0b4.