Line data Source code
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 : #ifdef _POSIX_C_SOURCE
15 : #undef _POSIX_C_SOURCE
16 : #endif
17 :
18 : #include "hdf5_api.h"
19 :
20 : #include "cpl_float.h"
21 : #include "cpl_string.h"
22 : #include "gdal_frmts.h"
23 : #include "gdal_pam.h"
24 : #include "gdal_priv.h"
25 : #include "gh5_convenience.h"
26 : #include "hdf5dataset.h"
27 : #include "hdf5drivercore.h"
28 : #include "ogr_spatialref.h"
29 : #include "memdataset.h"
30 :
31 : #include <algorithm>
32 : #include <cmath>
33 : #include <cstring>
34 : #include <limits>
35 :
36 : class HDF5ImageDataset final : public HDF5Dataset
37 : {
38 : typedef enum
39 : {
40 : UNKNOWN_PRODUCT = 0,
41 : CSK_PRODUCT
42 : } Hdf5ProductType;
43 :
44 : typedef enum
45 : {
46 : PROD_UNKNOWN = 0,
47 : PROD_CSK_L0,
48 : PROD_CSK_L1A,
49 : PROD_CSK_L1B,
50 : PROD_CSK_L1C,
51 : PROD_CSK_L1D
52 : } HDF5CSKProductEnum;
53 :
54 : friend class HDF5ImageRasterBand;
55 :
56 : OGRSpatialReference m_oSRS{};
57 : OGRSpatialReference m_oGCPSRS{};
58 : std::vector<gdal::GCP> m_aoGCPs{};
59 :
60 : hsize_t *dims;
61 : hsize_t *maxdims;
62 : int ndims;
63 : std::string m_osPath{};
64 : hid_t dataset_id;
65 : hid_t dataspace_id;
66 : hid_t native;
67 : #ifdef HDF5_HAVE_FLOAT16
68 : bool m_bConvertFromFloat16 = false;
69 : #endif
70 : Hdf5ProductType iSubdatasetType;
71 : HDF5CSKProductEnum iCSKProductType;
72 : GDALGeoTransform m_gt{};
73 : bool bHasGeoTransform;
74 : int m_nXIndex = -1;
75 : int m_nYIndex = -1;
76 : int m_nOtherDimIndex = -1;
77 :
78 : int m_nBlockXSize = 0;
79 : int m_nBlockYSize = 0;
80 : int m_nBandChunkSize = 1; //! Number of bands in a chunk
81 :
82 : enum WholeBandChunkOptim
83 : {
84 : WBC_DETECTION_IN_PROGRESS,
85 : WBC_DISABLED,
86 : WBC_ENABLED,
87 : };
88 :
89 : //! Flag to detect if the read pattern of HDF5ImageRasterBand::IRasterIO()
90 : // is whole band after whole band.
91 : WholeBandChunkOptim m_eWholeBandChunkOptim = WBC_DETECTION_IN_PROGRESS;
92 : //! Value of nBand during last HDF5ImageRasterBand::IRasterIO() call
93 : int m_nLastRasterIOBand = -1;
94 : //! Value of nXOff during last HDF5ImageRasterBand::IRasterIO() call
95 : int m_nLastRasterIOXOff = -1;
96 : //! Value of nYOff during last HDF5ImageRasterBand::IRasterIO() call
97 : int m_nLastRasterIOYOff = -1;
98 : //! Value of nXSize during last HDF5ImageRasterBand::IRasterIO() call
99 : int m_nLastRasterIOXSize = -1;
100 : //! Value of nYSize during last HDF5ImageRasterBand::IRasterIO() call
101 : int m_nLastRasterIOYSize = -1;
102 : //! Value such that m_abyBandChunk represent band data in the range
103 : // [m_iCurrentBandChunk * m_nBandChunkSize, (m_iCurrentBandChunk+1) * m_nBandChunkSize[
104 : int m_iCurrentBandChunk = -1;
105 : //! Cached values (in native data type) for bands in the range
106 : // [m_iCurrentBandChunk * m_nBandChunkSize, (m_iCurrentBandChunk+1) * m_nBandChunkSize[
107 : std::vector<GByte> m_abyBandChunk{};
108 :
109 : CPLErr CreateODIMH5Projection();
110 :
111 : CPL_DISALLOW_COPY_ASSIGN(HDF5ImageDataset)
112 :
113 : public:
114 : HDF5ImageDataset();
115 : ~HDF5ImageDataset() override;
116 :
117 : CPLErr CreateProjections();
118 : static GDALDataset *Open(GDALOpenInfo *);
119 : static int Identify(GDALOpenInfo *);
120 :
121 : const OGRSpatialReference *GetSpatialRef() const override;
122 : int GetGCPCount() override;
123 : const OGRSpatialReference *GetGCPSpatialRef() const override;
124 : const GDAL_GCP *GetGCPs() override;
125 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
126 :
127 : CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,
128 : int nYSize, void *pData, int nBufXSize, int nBufYSize,
129 : GDALDataType eBufType, int nBandCount,
130 : BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
131 : GSpacing nLineSpace, GSpacing nBandSpace,
132 : GDALRasterIOExtraArg *psExtraArg) override;
133 :
134 : const char *GetMetadataItem(const char *pszName,
135 : const char *pszDomain = "") override;
136 :
137 186 : Hdf5ProductType GetSubdatasetType() const
138 : {
139 186 : return iSubdatasetType;
140 : }
141 :
142 6 : HDF5CSKProductEnum GetCSKProductType() const
143 : {
144 6 : return iCSKProductType;
145 : }
146 :
147 186 : int IsComplexCSKL1A() const
148 : {
149 192 : return GetSubdatasetType() == CSK_PRODUCT &&
150 192 : GetCSKProductType() == PROD_CSK_L1A && ndims == 3;
151 : }
152 :
153 572 : int GetYIndex() const
154 : {
155 572 : return m_nYIndex;
156 : }
157 :
158 703 : int GetXIndex() const
159 : {
160 703 : return m_nXIndex;
161 : }
162 :
163 : /**
164 : * Identify if the subdataset has a known product format
165 : * It stores a product identifier in iSubdatasetType,
166 : * UNKNOWN_PRODUCT, if it isn't a recognizable format.
167 : */
168 : void IdentifyProductType();
169 :
170 : /**
171 : * Captures Geolocation information from a COSMO-SKYMED
172 : * file.
173 : * The geoid will always be WGS84
174 : * The projection type may be UTM or UPS, depending on the
175 : * latitude from the center of the image.
176 : * @param iProductType type of HDF5 subproduct, see HDF5CSKProduct
177 : */
178 : void CaptureCSKGeolocation(int iProductType);
179 :
180 : /**
181 : * Get Geotransform information for COSMO-SKYMED files
182 : * In case of success it stores the transformation
183 : * in m_gt. In case of failure it doesn't
184 : * modify m_gt
185 : * @param iProductType type of HDF5 subproduct, see HDF5CSKProduct
186 : */
187 : void CaptureCSKGeoTransform(int iProductType);
188 :
189 : /**
190 : * @param iProductType type of HDF5 subproduct, see HDF5CSKProduct
191 : */
192 : void CaptureCSKGCPs(int iProductType);
193 : };
194 :
195 : /************************************************************************/
196 : /* ==================================================================== */
197 : /* HDF5ImageDataset */
198 : /* ==================================================================== */
199 : /************************************************************************/
200 :
201 : /************************************************************************/
202 : /* HDF5ImageDataset() */
203 : /************************************************************************/
204 63 : HDF5ImageDataset::HDF5ImageDataset()
205 : : dims(nullptr), maxdims(nullptr), ndims(0), dataset_id(-1),
206 : dataspace_id(-1), native(-1), iSubdatasetType(UNKNOWN_PRODUCT),
207 63 : iCSKProductType(PROD_UNKNOWN), bHasGeoTransform(false)
208 : {
209 63 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
210 63 : m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
211 63 : }
212 :
213 : /************************************************************************/
214 : /* ~HDF5ImageDataset() */
215 : /************************************************************************/
216 126 : HDF5ImageDataset::~HDF5ImageDataset()
217 : {
218 : HDF5_GLOBAL_LOCK();
219 :
220 63 : FlushCache(true);
221 :
222 63 : if (dataset_id > 0)
223 63 : H5Dclose(dataset_id);
224 63 : if (dataspace_id > 0)
225 63 : H5Sclose(dataspace_id);
226 63 : if (native > 0)
227 62 : H5Tclose(native);
228 :
229 63 : CPLFree(dims);
230 63 : CPLFree(maxdims);
231 126 : }
232 :
233 : /************************************************************************/
234 : /* ==================================================================== */
235 : /* Hdf5imagerasterband */
236 : /* ==================================================================== */
237 : /************************************************************************/
238 : class HDF5ImageRasterBand final : public GDALPamRasterBand
239 : {
240 : friend class HDF5ImageDataset;
241 :
242 : bool bNoDataSet = false;
243 : double dfNoDataValue = -9999.0;
244 : bool m_bHasOffset = false;
245 : double m_dfOffset = 0.0;
246 : bool m_bHasScale = false;
247 : double m_dfScale = 1.0;
248 : int m_nIRasterIORecCounter = 0;
249 :
250 : public:
251 : HDF5ImageRasterBand(HDF5ImageDataset *, int, GDALDataType);
252 : ~HDF5ImageRasterBand() override;
253 :
254 : CPLErr IReadBlock(int, int, void *) override;
255 : double GetNoDataValue(int *) override;
256 : double GetOffset(int *) override;
257 : double GetScale(int *) override;
258 : // virtual CPLErr IWriteBlock( int, int, void * );
259 :
260 : CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,
261 : int nYSize, void *pData, int nBufXSize, int nBufYSize,
262 : GDALDataType eBufType, GSpacing nPixelSpace,
263 : GSpacing nLineSpace,
264 : GDALRasterIOExtraArg *psExtraArg) override;
265 : };
266 :
267 : /************************************************************************/
268 : /* ~HDF5ImageRasterBand() */
269 : /************************************************************************/
270 :
271 258 : HDF5ImageRasterBand::~HDF5ImageRasterBand()
272 : {
273 258 : }
274 :
275 : /************************************************************************/
276 : /* HDF5ImageRasterBand() */
277 : /************************************************************************/
278 129 : HDF5ImageRasterBand::HDF5ImageRasterBand(HDF5ImageDataset *poDSIn, int nBandIn,
279 129 : GDALDataType eType)
280 : {
281 129 : poDS = poDSIn;
282 129 : nBand = nBandIn;
283 129 : eDataType = eType;
284 129 : nBlockXSize = poDSIn->m_nBlockXSize;
285 129 : nBlockYSize = poDSIn->m_nBlockYSize;
286 :
287 : // netCDF convention for nodata
288 129 : bNoDataSet =
289 129 : GH5_FetchAttribute(poDSIn->dataset_id, "_FillValue", dfNoDataValue);
290 129 : if (!bNoDataSet)
291 120 : dfNoDataValue = -9999.0;
292 :
293 : // netCDF conventions for scale and offset
294 129 : m_bHasOffset =
295 129 : GH5_FetchAttribute(poDSIn->dataset_id, "add_offset", m_dfOffset);
296 129 : if (!m_bHasOffset)
297 128 : m_dfOffset = 0.0;
298 129 : m_bHasScale =
299 129 : GH5_FetchAttribute(poDSIn->dataset_id, "scale_factor", m_dfScale);
300 129 : if (!m_bHasScale)
301 128 : m_dfScale = 1.0;
302 129 : }
303 :
304 : /************************************************************************/
305 : /* GetNoDataValue() */
306 : /************************************************************************/
307 49 : double HDF5ImageRasterBand::GetNoDataValue(int *pbSuccess)
308 :
309 : {
310 49 : if (bNoDataSet)
311 : {
312 3 : if (pbSuccess)
313 3 : *pbSuccess = bNoDataSet;
314 :
315 3 : return dfNoDataValue;
316 : }
317 :
318 46 : return GDALPamRasterBand::GetNoDataValue(pbSuccess);
319 : }
320 :
321 : /************************************************************************/
322 : /* GetOffset() */
323 : /************************************************************************/
324 :
325 22 : double HDF5ImageRasterBand::GetOffset(int *pbSuccess)
326 :
327 : {
328 22 : if (m_bHasOffset)
329 : {
330 1 : if (pbSuccess)
331 1 : *pbSuccess = m_bHasOffset;
332 :
333 1 : return m_dfOffset;
334 : }
335 :
336 21 : return GDALPamRasterBand::GetOffset(pbSuccess);
337 : }
338 :
339 : /************************************************************************/
340 : /* GetScale() */
341 : /************************************************************************/
342 :
343 22 : double HDF5ImageRasterBand::GetScale(int *pbSuccess)
344 :
345 : {
346 22 : if (m_bHasScale)
347 : {
348 1 : if (pbSuccess)
349 1 : *pbSuccess = m_bHasScale;
350 :
351 1 : return m_dfScale;
352 : }
353 :
354 21 : return GDALPamRasterBand::GetScale(pbSuccess);
355 : }
356 :
357 : /************************************************************************/
358 : /* IReadBlock() */
359 : /************************************************************************/
360 103 : CPLErr HDF5ImageRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
361 : void *pImage)
362 : {
363 103 : HDF5ImageDataset *poGDS = static_cast<HDF5ImageDataset *>(poDS);
364 :
365 103 : memset(pImage, 0,
366 103 : static_cast<size_t>(nBlockXSize) * nBlockYSize *
367 103 : GDALGetDataTypeSizeBytes(eDataType));
368 :
369 103 : if (poGDS->eAccess == GA_Update)
370 : {
371 0 : return CE_None;
372 : }
373 :
374 103 : const int nXOff = nBlockXOff * nBlockXSize;
375 103 : const int nYOff = nBlockYOff * nBlockYSize;
376 103 : const int nXSize = std::min(nBlockXSize, nRasterXSize - nXOff);
377 103 : const int nYSize = std::min(nBlockYSize, nRasterYSize - nYOff);
378 103 : if (poGDS->m_eWholeBandChunkOptim == HDF5ImageDataset::WBC_ENABLED)
379 : {
380 : const bool bIsBandInterleavedData =
381 2 : poGDS->ndims == 3 && poGDS->m_nOtherDimIndex == 0 &&
382 3 : poGDS->GetYIndex() == 1 && poGDS->GetXIndex() == 2;
383 1 : if (poGDS->nBands == 1 || bIsBandInterleavedData)
384 : {
385 : GDALRasterIOExtraArg sExtraArg;
386 1 : INIT_RASTERIO_EXTRA_ARG(sExtraArg);
387 1 : const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
388 2 : return IRasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize, pImage,
389 : nXSize, nYSize, eDataType, nDTSize,
390 1 : static_cast<GSpacing>(nDTSize) * nBlockXSize,
391 1 : &sExtraArg);
392 : }
393 : }
394 :
395 : HDF5_GLOBAL_LOCK();
396 :
397 102 : hsize_t count[3] = {0, 0, 0};
398 102 : H5OFFSET_TYPE offset[3] = {0, 0, 0};
399 102 : hsize_t col_dims[3] = {0, 0, 0};
400 102 : hsize_t rank = std::min(poGDS->ndims, 2);
401 :
402 102 : if (poGDS->ndims == 3)
403 : {
404 1 : rank = 3;
405 1 : offset[poGDS->m_nOtherDimIndex] = nBand - 1;
406 1 : count[poGDS->m_nOtherDimIndex] = 1;
407 1 : col_dims[poGDS->m_nOtherDimIndex] = 1;
408 : }
409 :
410 102 : const int nYIndex = poGDS->GetYIndex();
411 : // Blocksize may not be a multiple of imagesize.
412 102 : if (nYIndex >= 0)
413 : {
414 101 : offset[nYIndex] = nYOff;
415 101 : count[nYIndex] = nYSize;
416 : }
417 102 : offset[poGDS->GetXIndex()] = nXOff;
418 102 : count[poGDS->GetXIndex()] = nXSize;
419 :
420 : // Select block from file space.
421 102 : herr_t status = H5Sselect_hyperslab(poGDS->dataspace_id, H5S_SELECT_SET,
422 : offset, nullptr, count, nullptr);
423 102 : if (status < 0)
424 0 : return CE_Failure;
425 :
426 : // Create memory space to receive the data.
427 102 : if (nYIndex >= 0)
428 101 : col_dims[nYIndex] = nBlockYSize;
429 102 : col_dims[poGDS->GetXIndex()] = nBlockXSize;
430 :
431 : const hid_t memspace =
432 102 : H5Screate_simple(static_cast<int>(rank), col_dims, nullptr);
433 102 : H5OFFSET_TYPE mem_offset[3] = {0, 0, 0};
434 102 : status = H5Sselect_hyperslab(memspace, H5S_SELECT_SET, mem_offset, nullptr,
435 : count, nullptr);
436 102 : if (status < 0)
437 : {
438 0 : H5Sclose(memspace);
439 0 : return CE_Failure;
440 : }
441 :
442 102 : status = H5Dread(poGDS->dataset_id, poGDS->native, memspace,
443 : poGDS->dataspace_id, H5P_DEFAULT, pImage);
444 :
445 102 : H5Sclose(memspace);
446 :
447 102 : if (status < 0)
448 : {
449 0 : CPLError(CE_Failure, CPLE_AppDefined, "H5Dread() failed for block.");
450 0 : return CE_Failure;
451 : }
452 :
453 : #ifdef HDF5_HAVE_FLOAT16
454 : if (eDataType == GDT_Float32 && poGDS->m_bConvertFromFloat16)
455 : {
456 : for (size_t i = static_cast<size_t>(nBlockXSize) * nBlockYSize; i > 0;
457 : /* do nothing */)
458 : {
459 : --i;
460 : uint16_t nVal16;
461 : memcpy(&nVal16, static_cast<uint16_t *>(pImage) + i,
462 : sizeof(nVal16));
463 : const uint32_t nVal32 = CPLHalfToFloat(nVal16);
464 : float fVal;
465 : memcpy(&fVal, &nVal32, sizeof(fVal));
466 : *(static_cast<float *>(pImage) + i) = fVal;
467 : }
468 : }
469 : else if (eDataType == GDT_CFloat32 && poGDS->m_bConvertFromFloat16)
470 : {
471 : for (size_t i = static_cast<size_t>(nBlockXSize) * nBlockYSize; i > 0;
472 : /* do nothing */)
473 : {
474 : --i;
475 : for (int j = 1; j >= 0; --j)
476 : {
477 : uint16_t nVal16;
478 : memcpy(&nVal16, static_cast<uint16_t *>(pImage) + 2 * i + j,
479 : sizeof(nVal16));
480 : const uint32_t nVal32 = CPLHalfToFloat(nVal16);
481 : float fVal;
482 : memcpy(&fVal, &nVal32, sizeof(fVal));
483 : *(static_cast<float *>(pImage) + 2 * i + j) = fVal;
484 : }
485 : }
486 : }
487 : #endif
488 :
489 102 : return CE_None;
490 : }
491 :
492 : /************************************************************************/
493 : /* IRasterIO() */
494 : /************************************************************************/
495 :
496 311 : CPLErr HDF5ImageRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
497 : int nXSize, int nYSize, void *pData,
498 : int nBufXSize, int nBufYSize,
499 : GDALDataType eBufType,
500 : GSpacing nPixelSpace, GSpacing nLineSpace,
501 : GDALRasterIOExtraArg *psExtraArg)
502 :
503 : {
504 311 : HDF5ImageDataset *poGDS = static_cast<HDF5ImageDataset *>(poDS);
505 :
506 : #ifdef HDF5_HAVE_FLOAT16
507 : if (poGDS->m_bConvertFromFloat16)
508 : {
509 : return GDALPamRasterBand::IRasterIO(
510 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
511 : eBufType, nPixelSpace, nLineSpace, psExtraArg);
512 : }
513 : #endif
514 :
515 : const bool bIsBandInterleavedData =
516 450 : poGDS->ndims == 3 && poGDS->m_nOtherDimIndex == 0 &&
517 761 : poGDS->GetYIndex() == 1 && poGDS->GetXIndex() == 2;
518 :
519 311 : const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
520 :
521 : // Try to detect if we read whole bands by chunks of whole lines
522 : // If so, then read and cache whole band (or group of m_nBandChunkSize bands)
523 : // to save HDF5 decompression.
524 311 : if (m_nIRasterIORecCounter == 0)
525 : {
526 263 : bool bInvalidateWholeBandChunkOptim = false;
527 263 : if (!(nXSize == nBufXSize && nYSize == nBufYSize))
528 : {
529 1 : bInvalidateWholeBandChunkOptim = true;
530 : }
531 : // Is the first request on band 1, line 0 and one or several full lines?
532 262 : else if (poGDS->m_eWholeBandChunkOptim !=
533 70 : HDF5ImageDataset::WBC_ENABLED &&
534 70 : nBand == 1 && nXOff == 0 && nYOff == 0 &&
535 20 : nXSize == nRasterXSize)
536 : {
537 19 : poGDS->m_eWholeBandChunkOptim =
538 : HDF5ImageDataset::WBC_DETECTION_IN_PROGRESS;
539 19 : poGDS->m_nLastRasterIOBand = 1;
540 19 : poGDS->m_nLastRasterIOXOff = nXOff;
541 19 : poGDS->m_nLastRasterIOYOff = nYOff;
542 19 : poGDS->m_nLastRasterIOXSize = nXSize;
543 19 : poGDS->m_nLastRasterIOYSize = nYSize;
544 : }
545 243 : else if (poGDS->m_eWholeBandChunkOptim ==
546 : HDF5ImageDataset::WBC_DETECTION_IN_PROGRESS)
547 : {
548 50 : if (poGDS->m_nLastRasterIOBand == 1 && nBand == 1)
549 : {
550 : // Is this request a continuation of the previous one?
551 44 : if (nXOff == 0 && poGDS->m_nLastRasterIOXOff == 0 &&
552 44 : nYOff == poGDS->m_nLastRasterIOYOff +
553 44 : poGDS->m_nLastRasterIOYSize &&
554 43 : poGDS->m_nLastRasterIOXSize == nRasterXSize &&
555 43 : nXSize == nRasterXSize)
556 : {
557 43 : poGDS->m_nLastRasterIOXOff = nXOff;
558 43 : poGDS->m_nLastRasterIOYOff = nYOff;
559 43 : poGDS->m_nLastRasterIOXSize = nXSize;
560 43 : poGDS->m_nLastRasterIOYSize = nYSize;
561 : }
562 : else
563 : {
564 1 : bInvalidateWholeBandChunkOptim = true;
565 : }
566 : }
567 6 : else if (poGDS->m_nLastRasterIOBand == 1 && nBand == 2)
568 : {
569 : // Are we switching to band 2 while having fully read band 1?
570 5 : if (nXOff == 0 && nYOff == 0 && nXSize == nRasterXSize &&
571 5 : poGDS->m_nLastRasterIOXOff == 0 &&
572 5 : poGDS->m_nLastRasterIOXSize == nRasterXSize &&
573 5 : poGDS->m_nLastRasterIOYOff + poGDS->m_nLastRasterIOYSize ==
574 5 : nRasterYSize)
575 : {
576 11 : if ((poGDS->m_nBandChunkSize > 1 ||
577 5 : nBufYSize < nRasterYSize) &&
578 1 : static_cast<int64_t>(poGDS->m_nBandChunkSize) *
579 1 : nRasterXSize * nRasterYSize * nDTSize <
580 1 : CPLGetUsablePhysicalRAM() / 10)
581 : {
582 1 : poGDS->m_eWholeBandChunkOptim =
583 : HDF5ImageDataset::WBC_ENABLED;
584 : }
585 : else
586 : {
587 3 : bInvalidateWholeBandChunkOptim = true;
588 : }
589 : }
590 : else
591 : {
592 1 : bInvalidateWholeBandChunkOptim = true;
593 : }
594 : }
595 : else
596 : {
597 1 : bInvalidateWholeBandChunkOptim = true;
598 : }
599 : }
600 263 : if (bInvalidateWholeBandChunkOptim)
601 : {
602 7 : poGDS->m_eWholeBandChunkOptim = HDF5ImageDataset::WBC_DISABLED;
603 7 : poGDS->m_nLastRasterIOBand = -1;
604 7 : poGDS->m_nLastRasterIOXOff = -1;
605 7 : poGDS->m_nLastRasterIOYOff = -1;
606 7 : poGDS->m_nLastRasterIOXSize = -1;
607 7 : poGDS->m_nLastRasterIOYSize = -1;
608 : }
609 : }
610 :
611 311 : if (poGDS->m_eWholeBandChunkOptim == HDF5ImageDataset::WBC_ENABLED &&
612 193 : nXSize == nBufXSize && nYSize == nBufYSize)
613 : {
614 193 : if (poGDS->nBands == 1 || bIsBandInterleavedData)
615 : {
616 193 : if (poGDS->m_iCurrentBandChunk < 0)
617 1 : CPLDebug("HDF5", "Using whole band chunk caching");
618 193 : const int iBandChunk = (nBand - 1) / poGDS->m_nBandChunkSize;
619 193 : if (iBandChunk != poGDS->m_iCurrentBandChunk)
620 : {
621 15 : poGDS->m_abyBandChunk.resize(
622 15 : static_cast<size_t>(poGDS->m_nBandChunkSize) *
623 15 : nRasterXSize * nRasterYSize * nDTSize);
624 :
625 : HDF5_GLOBAL_LOCK();
626 :
627 : hsize_t count[3] = {
628 30 : std::min(static_cast<hsize_t>(poGDS->nBands),
629 30 : static_cast<hsize_t>(iBandChunk + 1) *
630 15 : poGDS->m_nBandChunkSize) -
631 15 : static_cast<hsize_t>(iBandChunk) *
632 15 : poGDS->m_nBandChunkSize,
633 15 : static_cast<hsize_t>(nRasterYSize),
634 15 : static_cast<hsize_t>(nRasterXSize)};
635 15 : H5OFFSET_TYPE offset[3] = {
636 15 : static_cast<H5OFFSET_TYPE>(iBandChunk) *
637 15 : poGDS->m_nBandChunkSize,
638 : static_cast<H5OFFSET_TYPE>(0),
639 15 : static_cast<H5OFFSET_TYPE>(0)};
640 : herr_t status =
641 15 : H5Sselect_hyperslab(poGDS->dataspace_id, H5S_SELECT_SET,
642 : offset, nullptr, count, nullptr);
643 15 : if (status < 0)
644 0 : return CE_Failure;
645 :
646 : const hid_t memspace =
647 15 : H5Screate_simple(poGDS->ndims, count, nullptr);
648 15 : H5OFFSET_TYPE mem_offset[3] = {0, 0, 0};
649 : status =
650 15 : H5Sselect_hyperslab(memspace, H5S_SELECT_SET, mem_offset,
651 : nullptr, count, nullptr);
652 15 : if (status < 0)
653 : {
654 0 : H5Sclose(memspace);
655 0 : return CE_Failure;
656 : }
657 :
658 15 : status = H5Dread(poGDS->dataset_id, poGDS->native, memspace,
659 : poGDS->dataspace_id, H5P_DEFAULT,
660 15 : poGDS->m_abyBandChunk.data());
661 :
662 15 : H5Sclose(memspace);
663 :
664 15 : if (status < 0)
665 : {
666 0 : CPLError(
667 : CE_Failure, CPLE_AppDefined,
668 : "HDF5ImageRasterBand::IRasterIO(): H5Dread() failed");
669 0 : return CE_Failure;
670 : }
671 :
672 15 : poGDS->m_iCurrentBandChunk = iBandChunk;
673 : }
674 :
675 1447 : for (int iY = 0; iY < nYSize; ++iY)
676 : {
677 2508 : GDALCopyWords(poGDS->m_abyBandChunk.data() +
678 1254 : static_cast<size_t>((nBand - 1) %
679 1254 : poGDS->m_nBandChunkSize) *
680 1254 : nRasterYSize * nRasterXSize * nDTSize +
681 1254 : static_cast<size_t>(nYOff + iY) *
682 1254 : nRasterXSize * nDTSize +
683 1254 : nXOff * nDTSize,
684 : eDataType, nDTSize,
685 : static_cast<GByte *>(pData) +
686 1254 : static_cast<size_t>(iY) * nLineSpace,
687 : eBufType, static_cast<int>(nPixelSpace), nXSize);
688 : }
689 193 : return CE_None;
690 : }
691 : }
692 :
693 : const bool bIsExpectedLayout =
694 204 : (bIsBandInterleavedData ||
695 171 : (poGDS->ndims == 2 && poGDS->GetYIndex() == 0 &&
696 85 : poGDS->GetXIndex() == 1));
697 118 : if (eRWFlag == GF_Read && bIsExpectedLayout && nXSize == nBufXSize &&
698 116 : nYSize == nBufYSize && eBufType == eDataType &&
699 68 : nPixelSpace == nDTSize && nLineSpace == nXSize * nPixelSpace)
700 : {
701 : HDF5_GLOBAL_LOCK();
702 :
703 68 : hsize_t count[3] = {1, static_cast<hsize_t>(nYSize),
704 68 : static_cast<hsize_t>(nXSize)};
705 68 : H5OFFSET_TYPE offset[3] = {static_cast<H5OFFSET_TYPE>(nBand - 1),
706 68 : static_cast<H5OFFSET_TYPE>(nYOff),
707 68 : static_cast<H5OFFSET_TYPE>(nXOff)};
708 68 : if (poGDS->ndims == 2)
709 : {
710 47 : count[0] = count[1];
711 47 : count[1] = count[2];
712 :
713 47 : offset[0] = offset[1];
714 47 : offset[1] = offset[2];
715 : }
716 68 : herr_t status = H5Sselect_hyperslab(poGDS->dataspace_id, H5S_SELECT_SET,
717 : offset, nullptr, count, nullptr);
718 68 : if (status < 0)
719 0 : return CE_Failure;
720 :
721 68 : const hid_t memspace = H5Screate_simple(poGDS->ndims, count, nullptr);
722 68 : H5OFFSET_TYPE mem_offset[3] = {0, 0, 0};
723 68 : status = H5Sselect_hyperslab(memspace, H5S_SELECT_SET, mem_offset,
724 : nullptr, count, nullptr);
725 68 : if (status < 0)
726 : {
727 0 : H5Sclose(memspace);
728 0 : return CE_Failure;
729 : }
730 :
731 68 : status = H5Dread(poGDS->dataset_id, poGDS->native, memspace,
732 : poGDS->dataspace_id, H5P_DEFAULT, pData);
733 :
734 68 : H5Sclose(memspace);
735 :
736 68 : if (status < 0)
737 : {
738 0 : CPLError(CE_Failure, CPLE_AppDefined,
739 : "HDF5ImageRasterBand::IRasterIO(): H5Dread() failed");
740 0 : return CE_Failure;
741 : }
742 :
743 68 : return CE_None;
744 : }
745 :
746 : // If the request is still small enough, try to read from libhdf5 with
747 : // the natural interleaving into a temporary MEMDataset, and then read
748 : // from it with the requested interleaving and data type.
749 50 : if (eRWFlag == GF_Read && bIsExpectedLayout && nXSize == nBufXSize &&
750 100 : nYSize == nBufYSize &&
751 48 : static_cast<GIntBig>(nXSize) * nYSize < CPLGetUsablePhysicalRAM() / 10)
752 : {
753 : auto poMemDS = std::unique_ptr<GDALDataset>(
754 48 : MEMDataset::Create("", nXSize, nYSize, 1, eDataType, nullptr));
755 48 : if (poMemDS)
756 : {
757 48 : void *pMemData = poMemDS->GetInternalHandle("MEMORY1");
758 48 : CPLAssert(pMemData);
759 : // Read from HDF5 into the temporary MEMDataset using the
760 : // natural interleaving of the HDF5 dataset
761 48 : ++m_nIRasterIORecCounter;
762 : CPLErr eErr =
763 96 : IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pMemData,
764 : nXSize, nYSize, eDataType, nDTSize,
765 48 : static_cast<GSpacing>(nXSize) * nDTSize, psExtraArg);
766 48 : --m_nIRasterIORecCounter;
767 48 : if (eErr != CE_None)
768 : {
769 0 : return CE_Failure;
770 : }
771 : // Copy to the final buffer using requested data type and spacings.
772 48 : return poMemDS->GetRasterBand(1)->RasterIO(
773 : GF_Read, 0, 0, nXSize, nYSize, pData, nXSize, nYSize, eBufType,
774 48 : nPixelSpace, nLineSpace, nullptr);
775 : }
776 : }
777 :
778 2 : return GDALPamRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
779 : pData, nBufXSize, nBufYSize, eBufType,
780 2 : nPixelSpace, nLineSpace, psExtraArg);
781 : }
782 :
783 : /************************************************************************/
784 : /* IRasterIO() */
785 : /************************************************************************/
786 :
787 12 : CPLErr HDF5ImageDataset::IRasterIO(
788 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
789 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
790 : int nBandCount, BANDMAP_TYPE panBandMap, GSpacing nPixelSpace,
791 : GSpacing nLineSpace, GSpacing nBandSpace, GDALRasterIOExtraArg *psExtraArg)
792 :
793 : {
794 : #ifdef HDF5_HAVE_FLOAT16
795 : if (m_bConvertFromFloat16)
796 : {
797 : return HDF5Dataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
798 : pData, nBufXSize, nBufYSize, eBufType,
799 : nBandCount, panBandMap, nPixelSpace,
800 : nLineSpace, nBandSpace, psExtraArg);
801 : }
802 : #endif
803 :
804 16 : const auto IsConsecutiveBands = [](const int *panVals, int nCount)
805 : {
806 71 : for (int i = 1; i < nCount; ++i)
807 : {
808 57 : if (panVals[i] != panVals[i - 1] + 1)
809 2 : return false;
810 : }
811 14 : return true;
812 : };
813 :
814 12 : const auto eDT = GetRasterBand(1)->GetRasterDataType();
815 12 : const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
816 :
817 : // Band-interleaved data and request
818 20 : const bool bIsBandInterleavedData = ndims == 3 && m_nOtherDimIndex == 0 &&
819 32 : GetYIndex() == 1 && GetXIndex() == 2;
820 12 : if (eRWFlag == GF_Read && bIsBandInterleavedData && nXSize == nBufXSize &&
821 8 : nYSize == nBufYSize && IsConsecutiveBands(panBandMap, nBandCount) &&
822 6 : eBufType == eDT && nPixelSpace == nDTSize &&
823 24 : nLineSpace == nXSize * nPixelSpace && nBandSpace == nYSize * nLineSpace)
824 : {
825 : HDF5_GLOBAL_LOCK();
826 :
827 5 : hsize_t count[3] = {static_cast<hsize_t>(nBandCount),
828 5 : static_cast<hsize_t>(nYSize),
829 5 : static_cast<hsize_t>(nXSize)};
830 : H5OFFSET_TYPE offset[3] = {
831 5 : static_cast<H5OFFSET_TYPE>(panBandMap[0] - 1),
832 5 : static_cast<H5OFFSET_TYPE>(nYOff),
833 5 : static_cast<H5OFFSET_TYPE>(nXOff)};
834 5 : herr_t status = H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET,
835 : offset, nullptr, count, nullptr);
836 5 : if (status < 0)
837 0 : return CE_Failure;
838 :
839 5 : const hid_t memspace = H5Screate_simple(ndims, count, nullptr);
840 5 : H5OFFSET_TYPE mem_offset[3] = {0, 0, 0};
841 5 : status = H5Sselect_hyperslab(memspace, H5S_SELECT_SET, mem_offset,
842 : nullptr, count, nullptr);
843 5 : if (status < 0)
844 : {
845 0 : H5Sclose(memspace);
846 0 : return CE_Failure;
847 : }
848 :
849 5 : status = H5Dread(dataset_id, native, memspace, dataspace_id,
850 : H5P_DEFAULT, pData);
851 :
852 5 : H5Sclose(memspace);
853 :
854 5 : if (status < 0)
855 : {
856 0 : CPLError(CE_Failure, CPLE_AppDefined,
857 : "HDF5ImageDataset::IRasterIO(): H5Dread() failed");
858 0 : return CE_Failure;
859 : }
860 :
861 5 : return CE_None;
862 : }
863 :
864 : // Pixel-interleaved data and request
865 :
866 11 : const bool bIsPixelInterleaveData = ndims == 3 && m_nOtherDimIndex == 2 &&
867 18 : GetYIndex() == 0 && GetXIndex() == 1;
868 7 : if (eRWFlag == GF_Read && bIsPixelInterleaveData && nXSize == nBufXSize &&
869 4 : nYSize == nBufYSize && IsConsecutiveBands(panBandMap, nBandCount) &&
870 4 : eBufType == eDT && nBandSpace == nDTSize &&
871 17 : nPixelSpace == nBandCount * nBandSpace &&
872 3 : nLineSpace == nXSize * nPixelSpace)
873 : {
874 : HDF5_GLOBAL_LOCK();
875 :
876 3 : hsize_t count[3] = {static_cast<hsize_t>(nYSize),
877 3 : static_cast<hsize_t>(nXSize),
878 3 : static_cast<hsize_t>(nBandCount)};
879 : H5OFFSET_TYPE offset[3] = {
880 3 : static_cast<H5OFFSET_TYPE>(nYOff),
881 3 : static_cast<H5OFFSET_TYPE>(nXOff),
882 3 : static_cast<H5OFFSET_TYPE>(panBandMap[0] - 1)};
883 3 : herr_t status = H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET,
884 : offset, nullptr, count, nullptr);
885 3 : if (status < 0)
886 0 : return CE_Failure;
887 :
888 3 : const hid_t memspace = H5Screate_simple(ndims, count, nullptr);
889 3 : H5OFFSET_TYPE mem_offset[3] = {0, 0, 0};
890 3 : status = H5Sselect_hyperslab(memspace, H5S_SELECT_SET, mem_offset,
891 : nullptr, count, nullptr);
892 3 : if (status < 0)
893 : {
894 0 : H5Sclose(memspace);
895 0 : return CE_Failure;
896 : }
897 :
898 3 : status = H5Dread(dataset_id, native, memspace, dataspace_id,
899 : H5P_DEFAULT, pData);
900 :
901 3 : H5Sclose(memspace);
902 :
903 3 : if (status < 0)
904 : {
905 0 : CPLError(CE_Failure, CPLE_AppDefined,
906 : "HDF5ImageDataset::IRasterIO(): H5Dread() failed");
907 0 : return CE_Failure;
908 : }
909 :
910 3 : return CE_None;
911 : }
912 :
913 : // If the request is still small enough, try to read from libhdf5 with
914 : // the natural interleaving into a temporary MEMDataset, and then read
915 : // from it with the requested interleaving and data type.
916 4 : if (eRWFlag == GF_Read &&
917 4 : (bIsBandInterleavedData || bIsPixelInterleaveData) &&
918 8 : nXSize == nBufXSize && nYSize == nBufYSize &&
919 12 : IsConsecutiveBands(panBandMap, nBandCount) &&
920 3 : static_cast<GIntBig>(nXSize) * nYSize <
921 3 : CPLGetUsablePhysicalRAM() / 10 / nBandCount)
922 : {
923 3 : const char *const apszOptions[] = {
924 3 : bIsPixelInterleaveData ? "INTERLEAVE=PIXEL" : nullptr, nullptr};
925 : auto poMemDS = std::unique_ptr<GDALDataset>(
926 3 : MEMDataset::Create("", nXSize, nYSize, nBandCount, eDT,
927 3 : const_cast<char **>(apszOptions)));
928 3 : if (poMemDS)
929 : {
930 3 : void *pMemData = poMemDS->GetInternalHandle("MEMORY1");
931 3 : CPLAssert(pMemData);
932 : // Read from HDF5 into the temporary MEMDataset using the
933 : // natural interleaving of the HDF5 dataset
934 9 : if (IRasterIO(
935 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pMemData, nXSize,
936 : nYSize, eDT, nBandCount, panBandMap,
937 1 : bIsBandInterleavedData ? nDTSize : nDTSize * nBandCount,
938 : bIsBandInterleavedData
939 2 : ? static_cast<GSpacing>(nXSize) * nDTSize
940 1 : : static_cast<GSpacing>(nXSize) * nDTSize * nBandCount,
941 : bIsBandInterleavedData
942 2 : ? static_cast<GSpacing>(nYSize) * nXSize * nDTSize
943 : : nDTSize,
944 3 : psExtraArg) != CE_None)
945 : {
946 0 : return CE_Failure;
947 : }
948 : // Copy to the final buffer using requested data type and spacings.
949 3 : return poMemDS->RasterIO(GF_Read, 0, 0, nXSize, nYSize, pData,
950 : nXSize, nYSize, eBufType, nBandCount,
951 : nullptr, nPixelSpace, nLineSpace,
952 3 : nBandSpace, nullptr);
953 : }
954 : }
955 :
956 1 : return HDF5Dataset::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize, pData,
957 : nBufXSize, nBufYSize, eBufType, nBandCount,
958 : panBandMap, nPixelSpace, nLineSpace,
959 1 : nBandSpace, psExtraArg);
960 : }
961 :
962 : /************************************************************************/
963 : /* Open() */
964 : /************************************************************************/
965 63 : GDALDataset *HDF5ImageDataset::Open(GDALOpenInfo *poOpenInfo)
966 : {
967 63 : if (!STARTS_WITH_CI(poOpenInfo->pszFilename, "HDF5:"))
968 0 : return nullptr;
969 :
970 : HDF5_GLOBAL_LOCK();
971 :
972 : // Confirm the requested access is supported.
973 63 : if (poOpenInfo->eAccess == GA_Update)
974 : {
975 0 : ReportUpdateNotSupportedByDriver("HDF5");
976 0 : return nullptr;
977 : }
978 :
979 126 : auto poDS = std::make_unique<HDF5ImageDataset>();
980 :
981 : // Create a corresponding GDALDataset.
982 : const CPLStringList aosFilenameTokens(
983 63 : CSLTokenizeString2(poOpenInfo->pszFilename, ":",
984 126 : CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES));
985 :
986 63 : if (aosFilenameTokens.size() != 3 && aosFilenameTokens.size() != 4)
987 : {
988 0 : return nullptr;
989 : }
990 :
991 63 : poDS->SetDescription(poOpenInfo->pszFilename);
992 :
993 : // Check for drive name in windows HDF5:"D:\...
994 126 : CPLString osSubdatasetName;
995 :
996 126 : CPLString osFilename(aosFilenameTokens[1]);
997 :
998 63 : if ((strlen(aosFilenameTokens[1]) == 1 &&
999 126 : aosFilenameTokens[3] != nullptr) ||
1000 63 : (STARTS_WITH(aosFilenameTokens[1], "/vsicurl/http") &&
1001 0 : aosFilenameTokens[3] != nullptr))
1002 : {
1003 0 : osFilename += ":";
1004 0 : osFilename += aosFilenameTokens[2];
1005 0 : osSubdatasetName = aosFilenameTokens[3];
1006 : }
1007 : else
1008 : {
1009 63 : osSubdatasetName = aosFilenameTokens[2];
1010 : }
1011 :
1012 63 : poDS->SetSubdatasetName(osSubdatasetName);
1013 :
1014 63 : poDS->SetPhysicalFilename(osFilename);
1015 :
1016 : // Try opening the dataset.
1017 63 : poDS->m_hHDF5 = GDAL_HDF5Open(osFilename);
1018 63 : if (poDS->m_hHDF5 < 0)
1019 : {
1020 0 : return nullptr;
1021 : }
1022 :
1023 63 : poDS->hGroupID = H5Gopen(poDS->m_hHDF5, "/");
1024 63 : if (poDS->hGroupID < 0)
1025 : {
1026 0 : return nullptr;
1027 : }
1028 :
1029 : // This is an HDF5 file.
1030 63 : poDS->ReadGlobalAttributes(FALSE);
1031 :
1032 : // Create HDF5 Data Hierarchy in a link list.
1033 126 : HDF5GroupObjects *poH5Objects = poDS->HDF5FindDatasetObjectsbyPath(
1034 63 : poDS->poH5RootGroup, osSubdatasetName);
1035 :
1036 63 : if (poH5Objects == nullptr)
1037 : {
1038 0 : return nullptr;
1039 : }
1040 :
1041 : // Retrieve HDF5 data information.
1042 63 : poDS->m_osPath = poH5Objects->pszPath;
1043 63 : poDS->dataset_id = H5Dopen(poDS->m_hHDF5, poH5Objects->pszPath);
1044 63 : poDS->dataspace_id = H5Dget_space(poDS->dataset_id);
1045 63 : poDS->ndims = H5Sget_simple_extent_ndims(poDS->dataspace_id);
1046 63 : if (poDS->ndims <= 0)
1047 : {
1048 0 : return nullptr;
1049 : }
1050 126 : poDS->dims =
1051 63 : static_cast<hsize_t *>(CPLCalloc(poDS->ndims, sizeof(hsize_t)));
1052 126 : poDS->maxdims =
1053 63 : static_cast<hsize_t *>(CPLCalloc(poDS->ndims, sizeof(hsize_t)));
1054 63 : if (H5Sget_simple_extent_dims(poDS->dataspace_id, poDS->dims,
1055 126 : poDS->maxdims) != poDS->ndims)
1056 : {
1057 0 : return nullptr;
1058 : }
1059 201 : for (int i = 0; i < poDS->ndims; ++i)
1060 : {
1061 139 : if (poDS->dims[i] >
1062 139 : static_cast<hsize_t>(std::numeric_limits<int>::max()))
1063 : {
1064 1 : CPLError(CE_Failure, CPLE_NotSupported,
1065 : "At least one dimension size exceeds INT_MAX !");
1066 1 : return nullptr;
1067 : }
1068 : }
1069 :
1070 62 : auto datatype = H5Dget_type(poDS->dataset_id);
1071 62 : poDS->native = H5Tget_native_type(datatype, H5T_DIR_ASCEND);
1072 62 : H5Tclose(datatype);
1073 :
1074 62 : const auto eGDALDataType = poDS->GetDataType(poDS->native);
1075 62 : if (eGDALDataType == GDT_Unknown)
1076 : {
1077 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unhandled HDF5 data type");
1078 0 : return nullptr;
1079 : }
1080 :
1081 : #ifdef HDF5_HAVE_FLOAT16
1082 : if (H5Tequal(H5T_NATIVE_FLOAT16, poDS->native) ||
1083 : IsNativeCFloat16(poDS->native))
1084 : {
1085 : poDS->m_bConvertFromFloat16 = true;
1086 : }
1087 : #endif
1088 :
1089 : // CSK code in IdentifyProductType() and CreateProjections()
1090 : // uses dataset metadata.
1091 62 : poDS->SetMetadata(poDS->m_aosMetadata.List());
1092 :
1093 : // Check if the hdf5 is a well known product type
1094 62 : poDS->IdentifyProductType();
1095 :
1096 62 : poDS->m_nYIndex = poDS->IsComplexCSKL1A() ? 0 : poDS->ndims - 2;
1097 62 : poDS->m_nXIndex = poDS->IsComplexCSKL1A() ? 1 : poDS->ndims - 1;
1098 :
1099 62 : if (poDS->IsComplexCSKL1A())
1100 : {
1101 0 : poDS->m_nOtherDimIndex = 2;
1102 : }
1103 62 : else if (poDS->ndims == 3)
1104 : {
1105 15 : poDS->m_nOtherDimIndex = 0;
1106 : }
1107 :
1108 62 : if (HDF5EOSParser::HasHDFEOS(poDS->hGroupID))
1109 : {
1110 40 : HDF5EOSParser oHDFEOSParser;
1111 20 : if (oHDFEOSParser.Parse(poDS->hGroupID))
1112 : {
1113 20 : CPLDebug("HDF5", "Successfully parsed HDFEOS metadata");
1114 40 : HDF5EOSParser::GridDataFieldMetadata oGridDataFieldMetadata;
1115 40 : HDF5EOSParser::SwathFieldMetadata oSwathFieldMetadata;
1116 20 : if (oHDFEOSParser.GetDataModel() ==
1117 6 : HDF5EOSParser::DataModel::GRID &&
1118 6 : oHDFEOSParser.GetGridDataFieldMetadata(
1119 26 : osSubdatasetName.c_str(), oGridDataFieldMetadata) &&
1120 6 : static_cast<int>(oGridDataFieldMetadata.aoDimensions.size()) ==
1121 6 : poDS->ndims)
1122 : {
1123 6 : int iDim = 0;
1124 20 : for (const auto &oDim : oGridDataFieldMetadata.aoDimensions)
1125 : {
1126 14 : if (oDim.osName == "XDim")
1127 6 : poDS->m_nXIndex = iDim;
1128 8 : else if (oDim.osName == "YDim")
1129 6 : poDS->m_nYIndex = iDim;
1130 : else
1131 2 : poDS->m_nOtherDimIndex = iDim;
1132 14 : ++iDim;
1133 : }
1134 :
1135 6 : if (oGridDataFieldMetadata.poGridMetadata->GetGeoTransform(
1136 6 : poDS->m_gt))
1137 6 : poDS->bHasGeoTransform = true;
1138 :
1139 12 : auto poSRS = oGridDataFieldMetadata.poGridMetadata->GetSRS();
1140 6 : if (poSRS)
1141 6 : poDS->m_oSRS = *(poSRS.get());
1142 : }
1143 14 : else if (oHDFEOSParser.GetDataModel() ==
1144 14 : HDF5EOSParser::DataModel::SWATH &&
1145 14 : oHDFEOSParser.GetSwathFieldMetadata(
1146 14 : osSubdatasetName.c_str(), oSwathFieldMetadata) &&
1147 : static_cast<int>(
1148 14 : oSwathFieldMetadata.aoDimensions.size()) ==
1149 14 : poDS->ndims &&
1150 42 : oSwathFieldMetadata.iXDim >= 0 &&
1151 14 : oSwathFieldMetadata.iYDim >= 0)
1152 : {
1153 14 : poDS->m_nXIndex = oSwathFieldMetadata.iXDim;
1154 14 : poDS->m_nYIndex = oSwathFieldMetadata.iYDim;
1155 14 : poDS->m_nOtherDimIndex = oSwathFieldMetadata.iOtherDim;
1156 14 : if (!oSwathFieldMetadata.osLongitudeSubdataset.empty())
1157 : {
1158 28 : CPLStringList aosGeolocation;
1159 :
1160 : // Arbitrary
1161 14 : aosGeolocation.AddNameValue("SRS", SRS_WKT_WGS84_LAT_LONG);
1162 : aosGeolocation.AddNameValue(
1163 28 : "X_DATASET", ("HDF5:\"" + osFilename + "\":" +
1164 : oSwathFieldMetadata.osLongitudeSubdataset)
1165 14 : .c_str());
1166 14 : aosGeolocation.AddNameValue("X_BAND", "1");
1167 : aosGeolocation.AddNameValue(
1168 28 : "Y_DATASET", ("HDF5:\"" + osFilename + "\":" +
1169 : oSwathFieldMetadata.osLatitudeSubdataset)
1170 14 : .c_str());
1171 14 : aosGeolocation.AddNameValue("Y_BAND", "1");
1172 : aosGeolocation.AddNameValue(
1173 : "PIXEL_OFFSET",
1174 14 : CPLSPrintf("%d", oSwathFieldMetadata.nPixelOffset));
1175 : aosGeolocation.AddNameValue(
1176 : "PIXEL_STEP",
1177 14 : CPLSPrintf("%d", oSwathFieldMetadata.nPixelStep));
1178 : aosGeolocation.AddNameValue(
1179 : "LINE_OFFSET",
1180 14 : CPLSPrintf("%d", oSwathFieldMetadata.nLineOffset));
1181 : aosGeolocation.AddNameValue(
1182 : "LINE_STEP",
1183 14 : CPLSPrintf("%d", oSwathFieldMetadata.nLineStep));
1184 : // Not totally sure about that
1185 : aosGeolocation.AddNameValue("GEOREFERENCING_CONVENTION",
1186 14 : "PIXEL_CENTER");
1187 28 : poDS->SetMetadata(aosGeolocation.List(),
1188 14 : GDAL_MDD_GEOLOCATION);
1189 : }
1190 : }
1191 : }
1192 : }
1193 :
1194 62 : poDS->nRasterYSize =
1195 62 : poDS->GetYIndex() < 0
1196 62 : ? 1
1197 61 : : static_cast<int>(poDS->dims[poDS->GetYIndex()]); // nRows
1198 124 : poDS->nRasterXSize =
1199 62 : static_cast<int>(poDS->dims[poDS->GetXIndex()]); // nCols
1200 62 : int nBands = 1;
1201 62 : if (poDS->m_nOtherDimIndex >= 0)
1202 : {
1203 15 : nBands = static_cast<int>(poDS->dims[poDS->m_nOtherDimIndex]);
1204 : }
1205 :
1206 124 : CPLStringList aosMetadata;
1207 124 : std::map<std::string, CPLStringList> oMapBandSpecificMetadata;
1208 62 : if (poH5Objects->nType == H5G_DATASET)
1209 : {
1210 62 : HDF5Dataset::CreateMetadata(poDS->m_hHDF5, poH5Objects, H5G_DATASET,
1211 : false, aosMetadata);
1212 76 : if (nBands > 1 && poDS->nRasterXSize != nBands &&
1213 14 : poDS->nRasterYSize != nBands)
1214 : {
1215 : // Heuristics to detect non-scalar attributes, that are intended
1216 : // to be attached to a specific band.
1217 28 : const CPLStringList aosMetadataDup(aosMetadata);
1218 52 : for (const auto &[pszKey, pszValue] :
1219 66 : cpl::IterateNameValue(aosMetadataDup))
1220 : {
1221 26 : const hid_t hAttrID = H5Aopen_name(poDS->dataset_id, pszKey);
1222 26 : const hid_t hAttrSpace = H5Aget_space(hAttrID);
1223 47 : if (H5Sget_simple_extent_ndims(hAttrSpace) == 1 &&
1224 21 : H5Sget_simple_extent_npoints(hAttrSpace) == nBands)
1225 : {
1226 : CPLStringList aosTokens(
1227 40 : CSLTokenizeString2(pszValue, " ", 0));
1228 20 : if (aosTokens.size() == nBands)
1229 : {
1230 40 : std::string osAttrName(pszKey);
1231 30 : if (osAttrName.size() > strlen("_coefficients") &&
1232 30 : osAttrName.substr(osAttrName.size() -
1233 10 : strlen("_coefficients")) ==
1234 : "_coefficients")
1235 : {
1236 5 : osAttrName.pop_back();
1237 : }
1238 25 : else if (osAttrName.size() > strlen("_wavelengths") &&
1239 25 : osAttrName.substr(osAttrName.size() -
1240 10 : strlen("_wavelengths")) ==
1241 : "_wavelengths")
1242 : {
1243 5 : osAttrName.pop_back();
1244 : }
1245 15 : else if (osAttrName.size() > strlen("_list") &&
1246 15 : osAttrName.substr(osAttrName.size() -
1247 5 : strlen("_list")) == "_list")
1248 : {
1249 5 : osAttrName.resize(osAttrName.size() -
1250 : strlen("_list"));
1251 : }
1252 20 : oMapBandSpecificMetadata[osAttrName] =
1253 40 : std::move(aosTokens);
1254 20 : aosMetadata.SetNameValue(pszKey, nullptr);
1255 : }
1256 : }
1257 26 : H5Sclose(hAttrSpace);
1258 26 : H5Aclose(hAttrID);
1259 : }
1260 : }
1261 : }
1262 :
1263 62 : poDS->m_nBlockXSize = poDS->GetRasterXSize();
1264 62 : poDS->m_nBlockYSize = 1;
1265 62 : poDS->m_nBandChunkSize = 1;
1266 :
1267 : // Check for chunksize and set it as the blocksize (optimizes read).
1268 62 : const hid_t listid = H5Dget_create_plist(poDS->dataset_id);
1269 62 : if (listid > 0)
1270 : {
1271 62 : if (H5Pget_layout(listid) == H5D_CHUNKED)
1272 : {
1273 12 : hsize_t panChunkDims[3] = {0, 0, 0};
1274 12 : const int nDimSize = H5Pget_chunk(listid, 3, panChunkDims);
1275 12 : CPL_IGNORE_RET_VAL(nDimSize);
1276 12 : CPLAssert(nDimSize == poDS->ndims);
1277 24 : poDS->m_nBlockXSize =
1278 12 : static_cast<int>(panChunkDims[poDS->GetXIndex()]);
1279 12 : if (poDS->GetYIndex() >= 0)
1280 12 : poDS->m_nBlockYSize =
1281 12 : static_cast<int>(panChunkDims[poDS->GetYIndex()]);
1282 12 : if (nBands > 1)
1283 : {
1284 6 : poDS->m_nBandChunkSize =
1285 3 : static_cast<int>(panChunkDims[poDS->m_nOtherDimIndex]);
1286 :
1287 3 : if (poDS->m_nBandChunkSize > 1)
1288 3 : poDS->SetMetadataItem(GDALMD_INTERLEAVE, "PIXEL",
1289 3 : GDAL_MDD_IMAGE_STRUCTURE);
1290 :
1291 6 : poDS->SetMetadataItem("BAND_CHUNK_SIZE",
1292 3 : CPLSPrintf("%d", poDS->m_nBandChunkSize),
1293 3 : GDAL_MDD_IMAGE_STRUCTURE);
1294 : }
1295 : }
1296 :
1297 62 : const int nFilters = H5Pget_nfilters(listid);
1298 68 : for (int i = 0; i < nFilters; ++i)
1299 : {
1300 6 : unsigned int flags = 0;
1301 6 : size_t cd_nelmts = 0;
1302 6 : char szName[64 + 1] = {0};
1303 6 : const auto eFilter = H5Pget_filter(listid, i, &flags, &cd_nelmts,
1304 : nullptr, 64, szName);
1305 6 : if (eFilter == H5Z_FILTER_DEFLATE)
1306 : {
1307 5 : poDS->SetMetadataItem(GDALMD_COMPRESSION, "DEFLATE",
1308 5 : GDAL_MDD_IMAGE_STRUCTURE);
1309 : }
1310 1 : else if (eFilter == H5Z_FILTER_SZIP)
1311 : {
1312 0 : poDS->SetMetadataItem(GDALMD_COMPRESSION, "SZIP",
1313 0 : GDAL_MDD_IMAGE_STRUCTURE);
1314 : }
1315 : }
1316 :
1317 62 : H5Pclose(listid);
1318 : }
1319 :
1320 191 : for (int i = 0; i < nBands; i++)
1321 : {
1322 129 : auto poBand = std::make_unique<HDF5ImageRasterBand>(poDS.get(), i + 1,
1323 129 : eGDALDataType);
1324 :
1325 129 : if (poH5Objects->nType == H5G_DATASET)
1326 : {
1327 129 : poBand->SetMetadata(aosMetadata.List());
1328 169 : for (const auto &oIter : oMapBandSpecificMetadata)
1329 : {
1330 40 : poBand->SetMetadataItem(oIter.first.c_str(), oIter.second[i]);
1331 : }
1332 : }
1333 :
1334 129 : poDS->SetBand(i + 1, std::move(poBand));
1335 : }
1336 :
1337 62 : if (!poDS->GetMetadata(GDAL_MDD_GEOLOCATION))
1338 48 : poDS->CreateProjections();
1339 :
1340 : // Setup/check for pam .aux.xml.
1341 62 : poDS->TryLoadXML();
1342 :
1343 : // If the PAM .aux.xml file contains the serialized GEOLOCATION metadata
1344 : // domain, make sure to patch the X_DATASET and Y_DATASET values, when
1345 : // they point to the current file, to use the filename with which we have
1346 : // opened the dataset. Helps in scenarios where the .aux.xml file has
1347 : // been saved with a relative filename for example.
1348 : // Or scenarios like https://github.com/OSGeo/gdal/issues/12824 mixing
1349 : // native Windows and WSL use.
1350 : CSLConstList papszGeoLocationAfterLoadXML =
1351 62 : poDS->GetMetadata(GDAL_MDD_GEOLOCATION);
1352 62 : if (papszGeoLocationAfterLoadXML)
1353 : {
1354 42 : for (const char *pszItem : {"X_DATASET", "Y_DATASET"})
1355 : {
1356 : const char *pszSubdataset =
1357 28 : CSLFetchNameValue(papszGeoLocationAfterLoadXML, pszItem);
1358 56 : if (pszSubdataset && STARTS_WITH(pszSubdataset, "HDF5:\"") &&
1359 28 : strstr(pszSubdataset, CPLGetFilename(osFilename.c_str())))
1360 : {
1361 28 : auto hSubDSInfo = GDALGetSubdatasetInfo(pszSubdataset);
1362 28 : if (hSubDSInfo)
1363 : {
1364 : char *pszOriPath =
1365 28 : GDALSubdatasetInfoGetPathComponent(hSubDSInfo);
1366 28 : if (EQUAL(CPLGetFilename(pszOriPath),
1367 : CPLGetFilename(osFilename.c_str())))
1368 : {
1369 28 : char *pszNewVal = GDALSubdatasetInfoModifyPathComponent(
1370 : hSubDSInfo, osFilename.c_str());
1371 28 : poDS->SetMetadataItem(pszItem, pszNewVal,
1372 28 : GDAL_MDD_GEOLOCATION);
1373 28 : CPLFree(pszNewVal);
1374 : }
1375 28 : CPLFree(pszOriPath);
1376 28 : GDALDestroySubdatasetInfo(hSubDSInfo);
1377 : }
1378 : }
1379 : }
1380 : }
1381 :
1382 62 : poDS->SetPamFlags(poDS->GetPamFlags() & ~GPF_DIRTY);
1383 :
1384 : // Setup overviews.
1385 62 : poDS->oOvManager.Initialize(poDS.get(), ":::VIRTUAL:::");
1386 :
1387 62 : return poDS.release();
1388 : }
1389 :
1390 : /************************************************************************/
1391 : /* HDF5ImageDatasetDriverUnload() */
1392 : /************************************************************************/
1393 :
1394 9 : static void HDF5ImageDatasetDriverUnload(GDALDriver *)
1395 : {
1396 9 : HDF5UnloadFileDriver();
1397 9 : }
1398 :
1399 : /************************************************************************/
1400 : /* GDALRegister_HDF5Image() */
1401 : /************************************************************************/
1402 14 : void GDALRegister_HDF5Image()
1403 :
1404 : {
1405 14 : if (!GDAL_CHECK_VERSION("HDF5Image driver"))
1406 0 : return;
1407 :
1408 14 : if (GDALGetDriverByName(HDF5_IMAGE_DRIVER_NAME) != nullptr)
1409 0 : return;
1410 :
1411 14 : GDALDriver *poDriver = new GDALDriver();
1412 :
1413 14 : HDF5ImageDriverSetCommonMetadata(poDriver);
1414 :
1415 14 : poDriver->pfnOpen = HDF5ImageDataset::Open;
1416 14 : poDriver->pfnUnloadDriver = HDF5ImageDatasetDriverUnload;
1417 :
1418 14 : GetGDALDriverManager()->RegisterDriver(poDriver);
1419 : }
1420 :
1421 : /************************************************************************/
1422 : /* CreateODIMH5Projection() */
1423 : /************************************************************************/
1424 :
1425 : // Reference:
1426 : // http://www.knmi.nl/opera/opera3/OPERA_2008_03_WP2.1b_ODIM_H5_v2.1.pdf
1427 : //
1428 : // 4.3.2 where for geographically referenced image Groups
1429 : // We don't use the where_xscale and where_yscale parameters, but recompute them
1430 : // from the lower-left and upper-right coordinates. There's some difference.
1431 : // As all those parameters are linked together, I'm not sure which one should be
1432 : // considered as the reference.
1433 :
1434 0 : CPLErr HDF5ImageDataset::CreateODIMH5Projection()
1435 : {
1436 0 : const char *const pszProj4String = GetMetadataItem("where_projdef");
1437 0 : const char *const pszLL_lon = GetMetadataItem("where_LL_lon");
1438 0 : const char *const pszLL_lat = GetMetadataItem("where_LL_lat");
1439 0 : const char *const pszUR_lon = GetMetadataItem("where_UR_lon");
1440 0 : const char *const pszUR_lat = GetMetadataItem("where_UR_lat");
1441 0 : if (pszProj4String == nullptr || pszLL_lon == nullptr ||
1442 0 : pszLL_lat == nullptr || pszUR_lon == nullptr || pszUR_lat == nullptr)
1443 0 : return CE_Failure;
1444 :
1445 0 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1446 0 : if (m_oSRS.importFromProj4(pszProj4String) != OGRERR_NONE)
1447 0 : return CE_Failure;
1448 :
1449 0 : OGRSpatialReference oSRSWGS84;
1450 0 : oSRSWGS84.SetWellKnownGeogCS("WGS84");
1451 0 : oSRSWGS84.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1452 :
1453 : OGRCoordinateTransformation *poCT =
1454 0 : OGRCreateCoordinateTransformation(&oSRSWGS84, &m_oSRS);
1455 0 : if (poCT == nullptr)
1456 0 : return CE_Failure;
1457 :
1458 : // Reproject corners from long,lat WGS84 to the target SRS.
1459 0 : double dfLLX = CPLAtof(pszLL_lon);
1460 0 : double dfLLY = CPLAtof(pszLL_lat);
1461 0 : double dfURX = CPLAtof(pszUR_lon);
1462 0 : double dfURY = CPLAtof(pszUR_lat);
1463 0 : if (!poCT->Transform(1, &dfLLX, &dfLLY) ||
1464 0 : !poCT->Transform(1, &dfURX, &dfURY))
1465 : {
1466 0 : delete poCT;
1467 0 : return CE_Failure;
1468 : }
1469 0 : delete poCT;
1470 :
1471 : // Compute the geotransform now.
1472 0 : const double dfPixelX = (dfURX - dfLLX) / nRasterXSize;
1473 0 : const double dfPixelY = (dfURY - dfLLY) / nRasterYSize;
1474 :
1475 0 : bHasGeoTransform = true;
1476 0 : m_gt.xorig = dfLLX;
1477 0 : m_gt.xscale = dfPixelX;
1478 0 : m_gt.xrot = 0;
1479 0 : m_gt.yorig = dfURY;
1480 0 : m_gt.yrot = 0;
1481 0 : m_gt.yscale = -dfPixelY;
1482 :
1483 0 : return CE_None;
1484 : }
1485 :
1486 : /************************************************************************/
1487 : /* CreateProjections() */
1488 : /************************************************************************/
1489 48 : CPLErr HDF5ImageDataset::CreateProjections()
1490 : {
1491 48 : switch (iSubdatasetType)
1492 : {
1493 2 : case CSK_PRODUCT:
1494 : {
1495 2 : int productType = PROD_UNKNOWN;
1496 :
1497 2 : if (GetMetadataItem("Product_Type") != nullptr)
1498 : {
1499 : // Get the format's level.
1500 : const char *osMissionLevel =
1501 2 : HDF5Dataset::GetMetadataItem("Product_Type");
1502 :
1503 2 : if (STARTS_WITH_CI(osMissionLevel, "RAW"))
1504 0 : productType = PROD_CSK_L0;
1505 :
1506 2 : if (STARTS_WITH_CI(osMissionLevel, "SSC"))
1507 0 : productType = PROD_CSK_L1A;
1508 :
1509 2 : if (STARTS_WITH_CI(osMissionLevel, "DGM"))
1510 1 : productType = PROD_CSK_L1B;
1511 :
1512 2 : if (STARTS_WITH_CI(osMissionLevel, "GEC"))
1513 1 : productType = PROD_CSK_L1C;
1514 :
1515 2 : if (STARTS_WITH_CI(osMissionLevel, "GTC"))
1516 0 : productType = PROD_CSK_L1D;
1517 : }
1518 :
1519 2 : CaptureCSKGeoTransform(productType);
1520 2 : CaptureCSKGeolocation(productType);
1521 2 : CaptureCSKGCPs(productType);
1522 :
1523 2 : break;
1524 : }
1525 46 : case UNKNOWN_PRODUCT:
1526 : {
1527 46 : if (ndims == 2)
1528 : {
1529 : // Check if there's a CF conventions grid_mapping attribute.
1530 : const char *pszGridMapping =
1531 36 : GetRasterBand(1)->GetMetadataItem("grid_mapping");
1532 36 : if (pszGridMapping)
1533 : {
1534 : const std::string osGroupName =
1535 18 : CPLGetDirnameSafe(m_osPath.c_str());
1536 : HDF5GroupObjects *poH5Objects =
1537 9 : HDF5FindDatasetObjectsbyPath(
1538 : poH5RootGroup,
1539 18 : (osGroupName + "/" + pszGridMapping).c_str());
1540 9 : if (poH5Objects && poH5Objects->nType == H5G_DATASET)
1541 : {
1542 18 : CPLStringList aosMetadata;
1543 9 : CreateMetadata(m_hHDF5, poH5Objects, H5G_DATASET, false,
1544 : aosMetadata);
1545 :
1546 : // NISAR specific
1547 9 : const int nEPSGCode = atoi(
1548 : aosMetadata.FetchNameValueDef("epsg_code", "0"));
1549 9 : if (nEPSGCode > 0)
1550 : {
1551 1 : if (m_oSRS.importFromEPSG(nEPSGCode) != OGRERR_NONE)
1552 : {
1553 0 : m_oSRS.Clear();
1554 : }
1555 : }
1556 : else
1557 : {
1558 8 : if (m_oSRS.importFromCF1(
1559 8 : aosMetadata.List(),
1560 8 : /* pszUnits = */ nullptr) != OGRERR_NONE)
1561 : {
1562 0 : m_oSRS.Clear();
1563 : }
1564 : }
1565 : }
1566 : }
1567 :
1568 : // NISAR Level 2 specific
1569 36 : const char *pszMissionName = GetMetadataItem("mission_name");
1570 36 : if (pszMissionName && EQUAL(pszMissionName, "NISAR"))
1571 : {
1572 : const std::string osGroupName =
1573 4 : CPLGetDirnameSafe(m_osPath.c_str());
1574 : HDF5GroupObjects *poXCoordinates =
1575 2 : HDF5FindDatasetObjectsbyPath(
1576 : poH5RootGroup,
1577 4 : (osGroupName + "/xCoordinates").c_str());
1578 : HDF5GroupObjects *poYCoordinates =
1579 2 : HDF5FindDatasetObjectsbyPath(
1580 : poH5RootGroup,
1581 4 : (osGroupName + "/yCoordinates").c_str());
1582 2 : if (poXCoordinates && poYCoordinates &&
1583 2 : poXCoordinates->nType == H5G_DATASET &&
1584 2 : poYCoordinates->nType == H5G_DATASET &&
1585 2 : poXCoordinates->nRank == 1 &&
1586 2 : poYCoordinates->nRank == 1 &&
1587 2 : H5Tequal(H5T_NATIVE_DOUBLE, poXCoordinates->native) &&
1588 2 : H5Tequal(H5T_NATIVE_DOUBLE, poYCoordinates->native) &&
1589 6 : poXCoordinates->paDims[0] == dims[1] &&
1590 2 : poYCoordinates->paDims[0] == dims[0])
1591 : {
1592 2 : hid_t hDSX = H5Dopen(m_hHDF5, poXCoordinates->pszPath);
1593 2 : hid_t hDSY = H5Dopen(m_hHDF5, poYCoordinates->pszPath);
1594 :
1595 2 : if (hDSX >= 0 && hDSY >= 0)
1596 : {
1597 : const auto Read1DDoubleArray =
1598 4 : [](hid_t hDS, hsize_t count,
1599 : std::vector<double> &v)
1600 : {
1601 4 : hsize_t anCount[1] = {count};
1602 4 : hsize_t anStart[1] = {0};
1603 4 : v.resize(static_cast<size_t>(count));
1604 : hid_t hMem =
1605 4 : H5Screate_simple(1, anCount, nullptr);
1606 4 : const hid_t hDataspace = H5Dget_space(hDS);
1607 4 : H5Sselect_hyperslab(hDataspace, H5S_SELECT_SET,
1608 : anStart, nullptr, anCount,
1609 : nullptr);
1610 4 : bool ret = H5Dread(hDS, H5T_NATIVE_DOUBLE, hMem,
1611 : hDataspace, H5P_DEFAULT,
1612 4 : v.data()) >= 0;
1613 4 : H5Sclose(hDataspace);
1614 4 : H5Sclose(hMem);
1615 4 : return ret;
1616 : };
1617 :
1618 4 : std::vector<double> adfX, adfY;
1619 6 : if (Read1DDoubleArray(
1620 2 : hDSX, poXCoordinates->paDims[0], adfX) &&
1621 2 : Read1DDoubleArray(
1622 4 : hDSY, poYCoordinates->paDims[0], adfY) &&
1623 6 : adfX.size() >= 3 && adfY.size() >= 3)
1624 : {
1625 2 : constexpr double EPSILON = 1e-3;
1626 2 : double dfXSpacing = adfX[1] - adfX[0];
1627 6 : for (size_t i = 1; i + 1 < adfX.size(); ++i)
1628 : {
1629 : const double dfDelta =
1630 4 : adfX[i + 1] - adfX[i];
1631 4 : if (!(std::fabs(dfDelta - dfXSpacing) <=
1632 4 : EPSILON * std::fabs(dfXSpacing)))
1633 : {
1634 : dfXSpacing = std::numeric_limits<
1635 0 : double>::quiet_NaN();
1636 0 : break;
1637 : }
1638 : }
1639 :
1640 2 : double dfYSpacing = adfY[1] - adfY[0];
1641 8 : for (size_t i = 1; i + 1 < adfY.size(); ++i)
1642 : {
1643 : const double dfDelta =
1644 6 : adfY[i + 1] - adfY[i];
1645 6 : if (!(std::fabs(dfDelta - dfYSpacing) <=
1646 6 : EPSILON * std::fabs(dfYSpacing)))
1647 : {
1648 : dfYSpacing = std::numeric_limits<
1649 0 : double>::quiet_NaN();
1650 0 : break;
1651 : }
1652 : }
1653 :
1654 4 : if (!std::isnan(dfXSpacing) &&
1655 2 : !std::isnan(dfYSpacing))
1656 : {
1657 2 : bHasGeoTransform = true;
1658 2 : m_gt.xorig = adfX[0] - dfXSpacing / 2;
1659 2 : m_gt.xscale = dfXSpacing;
1660 2 : m_gt.yorig = adfY[0] - dfYSpacing / 2;
1661 2 : m_gt.yscale = dfYSpacing;
1662 : }
1663 : }
1664 : }
1665 2 : if (hDSX >= 0)
1666 2 : H5Dclose(hDSX);
1667 2 : if (hDSY >= 0)
1668 2 : H5Dclose(hDSY);
1669 : }
1670 : }
1671 : }
1672 :
1673 46 : constexpr int NBGCPLAT = 100;
1674 46 : constexpr int NBGCPLON = 30;
1675 :
1676 46 : const int nDeltaLat = nRasterYSize / NBGCPLAT;
1677 46 : const int nDeltaLon = nRasterXSize / NBGCPLON;
1678 :
1679 46 : if (nDeltaLat == 0 || nDeltaLon == 0)
1680 38 : return CE_None;
1681 :
1682 : // Create HDF5 Data Hierarchy in a link list.
1683 : HDF5GroupObjects *poH5Objects =
1684 8 : HDF5FindDatasetObjects(poH5RootGroup, "Latitude");
1685 8 : if (!poH5Objects)
1686 : {
1687 8 : if (GetMetadataItem("where_projdef") != nullptr)
1688 0 : return CreateODIMH5Projection();
1689 8 : return CE_None;
1690 : }
1691 :
1692 : // The Latitude and Longitude arrays must have a rank of 2 to
1693 : // retrieve GCPs.
1694 0 : if (poH5Objects->nRank != 2 ||
1695 0 : poH5Objects->paDims[0] != static_cast<size_t>(nRasterYSize) ||
1696 0 : poH5Objects->paDims[1] != static_cast<size_t>(nRasterXSize))
1697 : {
1698 0 : return CE_None;
1699 : }
1700 :
1701 : // Retrieve HDF5 data information.
1702 : const hid_t LatitudeDatasetID =
1703 0 : H5Dopen(m_hHDF5, poH5Objects->pszPath);
1704 : // LatitudeDataspaceID = H5Dget_space(dataset_id);
1705 :
1706 0 : poH5Objects = HDF5FindDatasetObjects(poH5RootGroup, "Longitude");
1707 : // GCPs.
1708 0 : if (poH5Objects == nullptr || poH5Objects->nRank != 2 ||
1709 0 : poH5Objects->paDims[0] != static_cast<size_t>(nRasterYSize) ||
1710 0 : poH5Objects->paDims[1] != static_cast<size_t>(nRasterXSize))
1711 : {
1712 0 : if (LatitudeDatasetID > 0)
1713 0 : H5Dclose(LatitudeDatasetID);
1714 0 : return CE_None;
1715 : }
1716 :
1717 : const hid_t LongitudeDatasetID =
1718 0 : H5Dopen(m_hHDF5, poH5Objects->pszPath);
1719 : // LongitudeDataspaceID = H5Dget_space(dataset_id);
1720 :
1721 0 : if (LatitudeDatasetID > 0 && LongitudeDatasetID > 0)
1722 : {
1723 : float *const Latitude =
1724 0 : static_cast<float *>(VSI_MALLOC3_VERBOSE(
1725 : nRasterYSize, nRasterXSize, sizeof(float)));
1726 : float *const Longitude =
1727 0 : static_cast<float *>(VSI_MALLOC3_VERBOSE(
1728 : nRasterYSize, nRasterXSize, sizeof(float)));
1729 0 : if (!Latitude || !Longitude)
1730 : {
1731 0 : CPLFree(Latitude);
1732 0 : CPLFree(Longitude);
1733 0 : H5Dclose(LatitudeDatasetID);
1734 0 : H5Dclose(LongitudeDatasetID);
1735 0 : return CE_Failure;
1736 : }
1737 0 : memset(Latitude, 0,
1738 0 : static_cast<size_t>(nRasterXSize) * nRasterYSize *
1739 : sizeof(float));
1740 0 : memset(Longitude, 0,
1741 0 : static_cast<size_t>(nRasterXSize) * nRasterYSize *
1742 : sizeof(float));
1743 :
1744 : // netCDF convention for nodata
1745 0 : double dfLatNoData = 0;
1746 0 : bool bHasLatNoData = GH5_FetchAttribute(
1747 : LatitudeDatasetID, "_FillValue", dfLatNoData);
1748 :
1749 0 : double dfLongNoData = 0;
1750 0 : bool bHasLongNoData = GH5_FetchAttribute(
1751 : LongitudeDatasetID, "_FillValue", dfLongNoData);
1752 :
1753 0 : H5Dread(LatitudeDatasetID, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
1754 : H5P_DEFAULT, Latitude);
1755 :
1756 0 : H5Dread(LongitudeDatasetID, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
1757 : H5P_DEFAULT, Longitude);
1758 :
1759 0 : m_oSRS.Clear();
1760 0 : m_oGCPSRS.SetWellKnownGeogCS("WGS84");
1761 :
1762 0 : const int nYLimit =
1763 0 : (static_cast<int>(nRasterYSize) / nDeltaLat) * nDeltaLat;
1764 0 : const int nXLimit =
1765 0 : (static_cast<int>(nRasterXSize) / nDeltaLon) * nDeltaLon;
1766 :
1767 : // The original code in
1768 : // https://trac.osgeo.org/gdal/changeset/8066 always add +180 to
1769 : // the longitudes, but without justification I suspect this
1770 : // might be due to handling products crossing the antimeridian.
1771 : // Trying to do it just when needed through a heuristics.
1772 0 : bool bHasLonNearMinus180 = false;
1773 0 : bool bHasLonNearPlus180 = false;
1774 0 : bool bHasLonNearZero = false;
1775 0 : for (int j = 0; j < nYLimit; j += nDeltaLat)
1776 : {
1777 0 : for (int i = 0; i < nXLimit; i += nDeltaLon)
1778 : {
1779 0 : const int iGCP = j * nRasterXSize + i;
1780 0 : if ((bHasLatNoData && static_cast<float>(dfLatNoData) ==
1781 0 : Latitude[iGCP]) ||
1782 0 : (bHasLongNoData &&
1783 0 : static_cast<float>(dfLongNoData) ==
1784 0 : Longitude[iGCP]))
1785 0 : continue;
1786 0 : if (Longitude[iGCP] > 170 && Longitude[iGCP] <= 180)
1787 0 : bHasLonNearPlus180 = true;
1788 0 : if (Longitude[iGCP] < -170 && Longitude[iGCP] >= -180)
1789 0 : bHasLonNearMinus180 = true;
1790 0 : if (fabs(Longitude[iGCP]) < 90)
1791 0 : bHasLonNearZero = true;
1792 : }
1793 : }
1794 :
1795 : // Fill the GCPs list.
1796 : const char *pszShiftGCP =
1797 0 : CPLGetConfigOption("HDF5_SHIFT_GCPX_BY_180", nullptr);
1798 : const bool bAdd180 =
1799 0 : (bHasLonNearPlus180 && bHasLonNearMinus180 &&
1800 0 : !bHasLonNearZero && pszShiftGCP == nullptr) ||
1801 0 : (pszShiftGCP != nullptr && CPLTestBool(pszShiftGCP));
1802 :
1803 0 : for (int j = 0; j < nYLimit; j += nDeltaLat)
1804 : {
1805 0 : for (int i = 0; i < nXLimit; i += nDeltaLon)
1806 : {
1807 0 : const int iGCP = j * nRasterXSize + i;
1808 0 : if ((bHasLatNoData && static_cast<float>(dfLatNoData) ==
1809 0 : Latitude[iGCP]) ||
1810 0 : (bHasLongNoData &&
1811 0 : static_cast<float>(dfLongNoData) ==
1812 0 : Longitude[iGCP]))
1813 0 : continue;
1814 0 : double dfGCPX = static_cast<double>(Longitude[iGCP]);
1815 0 : if (bAdd180)
1816 0 : dfGCPX += 180.0;
1817 0 : const double dfGCPY =
1818 0 : static_cast<double>(Latitude[iGCP]);
1819 :
1820 0 : m_aoGCPs.emplace_back("", "", i + 0.5, j + 0.5, dfGCPX,
1821 0 : dfGCPY);
1822 : }
1823 : }
1824 :
1825 0 : CPLFree(Latitude);
1826 0 : CPLFree(Longitude);
1827 : }
1828 :
1829 0 : if (LatitudeDatasetID > 0)
1830 0 : H5Dclose(LatitudeDatasetID);
1831 0 : if (LongitudeDatasetID > 0)
1832 0 : H5Dclose(LongitudeDatasetID);
1833 :
1834 0 : break;
1835 : }
1836 : }
1837 :
1838 2 : return CE_None;
1839 : }
1840 :
1841 : /************************************************************************/
1842 : /* GetMetadataItem() */
1843 : /************************************************************************/
1844 :
1845 160 : const char *HDF5ImageDataset::GetMetadataItem(const char *pszName,
1846 : const char *pszDomain)
1847 : {
1848 160 : if (pszDomain && EQUAL(pszDomain, "__DEBUG__") &&
1849 66 : EQUAL(pszName, "WholeBandChunkOptim"))
1850 : {
1851 66 : switch (m_eWholeBandChunkOptim)
1852 : {
1853 4 : case WBC_DETECTION_IN_PROGRESS:
1854 4 : return "DETECTION_IN_PROGRESS";
1855 3 : case WBC_DISABLED:
1856 3 : return "DISABLED";
1857 59 : case WBC_ENABLED:
1858 59 : return "ENABLED";
1859 : }
1860 : }
1861 94 : return GDALPamDataset::GetMetadataItem(pszName, pszDomain);
1862 : }
1863 :
1864 : /************************************************************************/
1865 : /* GetSpatialRef() */
1866 : /************************************************************************/
1867 :
1868 14 : const OGRSpatialReference *HDF5ImageDataset::GetSpatialRef() const
1869 : {
1870 14 : if (!m_oSRS.IsEmpty())
1871 13 : return &m_oSRS;
1872 1 : return GDALPamDataset::GetSpatialRef();
1873 : }
1874 :
1875 : /************************************************************************/
1876 : /* GetGCPCount() */
1877 : /************************************************************************/
1878 :
1879 4 : int HDF5ImageDataset::GetGCPCount()
1880 :
1881 : {
1882 4 : if (!m_aoGCPs.empty())
1883 1 : return static_cast<int>(m_aoGCPs.size());
1884 :
1885 3 : return GDALPamDataset::GetGCPCount();
1886 : }
1887 :
1888 : /************************************************************************/
1889 : /* GetGCPSpatialRef() */
1890 : /************************************************************************/
1891 :
1892 1 : const OGRSpatialReference *HDF5ImageDataset::GetGCPSpatialRef() const
1893 :
1894 : {
1895 1 : if (!m_aoGCPs.empty() && !m_oGCPSRS.IsEmpty())
1896 1 : return &m_oGCPSRS;
1897 :
1898 0 : return GDALPamDataset::GetGCPSpatialRef();
1899 : }
1900 :
1901 : /************************************************************************/
1902 : /* GetGCPs() */
1903 : /************************************************************************/
1904 :
1905 3 : const GDAL_GCP *HDF5ImageDataset::GetGCPs()
1906 : {
1907 3 : if (!m_aoGCPs.empty())
1908 1 : return gdal::GCP::c_ptr(m_aoGCPs);
1909 :
1910 2 : return GDALPamDataset::GetGCPs();
1911 : }
1912 :
1913 : /************************************************************************/
1914 : /* GetGeoTransform() */
1915 : /************************************************************************/
1916 :
1917 9 : CPLErr HDF5ImageDataset::GetGeoTransform(GDALGeoTransform >) const
1918 : {
1919 9 : if (bHasGeoTransform)
1920 : {
1921 8 : gt = m_gt;
1922 8 : return CE_None;
1923 : }
1924 :
1925 1 : return GDALPamDataset::GetGeoTransform(gt);
1926 : }
1927 :
1928 : /************************************************************************/
1929 : /* IdentifyProductType() */
1930 : /************************************************************************/
1931 :
1932 : /**
1933 : * Identify if the subdataset has a known product format
1934 : * It stores a product identifier in iSubdatasetType,
1935 : * UNKNOWN_PRODUCT, if it isn't a recognizable format.
1936 : */
1937 62 : void HDF5ImageDataset::IdentifyProductType()
1938 : {
1939 62 : iSubdatasetType = UNKNOWN_PRODUCT;
1940 :
1941 : // COSMO-SKYMED
1942 :
1943 : // Get the Mission Id as a char *, because the field may not exist.
1944 62 : const char *const pszMissionId = HDF5Dataset::GetMetadataItem("Mission_ID");
1945 :
1946 : // If there is a Mission_ID field.
1947 62 : if (pszMissionId != nullptr && strstr(GetDescription(), "QLK") == nullptr)
1948 : {
1949 : // Check if the mission type is CSK, KMPS or CSG.
1950 : // KMPS: Komsat-5 is Korean mission with a SAR instrument.
1951 : // CSG: Cosmo Skymed 2nd Generation
1952 2 : if (EQUAL(pszMissionId, "CSK") || EQUAL(pszMissionId, "KMPS") ||
1953 0 : EQUAL(pszMissionId, "CSG"))
1954 : {
1955 2 : iSubdatasetType = CSK_PRODUCT;
1956 :
1957 2 : if (GetMetadataItem("Product_Type") != nullptr)
1958 : {
1959 : // Get the format's level.
1960 : const char *osMissionLevel =
1961 2 : HDF5Dataset::GetMetadataItem("Product_Type");
1962 :
1963 2 : if (STARTS_WITH_CI(osMissionLevel, "RAW"))
1964 0 : iCSKProductType = PROD_CSK_L0;
1965 :
1966 2 : if (STARTS_WITH_CI(osMissionLevel, "SCS"))
1967 0 : iCSKProductType = PROD_CSK_L1A;
1968 :
1969 2 : if (STARTS_WITH_CI(osMissionLevel, "DGM"))
1970 1 : iCSKProductType = PROD_CSK_L1B;
1971 :
1972 2 : if (STARTS_WITH_CI(osMissionLevel, "GEC"))
1973 1 : iCSKProductType = PROD_CSK_L1C;
1974 :
1975 2 : if (STARTS_WITH_CI(osMissionLevel, "GTC"))
1976 0 : iCSKProductType = PROD_CSK_L1D;
1977 : }
1978 : }
1979 : }
1980 62 : }
1981 :
1982 : /************************************************************************/
1983 : /* CaptureCSKGeolocation() */
1984 : /************************************************************************/
1985 :
1986 : /**
1987 : * Captures Geolocation information from a COSMO-SKYMED
1988 : * file.
1989 : * The geoid will always be WGS84
1990 : * The projection type may be UTM or UPS, depending on the
1991 : * latitude from the center of the image.
1992 : * @param iProductType type of CSK subproduct, see HDF5CSKProduct
1993 : */
1994 2 : void HDF5ImageDataset::CaptureCSKGeolocation(int iProductType)
1995 : {
1996 : // Set the ellipsoid to WGS84.
1997 2 : m_oSRS.SetWellKnownGeogCS("WGS84");
1998 :
1999 2 : if (iProductType == PROD_CSK_L1C || iProductType == PROD_CSK_L1D)
2000 : {
2001 1 : double *dfProjFalseEastNorth = nullptr;
2002 1 : double *dfProjScaleFactor = nullptr;
2003 1 : double *dfCenterCoord = nullptr;
2004 :
2005 : // Check if all the metadata attributes are present.
2006 1 : if (HDF5ReadDoubleAttr("Map Projection False East-North",
2007 1 : &dfProjFalseEastNorth) == CE_Failure ||
2008 1 : HDF5ReadDoubleAttr("Map Projection Scale Factor",
2009 1 : &dfProjScaleFactor) == CE_Failure ||
2010 1 : HDF5ReadDoubleAttr("Map Projection Centre", &dfCenterCoord) ==
2011 2 : CE_Failure ||
2012 1 : GetMetadataItem("Projection_ID") == nullptr)
2013 : {
2014 0 : m_oSRS.Clear();
2015 0 : m_oGCPSRS.Clear();
2016 0 : CPLError(CE_Failure, CPLE_OpenFailed,
2017 : "The CSK hdf5 file geolocation information is "
2018 : "malformed");
2019 : }
2020 : else
2021 : {
2022 : // Fetch projection Type.
2023 2 : CPLString osProjectionID = GetMetadataItem("Projection_ID");
2024 :
2025 : // If the projection is UTM.
2026 1 : if (EQUAL(osProjectionID, "UTM"))
2027 : {
2028 : // @TODO: use SetUTM
2029 1 : m_oSRS.SetProjCS(SRS_PT_TRANSVERSE_MERCATOR);
2030 1 : m_oSRS.SetTM(dfCenterCoord[0], dfCenterCoord[1],
2031 : dfProjScaleFactor[0], dfProjFalseEastNorth[0],
2032 1 : dfProjFalseEastNorth[1]);
2033 : }
2034 : else
2035 : {
2036 : // TODO: No UPS projected files to test!
2037 : // If the projection is UPS.
2038 0 : if (EQUAL(osProjectionID, "UPS"))
2039 : {
2040 0 : m_oSRS.SetProjCS(SRS_PT_POLAR_STEREOGRAPHIC);
2041 0 : m_oSRS.SetPS(dfCenterCoord[0], dfCenterCoord[1],
2042 : dfProjScaleFactor[0], dfProjFalseEastNorth[0],
2043 0 : dfProjFalseEastNorth[1]);
2044 : }
2045 : }
2046 :
2047 1 : CPLFree(dfCenterCoord);
2048 1 : CPLFree(dfProjScaleFactor);
2049 1 : CPLFree(dfProjFalseEastNorth);
2050 1 : }
2051 : }
2052 : else
2053 : {
2054 1 : m_oGCPSRS = m_oSRS;
2055 : }
2056 2 : }
2057 :
2058 : /************************************************************************/
2059 : /* CaptureCSKGeoTransform() */
2060 : /************************************************************************/
2061 :
2062 : /**
2063 : * Get Geotransform information for COSMO-SKYMED files
2064 : * In case of success it stores the transformation
2065 : * in m_gt. In case of failure it doesn't
2066 : * modify m_gt
2067 : * @param iProductType type of CSK subproduct, see HDF5CSKProduct
2068 : */
2069 2 : void HDF5ImageDataset::CaptureCSKGeoTransform(int iProductType)
2070 : {
2071 2 : const char *pszSubdatasetName = GetSubdatasetName();
2072 :
2073 2 : bHasGeoTransform = false;
2074 : // If the product level is not L1C or L1D then
2075 : // it doesn't have a valid projection.
2076 2 : if (iProductType == PROD_CSK_L1C || iProductType == PROD_CSK_L1D)
2077 : {
2078 : // If there is a subdataset.
2079 1 : if (pszSubdatasetName != nullptr)
2080 : {
2081 2 : CPLString osULPath = pszSubdatasetName;
2082 1 : osULPath += "/Top Left East-North";
2083 :
2084 2 : CPLString osLineSpacingPath = pszSubdatasetName;
2085 1 : osLineSpacingPath += "/Line Spacing";
2086 :
2087 2 : CPLString osColumnSpacingPath = pszSubdatasetName;
2088 1 : osColumnSpacingPath += "/Column Spacing";
2089 :
2090 1 : double *pdOutUL = nullptr;
2091 1 : double *pdLineSpacing = nullptr;
2092 1 : double *pdColumnSpacing = nullptr;
2093 :
2094 : // If it could find the attributes on the metadata.
2095 1 : if (HDF5ReadDoubleAttr(osULPath.c_str(), &pdOutUL) == CE_Failure ||
2096 1 : HDF5ReadDoubleAttr(osLineSpacingPath.c_str(), &pdLineSpacing) ==
2097 2 : CE_Failure ||
2098 1 : HDF5ReadDoubleAttr(osColumnSpacingPath.c_str(),
2099 : &pdColumnSpacing) == CE_Failure)
2100 : {
2101 0 : bHasGeoTransform = false;
2102 : }
2103 : else
2104 : {
2105 : // geotransform[1] : width of pixel
2106 : // geotransform[4] : rotational coefficient, zero for north up
2107 : // images.
2108 : // geotransform[2] : rotational coefficient, zero for north up
2109 : // images.
2110 : // geotransform[5] : height of pixel (but negative)
2111 :
2112 1 : m_gt.xorig = pdOutUL[0];
2113 1 : m_gt.xscale = pdLineSpacing[0];
2114 1 : m_gt.xrot = 0;
2115 1 : m_gt.yorig = pdOutUL[1];
2116 1 : m_gt.yrot = 0;
2117 1 : m_gt.yscale = -pdColumnSpacing[0];
2118 :
2119 1 : CPLFree(pdOutUL);
2120 1 : CPLFree(pdLineSpacing);
2121 1 : CPLFree(pdColumnSpacing);
2122 :
2123 1 : bHasGeoTransform = true;
2124 : }
2125 : }
2126 : }
2127 2 : }
2128 :
2129 : /************************************************************************/
2130 : /* CaptureCSKGCPs() */
2131 : /************************************************************************/
2132 :
2133 : /**
2134 : * Retrieves and stores the GCPs from a COSMO-SKYMED dataset
2135 : * It only retrieves the GCPs for L0, L1A and L1B products
2136 : * for L1C and L1D products, geotransform is provided.
2137 : * The GCPs provided will be the Image's corners.
2138 : * @param iProductType type of CSK product @see HDF5CSKProductEnum
2139 : */
2140 2 : void HDF5ImageDataset::CaptureCSKGCPs(int iProductType)
2141 : {
2142 : // Only retrieve GCPs for L0,L1A and L1B products.
2143 2 : if (iProductType == PROD_CSK_L0 || iProductType == PROD_CSK_L1A ||
2144 : iProductType == PROD_CSK_L1B)
2145 : {
2146 10 : CPLString osCornerName[4];
2147 1 : double pdCornerPixel[4] = {0.0, 0.0, 0.0, 0.0};
2148 1 : double pdCornerLine[4] = {0.0, 0.0, 0.0, 0.0};
2149 :
2150 1 : const char *const pszSubdatasetName = GetSubdatasetName();
2151 :
2152 : // Load the subdataset name first.
2153 5 : for (int i = 0; i < 4; i++)
2154 4 : osCornerName[i] = pszSubdatasetName;
2155 :
2156 : // Load the attribute name, and raster coordinates for
2157 : // all the corners.
2158 1 : osCornerName[0] += "/Top Left Geodetic Coordinates";
2159 1 : pdCornerPixel[0] = 0;
2160 1 : pdCornerLine[0] = 0;
2161 :
2162 1 : osCornerName[1] += "/Top Right Geodetic Coordinates";
2163 1 : pdCornerPixel[1] = GetRasterXSize();
2164 1 : pdCornerLine[1] = 0;
2165 :
2166 1 : osCornerName[2] += "/Bottom Left Geodetic Coordinates";
2167 1 : pdCornerPixel[2] = 0;
2168 1 : pdCornerLine[2] = GetRasterYSize();
2169 :
2170 1 : osCornerName[3] += "/Bottom Right Geodetic Coordinates";
2171 1 : pdCornerPixel[3] = GetRasterXSize();
2172 1 : pdCornerLine[3] = GetRasterYSize();
2173 :
2174 : // For all the image's corners.
2175 5 : for (int i = 0; i < 4; i++)
2176 : {
2177 4 : double *pdCornerCoordinates = nullptr;
2178 :
2179 : // Retrieve the attributes.
2180 4 : if (HDF5ReadDoubleAttr(osCornerName[i].c_str(),
2181 4 : &pdCornerCoordinates) == CE_Failure)
2182 : {
2183 0 : CPLError(CE_Failure, CPLE_OpenFailed,
2184 : "Error retrieving CSK GCPs");
2185 0 : m_aoGCPs.clear();
2186 0 : break;
2187 : }
2188 :
2189 0 : m_aoGCPs.emplace_back(osCornerName[i].c_str(), "", pdCornerPixel[i],
2190 4 : pdCornerLine[i],
2191 4 : /* X = */ pdCornerCoordinates[1],
2192 : /* Y = */ pdCornerCoordinates[0],
2193 4 : /* Z = */ pdCornerCoordinates[2]);
2194 :
2195 : // Free the returned coordinates.
2196 4 : CPLFree(pdCornerCoordinates);
2197 : }
2198 : }
2199 2 : }
|
