LCOV - code coverage report
Current view: top level - apps - gdalalg_raster_reclassify.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 78 81 96.3 %
Date: 2025-06-19 12:30:01 Functions: 4 4 100.0 %

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

Generated by: LCOV version 1.14