LCOV - code coverage report
Current view: top level - frmts/vrt - vrtwarped.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 598 714 83.8 %
Date: 2024-05-07 17:03:27 Functions: 37 41 90.2 %

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

Generated by: LCOV version 1.14