From d56b7d89462fa9978c7f93c2491594e86dc33f9f Mon Sep 17 00:00:00 2001
From: vincentsarago <vincent.sarago@gmail.com>
Date: Thu, 16 Nov 2023 13:14:36 +0100
Subject: [PATCH] cheap transform calculation

---
 rio_tiler/utils.py | 38 ++++++++++++++++++++++++++++++++++----
 1 file changed, 34 insertions(+), 4 deletions(-)

diff --git a/rio_tiler/utils.py b/rio_tiler/utils.py
index 281a0b58..265cd6e9 100644
--- a/rio_tiler/utils.py
+++ b/rio_tiler/utils.py
@@ -13,12 +13,13 @@
 from rasterio.dtypes import _gdal_typename
 from rasterio.enums import ColorInterp, MaskFlags, Resampling
 from rasterio.errors import NotGeoreferencedWarning
+from rasterio.features import bounds as featureBounds
 from rasterio.features import is_valid_geom
 from rasterio.io import DatasetReader, DatasetWriter, MemoryFile
 from rasterio.rio.helpers import coords
 from rasterio.transform import from_bounds, rowcol
 from rasterio.vrt import WarpedVRT
-from rasterio.warp import calculate_default_transform, transform_geom
+from rasterio.warp import calculate_default_transform, transform_bounds, transform_geom
 
 from rio_tiler.colormap import apply_cmap
 from rio_tiler.constants import WEB_MERCATOR_CRS
@@ -222,6 +223,33 @@ def get_overview_level(
     return ovr_idx
 
 
+def bbox_to_feature(west: float, south: float, east: float, north: float) -> Dict:
+    """Create a GeoJSON feature from a bbox."""
+    return {
+        "type": "Polygon",
+        "coordinates": [
+            [[west, south], [west, north], [east, north], [east, south], [west, south]]
+        ],
+    }
+
+
+def _cheap_transform_calculation(src, dst_crs: CRS) -> Affine:
+    center = ((src.bounds[1] + src.bounds[3]) / 2, (src.bounds[0] + src.bounds[2]) / 2)
+    xres, yres = src.res[0] / 2, src.res[1] / 2
+    geom = bbox_to_feature(
+        center[0] - xres,
+        center[1] - yres,
+        center[0] + xres,
+        center[1] + xres,
+    )
+    out = transform_geom(src.crs, dst_crs, geom)
+    bbox = featureBounds(out)
+    xres = bbox[2] - bbox[0]
+    yres = bbox[3] - bbox[1]
+    left, _, _, top = transform_bounds(src.crs, dst_crs, *src.bounds)
+    return Affine(xres, 0.0, left, 0.0, -yres, top)
+
+
 def get_vrt_transform(
     src_dst: Union[DatasetReader, DatasetWriter, WarpedVRT],
     bounds: BBox,
@@ -244,9 +272,11 @@ def get_vrt_transform(
 
     """
     if src_dst.crs != dst_crs:
-        dst_transform, _, _ = calculate_default_transform(
-            src_dst.crs, dst_crs, src_dst.width, src_dst.height, *src_dst.bounds
-        )
+        # dst_transform, _, _ = calculate_default_transform(
+        #     src_dst.crs, dst_crs, src_dst.width, src_dst.height, *src_dst.bounds
+        # )
+        dst_transform = _cheap_transform_calculation(src_dst, dst_crs)
+
     else:
         dst_transform = src_dst.transform