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 292 : }
409 :
410 : /************************************************************************/
411 : /* GDALComputedDataset() */
412 : /************************************************************************/
413 :
414 37 : GDALComputedDataset::GDALComputedDataset(
415 : GDALComputedRasterBand *poBand, int nXSize, int nYSize, GDALDataType eDT,
416 : int nBlockXSize, int nBlockYSize, GDALComputedRasterBand::Operation op,
417 37 : const std::vector<const GDALRasterBand *> &bands, double constant)
418 37 : : m_oVRTDS(nXSize, nYSize, nBlockXSize, nBlockYSize)
419 : {
420 37 : bool bCanUseBuiltin = true;
421 74 : std::vector<std::string> aosNames;
422 132 : for (const GDALRasterBand *poIterBand : bands)
423 : {
424 190 : const auto [bIsExprBand, name] = AddSourceBand(poIterBand);
425 95 : bCanUseBuiltin = bCanUseBuiltin && !bIsExprBand;
426 95 : if (!bIsExprBand)
427 66 : aosNames.push_back(name);
428 : else
429 29 : aosNames.push_back(std::string("(").append(name).append(")"));
430 : }
431 :
432 37 : nRasterXSize = nXSize;
433 37 : nRasterYSize = nYSize;
434 :
435 37 : if (auto poSrcDS = m_apoBands.front().second->GetDataset())
436 : {
437 37 : GDALGeoTransform gt;
438 37 : if (poSrcDS->GetGeoTransform(gt) == CE_None)
439 : {
440 16 : m_oVRTDS.SetGeoTransform(gt);
441 : }
442 :
443 37 : if (const auto *poSRS = poSrcDS->GetSpatialRef())
444 : {
445 16 : m_oVRTDS.SetSpatialRef(poSRS);
446 : }
447 : }
448 :
449 37 : m_aosOptions.SetNameValue("subclass", "VRTDerivedRasterBand");
450 37 : if (op == GDALComputedRasterBand::Operation::OP_TERNARY)
451 : {
452 18 : m_expr = aosNames[0];
453 18 : m_expr += " ? ";
454 18 : m_expr += aosNames[1];
455 18 : m_expr += " : ";
456 18 : m_expr += aosNames[2];
457 :
458 18 : m_aosOptions.SetNameValue("PixelFunctionType", "expression");
459 18 : m_aosOptions.SetNameValue("_PIXELFN_ARG_expression", m_expr.c_str());
460 : }
461 : else
462 : {
463 19 : m_expr = OperationToFunctionName(op);
464 19 : m_expr += '(';
465 19 : bool first = true;
466 60 : for (const auto &name : aosNames)
467 : {
468 41 : if (!first)
469 22 : m_expr += ", ";
470 41 : m_expr += name;
471 41 : first = false;
472 : }
473 19 : if (!std::isnan(constant))
474 : {
475 8 : if (!first)
476 8 : m_expr += ", ";
477 8 : m_expr += CPLSPrintf("%.17g", constant);
478 : }
479 19 : m_expr += ')';
480 :
481 19 : if (bCanUseBuiltin)
482 : {
483 : m_aosOptions.SetNameValue("PixelFunctionType",
484 14 : OperationToFunctionName(op));
485 14 : if (!std::isnan(constant))
486 : {
487 : m_aosOptions.SetNameValue("_PIXELFN_ARG_k",
488 4 : CPLSPrintf("%.17g", constant));
489 : }
490 14 : m_aosOptions.SetNameValue("_PIXELFN_ARG_propagateNoData", "true");
491 : }
492 : else
493 : {
494 5 : m_aosOptions.SetNameValue("PixelFunctionType", "expression");
495 : m_aosOptions.SetNameValue("_PIXELFN_ARG_expression",
496 5 : m_expr.c_str());
497 : }
498 : }
499 37 : m_oVRTDS.AddBand(eDT, m_aosOptions.List());
500 :
501 37 : SetBand(1, poBand);
502 :
503 37 : AddSources(poBand);
504 37 : }
505 :
506 : /************************************************************************/
507 : /* ~GDALComputedDataset() */
508 : /************************************************************************/
509 :
510 : GDALComputedDataset::~GDALComputedDataset() = default;
511 :
512 : /************************************************************************/
513 : /* HaveAllBandsSameNoDataValue() */
514 : /************************************************************************/
515 :
516 462 : static bool HaveAllBandsSameNoDataValue(GDALRasterBand **apoBands,
517 : size_t nBands, bool &hasAtLeastOneNDV,
518 : double &singleNDV)
519 : {
520 462 : hasAtLeastOneNDV = false;
521 462 : singleNDV = 0;
522 :
523 462 : int bFirstBandHasNoData = false;
524 1867 : for (size_t i = 0; i < nBands; ++i)
525 : {
526 1411 : int bHasNoData = false;
527 1411 : const double dfNoData = apoBands[i]->GetNoDataValue(&bHasNoData);
528 1411 : if (bHasNoData)
529 18 : hasAtLeastOneNDV = true;
530 1411 : if (i == 0)
531 : {
532 462 : bFirstBandHasNoData = bHasNoData;
533 462 : singleNDV = dfNoData;
534 : }
535 949 : else if (bHasNoData != bFirstBandHasNoData)
536 : {
537 6 : return false;
538 : }
539 955 : else if (bFirstBandHasNoData &&
540 8 : !((std::isnan(singleNDV) && std::isnan(dfNoData)) ||
541 6 : (singleNDV == dfNoData)))
542 : {
543 4 : return false;
544 : }
545 : }
546 456 : return true;
547 : }
548 :
549 : /************************************************************************/
550 : /* GDALComputedDataset::AddSources() */
551 : /************************************************************************/
552 :
553 329 : void GDALComputedDataset::AddSources(GDALComputedRasterBand *poBand)
554 : {
555 : auto poSourcedRasterBand =
556 329 : cpl::down_cast<VRTSourcedRasterBand *>(m_oVRTDS.GetRasterBand(1));
557 :
558 329 : bool hasAtLeastOneNDV = false;
559 329 : double singleNDV = 0;
560 658 : std::vector<GDALRasterBand *> apoSrcBands;
561 1453 : for (auto &[_, band] : m_apoBands)
562 : {
563 1124 : apoSrcBands.push_back(band);
564 : }
565 329 : const bool bSameNDV = HaveAllBandsSameNoDataValue(
566 : apoSrcBands.data(), m_apoBands.size(), hasAtLeastOneNDV, singleNDV);
567 :
568 658 : std::set<std::string> alreadyAdded;
569 1453 : for (auto &[name, band] : m_apoBands)
570 : {
571 1124 : if (alreadyAdded.insert(name).second)
572 : {
573 480 : int bHasNoData = false;
574 480 : const double dfNoData = band->GetNoDataValue(&bHasNoData);
575 480 : if (bHasNoData)
576 : {
577 9 : poSourcedRasterBand->AddComplexSource(
578 : band, -1, -1, -1, -1, -1, -1, -1, -1, 0, 1, dfNoData);
579 : }
580 : else
581 : {
582 471 : poSourcedRasterBand->AddSimpleSource(band);
583 : }
584 480 : poSourcedRasterBand->m_papoSources.back()->SetName(name.c_str());
585 : }
586 : }
587 329 : if (hasAtLeastOneNDV)
588 : {
589 5 : poBand->m_bHasNoData = true;
590 5 : if (bSameNDV)
591 : {
592 2 : poBand->m_dfNoDataValue = singleNDV;
593 : }
594 : else
595 : {
596 3 : poBand->m_dfNoDataValue = std::numeric_limits<double>::quiet_NaN();
597 : }
598 5 : poSourcedRasterBand->SetNoDataValue(poBand->m_dfNoDataValue);
599 : }
600 329 : }
601 :
602 : /************************************************************************/
603 : /* OperationToFunctionName() */
604 : /************************************************************************/
605 :
606 360 : /* static */ const char *GDALComputedDataset::OperationToFunctionName(
607 : GDALComputedRasterBand::Operation op, bool bForceMuparser)
608 : {
609 360 : const char *ret = "";
610 360 : switch (op)
611 : {
612 64 : case GDALComputedRasterBand::Operation::OP_ADD:
613 64 : ret = bForceMuparser ? "+" : "sum";
614 64 : break;
615 14 : case GDALComputedRasterBand::Operation::OP_SUBTRACT:
616 14 : ret = bForceMuparser ? "-" : "diff";
617 14 : break;
618 76 : case GDALComputedRasterBand::Operation::OP_MULTIPLY:
619 76 : ret = bForceMuparser ? "*" : "mul";
620 76 : break;
621 0 : case GDALComputedRasterBand::Operation::OP_DIVIDE:
622 0 : ret = bForceMuparser ? "/" : "div";
623 0 : break;
624 15 : case GDALComputedRasterBand::Operation::OP_MIN:
625 15 : ret = "min";
626 15 : break;
627 16 : case GDALComputedRasterBand::Operation::OP_MAX:
628 16 : ret = "max";
629 16 : break;
630 4 : case GDALComputedRasterBand::Operation::OP_MEAN:
631 4 : ret = "mean";
632 4 : break;
633 13 : case GDALComputedRasterBand::Operation::OP_GT:
634 13 : ret = ">";
635 13 : break;
636 16 : case GDALComputedRasterBand::Operation::OP_GE:
637 16 : ret = ">=";
638 16 : break;
639 13 : case GDALComputedRasterBand::Operation::OP_LT:
640 13 : ret = "<";
641 13 : break;
642 14 : case GDALComputedRasterBand::Operation::OP_LE:
643 14 : ret = "<=";
644 14 : break;
645 18 : case GDALComputedRasterBand::Operation::OP_EQ:
646 18 : ret = "==";
647 18 : break;
648 20 : case GDALComputedRasterBand::Operation::OP_NE:
649 20 : ret = "!=";
650 20 : break;
651 18 : case GDALComputedRasterBand::Operation::OP_LOGICAL_AND:
652 18 : ret = "&&";
653 18 : break;
654 23 : case GDALComputedRasterBand::Operation::OP_LOGICAL_OR:
655 23 : ret = "||";
656 23 : break;
657 24 : case GDALComputedRasterBand::Operation::OP_CAST:
658 : case GDALComputedRasterBand::Operation::OP_TERNARY:
659 24 : break;
660 4 : case GDALComputedRasterBand::Operation::OP_ABS:
661 4 : ret = bForceMuparser ? "abs" : "mod";
662 4 : break;
663 4 : case GDALComputedRasterBand::Operation::OP_SQRT:
664 4 : ret = "sqrt";
665 4 : break;
666 0 : case GDALComputedRasterBand::Operation::OP_LOG:
667 0 : ret = "log";
668 0 : break;
669 4 : case GDALComputedRasterBand::Operation::OP_LOG10:
670 4 : ret = "log10";
671 4 : break;
672 0 : case GDALComputedRasterBand::Operation::OP_POW:
673 0 : ret = "pow";
674 0 : break;
675 : }
676 360 : return ret;
677 : }
678 :
679 : /************************************************************************/
680 : /* GDALComputedRasterBand() */
681 : /************************************************************************/
682 :
683 0 : GDALComputedRasterBand::GDALComputedRasterBand(
684 0 : const GDALComputedRasterBand &other, bool)
685 0 : : GDALRasterBand()
686 : {
687 0 : nRasterXSize = other.nRasterXSize;
688 0 : nRasterYSize = other.nRasterYSize;
689 0 : eDataType = other.eDataType;
690 0 : nBlockXSize = other.nBlockXSize;
691 0 : nBlockYSize = other.nBlockYSize;
692 0 : }
693 :
694 : //! @cond Doxygen_Suppress
695 :
696 : /************************************************************************/
697 : /* GDALComputedRasterBand() */
698 : /************************************************************************/
699 :
700 37 : GDALComputedRasterBand::GDALComputedRasterBand(
701 : Operation op, const std::vector<const GDALRasterBand *> &bands,
702 37 : double constant)
703 : {
704 37 : CPLAssert(op == Operation::OP_ADD || op == Operation::OP_MIN ||
705 : op == Operation::OP_MAX || op == Operation::OP_MEAN ||
706 : op == Operation::OP_TERNARY);
707 :
708 37 : CPLAssert(!bands.empty());
709 37 : nRasterXSize = bands[0]->GetXSize();
710 37 : nRasterYSize = bands[0]->GetYSize();
711 37 : eDataType = bands[0]->GetRasterDataType();
712 95 : for (size_t i = 1; i < bands.size(); ++i)
713 : {
714 58 : eDataType = GDALDataTypeUnion(eDataType, bands[i]->GetRasterDataType());
715 : }
716 :
717 37 : bool hasAtLeastOneNDV = false;
718 37 : double singleNDV = 0;
719 : const bool bSameNDV =
720 37 : HaveAllBandsSameNoDataValue(const_cast<GDALRasterBand **>(bands.data()),
721 : bands.size(), hasAtLeastOneNDV, singleNDV);
722 :
723 37 : if (!bSameNDV)
724 : {
725 0 : eDataType = eDataType == GDT_Float64 ? GDT_Float64 : GDT_Float32;
726 : }
727 37 : else if (op == Operation::OP_TERNARY)
728 : {
729 18 : CPLAssert(bands.size() == 3);
730 18 : eDataType = GDALDataTypeUnion(bands[1]->GetRasterDataType(),
731 18 : bands[2]->GetRasterDataType());
732 : }
733 19 : else if (!std::isnan(constant) && eDataType != GDT_Float64)
734 : {
735 4 : if (op == Operation::OP_MIN || op == Operation::OP_MAX)
736 : {
737 4 : eDataType = GDALDataTypeUnionWithValue(eDataType, constant, false);
738 : }
739 : else
740 : {
741 0 : eDataType =
742 0 : (static_cast<double>(static_cast<float>(constant)) == constant)
743 0 : ? GDT_Float32
744 : : GDT_Float64;
745 : }
746 : }
747 37 : bands[0]->GetBlockSize(&nBlockXSize, &nBlockYSize);
748 : auto l_poDS = std::make_unique<GDALComputedDataset>(
749 37 : this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
750 74 : op, bands, constant);
751 37 : m_poOwningDS.reset(l_poDS.release());
752 37 : }
753 :
754 : /************************************************************************/
755 : /* GDALComputedRasterBand() */
756 : /************************************************************************/
757 :
758 96 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
759 : const GDALRasterBand &firstBand,
760 96 : const GDALRasterBand &secondBand)
761 : {
762 96 : nRasterXSize = firstBand.GetXSize();
763 96 : nRasterYSize = firstBand.GetYSize();
764 :
765 96 : bool hasAtLeastOneNDV = false;
766 96 : double singleNDV = 0;
767 : GDALRasterBand *apoBands[] = {const_cast<GDALRasterBand *>(&firstBand),
768 96 : const_cast<GDALRasterBand *>(&secondBand)};
769 : const bool bSameNDV =
770 96 : HaveAllBandsSameNoDataValue(apoBands, 2, hasAtLeastOneNDV, singleNDV);
771 :
772 96 : const auto firstDT = firstBand.GetRasterDataType();
773 96 : const auto secondDT = secondBand.GetRasterDataType();
774 96 : if (!bSameNDV)
775 6 : eDataType = (firstDT == GDT_Float64 || secondDT == GDT_Float64)
776 6 : ? GDT_Float64
777 : : GDT_Float32;
778 93 : else if (IsComparisonOperator(op))
779 43 : eDataType = GDT_UInt8;
780 50 : else if (op == Operation::OP_ADD && firstDT == GDT_UInt8 &&
781 : secondDT == GDT_UInt8)
782 3 : eDataType = GDT_UInt16;
783 47 : else if (firstDT == GDT_Float32 && secondDT == GDT_Float32)
784 5 : eDataType = GDT_Float32;
785 42 : else if ((op == Operation::OP_MIN || op == Operation::OP_MAX) &&
786 : firstDT == secondDT)
787 1 : eDataType = firstDT;
788 : else
789 41 : eDataType = GDT_Float64;
790 96 : firstBand.GetBlockSize(&nBlockXSize, &nBlockYSize);
791 : auto l_poDS = std::make_unique<GDALComputedDataset>(
792 96 : this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
793 192 : op, &firstBand, nullptr, &secondBand, nullptr);
794 96 : m_poOwningDS.reset(l_poDS.release());
795 96 : }
796 :
797 : /************************************************************************/
798 : /* GDALComputedRasterBand() */
799 : /************************************************************************/
800 :
801 22 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op, double constant,
802 22 : const GDALRasterBand &band)
803 : {
804 22 : CPLAssert(op == Operation::OP_DIVIDE || IsComparisonOperator(op) ||
805 : op == Operation::OP_POW);
806 :
807 22 : nRasterXSize = band.GetXSize();
808 22 : nRasterYSize = band.GetYSize();
809 22 : const auto firstDT = band.GetRasterDataType();
810 22 : if (IsComparisonOperator(op))
811 18 : eDataType = GDT_UInt8;
812 4 : else if (firstDT == GDT_Float32 &&
813 0 : static_cast<double>(static_cast<float>(constant)) == constant)
814 0 : eDataType = GDT_Float32;
815 : else
816 4 : eDataType = GDT_Float64;
817 22 : band.GetBlockSize(&nBlockXSize, &nBlockYSize);
818 : auto l_poDS = std::make_unique<GDALComputedDataset>(
819 22 : this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
820 44 : op, nullptr, &constant, &band, nullptr);
821 22 : m_poOwningDS.reset(l_poDS.release());
822 22 : }
823 :
824 : /************************************************************************/
825 : /* GDALComputedRasterBand() */
826 : /************************************************************************/
827 :
828 141 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
829 : const GDALRasterBand &band,
830 141 : double constant)
831 : {
832 141 : nRasterXSize = band.GetXSize();
833 141 : nRasterYSize = band.GetYSize();
834 141 : const auto firstDT = band.GetRasterDataType();
835 141 : if (IsComparisonOperator(op))
836 74 : eDataType = GDT_UInt8;
837 67 : else if (op == Operation::OP_ADD && firstDT == GDT_UInt8 &&
838 10 : constant >= -128 && constant <= 127 &&
839 10 : std::floor(constant) == constant)
840 10 : eDataType = GDT_UInt8;
841 57 : else if (firstDT == GDT_Float32 &&
842 8 : static_cast<double>(static_cast<float>(constant)) == constant)
843 8 : eDataType = GDT_Float32;
844 : else
845 49 : eDataType = GDT_Float64;
846 141 : band.GetBlockSize(&nBlockXSize, &nBlockYSize);
847 : auto l_poDS = std::make_unique<GDALComputedDataset>(
848 141 : this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
849 282 : op, &band, nullptr, nullptr, &constant);
850 141 : m_poOwningDS.reset(l_poDS.release());
851 141 : }
852 :
853 : /************************************************************************/
854 : /* GDALComputedRasterBand() */
855 : /************************************************************************/
856 :
857 9 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
858 9 : const GDALRasterBand &band)
859 : {
860 9 : CPLAssert(op == Operation::OP_ABS || op == Operation::OP_SQRT ||
861 : op == Operation::OP_LOG || op == Operation::OP_LOG10);
862 9 : nRasterXSize = band.GetXSize();
863 9 : nRasterYSize = band.GetYSize();
864 9 : eDataType =
865 9 : band.GetRasterDataType() == GDT_Float32 ? GDT_Float32 : GDT_Float64;
866 9 : band.GetBlockSize(&nBlockXSize, &nBlockYSize);
867 : auto l_poDS = std::make_unique<GDALComputedDataset>(
868 9 : this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
869 18 : op, &band, nullptr, nullptr, nullptr);
870 9 : m_poOwningDS.reset(l_poDS.release());
871 9 : }
872 :
873 : /************************************************************************/
874 : /* GDALComputedRasterBand() */
875 : /************************************************************************/
876 :
877 24 : GDALComputedRasterBand::GDALComputedRasterBand(Operation op,
878 : const GDALRasterBand &band,
879 24 : GDALDataType dt)
880 : {
881 24 : CPLAssert(op == Operation::OP_CAST);
882 24 : nRasterXSize = band.GetXSize();
883 24 : nRasterYSize = band.GetYSize();
884 24 : eDataType = dt;
885 24 : band.GetBlockSize(&nBlockXSize, &nBlockYSize);
886 : auto l_poDS = std::make_unique<GDALComputedDataset>(
887 24 : this, nRasterXSize, nRasterYSize, eDataType, nBlockXSize, nBlockYSize,
888 48 : op, &band, nullptr, nullptr, nullptr);
889 24 : m_poOwningDS.reset(l_poDS.release());
890 24 : }
891 :
892 : //! @endcond
893 :
894 : /************************************************************************/
895 : /* ~GDALComputedRasterBand() */
896 : /************************************************************************/
897 :
898 500 : GDALComputedRasterBand::~GDALComputedRasterBand()
899 : {
900 329 : if (m_poOwningDS)
901 329 : cpl::down_cast<GDALComputedDataset *>(m_poOwningDS.get())->nBands = 0;
902 329 : poDS = nullptr;
903 500 : }
904 :
905 : /************************************************************************/
906 : /* GDALComputedRasterBand::GetNoDataValue() */
907 : /************************************************************************/
908 :
909 480 : double GDALComputedRasterBand::GetNoDataValue(int *pbHasNoData)
910 : {
911 480 : if (pbHasNoData)
912 480 : *pbHasNoData = m_bHasNoData;
913 480 : return m_dfNoDataValue;
914 : }
915 :
916 : /************************************************************************/
917 : /* GDALComputedRasterBandRelease() */
918 : /************************************************************************/
919 :
920 : /** Release a GDALComputedRasterBandH
921 : *
922 : * @since 3.12
923 : */
924 171 : void GDALComputedRasterBandRelease(GDALComputedRasterBandH hBand)
925 : {
926 171 : delete GDALComputedRasterBand::FromHandle(hBand);
927 171 : }
928 :
929 : /************************************************************************/
930 : /* IReadBlock() */
931 : /************************************************************************/
932 :
933 247 : CPLErr GDALComputedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
934 : void *pData)
935 : {
936 247 : auto l_poDS = cpl::down_cast<GDALComputedDataset *>(poDS);
937 247 : return l_poDS->m_oVRTDS.GetRasterBand(1)->ReadBlock(nBlockXOff, nBlockYOff,
938 247 : pData);
939 : }
940 :
941 : /************************************************************************/
942 : /* IRasterIO() */
943 : /************************************************************************/
944 :
945 7 : CPLErr GDALComputedRasterBand::IRasterIO(
946 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
947 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
948 : GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg)
949 : {
950 7 : auto l_poDS = cpl::down_cast<GDALComputedDataset *>(poDS);
951 7 : return l_poDS->m_oVRTDS.GetRasterBand(1)->RasterIO(
952 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
953 7 : eBufType, nPixelSpace, nLineSpace, psExtraArg);
954 : }
|
