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          97 : GDALRasterReclassifyAlgorithm::GDALRasterReclassifyAlgorithm(
      34          97 :     bool standaloneStep)
      35             :     : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
      36          97 :                                       standaloneStep)
      37             : {
      38             :     AddArg("mapping", 'm',
      39             :            _("Reclassification mappings (or specify a @<filename> to point to "
      40             :              "a file containing mappings"),
      41         194 :            &m_mapping)
      42          97 :         .SetRequired();
      43             :     AddArg("keep-color-table", 0, _("Preserve the input color table"),
      44          97 :            &m_keepColorTable);
      45          97 :     AddOutputDataTypeArg(&m_type);
      46          97 : }
      47             : 
      48             : /************************************************************************/
      49             : /*                 GDALRasterReclassifyValidateMappings                 */
      50             : /************************************************************************/
      51             : 
      52          12 : static bool GDALReclassifyValidateMappings(GDALDataset &input,
      53             :                                            const std::string &mappings,
      54             :                                            GDALDataType eDstType)
      55             : {
      56             :     int hasNoData;
      57             :     std::optional<double> noData =
      58          12 :         input.GetRasterBand(1)->GetNoDataValue(&hasNoData);
      59          12 :     if (!hasNoData)
      60             :     {
      61           8 :         noData.reset();
      62             :     }
      63             : 
      64          12 :     gdal::Reclassifier reclassifier;
      65          24 :     return reclassifier.Init(mappings.c_str(), noData, eDstType) == CE_None;
      66             : }
      67             : 
      68             : /************************************************************************/
      69             : /*                 GDALRasterReclassifyCreateVRTDerived                 */
      70             : /************************************************************************/
      71             : 
      72             : static std::unique_ptr<GDALDataset>
      73          11 : GDALReclassifyCreateVRTDerived(GDALDataset &input, const std::string &mappings,
      74             :                                GDALDataType eDstType, bool keepColorTable)
      75             : {
      76          22 :     CPLXMLTreeCloser root(CPLCreateXMLNode(nullptr, CXT_Element, "VRTDataset"));
      77             : 
      78          11 :     const auto nX = input.GetRasterXSize();
      79          11 :     const auto nY = input.GetRasterYSize();
      80             : 
      81          24 :     for (int iBand = 1; iBand <= input.GetRasterCount(); ++iBand)
      82             :     {
      83          13 :         GDALRasterBand *poSrcBand = input.GetRasterBand(iBand);
      84          13 :         const GDALDataType srcType = poSrcBand->GetRasterDataType();
      85          13 :         const GDALDataType bandType =
      86          13 :             eDstType == GDT_Unknown ? srcType : eDstType;
      87          13 :         const GDALDataType xferType = GDALDataTypeUnion(srcType, bandType);
      88             : 
      89             :         CPLXMLNode *band =
      90          13 :             CPLCreateXMLNode(root.get(), CXT_Element, "VRTRasterBand");
      91          13 :         CPLAddXMLAttributeAndValue(band, "subClass", "VRTDerivedRasterBand");
      92          13 :         CPLAddXMLAttributeAndValue(band, "dataType",
      93             :                                    GDALGetDataTypeName(bandType));
      94             : 
      95             :         CPLXMLNode *pixelFunctionType =
      96          13 :             CPLCreateXMLNode(band, CXT_Element, "PixelFunctionType");
      97          13 :         CPLCreateXMLNode(pixelFunctionType, CXT_Text, "reclassify");
      98             :         CPLXMLNode *arguments =
      99          13 :             CPLCreateXMLNode(band, CXT_Element, "PixelFunctionArguments");
     100          13 :         CPLAddXMLAttributeAndValue(arguments, "mapping", mappings.c_str());
     101             : 
     102             :         CPLXMLNode *sourceTransferType =
     103          13 :             CPLCreateXMLNode(band, CXT_Element, "SourceTransferType");
     104          13 :         CPLCreateXMLNode(sourceTransferType, CXT_Text,
     105             :                          GDALGetDataTypeName(xferType));
     106             : 
     107          13 :         if (const GDALColorTable *poTable =
     108          13 :                 keepColorTable ? poSrcBand->GetColorTable() : nullptr)
     109             :         {
     110             :             CPLXMLNode *colorTable =
     111           1 :                 CPLCreateXMLNode(band, CXT_Element, "ColorTable");
     112           5 :             for (int i = 0; i < poTable->GetColorEntryCount(); ++i)
     113             :             {
     114           4 :                 const GDALColorEntry *poSrcEntry = poTable->GetColorEntry(i);
     115             : 
     116             :                 CPLXMLNode *entry =
     117           4 :                     CPLCreateXMLNode(colorTable, CXT_Element, "Entry");
     118           4 :                 CPLAddXMLAttributeAndValue(
     119           8 :                     entry, "c1", std::to_string(poSrcEntry->c1).c_str());
     120           4 :                 CPLAddXMLAttributeAndValue(
     121           8 :                     entry, "c2", std::to_string(poSrcEntry->c2).c_str());
     122           4 :                 CPLAddXMLAttributeAndValue(
     123           8 :                     entry, "c3", std::to_string(poSrcEntry->c3).c_str());
     124           4 :                 CPLAddXMLAttributeAndValue(
     125           8 :                     entry, "c4", std::to_string(poSrcEntry->c4).c_str());
     126             :             }
     127             :         }
     128             :     }
     129             : 
     130          22 :     auto ds = std::make_unique<VRTDataset>(nX, nY);
     131          11 :     if (ds->XMLInit(root.get(), "") != CE_None)
     132             :     {
     133           0 :         return nullptr;
     134             :     };
     135             : 
     136          24 :     for (int iBand = 1; iBand <= input.GetRasterCount(); ++iBand)
     137             :     {
     138          13 :         auto poSrcBand = input.GetRasterBand(iBand);
     139             :         auto poDstBand =
     140          13 :             cpl::down_cast<VRTDerivedRasterBand *>(ds->GetRasterBand(iBand));
     141          13 :         GDALCopyNoDataValue(poDstBand, poSrcBand);
     142          13 :         poDstBand->AddSimpleSource(poSrcBand);
     143             :     }
     144             : 
     145          11 :     GDALGeoTransform gt;
     146          11 :     if (input.GetGeoTransform(gt) == CE_None)
     147           7 :         ds->SetGeoTransform(gt);
     148          11 :     ds->SetSpatialRef(input.GetSpatialRef());
     149             : 
     150          11 :     return ds;
     151             : }
     152             : 
     153             : /************************************************************************/
     154             : /*               GDALRasterReclassifyAlgorithm::RunStep()               */
     155             : /************************************************************************/
     156             : 
     157          15 : bool GDALRasterReclassifyAlgorithm::RunStep(GDALPipelineStepRunContext &)
     158             : {
     159          15 :     const auto poSrcDS = m_inputDataset[0].GetDatasetRef();
     160          15 :     CPLAssert(poSrcDS);
     161          15 :     CPLAssert(m_outputDataset.GetName().empty());
     162          15 :     CPLAssert(!m_outputDataset.GetDatasetRef());
     163             : 
     164             :     // Already validated by argument parser
     165             :     const GDALDataType eDstType =
     166          15 :         m_type.empty() ? GDT_Unknown : GDALGetDataTypeByName(m_type.c_str());
     167             : 
     168          15 :     const auto nErrorCount = CPLGetErrorCounter();
     169          15 :     if (!m_mapping.empty() && m_mapping[0] == '@')
     170             :     {
     171             :         auto f =
     172           6 :             VSIVirtualHandleUniquePtr(VSIFOpenL(m_mapping.c_str() + 1, "r"));
     173           6 :         if (!f)
     174             :         {
     175           1 :             ReportError(CE_Failure, CPLE_FileIO, "Cannot open %s",
     176           1 :                         m_mapping.c_str() + 1);
     177           1 :             return false;
     178             :         }
     179             : 
     180           5 :         m_mapping.clear();
     181             :         try
     182             :         {
     183           5 :             constexpr int MAX_CHARS_PER_LINE = 1000 * 1000;
     184           5 :             constexpr size_t MAX_MAPPING_SIZE = 10 * 1000 * 1000;
     185             :             while (const char *line =
     186        1034 :                        CPLReadLine2L(f.get(), MAX_CHARS_PER_LINE, nullptr))
     187             :             {
     188        1200 :                 while (isspace(*line))
     189             :                 {
     190         170 :                     line++;
     191             :                 }
     192             : 
     193        1030 :                 if (line[0])
     194             :                 {
     195        1024 :                     if (!m_mapping.empty())
     196             :                     {
     197        1018 :                         m_mapping.append(";");
     198             :                     }
     199             : 
     200        1024 :                     const char *comment = strchr(line, '#');
     201        1024 :                     if (!comment)
     202             :                     {
     203        1018 :                         m_mapping.append(line);
     204             :                     }
     205             :                     else
     206             :                     {
     207             :                         m_mapping.append(line,
     208           6 :                                          static_cast<size_t>(comment - line));
     209             :                     }
     210        1024 :                     if (m_mapping.size() > MAX_MAPPING_SIZE)
     211             :                     {
     212           1 :                         ReportError(CE_Failure, CPLE_AppDefined,
     213             :                                     "Too large mapping size");
     214           1 :                         return false;
     215             :                     }
     216             :                 }
     217        1029 :             }
     218             :         }
     219           0 :         catch (const std::exception &)
     220             :         {
     221           0 :             ReportError(CE_Failure, CPLE_OutOfMemory,
     222             :                         "Out of memory while ingesting mapping file");
     223             :         }
     224             :     }
     225          13 :     if (nErrorCount == CPLGetErrorCounter())
     226             :     {
     227          12 :         if (!GDALReclassifyValidateMappings(*poSrcDS, m_mapping, eDstType))
     228             :         {
     229           1 :             return false;
     230             :         }
     231             : 
     232          11 :         m_outputDataset.Set(GDALReclassifyCreateVRTDerived(
     233          11 :             *poSrcDS, m_mapping, eDstType, m_keepColorTable));
     234             :     }
     235          12 :     return m_outputDataset.GetDatasetRef() != nullptr;
     236             : }
     237             : 
     238             : GDALRasterReclassifyAlgorithmStandalone::
     239             :     ~GDALRasterReclassifyAlgorithmStandalone() = default;
     240             : 
     241             : //! @endcond