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