• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

OSGeo / gdal / 15885686134

25 Jun 2025 07:44PM UTC coverage: 71.084%. Remained the same
15885686134

push

github

rouault
gdal_priv.h: fix C++11 compatibility

573814 of 807237 relevant lines covered (71.08%)

250621.56 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

78.17
/frmts/hdf5/hdf5imagedataset.cpp
1
/******************************************************************************
2
 *
3
 * Project:  Hierarchical Data Format Release 5 (HDF5)
4
 * Purpose:  Read subdatasets of HDF5 file.
5
 * Author:   Denis Nadeau <denis.nadeau@gmail.com>
6
 *
7
 ******************************************************************************
8
 * Copyright (c) 2005, Frank Warmerdam <warmerdam@pobox.com>
9
 * Copyright (c) 2007-2013, Even Rouault <even dot rouault at spatialys.com>
10
 *
11
 * SPDX-License-Identifier: MIT
12
 ****************************************************************************/
13

14
#include "hdf5_api.h"
15

16
#include "cpl_float.h"
17
#include "cpl_string.h"
18
#include "gdal_frmts.h"
19
#include "gdal_pam.h"
20
#include "gdal_priv.h"
21
#include "gh5_convenience.h"
22
#include "hdf5dataset.h"
23
#include "hdf5drivercore.h"
24
#include "ogr_spatialref.h"
25
#include "memdataset.h"
26

27
#include <algorithm>
28
#include <limits>
29

30
class HDF5ImageDataset final : public HDF5Dataset
31
{
32
    typedef enum
33
    {
34
        UNKNOWN_PRODUCT = 0,
35
        CSK_PRODUCT
36
    } Hdf5ProductType;
37

38
    typedef enum
39
    {
40
        PROD_UNKNOWN = 0,
41
        PROD_CSK_L0,
42
        PROD_CSK_L1A,
43
        PROD_CSK_L1B,
44
        PROD_CSK_L1C,
45
        PROD_CSK_L1D
46
    } HDF5CSKProductEnum;
47

48
    friend class HDF5ImageRasterBand;
49

50
    OGRSpatialReference m_oSRS{};
51
    OGRSpatialReference m_oGCPSRS{};
52
    std::vector<gdal::GCP> m_aoGCPs{};
53

54
    hsize_t *dims;
55
    hsize_t *maxdims;
56
    HDF5GroupObjects *poH5Objects;
57
    int ndims;
58
    int dimensions;
59
    hid_t dataset_id;
60
    hid_t dataspace_id;
61
    hid_t native;
62
#ifdef HDF5_HAVE_FLOAT16
63
    bool m_bConvertFromFloat16 = false;
64
#endif
65
    Hdf5ProductType iSubdatasetType;
66
    HDF5CSKProductEnum iCSKProductType;
67
    GDALGeoTransform m_gt{};
68
    bool bHasGeoTransform;
69
    int m_nXIndex = -1;
70
    int m_nYIndex = -1;
71
    int m_nOtherDimIndex = -1;
72

73
    int m_nBlockXSize = 0;
74
    int m_nBlockYSize = 0;
75
    int m_nBandChunkSize = 1;  //! Number of bands in a chunk
76

77
    enum WholeBandChunkOptim
78
    {
79
        WBC_DETECTION_IN_PROGRESS,
80
        WBC_DISABLED,
81
        WBC_ENABLED,
82
    };
83

84
    //! Flag to detect if the read pattern of HDF5ImageRasterBand::IRasterIO()
85
    // is whole band after whole band.
86
    WholeBandChunkOptim m_eWholeBandChunkOptim = WBC_DETECTION_IN_PROGRESS;
87
    //! Value of nBand during last HDF5ImageRasterBand::IRasterIO() call
88
    int m_nLastRasterIOBand = -1;
89
    //! Value of nXOff during last HDF5ImageRasterBand::IRasterIO() call
90
    int m_nLastRasterIOXOff = -1;
91
    //! Value of nYOff during last HDF5ImageRasterBand::IRasterIO() call
92
    int m_nLastRasterIOYOff = -1;
93
    //! Value of nXSize during last HDF5ImageRasterBand::IRasterIO() call
94
    int m_nLastRasterIOXSize = -1;
95
    //! Value of nYSize during last HDF5ImageRasterBand::IRasterIO() call
96
    int m_nLastRasterIOYSize = -1;
97
    //! Value such that m_abyBandChunk represent band data in the range
98
    // [m_iCurrentBandChunk * m_nBandChunkSize, (m_iCurrentBandChunk+1) * m_nBandChunkSize[
99
    int m_iCurrentBandChunk = -1;
100
    //! Cached values (in native data type) for bands in the range
101
    // [m_iCurrentBandChunk * m_nBandChunkSize, (m_iCurrentBandChunk+1) * m_nBandChunkSize[
102
    std::vector<GByte> m_abyBandChunk{};
103

104
    CPLErr CreateODIMH5Projection();
105

106
    CPL_DISALLOW_COPY_ASSIGN(HDF5ImageDataset)
107

108
  public:
109
    HDF5ImageDataset();
110
    virtual ~HDF5ImageDataset();
111

112
    CPLErr CreateProjections();
113
    static GDALDataset *Open(GDALOpenInfo *);
114
    static int Identify(GDALOpenInfo *);
115

116
    const OGRSpatialReference *GetSpatialRef() const override;
117
    virtual int GetGCPCount() override;
118
    const OGRSpatialReference *GetGCPSpatialRef() const override;
119
    virtual const GDAL_GCP *GetGCPs() override;
120
    virtual CPLErr GetGeoTransform(GDALGeoTransform &gt) const override;
121

122
    CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,
123
                     int nYSize, void *pData, int nBufXSize, int nBufYSize,
124
                     GDALDataType eBufType, int nBandCount,
125
                     BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
126
                     GSpacing nLineSpace, GSpacing nBandSpace,
127
                     GDALRasterIOExtraArg *psExtraArg) override;
128

129
    const char *GetMetadataItem(const char *pszName,
130
                                const char *pszDomain = "") override;
131

132
    Hdf5ProductType GetSubdatasetType() const
171✔
133
    {
134
        return iSubdatasetType;
171✔
135
    }
136

137
    HDF5CSKProductEnum GetCSKProductType() const
6✔
138
    {
139
        return iCSKProductType;
6✔
140
    }
141

142
    int IsComplexCSKL1A() const
171✔
143
    {
144
        return GetSubdatasetType() == CSK_PRODUCT &&
177✔
145
               GetCSKProductType() == PROD_CSK_L1A && ndims == 3;
177✔
146
    }
147

148
    int GetYIndex() const
640✔
149
    {
150
        return m_nYIndex;
640✔
151
    }
152

153
    int GetXIndex() const
776✔
154
    {
155
        return m_nXIndex;
776✔
156
    }
157

158
    /**
159
     * Identify if the subdataset has a known product format
160
     * It stores a product identifier in iSubdatasetType,
161
     * UNKNOWN_PRODUCT, if it isn't a recognizable format.
162
     */
163
    void IdentifyProductType();
164

165
    /**
166
     * Captures Geolocation information from a COSMO-SKYMED
167
     * file.
168
     * The geoid will always be WGS84
169
     * The projection type may be UTM or UPS, depending on the
170
     * latitude from the center of the image.
171
     * @param iProductType type of HDF5 subproduct, see HDF5CSKProduct
172
     */
173
    void CaptureCSKGeolocation(int iProductType);
174

175
    /**
176
     * Get Geotransform information for COSMO-SKYMED files
177
     * In case of success it stores the transformation
178
     * in m_gt. In case of failure it doesn't
179
     * modify m_gt
180
     * @param iProductType type of HDF5 subproduct, see HDF5CSKProduct
181
     */
182
    void CaptureCSKGeoTransform(int iProductType);
183

184
    /**
185
     * @param iProductType type of HDF5 subproduct, see HDF5CSKProduct
186
     */
187
    void CaptureCSKGCPs(int iProductType);
188
};
189

190
/************************************************************************/
191
/* ==================================================================== */
192
/*                              HDF5ImageDataset                        */
193
/* ==================================================================== */
194
/************************************************************************/
195

196
/************************************************************************/
197
/*                           HDF5ImageDataset()                         */
198
/************************************************************************/
199
HDF5ImageDataset::HDF5ImageDataset()
58✔
200
    : dims(nullptr), maxdims(nullptr), poH5Objects(nullptr), ndims(0),
201
      dimensions(0), dataset_id(-1), dataspace_id(-1), native(-1),
202
      iSubdatasetType(UNKNOWN_PRODUCT), iCSKProductType(PROD_UNKNOWN),
203
      bHasGeoTransform(false)
58✔
204
{
205
    m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
58✔
206
    m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
58✔
207
}
58✔
208

209
/************************************************************************/
210
/*                            ~HDF5ImageDataset()                       */
211
/************************************************************************/
212
HDF5ImageDataset::~HDF5ImageDataset()
116✔
213
{
214
    HDF5_GLOBAL_LOCK();
215

216
    FlushCache(true);
58✔
217

218
    if (dataset_id > 0)
58✔
219
        H5Dclose(dataset_id);
58✔
220
    if (dataspace_id > 0)
58✔
221
        H5Sclose(dataspace_id);
58✔
222
    if (native > 0)
58✔
223
        H5Tclose(native);
57✔
224

225
    CPLFree(dims);
58✔
226
    CPLFree(maxdims);
58✔
227
}
116✔
228

229
/************************************************************************/
230
/* ==================================================================== */
231
/*                            Hdf5imagerasterband                       */
232
/* ==================================================================== */
233
/************************************************************************/
234
class HDF5ImageRasterBand final : public GDALPamRasterBand
235
{
236
    friend class HDF5ImageDataset;
237

238
    bool bNoDataSet = false;
239
    double dfNoDataValue = -9999.0;
240
    bool m_bHasOffset = false;
241
    double m_dfOffset = 0.0;
242
    bool m_bHasScale = false;
243
    double m_dfScale = 1.0;
244
    int m_nIRasterIORecCounter = 0;
245

246
  public:
247
    HDF5ImageRasterBand(HDF5ImageDataset *, int, GDALDataType);
248
    virtual ~HDF5ImageRasterBand();
249

250
    virtual CPLErr IReadBlock(int, int, void *) override;
251
    virtual double GetNoDataValue(int *) override;
252
    virtual double GetOffset(int *) override;
253
    virtual double GetScale(int *) override;
254
    // virtual CPLErr IWriteBlock( int, int, void * );
255

256
    CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,
257
                     int nYSize, void *pData, int nBufXSize, int nBufYSize,
258
                     GDALDataType eBufType, GSpacing nPixelSpace,
259
                     GSpacing nLineSpace,
260
                     GDALRasterIOExtraArg *psExtraArg) override;
261
};
262

263
/************************************************************************/
264
/*                        ~HDF5ImageRasterBand()                        */
265
/************************************************************************/
266

267
HDF5ImageRasterBand::~HDF5ImageRasterBand()
246✔
268
{
269
}
246✔
270

271
/************************************************************************/
272
/*                           HDF5ImageRasterBand()                      */
273
/************************************************************************/
274
HDF5ImageRasterBand::HDF5ImageRasterBand(HDF5ImageDataset *poDSIn, int nBandIn,
123✔
275
                                         GDALDataType eType)
123✔
276
{
277
    poDS = poDSIn;
123✔
278
    nBand = nBandIn;
123✔
279
    eDataType = eType;
123✔
280
    nBlockXSize = poDSIn->m_nBlockXSize;
123✔
281
    nBlockYSize = poDSIn->m_nBlockYSize;
123✔
282

283
    // netCDF convention for nodata
284
    bNoDataSet =
123✔
285
        GH5_FetchAttribute(poDSIn->dataset_id, "_FillValue", dfNoDataValue);
123✔
286
    if (!bNoDataSet)
123✔
287
        dfNoDataValue = -9999.0;
114✔
288

289
    // netCDF conventions for scale and offset
290
    m_bHasOffset =
123✔
291
        GH5_FetchAttribute(poDSIn->dataset_id, "add_offset", m_dfOffset);
123✔
292
    if (!m_bHasOffset)
123✔
293
        m_dfOffset = 0.0;
122✔
294
    m_bHasScale =
123✔
295
        GH5_FetchAttribute(poDSIn->dataset_id, "scale_factor", m_dfScale);
123✔
296
    if (!m_bHasScale)
123✔
297
        m_dfScale = 1.0;
122✔
298
}
123✔
299

300
/************************************************************************/
301
/*                           GetNoDataValue()                           */
302
/************************************************************************/
303
double HDF5ImageRasterBand::GetNoDataValue(int *pbSuccess)
48✔
304

305
{
306
    if (bNoDataSet)
48✔
307
    {
308
        if (pbSuccess)
3✔
309
            *pbSuccess = bNoDataSet;
3✔
310

311
        return dfNoDataValue;
3✔
312
    }
313

314
    return GDALPamRasterBand::GetNoDataValue(pbSuccess);
45✔
315
}
316

317
/************************************************************************/
318
/*                             GetOffset()                              */
319
/************************************************************************/
320

321
double HDF5ImageRasterBand::GetOffset(int *pbSuccess)
22✔
322

323
{
324
    if (m_bHasOffset)
22✔
325
    {
326
        if (pbSuccess)
1✔
327
            *pbSuccess = m_bHasOffset;
1✔
328

329
        return m_dfOffset;
1✔
330
    }
331

332
    return GDALPamRasterBand::GetOffset(pbSuccess);
21✔
333
}
334

335
/************************************************************************/
336
/*                             GetScale()                               */
337
/************************************************************************/
338

339
double HDF5ImageRasterBand::GetScale(int *pbSuccess)
22✔
340

341
{
342
    if (m_bHasScale)
22✔
343
    {
344
        if (pbSuccess)
1✔
345
            *pbSuccess = m_bHasScale;
1✔
346

347
        return m_dfScale;
1✔
348
    }
349

350
    return GDALPamRasterBand::GetScale(pbSuccess);
21✔
351
}
352

353
/************************************************************************/
354
/*                             IReadBlock()                             */
355
/************************************************************************/
356
CPLErr HDF5ImageRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
103✔
357
                                       void *pImage)
358
{
359
    HDF5ImageDataset *poGDS = static_cast<HDF5ImageDataset *>(poDS);
103✔
360

361
    memset(pImage, 0,
103✔
362
           static_cast<size_t>(nBlockXSize) * nBlockYSize *
103✔
363
               GDALGetDataTypeSizeBytes(eDataType));
103✔
364

365
    if (poGDS->eAccess == GA_Update)
103✔
366
    {
367
        return CE_None;
×
368
    }
369

370
    const int nXOff = nBlockXOff * nBlockXSize;
103✔
371
    const int nYOff = nBlockYOff * nBlockYSize;
103✔
372
    const int nXSize = std::min(nBlockXSize, nRasterXSize - nXOff);
103✔
373
    const int nYSize = std::min(nBlockYSize, nRasterYSize - nYOff);
103✔
374
    if (poGDS->m_eWholeBandChunkOptim == HDF5ImageDataset::WBC_ENABLED)
103✔
375
    {
376
        const bool bIsBandInterleavedData =
377
            poGDS->ndims == 3 && poGDS->m_nOtherDimIndex == 0 &&
2✔
378
            poGDS->GetYIndex() == 1 && poGDS->GetXIndex() == 2;
3✔
379
        if (poGDS->nBands == 1 || bIsBandInterleavedData)
1✔
380
        {
381
            GDALRasterIOExtraArg sExtraArg;
382
            INIT_RASTERIO_EXTRA_ARG(sExtraArg);
1✔
383
            const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
1✔
384
            return IRasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize, pImage,
2✔
385
                             nXSize, nYSize, eDataType, nDTSize,
386
                             static_cast<GSpacing>(nDTSize) * nBlockXSize,
1✔
387
                             &sExtraArg);
1✔
388
        }
389
    }
390

391
    HDF5_GLOBAL_LOCK();
392

393
    hsize_t count[3] = {0, 0, 0};
102✔
394
    H5OFFSET_TYPE offset[3] = {0, 0, 0};
102✔
395
    hsize_t col_dims[3] = {0, 0, 0};
102✔
396
    hsize_t rank = std::min(poGDS->ndims, 2);
102✔
397

398
    if (poGDS->ndims == 3)
102✔
399
    {
400
        rank = 3;
1✔
401
        offset[poGDS->m_nOtherDimIndex] = nBand - 1;
1✔
402
        count[poGDS->m_nOtherDimIndex] = 1;
1✔
403
        col_dims[poGDS->m_nOtherDimIndex] = 1;
1✔
404
    }
405

406
    const int nYIndex = poGDS->GetYIndex();
102✔
407
    // Blocksize may not be a multiple of imagesize.
408
    if (nYIndex >= 0)
102✔
409
    {
410
        offset[nYIndex] = nYOff;
101✔
411
        count[nYIndex] = nYSize;
101✔
412
    }
413
    offset[poGDS->GetXIndex()] = nXOff;
102✔
414
    count[poGDS->GetXIndex()] = nXSize;
102✔
415

416
    // Select block from file space.
417
    herr_t status = H5Sselect_hyperslab(poGDS->dataspace_id, H5S_SELECT_SET,
102✔
418
                                        offset, nullptr, count, nullptr);
419
    if (status < 0)
102✔
420
        return CE_Failure;
×
421

422
    // Create memory space to receive the data.
423
    if (nYIndex >= 0)
102✔
424
        col_dims[nYIndex] = nBlockYSize;
101✔
425
    col_dims[poGDS->GetXIndex()] = nBlockXSize;
102✔
426

427
    const hid_t memspace =
428
        H5Screate_simple(static_cast<int>(rank), col_dims, nullptr);
102✔
429
    H5OFFSET_TYPE mem_offset[3] = {0, 0, 0};
102✔
430
    status = H5Sselect_hyperslab(memspace, H5S_SELECT_SET, mem_offset, nullptr,
102✔
431
                                 count, nullptr);
432
    if (status < 0)
102✔
433
    {
434
        H5Sclose(memspace);
×
435
        return CE_Failure;
×
436
    }
437

438
    status = H5Dread(poGDS->dataset_id, poGDS->native, memspace,
102✔
439
                     poGDS->dataspace_id, H5P_DEFAULT, pImage);
440

441
    H5Sclose(memspace);
102✔
442

443
    if (status < 0)
102✔
444
    {
445
        CPLError(CE_Failure, CPLE_AppDefined, "H5Dread() failed for block.");
×
446
        return CE_Failure;
×
447
    }
448

449
#ifdef HDF5_HAVE_FLOAT16
450
    if (eDataType == GDT_Float32 && poGDS->m_bConvertFromFloat16)
451
    {
452
        for (size_t i = static_cast<size_t>(nBlockXSize) * nBlockYSize; i > 0;
453
             /* do nothing */)
454
        {
455
            --i;
456
            uint16_t nVal16;
457
            memcpy(&nVal16, static_cast<uint16_t *>(pImage) + i,
458
                   sizeof(nVal16));
459
            const uint32_t nVal32 = CPLHalfToFloat(nVal16);
460
            float fVal;
461
            memcpy(&fVal, &nVal32, sizeof(fVal));
462
            *(static_cast<float *>(pImage) + i) = fVal;
463
        }
464
    }
465
    else if (eDataType == GDT_CFloat32 && poGDS->m_bConvertFromFloat16)
466
    {
467
        for (size_t i = static_cast<size_t>(nBlockXSize) * nBlockYSize; i > 0;
468
             /* do nothing */)
469
        {
470
            --i;
471
            for (int j = 1; j >= 0; --j)
472
            {
473
                uint16_t nVal16;
474
                memcpy(&nVal16, static_cast<uint16_t *>(pImage) + 2 * i + j,
475
                       sizeof(nVal16));
476
                const uint32_t nVal32 = CPLHalfToFloat(nVal16);
477
                float fVal;
478
                memcpy(&fVal, &nVal32, sizeof(fVal));
479
                *(static_cast<float *>(pImage) + 2 * i + j) = fVal;
480
            }
481
        }
482
    }
483
#endif
484

485
    return CE_None;
102✔
486
}
487

488
/************************************************************************/
489
/*                             IRasterIO()                              */
490
/************************************************************************/
491

492
CPLErr HDF5ImageRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
311✔
493
                                      int nXSize, int nYSize, void *pData,
494
                                      int nBufXSize, int nBufYSize,
495
                                      GDALDataType eBufType,
496
                                      GSpacing nPixelSpace, GSpacing nLineSpace,
497
                                      GDALRasterIOExtraArg *psExtraArg)
498

499
{
500
    HDF5ImageDataset *poGDS = static_cast<HDF5ImageDataset *>(poDS);
311✔
501

502
#ifdef HDF5_HAVE_FLOAT16
503
    if (poGDS->m_bConvertFromFloat16)
504
    {
505
        return GDALPamRasterBand::IRasterIO(
506
            eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
507
            eBufType, nPixelSpace, nLineSpace, psExtraArg);
508
    }
509
#endif
510

511
    const bool bIsBandInterleavedData =
512
        poGDS->ndims == 3 && poGDS->m_nOtherDimIndex == 0 &&
450✔
513
        poGDS->GetYIndex() == 1 && poGDS->GetXIndex() == 2;
761✔
514

515
    const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
311✔
516

517
    // Try to detect if we read whole bands by chunks of whole lines
518
    // If so, then read and cache whole band (or group of m_nBandChunkSize bands)
519
    // to save HDF5 decompression.
520
    if (m_nIRasterIORecCounter == 0)
311✔
521
    {
522
        bool bInvalidateWholeBandChunkOptim = false;
263✔
523
        if (!(nXSize == nBufXSize && nYSize == nBufYSize))
263✔
524
        {
525
            bInvalidateWholeBandChunkOptim = true;
1✔
526
        }
527
        // Is the first request on band 1, line 0 and one or several full lines?
528
        else if (poGDS->m_eWholeBandChunkOptim !=
262✔
529
                     HDF5ImageDataset::WBC_ENABLED &&
70✔
530
                 nBand == 1 && nXOff == 0 && nYOff == 0 &&
70✔
531
                 nXSize == nRasterXSize)
20✔
532
        {
533
            poGDS->m_eWholeBandChunkOptim =
19✔
534
                HDF5ImageDataset::WBC_DETECTION_IN_PROGRESS;
535
            poGDS->m_nLastRasterIOBand = 1;
19✔
536
            poGDS->m_nLastRasterIOXOff = nXOff;
19✔
537
            poGDS->m_nLastRasterIOYOff = nYOff;
19✔
538
            poGDS->m_nLastRasterIOXSize = nXSize;
19✔
539
            poGDS->m_nLastRasterIOYSize = nYSize;
19✔
540
        }
541
        else if (poGDS->m_eWholeBandChunkOptim ==
243✔
542
                 HDF5ImageDataset::WBC_DETECTION_IN_PROGRESS)
543
        {
544
            if (poGDS->m_nLastRasterIOBand == 1 && nBand == 1)
50✔
545
            {
546
                // Is this request a continuation of the previous one?
547
                if (nXOff == 0 && poGDS->m_nLastRasterIOXOff == 0 &&
44✔
548
                    nYOff == poGDS->m_nLastRasterIOYOff +
44✔
549
                                 poGDS->m_nLastRasterIOYSize &&
44✔
550
                    poGDS->m_nLastRasterIOXSize == nRasterXSize &&
43✔
551
                    nXSize == nRasterXSize)
43✔
552
                {
553
                    poGDS->m_nLastRasterIOXOff = nXOff;
43✔
554
                    poGDS->m_nLastRasterIOYOff = nYOff;
43✔
555
                    poGDS->m_nLastRasterIOXSize = nXSize;
43✔
556
                    poGDS->m_nLastRasterIOYSize = nYSize;
43✔
557
                }
558
                else
559
                {
560
                    bInvalidateWholeBandChunkOptim = true;
1✔
561
                }
562
            }
563
            else if (poGDS->m_nLastRasterIOBand == 1 && nBand == 2)
6✔
564
            {
565
                // Are we switching to band 2 while having fully read band 1?
566
                if (nXOff == 0 && nYOff == 0 && nXSize == nRasterXSize &&
5✔
567
                    poGDS->m_nLastRasterIOXOff == 0 &&
5✔
568
                    poGDS->m_nLastRasterIOXSize == nRasterXSize &&
5✔
569
                    poGDS->m_nLastRasterIOYOff + poGDS->m_nLastRasterIOYSize ==
5✔
570
                        nRasterYSize)
5✔
571
                {
572
                    if ((poGDS->m_nBandChunkSize > 1 ||
11✔
573
                         nBufYSize < nRasterYSize) &&
5✔
574
                        static_cast<int64_t>(poGDS->m_nBandChunkSize) *
1✔
575
                                nRasterXSize * nRasterYSize * nDTSize <
1✔
576
                            CPLGetUsablePhysicalRAM() / 10)
1✔
577
                    {
578
                        poGDS->m_eWholeBandChunkOptim =
1✔
579
                            HDF5ImageDataset::WBC_ENABLED;
580
                    }
581
                    else
582
                    {
583
                        bInvalidateWholeBandChunkOptim = true;
3✔
584
                    }
585
                }
586
                else
587
                {
588
                    bInvalidateWholeBandChunkOptim = true;
1✔
589
                }
590
            }
591
            else
592
            {
593
                bInvalidateWholeBandChunkOptim = true;
1✔
594
            }
595
        }
596
        if (bInvalidateWholeBandChunkOptim)
263✔
597
        {
598
            poGDS->m_eWholeBandChunkOptim = HDF5ImageDataset::WBC_DISABLED;
7✔
599
            poGDS->m_nLastRasterIOBand = -1;
7✔
600
            poGDS->m_nLastRasterIOXOff = -1;
7✔
601
            poGDS->m_nLastRasterIOYOff = -1;
7✔
602
            poGDS->m_nLastRasterIOXSize = -1;
7✔
603
            poGDS->m_nLastRasterIOYSize = -1;
7✔
604
        }
605
    }
606

607
    if (poGDS->m_eWholeBandChunkOptim == HDF5ImageDataset::WBC_ENABLED &&
311✔
608
        nXSize == nBufXSize && nYSize == nBufYSize)
193✔
609
    {
610
        if (poGDS->nBands == 1 || bIsBandInterleavedData)
193✔
611
        {
612
            if (poGDS->m_iCurrentBandChunk < 0)
193✔
613
                CPLDebug("HDF5", "Using whole band chunk caching");
1✔
614
            const int iBandChunk = (nBand - 1) / poGDS->m_nBandChunkSize;
193✔
615
            if (iBandChunk != poGDS->m_iCurrentBandChunk)
193✔
616
            {
617
                poGDS->m_abyBandChunk.resize(
15✔
618
                    static_cast<size_t>(poGDS->m_nBandChunkSize) *
15✔
619
                    nRasterXSize * nRasterYSize * nDTSize);
15✔
620

621
                HDF5_GLOBAL_LOCK();
622

623
                hsize_t count[3] = {
624
                    std::min(static_cast<hsize_t>(poGDS->nBands),
30✔
625
                             static_cast<hsize_t>(iBandChunk + 1) *
30✔
626
                                 poGDS->m_nBandChunkSize) -
15✔
627
                        static_cast<hsize_t>(iBandChunk) *
15✔
628
                            poGDS->m_nBandChunkSize,
15✔
629
                    static_cast<hsize_t>(nRasterYSize),
15✔
630
                    static_cast<hsize_t>(nRasterXSize)};
15✔
631
                H5OFFSET_TYPE offset[3] = {
15✔
632
                    static_cast<H5OFFSET_TYPE>(iBandChunk) *
15✔
633
                        poGDS->m_nBandChunkSize,
15✔
634
                    static_cast<H5OFFSET_TYPE>(0),
635
                    static_cast<H5OFFSET_TYPE>(0)};
15✔
636
                herr_t status =
637
                    H5Sselect_hyperslab(poGDS->dataspace_id, H5S_SELECT_SET,
15✔
638
                                        offset, nullptr, count, nullptr);
639
                if (status < 0)
15✔
640
                    return CE_Failure;
×
641

642
                const hid_t memspace =
643
                    H5Screate_simple(poGDS->ndims, count, nullptr);
15✔
644
                H5OFFSET_TYPE mem_offset[3] = {0, 0, 0};
15✔
645
                status =
646
                    H5Sselect_hyperslab(memspace, H5S_SELECT_SET, mem_offset,
15✔
647
                                        nullptr, count, nullptr);
648
                if (status < 0)
15✔
649
                {
650
                    H5Sclose(memspace);
×
651
                    return CE_Failure;
×
652
                }
653

654
                status = H5Dread(poGDS->dataset_id, poGDS->native, memspace,
15✔
655
                                 poGDS->dataspace_id, H5P_DEFAULT,
656
                                 poGDS->m_abyBandChunk.data());
15✔
657

658
                H5Sclose(memspace);
15✔
659

660
                if (status < 0)
15✔
661
                {
662
                    CPLError(
×
663
                        CE_Failure, CPLE_AppDefined,
664
                        "HDF5ImageRasterBand::IRasterIO(): H5Dread() failed");
665
                    return CE_Failure;
×
666
                }
667

668
                poGDS->m_iCurrentBandChunk = iBandChunk;
15✔
669
            }
670

671
            for (int iY = 0; iY < nYSize; ++iY)
1,447✔
672
            {
673
                GDALCopyWords(poGDS->m_abyBandChunk.data() +
2,508✔
674
                                  static_cast<size_t>((nBand - 1) %
1,254✔
675
                                                      poGDS->m_nBandChunkSize) *
1,254✔
676
                                      nRasterYSize * nRasterXSize * nDTSize +
1,254✔
677
                                  static_cast<size_t>(nYOff + iY) *
1,254✔
678
                                      nRasterXSize * nDTSize +
1,254✔
679
                                  nXOff * nDTSize,
1,254✔
680
                              eDataType, nDTSize,
681
                              static_cast<GByte *>(pData) +
682
                                  static_cast<size_t>(iY) * nLineSpace,
1,254✔
683
                              eBufType, static_cast<int>(nPixelSpace), nXSize);
684
            }
685
            return CE_None;
193✔
686
        }
687
    }
688

689
    const bool bIsExpectedLayout =
690
        (bIsBandInterleavedData ||
204✔
691
         (poGDS->ndims == 2 && poGDS->GetYIndex() == 0 &&
171✔
692
          poGDS->GetXIndex() == 1));
85✔
693
    if (eRWFlag == GF_Read && bIsExpectedLayout && nXSize == nBufXSize &&
118✔
694
        nYSize == nBufYSize && eBufType == eDataType &&
116✔
695
        nPixelSpace == nDTSize && nLineSpace == nXSize * nPixelSpace)
68✔
696
    {
697
        HDF5_GLOBAL_LOCK();
698

699
        hsize_t count[3] = {1, static_cast<hsize_t>(nYSize),
68✔
700
                            static_cast<hsize_t>(nXSize)};
68✔
701
        H5OFFSET_TYPE offset[3] = {static_cast<H5OFFSET_TYPE>(nBand - 1),
68✔
702
                                   static_cast<H5OFFSET_TYPE>(nYOff),
68✔
703
                                   static_cast<H5OFFSET_TYPE>(nXOff)};
68✔
704
        if (poGDS->ndims == 2)
68✔
705
        {
706
            count[0] = count[1];
47✔
707
            count[1] = count[2];
47✔
708

709
            offset[0] = offset[1];
47✔
710
            offset[1] = offset[2];
47✔
711
        }
712
        herr_t status = H5Sselect_hyperslab(poGDS->dataspace_id, H5S_SELECT_SET,
68✔
713
                                            offset, nullptr, count, nullptr);
714
        if (status < 0)
68✔
715
            return CE_Failure;
×
716

717
        const hid_t memspace = H5Screate_simple(poGDS->ndims, count, nullptr);
68✔
718
        H5OFFSET_TYPE mem_offset[3] = {0, 0, 0};
68✔
719
        status = H5Sselect_hyperslab(memspace, H5S_SELECT_SET, mem_offset,
68✔
720
                                     nullptr, count, nullptr);
721
        if (status < 0)
68✔
722
        {
723
            H5Sclose(memspace);
×
724
            return CE_Failure;
×
725
        }
726

727
        status = H5Dread(poGDS->dataset_id, poGDS->native, memspace,
68✔
728
                         poGDS->dataspace_id, H5P_DEFAULT, pData);
729

730
        H5Sclose(memspace);
68✔
731

732
        if (status < 0)
68✔
733
        {
734
            CPLError(CE_Failure, CPLE_AppDefined,
×
735
                     "HDF5ImageRasterBand::IRasterIO(): H5Dread() failed");
736
            return CE_Failure;
×
737
        }
738

739
        return CE_None;
68✔
740
    }
741

742
    // If the request is still small enough, try to read from libhdf5 with
743
    // the natural interleaving into a temporary MEMDataset, and then read
744
    // from it with the requested interleaving and data type.
745
    if (eRWFlag == GF_Read && bIsExpectedLayout && nXSize == nBufXSize &&
50✔
746
        nYSize == nBufYSize &&
100✔
747
        static_cast<GIntBig>(nXSize) * nYSize < CPLGetUsablePhysicalRAM() / 10)
48✔
748
    {
749
        auto poMemDS = std::unique_ptr<GDALDataset>(
750
            MEMDataset::Create("", nXSize, nYSize, 1, eDataType, nullptr));
48✔
751
        if (poMemDS)
48✔
752
        {
753
            void *pMemData = poMemDS->GetInternalHandle("MEMORY1");
48✔
754
            CPLAssert(pMemData);
48✔
755
            // Read from HDF5 into the temporary MEMDataset using the
756
            // natural interleaving of the HDF5 dataset
757
            ++m_nIRasterIORecCounter;
48✔
758
            CPLErr eErr =
759
                IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pMemData,
96✔
760
                          nXSize, nYSize, eDataType, nDTSize,
761
                          static_cast<GSpacing>(nXSize) * nDTSize, psExtraArg);
48✔
762
            --m_nIRasterIORecCounter;
48✔
763
            if (eErr != CE_None)
48✔
764
            {
765
                return CE_Failure;
×
766
            }
767
            // Copy to the final buffer using requested data type and spacings.
768
            return poMemDS->GetRasterBand(1)->RasterIO(
48✔
769
                GF_Read, 0, 0, nXSize, nYSize, pData, nXSize, nYSize, eBufType,
770
                nPixelSpace, nLineSpace, nullptr);
48✔
771
        }
772
    }
773

774
    return GDALPamRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
2✔
775
                                        pData, nBufXSize, nBufYSize, eBufType,
776
                                        nPixelSpace, nLineSpace, psExtraArg);
2✔
777
}
778

779
/************************************************************************/
780
/*                             IRasterIO()                              */
781
/************************************************************************/
782

783
CPLErr HDF5ImageDataset::IRasterIO(
90✔
784
    GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
785
    void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
786
    int nBandCount, BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
787
    GSpacing nLineSpace, GSpacing nBandSpace, GDALRasterIOExtraArg *psExtraArg)
788

789
{
790
#ifdef HDF5_HAVE_FLOAT16
791
    if (m_bConvertFromFloat16)
792
    {
793
        return HDF5Dataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
794
                                      pData, nBufXSize, nBufYSize, eBufType,
795
                                      nBandCount, panBandMap, nPixelSpace,
796
                                      nLineSpace, nBandSpace, psExtraArg);
797
    }
798
#endif
799

800
    const auto IsConsecutiveBands = [](const int *panVals, int nCount)
134✔
801
    {
802
        for (int i = 1; i < nCount; ++i)
151✔
803
        {
804
            if (panVals[i] != panVals[i - 1] + 1)
19✔
805
                return false;
2✔
806
        }
807
        return true;
132✔
808
    };
809

810
    const auto eDT = GetRasterBand(1)->GetRasterDataType();
90✔
811
    const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
90✔
812

813
    // Band-interleaved data and request
814
    const bool bIsBandInterleavedData = ndims == 3 && m_nOtherDimIndex == 0 &&
176✔
815
                                        GetYIndex() == 1 && GetXIndex() == 2;
266✔
816
    if (eRWFlag == GF_Read && bIsBandInterleavedData && nXSize == nBufXSize &&
90✔
817
        nYSize == nBufYSize && IsConsecutiveBands(panBandMap, nBandCount) &&
86✔
818
        eBufType == eDT && nPixelSpace == nDTSize &&
84✔
819
        nLineSpace == nXSize * nPixelSpace && nBandSpace == nYSize * nLineSpace)
180✔
820
    {
821
        HDF5_GLOBAL_LOCK();
822

823
        hsize_t count[3] = {static_cast<hsize_t>(nBandCount),
43✔
824
                            static_cast<hsize_t>(nYSize),
43✔
825
                            static_cast<hsize_t>(nXSize)};
43✔
826
        H5OFFSET_TYPE offset[3] = {
827
            static_cast<H5OFFSET_TYPE>(panBandMap[0] - 1),
43✔
828
            static_cast<H5OFFSET_TYPE>(nYOff),
43✔
829
            static_cast<H5OFFSET_TYPE>(nXOff)};
43✔
830
        herr_t status = H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET,
43✔
831
                                            offset, nullptr, count, nullptr);
832
        if (status < 0)
43✔
833
            return CE_Failure;
×
834

835
        const hid_t memspace = H5Screate_simple(ndims, count, nullptr);
43✔
836
        H5OFFSET_TYPE mem_offset[3] = {0, 0, 0};
43✔
837
        status = H5Sselect_hyperslab(memspace, H5S_SELECT_SET, mem_offset,
43✔
838
                                     nullptr, count, nullptr);
839
        if (status < 0)
43✔
840
        {
841
            H5Sclose(memspace);
×
842
            return CE_Failure;
×
843
        }
844

845
        status = H5Dread(dataset_id, native, memspace, dataspace_id,
43✔
846
                         H5P_DEFAULT, pData);
847

848
        H5Sclose(memspace);
43✔
849

850
        if (status < 0)
43✔
851
        {
852
            CPLError(CE_Failure, CPLE_AppDefined,
×
853
                     "HDF5ImageDataset::IRasterIO(): H5Dread() failed");
854
            return CE_Failure;
×
855
        }
856

857
        return CE_None;
43✔
858
    }
859

860
    // Pixel-interleaved data and request
861

862
    const bool bIsPixelInterleaveData = ndims == 3 && m_nOtherDimIndex == 2 &&
51✔
863
                                        GetYIndex() == 0 && GetXIndex() == 1;
98✔
864
    if (eRWFlag == GF_Read && bIsPixelInterleaveData && nXSize == nBufXSize &&
47✔
865
        nYSize == nBufYSize && IsConsecutiveBands(panBandMap, nBandCount) &&
4✔
866
        eBufType == eDT && nBandSpace == nDTSize &&
4✔
867
        nPixelSpace == nBandCount * nBandSpace &&
97✔
868
        nLineSpace == nXSize * nPixelSpace)
3✔
869
    {
870
        HDF5_GLOBAL_LOCK();
871

872
        hsize_t count[3] = {static_cast<hsize_t>(nYSize),
3✔
873
                            static_cast<hsize_t>(nXSize),
3✔
874
                            static_cast<hsize_t>(nBandCount)};
3✔
875
        H5OFFSET_TYPE offset[3] = {
876
            static_cast<H5OFFSET_TYPE>(nYOff),
3✔
877
            static_cast<H5OFFSET_TYPE>(nXOff),
3✔
878
            static_cast<H5OFFSET_TYPE>(panBandMap[0] - 1)};
3✔
879
        herr_t status = H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET,
3✔
880
                                            offset, nullptr, count, nullptr);
881
        if (status < 0)
3✔
882
            return CE_Failure;
×
883

884
        const hid_t memspace = H5Screate_simple(ndims, count, nullptr);
3✔
885
        H5OFFSET_TYPE mem_offset[3] = {0, 0, 0};
3✔
886
        status = H5Sselect_hyperslab(memspace, H5S_SELECT_SET, mem_offset,
3✔
887
                                     nullptr, count, nullptr);
888
        if (status < 0)
3✔
889
        {
890
            H5Sclose(memspace);
×
891
            return CE_Failure;
×
892
        }
893

894
        status = H5Dread(dataset_id, native, memspace, dataspace_id,
3✔
895
                         H5P_DEFAULT, pData);
896

897
        H5Sclose(memspace);
3✔
898

899
        if (status < 0)
3✔
900
        {
901
            CPLError(CE_Failure, CPLE_AppDefined,
×
902
                     "HDF5ImageDataset::IRasterIO(): H5Dread() failed");
903
            return CE_Failure;
×
904
        }
905

906
        return CE_None;
3✔
907
    }
908

909
    // If the request is still small enough, try to read from libhdf5 with
910
    // the natural interleaving into a temporary MEMDataset, and then read
911
    // from it with the requested interleaving and data type.
912
    if (eRWFlag == GF_Read &&
44✔
913
        (bIsBandInterleavedData || bIsPixelInterleaveData) &&
44✔
914
        nXSize == nBufXSize && nYSize == nBufYSize &&
88✔
915
        IsConsecutiveBands(panBandMap, nBandCount) &&
132✔
916
        static_cast<GIntBig>(nXSize) * nYSize <
43✔
917
            CPLGetUsablePhysicalRAM() / 10 / nBandCount)
43✔
918
    {
919
        const char *const apszOptions[] = {
43✔
920
            bIsPixelInterleaveData ? "INTERLEAVE=PIXEL" : nullptr, nullptr};
43✔
921
        auto poMemDS = std::unique_ptr<GDALDataset>(
922
            MEMDataset::Create("", nXSize, nYSize, nBandCount, eDT,
43✔
923
                               const_cast<char **>(apszOptions)));
43✔
924
        if (poMemDS)
43✔
925
        {
926
            void *pMemData = poMemDS->GetInternalHandle("MEMORY1");
43✔
927
            CPLAssert(pMemData);
43✔
928
            // Read from HDF5 into the temporary MEMDataset using the
929
            // natural interleaving of the HDF5 dataset
930
            if (IRasterIO(
129✔
931
                    eRWFlag, nXOff, nYOff, nXSize, nYSize, pMemData, nXSize,
932
                    nYSize, eDT, nBandCount, panBandMap,
933
                    bIsBandInterleavedData ? nDTSize : nDTSize * nBandCount,
1✔
934
                    bIsBandInterleavedData
935
                        ? static_cast<GSpacing>(nXSize) * nDTSize
42✔
936
                        : static_cast<GSpacing>(nXSize) * nDTSize * nBandCount,
1✔
937
                    bIsBandInterleavedData
938
                        ? static_cast<GSpacing>(nYSize) * nXSize * nDTSize
42✔
939
                        : nDTSize,
940
                    psExtraArg) != CE_None)
43✔
941
            {
942
                return CE_Failure;
×
943
            }
944
            // Copy to the final buffer using requested data type and spacings.
945
            return poMemDS->RasterIO(GF_Read, 0, 0, nXSize, nYSize, pData,
43✔
946
                                     nXSize, nYSize, eBufType, nBandCount,
947
                                     nullptr, nPixelSpace, nLineSpace,
948
                                     nBandSpace, nullptr);
43✔
949
        }
950
    }
951

952
    return HDF5Dataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
1✔
953
                                  nBufXSize, nBufYSize, eBufType, nBandCount,
954
                                  panBandMap, nPixelSpace, nLineSpace,
955
                                  nBandSpace, psExtraArg);
1✔
956
}
957

958
/************************************************************************/
959
/*                                Open()                                */
960
/************************************************************************/
961
GDALDataset *HDF5ImageDataset::Open(GDALOpenInfo *poOpenInfo)
58✔
962
{
963
    if (!STARTS_WITH_CI(poOpenInfo->pszFilename, "HDF5:"))
58✔
964
        return nullptr;
×
965

966
    HDF5_GLOBAL_LOCK();
967

968
    // Confirm the requested access is supported.
969
    if (poOpenInfo->eAccess == GA_Update)
58✔
970
    {
971
        ReportUpdateNotSupportedByDriver("HDF5");
×
972
        return nullptr;
×
973
    }
974

975
    HDF5ImageDataset *poDS = new HDF5ImageDataset();
58✔
976

977
    // Create a corresponding GDALDataset.
978
    char **papszName =
979
        CSLTokenizeString2(poOpenInfo->pszFilename, ":",
58✔
980
                           CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES);
981

982
    if (!(CSLCount(papszName) == 3 || CSLCount(papszName) == 4))
58✔
983
    {
984
        CSLDestroy(papszName);
×
985
        delete poDS;
×
986
        return nullptr;
×
987
    }
988

989
    poDS->SetDescription(poOpenInfo->pszFilename);
58✔
990

991
    // Check for drive name in windows HDF5:"D:\...
992
    CPLString osSubdatasetName;
116✔
993

994
    CPLString osFilename(papszName[1]);
116✔
995

996
    if ((strlen(papszName[1]) == 1 && papszName[3] != nullptr) ||
58✔
997
        (STARTS_WITH(papszName[1], "/vsicurl/http") && papszName[3] != nullptr))
58✔
998
    {
999
        osFilename += ":";
×
1000
        osFilename += papszName[2];
×
1001
        osSubdatasetName = papszName[3];
×
1002
    }
1003
    else
1004
    {
1005
        osSubdatasetName = papszName[2];
58✔
1006
    }
1007

1008
    poDS->SetSubdatasetName(osSubdatasetName);
58✔
1009

1010
    CSLDestroy(papszName);
58✔
1011
    papszName = nullptr;
58✔
1012

1013
    poDS->SetPhysicalFilename(osFilename);
58✔
1014

1015
    // Try opening the dataset.
1016
    poDS->m_hHDF5 = GDAL_HDF5Open(osFilename);
58✔
1017
    if (poDS->m_hHDF5 < 0)
58✔
1018
    {
1019
        delete poDS;
×
1020
        return nullptr;
×
1021
    }
1022

1023
    poDS->hGroupID = H5Gopen(poDS->m_hHDF5, "/");
58✔
1024
    if (poDS->hGroupID < 0)
58✔
1025
    {
1026
        delete poDS;
×
1027
        return nullptr;
×
1028
    }
1029

1030
    // This is an HDF5 file.
1031
    poDS->ReadGlobalAttributes(FALSE);
58✔
1032

1033
    // Create HDF5 Data Hierarchy in a link list.
1034
    poDS->poH5Objects = poDS->HDF5FindDatasetObjectsbyPath(poDS->poH5RootGroup,
58✔
1035
                                                           osSubdatasetName);
1036

1037
    if (poDS->poH5Objects == nullptr)
58✔
1038
    {
1039
        delete poDS;
×
1040
        return nullptr;
×
1041
    }
1042

1043
    // Retrieve HDF5 data information.
1044
    poDS->dataset_id = H5Dopen(poDS->m_hHDF5, poDS->poH5Objects->pszPath);
58✔
1045
    poDS->dataspace_id = H5Dget_space(poDS->dataset_id);
58✔
1046
    poDS->ndims = H5Sget_simple_extent_ndims(poDS->dataspace_id);
58✔
1047
    if (poDS->ndims <= 0)
58✔
1048
    {
1049
        delete poDS;
×
1050
        return nullptr;
×
1051
    }
1052
    poDS->dims =
58✔
1053
        static_cast<hsize_t *>(CPLCalloc(poDS->ndims, sizeof(hsize_t)));
58✔
1054
    poDS->maxdims =
58✔
1055
        static_cast<hsize_t *>(CPLCalloc(poDS->ndims, sizeof(hsize_t)));
58✔
1056
    poDS->dimensions = H5Sget_simple_extent_dims(poDS->dataspace_id, poDS->dims,
58✔
1057
                                                 poDS->maxdims);
1058
    for (int i = 0; i < poDS->dimensions; ++i)
185✔
1059
    {
1060
        if (poDS->dims[i] >
256✔
1061
            static_cast<hsize_t>(std::numeric_limits<int>::max()))
128✔
1062
        {
1063
            CPLError(CE_Failure, CPLE_NotSupported,
1✔
1064
                     "At least one dimension size exceeds INT_MAX !");
1065
            delete poDS;
1✔
1066
            return nullptr;
1✔
1067
        }
1068
    }
1069

1070
    auto datatype = H5Dget_type(poDS->dataset_id);
57✔
1071
    poDS->native = H5Tget_native_type(datatype, H5T_DIR_ASCEND);
57✔
1072
    H5Tclose(datatype);
57✔
1073

1074
    const auto eGDALDataType = poDS->GetDataType(poDS->native);
57✔
1075
    if (eGDALDataType == GDT_Unknown)
57✔
1076
    {
1077
        CPLError(CE_Failure, CPLE_AppDefined, "Unhandled HDF5 data type");
×
1078
        delete poDS;
×
1079
        return nullptr;
×
1080
    }
1081

1082
#ifdef HDF5_HAVE_FLOAT16
1083
    if (H5Tequal(H5T_NATIVE_FLOAT16, poDS->native) ||
1084
        IsNativeCFloat16(poDS->native))
1085
    {
1086
        poDS->m_bConvertFromFloat16 = true;
1087
    }
1088
#endif
1089

1090
    // CSK code in IdentifyProductType() and CreateProjections()
1091
    // uses dataset metadata.
1092
    poDS->SetMetadata(poDS->m_aosMetadata.List());
57✔
1093

1094
    // Check if the hdf5 is a well known product type
1095
    poDS->IdentifyProductType();
57✔
1096

1097
    poDS->m_nYIndex = poDS->IsComplexCSKL1A() ? 0 : poDS->ndims - 2;
57✔
1098
    poDS->m_nXIndex = poDS->IsComplexCSKL1A() ? 1 : poDS->ndims - 1;
57✔
1099

1100
    if (poDS->IsComplexCSKL1A())
57✔
1101
    {
1102
        poDS->m_nOtherDimIndex = 2;
×
1103
    }
1104
    else if (poDS->ndims == 3)
57✔
1105
    {
1106
        poDS->m_nOtherDimIndex = 0;
14✔
1107
    }
1108

1109
    if (HDF5EOSParser::HasHDFEOS(poDS->hGroupID))
57✔
1110
    {
1111
        HDF5EOSParser oHDFEOSParser;
34✔
1112
        if (oHDFEOSParser.Parse(poDS->hGroupID))
17✔
1113
        {
1114
            CPLDebug("HDF5", "Successfully parsed HDFEOS metadata");
17✔
1115
            HDF5EOSParser::GridDataFieldMetadata oGridDataFieldMetadata;
34✔
1116
            HDF5EOSParser::SwathDataFieldMetadata oSwathDataFieldMetadata;
34✔
1117
            if (oHDFEOSParser.GetDataModel() ==
17✔
1118
                    HDF5EOSParser::DataModel::GRID &&
5✔
1119
                oHDFEOSParser.GetGridDataFieldMetadata(
5✔
1120
                    osSubdatasetName.c_str(), oGridDataFieldMetadata) &&
22✔
1121
                static_cast<int>(oGridDataFieldMetadata.aoDimensions.size()) ==
5✔
1122
                    poDS->ndims)
5✔
1123
            {
1124
                int iDim = 0;
5✔
1125
                for (const auto &oDim : oGridDataFieldMetadata.aoDimensions)
17✔
1126
                {
1127
                    if (oDim.osName == "XDim")
12✔
1128
                        poDS->m_nXIndex = iDim;
5✔
1129
                    else if (oDim.osName == "YDim")
7✔
1130
                        poDS->m_nYIndex = iDim;
5✔
1131
                    else
1132
                        poDS->m_nOtherDimIndex = iDim;
2✔
1133
                    ++iDim;
12✔
1134
                }
1135

1136
                if (oGridDataFieldMetadata.poGridMetadata->GetGeoTransform(
5✔
1137
                        poDS->m_gt))
5✔
1138
                    poDS->bHasGeoTransform = true;
5✔
1139

1140
                auto poSRS = oGridDataFieldMetadata.poGridMetadata->GetSRS();
10✔
1141
                if (poSRS)
5✔
1142
                    poDS->m_oSRS = *(poSRS.get());
5✔
1143
            }
1144
            else if (oHDFEOSParser.GetDataModel() ==
12✔
1145
                         HDF5EOSParser::DataModel::SWATH &&
12✔
1146
                     oHDFEOSParser.GetSwathDataFieldMetadata(
12✔
1147
                         osSubdatasetName.c_str(), oSwathDataFieldMetadata) &&
6✔
1148
                     static_cast<int>(
1149
                         oSwathDataFieldMetadata.aoDimensions.size()) ==
6✔
1150
                         poDS->ndims &&
6✔
1151
                     oSwathDataFieldMetadata.iXDim >= 0 &&
30✔
1152
                     oSwathDataFieldMetadata.iYDim >= 0)
6✔
1153
            {
1154
                poDS->m_nXIndex = oSwathDataFieldMetadata.iXDim;
6✔
1155
                poDS->m_nYIndex = oSwathDataFieldMetadata.iYDim;
6✔
1156
                poDS->m_nOtherDimIndex = oSwathDataFieldMetadata.iOtherDim;
6✔
1157
                if (!oSwathDataFieldMetadata.osLongitudeSubdataset.empty())
6✔
1158
                {
1159
                    // Arbitrary
1160
                    poDS->SetMetadataItem("SRS", SRS_WKT_WGS84_LAT_LONG,
6✔
1161
                                          "GEOLOCATION");
6✔
1162
                    poDS->SetMetadataItem(
6✔
1163
                        "X_DATASET",
1164
                        ("HDF5:\"" + osFilename +
12✔
1165
                         "\":" + oSwathDataFieldMetadata.osLongitudeSubdataset)
12✔
1166
                            .c_str(),
1167
                        "GEOLOCATION");
6✔
1168
                    poDS->SetMetadataItem("X_BAND", "1", "GEOLOCATION");
6✔
1169
                    poDS->SetMetadataItem(
6✔
1170
                        "Y_DATASET",
1171
                        ("HDF5:\"" + osFilename +
12✔
1172
                         "\":" + oSwathDataFieldMetadata.osLatitudeSubdataset)
12✔
1173
                            .c_str(),
1174
                        "GEOLOCATION");
6✔
1175
                    poDS->SetMetadataItem("Y_BAND", "1", "GEOLOCATION");
6✔
1176
                    poDS->SetMetadataItem(
6✔
1177
                        "PIXEL_OFFSET",
1178
                        CPLSPrintf("%d", oSwathDataFieldMetadata.nPixelOffset),
1179
                        "GEOLOCATION");
6✔
1180
                    poDS->SetMetadataItem(
6✔
1181
                        "PIXEL_STEP",
1182
                        CPLSPrintf("%d", oSwathDataFieldMetadata.nPixelStep),
1183
                        "GEOLOCATION");
6✔
1184
                    poDS->SetMetadataItem(
6✔
1185
                        "LINE_OFFSET",
1186
                        CPLSPrintf("%d", oSwathDataFieldMetadata.nLineOffset),
1187
                        "GEOLOCATION");
6✔
1188
                    poDS->SetMetadataItem(
6✔
1189
                        "LINE_STEP",
1190
                        CPLSPrintf("%d", oSwathDataFieldMetadata.nLineStep),
1191
                        "GEOLOCATION");
6✔
1192
                    // Not totally sure about that
1193
                    poDS->SetMetadataItem("GEOREFERENCING_CONVENTION",
6✔
1194
                                          "PIXEL_CENTER", "GEOLOCATION");
6✔
1195
                }
1196
            }
1197
        }
1198
    }
1199

1200
    poDS->nRasterYSize =
57✔
1201
        poDS->GetYIndex() < 0
57✔
1202
            ? 1
57✔
1203
            : static_cast<int>(poDS->dims[poDS->GetYIndex()]);  // nRows
56✔
1204
    poDS->nRasterXSize =
57✔
1205
        static_cast<int>(poDS->dims[poDS->GetXIndex()]);  // nCols
57✔
1206
    int nBands = 1;
57✔
1207
    if (poDS->m_nOtherDimIndex >= 0)
57✔
1208
    {
1209
        nBands = static_cast<int>(poDS->dims[poDS->m_nOtherDimIndex]);
14✔
1210
    }
1211

1212
    CPLStringList aosMetadata;
114✔
1213
    std::map<std::string, CPLStringList> oMapBandSpecificMetadata;
57✔
1214
    if (poDS->poH5Objects->nType == H5G_DATASET)
57✔
1215
    {
1216
        HDF5Dataset::CreateMetadata(poDS->m_hHDF5, poDS->poH5Objects,
57✔
1217
                                    H5G_DATASET, false, aosMetadata);
1218
        if (nBands > 1 && poDS->nRasterXSize != nBands &&
57✔
1219
            poDS->nRasterYSize != nBands)
13✔
1220
        {
1221
            // Heuristics to detect non-scalar attributes, that are intended
1222
            // to be attached to a specific band.
1223
            const CPLStringList aosMetadataDup(aosMetadata);
26✔
1224
            for (const auto &[pszKey, pszValue] :
52✔
1225
                 cpl::IterateNameValue(aosMetadataDup))
65✔
1226
            {
1227
                const hid_t hAttrID = H5Aopen_name(poDS->dataset_id, pszKey);
26✔
1228
                const hid_t hAttrSpace = H5Aget_space(hAttrID);
26✔
1229
                if (H5Sget_simple_extent_ndims(hAttrSpace) == 1 &&
47✔
1230
                    H5Sget_simple_extent_npoints(hAttrSpace) == nBands)
21✔
1231
                {
1232
                    CPLStringList aosTokens(
1233
                        CSLTokenizeString2(pszValue, " ", 0));
40✔
1234
                    if (aosTokens.size() == nBands)
20✔
1235
                    {
1236
                        std::string osAttrName(pszKey);
40✔
1237
                        if (osAttrName.size() > strlen("_coefficients") &&
30✔
1238
                            osAttrName.substr(osAttrName.size() -
30✔
1239
                                              strlen("_coefficients")) ==
10✔
1240
                                "_coefficients")
1241
                        {
1242
                            osAttrName.pop_back();
5✔
1243
                        }
1244
                        else if (osAttrName.size() > strlen("_wavelengths") &&
25✔
1245
                                 osAttrName.substr(osAttrName.size() -
25✔
1246
                                                   strlen("_wavelengths")) ==
10✔
1247
                                     "_wavelengths")
1248
                        {
1249
                            osAttrName.pop_back();
5✔
1250
                        }
1251
                        else if (osAttrName.size() > strlen("_list") &&
15✔
1252
                                 osAttrName.substr(osAttrName.size() -
15✔
1253
                                                   strlen("_list")) == "_list")
5✔
1254
                        {
1255
                            osAttrName.resize(osAttrName.size() -
5✔
1256
                                              strlen("_list"));
1257
                        }
1258
                        oMapBandSpecificMetadata[osAttrName] =
20✔
1259
                            std::move(aosTokens);
40✔
1260
                        aosMetadata.SetNameValue(pszKey, nullptr);
20✔
1261
                    }
1262
                }
1263
                H5Sclose(hAttrSpace);
26✔
1264
                H5Aclose(hAttrID);
26✔
1265
            }
1266
        }
1267
    }
1268

1269
    poDS->m_nBlockXSize = poDS->GetRasterXSize();
57✔
1270
    poDS->m_nBlockYSize = 1;
57✔
1271
    poDS->m_nBandChunkSize = 1;
57✔
1272

1273
    // Check for chunksize and set it as the blocksize (optimizes read).
1274
    const hid_t listid = H5Dget_create_plist(poDS->dataset_id);
57✔
1275
    if (listid > 0)
57✔
1276
    {
1277
        if (H5Pget_layout(listid) == H5D_CHUNKED)
57✔
1278
        {
1279
            hsize_t panChunkDims[3] = {0, 0, 0};
12✔
1280
            const int nDimSize = H5Pget_chunk(listid, 3, panChunkDims);
12✔
1281
            CPL_IGNORE_RET_VAL(nDimSize);
12✔
1282
            CPLAssert(nDimSize == poDS->ndims);
12✔
1283
            poDS->m_nBlockXSize =
12✔
1284
                static_cast<int>(panChunkDims[poDS->GetXIndex()]);
12✔
1285
            if (poDS->GetYIndex() >= 0)
12✔
1286
                poDS->m_nBlockYSize =
12✔
1287
                    static_cast<int>(panChunkDims[poDS->GetYIndex()]);
12✔
1288
            if (nBands > 1)
12✔
1289
            {
1290
                poDS->m_nBandChunkSize =
3✔
1291
                    static_cast<int>(panChunkDims[poDS->m_nOtherDimIndex]);
3✔
1292

1293
                poDS->SetMetadataItem("BAND_CHUNK_SIZE",
3✔
1294
                                      CPLSPrintf("%d", poDS->m_nBandChunkSize),
1295
                                      "IMAGE_STRUCTURE");
3✔
1296
            }
1297
        }
1298

1299
        const int nFilters = H5Pget_nfilters(listid);
57✔
1300
        for (int i = 0; i < nFilters; ++i)
63✔
1301
        {
1302
            unsigned int flags = 0;
6✔
1303
            size_t cd_nelmts = 0;
6✔
1304
            char szName[64 + 1] = {0};
6✔
1305
            const auto eFilter = H5Pget_filter(listid, i, &flags, &cd_nelmts,
6✔
1306
                                               nullptr, 64, szName);
1307
            if (eFilter == H5Z_FILTER_DEFLATE)
6✔
1308
            {
1309
                poDS->SetMetadataItem("COMPRESSION", "DEFLATE",
5✔
1310
                                      "IMAGE_STRUCTURE");
5✔
1311
            }
1312
            else if (eFilter == H5Z_FILTER_SZIP)
1✔
1313
            {
1314
                poDS->SetMetadataItem("COMPRESSION", "SZIP", "IMAGE_STRUCTURE");
×
1315
            }
1316
        }
1317

1318
        H5Pclose(listid);
57✔
1319
    }
1320

1321
    for (int i = 0; i < nBands; i++)
180✔
1322
    {
1323
        HDF5ImageRasterBand *const poBand =
1324
            new HDF5ImageRasterBand(poDS, i + 1, eGDALDataType);
123✔
1325

1326
        poDS->SetBand(i + 1, poBand);
123✔
1327

1328
        if (poDS->poH5Objects->nType == H5G_DATASET)
123✔
1329
        {
1330
            poBand->SetMetadata(aosMetadata.List());
123✔
1331
            for (const auto &oIter : oMapBandSpecificMetadata)
163✔
1332
            {
1333
                poBand->SetMetadataItem(oIter.first.c_str(), oIter.second[i]);
40✔
1334
            }
1335
        }
1336
    }
1337

1338
    if (!poDS->GetMetadata("GEOLOCATION"))
57✔
1339
        poDS->CreateProjections();
51✔
1340

1341
    // Setup/check for pam .aux.xml.
1342
    poDS->TryLoadXML();
57✔
1343

1344
    // Setup overviews.
1345
    poDS->oOvManager.Initialize(poDS, ":::VIRTUAL:::");
57✔
1346

1347
    return poDS;
57✔
1348
}
1349

1350
/************************************************************************/
1351
/*                   HDF5ImageDatasetDriverUnload()                     */
1352
/************************************************************************/
1353

1354
static void HDF5ImageDatasetDriverUnload(GDALDriver *)
6✔
1355
{
1356
    HDF5UnloadFileDriver();
6✔
1357
}
6✔
1358

1359
/************************************************************************/
1360
/*                        GDALRegister_HDF5Image()                      */
1361
/************************************************************************/
1362
void GDALRegister_HDF5Image()
11✔
1363

1364
{
1365
    if (!GDAL_CHECK_VERSION("HDF5Image driver"))
11✔
1366
        return;
×
1367

1368
    if (GDALGetDriverByName(HDF5_IMAGE_DRIVER_NAME) != nullptr)
11✔
1369
        return;
×
1370

1371
    GDALDriver *poDriver = new GDALDriver();
11✔
1372

1373
    HDF5ImageDriverSetCommonMetadata(poDriver);
11✔
1374

1375
    poDriver->pfnOpen = HDF5ImageDataset::Open;
11✔
1376
    poDriver->pfnUnloadDriver = HDF5ImageDatasetDriverUnload;
11✔
1377

1378
    GetGDALDriverManager()->RegisterDriver(poDriver);
11✔
1379
}
1380

1381
/************************************************************************/
1382
/*                       CreateODIMH5Projection()                       */
1383
/************************************************************************/
1384

1385
// Reference:
1386
//   http://www.knmi.nl/opera/opera3/OPERA_2008_03_WP2.1b_ODIM_H5_v2.1.pdf
1387
//
1388
// 4.3.2 where for geographically referenced image Groups
1389
// We don't use the where_xscale and where_yscale parameters, but recompute them
1390
// from the lower-left and upper-right coordinates. There's some difference.
1391
// As all those parameters are linked together, I'm not sure which one should be
1392
// considered as the reference.
1393

1394
CPLErr HDF5ImageDataset::CreateODIMH5Projection()
×
1395
{
1396
    const char *const pszProj4String = GetMetadataItem("where_projdef");
×
1397
    const char *const pszLL_lon = GetMetadataItem("where_LL_lon");
×
1398
    const char *const pszLL_lat = GetMetadataItem("where_LL_lat");
×
1399
    const char *const pszUR_lon = GetMetadataItem("where_UR_lon");
×
1400
    const char *const pszUR_lat = GetMetadataItem("where_UR_lat");
×
1401
    if (pszProj4String == nullptr || pszLL_lon == nullptr ||
×
1402
        pszLL_lat == nullptr || pszUR_lon == nullptr || pszUR_lat == nullptr)
×
1403
        return CE_Failure;
×
1404

1405
    m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
×
1406
    if (m_oSRS.importFromProj4(pszProj4String) != OGRERR_NONE)
×
1407
        return CE_Failure;
×
1408

1409
    OGRSpatialReference oSRSWGS84;
×
1410
    oSRSWGS84.SetWellKnownGeogCS("WGS84");
×
1411
    oSRSWGS84.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
×
1412

1413
    OGRCoordinateTransformation *poCT =
1414
        OGRCreateCoordinateTransformation(&oSRSWGS84, &m_oSRS);
×
1415
    if (poCT == nullptr)
×
1416
        return CE_Failure;
×
1417

1418
    // Reproject corners from long,lat WGS84 to the target SRS.
1419
    double dfLLX = CPLAtof(pszLL_lon);
×
1420
    double dfLLY = CPLAtof(pszLL_lat);
×
1421
    double dfURX = CPLAtof(pszUR_lon);
×
1422
    double dfURY = CPLAtof(pszUR_lat);
×
1423
    if (!poCT->Transform(1, &dfLLX, &dfLLY) ||
×
1424
        !poCT->Transform(1, &dfURX, &dfURY))
×
1425
    {
1426
        delete poCT;
×
1427
        return CE_Failure;
×
1428
    }
1429
    delete poCT;
×
1430

1431
    // Compute the geotransform now.
1432
    const double dfPixelX = (dfURX - dfLLX) / nRasterXSize;
×
1433
    const double dfPixelY = (dfURY - dfLLY) / nRasterYSize;
×
1434

1435
    bHasGeoTransform = true;
×
1436
    m_gt[0] = dfLLX;
×
1437
    m_gt[1] = dfPixelX;
×
1438
    m_gt[2] = 0;
×
1439
    m_gt[3] = dfURY;
×
1440
    m_gt[4] = 0;
×
1441
    m_gt[5] = -dfPixelY;
×
1442

1443
    return CE_None;
×
1444
}
1445

1446
/************************************************************************/
1447
/*                         CreateProjections()                          */
1448
/************************************************************************/
1449
CPLErr HDF5ImageDataset::CreateProjections()
51✔
1450
{
1451
    switch (iSubdatasetType)
51✔
1452
    {
1453
        case CSK_PRODUCT:
2✔
1454
        {
1455
            int productType = PROD_UNKNOWN;
2✔
1456

1457
            if (GetMetadataItem("Product_Type") != nullptr)
2✔
1458
            {
1459
                // Get the format's level.
1460
                const char *osMissionLevel =
1461
                    HDF5Dataset::GetMetadataItem("Product_Type");
2✔
1462

1463
                if (STARTS_WITH_CI(osMissionLevel, "RAW"))
2✔
1464
                    productType = PROD_CSK_L0;
×
1465

1466
                if (STARTS_WITH_CI(osMissionLevel, "SSC"))
2✔
1467
                    productType = PROD_CSK_L1A;
×
1468

1469
                if (STARTS_WITH_CI(osMissionLevel, "DGM"))
2✔
1470
                    productType = PROD_CSK_L1B;
1✔
1471

1472
                if (STARTS_WITH_CI(osMissionLevel, "GEC"))
2✔
1473
                    productType = PROD_CSK_L1C;
1✔
1474

1475
                if (STARTS_WITH_CI(osMissionLevel, "GTC"))
2✔
1476
                    productType = PROD_CSK_L1D;
×
1477
            }
1478

1479
            CaptureCSKGeoTransform(productType);
2✔
1480
            CaptureCSKGeolocation(productType);
2✔
1481
            CaptureCSKGCPs(productType);
2✔
1482

1483
            break;
2✔
1484
        }
1485
        case UNKNOWN_PRODUCT:
49✔
1486
        {
1487
            constexpr int NBGCPLAT = 100;
49✔
1488
            constexpr int NBGCPLON = 30;
49✔
1489

1490
            const int nDeltaLat = nRasterYSize / NBGCPLAT;
49✔
1491
            const int nDeltaLon = nRasterXSize / NBGCPLON;
49✔
1492

1493
            if (nDeltaLat == 0 || nDeltaLon == 0)
49✔
1494
                return CE_None;
41✔
1495

1496
            // Create HDF5 Data Hierarchy in a link list.
1497
            poH5Objects = HDF5FindDatasetObjects(poH5RootGroup, "Latitude");
8✔
1498
            if (!poH5Objects)
8✔
1499
            {
1500
                if (GetMetadataItem("where_projdef") != nullptr)
8✔
1501
                    return CreateODIMH5Projection();
×
1502
                return CE_None;
8✔
1503
            }
1504

1505
            // The Latitude and Longitude arrays must have a rank of 2 to
1506
            // retrieve GCPs.
1507
            if (poH5Objects->nRank != 2 ||
×
1508
                poH5Objects->paDims[0] != static_cast<size_t>(nRasterYSize) ||
×
1509
                poH5Objects->paDims[1] != static_cast<size_t>(nRasterXSize))
×
1510
            {
1511
                return CE_None;
×
1512
            }
1513

1514
            // Retrieve HDF5 data information.
1515
            const hid_t LatitudeDatasetID =
1516
                H5Dopen(m_hHDF5, poH5Objects->pszPath);
×
1517
            // LatitudeDataspaceID = H5Dget_space(dataset_id);
1518

1519
            poH5Objects = HDF5FindDatasetObjects(poH5RootGroup, "Longitude");
×
1520
            // GCPs.
1521
            if (poH5Objects == nullptr || poH5Objects->nRank != 2 ||
×
1522
                poH5Objects->paDims[0] != static_cast<size_t>(nRasterYSize) ||
×
1523
                poH5Objects->paDims[1] != static_cast<size_t>(nRasterXSize))
×
1524
            {
1525
                if (LatitudeDatasetID > 0)
×
1526
                    H5Dclose(LatitudeDatasetID);
×
1527
                return CE_None;
×
1528
            }
1529

1530
            const hid_t LongitudeDatasetID =
1531
                H5Dopen(m_hHDF5, poH5Objects->pszPath);
×
1532
            // LongitudeDataspaceID = H5Dget_space(dataset_id);
1533

1534
            if (LatitudeDatasetID > 0 && LongitudeDatasetID > 0)
×
1535
            {
1536
                float *const Latitude =
1537
                    static_cast<float *>(VSI_MALLOC3_VERBOSE(
×
1538
                        nRasterYSize, nRasterXSize, sizeof(float)));
1539
                float *const Longitude =
1540
                    static_cast<float *>(VSI_MALLOC3_VERBOSE(
×
1541
                        nRasterYSize, nRasterXSize, sizeof(float)));
1542
                if (!Latitude || !Longitude)
×
1543
                {
1544
                    CPLFree(Latitude);
×
1545
                    CPLFree(Longitude);
×
1546
                    H5Dclose(LatitudeDatasetID);
×
1547
                    H5Dclose(LongitudeDatasetID);
×
1548
                    return CE_Failure;
×
1549
                }
1550
                memset(Latitude, 0,
×
1551
                       static_cast<size_t>(nRasterXSize) * nRasterYSize *
×
1552
                           sizeof(float));
1553
                memset(Longitude, 0,
×
1554
                       static_cast<size_t>(nRasterXSize) * nRasterYSize *
×
1555
                           sizeof(float));
1556

1557
                // netCDF convention for nodata
1558
                double dfLatNoData = 0;
×
1559
                bool bHasLatNoData = GH5_FetchAttribute(
×
1560
                    LatitudeDatasetID, "_FillValue", dfLatNoData);
1561

1562
                double dfLongNoData = 0;
×
1563
                bool bHasLongNoData = GH5_FetchAttribute(
×
1564
                    LongitudeDatasetID, "_FillValue", dfLongNoData);
1565

1566
                H5Dread(LatitudeDatasetID, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
×
1567
                        H5P_DEFAULT, Latitude);
1568

1569
                H5Dread(LongitudeDatasetID, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
×
1570
                        H5P_DEFAULT, Longitude);
1571

1572
                m_oSRS.Clear();
×
1573
                m_oGCPSRS.SetWellKnownGeogCS("WGS84");
×
1574

1575
                const int nYLimit =
×
1576
                    (static_cast<int>(nRasterYSize) / nDeltaLat) * nDeltaLat;
×
1577
                const int nXLimit =
×
1578
                    (static_cast<int>(nRasterXSize) / nDeltaLon) * nDeltaLon;
×
1579

1580
                // The original code in
1581
                // https://trac.osgeo.org/gdal/changeset/8066 always add +180 to
1582
                // the longitudes, but without justification I suspect this
1583
                // might be due to handling products crossing the antimeridian.
1584
                // Trying to do it just when needed through a heuristics.
1585
                bool bHasLonNearMinus180 = false;
×
1586
                bool bHasLonNearPlus180 = false;
×
1587
                bool bHasLonNearZero = false;
×
1588
                for (int j = 0; j < nYLimit; j += nDeltaLat)
×
1589
                {
1590
                    for (int i = 0; i < nXLimit; i += nDeltaLon)
×
1591
                    {
1592
                        const int iGCP = j * nRasterXSize + i;
×
1593
                        if ((bHasLatNoData && static_cast<float>(dfLatNoData) ==
×
1594
                                                  Latitude[iGCP]) ||
×
1595
                            (bHasLongNoData &&
×
1596
                             static_cast<float>(dfLongNoData) ==
×
1597
                                 Longitude[iGCP]))
×
1598
                            continue;
×
1599
                        if (Longitude[iGCP] > 170 && Longitude[iGCP] <= 180)
×
1600
                            bHasLonNearPlus180 = true;
×
1601
                        if (Longitude[iGCP] < -170 && Longitude[iGCP] >= -180)
×
1602
                            bHasLonNearMinus180 = true;
×
1603
                        if (fabs(Longitude[iGCP]) < 90)
×
1604
                            bHasLonNearZero = true;
×
1605
                    }
1606
                }
1607

1608
                // Fill the GCPs list.
1609
                const char *pszShiftGCP =
1610
                    CPLGetConfigOption("HDF5_SHIFT_GCPX_BY_180", nullptr);
×
1611
                const bool bAdd180 =
1612
                    (bHasLonNearPlus180 && bHasLonNearMinus180 &&
×
1613
                     !bHasLonNearZero && pszShiftGCP == nullptr) ||
×
1614
                    (pszShiftGCP != nullptr && CPLTestBool(pszShiftGCP));
×
1615

1616
                for (int j = 0; j < nYLimit; j += nDeltaLat)
×
1617
                {
1618
                    for (int i = 0; i < nXLimit; i += nDeltaLon)
×
1619
                    {
1620
                        const int iGCP = j * nRasterXSize + i;
×
1621
                        if ((bHasLatNoData && static_cast<float>(dfLatNoData) ==
×
1622
                                                  Latitude[iGCP]) ||
×
1623
                            (bHasLongNoData &&
×
1624
                             static_cast<float>(dfLongNoData) ==
×
1625
                                 Longitude[iGCP]))
×
1626
                            continue;
×
1627
                        double dfGCPX = static_cast<double>(Longitude[iGCP]);
×
1628
                        if (bAdd180)
×
1629
                            dfGCPX += 180.0;
×
1630
                        const double dfGCPY =
×
1631
                            static_cast<double>(Latitude[iGCP]);
×
1632

1633
                        m_aoGCPs.emplace_back("", "", i + 0.5, j + 0.5, dfGCPX,
×
1634
                                              dfGCPY);
×
1635
                    }
1636
                }
1637

1638
                CPLFree(Latitude);
×
1639
                CPLFree(Longitude);
×
1640
            }
1641

1642
            if (LatitudeDatasetID > 0)
×
1643
                H5Dclose(LatitudeDatasetID);
×
1644
            if (LongitudeDatasetID > 0)
×
1645
                H5Dclose(LongitudeDatasetID);
×
1646

1647
            break;
×
1648
        }
1649
    }
1650

1651
    return CE_None;
2✔
1652
}
1653

1654
/************************************************************************/
1655
/*                         GetMetadataItem()                            */
1656
/************************************************************************/
1657

1658
const char *HDF5ImageDataset::GetMetadataItem(const char *pszName,
120✔
1659
                                              const char *pszDomain)
1660
{
1661
    if (pszDomain && EQUAL(pszDomain, "__DEBUG__") &&
120✔
1662
        EQUAL(pszName, "WholeBandChunkOptim"))
66✔
1663
    {
1664
        switch (m_eWholeBandChunkOptim)
66✔
1665
        {
1666
            case WBC_DETECTION_IN_PROGRESS:
4✔
1667
                return "DETECTION_IN_PROGRESS";
4✔
1668
            case WBC_DISABLED:
3✔
1669
                return "DISABLED";
3✔
1670
            case WBC_ENABLED:
59✔
1671
                return "ENABLED";
59✔
1672
        }
1673
    }
1674
    return GDALPamDataset::GetMetadataItem(pszName, pszDomain);
54✔
1675
}
1676

1677
/************************************************************************/
1678
/*                         GetSpatialRef()                              */
1679
/************************************************************************/
1680

1681
const OGRSpatialReference *HDF5ImageDataset::GetSpatialRef() const
12✔
1682
{
1683
    if (!m_oSRS.IsEmpty())
12✔
1684
        return &m_oSRS;
11✔
1685
    return GDALPamDataset::GetSpatialRef();
1✔
1686
}
1687

1688
/************************************************************************/
1689
/*                            GetGCPCount()                             */
1690
/************************************************************************/
1691

1692
int HDF5ImageDataset::GetGCPCount()
3✔
1693

1694
{
1695
    if (!m_aoGCPs.empty())
3✔
1696
        return static_cast<int>(m_aoGCPs.size());
1✔
1697

1698
    return GDALPamDataset::GetGCPCount();
2✔
1699
}
1700

1701
/************************************************************************/
1702
/*                        GetGCPSpatialRef()                            */
1703
/************************************************************************/
1704

1705
const OGRSpatialReference *HDF5ImageDataset::GetGCPSpatialRef() const
1✔
1706

1707
{
1708
    if (!m_aoGCPs.empty() && !m_oGCPSRS.IsEmpty())
1✔
1709
        return &m_oGCPSRS;
1✔
1710

1711
    return GDALPamDataset::GetGCPSpatialRef();
×
1712
}
1713

1714
/************************************************************************/
1715
/*                               GetGCPs()                              */
1716
/************************************************************************/
1717

1718
const GDAL_GCP *HDF5ImageDataset::GetGCPs()
2✔
1719
{
1720
    if (!m_aoGCPs.empty())
2✔
1721
        return gdal::GCP::c_ptr(m_aoGCPs);
1✔
1722

1723
    return GDALPamDataset::GetGCPs();
1✔
1724
}
1725

1726
/************************************************************************/
1727
/*                         GetGeoTransform()                            */
1728
/************************************************************************/
1729

1730
CPLErr HDF5ImageDataset::GetGeoTransform(GDALGeoTransform &gt) const
6✔
1731
{
1732
    if (bHasGeoTransform)
6✔
1733
    {
1734
        gt = m_gt;
5✔
1735
        return CE_None;
5✔
1736
    }
1737

1738
    return GDALPamDataset::GetGeoTransform(gt);
1✔
1739
}
1740

1741
/************************************************************************/
1742
/*                       IdentifyProductType()                          */
1743
/************************************************************************/
1744

1745
/**
1746
 * Identify if the subdataset has a known product format
1747
 * It stores a product identifier in iSubdatasetType,
1748
 * UNKNOWN_PRODUCT, if it isn't a recognizable format.
1749
 */
1750
void HDF5ImageDataset::IdentifyProductType()
57✔
1751
{
1752
    iSubdatasetType = UNKNOWN_PRODUCT;
57✔
1753

1754
    // COSMO-SKYMED
1755

1756
    // Get the Mission Id as a char *, because the field may not exist.
1757
    const char *const pszMissionId = HDF5Dataset::GetMetadataItem("Mission_ID");
57✔
1758

1759
    // If there is a Mission_ID field.
1760
    if (pszMissionId != nullptr && strstr(GetDescription(), "QLK") == nullptr)
57✔
1761
    {
1762
        // Check if the mission type is CSK, KMPS or CSG.
1763
        // KMPS: Komsat-5 is Korean mission with a SAR instrument.
1764
        // CSG: Cosmo Skymed 2nd Generation
1765
        if (EQUAL(pszMissionId, "CSK") || EQUAL(pszMissionId, "KMPS") ||
2✔
1766
            EQUAL(pszMissionId, "CSG"))
×
1767
        {
1768
            iSubdatasetType = CSK_PRODUCT;
2✔
1769

1770
            if (GetMetadataItem("Product_Type") != nullptr)
2✔
1771
            {
1772
                // Get the format's level.
1773
                const char *osMissionLevel =
1774
                    HDF5Dataset::GetMetadataItem("Product_Type");
2✔
1775

1776
                if (STARTS_WITH_CI(osMissionLevel, "RAW"))
2✔
1777
                    iCSKProductType = PROD_CSK_L0;
×
1778

1779
                if (STARTS_WITH_CI(osMissionLevel, "SCS"))
2✔
1780
                    iCSKProductType = PROD_CSK_L1A;
×
1781

1782
                if (STARTS_WITH_CI(osMissionLevel, "DGM"))
2✔
1783
                    iCSKProductType = PROD_CSK_L1B;
1✔
1784

1785
                if (STARTS_WITH_CI(osMissionLevel, "GEC"))
2✔
1786
                    iCSKProductType = PROD_CSK_L1C;
1✔
1787

1788
                if (STARTS_WITH_CI(osMissionLevel, "GTC"))
2✔
1789
                    iCSKProductType = PROD_CSK_L1D;
×
1790
            }
1791
        }
1792
    }
1793
}
57✔
1794

1795
/************************************************************************/
1796
/*                       CaptureCSKGeolocation()                        */
1797
/************************************************************************/
1798

1799
/**
1800
 * Captures Geolocation information from a COSMO-SKYMED
1801
 * file.
1802
 * The geoid will always be WGS84
1803
 * The projection type may be UTM or UPS, depending on the
1804
 * latitude from the center of the image.
1805
 * @param iProductType type of CSK subproduct, see HDF5CSKProduct
1806
 */
1807
void HDF5ImageDataset::CaptureCSKGeolocation(int iProductType)
2✔
1808
{
1809
    // Set the ellipsoid to WGS84.
1810
    m_oSRS.SetWellKnownGeogCS("WGS84");
2✔
1811

1812
    if (iProductType == PROD_CSK_L1C || iProductType == PROD_CSK_L1D)
2✔
1813
    {
1814
        double *dfProjFalseEastNorth = nullptr;
1✔
1815
        double *dfProjScaleFactor = nullptr;
1✔
1816
        double *dfCenterCoord = nullptr;
1✔
1817

1818
        // Check if all the metadata attributes are present.
1819
        if (HDF5ReadDoubleAttr("Map Projection False East-North",
1✔
1820
                               &dfProjFalseEastNorth) == CE_Failure ||
1✔
1821
            HDF5ReadDoubleAttr("Map Projection Scale Factor",
1✔
1822
                               &dfProjScaleFactor) == CE_Failure ||
1✔
1823
            HDF5ReadDoubleAttr("Map Projection Centre", &dfCenterCoord) ==
1✔
1824
                CE_Failure ||
2✔
1825
            GetMetadataItem("Projection_ID") == nullptr)
1✔
1826
        {
1827
            m_oSRS.Clear();
×
1828
            m_oGCPSRS.Clear();
×
1829
            CPLError(CE_Failure, CPLE_OpenFailed,
×
1830
                     "The CSK hdf5 file geolocation information is "
1831
                     "malformed");
1832
        }
1833
        else
1834
        {
1835
            // Fetch projection Type.
1836
            CPLString osProjectionID = GetMetadataItem("Projection_ID");
2✔
1837

1838
            // If the projection is UTM.
1839
            if (EQUAL(osProjectionID, "UTM"))
1✔
1840
            {
1841
                // @TODO: use SetUTM
1842
                m_oSRS.SetProjCS(SRS_PT_TRANSVERSE_MERCATOR);
1✔
1843
                m_oSRS.SetTM(dfCenterCoord[0], dfCenterCoord[1],
1✔
1844
                             dfProjScaleFactor[0], dfProjFalseEastNorth[0],
1845
                             dfProjFalseEastNorth[1]);
1✔
1846
            }
1847
            else
1848
            {
1849
                // TODO: No UPS projected files to test!
1850
                // If the projection is UPS.
1851
                if (EQUAL(osProjectionID, "UPS"))
×
1852
                {
1853
                    m_oSRS.SetProjCS(SRS_PT_POLAR_STEREOGRAPHIC);
×
1854
                    m_oSRS.SetPS(dfCenterCoord[0], dfCenterCoord[1],
×
1855
                                 dfProjScaleFactor[0], dfProjFalseEastNorth[0],
1856
                                 dfProjFalseEastNorth[1]);
×
1857
                }
1858
            }
1859

1860
            CPLFree(dfCenterCoord);
1✔
1861
            CPLFree(dfProjScaleFactor);
1✔
1862
            CPLFree(dfProjFalseEastNorth);
1✔
1863
        }
1✔
1864
    }
1865
    else
1866
    {
1867
        m_oGCPSRS = m_oSRS;
1✔
1868
    }
1869
}
2✔
1870

1871
/************************************************************************/
1872
/*                       CaptureCSKGeoTransform()                       */
1873
/************************************************************************/
1874

1875
/**
1876
 * Get Geotransform information for COSMO-SKYMED files
1877
 * In case of success it stores the transformation
1878
 * in m_gt. In case of failure it doesn't
1879
 * modify m_gt
1880
 * @param iProductType type of CSK subproduct, see HDF5CSKProduct
1881
 */
1882
void HDF5ImageDataset::CaptureCSKGeoTransform(int iProductType)
2✔
1883
{
1884
    const char *pszSubdatasetName = GetSubdatasetName();
2✔
1885

1886
    bHasGeoTransform = false;
2✔
1887
    // If the product level is not L1C or L1D then
1888
    // it doesn't have a valid projection.
1889
    if (iProductType == PROD_CSK_L1C || iProductType == PROD_CSK_L1D)
2✔
1890
    {
1891
        // If there is a subdataset.
1892
        if (pszSubdatasetName != nullptr)
1✔
1893
        {
1894
            CPLString osULPath = pszSubdatasetName;
2✔
1895
            osULPath += "/Top Left East-North";
1✔
1896

1897
            CPLString osLineSpacingPath = pszSubdatasetName;
2✔
1898
            osLineSpacingPath += "/Line Spacing";
1✔
1899

1900
            CPLString osColumnSpacingPath = pszSubdatasetName;
2✔
1901
            osColumnSpacingPath += "/Column Spacing";
1✔
1902

1903
            double *pdOutUL = nullptr;
1✔
1904
            double *pdLineSpacing = nullptr;
1✔
1905
            double *pdColumnSpacing = nullptr;
1✔
1906

1907
            // If it could find the attributes on the metadata.
1908
            if (HDF5ReadDoubleAttr(osULPath.c_str(), &pdOutUL) == CE_Failure ||
1✔
1909
                HDF5ReadDoubleAttr(osLineSpacingPath.c_str(), &pdLineSpacing) ==
1✔
1910
                    CE_Failure ||
2✔
1911
                HDF5ReadDoubleAttr(osColumnSpacingPath.c_str(),
1✔
1912
                                   &pdColumnSpacing) == CE_Failure)
1913
            {
1914
                bHasGeoTransform = false;
×
1915
            }
1916
            else
1917
            {
1918
                // geotransform[1] : width of pixel
1919
                // geotransform[4] : rotational coefficient, zero for north up
1920
                // images.
1921
                // geotransform[2] : rotational coefficient, zero for north up
1922
                // images.
1923
                // geotransform[5] : height of pixel (but negative)
1924

1925
                m_gt[0] = pdOutUL[0];
1✔
1926
                m_gt[1] = pdLineSpacing[0];
1✔
1927
                m_gt[2] = 0;
1✔
1928
                m_gt[3] = pdOutUL[1];
1✔
1929
                m_gt[4] = 0;
1✔
1930
                m_gt[5] = -pdColumnSpacing[0];
1✔
1931

1932
                CPLFree(pdOutUL);
1✔
1933
                CPLFree(pdLineSpacing);
1✔
1934
                CPLFree(pdColumnSpacing);
1✔
1935

1936
                bHasGeoTransform = true;
1✔
1937
            }
1938
        }
1939
    }
1940
}
2✔
1941

1942
/************************************************************************/
1943
/*                          CaptureCSKGCPs()                            */
1944
/************************************************************************/
1945

1946
/**
1947
 * Retrieves and stores the GCPs from a COSMO-SKYMED dataset
1948
 * It only retrieves the GCPs for L0, L1A and L1B products
1949
 * for L1C and L1D products, geotransform is provided.
1950
 * The GCPs provided will be the Image's corners.
1951
 * @param iProductType type of CSK product @see HDF5CSKProductEnum
1952
 */
1953
void HDF5ImageDataset::CaptureCSKGCPs(int iProductType)
2✔
1954
{
1955
    // Only retrieve GCPs for L0,L1A and L1B products.
1956
    if (iProductType == PROD_CSK_L0 || iProductType == PROD_CSK_L1A ||
2✔
1957
        iProductType == PROD_CSK_L1B)
1958
    {
1959
        CPLString osCornerName[4];
10✔
1960
        double pdCornerPixel[4] = {0.0, 0.0, 0.0, 0.0};
1✔
1961
        double pdCornerLine[4] = {0.0, 0.0, 0.0, 0.0};
1✔
1962

1963
        const char *const pszSubdatasetName = GetSubdatasetName();
1✔
1964

1965
        // Load the subdataset name first.
1966
        for (int i = 0; i < 4; i++)
5✔
1967
            osCornerName[i] = pszSubdatasetName;
4✔
1968

1969
        // Load the attribute name, and raster coordinates for
1970
        // all the corners.
1971
        osCornerName[0] += "/Top Left Geodetic Coordinates";
1✔
1972
        pdCornerPixel[0] = 0;
1✔
1973
        pdCornerLine[0] = 0;
1✔
1974

1975
        osCornerName[1] += "/Top Right Geodetic Coordinates";
1✔
1976
        pdCornerPixel[1] = GetRasterXSize();
1✔
1977
        pdCornerLine[1] = 0;
1✔
1978

1979
        osCornerName[2] += "/Bottom Left Geodetic Coordinates";
1✔
1980
        pdCornerPixel[2] = 0;
1✔
1981
        pdCornerLine[2] = GetRasterYSize();
1✔
1982

1983
        osCornerName[3] += "/Bottom Right Geodetic Coordinates";
1✔
1984
        pdCornerPixel[3] = GetRasterXSize();
1✔
1985
        pdCornerLine[3] = GetRasterYSize();
1✔
1986

1987
        // For all the image's corners.
1988
        for (int i = 0; i < 4; i++)
5✔
1989
        {
1990
            double *pdCornerCoordinates = nullptr;
4✔
1991

1992
            // Retrieve the attributes.
1993
            if (HDF5ReadDoubleAttr(osCornerName[i].c_str(),
4✔
1994
                                   &pdCornerCoordinates) == CE_Failure)
4✔
1995
            {
1996
                CPLError(CE_Failure, CPLE_OpenFailed,
×
1997
                         "Error retrieving CSK GCPs");
1998
                m_aoGCPs.clear();
×
1999
                break;
×
2000
            }
2001

2002
            m_aoGCPs.emplace_back(osCornerName[i].c_str(), "", pdCornerPixel[i],
×
2003
                                  pdCornerLine[i],
4✔
2004
                                  /* X = */ pdCornerCoordinates[1],
4✔
2005
                                  /* Y = */ pdCornerCoordinates[0],
2006
                                  /* Z = */ pdCornerCoordinates[2]);
4✔
2007

2008
            // Free the returned coordinates.
2009
            CPLFree(pdCornerCoordinates);
4✔
2010
        }
2011
    }
2012
}
2✔
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2025 Coveralls, Inc