LCOV - code coverage report
Current view: top level - frmts/vrt - vrtwarped.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 701 905 77.5 %
Date: 2025-07-11 10:11:13 Functions: 42 49 85.7 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Virtual GDAL Datasets
       4             :  * Purpose:  Implementation of VRTWarpedRasterBand *and VRTWarpedDataset.
       5             :  * Author:   Frank Warmerdam <warmerdam@pobox.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2004, Frank Warmerdam <warmerdam@pobox.com>
       9             :  * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "cpl_port.h"
      15             : #include "vrtdataset.h"
      16             : 
      17             : #include <stddef.h>
      18             : #include <stdio.h>
      19             : #include <stdlib.h>
      20             : #include <string.h>
      21             : #include <algorithm>
      22             : #include <map>
      23             : 
      24             : // Suppress deprecation warning for GDALOpenVerticalShiftGrid and
      25             : // GDALApplyVerticalShiftGrid
      26             : #ifndef CPL_WARN_DEPRECATED_GDALOpenVerticalShiftGrid
      27             : #define CPL_WARN_DEPRECATED_GDALOpenVerticalShiftGrid(x)
      28             : #define CPL_WARN_DEPRECATED_GDALApplyVerticalShiftGrid(x)
      29             : #endif
      30             : 
      31             : #include "cpl_conv.h"
      32             : #include "cpl_error.h"
      33             : #include "cpl_minixml.h"
      34             : #include "cpl_progress.h"
      35             : #include "cpl_string.h"
      36             : #include "cpl_vsi.h"
      37             : #include "gdal.h"
      38             : #include "gdal_alg.h"
      39             : #include "gdal_alg_priv.h"
      40             : #include "gdal_priv.h"
      41             : #include "gdalwarper.h"
      42             : #include "ogr_geometry.h"
      43             : 
      44             : /************************************************************************/
      45             : /*                      GDALAutoCreateWarpedVRT()                       */
      46             : /************************************************************************/
      47             : 
      48             : /**
      49             :  * Create virtual warped dataset automatically.
      50             :  *
      51             :  * This function will create a warped virtual file representing the
      52             :  * input image warped into the target coordinate system.  A GenImgProj
      53             :  * transformation is created to accomplish any required GCP/Geotransform
      54             :  * warp and reprojection to the target coordinate system.  The output virtual
      55             :  * dataset will be "northup" in the target coordinate system.   The
      56             :  * GDALSuggestedWarpOutput() function is used to determine the bounds and
      57             :  * resolution of the output virtual file which should be large enough to
      58             :  * include all the input image
      59             :  *
      60             :  * If you want to create an alpha band if the source dataset has none, set
      61             :  * psOptionsIn->nDstAlphaBand = GDALGetRasterCount(hSrcDS) + 1.
      62             :  *
      63             :  * Note that the constructed GDALDatasetH will acquire one or more references
      64             :  * to the passed in hSrcDS.  Reference counting semantics on the source
      65             :  * dataset should be honoured.  That is, don't just GDALClose() it unless it
      66             :  * was opened with GDALOpenShared().
      67             :  *
      68             :  * It is possible to "transfer" the ownership of the source dataset
      69             :  * to the warped dataset in the following way:
      70             :  *
      71             :  * \code{.c}
      72             :  *      GDALDatasetH src_ds = GDALOpen("source.tif");
      73             :  *      GDALDatasetH warped_ds = GDALAutoCreateWarpedVRT( src_ds, ... );
      74             :  *      GDALReleaseDataset(src_ds); // src_ds is not "owned" fully by warped_ds.
      75             :  * Do NOT use GDALClose(src_ds) here
      76             :  *      ...
      77             :  *      ...
      78             :  *      GDALReleaseDataset(warped_ds); // or GDALClose(warped_ds);
      79             :  * \endcode
      80             :  *
      81             :  * Traditonal nested calls are also possible of course:
      82             :  *
      83             :  * \code{.c}
      84             :  *      GDALDatasetH src_ds = GDALOpen("source.tif");
      85             :  *      GDALDatasetH warped_ds = GDALAutoCreateWarpedVRT( src_ds, ... );
      86             :  *      ...
      87             :  *      ...
      88             :  *      GDALReleaseDataset(warped_ds); // or GDALClose(warped_ds);
      89             :  *      GDALReleaseDataset(src_ds); // or GDALClose(src_ds);
      90             :  * \endcode
      91             :  *
      92             :  * The returned dataset will have no associated filename for itself.  If you
      93             :  * want to write the virtual dataset description to a file, use the
      94             :  * GDALSetDescription() function (or SetDescription() method) on the dataset
      95             :  * to assign a filename before it is closed.
      96             :  *
      97             :  * @param hSrcDS The source dataset.
      98             :  *
      99             :  * @param pszSrcWKT The coordinate system of the source image.  If NULL, it
     100             :  * will be read from the source image.
     101             :  *
     102             :  * @param pszDstWKT The coordinate system to convert to.  If NULL no change
     103             :  * of coordinate system will take place.
     104             :  *
     105             :  * @param eResampleAlg One of GRA_NearestNeighbour, GRA_Bilinear, GRA_Cubic,
     106             :  * GRA_CubicSpline, GRA_Lanczos, GRA_Average, GRA_RMS or GRA_Mode.
     107             :  * Controls the sampling method used.
     108             :  *
     109             :  * @param dfMaxError Maximum error measured in input pixels that is allowed in
     110             :  * approximating the transformation (0.0 for exact calculations).
     111             :  *
     112             :  * @param psOptionsIn Additional warp options, normally NULL.
     113             :  *
     114             :  * @return NULL on failure, or a new virtual dataset handle on success.
     115             :  */
     116             : 
     117             : GDALDatasetH CPL_STDCALL
     118          25 : GDALAutoCreateWarpedVRT(GDALDatasetH hSrcDS, const char *pszSrcWKT,
     119             :                         const char *pszDstWKT, GDALResampleAlg eResampleAlg,
     120             :                         double dfMaxError, const GDALWarpOptions *psOptionsIn)
     121             : {
     122          25 :     return GDALAutoCreateWarpedVRTEx(hSrcDS, pszSrcWKT, pszDstWKT, eResampleAlg,
     123          25 :                                      dfMaxError, psOptionsIn, nullptr);
     124             : }
     125             : 
     126             : /************************************************************************/
     127             : /*                     GDALAutoCreateWarpedVRTEx()                      */
     128             : /************************************************************************/
     129             : 
     130             : /**
     131             :  * Create virtual warped dataset automatically.
     132             :  *
     133             :  * Compared to GDALAutoCreateWarpedVRT() this function adds one extra
     134             :  * argument: options to be passed to GDALCreateGenImgProjTransformer2().
     135             :  *
     136             :  * @since 3.2
     137             :  */
     138             : 
     139          25 : GDALDatasetH CPL_STDCALL GDALAutoCreateWarpedVRTEx(
     140             :     GDALDatasetH hSrcDS, const char *pszSrcWKT, const char *pszDstWKT,
     141             :     GDALResampleAlg eResampleAlg, double dfMaxError,
     142             :     const GDALWarpOptions *psOptionsIn, CSLConstList papszTransformerOptions)
     143             : {
     144          25 :     VALIDATE_POINTER1(hSrcDS, "GDALAutoCreateWarpedVRT", nullptr);
     145             : 
     146             :     /* -------------------------------------------------------------------- */
     147             :     /*      Populate the warp options.                                      */
     148             :     /* -------------------------------------------------------------------- */
     149          25 :     GDALWarpOptions *psWO = nullptr;
     150          25 :     if (psOptionsIn != nullptr)
     151           1 :         psWO = GDALCloneWarpOptions(psOptionsIn);
     152             :     else
     153          24 :         psWO = GDALCreateWarpOptions();
     154             : 
     155          25 :     psWO->eResampleAlg = eResampleAlg;
     156             : 
     157          25 :     psWO->hSrcDS = hSrcDS;
     158             : 
     159          25 :     GDALWarpInitDefaultBandMapping(psWO, GDALGetRasterCount(hSrcDS));
     160             : 
     161             :     /* -------------------------------------------------------------------- */
     162             :     /*      Setup no data values (if not done in psOptionsIn)               */
     163             :     /* -------------------------------------------------------------------- */
     164          25 :     if (psWO->padfSrcNoDataReal == nullptr &&
     165          25 :         psWO->padfDstNoDataReal == nullptr && psWO->nSrcAlphaBand == 0)
     166             :     {
     167             :         // If none of the provided input nodata values can be represented in the
     168             :         // data type of the corresponding source band, ignore them.
     169          25 :         int nCountInvalidSrcNoDataReal = 0;
     170          55 :         for (int i = 0; i < psWO->nBandCount; i++)
     171             :         {
     172             :             GDALRasterBandH rasterBand =
     173          30 :                 GDALGetRasterBand(psWO->hSrcDS, psWO->panSrcBands[i]);
     174             : 
     175             :             int hasNoDataValue;
     176             :             double noDataValue =
     177          30 :                 GDALGetRasterNoDataValue(rasterBand, &hasNoDataValue);
     178             : 
     179          30 :             if (hasNoDataValue &&
     180           7 :                 !GDALIsValueExactAs(noDataValue,
     181          30 :                                     GDALGetRasterDataType(rasterBand)))
     182             :             {
     183           6 :                 nCountInvalidSrcNoDataReal++;
     184             :             }
     185             :         }
     186             : 
     187          25 :         if (nCountInvalidSrcNoDataReal != psWO->nBandCount)
     188             :         {
     189          46 :             for (int i = 0; i < psWO->nBandCount; i++)
     190             :             {
     191             :                 GDALRasterBandH rasterBand =
     192          24 :                     GDALGetRasterBand(psWO->hSrcDS, psWO->panSrcBands[i]);
     193             : 
     194             :                 int hasNoDataValue;
     195             :                 double noDataValue =
     196          24 :                     GDALGetRasterNoDataValue(rasterBand, &hasNoDataValue);
     197             : 
     198          24 :                 if (hasNoDataValue)
     199             :                 {
     200             :                     // Check if the nodata value is out of range
     201           1 :                     int bClamped = FALSE;
     202           1 :                     int bRounded = FALSE;
     203           1 :                     CPL_IGNORE_RET_VAL(GDALAdjustValueToDataType(
     204             :                         GDALGetRasterDataType(rasterBand), noDataValue,
     205             :                         &bClamped, &bRounded));
     206           1 :                     if (!bClamped)
     207             :                     {
     208           1 :                         GDALWarpInitNoDataReal(psWO, -1e10);
     209           1 :                         if (psWO->padfSrcNoDataReal != nullptr &&
     210           1 :                             psWO->padfDstNoDataReal != nullptr)
     211             :                         {
     212           1 :                             psWO->padfSrcNoDataReal[i] = noDataValue;
     213           1 :                             psWO->padfDstNoDataReal[i] = noDataValue;
     214             :                         }
     215             :                     }
     216             :                 }
     217             :             }
     218             :         }
     219             : 
     220          25 :         if (psWO->padfDstNoDataReal != nullptr)
     221             :         {
     222           1 :             if (CSLFetchNameValue(psWO->papszWarpOptions, "INIT_DEST") ==
     223             :                 nullptr)
     224             :             {
     225           1 :                 psWO->papszWarpOptions = CSLSetNameValue(
     226             :                     psWO->papszWarpOptions, "INIT_DEST", "NO_DATA");
     227             :             }
     228             :         }
     229             :     }
     230             : 
     231             :     /* -------------------------------------------------------------------- */
     232             :     /*      Create the transformer.                                         */
     233             :     /* -------------------------------------------------------------------- */
     234          25 :     psWO->pfnTransformer = GDALGenImgProjTransform;
     235             : 
     236          25 :     char **papszOptions = nullptr;
     237          25 :     if (pszSrcWKT != nullptr)
     238           5 :         papszOptions = CSLSetNameValue(papszOptions, "SRC_SRS", pszSrcWKT);
     239          25 :     if (pszDstWKT != nullptr)
     240           8 :         papszOptions = CSLSetNameValue(papszOptions, "DST_SRS", pszDstWKT);
     241          25 :     papszOptions = CSLMerge(papszOptions, papszTransformerOptions);
     242          25 :     psWO->pTransformerArg =
     243          25 :         GDALCreateGenImgProjTransformer2(psWO->hSrcDS, nullptr, papszOptions);
     244          25 :     CSLDestroy(papszOptions);
     245             : 
     246          25 :     if (psWO->pTransformerArg == nullptr)
     247             :     {
     248           1 :         GDALDestroyWarpOptions(psWO);
     249           1 :         return nullptr;
     250             :     }
     251             : 
     252             :     /* -------------------------------------------------------------------- */
     253             :     /*      Figure out the desired output bounds and resolution.            */
     254             :     /* -------------------------------------------------------------------- */
     255          24 :     double adfDstGeoTransform[6] = {0.0};
     256          24 :     int nDstPixels = 0;
     257          24 :     int nDstLines = 0;
     258          24 :     CPLErr eErr = GDALSuggestedWarpOutput(
     259             :         hSrcDS, psWO->pfnTransformer, psWO->pTransformerArg, adfDstGeoTransform,
     260             :         &nDstPixels, &nDstLines);
     261          24 :     if (eErr != CE_None)
     262             :     {
     263           0 :         GDALDestroyTransformer(psWO->pTransformerArg);
     264           0 :         GDALDestroyWarpOptions(psWO);
     265           0 :         return nullptr;
     266             :     }
     267             : 
     268             :     /* -------------------------------------------------------------------- */
     269             :     /*      Update the transformer to include an output geotransform        */
     270             :     /*      back to pixel/line coordinates.                                 */
     271             :     /*                                                                      */
     272             :     /* -------------------------------------------------------------------- */
     273          24 :     GDALSetGenImgProjTransformerDstGeoTransform(psWO->pTransformerArg,
     274             :                                                 adfDstGeoTransform);
     275             : 
     276             :     /* -------------------------------------------------------------------- */
     277             :     /*      Do we want to apply an approximating transformation?            */
     278             :     /* -------------------------------------------------------------------- */
     279          24 :     if (dfMaxError > 0.0)
     280             :     {
     281           3 :         psWO->pTransformerArg = GDALCreateApproxTransformer(
     282             :             psWO->pfnTransformer, psWO->pTransformerArg, dfMaxError);
     283           3 :         psWO->pfnTransformer = GDALApproxTransform;
     284           3 :         GDALApproxTransformerOwnsSubtransformer(psWO->pTransformerArg, TRUE);
     285             :     }
     286             : 
     287             :     /* -------------------------------------------------------------------- */
     288             :     /*      Create the VRT file.                                            */
     289             :     /* -------------------------------------------------------------------- */
     290          24 :     GDALDatasetH hDstDS = GDALCreateWarpedVRT(hSrcDS, nDstPixels, nDstLines,
     291             :                                               adfDstGeoTransform, psWO);
     292             : 
     293          24 :     GDALDestroyWarpOptions(psWO);
     294             : 
     295          24 :     if (hDstDS != nullptr)
     296             :     {
     297          24 :         if (pszDstWKT != nullptr)
     298           7 :             GDALSetProjection(hDstDS, pszDstWKT);
     299          17 :         else if (pszSrcWKT != nullptr)
     300           0 :             GDALSetProjection(hDstDS, pszSrcWKT);
     301          17 :         else if (GDALGetGCPCount(hSrcDS) > 0)
     302           2 :             GDALSetProjection(hDstDS, GDALGetGCPProjection(hSrcDS));
     303             :         else
     304          15 :             GDALSetProjection(hDstDS, GDALGetProjectionRef(hSrcDS));
     305             :     }
     306             : 
     307          24 :     return hDstDS;
     308             : }
     309             : 
     310             : /************************************************************************/
     311             : /*                        GDALCreateWarpedVRT()                         */
     312             : /************************************************************************/
     313             : 
     314             : /**
     315             :  * Create virtual warped dataset.
     316             :  *
     317             :  * This function will create a warped virtual file representing the
     318             :  * input image warped based on a provided transformation.  Output bounds
     319             :  * and resolution are provided explicitly.
     320             :  *
     321             :  * If you want to create an alpha band if the source dataset has none, set
     322             :  * psOptions->nDstAlphaBand = GDALGetRasterCount(hSrcDS) + 1.
     323             :  *
     324             :  * Note that the constructed GDALDatasetH will acquire one or more references
     325             :  * to the passed in hSrcDS.  Reference counting semantics on the source
     326             :  * dataset should be honoured.  That is, don't just GDALClose() it unless it
     327             :  * was opened with GDALOpenShared().
     328             :  *
     329             :  * It is possible to "transfer" the ownership of the source dataset
     330             :  * to the warped dataset in the following way:
     331             :  *
     332             :  * \code{.c}
     333             :  *      GDALDatasetH src_ds = GDALOpen("source.tif");
     334             :  *      GDALDatasetH warped_ds = GDALAutoCreateWarpedVRT( src_ds, ... );
     335             :  *      GDALReleaseDataset(src_ds); // src_ds is not "owned" fully by warped_ds.
     336             :  * Do NOT use GDALClose(src_ds) here
     337             :  *      ...
     338             :  *      ...
     339             :  *      GDALReleaseDataset(warped_ds); // or GDALClose(warped_ds);
     340             :  * \endcode
     341             :  *
     342             :  * Traditonal nested calls are also possible of course:
     343             :  *
     344             :  * \code{.c}
     345             :  *      GDALDatasetH src_ds = GDALOpen("source.tif");
     346             :  *      GDALDatasetH warped_ds = GDALAutoCreateWarpedVRT( src_ds, ... );
     347             :  *      ...
     348             :  *      ...
     349             :  *      GDALReleaseDataset(warped_ds); // or GDALClose(warped_ds);
     350             :  *      GDALReleaseDataset(src_ds); // or GDALClose(src_ds);
     351             :  * \endcode
     352             :  *
     353             :  * The returned dataset will have no associated filename for itself.  If you
     354             :  * want to write the virtual dataset description to a file, use the
     355             :  * GDALSetDescription() function (or SetDescription() method) on the dataset
     356             :  * to assign a filename before it is closed.
     357             :  *
     358             :  * @param hSrcDS The source dataset.
     359             :  *
     360             :  * @param nPixels Width of the virtual warped dataset to create
     361             :  *
     362             :  * @param nLines Height of the virtual warped dataset to create
     363             :  *
     364             :  * @param padfGeoTransform Geotransform matrix of the virtual warped dataset
     365             :  * to create
     366             :  *
     367             :  * @param psOptions Warp options. Must be different from NULL.
     368             :  *
     369             :  * @return NULL on failure, or a new virtual dataset handle on success.
     370             :  */
     371             : 
     372          46 : GDALDatasetH CPL_STDCALL GDALCreateWarpedVRT(GDALDatasetH hSrcDS, int nPixels,
     373             :                                              int nLines,
     374             :                                              const double *padfGeoTransform,
     375             :                                              GDALWarpOptions *psOptions)
     376             : 
     377             : {
     378          46 :     VALIDATE_POINTER1(hSrcDS, "GDALCreateWarpedVRT", nullptr);
     379          46 :     VALIDATE_POINTER1(psOptions, "GDALCreateWarpedVRT", nullptr);
     380             : 
     381             :     /* -------------------------------------------------------------------- */
     382             :     /*      Create the VRTDataset and populate it with bands.               */
     383             :     /* -------------------------------------------------------------------- */
     384          46 :     VRTWarpedDataset *poDS = new VRTWarpedDataset(nPixels, nLines);
     385             : 
     386             :     // Call this before assigning hDstDS
     387          46 :     GDALWarpResolveWorkingDataType(psOptions);
     388             : 
     389          46 :     psOptions->hDstDS = poDS;
     390          46 :     poDS->SetGeoTransform(GDALGeoTransform(padfGeoTransform));
     391             : 
     392         101 :     for (int i = 0; i < psOptions->nBandCount; i++)
     393             :     {
     394          55 :         int nDstBand = psOptions->panDstBands[i];
     395         110 :         while (poDS->GetRasterCount() < nDstBand)
     396             :         {
     397          55 :             poDS->AddBand(psOptions->eWorkingDataType, nullptr);
     398             :         }
     399             : 
     400             :         VRTWarpedRasterBand *poBand =
     401          55 :             static_cast<VRTWarpedRasterBand *>(poDS->GetRasterBand(nDstBand));
     402             :         GDALRasterBand *poSrcBand = static_cast<GDALRasterBand *>(
     403          55 :             GDALGetRasterBand(hSrcDS, psOptions->panSrcBands[i]));
     404             : 
     405          55 :         poBand->CopyCommonInfoFrom(poSrcBand);
     406             :     }
     407             : 
     408          47 :     while (poDS->GetRasterCount() < psOptions->nDstAlphaBand)
     409             :     {
     410           1 :         poDS->AddBand(psOptions->eWorkingDataType, nullptr);
     411             :     }
     412          46 :     if (psOptions->nDstAlphaBand)
     413             :     {
     414           1 :         poDS->GetRasterBand(psOptions->nDstAlphaBand)
     415           1 :             ->SetColorInterpretation(GCI_AlphaBand);
     416             :     }
     417             : 
     418             :     /* -------------------------------------------------------------------- */
     419             :     /*      Initialize the warp on the VRTWarpedDataset.                    */
     420             :     /* -------------------------------------------------------------------- */
     421          46 :     const CPLErr eErr = poDS->Initialize(psOptions);
     422          46 :     if (eErr == CE_Failure)
     423             :     {
     424           0 :         psOptions->hDstDS = nullptr;
     425           0 :         delete poDS;
     426           0 :         return nullptr;
     427             :     }
     428             : 
     429          46 :     return poDS;
     430             : }
     431             : 
     432             : /*! @cond Doxygen_Suppress */
     433             : 
     434             : /************************************************************************/
     435             : /* ==================================================================== */
     436             : /*                          VRTWarpedDataset                            */
     437             : /* ==================================================================== */
     438             : /************************************************************************/
     439             : 
     440             : /************************************************************************/
     441             : /*                          VRTWarpedDataset()                          */
     442             : /************************************************************************/
     443             : 
     444         358 : VRTWarpedDataset::VRTWarpedDataset(int nXSize, int nYSize, int nBlockXSize,
     445         358 :                                    int nBlockYSize)
     446             :     : VRTDataset(nXSize, nYSize,
     447         357 :                  nBlockXSize > 0 ? nBlockXSize : std::min(nXSize, 512),
     448         357 :                  nBlockYSize > 0 ? nBlockYSize : std::min(nYSize, 128)),
     449        1072 :       m_poWarper(nullptr), m_nSrcOvrLevel(-2)
     450             : {
     451         358 :     eAccess = GA_Update;
     452         358 :     DisableReadWriteMutex();
     453         358 : }
     454             : 
     455             : /************************************************************************/
     456             : /*                         ~VRTWarpedDataset()                          */
     457             : /************************************************************************/
     458             : 
     459         716 : VRTWarpedDataset::~VRTWarpedDataset()
     460             : 
     461             : {
     462         358 :     VRTWarpedDataset::FlushCache(true);
     463         358 :     VRTWarpedDataset::CloseDependentDatasets();
     464         716 : }
     465             : 
     466             : /************************************************************************/
     467             : /*                        CloseDependentDatasets()                      */
     468             : /************************************************************************/
     469             : 
     470         370 : int VRTWarpedDataset::CloseDependentDatasets()
     471             : {
     472         370 :     bool bHasDroppedRef = CPL_TO_BOOL(VRTDataset::CloseDependentDatasets());
     473             : 
     474             :     /* -------------------------------------------------------------------- */
     475             :     /*      Cleanup overviews.                                              */
     476             :     /* -------------------------------------------------------------------- */
     477         394 :     for (auto &poDS : m_apoOverviews)
     478             :     {
     479          24 :         if (poDS && poDS->Release())
     480             :         {
     481           0 :             bHasDroppedRef = true;
     482             :         }
     483             :     }
     484             : 
     485         370 :     m_apoOverviews.clear();
     486             : 
     487             :     /* -------------------------------------------------------------------- */
     488             :     /*      Cleanup warper if one is in effect.                             */
     489             :     /* -------------------------------------------------------------------- */
     490         370 :     if (m_poWarper != nullptr)
     491             :     {
     492         356 :         const GDALWarpOptions *psWO = m_poWarper->GetOptions();
     493             : 
     494             :         /* --------------------------------------------------------------------
     495             :          */
     496             :         /*      We take care to only call GDALClose() on psWO->hSrcDS if the */
     497             :         /*      reference count drops to zero.  This is makes it so that we */
     498             :         /*      can operate reference counting semantics more-or-less */
     499             :         /*      properly even if the dataset isn't open in shared mode, */
     500             :         /*      though we require that the caller also honour the reference */
     501             :         /*      counting semantics even though it isn't a shared dataset. */
     502             :         /* --------------------------------------------------------------------
     503             :          */
     504         356 :         if (psWO != nullptr && psWO->hSrcDS != nullptr)
     505             :         {
     506         355 :             if (GDALReleaseDataset(psWO->hSrcDS))
     507             :             {
     508         278 :                 bHasDroppedRef = true;
     509             :             }
     510             :         }
     511             : 
     512             :         /* --------------------------------------------------------------------
     513             :          */
     514             :         /*      We are responsible for cleaning up the transformer ourselves. */
     515             :         /* --------------------------------------------------------------------
     516             :          */
     517         356 :         if (psWO != nullptr && psWO->pTransformerArg != nullptr)
     518         355 :             GDALDestroyTransformer(psWO->pTransformerArg);
     519             : 
     520         356 :         delete m_poWarper;
     521         356 :         m_poWarper = nullptr;
     522             :     }
     523             : 
     524             :     /* -------------------------------------------------------------------- */
     525             :     /*      Destroy the raster bands if they exist.                         */
     526             :     /* -------------------------------------------------------------------- */
     527        1053 :     for (int iBand = 0; iBand < nBands; iBand++)
     528             :     {
     529         683 :         delete papoBands[iBand];
     530             :     }
     531         370 :     nBands = 0;
     532             : 
     533         370 :     return bHasDroppedRef;
     534             : }
     535             : 
     536             : /************************************************************************/
     537             : /*                         SetSrcOverviewLevel()                        */
     538             : /************************************************************************/
     539             : 
     540          97 : CPLErr VRTWarpedDataset::SetMetadataItem(const char *pszName,
     541             :                                          const char *pszValue,
     542             :                                          const char *pszDomain)
     543             : 
     544             : {
     545          97 :     if ((pszDomain == nullptr || EQUAL(pszDomain, "")) &&
     546          97 :         EQUAL(pszName, "SrcOvrLevel"))
     547             :     {
     548          97 :         const int nOldValue = m_nSrcOvrLevel;
     549          97 :         if (pszValue == nullptr || EQUAL(pszValue, "AUTO"))
     550           1 :             m_nSrcOvrLevel = -2;
     551          96 :         else if (STARTS_WITH_CI(pszValue, "AUTO-"))
     552           1 :             m_nSrcOvrLevel = -2 - atoi(pszValue + 5);
     553          95 :         else if (EQUAL(pszValue, "NONE"))
     554           1 :             m_nSrcOvrLevel = -1;
     555          94 :         else if (CPLGetValueType(pszValue) == CPL_VALUE_INTEGER)
     556          94 :             m_nSrcOvrLevel = atoi(pszValue);
     557          97 :         if (m_nSrcOvrLevel != nOldValue)
     558           2 :             SetNeedsFlush();
     559          97 :         return CE_None;
     560             :     }
     561           0 :     return VRTDataset::SetMetadataItem(pszName, pszValue, pszDomain);
     562             : }
     563             : 
     564             : /************************************************************************/
     565             : /*                        VRTWarpedAddOptions()                         */
     566             : /************************************************************************/
     567             : 
     568         356 : static char **VRTWarpedAddOptions(char **papszWarpOptions)
     569             : {
     570             :     /* Avoid errors when adding an alpha band, but source dataset has */
     571             :     /* no alpha band (#4571), and generally don't leave our buffer uninitialized
     572             :      */
     573         356 :     if (CSLFetchNameValue(papszWarpOptions, "INIT_DEST") == nullptr)
     574          52 :         papszWarpOptions = CSLSetNameValue(papszWarpOptions, "INIT_DEST", "0");
     575             : 
     576             :     /* For https://github.com/OSGeo/gdal/issues/1985 */
     577         356 :     if (CSLFetchNameValue(papszWarpOptions,
     578         356 :                           "ERROR_OUT_IF_EMPTY_SOURCE_WINDOW") == nullptr)
     579             :     {
     580         226 :         papszWarpOptions = CSLSetNameValue(
     581             :             papszWarpOptions, "ERROR_OUT_IF_EMPTY_SOURCE_WINDOW", "FALSE");
     582             :     }
     583         356 :     return papszWarpOptions;
     584             : }
     585             : 
     586             : /************************************************************************/
     587             : /*                             Initialize()                             */
     588             : /*                                                                      */
     589             : /*      Initialize a dataset from passed in warp options.               */
     590             : /************************************************************************/
     591             : 
     592         157 : CPLErr VRTWarpedDataset::Initialize(void *psWO)
     593             : 
     594             : {
     595         157 :     if (m_poWarper != nullptr)
     596           0 :         delete m_poWarper;
     597             : 
     598         157 :     m_poWarper = new GDALWarpOperation();
     599             : 
     600             :     GDALWarpOptions *psWO_Dup =
     601         157 :         GDALCloneWarpOptions(static_cast<GDALWarpOptions *>(psWO));
     602             : 
     603         157 :     psWO_Dup->papszWarpOptions =
     604         157 :         VRTWarpedAddOptions(psWO_Dup->papszWarpOptions);
     605             : 
     606         157 :     CPLErr eErr = m_poWarper->Initialize(psWO_Dup);
     607             : 
     608             :     // The act of initializing this warped dataset with this warp options
     609             :     // will result in our assuming ownership of a reference to the
     610             :     // hSrcDS.
     611             : 
     612         157 :     if (eErr == CE_None &&
     613         156 :         static_cast<GDALWarpOptions *>(psWO)->hSrcDS != nullptr)
     614             :     {
     615         156 :         GDALReferenceDataset(psWO_Dup->hSrcDS);
     616             :     }
     617             : 
     618         157 :     GDALDestroyWarpOptions(psWO_Dup);
     619             : 
     620         157 :     if (nBands > 1)
     621             :     {
     622          35 :         GDALDataset::SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
     623             :     }
     624             : 
     625         157 :     return eErr;
     626             : }
     627             : 
     628             : /************************************************************************/
     629             : /*                        GDALWarpCoordRescaler                         */
     630             : /************************************************************************/
     631             : 
     632           1 : class GDALWarpCoordRescaler : public OGRCoordinateTransformation
     633             : {
     634             :     const double m_dfRatioX;
     635             :     const double m_dfRatioY;
     636             : 
     637             :     GDALWarpCoordRescaler &operator=(const GDALWarpCoordRescaler &) = delete;
     638             :     GDALWarpCoordRescaler(GDALWarpCoordRescaler &&) = delete;
     639             :     GDALWarpCoordRescaler &operator=(GDALWarpCoordRescaler &&) = delete;
     640             : 
     641             :   public:
     642           1 :     GDALWarpCoordRescaler(double dfRatioX, double dfRatioY)
     643           1 :         : m_dfRatioX(dfRatioX), m_dfRatioY(dfRatioY)
     644             :     {
     645           1 :     }
     646             : 
     647           0 :     GDALWarpCoordRescaler(const GDALWarpCoordRescaler &) = default;
     648             : 
     649             :     ~GDALWarpCoordRescaler() override;
     650             : 
     651           0 :     virtual const OGRSpatialReference *GetSourceCS() const override
     652             :     {
     653           0 :         return nullptr;
     654             :     }
     655             : 
     656           6 :     virtual const OGRSpatialReference *GetTargetCS() const override
     657             :     {
     658           6 :         return nullptr;
     659             :     }
     660             : 
     661           3 :     virtual int Transform(size_t nCount, double *x, double *y, double * /*z*/,
     662             :                           double * /*t*/, int *pabSuccess) override
     663             :     {
     664          17 :         for (size_t i = 0; i < nCount; i++)
     665             :         {
     666          14 :             x[i] *= m_dfRatioX;
     667          14 :             y[i] *= m_dfRatioY;
     668          14 :             if (pabSuccess)
     669          14 :                 pabSuccess[i] = TRUE;
     670             :         }
     671           3 :         return TRUE;
     672             :     }
     673             : 
     674           0 :     virtual OGRCoordinateTransformation *Clone() const override
     675             :     {
     676           0 :         return new GDALWarpCoordRescaler(*this);
     677             :     }
     678             : 
     679           0 :     virtual OGRCoordinateTransformation *GetInverse() const override
     680             :     {
     681           0 :         return nullptr;
     682             :     }
     683             : };
     684             : 
     685             : GDALWarpCoordRescaler::~GDALWarpCoordRescaler() = default;
     686             : 
     687             : /************************************************************************/
     688             : /*                        RescaleDstGeoTransform()                      */
     689             : /************************************************************************/
     690             : 
     691          22 : static void RescaleDstGeoTransform(double adfDstGeoTransform[6],
     692             :                                    int nRasterXSize, int nDstPixels,
     693             :                                    int nRasterYSize, int nDstLines)
     694             : {
     695          22 :     adfDstGeoTransform[1] *= static_cast<double>(nRasterXSize) / nDstPixels;
     696          22 :     adfDstGeoTransform[2] *= static_cast<double>(nRasterXSize) / nDstPixels;
     697          22 :     adfDstGeoTransform[4] *= static_cast<double>(nRasterYSize) / nDstLines;
     698          22 :     adfDstGeoTransform[5] *= static_cast<double>(nRasterYSize) / nDstLines;
     699          22 : }
     700             : 
     701             : /************************************************************************/
     702             : /*                        GetSrcOverviewLevel()                         */
     703             : /************************************************************************/
     704             : 
     705          79 : int VRTWarpedDataset::GetSrcOverviewLevel(int iOvr,
     706             :                                           bool &bThisLevelOnlyOut) const
     707             : {
     708          79 :     bThisLevelOnlyOut = false;
     709          79 :     if (m_nSrcOvrLevel < -2)
     710             :     {
     711           6 :         if (iOvr + m_nSrcOvrLevel + 2 >= 0)
     712             :         {
     713           3 :             return iOvr + m_nSrcOvrLevel + 2;
     714             :         }
     715             :     }
     716          73 :     else if (m_nSrcOvrLevel == -2)
     717             :     {
     718          67 :         return iOvr;
     719             :     }
     720           6 :     else if (m_nSrcOvrLevel >= 0)
     721             :     {
     722           0 :         bThisLevelOnlyOut = true;
     723           0 :         return m_nSrcOvrLevel;
     724             :     }
     725           9 :     return -1;
     726             : }
     727             : 
     728             : /************************************************************************/
     729             : /*                            GetOverviewSize()                         */
     730             : /************************************************************************/
     731             : 
     732          73 : bool VRTWarpedDataset::GetOverviewSize(GDALDataset *poSrcDS, int iOvr,
     733             :                                        int iSrcOvr, int &nOvrXSize,
     734             :                                        int &nOvrYSize, double &dfSrcRatioX,
     735             :                                        double &dfSrcRatioY) const
     736             : {
     737             :     auto poSrcOvrBand = iSrcOvr >= 0
     738          73 :                             ? poSrcDS->GetRasterBand(1)->GetOverview(iSrcOvr)
     739           3 :                             : poSrcDS->GetRasterBand(1);
     740          73 :     if (!poSrcOvrBand)
     741             :     {
     742           0 :         return false;
     743             :     }
     744          73 :     dfSrcRatioX = static_cast<double>(poSrcDS->GetRasterXSize()) /
     745          73 :                   poSrcOvrBand->GetXSize();
     746          73 :     dfSrcRatioY = static_cast<double>(poSrcDS->GetRasterYSize()) /
     747          73 :                   poSrcOvrBand->GetYSize();
     748             :     const double dfTargetRatio =
     749          73 :         static_cast<double>(poSrcDS->GetRasterXSize()) /
     750          73 :         poSrcDS->GetRasterBand(1)->GetOverview(iOvr)->GetXSize();
     751             : 
     752          73 :     nOvrXSize = static_cast<int>(nRasterXSize / dfTargetRatio + 0.5);
     753          73 :     nOvrYSize = static_cast<int>(nRasterYSize / dfTargetRatio + 0.5);
     754          73 :     return nOvrXSize >= 1 && nOvrYSize >= 1;
     755             : }
     756             : 
     757             : /************************************************************************/
     758             : /*                        CreateImplicitOverview()                      */
     759             : /************************************************************************/
     760             : 
     761          20 : VRTWarpedDataset *VRTWarpedDataset::CreateImplicitOverview(int iOvr) const
     762             : {
     763          20 :     if (!m_poWarper)
     764           0 :         return nullptr;
     765          20 :     const GDALWarpOptions *psWO = m_poWarper->GetOptions();
     766          20 :     if (!psWO->hSrcDS || GDALGetRasterCount(psWO->hSrcDS) == 0)
     767           0 :         return nullptr;
     768          20 :     GDALDataset *poSrcDS = GDALDataset::FromHandle(psWO->hSrcDS);
     769          20 :     GDALDataset *poSrcOvrDS = poSrcDS;
     770          20 :     bool bThisLevelOnly = false;
     771          20 :     const int iSrcOvr = GetSrcOverviewLevel(iOvr, bThisLevelOnly);
     772          20 :     if (iSrcOvr >= 0)
     773             :     {
     774             :         poSrcOvrDS =
     775          17 :             GDALCreateOverviewDataset(poSrcDS, iSrcOvr, bThisLevelOnly);
     776             :     }
     777          20 :     if (poSrcOvrDS == nullptr)
     778           0 :         return nullptr;
     779          20 :     if (poSrcOvrDS == poSrcDS)
     780           3 :         poSrcOvrDS->Reference();
     781             : 
     782          20 :     int nDstPixels = 0;
     783          20 :     int nDstLines = 0;
     784          20 :     double dfSrcRatioX = 0;
     785          20 :     double dfSrcRatioY = 0;
     786             :     // Figure out the desired output bounds and resolution.
     787          20 :     if (!GetOverviewSize(poSrcDS, iOvr, iSrcOvr, nDstPixels, nDstLines,
     788             :                          dfSrcRatioX, dfSrcRatioY))
     789             :     {
     790           0 :         poSrcOvrDS->ReleaseRef();
     791           0 :         return nullptr;
     792             :     }
     793             : 
     794             :     /* --------------------------------------------------------------------
     795             :      */
     796             :     /*      Create transformer and warping options. */
     797             :     /* --------------------------------------------------------------------
     798             :      */
     799          40 :     void *pTransformerArg = GDALCreateSimilarTransformer(
     800          20 :         psWO->pTransformerArg, dfSrcRatioX, dfSrcRatioY);
     801          20 :     if (pTransformerArg == nullptr)
     802             :     {
     803           0 :         poSrcOvrDS->ReleaseRef();
     804           0 :         return nullptr;
     805             :     }
     806             : 
     807          20 :     GDALWarpOptions *psWOOvr = GDALCloneWarpOptions(psWO);
     808          20 :     psWOOvr->hSrcDS = poSrcOvrDS;
     809          20 :     psWOOvr->pfnTransformer = psWO->pfnTransformer;
     810          20 :     psWOOvr->pTransformerArg = pTransformerArg;
     811             : 
     812             :     /* --------------------------------------------------------------------
     813             :      */
     814             :     /*      We need to rescale the potential CUTLINE */
     815             :     /* --------------------------------------------------------------------
     816             :      */
     817          20 :     if (psWOOvr->hCutline)
     818             :     {
     819           2 :         GDALWarpCoordRescaler oRescaler(1.0 / dfSrcRatioX, 1.0 / dfSrcRatioY);
     820           1 :         static_cast<OGRGeometry *>(psWOOvr->hCutline)->transform(&oRescaler);
     821             :     }
     822             : 
     823             :     /* --------------------------------------------------------------------
     824             :      */
     825             :     /*      Rescale the output geotransform on the transformer. */
     826             :     /* --------------------------------------------------------------------
     827             :      */
     828          20 :     double adfDstGeoTransform[6] = {0.0};
     829          20 :     GDALGetTransformerDstGeoTransform(psWOOvr->pTransformerArg,
     830             :                                       adfDstGeoTransform);
     831          20 :     RescaleDstGeoTransform(adfDstGeoTransform, nRasterXSize, nDstPixels,
     832          20 :                            nRasterYSize, nDstLines);
     833          20 :     GDALSetTransformerDstGeoTransform(psWOOvr->pTransformerArg,
     834             :                                       adfDstGeoTransform);
     835             : 
     836             :     /* --------------------------------------------------------------------
     837             :      */
     838             :     /*      Create the VRT file. */
     839             :     /* --------------------------------------------------------------------
     840             :      */
     841          20 :     GDALDatasetH hDstDS = GDALCreateWarpedVRT(poSrcOvrDS, nDstPixels, nDstLines,
     842             :                                               adfDstGeoTransform, psWOOvr);
     843             : 
     844          20 :     poSrcOvrDS->ReleaseRef();
     845             : 
     846          20 :     GDALDestroyWarpOptions(psWOOvr);
     847             : 
     848          20 :     if (hDstDS == nullptr)
     849             :     {
     850           0 :         GDALDestroyTransformer(pTransformerArg);
     851           0 :         return nullptr;
     852             :     }
     853             : 
     854          20 :     auto poOvrDS = static_cast<VRTWarpedDataset *>(hDstDS);
     855          20 :     poOvrDS->m_bIsOverview = true;
     856          20 :     return poOvrDS;
     857             : }
     858             : 
     859             : /************************************************************************/
     860             : /*                           GetOverviewCount()                         */
     861             : /************************************************************************/
     862             : 
     863          65 : int VRTWarpedDataset::GetOverviewCount() const
     864             : {
     865          65 :     if (m_poWarper)
     866             :     {
     867          65 :         const GDALWarpOptions *psWO = m_poWarper->GetOptions();
     868          65 :         if (!m_bIsOverview && psWO->hSrcDS && GDALGetRasterCount(psWO->hSrcDS))
     869             :         {
     870          65 :             GDALDataset *poSrcDS = GDALDataset::FromHandle(psWO->hSrcDS);
     871             :             int nSrcOverviewCount =
     872          65 :                 poSrcDS->GetRasterBand(1)->GetOverviewCount();
     873          65 :             int nCount = 0;
     874         124 :             for (int i = 0; i < nSrcOverviewCount; ++i)
     875             :             {
     876          59 :                 bool bThisLevelOnly = false;
     877          59 :                 const int iSrcOvr = GetSrcOverviewLevel(i, bThisLevelOnly);
     878          59 :                 if (iSrcOvr >= 0)
     879             :                 {
     880          53 :                     int nDstPixels = 0;
     881          53 :                     int nDstLines = 0;
     882          53 :                     double dfSrcRatioX = 0;
     883          53 :                     double dfSrcRatioY = 0;
     884          53 :                     if (!GetOverviewSize(poSrcDS, i, iSrcOvr, nDstPixels,
     885             :                                          nDstLines, dfSrcRatioX, dfSrcRatioY))
     886             :                     {
     887           0 :                         break;
     888             :                     }
     889             :                 }
     890          59 :                 ++nCount;
     891             :             }
     892          65 :             return nCount;
     893             :         }
     894             :     }
     895           0 :     return 0;
     896             : }
     897             : 
     898             : /************************************************************************/
     899             : /*                        CreateImplicitOverviews()                     */
     900             : /*                                                                      */
     901             : /*      For each overview of the source dataset, create an overview     */
     902             : /*      in the warped VRT dataset.                                      */
     903             : /************************************************************************/
     904             : 
     905           8 : void VRTWarpedDataset::CreateImplicitOverviews()
     906             : {
     907           8 :     if (m_bIsOverview)
     908           0 :         return;
     909           8 :     const int nOvrCount = GetOverviewCount();
     910           8 :     if (m_apoOverviews.empty())
     911           4 :         m_apoOverviews.resize(nOvrCount);
     912          18 :     for (int iOvr = 0; iOvr < nOvrCount; iOvr++)
     913             :     {
     914          10 :         if (!m_apoOverviews[iOvr])
     915             :         {
     916           2 :             m_apoOverviews[iOvr] = CreateImplicitOverview(iOvr);
     917             :         }
     918             :     }
     919             : }
     920             : 
     921             : /************************************************************************/
     922             : /*                            GetFileList()                             */
     923             : /************************************************************************/
     924             : 
     925          20 : char **VRTWarpedDataset::GetFileList()
     926             : {
     927          20 :     char **papszFileList = GDALDataset::GetFileList();
     928             : 
     929          20 :     if (m_poWarper != nullptr)
     930             :     {
     931          20 :         const GDALWarpOptions *psWO = m_poWarper->GetOptions();
     932          20 :         const char *pszFilename = nullptr;
     933             : 
     934          40 :         if (psWO->hSrcDS != nullptr &&
     935             :             (pszFilename =
     936          20 :                  static_cast<GDALDataset *>(psWO->hSrcDS)->GetDescription()) !=
     937             :                 nullptr)
     938             :         {
     939             :             VSIStatBufL sStat;
     940          20 :             if (VSIStatL(pszFilename, &sStat) == 0)
     941             :             {
     942          20 :                 papszFileList = CSLAddString(papszFileList, pszFilename);
     943             :             }
     944             :         }
     945             :     }
     946             : 
     947          20 :     return papszFileList;
     948             : }
     949             : 
     950             : /************************************************************************/
     951             : /* ==================================================================== */
     952             : /*                    VRTWarpedOverviewTransformer                      */
     953             : /* ==================================================================== */
     954             : /************************************************************************/
     955             : 
     956             : typedef struct
     957             : {
     958             :     GDALTransformerInfo sTI;
     959             : 
     960             :     GDALTransformerFunc pfnBaseTransformer;
     961             :     void *pBaseTransformerArg;
     962             :     bool bOwnSubtransformer;
     963             : 
     964             :     double dfXOverviewFactor;
     965             :     double dfYOverviewFactor;
     966             : } VWOTInfo;
     967             : 
     968             : static void *VRTCreateWarpedOverviewTransformer(
     969             :     GDALTransformerFunc pfnBaseTransformer, void *pBaseTransformArg,
     970             :     double dfXOverviewFactor, double dfYOverviewFactor);
     971             : static void VRTDestroyWarpedOverviewTransformer(void *pTransformArg);
     972             : 
     973             : static int VRTWarpedOverviewTransform(void *pTransformArg, int bDstToSrc,
     974             :                                       int nPointCount, double *padfX,
     975             :                                       double *padfY, double *padfZ,
     976             :                                       int *panSuccess);
     977             : 
     978             : #if 0   // TODO: Why?
     979             : /************************************************************************/
     980             : /*                VRTSerializeWarpedOverviewTransformer()               */
     981             : /************************************************************************/
     982             : 
     983             : static CPLXMLNode *
     984             : VRTSerializeWarpedOverviewTransformer( void *pTransformArg )
     985             : 
     986             : {
     987             :     VWOTInfo *psInfo = static_cast<VWOTInfo *>( pTransformArg );
     988             : 
     989             :     CPLXMLNode *psTree
     990             :         = CPLCreateXMLNode( NULL, CXT_Element, "WarpedOverviewTransformer" );
     991             : 
     992             :     CPLCreateXMLElementAndValue(
     993             :         psTree, "XFactor",
     994             :         CPLString().Printf("%g",psInfo->dfXOverviewFactor) );
     995             :     CPLCreateXMLElementAndValue(
     996             :         psTree, "YFactor",
     997             :         CPLString().Printf("%g",psInfo->dfYOverviewFactor) );
     998             : 
     999             : /* -------------------------------------------------------------------- */
    1000             : /*      Capture underlying transformer.                                 */
    1001             : /* -------------------------------------------------------------------- */
    1002             :     CPLXMLNode *psTransformerContainer
    1003             :         = CPLCreateXMLNode( psTree, CXT_Element, "BaseTransformer" );
    1004             : 
    1005             :     CPLXMLNode *psTransformer
    1006             :         = GDALSerializeTransformer( psInfo->pfnBaseTransformer,
    1007             :                                     psInfo->pBaseTransformerArg );
    1008             :     if( psTransformer != NULL )
    1009             :         CPLAddXMLChild( psTransformerContainer, psTransformer );
    1010             : 
    1011             :     return psTree;
    1012             : }
    1013             : 
    1014             : /************************************************************************/
    1015             : /*           VRTWarpedOverviewTransformerOwnsSubtransformer()           */
    1016             : /************************************************************************/
    1017             : 
    1018             : static void VRTWarpedOverviewTransformerOwnsSubtransformer( void *pTransformArg,
    1019             :                                                             bool bOwnFlag )
    1020             : {
    1021             :     VWOTInfo *psInfo = static_cast<VWOTInfo *>( pTransformArg );
    1022             : 
    1023             :     psInfo->bOwnSubtransformer = bOwnFlag;
    1024             : }
    1025             : 
    1026             : /************************************************************************/
    1027             : /*            VRTDeserializeWarpedOverviewTransformer()                 */
    1028             : /************************************************************************/
    1029             : 
    1030             : void* VRTDeserializeWarpedOverviewTransformer( CPLXMLNode *psTree )
    1031             : 
    1032             : {
    1033             :     const double dfXOverviewFactor =
    1034             :         CPLAtof(CPLGetXMLValue( psTree, "XFactor",  "1" ));
    1035             :     const double dfYOverviewFactor =
    1036             :         CPLAtof(CPLGetXMLValue( psTree, "YFactor",  "1" ));
    1037             :     GDALTransformerFunc pfnBaseTransform = NULL;
    1038             :     void *pBaseTransformerArg = NULL;
    1039             : 
    1040             :     CPLXMLNode *psContainer = CPLGetXMLNode( psTree, "BaseTransformer" );
    1041             : 
    1042             :     if( psContainer != NULL && psContainer->psChild != NULL )
    1043             :     {
    1044             :         GDALDeserializeTransformer( psContainer->psChild,
    1045             :                                     &pfnBaseTransform,
    1046             :                                     &pBaseTransformerArg );
    1047             :     }
    1048             : 
    1049             :     if( pfnBaseTransform == NULL )
    1050             :     {
    1051             :         CPLError( CE_Failure, CPLE_AppDefined,
    1052             :                   "Cannot get base transform for scaled coord transformer." );
    1053             :         return NULL;
    1054             :     }
    1055             :     else
    1056             :     {
    1057             :         void *pApproxCBData =
    1058             :                        VRTCreateWarpedOverviewTransformer( pfnBaseTransform,
    1059             :                                                            pBaseTransformerArg,
    1060             :                                                            dfXOverviewFactor,
    1061             :                                                            dfYOverviewFactor );
    1062             :         VRTWarpedOverviewTransformerOwnsSubtransformer( pApproxCBData, true );
    1063             : 
    1064             :         return pApproxCBData;
    1065             :     }
    1066             : }
    1067             : #endif  // TODO: Why disabled?
    1068             : 
    1069             : /************************************************************************/
    1070             : /*                   VRTCreateWarpedOverviewTransformer()               */
    1071             : /************************************************************************/
    1072             : 
    1073           4 : static void *VRTCreateWarpedOverviewTransformer(
    1074             :     GDALTransformerFunc pfnBaseTransformer, void *pBaseTransformerArg,
    1075             :     double dfXOverviewFactor, double dfYOverviewFactor)
    1076             : 
    1077             : {
    1078           4 :     if (pfnBaseTransformer == nullptr)
    1079           0 :         return nullptr;
    1080             : 
    1081           4 :     VWOTInfo *psSCTInfo = static_cast<VWOTInfo *>(CPLMalloc(sizeof(VWOTInfo)));
    1082           4 :     psSCTInfo->pfnBaseTransformer = pfnBaseTransformer;
    1083           4 :     psSCTInfo->pBaseTransformerArg = pBaseTransformerArg;
    1084           4 :     psSCTInfo->dfXOverviewFactor = dfXOverviewFactor;
    1085           4 :     psSCTInfo->dfYOverviewFactor = dfYOverviewFactor;
    1086           4 :     psSCTInfo->bOwnSubtransformer = false;
    1087             : 
    1088           4 :     memcpy(psSCTInfo->sTI.abySignature, GDAL_GTI2_SIGNATURE,
    1089             :            strlen(GDAL_GTI2_SIGNATURE));
    1090           4 :     psSCTInfo->sTI.pszClassName = "VRTWarpedOverviewTransformer";
    1091           4 :     psSCTInfo->sTI.pfnTransform = VRTWarpedOverviewTransform;
    1092           4 :     psSCTInfo->sTI.pfnCleanup = VRTDestroyWarpedOverviewTransformer;
    1093             : #if 0
    1094             :     psSCTInfo->sTI.pfnSerialize = VRTSerializeWarpedOverviewTransformer;
    1095             : #endif
    1096           4 :     return psSCTInfo;
    1097             : }
    1098             : 
    1099             : /************************************************************************/
    1100             : /*               VRTDestroyWarpedOverviewTransformer()                  */
    1101             : /************************************************************************/
    1102             : 
    1103           4 : static void VRTDestroyWarpedOverviewTransformer(void *pTransformArg)
    1104             : {
    1105           4 :     VWOTInfo *psInfo = static_cast<VWOTInfo *>(pTransformArg);
    1106             : 
    1107           4 :     if (psInfo->bOwnSubtransformer)
    1108           0 :         GDALDestroyTransformer(psInfo->pBaseTransformerArg);
    1109             : 
    1110           4 :     CPLFree(psInfo);
    1111           4 : }
    1112             : 
    1113             : /************************************************************************/
    1114             : /*                     VRTWarpedOverviewTransform()                     */
    1115             : /************************************************************************/
    1116             : 
    1117         259 : static int VRTWarpedOverviewTransform(void *pTransformArg, int bDstToSrc,
    1118             :                                       int nPointCount, double *padfX,
    1119             :                                       double *padfY, double *padfZ,
    1120             :                                       int *panSuccess)
    1121             : 
    1122             : {
    1123         259 :     VWOTInfo *psInfo = static_cast<VWOTInfo *>(pTransformArg);
    1124             : 
    1125         259 :     if (bDstToSrc)
    1126             :     {
    1127       63032 :         for (int i = 0; i < nPointCount; i++)
    1128             :         {
    1129       62773 :             padfX[i] *= psInfo->dfXOverviewFactor;
    1130       62773 :             padfY[i] *= psInfo->dfYOverviewFactor;
    1131             :         }
    1132             :     }
    1133             : 
    1134         259 :     const int bSuccess = psInfo->pfnBaseTransformer(
    1135             :         psInfo->pBaseTransformerArg, bDstToSrc, nPointCount, padfX, padfY,
    1136             :         padfZ, panSuccess);
    1137             : 
    1138         259 :     if (!bDstToSrc)
    1139             :     {
    1140           0 :         for (int i = 0; i < nPointCount; i++)
    1141             :         {
    1142           0 :             padfX[i] /= psInfo->dfXOverviewFactor;
    1143           0 :             padfY[i] /= psInfo->dfYOverviewFactor;
    1144             :         }
    1145             :     }
    1146             : 
    1147         259 :     return bSuccess;
    1148             : }
    1149             : 
    1150             : /************************************************************************/
    1151             : /*                           BuildOverviews()                           */
    1152             : /*                                                                      */
    1153             : /*      For overviews, we actually just build a whole new dataset       */
    1154             : /*      with an extra layer of transformation on the warper used to     */
    1155             : /*      accomplish downsampling by the desired factor.                  */
    1156             : /************************************************************************/
    1157             : 
    1158           6 : CPLErr VRTWarpedDataset::IBuildOverviews(
    1159             :     const char * /* pszResampling */, int nOverviews,
    1160             :     const int *panOverviewList, int /* nListBands */,
    1161             :     const int * /* panBandList */, GDALProgressFunc pfnProgress,
    1162             :     void *pProgressData, CSLConstList /*papszOptions*/)
    1163             : {
    1164           6 :     if (m_poWarper == nullptr || m_bIsOverview)
    1165           0 :         return CE_Failure;
    1166             : 
    1167             :     /* -------------------------------------------------------------------- */
    1168             :     /*      Initial progress result.                                        */
    1169             :     /* -------------------------------------------------------------------- */
    1170           6 :     if (!pfnProgress(0.0, nullptr, pProgressData))
    1171             :     {
    1172           0 :         CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
    1173           0 :         return CE_Failure;
    1174             :     }
    1175             : 
    1176           6 :     CreateImplicitOverviews();
    1177             : 
    1178             :     /* -------------------------------------------------------------------- */
    1179             :     /*      Establish which of the overview levels we already have, and     */
    1180             :     /*      which are new.                                                  */
    1181             :     /* -------------------------------------------------------------------- */
    1182           6 :     int nNewOverviews = 0;
    1183             :     int *panNewOverviewList =
    1184           6 :         static_cast<int *>(CPLCalloc(sizeof(int), nOverviews));
    1185           6 :     std::vector<bool> abFoundOverviewFactor(nOverviews);
    1186          14 :     for (int i = 0; i < nOverviews; i++)
    1187             :     {
    1188          20 :         for (GDALDataset *const poOverview : m_apoOverviews)
    1189             :         {
    1190          12 :             if (poOverview)
    1191             :             {
    1192          12 :                 const int nOvFactor = GDALComputeOvFactor(
    1193             :                     poOverview->GetRasterXSize(), GetRasterXSize(),
    1194             :                     poOverview->GetRasterYSize(), GetRasterYSize());
    1195             : 
    1196          20 :                 if (nOvFactor == panOverviewList[i] ||
    1197           8 :                     nOvFactor == GDALOvLevelAdjust2(panOverviewList[i],
    1198             :                                                     GetRasterXSize(),
    1199             :                                                     GetRasterYSize()))
    1200           4 :                     abFoundOverviewFactor[i] = true;
    1201             :             }
    1202             :         }
    1203             : 
    1204           8 :         if (!abFoundOverviewFactor[i])
    1205           4 :             panNewOverviewList[nNewOverviews++] = panOverviewList[i];
    1206             :     }
    1207             : 
    1208             :     /* -------------------------------------------------------------------- */
    1209             :     /*      Create each missing overview (we don't need to do anything      */
    1210             :     /*      to update existing overviews).                                  */
    1211             :     /* -------------------------------------------------------------------- */
    1212           6 :     CPLErr eErr = CE_None;
    1213          10 :     for (int i = 0; i < nNewOverviews; i++)
    1214             :     {
    1215             :         /* --------------------------------------------------------------------
    1216             :          */
    1217             :         /*      What size should this overview be. */
    1218             :         /* --------------------------------------------------------------------
    1219             :          */
    1220           4 :         const int nOXSize = (GetRasterXSize() + panNewOverviewList[i] - 1) /
    1221           4 :                             panNewOverviewList[i];
    1222             : 
    1223           4 :         const int nOYSize = (GetRasterYSize() + panNewOverviewList[i] - 1) /
    1224           4 :                             panNewOverviewList[i];
    1225             : 
    1226             :         /* --------------------------------------------------------------------
    1227             :          */
    1228             :         /*      Find the most appropriate base dataset onto which to build the
    1229             :          */
    1230             :         /*      new one. The preference will be an overview dataset with a
    1231             :          * ratio*/
    1232             :         /*      greater than ours, and which is not using */
    1233             :         /*      VRTWarpedOverviewTransform, since those ones are slow. The
    1234             :          * other*/
    1235             :         /*      ones are based on overviews of the source dataset. */
    1236             :         /* --------------------------------------------------------------------
    1237             :          */
    1238           4 :         VRTWarpedDataset *poBaseDataset = this;
    1239           8 :         for (auto *poOverview : m_apoOverviews)
    1240             :         {
    1241           4 :             if (poOverview && poOverview->GetRasterXSize() > nOXSize &&
    1242           4 :                 poOverview->m_poWarper->GetOptions()->pfnTransformer !=
    1243           8 :                     VRTWarpedOverviewTransform &&
    1244           4 :                 poOverview->GetRasterXSize() < poBaseDataset->GetRasterXSize())
    1245             :             {
    1246           4 :                 poBaseDataset = poOverview;
    1247             :             }
    1248             :         }
    1249             : 
    1250             :         /* --------------------------------------------------------------------
    1251             :          */
    1252             :         /*      Create the overview dataset. */
    1253             :         /* --------------------------------------------------------------------
    1254             :          */
    1255           4 :         VRTWarpedDataset *poOverviewDS = new VRTWarpedDataset(nOXSize, nOYSize);
    1256             : 
    1257           8 :         for (int iBand = 0; iBand < GetRasterCount(); iBand++)
    1258             :         {
    1259           4 :             GDALRasterBand *const poOldBand = GetRasterBand(iBand + 1);
    1260             :             VRTWarpedRasterBand *const poNewBand = new VRTWarpedRasterBand(
    1261           4 :                 poOverviewDS, iBand + 1, poOldBand->GetRasterDataType());
    1262             : 
    1263           4 :             poNewBand->CopyCommonInfoFrom(poOldBand);
    1264           4 :             poOverviewDS->SetBand(iBand + 1, poNewBand);
    1265             :         }
    1266             : 
    1267             :         /* --------------------------------------------------------------------
    1268             :          */
    1269             :         /*      Prepare update transformation information that will apply */
    1270             :         /*      the overview decimation. */
    1271             :         /* --------------------------------------------------------------------
    1272             :          */
    1273             :         GDALWarpOptions *psWO = const_cast<GDALWarpOptions *>(
    1274           4 :             poBaseDataset->m_poWarper->GetOptions());
    1275             : 
    1276             :         /* --------------------------------------------------------------------
    1277             :          */
    1278             :         /*      Initialize the new dataset with adjusted warp options, and */
    1279             :         /*      then restore to original condition. */
    1280             :         /* --------------------------------------------------------------------
    1281             :          */
    1282           4 :         GDALTransformerFunc pfnTransformerBase = psWO->pfnTransformer;
    1283           4 :         void *pTransformerBaseArg = psWO->pTransformerArg;
    1284             : 
    1285           4 :         psWO->pfnTransformer = VRTWarpedOverviewTransform;
    1286          12 :         psWO->pTransformerArg = VRTCreateWarpedOverviewTransformer(
    1287             :             pfnTransformerBase, pTransformerBaseArg,
    1288           4 :             poBaseDataset->GetRasterXSize() / static_cast<double>(nOXSize),
    1289           4 :             poBaseDataset->GetRasterYSize() / static_cast<double>(nOYSize));
    1290             : 
    1291           4 :         eErr = poOverviewDS->Initialize(psWO);
    1292             : 
    1293           4 :         psWO->pfnTransformer = pfnTransformerBase;
    1294           4 :         psWO->pTransformerArg = pTransformerBaseArg;
    1295             : 
    1296           4 :         if (eErr != CE_None)
    1297             :         {
    1298           0 :             delete poOverviewDS;
    1299           0 :             break;
    1300             :         }
    1301             : 
    1302           4 :         m_apoOverviews.push_back(poOverviewDS);
    1303             :     }
    1304             : 
    1305           6 :     CPLFree(panNewOverviewList);
    1306             : 
    1307             :     /* -------------------------------------------------------------------- */
    1308             :     /*      Progress finished.                                              */
    1309             :     /* -------------------------------------------------------------------- */
    1310           6 :     pfnProgress(1.0, nullptr, pProgressData);
    1311             : 
    1312           6 :     SetNeedsFlush();
    1313             : 
    1314           6 :     return eErr;
    1315             : }
    1316             : 
    1317             : /*! @endcond */
    1318             : 
    1319             : /************************************************************************/
    1320             : /*                      GDALInitializeWarpedVRT()                       */
    1321             : /************************************************************************/
    1322             : 
    1323             : /**
    1324             :  * Set warp info on virtual warped dataset.
    1325             :  *
    1326             :  * Initializes all the warping information for a virtual warped dataset.
    1327             :  *
    1328             :  * This method is the same as the C++ method VRTWarpedDataset::Initialize().
    1329             :  *
    1330             :  * @param hDS dataset previously created with the VRT driver, and a
    1331             :  * SUBCLASS of "VRTWarpedDataset".
    1332             :  *
    1333             :  * @param psWO the warp options to apply.  Note that ownership of the
    1334             :  * transformation information is taken over by the function though everything
    1335             :  * else remains the property of the caller.
    1336             :  *
    1337             :  * @return CE_None on success or CE_Failure if an error occurs.
    1338             :  */
    1339             : 
    1340          94 : CPLErr CPL_STDCALL GDALInitializeWarpedVRT(GDALDatasetH hDS,
    1341             :                                            GDALWarpOptions *psWO)
    1342             : 
    1343             : {
    1344          94 :     VALIDATE_POINTER1(hDS, "GDALInitializeWarpedVRT", CE_Failure);
    1345             : 
    1346          94 :     return static_cast<VRTWarpedDataset *>(GDALDataset::FromHandle(hDS))
    1347          94 :         ->Initialize(psWO);
    1348             : }
    1349             : 
    1350             : /*! @cond Doxygen_Suppress */
    1351             : 
    1352             : /************************************************************************/
    1353             : /*                              XMLInit()                               */
    1354             : /************************************************************************/
    1355             : 
    1356         201 : CPLErr VRTWarpedDataset::XMLInit(const CPLXMLNode *psTree,
    1357             :                                  const char *pszVRTPathIn)
    1358             : 
    1359             : {
    1360             : 
    1361             :     /* -------------------------------------------------------------------- */
    1362             :     /*      Initialize blocksize before calling sub-init so that the        */
    1363             :     /*      band initializers can get it from the dataset object when       */
    1364             :     /*      they are created.                                               */
    1365             :     /* -------------------------------------------------------------------- */
    1366         201 :     m_nBlockXSize = atoi(CPLGetXMLValue(psTree, "BlockXSize", "512"));
    1367         201 :     m_nBlockYSize = atoi(CPLGetXMLValue(psTree, "BlockYSize", "128"));
    1368             : 
    1369             :     /* -------------------------------------------------------------------- */
    1370             :     /*      Initialize all the general VRT stuff.  This will even           */
    1371             :     /*      create the VRTWarpedRasterBands and initialize them.            */
    1372             :     /* -------------------------------------------------------------------- */
    1373             :     {
    1374         201 :         const CPLErr eErr = VRTDataset::XMLInit(psTree, pszVRTPathIn);
    1375             : 
    1376         201 :         if (eErr != CE_None)
    1377           0 :             return eErr;
    1378             :     }
    1379             : 
    1380             :     // Check that band block sizes didn't change the dataset block size.
    1381         650 :     for (int i = 1; i <= nBands; i++)
    1382             :     {
    1383         451 :         int nBlockXSize = 0;
    1384         451 :         int nBlockYSize = 0;
    1385         451 :         GetRasterBand(i)->GetBlockSize(&nBlockXSize, &nBlockYSize);
    1386         451 :         if (nBlockXSize != m_nBlockXSize || nBlockYSize != m_nBlockYSize)
    1387             :         {
    1388           2 :             CPLError(CE_Failure, CPLE_AppDefined,
    1389             :                      "Block size specified on band %d not consistent with "
    1390             :                      "dataset block size",
    1391             :                      i);
    1392           2 :             return CE_Failure;
    1393             :         }
    1394             :     }
    1395             : 
    1396         199 :     if (nBands > 1)
    1397             :     {
    1398          93 :         GDALDataset::SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
    1399             :     }
    1400             : 
    1401             :     /* -------------------------------------------------------------------- */
    1402             :     /*      Find the GDALWarpOptions XML tree.                              */
    1403             :     /* -------------------------------------------------------------------- */
    1404             :     const CPLXMLNode *const psOptionsTree =
    1405         199 :         CPLGetXMLNode(psTree, "GDALWarpOptions");
    1406         199 :     if (psOptionsTree == nullptr)
    1407             :     {
    1408           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    1409             :                  "Count not find required GDALWarpOptions in XML.");
    1410           0 :         return CE_Failure;
    1411             :     }
    1412             : 
    1413             :     /* -------------------------------------------------------------------- */
    1414             :     /*      Adjust the SourceDataset in the warp options to take into       */
    1415             :     /*      account that it is relative to the VRT if appropriate.          */
    1416             :     /* -------------------------------------------------------------------- */
    1417         199 :     const bool bRelativeToVRT = CPL_TO_BOOL(atoi(
    1418             :         CPLGetXMLValue(psOptionsTree, "SourceDataset.relativeToVRT", "0")));
    1419             : 
    1420             :     const char *pszRelativePath =
    1421         199 :         CPLGetXMLValue(psOptionsTree, "SourceDataset", "");
    1422         199 :     char *pszAbsolutePath = nullptr;
    1423             : 
    1424         199 :     if (bRelativeToVRT)
    1425         132 :         pszAbsolutePath = CPLStrdup(
    1426         264 :             CPLProjectRelativeFilenameSafe(pszVRTPathIn, pszRelativePath)
    1427             :                 .c_str());
    1428             :     else
    1429          67 :         pszAbsolutePath = CPLStrdup(pszRelativePath);
    1430             : 
    1431         199 :     CPLXMLNode *psOptionsTreeCloned = CPLCloneXMLTree(psOptionsTree);
    1432         199 :     CPLSetXMLValue(psOptionsTreeCloned, "SourceDataset", pszAbsolutePath);
    1433         199 :     CPLFree(pszAbsolutePath);
    1434             : 
    1435             :     /* -------------------------------------------------------------------- */
    1436             :     /*      And instantiate the warp options, and corresponding warp        */
    1437             :     /*      operation.                                                      */
    1438             :     /* -------------------------------------------------------------------- */
    1439         199 :     GDALWarpOptions *psWO = GDALDeserializeWarpOptions(psOptionsTreeCloned);
    1440         199 :     CPLDestroyXMLNode(psOptionsTreeCloned);
    1441         199 :     if (psWO == nullptr)
    1442           0 :         return CE_Failure;
    1443             : 
    1444         199 :     psWO->papszWarpOptions = VRTWarpedAddOptions(psWO->papszWarpOptions);
    1445             : 
    1446         199 :     eAccess = GA_Update;
    1447             : 
    1448         199 :     if (psWO->hDstDS != nullptr)
    1449             :     {
    1450           0 :         GDALClose(psWO->hDstDS);
    1451           0 :         psWO->hDstDS = nullptr;
    1452             :     }
    1453             : 
    1454         199 :     psWO->hDstDS = this;
    1455             : 
    1456             :     /* -------------------------------------------------------------------- */
    1457             :     /*      Deserialize vertical shift grids.                               */
    1458             :     /* -------------------------------------------------------------------- */
    1459         199 :     CPLXMLNode *psIter = psTree->psChild;
    1460        2406 :     for (; psWO->hSrcDS != nullptr && psIter != nullptr;
    1461        2207 :          psIter = psIter->psNext)
    1462             :     {
    1463        2207 :         if (psIter->eType != CXT_Element ||
    1464        1610 :             !EQUAL(psIter->pszValue, "VerticalShiftGrids"))
    1465             :         {
    1466        2207 :             continue;
    1467             :         }
    1468             : 
    1469           0 :         CPLError(CE_Warning, CPLE_AppDefined,
    1470             :                  "The VerticalShiftGrids in a warped VRT is now deprecated, "
    1471             :                  "and will no longer be handled in GDAL 4.0");
    1472             : 
    1473           0 :         const char *pszVGrids = CPLGetXMLValue(psIter, "Grids", nullptr);
    1474           0 :         if (pszVGrids)
    1475             :         {
    1476             :             int bInverse =
    1477           0 :                 CSLTestBoolean(CPLGetXMLValue(psIter, "Inverse", "FALSE"));
    1478             :             double dfToMeterSrc =
    1479           0 :                 CPLAtof(CPLGetXMLValue(psIter, "ToMeterSrc", "1.0"));
    1480             :             double dfToMeterDest =
    1481           0 :                 CPLAtof(CPLGetXMLValue(psIter, "ToMeterDest", "1.0"));
    1482           0 :             char **papszOptions = nullptr;
    1483           0 :             CPLXMLNode *psIter2 = psIter->psChild;
    1484           0 :             for (; psIter2 != nullptr; psIter2 = psIter2->psNext)
    1485             :             {
    1486           0 :                 if (psIter2->eType != CXT_Element ||
    1487           0 :                     !EQUAL(psIter2->pszValue, "Option"))
    1488             :                 {
    1489           0 :                     continue;
    1490             :                 }
    1491           0 :                 const char *pszName = CPLGetXMLValue(psIter2, "name", nullptr);
    1492             :                 const char *pszValue =
    1493           0 :                     CPLGetXMLValue(psIter2, nullptr, nullptr);
    1494           0 :                 if (pszName && pszValue)
    1495             :                 {
    1496             :                     papszOptions =
    1497           0 :                         CSLSetNameValue(papszOptions, pszName, pszValue);
    1498             :                 }
    1499             :             }
    1500             : 
    1501           0 :             int bError = FALSE;
    1502             :             GDALDatasetH hGridDataset =
    1503           0 :                 GDALOpenVerticalShiftGrid(pszVGrids, &bError);
    1504           0 :             if (bError && hGridDataset == nullptr)
    1505             :             {
    1506           0 :                 CPLError(CE_Warning, CPLE_AppDefined,
    1507             :                          "Cannot open %s. Source dataset will no "
    1508             :                          "be vertically adjusted regarding "
    1509             :                          "vertical datum",
    1510             :                          pszVGrids);
    1511             :             }
    1512           0 :             else if (hGridDataset != nullptr)
    1513             :             {
    1514             :                 // Transform from source vertical datum to WGS84
    1515           0 :                 GDALDatasetH hTmpDS = GDALApplyVerticalShiftGrid(
    1516             :                     psWO->hSrcDS, hGridDataset, bInverse, dfToMeterSrc,
    1517             :                     dfToMeterDest, papszOptions);
    1518           0 :                 GDALReleaseDataset(hGridDataset);
    1519           0 :                 if (hTmpDS == nullptr)
    1520             :                 {
    1521           0 :                     CPLError(CE_Warning, CPLE_AppDefined,
    1522             :                              "Source dataset will no "
    1523             :                              "be vertically adjusted regarding "
    1524             :                              "vertical datum %s",
    1525             :                              pszVGrids);
    1526             :                 }
    1527             :                 else
    1528             :                 {
    1529           0 :                     CPLDebug("GDALWARP",
    1530             :                              "Adjusting source dataset "
    1531             :                              "with vertical datum using %s",
    1532             :                              pszVGrids);
    1533           0 :                     GDALReleaseDataset(psWO->hSrcDS);
    1534           0 :                     psWO->hSrcDS = hTmpDS;
    1535             :                 }
    1536             :             }
    1537             : 
    1538           0 :             CSLDestroy(papszOptions);
    1539             :         }
    1540             :     }
    1541             : 
    1542             :     /* -------------------------------------------------------------------- */
    1543             :     /*      Instantiate the warp operation.                                 */
    1544             :     /* -------------------------------------------------------------------- */
    1545         199 :     m_poWarper = new GDALWarpOperation();
    1546             : 
    1547         199 :     const CPLErr eErr = m_poWarper->Initialize(psWO);
    1548         199 :     if (eErr != CE_None)
    1549             :     {
    1550             :         /* --------------------------------------------------------------------
    1551             :          */
    1552             :         /*      We are responsible for cleaning up the transformer ourselves. */
    1553             :         /* --------------------------------------------------------------------
    1554             :          */
    1555           0 :         if (psWO->pTransformerArg != nullptr)
    1556             :         {
    1557           0 :             GDALDestroyTransformer(psWO->pTransformerArg);
    1558           0 :             psWO->pTransformerArg = nullptr;
    1559             :         }
    1560             : 
    1561           0 :         if (psWO->hSrcDS != nullptr)
    1562             :         {
    1563           0 :             GDALClose(psWO->hSrcDS);
    1564           0 :             psWO->hSrcDS = nullptr;
    1565             :         }
    1566             :     }
    1567             : 
    1568         199 :     GDALDestroyWarpOptions(psWO);
    1569         199 :     if (eErr != CE_None)
    1570             :     {
    1571           0 :         delete m_poWarper;
    1572           0 :         m_poWarper = nullptr;
    1573             :     }
    1574             : 
    1575             :     /* -------------------------------------------------------------------- */
    1576             :     /*      Deserialize SrcOvrLevel                                         */
    1577             :     /* -------------------------------------------------------------------- */
    1578         199 :     const char *pszSrcOvrLevel = CPLGetXMLValue(psTree, "SrcOvrLevel", nullptr);
    1579         199 :     if (pszSrcOvrLevel != nullptr)
    1580             :     {
    1581           0 :         SetMetadataItem("SrcOvrLevel", pszSrcOvrLevel);
    1582             :     }
    1583             : 
    1584             :     /* -------------------------------------------------------------------- */
    1585             :     /*      Generate overviews, if appropriate.                             */
    1586             :     /* -------------------------------------------------------------------- */
    1587             : 
    1588             :     // OverviewList is historical, and quite inefficient, since it uses
    1589             :     // the full resolution source dataset, so only build it afterwards.
    1590             :     const CPLStringList aosOverviews(
    1591         199 :         CSLTokenizeString(CPLGetXMLValue(psTree, "OverviewList", "")));
    1592         199 :     if (!aosOverviews.empty())
    1593           2 :         CreateImplicitOverviews();
    1594             : 
    1595         203 :     for (int iOverview = 0; iOverview < aosOverviews.size(); ++iOverview)
    1596             :     {
    1597           4 :         int nOvFactor = atoi(aosOverviews[iOverview]);
    1598             : 
    1599           4 :         if (nOvFactor > 0)
    1600           4 :             BuildOverviews("NEAREST", 1, &nOvFactor, 0, nullptr, nullptr,
    1601             :                            nullptr,
    1602             :                            /*papszOptions=*/nullptr);
    1603             :         else
    1604           0 :             CPLError(CE_Failure, CPLE_AppDefined,
    1605             :                      "Bad value for overview factor : %s",
    1606             :                      aosOverviews[iOverview]);
    1607             :     }
    1608             : 
    1609         199 :     return eErr;
    1610             : }
    1611             : 
    1612             : /************************************************************************/
    1613             : /*                           SerializeToXML()                           */
    1614             : /************************************************************************/
    1615             : 
    1616          78 : CPLXMLNode *VRTWarpedDataset::SerializeToXML(const char *pszVRTPathIn)
    1617             : 
    1618             : {
    1619          78 :     CPLXMLNode *psTree = VRTDataset::SerializeToXML(pszVRTPathIn);
    1620             : 
    1621          78 :     if (psTree == nullptr)
    1622           0 :         return psTree;
    1623             : 
    1624             :     /* -------------------------------------------------------------------- */
    1625             :     /*      Set subclass.                                                   */
    1626             :     /* -------------------------------------------------------------------- */
    1627          78 :     CPLCreateXMLNode(CPLCreateXMLNode(psTree, CXT_Attribute, "subClass"),
    1628             :                      CXT_Text, "VRTWarpedDataset");
    1629             : 
    1630             :     /* -------------------------------------------------------------------- */
    1631             :     /*      Serialize the block size.                                       */
    1632             :     /* -------------------------------------------------------------------- */
    1633          78 :     CPLCreateXMLElementAndValue(psTree, "BlockXSize",
    1634             :                                 CPLSPrintf("%d", m_nBlockXSize));
    1635          78 :     CPLCreateXMLElementAndValue(psTree, "BlockYSize",
    1636             :                                 CPLSPrintf("%d", m_nBlockYSize));
    1637             : 
    1638             :     /* -------------------------------------------------------------------- */
    1639             :     /*      Serialize the overview list (only for non implicit overviews)   */
    1640             :     /* -------------------------------------------------------------------- */
    1641          78 :     if (!m_apoOverviews.empty())
    1642             :     {
    1643           3 :         int nSrcDSOvrCount = 0;
    1644           3 :         if (m_poWarper != nullptr && m_poWarper->GetOptions() != nullptr &&
    1645           9 :             m_poWarper->GetOptions()->hSrcDS != nullptr &&
    1646           3 :             GDALGetRasterCount(m_poWarper->GetOptions()->hSrcDS) > 0)
    1647             :         {
    1648             :             nSrcDSOvrCount =
    1649           3 :                 static_cast<GDALDataset *>(m_poWarper->GetOptions()->hSrcDS)
    1650           3 :                     ->GetRasterBand(1)
    1651           3 :                     ->GetOverviewCount();
    1652             :         }
    1653             : 
    1654           3 :         if (static_cast<int>(m_apoOverviews.size()) != nSrcDSOvrCount)
    1655             :         {
    1656           2 :             const size_t nLen = m_apoOverviews.size() * 8 + 10;
    1657           2 :             char *pszOverviewList = static_cast<char *>(CPLMalloc(nLen));
    1658           2 :             pszOverviewList[0] = '\0';
    1659           6 :             for (auto *poOverviewDS : m_apoOverviews)
    1660             :             {
    1661           4 :                 if (poOverviewDS)
    1662             :                 {
    1663             :                     const int nOvFactor = static_cast<int>(
    1664           4 :                         0.5 +
    1665           4 :                         GetRasterXSize() / static_cast<double>(
    1666           4 :                                                poOverviewDS->GetRasterXSize()));
    1667             : 
    1668           4 :                     snprintf(pszOverviewList + strlen(pszOverviewList),
    1669           4 :                              nLen - strlen(pszOverviewList), "%d ", nOvFactor);
    1670             :                 }
    1671             :             }
    1672             : 
    1673           2 :             CPLCreateXMLElementAndValue(psTree, "OverviewList",
    1674             :                                         pszOverviewList);
    1675             : 
    1676           2 :             CPLFree(pszOverviewList);
    1677             :         }
    1678             :     }
    1679             : 
    1680             :     /* -------------------------------------------------------------------- */
    1681             :     /*      Serialize source overview level.                                */
    1682             :     /* -------------------------------------------------------------------- */
    1683          78 :     if (m_nSrcOvrLevel != -2)
    1684             :     {
    1685           0 :         if (m_nSrcOvrLevel < -2)
    1686           0 :             CPLCreateXMLElementAndValue(
    1687             :                 psTree, "SrcOvrLevel",
    1688           0 :                 CPLSPrintf("AUTO%d", m_nSrcOvrLevel + 2));
    1689           0 :         else if (m_nSrcOvrLevel == -1)
    1690           0 :             CPLCreateXMLElementAndValue(psTree, "SrcOvrLevel", "NONE");
    1691             :         else
    1692           0 :             CPLCreateXMLElementAndValue(psTree, "SrcOvrLevel",
    1693             :                                         CPLSPrintf("%d", m_nSrcOvrLevel));
    1694             :     }
    1695             : 
    1696             :     /* ==================================================================== */
    1697             :     /*      Serialize the warp options.                                     */
    1698             :     /* ==================================================================== */
    1699          78 :     if (m_poWarper != nullptr)
    1700             :     {
    1701             :         /* --------------------------------------------------------------------
    1702             :          */
    1703             :         /*      We reset the destination dataset name so it doesn't get */
    1704             :         /*      written out in the serialize warp options. */
    1705             :         /* --------------------------------------------------------------------
    1706             :          */
    1707          78 :         char *const pszSavedName = CPLStrdup(GetDescription());
    1708          78 :         SetDescription("");
    1709             : 
    1710             :         CPLXMLNode *const psWOTree =
    1711          78 :             GDALSerializeWarpOptions(m_poWarper->GetOptions());
    1712          78 :         CPLAddXMLChild(psTree, psWOTree);
    1713             : 
    1714          78 :         SetDescription(pszSavedName);
    1715          78 :         CPLFree(pszSavedName);
    1716             : 
    1717             :         /* --------------------------------------------------------------------
    1718             :          */
    1719             :         /*      We need to consider making the source dataset relative to */
    1720             :         /*      the VRT file if possible.  Adjust accordingly. */
    1721             :         /* --------------------------------------------------------------------
    1722             :          */
    1723          78 :         CPLXMLNode *psSDS = CPLGetXMLNode(psWOTree, "SourceDataset");
    1724          78 :         int bRelativeToVRT = FALSE;
    1725             :         VSIStatBufL sStat;
    1726             : 
    1727          78 :         if (VSIStatExL(psSDS->psChild->pszValue, &sStat,
    1728          78 :                        VSI_STAT_EXISTS_FLAG) == 0)
    1729             :         {
    1730         152 :             std::string osVRTFilename = pszVRTPathIn;
    1731          76 :             std::string osSourceDataset = psSDS->psChild->pszValue;
    1732          76 :             char *pszCurDir = CPLGetCurrentDir();
    1733          76 :             if (CPLIsFilenameRelative(osSourceDataset.c_str()) &&
    1734          76 :                 !CPLIsFilenameRelative(osVRTFilename.c_str()) &&
    1735             :                 pszCurDir != nullptr)
    1736             :             {
    1737          48 :                 osSourceDataset = CPLFormFilenameSafe(
    1738          24 :                     pszCurDir, osSourceDataset.c_str(), nullptr);
    1739             :             }
    1740          52 :             else if (!CPLIsFilenameRelative(osSourceDataset.c_str()) &&
    1741          52 :                      CPLIsFilenameRelative(osVRTFilename.c_str()) &&
    1742             :                      pszCurDir != nullptr)
    1743             :             {
    1744           8 :                 osVRTFilename = CPLFormFilenameSafe(
    1745           4 :                     pszCurDir, osVRTFilename.c_str(), nullptr);
    1746             :             }
    1747          76 :             CPLFree(pszCurDir);
    1748          76 :             char *pszRelativePath = CPLStrdup(CPLExtractRelativePath(
    1749             :                 osVRTFilename.c_str(), osSourceDataset.c_str(),
    1750             :                 &bRelativeToVRT));
    1751             : 
    1752          76 :             CPLFree(psSDS->psChild->pszValue);
    1753          76 :             psSDS->psChild->pszValue = pszRelativePath;
    1754             :         }
    1755             : 
    1756          78 :         CPLCreateXMLNode(
    1757             :             CPLCreateXMLNode(psSDS, CXT_Attribute, "relativeToVRT"), CXT_Text,
    1758          78 :             bRelativeToVRT ? "1" : "0");
    1759             :     }
    1760             : 
    1761          78 :     return psTree;
    1762             : }
    1763             : 
    1764             : /************************************************************************/
    1765             : /*                            GetBlockSize()                            */
    1766             : /************************************************************************/
    1767             : 
    1768         683 : void VRTWarpedDataset::GetBlockSize(int *pnBlockXSize, int *pnBlockYSize) const
    1769             : 
    1770             : {
    1771         683 :     CPLAssert(nullptr != pnBlockXSize);
    1772         683 :     CPLAssert(nullptr != pnBlockYSize);
    1773             : 
    1774         683 :     *pnBlockXSize = m_nBlockXSize;
    1775         683 :     *pnBlockYSize = m_nBlockYSize;
    1776         683 : }
    1777             : 
    1778             : /************************************************************************/
    1779             : /*                            ProcessBlock()                            */
    1780             : /*                                                                      */
    1781             : /*      Warp a single requested block, and then push each band of       */
    1782             : /*      the result into the block cache.                                */
    1783             : /************************************************************************/
    1784             : 
    1785         243 : CPLErr VRTWarpedDataset::ProcessBlock(int iBlockX, int iBlockY)
    1786             : 
    1787             : {
    1788         243 :     if (m_poWarper == nullptr)
    1789           0 :         return CE_Failure;
    1790             : 
    1791         243 :     int nReqXSize = m_nBlockXSize;
    1792         243 :     if (iBlockX * m_nBlockXSize + nReqXSize > nRasterXSize)
    1793          62 :         nReqXSize = nRasterXSize - iBlockX * m_nBlockXSize;
    1794         243 :     int nReqYSize = m_nBlockYSize;
    1795         243 :     if (iBlockY * m_nBlockYSize + nReqYSize > nRasterYSize)
    1796          41 :         nReqYSize = nRasterYSize - iBlockY * m_nBlockYSize;
    1797             : 
    1798             :     GByte *pabyDstBuffer = static_cast<GByte *>(
    1799         243 :         m_poWarper->CreateDestinationBuffer(nReqXSize, nReqYSize));
    1800             : 
    1801         243 :     if (pabyDstBuffer == nullptr)
    1802             :     {
    1803           0 :         return CE_Failure;
    1804             :     }
    1805             : 
    1806             :     /* -------------------------------------------------------------------- */
    1807             :     /*      Warp into this buffer.                                          */
    1808             :     /* -------------------------------------------------------------------- */
    1809             : 
    1810         243 :     const GDALWarpOptions *psWO = m_poWarper->GetOptions();
    1811         486 :     const CPLErr eErr = m_poWarper->WarpRegionToBuffer(
    1812         243 :         iBlockX * m_nBlockXSize, iBlockY * m_nBlockYSize, nReqXSize, nReqYSize,
    1813         243 :         pabyDstBuffer, psWO->eWorkingDataType);
    1814             : 
    1815         243 :     if (eErr != CE_None)
    1816             :     {
    1817           0 :         m_poWarper->DestroyDestinationBuffer(pabyDstBuffer);
    1818           0 :         return eErr;
    1819             :     }
    1820             : 
    1821             :     /* -------------------------------------------------------------------- */
    1822             :     /*      Copy out into cache blocks for each band.                       */
    1823             :     /* -------------------------------------------------------------------- */
    1824         243 :     const int nWordSize = GDALGetDataTypeSizeBytes(psWO->eWorkingDataType);
    1825         879 :     for (int i = 0; i < psWO->nBandCount; i++)
    1826             :     {
    1827         636 :         int nDstBand = psWO->panDstBands[i];
    1828         636 :         if (GetRasterCount() < nDstBand)
    1829             :         {
    1830           0 :             continue;
    1831             :         }
    1832             : 
    1833         636 :         GDALRasterBand *poBand = GetRasterBand(nDstBand);
    1834             :         GDALRasterBlock *poBlock =
    1835         636 :             poBand->GetLockedBlockRef(iBlockX, iBlockY, TRUE);
    1836             : 
    1837         636 :         const GByte *pabyDstBandBuffer =
    1838             :             pabyDstBuffer +
    1839         636 :             static_cast<GPtrDiff_t>(i) * nReqXSize * nReqYSize * nWordSize;
    1840             : 
    1841         636 :         if (poBlock != nullptr)
    1842             :         {
    1843         636 :             if (poBlock->GetDataRef() != nullptr)
    1844             :             {
    1845         636 :                 if (nReqXSize == m_nBlockXSize && nReqYSize == m_nBlockYSize)
    1846             :                 {
    1847         413 :                     GDALCopyWords64(
    1848         413 :                         pabyDstBandBuffer, psWO->eWorkingDataType, nWordSize,
    1849             :                         poBlock->GetDataRef(), poBlock->GetDataType(),
    1850             :                         GDALGetDataTypeSizeBytes(poBlock->GetDataType()),
    1851         413 :                         static_cast<GPtrDiff_t>(m_nBlockXSize) * m_nBlockYSize);
    1852             :                 }
    1853             :                 else
    1854             :                 {
    1855             :                     GByte *pabyBlock =
    1856         223 :                         static_cast<GByte *>(poBlock->GetDataRef());
    1857             :                     const int nDTSize =
    1858         223 :                         GDALGetDataTypeSizeBytes(poBlock->GetDataType());
    1859       18716 :                     for (int iY = 0; iY < nReqYSize; iY++)
    1860             :                     {
    1861       18493 :                         GDALCopyWords(
    1862       18493 :                             pabyDstBandBuffer + static_cast<GPtrDiff_t>(iY) *
    1863       18493 :                                                     nReqXSize * nWordSize,
    1864       18493 :                             psWO->eWorkingDataType, nWordSize,
    1865       18493 :                             pabyBlock + static_cast<GPtrDiff_t>(iY) *
    1866       18493 :                                             m_nBlockXSize * nDTSize,
    1867             :                             poBlock->GetDataType(), nDTSize, nReqXSize);
    1868             :                     }
    1869             :                 }
    1870             :             }
    1871             : 
    1872         636 :             poBlock->DropLock();
    1873             :         }
    1874             :     }
    1875             : 
    1876         243 :     m_poWarper->DestroyDestinationBuffer(pabyDstBuffer);
    1877             : 
    1878         243 :     return CE_None;
    1879             : }
    1880             : 
    1881             : /************************************************************************/
    1882             : /*                              IRasterIO()                             */
    1883             : /************************************************************************/
    1884             : 
    1885             : // Specialized implementation of IRasterIO() that will be faster than
    1886             : // using the VRTWarpedRasterBand::IReadBlock() method in situations where
    1887             : // - a large enough chunk of data is requested at once
    1888             : // - and multi-threaded warping is enabled (it only kicks in if the warped
    1889             : //   chunk is large enough) and/or when reading the source dataset is
    1890             : //   multi-threaded (e.g JP2KAK or JP2OpenJPEG driver).
    1891        1018 : CPLErr VRTWarpedDataset::IRasterIO(
    1892             :     GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
    1893             :     void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
    1894             :     int nBandCount, BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
    1895             :     GSpacing nLineSpace, GSpacing nBandSpace, GDALRasterIOExtraArg *psExtraArg)
    1896             : {
    1897         995 :     const bool bWholeImage = nXOff == 0 && nYOff == 0 &&
    1898        2013 :                              nXSize == nRasterXSize && nYSize == nRasterYSize;
    1899             : 
    1900        2036 :     if (eRWFlag == GF_Write ||
    1901             :         // For too small request fall back to the block-based approach to
    1902             :         // benefit from caching
    1903        1018 :         (!bWholeImage &&
    1904         813 :          (nBufXSize <= m_nBlockXSize || nBufYSize <= m_nBlockYSize)) ||
    1905             :         // Or if we don't request all bands at once
    1906        2226 :         nBandCount < nBands ||
    1907         190 :         !CPLTestBool(
    1908             :             CPLGetConfigOption("GDAL_VRT_WARP_USE_DATASET_RASTERIO", "YES")))
    1909             :     {
    1910         836 :         return GDALDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
    1911             :                                       pData, nBufXSize, nBufYSize, eBufType,
    1912             :                                       nBandCount, panBandMap, nPixelSpace,
    1913         836 :                                       nLineSpace, nBandSpace, psExtraArg);
    1914             :     }
    1915             : 
    1916             :     // Try overviews for sub-sampled requests
    1917         182 :     if (nBufXSize < nXSize || nBufYSize < nYSize)
    1918             :     {
    1919           5 :         int bTried = FALSE;
    1920           5 :         const CPLErr eErr = TryOverviewRasterIO(
    1921             :             eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
    1922             :             eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
    1923             :             nBandSpace, psExtraArg, &bTried);
    1924             : 
    1925           5 :         if (bTried)
    1926             :         {
    1927           3 :             return eErr;
    1928             :         }
    1929             :     }
    1930             : 
    1931         179 :     if (m_poWarper == nullptr)
    1932           0 :         return CE_Failure;
    1933             : 
    1934         179 :     const GDALWarpOptions *psWO = m_poWarper->GetOptions();
    1935             : 
    1936         179 :     if (nBufXSize != nXSize || nBufYSize != nYSize)
    1937             :     {
    1938           2 :         if (!bWholeImage || !GDALTransformHasFastClone(psWO->pTransformerArg))
    1939             :         {
    1940           0 :             return GDALDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
    1941             :                                           pData, nBufXSize, nBufYSize, eBufType,
    1942             :                                           nBandCount, panBandMap, nPixelSpace,
    1943           0 :                                           nLineSpace, nBandSpace, psExtraArg);
    1944             :         }
    1945             : 
    1946             :         // Build a temporary dataset taking into account the rescaling
    1947           2 :         void *pTransformerArg = GDALCloneTransformer(psWO->pTransformerArg);
    1948           2 :         if (pTransformerArg == nullptr)
    1949             :         {
    1950           0 :             return GDALDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
    1951             :                                           pData, nBufXSize, nBufYSize, eBufType,
    1952             :                                           nBandCount, panBandMap, nPixelSpace,
    1953           0 :                                           nLineSpace, nBandSpace, psExtraArg);
    1954             :         }
    1955             : 
    1956           2 :         GDALWarpOptions *psRescaledWO = GDALCloneWarpOptions(psWO);
    1957           2 :         psRescaledWO->hSrcDS = psWO->hSrcDS;
    1958           2 :         psRescaledWO->pfnTransformer = psWO->pfnTransformer;
    1959           2 :         psRescaledWO->pTransformerArg = pTransformerArg;
    1960             : 
    1961             :         // Rescale the output geotransform on the transformer.
    1962           2 :         double adfDstGeoTransform[6] = {0.0};
    1963           2 :         GDALGetTransformerDstGeoTransform(psRescaledWO->pTransformerArg,
    1964             :                                           adfDstGeoTransform);
    1965           2 :         RescaleDstGeoTransform(adfDstGeoTransform, nRasterXSize, nBufXSize,
    1966             :                                nRasterYSize, nBufYSize);
    1967           2 :         GDALSetTransformerDstGeoTransform(psRescaledWO->pTransformerArg,
    1968             :                                           adfDstGeoTransform);
    1969             : 
    1970             :         GDALDatasetH hDstDS =
    1971           2 :             GDALCreateWarpedVRT(psWO->hSrcDS, nBufXSize, nBufYSize,
    1972             :                                 adfDstGeoTransform, psRescaledWO);
    1973             : 
    1974           2 :         GDALDestroyWarpOptions(psRescaledWO);
    1975             : 
    1976           2 :         if (hDstDS == nullptr)
    1977             :         {
    1978             :             // Not supposed to happen in nominal circumstances. Could perhaps
    1979             :             // happen if some memory allocation error occurred in code called
    1980             :             // by GDALCreateWarpedVRT()
    1981           0 :             GDALDestroyTransformer(pTransformerArg);
    1982           0 :             return GDALDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
    1983             :                                           pData, nBufXSize, nBufYSize, eBufType,
    1984             :                                           nBandCount, panBandMap, nPixelSpace,
    1985           0 :                                           nLineSpace, nBandSpace, psExtraArg);
    1986             :         }
    1987             : 
    1988           2 :         auto poOvrDS = static_cast<VRTWarpedDataset *>(hDstDS);
    1989           2 :         poOvrDS->m_bIsOverview = true;
    1990             : 
    1991             :         GDALRasterIOExtraArg sExtraArg;
    1992           2 :         INIT_RASTERIO_EXTRA_ARG(sExtraArg);
    1993           2 :         CPLErr eErr = poOvrDS->IRasterIO(GF_Read, 0, 0, nBufXSize, nBufYSize,
    1994             :                                          pData, nBufXSize, nBufYSize, eBufType,
    1995             :                                          nBandCount, panBandMap, nPixelSpace,
    1996             :                                          nLineSpace, nBandSpace, &sExtraArg);
    1997             : 
    1998           2 :         poOvrDS->ReleaseRef();
    1999           2 :         return eErr;
    2000             :     }
    2001             : 
    2002             :     // Build a map from warped output bands to their index
    2003         354 :     std::map<int, int> oMapBandToWarpingBandIndex;
    2004         177 :     bool bAllBandsIncreasingOrder =
    2005         177 :         (psWO->nBandCount == nBands && nBands == nBandCount);
    2006         374 :     for (int i = 0; i < psWO->nBandCount; ++i)
    2007             :     {
    2008         197 :         oMapBandToWarpingBandIndex[psWO->panDstBands[i]] = i;
    2009         197 :         if (psWO->panDstBands[i] != i + 1 || panBandMap[i] != i + 1)
    2010             :         {
    2011           3 :             bAllBandsIncreasingOrder = false;
    2012             :         }
    2013             :     }
    2014             : 
    2015             :     // Check that all requested bands are actually warped output bands.
    2016         374 :     for (int i = 0; i < nBandCount; ++i)
    2017             :     {
    2018         200 :         const int nRasterIOBand = panBandMap[i];
    2019         200 :         if (oMapBandToWarpingBandIndex.find(nRasterIOBand) ==
    2020         400 :             oMapBandToWarpingBandIndex.end())
    2021             :         {
    2022             :             // Not sure if that can happen...
    2023             :             // but if that does, that will likely later fail in ProcessBlock()
    2024           3 :             return GDALDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
    2025             :                                           pData, nBufXSize, nBufYSize, eBufType,
    2026             :                                           nBandCount, panBandMap, nPixelSpace,
    2027           3 :                                           nLineSpace, nBandSpace, psExtraArg);
    2028             :         }
    2029             :     }
    2030             : 
    2031         174 :     int nSrcXOff = 0;
    2032         174 :     int nSrcYOff = 0;
    2033         174 :     int nSrcXSize = 0;
    2034         174 :     int nSrcYSize = 0;
    2035         174 :     double dfSrcXExtraSize = 0;
    2036         174 :     double dfSrcYExtraSize = 0;
    2037         174 :     double dfSrcFillRatio = 0;
    2038             :     // Find the source window that corresponds to our target window
    2039         174 :     if (m_poWarper->ComputeSourceWindow(nXOff, nYOff, nXSize, nYSize, &nSrcXOff,
    2040             :                                         &nSrcYOff, &nSrcXSize, &nSrcYSize,
    2041             :                                         &dfSrcXExtraSize, &dfSrcYExtraSize,
    2042         174 :                                         &dfSrcFillRatio) != CE_None)
    2043             :     {
    2044           0 :         return CE_Failure;
    2045             :     }
    2046             : 
    2047         174 :     GByte *const pabyDst = static_cast<GByte *>(pData);
    2048         174 :     const int nWarpDTSize = GDALGetDataTypeSizeBytes(psWO->eWorkingDataType);
    2049             : 
    2050         174 :     const double dfMemRequired = m_poWarper->GetWorkingMemoryForWindow(
    2051             :         nSrcXSize, nSrcYSize, nXSize, nYSize);
    2052             :     // If we need more warp working memory than allowed, we have to use a
    2053             :     // splitting strategy until we get below the limit.
    2054         174 :     if (dfMemRequired > psWO->dfWarpMemoryLimit && nXSize >= 2 && nYSize >= 2)
    2055             :     {
    2056           3 :         CPLDebugOnly("VRT", "VRTWarpedDataset::IRasterIO(): exceeding warp "
    2057             :                             "memory. Splitting region");
    2058             : 
    2059             :         GDALRasterIOExtraArg sExtraArg;
    2060           3 :         INIT_RASTERIO_EXTRA_ARG(sExtraArg);
    2061             : 
    2062             :         bool bOK;
    2063             :         // Split along the longest dimension
    2064           3 :         if (nXSize >= nYSize)
    2065             :         {
    2066           2 :             const int nHalfXSize = nXSize / 2;
    2067           4 :             bOK = IRasterIO(GF_Read, nXOff, nYOff, nHalfXSize, nYSize, pabyDst,
    2068             :                             nHalfXSize, nYSize, eBufType, nBandCount,
    2069             :                             panBandMap, nPixelSpace, nLineSpace, nBandSpace,
    2070           4 :                             &sExtraArg) == CE_None &&
    2071           2 :                   IRasterIO(GF_Read, nXOff + nHalfXSize, nYOff,
    2072             :                             nXSize - nHalfXSize, nYSize,
    2073           2 :                             pabyDst + nHalfXSize * nPixelSpace,
    2074             :                             nXSize - nHalfXSize, nYSize, eBufType, nBandCount,
    2075             :                             panBandMap, nPixelSpace, nLineSpace, nBandSpace,
    2076             :                             &sExtraArg) == CE_None;
    2077             :         }
    2078             :         else
    2079             :         {
    2080           1 :             const int nHalfYSize = nYSize / 2;
    2081           2 :             bOK = IRasterIO(GF_Read, nXOff, nYOff, nXSize, nHalfYSize, pabyDst,
    2082             :                             nXSize, nHalfYSize, eBufType, nBandCount,
    2083             :                             panBandMap, nPixelSpace, nLineSpace, nBandSpace,
    2084           2 :                             &sExtraArg) == CE_None &&
    2085           1 :                   IRasterIO(GF_Read, nXOff, nYOff + nHalfYSize, nXSize,
    2086             :                             nYSize - nHalfYSize,
    2087           1 :                             pabyDst + nHalfYSize * nLineSpace, nXSize,
    2088             :                             nYSize - nHalfYSize, eBufType, nBandCount,
    2089             :                             panBandMap, nPixelSpace, nLineSpace, nBandSpace,
    2090             :                             &sExtraArg) == CE_None;
    2091             :         }
    2092           3 :         return bOK ? CE_None : CE_Failure;
    2093             :     }
    2094             : 
    2095         171 :     CPLDebugOnly("VRT",
    2096             :                  "Using optimized VRTWarpedDataset::IRasterIO() code path");
    2097             : 
    2098             :     // Allocate a warping destination buffer if needed.
    2099             :     // We can use directly the output buffer pData if:
    2100             :     // - we request exactly all warped bands, and that there are as many
    2101             :     //   warped bands as dataset bands (no alpha)
    2102             :     // - the output buffer data atype is the warping working data type
    2103             :     // - the output buffer has a band-sequential layout.
    2104             :     GByte *pabyWarpBuffer;
    2105             : 
    2106         169 :     if (bAllBandsIncreasingOrder && psWO->eWorkingDataType == eBufType &&
    2107          85 :         nPixelSpace == GDALGetDataTypeSizeBytes(eBufType) &&
    2108         413 :         nLineSpace == nPixelSpace * nXSize &&
    2109          73 :         (nBands == 1 || nBandSpace == nLineSpace * nYSize))
    2110             :     {
    2111          73 :         pabyWarpBuffer = static_cast<GByte *>(pData);
    2112          73 :         m_poWarper->InitializeDestinationBuffer(pabyWarpBuffer, nXSize, nYSize);
    2113             :     }
    2114             :     else
    2115             :     {
    2116             :         pabyWarpBuffer = static_cast<GByte *>(
    2117          98 :             m_poWarper->CreateDestinationBuffer(nXSize, nYSize));
    2118             : 
    2119          98 :         if (pabyWarpBuffer == nullptr)
    2120             :         {
    2121           0 :             return CE_Failure;
    2122             :         }
    2123             :     }
    2124             : 
    2125         342 :     const CPLErr eErr = m_poWarper->WarpRegionToBuffer(
    2126         171 :         nXOff, nYOff, nXSize, nYSize, pabyWarpBuffer, psWO->eWorkingDataType,
    2127             :         nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize, dfSrcXExtraSize,
    2128             :         dfSrcYExtraSize);
    2129             : 
    2130         171 :     if (pabyWarpBuffer != pData)
    2131             :     {
    2132          98 :         if (eErr == CE_None)
    2133             :         {
    2134             :             // Copy warping buffer into user destination buffer
    2135         204 :             for (int i = 0; i < nBandCount; i++)
    2136             :             {
    2137         106 :                 const int nRasterIOBand = panBandMap[i];
    2138             :                 const auto oIterToWarpingBandIndex =
    2139         106 :                     oMapBandToWarpingBandIndex.find(nRasterIOBand);
    2140             :                 // cannot happen due to earlier check
    2141         106 :                 CPLAssert(oIterToWarpingBandIndex !=
    2142             :                           oMapBandToWarpingBandIndex.end());
    2143             : 
    2144             :                 const GByte *const pabyWarpBandBuffer =
    2145             :                     pabyWarpBuffer +
    2146         106 :                     static_cast<GPtrDiff_t>(oIterToWarpingBandIndex->second) *
    2147         106 :                         nXSize * nYSize * nWarpDTSize;
    2148         106 :                 GByte *const pabyDstBand = pabyDst + i * nBandSpace;
    2149             : 
    2150       17277 :                 for (int iY = 0; iY < nYSize; iY++)
    2151             :                 {
    2152       17171 :                     GDALCopyWords(pabyWarpBandBuffer +
    2153       17171 :                                       static_cast<GPtrDiff_t>(iY) * nXSize *
    2154       17171 :                                           nWarpDTSize,
    2155       17171 :                                   psWO->eWorkingDataType, nWarpDTSize,
    2156       17171 :                                   pabyDstBand + iY * nLineSpace, eBufType,
    2157             :                                   static_cast<int>(nPixelSpace), nXSize);
    2158             :                 }
    2159             :             }
    2160             :         }
    2161             : 
    2162          98 :         m_poWarper->DestroyDestinationBuffer(pabyWarpBuffer);
    2163             :     }
    2164             : 
    2165         171 :     return eErr;
    2166             : }
    2167             : 
    2168             : /************************************************************************/
    2169             : /*                              AddBand()                               */
    2170             : /************************************************************************/
    2171             : 
    2172         229 : CPLErr VRTWarpedDataset::AddBand(GDALDataType eType, char ** /* papszOptions */)
    2173             : 
    2174             : {
    2175         229 :     if (eType == GDT_Unknown || eType == GDT_TypeCount)
    2176             :     {
    2177           1 :         ReportError(CE_Failure, CPLE_IllegalArg,
    2178             :                     "Illegal GDT_Unknown/GDT_TypeCount argument");
    2179           1 :         return CE_Failure;
    2180             :     }
    2181             : 
    2182         228 :     SetBand(GetRasterCount() + 1,
    2183         228 :             new VRTWarpedRasterBand(this, GetRasterCount() + 1, eType));
    2184             : 
    2185         228 :     return CE_None;
    2186             : }
    2187             : 
    2188             : /************************************************************************/
    2189             : /* ==================================================================== */
    2190             : /*                        VRTWarpedRasterBand                           */
    2191             : /* ==================================================================== */
    2192             : /************************************************************************/
    2193             : 
    2194             : /************************************************************************/
    2195             : /*                        VRTWarpedRasterBand()                         */
    2196             : /************************************************************************/
    2197             : 
    2198         683 : VRTWarpedRasterBand::VRTWarpedRasterBand(GDALDataset *poDSIn, int nBandIn,
    2199         683 :                                          GDALDataType eType)
    2200             : 
    2201             : {
    2202         683 :     Initialize(poDSIn->GetRasterXSize(), poDSIn->GetRasterYSize());
    2203             : 
    2204         683 :     poDS = poDSIn;
    2205         683 :     nBand = nBandIn;
    2206         683 :     eAccess = GA_Update;
    2207             : 
    2208         683 :     static_cast<VRTWarpedDataset *>(poDS)->GetBlockSize(&nBlockXSize,
    2209             :                                                         &nBlockYSize);
    2210             : 
    2211         683 :     if (eType != GDT_Unknown)
    2212         232 :         eDataType = eType;
    2213         683 : }
    2214             : 
    2215             : /************************************************************************/
    2216             : /*                        ~VRTWarpedRasterBand()                        */
    2217             : /************************************************************************/
    2218             : 
    2219        1366 : VRTWarpedRasterBand::~VRTWarpedRasterBand()
    2220             : 
    2221             : {
    2222         683 :     FlushCache(true);
    2223        1366 : }
    2224             : 
    2225             : /************************************************************************/
    2226             : /*                             IReadBlock()                             */
    2227             : /************************************************************************/
    2228             : 
    2229         243 : CPLErr VRTWarpedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
    2230             :                                        void *pImage)
    2231             : 
    2232             : {
    2233         243 :     VRTWarpedDataset *poWDS = static_cast<VRTWarpedDataset *>(poDS);
    2234             :     const GPtrDiff_t nDataBytes =
    2235         243 :         static_cast<GPtrDiff_t>(GDALGetDataTypeSizeBytes(eDataType)) *
    2236         243 :         nBlockXSize * nBlockYSize;
    2237             : 
    2238         243 :     GDALRasterBlock *poBlock = GetLockedBlockRef(nBlockXOff, nBlockYOff, TRUE);
    2239         243 :     if (poBlock == nullptr)
    2240           0 :         return CE_Failure;
    2241             : 
    2242         243 :     if (poWDS->m_poWarper)
    2243             :     {
    2244         243 :         const GDALWarpOptions *psWO = poWDS->m_poWarper->GetOptions();
    2245         243 :         if (nBand == psWO->nDstAlphaBand)
    2246             :         {
    2247             :             // For a reader starting by asking on band 1, we should normally
    2248             :             // not reach here, because ProcessBlock() on band 1 will have
    2249             :             // populated the block cache for the regular bands and the alpha
    2250             :             // band.
    2251             :             // But if there's no source window corresponding to the block,
    2252             :             // the alpha band block will not be written through RasterIO(),
    2253             :             // so we nee to initialize it.
    2254          88 :             memset(poBlock->GetDataRef(), 0, nDataBytes);
    2255             :         }
    2256             :     }
    2257             : 
    2258         243 :     const CPLErr eErr = poWDS->ProcessBlock(nBlockXOff, nBlockYOff);
    2259             : 
    2260         243 :     if (eErr == CE_None && pImage != poBlock->GetDataRef())
    2261             :     {
    2262           0 :         memcpy(pImage, poBlock->GetDataRef(), nDataBytes);
    2263             :     }
    2264             : 
    2265         243 :     poBlock->DropLock();
    2266             : 
    2267         243 :     return eErr;
    2268             : }
    2269             : 
    2270             : /************************************************************************/
    2271             : /*                            IWriteBlock()                             */
    2272             : /************************************************************************/
    2273             : 
    2274          89 : CPLErr VRTWarpedRasterBand::IWriteBlock(int nBlockXOff, int nBlockYOff,
    2275             :                                         void *pImage)
    2276             : 
    2277             : {
    2278          89 :     VRTWarpedDataset *poWDS = static_cast<VRTWarpedDataset *>(poDS);
    2279             : 
    2280             :     // This is a bit tricky. In the case we are warping a VRTWarpedDataset
    2281             :     // with a destination alpha band, IWriteBlock can be called on that alpha
    2282             :     // band by GDALWarpDstAlphaMasker
    2283             :     // We don't need to do anything since the data will have hopefully been
    2284             :     // read from the block cache before if the reader processes all the bands
    2285             :     // of a same block.
    2286          89 :     if (poWDS->m_poWarper->GetOptions()->nDstAlphaBand != nBand)
    2287             :     {
    2288             :         /* Otherwise, call the superclass method, that will fail of course */
    2289           0 :         return VRTRasterBand::IWriteBlock(nBlockXOff, nBlockYOff, pImage);
    2290             :     }
    2291             : 
    2292          89 :     return CE_None;
    2293             : }
    2294             : 
    2295             : /************************************************************************/
    2296             : /*                   EmitErrorMessageIfWriteNotSupported()              */
    2297             : /************************************************************************/
    2298             : 
    2299          91 : bool VRTWarpedRasterBand::EmitErrorMessageIfWriteNotSupported(
    2300             :     const char *pszCaller) const
    2301             : {
    2302          91 :     VRTWarpedDataset *poWDS = static_cast<VRTWarpedDataset *>(poDS);
    2303             :     // Cf comment in IWriteBlock()
    2304          91 :     if (poWDS->m_poWarper->GetOptions()->nDstAlphaBand != nBand)
    2305             :     {
    2306           2 :         ReportError(CE_Failure, CPLE_NoWriteAccess,
    2307             :                     "%s: attempt to write to a VRTWarpedRasterBand.",
    2308             :                     pszCaller);
    2309             : 
    2310           2 :         return true;
    2311             :     }
    2312          89 :     return false;
    2313             : }
    2314             : 
    2315             : /************************************************************************/
    2316             : /*                       GetBestOverviewLevel()                         */
    2317             : /************************************************************************/
    2318             : 
    2319           0 : int VRTWarpedRasterBand::GetBestOverviewLevel(
    2320             :     int &nXOff, int &nYOff, int &nXSize, int &nYSize, int nBufXSize,
    2321             :     int nBufYSize, GDALRasterIOExtraArg *psExtraArg) const
    2322             : {
    2323           0 :     VRTWarpedDataset *poWDS = static_cast<VRTWarpedDataset *>(poDS);
    2324             : 
    2325             :     /* -------------------------------------------------------------------- */
    2326             :     /*      Compute the desired downsampling factor.  It is                 */
    2327             :     /*      based on the least reduced axis, and represents the number      */
    2328             :     /*      of source pixels to one destination pixel.                      */
    2329             :     /* -------------------------------------------------------------------- */
    2330           0 :     const double dfDesiredDownsamplingFactor =
    2331           0 :         ((nXSize / static_cast<double>(nBufXSize)) <
    2332           0 :              (nYSize / static_cast<double>(nBufYSize)) ||
    2333             :          nBufYSize == 1)
    2334           0 :             ? nXSize / static_cast<double>(nBufXSize)
    2335           0 :             : nYSize / static_cast<double>(nBufYSize);
    2336             : 
    2337             :     /* -------------------------------------------------------------------- */
    2338             :     /*      Find the overview level that largest downsampling factor (most  */
    2339             :     /*      downsampled) that is still less than (or only a little more)    */
    2340             :     /*      downsampled than the request.                                   */
    2341             :     /* -------------------------------------------------------------------- */
    2342           0 :     const GDALWarpOptions *psWO = poWDS->m_poWarper->GetOptions();
    2343           0 :     GDALDataset *poSrcDS = GDALDataset::FromHandle(psWO->hSrcDS);
    2344           0 :     const int nOverviewCount = poSrcDS->GetRasterBand(1)->GetOverviewCount();
    2345             : 
    2346           0 :     int nBestOverviewXSize = 1;
    2347           0 :     int nBestOverviewYSize = 1;
    2348           0 :     double dfBestDownsamplingFactor = 0;
    2349           0 :     int nBestOverviewLevel = -1;
    2350             : 
    2351             :     const char *pszOversampligThreshold =
    2352           0 :         CPLGetConfigOption("GDAL_OVERVIEW_OVERSAMPLING_THRESHOLD", nullptr);
    2353             : 
    2354             :     // Cf https://github.com/OSGeo/gdal/pull/9040#issuecomment-1898524693
    2355             :     // Do not exactly use a oversampling threshold of 1.0 because of numerical
    2356             :     // instability.
    2357           0 :     const auto AdjustThreshold = [](double x)
    2358             :     {
    2359           0 :         constexpr double EPS = 1e-2;
    2360           0 :         return x == 1.0 ? x + EPS : x;
    2361             :     };
    2362           0 :     const double dfOversamplingThreshold = AdjustThreshold(
    2363           0 :         pszOversampligThreshold ? CPLAtof(pszOversampligThreshold)
    2364           0 :         : psExtraArg && psExtraArg->eResampleAlg != GRIORA_NearestNeighbour
    2365           0 :             ? 1.0
    2366             :             : 1.2);
    2367           0 :     for (int iOverview = 0; iOverview < nOverviewCount; iOverview++)
    2368             :     {
    2369           0 :         const GDALRasterBand *poSrcOvrBand = this;
    2370           0 :         bool bThisLevelOnly = false;
    2371             :         const int iSrcOvr =
    2372           0 :             poWDS->GetSrcOverviewLevel(iOverview, bThisLevelOnly);
    2373           0 :         if (iSrcOvr >= 0)
    2374             :         {
    2375           0 :             poSrcOvrBand = poSrcDS->GetRasterBand(1)->GetOverview(iSrcOvr);
    2376             :         }
    2377           0 :         if (poSrcOvrBand == nullptr)
    2378           0 :             break;
    2379             : 
    2380           0 :         int nDstPixels = 0;
    2381           0 :         int nDstLines = 0;
    2382           0 :         double dfSrcRatioX = 0;
    2383           0 :         double dfSrcRatioY = 0;
    2384           0 :         if (!poWDS->GetOverviewSize(poSrcDS, iOverview, iSrcOvr, nDstPixels,
    2385             :                                     nDstLines, dfSrcRatioX, dfSrcRatioY))
    2386             :         {
    2387           0 :             break;
    2388             :         }
    2389             : 
    2390             :         // Compute downsampling factor of this overview
    2391             :         const double dfDownsamplingFactor =
    2392           0 :             std::min(nRasterXSize / static_cast<double>(nDstPixels),
    2393           0 :                      nRasterYSize / static_cast<double>(nDstLines));
    2394             : 
    2395             :         // Is it nearly the requested factor and better (lower) than
    2396             :         // the current best factor?
    2397           0 :         if (dfDownsamplingFactor >=
    2398           0 :                 dfDesiredDownsamplingFactor * dfOversamplingThreshold ||
    2399             :             dfDownsamplingFactor <= dfBestDownsamplingFactor)
    2400             :         {
    2401           0 :             continue;
    2402             :         }
    2403             : 
    2404             :         // Ignore AVERAGE_BIT2GRAYSCALE overviews for RasterIO purposes.
    2405             :         const char *pszResampling = const_cast<GDALRasterBand *>(poSrcOvrBand)
    2406           0 :                                         ->GetMetadataItem("RESAMPLING");
    2407             : 
    2408           0 :         if (pszResampling != nullptr &&
    2409           0 :             STARTS_WITH_CI(pszResampling, "AVERAGE_BIT2"))
    2410           0 :             continue;
    2411             : 
    2412             :         // OK, this is our new best overview.
    2413           0 :         nBestOverviewXSize = nDstPixels;
    2414           0 :         nBestOverviewYSize = nDstLines;
    2415           0 :         nBestOverviewLevel = iOverview;
    2416           0 :         dfBestDownsamplingFactor = dfDownsamplingFactor;
    2417             :     }
    2418             : 
    2419             :     /* -------------------------------------------------------------------- */
    2420             :     /*      If we didn't find an overview that helps us, just return        */
    2421             :     /*      indicating failure and the full resolution image will be used.  */
    2422             :     /* -------------------------------------------------------------------- */
    2423           0 :     if (nBestOverviewLevel < 0)
    2424           0 :         return -1;
    2425             : 
    2426             :     /* -------------------------------------------------------------------- */
    2427             :     /*      Recompute the source window in terms of the selected            */
    2428             :     /*      overview.                                                       */
    2429             :     /* -------------------------------------------------------------------- */
    2430           0 :     const double dfXFactor =
    2431           0 :         nRasterXSize / static_cast<double>(nBestOverviewXSize);
    2432           0 :     const double dfYFactor =
    2433           0 :         nRasterYSize / static_cast<double>(nBestOverviewYSize);
    2434           0 :     CPLDebug("GDAL", "Selecting overview %d x %d", nBestOverviewXSize,
    2435             :              nBestOverviewYSize);
    2436             : 
    2437           0 :     const int nOXOff = std::min(nBestOverviewXSize - 1,
    2438           0 :                                 static_cast<int>(nXOff / dfXFactor + 0.5));
    2439           0 :     const int nOYOff = std::min(nBestOverviewYSize - 1,
    2440           0 :                                 static_cast<int>(nYOff / dfYFactor + 0.5));
    2441           0 :     int nOXSize = std::max(1, static_cast<int>(nXSize / dfXFactor + 0.5));
    2442           0 :     int nOYSize = std::max(1, static_cast<int>(nYSize / dfYFactor + 0.5));
    2443           0 :     if (nOXOff + nOXSize > nBestOverviewXSize)
    2444           0 :         nOXSize = nBestOverviewXSize - nOXOff;
    2445           0 :     if (nOYOff + nOYSize > nBestOverviewYSize)
    2446           0 :         nOYSize = nBestOverviewYSize - nOYOff;
    2447             : 
    2448           0 :     if (psExtraArg)
    2449             :     {
    2450           0 :         if (psExtraArg->bFloatingPointWindowValidity)
    2451             :         {
    2452           0 :             psExtraArg->dfXOff /= dfXFactor;
    2453           0 :             psExtraArg->dfXSize /= dfXFactor;
    2454           0 :             psExtraArg->dfYOff /= dfYFactor;
    2455           0 :             psExtraArg->dfYSize /= dfYFactor;
    2456             :         }
    2457           0 :         else if (psExtraArg->eResampleAlg != GRIORA_NearestNeighbour)
    2458             :         {
    2459           0 :             psExtraArg->bFloatingPointWindowValidity = true;
    2460           0 :             psExtraArg->dfXOff = nXOff / dfXFactor;
    2461           0 :             psExtraArg->dfXSize = nXSize / dfXFactor;
    2462           0 :             psExtraArg->dfYOff = nYOff / dfYFactor;
    2463           0 :             psExtraArg->dfYSize = nYSize / dfYFactor;
    2464             :         }
    2465             :     }
    2466             : 
    2467           0 :     nXOff = nOXOff;
    2468           0 :     nYOff = nOYOff;
    2469           0 :     nXSize = nOXSize;
    2470           0 :     nYSize = nOYSize;
    2471             : 
    2472           0 :     return nBestOverviewLevel;
    2473             : }
    2474             : 
    2475             : /************************************************************************/
    2476             : /*                              IRasterIO()                             */
    2477             : /************************************************************************/
    2478             : 
    2479        4127 : CPLErr VRTWarpedRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
    2480             :                                       int nXSize, int nYSize, void *pData,
    2481             :                                       int nBufXSize, int nBufYSize,
    2482             :                                       GDALDataType eBufType,
    2483             :                                       GSpacing nPixelSpace, GSpacing nLineSpace,
    2484             :                                       GDALRasterIOExtraArg *psExtraArg)
    2485             : {
    2486        4127 :     VRTWarpedDataset *poWDS = static_cast<VRTWarpedDataset *>(poDS);
    2487        4127 :     if (m_nIRasterIOCounter == 0 && poWDS->GetRasterCount() == 1)
    2488             :     {
    2489         174 :         int anBandMap[] = {nBand};
    2490         174 :         ++m_nIRasterIOCounter;
    2491         174 :         const CPLErr eErr = poWDS->IRasterIO(
    2492             :             eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
    2493             :             eBufType, 1, anBandMap, nPixelSpace, nLineSpace, 0, psExtraArg);
    2494         174 :         --m_nIRasterIOCounter;
    2495         174 :         return eErr;
    2496             :     }
    2497             : 
    2498             :     /* ==================================================================== */
    2499             :     /*      Do we have overviews that would be appropriate to satisfy       */
    2500             :     /*      this request?                                                   */
    2501             :     /* ==================================================================== */
    2502        3953 :     if ((nBufXSize < nXSize || nBufYSize < nYSize) && GetOverviewCount() &&
    2503             :         eRWFlag == GF_Read)
    2504             :     {
    2505             :         GDALRasterIOExtraArg sExtraArg;
    2506           0 :         GDALCopyRasterIOExtraArg(&sExtraArg, psExtraArg);
    2507             : 
    2508           0 :         const int nOverview = GetBestOverviewLevel(
    2509             :             nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize, &sExtraArg);
    2510           0 :         if (nOverview >= 0)
    2511             :         {
    2512           0 :             auto poOvrBand = GetOverview(nOverview);
    2513           0 :             if (!poOvrBand)
    2514           0 :                 return CE_Failure;
    2515             : 
    2516           0 :             return poOvrBand->RasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
    2517             :                                        pData, nBufXSize, nBufYSize, eBufType,
    2518           0 :                                        nPixelSpace, nLineSpace, &sExtraArg);
    2519             :         }
    2520             :     }
    2521             : 
    2522        3953 :     return GDALRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
    2523             :                                      pData, nBufXSize, nBufYSize, eBufType,
    2524        3953 :                                      nPixelSpace, nLineSpace, psExtraArg);
    2525             : }
    2526             : 
    2527             : /************************************************************************/
    2528             : /*                           SerializeToXML()                           */
    2529             : /************************************************************************/
    2530             : 
    2531         187 : CPLXMLNode *VRTWarpedRasterBand::SerializeToXML(const char *pszVRTPathIn,
    2532             :                                                 bool &bHasWarnedAboutRAMUsage,
    2533             :                                                 size_t &nAccRAMUsage)
    2534             : 
    2535             : {
    2536         187 :     CPLXMLNode *const psTree = VRTRasterBand::SerializeToXML(
    2537             :         pszVRTPathIn, bHasWarnedAboutRAMUsage, nAccRAMUsage);
    2538             : 
    2539             :     /* -------------------------------------------------------------------- */
    2540             :     /*      Set subclass.                                                   */
    2541             :     /* -------------------------------------------------------------------- */
    2542         187 :     CPLCreateXMLNode(CPLCreateXMLNode(psTree, CXT_Attribute, "subClass"),
    2543             :                      CXT_Text, "VRTWarpedRasterBand");
    2544             : 
    2545         187 :     return psTree;
    2546             : }
    2547             : 
    2548             : /************************************************************************/
    2549             : /*                          GetOverviewCount()                          */
    2550             : /************************************************************************/
    2551             : 
    2552          86 : int VRTWarpedRasterBand::GetOverviewCount()
    2553             : 
    2554             : {
    2555          86 :     VRTWarpedDataset *const poWDS = static_cast<VRTWarpedDataset *>(poDS);
    2556          86 :     if (poWDS->m_bIsOverview)
    2557           1 :         return 0;
    2558             : 
    2559          85 :     if (poWDS->m_apoOverviews.empty())
    2560             :     {
    2561          57 :         return poWDS->GetOverviewCount();
    2562             :     }
    2563             : 
    2564          28 :     return static_cast<int>(poWDS->m_apoOverviews.size());
    2565             : }
    2566             : 
    2567             : /************************************************************************/
    2568             : /*                            GetOverview()                             */
    2569             : /************************************************************************/
    2570             : 
    2571          41 : GDALRasterBand *VRTWarpedRasterBand::GetOverview(int iOverview)
    2572             : 
    2573             : {
    2574          41 :     VRTWarpedDataset *const poWDS = static_cast<VRTWarpedDataset *>(poDS);
    2575             : 
    2576          41 :     const int nOvrCount = GetOverviewCount();
    2577          41 :     if (iOverview < 0 || iOverview >= nOvrCount)
    2578           6 :         return nullptr;
    2579             : 
    2580          35 :     if (poWDS->m_apoOverviews.empty())
    2581          11 :         poWDS->m_apoOverviews.resize(nOvrCount);
    2582          35 :     if (!poWDS->m_apoOverviews[iOverview])
    2583          18 :         poWDS->m_apoOverviews[iOverview] =
    2584          18 :             poWDS->CreateImplicitOverview(iOverview);
    2585          35 :     if (!poWDS->m_apoOverviews[iOverview])
    2586           0 :         return nullptr;
    2587          35 :     return poWDS->m_apoOverviews[iOverview]->GetRasterBand(nBand);
    2588             : }
    2589             : 
    2590             : /*! @endcond */

Generated by: LCOV version 1.14