Investigate weird output resolution after #648

[vincentsarago]

with rio-tiler 6.2.4

with rio-tiler 6.2.3

data: cog_nasa.tif.zip

Pretty sure this is due to the change in #648

[vincentsarago]

Ok I think I know what's going on!

The file has bounds: BoundingBox(left=-180.0, bottom=-90.0, right=180.0, top=90.0), when we want tiles in WebMercator we call rasterio.warp.calculate_default_transform to get the transform of the reprojected input file... but our file bounds is much bigger than the webmercator bounds so the transform is 💩

rio-tiler/rio_tiler/utils.py

Lines 247 to 249 in dc7e87d

	dst_transform, _, _ = calculate_default_transform(
	src_dst.crs, dst_crs, src_dst.width, src_dst.height, *src_dst.bounds
	)

files looks good in EPSG:4326 projection

import rasterio
from rasterio.warp import calculate_default_transform
from rasterio.warp import transform_bounds

bbox = (-160, -70, 160, 70)
print("bbox read: ", bbox)
b = transform_bounds("epsg:4326", "epsg:3857", *bbox)
print("bbox read (epsg: 3857): ", b)
w, s, e, n = b

print()
with rasterio.open("tests/fixtures/cog_nasa.tif") as src_dst:
    print("dataset bbox :", src_dst.bounds)
    dst_transform, _, _ = calculate_default_transform(
        src_dst.crs, "epsg:3857", src_dst.width, src_dst.height, *src_dst.bounds
    )

    print("res: ", dst_transform.a, dst_transform.e)
    vrt_width = round((e - w) / dst_transform.a)
    vrt_height = round((s - n) / dst_transform.e)
    print("vrt_size: ", vrt_height, vrt_width)

print()
# Validate Theory using a file which do not cross the WebMercator Bounds
with rasterio.open("tests/fixtures/cog_nasa_crop.tif") as src_dst:
    print("dataset bbox :", src_dst.bounds)
    dst_transform, _, _ = calculate_default_transform(
        src_dst.crs, "epsg:3857", src_dst.width, src_dst.height, *src_dst.bounds
    )
    
    print("res: ", dst_transform.a, dst_transform.e)
    vrt_width = round((e - w) / dst_transform.a)
    vrt_height = round((s - n) / dst_transform.e)
    print("vrt_size: ", vrt_height, vrt_width)

bbox read:  (-160, -70, 160, 70)
bbox read (epsg: 3857):  (-17811118.526923772, -11068715.659379492, 17811118.526923772, 11068715.659379492)

dataset bbox : BoundingBox(left=-180.0, bottom=-90.0, right=180.0, top=90.0)
res:  604620.393207253 -604620.393207253
vrt_size:  37 59

dataset bbox : BoundingBox(left=-170.0, bottom=-80.0, right=170.0, top=80.0)
res:  65163.8369659509 -65163.8369659509
vrt_size:  340 547

[vincentsarago]

Note: this is a global issue, not a rio-tiler nor rasterio one.

When using GDAL to warp the data to EPSG:3857 it gives the exact same result than rio-tiler

$ gdalwarp -of VRT -t_srs epsg:3857 cog_nasa.tif cog_nasa.vrt
Creating output file that is 66P x 802L.
Processing cog_nasa.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

$ gdalinfo cog_nasa.vrt
Driver: VRT/Virtual Raster
Files: cog_nasa.vrt
       cog_nasa.tif
Size is 66, 802
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["unnamed",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (-20037508.342789243906736,242528680.943742722272873)
Pixel Size = (604620.393207252956927,-604620.393207252956927)
Metadata:
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (-20037508.343,242528680.944) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Left  (-20037508.343,-242376874.408) (180d 0' 0.00"W, 90d 0' 0.00"S)
Upper Right (19867437.609,242528680.944) (178d28'20.02"E, 90d 0' 0.00"N)
Lower Right (19867437.609,-242376874.408) (178d28'20.02"E, 90d 0' 0.00"S)
Center      (  -85035.367,   75903.268) (  0d45'49.99"W,  0d40'54.60"N)
Band 1 Block=66x128 Type=Float32, ColorInterp=Gray
  Description = ch4_wl
  Overviews: 33x401

This confirms that rio-tiler does the right thing ;-)

[Plantain]

This is hitting us as well.
We store space-projected imagery and access that with rio-tiler to cut WebMercator tiles. 6.2.4+ gives us terrible picture quality.

Reproduction case:
Data: http://static.skysight.io/rio-bug/orig.tiff

from rio_tiler.io import Reader
with Reader("orig.tiff") as src:
    img = src.tile(230,152,8, padding=2)
    f = open("orig.png", "wb")
    f.write(img.render(resampling_method="lanczos"))

6.2.3:

6.2.7: