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 9 : GDALGeoTransform gt;
122 9 : if (input.GetGeoTransform(gt) == CE_None)
123 7 : ds->SetGeoTransform(gt);
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
|
