VRTWarpedDataset::IRasterIO(): optimize I/O when requesting whole ima…

…ge at a resolution that doesn't match an overview

Fixes rasterio/rasterio#3203

With that optimization:

```
$ gdalwarp /vsicurl/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 in.vrt

$ CPL_DEBUG=VSICURL GDAL_INGESTED_BYTES_AT_OPEN=32768 GDAL_HTTP_MULTIRANGE=YES GDAL_HTTP_MERGE_CONSECUTIVE_RANGES=YES GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR python -c "from osgeo import gdal; ds = gdal.Open('in.vrt'); ds.ReadRaster(0, 0, ds.RasterXSize, ds.RasterYSize, ds.RasterXSize // 5, ds.RasterYSize // 5)"
/home/even/gdal/gdal/build_cmake/swig/python/osgeo/gdal.py:311: FutureWarning: Neither gdal.UseExceptions() nor gdal.DontUseExceptions() has been explicitly called. In GDAL 4.0, exceptions will be enabled by default.
  warnings.warn(
VSICURL: GetFileSize(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)=200959614  response_code=200
VSICURL: Downloading 0-32767 (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)...
VSICURL: Got response_code=206
VSICURL: GetFileList(/vsicurl/https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/15/T/VK/2023/10/S2B_15TVK_20231008_0_L2A)
VSICURL: GetFileSize(https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/15/T/VK/2023/10/S2B_15TVK_20231008_0_L2A/TCI.xml)=0  response_code=404
VSICURL: GetFileSize(https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/15/T/VK/2023/10/S2B_15TVK_20231008_0_L2A/TCI.XML)=0  response_code=404
VSICURL: Downloading 3096576-13139967 (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)...
VSICURL: Got response_code=206
```

autotest/gdrivers/vrtwarp.py

@@ -734,8 +734,15 @@ def test_vrtwarp_irasterio_optim_three_band():
     assert warped_vrt_ds.ReadRaster(buf_type=gdal.GDT_UInt16) == expected_data
 
     with gdaltest.config_option("GDAL_VRT_WARP_USE_DATASET_RASTERIO", "NO"):
-        expected_data = warped_vrt_ds.ReadRaster(buf_xsize=20, buf_ysize=20)
-    assert warped_vrt_ds.ReadRaster(buf_xsize=20, buf_ysize=20) == expected_data
+        expected_data = warped_vrt_ds.ReadRaster(buf_xsize=20, buf_ysize=40)
+    assert warped_vrt_ds.ReadRaster(buf_xsize=20, buf_ysize=40) == expected_data
+
+    with gdaltest.config_option("GDAL_VRT_WARP_USE_DATASET_RASTERIO", "NO"):
+        expected_data = warped_vrt_ds.ReadRaster(1, 2, 3, 4, buf_xsize=20, buf_ysize=40)
+    assert (
+        warped_vrt_ds.ReadRaster(1, 2, 3, 4, buf_xsize=20, buf_ysize=40)
+        == expected_data
+    )
 
 
 ###############################################################################

frmts/vrt/vrtwarped.cpp

@@ -1903,10 +1903,14 @@ CPLErr VRTWarpedDataset::IRasterIO(
     int nBandCount, BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
     GSpacing nLineSpace, GSpacing nBandSpace, GDALRasterIOExtraArg *psExtraArg)
 {
+    const bool bWholeImage = nXOff == 0 && nYOff == 0 &&
+                             nXSize == nRasterXSize && nYSize == nRasterYSize;
+
     if (eRWFlag == GF_Write ||
         // For too small request fall back to the block-based approach to
         // benefit from caching
-        nBufXSize <= m_nBlockXSize || nBufYSize <= m_nBlockYSize ||
+        (!bWholeImage &&
+         (nBufXSize <= m_nBlockXSize || nBufYSize <= m_nBlockYSize)) ||
         // Or if we don't request all bands at once
         nBandCount < nBands ||
         !CPLTestBool(
@@ -1933,23 +1937,74 @@ CPLErr VRTWarpedDataset::IRasterIO(
         }
     }
 
-    // Fallback to default block-based implementation when not requesting at
-    // the nominal resolution (in theory we could construct a warper taking
-    // into account that, like we do for virtual warped overviews, but that
-    // would be a complication).
-    if (nBufXSize != nXSize || nBufYSize != nYSize)
-    {
-        return GDALDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
-                                      pData, nBufXSize, nBufYSize, eBufType,
-                                      nBandCount, panBandMap, nPixelSpace,
-                                      nLineSpace, nBandSpace, psExtraArg);
-    }
-
     if (m_poWarper == nullptr)
         return CE_Failure;
 
     const GDALWarpOptions *psWO = m_poWarper->GetOptions();
 
+    if (nBufXSize != nXSize || nBufYSize != nYSize)
+    {
+        if (!bWholeImage || !GDALTransformHasFastClone(psWO->pTransformerArg))
+        {
+            return GDALDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
+                                          pData, nBufXSize, nBufYSize, eBufType,
+                                          nBandCount, panBandMap, nPixelSpace,
+                                          nLineSpace, nBandSpace, psExtraArg);
+        }
+
+        // Build a temporary dataset taking into account the rescaling
+        void *pTransformerArg = GDALCloneTransformer(psWO->pTransformerArg);
+        if (pTransformerArg == nullptr)
+        {
+            return GDALDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
+                                          pData, nBufXSize, nBufYSize, eBufType,
+                                          nBandCount, panBandMap, nPixelSpace,
+                                          nLineSpace, nBandSpace, psExtraArg);
+        }
+
+        GDALWarpOptions *psRescaledWO = GDALCloneWarpOptions(psWO);
+        psRescaledWO->hSrcDS = psWO->hSrcDS;
+        psRescaledWO->pfnTransformer = psWO->pfnTransformer;
+        psRescaledWO->pTransformerArg = pTransformerArg;
+
+        // Rescale the output geotransform on the transformer.
+        double adfDstGeoTransform[6] = {0.0};
+        GDALGetTransformerDstGeoTransform(psRescaledWO->pTransformerArg,
+                                          adfDstGeoTransform);
+        RescaleDstGeoTransform(adfDstGeoTransform, nRasterXSize, nBufXSize,
+                               nRasterYSize, nBufYSize);
+        GDALSetTransformerDstGeoTransform(psRescaledWO->pTransformerArg,
+                                          adfDstGeoTransform);
+
+        GDALDatasetH hDstDS =
+            GDALCreateWarpedVRT(psWO->hSrcDS, nBufXSize, nBufYSize,
+                                adfDstGeoTransform, psRescaledWO);
+
+        GDALDestroyWarpOptions(psRescaledWO);
+
+        if (hDstDS == nullptr)
+        {
+            GDALDestroyTransformer(pTransformerArg);
+            return GDALDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
+                                          pData, nBufXSize, nBufYSize, eBufType,
+                                          nBandCount, panBandMap, nPixelSpace,
+                                          nLineSpace, nBandSpace, psExtraArg);
+        }
+
+        auto poOvrDS = static_cast<VRTWarpedDataset *>(hDstDS);
+        poOvrDS->m_bIsOverview = true;
+
+        GDALRasterIOExtraArg sExtraArg;
+        INIT_RASTERIO_EXTRA_ARG(sExtraArg);
+        CPLErr eErr = poOvrDS->IRasterIO(GF_Read, 0, 0, nBufXSize, nBufYSize,
+                                         pData, nBufXSize, nBufYSize, eBufType,
+                                         nBandCount, panBandMap, nPixelSpace,
+                                         nLineSpace, nBandSpace, &sExtraArg);
+
+        poOvrDS->ReleaseRef();
+        return eErr;
+    }
+
     // Build a map from warped output bands to their index
     std::map<int, int> oMapBandToWarpingBandIndex;
     bool bAllBandsIncreasingOrder =