/src/gdal/frmts/vrt/vrtpansharpened.cpp

Source
/******************************************************************************
 *
 * Project:  Virtual GDAL Datasets
 * Purpose:  Implementation of VRTPansharpenedRasterBand and
 *VRTPansharpenedDataset. Author:   Even Rouault <even.rouault at spatialys.com>
 *
 ******************************************************************************
 * Copyright (c) 2015, Even Rouault <even.rouault at spatialys.com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "cpl_port.h"
#include "gdal_vrt.h"
#include "vrtdataset.h"

#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdlib>
#include <cstring>

#include <algorithm>
#include <limits>
#include <map>
#include <set>
#include <string>
#include <utility>
#include <vector>

#include "cpl_conv.h"
#include "cpl_error.h"
#include "cpl_minixml.h"
#include "cpl_string.h"
#include "cpl_vsi.h"
#include "gdal.h"
#include "gdal_priv.h"
#include "gdalpansharpen.h"
#include "ogr_core.h"
#include "ogr_spatialref.h"

/************************************************************************/
/*                     GDALCreatePansharpenedVRT()                      */
/************************************************************************/

/**
 * Create a virtual pansharpened dataset.
 *
 * This function will create a virtual pansharpened dataset.
 *
 * Starting with GDAL 3.12, a reference will be taken on the dataset owing
 * each input bands.
 * Before 3.12, no reference were taken on the passed bands. Consequently,
 * they or their dataset to which they belong had to be kept open until
 * this virtual pansharpened dataset is closed.
 *
 * The returned dataset will have no associated filename for itself.  If you
 * want to write the virtual dataset description to a file, use the
 * GDALSetDescription() function (or SetDescription() method) on the dataset
 * to assign a filename before it is closed.
 *
 * @param pszXML Pansharpened VRT XML where &lt;SpectralBand&gt; elements have
 * no explicit SourceFilename and SourceBand. The spectral bands in the XML will
 * be assigned the successive values of the pahInputSpectralBands array. Must
 * not be NULL.
 *
 * @param hPanchroBand Panchromatic band. Must not be NULL.
 *
 * @param nInputSpectralBands Number of input spectral bands.  Must be greater
 * than zero.
 *
 * @param pahInputSpectralBands Array of nInputSpectralBands spectral bands.
 *
 * @return NULL on failure, or a new virtual dataset handle on success to be
 * closed with GDALClose().
 *
 * @since GDAL 2.1
 */

GDALDatasetH GDALCreatePansharpenedVRT(const char *pszXML,
                                       GDALRasterBandH hPanchroBand,
                                       int nInputSpectralBands,
                                       GDALRasterBandH *pahInputSpectralBands)
{
    VALIDATE_POINTER1(pszXML, "GDALCreatePansharpenedVRT", nullptr);
    VALIDATE_POINTER1(hPanchroBand, "GDALCreatePansharpenedVRT", nullptr);
    VALIDATE_POINTER1(pahInputSpectralBands, "GDALCreatePansharpenedVRT",
                      nullptr);

    CPLXMLNode *psTree = CPLParseXMLString(pszXML);
    if (psTree == nullptr)
        return nullptr;
    VRTPansharpenedDataset *poDS = new VRTPansharpenedDataset(0, 0);
    CPLErr eErr = poDS->XMLInit(psTree, nullptr, hPanchroBand,
                                nInputSpectralBands, pahInputSpectralBands);
    CPLDestroyXMLNode(psTree);
    if (eErr != CE_None)
    {
        delete poDS;
        return nullptr;
    }
    return GDALDataset::ToHandle(poDS);
}

/*! @cond Doxygen_Suppress */

/************************************************************************/
/* ==================================================================== */
/*                        VRTPansharpenedDataset                        */
/* ==================================================================== */
/************************************************************************/

/************************************************************************/
/*                       VRTPansharpenedDataset()                       */
/************************************************************************/

VRTPansharpenedDataset::VRTPansharpenedDataset(int nXSize, int nYSize,
                                               int nBlockXSize, int nBlockYSize)
    : VRTDataset(nXSize, nYSize,
                 nBlockXSize > 0 ? nBlockXSize : std::min(nXSize, 512),
                 nBlockYSize > 0 ? nBlockYSize : std::min(nYSize, 512)),
      m_poMainDataset(nullptr), m_bLoadingOtherBands(FALSE),
      m_pabyLastBufferBandRasterIO(nullptr), m_nLastBandRasterIOXOff(0),
      m_nLastBandRasterIOYOff(0), m_nLastBandRasterIOXSize(0),
      m_nLastBandRasterIOYSize(0), m_eLastBandRasterIODataType(GDT_Unknown),
      m_eGTAdjustment(GTAdjust_Union), m_bNoDataDisabled(FALSE)
{
    eAccess = GA_Update;
    m_poMainDataset = this;
}

/************************************************************************/
/*                      ~VRTPansharpenedDataset()                       */
/************************************************************************/

VRTPansharpenedDataset::~VRTPansharpenedDataset()

{
    VRTPansharpenedDataset::FlushCache(true);
    VRTPansharpenedDataset::CloseDependentDatasets();
    CPLFree(m_pabyLastBufferBandRasterIO);
}

/************************************************************************/
/*                       CloseDependentDatasets()                       */
/************************************************************************/

int VRTPansharpenedDataset::CloseDependentDatasets()
{
    int bHasDroppedRef = VRTDataset::CloseDependentDatasets();

    /* -------------------------------------------------------------------- */
    /*      Destroy the raster bands if they exist.                         */
    /* -------------------------------------------------------------------- */
    for (int iBand = 0; iBand < nBands; iBand++)
    {
        delete papoBands[iBand];
    }
    nBands = 0;

    // Destroy the overviews before m_poPansharpener as they might reference
    // files that are in m_apoDatasetsToReleaseRef.
    bHasDroppedRef |= !m_apoOverviewDatasets.empty();
    m_apoOverviewDatasets.clear();

    if (m_poPansharpener != nullptr)
    {
        // Delete the pansharper object before closing the dataset
        // because it may have warped the bands into an intermediate VRT
        m_poPansharpener.reset();

        // Close in reverse order (VRT firsts and real datasets after)
        for (int i = static_cast<int>(m_apoDatasetsToReleaseRef.size()) - 1;
             i >= 0; i--)
        {
            bHasDroppedRef = TRUE;
            m_apoDatasetsToReleaseRef[i].reset();
        }
        m_apoDatasetsToReleaseRef.clear();
    }

    return bHasDroppedRef;
}

/************************************************************************/
/*                            GetFileList()                             */
/************************************************************************/

char **VRTPansharpenedDataset::GetFileList()
{
    char **papszFileList = GDALDataset::GetFileList();

    if (m_poPansharpener != nullptr)
    {
        const GDALPansharpenOptions *psOptions = m_poPansharpener->GetOptions();
        if (psOptions != nullptr)
        {
            std::set<CPLString> oSetNames;
            if (psOptions->hPanchroBand != nullptr)
            {
                GDALDataset *poDS = GDALDataset::FromHandle(
                    GDALGetBandDataset(psOptions->hPanchroBand));
                if (poDS != nullptr)
                {
                    const char *pszName =
                        poDS->GetMetadataItem("UNDERLYING_DATASET");
                    if (!pszName)
                        pszName = poDS->GetDescription();
                    papszFileList = CSLAddString(papszFileList, pszName);
                    oSetNames.insert(pszName);
                }
            }
            for (int i = 0; i < psOptions->nInputSpectralBands; i++)
            {
                if (psOptions->pahInputSpectralBands[i] != nullptr)
                {
                    GDALDataset *poDS =
                        GDALDataset::FromHandle(GDALGetBandDataset(
                            psOptions->pahInputSpectralBands[i]));
                    if (poDS != nullptr)
                    {
                        const char *pszName =
                            poDS->GetMetadataItem("UNDERLYING_DATASET");
                        if (!pszName)
                            pszName = poDS->GetDescription();
                        if (!cpl::contains(oSetNames, pszName))
                        {
                            papszFileList =
                                CSLAddString(papszFileList, pszName);
                            oSetNames.insert(pszName);
                        }
                    }
                }
            }
        }
    }

    return papszFileList;
}

/************************************************************************/
/*                              XMLInit()                               */
/************************************************************************/

CPLErr VRTPansharpenedDataset::XMLInit(const CPLXMLNode *psTree,
                                       const char *pszVRTPathIn)

{
    return XMLInit(psTree, pszVRTPathIn, nullptr, 0, nullptr);
}

CPLErr VRTPansharpenedDataset::XMLInit(const CPLXMLNode *psTree,
                                       const char *pszVRTPathIn,
                                       GDALRasterBandH hPanchroBandIn,
                                       int nInputSpectralBandsIn,
                                       GDALRasterBandH *pahInputSpectralBandsIn)
{
    CPLErr eErr;
    std::unique_ptr<GDALPansharpenOptions,
                    decltype(&GDALDestroyPansharpenOptions)>
        psPanOptions(nullptr, GDALDestroyPansharpenOptions);

    /* -------------------------------------------------------------------- */
    /*      Initialize blocksize before calling sub-init so that the        */
    /*      band initializers can get it from the dataset object when       */
    /*      they are created.                                               */
    /* -------------------------------------------------------------------- */
    m_nBlockXSize = atoi(CPLGetXMLValue(psTree, "BlockXSize", "512"));
    m_nBlockYSize = atoi(CPLGetXMLValue(psTree, "BlockYSize", "512"));

    /* -------------------------------------------------------------------- */
    /*      Parse PansharpeningOptions                                      */
    /* -------------------------------------------------------------------- */

    const CPLXMLNode *psOptions = CPLGetXMLNode(psTree, "PansharpeningOptions");
    if (psOptions == nullptr)
    {
        CPLError(CE_Failure, CPLE_AppDefined, "Missing PansharpeningOptions");
        return CE_Failure;
    }

    if (CPLGetXMLValue(psOptions, "MSShiftX", nullptr) != nullptr ||
        CPLGetXMLValue(psOptions, "MSShiftY", nullptr) != nullptr)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Options MSShiftX and MSShiftY are no longer supported. "
                 "Spatial adjustment between panchromatic and multispectral "
                 "bands is done through their geotransform.");
        return CE_Failure;
    }

    CPLString osSourceFilename;
    GDALDataset *poPanDataset = nullptr;
    GDALRasterBand *poPanBand = nullptr;
    std::map<CPLString, GDALDataset *> oMapNamesToDataset;
    int nPanBand;

    if (hPanchroBandIn == nullptr)
    {
        const CPLXMLNode *psPanchroBand =
            CPLGetXMLNode(psOptions, "PanchroBand");
        if (psPanchroBand == nullptr)
        {
            CPLError(CE_Failure, CPLE_AppDefined, "PanchroBand missing");
            return CE_Failure;
        }

        osSourceFilename = CPLGetXMLValue(psPanchroBand, "SourceFilename", "");
        if (osSourceFilename.empty())
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "PanchroBand.SourceFilename missing");
            return CE_Failure;
        }
        const bool bRelativeToVRT = CPL_TO_BOOL(atoi(CPLGetXMLValue(
            psPanchroBand, "SourceFilename.relativetoVRT", "0")));
        if (bRelativeToVRT)
        {
            const std::string osAbs = CPLProjectRelativeFilenameSafe(
                pszVRTPathIn, osSourceFilename.c_str());
            m_oMapToRelativeFilenames[osAbs] = osSourceFilename;
            osSourceFilename = osAbs;
        }

        const CPLStringList aosOpenOptions(
            GDALDeserializeOpenOptionsFromXML(psPanchroBand));

        poPanDataset = GDALDataset::Open(osSourceFilename,
                                         GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR,
                                         nullptr, aosOpenOptions.List());
        if (poPanDataset == nullptr)
        {
            return CE_Failure;
        }

        const char *pszSourceBand =
            CPLGetXMLValue(psPanchroBand, "SourceBand", "1");
        nPanBand = atoi(pszSourceBand);
        if (poPanBand == nullptr)
            poPanBand = poPanDataset->GetRasterBand(nPanBand);
        if (poPanBand == nullptr)
        {
            CPLError(CE_Failure, CPLE_AppDefined, "%s invalid band of %s",
                     pszSourceBand, osSourceFilename.c_str());
            GDALClose(poPanDataset);
            return CE_Failure;
        }
        oMapNamesToDataset[osSourceFilename] = poPanDataset;
    }
    else
    {
        poPanBand = GDALRasterBand::FromHandle(hPanchroBandIn);
        nPanBand = poPanBand->GetBand();
        poPanDataset = poPanBand->GetDataset();
        if (poPanDataset == nullptr)
        {
            CPLError(
                CE_Failure, CPLE_AppDefined,
                "Cannot retrieve dataset associated with panchromatic band");
            return CE_Failure;
        }
        oMapNamesToDataset[CPLSPrintf("%p", poPanDataset)] = poPanDataset;
        poPanDataset->Reference();
    }
    m_apoDatasetsToReleaseRef.push_back(
        std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>(
            poPanDataset));
    poPanDataset = m_apoDatasetsToReleaseRef.back().get();
    // Figure out which kind of adjustment we should do if the pan and spectral
    // bands do not share the same geotransform
    const char *pszGTAdjustment =
        CPLGetXMLValue(psOptions, "SpatialExtentAdjustment", "Union");
    if (EQUAL(pszGTAdjustment, "Union"))
        m_eGTAdjustment = GTAdjust_Union;
    else if (EQUAL(pszGTAdjustment, "Intersection"))
        m_eGTAdjustment = GTAdjust_Intersection;
    else if (EQUAL(pszGTAdjustment, "None"))
        m_eGTAdjustment = GTAdjust_None;
    else if (EQUAL(pszGTAdjustment, "NoneWithoutWarning"))
        m_eGTAdjustment = GTAdjust_NoneWithoutWarning;
    else
    {
        m_eGTAdjustment = GTAdjust_Union;
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Unsupported value for GeoTransformAdjustment. Defaulting to "
                 "Union");
    }

    const char *pszNumThreads =
        CPLGetXMLValue(psOptions, "NumThreads", nullptr);
    int nThreads = 0;
    if (pszNumThreads != nullptr)
    {
        if (EQUAL(pszNumThreads, "ALL_CPUS"))
            nThreads = -1;
        else
            nThreads = atoi(pszNumThreads);
    }

    const char *pszAlgorithm =
        CPLGetXMLValue(psOptions, "Algorithm", "WeightedBrovey");
    if (!EQUAL(pszAlgorithm, "WeightedBrovey"))
    {
        CPLError(CE_Failure, CPLE_AppDefined, "Algorithm %s unsupported",
                 pszAlgorithm);
        return CE_Failure;
    }

    std::vector<double> adfWeights;
    if (const CPLXMLNode *psAlgOptions =
            CPLGetXMLNode(psOptions, "AlgorithmOptions"))
    {
        const char *pszWeights =
            CPLGetXMLValue(psAlgOptions, "Weights", nullptr);
        if (pszWeights != nullptr)
        {
            char **papszTokens = CSLTokenizeString2(pszWeights, " ,", 0);
            for (int i = 0; papszTokens && papszTokens[i]; i++)
                adfWeights.push_back(CPLAtof(papszTokens[i]));
            CSLDestroy(papszTokens);
        }
    }

    GDALRIOResampleAlg eResampleAlg = GDALRasterIOGetResampleAlg(
        CPLGetXMLValue(psOptions, "Resampling", "Cubic"));

    std::vector<GDALRasterBand *> ahSpectralBands;
    std::map<int, int> aMapDstBandToSpectralBand;
    std::map<int, int>::iterator aMapDstBandToSpectralBandIter;
    int nBitDepth = 0;
    bool bFoundNonMatchingGT = false;
    GDALGeoTransform panGT;
    const bool bPanGeoTransformValid =
        (poPanDataset->GetGeoTransform(panGT) == CE_None);
    if (!bPanGeoTransformValid)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Panchromatic band has no associated geotransform");
        return CE_Failure;
    }
    int nPanXSize = poPanBand->GetXSize();
    int nPanYSize = poPanBand->GetYSize();
    int bHasNoData = FALSE;
    double dfNoData = poPanBand->GetNoDataValue(&bHasNoData);
    double dfLRPanX = panGT[0] + nPanXSize * panGT[1] + nPanYSize * panGT[2];
    double dfLRPanY = panGT[3] + nPanXSize * panGT[4] + nPanYSize * panGT[5];
    bool bFoundRotatingTerms = (panGT[2] != 0.0 || panGT[4] != 0.0);
    double dfMinX = panGT[0];
    double dfMaxX = dfLRPanX;
    double dfMaxY = panGT[3];
    double dfMinY = dfLRPanY;
    if (dfMinY > dfMaxY)
        std::swap(dfMinY, dfMaxY);
    const double dfPanMinY = dfMinY;
    const double dfPanMaxY = dfMaxY;

    const auto poPanSRS = poPanDataset->GetSpatialRef();

    /* -------------------------------------------------------------------- */
    /*      First pass on spectral datasets to check their georeferencing.  */
    /* -------------------------------------------------------------------- */
    int iSpectralBand = 0;
    for (const CPLXMLNode *psIter = psOptions->psChild; psIter;
         psIter = psIter->psNext)
    {
        GDALDataset *poDataset;
        if (psIter->eType != CXT_Element ||
            !EQUAL(psIter->pszValue, "SpectralBand"))
            continue;

        if (nInputSpectralBandsIn && pahInputSpectralBandsIn != nullptr)
        {
            if (iSpectralBand == nInputSpectralBandsIn)
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "More SpectralBand elements than in source array");
                goto error;
            }
            poDataset = GDALRasterBand::FromHandle(
                            pahInputSpectralBandsIn[iSpectralBand])
                            ->GetDataset();
            if (poDataset == nullptr)
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "SpectralBand has no associated dataset");
                goto error;
            }
            osSourceFilename = poDataset->GetDescription();

            oMapNamesToDataset[CPLSPrintf("%p", poDataset)] = poDataset;
            poDataset->Reference();
        }
        else
        {
            osSourceFilename = CPLGetXMLValue(psIter, "SourceFilename", "");
            if (osSourceFilename.empty())
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "SpectralBand.SourceFilename missing");
                goto error;
            }
            int bRelativeToVRT = atoi(
                CPLGetXMLValue(psIter, "SourceFilename.relativetoVRT", "0"));
            if (bRelativeToVRT)
            {
                const std::string osAbs = CPLProjectRelativeFilenameSafe(
                    pszVRTPathIn, osSourceFilename.c_str());
                m_oMapToRelativeFilenames[osAbs] = osSourceFilename;
                osSourceFilename = osAbs;
            }
            poDataset = oMapNamesToDataset[osSourceFilename];
            if (poDataset == nullptr)
            {
                const CPLStringList aosOpenOptions(
                    GDALDeserializeOpenOptionsFromXML(psIter));

                poDataset = GDALDataset::Open(
                    osSourceFilename, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR,
                    nullptr, aosOpenOptions.List());
                if (poDataset == nullptr)
                {
                    goto error;
                }
                oMapNamesToDataset[osSourceFilename] = poDataset;
            }
            else
            {
                poDataset->Reference();
            }
        }
        m_apoDatasetsToReleaseRef.push_back(
            std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>(
                poDataset));
        poDataset = m_apoDatasetsToReleaseRef.back().get();

        // Check that the spectral band has a georeferencing consistent
        // of the pan band. Allow an error of at most the size of one pixel
        // of the spectral band.
        const auto poSpectralSRS = poDataset->GetSpatialRef();
        if (poPanSRS)
        {
            if (poSpectralSRS)
            {
                if (!poPanSRS->IsSame(poSpectralSRS))
                {
                    CPLError(CE_Warning, CPLE_AppDefined,
                             "Pan dataset and %s do not seem to "
                             "have same projection. Results might "
                             "be incorrect",
                             osSourceFilename.c_str());
                }
            }
            else
            {
                CPLError(CE_Warning, CPLE_AppDefined,
                         "Pan dataset has a projection, whereas %s "
                         "not. Results might be incorrect",
                         osSourceFilename.c_str());
            }
        }
        else if (poSpectralSRS)
        {
            CPLError(CE_Warning, CPLE_AppDefined,
                     "Pan dataset has no projection, whereas %s has "
                     "one. Results might be incorrect",
                     osSourceFilename.c_str());
        }

        GDALGeoTransform spectralGT;
        if (poDataset->GetGeoTransform(spectralGT) != CE_None)
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "Spectral band has no associated geotransform");
            goto error;
        }
        if (spectralGT[5] * panGT[5] < 0)
        {
            CPLError(CE_Failure, CPLE_NotSupported,
                     "Spectral band vertical orientation is "
                     "different from pan one");
            goto error;
        }

        int bIsThisOneNonMatching = FALSE;
        double dfPixelSize = std::max(spectralGT[1], std::abs(spectralGT[5]));
        if (std::abs(panGT[0] - spectralGT[0]) > dfPixelSize ||
            std::abs(panGT[3] - spectralGT[3]) > dfPixelSize)
        {
            bIsThisOneNonMatching = TRUE;
            if (m_eGTAdjustment == GTAdjust_None)
            {
                CPLError(CE_Warning, CPLE_AppDefined,
                         "Georeferencing of top-left corner of pan "
                         "dataset and %s do not match",
                         osSourceFilename.c_str());
            }
        }
        bFoundRotatingTerms = bFoundRotatingTerms ||
                              (spectralGT[2] != 0.0 || spectralGT[4] != 0.0);
        double dfLRSpectralX = spectralGT[0] +
                               poDataset->GetRasterXSize() * spectralGT[1] +
                               poDataset->GetRasterYSize() * spectralGT[2];
        double dfLRSpectralY = spectralGT[3] +
                               poDataset->GetRasterXSize() * spectralGT[4] +
                               poDataset->GetRasterYSize() * spectralGT[5];
        if (std::abs(dfLRPanX - dfLRSpectralX) > dfPixelSize ||
            std::abs(dfLRPanY - dfLRSpectralY) > dfPixelSize)
        {
            bIsThisOneNonMatching = TRUE;
            if (m_eGTAdjustment == GTAdjust_None)
            {
                CPLError(CE_Warning, CPLE_AppDefined,
                         "Georeferencing of bottom-right corner of "
                         "pan dataset and %s do not match",
                         osSourceFilename.c_str());
            }
        }

        double dfThisMinY = dfLRSpectralY;
        double dfThisMaxY = spectralGT[3];
        if (dfThisMinY > dfThisMaxY)
            std::swap(dfThisMinY, dfThisMaxY);
        if (bIsThisOneNonMatching && m_eGTAdjustment == GTAdjust_Union)
        {
            dfMinX = std::min(dfMinX, spectralGT[0]);
            dfMinY = std::min(dfMinY, dfThisMinY);
            dfMaxX = std::max(dfMaxX, dfLRSpectralX);
            dfMaxY = std::max(dfMaxY, dfThisMaxY);
        }
        else if (bIsThisOneNonMatching &&
                 m_eGTAdjustment == GTAdjust_Intersection)
        {
            dfMinX = std::max(dfMinX, spectralGT[0]);
            dfMinY = std::max(dfMinY, dfThisMinY);
            dfMaxX = std::min(dfMaxX, dfLRSpectralX);
            dfMaxY = std::min(dfMaxY, dfThisMaxY);
        }

        bFoundNonMatchingGT = bFoundNonMatchingGT || bIsThisOneNonMatching;

        iSpectralBand++;
    }

    if (nInputSpectralBandsIn && iSpectralBand != nInputSpectralBandsIn)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Less SpectralBand elements than in source array");
        goto error;
    }

    /* -------------------------------------------------------------------- */
    /*      On-the-fly spatial extent adjustment if needed and asked.       */
    /* -------------------------------------------------------------------- */

    if (bFoundNonMatchingGT && (m_eGTAdjustment == GTAdjust_Union ||
                                m_eGTAdjustment == GTAdjust_Intersection))
    {
        if (bFoundRotatingTerms)
        {
            CPLError(
                CE_Failure, CPLE_NotSupported,
                "One of the panchromatic or spectral datasets has rotating "
                "terms in their geotransform matrix. Adjustment not possible");
            goto error;
        }
        if (m_eGTAdjustment == GTAdjust_Intersection &&
            (dfMinX >= dfMaxX || dfMinY >= dfMaxY))
        {
            CPLError(
                CE_Failure, CPLE_NotSupported,
                "One of the panchromatic or spectral datasets has rotating "
                "terms in their geotransform matrix. Adjustment not possible");
            goto error;
        }
        if (m_eGTAdjustment == GTAdjust_Union)
            CPLDebug("VRT", "Do union of bounding box of panchromatic and "
                            "spectral datasets");
        else
            CPLDebug("VRT", "Do intersection of bounding box of panchromatic "
                            "and spectral datasets");

        // If the pandataset needs adjustments, make sure the coordinates of the
        // union/intersection properly align with the grid of the pandataset
        // to avoid annoying sub-pixel shifts on the panchro band.
        double dfPixelSize = std::max(panGT[1], std::abs(panGT[5]));
        if (std::abs(panGT[0] - dfMinX) > dfPixelSize ||
            std::abs(dfPanMaxY - dfMaxY) > dfPixelSize ||
            std::abs(dfLRPanX - dfMaxX) > dfPixelSize ||
            std::abs(dfPanMinY - dfMinY) > dfPixelSize)
        {
            dfMinX =
                panGT[0] +
                std::floor((dfMinX - panGT[0]) / panGT[1] + 0.5) * panGT[1];
            dfMaxX =
                dfLRPanX +
                std::floor((dfMaxX - dfLRPanX) / panGT[1] + 0.5) * panGT[1];
            if (panGT[5] > 0)
            {
                dfMinY = dfPanMinY +
                         std::floor((dfMinY - dfPanMinY) / panGT[5] + 0.5) *
                             panGT[5];
                dfMaxY = dfPanMinY +
                         std::floor((dfMaxY - dfPanMinY) / panGT[5] + 0.5) *
                             panGT[5];
            }
            else
            {
                dfMinY = dfPanMaxY +
                         std::floor((dfMinY - dfPanMaxY) / panGT[5] + 0.5) *
                             panGT[5];
                dfMaxY = dfPanMaxY +
                         std::floor((dfMaxY - dfPanMaxY) / panGT[5] + 0.5) *
                             panGT[5];
            }
        }

        std::map<CPLString, GDALDataset *>::iterator oIter =
            oMapNamesToDataset.begin();
        for (; oIter != oMapNamesToDataset.end(); ++oIter)
        {
            GDALDataset *poSrcDS = oIter->second;
            GDALGeoTransform gt;
            if (poSrcDS->GetGeoTransform(gt) != CE_None)
                continue;

            // Check if this dataset needs adjustments
            dfPixelSize = std::max(gt.xscale, std::abs(gt.yscale));
            dfPixelSize = std::max(panGT[1], dfPixelSize);
            dfPixelSize = std::max(std::abs(panGT[5]), dfPixelSize);
            double dfThisMinX = gt.xorig;
            double dfThisMaxX =
                gt.xorig + poSrcDS->GetRasterXSize() * gt.xscale;
            double dfThisMaxY = gt.yorig;
            double dfThisMinY =
                gt.yorig + poSrcDS->GetRasterYSize() * gt.yscale;
            if (dfThisMinY > dfThisMaxY)
                std::swap(dfThisMinY, dfThisMaxY);
            if (std::abs(dfThisMinX - dfMinX) <= dfPixelSize &&
                std::abs(dfThisMaxY - dfMaxY) <= dfPixelSize &&
                std::abs(dfThisMaxX - dfMaxX) <= dfPixelSize &&
                std::abs(dfThisMinY - dfMinY) <= dfPixelSize)
            {
                continue;
            }

            GDALGeoTransform adjustedGT;
            adjustedGT[0] = dfMinX;
            adjustedGT[1] = gt.xscale;
            adjustedGT[2] = 0;
            adjustedGT[3] = gt.yscale > 0 ? dfMinY : dfMaxY;
            adjustedGT[4] = 0;
            adjustedGT[5] = gt.yscale;
            const double dfAdjustedRasterXSize =
                0.5 + (dfMaxX - dfMinX) / adjustedGT[1];
            const double dfAdjustedRasterYSize =
                0.5 + (dfMaxY - dfMinY) / std::abs(adjustedGT[5]);
            if (!(dfAdjustedRasterXSize >= 1 &&
                  dfAdjustedRasterXSize < INT_MAX &&
                  dfAdjustedRasterYSize >= 1 &&
                  dfAdjustedRasterYSize < INT_MAX))
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Datasets are too disjoint");
                goto error;
            }
            const int nAdjustRasterXSize =
                static_cast<int>(dfAdjustedRasterXSize);
            const int nAdjustRasterYSize =
                static_cast<int>(dfAdjustedRasterYSize);

            VRTDataset *poVDS =
                new VRTDataset(nAdjustRasterXSize, nAdjustRasterYSize);
            poVDS->SetWritable(FALSE);
            poVDS->SetMetadataItem("UNDERLYING_DATASET",
                                   poSrcDS->GetDescription());
            poVDS->SetGeoTransform(adjustedGT);
            poVDS->SetProjection(poPanDataset->GetProjectionRef());

            for (int i = 0; i < poSrcDS->GetRasterCount(); i++)
            {
                GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(i + 1);
                poVDS->AddBand(poSrcBand->GetRasterDataType(), nullptr);
                VRTSourcedRasterBand *poVRTBand =
                    static_cast<VRTSourcedRasterBand *>(
                        poVDS->GetRasterBand(i + 1));

                const char *pszNBITS =
                    poSrcBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
                if (pszNBITS)
                    poVRTBand->SetMetadataItem("NBITS", pszNBITS,
                                               "IMAGE_STRUCTURE");

                VRTSimpleSource *poSimpleSource = new VRTSimpleSource();
                poVRTBand->ConfigureSource(
                    poSimpleSource, poSrcBand, FALSE,
                    static_cast<int>(
                        std::floor((dfMinX - gt.xorig) / gt.xscale + 0.001)),
                    gt.yscale < 0
                        ? static_cast<int>(std::floor(
                              (dfMaxY - dfThisMaxY) / gt.yscale + 0.001))
                        : static_cast<int>(std::floor(
                              (dfMinY - dfThisMinY) / gt.yscale + 0.001)),
                    static_cast<int>(0.5 + (dfMaxX - dfMinX) / gt.xscale),
                    static_cast<int>(0.5 +
                                     (dfMaxY - dfMinY) / std::abs(gt.yscale)),
                    0, 0, nAdjustRasterXSize, nAdjustRasterYSize);

                poVRTBand->AddSource(poSimpleSource);
            }

            oIter->second = poVDS;
            if (poSrcDS == poPanDataset)
            {
                panGT = adjustedGT;
                poPanDataset = poVDS;
                poPanBand = poPanDataset->GetRasterBand(nPanBand);
                nPanXSize = poPanDataset->GetRasterXSize();
                nPanYSize = poPanDataset->GetRasterYSize();
            }

            m_apoDatasetsToReleaseRef.push_back(
                std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>(
                    poVDS));
        }
    }

    if (nRasterXSize == 0 && nRasterYSize == 0)
    {
        nRasterXSize = nPanXSize;
        nRasterYSize = nPanYSize;
    }
    else if (nRasterXSize != nPanXSize || nRasterYSize != nPanYSize)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Inconsistent declared VRT dimensions with panchro dataset");
        goto error;
    }

    /* -------------------------------------------------------------------- */
    /*      Initialize all the general VRT stuff.  This will even           */
    /*      create the VRTPansharpenedRasterBands and initialize them.      */
    /* -------------------------------------------------------------------- */
    eErr = VRTDataset::XMLInit(psTree, pszVRTPathIn);

    if (eErr != CE_None)
    {
        goto error;
    }

    /* -------------------------------------------------------------------- */
    /*      Inherit georeferencing info from panchro band if not defined    */
    /*      in VRT.                                                         */
    /* -------------------------------------------------------------------- */

    {
        GDALGeoTransform outGT;
        if (GetGeoTransform(outGT) != CE_None && GetGCPCount() == 0 &&
            GetSpatialRef() == nullptr)
        {
            if (bPanGeoTransformValid)
            {
                SetGeoTransform(panGT);
            }
            if (poPanSRS)
            {
                SetSpatialRef(poPanSRS);
            }
        }
    }

    /* -------------------------------------------------------------------- */
    /*      Parse rest of PansharpeningOptions                              */
    /* -------------------------------------------------------------------- */
    iSpectralBand = 0;
    for (const CPLXMLNode *psIter = psOptions->psChild; psIter;
         psIter = psIter->psNext)
    {
        if (psIter->eType != CXT_Element ||
            !EQUAL(psIter->pszValue, "SpectralBand"))
            continue;

        GDALDataset *poDataset;
        GDALRasterBand *poBand;

        const char *pszDstBand = CPLGetXMLValue(psIter, "dstBand", nullptr);
        int nDstBand = -1;
        if (pszDstBand != nullptr)
        {
            nDstBand = atoi(pszDstBand);
            if (nDstBand <= 0)
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "SpectralBand.dstBand = '%s' invalid", pszDstBand);
                goto error;
            }
        }

        if (nInputSpectralBandsIn && pahInputSpectralBandsIn != nullptr)
        {
            poBand = GDALRasterBand::FromHandle(
                pahInputSpectralBandsIn[iSpectralBand]);
            poDataset = poBand->GetDataset();
            if (poDataset)
            {
                poDataset = oMapNamesToDataset[CPLSPrintf("%p", poDataset)];
                CPLAssert(poDataset);
                if (poDataset)
                    poBand = poDataset->GetRasterBand(poBand->GetBand());
            }
        }
        else
        {
            osSourceFilename = CPLGetXMLValue(psIter, "SourceFilename", "");
            const bool bRelativeToVRT = CPL_TO_BOOL(atoi(
                CPLGetXMLValue(psIter, "SourceFilename.relativetoVRT", "0")));
            if (bRelativeToVRT)
                osSourceFilename = CPLProjectRelativeFilenameSafe(
                    pszVRTPathIn, osSourceFilename.c_str());
            poDataset = oMapNamesToDataset[osSourceFilename];
            CPLAssert(poDataset);
            const char *pszSourceBand =
                CPLGetXMLValue(psIter, "SourceBand", "1");
            const int nBand = atoi(pszSourceBand);
            poBand = poDataset->GetRasterBand(nBand);
            if (poBand == nullptr)
            {
                CPLError(CE_Failure, CPLE_AppDefined, "%s invalid band of %s",
                         pszSourceBand, osSourceFilename.c_str());
                goto error;
            }
        }

        if (bHasNoData)
        {
            double dfSpectralNoData = poPanBand->GetNoDataValue(&bHasNoData);
            if (bHasNoData && dfSpectralNoData != dfNoData)
                bHasNoData = FALSE;
        }

        ahSpectralBands.push_back(poBand);
        if (nDstBand >= 1)
        {
            if (aMapDstBandToSpectralBand.find(nDstBand - 1) !=
                aMapDstBandToSpectralBand.end())
            {
                CPLError(
                    CE_Failure, CPLE_AppDefined,
                    "Another spectral band is already mapped to output band %d",
                    nDstBand);
                goto error;
            }
            aMapDstBandToSpectralBand[nDstBand - 1] =
                static_cast<int>(ahSpectralBands.size() - 1);
        }

        iSpectralBand++;
    }

    if (ahSpectralBands.empty())
    {
        CPLError(CE_Failure, CPLE_AppDefined, "No spectral band defined");
        goto error;
    }

    {
        const char *pszNoData = CPLGetXMLValue(psOptions, "NoData", nullptr);
        if (pszNoData)
        {
            if (EQUAL(pszNoData, "NONE"))
            {
                m_bNoDataDisabled = TRUE;
                bHasNoData = FALSE;
            }
            else
            {
                bHasNoData = TRUE;
                dfNoData = CPLAtof(pszNoData);
            }
        }
    }

    if (GetRasterCount() == 0)
    {
        for (int i = 0; i < static_cast<int>(aMapDstBandToSpectralBand.size());
             i++)
        {
            if (aMapDstBandToSpectralBand.find(i) ==
                aMapDstBandToSpectralBand.end())
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Hole in SpectralBand.dstBand numbering");
                goto error;
            }
            GDALRasterBand *poInputBand = GDALRasterBand::FromHandle(
                ahSpectralBands[aMapDstBandToSpectralBand[i]]);
            GDALRasterBand *poBand = new VRTPansharpenedRasterBand(
                this, i + 1, poInputBand->GetRasterDataType());
            poBand->SetColorInterpretation(
                poInputBand->GetColorInterpretation());
            if (bHasNoData)
                poBand->SetNoDataValue(dfNoData);
            SetBand(i + 1, poBand);
        }
    }
    else
    {
        int nIdxAsPansharpenedBand = 0;
        for (int i = 0; i < nBands; i++)
        {
            if (static_cast<VRTRasterBand *>(GetRasterBand(i + 1))
                    ->IsPansharpenRasterBand())
            {
                if (aMapDstBandToSpectralBand.find(i) ==
                    aMapDstBandToSpectralBand.end())
                {
                    CPLError(CE_Failure, CPLE_AppDefined,
                             "Band %d of type VRTPansharpenedRasterBand, but "
                             "no corresponding SpectralBand",
                             i + 1);
                    goto error;
                }
                else
                {
                    static_cast<VRTPansharpenedRasterBand *>(
                        GetRasterBand(i + 1))
                        ->SetIndexAsPansharpenedBand(nIdxAsPansharpenedBand);
                    nIdxAsPansharpenedBand++;
                }
            }
        }
    }

    // Figure out bit depth
    {
        const char *pszBitDepth =
            CPLGetXMLValue(psOptions, "BitDepth", nullptr);
        if (pszBitDepth == nullptr)
            pszBitDepth = GDALRasterBand::FromHandle(ahSpectralBands[0])
                              ->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
        if (pszBitDepth)
            nBitDepth = atoi(pszBitDepth);
        if (nBitDepth)
        {
            for (int i = 0; i < nBands; i++)
            {
                if (!static_cast<VRTRasterBand *>(GetRasterBand(i + 1))
                         ->IsPansharpenRasterBand())
                    continue;
                if (GetRasterBand(i + 1)->GetMetadataItem(
                        "NBITS", "IMAGE_STRUCTURE") == nullptr)
                {
                    if (nBitDepth != 8 && nBitDepth != 16 && nBitDepth != 32)
                    {
                        GetRasterBand(i + 1)->SetMetadataItem(
                            "NBITS", CPLSPrintf("%d", nBitDepth),
                            "IMAGE_STRUCTURE");
                    }
                }
                else if (nBitDepth == 8 || nBitDepth == 16 || nBitDepth == 32)
                {
                    GetRasterBand(i + 1)->SetMetadataItem("NBITS", nullptr,
                                                          "IMAGE_STRUCTURE");
                }
            }
        }
    }

    if (GDALGetRasterBandXSize(ahSpectralBands[0]) >
            GDALGetRasterBandXSize(poPanBand) ||
        GDALGetRasterBandYSize(ahSpectralBands[0]) >
            GDALGetRasterBandYSize(poPanBand))
    {
        CPLError(CE_Warning, CPLE_AppDefined,
                 "Dimensions of spectral band larger than panchro band");
    }

    aMapDstBandToSpectralBandIter = aMapDstBandToSpectralBand.begin();
    for (; aMapDstBandToSpectralBandIter != aMapDstBandToSpectralBand.end();
         ++aMapDstBandToSpectralBandIter)
    {
        const int nDstBand = 1 + aMapDstBandToSpectralBandIter->first;
        if (nDstBand > nBands ||
            !static_cast<VRTRasterBand *>(GetRasterBand(nDstBand))
                 ->IsPansharpenRasterBand())
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "SpectralBand.dstBand = '%d' invalid", nDstBand);
            goto error;
        }
    }

    if (adfWeights.empty())
    {
        for (int i = 0; i < static_cast<int>(ahSpectralBands.size()); i++)
        {
            adfWeights.push_back(1.0 / ahSpectralBands.size());
        }
    }
    else if (adfWeights.size() != ahSpectralBands.size())
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "%d weights defined, but %d input spectral bands",
                 static_cast<int>(adfWeights.size()),
                 static_cast<int>(ahSpectralBands.size()));
        goto error;
    }

    if (aMapDstBandToSpectralBand.empty())
    {
        CPLError(CE_Warning, CPLE_AppDefined,
                 "No spectral band is mapped to an output band");
    }

    /* -------------------------------------------------------------------- */
    /*      Instantiate poPansharpener                                      */
    /* -------------------------------------------------------------------- */
    psPanOptions.reset(GDALCreatePansharpenOptions());
    psPanOptions->ePansharpenAlg = GDAL_PSH_WEIGHTED_BROVEY;
    psPanOptions->eResampleAlg = eResampleAlg;
    psPanOptions->nBitDepth = nBitDepth;
    psPanOptions->nWeightCount = static_cast<int>(adfWeights.size());
    psPanOptions->padfWeights =
        static_cast<double *>(CPLMalloc(sizeof(double) * adfWeights.size()));
    memcpy(psPanOptions->padfWeights, &adfWeights[0],
           sizeof(double) * adfWeights.size());
    psPanOptions->hPanchroBand = poPanBand;
    psPanOptions->nInputSpectralBands =
        static_cast<int>(ahSpectralBands.size());
    psPanOptions->pahInputSpectralBands = static_cast<GDALRasterBandH *>(
        CPLMalloc(sizeof(GDALRasterBandH) * ahSpectralBands.size()));
    memcpy(psPanOptions->pahInputSpectralBands, &ahSpectralBands[0],
           sizeof(GDALRasterBandH) * ahSpectralBands.size());
    psPanOptions->nOutPansharpenedBands =
        static_cast<int>(aMapDstBandToSpectralBand.size());
    psPanOptions->panOutPansharpenedBands = static_cast<int *>(
        CPLMalloc(sizeof(int) * aMapDstBandToSpectralBand.size()));
    aMapDstBandToSpectralBandIter = aMapDstBandToSpectralBand.begin();
    for (int i = 0;
         aMapDstBandToSpectralBandIter != aMapDstBandToSpectralBand.end();
         ++aMapDstBandToSpectralBandIter, ++i)
    {
        psPanOptions->panOutPansharpenedBands[i] =
            aMapDstBandToSpectralBandIter->second;
    }
    psPanOptions->bHasNoData = bHasNoData;
    psPanOptions->dfNoData = dfNoData;
    psPanOptions->nThreads = nThreads;

    if (nBands == psPanOptions->nOutPansharpenedBands)
        SetMetadataItem("INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE");

    m_poPansharpener = std::make_unique<GDALPansharpenOperation>();
    eErr = m_poPansharpener->Initialize(psPanOptions.get());
    if (eErr != CE_None)
    {
        // Delete the pansharper object before closing the dataset
        // because it may have warped the bands into an intermediate VRT
        m_poPansharpener.reset();

        // Close in reverse order (VRT firsts and real datasets after)
        for (int i = static_cast<int>(m_apoDatasetsToReleaseRef.size()) - 1;
             i >= 0; i--)
        {
            m_apoDatasetsToReleaseRef[i].reset();
        }
        m_apoDatasetsToReleaseRef.clear();
    }

    return eErr;

error:
    // Close in reverse order (VRT firsts and real datasets after)
    for (int i = static_cast<int>(m_apoDatasetsToReleaseRef.size()) - 1; i >= 0;
         i--)
    {
        m_apoDatasetsToReleaseRef[i].reset();
    }
    m_apoDatasetsToReleaseRef.clear();
    return CE_Failure;
}

/************************************************************************/
/*                           SerializeToXML()                           */
/************************************************************************/

CPLXMLNode *VRTPansharpenedDataset::SerializeToXML(const char *pszVRTPathIn)

{
    CPLXMLNode *psTree = VRTDataset::SerializeToXML(pszVRTPathIn);

    if (psTree == nullptr)
        return psTree;

    /* -------------------------------------------------------------------- */
    /*      Set subclass.                                                   */
    /* -------------------------------------------------------------------- */
    CPLCreateXMLNode(CPLCreateXMLNode(psTree, CXT_Attribute, "subClass"),
                     CXT_Text, "VRTPansharpenedDataset");

    /* -------------------------------------------------------------------- */
    /*      Serialize the block size.                                       */
    /* -------------------------------------------------------------------- */
    CPLCreateXMLElementAndValue(psTree, "BlockXSize",
                                CPLSPrintf("%d", m_nBlockXSize));
    CPLCreateXMLElementAndValue(psTree, "BlockYSize",
                                CPLSPrintf("%d", m_nBlockYSize));

    /* -------------------------------------------------------------------- */
    /*      Serialize the options.                                          */
    /* -------------------------------------------------------------------- */
    if (m_poPansharpener == nullptr)
        return psTree;
    const GDALPansharpenOptions *psOptions = m_poPansharpener->GetOptions();
    if (psOptions == nullptr)
        return psTree;

    CPLXMLNode *psOptionsNode =
        CPLCreateXMLNode(psTree, CXT_Element, "PansharpeningOptions");

    if (psOptions->ePansharpenAlg == GDAL_PSH_WEIGHTED_BROVEY)
    {
        CPLCreateXMLElementAndValue(psOptionsNode, "Algorithm",
                                    "WeightedBrovey");
    }
    else
    {
        CPLAssert(false);
    }
    if (psOptions->nWeightCount)
    {
        CPLString osWeights;
        for (int i = 0; i < psOptions->nWeightCount; i++)
        {
            if (i)
                osWeights += ",";
            osWeights += CPLSPrintf("%.16g", psOptions->padfWeights[i]);
        }
        CPLCreateXMLElementAndValue(
            CPLCreateXMLNode(psOptionsNode, CXT_Element, "AlgorithmOptions"),
            "Weights", osWeights.c_str());
    }
    CPLCreateXMLElementAndValue(
        psOptionsNode, "Resampling",
        GDALRasterIOGetResampleAlg(psOptions->eResampleAlg));

    if (psOptions->nThreads == -1)
    {
        CPLCreateXMLElementAndValue(psOptionsNode, "NumThreads", "ALL_CPUS");
    }
    else if (psOptions->nThreads > 1)
    {
        CPLCreateXMLElementAndValue(psOptionsNode, "NumThreads",
                                    CPLSPrintf("%d", psOptions->nThreads));
    }

    if (psOptions->nBitDepth)
        CPLCreateXMLElementAndValue(psOptionsNode, "BitDepth",
                                    CPLSPrintf("%d", psOptions->nBitDepth));

    const char *pszAdjust = nullptr;
    switch (m_eGTAdjustment)
    {
        case GTAdjust_Union:
            pszAdjust = "Union";
            break;
        case GTAdjust_Intersection:
            pszAdjust = "Intersection";
            break;
        case GTAdjust_None:
            pszAdjust = "None";
            break;
        case GTAdjust_NoneWithoutWarning:
            pszAdjust = "NoneWithoutWarning";
            break;
        default:
            break;
    }

    if (psOptions->bHasNoData)
    {
        CPLCreateXMLElementAndValue(psOptionsNode, "NoData",
                                    CPLSPrintf("%.16g", psOptions->dfNoData));
    }
    else if (m_bNoDataDisabled)
    {
        CPLCreateXMLElementAndValue(psOptionsNode, "NoData", "None");
    }

    if (pszAdjust)
        CPLCreateXMLElementAndValue(psOptionsNode, "SpatialExtentAdjustment",
                                    pszAdjust);

    if (psOptions->hPanchroBand)
    {
        CPLXMLNode *psBand =
            CPLCreateXMLNode(psOptionsNode, CXT_Element, "PanchroBand");
        GDALRasterBand *poBand =
            GDALRasterBand::FromHandle(psOptions->hPanchroBand);
        if (poBand->GetDataset())
        {
            auto poPanchroDS = poBand->GetDataset();
            std::map<CPLString, CPLString>::iterator oIter =
                m_oMapToRelativeFilenames.find(poPanchroDS->GetDescription());
            if (oIter == m_oMapToRelativeFilenames.end())
            {
                const char *pszUnderlyingDS =
                    poPanchroDS->GetMetadataItem("UNDERLYING_DATASET");
                if (!pszUnderlyingDS)
                    pszUnderlyingDS = poPanchroDS->GetDescription();
                CPLCreateXMLElementAndValue(psBand, "SourceFilename",
                                            pszUnderlyingDS);
            }
            else
            {
                CPLXMLNode *psSourceFilename = CPLCreateXMLElementAndValue(
                    psBand, "SourceFilename", oIter->second);
                CPLCreateXMLNode(CPLCreateXMLNode(psSourceFilename,
                                                  CXT_Attribute,
                                                  "relativeToVRT"),
                                 CXT_Text, "1");
            }

            GDALSerializeOpenOptionsToXML(psBand,
                                          poPanchroDS->GetOpenOptions());

            CPLCreateXMLElementAndValue(psBand, "SourceBand",
                                        CPLSPrintf("%d", poBand->GetBand()));
        }
    }
    for (int i = 0; i < psOptions->nInputSpectralBands; i++)
    {
        CPLXMLNode *psBand =
            CPLCreateXMLNode(psOptionsNode, CXT_Element, "SpectralBand");

        for (int j = 0; j < psOptions->nOutPansharpenedBands; j++)
        {
            if (psOptions->panOutPansharpenedBands[j] == i)
            {
                for (int k = 0; k < nBands; k++)
                {
                    if (static_cast<VRTRasterBand *>(GetRasterBand(k + 1))
                            ->IsPansharpenRasterBand())
                    {
                        if (static_cast<VRTPansharpenedRasterBand *>(
                                GetRasterBand(k + 1))
                                ->GetIndexAsPansharpenedBand() == j)
                        {
                            CPLCreateXMLNode(CPLCreateXMLNode(psBand,
                                                              CXT_Attribute,
                                                              "dstBand"),
                                             CXT_Text, CPLSPrintf("%d", k + 1));
                            break;
                        }
                    }
                }
                break;
            }
        }

        GDALRasterBand *poBand =
            GDALRasterBand::FromHandle(psOptions->pahInputSpectralBands[i]);
        if (poBand->GetDataset())
        {
            auto poSpectralDS = poBand->GetDataset();
            std::map<CPLString, CPLString>::iterator oIter =
                m_oMapToRelativeFilenames.find(poSpectralDS->GetDescription());
            if (oIter == m_oMapToRelativeFilenames.end())
            {
                const char *pszUnderlyingDS =
                    poSpectralDS->GetMetadataItem("UNDERLYING_DATASET");
                if (!pszUnderlyingDS)
                    pszUnderlyingDS = poSpectralDS->GetDescription();
                CPLCreateXMLElementAndValue(psBand, "SourceFilename",
                                            pszUnderlyingDS);
            }
            else
            {
                CPLXMLNode *psSourceFilename = CPLCreateXMLElementAndValue(
                    psBand, "SourceFilename", oIter->second);
                CPLCreateXMLNode(CPLCreateXMLNode(psSourceFilename,
                                                  CXT_Attribute,
                                                  "relativeToVRT"),
                                 CXT_Text, "1");
            }

            GDALSerializeOpenOptionsToXML(psBand,
                                          poSpectralDS->GetOpenOptions());

            CPLCreateXMLElementAndValue(psBand, "SourceBand",
                                        CPLSPrintf("%d", poBand->GetBand()));
        }
    }

    return psTree;
}

/************************************************************************/
/*                            GetBlockSize()                            */
/************************************************************************/

void VRTPansharpenedDataset::GetBlockSize(int *pnBlockXSize,
                                          int *pnBlockYSize) const

{
    assert(nullptr != pnBlockXSize);
    assert(nullptr != pnBlockYSize);

    *pnBlockXSize = m_nBlockXSize;
    *pnBlockYSize = m_nBlockYSize;
}

/************************************************************************/
/*                              AddBand()                               */
/************************************************************************/

CPLErr VRTPansharpenedDataset::AddBand(CPL_UNUSED GDALDataType eType,
                                       CPL_UNUSED CSLConstList papszOptions)

{
    CPLError(CE_Failure, CPLE_NotSupported, "AddBand() not supported");

    return CE_Failure;
}

/************************************************************************/
/*                             IRasterIO()                              */
/************************************************************************/

CPLErr VRTPansharpenedDataset::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)
{
    if (eRWFlag == GF_Write)
        return CE_Failure;

    /* Try to pass the request to the most appropriate overview dataset */
    if (nBufXSize < nXSize && nBufYSize < nYSize)
    {
        int bTried;
        CPLErr eErr = TryOverviewRasterIO(
            eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
            eBufType, nBandCount, panBandMap, nPixelSpace, nLineSpace,
            nBandSpace, psExtraArg, &bTried);
        if (bTried)
            return eErr;
    }

    const int nDataTypeSize = GDALGetDataTypeSizeBytes(eBufType);
    if (nXSize == nBufXSize && nYSize == nBufYSize &&
        nDataTypeSize == nPixelSpace && nLineSpace == nPixelSpace * nBufXSize &&
        nBandSpace == nLineSpace * nBufYSize && nBandCount == nBands)
    {
        for (int i = 0; i < nBands; i++)
        {
            if (panBandMap[i] != i + 1 ||
                !static_cast<VRTRasterBand *>(GetRasterBand(i + 1))
                     ->IsPansharpenRasterBand())
            {
                goto default_path;
            }
        }

        //{static int bDone = 0; if (!bDone) printf("(2)\n"); bDone = 1; }
        return m_poPansharpener->ProcessRegion(nXOff, nYOff, nXSize, nYSize,
                                               pData, eBufType);
    }

default_path:
    //{static int bDone = 0; if (!bDone) printf("(3)\n"); bDone = 1; }
    return VRTDataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
                                 nBufXSize, nBufYSize, eBufType, nBandCount,
                                 panBandMap, nPixelSpace, nLineSpace,
                                 nBandSpace, psExtraArg);
}

/************************************************************************/
/* ==================================================================== */
/*                        VRTPansharpenedRasterBand                     */
/* ==================================================================== */
/************************************************************************/

/************************************************************************/
/*                     VRTPansharpenedRasterBand()                      */
/************************************************************************/

VRTPansharpenedRasterBand::VRTPansharpenedRasterBand(GDALDataset *poDSIn,
                                                     int nBandIn,
                                                     GDALDataType eDataTypeIn)
    : m_nIndexAsPansharpenedBand(nBandIn - 1)
{
    Initialize(poDSIn->GetRasterXSize(), poDSIn->GetRasterYSize());

    poDS = poDSIn;
    nBand = nBandIn;
    eAccess = GA_Update;
    eDataType = eDataTypeIn;

    static_cast<VRTPansharpenedDataset *>(poDS)->GetBlockSize(&nBlockXSize,
                                                              &nBlockYSize);
}

/************************************************************************/
/*                     ~VRTPansharpenedRasterBand()                     */
/************************************************************************/

VRTPansharpenedRasterBand::~VRTPansharpenedRasterBand()

{
    FlushCache(true);
}

/************************************************************************/
/*                             IReadBlock()                             */
/************************************************************************/

CPLErr VRTPansharpenedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
                                             void *pImage)

{
    const int nReqXOff = nBlockXOff * nBlockXSize;
    const int nReqYOff = nBlockYOff * nBlockYSize;
    const int nReqXSize = std::min(nBlockXSize, nRasterXSize - nReqXOff);
    const int nReqYSize = std::min(nBlockYSize, nRasterYSize - nReqYOff);

    //{static int bDone = 0; if (!bDone) printf("(4)\n"); bDone = 1; }
    const int nDataTypeSize = GDALGetDataTypeSizeBytes(eDataType);
    GDALRasterIOExtraArg sExtraArg;
    INIT_RASTERIO_EXTRA_ARG(sExtraArg);
    if (IRasterIO(GF_Read, nReqXOff, nReqYOff, nReqXSize, nReqYSize, pImage,
                  nReqXSize, nReqYSize, eDataType, nDataTypeSize,
                  static_cast<GSpacing>(nDataTypeSize) * nReqXSize,
                  &sExtraArg) != CE_None)
    {
        return CE_Failure;
    }

    if (nReqXSize < nBlockXSize)
    {
        for (int j = nReqYSize - 1; j >= 0; j--)
        {
            memmove(static_cast<GByte *>(pImage) +
                        static_cast<size_t>(j) * nDataTypeSize * nBlockXSize,
                    static_cast<GByte *>(pImage) +
                        static_cast<size_t>(j) * nDataTypeSize * nReqXSize,
                    static_cast<size_t>(nReqXSize) * nDataTypeSize);
            memset(static_cast<GByte *>(pImage) +
                       (static_cast<size_t>(j) * nBlockXSize + nReqXSize) *
                           nDataTypeSize,
                   0,
                   static_cast<size_t>(nBlockXSize - nReqXSize) *
                       nDataTypeSize);
        }
    }
    if (nReqYSize < nBlockYSize)
    {
        memset(static_cast<GByte *>(pImage) +
                   static_cast<size_t>(nReqYSize) * nBlockXSize * nDataTypeSize,
               0,
               static_cast<size_t>(nBlockYSize - nReqYSize) * nBlockXSize *
                   nDataTypeSize);
    }

    // Cache other bands
    CPLErr eErr = CE_None;
    VRTPansharpenedDataset *poGDS = static_cast<VRTPansharpenedDataset *>(poDS);
    if (poGDS->nBands != 1 && !poGDS->m_bLoadingOtherBands)
    {
        poGDS->m_bLoadingOtherBands = TRUE;

        for (int iOtherBand = 1; iOtherBand <= poGDS->nBands; iOtherBand++)
        {
            if (iOtherBand == nBand)
                continue;

            GDALRasterBlock *poBlock =
                poGDS->GetRasterBand(iOtherBand)
                    ->GetLockedBlockRef(nBlockXOff, nBlockYOff);
            if (poBlock == nullptr)
            {
                eErr = CE_Failure;
                break;
            }
            poBlock->DropLock();
        }

        poGDS->m_bLoadingOtherBands = FALSE;
    }

    return eErr;
}

/************************************************************************/
/*                             IRasterIO()                              */
/************************************************************************/

CPLErr VRTPansharpenedRasterBand::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)
{
    if (eRWFlag == GF_Write)
        return CE_Failure;

    VRTPansharpenedDataset *poGDS = static_cast<VRTPansharpenedDataset *>(poDS);

    /* Try to pass the request to the most appropriate overview dataset */
    if (nBufXSize < nXSize && nBufYSize < nYSize)
    {
        int bTried;
        CPLErr eErr = TryOverviewRasterIO(
            eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
            eBufType, nPixelSpace, nLineSpace, psExtraArg, &bTried);
        if (bTried)
            return eErr;
    }

    const int nDataTypeSize = GDALGetDataTypeSizeBytes(eBufType);
    if (nDataTypeSize > 0 && nXSize == nBufXSize && nYSize == nBufYSize &&
        nDataTypeSize == nPixelSpace && nLineSpace == nPixelSpace * nBufXSize)
    {
        const GDALPansharpenOptions *psOptions =
            poGDS->m_poPansharpener->GetOptions();

        // Have we already done this request for another band ?
        // If so use the cached result
        const size_t nBufferSizePerBand =
            static_cast<size_t>(nXSize) * nYSize * nDataTypeSize;
        if (nXOff == poGDS->m_nLastBandRasterIOXOff &&
            nYOff >= poGDS->m_nLastBandRasterIOYOff &&
            nXSize == poGDS->m_nLastBandRasterIOXSize &&
            nYOff + nYSize <= poGDS->m_nLastBandRasterIOYOff +
                                  poGDS->m_nLastBandRasterIOYSize &&
            eBufType == poGDS->m_eLastBandRasterIODataType)
        {
            //{static int bDone = 0; if (!bDone) printf("(6)\n"); bDone = 1; }
            if (poGDS->m_pabyLastBufferBandRasterIO == nullptr)
                return CE_Failure;
            const size_t nBufferSizePerBandCached =
                static_cast<size_t>(nXSize) * poGDS->m_nLastBandRasterIOYSize *
                nDataTypeSize;
            memcpy(pData,
                   poGDS->m_pabyLastBufferBandRasterIO +
                       nBufferSizePerBandCached * m_nIndexAsPansharpenedBand +
                       static_cast<size_t>(nYOff -
                                           poGDS->m_nLastBandRasterIOYOff) *
                           nXSize * nDataTypeSize,
                   nBufferSizePerBand);
            return CE_None;
        }

        int nYSizeToCache = nYSize;
        if (nYSize == 1 && nXSize == nRasterXSize)
        {
            //{static int bDone = 0; if (!bDone) printf("(7)\n"); bDone = 1; }
            // For efficiency, try to cache at leak 256 K
            nYSizeToCache = (256 * 1024) / nXSize / nDataTypeSize;
            if (nYSizeToCache == 0)
                nYSizeToCache = 1;
            else if (nYOff + nYSizeToCache > nRasterYSize)
                nYSizeToCache = nRasterYSize - nYOff;
        }
        const GUIntBig nBufferSize = static_cast<GUIntBig>(nXSize) *
                                     nYSizeToCache * nDataTypeSize *
                                     psOptions->nOutPansharpenedBands;
        // Check the we don't overflow (for 32 bit platforms)
        if (nBufferSize > std::numeric_limits<size_t>::max() / 2)
        {
            CPLError(CE_Failure, CPLE_OutOfMemory,
                     "Out of memory error while allocating working buffers");
            return CE_Failure;
        }
        GByte *pabyTemp = static_cast<GByte *>(
            VSI_REALLOC_VERBOSE(poGDS->m_pabyLastBufferBandRasterIO,
                                static_cast<size_t>(nBufferSize)));
        if (pabyTemp == nullptr)
        {
            return CE_Failure;
        }
        poGDS->m_nLastBandRasterIOXOff = nXOff;
        poGDS->m_nLastBandRasterIOYOff = nYOff;
        poGDS->m_nLastBandRasterIOXSize = nXSize;
        poGDS->m_nLastBandRasterIOYSize = nYSizeToCache;
        poGDS->m_eLastBandRasterIODataType = eBufType;
        poGDS->m_pabyLastBufferBandRasterIO = pabyTemp;

        CPLErr eErr = poGDS->m_poPansharpener->ProcessRegion(
            nXOff, nYOff, nXSize, nYSizeToCache,
            poGDS->m_pabyLastBufferBandRasterIO, eBufType);
        if (eErr == CE_None)
        {
            //{static int bDone = 0; if (!bDone) printf("(8)\n"); bDone = 1; }
            size_t nBufferSizePerBandCached = static_cast<size_t>(nXSize) *
                                              poGDS->m_nLastBandRasterIOYSize *
                                              nDataTypeSize;
            memcpy(pData,
                   poGDS->m_pabyLastBufferBandRasterIO +
                       nBufferSizePerBandCached * m_nIndexAsPansharpenedBand,
                   nBufferSizePerBand);
        }
        else
        {
            VSIFree(poGDS->m_pabyLastBufferBandRasterIO);
            poGDS->m_pabyLastBufferBandRasterIO = nullptr;
        }
        return eErr;
    }

    //{static int bDone = 0; if (!bDone) printf("(9)\n"); bDone = 1; }
    CPLErr eErr = VRTRasterBand::IRasterIO(
        eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
        eBufType, nPixelSpace, nLineSpace, psExtraArg);

    return eErr;
}

/************************************************************************/
/*                           SerializeToXML()                           */
/************************************************************************/

CPLXMLNode *
VRTPansharpenedRasterBand::SerializeToXML(const char *pszVRTPathIn,
                                          bool &bHasWarnedAboutRAMUsage,
                                          size_t &nAccRAMUsage)

{
    CPLXMLNode *psTree = VRTRasterBand::SerializeToXML(
        pszVRTPathIn, bHasWarnedAboutRAMUsage, nAccRAMUsage);

    /* -------------------------------------------------------------------- */
    /*      Set subclass.                                                   */
    /* -------------------------------------------------------------------- */
    CPLCreateXMLNode(CPLCreateXMLNode(psTree, CXT_Attribute, "subClass"),
                     CXT_Text, "VRTPansharpenedRasterBand");

    return psTree;
}

/************************************************************************/
/*                          GetOverviewCount()                          */
/************************************************************************/

int VRTPansharpenedRasterBand::GetOverviewCount()
{
    VRTPansharpenedDataset *poGDS = static_cast<VRTPansharpenedDataset *>(poDS);

    // Build on-the-fly overviews from overviews of pan and spectral bands
    if (poGDS->m_poPansharpener != nullptr &&
        poGDS->m_apoOverviewDatasets.empty() && poGDS->m_poMainDataset == poGDS)
    {
        const GDALPansharpenOptions *psOptions =
            poGDS->m_poPansharpener->GetOptions();

        GDALRasterBand *poPanBand =
            GDALRasterBand::FromHandle(psOptions->hPanchroBand);
        int nOvrCount = poPanBand->GetOverviewCount();
        if (nOvrCount > 0)
        {
            for (int i = 0; i < poGDS->GetRasterCount(); i++)
            {
                if (!static_cast<VRTRasterBand *>(poGDS->GetRasterBand(i + 1))
                         ->IsPansharpenRasterBand())
                {
                    return 0;
                }
            }

            // Limit number of overviews of the VRTPansharpenedRasterBand to
            // the minimum number of overviews of the pan and spectral source
            // bands, and also make sure all spectral bands have overviews
            // of the same dimension for a given level.
            std::vector<std::pair<int, int>> sizeSpectralOverviews;
            for (int i = 0; i < psOptions->nInputSpectralBands; i++)
            {
                auto poSpectralBand = GDALRasterBand::FromHandle(
                    psOptions->pahInputSpectralBands[i]);
                nOvrCount =
                    std::min(nOvrCount, poSpectralBand->GetOverviewCount());
                if (i == 0)
                {
                    for (int iOvr = 0; iOvr < nOvrCount; ++iOvr)
                    {
                        auto poOvrBand = poSpectralBand->GetOverview(iOvr);
                        sizeSpectralOverviews.emplace_back(
                            poOvrBand->GetXSize(), poOvrBand->GetYSize());
                    }
                }
                else
                {
                    for (int iOvr = 0; iOvr < nOvrCount; ++iOvr)
                    {
                        auto poOvrBand = poSpectralBand->GetOverview(iOvr);
                        if (sizeSpectralOverviews[iOvr].first !=
                                poOvrBand->GetXSize() ||
                            sizeSpectralOverviews[iOvr].second !=
                                poOvrBand->GetYSize())
                        {
                            nOvrCount = iOvr;
                            break;
                        }
                    }
                }
            }

            auto poPanBandDS = poPanBand->GetDataset();
            for (int j = 0; j < nOvrCount; j++)
            {
                std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>
                    poPanOvrDS(GDALCreateOverviewDataset(poPanBandDS, j, true));
                if (!poPanOvrDS)
                {
                    CPLError(CE_Warning, CPLE_AppDefined,
                             "GDALCreateOverviewDataset(poPanBandDS, %d, true) "
                             "failed",
                             j);
                    break;
                }
                GDALRasterBand *poPanOvrBand =
                    poPanOvrDS->GetRasterBand(poPanBand->GetBand());
                auto poOvrDS = std::make_unique<VRTPansharpenedDataset>(
                    poPanOvrBand->GetXSize(), poPanOvrBand->GetYSize());
                poOvrDS->m_apoDatasetsToReleaseRef.push_back(
                    std::move(poPanOvrDS));
                for (int i = 0; i < poGDS->GetRasterCount(); i++)
                {
                    GDALRasterBand *poSrcBand = poGDS->GetRasterBand(i + 1);
                    auto poBand = std::make_unique<VRTPansharpenedRasterBand>(
                        poOvrDS.get(), i + 1, poSrcBand->GetRasterDataType());
                    const char *pszNBITS =
                        poSrcBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
                    if (pszNBITS)
                        poBand->SetMetadataItem("NBITS", pszNBITS,
                                                "IMAGE_STRUCTURE");
                    int bHasNoData = FALSE;
                    const double dfNoData =
                        poSrcBand->GetNoDataValue(&bHasNoData);
                    if (bHasNoData)
                        poBand->SetNoDataValue(dfNoData);
                    poBand->SetColorInterpretation(
                        poSrcBand->GetColorInterpretation());
                    poOvrDS->SetBand(i + 1, std::move(poBand));
                }

                std::unique_ptr<GDALPansharpenOptions,
                                decltype(&GDALDestroyPansharpenOptions)>
                    psPanOvrOptions(GDALClonePansharpenOptions(psOptions),
                                    GDALDestroyPansharpenOptions);
                psPanOvrOptions->hPanchroBand = poPanOvrBand;
                for (int i = 0; i < psOptions->nInputSpectralBands; i++)
                {
                    auto poSpectralBand = GDALRasterBand::FromHandle(
                        psOptions->pahInputSpectralBands[i]);
                    std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>
                        poSpectralOvrDS(GDALCreateOverviewDataset(
                            poSpectralBand->GetDataset(), j, true));
                    if (!poSpectralOvrDS)
                    {
                        CPLError(CE_Warning, CPLE_AppDefined,
                                 "GDALCreateOverviewDataset(poSpectralBand->"
                                 "GetDataset(), %d, true) failed",
                                 j);
                        return 0;
                    }
                    psPanOvrOptions->pahInputSpectralBands[i] =
                        poSpectralOvrDS->GetRasterBand(
                            poSpectralBand->GetBand());
                    poOvrDS->m_apoDatasetsToReleaseRef.push_back(
                        std::move(poSpectralOvrDS));
                }
                poOvrDS->m_poPansharpener =
                    std::make_unique<GDALPansharpenOperation>();
                if (poOvrDS->m_poPansharpener->Initialize(
                        psPanOvrOptions.get()) != CE_None)
                {
                    CPLError(CE_Warning, CPLE_AppDefined,
                             "Unable to initialize pansharpener.");
                    return 0;
                }

                poOvrDS->m_poMainDataset = poGDS;
                poOvrDS->SetMetadataItem("INTERLEAVE", "PIXEL",
                                         "IMAGE_STRUCTURE");

                poGDS->m_apoOverviewDatasets.push_back(std::move(poOvrDS));
            }
        }
    }
    return static_cast<int>(poGDS->m_apoOverviewDatasets.size());
}

/************************************************************************/
/*                            GetOverview()                             */
/************************************************************************/

GDALRasterBand *VRTPansharpenedRasterBand::GetOverview(int iOvr)
{
    if (iOvr < 0 || iOvr >= GetOverviewCount())
        return nullptr;

    VRTPansharpenedDataset *poGDS = static_cast<VRTPansharpenedDataset *>(poDS);

    return poGDS->m_apoOverviewDatasets[iOvr]->GetRasterBand(nBand);
}

/*! @endcond */

Coverage Report

Created: 2026-02-14 06:52