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 : #include <set>
19 : #include <utility>
20 :
21 : /************************************************************************/
22 : /* GDALComputedDataset */
23 : /************************************************************************/
24 :
25 658 : class GDALComputedDataset final : public GDALDataset
26 : {
27 : friend class GDALComputedRasterBand;
28 :
29 : std::string m_expr{};
30 : CPLStringList m_aosOptions{};
31 : std::vector<std::pair<std::string, GDALRasterBand *>> m_apoBands{};
32 : VRTDataset m_oVRTDS;
33 :
34 : std::pair<bool, std::string> AddSourceBand(const GDALRasterBand *band);
35 :
36 : void AddSources(GDALComputedRasterBand *poBand);
37 :
38 : static const char *
39 : OperationToFunctionName(GDALComputedRasterBand::Operation op,
40 : bool bForceMuparser = false);
41 :
42 : GDALComputedDataset &operator=(const GDALComputedDataset &) = delete;
43 : GDALComputedDataset(GDALComputedDataset &&) = delete;
44 : GDALComputedDataset &operator=(GDALComputedDataset &&) = delete;
45 :
46 : public:
47 : GDALComputedDataset(GDALComputedRasterBand *poBand, int nXSize, int nYSize,
48 : GDALDataType eDT, int nBlockXSize, int nBlockYSize,
49 : GDALComputedRasterBand::Operation op,
50 : const GDALRasterBand *firstBand, double *pFirstConstant,
51 : const GDALRasterBand *secondBand,
52 : double *pSecondConstant);
53 :
54 : GDALComputedDataset(GDALComputedRasterBand *poBand, int nXSize, int nYSize,
55 : GDALDataType eDT, int nBlockXSize, int nBlockYSize,
56 : GDALComputedRasterBand::Operation op,
57 : const std::vector<const GDALRasterBand *> &bands,
58 : double constant);
59 :
60 : ~GDALComputedDataset() override;
61 :
62 1 : CPLErr GetGeoTransform(GDALGeoTransform >) const override
63 : {
64 1 : return m_oVRTDS.GetGeoTransform(gt);
65 : }
66 :
67 1 : const OGRSpatialReference *GetSpatialRef() const override
68 : {
69 1 : return m_oVRTDS.GetSpatialRef();
70 : }
71 :
72 1 : CSLConstList GetMetadata(const char *pszDomain) override
73 : {
74 1 : return m_oVRTDS.GetMetadata(pszDomain);
75 : }
76 :
77 176 : const char *GetMetadataItem(const char *pszName,
78 : const char *pszDomain) override
79 : {
80 176 : return m_oVRTDS.GetMetadataItem(pszName, pszDomain);
81 : }
82 :
83 6 : void *GetInternalHandle(const char *pszHandleName) override
84 : {
85 6 : if (pszHandleName && EQUAL(pszHandleName, "VRT_DATASET"))
86 3 : return &m_oVRTDS;
87 3 : return nullptr;
88 : }
89 : };
90 :
91 : /************************************************************************/
92 : /* IsComparisonOperator() */
93 : /************************************************************************/
94 :
95 544 : static bool IsComparisonOperator(GDALComputedRasterBand::Operation op)
96 : {
97 544 : switch (op)
98 : {
99 288 : case GDALComputedRasterBand::Operation::OP_GT:
100 : case GDALComputedRasterBand::Operation::OP_GE:
101 : case GDALComputedRasterBand::Operation::OP_LT:
102 : case GDALComputedRasterBand::Operation::OP_LE:
103 : case GDALComputedRasterBand::Operation::OP_EQ:
104 : case GDALComputedRasterBand::Operation::OP_NE:
105 : case GDALComputedRasterBand::Operation::OP_LOGICAL_AND:
106 : case GDALComputedRasterBand::Operation::OP_LOGICAL_OR:
107 288 : return true;
108 256 : default:
109 256 : break;
110 : }
111 256 : return false;
112 : }
113 :
114 : /************************************************************************/
115 : /* AddSourceBand() */
116 : /************************************************************************/
117 :
118 : /** Returns a pair (bIsExprBand, osBandName) */
119 : std::pair<bool, std::string>
120 483 : GDALComputedDataset::AddSourceBand(const GDALRasterBand *band)
121 : {
122 483 : auto bandDS = band->GetDataset();
123 483 : if (auto poComputedDS = dynamic_cast<GDALComputedDataset *>(bandDS))
124 : {
125 151 : m_apoBands.insert(m_apoBands.end(), poComputedDS->m_apoBands.begin(),
126 302 : poComputedDS->m_apoBands.end());
127 151 : return {true, poComputedDS->m_expr};
128 : }
129 : else
130 : {
131 332 : std::string osName = CPLSPrintf("band_%p", band);
132 332 : m_apoBands.emplace_back(osName, const_cast<GDALRasterBand *>(band));
133 332 : return {false, osName};
134 : }
135 : }
136 :
137 : /************************************************************************/
138 : /* GDALComputedDataset() */
139 : /************************************************************************/
140 :
141 292 : GDALComputedDataset::GDALComputedDataset(
142 : GDALComputedRasterBand *poBand, int nXSize, int nYSize, GDALDataType eDT,
143 : int nBlockXSize, int nBlockYSize, GDALComputedRasterBand::Operation op,
144 : const GDALRasterBand *firstBand, double *pFirstConstant,
145 292 : const GDALRasterBand *secondBand, double *pSecondConstant)
146 292 : : m_oVRTDS(nXSize, nYSize, nBlockXSize, nBlockYSize)
147 : {
148 292 : CPLAssert(firstBand != nullptr || secondBand != nullptr);
149 584 : std::string osFirstBand;
150 292 : bool bCanUseBuiltin = true;
151 292 : if (firstBand)
152 : {
153 540 : const auto [bIsExprBand, name] = AddSourceBand(firstBand);
154 270 : bCanUseBuiltin = bCanUseBuiltin && !bIsExprBand;
155 270 : if (bIsExprBand)
156 89 : osFirstBand = "(";
157 270 : osFirstBand += name;
158 270 : if (bIsExprBand)
159 89 : osFirstBand += ')';
160 : }
161 584 : std::string osSecondBand;
162 292 : if (secondBand)
163 : {
164 236 : const auto [bIsExprBand, name] = AddSourceBand(secondBand);
165 118 : bCanUseBuiltin = bCanUseBuiltin && !bIsExprBand;
166 118 : if (bIsExprBand)
167 33 : osSecondBand = "(";
168 118 : osSecondBand += name;
169 118 : if (bIsExprBand)
170 33 : osSecondBand += ')';
171 : }
172 :
173 292 : nRasterXSize = nXSize;
174 292 : nRasterYSize = nYSize;
175 :
176 292 : if (auto poSrcDS = m_apoBands.front().second->GetDataset())
177 : {
178 292 : GDALGeoTransform gt;
179 292 : if (poSrcDS->GetGeoTransform(gt) == CE_None)
180 : {
181 140 : m_oVRTDS.SetGeoTransform(gt);
182 : }
183 :
184 292 : if (const auto *poSRS = poSrcDS->GetSpatialRef())
185 : {
186 140 : m_oVRTDS.SetSpatialRef(poSRS);
187 : }
188 : }
189 :
190 292 : if (op == GDALComputedRasterBand::Operation::OP_CAST)
191 : {
192 : #ifdef DEBUG
193 : // Just for code coverage...
194 24 : CPL_IGNORE_RET_VAL(GDALComputedDataset::OperationToFunctionName(op));
195 : #endif
196 :
197 24 : m_expr = osFirstBand;
198 24 : if (m_apoBands.size() > 1)
199 : {
200 5 : m_aosOptions.SetNameValue("subclass", "VRTDerivedRasterBand");
201 5 : m_aosOptions.SetNameValue("PixelFunctionType", "expression");
202 : m_aosOptions.SetNameValue("_PIXELFN_ARG_expression",
203 5 : m_expr.c_str());
204 : }
205 : else
206 : {
207 19 : m_aosOptions.SetNameValue("subclass", "VRTSourcedRasterBand");
208 : }
209 : }
210 : else
211 : {
212 268 : m_aosOptions.SetNameValue("subclass", "VRTDerivedRasterBand");
213 268 : if (IsComparisonOperator(op))
214 : {
215 135 : m_aosOptions.SetNameValue("PixelFunctionType", "expression");
216 135 : if (firstBand && secondBand)
217 : {
218 43 : m_expr = osFirstBand;
219 43 : m_expr += ' ';
220 43 : m_expr += GDALComputedDataset::OperationToFunctionName(op);
221 43 : m_expr += ' ';
222 43 : m_expr += osSecondBand;
223 : }
224 92 : else if (firstBand && pSecondConstant)
225 : {
226 74 : m_expr = osFirstBand;
227 74 : m_expr += ' ';
228 74 : m_expr += GDALComputedDataset::OperationToFunctionName(op);
229 74 : m_expr += ' ';
230 74 : m_expr += CPLSPrintf("%.17g", *pSecondConstant);
231 : }
232 18 : else if (pFirstConstant && secondBand)
233 : {
234 18 : m_expr = CPLSPrintf("%.17g", *pFirstConstant);
235 18 : m_expr += ' ';
236 18 : m_expr += GDALComputedDataset::OperationToFunctionName(op);
237 18 : m_expr += ' ';
238 18 : m_expr += osSecondBand;
239 : }
240 : else
241 : {
242 0 : CPLAssert(false);
243 : }
244 : m_aosOptions.SetNameValue("_PIXELFN_ARG_expression",
245 135 : m_expr.c_str());
246 : }
247 133 : else if (op == GDALComputedRasterBand::Operation::OP_SUBTRACT &&
248 : pSecondConstant)
249 : {
250 2 : m_aosOptions.SetNameValue("PixelFunctionType", "sum");
251 : m_aosOptions.SetNameValue("_PIXELFN_ARG_k",
252 2 : CPLSPrintf("%.17g", -(*pSecondConstant)));
253 2 : m_expr = osFirstBand;
254 2 : m_expr += CPLSPrintf(" - %.17g", *pSecondConstant);
255 : }
256 131 : else if (op == GDALComputedRasterBand::Operation::OP_DIVIDE)
257 : {
258 6 : if (pSecondConstant)
259 : {
260 1 : m_aosOptions.SetNameValue("PixelFunctionType", "mul");
261 : m_aosOptions.SetNameValue(
262 : "_PIXELFN_ARG_k",
263 1 : CPLSPrintf("%.17g", 1.0 / (*pSecondConstant)));
264 1 : m_expr = osFirstBand;
265 1 : m_expr += CPLSPrintf(" * %.17g", 1.0 / (*pSecondConstant));
266 : }
267 5 : else if (pFirstConstant)
268 : {
269 2 : m_aosOptions.SetNameValue("PixelFunctionType", "inv");
270 : m_aosOptions.SetNameValue("_PIXELFN_ARG_k",
271 2 : CPLSPrintf("%.17g", *pFirstConstant));
272 2 : m_expr = CPLSPrintf("%.17g / ", *pFirstConstant);
273 2 : m_expr += osSecondBand;
274 : }
275 : else
276 : {
277 3 : m_aosOptions.SetNameValue("PixelFunctionType", "div");
278 3 : m_expr = osFirstBand;
279 3 : m_expr += " / ";
280 3 : m_expr += osSecondBand;
281 : }
282 : }
283 125 : else if (op == GDALComputedRasterBand::Operation::OP_LOG)
284 : {
285 2 : CPLAssert(firstBand);
286 2 : CPLAssert(!secondBand);
287 2 : CPLAssert(!pFirstConstant);
288 2 : CPLAssert(!pSecondConstant);
289 2 : m_aosOptions.SetNameValue("PixelFunctionType", "expression");
290 2 : m_expr = "log(";
291 2 : m_expr += osFirstBand;
292 2 : m_expr += ')';
293 : m_aosOptions.SetNameValue("_PIXELFN_ARG_expression",
294 2 : m_expr.c_str());
295 : }
296 123 : else if (op == GDALComputedRasterBand::Operation::OP_POW)
297 : {
298 6 : if (firstBand && secondBand)
299 : {
300 2 : m_aosOptions.SetNameValue("PixelFunctionType", "expression");
301 2 : m_expr = osFirstBand;
302 2 : m_expr += " ^ ";
303 2 : m_expr += osSecondBand;
304 2 : m_aosOptions.SetNameValue("_PIXELFN_ARG_expression",
305 2 : m_expr.c_str());
306 : }
307 4 : else if (firstBand && pSecondConstant)
308 : {
309 2 : m_aosOptions.SetNameValue("PixelFunctionType", "pow");
310 : m_aosOptions.SetNameValue(
311 : "_PIXELFN_ARG_power",
312 2 : CPLSPrintf("%.17g", *pSecondConstant));
313 2 : m_expr = osFirstBand;
314 2 : m_expr += " ^ ";
315 2 : m_expr += CPLSPrintf("%.17g", *pSecondConstant);
316 : }
317 2 : else if (pFirstConstant && secondBand)
318 : {
319 2 : m_aosOptions.SetNameValue("PixelFunctionType", "exp");
320 : m_aosOptions.SetNameValue("_PIXELFN_ARG_base",
321 2 : CPLSPrintf("%.17g", *pFirstConstant));
322 2 : m_expr = CPLSPrintf("%.17g", *pFirstConstant);
323 2 : m_expr += " ^ ";
324 2 : m_expr += osSecondBand;
325 : }
326 : else
327 : {
328 0 : CPLAssert(false);
329 : }
330 : }
331 : else
332 : {
333 117 : if (firstBand && secondBand)
334 : {
335 48 : if (op == GDALComputedRasterBand::Operation::OP_MIN ||
336 47 : op == GDALComputedRasterBand::Operation::OP_MAX ||
337 : op == GDALComputedRasterBand::Operation::OP_MEAN)
338 : {
339 : m_expr +=
340 1 : GDALComputedDataset::OperationToFunctionName(op, true);
341 1 : m_expr += '(';
342 1 : m_expr += osFirstBand;
343 1 : m_expr += ',';
344 1 : m_expr += osSecondBand;
345 1 : m_expr += ')';
346 : }
347 : else
348 : {
349 47 : m_expr = osFirstBand;
350 47 : m_expr += ' ';
351 : m_expr +=
352 47 : GDALComputedDataset::OperationToFunctionName(op, true);
353 47 : m_expr += ' ';
354 47 : m_expr += osSecondBand;
355 : }
356 : }
357 69 : else if (firstBand && pSecondConstant)
358 : {
359 62 : m_expr = osFirstBand;
360 62 : m_expr += ' ';
361 : m_expr +=
362 62 : GDALComputedDataset::OperationToFunctionName(op, true);
363 62 : m_expr += ' ';
364 62 : m_expr += CPLSPrintf("%.17g", *pSecondConstant);
365 : }
366 7 : else if (pFirstConstant && secondBand)
367 : {
368 0 : m_expr = CPLSPrintf("%.17g", *pFirstConstant);
369 0 : m_expr += ' ';
370 : m_expr +=
371 0 : GDALComputedDataset::OperationToFunctionName(op, true);
372 0 : m_expr += ' ';
373 0 : m_expr += osSecondBand;
374 : }
375 : else
376 : {
377 7 : m_expr = GDALComputedDataset::OperationToFunctionName(op, true);
378 7 : m_expr += '(';
379 7 : m_expr += osFirstBand;
380 7 : m_expr += ')';
381 : }
382 :
383 117 : if (bCanUseBuiltin)
384 : {
385 : m_aosOptions.SetNameValue("PixelFunctionType",
386 51 : OperationToFunctionName(op));
387 51 : if (pSecondConstant)
388 : {
389 : m_aosOptions.SetNameValue(
390 : "_PIXELFN_ARG_k",
391 31 : CPLSPrintf("%.17g", *pSecondConstant));
392 : }
393 : }
394 : else
395 : {
396 66 : m_aosOptions.SetNameValue("PixelFunctionType", "expression");
397 : m_aosOptions.SetNameValue("_PIXELFN_ARG_expression",
398 66 : m_expr.c_str());
399 : }
400 : }
401 : }
402 292 : m_aosOptions.SetNameValue("_PIXELFN_ARG_propagateNoData", "true");
403 292 : m_oVRTDS.AddBand(eDT, m_aosOptions.List());
404 :
405 292 : SetBand(1, poBand);
406 :
407 292 : AddSources(poBand);
408 :
409 : // Otherwise coverity scan suggets adding annoying std::move()s
410 292 : CPL_IGNORE_RET_VAL(osFirstBand);
411 292 : CPL_IGNORE_RET_VAL(osSecondBand);
412 292 : }
413 :
414 : /************************************************************************/
415 : /* GDALComputedDataset() */
416 : /************************************************************************/
417 :
418 37 : GDALComputedDataset::GDALComputedDataset(
419 : GDALComputedRasterBand *poBand, int nXSize, int nYSize, GDALDataType eDT,
420 : int nBlockXSize, int nBlockYSize, GDALComputedRasterBand::Operation op,
421 37 : const std::vector<const GDALRasterBand *> &bands, double constant)
422 37 : : m_oVRTDS(nXSize, nYSize, nBlockXSize, nBlockYSize)
423 : {
424 37 : bool bCanUseBuiltin = true;
425 74 : std::vector<std::string> aosNames;
426 132 : for (const GDALRasterBand *poIterBand : bands)
427 : {
428 190 : const auto [bIsExprBand, name] = AddSourceBand(poIterBand);
429 95 : bCanUseBuiltin = bCanUseBuiltin && !bIsExprBand;
430 95 : if (!bIsExprBand)
431 66 : aosNames.push_back(std::move(name));
432 : else
433 29 : aosNames.push_back(std::string("(").append(name).append(")"));
434 : }
435 :
436 37 : nRasterXSize = nXSize;
437 37 : nRasterYSize = nYSize;
438 :
439 37 : if (auto poSrcDS = m_apoBands.front().second->GetDataset())
440 : {
441 37 : GDALGeoTransform gt;
442 37 : if (poSrcDS->GetGeoTransform(gt) == CE_None)
443 : {
444 16 : m_oVRTDS.SetGeoTransform(gt);
445 : }
446 :
447 37 : if (const auto *poSRS = poSrcDS->GetSpatialRef())
448 : {
449 16 : m_oVRTDS.SetSpatialRef(poSRS);
450 : }
451 : }
452 :
453 37 : m_aosOptions.SetNameValue("subclass", "VRTDerivedRasterBand");
454 37 : if (op == GDALComputedRasterBand::Operation::OP_TERNARY)
455 : {
456 18 : m_expr = aosNames[0];
457 18 : m_expr += " ? ";
458 18 : m_expr += aosNames[1];
459 18 : m_expr += " : ";
460 18 : m_expr += aosNames[2];
461 :
462 18 : m_aosOptions.SetNameValue("PixelFunctionType", "expression");
463 18 : m_aosOptions.SetNameValue("_PIXELFN_ARG_expression", m_expr.c_str());
464 : }
465 : else
466 : {
467 19 : m_expr = OperationToFunctionName(op);
468 19 : m_expr += '(';
469 19 : bool first = true;
470 60 : for (const auto &name : aosNames)
471 : {
472 41 : if (!first)
473 22 : m_expr += ", ";
474 41 : m_expr += name;
475 41 : first = false;
476 : }
477 19 : if (!std::isnan(constant))
478 : {
479 8 : if (!first)
480 8 : m_expr += ", ";
481 8 : m_expr += CPLSPrintf("%.17g", constant);
482 : }
483 19 : m_expr += ')';
484 :
485 19 : if (bCanUseBuiltin)
486 : {
487 : m_aosOptions.SetNameValue("PixelFunctionType",
488 14 : OperationToFunctionName(op));
489 14 : if (!std::isnan(constant))
490 : {
491 : m_aosOptions.SetNameValue("_PIXELFN_ARG_k",
492 4 : CPLSPrintf("%.17g", constant));
493 : }
494 14 : m_aosOptions.SetNameValue("_PIXELFN_ARG_propagateNoData", "true");
495 : }
496 : else
497 : {
498 5 : m_aosOptions.SetNameValue("PixelFunctionType", "expression");
499 : m_aosOptions.SetNameValue("_PIXELFN_ARG_expression",
500 5 : m_expr.c_str());
501 : }
502 : }
503 37 : m_oVRTDS.AddBand(eDT, m_aosOptions.List());
504 :
505 37 : SetBand(1, poBand);
506 :
507 37 : AddSources(poBand);
508 37 : }
509 :
510 : /************************************************************************/
511 : /* ~GDALComputedDataset() */
512 : /************************************************************************/
513 :
514 : GDALComputedDataset::~GDALComputedDataset() = default;
515 :
516 : /************************************************************************/
517 : /* HaveAllBandsSameNoDataValue() */
518 : /************************************************************************/
519 :
520 462 : static bool HaveAllBandsSameNoDataValue(GDALRasterBand **apoBands,
521 : size_t nBands, bool &hasAtLeastOneNDV,
522 : double &singleNDV)
523 : {
524 462 : hasAtLeastOneNDV = false;
525 462 : singleNDV = 0;
526 :
527 462 : int bFirstBandHasNoData = false;
528 1867 : for (size_t i = 0; i < nBands; ++i)
529 : {
530 1411 : int bHasNoData = false;
531 1411 : const double dfNoData = apoBands[i]->GetNoDataValue(&bHasNoData);
532 1411 : if (bHasNoData)
533 18 : hasAtLeastOneNDV = true;
534 1411 : if (i == 0)
535 : {
536 462 : bFirstBandHasNoData = bHasNoData;
537 462 : singleNDV = dfNoData;
538 : }
539 949 : else if (bHasNoData != bFirstBandHasNoData)
540 : {
541 6 : return false;
542 : }
543 955 : else if (bFirstBandHasNoData &&
544 8 : !((std::isnan(singleNDV) && std::isnan(dfNoData)) ||
545 6 : (singleNDV == dfNoData)))
546 : {
547 4 : return false;
548 : }
549 : }
550 456 : return true;
551 : }
552 :
553 : /************************************************************************/
554 : /* GDALComputedDataset::AddSources() */
555 : /************************************************************************/
556 :
557 329 : void GDALComputedDataset::AddSources(GDALComputedRasterBand *poBand)
558 : {
559 : auto poSourcedRasterBand =
560 329 : cpl::down_cast<VRTSourcedRasterBand *>(m_oVRTDS.GetRasterBand(1));
561 :
562 329 : bool hasAtLeastOneNDV = false;
563 329 : double singleNDV = 0;
564 658 : std::vector<GDALRasterBand *> apoSrcBands;
565 1453 : for (auto &[_, band] : m_apoBands)
566 : {
567 1124 : apoSrcBands.push_back(band);
568 : }
569 329 : const bool bSameNDV = HaveAllBandsSameNoDataValue(
570 : apoSrcBands.data(), m_apoBands.size(), hasAtLeastOneNDV, singleNDV);
571 :
572 658 : std::set<std::string> alreadyAdded;
573 1453 : for (auto &[name, band] : m_apoBands)
574 : {
575 1124 : if (alreadyAdded.insert(name).second)
576 : {
577 480 : int bHasNoData = false;
578 480 : const double dfNoData = band->GetNoDataValue(&bHasNoData);
579 480 : if (bHasNoData)
580 : {
581 9 : poSourcedRasterBand->AddComplexSource(
582 : band, -1, -1, -1, -1, -1, -1, -1, -1, 0, 1, dfNoData);
583 : }
584 : else
585 : {
586 471 : poSourcedRasterBand->AddSimpleSource(band);
587 : }
588 480 : poSourcedRasterBand->m_papoSources.back()->SetName(name.c_str());
589 : }
590 : }
591 329 : if (hasAtLeastOneNDV)
592 : {
593 5 : poBand->m_bHasNoData = true;
594 5 : if (bSameNDV)
595 : {
596 2 : poBand->m_dfNoDataValue = singleNDV;
597 : }
598 : else
599 : {
600 3 : poBand->m_dfNoDataValue = std::numeric_limits<double>::quiet_NaN();
601 : }
602 5 : poSourcedRasterBand->SetNoDataValue(poBand->m_dfNoDataValue);
603 : }
604 329 : }
605 :
606 : /************************************************************************/
607 : /* OperationToFunctionName() */
608 : /************************************************************************/
609 :
610 360 : /* static */ const char *GDALComputedDataset::OperationToFunctionName(
611 : GDALComputedRasterBand::Operation op, bool bForceMuparser)
612 : {
613 360 : const char *ret = "";
614 360 : switch (op)
615 : {
616 64 : case GDALComputedRasterBand::Operation::OP_ADD:
617 64 : ret = bForceMuparser ? "+" : "sum";
618 64 : break;
619 14 : case GDALComputedRasterBand::Operation::OP_SUBTRACT:
620 14 : ret = bForceMuparser ? "-" : "diff";
621 14 : break;
622 76 : case GDALComputedRasterBand::Operation::OP_MULTIPLY:
623 76 : ret = bForceMuparser ? "*" : "mul";
624 76 : break;
625 0 : case GDALComputedRasterBand::Operation::OP_DIVIDE:
626 0 : ret = bForceMuparser ? "/" : "div";
627 0 : break;
628 15 : case GDALComputedRasterBand::Operation::OP_MIN:
629 15 : ret = "min";
630 15 : break;
631 16 : case GDALComputedRasterBand::Operation::OP_MAX:
632 16 : ret = "max";
633 16 : break;
634 4 : case GDALComputedRasterBand::Operation::OP_MEAN:
635 4 : ret = "mean";
636 4 : break;
637 13 : case GDALComputedRasterBand::Operation::OP_GT:
638 13 : ret = ">";
639 13 : break;
640 16 : case GDALComputedRasterBand::Operation::OP_GE:
641 16 : ret = ">=";
642 16 : break;
643 13 : case GDALComputedRasterBand::Operation::OP_LT:
644 13 : ret = "<";
645 13 : break;
646 14 : case GDALComputedRasterBand::Operation::OP_LE:
647 14 : ret = "<=";
648 14 : break;
649 18 : case GDALComputedRasterBand::Operation::OP_EQ:
650 18 : ret = "==";
651 18 : break;
652 20 : case GDALComputedRasterBand::Operation::OP_NE:
653 20 : ret = "!=";
654 20 : break;
655 18 : case GDALComputedRasterBand::Operation::OP_LOGICAL_AND:
656 18 : ret = "&&";
657 18 : break;
658 23 : case GDALComputedRasterBand::Operation::OP_LOGICAL_OR:
659 23 : ret = "||";
660 23 : break;
661 24 : case GDALComputedRasterBand::Operation::OP_CAST:
662 : case GDALComputedRasterBand::Operation::OP_TERNARY:
663 24 : break;
664 4 : case GDALComputedRasterBand::Operation::OP_ABS:
665 4 : ret = bForceMuparser ? "abs" : "mod";
666 4 : break;
667 4 : case GDALComputedRasterBand::Operation::OP_SQRT:
668 4 : ret = "sqrt";
669 4 : break;
670 0 : case GDALComputedRasterBand::Operation::OP_LOG:
671 0 : ret = "log";
672 0 : break;
673 4 : case GDALComputedRasterBand::Operation::OP_LOG10:
674 4 : ret = "log10";
675 4 : break;
676 0 : case GDALComputedRasterBand::Operation::OP_POW:
677 0 : ret = "pow";
678 0 : break;
679 : }
680 360 : return ret;
681 : }
682 :
683 : /************************************************************************/
684 : /* GDALComputedRasterBand() */
685 : /************************************************************************/
686 :
687 0 : GDALComputedRasterBand::GDALComputedRasterBand(
688 0 : const GDALComputedRasterBand &other, bool)
689 0 : : GDALRasterBand()
690 : {
691 0 : nRasterXSize = other.nRasterXSize;
692 0 : nRasterYSize = other.nRasterYSize;
693 0 : eDataType = other.eDataType;
694 0 : nBlockXSize = other.nBlockXSize;
695 0 : nBlockYSize = other.nBlockYSize;
696 0 : }
697 :
698 : //! @cond Doxygen_Suppress
699 :
700 : /************************************************************************/
701 : /* GDALComputedRasterBand() */
702 : /************************************************************************/
703 :
704 37 : GDALComputedRasterBand::GDALComputedRasterBand(
705 : Operation op, const std::vector<const GDALRasterBand *> &bands,
706 37 : double constant)
707 : {
708 37 : CPLAssert(op == Operation::OP_ADD || op == Operation::OP_MIN ||
709 : op == Operation::OP_MAX || op == Operation::OP_MEAN ||
710 : op == Operation::OP_TERNARY);
711 :
712 37 : CPLAssert(!bands.empty());
713 37 : nRasterXSize = bands[0]->GetXSize();
714 37 : nRasterYSize = bands[0]->GetYSize();
715 37 : eDataType = bands[0]->GetRasterDataType();
716 95 : for (size_t i = 1; i < bands.size(); ++i)
717 : {
718 58 : eDataType = GDALDataTypeUnion(eDataType, bands[i]->GetRasterDataType());
719 : }
720 :
721 37 : bool hasAtLeastOneNDV = false;
722 37 : double singleNDV = 0;
723 : const bool bSameNDV =
724 37 : HaveAllBandsSameNoDataValue(const_cast<GDALRasterBand **>(bands.data()),
725 : bands.size(), hasAtLeastOneNDV, singleNDV);
726 :
727 37 : if (!bSameNDV)
728 : {
729 0 : eDataType = eDataType == GDT_Float64 ? GDT_Float64 : GDT_Float32;
730 : }
731 37 : else if (op == Operation::OP_TERNARY)
732 : {
733 18 : CPLAssert(bands.size() == 3);
734 18 : eDataType = GDALDataTypeUnion(bands[1]->GetRasterDataType(),
735 18 : bands[2]->GetRasterDataType());
736 : }
737 19 : else if (!std::isnan(constant) && eDataType != GDT_Float64)
738 : {
739 4 : if (op == Operation::OP_MIN || op == Operation::OP_MAX)
740 : {
741 4 : eDataType = GDALDataTypeUnionWithValue(eDataType, constant, false);
742 : }
743 : else
744 : {
745 0 : eDataType =
746 0 : (static_cast<double>(static_cast<float>(constant)) == constant)
747 0 : ? GDT_Float32
748 : : GDT_Float64;
749 : }
750 : }
751 37 : bands[0]->GetBlockSize(&nBlockXSize, &nBlockYSize);
752 : auto l_poDS = std::make_unique<GDALComputedDataset>(
753 37 : this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
754 74 : op, bands, constant);
755 37 : m_poOwningDS.reset(l_poDS.release());
756 37 : }
757 :
758 : /************************************************************************/
759 : /* GDALComputedRasterBand() */
760 : /************************************************************************/
761 :
762 96 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
763 : const GDALRasterBand &firstBand,
764 96 : const GDALRasterBand &secondBand)
765 : {
766 96 : nRasterXSize = firstBand.GetXSize();
767 96 : nRasterYSize = firstBand.GetYSize();
768 :
769 96 : bool hasAtLeastOneNDV = false;
770 96 : double singleNDV = 0;
771 : GDALRasterBand *apoBands[] = {const_cast<GDALRasterBand *>(&firstBand),
772 96 : const_cast<GDALRasterBand *>(&secondBand)};
773 : const bool bSameNDV =
774 96 : HaveAllBandsSameNoDataValue(apoBands, 2, hasAtLeastOneNDV, singleNDV);
775 :
776 96 : const auto firstDT = firstBand.GetRasterDataType();
777 96 : const auto secondDT = secondBand.GetRasterDataType();
778 96 : if (!bSameNDV)
779 6 : eDataType = (firstDT == GDT_Float64 || secondDT == GDT_Float64)
780 6 : ? GDT_Float64
781 : : GDT_Float32;
782 93 : else if (IsComparisonOperator(op))
783 43 : eDataType = GDT_UInt8;
784 50 : else if (op == Operation::OP_ADD && firstDT == GDT_UInt8 &&
785 : secondDT == GDT_UInt8)
786 3 : eDataType = GDT_UInt16;
787 47 : else if (firstDT == GDT_Float32 && secondDT == GDT_Float32)
788 5 : eDataType = GDT_Float32;
789 42 : else if ((op == Operation::OP_MIN || op == Operation::OP_MAX) &&
790 : firstDT == secondDT)
791 1 : eDataType = firstDT;
792 : else
793 41 : eDataType = GDT_Float64;
794 96 : firstBand.GetBlockSize(&nBlockXSize, &nBlockYSize);
795 : auto l_poDS = std::make_unique<GDALComputedDataset>(
796 96 : this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
797 192 : op, &firstBand, nullptr, &secondBand, nullptr);
798 96 : m_poOwningDS.reset(l_poDS.release());
799 96 : }
800 :
801 : /************************************************************************/
802 : /* GDALComputedRasterBand() */
803 : /************************************************************************/
804 :
805 22 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op, double constant,
806 22 : const GDALRasterBand &band)
807 : {
808 22 : CPLAssert(op == Operation::OP_DIVIDE || IsComparisonOperator(op) ||
809 : op == Operation::OP_POW);
810 :
811 22 : nRasterXSize = band.GetXSize();
812 22 : nRasterYSize = band.GetYSize();
813 22 : const auto firstDT = band.GetRasterDataType();
814 22 : if (IsComparisonOperator(op))
815 18 : eDataType = GDT_UInt8;
816 4 : else if (firstDT == GDT_Float32 &&
817 0 : static_cast<double>(static_cast<float>(constant)) == constant)
818 0 : eDataType = GDT_Float32;
819 : else
820 4 : eDataType = GDT_Float64;
821 22 : band.GetBlockSize(&nBlockXSize, &nBlockYSize);
822 : auto l_poDS = std::make_unique<GDALComputedDataset>(
823 22 : this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
824 44 : op, nullptr, &constant, &band, nullptr);
825 22 : m_poOwningDS.reset(l_poDS.release());
826 22 : }
827 :
828 : /************************************************************************/
829 : /* GDALComputedRasterBand() */
830 : /************************************************************************/
831 :
832 141 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
833 : const GDALRasterBand &band,
834 141 : double constant)
835 : {
836 141 : nRasterXSize = band.GetXSize();
837 141 : nRasterYSize = band.GetYSize();
838 141 : const auto firstDT = band.GetRasterDataType();
839 141 : if (IsComparisonOperator(op))
840 74 : eDataType = GDT_UInt8;
841 67 : else if (op == Operation::OP_ADD && firstDT == GDT_UInt8 &&
842 10 : constant >= -128 && constant <= 127 &&
843 10 : std::floor(constant) == constant)
844 10 : eDataType = GDT_UInt8;
845 57 : else if (firstDT == GDT_Float32 &&
846 8 : static_cast<double>(static_cast<float>(constant)) == constant)
847 8 : eDataType = GDT_Float32;
848 : else
849 49 : eDataType = GDT_Float64;
850 141 : band.GetBlockSize(&nBlockXSize, &nBlockYSize);
851 : auto l_poDS = std::make_unique<GDALComputedDataset>(
852 141 : this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
853 282 : op, &band, nullptr, nullptr, &constant);
854 141 : m_poOwningDS.reset(l_poDS.release());
855 141 : }
856 :
857 : /************************************************************************/
858 : /* GDALComputedRasterBand() */
859 : /************************************************************************/
860 :
861 9 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
862 9 : const GDALRasterBand &band)
863 : {
864 9 : CPLAssert(op == Operation::OP_ABS || op == Operation::OP_SQRT ||
865 : op == Operation::OP_LOG || op == Operation::OP_LOG10);
866 9 : nRasterXSize = band.GetXSize();
867 9 : nRasterYSize = band.GetYSize();
868 9 : eDataType =
869 9 : band.GetRasterDataType() == GDT_Float32 ? GDT_Float32 : GDT_Float64;
870 9 : band.GetBlockSize(&nBlockXSize, &nBlockYSize);
871 : auto l_poDS = std::make_unique<GDALComputedDataset>(
872 9 : this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
873 18 : op, &band, nullptr, nullptr, nullptr);
874 9 : m_poOwningDS.reset(l_poDS.release());
875 9 : }
876 :
877 : /************************************************************************/
878 : /* GDALComputedRasterBand() */
879 : /************************************************************************/
880 :
881 24 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
882 : const GDALRasterBand &band,
883 24 : GDALDataType dt)
884 : {
885 24 : CPLAssert(op == Operation::OP_CAST);
886 24 : nRasterXSize = band.GetXSize();
887 24 : nRasterYSize = band.GetYSize();
888 24 : eDataType = dt;
889 24 : band.GetBlockSize(&nBlockXSize, &nBlockYSize);
890 : auto l_poDS = std::make_unique<GDALComputedDataset>(
891 24 : this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
892 48 : op, &band, nullptr, nullptr, nullptr);
893 24 : m_poOwningDS.reset(l_poDS.release());
894 24 : }
895 :
896 : //! @endcond
897 :
898 : /************************************************************************/
899 : /* ~GDALComputedRasterBand() */
900 : /************************************************************************/
901 :
902 500 : GDALComputedRasterBand::~GDALComputedRasterBand()
903 : {
904 329 : if (m_poOwningDS)
905 329 : cpl::down_cast<GDALComputedDataset *>(m_poOwningDS.get())->nBands = 0;
906 329 : poDS = nullptr;
907 500 : }
908 :
909 : /************************************************************************/
910 : /* GDALComputedRasterBand::GetNoDataValue() */
911 : /************************************************************************/
912 :
913 480 : double GDALComputedRasterBand::GetNoDataValue(int *pbHasNoData)
914 : {
915 480 : if (pbHasNoData)
916 480 : *pbHasNoData = m_bHasNoData;
917 480 : return m_dfNoDataValue;
918 : }
919 :
920 : /************************************************************************/
921 : /* GDALComputedRasterBandRelease() */
922 : /************************************************************************/
923 :
924 : /** Release a GDALComputedRasterBandH
925 : *
926 : * @since 3.12
927 : */
928 171 : void GDALComputedRasterBandRelease(GDALComputedRasterBandH hBand)
929 : {
930 171 : delete GDALComputedRasterBand::FromHandle(hBand);
931 171 : }
932 :
933 : /************************************************************************/
934 : /* IReadBlock() */
935 : /************************************************************************/
936 :
937 247 : CPLErr GDALComputedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
938 : void *pData)
939 : {
940 247 : auto l_poDS = cpl::down_cast<GDALComputedDataset *>(poDS);
941 247 : return l_poDS->m_oVRTDS.GetRasterBand(1)->ReadBlock(nBlockXOff, nBlockYOff,
942 247 : pData);
943 : }
944 :
945 : /************************************************************************/
946 : /* IRasterIO() */
947 : /************************************************************************/
948 :
949 7 : CPLErr GDALComputedRasterBand::IRasterIO(
950 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
951 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
952 : GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg)
953 : {
954 7 : auto l_poDS = cpl::down_cast<GDALComputedDataset *>(poDS);
955 7 : return l_poDS->m_oVRTDS.GetRasterBand(1)->RasterIO(
956 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
957 7 : eBufType, nPixelSpace, nLineSpace, psExtraArg);
958 : }
|
