/src/gdal/gcore/gdalcomputedrasterband.cpp
Line | Count | Source (jump to first uncovered line) |
1 | | /****************************************************************************** |
2 | | * |
3 | | * Project: GDAL |
4 | | * Purpose: GDALComputedDataset and GDALComputedRasterBand |
5 | | * Author: Even Rouault <even dot rouault at spatialys.com> |
6 | | * |
7 | | ****************************************************************************** |
8 | | * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com> |
9 | | * |
10 | | * SPDX-License-Identifier: MIT |
11 | | ****************************************************************************/ |
12 | | |
13 | | #include "gdal_priv.h" |
14 | | #include "vrtdataset.h" |
15 | | |
16 | | #include <cmath> |
17 | | #include <limits> |
18 | | |
19 | | /************************************************************************/ |
20 | | /* GDALComputedDataset */ |
21 | | /************************************************************************/ |
22 | | |
23 | | class GDALComputedDataset final : public GDALDataset |
24 | | { |
25 | | friend class GDALComputedRasterBand; |
26 | | |
27 | | const GDALComputedRasterBand::Operation m_op; |
28 | | CPLStringList m_aosOptions{}; |
29 | | std::vector<std::unique_ptr<GDALDataset, GDALDatasetUniquePtrReleaser>> |
30 | | m_bandDS{}; |
31 | | std::vector<GDALRasterBand *> m_poBands{}; |
32 | | VRTDataset m_oVRTDS; |
33 | | |
34 | | void AddSources(GDALComputedRasterBand *poBand); |
35 | | |
36 | | static const char * |
37 | | OperationToFunctionName(GDALComputedRasterBand::Operation op); |
38 | | |
39 | | GDALComputedDataset &operator=(const GDALComputedDataset &) = delete; |
40 | | GDALComputedDataset(GDALComputedDataset &&) = delete; |
41 | | GDALComputedDataset &operator=(GDALComputedDataset &&) = delete; |
42 | | |
43 | | public: |
44 | | GDALComputedDataset(const GDALComputedDataset &other); |
45 | | |
46 | | GDALComputedDataset(GDALComputedRasterBand *poBand, int nXSize, int nYSize, |
47 | | GDALDataType eDT, int nBlockXSize, int nBlockYSize, |
48 | | GDALComputedRasterBand::Operation op, |
49 | | const GDALRasterBand *firstBand, double *pFirstConstant, |
50 | | const GDALRasterBand *secondBand, |
51 | | double *pSecondConstant); |
52 | | |
53 | | GDALComputedDataset(GDALComputedRasterBand *poBand, int nXSize, int nYSize, |
54 | | GDALDataType eDT, int nBlockXSize, int nBlockYSize, |
55 | | GDALComputedRasterBand::Operation op, |
56 | | const std::vector<const GDALRasterBand *> &bands, |
57 | | double constant); |
58 | | |
59 | | ~GDALComputedDataset() override; |
60 | | |
61 | | CPLErr GetGeoTransform(double *padfGeoTransform) override |
62 | 0 | { |
63 | 0 | return m_oVRTDS.GetGeoTransform(padfGeoTransform); |
64 | 0 | } |
65 | | |
66 | | const OGRSpatialReference *GetSpatialRef() const override |
67 | 0 | { |
68 | 0 | return m_oVRTDS.GetSpatialRef(); |
69 | 0 | } |
70 | | |
71 | | char **GetMetadata(const char *pszDomain) override |
72 | 0 | { |
73 | 0 | return m_oVRTDS.GetMetadata(pszDomain); |
74 | 0 | } |
75 | | |
76 | | const char *GetMetadataItem(const char *pszName, |
77 | | const char *pszDomain) override |
78 | 0 | { |
79 | 0 | return m_oVRTDS.GetMetadataItem(pszName, pszDomain); |
80 | 0 | } |
81 | | |
82 | | void *GetInternalHandle(const char *pszHandleName) override |
83 | 0 | { |
84 | 0 | if (pszHandleName && EQUAL(pszHandleName, "VRT_DATASET")) |
85 | 0 | return &m_oVRTDS; |
86 | 0 | return nullptr; |
87 | 0 | } |
88 | | }; |
89 | | |
90 | | /************************************************************************/ |
91 | | /* IsComparisonOperator() */ |
92 | | /************************************************************************/ |
93 | | |
94 | | static bool IsComparisonOperator(GDALComputedRasterBand::Operation op) |
95 | 0 | { |
96 | 0 | switch (op) |
97 | 0 | { |
98 | 0 | case GDALComputedRasterBand::Operation::OP_GT: |
99 | 0 | case GDALComputedRasterBand::Operation::OP_GE: |
100 | 0 | case GDALComputedRasterBand::Operation::OP_LT: |
101 | 0 | case GDALComputedRasterBand::Operation::OP_LE: |
102 | 0 | case GDALComputedRasterBand::Operation::OP_EQ: |
103 | 0 | case GDALComputedRasterBand::Operation::OP_NE: |
104 | 0 | case GDALComputedRasterBand::Operation::OP_LOGICAL_AND: |
105 | 0 | case GDALComputedRasterBand::Operation::OP_LOGICAL_OR: |
106 | 0 | return true; |
107 | 0 | default: |
108 | 0 | break; |
109 | 0 | } |
110 | 0 | return false; |
111 | 0 | } |
112 | | |
113 | | /************************************************************************/ |
114 | | /* GDALComputedDataset() */ |
115 | | /************************************************************************/ |
116 | | |
117 | | GDALComputedDataset::GDALComputedDataset(const GDALComputedDataset &other) |
118 | 0 | : GDALDataset(), m_op(other.m_op), m_aosOptions(other.m_aosOptions), |
119 | 0 | m_poBands(other.m_poBands), |
120 | 0 | m_oVRTDS(other.GetRasterXSize(), other.GetRasterYSize(), |
121 | 0 | other.m_oVRTDS.GetBlockXSize(), other.m_oVRTDS.GetBlockYSize()) |
122 | 0 | { |
123 | 0 | nRasterXSize = other.nRasterXSize; |
124 | 0 | nRasterYSize = other.nRasterYSize; |
125 | |
|
126 | 0 | auto poBand = new GDALComputedRasterBand( |
127 | 0 | const_cast<const GDALComputedRasterBand &>( |
128 | 0 | *cpl::down_cast<GDALComputedRasterBand *>( |
129 | 0 | const_cast<GDALComputedDataset &>(other).GetRasterBand(1))), |
130 | 0 | true); |
131 | 0 | SetBand(1, poBand); |
132 | |
|
133 | 0 | double adfGT[6]; |
134 | 0 | if (const_cast<VRTDataset &>(other.m_oVRTDS).GetGeoTransform(adfGT) == |
135 | 0 | CE_None) |
136 | 0 | { |
137 | 0 | m_oVRTDS.SetGeoTransform(adfGT); |
138 | 0 | } |
139 | |
|
140 | 0 | if (const auto *poSRS = |
141 | 0 | const_cast<VRTDataset &>(other.m_oVRTDS).GetSpatialRef()) |
142 | 0 | { |
143 | 0 | m_oVRTDS.SetSpatialRef(poSRS); |
144 | 0 | } |
145 | |
|
146 | 0 | m_oVRTDS.AddBand(other.m_oVRTDS.GetRasterBand(1)->GetRasterDataType(), |
147 | 0 | m_aosOptions.List()); |
148 | |
|
149 | 0 | AddSources(poBand); |
150 | 0 | } |
151 | | |
152 | | /************************************************************************/ |
153 | | /* GDALComputedDataset() */ |
154 | | /************************************************************************/ |
155 | | |
156 | | GDALComputedDataset::GDALComputedDataset( |
157 | | GDALComputedRasterBand *poBand, int nXSize, int nYSize, GDALDataType eDT, |
158 | | int nBlockXSize, int nBlockYSize, GDALComputedRasterBand::Operation op, |
159 | | const GDALRasterBand *firstBand, double *pFirstConstant, |
160 | | const GDALRasterBand *secondBand, double *pSecondConstant) |
161 | 0 | : m_op(op), m_oVRTDS(nXSize, nYSize, nBlockXSize, nBlockYSize) |
162 | 0 | { |
163 | 0 | CPLAssert(firstBand != nullptr || secondBand != nullptr); |
164 | 0 | if (firstBand) |
165 | 0 | m_poBands.push_back(const_cast<GDALRasterBand *>(firstBand)); |
166 | 0 | if (secondBand) |
167 | 0 | m_poBands.push_back(const_cast<GDALRasterBand *>(secondBand)); |
168 | |
|
169 | 0 | nRasterXSize = nXSize; |
170 | 0 | nRasterYSize = nYSize; |
171 | |
|
172 | 0 | if (auto poSrcDS = m_poBands.front()->GetDataset()) |
173 | 0 | { |
174 | 0 | double adfGT[6]; |
175 | 0 | if (poSrcDS->GetGeoTransform(adfGT) == CE_None) |
176 | 0 | { |
177 | 0 | m_oVRTDS.SetGeoTransform(adfGT); |
178 | 0 | } |
179 | |
|
180 | 0 | if (const auto *poSRS = poSrcDS->GetSpatialRef()) |
181 | 0 | { |
182 | 0 | m_oVRTDS.SetSpatialRef(poSRS); |
183 | 0 | } |
184 | 0 | } |
185 | |
|
186 | 0 | if (op == GDALComputedRasterBand::Operation::OP_CAST) |
187 | 0 | { |
188 | 0 | #ifdef DEBUG |
189 | | // Just for code coverage... |
190 | 0 | CPL_IGNORE_RET_VAL(GDALComputedDataset::OperationToFunctionName(op)); |
191 | 0 | #endif |
192 | 0 | m_aosOptions.SetNameValue("subclass", "VRTSourcedRasterBand"); |
193 | 0 | } |
194 | 0 | else |
195 | 0 | { |
196 | 0 | m_aosOptions.SetNameValue("subclass", "VRTDerivedRasterBand"); |
197 | 0 | if (IsComparisonOperator(op)) |
198 | 0 | { |
199 | 0 | m_aosOptions.SetNameValue("PixelFunctionType", "expression"); |
200 | 0 | if (firstBand && secondBand) |
201 | 0 | { |
202 | 0 | m_aosOptions.SetNameValue( |
203 | 0 | "_PIXELFN_ARG_expression", |
204 | 0 | CPLSPrintf( |
205 | 0 | "source1 %s source2", |
206 | 0 | GDALComputedDataset::OperationToFunctionName(op))); |
207 | 0 | } |
208 | 0 | else if (firstBand && pSecondConstant) |
209 | 0 | { |
210 | 0 | m_aosOptions.SetNameValue( |
211 | 0 | "_PIXELFN_ARG_expression", |
212 | 0 | CPLSPrintf("source1 %s %.17g", |
213 | 0 | GDALComputedDataset::OperationToFunctionName(op), |
214 | 0 | *pSecondConstant)); |
215 | 0 | } |
216 | 0 | else if (pFirstConstant && secondBand) |
217 | 0 | { |
218 | 0 | m_aosOptions.SetNameValue( |
219 | 0 | "_PIXELFN_ARG_expression", |
220 | 0 | CPLSPrintf( |
221 | 0 | "%.17g %s source1", *pFirstConstant, |
222 | 0 | GDALComputedDataset::OperationToFunctionName(op))); |
223 | 0 | } |
224 | 0 | else |
225 | 0 | { |
226 | 0 | CPLAssert(false); |
227 | 0 | } |
228 | 0 | } |
229 | 0 | else if (op == GDALComputedRasterBand::Operation::OP_SUBTRACT && |
230 | 0 | pSecondConstant) |
231 | 0 | { |
232 | 0 | m_aosOptions.SetNameValue("PixelFunctionType", "sum"); |
233 | 0 | m_aosOptions.SetNameValue("_PIXELFN_ARG_k", |
234 | 0 | CPLSPrintf("%.17g", -(*pSecondConstant))); |
235 | 0 | } |
236 | 0 | else if (op == GDALComputedRasterBand::Operation::OP_DIVIDE) |
237 | 0 | { |
238 | 0 | if (pSecondConstant) |
239 | 0 | { |
240 | 0 | m_aosOptions.SetNameValue("PixelFunctionType", "mul"); |
241 | 0 | m_aosOptions.SetNameValue( |
242 | 0 | "_PIXELFN_ARG_k", |
243 | 0 | CPLSPrintf("%.17g", 1.0 / (*pSecondConstant))); |
244 | 0 | } |
245 | 0 | else if (pFirstConstant) |
246 | 0 | { |
247 | 0 | m_aosOptions.SetNameValue("PixelFunctionType", "inv"); |
248 | 0 | m_aosOptions.SetNameValue("_PIXELFN_ARG_k", |
249 | 0 | CPLSPrintf("%.17g", *pFirstConstant)); |
250 | 0 | } |
251 | 0 | else |
252 | 0 | { |
253 | 0 | m_aosOptions.SetNameValue("PixelFunctionType", "div"); |
254 | 0 | } |
255 | 0 | } |
256 | 0 | else |
257 | 0 | { |
258 | 0 | m_aosOptions.SetNameValue("PixelFunctionType", |
259 | 0 | OperationToFunctionName(op)); |
260 | 0 | if (pSecondConstant) |
261 | 0 | m_aosOptions.SetNameValue( |
262 | 0 | "_PIXELFN_ARG_k", CPLSPrintf("%.17g", *pSecondConstant)); |
263 | 0 | } |
264 | 0 | } |
265 | 0 | m_aosOptions.SetNameValue("_PIXELFN_ARG_propagateNoData", "true"); |
266 | 0 | m_oVRTDS.AddBand(eDT, m_aosOptions.List()); |
267 | |
|
268 | 0 | SetBand(1, poBand); |
269 | |
|
270 | 0 | AddSources(poBand); |
271 | 0 | } |
272 | | |
273 | | /************************************************************************/ |
274 | | /* GDALComputedDataset() */ |
275 | | /************************************************************************/ |
276 | | |
277 | | GDALComputedDataset::GDALComputedDataset( |
278 | | GDALComputedRasterBand *poBand, int nXSize, int nYSize, GDALDataType eDT, |
279 | | int nBlockXSize, int nBlockYSize, GDALComputedRasterBand::Operation op, |
280 | | const std::vector<const GDALRasterBand *> &bands, double constant) |
281 | 0 | : m_op(op), m_oVRTDS(nXSize, nYSize, nBlockXSize, nBlockYSize) |
282 | 0 | { |
283 | 0 | for (const GDALRasterBand *poIterBand : bands) |
284 | 0 | m_poBands.push_back(const_cast<GDALRasterBand *>(poIterBand)); |
285 | |
|
286 | 0 | nRasterXSize = nXSize; |
287 | 0 | nRasterYSize = nYSize; |
288 | |
|
289 | 0 | if (auto poSrcDS = m_poBands.front()->GetDataset()) |
290 | 0 | { |
291 | 0 | double adfGT[6]; |
292 | 0 | if (poSrcDS->GetGeoTransform(adfGT) == CE_None) |
293 | 0 | { |
294 | 0 | m_oVRTDS.SetGeoTransform(adfGT); |
295 | 0 | } |
296 | |
|
297 | 0 | if (const auto *poSRS = poSrcDS->GetSpatialRef()) |
298 | 0 | { |
299 | 0 | m_oVRTDS.SetSpatialRef(poSRS); |
300 | 0 | } |
301 | 0 | } |
302 | |
|
303 | 0 | m_aosOptions.SetNameValue("subclass", "VRTDerivedRasterBand"); |
304 | 0 | if (op == GDALComputedRasterBand::Operation::OP_TERNARY) |
305 | 0 | { |
306 | 0 | m_aosOptions.SetNameValue("PixelFunctionType", "expression"); |
307 | 0 | m_aosOptions.SetNameValue("_PIXELFN_ARG_expression", |
308 | 0 | "source1 ? source2 : source3"); |
309 | 0 | } |
310 | 0 | else |
311 | 0 | { |
312 | 0 | m_aosOptions.SetNameValue("PixelFunctionType", |
313 | 0 | OperationToFunctionName(op)); |
314 | 0 | if (!std::isnan(constant)) |
315 | 0 | { |
316 | 0 | m_aosOptions.SetNameValue("_PIXELFN_ARG_k", |
317 | 0 | CPLSPrintf("%.17g", constant)); |
318 | 0 | } |
319 | 0 | m_aosOptions.SetNameValue("_PIXELFN_ARG_propagateNoData", "true"); |
320 | 0 | } |
321 | 0 | m_oVRTDS.AddBand(eDT, m_aosOptions.List()); |
322 | |
|
323 | 0 | SetBand(1, poBand); |
324 | |
|
325 | 0 | AddSources(poBand); |
326 | 0 | } |
327 | | |
328 | | /************************************************************************/ |
329 | | /* ~GDALComputedDataset() */ |
330 | | /************************************************************************/ |
331 | | |
332 | 0 | GDALComputedDataset::~GDALComputedDataset() = default; |
333 | | |
334 | | /************************************************************************/ |
335 | | /* HaveAllBandsSameNoDataValue() */ |
336 | | /************************************************************************/ |
337 | | |
338 | | static bool HaveAllBandsSameNoDataValue(GDALRasterBand **apoBands, |
339 | | size_t nBands, bool &hasAtLeastOneNDV, |
340 | | double &singleNDV) |
341 | 0 | { |
342 | 0 | hasAtLeastOneNDV = false; |
343 | 0 | singleNDV = 0; |
344 | |
|
345 | 0 | int bFirstBandHasNoData = false; |
346 | 0 | for (size_t i = 0; i < nBands; ++i) |
347 | 0 | { |
348 | 0 | int bHasNoData = false; |
349 | 0 | const double dfNoData = apoBands[i]->GetNoDataValue(&bHasNoData); |
350 | 0 | if (bHasNoData) |
351 | 0 | hasAtLeastOneNDV = true; |
352 | 0 | if (i == 0) |
353 | 0 | { |
354 | 0 | bFirstBandHasNoData = bHasNoData; |
355 | 0 | singleNDV = dfNoData; |
356 | 0 | } |
357 | 0 | else if (bHasNoData != bFirstBandHasNoData) |
358 | 0 | { |
359 | 0 | return false; |
360 | 0 | } |
361 | 0 | else if (bFirstBandHasNoData && |
362 | 0 | !((std::isnan(singleNDV) && std::isnan(dfNoData)) || |
363 | 0 | (singleNDV == dfNoData))) |
364 | 0 | { |
365 | 0 | return false; |
366 | 0 | } |
367 | 0 | } |
368 | 0 | return true; |
369 | 0 | } |
370 | | |
371 | | /************************************************************************/ |
372 | | /* GDALComputedDataset::AddSources() */ |
373 | | /************************************************************************/ |
374 | | |
375 | | void GDALComputedDataset::AddSources(GDALComputedRasterBand *poBand) |
376 | 0 | { |
377 | 0 | auto poSourcedRasterBand = |
378 | 0 | cpl::down_cast<VRTSourcedRasterBand *>(m_oVRTDS.GetRasterBand(1)); |
379 | |
|
380 | 0 | bool hasAtLeastOneNDV = false; |
381 | 0 | double singleNDV = 0; |
382 | 0 | const bool bSameNDV = HaveAllBandsSameNoDataValue( |
383 | 0 | m_poBands.data(), m_poBands.size(), hasAtLeastOneNDV, singleNDV); |
384 | | |
385 | | // For inputs that are instances of GDALComputedDataset, clone them |
386 | | // to make sure we do not depend on temporary instances, |
387 | | // such as "a + b + c", which is evaluated as "(a + b) + c", and the |
388 | | // temporary band/dataset corresponding to a + b will go out of scope |
389 | | // quickly. |
390 | 0 | for (GDALRasterBand *&band : m_poBands) |
391 | 0 | { |
392 | 0 | auto poDS = band->GetDataset(); |
393 | 0 | if (auto poComputedDS = dynamic_cast<GDALComputedDataset *>(poDS)) |
394 | 0 | { |
395 | 0 | auto poComputedDSNew = |
396 | 0 | std::make_unique<GDALComputedDataset>(*poComputedDS); |
397 | 0 | band = poComputedDSNew->GetRasterBand(1); |
398 | 0 | m_bandDS.emplace_back(poComputedDSNew.release()); |
399 | 0 | } |
400 | |
|
401 | 0 | int bHasNoData = false; |
402 | 0 | const double dfNoData = band->GetNoDataValue(&bHasNoData); |
403 | 0 | if (bHasNoData) |
404 | 0 | { |
405 | 0 | poSourcedRasterBand->AddComplexSource(band, -1, -1, -1, -1, -1, -1, |
406 | 0 | -1, -1, 0, 1, dfNoData); |
407 | 0 | } |
408 | 0 | else |
409 | 0 | { |
410 | 0 | poSourcedRasterBand->AddSimpleSource(band); |
411 | 0 | } |
412 | 0 | poSourcedRasterBand->papoSources[poSourcedRasterBand->nSources - 1] |
413 | 0 | ->SetName(CPLSPrintf("source%d", poSourcedRasterBand->nSources)); |
414 | 0 | } |
415 | 0 | if (hasAtLeastOneNDV) |
416 | 0 | { |
417 | 0 | poBand->m_bHasNoData = true; |
418 | 0 | if (bSameNDV) |
419 | 0 | { |
420 | 0 | poBand->m_dfNoDataValue = singleNDV; |
421 | 0 | } |
422 | 0 | else |
423 | 0 | { |
424 | 0 | poBand->m_dfNoDataValue = std::numeric_limits<double>::quiet_NaN(); |
425 | 0 | } |
426 | 0 | poSourcedRasterBand->SetNoDataValue(poBand->m_dfNoDataValue); |
427 | 0 | } |
428 | 0 | } |
429 | | |
430 | | /************************************************************************/ |
431 | | /* OperationToFunctionName() */ |
432 | | /************************************************************************/ |
433 | | |
434 | | /* static */ const char *GDALComputedDataset::OperationToFunctionName( |
435 | | GDALComputedRasterBand::Operation op) |
436 | 0 | { |
437 | 0 | const char *ret = ""; |
438 | 0 | switch (op) |
439 | 0 | { |
440 | 0 | case GDALComputedRasterBand::Operation::OP_ADD: |
441 | 0 | ret = "sum"; |
442 | 0 | break; |
443 | 0 | case GDALComputedRasterBand::Operation::OP_SUBTRACT: |
444 | 0 | ret = "diff"; |
445 | 0 | break; |
446 | 0 | case GDALComputedRasterBand::Operation::OP_MULTIPLY: |
447 | 0 | ret = "mul"; |
448 | 0 | break; |
449 | 0 | case GDALComputedRasterBand::Operation::OP_DIVIDE: |
450 | 0 | ret = "div"; |
451 | 0 | break; |
452 | 0 | case GDALComputedRasterBand::Operation::OP_MIN: |
453 | 0 | ret = "min"; |
454 | 0 | break; |
455 | 0 | case GDALComputedRasterBand::Operation::OP_MAX: |
456 | 0 | ret = "max"; |
457 | 0 | break; |
458 | 0 | case GDALComputedRasterBand::Operation::OP_MEAN: |
459 | 0 | ret = "mean"; |
460 | 0 | break; |
461 | 0 | case GDALComputedRasterBand::Operation::OP_GT: |
462 | 0 | ret = ">"; |
463 | 0 | break; |
464 | 0 | case GDALComputedRasterBand::Operation::OP_GE: |
465 | 0 | ret = ">="; |
466 | 0 | break; |
467 | 0 | case GDALComputedRasterBand::Operation::OP_LT: |
468 | 0 | ret = "<"; |
469 | 0 | break; |
470 | 0 | case GDALComputedRasterBand::Operation::OP_LE: |
471 | 0 | ret = "<="; |
472 | 0 | break; |
473 | 0 | case GDALComputedRasterBand::Operation::OP_EQ: |
474 | 0 | ret = "=="; |
475 | 0 | break; |
476 | 0 | case GDALComputedRasterBand::Operation::OP_NE: |
477 | 0 | ret = "!="; |
478 | 0 | break; |
479 | 0 | case GDALComputedRasterBand::Operation::OP_LOGICAL_AND: |
480 | 0 | ret = "&&"; |
481 | 0 | break; |
482 | 0 | case GDALComputedRasterBand::Operation::OP_LOGICAL_OR: |
483 | 0 | ret = "||"; |
484 | 0 | break; |
485 | 0 | case GDALComputedRasterBand::Operation::OP_CAST: |
486 | 0 | case GDALComputedRasterBand::Operation::OP_TERNARY: |
487 | 0 | break; |
488 | 0 | } |
489 | 0 | return ret; |
490 | 0 | } |
491 | | |
492 | | /************************************************************************/ |
493 | | /* GDALComputedRasterBand() */ |
494 | | /************************************************************************/ |
495 | | |
496 | | GDALComputedRasterBand::GDALComputedRasterBand( |
497 | | const GDALComputedRasterBand &other, bool) |
498 | 0 | : GDALRasterBand() |
499 | 0 | { |
500 | 0 | nRasterXSize = other.nRasterXSize; |
501 | 0 | nRasterYSize = other.nRasterYSize; |
502 | 0 | eDataType = other.eDataType; |
503 | 0 | nBlockXSize = other.nBlockXSize; |
504 | 0 | nBlockYSize = other.nBlockYSize; |
505 | 0 | } |
506 | | |
507 | | //! @cond Doxygen_Suppress |
508 | | |
509 | | /************************************************************************/ |
510 | | /* GDALComputedRasterBand() */ |
511 | | /************************************************************************/ |
512 | | |
513 | | GDALComputedRasterBand::GDALComputedRasterBand( |
514 | | Operation op, const std::vector<const GDALRasterBand *> &bands, |
515 | | double constant) |
516 | 0 | { |
517 | 0 | CPLAssert(op == Operation::OP_ADD || op == Operation::OP_MIN || |
518 | 0 | op == Operation::OP_MAX || op == Operation::OP_MEAN || |
519 | 0 | op == Operation::OP_TERNARY); |
520 | | |
521 | 0 | CPLAssert(!bands.empty()); |
522 | 0 | nRasterXSize = bands[0]->GetXSize(); |
523 | 0 | nRasterYSize = bands[0]->GetYSize(); |
524 | 0 | eDataType = bands[0]->GetRasterDataType(); |
525 | 0 | for (size_t i = 1; i < bands.size(); ++i) |
526 | 0 | { |
527 | 0 | eDataType = GDALDataTypeUnion(eDataType, bands[i]->GetRasterDataType()); |
528 | 0 | } |
529 | |
|
530 | 0 | bool hasAtLeastOneNDV = false; |
531 | 0 | double singleNDV = 0; |
532 | 0 | const bool bSameNDV = |
533 | 0 | HaveAllBandsSameNoDataValue(const_cast<GDALRasterBand **>(bands.data()), |
534 | 0 | bands.size(), hasAtLeastOneNDV, singleNDV); |
535 | |
|
536 | 0 | if (!bSameNDV) |
537 | 0 | { |
538 | 0 | eDataType = eDataType == GDT_Float64 ? GDT_Float64 : GDT_Float32; |
539 | 0 | } |
540 | 0 | else if (op == Operation::OP_TERNARY) |
541 | 0 | { |
542 | 0 | CPLAssert(bands.size() == 3); |
543 | 0 | eDataType = GDALDataTypeUnion(bands[1]->GetRasterDataType(), |
544 | 0 | bands[2]->GetRasterDataType()); |
545 | 0 | } |
546 | 0 | else if (!std::isnan(constant) && eDataType != GDT_Float64) |
547 | 0 | { |
548 | 0 | if (op == Operation::OP_MIN || op == Operation::OP_MAX) |
549 | 0 | { |
550 | 0 | eDataType = GDALDataTypeUnionWithValue(eDataType, constant, false); |
551 | 0 | } |
552 | 0 | else |
553 | 0 | { |
554 | 0 | eDataType = (static_cast<float>(constant) == constant) |
555 | 0 | ? GDT_Float32 |
556 | 0 | : GDT_Float64; |
557 | 0 | } |
558 | 0 | } |
559 | 0 | bands[0]->GetBlockSize(&nBlockXSize, &nBlockYSize); |
560 | 0 | auto l_poDS = std::make_unique<GDALComputedDataset>( |
561 | 0 | this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize, |
562 | 0 | op, bands, constant); |
563 | 0 | m_poOwningDS.reset(l_poDS.release()); |
564 | 0 | } |
565 | | |
566 | | /************************************************************************/ |
567 | | /* GDALComputedRasterBand() */ |
568 | | /************************************************************************/ |
569 | | |
570 | | GDALComputedRasterBand::GDALComputedRasterBand(Operation op, |
571 | | const GDALRasterBand &firstBand, |
572 | | const GDALRasterBand &secondBand) |
573 | 0 | { |
574 | 0 | nRasterXSize = firstBand.GetXSize(); |
575 | 0 | nRasterYSize = firstBand.GetYSize(); |
576 | |
|
577 | 0 | bool hasAtLeastOneNDV = false; |
578 | 0 | double singleNDV = 0; |
579 | 0 | GDALRasterBand *apoBands[] = {const_cast<GDALRasterBand *>(&firstBand), |
580 | 0 | const_cast<GDALRasterBand *>(&secondBand)}; |
581 | 0 | const bool bSameNDV = |
582 | 0 | HaveAllBandsSameNoDataValue(apoBands, 2, hasAtLeastOneNDV, singleNDV); |
583 | |
|
584 | 0 | const auto firstDT = firstBand.GetRasterDataType(); |
585 | 0 | const auto secondDT = secondBand.GetRasterDataType(); |
586 | 0 | if (!bSameNDV) |
587 | 0 | eDataType = (firstDT == GDT_Float64 || secondDT == GDT_Float64) |
588 | 0 | ? GDT_Float64 |
589 | 0 | : GDT_Float32; |
590 | 0 | else if (IsComparisonOperator(op)) |
591 | 0 | eDataType = GDT_Byte; |
592 | 0 | else if (op == Operation::OP_ADD && firstDT == GDT_Byte && |
593 | 0 | secondDT == GDT_Byte) |
594 | 0 | eDataType = GDT_UInt16; |
595 | 0 | else if (firstDT == GDT_Float32 && secondDT == GDT_Float32) |
596 | 0 | eDataType = GDT_Float32; |
597 | 0 | else if ((op == Operation::OP_MIN || op == Operation::OP_MAX) && |
598 | 0 | firstDT == secondDT) |
599 | 0 | eDataType = firstDT; |
600 | 0 | else |
601 | 0 | eDataType = GDT_Float64; |
602 | 0 | firstBand.GetBlockSize(&nBlockXSize, &nBlockYSize); |
603 | 0 | auto l_poDS = std::make_unique<GDALComputedDataset>( |
604 | 0 | this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize, |
605 | 0 | op, &firstBand, nullptr, &secondBand, nullptr); |
606 | 0 | m_poOwningDS.reset(l_poDS.release()); |
607 | 0 | } |
608 | | |
609 | | /************************************************************************/ |
610 | | /* GDALComputedRasterBand() */ |
611 | | /************************************************************************/ |
612 | | |
613 | | GDALComputedRasterBand::GDALComputedRasterBand(Operation op, double constant, |
614 | | const GDALRasterBand &band) |
615 | 0 | { |
616 | 0 | CPLAssert(op == Operation::OP_DIVIDE || IsComparisonOperator(op)); |
617 | | |
618 | 0 | nRasterXSize = band.GetXSize(); |
619 | 0 | nRasterYSize = band.GetYSize(); |
620 | 0 | const auto firstDT = band.GetRasterDataType(); |
621 | 0 | if (IsComparisonOperator(op)) |
622 | 0 | eDataType = GDT_Byte; |
623 | 0 | else if (firstDT == GDT_Float32 && static_cast<float>(constant) == constant) |
624 | 0 | eDataType = GDT_Float32; |
625 | 0 | else |
626 | 0 | eDataType = GDT_Float64; |
627 | 0 | band.GetBlockSize(&nBlockXSize, &nBlockYSize); |
628 | 0 | auto l_poDS = std::make_unique<GDALComputedDataset>( |
629 | 0 | this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize, |
630 | 0 | op, nullptr, &constant, &band, nullptr); |
631 | 0 | m_poOwningDS.reset(l_poDS.release()); |
632 | 0 | } |
633 | | |
634 | | /************************************************************************/ |
635 | | /* GDALComputedRasterBand() */ |
636 | | /************************************************************************/ |
637 | | |
638 | | GDALComputedRasterBand::GDALComputedRasterBand(Operation op, |
639 | | const GDALRasterBand &band, |
640 | | double constant) |
641 | 0 | { |
642 | 0 | nRasterXSize = band.GetXSize(); |
643 | 0 | nRasterYSize = band.GetYSize(); |
644 | 0 | const auto firstDT = band.GetRasterDataType(); |
645 | 0 | if (IsComparisonOperator(op)) |
646 | 0 | eDataType = GDT_Byte; |
647 | 0 | else if (op == Operation::OP_ADD && firstDT == GDT_Byte && |
648 | 0 | constant >= -128 && constant <= 127 && |
649 | 0 | std::floor(constant) == constant) |
650 | 0 | eDataType = GDT_Byte; |
651 | 0 | else if (firstDT == GDT_Float32 && static_cast<float>(constant) == constant) |
652 | 0 | eDataType = GDT_Float32; |
653 | 0 | else |
654 | 0 | eDataType = GDT_Float64; |
655 | 0 | band.GetBlockSize(&nBlockXSize, &nBlockYSize); |
656 | 0 | auto l_poDS = std::make_unique<GDALComputedDataset>( |
657 | 0 | this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize, |
658 | 0 | op, &band, nullptr, nullptr, &constant); |
659 | 0 | m_poOwningDS.reset(l_poDS.release()); |
660 | 0 | } |
661 | | |
662 | | /************************************************************************/ |
663 | | /* GDALComputedRasterBand() */ |
664 | | /************************************************************************/ |
665 | | |
666 | | GDALComputedRasterBand::GDALComputedRasterBand(Operation op, |
667 | | const GDALRasterBand &band, |
668 | | GDALDataType dt) |
669 | 0 | { |
670 | 0 | CPLAssert(op == Operation::OP_CAST); |
671 | 0 | nRasterXSize = band.GetXSize(); |
672 | 0 | nRasterYSize = band.GetYSize(); |
673 | 0 | eDataType = dt; |
674 | 0 | band.GetBlockSize(&nBlockXSize, &nBlockYSize); |
675 | 0 | auto l_poDS = std::make_unique<GDALComputedDataset>( |
676 | 0 | this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize, |
677 | 0 | op, &band, nullptr, nullptr, nullptr); |
678 | 0 | m_poOwningDS.reset(l_poDS.release()); |
679 | 0 | } |
680 | | |
681 | | //! @endcond |
682 | | |
683 | | /************************************************************************/ |
684 | | /* ~GDALComputedRasterBand() */ |
685 | | /************************************************************************/ |
686 | | |
687 | | GDALComputedRasterBand::~GDALComputedRasterBand() |
688 | 0 | { |
689 | 0 | if (m_poOwningDS) |
690 | 0 | cpl::down_cast<GDALComputedDataset *>(m_poOwningDS.get())->nBands = 0; |
691 | 0 | poDS = nullptr; |
692 | 0 | } |
693 | | |
694 | | /************************************************************************/ |
695 | | /* GDALComputedRasterBand::GetNoDataValue() */ |
696 | | /************************************************************************/ |
697 | | |
698 | | double GDALComputedRasterBand::GetNoDataValue(int *pbHasNoData) |
699 | 0 | { |
700 | 0 | if (pbHasNoData) |
701 | 0 | *pbHasNoData = m_bHasNoData; |
702 | 0 | return m_dfNoDataValue; |
703 | 0 | } |
704 | | |
705 | | /************************************************************************/ |
706 | | /* GDALComputedRasterBandRelease() */ |
707 | | /************************************************************************/ |
708 | | |
709 | | /** Release a GDALComputedRasterBandH |
710 | | * |
711 | | * @since 3.12 |
712 | | */ |
713 | | void GDALComputedRasterBandRelease(GDALComputedRasterBandH hBand) |
714 | 0 | { |
715 | 0 | delete GDALComputedRasterBand::FromHandle(hBand); |
716 | 0 | } |
717 | | |
718 | | /************************************************************************/ |
719 | | /* IReadBlock() */ |
720 | | /************************************************************************/ |
721 | | |
722 | | CPLErr GDALComputedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, |
723 | | void *pData) |
724 | 0 | { |
725 | 0 | auto l_poDS = cpl::down_cast<GDALComputedDataset *>(poDS); |
726 | 0 | return l_poDS->m_oVRTDS.GetRasterBand(1)->ReadBlock(nBlockXOff, nBlockYOff, |
727 | 0 | pData); |
728 | 0 | } |
729 | | |
730 | | /************************************************************************/ |
731 | | /* IRasterIO() */ |
732 | | /************************************************************************/ |
733 | | |
734 | | CPLErr GDALComputedRasterBand::IRasterIO( |
735 | | GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize, |
736 | | void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType, |
737 | | GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg) |
738 | 0 | { |
739 | 0 | auto l_poDS = cpl::down_cast<GDALComputedDataset *>(poDS); |
740 | 0 | return l_poDS->m_oVRTDS.GetRasterBand(1)->RasterIO( |
741 | 0 | eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize, |
742 | 0 | eBufType, nPixelSpace, nLineSpace, psExtraArg); |
743 | 0 | } |