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