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          46 : GDALMdimReprojectAlgorithm::GDALMdimReprojectAlgorithm(bool standaloneStep)
      33             :     : GDALMdimPipelineStepAlgorithm(
      34             :           NAME, DESCRIPTION, HELP_URL,
      35          46 :           ConstructorOptions().SetStandaloneStep(standaloneStep))
      36             : {
      37          92 :     AddArg(GDAL_ARG_NAME_OUTPUT_CRS, 'd', _("Output CRS"), &m_dstCrs)
      38          92 :         .SetIsCRSArg()
      39          92 :         .AddHiddenAlias("t_srs")
      40          46 :         .AddHiddenAlias("dst-crs");
      41             : 
      42          46 :     GDALRasterReprojectUtils::AddResamplingArg(this, m_resampling);
      43          46 : }
      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           1 :     GetAttributes(CSLConstList papszOptions = nullptr) const override
      86             :     {
      87           1 :         return m_poSrcGroup->GetAttributes(papszOptions);
      88             :     }
      89             : 
      90             :     std::vector<std::string>
      91           1 :     GetGroupNames(CSLConstList papszOptions = nullptr) const override
      92             :     {
      93           1 :         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           1 :     std::vector<std::string> GetMDArrayNames(CSLConstList) const override
     108             :     {
     109           1 :         return m_aosArrayNames;
     110             :     }
     111             : 
     112           5 :     std::shared_ptr<GDALMDArray> OpenMDArray(const std::string &osName,
     113             :                                              CSLConstList) const override
     114             :     {
     115           5 :         const auto oIter = m_oMapArrays.find(osName);
     116           5 :         if (oIter != m_oMapArrays.end())
     117           3 :             return oIter->second;
     118           2 :         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           3 : GDALMdimReprojectGroup::GDALMdimReprojectGroup(
     130             :     const std::string &osParentName,
     131             :     const std::shared_ptr<GDALGroup> &poSrcGroup,
     132           3 :     const GDALMdimReprojectParams &sParams)
     133             :     : GDALGroup(osParentName, poSrcGroup->GetName()), m_poSrcGroup(poSrcGroup),
     134           3 :       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           6 :     std::map<std::string, std::pair<bool, bool>> oMapArrayDims;
     142           6 :     std::map<std::string, std::shared_ptr<GDALDimension>> oMapNewDims;
     143           6 :     std::vector<std::shared_ptr<GDALMDArray>> apoIndexingVariables;
     144             : 
     145          12 :     for (const std::string &osName : m_poSrcGroup->GetMDArrayNames())
     146             :     {
     147           9 :         auto poSrcArray = m_poSrcGroup->OpenMDArray(osName);
     148           9 :         if (!poSrcArray)
     149           0 :             continue;
     150           9 :         if (poSrcArray->GetDimensionCount() == 0)
     151             :         {
     152           0 :             m_aosArrayNames.push_back(osName);
     153           0 :             m_oMapArrays[osName] = std::move(poSrcArray);
     154             :         }
     155           9 :         else if (poSrcArray->GetDimensionCount() == 1)
     156             :         {
     157           6 :             m_aosArrayNames.push_back(osName);
     158           6 :             const auto &poDim = poSrcArray->GetDimensions()[0];
     159           6 :             if (poDim->GetName() == osName)
     160             :             {
     161           6 :                 apoIndexingVariables.push_back(poSrcArray);
     162             :             }
     163             :             else
     164             :             {
     165           0 :                 oMapArrayDims[poDim->GetName()].first = true;
     166           0 :                 m_oMapArrays[osName] = std::move(poSrcArray);
     167             :             }
     168             :         }
     169             :         else
     170             :         {
     171             :             std::vector<std::shared_ptr<GDALDimension>> apoNewDims(
     172           6 :                 poSrcArray->GetDimensionCount());
     173           6 :             CPLStringList aosOptions;
     174           3 :             aosOptions.SetNameValue("PARENT_PATH", GetFullName().c_str());
     175           3 :             aosOptions.SetNameValue("NAME", osName.c_str());
     176           3 :             auto poDstArray = poSrcArray->GetResampled(
     177           3 :                 apoNewDims, m_sParams.eResampleAlg, m_sParams.poTargetSRS.get(),
     178           9 :                 aosOptions.List());
     179           3 :             if (poDstArray)
     180             :             {
     181           3 :                 m_aosArrayNames.push_back(osName);
     182           3 :                 CPLAssert(poDstArray->GetDimensionCount() ==
     183             :                           poSrcArray->GetDimensionCount());
     184           9 :                 for (size_t i = 0; i < poDstArray->GetDimensionCount(); ++i)
     185             :                 {
     186           6 :                     const auto &poSrcDim = poSrcArray->GetDimensions()[i];
     187           6 :                     const auto &poDstDim = poDstArray->GetDimensions()[i];
     188           6 :                     if (poSrcDim.get() != poDstDim.get())
     189             :                     {
     190           6 :                         oMapArrayDims[poSrcDim->GetName()].second = true;
     191           6 :                         oMapNewDims[poDstDim->GetName()] = poDstDim;
     192             : 
     193          12 :                         auto poVar = poDstDim->GetIndexingVariable();
     194           6 :                         if (poVar)
     195             :                         {
     196           6 :                             m_aosArrayNames.push_back(poVar->GetName());
     197           6 :                             m_oMapArrays[poVar->GetName()] = poVar;
     198             :                         }
     199             :                     }
     200             :                 }
     201           3 :                 m_oMapArrays[osName] = std::move(poDstArray);
     202             :             }
     203             :         }
     204             :     }
     205             : 
     206           9 :     for (auto &poArray : apoIndexingVariables)
     207             :     {
     208           6 :         auto oIter = oMapArrayDims.find(poArray->GetName());
     209           6 :         if (oIter == oMapArrayDims.end() || !oIter->second.second)
     210             :         {
     211           0 :             m_aosArrayNames.push_back(poArray->GetName());
     212           0 :             m_oMapArrays[poArray->GetName()] = poArray;
     213             :         }
     214             :     }
     215             : 
     216           9 :     for (const auto &poDim : m_poSrcGroup->GetDimensions())
     217             :     {
     218           6 :         auto oIter = oMapArrayDims.find(poDim->GetName());
     219           6 :         if (oIter == oMapArrayDims.end() || !oIter->second.second)
     220             :         {
     221           0 :             m_apoDims.push_back(poDim);
     222             :         }
     223             :     }
     224           9 :     for (const auto &[_, poDim] : oMapNewDims)
     225             :     {
     226           6 :         m_apoDims.push_back(poDim);
     227             :     }
     228           3 : }
     229             : 
     230             : /************************************************************************/
     231             : /*               GDALMdimReprojectGroup::GetDimensions()                */
     232             : /************************************************************************/
     233             : 
     234             : std::vector<std::shared_ptr<GDALDimension>>
     235           1 : GDALMdimReprojectGroup::GetDimensions(CSLConstList) const
     236             : {
     237           1 :     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           3 :     GDALMdimReprojectDataset(GDALDataset *poSrcDS,
     250             :                              const GDALMdimReprojectParams &params)
     251           3 :     {
     252           6 :         auto poSrcRootGroup = poSrcDS->GetRootGroup();
     253           3 :         if (poSrcRootGroup)
     254             :         {
     255           6 :             m_poRootGroup = std::make_shared<GDALMdimReprojectGroup>(
     256           9 :                 std::string(), std::move(poSrcRootGroup), params);
     257             :         }
     258           3 :     }
     259             : 
     260             :     std::shared_ptr<GDALGroup> GetRootGroup() const override;
     261             : };
     262             : 
     263           1 : std::shared_ptr<GDALGroup> GDALMdimReprojectDataset::GetRootGroup() const
     264             : {
     265           1 :     return m_poRootGroup;
     266             : }
     267             : 
     268             : }  // namespace
     269             : 
     270             : /************************************************************************/
     271             : /*                GDALMdimReprojectAlgorithm::RunStep()                 */
     272             : /************************************************************************/
     273             : 
     274           3 : bool GDALMdimReprojectAlgorithm::RunStep(GDALPipelineStepRunContext &)
     275             : {
     276           3 :     auto poSrcDS = m_inputDataset[0].GetDatasetRef();
     277           3 :     CPLAssert(poSrcDS);
     278           3 :     CPLAssert(m_outputDataset.GetName().empty());
     279           3 :     CPLAssert(!m_outputDataset.GetDatasetRef());
     280             : 
     281           3 :     GDALMdimReprojectParams sParams;
     282           3 :     sParams.eResampleAlg = GDALRasterIOGetResampleAlg(m_resampling.c_str());
     283           3 :     if (!m_dstCrs.empty())
     284             :     {
     285           2 :         sParams.poTargetSRS = OGRSpatialReferenceRefCountedPtr::makeInstance();
     286           2 :         sParams.poTargetSRS->SetFromUserInput(m_dstCrs.c_str());
     287             :     }
     288           3 :     m_outputDataset.Set(
     289           6 :         std::make_unique<GDALMdimReprojectDataset>(poSrcDS, sParams));
     290             : 
     291           6 :     return true;
     292             : }
     293             : 
     294             : GDALMdimReprojectAlgorithmStandalone::~GDALMdimReprojectAlgorithmStandalone() =
     295             :     default;
     296             : 
     297             : //! @endcond