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

 *

 * Project:  GDAL

 * Purpose:  gdal "mdim mosaic" subcommand

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

 *

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

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

 *

 * SPDX-License-Identifier: MIT

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

#include "gdalalg_mdim_mosaic.h"

#include "cpl_conv.h"

#include "cpl_vsi_virtual.h"

#include "gdal_priv.h"

#include "vrtdataset.h"

#include <algorithm>

#include <cmath>

#include <optional>

//! @cond Doxygen_Suppress

#ifndef _

#define _(x) (x)

#endif

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

/*          GDALMdimMosaicAlgorithm::GDALMdimMosaicAlgorithm()          */

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

GDALMdimMosaicAlgorithm::GDALMdimMosaicAlgorithm()

    : GDALAlgorithm(NAME, DESCRIPTION, HELP_URL)

{

    AddProgressArg();

    AddOutputFormatArg(&m_outputFormat)

        .AddMetadataItem(GAAMDI_REQUIRED_CAPABILITIES,

                         {GDAL_DCAP_CREATE_MULTIDIMENSIONAL});

    AddOpenOptionsArg(&m_openOptions);

    AddInputFormatsArg(&m_inputFormats)

        .AddMetadataItem(GAAMDI_REQUIRED_CAPABILITIES,

                         {GDAL_ALG_DCAP_RASTER_OR_MULTIDIM_RASTER});

    AddInputDatasetArg(&m_inputDatasets, GDAL_OF_MULTIDIM_RASTER)

        .SetDatasetInputFlags(GADV_NAME)

        .SetDatasetOutputFlags(0)

        .SetAutoOpenDataset(false)

        .SetMinCount(1);

    AddOutputDatasetArg(&m_outputDataset, GDAL_OF_MULTIDIM_RASTER);

    AddCreationOptionsArg(&m_creationOptions);

    AddOverwriteArg(&m_overwrite);

    AddArrayNameArg(&m_array, _("Name of array(s) to mosaic."));

}

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

/*                          GetDimensionDesc()                          */

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

std::optional<GDALMdimMosaicAlgorithm::DimensionDesc>

GDALMdimMosaicAlgorithm::GetDimensionDesc(

    const std::string &osDSName,

    const std::shared_ptr<GDALDimension> &poDim) const

{

    auto poVar = poDim->GetIndexingVariable();

    if (poVar)

    {

        if (poVar->GetDimensionCount() != 1)

        {

            ReportError(

                CE_Failure, CPLE_AppDefined,

                "Dataset %s: indexing variable %s of dimension %s is not 1D",

                osDSName.c_str(), poVar->GetName().c_str(),

                poDim->GetName().c_str());

            return std::nullopt;

        }

        if (poVar->GetDataType().GetClass() != GEDTC_NUMERIC)

        {

            ReportError(

                CE_Failure, CPLE_AppDefined,

                "Dataset %s: indexing variable %s of dimension %s has a "

                "non-numeric data type",

                osDSName.c_str(), poVar->GetName().c_str(),

                poDim->GetName().c_str());

            return std::nullopt;

        }

    }

    DimensionDesc desc;

    desc.osName = poDim->GetName();

    desc.osType = poDim->GetType();

    desc.osDirection = poDim->GetDirection();

    const auto nSize = poDim->GetSize();

    if (poVar)

        desc.attributes = poVar->GetAttributes();

    CPLAssert(nSize > 0);

    desc.nSize = nSize;

    desc.bHasIndexingVar = poVar != nullptr;

    if (!poVar)

    {

        desc.dfStart = 0;

        desc.dfIncrement = 1;

    }

    else if (nSize <= 2 ||

             !poVar->IsRegularlySpaced(desc.dfStart, desc.dfIncrement))

    {

        constexpr uint64_t LIMIT = 100 * 1000 * 1000;

        if (nSize > LIMIT)

        {

            ReportError(CE_Failure, CPLE_AppDefined,

                        "Dataset %s: indexing variable %s of dimension %s has "

                        "too large size",

                        osDSName.c_str(), poVar->GetName().c_str(),

                        desc.osName.c_str());

            return std::nullopt;

        }

        std::vector<double> adfValues(static_cast<size_t>(nSize));

        GUInt64 arrayStartIdx[] = {0};

        size_t anCount[] = {adfValues.size()};

        GInt64 arrayStep[] = {1};

        GPtrDiff_t bufferStride[] = {1};

        if (!poVar->Read(arrayStartIdx, anCount, arrayStep, bufferStride,

                         GDALExtendedDataType::Create(GDT_Float64),

                         adfValues.data()))

        {

            return std::nullopt;

        }

        if (std::isnan(adfValues[0]))

        {

            ReportError(CE_Failure, CPLE_AppDefined,

                        "Dataset %s: indexing variable %s of dimension %s has "

                        "NaN values",

                        osDSName.c_str(), poVar->GetName().c_str(),

                        desc.osName.c_str());

            return std::nullopt;

        }

        if (nSize > 1)

        {

            const int nSign = adfValues[1] > adfValues[0] ? 1 : -1;

            if (std::isnan(adfValues[1]))

            {

                ReportError(CE_Failure, CPLE_AppDefined,

                            "Dataset %s: indexing variable %s of dimension %s "

                            "has NaN values",

                            osDSName.c_str(), poVar->GetName().c_str(),

                            desc.osName.c_str());

                return std::nullopt;

            }

            for (size_t i = 2; i < nSize; ++i)

            {

                if (std::isnan(adfValues[i]))

                {

                    ReportError(CE_Failure, CPLE_AppDefined,

                                "Dataset %s: indexing variable %s of dimension "

                                "%s has NaN values",

                                osDSName.c_str(), poVar->GetName().c_str(),

                                desc.osName.c_str());

                    return std::nullopt;

                }

                const int nSign2 = adfValues[i] > adfValues[i - 1] ? 1 : -1;

                if (nSign != nSign2)

                {

                    ReportError(CE_Failure, CPLE_AppDefined,

                                "Dataset %s: indexing variable %s of dimension "

                                "%s is not strictly increasing or decreasing",

                                osDSName.c_str(), poVar->GetName().c_str(),

                                desc.osName.c_str());

                    return std::nullopt;

                }

            }

            desc.nProgressionSign = nSign;

        }

        std::sort(adfValues.begin(), adfValues.end());

        desc.aaValues.push_back(std::move(adfValues));

    }

    return desc;

}

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

/*           GDALMdimMosaicAlgorithm::BuildArrayParameters()            */

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

bool GDALMdimMosaicAlgorithm::BuildArrayParameters(

    const CPLStringList &aosInputDatasetNames,

    std::vector<ArrayParameters> &aoArrayParameters)

{

    for (const char *pszDatasetName : cpl::Iterate(aosInputDatasetNames))

    {

        auto poDS = std::unique_ptr<GDALDataset>(GDALDataset::Open(

            pszDatasetName, GDAL_OF_MULTIDIM_RASTER | GDAL_OF_VERBOSE_ERROR,

            CPLStringList(m_inputFormats).List(),

            CPLStringList(m_openOptions).List(), nullptr));

        if (!poDS)

            return false;

        auto poRG = poDS->GetRootGroup();

        if (!poRG)

        {

            ReportError(CE_Failure, CPLE_AppDefined,

                        "Cannot get root group for dataset %s", pszDatasetName);

            return false;

        }

        std::vector<std::shared_ptr<GDALMDArray>> apoArrays;

        if (!m_array.empty())

        {

            for (const auto &array : m_array)

            {

                auto poArray = poRG->OpenMDArrayFromFullname(array);

                if (!poArray)

                {

                    ReportError(CE_Failure, CPLE_AppDefined,

                                "Cannot find array %s in dataset %s",

                                array.c_str(), pszDatasetName);

                    return false;

                }

                apoArrays.push_back(std::move(poArray));

            }

        }

        else

        {

            for (const std::string &arrayName :

                 poRG->GetMDArrayFullNamesRecursive())

            {

                auto poArray = poRG->OpenMDArrayFromFullname(arrayName);

                if (!poArray)

                {

                    ReportError(CE_Failure, CPLE_AppDefined,

                                "Cannot open array %s of dataset %s",

                                arrayName.c_str(), pszDatasetName);

                    return false;

                }

                if (poArray->GetDimensionCount() < 2)

                    continue;

                m_array.push_back(arrayName);

                apoArrays.push_back(std::move(poArray));

            }

            if (apoArrays.empty())

            {

                ReportError(

                    CE_Failure, CPLE_AppDefined,

                    "No array of dimension count >= 2 found in dataset %s",

                    pszDatasetName);

                return false;

            }

        }

        if (aoArrayParameters.empty())

            aoArrayParameters.resize(apoArrays.size());

        for (size_t iArray = 0; iArray < apoArrays.size(); ++iArray)

        {

            auto &arrayParameters = aoArrayParameters[iArray];

            auto &poArray = apoArrays[iArray];

            if (poArray->GetDimensionCount() == 0)

            {

                ReportError(CE_Failure, CPLE_AppDefined,

                            "Array %s of dataset %s has no dimensions",

                            poArray->GetName().c_str(), pszDatasetName);

                return false;

            }

            std::vector<SourceShortDimDesc> aoSourceShortDimDesc;

            const auto AddToSourceShortDimDesc =

                [&aoSourceShortDimDesc](const DimensionDesc &dimDesc,

                                        uint64_t nSize)

            {

                SourceShortDimDesc sourceDesc;

                sourceDesc.nSize = nSize;

                sourceDesc.bIsRegularlySpaced = dimDesc.aaValues.empty();

                if (sourceDesc.bIsRegularlySpaced)

                    sourceDesc.dfStart = dimDesc.dfStart;

                else

                    sourceDesc.dfStart = dimDesc.aaValues[0][0];

                aoSourceShortDimDesc.push_back(std::move(sourceDesc));

            };

            const auto anBlockSize = poArray->GetBlockSize();

            CPLAssert(anBlockSize.size() == poArray->GetDimensionCount());

            if (!arrayParameters.poFirstSourceArray)

            {

                arrayParameters.poFirstSourceArray = poArray;

                CPLAssert(arrayParameters.mosaicDimensions.empty());

                size_t iDim = 0;

                for (auto &poDim : poArray->GetDimensions())

                {

                    auto descOpt = GetDimensionDesc(pszDatasetName, poDim);

                    if (!descOpt.has_value())

                        return false;

                    auto &desc = descOpt.value();

                    AddToSourceShortDimDesc(desc, poDim->GetSize());

                    desc.nBlockSize = anBlockSize[iDim];

                    arrayParameters.mosaicDimensions.push_back(std::move(desc));

                    ++iDim;

                }

            }

            else

            {

                if (poArray->GetDimensionCount() !=

                    arrayParameters.mosaicDimensions.size())

                {

                    ReportError(CE_Failure, CPLE_AppDefined,

                                "Array %s of dataset %s does not have the same "

                                "number of dimensions as in other datasets",

                                poArray->GetName().c_str(), pszDatasetName);

                    return false;

                }

                if (poArray->GetDataType() !=

                    arrayParameters.poFirstSourceArray->GetDataType())

                {

                    ReportError(CE_Failure, CPLE_AppDefined,

                                "Array %s of dataset %s does not have the same "

                                "data type as in other datasets",

                                poArray->GetName().c_str(), pszDatasetName);

                    return false;

                }

                const void *poFirstRawNoData =

                    arrayParameters.poFirstSourceArray->GetRawNoDataValue();

                const void *poRawNoData = poArray->GetRawNoDataValue();

                if (!((!poFirstRawNoData && !poRawNoData) ||

                      (poFirstRawNoData && poRawNoData &&

                       memcmp(poFirstRawNoData, poRawNoData,

                              poArray->GetDataType().GetSize()) == 0)))

                {

                    ReportError(CE_Failure, CPLE_AppDefined,

                                "Array %s of dataset %s does not have the same "

                                "nodata value as in other datasets",

                                poArray->GetName().c_str(), pszDatasetName);

                    return false;

                }

                const std::vector<std::shared_ptr<GDALDimension>> apoDims =

                    poArray->GetDimensions();

                for (size_t iDim = 0;

                     iDim < arrayParameters.mosaicDimensions.size(); ++iDim)

                {

                    DimensionDesc &desc =

                        arrayParameters.mosaicDimensions[iDim];

                    auto &poDim = apoDims[iDim];

                    if (poDim->GetName() != desc.osName)

                    {

                        ReportError(

                            CE_Failure, CPLE_AppDefined,

                            "Dimension %d of array %s of dataset %s does "

                            "not have the same name as in other datasets",

                            static_cast<int>(iDim), poArray->GetName().c_str(),

                            pszDatasetName);

                        return false;

                    }

                    if (desc.nBlockSize != anBlockSize[iDim])

                        desc.nBlockSize = 0;

                    auto descThisDatasetOpt =

                        GetDimensionDesc(pszDatasetName, poDim);

                    if (!descThisDatasetOpt.has_value())

                        return false;

                    auto &descThisDataset = descThisDatasetOpt.value();

                    AddToSourceShortDimDesc(descThisDataset, poDim->GetSize());

                    if (descThisDataset.bHasIndexingVar &&

                        !desc.bHasIndexingVar)

                    {

                        ReportError(CE_Failure, CPLE_AppDefined,

                                    "Dimension %s of array %s of dataset %s "

                                    "has an indexing variable, contrary "

                                    "to other datasets",

                                    poDim->GetName().c_str(),

                                    poArray->GetName().c_str(), pszDatasetName);

                        return false;

                    }

                    else if (!descThisDataset.bHasIndexingVar &&

                             desc.bHasIndexingVar)

                    {

                        ReportError(

                            CE_Failure, CPLE_AppDefined,

                            "Dimension %s of array %s of dataset %s "

                            "does not have an indexing variable, contrary "

                            "to other datasets",

                            poDim->GetName().c_str(),

                            poArray->GetName().c_str(), pszDatasetName);

                        return false;

                    }

                    else if (!desc.bHasIndexingVar &&

                             descThisDataset.nSize != desc.nSize)

                    {

                        ReportError(

                            CE_Failure, CPLE_AppDefined,

                            "Dimension %s of array %s of dataset %s "

                            "does not have the same size as in other datasets",

                            poDim->GetName().c_str(),

                            poArray->GetName().c_str(), pszDatasetName);

                        return false;

                    }

                    else if (desc.bHasIndexingVar && desc.aaValues.empty())

                    {

                        if (!descThisDataset.aaValues.empty())

                        {

                            ReportError(

                                CE_Failure, CPLE_AppDefined,

                                "Dimension %s of array %s of dataset %s "

                                "has irregularly-spaced values, contrary "

                                "to other datasets",

                                poDim->GetName().c_str(),

                                poArray->GetName().c_str(), pszDatasetName);

                            return false;

                        }

                        if (std::fabs(descThisDataset.dfIncrement -

                                      desc.dfIncrement) >

                            1e-10 * std::fabs(desc.dfIncrement))

                        {

                            ReportError(

                                CE_Failure, CPLE_AppDefined,

                                "Dimension %s of array %s of dataset %s is "

                                "indexed by a variable with spacing %g, "

                                "whereas it is %g in other datasets",

                                poDim->GetName().c_str(),

                                poArray->GetName().c_str(), pszDatasetName,

                                descThisDataset.dfIncrement, desc.dfIncrement);

                            return false;

                        }

                        const double dfPos =

                            (descThisDataset.dfStart - desc.dfStart) /

                            desc.dfIncrement;

                        if (std::fabs(std::round(dfPos) - dfPos) > 1e-3)

                        {

                            ReportError(CE_Failure, CPLE_AppDefined,

                                        "Dimension %s of array %s of dataset "

                                        "%s is indexed "

                                        "by a variable whose start value is "

                                        "not aligned "

                                        "with the one of other datasets",

                                        poDim->GetName().c_str(),

                                        poArray->GetName().c_str(),

                                        pszDatasetName);

                            return false;

                        }

                        const double dfEnd = std::max(

                            desc.dfStart + static_cast<double>(desc.nSize) *

                                               desc.dfIncrement,

                            descThisDataset.dfStart +

                                static_cast<double>(descThisDataset.nSize) *

                                    descThisDataset.dfIncrement);

                        desc.dfStart =

                            std::min(desc.dfStart, descThisDataset.dfStart);

                        const double dfSize =

                            (dfEnd - desc.dfStart) / desc.dfIncrement;

                        constexpr double MAX_INTEGER_REPRESENTABLE =

                            static_cast<double>(1ULL << 53);

                        if (dfSize > MAX_INTEGER_REPRESENTABLE)

                        {

                            ReportError(

                                CE_Failure, CPLE_AppDefined,

                                "Dimension %s of array %s of dataset %s "

                                "would be too large if merged.",

                                poDim->GetName().c_str(),

                                poArray->GetName().c_str(), pszDatasetName);

                            return false;

                        }

                        desc.nSize = static_cast<uint64_t>(dfSize + 0.5);

                    }

                    else if (desc.bHasIndexingVar)

                    {

                        if (descThisDataset.aaValues.empty())

                        {

                            ReportError(

                                CE_Failure, CPLE_AppDefined,

                                "Dimension %s of array %s of dataset %s "

                                "has regularly spaced labels, contrary to "

                                "other datasets",

                                poDim->GetName().c_str(),

                                poArray->GetName().c_str(), pszDatasetName);

                            return false;

                        }

                        if (descThisDataset.nProgressionSign !=

                            desc.nProgressionSign)

                        {

                            ReportError(

                                CE_Failure, CPLE_AppDefined,

                                "Dataset %s: values in indexing variable %s of "

                                "dimension %s must be either increasing or "

                                "decreasing in all input datasets.",

                                pszDatasetName,

                                poDim->GetIndexingVariable()->GetName().c_str(),

                                desc.osName.c_str());

                            return false;

                        }

                        CPLAssert(descThisDataset.aaValues.size() == 1);

                        if (descThisDataset.aaValues[0][0] <

                            desc.aaValues[0][0])

                        {

                            if (descThisDataset.aaValues[0].back() >=

                                desc.aaValues[0][0])

                            {

                                ReportError(

                                    CE_Failure, CPLE_AppDefined,

                                    "Dataset %s: values in indexing variable "

                                    "%s of "

                                    "dimension %s are not the same as in other "

                                    "datasets",

                                    pszDatasetName,

                                    poDim->GetIndexingVariable()

                                        ->GetName()

                                        .c_str(),

                                    desc.osName.c_str());

                                return false;

                            }

                            desc.aaValues.insert(

                                desc.aaValues.begin(),

                                std::move(descThisDataset.aaValues[0]));

                        }

                        else

                        {

                            for (size_t i = 0; i < desc.aaValues.size(); ++i)

                            {

                                if (descThisDataset.aaValues[0][0] ==

                                    desc.aaValues[i][0])

                                {

                                    if (descThisDataset.aaValues[0] !=

                                        desc.aaValues[i])

                                    {

                                        ReportError(

                                            CE_Failure, CPLE_AppDefined,

                                            "Dataset %s: values in indexing "

                                            "variable %s of dimension %s are "

                                            "not "

                                            "the same as in other datasets",

                                            pszDatasetName,

                                            poDim->GetIndexingVariable()

                                                ->GetName()

                                                .c_str(),

                                            desc.osName.c_str());

                                        return false;

                                    }

                                }

                                else if (descThisDataset.aaValues[0][0] >

                                             desc.aaValues[i][0] &&

                                         (i + 1 == desc.aaValues.size() ||

                                          descThisDataset.aaValues[0][0] <

                                              desc.aaValues[i + 1][0]))

                                {

                                    if (descThisDataset.aaValues[0][0] <=

                                        desc.aaValues[i].back())

                                    {

                                        ReportError(

                                            CE_Failure, CPLE_AppDefined,

                                            "Dataset %s: values in indexing "

                                            "variable %s of dimension %s are "

                                            "overlapping with the ones of "

                                            "other "

                                            "datasets",

                                            pszDatasetName,

                                            poDim->GetIndexingVariable()

                                                ->GetName()

                                                .c_str(),

                                            desc.osName.c_str());

                                        return false;

                                    }

                                    if (i + 1 < desc.aaValues.size() &&

                                        descThisDataset.aaValues[0].back() >=

                                            desc.aaValues[i + 1][0])

                                    {

                                        ReportError(

                                            CE_Failure, CPLE_AppDefined,

                                            "Dataset %s: values in indexing "

                                            "variable %s of dimension %s are "

                                            "overlapping with the ones of "

                                            "other "

                                            "datasets",

                                            pszDatasetName,

                                            poDim->GetIndexingVariable()

                                                ->GetName()

                                                .c_str(),

                                            desc.osName.c_str());

                                        return false;

                                    }

                                    desc.aaValues.insert(

                                        desc.aaValues.begin() + i + 1,

                                        std::move(descThisDataset.aaValues[0]));

                                    break;

                                }

                            }

                        }

                    }

                }

            }

            arrayParameters.aaoSourceShortDimDesc.push_back(

                std::move(aoSourceShortDimDesc));

        }

    }

    return true;

}

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

/*           GDALMdimMosaicAlgorithm::GetInputDatasetNames()            */

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

bool GDALMdimMosaicAlgorithm::GetInputDatasetNames(

    GDALProgressFunc pfnProgress, void *pProgressData,

    CPLStringList &aosInputDatasetNames) const

{

    for (auto &ds : m_inputDatasets)

    {

        if (ds.GetName()[0] == '@')

        {

            auto f = VSIVirtualHandleUniquePtr(

                VSIFOpenL(ds.GetName().c_str() + 1, "r"));

            if (!f)

            {

                ReportError(CE_Failure, CPLE_FileIO, "Cannot open %s",

                            ds.GetName().c_str() + 1);

                return false;

            }

            while (const char *filename = CPLReadLineL(f.get()))

            {

                aosInputDatasetNames.push_back(filename);

            }

        }

        else if (ds.GetName().find_first_of("*?[") != std::string::npos)

        {

            CPLStringList aosMatches(VSIGlob(ds.GetName().c_str(), nullptr,

                                             pfnProgress, pProgressData));

            for (const char *pszStr : aosMatches)

            {

                aosInputDatasetNames.push_back(pszStr);

            }

        }

        else

        {

            std::string osDatasetName = ds.GetName();

            if (!GetReferencePathForRelativePaths().empty())

            {

                osDatasetName = GDALDataset::BuildFilename(

                    osDatasetName.c_str(),

                    GetReferencePathForRelativePaths().c_str(), true);

            }

            aosInputDatasetNames.push_back(osDatasetName.c_str());

        }

    }

    return true;

}

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

/*                  GDALMdimMosaicAlgorithm::RunImpl()                  */

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

bool GDALMdimMosaicAlgorithm::RunImpl(GDALProgressFunc pfnProgress,

                                      void *pProgressData)

{

    CPLAssert(!m_outputDataset.GetDatasetRef());

    if (m_outputFormat.empty())

    {

        const auto aosFormats =

            CPLStringList(GDALGetOutputDriversForDatasetName(

                m_outputDataset.GetName().c_str(), GDAL_OF_MULTIDIM_RASTER,

                /* bSingleMatch = */ true,

                /* bWarn = */ true));

        if (aosFormats.size() != 1)

        {

            ReportError(CE_Failure, CPLE_AppDefined,

                        "Cannot guess driver for %s",

                        m_outputDataset.GetName().c_str());

            return false;

        }

        m_outputFormat = aosFormats[0];

    }

    auto poOutDrv =

        GetGDALDriverManager()->GetDriverByName(m_outputFormat.c_str());

    if (!poOutDrv)

    {

        ReportError(CE_Failure, CPLE_AppDefined, "Driver %s does not exist",

                    m_outputFormat.c_str());

        return false;

    }

    const bool bIsVRT = EQUAL(m_outputFormat.c_str(), "VRT");

    CPLStringList aosInputDatasetNames;

    const double dfIntermediatePercentage = bIsVRT ? 1.0 : 0.1;

    std::unique_ptr<void, decltype(&GDALDestroyScaledProgress)> pScaledData(

        GDALCreateScaledProgress(0.0, dfIntermediatePercentage, pfnProgress,

                                 pProgressData),

        GDALDestroyScaledProgress);

    if (!GetInputDatasetNames(GDALScaledProgress, pScaledData.get(),

                              aosInputDatasetNames))

        return false;

    for (std::string &name : m_array)

    {

        if (!name.empty() && name[0] != '/')

            name = "/" + name;

    }

    std::vector<ArrayParameters> aoArrayParameters;

    if (!BuildArrayParameters(aosInputDatasetNames, aoArrayParameters))

    {

        return false;

    }

    auto poVRTDS = VRTDataset::CreateVRTMultiDimensional("", nullptr, nullptr);

    CPLAssert(poVRTDS);

    auto poDstGroup = poVRTDS->GetRootVRTGroup();

    CPLAssert(poDstGroup);

    std::map<std::string, std::shared_ptr<GDALDimension>>

        oMapAlreadyCreatedDims;

    // Iterate over arrays

    for (auto &arrayParameters : aoArrayParameters)

    {

        auto &poFirstSourceArray = arrayParameters.poFirstSourceArray;

        CPLAssert(poFirstSourceArray);

        auto &aaoSourceShortDimDesc = arrayParameters.aaoSourceShortDimDesc;

        CPLAssert(aaoSourceShortDimDesc.size() ==

                  static_cast<size_t>(aosInputDatasetNames.size()));

        // Create mosaic array dimensions

        std::vector<std::shared_ptr<GDALDimension>> apoDstDims;

        for (auto &desc : arrayParameters.mosaicDimensions)

        {

            auto oIterCreatedDims = oMapAlreadyCreatedDims.find(desc.osName);

            if (oIterCreatedDims != oMapAlreadyCreatedDims.end())

            {

                apoDstDims.push_back(oIterCreatedDims->second);

            }

            else

            {

                uint64_t nDimSize64 = desc.nSize;

                if (!desc.aaValues.empty())

                {

                    nDimSize64 = 0;

                    for (const auto &aValues : desc.aaValues)

                        nDimSize64 += aValues.size();

                }

                auto dstDim = poDstGroup->CreateDimension(

                    desc.osName, desc.osType, desc.osDirection, nDimSize64);

                if (!dstDim)

                    return false;

                if (desc.bHasIndexingVar)

                {

                    auto var = poDstGroup->CreateVRTMDArray(

                        desc.osName, {dstDim},

                        GDALExtendedDataType::Create(GDT_Float64));

                    if (!var)

                        return false;

                    for (const auto &attr : desc.attributes)

                    {

                        CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);

                        auto dstAttr = var->CreateAttribute(

                            attr->GetName(), attr->GetDimensionsSize(),

                            attr->GetDataType());

                        if (dstAttr)

                        {

                            auto raw(attr->ReadAsRaw());

                            CPL_IGNORE_RET_VAL(

                                dstAttr->Write(raw.data(), raw.size()));

                        }

                    }

                    if (desc.aaValues.empty())

                    {

                        auto poSource =

                            std::make_unique<VRTMDArraySourceRegularlySpaced>(

                                desc.dfStart, desc.dfIncrement);

                        var->AddSource(std::move(poSource));

                    }

                    else

                    {

                        const size_t nDimSize = static_cast<size_t>(nDimSize64);

                        std::vector<GUInt64> anOffset = {0};

                        std::vector<size_t> anCount = {nDimSize};

                        std::vector<double> adfValues;

                        adfValues.reserve(nDimSize);

                        if (desc.nProgressionSign >= 0)

                        {

                            for (const auto &aValues : desc.aaValues)

                                adfValues.insert(adfValues.end(),

                                                 aValues.begin(),

                                                 aValues.end());

                        }

                        else

                        {

                            for (auto oIter = desc.aaValues.rbegin();

                                 oIter != desc.aaValues.rend(); ++oIter)

                            {

                                adfValues.insert(adfValues.end(),

                                                 oIter->rbegin(),

                                                 oIter->rend());

                            }

                        }

                        std::vector<GByte> abyValues(nDimSize * sizeof(double));

                        memcpy(abyValues.data(), adfValues.data(),

                               nDimSize * sizeof(double));

                        auto poSource =

                            std::make_unique<VRTMDArraySourceInlinedValues>(

                                var.get(),

                                /* bIsConstantValue = */ false,

                                std::move(anOffset), std::move(anCount),

                                std::move(abyValues));

                        var->AddSource(std::move(poSource));

                    }

                    dstDim->SetIndexingVariable(std::move(var));

                }

                oMapAlreadyCreatedDims[dstDim->GetName()] = dstDim;

                apoDstDims.push_back(std::move(dstDim));

            }

        }

        // Create mosaic array

        CPLStringList aosArrayCO;

        std::string osBlockSize;

        for (size_t i = 0; i < apoDstDims.size(); ++i)

        {

            if (i > 0)

                osBlockSize += ',';

            osBlockSize +=

                std::to_string(arrayParameters.mosaicDimensions[i].nBlockSize);

        }

        if (!osBlockSize.empty())

            aosArrayCO.SetNameValue("BLOCKSIZE", osBlockSize.c_str());

        auto poDstArray = poDstGroup->CreateVRTMDArray(

            CPLGetFilename(poFirstSourceArray->GetName().c_str()), apoDstDims,

            poFirstSourceArray->GetDataType(), aosArrayCO);

        if (!poDstArray)

            return false;

        GUInt64 nCurCost = 0;

        poDstArray->CopyFromAllExceptValues(poFirstSourceArray.get(), false,

                                            nCurCost, 0, nullptr, nullptr);

        // Add sources to mosaic array

        for (int iDS = 0; iDS < aosInputDatasetNames.size(); ++iDS)

        {

            const auto &aoSourceShortDimDesc = aaoSourceShortDimDesc[iDS];

            const auto dimCount = arrayParameters.mosaicDimensions.size();

            std::vector<GUInt64> anSrcOffset(dimCount);

            std::vector<GUInt64> anCount(dimCount);

            std::vector<GUInt64> anDstOffset;

            CPLAssert(aoSourceShortDimDesc.size() == dimCount);

            for (size_t iDim = 0; iDim < dimCount; ++iDim)

            {

                const DimensionDesc &desc =

                    arrayParameters.mosaicDimensions[iDim];

                const SourceShortDimDesc &sourceDesc =

                    aoSourceShortDimDesc[iDim];

                if (sourceDesc.bIsRegularlySpaced)

                {