adjust latitude in WarpedVRT parameter calculation (#660)

.gitignore

@@ -109,3 +109,4 @@ tests/fixtures/mask*
 .vscode/settings.json
 
 docs/src/api/*
+dev_notebooks/*

rio_tiler/utils.py

@@ -21,7 +21,7 @@
 from rasterio.warp import calculate_default_transform, transform_geom
 
 from rio_tiler.colormap import apply_cmap
-from rio_tiler.constants import WEB_MERCATOR_CRS
+from rio_tiler.constants import WEB_MERCATOR_CRS, WGS84_CRS
 from rio_tiler.errors import RioTilerError
 from rio_tiler.types import BBox, ColorMapType, IntervalTuple, RIOResampling
 
@@ -244,9 +244,38 @@ def get_vrt_transform(
 
     """
     if src_dst.crs != dst_crs:
+        src_width = src_dst.width
+        src_height = src_dst.height
+        src_bounds = list(src_dst.bounds)
+
+        # Fix for https://github.com/cogeotiff/rio-tiler/issues/654
+        #
+        # When using `calculate_default_transform` with dataset
+        # which span at high/low latitude outside the area_of_use
+        # of the WebMercator projection, we `crop` the dataset
+        # to get the transform (resolution).
+        #
+        # Note: Should be handled in gdal 3.8 directly
+        # https://github.com/OSGeo/gdal/pull/8775
+        if (
+            src_dst.crs == WGS84_CRS
+            and dst_crs == WEB_MERCATOR_CRS
+            and (src_bounds[1] < -85.06 or src_bounds[3] > 85.06)
+        ):
+            warnings.warn(
+                "Adjusting dataset latitudes to avoid re-projection overflow",
+                UserWarning,
+            )
+            src_bounds[1] = max(src_bounds[1], -85.06)
+            src_bounds[3] = min(src_bounds[3], 85.06)
+            w = windows.from_bounds(*src_bounds, transform=src_dst.transform)
+            src_height = round(w.height)
+            src_width = round(w.width)
+
         dst_transform, _, _ = calculate_default_transform(
-            src_dst.crs, dst_crs, src_dst.width, src_dst.height, *src_dst.bounds
+            src_dst.crs, dst_crs, src_width, src_height, *src_bounds
         )
+
     else:
         dst_transform = src_dst.transform
 

tests/test_utils.py

@@ -8,6 +8,7 @@
 import numpy as np
 import pytest
 import rasterio
+from rasterio.crs import CRS
 from rasterio.enums import ColorInterp
 from rasterio.errors import NotGeoreferencedWarning
 from rasterio.features import bounds as featureBounds
@@ -588,3 +589,36 @@ def test_get_array_statistics_coverage():
     assert stats[0]["max"] == 4
     assert stats[0]["mean"] == 2.5
     assert stats[0]["count"] == 4
+
+
+def test_get_vrt_transform_world_file(dataset_fixture):
+    """Should return correct transform and size."""
+    bounds = (
+        -17811118.526923772,
+        -6446275.841017159,
+        17811118.526923772,
+        6446275.841017159,
+    )
+    with MemoryFile(
+        dataset_fixture(
+            crs=CRS.from_epsg(4326),
+            bounds=(-180.0, -90, 180.0, 90.0),
+            dtype="uint8",
+            nodata_type="nodata",
+            width=720,
+            height=360,
+        )
+    ) as memfile:
+        with memfile.open() as src_dst:
+            # adjusting latitudes
+            # with pytest.warns(UserWarning):
+            vrt_transform, vrt_width, vrt_height = utils.get_vrt_transform(
+                src_dst,
+                bounds,
+                dst_crs="epsg:3857",
+            )
+
+    assert vrt_transform[2] == -17811118.526923772
+    assert vrt_transform[5] == 6446275.841017159
+    assert vrt_width == 501  # 59 without the latitude adjust patch
+    assert vrt_height == 181  # 21 without the latitude adjust patch