Rasterio: WarpedVRT transform and results are not as expected

[brendan-ward]

I am comparing the prior method of using reproject for cutting out spherical mercator tiles to the new WarpedVRT method, and I'm not seeing the results I expect.

Where:
tile = mercantile.Tile(70, 101, 8)

Using the reproject method I first derive a transform for the output (using the new mercantile.xy_bounds() from mercantile 0.10.0):

transform = rasterio.transform.from_bounds(*mercantile.xy_bounds(*tile), 256, 256)

out = np.empty(shape=(256, 256), dtype=src.profile['dtype'])
reproject(
    rasterio.band(src, src.indexes),
    out,
    dst_crs='EPSG:3857',
    dst_transform=transform,
    src_nodata=src.nodata,
    dst_nodata=src.nodata,
    num_threads=1,
    resampling=Resampling.nearest
)

transform is

| 611.50, 0.00,-9079495.97|
| 0.00,-611.50, 4226661.92|
| 0.00, 0.00, 1.00|

for bounds in Spherical Mercator:
Bbox(left=-9079495.967826376, bottom=4070118.882129065, right=-8922952.933898335, top=4226661.916057106)

Which looks correct: the transform of the extracted tile matches the bounds of the tile.

Whereas when I use the WarpedVRT method derived from the method in the docs:

with rasterio.open('/tmp/test.tif') as src:
    with WarpedVRT(src, dst_crs='EPSG:3857',
        src_nodata=src.nodata,
        dst_nodata=src.nodata,
        resampling=Resampling.nearest) as vrt:

        dst_window = vrt.window(*mercantile.xy_bounds(*tile), boundless=True)
        vrt_data = vrt.read(window=dst_window, out_shape=(1, 256, 256))

        dst_transform = vrt.window_transform(dst_window)
        scaling = Affine.scale(dst_window.num_cols / 256,
                               dst_window.num_rows / 256)
        dst_transform *= scaling

my dst_transform is:

| 612.50, 0.00,-9079636.67|
| 0.00,-612.50, 4226764.45|
| 0.00, 0.00, 1.00|

The output arrays are also not the same between the two approaches, so it isn't just that the calculation of the dst_transform is wrong.

Perhaps I am interpreting the docs wrong? They imply that using the WarpedVRT method should get to the same point: a tile in spherical mercator based on warping and extracting data for the bounds of that tile. But this doesn't seem to be the case?

Is there something I should be doing differently?

[sgillies]

@brendan-ward It's my intent that the results should be the same and that WarpedVRT is the convenient way to do this kind of thing for GDAL sources. I'll look into this for a bit before the 1.0a9 release, but I'm inclined to release with this bug since the fractional transform differences are ~= 1/600.

[sgillies]

@brendan-ward my current hypothesis: reproject() is correct and reads from WarpedVRT would be correct if we had full support for floating point values our our Window class.

[sgillies]

Yes, this was it. I've got a branch that supports Windows with float offsets and height/width and the transforms (from your example above) come out the same. Integer indexes are assumed throughout the code and it'll take a little more work to relax this and get all the tests to pass. I expect to be done tomorrow and then will make the 1.0a9 release.

[sgillies]

Putting this one off post-1.0a9

[mojodna]

@sgillies I'd love to test your patch even if it's partially-baked. I'm running into floating point-related rendering issues that are painful to work around:

(here's an attempt that works in places but still isn't quite right: https://github.com/mojodna/marblecutter/commit/7f2038ac7265822cf933711165642f2da509c6a6)

[sgillies]

@mojodna it's the float-windows branch.

[mojodna]

float-windows is working well for me. It's the difference between

and

I haven't run into anything unexpected. Thanks!

[sgillies]

Great. I've made a lot of progress on the branch today. All the tests are passing except, not surprisingly, the ones where I'm using merge to make a mosaic from single pixels windows.

[sgillies]

@brendan-ward I'm happy to say that #1074 will solve the discrepancy you've reported here.

import rasterio
from rasterio.enums import Resampling
from rasterio.transform import Affine, from_bounds
from rasterio.vrt import WarpedVRT
from rasterio.windows import Window, rowcol, crop


with rasterio.open('tests/data/RGB.byte.tif') as src:
    with WarpedVRT(src, dst_crs='EPSG:3857', src_nodata=src.nodata,
                   dst_nodata=src.nodata, resampling=Resampling.nearest) as vrt:

        tile = mercantile.tile(*vrt.lnglat(), 8)

        dst_window = vrt.window(*mercantile.xy_bounds(*tile), boundless=True)

        vrt_data = vrt.read(window=dst_window, out_shape=(1, 256, 256), boundless=Tru
e)

        dst_transform = vrt.window_transform(dst_window)
        scaling = Affine.scale(dst_window.num_cols / 256,
                               dst_window.num_rows / 256)
        dst_transform_scaled = dst_transform * scaling
        print("Warped VRT transform:")
        print(dst_transform_scaled)

transform = rasterio.transform.from_bounds(*mercantile.xy_bounds(*tile), 256, 256)
print("Transform from bounds:")
print(transform)

Output:

$ python tensixtythree.py
Warped VRT transform:
| 611.50, 0.00,-8766409.90|
| 0.00,-611.50, 2974317.64|
| 0.00, 0.00, 1.00|
Transform from bounds:
| 611.50, 0.00,-8766409.90|
| 0.00,-611.50, 2974317.64|
| 0.00, 0.00, 1.00|