Reader resampling_method has no effect

[DanSchoppe]

On the latest release of rio-tiler (as well as several earlier 6.x releases I've tried), the resampling_method argument to Reader.tile() seems to have no effect. I haven't tested whether it affects Reader methods beyond tile().

Here's a simple recreation:

from rio_tiler.io import Reader

with Reader('https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/15/T/VK/2023/10/S2B_15TVK_20231008_0_L2A/TCI.tif') as cog:
    img = cog.tile(
        tile_x=31588,
        tile_y=47192,
        tile_z=17,
        resampling_method="cubic",
    )

with open('tile.png', 'wb') as f:
    f.write(img.render())

The resulting tile looks like:

... regardless of the value provided in resampling_method.

Earlier versions of rio-tiler would have produced something like:

[vincentsarago]

@DanSchoppe there are now two resampling options

rio-tiler/rio_tiler/reader.py

Lines 91 to 92 in 775ecfd

	resampling_method: RIOResampling = "nearest",
	reproject_method: WarpResampling = "nearest",

resampling_method: upscaling / downscaling the data
reproject_method: wrap resampling

At zoom 17 there are no downscaling which is why the cubic method has no effect.

[DanSchoppe]

Lovely, thanks for the response @vincentsarago !

@DanSchoppe there are now two resampling options

I saw that in the v5 migration guide. I had interpreted this to mean that reproject_method would inherit the value from resampling_method if not set explicitly, but am not seeing any evidence of that in the documentation, so I must have misunderstood.

I can confirm that adding a reproject_method argument gets me back to my expected result:

At zoom 17 there are no downscaling which is why the cubic method has no effect.

I think I understand. At zoom 17, each tile pixel is finer resolution than the native resolution of the dataset, such that "resampling" is not applicable? I'll have to learn about rio-tiler's definitions of resampling and reprojection so I can make intelligent decisions on when I should be using each.

One last question: I originally bumped into this resampling behavior change via TiTiler. Is it true that there's currently no way to get reproject_method from a TiTiler query string argument all the way through to rio-tiler?

[DanSchoppe]

@vincentsarago Just reading up on rasterio's resampling documentation. Resampling refers to both upsampling and downsampling, and in either case it seems like the resampling argument to .read() should be meaningful. In light of this, I guess it's surprising that cubic resampling doesn't have an effect on my original example.

Is it true that earlier versions (<5) of rio-tiler were passing the resampling_method value to rasterio's resampling and reproject arguments, and that's why earlier versions of TiTiler (like 0.11) would react to changes to the resampling query param in both the upscale and downscale cases?

Put differently, are you certain that this separation of resampling_method and reproject_method is the root cause of the change to the above tile request, and not some other regression? If so, I'm happy to forge ahead with adding reproject to TiTiler and adding a reproject query parameter. Just wanting to make sure I'm understanding this correctly.

[vincentsarago]

@DanSchoppe Hope those example will help a bit more

Something we need to remember is that when we are doing reprojection (e.g tile read) we are using WarpedVRT where we set the resampling

rio-tiler/rio_tiler/reader.py

Lines 136 to 140 in 775ecfd

	vrt_params = {
	"crs": dst_crs,
	"add_alpha": True,
	"resampling": warp_resampling,
	}

which is the WarpResampling https://gdal.org/api/gdalwarp_cpp.html#_CPPv415GDALResampleAlg

Let's check the influence of both reproject_method and resampling_method:

import numpy
from rio_tiler.io import Reader

with Reader('https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/15/T/VK/2023/10/S2B_15TVK_20231008_0_L2A/TCI.tif') as cog:
    img_cubic = cog.tile(
        tile_x=31588,
        tile_y=47192,
        tile_z=17,
        reproject_method="cubic",
    )
    img_near = cog.tile(
        tile_x=31588,
        tile_y=47192,
        tile_z=17,
        reproject_method="nearest",
    )
    print("Reproject cubic vs nearest (WarpedVRT)")
    print(numpy.array_equal(img_cubic.array, img_near.array))
    print()

    img_cubic = cog.tile(
        tile_x=31588,
        tile_y=47192,
        tile_z=17,
        resampling_method="cubic",
    )
    img_near = cog.tile(
        tile_x=31588,
        tile_y=47192,
        tile_z=17,
        resampling_method="nearest",
    )
    print("Resampling cubic vs nearest (WarpedVRT)")
    print(numpy.array_equal(img_cubic.array, img_near.array))
    print()
    
    img_mode = cog.tile(
        tile_x=31588,
        tile_y=47192,
        tile_z=17,
        resampling_method="mode",
    )
    img_near = cog.tile(
        tile_x=31588,
        tile_y=47192,
        tile_z=17,
        resampling_method="nearest",
    )
    print("Resampling mode vs nearest (WarpedVRT)")
    print(numpy.array_equal(img_mode.array, img_near.array))
    print()

Reproject cubic vs nearest (WarpedVRT)
False

Resampling cubic vs nearest (WarpedVRT)
True

Resampling mode vs nearest (WarpedVRT)
True

☝️ When we are passing different values for resampling the data looks the same, which is weird. Clearly the resampling method has no effects!

Let's check when we are not using WarpedVRT (so no reprojection)

import numpy
from rio_tiler.io import Reader

with Reader('https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/15/T/VK/2023/10/S2B_15TVK_20231008_0_L2A/TCI.tif') as cog:
    print("Check 'resampling' when no-projection (no WarpedVRT)")
    bounds = cog.tms.xy_bounds(31588, 47192, 17)
    img_cubic = cog.part(
        bounds,
        bounds_crs=cog.tms.rasterio_crs,
        dst_crs=cog.crs,
        resampling_method="cubic",
        width=512,
        height=512,
    )
    img_near = cog.part(
        bounds,
        bounds_crs=cog.tms.rasterio_crs,
        dst_crs=cog.crs,
        resampling_method="nearest",
        width=512,
        height=512,
    )
    print("Resampling cubic vs nearest")
    print(numpy.array_equal(img_cubic.array, img_near.array))
    print()

Check 'resampling' when no-projection (no WarpedVRT)
Resampling cubic vs nearest
False

[vincentsarago]

Ooo I think I know what's going on 😬

when we are doing reprojection (using WarpedVRT) we are creating a VRT which is aligned with the bounds (This is to make sure we fetch the overviews)

rio-tiler/rio_tiler/reader.py

Lines 360 to 366 in 775ecfd

	vrt_transform, vrt_width, vrt_height = get_vrt_transform(
	src_dst,
	bounds,
	height=height,
	width=width,
	dst_crs=dst_crs,
	)

Because the VRT will be of the same size as the desired output, there will be no downscaling/upscaling happening there

rio-tiler/rio_tiler/reader.py

Line 226 in 775ecfd

resampling=io_resampling,

which is why the resampling_method will have no effect 😭

[vincentsarago]

I think I have a fix. BUT I will have to do extensive tests 😓

[DanSchoppe]

Ok. I'm not sure what to say other than I really appreciate you and your effort on this 😬

Happy to contribute some testing to a feature branch when you're ready. Also happy to keep pushing on developmentseed/titiler#714 if you think it's a useful contribution, regardless of whether it turns out to be helpful for this particular issue.

[vincentsarago]

Ok. I'm not sure what to say other than I really appreciate you and your effort on this 😬

🙏

let's wait a bit in titiler and first fix things here

[DanSchoppe]

For the benefit of anyone finding this in the future: this conversation picks up in #647 and is resolved in #648