Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  "pansharpen" step of "raster pipeline"
       5             :  * Author:   Even Rouault <even dot rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "gdalalg_raster_pansharpen.h"
      14             : 
      15             : #include "gdal_priv.h"
      16             : #include "cpl_minixml.h"
      17             : 
      18             : //! @cond Doxygen_Suppress
      19             : 
      20             : #ifndef _
      21             : #define _(x) (x)
      22             : #endif
      23             : 
      24             : /************************************************************************/
      25             : /*       GDALRasterPansharpenAlgorithm::GetConstructorOptions()         */
      26             : /************************************************************************/
      27             : 
      28             : /* static */ GDALRasterPansharpenAlgorithm::ConstructorOptions
      29          33 : GDALRasterPansharpenAlgorithm::GetConstructorOptions(bool standaloneStep)
      30             : {
      31          33 :     ConstructorOptions opts;
      32          33 :     opts.SetStandaloneStep(standaloneStep);
      33          33 :     opts.SetAddDefaultArguments(false);
      34          33 :     opts.SetInputDatasetAlias("panchromatic");
      35          33 :     opts.SetInputDatasetHelpMsg(_("Input panchromatic raster dataset"));
      36          33 :     return opts;
      37             : }
      38             : 
      39             : /************************************************************************/
      40             : /*     GDALRasterPansharpenAlgorithm::GDALRasterPansharpenAlgorithm()   */
      41             : /************************************************************************/
      42             : 
      43          33 : GDALRasterPansharpenAlgorithm::GDALRasterPansharpenAlgorithm(
      44          33 :     bool standaloneStep)
      45             :     : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
      46          33 :                                       GetConstructorOptions(standaloneStep))
      47             : {
      48          33 :     const auto AddSpectralDatasetArg = [this]()
      49             :     {
      50             :         auto &arg = AddArg("spectral", 0, _("Input spectral band dataset"),
      51          66 :                            &m_spectralDatasets)
      52          33 :                         .SetPositional()
      53          33 :                         .SetRequired()
      54          33 :                         .SetMinCount(1)
      55             :                         // due to ",band=" comma syntax
      56          33 :                         .SetAutoOpenDataset(false)
      57             :                         // due to ",band=" comma syntax
      58          33 :                         .SetPackedValuesAllowed(false)
      59          33 :                         .SetMetaVar("SPECTRAL");
      60             : 
      61          33 :         SetAutoCompleteFunctionForFilename(arg, GDAL_OF_RASTER);
      62          33 :     };
      63             : 
      64          33 :     if (standaloneStep)
      65             :     {
      66          13 :         AddRasterInputArgs(false, false);
      67          13 :         AddSpectralDatasetArg();
      68          13 :         AddProgressArg();
      69          13 :         AddRasterOutputArgs(false);
      70             :     }
      71             :     else
      72             :     {
      73          20 :         AddRasterHiddenInputDatasetArg();
      74          20 :         AddSpectralDatasetArg();
      75             :     }
      76             : 
      77          66 :     AddArg("resampling", 'r', _("Resampling algorithm"), &m_resampling)
      78          33 :         .SetDefault(m_resampling)
      79             :         .SetChoices("nearest", "bilinear", "cubic", "cubicspline", "lanczos",
      80          33 :                     "average");
      81          33 :     AddArg("weights", 0, _("Weight for each input spectral band"), &m_weights);
      82          33 :     AddArg("nodata", 0, _("Override nodata value of input bands"), &m_nodata);
      83          66 :     AddArg("bit-depth", 0, _("Override bit depth of input bands"), &m_bitDepth)
      84          33 :         .SetMinValueIncluded(8);
      85             :     AddArg("spatial-extent-adjustment", 0,
      86             :            _("Select behavior when bands have not the same extent"),
      87          66 :            &m_spatialExtentAdjustment)
      88          33 :         .SetDefault(m_spatialExtentAdjustment)
      89          33 :         .SetChoices("union", "intersection", "none", "none-without-warning");
      90          33 :     AddNumThreadsArg(&m_numThreads, &m_numThreadsStr);
      91          33 : }
      92             : 
      93             : /************************************************************************/
      94             : /*            GDALRasterPansharpenAlgorithm::RunStep()                  */
      95             : /************************************************************************/
      96             : 
      97          11 : bool GDALRasterPansharpenAlgorithm::RunStep(GDALPipelineStepRunContext &)
      98             : {
      99          11 :     auto poPanDS = m_inputDataset[0].GetDatasetRef();
     100          11 :     CPLAssert(poPanDS);
     101          11 :     CPLAssert(m_outputDataset.GetName().empty());
     102          11 :     CPLAssert(!m_outputDataset.GetDatasetRef());
     103             : 
     104          11 :     if (poPanDS->GetRasterCount() != 1)
     105             :     {
     106           1 :         ReportError(CE_Failure, CPLE_AppDefined,
     107             :                     "Input panchromatic dataset must have a single band");
     108           1 :         return false;
     109             :     }
     110             : 
     111             :     // to keep in this scope to keep datasets of spectral bands open until
     112             :     // GDALCreatePansharpenedVRT() runs
     113             :     std::vector<std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>>
     114          20 :         apoDatasetsToReleaseRef;
     115          20 :     std::vector<GDALRasterBandH> ahSpectralBands;
     116             : 
     117          20 :     for (auto &spectralDataset : m_spectralDatasets)
     118             :     {
     119          12 :         if (auto poSpectralDS = spectralDataset.GetDatasetRef())
     120             :         {
     121          12 :             for (int i = 1; i <= poSpectralDS->GetRasterCount(); ++i)
     122             :             {
     123           9 :                 ahSpectralBands.push_back(
     124           9 :                     GDALRasterBand::ToHandle(poSpectralDS->GetRasterBand(i)));
     125             :             }
     126             :         }
     127             :         else
     128             :         {
     129           9 :             const auto &name = spectralDataset.GetName();
     130           9 :             std::string dsName(name);
     131           9 :             const auto pos = name.find(",band=");
     132           9 :             int iBand = 0;
     133           9 :             if (pos != std::string::npos)
     134             :             {
     135           4 :                 dsName.resize(pos);
     136           4 :                 iBand = atoi(name.c_str() + pos + strlen(",band="));
     137             :             }
     138             :             std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser> poDS(
     139             :                 GDALDataset::Open(dsName.c_str(),
     140           9 :                                   GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR));
     141           9 :             if (!poDS)
     142           1 :                 return false;
     143             : 
     144           8 :             if (iBand <= 0)
     145             :             {
     146          16 :                 for (int i = 1; i <= poDS->GetRasterCount(); ++i)
     147             :                 {
     148          12 :                     ahSpectralBands.push_back(
     149          12 :                         GDALRasterBand::ToHandle(poDS->GetRasterBand(i)));
     150             :                 }
     151             :             }
     152           4 :             else if (iBand > poDS->GetRasterCount())
     153             :             {
     154           1 :                 ReportError(CE_Failure, CPLE_IllegalArg, "Illegal band in '%s'",
     155             :                             name.c_str());
     156           1 :                 return false;
     157             :             }
     158             :             else
     159             :             {
     160           3 :                 ahSpectralBands.push_back(
     161           3 :                     GDALRasterBand::ToHandle(poDS->GetRasterBand(iBand)));
     162             :             }
     163             : 
     164           7 :             apoDatasetsToReleaseRef.push_back(std::move(poDS));
     165             :         }
     166             :     }
     167             : 
     168          16 :     CPLXMLTreeCloser root(CPLCreateXMLNode(nullptr, CXT_Element, "VRTDataset"));
     169          32 :     for (int i = 0; i < static_cast<int>(ahSpectralBands.size()); ++i)
     170             :     {
     171             :         auto psBandNode =
     172          24 :             CPLCreateXMLNode(root.get(), CXT_Element, "VRTRasterBand");
     173          24 :         CPLAddXMLAttributeAndValue(
     174             :             psBandNode, "dataType",
     175          24 :             GDALGetDataTypeName(GDALGetRasterDataType(ahSpectralBands[i])));
     176          24 :         CPLAddXMLAttributeAndValue(psBandNode, "band", CPLSPrintf("%d", i + 1));
     177          24 :         CPLAddXMLAttributeAndValue(psBandNode, "subClass",
     178             :                                    "VRTPansharpenedRasterBand");
     179             :     }
     180           8 :     CPLAddXMLAttributeAndValue(root.get(), "subClass",
     181             :                                "VRTPansharpenedDataset");
     182             :     auto psPansharpeningOptionsNode =
     183           8 :         CPLCreateXMLNode(root.get(), CXT_Element, "PansharpeningOptions");
     184           8 :     if (!m_weights.empty())
     185             :     {
     186           2 :         auto psAlgorithmOptionsNode = CPLCreateXMLNode(
     187             :             psPansharpeningOptionsNode, CXT_Element, "AlgorithmOptions");
     188           4 :         std::string osWeights;
     189           8 :         for (double w : m_weights)
     190             :         {
     191           6 :             if (!osWeights.empty())
     192           4 :                 osWeights += ',';
     193           6 :             osWeights += CPLSPrintf("%.17g", w);
     194             :         }
     195           2 :         CPLCreateXMLElementAndValue(psAlgorithmOptionsNode, "Weights",
     196             :                                     osWeights.c_str());
     197             :     }
     198           8 :     CPLCreateXMLElementAndValue(psPansharpeningOptionsNode, "Resampling",
     199             :                                 m_resampling.c_str());
     200           8 :     CPLCreateXMLElementAndValue(psPansharpeningOptionsNode, "NumThreads",
     201             :                                 m_numThreadsStr.c_str());
     202           8 :     if (m_bitDepth > 0)
     203             :     {
     204           1 :         CPLCreateXMLElementAndValue(psPansharpeningOptionsNode, "BitDepth",
     205             :                                     CPLSPrintf("%d", m_bitDepth));
     206             :     }
     207           8 :     if (GetArg("nodata")->IsExplicitlySet())
     208             :     {
     209           1 :         CPLCreateXMLElementAndValue(psPansharpeningOptionsNode, "NoData",
     210             :                                     CPLSPrintf("%.17g", m_nodata));
     211             :     }
     212           8 :     CPLCreateXMLElementAndValue(
     213             :         psPansharpeningOptionsNode, "SpatialExtentAdjustment",
     214          16 :         CPLString(m_spatialExtentAdjustment).replaceAll('-', 0).c_str());
     215          32 :     for (int i = 0; i < static_cast<int>(ahSpectralBands.size()); ++i)
     216             :     {
     217          24 :         auto psSpectraBandNode = CPLCreateXMLNode(psPansharpeningOptionsNode,
     218             :                                                   CXT_Element, "SpectralBand");
     219          24 :         CPLAddXMLAttributeAndValue(psSpectraBandNode, "dstBand",
     220             :                                    CPLSPrintf("%d", i + 1));
     221             :     }
     222          16 :     CPLCharUniquePtr pszXML(CPLSerializeXMLTree(root.get()));
     223             :     auto poVRTDS = std::unique_ptr<GDALDataset>(
     224             :         GDALDataset::FromHandle(GDALCreatePansharpenedVRT(
     225           8 :             pszXML.get(), GDALRasterBand::ToHandle(poPanDS->GetRasterBand(1)),
     226          16 :             static_cast<int>(ahSpectralBands.size()), ahSpectralBands.data())));
     227           8 :     const bool bRet = poVRTDS != nullptr;
     228           8 :     if (poVRTDS)
     229             :     {
     230           7 :         m_outputDataset.Set(std::move(poVRTDS));
     231             :     }
     232           8 :     return bRet;
     233             : }
     234             : 
     235             : GDALRasterPansharpenAlgorithmStandalone::
     236             :     ~GDALRasterPansharpenAlgorithmStandalone() = default;
     237             : 
     238             : //! @endcond