Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  "reclassify" step of "raster pipeline"
       5             :  * Author:   Daniel Baston
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2025, ISciences LLC
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "gdalalg_raster_reclassify.h"
      14             : 
      15             : #include "cpl_vsi_virtual.h"
      16             : #include "gdal_priv.h"
      17             : #include "gdal_utils.h"
      18             : #include "../frmts/vrt/vrtdataset.h"
      19             : #include "../frmts/vrt/vrtreclassifier.h"
      20             : 
      21             : #include <array>
      22             : 
      23             : //! @cond Doxygen_Suppress
      24             : 
      25             : #ifndef _
      26             : #define _(x) (x)
      27             : #endif
      28             : 
      29             : /************************************************************************/
      30             : /*    GDALRasterReclassifyAlgorithm::GDALRasterReclassifyAlgorithm()    */
      31             : /************************************************************************/
      32             : 
      33          99 : GDALRasterReclassifyAlgorithm::GDALRasterReclassifyAlgorithm(
      34          99 :     bool standaloneStep)
      35             :     : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
      36          99 :                                       standaloneStep)
      37             : {
      38             :     AddArg("mapping", 'm',
      39             :            _("Reclassification mappings (or specify a @<filename> to point to "
      40             :              "a file containing mappings"),
      41         198 :            &m_mapping)
      42          99 :         .SetRequired();
      43             :     AddArg("keep-color-table", 0, _("Preserve the input color table"),
      44          99 :            &m_keepColorTable);
      45          99 :     AddOutputDataTypeArg(&m_type);
      46          99 : }
      47             : 
      48             : /************************************************************************/
      49             : /*                 GDALRasterReclassifyValidateMappings                 */
      50             : /************************************************************************/
      51             : 
      52          14 : static bool GDALReclassifyValidateMappings(GDALDataset &input,
      53             :                                            const std::string &mappings,
      54             :                                            GDALDataType eDstType)
      55             : {
      56             :     int hasNoData;
      57             :     std::optional<double> noData =
      58          14 :         input.GetRasterBand(1)->GetNoDataValue(&hasNoData);
      59          14 :     if (!hasNoData)
      60             :     {
      61          10 :         noData.reset();
      62             :     }
      63             : 
      64          14 :     gdal::Reclassifier reclassifier;
      65          28 :     return reclassifier.Init(mappings.c_str(), noData, eDstType) == CE_None;
      66             : }
      67             : 
      68             : /************************************************************************/
      69             : /*                 GDALRasterReclassifyCreateVRTDerived                 */
      70             : /************************************************************************/
      71             : 
      72             : static std::unique_ptr<GDALDataset>
      73          12 : GDALReclassifyCreateVRTDerived(GDALDataset &input, const std::string &mappings,
      74             :                                GDALDataType eDstType, bool keepColorTable)
      75             : {
      76          12 :     const auto nX = input.GetRasterXSize();
      77          12 :     const auto nY = input.GetRasterYSize();
      78             : 
      79          24 :     auto poDS = VRTDataset::CreateVRTDataset("", nX, nY, 0, eDstType, nullptr);
      80             : 
      81          12 :     GDALGeoTransform gt;
      82          12 :     if (input.GetGeoTransform(gt) == CE_None)
      83           8 :         poDS->SetGeoTransform(gt);
      84          12 :     poDS->SetSpatialRef(input.GetSpatialRef());
      85             : 
      86          24 :     CPLStringList papszBandArgs;
      87          12 :     papszBandArgs.SetNameValue("subclass", "VRTDerivedRasterBand");
      88             : 
      89          26 :     for (int iBand = 1; iBand <= input.GetRasterCount(); ++iBand)
      90             :     {
      91          14 :         GDALRasterBand *poSrcBand = input.GetRasterBand(iBand);
      92          14 :         const GDALDataType srcType = poSrcBand->GetRasterDataType();
      93          14 :         const GDALDataType bandType =
      94          14 :             eDstType == GDT_Unknown ? srcType : eDstType;
      95          14 :         const GDALDataType xferType = GDALDataTypeUnion(srcType, bandType);
      96             : 
      97          14 :         if (poDS->AddBand(bandType, papszBandArgs) != CE_None)
      98             :         {
      99           0 :             return nullptr;
     100             :         }
     101             : 
     102             :         VRTDerivedRasterBand *poDstBand =
     103          14 :             cpl::down_cast<VRTDerivedRasterBand *>(poDS->GetRasterBand(iBand));
     104          14 :         poDstBand->SetSourceTransferType(xferType);
     105          14 :         poDstBand->SetPixelFunctionName("reclassify");
     106          14 :         poDstBand->AddPixelFunctionArgument("mapping", mappings.c_str());
     107             : 
     108          14 :         if (keepColorTable && poSrcBand->GetColorTable() != nullptr)
     109             :         {
     110           1 :             poDstBand->SetColorTable(poSrcBand->GetColorTable());
     111             :         }
     112             : 
     113          14 :         GDALCopyNoDataValue(poDstBand, poSrcBand);
     114          14 :         poDstBand->AddSimpleSource(poSrcBand);
     115             :     }
     116             : 
     117          12 :     return poDS;
     118             : }
     119             : 
     120             : /************************************************************************/
     121             : /*               GDALRasterReclassifyAlgorithm::RunStep()               */
     122             : /************************************************************************/
     123             : 
     124          17 : bool GDALRasterReclassifyAlgorithm::RunStep(GDALPipelineStepRunContext &)
     125             : {
     126          17 :     const auto poSrcDS = m_inputDataset[0].GetDatasetRef();
     127          17 :     CPLAssert(poSrcDS);
     128          17 :     CPLAssert(m_outputDataset.GetName().empty());
     129          17 :     CPLAssert(!m_outputDataset.GetDatasetRef());
     130             : 
     131             :     // Already validated by argument parser
     132             :     const GDALDataType eDstType =
     133          17 :         m_type.empty() ? GDT_Unknown : GDALGetDataTypeByName(m_type.c_str());
     134             : 
     135          17 :     const auto nErrorCount = CPLGetErrorCounter();
     136          17 :     if (!m_mapping.empty() && m_mapping[0] == '@')
     137             :     {
     138             :         auto f =
     139           6 :             VSIVirtualHandleUniquePtr(VSIFOpenL(m_mapping.c_str() + 1, "r"));
     140           6 :         if (!f)
     141             :         {
     142           1 :             ReportError(CE_Failure, CPLE_FileIO, "Cannot open %s",
     143           1 :                         m_mapping.c_str() + 1);
     144           1 :             return false;
     145             :         }
     146             : 
     147           5 :         m_mapping.clear();
     148             :         try
     149             :         {
     150           5 :             constexpr int MAX_CHARS_PER_LINE = 1000 * 1000;
     151           5 :             constexpr size_t MAX_MAPPING_SIZE = 10 * 1000 * 1000;
     152             :             while (const char *line =
     153        1034 :                        CPLReadLine2L(f.get(), MAX_CHARS_PER_LINE, nullptr))
     154             :             {
     155        1200 :                 while (isspace(*line))
     156             :                 {
     157         170 :                     line++;
     158             :                 }
     159             : 
     160        1030 :                 if (line[0])
     161             :                 {
     162        1024 :                     if (!m_mapping.empty())
     163             :                     {
     164        1018 :                         m_mapping.append(";");
     165             :                     }
     166             : 
     167        1024 :                     const char *comment = strchr(line, '#');
     168        1024 :                     if (!comment)
     169             :                     {
     170        1018 :                         m_mapping.append(line);
     171             :                     }
     172             :                     else
     173             :                     {
     174             :                         m_mapping.append(line,
     175           6 :                                          static_cast<size_t>(comment - line));
     176             :                     }
     177        1024 :                     if (m_mapping.size() > MAX_MAPPING_SIZE)
     178             :                     {
     179           1 :                         ReportError(CE_Failure, CPLE_AppDefined,
     180             :                                     "Too large mapping size");
     181           1 :                         return false;
     182             :                     }
     183             :                 }
     184        1029 :             }
     185             :         }
     186           0 :         catch (const std::exception &)
     187             :         {
     188           0 :             ReportError(CE_Failure, CPLE_OutOfMemory,
     189             :                         "Out of memory while ingesting mapping file");
     190             :         }
     191             :     }
     192          15 :     if (nErrorCount == CPLGetErrorCounter())
     193             :     {
     194          14 :         if (!GDALReclassifyValidateMappings(*poSrcDS, m_mapping, eDstType))
     195             :         {
     196           2 :             return false;
     197             :         }
     198             : 
     199          12 :         m_outputDataset.Set(GDALReclassifyCreateVRTDerived(
     200          12 :             *poSrcDS, m_mapping, eDstType, m_keepColorTable));
     201             :     }
     202          13 :     return m_outputDataset.GetDatasetRef() != nullptr;
     203             : }
     204             : 
     205             : GDALRasterReclassifyAlgorithmStandalone::
     206             :     ~GDALRasterReclassifyAlgorithmStandalone() = default;
     207             : 
     208             : //! @endcond