LCOV - code coverage report
Current view: top level - frmts/vrt - gdaltileindexdataset.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 1537 1674 91.8 %
Date: 2024-05-06 22:33:47 Functions: 54 54 100.0 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Virtual GDAL Datasets
       4             :  * Purpose:  Tile index based VRT
       5             :  * Author:   Even Rouault <even dot rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2023, Even Rouault <even dot rouault at spatialys.com>
       9             :  *
      10             :  * Permission is hereby granted, free of charge, to any person obtaining a
      11             :  * copy of this software and associated documentation files (the "Software"),
      12             :  * to deal in the Software without restriction, including without limitation
      13             :  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
      14             :  * and/or sell copies of the Software, and to permit persons to whom the
      15             :  * Software is furnished to do so, subject to the following conditions:
      16             :  *
      17             :  * The above copyright notice and this permission notice shall be included
      18             :  * in all copies or substantial portions of the Software.
      19             :  *
      20             :  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
      21             :  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      22             :  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
      23             :  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      24             :  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
      25             :  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
      26             :  * DEALINGS IN THE SOFTWARE.
      27             :  ****************************************************************************/
      28             : 
      29             : /*! @cond Doxygen_Suppress */
      30             : 
      31             : #include <array>
      32             : #include <algorithm>
      33             : #include <limits>
      34             : #include <set>
      35             : #include <tuple>
      36             : #include <utility>
      37             : #include <vector>
      38             : 
      39             : #include "cpl_port.h"
      40             : #include "cpl_mem_cache.h"
      41             : #include "cpl_minixml.h"
      42             : #include "vrtdataset.h"
      43             : #include "vrt_priv.h"
      44             : #include "ogrsf_frmts.h"
      45             : #include "gdal_proxy.h"
      46             : #include "gdal_utils.h"
      47             : 
      48             : #if defined(__SSE2__) || defined(_M_X64)
      49             : #define USE_SSE2_OPTIM
      50             : #include <emmintrin.h>
      51             : // MSVC doesn't define __SSE4_1__, but if -arch:AVX2 is enabled, we do have SSE4.1
      52             : #if defined(__SSE4_1__) || defined(__AVX2__)
      53             : #define USE_SSE41_OPTIM
      54             : #include <smmintrin.h>
      55             : #endif
      56             : #endif
      57             : 
      58             : // Semantincs of indices of a GeoTransform (double[6]) matrix
      59             : constexpr int GT_TOPLEFT_X = 0;
      60             : constexpr int GT_WE_RES = 1;
      61             : constexpr int GT_ROTATION_PARAM1 = 2;
      62             : constexpr int GT_TOPLEFT_Y = 3;
      63             : constexpr int GT_ROTATION_PARAM2 = 4;
      64             : constexpr int GT_NS_RES = 5;
      65             : 
      66             : constexpr const char *GTI_PREFIX = "GTI:";
      67             : 
      68             : constexpr const char *MD_DS_TILE_INDEX_LAYER = "TILE_INDEX_LAYER";
      69             : 
      70             : constexpr const char *MD_RESX = "RESX";
      71             : constexpr const char *MD_RESY = "RESY";
      72             : constexpr const char *MD_BAND_COUNT = "BAND_COUNT";
      73             : constexpr const char *MD_DATA_TYPE = "DATA_TYPE";
      74             : constexpr const char *MD_NODATA = "NODATA";
      75             : constexpr const char *MD_MINX = "MINX";
      76             : constexpr const char *MD_MINY = "MINY";
      77             : constexpr const char *MD_MAXX = "MAXX";
      78             : constexpr const char *MD_MAXY = "MAXY";
      79             : constexpr const char *MD_GEOTRANSFORM = "GEOTRANSFORM";
      80             : constexpr const char *MD_XSIZE = "XSIZE";
      81             : constexpr const char *MD_YSIZE = "YSIZE";
      82             : constexpr const char *MD_COLOR_INTERPRETATION = "COLOR_INTERPRETATION";
      83             : constexpr const char *MD_SRS = "SRS";
      84             : constexpr const char *MD_LOCATION_FIELD = "LOCATION_FIELD";
      85             : constexpr const char *MD_SORT_FIELD = "SORT_FIELD";
      86             : constexpr const char *MD_SORT_FIELD_ASC = "SORT_FIELD_ASC";
      87             : constexpr const char *MD_BLOCK_X_SIZE = "BLOCKXSIZE";
      88             : constexpr const char *MD_BLOCK_Y_SIZE = "BLOCKYSIZE";
      89             : constexpr const char *MD_MASK_BAND = "MASK_BAND";
      90             : constexpr const char *MD_RESAMPLING = "RESAMPLING";
      91             : 
      92             : constexpr const char *const apszTIOptions[] = {MD_RESX,
      93             :                                                MD_RESY,
      94             :                                                MD_BAND_COUNT,
      95             :                                                MD_DATA_TYPE,
      96             :                                                MD_NODATA,
      97             :                                                MD_MINX,
      98             :                                                MD_MINY,
      99             :                                                MD_MAXX,
     100             :                                                MD_MAXY,
     101             :                                                MD_GEOTRANSFORM,
     102             :                                                MD_XSIZE,
     103             :                                                MD_YSIZE,
     104             :                                                MD_COLOR_INTERPRETATION,
     105             :                                                MD_SRS,
     106             :                                                MD_LOCATION_FIELD,
     107             :                                                MD_SORT_FIELD,
     108             :                                                MD_SORT_FIELD_ASC,
     109             :                                                MD_BLOCK_X_SIZE,
     110             :                                                MD_BLOCK_Y_SIZE,
     111             :                                                MD_MASK_BAND,
     112             :                                                MD_RESAMPLING};
     113             : 
     114             : constexpr const char *const MD_BAND_OFFSET = "OFFSET";
     115             : constexpr const char *const MD_BAND_SCALE = "SCALE";
     116             : constexpr const char *const MD_BAND_UNITTYPE = "UNITTYPE";
     117             : constexpr const char *const apszReservedBandItems[] = {
     118             :     MD_BAND_OFFSET, MD_BAND_SCALE, MD_BAND_UNITTYPE};
     119             : 
     120             : constexpr const char *GTI_XML_BANDCOUNT = "BandCount";
     121             : constexpr const char *GTI_XML_DATATYPE = "DataType";
     122             : constexpr const char *GTI_XML_NODATAVALUE = "NoDataValue";
     123             : constexpr const char *GTI_XML_COLORINTERP = "ColorInterp";
     124             : constexpr const char *GTI_XML_LOCATIONFIELD = "LocationField";
     125             : constexpr const char *GTI_XML_SORTFIELD = "SortField";
     126             : constexpr const char *GTI_XML_SORTFIELDASC = "SortFieldAsc";
     127             : constexpr const char *GTI_XML_MASKBAND = "MaskBand";
     128             : constexpr const char *GTI_XML_OVERVIEW_ELEMENT = "Overview";
     129             : constexpr const char *GTI_XML_OVERVIEW_DATASET = "Dataset";
     130             : constexpr const char *GTI_XML_OVERVIEW_LAYER = "Layer";
     131             : constexpr const char *GTI_XML_OVERVIEW_FACTOR = "Factor";
     132             : 
     133             : constexpr const char *GTI_XML_BAND_ELEMENT = "Band";
     134             : constexpr const char *GTI_XML_BAND_NUMBER = "band";
     135             : constexpr const char *GTI_XML_BAND_DATATYPE = "dataType";
     136             : constexpr const char *GTI_XML_BAND_DESCRIPTION = "Description";
     137             : constexpr const char *GTI_XML_BAND_OFFSET = "Offset";
     138             : constexpr const char *GTI_XML_BAND_SCALE = "Scale";
     139             : constexpr const char *GTI_XML_BAND_NODATAVALUE = "NoDataValue";
     140             : constexpr const char *GTI_XML_BAND_UNITTYPE = "UnitType";
     141             : constexpr const char *GTI_XML_BAND_COLORINTERP = "ColorInterp";
     142             : constexpr const char *GTI_XML_CATEGORYNAMES = "CategoryNames";
     143             : constexpr const char *GTI_XML_COLORTABLE = "ColorTable";
     144             : constexpr const char *GTI_XML_RAT = "GDALRasterAttributeTable";
     145             : 
     146             : /************************************************************************/
     147             : /*                           ENDS_WITH_CI()                             */
     148             : /************************************************************************/
     149             : 
     150       47763 : static inline bool ENDS_WITH_CI(const char *a, const char *b)
     151             : {
     152       47763 :     return strlen(a) >= strlen(b) && EQUAL(a + strlen(a) - strlen(b), b);
     153             : }
     154             : 
     155             : /************************************************************************/
     156             : /*                       GDALTileIndexDataset                           */
     157             : /************************************************************************/
     158             : 
     159             : class GDALTileIndexBand;
     160             : 
     161             : class GDALTileIndexDataset final : public GDALPamDataset
     162             : {
     163             :   public:
     164             :     GDALTileIndexDataset();
     165             :     ~GDALTileIndexDataset() override;
     166             : 
     167             :     bool Open(GDALOpenInfo *poOpenInfo);
     168             : 
     169             :     CPLErr FlushCache(bool bAtClosing) override;
     170             : 
     171             :     CPLErr GetGeoTransform(double *padfGeoTransform) override;
     172             :     const OGRSpatialReference *GetSpatialRef() const override;
     173             : 
     174             :     CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,
     175             :                      int nYSize, void *pData, int nBufXSize, int nBufYSize,
     176             :                      GDALDataType eBufType, int nBandCount, int *panBandMap,
     177             :                      GSpacing nPixelSpace, GSpacing nLineSpace,
     178             :                      GSpacing nBandSpace,
     179             :                      GDALRasterIOExtraArg *psExtraArg) override;
     180             : 
     181             :     const char *GetMetadataItem(const char *pszName,
     182             :                                 const char *pszDomain) override;
     183             :     CPLErr SetMetadataItem(const char *pszName, const char *pszValue,
     184             :                            const char *pszDomain) override;
     185             :     CPLErr SetMetadata(char **papszMD, const char *pszDomain) override;
     186             : 
     187             :     void LoadOverviews();
     188             : 
     189             :     std::vector<GTISourceDesc> GetSourcesMoreRecentThan(int64_t mTime);
     190             : 
     191             :   private:
     192             :     friend class GDALTileIndexBand;
     193             : 
     194             :     //! Optional GTI XML
     195             :     CPLXMLTreeCloser m_psXMLTree{nullptr};
     196             : 
     197             :     //! Whether the GTI XML might be modified (by SetMetadata/SetMetadataItem)
     198             :     bool m_bXMLUpdatable = false;
     199             : 
     200             :     //! Whether the GTI XML has been modified (by SetMetadata/SetMetadataItem)
     201             :     bool m_bXMLModified = false;
     202             : 
     203             :     //! Unique string (without the process) for this tile index. Passed to
     204             :     //! GDALProxyPoolDataset to ensure that sources are unique for a given owner
     205             :     const std::string m_osUniqueHandle;
     206             : 
     207             :     //! Vector dataset with the sources
     208             :     std::unique_ptr<GDALDataset> m_poVectorDS{};
     209             : 
     210             :     //! Vector layer with the sources
     211             :     OGRLayer *m_poLayer = nullptr;
     212             : 
     213             :     //! Geotransform matrix of the tile index
     214             :     std::array<double, 6> m_adfGeoTransform{0, 0, 0, 0, 0, 0};
     215             : 
     216             :     //! Index of the "location" (or alternate name given by user) field
     217             :     //! (within m_poLayer->GetLayerDefn()), that contain source dataset names.
     218             :     int m_nLocationFieldIndex = -1;
     219             : 
     220             :     //! SRS of the tile index.
     221             :     OGRSpatialReference m_oSRS{};
     222             : 
     223             :     //! Cache from dataset name to dataset handle.
     224             :     //! Note that the dataset objects are ultimately GDALProxyPoolDataset,
     225             :     //! and that the GDALProxyPoolDataset limits the number of simultaneously
     226             :     //! opened real datasets (controlled by GDAL_MAX_DATASET_POOL_SIZE). Hence 500 is not too big.
     227             :     lru11::Cache<std::string, std::shared_ptr<GDALDataset>> m_oMapSharedSources{
     228             :         500};
     229             : 
     230             :     //! Mask band (e.g. for JPEG compressed + mask band)
     231             :     std::unique_ptr<GDALTileIndexBand> m_poMaskBand{};
     232             : 
     233             :     //! Whether all bands of the tile index have the same data type.
     234             :     bool m_bSameDataType = true;
     235             : 
     236             :     //! Whether all bands of the tile index have the same nodata value.
     237             :     bool m_bSameNoData = true;
     238             : 
     239             :     //! Minimum X of the current pixel request, in georeferenced units.
     240             :     double m_dfLastMinXFilter = std::numeric_limits<double>::quiet_NaN();
     241             : 
     242             :     //! Minimum Y of the current pixel request, in georeferenced units.
     243             :     double m_dfLastMinYFilter = std::numeric_limits<double>::quiet_NaN();
     244             : 
     245             :     //! Maximum X of the current pixel request, in georeferenced units.
     246             :     double m_dfLastMaxXFilter = std::numeric_limits<double>::quiet_NaN();
     247             : 
     248             :     //! Maximum Y of the current pixel request, in georeferenced units.
     249             :     double m_dfLastMaxYFilter = std::numeric_limits<double>::quiet_NaN();
     250             : 
     251             :     //! Index of the field (within m_poLayer->GetLayerDefn()) used to sort, or -1 if none.
     252             :     int m_nSortFieldIndex = -1;
     253             : 
     254             :     //! Whether sorting must be ascending (true) or descending (false).
     255             :     bool m_bSortFieldAsc = true;
     256             : 
     257             :     //! Resampling method by default for warping or when a source has not
     258             :     //! the same resolution as the tile index.
     259             :     std::string m_osResampling = "near";
     260             :     GDALRIOResampleAlg m_eResampling = GRIORA_NearestNeighbour;
     261             : 
     262             :     //! WKT2 representation of the tile index SRS (if needed, typically for on-the-fly warping).
     263             :     std::string m_osWKT{};
     264             : 
     265             :     //! Whether we had to open of the sources at tile index opening.
     266             :     bool m_bScannedOneFeatureAtOpening = false;
     267             : 
     268             :     //! Array of overview descriptors.
     269             :     //! Each descriptor is a tuple (dataset_name, concatenated_open_options, layer_name, overview_factor).
     270             :     std::vector<std::tuple<std::string, CPLStringList, std::string, double>>
     271             :         m_aoOverviewDescriptor{};
     272             : 
     273             :     //! Array of overview datasets.
     274             :     std::vector<std::unique_ptr<GDALDataset>> m_apoOverviews{};
     275             : 
     276             :     //! Cache of buffers used by VRTComplexSource to avoid memory reallocation.
     277             :     VRTSource::WorkingState m_oWorkingState{};
     278             : 
     279             :     //! Structure describing one of the source raster in the tile index.
     280             :     struct SourceDesc
     281             :     {
     282             :         //! Source dataset name.
     283             :         std::string osName{};
     284             : 
     285             :         //! Source dataset handle.
     286             :         std::shared_ptr<GDALDataset> poDS{};
     287             : 
     288             :         //! VRTSimpleSource or VRTComplexSource for the source.
     289             :         std::unique_ptr<VRTSimpleSource> poSource{};
     290             : 
     291             :         //! OGRFeature corresponding to the source in the tile index.
     292             :         std::unique_ptr<OGRFeature> poFeature{};
     293             : 
     294             :         //! Work buffer containing the value of the mask band for the current pixel query.
     295             :         std::vector<GByte> abyMask{};
     296             : 
     297             :         //! Whether the source covers the whole area of interest of the current pixel query.
     298             :         bool bCoversWholeAOI = false;
     299             : 
     300             :         //! Whether the source has a nodata value at least in one of its band.
     301             :         bool bHasNoData = false;
     302             : 
     303             :         //! Whether all bands of the source have the same nodata value.
     304             :         bool bSameNoData = false;
     305             : 
     306             :         //! Nodata value of all bands (when bSameNoData == true).
     307             :         double dfSameNoData = 0;
     308             : 
     309             :         //! Mask band of the source.
     310             :         GDALRasterBand *poMaskBand = nullptr;
     311             :     };
     312             : 
     313             :     //! Array of sources participating to the current pixel query.
     314             :     std::vector<SourceDesc> m_aoSourceDesc{};
     315             : 
     316             :     //! From a source dataset name, return its SourceDesc description structure.
     317             :     bool GetSourceDesc(const std::string &osTileName, SourceDesc &oSourceDesc);
     318             : 
     319             :     //! Collect sources corresponding to the georeferenced window of interest,
     320             :     //! and store them in m_aoSourceDesc[].
     321             :     bool CollectSources(double dfXOff, double dfYOff, double dfXSize,
     322             :                         double dfYSize);
     323             : 
     324             :     //! Sort sources according to m_nSortFieldIndex.
     325             :     void SortSourceDesc();
     326             : 
     327             :     //! Whether the output buffer needs to be nodata initialized, or if
     328             :     //! sources are fully covering it.
     329             :     bool NeedInitBuffer(int nBandCount, const int *panBandMap) const;
     330             : 
     331             :     //! Nodata initialize the output buffer.
     332             :     void InitBuffer(void *pData, int nBufXSize, int nBufYSize,
     333             :                     GDALDataType eBufType, int nBandCount,
     334             :                     const int *panBandMap, GSpacing nPixelSpace,
     335             :                     GSpacing nLineSpace, GSpacing nBandSpace) const;
     336             : 
     337             :     //! Whether m_poVectorDS supports SetMetadata()/SetMetadataItem()
     338             :     bool TileIndexSupportsEditingLayerMetadata() const;
     339             : 
     340             :     CPL_DISALLOW_COPY_ASSIGN(GDALTileIndexDataset)
     341             : };
     342             : 
     343             : /************************************************************************/
     344             : /*                            GDALTileIndexBand                          */
     345             : /************************************************************************/
     346             : 
     347             : class GDALTileIndexBand final : public GDALPamRasterBand
     348             : {
     349             :   public:
     350             :     GDALTileIndexBand(GDALTileIndexDataset *poDSIn, int nBandIn,
     351             :                       GDALDataType eDT, int nBlockXSizeIn, int nBlockYSizeIn);
     352             : 
     353          59 :     double GetNoDataValue(int *pbHasNoData) override
     354             :     {
     355          59 :         if (pbHasNoData)
     356          59 :             *pbHasNoData = m_bNoDataValueSet;
     357          59 :         return m_dfNoDataValue;
     358             :     }
     359             : 
     360          37 :     GDALColorInterp GetColorInterpretation() override
     361             :     {
     362          37 :         return m_eColorInterp;
     363             :     }
     364             : 
     365             :     CPLErr IReadBlock(int nBlockXOff, int nBlockYOff, void *pData) override;
     366             : 
     367             :     CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,
     368             :                      int nYSize, void *pData, int nBufXSize, int nBufYSize,
     369             :                      GDALDataType eBufType, GSpacing nPixelSpace,
     370             :                      GSpacing nLineSpace,
     371             :                      GDALRasterIOExtraArg *psExtraArg) override;
     372             : 
     373          12 :     int GetMaskFlags() override
     374             :     {
     375          12 :         if (m_poDS->m_poMaskBand && m_poDS->m_poMaskBand.get() != this)
     376           4 :             return GMF_PER_DATASET;
     377           8 :         return GDALPamRasterBand::GetMaskFlags();
     378             :     }
     379             : 
     380          18 :     GDALRasterBand *GetMaskBand() override
     381             :     {
     382          18 :         if (m_poDS->m_poMaskBand && m_poDS->m_poMaskBand.get() != this)
     383           7 :             return m_poDS->m_poMaskBand.get();
     384          11 :         return GDALPamRasterBand::GetMaskBand();
     385             :     }
     386             : 
     387           6 :     double GetOffset(int *pbHasValue) override
     388             :     {
     389           6 :         int bHasValue = FALSE;
     390           6 :         double dfVal = GDALPamRasterBand::GetOffset(&bHasValue);
     391           6 :         if (bHasValue)
     392             :         {
     393           0 :             if (pbHasValue)
     394           0 :                 *pbHasValue = true;
     395           0 :             return dfVal;
     396             :         }
     397           6 :         if (pbHasValue)
     398           6 :             *pbHasValue = !std::isnan(m_dfOffset);
     399           6 :         return std::isnan(m_dfOffset) ? 0.0 : m_dfOffset;
     400             :     }
     401             : 
     402           6 :     double GetScale(int *pbHasValue) override
     403             :     {
     404           6 :         int bHasValue = FALSE;
     405           6 :         double dfVal = GDALPamRasterBand::GetScale(&bHasValue);
     406           6 :         if (bHasValue)
     407             :         {
     408           0 :             if (pbHasValue)
     409           0 :                 *pbHasValue = true;
     410           0 :             return dfVal;
     411             :         }
     412           6 :         if (pbHasValue)
     413           6 :             *pbHasValue = !std::isnan(m_dfScale);
     414           6 :         return std::isnan(m_dfScale) ? 1.0 : m_dfScale;
     415             :     }
     416             : 
     417           6 :     const char *GetUnitType() override
     418             :     {
     419           6 :         const char *pszVal = GDALPamRasterBand::GetUnitType();
     420           6 :         if (pszVal && *pszVal)
     421           0 :             return pszVal;
     422           6 :         return m_osUnit.c_str();
     423             :     }
     424             : 
     425           2 :     char **GetCategoryNames() override
     426             :     {
     427           2 :         return m_aosCategoryNames.List();
     428             :     }
     429             : 
     430           6 :     GDALColorTable *GetColorTable() override
     431             :     {
     432           6 :         return m_poColorTable.get();
     433             :     }
     434             : 
     435           2 :     GDALRasterAttributeTable *GetDefaultRAT() override
     436             :     {
     437           2 :         return m_poRAT.get();
     438             :     }
     439             : 
     440             :     int GetOverviewCount() override;
     441             :     GDALRasterBand *GetOverview(int iOvr) override;
     442             : 
     443             :     char **GetMetadataDomainList() override;
     444             :     const char *GetMetadataItem(const char *pszName,
     445             :                                 const char *pszDomain) override;
     446             :     CPLErr SetMetadataItem(const char *pszName, const char *pszValue,
     447             :                            const char *pszDomain) override;
     448             :     CPLErr SetMetadata(char **papszMD, const char *pszDomain) override;
     449             : 
     450             :   private:
     451             :     friend class GDALTileIndexDataset;
     452             : 
     453             :     //! Dataset that owns this band.
     454             :     GDALTileIndexDataset *m_poDS = nullptr;
     455             : 
     456             :     //! Whether a nodata value is set to this band.
     457             :     bool m_bNoDataValueSet = false;
     458             : 
     459             :     //! Nodata value.
     460             :     double m_dfNoDataValue = 0;
     461             : 
     462             :     //! Color interpretation.
     463             :     GDALColorInterp m_eColorInterp = GCI_Undefined;
     464             : 
     465             :     //! Cached value for GetMetadataItem("Pixel_X_Y", "LocationInfo").
     466             :     std::string m_osLastLocationInfo{};
     467             : 
     468             :     //! Scale value (returned by GetScale())
     469             :     double m_dfScale = std::numeric_limits<double>::quiet_NaN();
     470             : 
     471             :     //! Offset value (returned by GetOffset())
     472             :     double m_dfOffset = std::numeric_limits<double>::quiet_NaN();
     473             : 
     474             :     //! Unit type (returned by GetUnitType()).
     475             :     std::string m_osUnit{};
     476             : 
     477             :     //! Category names (returned by GetCategoryNames()).
     478             :     CPLStringList m_aosCategoryNames{};
     479             : 
     480             :     //! Color table (returned by GetColorTable()).
     481             :     std::unique_ptr<GDALColorTable> m_poColorTable{};
     482             : 
     483             :     //! Raster attribute table (returned by GetDefaultRAT()).
     484             :     std::unique_ptr<GDALRasterAttributeTable> m_poRAT{};
     485             : 
     486             :     CPL_DISALLOW_COPY_ASSIGN(GDALTileIndexBand)
     487             : };
     488             : 
     489             : /************************************************************************/
     490             : /*                        IsSameNaNAware()                              */
     491             : /************************************************************************/
     492             : 
     493         151 : static inline bool IsSameNaNAware(double a, double b)
     494             : {
     495         151 :     return a == b || (std::isnan(a) && std::isnan(b));
     496             : }
     497             : 
     498             : /************************************************************************/
     499             : /*                         GDALTileIndexDataset()                        */
     500             : /************************************************************************/
     501             : 
     502         185 : GDALTileIndexDataset::GDALTileIndexDataset()
     503         185 :     : m_osUniqueHandle(CPLSPrintf("%p", this))
     504             : {
     505         185 : }
     506             : 
     507             : /************************************************************************/
     508             : /*                        GetAbsoluteFileName()                         */
     509             : /************************************************************************/
     510             : 
     511         278 : static std::string GetAbsoluteFileName(const char *pszTileName,
     512             :                                        const char *pszVRTName)
     513             : {
     514         278 :     if (CPLIsFilenameRelative(pszTileName) &&
     515         282 :         !STARTS_WITH(pszTileName, "<VRTDataset") &&
     516           4 :         !STARTS_WITH(pszVRTName, "<GDALTileIndexDataset"))
     517             :     {
     518           4 :         const auto oSubDSInfo(GDALGetSubdatasetInfo(pszTileName));
     519           4 :         if (oSubDSInfo && !oSubDSInfo->GetPathComponent().empty())
     520             :         {
     521           4 :             const std::string osPath(oSubDSInfo->GetPathComponent());
     522             :             const std::string osRet =
     523           2 :                 CPLIsFilenameRelative(osPath.c_str())
     524             :                     ? oSubDSInfo->ModifyPathComponent(
     525             :                           CPLProjectRelativeFilename(CPLGetPath(pszVRTName),
     526             :                                                      osPath.c_str()))
     527           6 :                     : std::string(pszTileName);
     528           2 :             GDALDestroySubdatasetInfo(oSubDSInfo);
     529           2 :             return osRet;
     530             :         }
     531             : 
     532             :         const std::string osRelativeMadeAbsolute =
     533           2 :             CPLProjectRelativeFilename(CPLGetPath(pszVRTName), pszTileName);
     534             :         VSIStatBufL sStat;
     535           2 :         if (VSIStatL(osRelativeMadeAbsolute.c_str(), &sStat) == 0)
     536           2 :             return osRelativeMadeAbsolute;
     537             :     }
     538         274 :     return pszTileName;
     539             : }
     540             : 
     541             : /************************************************************************/
     542             : /*                    GTIDoPaletteExpansionIfNeeded()                   */
     543             : /************************************************************************/
     544             : 
     545             : //! Do palette -> RGB(A) expansion
     546             : static bool
     547         217 : GTIDoPaletteExpansionIfNeeded(std::shared_ptr<GDALDataset> &poTileDS,
     548             :                               int nBandCount)
     549             : {
     550         408 :     if (poTileDS->GetRasterCount() == 1 &&
     551         410 :         (nBandCount == 3 || nBandCount == 4) &&
     552           4 :         poTileDS->GetRasterBand(1)->GetColorTable() != nullptr)
     553             :     {
     554             : 
     555           4 :         CPLStringList aosOptions;
     556           4 :         aosOptions.AddString("-of");
     557           4 :         aosOptions.AddString("VRT");
     558             : 
     559           4 :         aosOptions.AddString("-expand");
     560           4 :         aosOptions.AddString(nBandCount == 3 ? "rgb" : "rgba");
     561             : 
     562             :         GDALTranslateOptions *psOptions =
     563           4 :             GDALTranslateOptionsNew(aosOptions.List(), nullptr);
     564           4 :         int bUsageError = false;
     565             :         auto poRGBDS = std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(
     566             :             GDALTranslate("", GDALDataset::ToHandle(poTileDS.get()), psOptions,
     567           4 :                           &bUsageError)));
     568           4 :         GDALTranslateOptionsFree(psOptions);
     569           4 :         if (!poRGBDS)
     570             :         {
     571           0 :             return false;
     572             :         }
     573             : 
     574           4 :         poTileDS.reset(poRGBDS.release());
     575             :     }
     576         217 :     return true;
     577             : }
     578             : 
     579             : /************************************************************************/
     580             : /*                                Open()                                */
     581             : /************************************************************************/
     582             : 
     583         185 : bool GDALTileIndexDataset::Open(GDALOpenInfo *poOpenInfo)
     584             : {
     585         185 :     eAccess = poOpenInfo->eAccess;
     586             : 
     587         185 :     CPLXMLNode *psRoot = nullptr;
     588         185 :     const char *pszIndexDataset = poOpenInfo->pszFilename;
     589             : 
     590         185 :     if (STARTS_WITH(poOpenInfo->pszFilename, GTI_PREFIX))
     591             :     {
     592           3 :         pszIndexDataset = poOpenInfo->pszFilename + strlen(GTI_PREFIX);
     593             :     }
     594         182 :     else if (STARTS_WITH(poOpenInfo->pszFilename, "<GDALTileIndexDataset"))
     595             :     {
     596             :         // CPLParseXMLString() emits an error in case of failure
     597          22 :         m_psXMLTree.reset(CPLParseXMLString(poOpenInfo->pszFilename));
     598          22 :         if (m_psXMLTree == nullptr)
     599           1 :             return false;
     600             :     }
     601         160 :     else if (poOpenInfo->nHeaderBytes > 0 &&
     602         160 :              strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
     603             :                     "<GDALTileIndexDataset"))
     604             :     {
     605             :         // CPLParseXMLFile() emits an error in case of failure
     606           6 :         m_psXMLTree.reset(CPLParseXMLFile(poOpenInfo->pszFilename));
     607           6 :         if (m_psXMLTree == nullptr)
     608           1 :             return false;
     609           5 :         m_bXMLUpdatable = (poOpenInfo->eAccess == GA_Update);
     610             :     }
     611             : 
     612         183 :     if (m_psXMLTree)
     613             :     {
     614          26 :         psRoot = CPLGetXMLNode(m_psXMLTree.get(), "=GDALTileIndexDataset");
     615          26 :         if (psRoot == nullptr)
     616             :         {
     617           1 :             CPLError(CE_Failure, CPLE_AppDefined,
     618             :                      "Missing GDALTileIndexDataset root element.");
     619           1 :             return false;
     620             :         }
     621             : 
     622          25 :         pszIndexDataset = CPLGetXMLValue(psRoot, "IndexDataset", nullptr);
     623          25 :         if (!pszIndexDataset)
     624             :         {
     625           1 :             CPLError(CE_Failure, CPLE_AppDefined,
     626             :                      "Missing IndexDataset element.");
     627           1 :             return false;
     628             :         }
     629             :     }
     630             : 
     631         181 :     if (ENDS_WITH_CI(pszIndexDataset, ".gti.gpkg") &&
     632         339 :         poOpenInfo->nHeaderBytes >= 100 &&
     633         158 :         STARTS_WITH(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
     634             :                     "SQLite format 3"))
     635             :     {
     636         154 :         const char *const apszAllowedDrivers[] = {"GPKG", nullptr};
     637         154 :         m_poVectorDS.reset(GDALDataset::Open(
     638         308 :             std::string("GPKG:\"").append(pszIndexDataset).append("\"").c_str(),
     639         154 :             GDAL_OF_VECTOR | GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR |
     640         154 :                 ((poOpenInfo->nOpenFlags & GDAL_OF_UPDATE) ? GDAL_OF_UPDATE
     641             :                                                            : GDAL_OF_READONLY),
     642             :             apszAllowedDrivers));
     643         154 :         if (!m_poVectorDS)
     644             :         {
     645           1 :             return false;
     646             :         }
     647         156 :         if (m_poVectorDS->GetLayerCount() == 0 &&
     648           2 :             (m_poVectorDS->GetRasterCount() != 0 ||
     649           1 :              m_poVectorDS->GetMetadata("SUBDATASETS") != nullptr))
     650             :         {
     651           1 :             return false;
     652             :         }
     653             :     }
     654             :     else
     655             :     {
     656          27 :         m_poVectorDS.reset(GDALDataset::Open(
     657          27 :             pszIndexDataset, GDAL_OF_VECTOR | GDAL_OF_VERBOSE_ERROR |
     658          27 :                                  ((poOpenInfo->nOpenFlags & GDAL_OF_UPDATE)
     659          27 :                                       ? GDAL_OF_UPDATE
     660             :                                       : GDAL_OF_READONLY)));
     661          27 :         if (!m_poVectorDS)
     662             :         {
     663           1 :             return false;
     664             :         }
     665             :     }
     666             : 
     667         179 :     if (m_poVectorDS->GetLayerCount() == 0)
     668             :     {
     669           1 :         CPLError(CE_Failure, CPLE_AppDefined, "%s has no vector layer",
     670             :                  poOpenInfo->pszFilename);
     671           1 :         return false;
     672             :     }
     673             : 
     674         178 :     double dfOvrFactor = 1.0;
     675         178 :     if (const char *pszFactor =
     676         178 :             CSLFetchNameValue(poOpenInfo->papszOpenOptions, "FACTOR"))
     677             :     {
     678           5 :         dfOvrFactor = CPLAtof(pszFactor);
     679           5 :         if (!(dfOvrFactor > 1.0))
     680             :         {
     681           1 :             CPLError(CE_Failure, CPLE_AppDefined, "Wrong overview factor");
     682           1 :             return false;
     683             :         }
     684             :     }
     685             : 
     686             :     const char *pszLayerName;
     687             : 
     688         177 :     if ((pszLayerName = CSLFetchNameValue(poOpenInfo->papszOpenOptions,
     689         177 :                                           "LAYER")) != nullptr)
     690             :     {
     691           5 :         m_poLayer = m_poVectorDS->GetLayerByName(pszLayerName);
     692           5 :         if (!m_poLayer)
     693             :         {
     694           2 :             CPLError(CE_Failure, CPLE_AppDefined, "Layer %s does not exist",
     695             :                      pszLayerName);
     696           2 :             return false;
     697             :         }
     698             :     }
     699         172 :     else if (psRoot && (pszLayerName = CPLGetXMLValue(psRoot, "IndexLayer",
     700             :                                                       nullptr)) != nullptr)
     701             :     {
     702           2 :         m_poLayer = m_poVectorDS->GetLayerByName(pszLayerName);
     703           2 :         if (!m_poLayer)
     704             :         {
     705           1 :             CPLError(CE_Failure, CPLE_AppDefined, "Layer %s does not exist",
     706             :                      pszLayerName);
     707           1 :             return false;
     708             :         }
     709             :     }
     710         322 :     else if (!psRoot && (pszLayerName = m_poVectorDS->GetMetadataItem(
     711         152 :                              MD_DS_TILE_INDEX_LAYER)) != nullptr)
     712             :     {
     713           2 :         m_poLayer = m_poVectorDS->GetLayerByName(pszLayerName);
     714           2 :         if (!m_poLayer)
     715             :         {
     716           1 :             CPLError(CE_Failure, CPLE_AppDefined, "Layer %s does not exist",
     717             :                      pszLayerName);
     718           1 :             return false;
     719             :         }
     720             :     }
     721         168 :     else if (m_poVectorDS->GetLayerCount() == 1)
     722             :     {
     723         166 :         m_poLayer = m_poVectorDS->GetLayer(0);
     724         166 :         if (!m_poLayer)
     725             :         {
     726           0 :             CPLError(CE_Failure, CPLE_AppDefined, "Cannot open layer 0");
     727           0 :             return false;
     728             :         }
     729             :     }
     730             :     else
     731             :     {
     732           2 :         if (STARTS_WITH(poOpenInfo->pszFilename, GTI_PREFIX))
     733             :         {
     734           1 :             CPLError(CE_Failure, CPLE_AppDefined,
     735             :                      "%s has more than one layer. LAYER open option "
     736             :                      "must be defined to specify which one to "
     737             :                      "use as the tile index",
     738             :                      pszIndexDataset);
     739             :         }
     740           1 :         else if (psRoot)
     741             :         {
     742           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     743             :                      "%s has more than one layer. IndexLayer element must be "
     744             :                      "defined to specify which one to "
     745             :                      "use as the tile index",
     746             :                      pszIndexDataset);
     747             :         }
     748             :         else
     749             :         {
     750           1 :             CPLError(CE_Failure, CPLE_AppDefined,
     751             :                      "%s has more than one layer. %s "
     752             :                      "metadata item must be defined to specify which one to "
     753             :                      "use as the tile index",
     754             :                      pszIndexDataset, MD_DS_TILE_INDEX_LAYER);
     755             :         }
     756           2 :         return false;
     757             :     }
     758             : 
     759             :     // Try to get the metadata from an embedded xml:GTI domain
     760         171 :     if (!m_psXMLTree)
     761             :     {
     762         150 :         char **papszMD = m_poLayer->GetMetadata("xml:GTI");
     763         150 :         if (papszMD && papszMD[0])
     764             :         {
     765           1 :             m_psXMLTree.reset(CPLParseXMLString(papszMD[0]));
     766           1 :             if (m_psXMLTree == nullptr)
     767           0 :                 return false;
     768             : 
     769           1 :             psRoot = CPLGetXMLNode(m_psXMLTree.get(), "=GDALTileIndexDataset");
     770           1 :             if (psRoot == nullptr)
     771             :             {
     772           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
     773             :                          "Missing GDALTileIndexDataset root element.");
     774           0 :                 return false;
     775             :             }
     776             :         }
     777             :     }
     778             : 
     779             :     // Get the value of an option.
     780             :     // The order of lookup is the following one (first to last):
     781             :     // - open options
     782             :     // - XML file
     783             :     // - Layer metadata items.
     784       13349 :     const auto GetOption = [poOpenInfo, psRoot, this](const char *pszItem)
     785             :     {
     786             :         const char *pszVal =
     787        4138 :             CSLFetchNameValue(poOpenInfo->papszOpenOptions, pszItem);
     788        4138 :         if (pszVal)
     789           2 :             return pszVal;
     790             : 
     791        4136 :         if (psRoot)
     792             :         {
     793         492 :             pszVal = CPLGetXMLValue(psRoot, pszItem, nullptr);
     794         492 :             if (pszVal)
     795          19 :                 return pszVal;
     796             : 
     797         473 :             if (EQUAL(pszItem, MD_BAND_COUNT))
     798          19 :                 pszItem = GTI_XML_BANDCOUNT;
     799         454 :             else if (EQUAL(pszItem, MD_DATA_TYPE))
     800          22 :                 pszItem = GTI_XML_DATATYPE;
     801         432 :             else if (EQUAL(pszItem, MD_NODATA))
     802          17 :                 pszItem = GTI_XML_NODATAVALUE;
     803         415 :             else if (EQUAL(pszItem, MD_COLOR_INTERPRETATION))
     804          22 :                 pszItem = GTI_XML_COLORINTERP;
     805         393 :             else if (EQUAL(pszItem, MD_LOCATION_FIELD))
     806          22 :                 pszItem = GTI_XML_LOCATIONFIELD;
     807         371 :             else if (EQUAL(pszItem, MD_SORT_FIELD))
     808          22 :                 pszItem = GTI_XML_SORTFIELD;
     809         349 :             else if (EQUAL(pszItem, MD_SORT_FIELD_ASC))
     810           2 :                 pszItem = GTI_XML_SORTFIELDASC;
     811         347 :             else if (EQUAL(pszItem, MD_MASK_BAND))
     812          17 :                 pszItem = GTI_XML_MASKBAND;
     813         473 :             pszVal = CPLGetXMLValue(psRoot, pszItem, nullptr);
     814         473 :             if (pszVal)
     815           7 :                 return pszVal;
     816             :         }
     817             : 
     818        4110 :         return m_poLayer->GetMetadataItem(pszItem);
     819         171 :     };
     820             : 
     821         171 :     const char *pszFilter = GetOption("Filter");
     822         171 :     if (pszFilter)
     823             :     {
     824           1 :         if (m_poLayer->SetAttributeFilter(pszFilter) != OGRERR_NONE)
     825           0 :             return false;
     826             :     }
     827             : 
     828         171 :     const char *pszLocationFieldName = GetOption(MD_LOCATION_FIELD);
     829         171 :     if (!pszLocationFieldName)
     830             :     {
     831         168 :         constexpr const char *DEFAULT_LOCATION_FIELD_NAME = "location";
     832         168 :         pszLocationFieldName = DEFAULT_LOCATION_FIELD_NAME;
     833             :     }
     834         171 :     auto poLayerDefn = m_poLayer->GetLayerDefn();
     835         171 :     m_nLocationFieldIndex = poLayerDefn->GetFieldIndex(pszLocationFieldName);
     836         171 :     if (m_nLocationFieldIndex < 0)
     837             :     {
     838           1 :         CPLError(CE_Failure, CPLE_AppDefined, "Cannot find field %s",
     839             :                  pszLocationFieldName);
     840           1 :         return false;
     841             :     }
     842         170 :     if (poLayerDefn->GetFieldDefn(m_nLocationFieldIndex)->GetType() !=
     843             :         OFTString)
     844             :     {
     845           1 :         CPLError(CE_Failure, CPLE_AppDefined, "Field %s is not of type string",
     846             :                  pszLocationFieldName);
     847           1 :         return false;
     848             :     }
     849             : 
     850         169 :     const char *pszSortFieldName = GetOption(MD_SORT_FIELD);
     851         169 :     if (pszSortFieldName)
     852             :     {
     853          40 :         m_nSortFieldIndex = poLayerDefn->GetFieldIndex(pszSortFieldName);
     854          40 :         if (m_nSortFieldIndex < 0)
     855             :         {
     856           1 :             CPLError(CE_Failure, CPLE_AppDefined, "Cannot find field %s",
     857             :                      pszSortFieldName);
     858           1 :             return false;
     859             :         }
     860             : 
     861             :         const auto eFieldType =
     862          39 :             poLayerDefn->GetFieldDefn(m_nSortFieldIndex)->GetType();
     863          39 :         if (eFieldType != OFTString && eFieldType != OFTInteger &&
     864          21 :             eFieldType != OFTInteger64 && eFieldType != OFTReal &&
     865           7 :             eFieldType != OFTDate && eFieldType != OFTDateTime)
     866             :         {
     867           1 :             CPLError(CE_Failure, CPLE_AppDefined,
     868             :                      "Unsupported type for field %s", pszSortFieldName);
     869           1 :             return false;
     870             :         }
     871             : 
     872          38 :         const char *pszSortFieldAsc = GetOption(MD_SORT_FIELD_ASC);
     873          38 :         if (pszSortFieldAsc)
     874             :         {
     875           3 :             m_bSortFieldAsc = CPLTestBool(pszSortFieldAsc);
     876             :         }
     877             :     }
     878             : 
     879         167 :     const char *pszResX = GetOption(MD_RESX);
     880         167 :     const char *pszResY = GetOption(MD_RESY);
     881         167 :     if (pszResX && !pszResY)
     882             :     {
     883           1 :         CPLError(CE_Failure, CPLE_AppDefined,
     884             :                  "%s metadata item defined, but not %s", MD_RESX, MD_RESY);
     885           1 :         return false;
     886             :     }
     887         166 :     if (!pszResX && pszResY)
     888             :     {
     889           1 :         CPLError(CE_Failure, CPLE_AppDefined,
     890             :                  "%s metadata item defined, but not %s", MD_RESY, MD_RESX);
     891           1 :         return false;
     892             :     }
     893             : 
     894         165 :     const char *pszResampling = GetOption(MD_RESAMPLING);
     895         165 :     if (pszResampling)
     896             :     {
     897           6 :         const auto nErrorCountBefore = CPLGetErrorCounter();
     898           6 :         m_eResampling = GDALRasterIOGetResampleAlg(pszResampling);
     899           6 :         if (nErrorCountBefore != CPLGetErrorCounter())
     900             :         {
     901           0 :             return false;
     902             :         }
     903           6 :         m_osResampling = pszResampling;
     904             :     }
     905             : 
     906         165 :     const char *pszMinX = GetOption(MD_MINX);
     907         165 :     const char *pszMinY = GetOption(MD_MINY);
     908         165 :     const char *pszMaxX = GetOption(MD_MAXX);
     909         165 :     const char *pszMaxY = GetOption(MD_MAXY);
     910         165 :     const int nCountMinMaxXY = (pszMinX ? 1 : 0) + (pszMinY ? 1 : 0) +
     911         165 :                                (pszMaxX ? 1 : 0) + (pszMaxY ? 1 : 0);
     912         165 :     if (nCountMinMaxXY != 0 && nCountMinMaxXY != 4)
     913             :     {
     914           4 :         CPLError(CE_Failure, CPLE_AppDefined,
     915             :                  "None or all of %s, %s, %s and %s must be specified", MD_MINX,
     916             :                  MD_MINY, MD_MAXX, MD_MAXY);
     917           4 :         return false;
     918             :     }
     919             : 
     920         161 :     const char *pszXSize = GetOption(MD_XSIZE);
     921         161 :     const char *pszYSize = GetOption(MD_YSIZE);
     922         161 :     const char *pszGeoTransform = GetOption(MD_GEOTRANSFORM);
     923         161 :     const int nCountXSizeYSizeGT =
     924         161 :         (pszXSize ? 1 : 0) + (pszYSize ? 1 : 0) + (pszGeoTransform ? 1 : 0);
     925         161 :     if (nCountXSizeYSizeGT != 0 && nCountXSizeYSizeGT != 3)
     926             :     {
     927           3 :         CPLError(CE_Failure, CPLE_AppDefined,
     928             :                  "None or all of %s, %s, %s must be specified", MD_XSIZE,
     929             :                  MD_YSIZE, MD_GEOTRANSFORM);
     930           3 :         return false;
     931             :     }
     932             : 
     933         158 :     const char *pszDataType = GetOption(MD_DATA_TYPE);
     934         158 :     const char *pszColorInterp = GetOption(MD_COLOR_INTERPRETATION);
     935         158 :     int nBandCount = 0;
     936         316 :     std::vector<GDALDataType> aeDataTypes;
     937         316 :     std::vector<std::pair<bool, double>> aNoData;
     938         316 :     std::vector<GDALColorInterp> aeColorInterp;
     939             : 
     940         158 :     const char *pszSRS = GetOption(MD_SRS);
     941         158 :     m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     942         158 :     if (pszSRS)
     943             :     {
     944           2 :         if (m_oSRS.SetFromUserInput(
     945             :                 pszSRS,
     946           2 :                 OGRSpatialReference::SET_FROM_USER_INPUT_LIMITATIONS_get()) !=
     947             :             OGRERR_NONE)
     948             :         {
     949           1 :             CPLError(CE_Failure, CPLE_AppDefined, "Invalid %s", MD_SRS);
     950           1 :             return false;
     951             :         }
     952             :     }
     953             :     else
     954             :     {
     955         156 :         const auto poSRS = m_poLayer->GetSpatialRef();
     956             :         // Ignore GPKG "Undefined geographic SRS" and "Undefined Cartesian SRS"
     957         156 :         if (poSRS && !STARTS_WITH(poSRS->GetName(), "Undefined "))
     958         100 :             m_oSRS = *poSRS;
     959             :     }
     960             : 
     961         314 :     std::vector<const CPLXMLNode *> apoXMLNodeBands;
     962         157 :     if (psRoot)
     963             :     {
     964          22 :         int nExpectedBandNumber = 1;
     965          86 :         for (const CPLXMLNode *psIter = psRoot->psChild; psIter;
     966          64 :              psIter = psIter->psNext)
     967             :         {
     968          67 :             if (psIter->eType == CXT_Element &&
     969          67 :                 strcmp(psIter->pszValue, GTI_XML_BAND_ELEMENT) == 0)
     970             :             {
     971             :                 const char *pszBand =
     972           8 :                     CPLGetXMLValue(psIter, GTI_XML_BAND_NUMBER, nullptr);
     973           8 :                 if (!pszBand)
     974             :                 {
     975           1 :                     CPLError(CE_Failure, CPLE_AppDefined,
     976             :                              "%s attribute missing on %s element",
     977             :                              GTI_XML_BAND_NUMBER, GTI_XML_BAND_ELEMENT);
     978           3 :                     return false;
     979             :                 }
     980           7 :                 const int nBand = atoi(pszBand);
     981           7 :                 if (nBand <= 0)
     982             :                 {
     983           1 :                     CPLError(CE_Failure, CPLE_AppDefined,
     984             :                              "Invalid band number");
     985           1 :                     return false;
     986             :                 }
     987           6 :                 if (nBand != nExpectedBandNumber)
     988             :                 {
     989           1 :                     CPLError(CE_Failure, CPLE_AppDefined,
     990             :                              "Invalid band number: found %d, expected %d",
     991             :                              nBand, nExpectedBandNumber);
     992           1 :                     return false;
     993             :                 }
     994           5 :                 apoXMLNodeBands.push_back(psIter);
     995           5 :                 ++nExpectedBandNumber;
     996             :             }
     997             :         }
     998             :     }
     999             : 
    1000         154 :     const char *pszBandCount = GetOption(MD_BAND_COUNT);
    1001         154 :     if (pszBandCount)
    1002          22 :         nBandCount = atoi(pszBandCount);
    1003             : 
    1004         154 :     if (!apoXMLNodeBands.empty())
    1005             :     {
    1006           5 :         if (!pszBandCount)
    1007           4 :             nBandCount = static_cast<int>(apoXMLNodeBands.size());
    1008           1 :         else if (nBandCount != static_cast<int>(apoXMLNodeBands.size()))
    1009             :         {
    1010           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    1011             :                      "Inconsistent %s with actual number of %s elements",
    1012             :                      GTI_XML_BANDCOUNT, GTI_XML_BAND_ELEMENT);
    1013           1 :             return false;
    1014             :         }
    1015             :     }
    1016             : 
    1017         153 :     bool bHasMaskBand = false;
    1018         178 :     if ((!pszBandCount && apoXMLNodeBands.empty()) ||
    1019          25 :         (!(pszResX && pszResY) && nCountXSizeYSizeGT == 0))
    1020             :     {
    1021         142 :         CPLDebug("VRT", "Inspecting one feature due to missing metadata items");
    1022         142 :         m_bScannedOneFeatureAtOpening = true;
    1023             : 
    1024             :         auto poFeature =
    1025         142 :             std::unique_ptr<OGRFeature>(m_poLayer->GetNextFeature());
    1026         282 :         if (!poFeature ||
    1027         140 :             !poFeature->IsFieldSetAndNotNull(m_nLocationFieldIndex))
    1028             :         {
    1029           2 :             CPLError(
    1030             :                 CE_Failure, CPLE_AppDefined,
    1031             :                 "BAND_COUNT(+DATA_TYPE+COLOR_INTERPRETATION)+ (RESX+RESY or "
    1032             :                 "XSIZE+YSIZE+GEOTRANSFORM) metadata items "
    1033             :                 "missing");
    1034           2 :             return false;
    1035             :         }
    1036             : 
    1037             :         const char *pszTileName =
    1038         140 :             poFeature->GetFieldAsString(m_nLocationFieldIndex);
    1039             :         const std::string osTileName(
    1040         140 :             GetAbsoluteFileName(pszTileName, poOpenInfo->pszFilename));
    1041         140 :         pszTileName = osTileName.c_str();
    1042             : 
    1043             :         auto poTileDS = std::shared_ptr<GDALDataset>(
    1044             :             GDALDataset::Open(pszTileName,
    1045             :                               GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR),
    1046         140 :             GDALDatasetUniquePtrReleaser());
    1047         140 :         if (!poTileDS)
    1048             :         {
    1049           1 :             return false;
    1050             :         }
    1051             : 
    1052             :         // do palette -> RGB(A) expansion if needed
    1053         139 :         if (!GTIDoPaletteExpansionIfNeeded(poTileDS, nBandCount))
    1054           0 :             return false;
    1055             : 
    1056         139 :         const int nTileBandCount = poTileDS->GetRasterCount();
    1057         301 :         for (int i = 0; i < nTileBandCount; ++i)
    1058             :         {
    1059         162 :             auto poTileBand = poTileDS->GetRasterBand(i + 1);
    1060         162 :             aeDataTypes.push_back(poTileBand->GetRasterDataType());
    1061         162 :             int bHasNoData = FALSE;
    1062         162 :             const double dfNoData = poTileBand->GetNoDataValue(&bHasNoData);
    1063         162 :             aNoData.emplace_back(CPL_TO_BOOL(bHasNoData), dfNoData);
    1064         162 :             aeColorInterp.push_back(poTileBand->GetColorInterpretation());
    1065             : 
    1066         162 :             if (poTileBand->GetMaskFlags() == GMF_PER_DATASET)
    1067           1 :                 bHasMaskBand = true;
    1068             :         }
    1069         139 :         if (!pszBandCount && nBandCount == 0)
    1070         125 :             nBandCount = nTileBandCount;
    1071             : 
    1072         139 :         auto poTileSRS = poTileDS->GetSpatialRef();
    1073         139 :         if (!m_oSRS.IsEmpty() && poTileSRS && !m_oSRS.IsSame(poTileSRS))
    1074             :         {
    1075           3 :             CPLStringList aosOptions;
    1076           3 :             aosOptions.AddString("-of");
    1077           3 :             aosOptions.AddString("VRT");
    1078             : 
    1079           3 :             char *pszWKT = nullptr;
    1080           3 :             const char *const apszWKTOptions[] = {"FORMAT=WKT2_2019", nullptr};
    1081           3 :             m_oSRS.exportToWkt(&pszWKT, apszWKTOptions);
    1082           3 :             if (pszWKT)
    1083           3 :                 m_osWKT = pszWKT;
    1084           3 :             CPLFree(pszWKT);
    1085             : 
    1086           3 :             if (m_osWKT.empty())
    1087             :             {
    1088           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
    1089             :                          "Cannot export VRT SRS to WKT2");
    1090           0 :                 return false;
    1091             :             }
    1092             : 
    1093           3 :             aosOptions.AddString("-t_srs");
    1094           3 :             aosOptions.AddString(m_osWKT.c_str());
    1095             : 
    1096             :             GDALWarpAppOptions *psWarpOptions =
    1097           3 :                 GDALWarpAppOptionsNew(aosOptions.List(), nullptr);
    1098           3 :             GDALDatasetH ahSrcDS[] = {GDALDataset::ToHandle(poTileDS.get())};
    1099           3 :             int bUsageError = false;
    1100             :             auto poWarpDS =
    1101             :                 std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(GDALWarp(
    1102           3 :                     "", nullptr, 1, ahSrcDS, psWarpOptions, &bUsageError)));
    1103           3 :             GDALWarpAppOptionsFree(psWarpOptions);
    1104           3 :             if (!poWarpDS)
    1105             :             {
    1106           0 :                 return false;
    1107             :             }
    1108             : 
    1109           3 :             poTileDS.reset(poWarpDS.release());
    1110           3 :             poTileSRS = poTileDS->GetSpatialRef();
    1111           3 :             CPL_IGNORE_RET_VAL(poTileSRS);
    1112             :         }
    1113             : 
    1114             :         double adfGeoTransformTile[6];
    1115         139 :         if (poTileDS->GetGeoTransform(adfGeoTransformTile) != CE_None)
    1116             :         {
    1117           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    1118             :                      "Cannot find geotransform on %s", pszTileName);
    1119           1 :             return false;
    1120             :         }
    1121         138 :         if (!(adfGeoTransformTile[GT_ROTATION_PARAM1] == 0))
    1122             :         {
    1123           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    1124             :                      "3rd value of GeoTransform of %s must be 0", pszTileName);
    1125           1 :             return false;
    1126             :         }
    1127         137 :         if (!(adfGeoTransformTile[GT_ROTATION_PARAM2] == 0))
    1128             :         {
    1129           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    1130             :                      "5th value of GeoTransform of %s must be 0", pszTileName);
    1131           1 :             return false;
    1132             :         }
    1133         136 :         if (!(adfGeoTransformTile[GT_NS_RES] < 0))
    1134             :         {
    1135           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    1136             :                      "6th value of GeoTransform of %s must be < 0",
    1137             :                      pszTileName);
    1138           1 :             return false;
    1139             :         }
    1140             : 
    1141         135 :         const double dfResX = adfGeoTransformTile[GT_WE_RES];
    1142         135 :         const double dfResY = -adfGeoTransformTile[GT_NS_RES];
    1143             : 
    1144         135 :         OGREnvelope sEnvelope;
    1145         135 :         if (m_poLayer->GetExtent(&sEnvelope, /* bForce = */ false) ==
    1146             :             OGRERR_FAILURE)
    1147             :         {
    1148           1 :             if (m_poLayer->GetExtent(&sEnvelope, /* bForce = */ true) ==
    1149             :                 OGRERR_FAILURE)
    1150             :             {
    1151           1 :                 CPLError(CE_Failure, CPLE_AppDefined,
    1152             :                          "Cannot get layer extent");
    1153           1 :                 return false;
    1154             :             }
    1155           0 :             CPLError(CE_Warning, CPLE_AppDefined,
    1156             :                      "Could get layer extent, but using a slower method");
    1157             :         }
    1158             : 
    1159         134 :         const double dfXSize = (sEnvelope.MaxX - sEnvelope.MinX) / dfResX;
    1160         134 :         if (!(dfXSize >= 0 && dfXSize < INT_MAX))
    1161             :         {
    1162           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    1163             :                      "Too small %s, or wrong layer extent", MD_RESX);
    1164           1 :             return false;
    1165             :         }
    1166             : 
    1167         133 :         const double dfYSize = (sEnvelope.MaxY - sEnvelope.MinY) / dfResY;
    1168         133 :         if (!(dfYSize >= 0 && dfYSize < INT_MAX))
    1169             :         {
    1170           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    1171             :                      "Too small %s, or wrong layer extent", MD_RESY);
    1172           1 :             return false;
    1173             :         }
    1174             : 
    1175         132 :         m_adfGeoTransform[GT_TOPLEFT_X] = sEnvelope.MinX;
    1176         132 :         m_adfGeoTransform[GT_WE_RES] = dfResX;
    1177         132 :         m_adfGeoTransform[GT_ROTATION_PARAM1] = 0;
    1178         132 :         m_adfGeoTransform[GT_TOPLEFT_Y] = sEnvelope.MaxY;
    1179         132 :         m_adfGeoTransform[GT_ROTATION_PARAM2] = 0;
    1180         132 :         m_adfGeoTransform[GT_NS_RES] = -dfResY;
    1181         132 :         nRasterXSize = static_cast<int>(std::ceil(dfXSize));
    1182         132 :         nRasterYSize = static_cast<int>(std::ceil(dfYSize));
    1183             :     }
    1184             : 
    1185         143 :     if (pszXSize && pszYSize && pszGeoTransform)
    1186             :     {
    1187          12 :         const int nXSize = atoi(pszXSize);
    1188          12 :         if (nXSize <= 0)
    1189             :         {
    1190           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    1191             :                      "%s metadata item must be > 0", MD_XSIZE);
    1192           6 :             return false;
    1193             :         }
    1194             : 
    1195          11 :         const int nYSize = atoi(pszYSize);
    1196          11 :         if (nYSize <= 0)
    1197             :         {
    1198           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    1199             :                      "%s metadata item must be > 0", MD_YSIZE);
    1200           1 :             return false;
    1201             :         }
    1202             : 
    1203             :         const CPLStringList aosTokens(
    1204          10 :             CSLTokenizeString2(pszGeoTransform, ",", 0));
    1205          10 :         if (aosTokens.size() != 6)
    1206             :         {
    1207           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    1208             :                      "%s metadata item must be 6 numeric values "
    1209             :                      "separated with comma",
    1210             :                      MD_GEOTRANSFORM);
    1211           1 :             return false;
    1212             :         }
    1213          63 :         for (int i = 0; i < 6; ++i)
    1214             :         {
    1215          54 :             m_adfGeoTransform[i] = CPLAtof(aosTokens[i]);
    1216             :         }
    1217           9 :         if (!(m_adfGeoTransform[GT_ROTATION_PARAM1] == 0))
    1218             :         {
    1219           1 :             CPLError(CE_Failure, CPLE_AppDefined, "3rd value of %s must be 0",
    1220             :                      MD_GEOTRANSFORM);
    1221           1 :             return false;
    1222             :         }
    1223           8 :         if (!(m_adfGeoTransform[GT_ROTATION_PARAM2] == 0))
    1224             :         {
    1225           1 :             CPLError(CE_Failure, CPLE_AppDefined, "5th value of %s must be 0",
    1226             :                      MD_GEOTRANSFORM);
    1227           1 :             return false;
    1228             :         }
    1229           7 :         if (!(m_adfGeoTransform[GT_NS_RES] < 0))
    1230             :         {
    1231           1 :             CPLError(CE_Failure, CPLE_AppDefined, "6th value of %s must be < 0",
    1232             :                      MD_GEOTRANSFORM);
    1233           1 :             return false;
    1234             :         }
    1235             : 
    1236           6 :         nRasterXSize = nXSize;
    1237          12 :         nRasterYSize = nYSize;
    1238             :     }
    1239         131 :     else if (pszResX && pszResY)
    1240             :     {
    1241          16 :         const double dfResX = CPLAtof(pszResX);
    1242          16 :         if (!(dfResX > 0))
    1243             :         {
    1244           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    1245             :                      "RESX metadata item must be > 0");
    1246           6 :             return false;
    1247             :         }
    1248          15 :         const double dfResY = CPLAtof(pszResY);
    1249          15 :         if (!(dfResY > 0))
    1250             :         {
    1251           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    1252             :                      "RESY metadata item must be > 0");
    1253           1 :             return false;
    1254             :         }
    1255             : 
    1256          14 :         OGREnvelope sEnvelope;
    1257             : 
    1258          14 :         if (nCountMinMaxXY == 4)
    1259             :         {
    1260           7 :             if (pszXSize || pszYSize || pszGeoTransform)
    1261             :             {
    1262           0 :                 CPLError(CE_Warning, CPLE_AppDefined,
    1263             :                          "Ignoring %s, %s and %s when %s, "
    1264             :                          "%s, %s and %s are specified",
    1265             :                          MD_XSIZE, MD_YSIZE, MD_GEOTRANSFORM, MD_MINX, MD_MINY,
    1266             :                          MD_MAXX, MD_MAXY);
    1267             :             }
    1268           7 :             const double dfMinX = CPLAtof(pszMinX);
    1269           7 :             const double dfMinY = CPLAtof(pszMinY);
    1270           7 :             const double dfMaxX = CPLAtof(pszMaxX);
    1271           7 :             const double dfMaxY = CPLAtof(pszMaxY);
    1272           7 :             if (!(dfMaxX > dfMinX))
    1273             :             {
    1274           1 :                 CPLError(CE_Failure, CPLE_AppDefined,
    1275             :                          "%s metadata item must be > %s", MD_MAXX, MD_MINX);
    1276           1 :                 return false;
    1277             :             }
    1278           6 :             if (!(dfMaxY > dfMinY))
    1279             :             {
    1280           1 :                 CPLError(CE_Failure, CPLE_AppDefined,
    1281             :                          "%s metadata item must be > %s", MD_MAXY, MD_MINY);
    1282           1 :                 return false;
    1283             :             }
    1284           5 :             sEnvelope.MinX = dfMinX;
    1285           5 :             sEnvelope.MinY = dfMinY;
    1286           5 :             sEnvelope.MaxX = dfMaxX;
    1287           5 :             sEnvelope.MaxY = dfMaxY;
    1288             :         }
    1289           7 :         else if (m_poLayer->GetExtent(&sEnvelope, /* bForce = */ false) ==
    1290             :                  OGRERR_FAILURE)
    1291             :         {
    1292           0 :             if (m_poLayer->GetExtent(&sEnvelope, /* bForce = */ true) ==
    1293             :                 OGRERR_FAILURE)
    1294             :             {
    1295           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
    1296             :                          "Cannot get layer extent");
    1297           0 :                 return false;
    1298             :             }
    1299           0 :             CPLError(CE_Warning, CPLE_AppDefined,
    1300             :                      "Could get layer extent, but using a slower method");
    1301             :         }
    1302             : 
    1303          12 :         const double dfXSize = (sEnvelope.MaxX - sEnvelope.MinX) / dfResX;
    1304          12 :         if (!(dfXSize >= 0 && dfXSize < INT_MAX))
    1305             :         {
    1306           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    1307             :                      "Too small %s, or wrong layer extent", MD_RESX);
    1308           1 :             return false;
    1309             :         }
    1310             : 
    1311          11 :         const double dfYSize = (sEnvelope.MaxY - sEnvelope.MinY) / dfResY;
    1312          11 :         if (!(dfYSize >= 0 && dfYSize < INT_MAX))
    1313             :         {
    1314           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    1315             :                      "Too small %s, or wrong layer extent", MD_RESY);
    1316           1 :             return false;
    1317             :         }
    1318             : 
    1319          10 :         m_adfGeoTransform[GT_TOPLEFT_X] = sEnvelope.MinX;
    1320          10 :         m_adfGeoTransform[GT_WE_RES] = dfResX;
    1321          10 :         m_adfGeoTransform[GT_ROTATION_PARAM1] = 0;
    1322          10 :         m_adfGeoTransform[GT_TOPLEFT_Y] = sEnvelope.MaxY;
    1323          10 :         m_adfGeoTransform[GT_ROTATION_PARAM2] = 0;
    1324          10 :         m_adfGeoTransform[GT_NS_RES] = -dfResY;
    1325          10 :         nRasterXSize = static_cast<int>(std::ceil(dfXSize));
    1326          10 :         nRasterYSize = static_cast<int>(std::ceil(dfYSize));
    1327             :     }
    1328             : 
    1329         131 :     if (nBandCount == 0 && !pszBandCount)
    1330             :     {
    1331           0 :         CPLError(CE_Failure, CPLE_AppDefined, "%s metadata item missing",
    1332             :                  MD_BAND_COUNT);
    1333           0 :         return false;
    1334             :     }
    1335             : 
    1336         131 :     if (!GDALCheckBandCount(nBandCount, false))
    1337           1 :         return false;
    1338             : 
    1339         130 :     if (aeDataTypes.empty() && !pszDataType)
    1340             :     {
    1341           9 :         aeDataTypes.resize(nBandCount, GDT_Byte);
    1342             :     }
    1343         121 :     else if (pszDataType)
    1344             :     {
    1345           8 :         aeDataTypes.clear();
    1346           8 :         const CPLStringList aosTokens(CSLTokenizeString2(pszDataType, ", ", 0));
    1347           8 :         if (aosTokens.size() == 1)
    1348             :         {
    1349           6 :             const auto eDataType = GDALGetDataTypeByName(aosTokens[0]);
    1350           6 :             if (eDataType == GDT_Unknown)
    1351             :             {
    1352           1 :                 CPLError(CE_Failure, CPLE_AppDefined, "Invalid value for %s",
    1353             :                          MD_DATA_TYPE);
    1354           1 :                 return false;
    1355             :             }
    1356           5 :             aeDataTypes.resize(nBandCount, eDataType);
    1357             :         }
    1358           2 :         else if (aosTokens.size() == nBandCount)
    1359             :         {
    1360           2 :             for (int i = 0; i < nBandCount; ++i)
    1361             :             {
    1362           2 :                 const auto eDataType = GDALGetDataTypeByName(aosTokens[i]);
    1363           2 :                 if (eDataType == GDT_Unknown)
    1364             :                 {
    1365           1 :                     CPLError(CE_Failure, CPLE_AppDefined,
    1366             :                              "Invalid value for %s", MD_DATA_TYPE);
    1367           1 :                     return false;
    1368             :                 }
    1369           1 :                 aeDataTypes.push_back(eDataType);
    1370             :             }
    1371             :         }
    1372             :         else
    1373             :         {
    1374           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    1375             :                      "Number of values in %s must be 1 or %s", MD_DATA_TYPE,
    1376             :                      MD_BAND_COUNT);
    1377           1 :             return false;
    1378             :         }
    1379             :     }
    1380             : 
    1381         127 :     const char *pszNoData = GetOption(MD_NODATA);
    1382         127 :     if (pszNoData)
    1383             :     {
    1384          20 :         const auto IsValidNoDataStr = [](const char *pszStr)
    1385             :         {
    1386          20 :             if (EQUAL(pszStr, "inf") || EQUAL(pszStr, "-inf") ||
    1387          16 :                 EQUAL(pszStr, "nan"))
    1388           6 :                 return true;
    1389          14 :             const auto eType = CPLGetValueType(pszStr);
    1390          14 :             return eType == CPL_VALUE_INTEGER || eType == CPL_VALUE_REAL;
    1391             :         };
    1392             : 
    1393          18 :         aNoData.clear();
    1394          18 :         const CPLStringList aosTokens(CSLTokenizeString2(pszNoData, ", ", 0));
    1395          18 :         if (aosTokens.size() == 1)
    1396             :         {
    1397          14 :             if (!EQUAL(aosTokens[0], "NONE"))
    1398             :             {
    1399          11 :                 if (!IsValidNoDataStr(aosTokens[0]))
    1400             :                 {
    1401           1 :                     CPLError(CE_Failure, CPLE_AppDefined,
    1402             :                              "Invalid value for %s", MD_NODATA);
    1403           1 :                     return false;
    1404             :                 }
    1405          10 :                 aNoData.resize(nBandCount,
    1406          20 :                                std::pair(true, CPLAtof(aosTokens[0])));
    1407             :             }
    1408             :         }
    1409           4 :         else if (aosTokens.size() == nBandCount)
    1410             :         {
    1411          12 :             for (int i = 0; i < nBandCount; ++i)
    1412             :             {
    1413          10 :                 if (EQUAL(aosTokens[i], "NONE"))
    1414             :                 {
    1415           1 :                     aNoData.emplace_back(false, 0);
    1416             :                 }
    1417           9 :                 else if (IsValidNoDataStr(aosTokens[i]))
    1418             :                 {
    1419           8 :                     aNoData.emplace_back(true, CPLAtof(aosTokens[i]));
    1420             :                 }
    1421             :                 else
    1422             :                 {
    1423           1 :                     CPLError(CE_Failure, CPLE_AppDefined,
    1424             :                              "Invalid value for %s", MD_NODATA);
    1425           1 :                     return false;
    1426             :                 }
    1427             :             }
    1428             :         }
    1429             :         else
    1430             :         {
    1431           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    1432             :                      "Number of values in %s must be 1 or %s", MD_NODATA,
    1433             :                      MD_BAND_COUNT);
    1434           1 :             return false;
    1435             :         }
    1436             :     }
    1437             : 
    1438         124 :     if (pszColorInterp)
    1439             :     {
    1440          11 :         aeColorInterp.clear();
    1441             :         const CPLStringList aosTokens(
    1442          11 :             CSLTokenizeString2(pszColorInterp, ", ", 0));
    1443          11 :         if (aosTokens.size() == 1)
    1444             :         {
    1445           7 :             const auto eInterp = GDALGetColorInterpretationByName(aosTokens[0]);
    1446          12 :             if (eInterp == GCI_Undefined &&
    1447           5 :                 !EQUAL(aosTokens[0],
    1448             :                        GDALGetColorInterpretationName(GCI_Undefined)))
    1449             :             {
    1450           1 :                 CPLError(CE_Failure, CPLE_AppDefined, "Invalid value for %s",
    1451             :                          MD_COLOR_INTERPRETATION);
    1452           1 :                 return false;
    1453             :             }
    1454           6 :             aeColorInterp.resize(nBandCount, eInterp);
    1455             :         }
    1456           4 :         else if (aosTokens.size() == nBandCount)
    1457             :         {
    1458          11 :             for (int i = 0; i < nBandCount; ++i)
    1459             :             {
    1460             :                 const auto eInterp =
    1461           9 :                     GDALGetColorInterpretationByName(aosTokens[i]);
    1462          11 :                 if (eInterp == GCI_Undefined &&
    1463           2 :                     !EQUAL(aosTokens[i],
    1464             :                            GDALGetColorInterpretationName(GCI_Undefined)))
    1465             :                 {
    1466           1 :                     CPLError(CE_Failure, CPLE_AppDefined,
    1467             :                              "Invalid value for %s", MD_COLOR_INTERPRETATION);
    1468           1 :                     return false;
    1469             :                 }
    1470           8 :                 aeColorInterp.emplace_back(eInterp);
    1471             :             }
    1472             :         }
    1473             :         else
    1474             :         {
    1475           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    1476             :                      "Number of values in %s must be 1 or "
    1477             :                      "%s",
    1478             :                      MD_COLOR_INTERPRETATION, MD_BAND_COUNT);
    1479           1 :             return false;
    1480             :         }
    1481             :     }
    1482             : 
    1483             :     /* -------------------------------------------------------------------- */
    1484             :     /*      Create bands.                                                   */
    1485             :     /* -------------------------------------------------------------------- */
    1486         121 :     if (aeDataTypes.size() != static_cast<size_t>(nBandCount))
    1487             :     {
    1488           1 :         CPLError(
    1489             :             CE_Failure, CPLE_AppDefined,
    1490             :             "Number of data types values found not matching number of bands");
    1491           1 :         return false;
    1492             :     }
    1493         120 :     if (!aNoData.empty() && aNoData.size() != static_cast<size_t>(nBandCount))
    1494             :     {
    1495           1 :         CPLError(CE_Failure, CPLE_AppDefined,
    1496             :                  "Number of nodata values found not matching number of bands");
    1497           1 :         return false;
    1498             :     }
    1499         230 :     if (!aeColorInterp.empty() &&
    1500         111 :         aeColorInterp.size() != static_cast<size_t>(nBandCount))
    1501             :     {
    1502           1 :         CPLError(CE_Failure, CPLE_AppDefined,
    1503             :                  "Number of color interpretation values found not matching "
    1504             :                  "number of bands");
    1505           1 :         return false;
    1506             :     }
    1507             : 
    1508         118 :     int nBlockXSize = 256;
    1509         118 :     const char *pszBlockXSize = GetOption(MD_BLOCK_X_SIZE);
    1510         118 :     if (pszBlockXSize)
    1511             :     {
    1512           3 :         nBlockXSize = atoi(pszBlockXSize);
    1513           3 :         if (nBlockXSize <= 0)
    1514             :         {
    1515           1 :             CPLError(CE_Failure, CPLE_AppDefined, "Invalid %s",
    1516             :                      MD_BLOCK_X_SIZE);
    1517           1 :             return false;
    1518             :         }
    1519             :     }
    1520             : 
    1521         117 :     int nBlockYSize = 256;
    1522         117 :     const char *pszBlockYSize = GetOption(MD_BLOCK_Y_SIZE);
    1523         117 :     if (pszBlockYSize)
    1524             :     {
    1525           3 :         nBlockYSize = atoi(pszBlockYSize);
    1526           3 :         if (nBlockYSize <= 0)
    1527             :         {
    1528           1 :             CPLError(CE_Failure, CPLE_AppDefined, "Invalid %s",
    1529             :                      MD_BLOCK_Y_SIZE);
    1530           1 :             return false;
    1531             :         }
    1532             :     }
    1533             : 
    1534         116 :     if (nBlockXSize > INT_MAX / nBlockYSize)
    1535             :     {
    1536           1 :         CPLError(CE_Failure, CPLE_AppDefined, "Too big %s * %s",
    1537             :                  MD_BLOCK_X_SIZE, MD_BLOCK_Y_SIZE);
    1538           1 :         return false;
    1539             :     }
    1540             : 
    1541         115 :     if (dfOvrFactor > 1.0)
    1542             :     {
    1543           4 :         m_adfGeoTransform[GT_WE_RES] *= dfOvrFactor;
    1544           4 :         m_adfGeoTransform[GT_NS_RES] *= dfOvrFactor;
    1545           4 :         nRasterXSize = static_cast<int>(std::ceil(nRasterXSize / dfOvrFactor));
    1546           4 :         nRasterYSize = static_cast<int>(std::ceil(nRasterYSize / dfOvrFactor));
    1547             :     }
    1548             : 
    1549         115 :     GDALTileIndexBand *poFirstBand = nullptr;
    1550         261 :     for (int i = 0; i < nBandCount; ++i)
    1551             :     {
    1552         146 :         GDALDataType eDataType = aeDataTypes[i];
    1553         146 :         if (!apoXMLNodeBands.empty())
    1554             :         {
    1555           4 :             const char *pszVal = CPLGetXMLValue(apoXMLNodeBands[i],
    1556             :                                                 GTI_XML_BAND_DATATYPE, nullptr);
    1557           4 :             if (pszVal)
    1558             :             {
    1559           3 :                 eDataType = GDALGetDataTypeByName(pszVal);
    1560           3 :                 if (eDataType == GDT_Unknown)
    1561           0 :                     return false;
    1562             :             }
    1563             :         }
    1564             :         auto poBandUniquePtr = std::make_unique<GDALTileIndexBand>(
    1565         292 :             this, i + 1, eDataType, nBlockXSize, nBlockYSize);
    1566         146 :         auto poBand = poBandUniquePtr.get();
    1567         146 :         SetBand(i + 1, poBandUniquePtr.release());
    1568         146 :         if (!poFirstBand)
    1569         115 :             poFirstBand = poBand;
    1570         146 :         if (poBand->GetRasterDataType() != poFirstBand->GetRasterDataType())
    1571             :         {
    1572           0 :             m_bSameDataType = false;
    1573             :         }
    1574             : 
    1575         146 :         if (!apoXMLNodeBands.empty())
    1576             :         {
    1577           4 :             const char *pszVal = CPLGetXMLValue(
    1578           4 :                 apoXMLNodeBands[i], GTI_XML_BAND_DESCRIPTION, nullptr);
    1579           4 :             if (pszVal)
    1580             :             {
    1581           2 :                 poBand->GDALRasterBand::SetDescription(pszVal);
    1582             :             }
    1583             :         }
    1584             : 
    1585         146 :         if (!aNoData.empty() && aNoData[i].first)
    1586             :         {
    1587          21 :             poBand->m_bNoDataValueSet = true;
    1588          21 :             poBand->m_dfNoDataValue = aNoData[i].second;
    1589             :         }
    1590         146 :         if (!apoXMLNodeBands.empty())
    1591             :         {
    1592           4 :             const char *pszVal = CPLGetXMLValue(
    1593           4 :                 apoXMLNodeBands[i], GTI_XML_BAND_NODATAVALUE, nullptr);
    1594           4 :             if (pszVal)
    1595             :             {
    1596           3 :                 poBand->m_bNoDataValueSet = true;
    1597           3 :                 poBand->m_dfNoDataValue = CPLAtof(pszVal);
    1598             :             }
    1599             :         }
    1600         291 :         if (poBand->m_bNoDataValueSet != poFirstBand->m_bNoDataValueSet ||
    1601         145 :             !IsSameNaNAware(poBand->m_dfNoDataValue,
    1602             :                             poFirstBand->m_dfNoDataValue))
    1603             :         {
    1604           6 :             m_bSameNoData = false;
    1605             :         }
    1606             : 
    1607         146 :         if (!aeColorInterp.empty())
    1608             :         {
    1609         135 :             poBand->m_eColorInterp = aeColorInterp[i];
    1610             :         }
    1611         146 :         if (!apoXMLNodeBands.empty())
    1612             :         {
    1613           4 :             const char *pszVal = CPLGetXMLValue(
    1614           4 :                 apoXMLNodeBands[i], GTI_XML_BAND_COLORINTERP, nullptr);
    1615           4 :             if (pszVal)
    1616             :             {
    1617           4 :                 poBand->m_eColorInterp =
    1618           4 :                     GDALGetColorInterpretationByName(pszVal);
    1619             :             }
    1620             :         }
    1621             : 
    1622         146 :         if (const char *pszScale =
    1623         146 :                 GetOption(CPLSPrintf("BAND_%d_%s", i + 1, MD_BAND_SCALE)))
    1624             :         {
    1625           6 :             poBand->m_dfScale = CPLAtof(pszScale);
    1626             :         }
    1627         146 :         if (!apoXMLNodeBands.empty())
    1628             :         {
    1629             :             const char *pszVal =
    1630           4 :                 CPLGetXMLValue(apoXMLNodeBands[i], GTI_XML_BAND_SCALE, nullptr);
    1631           4 :             if (pszVal)
    1632             :             {
    1633           2 :                 poBand->m_dfScale = CPLAtof(pszVal);
    1634             :             }
    1635             :         }
    1636             : 
    1637         146 :         if (const char *pszOffset =
    1638         146 :                 GetOption(CPLSPrintf("BAND_%d_%s", i + 1, MD_BAND_OFFSET)))
    1639             :         {
    1640           6 :             poBand->m_dfOffset = CPLAtof(pszOffset);
    1641             :         }
    1642         146 :         if (!apoXMLNodeBands.empty())
    1643             :         {
    1644           4 :             const char *pszVal = CPLGetXMLValue(apoXMLNodeBands[i],
    1645             :                                                 GTI_XML_BAND_OFFSET, nullptr);
    1646           4 :             if (pszVal)
    1647             :             {
    1648           2 :                 poBand->m_dfOffset = CPLAtof(pszVal);
    1649             :             }
    1650             :         }
    1651             : 
    1652         146 :         if (const char *pszUnit =
    1653         146 :                 GetOption(CPLSPrintf("BAND_%d_%s", i + 1, MD_BAND_UNITTYPE)))
    1654             :         {
    1655           6 :             poBand->m_osUnit = pszUnit;
    1656             :         }
    1657         146 :         if (!apoXMLNodeBands.empty())
    1658             :         {
    1659           4 :             const char *pszVal = CPLGetXMLValue(apoXMLNodeBands[i],
    1660             :                                                 GTI_XML_BAND_UNITTYPE, nullptr);
    1661           4 :             if (pszVal)
    1662             :             {
    1663           2 :                 poBand->m_osUnit = pszVal;
    1664             :             }
    1665             :         }
    1666             : 
    1667         146 :         if (!apoXMLNodeBands.empty())
    1668             :         {
    1669           4 :             const CPLXMLNode *psBandNode = apoXMLNodeBands[i];
    1670           4 :             poBand->oMDMD.XMLInit(psBandNode, TRUE);
    1671             : 
    1672           4 :             if (const CPLXMLNode *psCategoryNames =
    1673           4 :                     CPLGetXMLNode(psBandNode, GTI_XML_CATEGORYNAMES))
    1674             :             {
    1675             :                 poBand->m_aosCategoryNames =
    1676           2 :                     VRTParseCategoryNames(psCategoryNames);
    1677             :             }
    1678             : 
    1679           4 :             if (const CPLXMLNode *psColorTable =
    1680           4 :                     CPLGetXMLNode(psBandNode, GTI_XML_COLORTABLE))
    1681             :             {
    1682           2 :                 poBand->m_poColorTable = VRTParseColorTable(psColorTable);
    1683             :             }
    1684             : 
    1685           4 :             if (const CPLXMLNode *psRAT =
    1686           4 :                     CPLGetXMLNode(psBandNode, GTI_XML_RAT))
    1687             :             {
    1688             :                 poBand->m_poRAT =
    1689           2 :                     std::make_unique<GDALDefaultRasterAttributeTable>();
    1690           2 :                 poBand->m_poRAT->XMLInit(psRAT, "");
    1691             :             }
    1692             :         }
    1693             :     }
    1694             : 
    1695         115 :     const char *pszMaskBand = GetOption(MD_MASK_BAND);
    1696         115 :     if (pszMaskBand)
    1697           7 :         bHasMaskBand = CPLTestBool(pszMaskBand);
    1698         115 :     if (bHasMaskBand)
    1699             :     {
    1700           8 :         m_poMaskBand = std::make_unique<GDALTileIndexBand>(
    1701          16 :             this, 0, GDT_Byte, nBlockXSize, nBlockYSize);
    1702             :     }
    1703             : 
    1704         115 :     if (dfOvrFactor == 1.0)
    1705             :     {
    1706         111 :         if (psRoot)
    1707             :         {
    1708          58 :             for (const CPLXMLNode *psIter = psRoot->psChild; psIter;
    1709          43 :                  psIter = psIter->psNext)
    1710             :             {
    1711          44 :                 if (psIter->eType == CXT_Element &&
    1712          44 :                     strcmp(psIter->pszValue, GTI_XML_OVERVIEW_ELEMENT) == 0)
    1713             :                 {
    1714           8 :                     const char *pszDataset = CPLGetXMLValue(
    1715             :                         psIter, GTI_XML_OVERVIEW_DATASET, nullptr);
    1716             :                     const char *pszLayer =
    1717           8 :                         CPLGetXMLValue(psIter, GTI_XML_OVERVIEW_LAYER, nullptr);
    1718           8 :                     const char *pszFactor = CPLGetXMLValue(
    1719             :                         psIter, GTI_XML_OVERVIEW_FACTOR, nullptr);
    1720           8 :                     if (!pszDataset && !pszLayer && !pszFactor)
    1721             :                     {
    1722           1 :                         CPLError(
    1723             :                             CE_Failure, CPLE_AppDefined,
    1724             :                             "At least one of %s, %s or %s element "
    1725             :                             "must be present as an %s child",
    1726             :                             GTI_XML_OVERVIEW_DATASET, GTI_XML_OVERVIEW_LAYER,
    1727             :                             GTI_XML_OVERVIEW_FACTOR, GTI_XML_OVERVIEW_ELEMENT);
    1728           1 :                         return false;
    1729             :                     }
    1730             :                     m_aoOverviewDescriptor.emplace_back(
    1731          14 :                         std::string(pszDataset ? pszDataset : ""),
    1732          14 :                         CPLStringList(
    1733             :                             GDALDeserializeOpenOptionsFromXML(psIter)),
    1734          14 :                         std::string(pszLayer ? pszLayer : ""),
    1735          21 :                         pszFactor ? CPLAtof(pszFactor) : 0.0);
    1736             :                 }
    1737             :             }
    1738             :         }
    1739             :         else
    1740             :         {
    1741          96 :             for (int iOvr = 1;; ++iOvr)
    1742             :             {
    1743             :                 const char *pszOvrDSName =
    1744         101 :                     GetOption(CPLSPrintf("OVERVIEW_%d_DATASET", iOvr));
    1745             :                 const char *pszOpenOptions =
    1746         101 :                     GetOption(CPLSPrintf("OVERVIEW_%d_OPEN_OPTIONS", iOvr));
    1747             :                 const char *pszOvrLayer =
    1748         101 :                     GetOption(CPLSPrintf("OVERVIEW_%d_LAYER", iOvr));
    1749             :                 const char *pszOvrFactor =
    1750         101 :                     GetOption(CPLSPrintf("OVERVIEW_%d_FACTOR", iOvr));
    1751         101 :                 if (!pszOvrDSName && !pszOvrLayer && !pszOvrFactor)
    1752          96 :                     break;
    1753             :                 m_aoOverviewDescriptor.emplace_back(
    1754          10 :                     std::string(pszOvrDSName ? pszOvrDSName : ""),
    1755          10 :                     pszOpenOptions ? CPLStringList(CSLTokenizeString2(
    1756             :                                          pszOpenOptions, ",", 0))
    1757             :                                    : CPLStringList(),
    1758          10 :                     std::string(pszOvrLayer ? pszOvrLayer : ""),
    1759          10 :                     pszOvrFactor ? CPLAtof(pszOvrFactor) : 0.0);
    1760           5 :             }
    1761             :         }
    1762             :     }
    1763             : 
    1764         114 :     if (psRoot)
    1765             :     {
    1766          16 :         oMDMD.XMLInit(psRoot, TRUE);
    1767             :     }
    1768             :     else
    1769             :     {
    1770             :         // Set on the dataset all metadata items from the index layer which are
    1771             :         // not "reserved" keywords.
    1772          98 :         CSLConstList papszLayerMD = m_poLayer->GetMetadata();
    1773         376 :         for (const auto &[pszKey, pszValue] :
    1774         474 :              cpl::IterateNameValue(papszLayerMD))
    1775             :         {
    1776         188 :             if (STARTS_WITH_CI(pszKey, "OVERVIEW_"))
    1777             :             {
    1778           7 :                 continue;
    1779             :             }
    1780             : 
    1781         181 :             bool bIsVRTItem = false;
    1782        2620 :             for (const char *pszTest : apszTIOptions)
    1783             :             {
    1784        2560 :                 if (EQUAL(pszKey, pszTest))
    1785             :                 {
    1786         121 :                     bIsVRTItem = true;
    1787         121 :                     break;
    1788             :                 }
    1789             :             }
    1790         181 :             if (!bIsVRTItem)
    1791             :             {
    1792          60 :                 if (STARTS_WITH_CI(pszKey, "BAND_"))
    1793             :                 {
    1794          52 :                     const int nBandNr = atoi(pszKey + strlen("BAND_"));
    1795             :                     const char *pszNextUnderscore =
    1796          52 :                         strchr(pszKey + strlen("BAND_"), '_');
    1797          52 :                     if (pszNextUnderscore && nBandNr >= 1 && nBandNr <= nBands)
    1798             :                     {
    1799          42 :                         const char *pszKeyWithoutBand = pszNextUnderscore + 1;
    1800          42 :                         bool bIsReservedBandItem = false;
    1801         132 :                         for (const char *pszItem : apszReservedBandItems)
    1802             :                         {
    1803         108 :                             if (EQUAL(pszKeyWithoutBand, pszItem))
    1804             :                             {
    1805          18 :                                 bIsReservedBandItem = true;
    1806          18 :                                 break;
    1807             :                             }
    1808             :                         }
    1809          42 :                         if (!bIsReservedBandItem)
    1810             :                         {
    1811          24 :                             GetRasterBand(nBandNr)
    1812          24 :                                 ->GDALRasterBand::SetMetadataItem(
    1813             :                                     pszKeyWithoutBand, pszValue);
    1814             :                         }
    1815             :                     }
    1816             :                 }
    1817             :                 else
    1818             :                 {
    1819           8 :                     GDALDataset::SetMetadataItem(pszKey, pszValue);
    1820             :                 }
    1821             :             }
    1822             :         }
    1823             :     }
    1824             : 
    1825         114 :     if (nBandCount > 1 && !GetMetadata("IMAGE_STRUCTURE"))
    1826             :     {
    1827          16 :         GDALDataset::SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");
    1828             :     }
    1829             : 
    1830             :     /* -------------------------------------------------------------------- */
    1831             :     /*      Initialize any PAM information.                                 */
    1832             :     /* -------------------------------------------------------------------- */
    1833         114 :     SetDescription(poOpenInfo->pszFilename);
    1834         114 :     TryLoadXML();
    1835             : 
    1836             :     /* -------------------------------------------------------------------- */
    1837             :     /*      Check for overviews.                                            */
    1838             :     /* -------------------------------------------------------------------- */
    1839         114 :     oOvManager.Initialize(this, poOpenInfo->pszFilename);
    1840             : 
    1841         114 :     return true;
    1842             : }
    1843             : 
    1844             : /************************************************************************/
    1845             : /*                        GetMetadataItem()                             */
    1846             : /************************************************************************/
    1847             : 
    1848          34 : const char *GDALTileIndexDataset::GetMetadataItem(const char *pszName,
    1849             :                                                   const char *pszDomain)
    1850             : {
    1851          34 :     if (pszName && pszDomain && EQUAL(pszDomain, "__DEBUG__"))
    1852             :     {
    1853           9 :         if (EQUAL(pszName, "SCANNED_ONE_FEATURE_AT_OPENING"))
    1854             :         {
    1855           4 :             return m_bScannedOneFeatureAtOpening ? "YES" : "NO";
    1856             :         }
    1857           5 :         else if (EQUAL(pszName, "NUMBER_OF_CONTRIBUTING_SOURCES"))
    1858             :         {
    1859           5 :             return CPLSPrintf("%d", static_cast<int>(m_aoSourceDesc.size()));
    1860             :         }
    1861             :     }
    1862          25 :     return GDALPamDataset::GetMetadataItem(pszName, pszDomain);
    1863             : }
    1864             : 
    1865             : /************************************************************************/
    1866             : /*                TileIndexSupportsEditingLayerMetadata()               */
    1867             : /************************************************************************/
    1868             : 
    1869          17 : bool GDALTileIndexDataset::TileIndexSupportsEditingLayerMetadata() const
    1870             : {
    1871          27 :     return eAccess == GA_Update && m_poVectorDS->GetDriver() &&
    1872          27 :            EQUAL(m_poVectorDS->GetDriver()->GetDescription(), "GPKG");
    1873             : }
    1874             : 
    1875             : /************************************************************************/
    1876             : /*                        SetMetadataItem()                             */
    1877             : /************************************************************************/
    1878             : 
    1879           3 : CPLErr GDALTileIndexDataset::SetMetadataItem(const char *pszName,
    1880             :                                              const char *pszValue,
    1881             :                                              const char *pszDomain)
    1882             : {
    1883           3 :     if (m_bXMLUpdatable)
    1884             :     {
    1885           1 :         m_bXMLModified = true;
    1886           1 :         return GDALDataset::SetMetadataItem(pszName, pszValue, pszDomain);
    1887             :     }
    1888           2 :     else if (TileIndexSupportsEditingLayerMetadata())
    1889             :     {
    1890           1 :         m_poLayer->SetMetadataItem(pszName, pszValue, pszDomain);
    1891           1 :         return GDALDataset::SetMetadataItem(pszName, pszValue, pszDomain);
    1892             :     }
    1893             :     else
    1894             :     {
    1895           1 :         return GDALPamDataset::SetMetadataItem(pszName, pszValue, pszDomain);
    1896             :     }
    1897             : }
    1898             : 
    1899             : /************************************************************************/
    1900             : /*                           SetMetadata()                              */
    1901             : /************************************************************************/
    1902             : 
    1903           3 : CPLErr GDALTileIndexDataset::SetMetadata(char **papszMD, const char *pszDomain)
    1904             : {
    1905           3 :     if (m_bXMLUpdatable)
    1906             :     {
    1907           1 :         m_bXMLModified = true;
    1908           1 :         return GDALDataset::SetMetadata(papszMD, pszDomain);
    1909             :     }
    1910           2 :     else if (TileIndexSupportsEditingLayerMetadata())
    1911             :     {
    1912           2 :         if (!pszDomain || pszDomain[0] == 0)
    1913             :         {
    1914           4 :             CPLStringList aosMD(CSLDuplicate(papszMD));
    1915             : 
    1916             :             // Reinject dataset reserved items
    1917          44 :             for (const char *pszItem : apszTIOptions)
    1918             :             {
    1919          42 :                 if (!aosMD.FetchNameValue(pszItem))
    1920             :                 {
    1921          42 :                     const char *pszValue = m_poLayer->GetMetadataItem(pszItem);
    1922          42 :                     if (pszValue)
    1923             :                     {
    1924           2 :                         aosMD.SetNameValue(pszItem, pszValue);
    1925             :                     }
    1926             :                 }
    1927             :             }
    1928             : 
    1929             :             // Reinject band metadata
    1930           2 :             char **papszExistingLayerMD = m_poLayer->GetMetadata();
    1931          17 :             for (int i = 0; papszExistingLayerMD && papszExistingLayerMD[i];
    1932             :                  ++i)
    1933             :             {
    1934          15 :                 if (STARTS_WITH_CI(papszExistingLayerMD[i], "BAND_"))
    1935             :                 {
    1936          12 :                     aosMD.AddString(papszExistingLayerMD[i]);
    1937             :                 }
    1938             :             }
    1939             : 
    1940           4 :             m_poLayer->SetMetadata(aosMD.List(), pszDomain);
    1941             :         }
    1942             :         else
    1943             :         {
    1944           0 :             m_poLayer->SetMetadata(papszMD, pszDomain);
    1945             :         }
    1946           2 :         return GDALDataset::SetMetadata(papszMD, pszDomain);
    1947             :     }
    1948             :     else
    1949             :     {
    1950           0 :         return GDALPamDataset::SetMetadata(papszMD, pszDomain);
    1951             :     }
    1952             : }
    1953             : 
    1954             : /************************************************************************/
    1955             : /*                     GDALTileIndexDatasetIdentify()                   */
    1956             : /************************************************************************/
    1957             : 
    1958       71352 : static int GDALTileIndexDatasetIdentify(GDALOpenInfo *poOpenInfo)
    1959             : {
    1960       71352 :     if (STARTS_WITH(poOpenInfo->pszFilename, GTI_PREFIX))
    1961           6 :         return true;
    1962             : 
    1963       71346 :     if (STARTS_WITH(poOpenInfo->pszFilename, "<GDALTileIndexDataset"))
    1964          44 :         return true;
    1965             : 
    1966      166708 :     if (poOpenInfo->nHeaderBytes >= 100 &&
    1967       24104 :         STARTS_WITH(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
    1968         726 :                     "SQLite format 3") &&
    1969       96132 :         ENDS_WITH_CI(poOpenInfo->pszFilename, ".gti.gpkg") &&
    1970         341 :         !STARTS_WITH(poOpenInfo->pszFilename, "GPKG:"))
    1971             :     {
    1972             :         // Most likely handled by GTI driver, but we can't be sure
    1973         341 :         return GDAL_IDENTIFY_UNKNOWN;
    1974             :     }
    1975             : 
    1976       95288 :     return poOpenInfo->nHeaderBytes > 0 &&
    1977       94391 :            (poOpenInfo->nOpenFlags & GDAL_OF_RASTER) != 0 &&
    1978       23441 :            (strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
    1979       23424 :                    "<GDALTileIndexDataset") ||
    1980       46857 :             ENDS_WITH_CI(poOpenInfo->pszFilename, ".gti.fgb") ||
    1981       94380 :             ENDS_WITH_CI(poOpenInfo->pszFilename, ".gti.parquet"));
    1982             : }
    1983             : 
    1984             : /************************************************************************/
    1985             : /*                      GDALTileIndexDatasetOpen()                       */
    1986             : /************************************************************************/
    1987             : 
    1988         185 : static GDALDataset *GDALTileIndexDatasetOpen(GDALOpenInfo *poOpenInfo)
    1989             : {
    1990         185 :     if (GDALTileIndexDatasetIdentify(poOpenInfo) == GDAL_IDENTIFY_FALSE)
    1991           0 :         return nullptr;
    1992         370 :     auto poDS = std::make_unique<GDALTileIndexDataset>();
    1993         185 :     if (!poDS->Open(poOpenInfo))
    1994          71 :         return nullptr;
    1995         114 :     return poDS.release();
    1996             : }
    1997             : 
    1998             : /************************************************************************/
    1999             : /*                          ~GDALTileIndexDataset()                      */
    2000             : /************************************************************************/
    2001             : 
    2002         370 : GDALTileIndexDataset::~GDALTileIndexDataset()
    2003             : {
    2004         185 :     GDALTileIndexDataset::FlushCache(true);
    2005         370 : }
    2006             : 
    2007             : /************************************************************************/
    2008             : /*                              FlushCache()                            */
    2009             : /************************************************************************/
    2010             : 
    2011         186 : CPLErr GDALTileIndexDataset::FlushCache(bool bAtClosing)
    2012             : {
    2013         186 :     CPLErr eErr = CE_None;
    2014         186 :     if (bAtClosing && m_bXMLModified)
    2015             :     {
    2016             :         CPLXMLNode *psRoot =
    2017           1 :             CPLGetXMLNode(m_psXMLTree.get(), "=GDALTileIndexDataset");
    2018             : 
    2019             :         // Suppress existing dataset metadata
    2020             :         while (true)
    2021             :         {
    2022           1 :             CPLXMLNode *psExistingMetadata = CPLGetXMLNode(psRoot, "Metadata");
    2023           1 :             if (!psExistingMetadata)
    2024           1 :                 break;
    2025           0 :             CPLRemoveXMLChild(psRoot, psExistingMetadata);
    2026           0 :         }
    2027             : 
    2028             :         // Serialize new dataset metadata
    2029           1 :         if (CPLXMLNode *psMD = oMDMD.Serialize())
    2030           1 :             CPLAddXMLChild(psRoot, psMD);
    2031             : 
    2032             :         // Update existing band metadata
    2033           1 :         if (CPLGetXMLNode(psRoot, GTI_XML_BAND_ELEMENT))
    2034             :         {
    2035           0 :             for (CPLXMLNode *psIter = psRoot->psChild; psIter;
    2036           0 :                  psIter = psIter->psNext)
    2037             :             {
    2038           0 :                 if (psIter->eType == CXT_Element &&
    2039           0 :                     strcmp(psIter->pszValue, GTI_XML_BAND_ELEMENT))
    2040             :                 {
    2041             :                     const char *pszBand =
    2042           0 :                         CPLGetXMLValue(psIter, GTI_XML_BAND_NUMBER, nullptr);
    2043           0 :                     if (pszBand)
    2044             :                     {
    2045           0 :                         const int nBand = atoi(pszBand);
    2046           0 :                         if (nBand >= 1 && nBand <= nBands)
    2047             :                         {
    2048             :                             while (true)
    2049             :                             {
    2050             :                                 CPLXMLNode *psExistingMetadata =
    2051           0 :                                     CPLGetXMLNode(psIter, "Metadata");
    2052           0 :                                 if (!psExistingMetadata)
    2053           0 :                                     break;
    2054           0 :                                 CPLRemoveXMLChild(psIter, psExistingMetadata);
    2055           0 :                             }
    2056             : 
    2057           0 :                             auto poBand = cpl::down_cast<GDALTileIndexBand *>(
    2058           0 :                                 papoBands[nBand - 1]);
    2059           0 :                             if (CPLXMLNode *psMD = poBand->oMDMD.Serialize())
    2060           0 :                                 CPLAddXMLChild(psIter, psMD);
    2061             :                         }
    2062             :                     }
    2063             :                 }
    2064             :             }
    2065             :         }
    2066             :         else
    2067             :         {
    2068             :             // Create new band objects if they have metadata
    2069           2 :             std::vector<CPLXMLTreeCloser> aoBandXML;
    2070           1 :             bool bHasBandMD = false;
    2071           2 :             for (int i = 1; i <= nBands; ++i)
    2072             :             {
    2073             :                 auto poBand =
    2074           1 :                     cpl::down_cast<GDALTileIndexBand *>(papoBands[i - 1]);
    2075           1 :                 auto psMD = poBand->oMDMD.Serialize();
    2076           1 :                 if (psMD)
    2077           1 :                     bHasBandMD = true;
    2078           1 :                 aoBandXML.emplace_back(CPLXMLTreeCloser(psMD));
    2079             :             }
    2080           1 :             if (bHasBandMD)
    2081             :             {
    2082           2 :                 for (int i = 1; i <= nBands; ++i)
    2083             :                 {
    2084             :                     auto poBand =
    2085           1 :                         cpl::down_cast<GDALTileIndexBand *>(papoBands[i - 1]);
    2086             : 
    2087           1 :                     CPLXMLNode *psBand = CPLCreateXMLNode(psRoot, CXT_Element,
    2088             :                                                           GTI_XML_BAND_ELEMENT);
    2089           1 :                     CPLAddXMLAttributeAndValue(psBand, GTI_XML_BAND_NUMBER,
    2090             :                                                CPLSPrintf("%d", i));
    2091           1 :                     CPLAddXMLAttributeAndValue(
    2092             :                         psBand, GTI_XML_BAND_DATATYPE,
    2093             :                         GDALGetDataTypeName(poBand->GetRasterDataType()));
    2094             : 
    2095           1 :                     const char *pszDescription = poBand->GetDescription();
    2096           1 :                     if (pszDescription && pszDescription[0])
    2097           0 :                         CPLSetXMLValue(psBand, GTI_XML_BAND_DESCRIPTION,
    2098             :                                        pszDescription);
    2099             : 
    2100           1 :                     const auto eColorInterp = poBand->GetColorInterpretation();
    2101           1 :                     if (eColorInterp != GCI_Undefined)
    2102           1 :                         CPLSetXMLValue(
    2103             :                             psBand, GTI_XML_BAND_COLORINTERP,
    2104             :                             GDALGetColorInterpretationName(eColorInterp));
    2105             : 
    2106           1 :                     if (!std::isnan(poBand->m_dfOffset))
    2107           0 :                         CPLSetXMLValue(psBand, GTI_XML_BAND_OFFSET,
    2108             :                                        CPLSPrintf("%.16g", poBand->m_dfOffset));
    2109             : 
    2110           1 :                     if (!std::isnan(poBand->m_dfScale))
    2111           0 :                         CPLSetXMLValue(psBand, GTI_XML_BAND_SCALE,
    2112             :                                        CPLSPrintf("%.16g", poBand->m_dfScale));
    2113             : 
    2114           1 :                     if (!poBand->m_osUnit.empty())
    2115           0 :                         CPLSetXMLValue(psBand, GTI_XML_BAND_UNITTYPE,
    2116             :                                        poBand->m_osUnit.c_str());
    2117             : 
    2118           1 :                     if (poBand->m_bNoDataValueSet)
    2119             :                     {
    2120           0 :                         CPLSetXMLValue(
    2121             :                             psBand, GTI_XML_BAND_NODATAVALUE,
    2122           0 :                             VRTSerializeNoData(poBand->m_dfNoDataValue,
    2123             :                                                poBand->GetRasterDataType(), 18)
    2124             :                                 .c_str());
    2125             :                     }
    2126           1 :                     if (aoBandXML[i - 1])
    2127             :                     {
    2128           1 :                         CPLAddXMLChild(psBand, aoBandXML[i - 1].release());
    2129             :                     }
    2130             :                 }
    2131             :             }
    2132             :         }
    2133             : 
    2134           1 :         if (!CPLSerializeXMLTreeToFile(m_psXMLTree.get(), GetDescription()))
    2135           0 :             eErr = CE_Failure;
    2136             :     }
    2137             : 
    2138             :     // We also clear the cache of opened sources, in case the user would
    2139             :     // change the content of a source and would want the GTI dataset to see
    2140             :     // the refreshed content.
    2141         186 :     m_oMapSharedSources.clear();
    2142         186 :     m_dfLastMinXFilter = std::numeric_limits<double>::quiet_NaN();
    2143         186 :     m_dfLastMinYFilter = std::numeric_limits<double>::quiet_NaN();
    2144         186 :     m_dfLastMaxXFilter = std::numeric_limits<double>::quiet_NaN();
    2145         186 :     m_dfLastMaxYFilter = std::numeric_limits<double>::quiet_NaN();
    2146         186 :     m_aoSourceDesc.clear();
    2147         186 :     if (GDALPamDataset::FlushCache(bAtClosing) != CE_None)
    2148           0 :         eErr = CE_Failure;
    2149         186 :     return eErr;
    2150             : }
    2151             : 
    2152             : /************************************************************************/
    2153             : /*                            LoadOverviews()                           */
    2154             : /************************************************************************/
    2155             : 
    2156          28 : void GDALTileIndexDataset::LoadOverviews()
    2157             : {
    2158          28 :     if (m_apoOverviews.empty() && !m_aoOverviewDescriptor.empty())
    2159             :     {
    2160          22 :         for (const auto &[osDSName, aosOpenOptions, osLyrName, dfFactor] :
    2161          33 :              m_aoOverviewDescriptor)
    2162             :         {
    2163          22 :             CPLStringList aosNewOpenOptions(aosOpenOptions);
    2164          11 :             if (dfFactor != 0)
    2165             :             {
    2166             :                 aosNewOpenOptions.SetNameValue("@FACTOR",
    2167           4 :                                                CPLSPrintf("%.18g", dfFactor));
    2168             :             }
    2169          11 :             if (!osLyrName.empty())
    2170             :             {
    2171           4 :                 aosNewOpenOptions.SetNameValue("@LAYER", osLyrName.c_str());
    2172             :             }
    2173             : 
    2174             :             std::unique_ptr<GDALDataset> poOvrDS(GDALDataset::Open(
    2175          22 :                 !osDSName.empty() ? osDSName.c_str() : GetDescription(),
    2176             :                 GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, nullptr,
    2177          44 :                 aosNewOpenOptions.List(), nullptr));
    2178          11 :             if (poOvrDS)
    2179             :             {
    2180           6 :                 if (poOvrDS->GetRasterCount() == GetRasterCount())
    2181             :                 {
    2182           6 :                     m_apoOverviews.emplace_back(std::move(poOvrDS));
    2183             :                 }
    2184             :                 else
    2185             :                 {
    2186           0 :                     CPLError(CE_Warning, CPLE_AppDefined,
    2187             :                              "%s has not the same number of bands as %s",
    2188           0 :                              poOvrDS->GetDescription(), GetDescription());
    2189             :                 }
    2190             :             }
    2191             :         }
    2192             :     }
    2193          28 : }
    2194             : 
    2195             : /************************************************************************/
    2196             : /*                          GetOverviewCount()                          */
    2197             : /************************************************************************/
    2198             : 
    2199          40 : int GDALTileIndexBand::GetOverviewCount()
    2200             : {
    2201          40 :     const int nPAMOverviews = GDALPamRasterBand::GetOverviewCount();
    2202          40 :     if (nPAMOverviews)
    2203          12 :         return nPAMOverviews;
    2204             : 
    2205          28 :     m_poDS->LoadOverviews();
    2206          28 :     return static_cast<int>(m_poDS->m_apoOverviews.size());
    2207             : }
    2208             : 
    2209             : /************************************************************************/
    2210             : /*                             GetOverview()                            */
    2211             : /************************************************************************/
    2212             : 
    2213          23 : GDALRasterBand *GDALTileIndexBand::GetOverview(int iOvr)
    2214             : {
    2215          23 :     if (iOvr < 0 || iOvr >= GetOverviewCount())
    2216           6 :         return nullptr;
    2217             : 
    2218          17 :     const int nPAMOverviews = GDALPamRasterBand::GetOverviewCount();
    2219          17 :     if (nPAMOverviews)
    2220           7 :         return GDALPamRasterBand::GetOverview(iOvr);
    2221             : 
    2222          10 :     if (nBand == 0)
    2223             :     {
    2224           1 :         auto poBand = m_poDS->m_apoOverviews[iOvr]->GetRasterBand(1);
    2225           1 :         if (!poBand)
    2226           0 :             return nullptr;
    2227           1 :         return poBand->GetMaskBand();
    2228             :     }
    2229             :     else
    2230             :     {
    2231           9 :         return m_poDS->m_apoOverviews[iOvr]->GetRasterBand(nBand);
    2232             :     }
    2233             : }
    2234             : 
    2235             : /************************************************************************/
    2236             : /*                           GetGeoTransform()                          */
    2237             : /************************************************************************/
    2238             : 
    2239          16 : CPLErr GDALTileIndexDataset::GetGeoTransform(double *padfGeoTransform)
    2240             : {
    2241          16 :     memcpy(padfGeoTransform, m_adfGeoTransform.data(), 6 * sizeof(double));
    2242          16 :     return CE_None;
    2243             : }
    2244             : 
    2245             : /************************************************************************/
    2246             : /*                            GetSpatialRef()                           */
    2247             : /************************************************************************/
    2248             : 
    2249          12 : const OGRSpatialReference *GDALTileIndexDataset::GetSpatialRef() const
    2250             : {
    2251          12 :     return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
    2252             : }
    2253             : 
    2254             : /************************************************************************/
    2255             : /*                           GDALTileIndexBand()                         */
    2256             : /************************************************************************/
    2257             : 
    2258         154 : GDALTileIndexBand::GDALTileIndexBand(GDALTileIndexDataset *poDSIn, int nBandIn,
    2259             :                                      GDALDataType eDT, int nBlockXSizeIn,
    2260         154 :                                      int nBlockYSizeIn)
    2261             : {
    2262         154 :     m_poDS = poDSIn;
    2263         154 :     nBand = nBandIn;
    2264         154 :     eDataType = eDT;
    2265         154 :     nRasterXSize = poDSIn->GetRasterXSize();
    2266         154 :     nRasterYSize = poDSIn->GetRasterYSize();
    2267         154 :     nBlockXSize = nBlockXSizeIn;
    2268         154 :     nBlockYSize = nBlockYSizeIn;
    2269         154 : }
    2270             : 
    2271             : /************************************************************************/
    2272             : /*                             IReadBlock()                             */
    2273             : /************************************************************************/
    2274             : 
    2275           3 : CPLErr GDALTileIndexBand::IReadBlock(int nBlockXOff, int nBlockYOff,
    2276             :                                      void *pImage)
    2277             : 
    2278             : {
    2279           3 :     const int nPixelSize = GDALGetDataTypeSizeBytes(eDataType);
    2280             : 
    2281           3 :     int nReadXSize = nBlockXSize;
    2282           3 :     int nReadYSize = nBlockYSize;
    2283           3 :     GetActualBlockSize(nBlockXOff, nBlockYOff, &nReadXSize, &nReadYSize);
    2284             : 
    2285             :     GDALRasterIOExtraArg sExtraArg;
    2286           3 :     INIT_RASTERIO_EXTRA_ARG(sExtraArg);
    2287             : 
    2288           6 :     return IRasterIO(
    2289           3 :         GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize, nReadXSize,
    2290             :         nReadYSize, pImage, nReadXSize, nReadYSize, eDataType, nPixelSize,
    2291           6 :         static_cast<GSpacing>(nPixelSize) * nBlockXSize, &sExtraArg);
    2292             : }
    2293             : 
    2294             : /************************************************************************/
    2295             : /*                             IRasterIO()                              */
    2296             : /************************************************************************/
    2297             : 
    2298         119 : CPLErr GDALTileIndexBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
    2299             :                                     int nXSize, int nYSize, void *pData,
    2300             :                                     int nBufXSize, int nBufYSize,
    2301             :                                     GDALDataType eBufType, GSpacing nPixelSpace,
    2302             :                                     GSpacing nLineSpace,
    2303             :                                     GDALRasterIOExtraArg *psExtraArg)
    2304             : {
    2305         119 :     int anBand[] = {nBand};
    2306             : 
    2307         119 :     return m_poDS->IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
    2308             :                              nBufXSize, nBufYSize, eBufType, 1, anBand,
    2309         238 :                              nPixelSpace, nLineSpace, 0, psExtraArg);
    2310             : }
    2311             : 
    2312             : /************************************************************************/
    2313             : /*                      GetMetadataDomainList()                         */
    2314             : /************************************************************************/
    2315             : 
    2316           1 : char **GDALTileIndexBand::GetMetadataDomainList()
    2317             : {
    2318           1 :     return CSLAddString(GDALRasterBand::GetMetadataDomainList(),
    2319           1 :                         "LocationInfo");
    2320             : }
    2321             : 
    2322             : /************************************************************************/
    2323             : /*                          GetMetadataItem()                           */
    2324             : /************************************************************************/
    2325             : 
    2326          27 : const char *GDALTileIndexBand::GetMetadataItem(const char *pszName,
    2327             :                                                const char *pszDomain)
    2328             : 
    2329             : {
    2330             :     /* ==================================================================== */
    2331             :     /*      LocationInfo handling.                                          */
    2332             :     /* ==================================================================== */
    2333          27 :     if (pszDomain != nullptr && EQUAL(pszDomain, "LocationInfo") &&
    2334          19 :         (STARTS_WITH_CI(pszName, "Pixel_") ||
    2335           6 :          STARTS_WITH_CI(pszName, "GeoPixel_")))
    2336             :     {
    2337             :         // What pixel are we aiming at?
    2338          18 :         int iPixel = 0;
    2339          18 :         int iLine = 0;
    2340             : 
    2341          18 :         if (STARTS_WITH_CI(pszName, "Pixel_"))
    2342             :         {
    2343          13 :             pszName += strlen("Pixel_");
    2344          13 :             iPixel = atoi(pszName);
    2345          13 :             const char *const pszUnderscore = strchr(pszName, '_');
    2346          13 :             if (!pszUnderscore)
    2347           2 :                 return nullptr;
    2348          11 :             iLine = atoi(pszUnderscore + 1);
    2349             :         }
    2350           5 :         else if (STARTS_WITH_CI(pszName, "GeoPixel_"))
    2351             :         {
    2352           5 :             pszName += strlen("GeoPixel_");
    2353           5 :             const double dfGeoX = CPLAtof(pszName);
    2354           5 :             const char *const pszUnderscore = strchr(pszName, '_');
    2355           5 :             if (!pszUnderscore)
    2356           2 :                 return nullptr;
    2357           3 :             const double dfGeoY = CPLAtof(pszUnderscore + 1);
    2358             : 
    2359           3 :             double adfInvGeoTransform[6] = {0.0};
    2360           3 :             if (!GDALInvGeoTransform(m_poDS->m_adfGeoTransform.data(),
    2361             :                                      adfInvGeoTransform))
    2362           0 :                 return nullptr;
    2363             : 
    2364           3 :             iPixel = static_cast<int>(floor(adfInvGeoTransform[0] +
    2365           3 :                                             adfInvGeoTransform[1] * dfGeoX +
    2366           3 :                                             adfInvGeoTransform[2] * dfGeoY));
    2367           3 :             iLine = static_cast<int>(floor(adfInvGeoTransform[3] +
    2368           3 :                                            adfInvGeoTransform[4] * dfGeoX +
    2369           3 :                                            adfInvGeoTransform[5] * dfGeoY));
    2370             :         }
    2371             :         else
    2372             :         {
    2373           0 :             return nullptr;
    2374             :         }
    2375             : 
    2376          23 :         if (iPixel < 0 || iLine < 0 || iPixel >= GetXSize() ||
    2377           9 :             iLine >= GetYSize())
    2378           6 :             return nullptr;
    2379             : 
    2380           8 :         if (!m_poDS->CollectSources(iPixel, iLine, 1, 1))
    2381           0 :             return nullptr;
    2382             : 
    2383             :         // Format into XML.
    2384           8 :         m_osLastLocationInfo = "<LocationInfo>";
    2385             : 
    2386           8 :         if (!m_poDS->m_aoSourceDesc.empty())
    2387             :         {
    2388             :             const auto AddSource =
    2389           6 :                 [&](const GDALTileIndexDataset::SourceDesc &oSourceDesc)
    2390             :             {
    2391           6 :                 m_osLastLocationInfo += "<File>";
    2392             :                 char *const pszXMLEscaped =
    2393           6 :                     CPLEscapeString(oSourceDesc.osName.c_str(), -1, CPLES_XML);
    2394           6 :                 m_osLastLocationInfo += pszXMLEscaped;
    2395           6 :                 CPLFree(pszXMLEscaped);
    2396           6 :                 m_osLastLocationInfo += "</File>";
    2397          11 :             };
    2398             : 
    2399           5 :             const int anBand[] = {nBand};
    2400           5 :             if (!m_poDS->NeedInitBuffer(1, anBand))
    2401             :             {
    2402           4 :                 AddSource(m_poDS->m_aoSourceDesc.back());
    2403             :             }
    2404             :             else
    2405             :             {
    2406           3 :                 for (const auto &oSourceDesc : m_poDS->m_aoSourceDesc)
    2407             :                 {
    2408           2 :                     if (oSourceDesc.poDS)
    2409           2 :                         AddSource(oSourceDesc);
    2410             :                 }
    2411             :             }
    2412             :         }
    2413             : 
    2414           8 :         m_osLastLocationInfo += "</LocationInfo>";
    2415             : 
    2416           8 :         return m_osLastLocationInfo.c_str();
    2417             :     }
    2418             : 
    2419           9 :     return GDALPamRasterBand::GetMetadataItem(pszName, pszDomain);
    2420             : }
    2421             : 
    2422             : /************************************************************************/
    2423             : /*                        SetMetadataItem()                             */
    2424             : /************************************************************************/
    2425             : 
    2426          13 : CPLErr GDALTileIndexBand::SetMetadataItem(const char *pszName,
    2427             :                                           const char *pszValue,
    2428             :                                           const char *pszDomain)
    2429             : {
    2430          13 :     if (nBand > 0 && m_poDS->m_bXMLUpdatable)
    2431             :     {
    2432           1 :         m_poDS->m_bXMLModified = true;
    2433           1 :         return GDALRasterBand::SetMetadataItem(pszName, pszValue, pszDomain);
    2434             :     }
    2435          12 :     else if (nBand > 0 && m_poDS->TileIndexSupportsEditingLayerMetadata())
    2436             :     {
    2437           6 :         m_poDS->m_poLayer->SetMetadataItem(
    2438           6 :             CPLSPrintf("BAND_%d_%s", nBand, pszName), pszValue, pszDomain);
    2439           6 :         return GDALRasterBand::SetMetadataItem(pszName, pszValue, pszDomain);
    2440             :     }
    2441             :     else
    2442             :     {
    2443           6 :         return GDALPamRasterBand::SetMetadataItem(pszName, pszValue, pszDomain);
    2444             :     }
    2445             : }
    2446             : 
    2447             : /************************************************************************/
    2448             : /*                           SetMetadata()                              */
    2449             : /************************************************************************/
    2450             : 
    2451           2 : CPLErr GDALTileIndexBand::SetMetadata(char **papszMD, const char *pszDomain)
    2452             : {
    2453           2 :     if (nBand > 0 && m_poDS->m_bXMLUpdatable)
    2454             :     {
    2455           1 :         m_poDS->m_bXMLModified = true;
    2456           1 :         return GDALRasterBand::SetMetadata(papszMD, pszDomain);
    2457             :     }
    2458           1 :     else if (nBand > 0 && m_poDS->TileIndexSupportsEditingLayerMetadata())
    2459             :     {
    2460           2 :         CPLStringList aosMD;
    2461             : 
    2462           1 :         if (!pszDomain || pszDomain[0] == 0)
    2463             :         {
    2464             :             // Reinject dataset metadata
    2465           1 :             char **papszLayerMD = m_poDS->m_poLayer->GetMetadata(pszDomain);
    2466          14 :             for (const char *const *papszIter = papszLayerMD;
    2467          14 :                  papszIter && *papszIter; ++papszIter)
    2468             :             {
    2469          13 :                 if (!STARTS_WITH(*papszIter, "BAND_") ||
    2470          12 :                     STARTS_WITH(*papszIter, MD_BAND_COUNT))
    2471           1 :                     aosMD.AddString(*papszIter);
    2472             :             }
    2473             :         }
    2474             : 
    2475           8 :         for (int i = 0; papszMD && papszMD[i]; ++i)
    2476             :         {
    2477           7 :             aosMD.AddString(CPLSPrintf("BAND_%d_%s", nBand, papszMD[i]));
    2478             :         }
    2479             : 
    2480           1 :         if (!pszDomain || pszDomain[0] == 0)
    2481             :         {
    2482           4 :             for (const char *pszItem : apszReservedBandItems)
    2483             :             {
    2484           3 :                 const char *pszKey = CPLSPrintf("BAND_%d_%s", nBand, pszItem);
    2485           3 :                 if (!aosMD.FetchNameValue(pszKey))
    2486             :                 {
    2487           3 :                     if (const char *pszVal =
    2488           3 :                             m_poDS->m_poLayer->GetMetadataItem(pszKey))
    2489             :                     {
    2490           3 :                         aosMD.SetNameValue(pszKey, pszVal);
    2491             :                     }
    2492             :                 }
    2493             :             }
    2494             :         }
    2495             : 
    2496           1 :         m_poDS->m_poLayer->SetMetadata(aosMD.List(), pszDomain);
    2497           1 :         return GDALRasterBand::SetMetadata(papszMD, pszDomain);
    2498             :     }
    2499             :     else
    2500             :     {
    2501           0 :         return GDALPamRasterBand::SetMetadata(papszMD, pszDomain);
    2502             :     }
    2503             : }
    2504             : 
    2505             : /************************************************************************/
    2506             : /*                         GetSrcDstWin()                               */
    2507             : /************************************************************************/
    2508             : 
    2509         135 : static bool GetSrcDstWin(const double adfTileGT[6], int nTileXSize,
    2510             :                          int nTileYSize, const double adfVRTGT[6],
    2511             :                          int nVRTXSize, int nVRTYSize, double *pdfSrcXOff,
    2512             :                          double *pdfSrcYOff, double *pdfSrcXSize,
    2513             :                          double *pdfSrcYSize, double *pdfDstXOff,
    2514             :                          double *pdfDstYOff, double *pdfDstXSize,
    2515             :                          double *pdfDstYSize)
    2516             : {
    2517         135 :     const double minX = adfVRTGT[GT_TOPLEFT_X];
    2518         135 :     const double we_res = adfVRTGT[GT_WE_RES];
    2519         135 :     const double maxX = minX + nVRTXSize * we_res;
    2520         135 :     const double maxY = adfVRTGT[GT_TOPLEFT_Y];
    2521         135 :     const double ns_res = adfVRTGT[GT_NS_RES];
    2522         135 :     const double minY = maxY + nVRTYSize * ns_res;
    2523             : 
    2524             :     /* Check that the destination bounding box intersects the source bounding
    2525             :      * box */
    2526         135 :     if (adfTileGT[GT_TOPLEFT_X] + nTileXSize * adfTileGT[GT_WE_RES] <= minX)
    2527           0 :         return false;
    2528         135 :     if (adfTileGT[GT_TOPLEFT_X] >= maxX)
    2529           0 :         return false;
    2530         135 :     if (adfTileGT[GT_TOPLEFT_Y] + nTileYSize * adfTileGT[GT_NS_RES] >= maxY)
    2531           0 :         return false;
    2532         135 :     if (adfTileGT[GT_TOPLEFT_Y] <= minY)
    2533           0 :         return false;
    2534             : 
    2535         135 :     if (adfTileGT[GT_TOPLEFT_X] < minX)
    2536             :     {
    2537           1 :         *pdfSrcXOff = (minX - adfTileGT[GT_TOPLEFT_X]) / adfTileGT[GT_WE_RES];
    2538           1 :         *pdfDstXOff = 0.0;
    2539             :     }
    2540             :     else
    2541             :     {
    2542         134 :         *pdfSrcXOff = 0.0;
    2543         134 :         *pdfDstXOff = ((adfTileGT[GT_TOPLEFT_X] - minX) / we_res);
    2544             :     }
    2545         135 :     if (maxY < adfTileGT[GT_TOPLEFT_Y])
    2546             :     {
    2547           1 :         *pdfSrcYOff = (adfTileGT[GT_TOPLEFT_Y] - maxY) / -adfTileGT[GT_NS_RES];
    2548           1 :         *pdfDstYOff = 0.0;
    2549             :     }
    2550             :     else
    2551             :     {
    2552         134 :         *pdfSrcYOff = 0.0;
    2553         134 :         *pdfDstYOff = ((maxY - adfTileGT[GT_TOPLEFT_Y]) / -ns_res);
    2554             :     }
    2555             : 
    2556         135 :     *pdfSrcXSize = nTileXSize;
    2557         135 :     *pdfSrcYSize = nTileYSize;
    2558         135 :     if (*pdfSrcXOff > 0)
    2559           1 :         *pdfSrcXSize -= *pdfSrcXOff;
    2560         135 :     if (*pdfSrcYOff > 0)
    2561           1 :         *pdfSrcYSize -= *pdfSrcYOff;
    2562             : 
    2563         135 :     const double dfSrcToDstXSize = adfTileGT[GT_WE_RES] / we_res;
    2564         135 :     *pdfDstXSize = *pdfSrcXSize * dfSrcToDstXSize;
    2565         135 :     const double dfSrcToDstYSize = adfTileGT[GT_NS_RES] / ns_res;
    2566         135 :     *pdfDstYSize = *pdfSrcYSize * dfSrcToDstYSize;
    2567             : 
    2568         135 :     if (*pdfDstXOff + *pdfDstXSize > nVRTXSize)
    2569             :     {
    2570           2 :         *pdfDstXSize = nVRTXSize - *pdfDstXOff;
    2571           2 :         *pdfSrcXSize = *pdfDstXSize / dfSrcToDstXSize;
    2572             :     }
    2573             : 
    2574         135 :     if (*pdfDstYOff + *pdfDstYSize > nVRTYSize)
    2575             :     {
    2576           1 :         *pdfDstYSize = nVRTYSize - *pdfDstYOff;
    2577           1 :         *pdfSrcYSize = *pdfDstYSize / dfSrcToDstYSize;
    2578             :     }
    2579             : 
    2580         270 :     return *pdfSrcXSize > 0 && *pdfDstXSize > 0 && *pdfSrcYSize > 0 &&
    2581         270 :            *pdfDstYSize > 0;
    2582             : }
    2583             : 
    2584             : /************************************************************************/
    2585             : /*                   GDALDatasetCastToGTIDataset()                    */
    2586             : /************************************************************************/
    2587             : 
    2588           1 : GDALTileIndexDataset *GDALDatasetCastToGTIDataset(GDALDataset *poDS)
    2589             : {
    2590           1 :     return dynamic_cast<GDALTileIndexDataset *>(poDS);
    2591             : }
    2592             : 
    2593             : /************************************************************************/
    2594             : /*                   GTIGetSourcesMoreRecentThan()                    */
    2595             : /************************************************************************/
    2596             : 
    2597             : std::vector<GTISourceDesc>
    2598           1 : GTIGetSourcesMoreRecentThan(GDALTileIndexDataset *poDS, int64_t mTime)
    2599             : {
    2600           1 :     return poDS->GetSourcesMoreRecentThan(mTime);
    2601             : }
    2602             : 
    2603             : /************************************************************************/
    2604             : /*                       GetSourcesMoreRecentThan()                     */
    2605             : /************************************************************************/
    2606             : 
    2607             : std::vector<GTISourceDesc>
    2608           1 : GDALTileIndexDataset::GetSourcesMoreRecentThan(int64_t mTime)
    2609             : {
    2610           1 :     std::vector<GTISourceDesc> oRes;
    2611             : 
    2612           1 :     m_poLayer->SetSpatialFilter(nullptr);
    2613           3 :     for (auto &&poFeature : m_poLayer)
    2614             :     {
    2615           2 :         if (!poFeature->IsFieldSetAndNotNull(m_nLocationFieldIndex))
    2616             :         {
    2617           1 :             continue;
    2618             :         }
    2619             : 
    2620           2 :         auto poGeom = poFeature->GetGeometryRef();
    2621           2 :         if (!poGeom || poGeom->IsEmpty())
    2622           0 :             continue;
    2623             : 
    2624           2 :         OGREnvelope sEnvelope;
    2625           2 :         poGeom->getEnvelope(&sEnvelope);
    2626             : 
    2627           2 :         double dfXOff = (sEnvelope.MinX - m_adfGeoTransform[GT_TOPLEFT_X]) /
    2628           2 :                         m_adfGeoTransform[GT_WE_RES];
    2629           2 :         if (dfXOff >= nRasterXSize)
    2630           0 :             continue;
    2631             : 
    2632           2 :         double dfYOff = (sEnvelope.MaxY - m_adfGeoTransform[GT_TOPLEFT_Y]) /
    2633           2 :                         m_adfGeoTransform[GT_NS_RES];
    2634           2 :         if (dfYOff >= nRasterYSize)
    2635           0 :             continue;
    2636             : 
    2637             :         double dfXSize =
    2638           2 :             (sEnvelope.MaxX - sEnvelope.MinX) / m_adfGeoTransform[GT_WE_RES];
    2639           2 :         if (dfXOff < 0)
    2640             :         {
    2641           0 :             dfXSize += dfXOff;
    2642           0 :             dfXOff = 0;
    2643           0 :             if (dfXSize <= 0)
    2644           0 :                 continue;
    2645             :         }
    2646             : 
    2647           2 :         double dfYSize = (sEnvelope.MaxY - sEnvelope.MinY) /
    2648           2 :                          std::fabs(m_adfGeoTransform[GT_NS_RES]);
    2649           2 :         if (dfYOff < 0)
    2650             :         {
    2651           0 :             dfYSize += dfYOff;
    2652           0 :             dfYOff = 0;
    2653           0 :             if (dfYSize <= 0)
    2654           0 :                 continue;
    2655             :         }
    2656             : 
    2657             :         const char *pszTileName =
    2658           2 :             poFeature->GetFieldAsString(m_nLocationFieldIndex);
    2659             :         const std::string osTileName(
    2660           2 :             GetAbsoluteFileName(pszTileName, GetDescription()));
    2661             :         VSIStatBufL sStatSource;
    2662           4 :         if (VSIStatL(osTileName.c_str(), &sStatSource) != 0 ||
    2663           2 :             sStatSource.st_mtime <= mTime)
    2664             :         {
    2665           1 :             continue;
    2666             :         }
    2667             : 
    2668           1 :         constexpr double EPS = 1e-8;
    2669           2 :         GTISourceDesc oSourceDesc;
    2670           1 :         oSourceDesc.osFilename = osTileName;
    2671           1 :         oSourceDesc.nDstXOff = static_cast<int>(dfXOff + EPS);
    2672           1 :         oSourceDesc.nDstYOff = static_cast<int>(dfYOff + EPS);
    2673           1 :         oSourceDesc.nDstXSize = static_cast<int>(dfXSize + 0.5);
    2674           1 :         oSourceDesc.nDstYSize = static_cast<int>(dfYSize + 0.5);
    2675           1 :         oRes.emplace_back(std::move(oSourceDesc));
    2676             :     }
    2677             : 
    2678           1 :     return oRes;
    2679             : }
    2680             : 
    2681             : /************************************************************************/
    2682             : /*                         GetSourceDesc()                              */
    2683             : /************************************************************************/
    2684             : 
    2685         136 : bool GDALTileIndexDataset::GetSourceDesc(const std::string &osTileName,
    2686             :                                          SourceDesc &oSourceDesc)
    2687             : {
    2688         136 :     std::shared_ptr<GDALDataset> poTileDS;
    2689         136 :     if (!m_oMapSharedSources.tryGet(osTileName, poTileDS))
    2690             :     {
    2691         158 :         poTileDS = std::shared_ptr<GDALDataset>(
    2692             :             GDALProxyPoolDataset::Create(
    2693             :                 osTileName.c_str(), nullptr, GA_ReadOnly,
    2694             :                 /* bShared = */ true, m_osUniqueHandle.c_str()),
    2695          79 :             GDALDatasetUniquePtrReleaser());
    2696          79 :         if (!poTileDS || poTileDS->GetRasterCount() == 0)
    2697             :         {
    2698           1 :             return false;
    2699             :         }
    2700             : 
    2701             :         // do palette -> RGB(A) expansion if needed
    2702          78 :         if (!GTIDoPaletteExpansionIfNeeded(poTileDS, nBands))
    2703           0 :             return false;
    2704             : 
    2705             :         const OGRSpatialReference *poTileSRS;
    2706          78 :         if (!m_oSRS.IsEmpty() &&
    2707         104 :             (poTileSRS = poTileDS->GetSpatialRef()) != nullptr &&
    2708          26 :             !m_oSRS.IsSame(poTileSRS))
    2709             :         {
    2710           2 :             CPLDebug("VRT",
    2711             :                      "Tile %s has not the same SRS as the VRT. "
    2712             :                      "Proceed to on-the-fly warping",
    2713             :                      osTileName.c_str());
    2714             : 
    2715           2 :             CPLStringList aosOptions;
    2716           2 :             aosOptions.AddString("-of");
    2717           2 :             aosOptions.AddString("VRT");
    2718             : 
    2719           2 :             if ((poTileDS->GetRasterBand(1)->GetColorTable() == nullptr &&
    2720           2 :                  poTileDS->GetRasterBand(1)->GetCategoryNames() == nullptr) ||
    2721           0 :                 m_eResampling == GRIORA_Mode)
    2722             :             {
    2723           2 :                 aosOptions.AddString("-r");
    2724           2 :                 aosOptions.AddString(m_osResampling.c_str());
    2725             :             }
    2726             : 
    2727           2 :             if (m_osWKT.empty())
    2728             :             {
    2729           0 :                 char *pszWKT = nullptr;
    2730           0 :                 const char *const apszWKTOptions[] = {"FORMAT=WKT2_2019",
    2731             :                                                       nullptr};
    2732           0 :                 m_oSRS.exportToWkt(&pszWKT, apszWKTOptions);
    2733           0 :                 if (pszWKT)
    2734           0 :                     m_osWKT = pszWKT;
    2735           0 :                 CPLFree(pszWKT);
    2736             :             }
    2737           2 :             if (m_osWKT.empty())
    2738             :             {
    2739           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
    2740             :                          "Cannot export VRT SRS to WKT2");
    2741           0 :                 return false;
    2742             :             }
    2743             : 
    2744           2 :             aosOptions.AddString("-t_srs");
    2745           2 :             aosOptions.AddString(m_osWKT.c_str());
    2746             : 
    2747             :             // First pass to get the extent of the tile in the
    2748             :             // target VRT SRS
    2749             :             GDALWarpAppOptions *psWarpOptions =
    2750           2 :                 GDALWarpAppOptionsNew(aosOptions.List(), nullptr);
    2751           2 :             GDALDatasetH ahSrcDS[] = {GDALDataset::ToHandle(poTileDS.get())};
    2752           2 :             int bUsageError = false;
    2753             :             auto poWarpDS =
    2754             :                 std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(GDALWarp(
    2755           2 :                     "", nullptr, 1, ahSrcDS, psWarpOptions, &bUsageError)));
    2756           2 :             GDALWarpAppOptionsFree(psWarpOptions);
    2757           2 :             if (!poWarpDS)
    2758             :             {
    2759           0 :                 return false;
    2760             :             }
    2761             : 
    2762             :             // Second pass to create a warped source VRT whose
    2763             :             // extent is aligned on the one of the target VRT
    2764             :             double adfWarpDSGeoTransform[6];
    2765           2 :             const auto eErr = poWarpDS->GetGeoTransform(adfWarpDSGeoTransform);
    2766           2 :             CPL_IGNORE_RET_VAL(eErr);
    2767           2 :             CPLAssert(eErr == CE_None);
    2768           2 :             const double dfVRTMinX = m_adfGeoTransform[GT_TOPLEFT_X];
    2769           2 :             const double dfVRTResX = m_adfGeoTransform[GT_WE_RES];
    2770           2 :             const double dfVRTMaxY = m_adfGeoTransform[GT_TOPLEFT_Y];
    2771           2 :             const double dfVRTResYAbs = -m_adfGeoTransform[GT_NS_RES];
    2772           2 :             const double dfWarpMinX =
    2773           2 :                 std::floor((adfWarpDSGeoTransform[GT_TOPLEFT_X] - dfVRTMinX) /
    2774           2 :                            dfVRTResX) *
    2775             :                     dfVRTResX +
    2776             :                 dfVRTMinX;
    2777             :             const double dfWarpMaxX =
    2778           4 :                 std::ceil((adfWarpDSGeoTransform[GT_TOPLEFT_X] +
    2779           4 :                            adfWarpDSGeoTransform[GT_WE_RES] *
    2780           2 :                                poWarpDS->GetRasterXSize() -
    2781             :                            dfVRTMinX) /
    2782           2 :                           dfVRTResX) *
    2783             :                     dfVRTResX +
    2784           2 :                 dfVRTMinX;
    2785           2 :             const double dfWarpMaxY =
    2786             :                 dfVRTMaxY -
    2787           2 :                 std::floor((dfVRTMaxY - adfWarpDSGeoTransform[GT_TOPLEFT_Y]) /
    2788           2 :                            dfVRTResYAbs) *
    2789             :                     dfVRTResYAbs;
    2790             :             const double dfWarpMinY =
    2791             :                 dfVRTMaxY -
    2792           4 :                 std::ceil((dfVRTMaxY - (adfWarpDSGeoTransform[GT_TOPLEFT_Y] +
    2793           4 :                                         adfWarpDSGeoTransform[GT_NS_RES] *
    2794           2 :                                             poWarpDS->GetRasterYSize())) /
    2795           2 :                           dfVRTResYAbs) *
    2796           2 :                     dfVRTResYAbs;
    2797             : 
    2798           2 :             aosOptions.AddString("-te");
    2799           2 :             aosOptions.AddString(CPLSPrintf("%.18g", dfWarpMinX));
    2800           2 :             aosOptions.AddString(CPLSPrintf("%.18g", dfWarpMinY));
    2801           2 :             aosOptions.AddString(CPLSPrintf("%.18g", dfWarpMaxX));
    2802           2 :             aosOptions.AddString(CPLSPrintf("%.18g", dfWarpMaxY));
    2803             : 
    2804           2 :             aosOptions.AddString("-tr");
    2805           2 :             aosOptions.AddString(CPLSPrintf("%.18g", dfVRTResX));
    2806           2 :             aosOptions.AddString(CPLSPrintf("%.18g", dfVRTResYAbs));
    2807             : 
    2808           2 :             aosOptions.AddString("-dstalpha");
    2809             : 
    2810           2 :             psWarpOptions = GDALWarpAppOptionsNew(aosOptions.List(), nullptr);
    2811           2 :             poWarpDS.reset(GDALDataset::FromHandle(GDALWarp(
    2812             :                 "", nullptr, 1, ahSrcDS, psWarpOptions, &bUsageError)));
    2813           2 :             GDALWarpAppOptionsFree(psWarpOptions);
    2814           2 :             if (!poWarpDS)
    2815             :             {
    2816           0 :                 return false;
    2817             :             }
    2818             : 
    2819           2 :             poTileDS.reset(poWarpDS.release());
    2820             :         }
    2821             : 
    2822          78 :         m_oMapSharedSources.insert(osTileName, poTileDS);
    2823             :     }
    2824             : 
    2825             :     double adfGeoTransformTile[6];
    2826         135 :     if (poTileDS->GetGeoTransform(adfGeoTransformTile) != CE_None)
    2827             :     {
    2828           0 :         CPLError(CE_Failure, CPLE_AppDefined, "%s lacks geotransform",
    2829             :                  osTileName.c_str());
    2830           0 :         return false;
    2831             :     }
    2832             : 
    2833         135 :     bool bHasNoData = false;
    2834         135 :     bool bSameNoData = true;
    2835         135 :     double dfNoDataValue = 0;
    2836         135 :     GDALRasterBand *poMaskBand = nullptr;
    2837         135 :     const int nBandCount = poTileDS->GetRasterCount();
    2838         361 :     for (int iBand = 0; iBand < nBandCount; ++iBand)
    2839             :     {
    2840         226 :         auto poTileBand = poTileDS->GetRasterBand(iBand + 1);
    2841         226 :         int bThisBandHasNoData = false;
    2842             :         const double dfThisBandNoDataValue =
    2843         226 :             poTileBand->GetNoDataValue(&bThisBandHasNoData);
    2844         226 :         if (bThisBandHasNoData)
    2845             :         {
    2846          10 :             bHasNoData = true;
    2847          10 :             dfNoDataValue = dfThisBandNoDataValue;
    2848             :         }
    2849         317 :         if (iBand > 0 &&
    2850          91 :             (static_cast<int>(bThisBandHasNoData) !=
    2851          91 :                  static_cast<int>(bHasNoData) ||
    2852           4 :              (bHasNoData &&
    2853           4 :               !IsSameNaNAware(dfNoDataValue, dfThisBandNoDataValue))))
    2854             :         {
    2855           0 :             bSameNoData = false;
    2856             :         }
    2857             : 
    2858         226 :         if (poTileBand->GetMaskFlags() == GMF_PER_DATASET)
    2859           2 :             poMaskBand = poTileBand->GetMaskBand();
    2860         224 :         else if (poTileBand->GetColorInterpretation() == GCI_AlphaBand)
    2861          31 :             poMaskBand = poTileBand;
    2862             :     }
    2863             : 
    2864         135 :     std::unique_ptr<VRTSimpleSource> poSource;
    2865         135 :     if (!bHasNoData)
    2866             :     {
    2867         129 :         poSource = std::make_unique<VRTSimpleSource>();
    2868             :     }
    2869             :     else
    2870             :     {
    2871          12 :         auto poComplexSource = std::make_unique<VRTComplexSource>();
    2872           6 :         poComplexSource->SetNoDataValue(dfNoDataValue);
    2873           6 :         poSource = std::move(poComplexSource);
    2874             :     }
    2875             : 
    2876         270 :     if (!GetSrcDstWin(adfGeoTransformTile, poTileDS->GetRasterXSize(),
    2877         135 :                       poTileDS->GetRasterYSize(), m_adfGeoTransform.data(),
    2878             :                       GetRasterXSize(), GetRasterYSize(),
    2879         135 :                       &poSource->m_dfSrcXOff, &poSource->m_dfSrcYOff,
    2880         135 :                       &poSource->m_dfSrcXSize, &poSource->m_dfSrcYSize,
    2881         135 :                       &poSource->m_dfDstXOff, &poSource->m_dfDstYOff,
    2882         135 :                       &poSource->m_dfDstXSize, &poSource->m_dfDstYSize))
    2883             :     {
    2884             :         // Should not happen on a consistent tile index
    2885           0 :         CPLDebug("VRT", "Tile %s does not actually intersect area of interest",
    2886             :                  osTileName.c_str());
    2887           0 :         return false;
    2888             :     }
    2889             : 
    2890         135 :     oSourceDesc.osName = osTileName;
    2891         135 :     oSourceDesc.poDS = std::move(poTileDS);
    2892         135 :     oSourceDesc.poSource = std::move(poSource);
    2893         135 :     oSourceDesc.bHasNoData = bHasNoData;
    2894         135 :     oSourceDesc.bSameNoData = bSameNoData;
    2895         135 :     if (bSameNoData)
    2896         135 :         oSourceDesc.dfSameNoData = dfNoDataValue;
    2897         135 :     oSourceDesc.poMaskBand = poMaskBand;
    2898         135 :     return true;
    2899             : }
    2900             : 
    2901             : /************************************************************************/
    2902             : /*                        CollectSources()                              */
    2903             : /************************************************************************/
    2904             : 
    2905         154 : bool GDALTileIndexDataset::CollectSources(double dfXOff, double dfYOff,
    2906             :                                           double dfXSize, double dfYSize)
    2907             : {
    2908             :     const double dfMinX =
    2909         154 :         m_adfGeoTransform[GT_TOPLEFT_X] + dfXOff * m_adfGeoTransform[GT_WE_RES];
    2910         154 :     const double dfMaxX = dfMinX + dfXSize * m_adfGeoTransform[GT_WE_RES];
    2911             :     const double dfMaxY =
    2912         154 :         m_adfGeoTransform[GT_TOPLEFT_Y] + dfYOff * m_adfGeoTransform[GT_NS_RES];
    2913         154 :     const double dfMinY = dfMaxY + dfYSize * m_adfGeoTransform[GT_NS_RES];
    2914             : 
    2915         154 :     if (dfMinX == m_dfLastMinXFilter && dfMinY == m_dfLastMinYFilter &&
    2916          43 :         dfMaxX == m_dfLastMaxXFilter && dfMaxY == m_dfLastMaxYFilter)
    2917             :     {
    2918          39 :         return true;
    2919             :     }
    2920             : 
    2921         115 :     m_dfLastMinXFilter = dfMinX;
    2922         115 :     m_dfLastMinYFilter = dfMinY;
    2923         115 :     m_dfLastMaxXFilter = dfMaxX;
    2924         115 :     m_dfLastMaxYFilter = dfMaxY;
    2925             : 
    2926         115 :     m_poLayer->SetSpatialFilterRect(dfMinX, dfMinY, dfMaxX, dfMaxY);
    2927         115 :     m_poLayer->ResetReading();
    2928             : 
    2929         115 :     m_aoSourceDesc.clear();
    2930             :     while (true)
    2931             :     {
    2932             :         auto poFeature =
    2933         284 :             std::unique_ptr<OGRFeature>(m_poLayer->GetNextFeature());
    2934         284 :         if (!poFeature)
    2935         115 :             break;
    2936         169 :         if (!poFeature->IsFieldSetAndNotNull(m_nLocationFieldIndex))
    2937             :         {
    2938           1 :             continue;
    2939             :         }
    2940             : 
    2941         168 :         SourceDesc oSourceDesc;
    2942         168 :         oSourceDesc.poFeature = std::move(poFeature);
    2943         168 :         m_aoSourceDesc.emplace_back(std::move(oSourceDesc));
    2944             : 
    2945         168 :         if (m_aoSourceDesc.size() > 10 * 1000 * 1000)
    2946             :         {
    2947             :             // Safety belt...
    2948           0 :             CPLError(CE_Failure, CPLE_AppDefined,
    2949             :                      "More than 10 million contributing sources to a "
    2950             :                      "single RasterIO() request is not supported");
    2951           0 :             return false;
    2952             :         }
    2953         169 :     }
    2954             : 
    2955         115 :     if (m_aoSourceDesc.size() > 1)
    2956             :     {
    2957          57 :         SortSourceDesc();
    2958             :     }
    2959             : 
    2960             :     // Try to find the last (most prioritary) fully opaque source covering
    2961             :     // the whole AOI. We only need to start rendering from it.
    2962         115 :     size_t i = m_aoSourceDesc.size();
    2963         170 :     while (i > 0)
    2964             :     {
    2965         136 :         --i;
    2966         136 :         auto &poFeature = m_aoSourceDesc[i].poFeature;
    2967             :         const char *pszTileName =
    2968         136 :             poFeature->GetFieldAsString(m_nLocationFieldIndex);
    2969             :         const std::string osTileName(
    2970         136 :             GetAbsoluteFileName(pszTileName, GetDescription()));
    2971             : 
    2972         136 :         SourceDesc oSourceDesc;
    2973         136 :         if (!GetSourceDesc(osTileName, oSourceDesc))
    2974           1 :             continue;
    2975             : 
    2976         135 :         const auto &poSource = oSourceDesc.poSource;
    2977         135 :         if (dfXOff >= poSource->m_dfDstXOff + poSource->m_dfDstXSize ||
    2978         135 :             dfYOff >= poSource->m_dfDstYOff + poSource->m_dfDstYSize ||
    2979         403 :             poSource->m_dfDstXOff >= dfXOff + dfXSize ||
    2980         133 :             poSource->m_dfDstYOff >= dfYOff + dfYSize)
    2981             :         {
    2982             :             // Can happen as some spatial filters select slightly more features
    2983             :             // than strictly needed.
    2984           2 :             continue;
    2985             :         }
    2986             : 
    2987             :         const bool bCoversWholeAOI =
    2988         133 :             (poSource->m_dfDstXOff <= dfXOff &&
    2989         126 :              poSource->m_dfDstYOff <= dfYOff &&
    2990         126 :              poSource->m_dfDstXOff + poSource->m_dfDstXSize >=
    2991         378 :                  dfXOff + dfXSize &&
    2992         119 :              poSource->m_dfDstYOff + poSource->m_dfDstYSize >=
    2993         119 :                  dfYOff + dfYSize);
    2994         133 :         oSourceDesc.bCoversWholeAOI = bCoversWholeAOI;
    2995             : 
    2996         133 :         m_aoSourceDesc[i] = std::move(oSourceDesc);
    2997             : 
    2998         133 :         if (m_aoSourceDesc[i].bCoversWholeAOI &&
    2999         133 :             !m_aoSourceDesc[i].bHasNoData && !m_aoSourceDesc[i].poMaskBand)
    3000             :         {
    3001          81 :             break;
    3002             :         }
    3003             :     }
    3004             : 
    3005         115 :     if (i > 0)
    3006             :     {
    3007             :         // Remove sources that will not be rendered
    3008          32 :         m_aoSourceDesc.erase(m_aoSourceDesc.begin(),
    3009          64 :                              m_aoSourceDesc.begin() + i);
    3010             :     }
    3011             : 
    3012             :     // Truncate the array when its last elements have no dataset
    3013         115 :     i = m_aoSourceDesc.size();
    3014         247 :     while (i > 0)
    3015             :     {
    3016         135 :         --i;
    3017         135 :         if (!m_aoSourceDesc[i].poDS)
    3018             :         {
    3019           3 :             m_aoSourceDesc.resize(i);
    3020           3 :             break;
    3021             :         }
    3022             :     }
    3023             : 
    3024         115 :     return true;
    3025             : }
    3026             : 
    3027             : /************************************************************************/
    3028             : /*                          SortSourceDesc()                            */
    3029             : /************************************************************************/
    3030             : 
    3031          57 : void GDALTileIndexDataset::SortSourceDesc()
    3032             : {
    3033          57 :     const auto eFieldType = m_nSortFieldIndex >= 0
    3034          57 :                                 ? m_poLayer->GetLayerDefn()
    3035          47 :                                       ->GetFieldDefn(m_nSortFieldIndex)
    3036          47 :                                       ->GetType()
    3037          57 :                                 : OFTMaxType;
    3038          57 :     std::sort(
    3039             :         m_aoSourceDesc.begin(), m_aoSourceDesc.end(),
    3040         871 :         [this, eFieldType](const SourceDesc &a, const SourceDesc &b)
    3041             :         {
    3042         100 :             const auto &poFeatureA = (m_bSortFieldAsc ? a : b).poFeature;
    3043         100 :             const auto &poFeatureB = (m_bSortFieldAsc ? b : a).poFeature;
    3044         280 :             if (m_nSortFieldIndex >= 0 &&
    3045         180 :                 poFeatureA->IsFieldSetAndNotNull(m_nSortFieldIndex) &&
    3046          80 :                 poFeatureB->IsFieldSetAndNotNull(m_nSortFieldIndex))
    3047             :             {
    3048          80 :                 if (eFieldType == OFTString)
    3049             :                 {
    3050             :                     const int nCmp =
    3051           5 :                         strcmp(poFeatureA->GetFieldAsString(m_nSortFieldIndex),
    3052             :                                poFeatureB->GetFieldAsString(m_nSortFieldIndex));
    3053           5 :                     if (nCmp < 0)
    3054           1 :                         return true;
    3055           4 :                     if (nCmp > 0)
    3056           2 :                         return false;
    3057             :                 }
    3058          75 :                 else if (eFieldType == OFTInteger || eFieldType == OFTInteger64)
    3059             :                 {
    3060             :                     const auto nA =
    3061          45 :                         poFeatureA->GetFieldAsInteger64(m_nSortFieldIndex);
    3062             :                     const auto nB =
    3063          45 :                         poFeatureB->GetFieldAsInteger64(m_nSortFieldIndex);
    3064          45 :                     if (nA < nB)
    3065           3 :                         return true;
    3066          42 :                     if (nA > nB)
    3067          42 :                         return false;
    3068             :                 }
    3069          30 :                 else if (eFieldType == OFTReal)
    3070             :                 {
    3071             :                     const auto dfA =
    3072           3 :                         poFeatureA->GetFieldAsDouble(m_nSortFieldIndex);
    3073             :                     const auto dfB =
    3074           3 :                         poFeatureB->GetFieldAsDouble(m_nSortFieldIndex);
    3075           3 :                     if (dfA < dfB)
    3076           1 :                         return true;
    3077           2 :                     if (dfA > dfB)
    3078           2 :                         return false;
    3079             :                 }
    3080          27 :                 else if (eFieldType == OFTDate || eFieldType == OFTDateTime)
    3081             :                 {
    3082             :                     const auto poFieldA =
    3083          27 :                         poFeatureA->GetRawFieldRef(m_nSortFieldIndex);
    3084             :                     const auto poFieldB =
    3085          27 :                         poFeatureB->GetRawFieldRef(m_nSortFieldIndex);
    3086             : 
    3087             : #define COMPARE_DATE_COMPONENT(comp)                                           \
    3088             :     do                                                                         \
    3089             :     {                                                                          \
    3090             :         if (poFieldA->Date.comp < poFieldB->Date.comp)                         \
    3091             :             return true;                                                       \
    3092             :         if (poFieldA->Date.comp > poFieldB->Date.comp)                         \
    3093             :             return false;                                                      \
    3094             :     } while (0)
    3095             : 
    3096          27 :                     COMPARE_DATE_COMPONENT(Year);
    3097          21 :                     COMPARE_DATE_COMPONENT(Month);
    3098          15 :                     COMPARE_DATE_COMPONENT(Day);
    3099           9 :                     COMPARE_DATE_COMPONENT(Hour);
    3100           8 :                     COMPARE_DATE_COMPONENT(Minute);
    3101           7 :                     COMPARE_DATE_COMPONENT(Second);
    3102             :                 }
    3103             :                 else
    3104             :                 {
    3105           0 :                     CPLAssert(false);
    3106             :                 }
    3107             :             }
    3108          28 :             return poFeatureA->GetFID() < poFeatureB->GetFID();
    3109             :         });
    3110          57 : }
    3111             : 
    3112             : /************************************************************************/
    3113             : /*                   CompositeSrcWithMaskIntoDest()                     */
    3114             : /************************************************************************/
    3115             : 
    3116             : static void
    3117          66 : CompositeSrcWithMaskIntoDest(const int nOutXSize, const int nOutYSize,
    3118             :                              const GDALDataType eBufType,
    3119             :                              const int nBufTypeSize, const GSpacing nPixelSpace,
    3120             :                              const GSpacing nLineSpace, const GByte *pabySrc,
    3121             :                              const GByte *const pabyMask, GByte *const pabyDest)
    3122             : {
    3123          66 :     size_t iMaskIdx = 0;
    3124          66 :     if (eBufType == GDT_Byte)
    3125             :     {
    3126             :         // Optimization for byte case
    3127         136 :         for (int iY = 0; iY < nOutYSize; iY++)
    3128             :         {
    3129          86 :             GByte *pabyDestLine =
    3130          86 :                 pabyDest + static_cast<GPtrDiff_t>(iY * nLineSpace);
    3131          86 :             int iX = 0;
    3132             : #ifdef USE_SSE2_OPTIM
    3133          86 :             if (nPixelSpace == 1)
    3134             :             {
    3135             :                 // SSE2 version up to 6 times faster than portable version
    3136          86 :                 const auto xmm_zero = _mm_setzero_si128();
    3137          86 :                 constexpr int SIZEOF_REG = static_cast<int>(sizeof(xmm_zero));
    3138         110 :                 for (; iX + SIZEOF_REG <= nOutXSize; iX += SIZEOF_REG)
    3139             :                 {
    3140          48 :                     auto xmm_mask = _mm_loadu_si128(
    3141             :                         reinterpret_cast<__m128i const *>(pabyMask + iMaskIdx));
    3142          24 :                     const auto xmm_src = _mm_loadu_si128(
    3143             :                         reinterpret_cast<__m128i const *>(pabySrc));
    3144          24 :                     auto xmm_dst = _mm_loadu_si128(
    3145             :                         reinterpret_cast<__m128i const *>(pabyDestLine));
    3146             : #ifdef USE_SSE41_OPTIM
    3147             :                     xmm_dst = _mm_blendv_epi8(xmm_dst, xmm_src, xmm_mask);
    3148             : #else
    3149             :                     // mask[i] = 0 becomes 255, and mask[i] != 0 becomes 0
    3150          24 :                     xmm_mask = _mm_cmpeq_epi8(xmm_mask, xmm_zero);
    3151             :                     // dst_data[i] = (mask[i] & dst_data[i]) |
    3152             :                     //               (~mask[i] & src_data[i])
    3153             :                     // That is:
    3154             :                     // dst_data[i] = dst_data[i] when mask[i] = 255
    3155             :                     // dst_data[i] = src_data[i] when mask[i] = 0
    3156          72 :                     xmm_dst = _mm_or_si128(_mm_and_si128(xmm_mask, xmm_dst),
    3157             :                                            _mm_andnot_si128(xmm_mask, xmm_src));
    3158             : #endif
    3159             :                     _mm_storeu_si128(reinterpret_cast<__m128i *>(pabyDestLine),
    3160             :                                      xmm_dst);
    3161          24 :                     pabyDestLine += SIZEOF_REG;
    3162          24 :                     pabySrc += SIZEOF_REG;
    3163          24 :                     iMaskIdx += SIZEOF_REG;
    3164             :                 }
    3165             :             }
    3166             : #endif
    3167         342 :             for (; iX < nOutXSize; iX++)
    3168             :             {
    3169         256 :                 if (pabyMask[iMaskIdx])
    3170             :                 {
    3171         218 :                     *pabyDestLine = *pabySrc;
    3172             :                 }
    3173         256 :                 pabyDestLine += static_cast<GPtrDiff_t>(nPixelSpace);
    3174         256 :                 pabySrc++;
    3175         256 :                 iMaskIdx++;
    3176             :             }
    3177             :         }
    3178             :     }
    3179             :     else
    3180             :     {
    3181          38 :         for (int iY = 0; iY < nOutYSize; iY++)
    3182             :         {
    3183          22 :             GByte *pabyDestLine =
    3184          22 :                 pabyDest + static_cast<GPtrDiff_t>(iY * nLineSpace);
    3185          54 :             for (int iX = 0; iX < nOutXSize; iX++)
    3186             :             {
    3187          32 :                 if (pabyMask[iMaskIdx])
    3188             :                 {
    3189          16 :                     memcpy(pabyDestLine, pabySrc, nBufTypeSize);
    3190             :                 }
    3191          32 :                 pabyDestLine += static_cast<GPtrDiff_t>(nPixelSpace);
    3192          32 :                 pabySrc += nBufTypeSize;
    3193          32 :                 iMaskIdx++;
    3194             :             }
    3195             :         }
    3196             :     }
    3197          66 : }
    3198             : 
    3199             : /************************************************************************/
    3200             : /*                         NeedInitBuffer()                             */
    3201             : /************************************************************************/
    3202             : 
    3203             : // Must be called after CollectSources()
    3204         151 : bool GDALTileIndexDataset::NeedInitBuffer(int nBandCount,
    3205             :                                           const int *panBandMap) const
    3206             : {
    3207         151 :     bool bNeedInitBuffer = true;
    3208             :     // If the last source (that is the most prioritary one) covers at least
    3209             :     // the window of interest and is fully opaque, then we don't need to
    3210             :     // initialize the buffer, and can directly render that source.
    3211         151 :     int bHasNoData = false;
    3212         299 :     if (!m_aoSourceDesc.empty() && m_aoSourceDesc.back().bCoversWholeAOI &&
    3213         145 :         (!m_aoSourceDesc.back().bHasNoData ||
    3214             :          // Also, if there's a single source and that the VRT bands and the
    3215             :          // source bands have the same nodata value, we can skip initialization.
    3216           6 :          (m_aoSourceDesc.size() == 1 && m_aoSourceDesc.back().bSameNoData &&
    3217           4 :           m_bSameNoData && m_bSameDataType &&
    3218           2 :           IsSameNaNAware(papoBands[0]->GetNoDataValue(&bHasNoData),
    3219           2 :                          m_aoSourceDesc.back().dfSameNoData) &&
    3220         301 :           bHasNoData)) &&
    3221         175 :         (!m_aoSourceDesc.back().poMaskBand ||
    3222             :          // Also, if there's a single source that has a mask band, and the VRT
    3223             :          // bands have no-nodata or a 0-nodata value, we can skip
    3224             :          // initialization.
    3225          43 :          (m_aoSourceDesc.size() == 1 && m_bSameDataType &&
    3226           7 :           !(nBandCount == 1 && panBandMap[0] == 0) && m_bSameNoData &&
    3227           7 :           papoBands[0]->GetNoDataValue(&bHasNoData) == 0)))
    3228             :     {
    3229         110 :         bNeedInitBuffer = false;
    3230             :     }
    3231         151 :     return bNeedInitBuffer;
    3232             : }
    3233             : 
    3234             : /************************************************************************/
    3235             : /*                            InitBuffer()                              */
    3236             : /************************************************************************/
    3237             : 
    3238          40 : void GDALTileIndexDataset::InitBuffer(void *pData, int nBufXSize, int nBufYSize,
    3239             :                                       GDALDataType eBufType, int nBandCount,
    3240             :                                       const int *panBandMap,
    3241             :                                       GSpacing nPixelSpace, GSpacing nLineSpace,
    3242             :                                       GSpacing nBandSpace) const
    3243             : {
    3244          40 :     const int nBufTypeSize = GDALGetDataTypeSizeBytes(eBufType);
    3245          40 :     if (m_bSameNoData && nBandCount > 1 &&
    3246           7 :         ((nPixelSpace == nBufTypeSize &&
    3247           7 :           nLineSpace == nBufXSize * nPixelSpace &&
    3248           7 :           nBandSpace == nBufYSize * nLineSpace) ||
    3249           0 :          (nBandSpace == nBufTypeSize &&
    3250           0 :           nPixelSpace == nBandCount * nBandSpace &&
    3251           0 :           nLineSpace == nBufXSize * nPixelSpace)))
    3252             :     {
    3253           7 :         const int nBandNr = panBandMap[0];
    3254             :         auto poVRTBand =
    3255             :             nBandNr == 0
    3256           7 :                 ? m_poMaskBand.get()
    3257           7 :                 : cpl::down_cast<GDALTileIndexBand *>(papoBands[nBandNr - 1]);
    3258           7 :         const double dfNoData = poVRTBand->m_dfNoDataValue;
    3259           7 :         if (dfNoData == 0.0)
    3260             :         {
    3261           5 :             memset(pData, 0,
    3262           5 :                    static_cast<size_t>(nBufXSize) * nBufYSize * nBandCount *
    3263           5 :                        nBufTypeSize);
    3264             :         }
    3265             :         else
    3266             :         {
    3267           2 :             GDALCopyWords64(
    3268             :                 &dfNoData, GDT_Float64, 0, pData, eBufType, nBufTypeSize,
    3269           2 :                 static_cast<size_t>(nBufXSize) * nBufYSize * nBandCount);
    3270           7 :         }
    3271             :     }
    3272             :     else
    3273             :     {
    3274          67 :         for (int i = 0; i < nBandCount; ++i)
    3275             :         {
    3276          34 :             const int nBandNr = panBandMap[i];
    3277          34 :             auto poVRTBand = nBandNr == 0 ? m_poMaskBand.get()
    3278          32 :                                           : cpl::down_cast<GDALTileIndexBand *>(
    3279          32 :                                                 papoBands[nBandNr - 1]);
    3280          34 :             GByte *pabyBandData = static_cast<GByte *>(pData) + i * nBandSpace;
    3281          34 :             if (nPixelSpace == nBufTypeSize &&
    3282          34 :                 poVRTBand->m_dfNoDataValue == 0.0)
    3283             :             {
    3284          30 :                 if (nLineSpace == nBufXSize * nPixelSpace)
    3285             :                 {
    3286          30 :                     memset(pabyBandData, 0,
    3287          30 :                            static_cast<size_t>(nBufYSize * nLineSpace));
    3288             :                 }
    3289             :                 else
    3290             :                 {
    3291           0 :                     for (int iLine = 0; iLine < nBufYSize; iLine++)
    3292             :                     {
    3293           0 :                         memset(static_cast<GByte *>(pabyBandData) +
    3294           0 :                                    static_cast<GIntBig>(iLine) * nLineSpace,
    3295           0 :                                0, static_cast<size_t>(nBufXSize * nPixelSpace));
    3296             :                     }
    3297          30 :                 }
    3298             :             }
    3299             :             else
    3300             :             {
    3301           4 :                 double dfWriteValue = poVRTBand->m_dfNoDataValue;
    3302             : 
    3303          12 :                 for (int iLine = 0; iLine < nBufYSize; iLine++)
    3304             :                 {
    3305           8 :                     GDALCopyWords(&dfWriteValue, GDT_Float64, 0,
    3306           8 :                                   static_cast<GByte *>(pabyBandData) +
    3307           8 :                                       static_cast<GIntBig>(nLineSpace) * iLine,
    3308             :                                   eBufType, static_cast<int>(nPixelSpace),
    3309             :                                   nBufXSize);
    3310             :                 }
    3311             :             }
    3312             :         }
    3313             :     }
    3314          40 : }
    3315             : 
    3316             : /************************************************************************/
    3317             : /*                             IRasterIO()                              */
    3318             : /************************************************************************/
    3319             : 
    3320         148 : CPLErr GDALTileIndexDataset::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
    3321             :                                        int nXSize, int nYSize, void *pData,
    3322             :                                        int nBufXSize, int nBufYSize,
    3323             :                                        GDALDataType eBufType, int nBandCount,
    3324             :                                        int *panBandMap, GSpacing nPixelSpace,
    3325             :                                        GSpacing nLineSpace, GSpacing nBandSpace,
    3326             :                                        GDALRasterIOExtraArg *psExtraArg)
    3327             : {
    3328         148 :     if (eRWFlag != GF_Read)
    3329           0 :         return CE_Failure;
    3330             : 
    3331         148 :     if (nBufXSize < nXSize && nBufYSize < nYSize && AreOverviewsEnabled())
    3332             :     {
    3333           2 :         int bTried = FALSE;
    3334           2 :         const CPLErr eErr = TryOverviewRasterIO(
    3335             :             eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
    3336             :             eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
    3337             :             nBandSpace, psExtraArg, &bTried);
    3338           2 :         if (bTried)
    3339           2 :             return eErr;
    3340             :     }
    3341             : 
    3342         146 :     double dfXOff = nXOff;
    3343         146 :     double dfYOff = nYOff;
    3344         146 :     double dfXSize = nXSize;
    3345         146 :     double dfYSize = nYSize;
    3346         146 :     if (psExtraArg->bFloatingPointWindowValidity)
    3347             :     {
    3348           0 :         dfXOff = psExtraArg->dfXOff;
    3349           0 :         dfYOff = psExtraArg->dfYOff;
    3350           0 :         dfXSize = psExtraArg->dfXSize;
    3351           0 :         dfYSize = psExtraArg->dfYSize;
    3352             :     }
    3353             : 
    3354         146 :     if (!CollectSources(dfXOff, dfYOff, dfXSize, dfYSize))
    3355             :     {
    3356           0 :         return CE_Failure;
    3357             :     }
    3358             : 
    3359             :     // We might be called with nBandCount == 1 && panBandMap[0] == 0
    3360             :     // to mean m_poMaskBand
    3361         146 :     int nBandNrMax = 0;
    3362         319 :     for (int i = 0; i < nBandCount; ++i)
    3363             :     {
    3364         173 :         const int nBandNr = panBandMap[i];
    3365         173 :         nBandNrMax = std::max(nBandNrMax, nBandNr);
    3366             :     }
    3367             : 
    3368         146 :     const bool bNeedInitBuffer = NeedInitBuffer(nBandCount, panBandMap);
    3369             : 
    3370         180 :     const auto RenderSource = [=](SourceDesc &oSourceDesc)
    3371             :     {
    3372         180 :         auto &poTileDS = oSourceDesc.poDS;
    3373         180 :         auto &poSource = oSourceDesc.poSource;
    3374         180 :         auto poComplexSource = dynamic_cast<VRTComplexSource *>(poSource.get());
    3375         180 :         CPLErr eErr = CE_None;
    3376             : 
    3377         180 :         if (poTileDS->GetRasterCount() + 1 == nBandNrMax &&
    3378           6 :             GetRasterBand(nBandNrMax)->GetColorInterpretation() ==
    3379         186 :                 GCI_AlphaBand &&
    3380           4 :             GetRasterBand(nBandNrMax)->GetRasterDataType() == GDT_Byte)
    3381             :         {
    3382             :             // Special case when there's typically a mix of RGB and RGBA source
    3383             :             // datasets and we read a RGB one.
    3384          14 :             for (int iBand = 0; iBand < nBandCount && eErr == CE_None; ++iBand)
    3385             :             {
    3386         339 :                 const int nBandNr = panBandMap[iBand];
    3387          10 :                 if (nBandNr == nBandNrMax)
    3388             :                 {
    3389             :                     // The window we will actually request from the source raster band.
    3390           4 :                     double dfReqXOff = 0.0;
    3391           4 :                     double dfReqYOff = 0.0;
    3392           4 :                     double dfReqXSize = 0.0;
    3393           4 :                     double dfReqYSize = 0.0;
    3394           4 :                     int nReqXOff = 0;
    3395           4 :                     int nReqYOff = 0;
    3396           4 :                     int nReqXSize = 0;
    3397           4 :                     int nReqYSize = 0;
    3398             : 
    3399             :                     // The window we will actual set _within_ the pData buffer.
    3400           4 :                     int nOutXOff = 0;
    3401           4 :                     int nOutYOff = 0;
    3402           4 :                     int nOutXSize = 0;
    3403           4 :                     int nOutYSize = 0;
    3404             : 
    3405           4 :                     bool bError = false;
    3406             : 
    3407           4 :                     auto poTileBand = poTileDS->GetRasterBand(1);
    3408           4 :                     poSource->SetRasterBand(poTileBand, false);
    3409           4 :                     if (poSource->GetSrcDstWindow(
    3410         184 :                             dfXOff, dfYOff, dfXSize, dfYSize, nBufXSize,
    3411           4 :                             nBufYSize, &dfReqXOff, &dfReqYOff, &dfReqXSize,
    3412             :                             &dfReqYSize, &nReqXOff, &nReqYOff, &nReqXSize,
    3413             :                             &nReqYSize, &nOutXOff, &nOutYOff, &nOutXSize,
    3414           4 :                             &nOutYSize, bError))
    3415             :                     {
    3416           4 :                         GByte *pabyOut =
    3417         184 :                             static_cast<GByte *>(pData) +
    3418         195 :                             static_cast<GPtrDiff_t>(iBand * nBandSpace +
    3419         254 :                                                     nOutXOff * nPixelSpace +
    3420         254 :                                                     nOutYOff * nLineSpace);
    3421             : 
    3422           4 :                         constexpr GByte n255 = 255;
    3423           8 :                         for (int iY = 0; iY < nOutYSize; iY++)
    3424             :                         {
    3425           4 :                             GDALCopyWords(&n255, GDT_Byte, 0,
    3426           4 :                                           pabyOut + static_cast<GPtrDiff_t>(
    3427           4 :                                                         iY * nLineSpace),
    3428           4 :                                           eBufType,
    3429             :                                           static_cast<int>(nPixelSpace),
    3430             :                                           nOutXSize);
    3431             :                         }
    3432             :                     }
    3433             :                 }
    3434             :                 else
    3435             :                 {
    3436           6 :                     auto poTileBand = poTileDS->GetRasterBand(nBandNr);
    3437           6 :                     if (poComplexSource)
    3438             :                     {
    3439           0 :                         int bHasNoData = false;
    3440             :                         const double dfNoDataValue =
    3441           0 :                             poTileBand->GetNoDataValue(&bHasNoData);
    3442           0 :                         poComplexSource->SetNoDataValue(
    3443           0 :                             bHasNoData ? dfNoDataValue : VRT_NODATA_UNSET);
    3444             :                     }
    3445           6 :                     poSource->SetRasterBand(poTileBand, false);
    3446             : 
    3447             :                     GDALRasterIOExtraArg sExtraArg;
    3448           6 :                     INIT_RASTERIO_EXTRA_ARG(sExtraArg);
    3449           6 :                     if (psExtraArg->eResampleAlg != GRIORA_NearestNeighbour)
    3450           0 :                         sExtraArg.eResampleAlg = psExtraArg->eResampleAlg;
    3451             :                     else
    3452           6 :                         sExtraArg.eResampleAlg = m_eResampling;
    3453             : 
    3454           6 :                     GByte *pabyBandData =
    3455           6 :                         static_cast<GByte *>(pData) + iBand * nBandSpace;
    3456          12 :                     eErr = poSource->RasterIO(
    3457         125 :                         poTileBand->GetRasterDataType(), nXOff, nYOff, nXSize,
    3458         125 :                         nYSize, pabyBandData, nBufXSize, nBufYSize, eBufType,
    3459           6 :                         nPixelSpace, nLineSpace, &sExtraArg, m_oWorkingState);
    3460             :                 }
    3461             :             }
    3462           4 :             return eErr;
    3463             :         }
    3464         176 :         else if (poTileDS->GetRasterCount() < nBandNrMax)
    3465             :         {
    3466           2 :             CPLError(CE_Failure, CPLE_AppDefined, "%s has not enough bands.",
    3467             :                      oSourceDesc.osName.c_str());
    3468           2 :             return CE_Failure;
    3469             :         }
    3470             : 
    3471         174 :         if ((oSourceDesc.poMaskBand && bNeedInitBuffer) || nBandNrMax == 0)
    3472             :         {
    3473             :             // The window we will actually request from the source raster band.
    3474          55 :             double dfReqXOff = 0.0;
    3475          55 :             double dfReqYOff = 0.0;
    3476          55 :             double dfReqXSize = 0.0;
    3477          55 :             double dfReqYSize = 0.0;
    3478          55 :             int nReqXOff = 0;
    3479          55 :             int nReqYOff = 0;
    3480          55 :             int nReqXSize = 0;
    3481          55 :             int nReqYSize = 0;
    3482             : 
    3483             :             // The window we will actual set _within_ the pData buffer.
    3484          55 :             int nOutXOff = 0;
    3485          55 :             int nOutYOff = 0;
    3486          55 :             int nOutXSize = 0;
    3487          55 :             int nOutYSize = 0;
    3488             : 
    3489          55 :             bool bError = false;
    3490             : 
    3491          55 :             auto poFirstTileBand = poTileDS->GetRasterBand(1);
    3492          55 :             poSource->SetRasterBand(poFirstTileBand, false);
    3493          55 :             if (poSource->GetSrcDstWindow(
    3494             :                     dfXOff, dfYOff, dfXSize, dfYSize, nBufXSize, nBufYSize,
    3495             :                     &dfReqXOff, &dfReqYOff, &dfReqXSize, &dfReqYSize, &nReqXOff,
    3496             :                     &nReqYOff, &nReqXSize, &nReqYSize, &nOutXOff, &nOutYOff,
    3497          55 :                     &nOutXSize, &nOutYSize, bError))
    3498             :             {
    3499          55 :                 int iMaskBandIdx = -1;
    3500          55 :                 if (eBufType == GDT_Byte && nBandNrMax == 0)
    3501             :                 {
    3502             :                     // when called from m_poMaskBand
    3503           4 :                     iMaskBandIdx = 0;
    3504             :                 }
    3505          51 :                 else if (oSourceDesc.poMaskBand)
    3506             :                 {
    3507             :                     // If we request a Byte buffer and the mask band is actually
    3508             :                     // one of the queried bands of this request, we can save
    3509             :                     // requesting it separately.
    3510          51 :                     const int nMaskBandNr = oSourceDesc.poMaskBand->GetBand();
    3511          39 :                     if (eBufType == GDT_Byte && nMaskBandNr >= 1 &&
    3512         129 :                         nMaskBandNr <= poTileDS->GetRasterCount() &&
    3513          39 :                         poTileDS->GetRasterBand(nMaskBandNr) ==
    3514          39 :                             oSourceDesc.poMaskBand)
    3515             :                     {
    3516          61 :                         for (int iBand = 0; iBand < nBandCount; ++iBand)
    3517             :                         {
    3518          44 :                             if (panBandMap[iBand] == nMaskBandNr)
    3519             :                             {
    3520          20 :                                 iMaskBandIdx = iBand;
    3521          20 :                                 break;
    3522             :                             }
    3523             :                         }
    3524             :                     }
    3525             :                 }
    3526             : 
    3527             :                 GDALRasterIOExtraArg sExtraArg;
    3528          55 :                 INIT_RASTERIO_EXTRA_ARG(sExtraArg);
    3529          55 :                 if (psExtraArg->eResampleAlg != GRIORA_NearestNeighbour)
    3530           0 :                     sExtraArg.eResampleAlg = psExtraArg->eResampleAlg;
    3531             :                 else
    3532          55 :                     sExtraArg.eResampleAlg = m_eResampling;
    3533          55 :                 sExtraArg.bFloatingPointWindowValidity = TRUE;
    3534          55 :                 sExtraArg.dfXOff = dfReqXOff;
    3535          55 :                 sExtraArg.dfYOff = dfReqYOff;
    3536          55 :                 sExtraArg.dfXSize = dfReqXSize;
    3537          55 :                 sExtraArg.dfYSize = dfReqYSize;
    3538             : 
    3539          76 :                 if (iMaskBandIdx < 0 && oSourceDesc.abyMask.empty() &&
    3540          21 :                     oSourceDesc.poMaskBand)
    3541             :                 {
    3542             :                     // Fetch the mask band
    3543             :                     try
    3544             :                     {
    3545          21 :                         oSourceDesc.abyMask.resize(
    3546          21 :                             static_cast<size_t>(nOutXSize) * nOutYSize);
    3547             :                     }
    3548           0 :                     catch (const std::bad_alloc &)
    3549             :                     {
    3550           0 :                         CPLError(CE_Failure, CPLE_OutOfMemory,
    3551             :                                  "Cannot allocate working buffer for mask");
    3552           0 :                         return CE_Failure;
    3553             :                     }
    3554             : 
    3555          21 :                     if (oSourceDesc.poMaskBand->RasterIO(
    3556             :                             GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize,
    3557          21 :                             oSourceDesc.abyMask.data(), nOutXSize, nOutYSize,
    3558          21 :                             GDT_Byte, 0, 0, &sExtraArg) != CE_None)
    3559             :                     {
    3560           0 :                         oSourceDesc.abyMask.clear();
    3561           0 :                         return CE_Failure;
    3562             :                     }
    3563             :                 }
    3564             : 
    3565             :                 // Allocate a temporary contiguous buffer to receive pixel data
    3566          55 :                 const int nBufTypeSize = GDALGetDataTypeSizeBytes(eBufType);
    3567          55 :                 const size_t nWorkBufferBandSize =
    3568          55 :                     static_cast<size_t>(nOutXSize) * nOutYSize * nBufTypeSize;
    3569          55 :                 std::vector<GByte> abyWorkBuffer;
    3570             :                 try
    3571             :                 {
    3572          55 :                     abyWorkBuffer.resize(nBandCount * nWorkBufferBandSize);
    3573             :                 }
    3574           0 :                 catch (const std::bad_alloc &)
    3575             :                 {
    3576           0 :                     CPLError(CE_Failure, CPLE_OutOfMemory,
    3577             :                              "Cannot allocate working buffer");
    3578           0 :                     return CE_Failure;
    3579             :                 }
    3580             : 
    3581             :                 const GByte *const pabyMask =
    3582          24 :                     iMaskBandIdx >= 0 ? abyWorkBuffer.data() +
    3583          24 :                                             iMaskBandIdx * nWorkBufferBandSize
    3584          79 :                                       : oSourceDesc.abyMask.data();
    3585             : 
    3586          55 :                 if (nBandNrMax == 0)
    3587             :                 {
    3588             :                     // Special case when called from m_poMaskBand
    3589          12 :                     if (poTileDS->GetRasterBand(1)->GetMaskBand()->RasterIO(
    3590             :                             GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize,
    3591           6 :                             abyWorkBuffer.data(), nOutXSize, nOutYSize,
    3592           6 :                             eBufType, 0, 0, &sExtraArg) != CE_None)
    3593             :                     {
    3594           0 :                         return CE_Failure;
    3595             :                     }
    3596             :                 }
    3597          98 :                 else if (poTileDS->RasterIO(
    3598             :                              GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize,
    3599          49 :                              abyWorkBuffer.data(), nOutXSize, nOutYSize,
    3600             :                              eBufType, nBandCount, panBandMap, 0, 0, 0,
    3601          49 :                              &sExtraArg) != CE_None)
    3602             :                 {
    3603           0 :                     return CE_Failure;
    3604             :                 }
    3605             : 
    3606             :                 // Compose the temporary contiguous buffer into the target
    3607             :                 // buffer, taking into account the mask
    3608          55 :                 GByte *pabyOut =
    3609             :                     static_cast<GByte *>(pData) +
    3610          55 :                     static_cast<GPtrDiff_t>(nOutXOff * nPixelSpace +
    3611          55 :                                             nOutYOff * nLineSpace);
    3612             : 
    3613         121 :                 for (int iBand = 0; iBand < nBandCount && eErr == CE_None;
    3614             :                      ++iBand)
    3615             :                 {
    3616          66 :                     GByte *pabyDestBand =
    3617          66 :                         pabyOut + static_cast<GPtrDiff_t>(iBand * nBandSpace);
    3618             :                     const GByte *pabySrc =
    3619          66 :                         abyWorkBuffer.data() + iBand * nWorkBufferBandSize;
    3620             : 
    3621          66 :                     CompositeSrcWithMaskIntoDest(nOutXSize, nOutYSize, eBufType,
    3622             :                                                  nBufTypeSize, nPixelSpace,
    3623             :                                                  nLineSpace, pabySrc, pabyMask,
    3624             :                                                  pabyDestBand);
    3625             :                 }
    3626          55 :             }
    3627             :         }
    3628         119 :         else if (m_bSameDataType && !bNeedInitBuffer && oSourceDesc.bHasNoData)
    3629             :         {
    3630             :             // We create a non-VRTComplexSource SimpleSource copy of poSource
    3631             :             // to be able to call DatasetRasterIO()
    3632           1 :             VRTSimpleSource oSimpleSource(poSource.get(), 1.0, 1.0);
    3633             : 
    3634             :             GDALRasterIOExtraArg sExtraArg;
    3635           1 :             INIT_RASTERIO_EXTRA_ARG(sExtraArg);
    3636           1 :             if (psExtraArg->eResampleAlg != GRIORA_NearestNeighbour)
    3637           0 :                 sExtraArg.eResampleAlg = psExtraArg->eResampleAlg;
    3638             :             else
    3639           1 :                 sExtraArg.eResampleAlg = m_eResampling;
    3640             : 
    3641           1 :             auto poTileBand = poTileDS->GetRasterBand(panBandMap[0]);
    3642           1 :             oSimpleSource.SetRasterBand(poTileBand, false);
    3643           1 :             eErr = oSimpleSource.DatasetRasterIO(
    3644           1 :                 papoBands[0]->GetRasterDataType(), nXOff, nYOff, nXSize, nYSize,
    3645             :                 pData, nBufXSize, nBufYSize, eBufType, nBandCount, panBandMap,
    3646           1 :                 nPixelSpace, nLineSpace, nBandSpace, &sExtraArg);
    3647             :         }
    3648         118 :         else if (m_bSameDataType && !poComplexSource)
    3649             :         {
    3650         116 :             auto poTileBand = poTileDS->GetRasterBand(panBandMap[0]);
    3651         116 :             poSource->SetRasterBand(poTileBand, false);
    3652             : 
    3653             :             GDALRasterIOExtraArg sExtraArg;
    3654         116 :             INIT_RASTERIO_EXTRA_ARG(sExtraArg);
    3655         116 :             if (poTileBand->GetColorTable())
    3656           0 :                 sExtraArg.eResampleAlg = GRIORA_NearestNeighbour;
    3657         116 :             else if (psExtraArg->eResampleAlg != GRIORA_NearestNeighbour)
    3658           0 :                 sExtraArg.eResampleAlg = psExtraArg->eResampleAlg;
    3659             :             else
    3660         116 :                 sExtraArg.eResampleAlg = m_eResampling;
    3661             : 
    3662         232 :             eErr = poSource->DatasetRasterIO(
    3663         116 :                 papoBands[0]->GetRasterDataType(), nXOff, nYOff, nXSize, nYSize,
    3664             :                 pData, nBufXSize, nBufYSize, eBufType, nBandCount, panBandMap,
    3665         116 :                 nPixelSpace, nLineSpace, nBandSpace, &sExtraArg);
    3666             :         }
    3667             :         else
    3668             :         {
    3669           4 :             for (int i = 0; i < nBandCount && eErr == CE_None; ++i)
    3670             :             {
    3671           2 :                 const int nBandNr = panBandMap[i];
    3672           2 :                 GByte *pabyBandData =
    3673           2 :                     static_cast<GByte *>(pData) + i * nBandSpace;
    3674           2 :                 auto poTileBand = poTileDS->GetRasterBand(nBandNr);
    3675           2 :                 if (poComplexSource)
    3676             :                 {
    3677           2 :                     int bHasNoData = false;
    3678             :                     const double dfNoDataValue =
    3679           2 :                         poTileBand->GetNoDataValue(&bHasNoData);
    3680           2 :                     poComplexSource->SetNoDataValue(
    3681           2 :                         bHasNoData ? dfNoDataValue : VRT_NODATA_UNSET);
    3682             :                 }
    3683           2 :                 poSource->SetRasterBand(poTileBand, false);
    3684             : 
    3685             :                 GDALRasterIOExtraArg sExtraArg;
    3686           2 :                 INIT_RASTERIO_EXTRA_ARG(sExtraArg);
    3687           2 :                 if (poTileBand->GetColorTable())
    3688           0 :                     sExtraArg.eResampleAlg = GRIORA_NearestNeighbour;
    3689           2 :                 else if (psExtraArg->eResampleAlg != GRIORA_NearestNeighbour)
    3690           0 :                     sExtraArg.eResampleAlg = psExtraArg->eResampleAlg;
    3691             :                 else
    3692           2 :                     sExtraArg.eResampleAlg = m_eResampling;
    3693             : 
    3694           4 :                 eErr = poSource->RasterIO(
    3695           2 :                     papoBands[nBandNr - 1]->GetRasterDataType(), nXOff, nYOff,
    3696             :                     nXSize, nYSize, pabyBandData, nBufXSize, nBufYSize,
    3697             :                     eBufType, nPixelSpace, nLineSpace, &sExtraArg,
    3698           2 :                     m_oWorkingState);
    3699             :             }
    3700             :         }
    3701         174 :         return eErr;
    3702         146 :     };
    3703             : 
    3704         146 :     if (!bNeedInitBuffer)
    3705             :     {
    3706         106 :         return RenderSource(m_aoSourceDesc.back());
    3707             :     }
    3708             :     else
    3709             :     {
    3710          40 :         InitBuffer(pData, nBufXSize, nBufYSize, eBufType, nBandCount,
    3711             :                    panBandMap, nPixelSpace, nLineSpace, nBandSpace);
    3712             : 
    3713             :         // Now render from bottom of the stack to top.
    3714         114 :         for (auto &oSourceDesc : m_aoSourceDesc)
    3715             :         {
    3716          74 :             if (oSourceDesc.poDS && RenderSource(oSourceDesc) != CE_None)
    3717           0 :                 return CE_Failure;
    3718             :         }
    3719             : 
    3720          40 :         return CE_None;
    3721             :     }
    3722             : }
    3723             : 
    3724             : /************************************************************************/
    3725             : /*                         GDALRegister_GTI()                           */
    3726             : /************************************************************************/
    3727             : 
    3728        1520 : void GDALRegister_GTI()
    3729             : {
    3730        1520 :     if (GDALGetDriverByName("GTI") != nullptr)
    3731         301 :         return;
    3732             : 
    3733        2438 :     auto poDriver = std::make_unique<VRTDriver>();
    3734             : 
    3735        1219 :     poDriver->SetDescription("GTI");
    3736        1219 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    3737        1219 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "GDAL Raster Tile Index");
    3738        1219 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSIONS, "gti.gpkg gti.fgb gti");
    3739        1219 :     poDriver->SetMetadataItem(GDAL_DMD_CONNECTION_PREFIX, GTI_PREFIX);
    3740        1219 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/gti.html");
    3741             : 
    3742        1219 :     poDriver->pfnOpen = GDALTileIndexDatasetOpen;
    3743        1219 :     poDriver->pfnIdentify = GDALTileIndexDatasetIdentify;
    3744             : 
    3745        1219 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
    3746             : 
    3747        1219 :     poDriver->SetMetadataItem(GDAL_DMD_OPENOPTIONLIST,
    3748             :                               "<OpenOptionList>"
    3749             :                               "  <Option name='LAYER' type='string'/>"
    3750             :                               "  <Option name='LOCATION_FIELD' type='string'/>"
    3751             :                               "  <Option name='SORT_FIELD' type='string'/>"
    3752             :                               "  <Option name='SORT_FIELD_ASC' type='boolean'/>"
    3753             :                               "  <Option name='FILTER' type='string'/>"
    3754             :                               "  <Option name='RESX' type='float'/>"
    3755             :                               "  <Option name='RESY' type='float'/>"
    3756             :                               "  <Option name='MINX' type='float'/>"
    3757             :                               "  <Option name='MINY' type='float'/>"
    3758             :                               "  <Option name='MAXX' type='float'/>"
    3759             :                               "  <Option name='MAXY' type='float'/>"
    3760        1219 :                               "</OpenOptionList>");
    3761             : 
    3762        1219 :     GetGDALDriverManager()->RegisterDriver(poDriver.release());
    3763             : }
    3764             : 
    3765             : /*! @endcond */

Generated by: LCOV version 1.14