Coverage Report

Created: 2025-06-22 06:59

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