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 1067 : poSourcedRasterBand->papoSources[poSourcedRasterBand->nSources - 1]
448 1067 : ->SetName(CPLSPrintf("source%d", poSourcedRasterBand->nSources));
449 : }
450 794 : if (hasAtLeastOneNDV)
451 : {
452 5 : poBand->m_bHasNoData = true;
453 5 : if (bSameNDV)
454 : {
455 2 : poBand->m_dfNoDataValue = singleNDV;
456 : }
457 : else
458 : {
459 3 : poBand->m_dfNoDataValue = std::numeric_limits<double>::quiet_NaN();
460 : }
461 5 : poSourcedRasterBand->SetNoDataValue(poBand->m_dfNoDataValue);
462 : }
463 794 : }
464 :
465 : /************************************************************************/
466 : /* OperationToFunctionName() */
467 : /************************************************************************/
468 :
469 277 : /* static */ const char *GDALComputedDataset::OperationToFunctionName(
470 : GDALComputedRasterBand::Operation op)
471 : {
472 277 : const char *ret = "";
473 277 : switch (op)
474 : {
475 47 : case GDALComputedRasterBand::Operation::OP_ADD:
476 47 : ret = "sum";
477 47 : break;
478 3 : case GDALComputedRasterBand::Operation::OP_SUBTRACT:
479 3 : ret = "diff";
480 3 : break;
481 41 : case GDALComputedRasterBand::Operation::OP_MULTIPLY:
482 41 : ret = "mul";
483 41 : break;
484 0 : case GDALComputedRasterBand::Operation::OP_DIVIDE:
485 0 : ret = "div";
486 0 : break;
487 9 : case GDALComputedRasterBand::Operation::OP_MIN:
488 9 : ret = "min";
489 9 : break;
490 9 : case GDALComputedRasterBand::Operation::OP_MAX:
491 9 : ret = "max";
492 9 : break;
493 2 : case GDALComputedRasterBand::Operation::OP_MEAN:
494 2 : ret = "mean";
495 2 : break;
496 13 : case GDALComputedRasterBand::Operation::OP_GT:
497 13 : ret = ">";
498 13 : break;
499 16 : case GDALComputedRasterBand::Operation::OP_GE:
500 16 : ret = ">=";
501 16 : break;
502 13 : case GDALComputedRasterBand::Operation::OP_LT:
503 13 : ret = "<";
504 13 : break;
505 14 : case GDALComputedRasterBand::Operation::OP_LE:
506 14 : ret = "<=";
507 14 : break;
508 18 : case GDALComputedRasterBand::Operation::OP_EQ:
509 18 : ret = "==";
510 18 : break;
511 20 : case GDALComputedRasterBand::Operation::OP_NE:
512 20 : ret = "!=";
513 20 : break;
514 18 : case GDALComputedRasterBand::Operation::OP_LOGICAL_AND:
515 18 : ret = "&&";
516 18 : break;
517 23 : case GDALComputedRasterBand::Operation::OP_LOGICAL_OR:
518 23 : ret = "||";
519 23 : break;
520 24 : case GDALComputedRasterBand::Operation::OP_CAST:
521 : case GDALComputedRasterBand::Operation::OP_TERNARY:
522 24 : break;
523 3 : case GDALComputedRasterBand::Operation::OP_ABS:
524 3 : ret = "mod";
525 3 : break;
526 2 : case GDALComputedRasterBand::Operation::OP_SQRT:
527 2 : ret = "sqrt";
528 2 : break;
529 0 : case GDALComputedRasterBand::Operation::OP_LOG:
530 0 : ret = "log";
531 0 : break;
532 2 : case GDALComputedRasterBand::Operation::OP_LOG10:
533 2 : ret = "log10";
534 2 : break;
535 0 : case GDALComputedRasterBand::Operation::OP_POW:
536 0 : ret = "pow";
537 0 : break;
538 : }
539 277 : return ret;
540 : }
541 :
542 : /************************************************************************/
543 : /* GDALComputedRasterBand() */
544 : /************************************************************************/
545 :
546 483 : GDALComputedRasterBand::GDALComputedRasterBand(
547 483 : const GDALComputedRasterBand &other, bool)
548 483 : : GDALRasterBand()
549 : {
550 483 : nRasterXSize = other.nRasterXSize;
551 483 : nRasterYSize = other.nRasterYSize;
552 483 : eDataType = other.eDataType;
553 483 : nBlockXSize = other.nBlockXSize;
554 483 : nBlockYSize = other.nBlockYSize;
555 483 : }
556 :
557 : //! @cond Doxygen_Suppress
558 :
559 : /************************************************************************/
560 : /* GDALComputedRasterBand() */
561 : /************************************************************************/
562 :
563 37 : GDALComputedRasterBand::GDALComputedRasterBand(
564 : Operation op, const std::vector<const GDALRasterBand *> &bands,
565 37 : double constant)
566 : {
567 37 : CPLAssert(op == Operation::OP_ADD || op == Operation::OP_MIN ||
568 : op == Operation::OP_MAX || op == Operation::OP_MEAN ||
569 : op == Operation::OP_TERNARY);
570 :
571 37 : CPLAssert(!bands.empty());
572 37 : nRasterXSize = bands[0]->GetXSize();
573 37 : nRasterYSize = bands[0]->GetYSize();
574 37 : eDataType = bands[0]->GetRasterDataType();
575 95 : for (size_t i = 1; i < bands.size(); ++i)
576 : {
577 58 : eDataType = GDALDataTypeUnion(eDataType, bands[i]->GetRasterDataType());
578 : }
579 :
580 37 : bool hasAtLeastOneNDV = false;
581 37 : double singleNDV = 0;
582 : const bool bSameNDV =
583 37 : HaveAllBandsSameNoDataValue(const_cast<GDALRasterBand **>(bands.data()),
584 : bands.size(), hasAtLeastOneNDV, singleNDV);
585 :
586 37 : if (!bSameNDV)
587 : {
588 0 : eDataType = eDataType == GDT_Float64 ? GDT_Float64 : GDT_Float32;
589 : }
590 37 : else if (op == Operation::OP_TERNARY)
591 : {
592 18 : CPLAssert(bands.size() == 3);
593 18 : eDataType = GDALDataTypeUnion(bands[1]->GetRasterDataType(),
594 18 : bands[2]->GetRasterDataType());
595 : }
596 19 : else if (!std::isnan(constant) && eDataType != GDT_Float64)
597 : {
598 4 : if (op == Operation::OP_MIN || op == Operation::OP_MAX)
599 : {
600 4 : eDataType = GDALDataTypeUnionWithValue(eDataType, constant, false);
601 : }
602 : else
603 : {
604 0 : eDataType = (static_cast<float>(constant) == constant)
605 0 : ? GDT_Float32
606 : : GDT_Float64;
607 : }
608 : }
609 37 : bands[0]->GetBlockSize(&nBlockXSize, &nBlockYSize);
610 : auto l_poDS = std::make_unique<GDALComputedDataset>(
611 37 : this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
612 74 : op, bands, constant);
613 37 : m_poOwningDS.reset(l_poDS.release());
614 37 : }
615 :
616 : /************************************************************************/
617 : /* GDALComputedRasterBand() */
618 : /************************************************************************/
619 :
620 78 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
621 : const GDALRasterBand &firstBand,
622 78 : const GDALRasterBand &secondBand)
623 : {
624 78 : nRasterXSize = firstBand.GetXSize();
625 78 : nRasterYSize = firstBand.GetYSize();
626 :
627 78 : bool hasAtLeastOneNDV = false;
628 78 : double singleNDV = 0;
629 : GDALRasterBand *apoBands[] = {const_cast<GDALRasterBand *>(&firstBand),
630 78 : const_cast<GDALRasterBand *>(&secondBand)};
631 : const bool bSameNDV =
632 78 : HaveAllBandsSameNoDataValue(apoBands, 2, hasAtLeastOneNDV, singleNDV);
633 :
634 78 : const auto firstDT = firstBand.GetRasterDataType();
635 78 : const auto secondDT = secondBand.GetRasterDataType();
636 78 : if (!bSameNDV)
637 6 : eDataType = (firstDT == GDT_Float64 || secondDT == GDT_Float64)
638 6 : ? GDT_Float64
639 : : GDT_Float32;
640 75 : else if (IsComparisonOperator(op))
641 43 : eDataType = GDT_Byte;
642 32 : else if (op == Operation::OP_ADD && firstDT == GDT_Byte &&
643 : secondDT == GDT_Byte)
644 2 : eDataType = GDT_UInt16;
645 30 : else if (firstDT == GDT_Float32 && secondDT == GDT_Float32)
646 1 : eDataType = GDT_Float32;
647 29 : else if ((op == Operation::OP_MIN || op == Operation::OP_MAX) &&
648 : firstDT == secondDT)
649 1 : eDataType = firstDT;
650 : else
651 28 : eDataType = GDT_Float64;
652 78 : firstBand.GetBlockSize(&nBlockXSize, &nBlockYSize);
653 : auto l_poDS = std::make_unique<GDALComputedDataset>(
654 78 : this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
655 156 : op, &firstBand, nullptr, &secondBand, nullptr);
656 78 : m_poOwningDS.reset(l_poDS.release());
657 78 : }
658 :
659 : /************************************************************************/
660 : /* GDALComputedRasterBand() */
661 : /************************************************************************/
662 :
663 22 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op, double constant,
664 22 : const GDALRasterBand &band)
665 : {
666 22 : CPLAssert(op == Operation::OP_DIVIDE || IsComparisonOperator(op) ||
667 : op == Operation::OP_POW);
668 :
669 22 : nRasterXSize = band.GetXSize();
670 22 : nRasterYSize = band.GetYSize();
671 22 : const auto firstDT = band.GetRasterDataType();
672 22 : if (IsComparisonOperator(op))
673 18 : eDataType = GDT_Byte;
674 4 : else if (firstDT == GDT_Float32 && static_cast<float>(constant) == constant)
675 0 : eDataType = GDT_Float32;
676 : else
677 4 : eDataType = GDT_Float64;
678 22 : band.GetBlockSize(&nBlockXSize, &nBlockYSize);
679 : auto l_poDS = std::make_unique<GDALComputedDataset>(
680 22 : this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
681 44 : op, nullptr, &constant, &band, nullptr);
682 22 : m_poOwningDS.reset(l_poDS.release());
683 22 : }
684 :
685 : /************************************************************************/
686 : /* GDALComputedRasterBand() */
687 : /************************************************************************/
688 :
689 141 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
690 : const GDALRasterBand &band,
691 141 : double constant)
692 : {
693 141 : nRasterXSize = band.GetXSize();
694 141 : nRasterYSize = band.GetYSize();
695 141 : const auto firstDT = band.GetRasterDataType();
696 141 : if (IsComparisonOperator(op))
697 74 : eDataType = GDT_Byte;
698 67 : else if (op == Operation::OP_ADD && firstDT == GDT_Byte &&
699 10 : constant >= -128 && constant <= 127 &&
700 10 : std::floor(constant) == constant)
701 10 : eDataType = GDT_Byte;
702 57 : else if (firstDT == GDT_Float32 && static_cast<float>(constant) == constant)
703 8 : eDataType = GDT_Float32;
704 : else
705 49 : eDataType = GDT_Float64;
706 141 : band.GetBlockSize(&nBlockXSize, &nBlockYSize);
707 : auto l_poDS = std::make_unique<GDALComputedDataset>(
708 141 : this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
709 282 : op, &band, nullptr, nullptr, &constant);
710 141 : m_poOwningDS.reset(l_poDS.release());
711 141 : }
712 :
713 : /************************************************************************/
714 : /* GDALComputedRasterBand() */
715 : /************************************************************************/
716 :
717 9 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
718 9 : const GDALRasterBand &band)
719 : {
720 9 : CPLAssert(op == Operation::OP_ABS || op == Operation::OP_SQRT ||
721 : op == Operation::OP_LOG || op == Operation::OP_LOG10);
722 9 : nRasterXSize = band.GetXSize();
723 9 : nRasterYSize = band.GetYSize();
724 9 : eDataType =
725 9 : band.GetRasterDataType() == GDT_Float32 ? GDT_Float32 : GDT_Float64;
726 9 : band.GetBlockSize(&nBlockXSize, &nBlockYSize);
727 : auto l_poDS = std::make_unique<GDALComputedDataset>(
728 9 : this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
729 18 : op, &band, nullptr, nullptr, nullptr);
730 9 : m_poOwningDS.reset(l_poDS.release());
731 9 : }
732 :
733 : /************************************************************************/
734 : /* GDALComputedRasterBand() */
735 : /************************************************************************/
736 :
737 24 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
738 : const GDALRasterBand &band,
739 24 : GDALDataType dt)
740 : {
741 24 : CPLAssert(op == Operation::OP_CAST);
742 24 : nRasterXSize = band.GetXSize();
743 24 : nRasterYSize = band.GetYSize();
744 24 : eDataType = dt;
745 24 : band.GetBlockSize(&nBlockXSize, &nBlockYSize);
746 : auto l_poDS = std::make_unique<GDALComputedDataset>(
747 24 : this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
748 48 : op, &band, nullptr, nullptr, nullptr);
749 24 : m_poOwningDS.reset(l_poDS.release());
750 24 : }
751 :
752 : //! @endcond
753 :
754 : /************************************************************************/
755 : /* ~GDALComputedRasterBand() */
756 : /************************************************************************/
757 :
758 1440 : GDALComputedRasterBand::~GDALComputedRasterBand()
759 : {
760 794 : if (m_poOwningDS)
761 311 : cpl::down_cast<GDALComputedDataset *>(m_poOwningDS.get())->nBands = 0;
762 794 : poDS = nullptr;
763 1440 : }
764 :
765 : /************************************************************************/
766 : /* GDALComputedRasterBand::GetNoDataValue() */
767 : /************************************************************************/
768 :
769 1377 : double GDALComputedRasterBand::GetNoDataValue(int *pbHasNoData)
770 : {
771 1377 : if (pbHasNoData)
772 1377 : *pbHasNoData = m_bHasNoData;
773 1377 : return m_dfNoDataValue;
774 : }
775 :
776 : /************************************************************************/
777 : /* GDALComputedRasterBandRelease() */
778 : /************************************************************************/
779 :
780 : /** Release a GDALComputedRasterBandH
781 : *
782 : * @since 3.12
783 : */
784 163 : void GDALComputedRasterBandRelease(GDALComputedRasterBandH hBand)
785 : {
786 163 : delete GDALComputedRasterBand::FromHandle(hBand);
787 163 : }
788 :
789 : /************************************************************************/
790 : /* IReadBlock() */
791 : /************************************************************************/
792 :
793 239 : CPLErr GDALComputedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
794 : void *pData)
795 : {
796 239 : auto l_poDS = cpl::down_cast<GDALComputedDataset *>(poDS);
797 239 : return l_poDS->m_oVRTDS.GetRasterBand(1)->ReadBlock(nBlockXOff, nBlockYOff,
798 239 : pData);
799 : }
800 :
801 : /************************************************************************/
802 : /* IRasterIO() */
803 : /************************************************************************/
804 :
805 208 : CPLErr GDALComputedRasterBand::IRasterIO(
806 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
807 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
808 : GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg)
809 : {
810 208 : auto l_poDS = cpl::down_cast<GDALComputedDataset *>(poDS);
811 208 : return l_poDS->m_oVRTDS.GetRasterBand(1)->RasterIO(
812 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
813 208 : eBufType, nPixelSpace, nLineSpace, psExtraArg);
814 : }
|
