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