Shift/zoom after updating to >=4

[fwfichtner]

Since the update from 3.5.1 to 4.0.1 we are observing an irregular location shift/zoom after translating files to Cloud Optimized GeoTiff. If we change nothing, but the rio-cogeo and morecantile (3.4.0-> 4.3.0) version, the results are completely different.

from pyproj import CRS
from rio_cogeo import cog_profiles, cog_translate
import rasterio

from rasterio.io import MemoryFile
from rasterio.warp import Resampling, calculate_default_transform, reproject

if __name__ == "__main__":
    # random Sentinel-2 file, for simplicity I just used https://github.com/recepsuluker/SatImage/blob/main/sentinel2Image20m.jp2
    p = "sentinel2Image20m.jp2"

    with rasterio.open(p) as src:
        reproj_crs = CRS.from_epsg(3857)

        reproj_transform, reproj_width, reproj_height = calculate_default_transform(
            src.crs,
            reproj_crs,
            src.width,
            src.height,
            *src.bounds,
        )
        reproj_meta = src.meta.copy()
        reproj_meta.update(
            {
                "driver": "GTiff",
                "crs": reproj_crs,
                "transform": reproj_transform,
                "width": reproj_width,
                "height": reproj_height,
                "count": 1,
                "dtype": "float32",
            }
        )

        with MemoryFile() as memfile:
            with memfile.open(**reproj_meta) as reproj_ds:
                print("Reprojecting...")
                reproject(
                    source=src.read(1),
                    destination=rasterio.band(reproj_ds, 1),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=reproj_transform,
                    dst_crs=reproj_crs,
                    resampling=Resampling.nearest,
                )

                dst_meta = reproj_meta.copy()
                dst_meta.update(**cog_profiles.get("deflate"))

                print("Generating COG file...")
                cog_translate(
                    reproj_ds,
                    "s2_cog_351.tif",  # TODO update this name
                    dst_meta,
                    web_optimized=True,
                    use_cog_driver=True,
                    forward_band_tags=True,
                )

The results are completely different (besides the tags I would not have expected any changes):

(base) ➜  cogeo_debug gdalinfo s2_cog_351.tif>351.txt
(base) ➜  cogeo_debug gdalinfo s2_cog_401.tif>401.txt
(base) ➜  cogeo_debug diff 351.txt 401.txt
2,3c2,3
< Files: s2_cog_351.tif
< Size is 5120, 5120
---
> Files: s2_cog_401.tif
> Size is 5539, 5553
49,50c49,50
< Origin = (997961.841291261138394,6555239.545736715197563)
< Pixel Size = (38.218514142588127,-38.218514142588127)
---
> Origin = (1001843.989920884137973,6542478.235257362946868)
> Pixel Size = (31.144913025519230,-31.144913025519230)
53a54,55
>   TILING_SCHEME_NAME=WebMercatorQuad
>   TILING_SCHEME_ZOOM_LEVEL=12
56c58
<   INTERLEAVE=PIXEL
---
>   INTERLEAVE=BAND
58,60d59
< Tiling Scheme:
<   NAME=GOOGLEMAPSCOMPATIBLE
<   ZOOM_LEVEL=11
62,66c61,65
< Upper Left  (  997961.841, 6555239.546) (  8d57'53.44"E, 50d37'30.26"N)
< Lower Left  (  997961.841, 6359560.753) (  8d57'53.44"E, 49d29'48.03"N)
< Upper Right ( 1193640.634, 6555239.546) ( 10d43'21.56"E, 50d37'30.26"N)
< Lower Right ( 1193640.634, 6359560.753) ( 10d43'21.56"E, 49d29'48.03"N)
< Center      ( 1095801.237, 6457400.150) (  9d50'37.50"E, 50d 3'51.09"N)
---
> Upper Left  ( 1001843.990, 6542478.235) (  8d59'58.98"E, 50d33' 8.25"N)
> Lower Left  ( 1001843.990, 6369530.533) (  8d59'58.98"E, 49d33'17.31"N)
> Upper Right ( 1174355.663, 6542478.235) ( 10d32'57.90"E, 50d33' 8.25"N)
> Lower Right ( 1174355.663, 6369530.533) ( 10d32'57.90"E, 49d33'17.31"N)
> Center      ( 1088099.827, 6456004.384) (  9d46'28.44"E, 50d 3'22.11"N)
68,70c67
<   Overviews: 2560x2560, 1280x1280, 640x640, 320x320
< Band 2 Block=512x512 Type=Float32, ColorInterp=Alpha
<   Overviews: 2560x2560, 1280x1280, 640x640, 320x320
---
>   Overviews: 2770x2777, 1385x1389, 693x695, 347x348

There are breaking changes in morecantile, but it is currently unclear to me how to adapt our code in order to reflect them.

Edit: Just after writing this I realized that this is connected to the web_optimized-flag.

[vincentsarago]

@fwfichtner yeah there was a change on the web-optimized part (especially when used in combination with the GDAL COG Driver). This was made to allow more TMS (because with GDAL we were limited to few ones and TMS 1.0 spec)

ref: https://github.com/cogeotiff/rio-cogeo/pull/261

The real question is to know which version give a more precise result

as you noted this case is mentioned in the breaking change in https://github.com/cogeotiff/rio-cogeo/blob/main/CHANGES.md#400-2023-05-31 but we could do a better job explaining what breaks 😓

[fwfichtner]

Got it, thanks! I am still not sure how it broke, but it does make more sense how it is now.