/******************************************************************************

 *

 * Project:  Virtual GDAL Datasets

 * Purpose:  Tile index based VRT

 * Author:   Even Rouault <even dot rouault at spatialys.com>

 *

 ******************************************************************************

 * Copyright (c) 2023, Even Rouault <even dot rouault at spatialys.com>

 *

 * SPDX-License-Identifier: MIT

 ****************************************************************************/

/*! @cond Doxygen_Suppress */

#include <array>

#include <algorithm>

#include <atomic>

#include <cmath>

#include <limits>

#include <mutex>

#include <set>

#include <tuple>

#include <utility>

#include <vector>

#include "cpl_port.h"

#include "cpl_error_internal.h"

#include "cpl_json.h"

#include "cpl_mem_cache.h"

#include "cpl_minixml.h"

#include "cpl_quad_tree.h"

#include "vrtdataset.h"

#include "vrt_priv.h"

#include "ogrsf_frmts.h"

#include "ogrwarpedlayer.h"

#include "gdal_frmts.h"

#include "gdal_proxy.h"

#include "gdalsubdatasetinfo.h"

#include "gdal_thread_pool.h"

#include "gdal_utils.h"

#include "gdalalg_raster_index.h"

#ifdef USE_NEON_OPTIMIZATIONS

#define USE_SSE2_OPTIM

#define USE_SSE41_OPTIM

#include "include_sse2neon.h"

#elif defined(__SSE2__) || defined(_M_X64)

#define USE_SSE2_OPTIM

#include <emmintrin.h>

// MSVC doesn't define __SSE4_1__, but if -arch:AVX2 is enabled, we do have SSE4.1

#if defined(__SSE4_1__) || defined(__AVX2__)

#define USE_SSE41_OPTIM

#include <smmintrin.h>

#endif

#endif

#ifndef _

#define _(x) (x)

#endif

// Semantincs of indices of a GeoTransform (double[6]) matrix

constexpr int GT_TOPLEFT_X = 0;

constexpr int GT_WE_RES = 1;

constexpr int GT_ROTATION_PARAM1 = 2;

constexpr int GT_TOPLEFT_Y = 3;

constexpr int GT_ROTATION_PARAM2 = 4;

constexpr int GT_NS_RES = 5;

constexpr const char *GTI_PREFIX = "GTI:";

constexpr const char *MD_DS_TILE_INDEX_LAYER = "TILE_INDEX_LAYER";

constexpr const char *MD_DS_TILE_INDEX_SQL = "TILE_INDEX_SQL";

constexpr const char *MD_DS_TILE_INDEX_SPATIAL_SQL = "TILE_INDEX_SPATIAL_SQL";

constexpr const char *MD_RESX = "RESX";

constexpr const char *MD_RESY = "RESY";

constexpr const char *MD_BAND_COUNT = "BAND_COUNT";

constexpr const char *MD_DATA_TYPE = "DATA_TYPE";

constexpr const char *MD_NODATA = "NODATA";

constexpr const char *MD_MINX = "MINX";

constexpr const char *MD_MINY = "MINY";

constexpr const char *MD_MAXX = "MAXX";

constexpr const char *MD_MAXY = "MAXY";

constexpr const char *MD_GEOTRANSFORM = "GEOTRANSFORM";

constexpr const char *MD_XSIZE = "XSIZE";

constexpr const char *MD_YSIZE = "YSIZE";

constexpr const char *MD_COLOR_INTERPRETATION = "COLOR_INTERPRETATION";

constexpr const char *MD_SRS = "SRS";

constexpr const char *MD_LOCATION_FIELD = "LOCATION_FIELD";

constexpr const char *MD_SORT_FIELD = "SORT_FIELD";

constexpr const char *MD_SORT_FIELD_ASC = "SORT_FIELD_ASC";

constexpr const char *MD_BLOCK_X_SIZE = "BLOCKXSIZE";

constexpr const char *MD_BLOCK_Y_SIZE = "BLOCKYSIZE";

constexpr const char *MD_MASK_BAND = "MASK_BAND";

constexpr const char *MD_RESAMPLING = "RESAMPLING";

constexpr const char *const apszTIOptions[] = {MD_RESX,

                                               MD_RESY,

                                               MD_BAND_COUNT,

                                               MD_DATA_TYPE,

                                               MD_NODATA,

                                               MD_MINX,

                                               MD_MINY,

                                               MD_MAXX,

                                               MD_MAXY,

                                               MD_GEOTRANSFORM,

                                               MD_XSIZE,

                                               MD_YSIZE,

                                               MD_COLOR_INTERPRETATION,

                                               MD_SRS,

                                               MD_LOCATION_FIELD,

                                               MD_SORT_FIELD,

                                               MD_SORT_FIELD_ASC,

                                               MD_BLOCK_X_SIZE,

                                               MD_BLOCK_Y_SIZE,

                                               MD_MASK_BAND,

                                               MD_RESAMPLING};

constexpr const char *const MD_BAND_OFFSET = "OFFSET";

constexpr const char *const MD_BAND_SCALE = "SCALE";

constexpr const char *const MD_BAND_UNITTYPE = "UNITTYPE";

constexpr const char *const apszReservedBandItems[] = {

    MD_BAND_OFFSET, MD_BAND_SCALE, MD_BAND_UNITTYPE};

constexpr const char *GTI_XML_BANDCOUNT = "BandCount";

constexpr const char *GTI_XML_DATATYPE = "DataType";

constexpr const char *GTI_XML_NODATAVALUE = "NoDataValue";

constexpr const char *GTI_XML_COLORINTERP = "ColorInterp";

constexpr const char *GTI_XML_LOCATIONFIELD = "LocationField";

constexpr const char *GTI_XML_SORTFIELD = "SortField";

constexpr const char *GTI_XML_SORTFIELDASC = "SortFieldAsc";

constexpr const char *GTI_XML_MASKBAND = "MaskBand";

constexpr const char *GTI_XML_OVERVIEW_ELEMENT = "Overview";

constexpr const char *GTI_XML_OVERVIEW_DATASET = "Dataset";

constexpr const char *GTI_XML_OVERVIEW_LAYER = "Layer";

constexpr const char *GTI_XML_OVERVIEW_FACTOR = "Factor";

constexpr const char *GTI_XML_BAND_ELEMENT = "Band";

constexpr const char *GTI_XML_BAND_NUMBER = "band";

constexpr const char *GTI_XML_BAND_DATATYPE = "dataType";

constexpr const char *GTI_XML_BAND_DESCRIPTION = "Description";

constexpr const char *GTI_XML_BAND_OFFSET = "Offset";

constexpr const char *GTI_XML_BAND_SCALE = "Scale";

constexpr const char *GTI_XML_BAND_NODATAVALUE = "NoDataValue";

constexpr const char *GTI_XML_BAND_UNITTYPE = "UnitType";

constexpr const char *GTI_XML_BAND_COLORINTERP = "ColorInterp";

constexpr const char *GTI_XML_CATEGORYNAMES = "CategoryNames";

constexpr const char *GTI_XML_COLORTABLE = "ColorTable";

constexpr const char *GTI_XML_RAT = "GDALRasterAttributeTable";

/************************************************************************/

/*                           ENDS_WITH_CI()                             */

/************************************************************************/

static inline bool ENDS_WITH_CI(const char *a, const char *b)

{

    return strlen(a) >= strlen(b) && EQUAL(a + strlen(a) - strlen(b), b);

}

/************************************************************************/

/*                       GDALTileIndexDataset                           */

/************************************************************************/

class GDALTileIndexBand;

class GDALTileIndexDataset final : public GDALPamDataset

{

  public:

    GDALTileIndexDataset();

    ~GDALTileIndexDataset() override;

    bool Open(GDALOpenInfo *poOpenInfo);

    CPLErr FlushCache(bool bAtClosing) override;

    CPLErr GetGeoTransform(GDALGeoTransform &gt) const override;

    const OGRSpatialReference *GetSpatialRef() const override;

    CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,

                     int nYSize, void *pData, int nBufXSize, int nBufYSize,

                     GDALDataType eBufType, int nBandCount,

                     BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,

                     GSpacing nLineSpace, GSpacing nBandSpace,

                     GDALRasterIOExtraArg *psExtraArg) override;

    const char *GetMetadataItem(const char *pszName,

                                const char *pszDomain) override;

    CPLErr SetMetadataItem(const char *pszName, const char *pszValue,

                           const char *pszDomain) override;

    CPLErr SetMetadata(char **papszMD, const char *pszDomain) override;

    void LoadOverviews();

    std::vector<GTISourceDesc> GetSourcesMoreRecentThan(int64_t mTime);

  private:

    friend class GDALTileIndexBand;

    //! Optional GTI XML

    CPLXMLTreeCloser m_psXMLTree{nullptr};

    //! Whether the GTI XML might be modified (by SetMetadata/SetMetadataItem)

    bool m_bXMLUpdatable = false;

    //! Whether the GTI XML has been modified (by SetMetadata/SetMetadataItem)

    bool m_bXMLModified = false;

    //! Unique string (without the process) for this tile index. Passed to

    //! GDALProxyPoolDataset to ensure that sources are unique for a given owner

    const std::string m_osUniqueHandle;

    //! Vector dataset with the sources

    std::unique_ptr<GDALDataset> m_poVectorDS{};

    //! Generic SQL request to return features. May be empty.

    std::string m_osSQL{};

    //! SQL request to return features with placeholders for spatial filtering. May be empty

    std::string m_osSpatialSQL{};

    //! Vector layer with the sources

    OGRLayer *m_poLayer = nullptr;

    //! Whether m_poLayer should be freed with m_poVectorDS->ReleaseResultSet()

    bool m_bIsSQLResultLayer = false;

    //! When the SRS of m_poLayer is not the one we expose

    std::unique_ptr<OGRLayer> m_poWarpedLayerKeeper{};

    //! Geotransform matrix of the tile index

    GDALGeoTransform m_gt{};

    //! Index of the "location" (or alternate name given by user) field

    //! (within m_poLayer->GetLayerDefn()), that contain source dataset names.

    int m_nLocationFieldIndex = -1;

    //! SRS of the tile index.

    OGRSpatialReference m_oSRS{};

    //! Cache from dataset name to dataset handle.

    //! Note that the dataset objects are ultimately GDALProxyPoolDataset,

    //! and that the GDALProxyPoolDataset limits the number of simultaneously

    //! opened real datasets (controlled by GDAL_MAX_DATASET_POOL_SIZE). Hence 500 is not too big.

    lru11::Cache<std::string, std::shared_ptr<GDALDataset>> m_oMapSharedSources{

        500};

    //! Mask band (e.g. for JPEG compressed + mask band)

    std::unique_ptr<GDALTileIndexBand> m_poMaskBand{};

    //! Whether all bands of the tile index have the same data type.

    bool m_bSameDataType = true;

    //! Whether all bands of the tile index have the same nodata value.

    bool m_bSameNoData = true;

    //! Minimum X of the current pixel request, in georeferenced units.

    double m_dfLastMinXFilter = std::numeric_limits<double>::quiet_NaN();

    //! Minimum Y of the current pixel request, in georeferenced units.

    double m_dfLastMinYFilter = std::numeric_limits<double>::quiet_NaN();

    //! Maximum X of the current pixel request, in georeferenced units.

    double m_dfLastMaxXFilter = std::numeric_limits<double>::quiet_NaN();

    //! Maximum Y of the current pixel request, in georeferenced units.

    double m_dfLastMaxYFilter = std::numeric_limits<double>::quiet_NaN();

    //! Index of the field (within m_poLayer->GetLayerDefn()) used to sort, or -1 if none.

    int m_nSortFieldIndex = -1;

    //! Whether sorting must be ascending (true) or descending (false).

    bool m_bSortFieldAsc = true;

    //! Resampling method by default for warping or when a source has not

    //! the same resolution as the tile index.

    std::string m_osResampling = "near";

    GDALRIOResampleAlg m_eResampling = GRIORA_NearestNeighbour;

    //! WKT2 representation of the tile index SRS (if needed, typically for on-the-fly warping).

    std::string m_osWKT{};

    //! Whether we had to open of the sources at tile index opening.

    bool m_bScannedOneFeatureAtOpening = false;

    //! Array of overview descriptors.

    //! Each descriptor is a tuple (dataset_name, concatenated_open_options, layer_name, overview_factor).

    std::vector<std::tuple<std::string, CPLStringList, std::string, double>>

        m_aoOverviewDescriptor{};

    //! Array of overview datasets.

    std::vector<std::unique_ptr<GDALDataset>> m_apoOverviews{};

    //! Cache of buffers used by VRTComplexSource to avoid memory reallocation.

    VRTSource::WorkingState m_oWorkingState{};

    //! Used by IRasterIO() when using multi-threading

    struct QueueWorkingStates

    {

        std::mutex oMutex{};

        std::vector<std::unique_ptr<VRTSource::WorkingState>> oStates{};

    };

    //! Used by IRasterIO() when using multi-threading

    QueueWorkingStates m_oQueueWorkingStates{};

    //! Structure describing one of the source raster in the tile index.

    struct SourceDesc

    {

        //! Source dataset name.

        std::string osName{};

        //! Source dataset handle.

        std::shared_ptr<GDALDataset> poDS{};

        //! VRTSimpleSource or VRTComplexSource for the source.

        std::unique_ptr<VRTSimpleSource> poSource{};

        //! OGRFeature corresponding to the source in the tile index.

        std::unique_ptr<OGRFeature> poFeature{};

        //! Work buffer containing the value of the mask band for the current pixel query.

        mutable std::vector<GByte> abyMask{};

        //! Whether the source covers the whole area of interest of the current pixel query.

        bool bCoversWholeAOI = false;

        //! Whether the source has a nodata value at least in one of its band.

        bool bHasNoData = false;

        //! Whether all bands of the source have the same nodata value.

        bool bSameNoData = false;

        //! Nodata value of all bands (when bSameNoData == true).

        double dfSameNoData = 0;

        //! Mask band of the source.

        GDALRasterBand *poMaskBand = nullptr;

    };

    //! Array of sources participating to the current pixel query.

    std::vector<SourceDesc> m_aoSourceDesc{};

    //! Maximum number of threads. Updated by CollectSources().

    int m_nNumThreads = -1;

    //! Whereas the multi-threading rendering code path must be used. Updated by CollectSources().

    bool m_bLastMustUseMultiThreading = false;

    //! From a source dataset name, return its SourceDesc description structure.

    bool GetSourceDesc(const std::string &osTileName, SourceDesc &oSourceDesc,

                       std::mutex *pMutex);

    //! Collect sources corresponding to the georeferenced window of interest,

    //! and store them in m_aoSourceDesc[].

    bool CollectSources(double dfXOff, double dfYOff, double dfXSize,

                        double dfYSize, bool bMultiThreadAllowed);

    //! Sort sources according to m_nSortFieldIndex.

    void SortSourceDesc();

    //! Whether the output buffer needs to be nodata initialized, or if

    //! sources are fully covering it.

    bool NeedInitBuffer(int nBandCount, const int *panBandMap) const;

    //! Nodata initialize the output buffer.

    void InitBuffer(void *pData, int nBufXSize, int nBufYSize,

                    GDALDataType eBufType, int nBandCount,

                    const int *panBandMap, GSpacing nPixelSpace,

                    GSpacing nLineSpace, GSpacing nBandSpace) const;

    //! Render one source. Used by IRasterIO()

    CPLErr RenderSource(const SourceDesc &oSourceDesc, bool bNeedInitBuffer,

                        int nBandNrMax, int nXOff, int nYOff, int nXSize,

                        int nYSize, double dfXOff, double dfYOff,

                        double dfXSize, double dfYSize, int nBufXSize,

                        int nBufYSize, void *pData, GDALDataType eBufType,

                        int nBandCount, BANDMAP_TYPE panBandMap,

                        GSpacing nPixelSpace, GSpacing nLineSpace,

                        GSpacing nBandSpace, GDALRasterIOExtraArg *psExtraArg,

                        VRTSource::WorkingState &oWorkingState) const;

    //! Whether m_poVectorDS supports SetMetadata()/SetMetadataItem()

    bool TileIndexSupportsEditingLayerMetadata() const;

    //! Return number of threads that can be used

    int GetNumThreads() const;

    /** Structure used to declare a threaded job to satisfy IRasterIO()

     * on a given source.

     */

    struct RasterIOJob

    {

        std::atomic<int> *pnCompletedJobs = nullptr;

        std::atomic<bool> *pbSuccess = nullptr;

        CPLErrorAccumulator *poErrorAccumulator = nullptr;

        GDALTileIndexDataset *poDS = nullptr;

        GDALTileIndexDataset::QueueWorkingStates *poQueueWorkingStates =

            nullptr;

        int nBandNrMax = 0;

        int nXOff = 0;

        int nYOff = 0;

        int nXSize = 0;

        int nYSize = 0;

        void *pData = nullptr;

        int nBufXSize = 0;

        int nBufYSize = 0;

        int nBandCount = 0;

        BANDMAP_TYPE panBandMap = nullptr;

        GDALDataType eBufType = GDT_Unknown;

        GSpacing nPixelSpace = 0;

        GSpacing nLineSpace = 0;

        GSpacing nBandSpace = 0;

        GDALRasterIOExtraArg *psExtraArg = nullptr;

        std::string osTileName{};

        static void Func(void *pData);

    };

    CPL_DISALLOW_COPY_ASSIGN(GDALTileIndexDataset)

};

/************************************************************************/

/*                            GDALTileIndexBand                          */

/************************************************************************/

class GDALTileIndexBand final : public GDALPamRasterBand

{

  public:

    GDALTileIndexBand(GDALTileIndexDataset *poDSIn, int nBandIn,

                      GDALDataType eDT, int nBlockXSizeIn, int nBlockYSizeIn);

    double GetNoDataValue(int *pbHasNoData) override

    {

        if (pbHasNoData)

            *pbHasNoData = m_bNoDataValueSet;

        return m_dfNoDataValue;

    }

    GDALColorInterp GetColorInterpretation() override

    {

        return m_eColorInterp;

    }

    CPLErr IReadBlock(int nBlockXOff, int nBlockYOff, void *pData) override;

    CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,

                     int nYSize, void *pData, int nBufXSize, int nBufYSize,

                     GDALDataType eBufType, GSpacing nPixelSpace,

                     GSpacing nLineSpace,

                     GDALRasterIOExtraArg *psExtraArg) override;

    int IGetDataCoverageStatus(int nXOff, int nYOff, int nXSize, int nYSize,

                               int nMaskFlagStop, double *pdfDataPct) override;

    int GetMaskFlags() override

    {

        if (m_poDS->m_poMaskBand && m_poDS->m_poMaskBand.get() != this)

            return GMF_PER_DATASET;

        return GDALPamRasterBand::GetMaskFlags();

    }

    GDALRasterBand *GetMaskBand() override

    {

        if (m_poDS->m_poMaskBand && m_poDS->m_poMaskBand.get() != this)

            return m_poDS->m_poMaskBand.get();

        return GDALPamRasterBand::GetMaskBand();

    }

    double GetOffset(int *pbHasValue) override

    {

        int bHasValue = FALSE;

        double dfVal = GDALPamRasterBand::GetOffset(&bHasValue);

        if (bHasValue)

        {

            if (pbHasValue)

                *pbHasValue = true;

            return dfVal;

        }

        if (pbHasValue)

            *pbHasValue = !std::isnan(m_dfOffset);

        return std::isnan(m_dfOffset) ? 0.0 : m_dfOffset;

    }

    double GetScale(int *pbHasValue) override

    {

        int bHasValue = FALSE;

        double dfVal = GDALPamRasterBand::GetScale(&bHasValue);

        if (bHasValue)

        {

            if (pbHasValue)

                *pbHasValue = true;

            return dfVal;

        }

        if (pbHasValue)

            *pbHasValue = !std::isnan(m_dfScale);

        return std::isnan(m_dfScale) ? 1.0 : m_dfScale;

    }

    const char *GetUnitType() override

    {

        const char *pszVal = GDALPamRasterBand::GetUnitType();

        if (pszVal && *pszVal)

            return pszVal;

        return m_osUnit.c_str();

    }

    char **GetCategoryNames() override

    {

        return m_aosCategoryNames.List();

    }

    GDALColorTable *GetColorTable() override

    {

        return m_poColorTable.get();

    }

    GDALRasterAttributeTable *GetDefaultRAT() override

    {

        return m_poRAT.get();

    }

    int GetOverviewCount() override;

    GDALRasterBand *GetOverview(int iOvr) override;

    char **GetMetadataDomainList() override;

    const char *GetMetadataItem(const char *pszName,

                                const char *pszDomain) override;

    CPLErr SetMetadataItem(const char *pszName, const char *pszValue,

                           const char *pszDomain) override;

    CPLErr SetMetadata(char **papszMD, const char *pszDomain) override;

  private:

    friend class GDALTileIndexDataset;

    //! Dataset that owns this band.

    GDALTileIndexDataset *m_poDS = nullptr;

    //! Whether a nodata value is set to this band.

    bool m_bNoDataValueSet = false;

    //! Nodata value.

    double m_dfNoDataValue = 0;

    //! Color interpretation.

    GDALColorInterp m_eColorInterp = GCI_Undefined;

    //! Cached value for GetMetadataItem("Pixel_X_Y", "LocationInfo").

    std::string m_osLastLocationInfo{};

    //! Scale value (returned by GetScale())

    double m_dfScale = std::numeric_limits<double>::quiet_NaN();

    //! Offset value (returned by GetOffset())

    double m_dfOffset = std::numeric_limits<double>::quiet_NaN();

    //! Unit type (returned by GetUnitType()).

    std::string m_osUnit{};

    //! Category names (returned by GetCategoryNames()).

    CPLStringList m_aosCategoryNames{};

    //! Color table (returned by GetColorTable()).

    std::unique_ptr<GDALColorTable> m_poColorTable{};

    //! Raster attribute table (returned by GetDefaultRAT()).

    std::unique_ptr<GDALRasterAttributeTable> m_poRAT{};

    CPL_DISALLOW_COPY_ASSIGN(GDALTileIndexBand)

};

/************************************************************************/

/*                        IsSameNaNAware()                              */

/************************************************************************/

static inline bool IsSameNaNAware(double a, double b)

{

    return a == b || (std::isnan(a) && std::isnan(b));

}

/************************************************************************/

/*                         GDALTileIndexDataset()                        */

/************************************************************************/

GDALTileIndexDataset::GDALTileIndexDataset()

    : m_osUniqueHandle(CPLSPrintf("%p", this))

{

}

/************************************************************************/

/*                        GetAbsoluteFileName()                         */

/************************************************************************/

static std::string GetAbsoluteFileName(const char *pszTileName,

                                       const char *pszVRTName)

{

    std::string osRet = VSIURIToVSIPath(pszTileName);

    if (osRet != pszTileName)

        return osRet;

    if (CPLIsFilenameRelative(pszTileName) &&

        !STARTS_WITH(pszTileName, "<VRTDataset") &&

        !STARTS_WITH(pszVRTName, "<GDALTileIndexDataset"))

    {

        const auto oSubDSInfo(GDALGetSubdatasetInfo(pszTileName));

        if (oSubDSInfo && !oSubDSInfo->GetPathComponent().empty())

        {

            const std::string osPath(oSubDSInfo->GetPathComponent());

            osRet = CPLIsFilenameRelative(osPath.c_str())

                        ? oSubDSInfo->ModifyPathComponent(

                              CPLProjectRelativeFilenameSafe(

                                  CPLGetPathSafe(pszVRTName).c_str(),

                                  osPath.c_str()))

                        : std::string(pszTileName);

            GDALDestroySubdatasetInfo(oSubDSInfo);

            return osRet;

        }

        std::string osRelativeMadeAbsolute = CPLProjectRelativeFilenameSafe(

            CPLGetPathSafe(pszVRTName).c_str(), pszTileName);

        VSIStatBufL sStat;

        if (VSIStatL(osRelativeMadeAbsolute.c_str(), &sStat) == 0)

            return osRelativeMadeAbsolute;

    }

    return pszTileName;

}

/************************************************************************/

/*                    GTIDoPaletteExpansionIfNeeded()                   */

/************************************************************************/

//! Do palette -> RGB(A) expansion

static bool

GTIDoPaletteExpansionIfNeeded(std::shared_ptr<GDALDataset> &poTileDS,

                              int nBandCount)

{

    if (poTileDS->GetRasterCount() == 1 &&

        (nBandCount == 3 || nBandCount == 4) &&

        poTileDS->GetRasterBand(1)->GetColorTable() != nullptr)

    {

        CPLStringList aosOptions;

        aosOptions.AddString("-of");

        aosOptions.AddString("VRT");

        aosOptions.AddString("-expand");

        aosOptions.AddString(nBandCount == 3 ? "rgb" : "rgba");

        GDALTranslateOptions *psOptions =

            GDALTranslateOptionsNew(aosOptions.List(), nullptr);

        int bUsageError = false;

        auto poRGBDS = std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(

            GDALTranslate("", GDALDataset::ToHandle(poTileDS.get()), psOptions,

                          &bUsageError)));

        GDALTranslateOptionsFree(psOptions);

        if (!poRGBDS)

        {

            return false;

        }

        poTileDS.reset(poRGBDS.release());

    }

    return true;

}

/************************************************************************/

/*                                Open()                                */

/************************************************************************/

bool GDALTileIndexDataset::Open(GDALOpenInfo *poOpenInfo)

{

    eAccess = poOpenInfo->eAccess;

    CPLXMLNode *psRoot = nullptr;

    const char *pszIndexDataset = poOpenInfo->pszFilename;

    if (STARTS_WITH(poOpenInfo->pszFilename, GTI_PREFIX))

    {

        pszIndexDataset = poOpenInfo->pszFilename + strlen(GTI_PREFIX);

    }

    else if (STARTS_WITH(poOpenInfo->pszFilename, "<GDALTileIndexDataset"))

    {

        // CPLParseXMLString() emits an error in case of failure

        m_psXMLTree.reset(CPLParseXMLString(poOpenInfo->pszFilename));

        if (m_psXMLTree == nullptr)

            return false;

    }

    else if (poOpenInfo->nHeaderBytes > 0 &&

             strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),

                    "<GDALTileIndexDataset"))

    {

        // CPLParseXMLFile() emits an error in case of failure

        m_psXMLTree.reset(CPLParseXMLFile(poOpenInfo->pszFilename));

        if (m_psXMLTree == nullptr)

            return false;

        m_bXMLUpdatable = (poOpenInfo->eAccess == GA_Update);

    }

    if (m_psXMLTree)

    {

        psRoot = CPLGetXMLNode(m_psXMLTree.get(), "=GDALTileIndexDataset");

        if (psRoot == nullptr)

        {

            CPLError(CE_Failure, CPLE_AppDefined,

                     "Missing GDALTileIndexDataset root element.");

            return false;

        }

        pszIndexDataset = CPLGetXMLValue(psRoot, "IndexDataset", nullptr);

        if (!pszIndexDataset)

        {

            CPLError(CE_Failure, CPLE_AppDefined,

                     "Missing IndexDataset element.");

            return false;

        }

    }

    if (ENDS_WITH_CI(pszIndexDataset, ".gti.gpkg") &&

        poOpenInfo->nHeaderBytes >= 100 &&

        STARTS_WITH(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),

                    "SQLite format 3"))

    {

        const char *const apszAllowedDrivers[] = {"GPKG", nullptr};

        m_poVectorDS.reset(GDALDataset::Open(

            std::string("GPKG:\"").append(pszIndexDataset).append("\"").c_str(),

            GDAL_OF_VECTOR | GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR |

                ((poOpenInfo->nOpenFlags & GDAL_OF_UPDATE) ? GDAL_OF_UPDATE

                                                           : GDAL_OF_READONLY),

            apszAllowedDrivers));

        if (!m_poVectorDS)

        {

            return false;

        }

        if (m_poVectorDS->GetLayerCount() == 0 &&

            (m_poVectorDS->GetRasterCount() != 0 ||

             m_poVectorDS->GetMetadata("SUBDATASETS") != nullptr))

        {

            return false;

        }

    }

    else

    {

        m_poVectorDS.reset(GDALDataset::Open(

            pszIndexDataset, GDAL_OF_VECTOR | GDAL_OF_VERBOSE_ERROR |

                                 ((poOpenInfo->nOpenFlags & GDAL_OF_UPDATE)

                                      ? GDAL_OF_UPDATE

                                      : GDAL_OF_READONLY)));

        if (!m_poVectorDS)

        {

            return false;

        }

    }

    if (m_poVectorDS->GetLayerCount() == 0)

    {

        CPLError(CE_Failure, CPLE_AppDefined, "%s has no vector layer",

                 poOpenInfo->pszFilename);

        return false;

    }

    double dfOvrFactor = 1.0;

    if (const char *pszFactor =

            CSLFetchNameValue(poOpenInfo->papszOpenOptions, "FACTOR"))

    {

        dfOvrFactor = CPLAtof(pszFactor);

        if (!(dfOvrFactor > 1.0))

        {

            CPLError(CE_Failure, CPLE_AppDefined, "Wrong overview factor");

            return false;

        }

    }

    m_osSQL = CSLFetchNameValueDef(poOpenInfo->papszOpenOptions, "SQL", "");

    if (m_osSQL.empty())

    {

        if (!psRoot)

        {

            if (const char *pszVal =

                    m_poVectorDS->GetMetadataItem(MD_DS_TILE_INDEX_SQL))

                m_osSQL = pszVal;

        }

        else

            m_osSQL = CPLGetXMLValue(psRoot, "SQL", "");

    }

    if (!m_osSQL.empty())

    {

        m_osSpatialSQL = CSLFetchNameValueDef(poOpenInfo->papszOpenOptions,

                                              "SPATIAL_SQL", "");

        if (m_osSpatialSQL.empty())

        {

            if (!psRoot)

            {

                if (const char *pszVal = m_poVectorDS->GetMetadataItem(

                        MD_DS_TILE_INDEX_SPATIAL_SQL))

                    m_osSpatialSQL = pszVal;

            }

            else

                m_osSpatialSQL = CPLGetXMLValue(psRoot, "SpatialSQL", "");

        }

    }

    const char *pszLayerName;

    if ((pszLayerName = CSLFetchNameValue(poOpenInfo->papszOpenOptions,

                                          "LAYER")) != nullptr)

    {

        m_poLayer = m_poVectorDS->GetLayerByName(pszLayerName);

        if (!m_poLayer)

        {

            CPLError(CE_Failure, CPLE_AppDefined, "Layer %s does not exist",

                     pszLayerName);

            return false;

        }

    }

    else if (psRoot && (pszLayerName = CPLGetXMLValue(psRoot, "IndexLayer",

                                                      nullptr)) != nullptr)

    {

        m_poLayer = m_poVectorDS->GetLayerByName(pszLayerName);

        if (!m_poLayer)

        {

            CPLError(CE_Failure, CPLE_AppDefined, "Layer %s does not exist",

                     pszLayerName);

            return false;

        }

    }

    else if (!psRoot && (pszLayerName = m_poVectorDS->GetMetadataItem(

                             MD_DS_TILE_INDEX_LAYER)) != nullptr)

    {

        m_poLayer = m_poVectorDS->GetLayerByName(pszLayerName);

        if (!m_poLayer)

        {

            CPLError(CE_Failure, CPLE_AppDefined, "Layer %s does not exist",

                     pszLayerName);

            return false;

        }

    }

    else if (!m_osSQL.empty())

    {

        m_poLayer = m_poVectorDS->ExecuteSQL(m_osSQL.c_str(), nullptr, nullptr);

        if (!m_poLayer)

        {

            CPLError(CE_Failure, CPLE_AppDefined, "SQL request %s failed",

                     m_osSQL.c_str());

            return false;

        }

        m_bIsSQLResultLayer = true;

    }

    else if (m_poVectorDS->GetLayerCount() == 1)

    {

        m_poLayer = m_poVectorDS->GetLayer(0);

        if (!m_poLayer)

        {

            CPLError(CE_Failure, CPLE_AppDefined, "Cannot open layer 0");

            return false;

        }

    }

    else

    {

        if (STARTS_WITH(poOpenInfo->pszFilename, GTI_PREFIX))

        {

            CPLError(CE_Failure, CPLE_AppDefined,

                     "%s has more than one layer. LAYER open option "

                     "must be defined to specify which one to "

                     "use as the tile index",

                     pszIndexDataset);

        }

        else if (psRoot)

        {

            CPLError(CE_Failure, CPLE_AppDefined,

                     "%s has more than one layer. IndexLayer element must be "

                     "defined to specify which one to "

                     "use as the tile index",

                     pszIndexDataset);

        }

        else

        {

            CPLError(CE_Failure, CPLE_AppDefined,

                     "%s has more than one layer. %s "

                     "metadata item must be defined to specify which one to "

                     "use as the tile index",

                     pszIndexDataset, MD_DS_TILE_INDEX_LAYER);

        }

        return false;

    }

    // Try to get the metadata from an embedded xml:GTI domain

    if (!m_psXMLTree)

    {

        char **papszMD = m_poLayer->GetMetadata("xml:GTI");

        if (papszMD && papszMD[0])

        {

            m_psXMLTree.reset(CPLParseXMLString(papszMD[0]));

            if (m_psXMLTree == nullptr)

                return false;

            psRoot = CPLGetXMLNode(m_psXMLTree.get(), "=GDALTileIndexDataset");

            if (psRoot == nullptr)

            {

                CPLError(CE_Failure, CPLE_AppDefined,

                         "Missing GDALTileIndexDataset root element.");

                return false;

            }

        }

    }

    // Get the value of an option.

    // The order of lookup is the following one (first to last):

    // - open options

    // - XML file

    // - Layer metadata items.

    const auto GetOption = [poOpenInfo, psRoot, this](const char *pszItem)

    {

        const char *pszVal =

            CSLFetchNameValue(poOpenInfo->papszOpenOptions, pszItem);

        if (pszVal)

            return pszVal;

        if (psRoot)

        {

            pszVal = CPLGetXMLValue(psRoot, pszItem, nullptr);

            if (pszVal)

                return pszVal;

            if (EQUAL(pszItem, MD_BAND_COUNT))

                pszItem = GTI_XML_BANDCOUNT;

            else if (EQUAL(pszItem, MD_DATA_TYPE))

                pszItem = GTI_XML_DATATYPE;

            else if (EQUAL(pszItem, MD_NODATA))

                pszItem = GTI_XML_NODATAVALUE;

            else if (EQUAL(pszItem, MD_COLOR_INTERPRETATION))

                pszItem = GTI_XML_COLORINTERP;

            else if (EQUAL(pszItem, MD_LOCATION_FIELD))

                pszItem = GTI_XML_LOCATIONFIELD;

            else if (EQUAL(pszItem, MD_SORT_FIELD))

                pszItem = GTI_XML_SORTFIELD;

            else if (EQUAL(pszItem, MD_SORT_FIELD_ASC))

                pszItem = GTI_XML_SORTFIELDASC;

            else if (EQUAL(pszItem, MD_MASK_BAND))

                pszItem = GTI_XML_MASKBAND;

            pszVal = CPLGetXMLValue(psRoot, pszItem, nullptr);

            if (pszVal)

                return pszVal;

        }

        return m_poLayer->GetMetadataItem(pszItem);

    };

    const char *pszFilter = GetOption("Filter");

    if (pszFilter)

    {

        if (m_poLayer->SetAttributeFilter(pszFilter) != OGRERR_NONE)

            return false;

    }

    const OGRFeatureDefn *poLayerDefn = m_poLayer->GetLayerDefn();

    std::string osLocationFieldName;

    {

        const char *pszLocationFieldName = GetOption(MD_LOCATION_FIELD);

        if (pszLocationFieldName)

        {

            osLocationFieldName = pszLocationFieldName;

        }

        else

        {

            // Is this a https://stac-utils.github.io/stac-geoparquet/latest/spec/stac-geoparquet-spec ?

            if (poLayerDefn->GetFieldIndex("assets.data.href") >= 0)

            {

                osLocationFieldName = "assets.data.href";

                CPLDebug("GTI", "Using %s as location field",

                         osLocationFieldName.c_str());

            }

            else if (poLayerDefn->GetFieldIndex("assets.image.href") >= 0)

            {

                osLocationFieldName = "assets.image.href";

                CPLDebug("GTI", "Using %s as location field",

                         osLocationFieldName.c_str());

            }

            else if (poLayerDefn->GetFieldIndex("stac_version") >= 0)

            {

                const int nFieldCount = poLayerDefn->GetFieldCount();

                // Look for "assets.xxxxx.href" fields

                int nAssetCount = 0;

                for (int i = 0; i < nFieldCount; ++i)

                {

                    const auto poFDefn = poLayerDefn->GetFieldDefn(i);

                    const char *pszFieldName = poFDefn->GetNameRef();

                    if (STARTS_WITH(pszFieldName, "assets.") &&

                        EQUAL(pszFieldName + strlen(pszFieldName) -

                                  strlen(".href"),

                              ".href") &&

                        // Assets with "metadata" in them are very much likely

                        // not rasters... We could potentially confirm that by

                        // inspecting the value of the assets.XXX.type or

                        // assets.XXX.roles fields of one feature

                        !strstr(pszFieldName, "metadata"))