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