Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  "reproject" step of "mdim pipeline"
       5             :  * Author:   Even Rouault <even dot rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2026, Even Rouault <even dot rouault at spatialys.com>
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "gdalalg_mdim_reproject.h"
      14             : #include "gdalalg_raster_reproject.h"
      15             : 
      16             : #include "gdal_priv.h"
      17             : 
      18             : #include <map>
      19             : #include <set>
      20             : #include <utility>
      21             : 
      22             : //! @cond Doxygen_Suppress
      23             : 
      24             : #ifndef _
      25             : #define _(x) (x)
      26             : #endif
      27             : 
      28             : /************************************************************************/
      29             : /*       GDALMdimReprojectAlgorithm::GDALMdimReprojectAlgorithm()       */
      30             : /************************************************************************/
      31             : 
      32          49 : GDALMdimReprojectAlgorithm::GDALMdimReprojectAlgorithm(bool standaloneStep)
      33             :     : GDALMdimPipelineStepAlgorithm(
      34             :           NAME, DESCRIPTION, HELP_URL,
      35          49 :           ConstructorOptions().SetStandaloneStep(standaloneStep))
      36             : {
      37          98 :     AddArg(GDAL_ARG_NAME_OUTPUT_CRS, 'd', _("Output CRS"), &m_dstCrs)
      38          98 :         .SetIsCRSArg()
      39          98 :         .AddHiddenAlias("t_srs")
      40          49 :         .AddHiddenAlias("dst-crs");
      41             : 
      42          49 :     GDALRasterReprojectUtils::AddResamplingArg(this, m_resampling);
      43          49 : }
      44             : 
      45             : namespace
      46             : {
      47             : 
      48             : /************************************************************************/
      49             : /*                       GDALMdimReprojectParams                        */
      50             : /************************************************************************/
      51             : 
      52             : struct GDALMdimReprojectParams
      53             : {
      54             :     GDALRIOResampleAlg eResampleAlg = GRIORA_NearestNeighbour;
      55             :     OGRSpatialReferenceRefCountedPtr poTargetSRS{};
      56             : };
      57             : 
      58             : /************************************************************************/
      59             : /*                        GDALMdimReprojectGroup                        */
      60             : /************************************************************************/
      61             : 
      62             : /** Wrapper around a source group object, that is essentially passthrough,
      63             :  * except with 2D+ arrays.
      64             :  */
      65             : class GDALMdimReprojectGroup final : public GDALGroup
      66             : {
      67             :     std::shared_ptr<GDALGroup> m_poSrcGroup{};
      68             :     const GDALMdimReprojectParams m_sParams;
      69             :     std::vector<std::string> m_aosArrayNames{};
      70             :     std::map<std::string, std::shared_ptr<GDALMDArray>> m_oMapArrays{};
      71             :     std::vector<std::shared_ptr<GDALDimension>> m_apoDims{};
      72             : 
      73             :   public:
      74             :     GDALMdimReprojectGroup(const std::string &osParentName,
      75             :                            const std::shared_ptr<GDALGroup> &poSrcGroup,
      76             :                            const GDALMdimReprojectParams &sParams);
      77             : 
      78             :     std::shared_ptr<GDALAttribute>
      79           0 :     GetAttribute(const std::string &osName) const override
      80             :     {
      81           0 :         return m_poSrcGroup->GetAttribute(osName);
      82             :     }
      83             : 
      84             :     std::vector<std::shared_ptr<GDALAttribute>>
      85           3 :     GetAttributes(CSLConstList papszOptions = nullptr) const override
      86             :     {
      87           3 :         return m_poSrcGroup->GetAttributes(papszOptions);
      88             :     }
      89             : 
      90             :     std::vector<std::string>
      91           3 :     GetGroupNames(CSLConstList papszOptions = nullptr) const override
      92             :     {
      93           3 :         return m_poSrcGroup->GetGroupNames(papszOptions);
      94             :     }
      95             : 
      96             :     std::shared_ptr<GDALGroup>
      97           0 :     OpenGroup(const std::string &osName,
      98             :               CSLConstList papszOptions = nullptr) const override
      99             :     {
     100           0 :         auto poSrcChildGroup = m_poSrcGroup->OpenGroup(osName, papszOptions);
     101           0 :         if (!poSrcChildGroup)
     102           0 :             return nullptr;
     103           0 :         return std::make_shared<GDALMdimReprojectGroup>(
     104           0 :             GetFullName(), std::move(poSrcChildGroup), m_sParams);
     105             :     }
     106             : 
     107           4 :     std::vector<std::string> GetMDArrayNames(CSLConstList) const override
     108             :     {
     109           4 :         return m_aosArrayNames;
     110             :     }
     111             : 
     112          16 :     std::shared_ptr<GDALMDArray> OpenMDArray(const std::string &osName,
     113             :                                              CSLConstList) const override
     114             :     {
     115          16 :         const auto oIter = m_oMapArrays.find(osName);
     116          16 :         if (oIter != m_oMapArrays.end())
     117          16 :             return oIter->second;
     118           0 :         return nullptr;
     119             :     }
     120             : 
     121             :     std::vector<std::shared_ptr<GDALDimension>>
     122             :         GetDimensions(CSLConstList) const override;
     123             : };
     124             : 
     125             : /************************************************************************/
     126             : /*           GDALMdimReprojectGroup::GDALMdimReprojectGroup()           */
     127             : /************************************************************************/
     128             : 
     129           5 : GDALMdimReprojectGroup::GDALMdimReprojectGroup(
     130             :     const std::string &osParentName,
     131             :     const std::shared_ptr<GDALGroup> &poSrcGroup,
     132           5 :     const GDALMdimReprojectParams &sParams)
     133             :     : GDALGroup(osParentName, poSrcGroup->GetName()), m_poSrcGroup(poSrcGroup),
     134           5 :       m_sParams(sParams)
     135             : {
     136             :     // We collect source arrays and dimensions, and remove dimensions that
     137             :     // are only referenced by the spatial dimensions of reprojected arrays.
     138             : 
     139             :     // First bool in pair: referenced by 1d variables (other than its own indexing variable)
     140             :     // Second bool in pair: referenced by >= 2d variables
     141          10 :     std::map<std::string, std::pair<bool, bool>> oMapArrayDims;
     142          10 :     std::map<std::string, std::shared_ptr<GDALDimension>> oMapNewDims;
     143          10 :     std::vector<std::shared_ptr<GDALMDArray>> apoIndexingVariables;
     144             : 
     145          21 :     for (const std::string &osName : m_poSrcGroup->GetMDArrayNames())
     146             :     {
     147          16 :         auto poSrcArray = m_poSrcGroup->OpenMDArray(osName);
     148          16 :         if (!poSrcArray)
     149           0 :             continue;
     150          16 :         if (poSrcArray->GetDimensionCount() == 0)
     151             :         {
     152           0 :             m_aosArrayNames.push_back(osName);
     153           0 :             m_oMapArrays[osName] = std::move(poSrcArray);
     154             :         }
     155          16 :         else if (poSrcArray->GetDimensionCount() == 1)
     156             :         {
     157          11 :             const auto &poDim = poSrcArray->GetDimensions()[0];
     158          11 :             if (poDim->GetName() == osName)
     159             :             {
     160          11 :                 apoIndexingVariables.push_back(std::move(poSrcArray));
     161             :             }
     162             :             else
     163             :             {
     164           0 :                 oMapArrayDims[poDim->GetName()].first = true;
     165           0 :                 m_aosArrayNames.push_back(osName);
     166           0 :                 m_oMapArrays[osName] = std::move(poSrcArray);
     167             :             }
     168             :         }
     169             :         else
     170             :         {
     171             :             std::vector<std::shared_ptr<GDALDimension>> apoNewDims(
     172          10 :                 poSrcArray->GetDimensionCount());
     173          10 :             CPLStringList aosOptions;
     174           5 :             aosOptions.SetNameValue("PARENT_PATH", GetFullName().c_str());
     175           5 :             aosOptions.SetNameValue("NAME", osName.c_str());
     176           5 :             auto poDstArray = poSrcArray->GetResampled(
     177           5 :                 apoNewDims, m_sParams.eResampleAlg, m_sParams.poTargetSRS.get(),
     178          15 :                 aosOptions.List());
     179           5 :             if (poDstArray)
     180             :             {
     181           5 :                 m_aosArrayNames.push_back(osName);
     182           5 :                 CPLAssert(poDstArray->GetDimensionCount() ==
     183             :                           poSrcArray->GetDimensionCount());
     184          16 :                 for (size_t i = 0; i < poDstArray->GetDimensionCount(); ++i)
     185             :                 {
     186          11 :                     const auto &poSrcDim = poSrcArray->GetDimensions()[i];
     187          11 :                     const auto &poDstDim = poDstArray->GetDimensions()[i];
     188          11 :                     if (poSrcDim.get() != poDstDim.get())
     189             :                     {
     190          10 :                         oMapArrayDims[poSrcDim->GetName()].second = true;
     191          10 :                         oMapNewDims[poDstDim->GetName()] = poDstDim;
     192             : 
     193          20 :                         auto poVar = poDstDim->GetIndexingVariable();
     194          10 :                         if (poVar)
     195             :                         {
     196          10 :                             m_aosArrayNames.push_back(poVar->GetName());
     197          10 :                             m_oMapArrays[poVar->GetName()] = poVar;
     198             :                         }
     199             :                     }
     200             :                 }
     201           5 :                 m_oMapArrays[osName] = std::move(poDstArray);
     202             :             }
     203             :         }
     204             :     }
     205             : 
     206          16 :     for (auto &poArray : apoIndexingVariables)
     207             :     {
     208          11 :         auto oIter = oMapArrayDims.find(poArray->GetName());
     209          11 :         if (oIter == oMapArrayDims.end() || !oIter->second.second)
     210             :         {
     211           1 :             m_aosArrayNames.push_back(poArray->GetName());
     212           1 :             m_oMapArrays[poArray->GetName()] = poArray;
     213             :         }
     214             :     }
     215             : 
     216          16 :     for (const auto &poDim : m_poSrcGroup->GetDimensions())
     217             :     {
     218          11 :         auto oIter = oMapArrayDims.find(poDim->GetName());
     219          11 :         if (oIter == oMapArrayDims.end() || !oIter->second.second)
     220             :         {
     221           1 :             m_apoDims.push_back(poDim);
     222             :         }
     223             :     }
     224          15 :     for (const auto &[_, poDim] : oMapNewDims)
     225             :     {
     226          10 :         m_apoDims.push_back(poDim);
     227             :     }
     228           5 : }
     229             : 
     230             : /************************************************************************/
     231             : /*               GDALMdimReprojectGroup::GetDimensions()                */
     232             : /************************************************************************/
     233             : 
     234             : std::vector<std::shared_ptr<GDALDimension>>
     235           2 : GDALMdimReprojectGroup::GetDimensions(CSLConstList) const
     236             : {
     237           2 :     return m_apoDims;
     238             : }
     239             : 
     240             : /************************************************************************/
     241             : /*                       GDALMdimReprojectDataset                       */
     242             : /************************************************************************/
     243             : 
     244             : class GDALMdimReprojectDataset final : public GDALDataset
     245             : {
     246             :     std::shared_ptr<GDALGroup> m_poRootGroup{};
     247             : 
     248             :   public:
     249           5 :     GDALMdimReprojectDataset(GDALDataset *poSrcDS,
     250             :                              const GDALMdimReprojectParams &params)
     251           5 :     {
     252          10 :         auto poSrcRootGroup = poSrcDS->GetRootGroup();
     253           5 :         if (poSrcRootGroup)
     254             :         {
     255          10 :             m_poRootGroup = std::make_shared<GDALMdimReprojectGroup>(
     256          15 :                 std::string(), std::move(poSrcRootGroup), params);
     257             :         }
     258           5 :     }
     259             : 
     260             :     std::shared_ptr<GDALGroup> GetRootGroup() const override;
     261             : };
     262             : 
     263           6 : std::shared_ptr<GDALGroup> GDALMdimReprojectDataset::GetRootGroup() const
     264             : {
     265           6 :     return m_poRootGroup;
     266             : }
     267             : 
     268             : }  // namespace
     269             : 
     270             : /************************************************************************/
     271             : /*                GDALMdimReprojectAlgorithm::RunStep()                 */
     272             : /************************************************************************/
     273             : 
     274           5 : bool GDALMdimReprojectAlgorithm::RunStep(GDALPipelineStepRunContext &)
     275             : {
     276           5 :     auto poSrcDS = m_inputDataset[0].GetDatasetRef();
     277           5 :     CPLAssert(poSrcDS);
     278           5 :     CPLAssert(m_outputDataset.GetName().empty());
     279           5 :     CPLAssert(!m_outputDataset.GetDatasetRef());
     280             : 
     281           5 :     GDALMdimReprojectParams sParams;
     282           5 :     sParams.eResampleAlg = GDALRasterIOGetResampleAlg(m_resampling.c_str());
     283           5 :     if (!m_dstCrs.empty())
     284             :     {
     285           4 :         sParams.poTargetSRS = OGRSpatialReferenceRefCountedPtr::makeInstance();
     286           4 :         sParams.poTargetSRS->SetFromUserInput(m_dstCrs.c_str());
     287             :     }
     288           5 :     m_outputDataset.Set(
     289          10 :         std::make_unique<GDALMdimReprojectDataset>(poSrcDS, sParams));
     290             : 
     291          10 :     return true;
     292             : }
     293             : 
     294             : GDALMdimReprojectAlgorithmStandalone::~GDALMdimReprojectAlgorithmStandalone() =
     295             :     default;
     296             : 
     297             : //! @endcond