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