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 644 : int GetYIndex() const
150 : {
151 644 : return m_nYIndex;
152 : }
153 :
154 778 : int GetXIndex() const
155 : {
156 778 : 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 90 : 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 134 : const auto IsConsecutiveBands = [](const int *panVals, int nCount)
802 : {
803 151 : for (int i = 1; i < nCount; ++i)
804 : {
805 19 : if (panVals[i] != panVals[i - 1] + 1)
806 2 : return false;
807 : }
808 132 : return true;
809 : };
810 :
811 90 : const auto eDT = GetRasterBand(1)->GetRasterDataType();
812 90 : const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
813 :
814 : // Band-interleaved data and request
815 176 : const bool bIsBandInterleavedData = ndims == 3 && m_nOtherDimIndex == 0 &&
816 266 : GetYIndex() == 1 && GetXIndex() == 2;
817 90 : if (eRWFlag == GF_Read && bIsBandInterleavedData && nXSize == nBufXSize &&
818 86 : nYSize == nBufYSize && IsConsecutiveBands(panBandMap, nBandCount) &&
819 84 : eBufType == eDT && nPixelSpace == nDTSize &&
820 180 : nLineSpace == nXSize * nPixelSpace && nBandSpace == nYSize * nLineSpace)
821 : {
822 : HDF5_GLOBAL_LOCK();
823 :
824 43 : hsize_t count[3] = {static_cast<hsize_t>(nBandCount),
825 43 : static_cast<hsize_t>(nYSize),
826 43 : static_cast<hsize_t>(nXSize)};
827 : H5OFFSET_TYPE offset[3] = {
828 43 : static_cast<H5OFFSET_TYPE>(panBandMap[0] - 1),
829 43 : static_cast<H5OFFSET_TYPE>(nYOff),
830 43 : static_cast<H5OFFSET_TYPE>(nXOff)};
831 43 : herr_t status = H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET,
832 : offset, nullptr, count, nullptr);
833 43 : if (status < 0)
834 0 : return CE_Failure;
835 :
836 43 : const hid_t memspace = H5Screate_simple(ndims, count, nullptr);
837 43 : H5OFFSET_TYPE mem_offset[3] = {0, 0, 0};
838 43 : status = H5Sselect_hyperslab(memspace, H5S_SELECT_SET, mem_offset,
839 : nullptr, count, nullptr);
840 43 : if (status < 0)
841 : {
842 0 : H5Sclose(memspace);
843 0 : return CE_Failure;
844 : }
845 :
846 43 : status = H5Dread(dataset_id, native, memspace, dataspace_id,
847 : H5P_DEFAULT, pData);
848 :
849 43 : H5Sclose(memspace);
850 :
851 43 : if (status < 0)
852 : {
853 0 : CPLError(CE_Failure, CPLE_AppDefined,
854 : "HDF5ImageDataset::IRasterIO(): H5Dread() failed");
855 0 : return CE_Failure;
856 : }
857 :
858 43 : return CE_None;
859 : }
860 :
861 : // Pixel-interleaved data and request
862 :
863 51 : const bool bIsPixelInterleaveData = ndims == 3 && m_nOtherDimIndex == 2 &&
864 98 : GetYIndex() == 0 && GetXIndex() == 1;
865 47 : if (eRWFlag == GF_Read && bIsPixelInterleaveData && nXSize == nBufXSize &&
866 4 : nYSize == nBufYSize && IsConsecutiveBands(panBandMap, nBandCount) &&
867 4 : eBufType == eDT && nBandSpace == nDTSize &&
868 97 : 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 44 : if (eRWFlag == GF_Read &&
914 44 : (bIsBandInterleavedData || bIsPixelInterleaveData) &&
915 88 : nXSize == nBufXSize && nYSize == nBufYSize &&
916 132 : IsConsecutiveBands(panBandMap, nBandCount) &&
917 43 : static_cast<GIntBig>(nXSize) * nYSize <
918 43 : CPLGetUsablePhysicalRAM() / 10 / nBandCount)
919 : {
920 43 : const char *const apszOptions[] = {
921 43 : bIsPixelInterleaveData ? "INTERLEAVE=PIXEL" : nullptr, nullptr};
922 : auto poMemDS = std::unique_ptr<GDALDataset>(
923 43 : MEMDataset::Create("", nXSize, nYSize, nBandCount, eDT,
924 43 : const_cast<char **>(apszOptions)));
925 43 : if (poMemDS)
926 : {
927 43 : void *pMemData = poMemDS->GetInternalHandle("MEMORY1");
928 43 : CPLAssert(pMemData);
929 : // Read from HDF5 into the temporary MEMDataset using the
930 : // natural interleaving of the HDF5 dataset
931 129 : if (IRasterIO(
932 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pMemData, nXSize,
933 : nYSize, eDT, nBandCount, panBandMap,
934 1 : bIsBandInterleavedData ? nDTSize : nDTSize * nBandCount,
935 : bIsBandInterleavedData
936 42 : ? static_cast<GSpacing>(nXSize) * nDTSize
937 1 : : static_cast<GSpacing>(nXSize) * nDTSize * nBandCount,
938 : bIsBandInterleavedData
939 42 : ? static_cast<GSpacing>(nYSize) * nXSize * nDTSize
940 : : nDTSize,
941 43 : psExtraArg) != CE_None)
942 : {
943 0 : return CE_Failure;
944 : }
945 : // Copy to the final buffer using requested data type and spacings.
946 43 : return poMemDS->RasterIO(GF_Read, 0, 0, nXSize, nYSize, pData,
947 : nXSize, nYSize, eBufType, nBandCount,
948 : nullptr, nPixelSpace, nLineSpace,
949 43 : 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 6 : std::shared_ptr<GDALMDArray> poXDim, poYDim;
1127 20 : for (const auto &oDim : oGridDataFieldMetadata.aoDimensions)
1128 : {
1129 14 : if (oDim.osName == "XDim")
1130 6 : poDS->m_nXIndex = iDim;
1131 8 : else if (oDim.osName == "YDim")
1132 6 : poDS->m_nYIndex = iDim;
1133 : else
1134 2 : poDS->m_nOtherDimIndex = iDim;
1135 14 : ++iDim;
1136 : }
1137 :
1138 6 : if (oGridDataFieldMetadata.poGridMetadata->GetGeoTransform(
1139 6 : poDS->m_gt))
1140 6 : poDS->bHasGeoTransform = true;
1141 :
1142 12 : auto poSRS = oGridDataFieldMetadata.poGridMetadata->GetSRS();
1143 6 : if (poSRS)
1144 6 : poDS->m_oSRS = *(poSRS.get());
1145 : }
1146 13 : else if (oHDFEOSParser.GetDataModel() ==
1147 13 : HDF5EOSParser::DataModel::SWATH &&
1148 13 : oHDFEOSParser.GetSwathDataFieldMetadata(
1149 7 : osSubdatasetName.c_str(), oSwathDataFieldMetadata) &&
1150 : static_cast<int>(
1151 7 : oSwathDataFieldMetadata.aoDimensions.size()) ==
1152 7 : poDS->ndims &&
1153 33 : oSwathDataFieldMetadata.iXDim >= 0 &&
1154 7 : oSwathDataFieldMetadata.iYDim >= 0)
1155 : {
1156 7 : poDS->m_nXIndex = oSwathDataFieldMetadata.iXDim;
1157 7 : poDS->m_nYIndex = oSwathDataFieldMetadata.iYDim;
1158 7 : poDS->m_nOtherDimIndex = oSwathDataFieldMetadata.iOtherDim;
1159 7 : if (!oSwathDataFieldMetadata.osLongitudeSubdataset.empty())
1160 : {
1161 14 : CPLStringList aosGeolocation;
1162 :
1163 : // Arbitrary
1164 7 : aosGeolocation.AddNameValue("SRS", SRS_WKT_WGS84_LAT_LONG);
1165 : aosGeolocation.AddNameValue(
1166 : "X_DATASET",
1167 14 : ("HDF5:\"" + osFilename +
1168 14 : "\":" + oSwathDataFieldMetadata.osLongitudeSubdataset)
1169 7 : .c_str());
1170 7 : aosGeolocation.AddNameValue("X_BAND", "1");
1171 : aosGeolocation.AddNameValue(
1172 : "Y_DATASET",
1173 14 : ("HDF5:\"" + osFilename +
1174 14 : "\":" + oSwathDataFieldMetadata.osLatitudeSubdataset)
1175 7 : .c_str());
1176 7 : aosGeolocation.AddNameValue("Y_BAND", "1");
1177 : aosGeolocation.AddNameValue(
1178 : "PIXEL_OFFSET",
1179 7 : CPLSPrintf("%d", oSwathDataFieldMetadata.nPixelOffset));
1180 : aosGeolocation.AddNameValue(
1181 : "PIXEL_STEP",
1182 7 : CPLSPrintf("%d", oSwathDataFieldMetadata.nPixelStep));
1183 : aosGeolocation.AddNameValue(
1184 : "LINE_OFFSET",
1185 7 : CPLSPrintf("%d", oSwathDataFieldMetadata.nLineOffset));
1186 : aosGeolocation.AddNameValue(
1187 : "LINE_STEP",
1188 7 : CPLSPrintf("%d", oSwathDataFieldMetadata.nLineStep));
1189 : // Not totally sure about that
1190 : aosGeolocation.AddNameValue("GEOREFERENCING_CONVENTION",
1191 7 : "PIXEL_CENTER");
1192 7 : poDS->SetMetadata(aosGeolocation.List(), "GEOLOCATION");
1193 : }
1194 : }
1195 : }
1196 : }
1197 :
1198 59 : poDS->nRasterYSize =
1199 59 : poDS->GetYIndex() < 0
1200 59 : ? 1
1201 58 : : static_cast<int>(poDS->dims[poDS->GetYIndex()]); // nRows
1202 59 : poDS->nRasterXSize =
1203 59 : static_cast<int>(poDS->dims[poDS->GetXIndex()]); // nCols
1204 59 : int nBands = 1;
1205 59 : if (poDS->m_nOtherDimIndex >= 0)
1206 : {
1207 15 : nBands = static_cast<int>(poDS->dims[poDS->m_nOtherDimIndex]);
1208 : }
1209 :
1210 118 : CPLStringList aosMetadata;
1211 59 : std::map<std::string, CPLStringList> oMapBandSpecificMetadata;
1212 59 : if (poDS->poH5Objects->nType == H5G_DATASET)
1213 : {
1214 59 : HDF5Dataset::CreateMetadata(poDS->m_hHDF5, poDS->poH5Objects,
1215 : H5G_DATASET, false, aosMetadata);
1216 59 : if (nBands > 1 && poDS->nRasterXSize != nBands &&
1217 14 : poDS->nRasterYSize != nBands)
1218 : {
1219 : // Heuristics to detect non-scalar attributes, that are intended
1220 : // to be attached to a specific band.
1221 28 : const CPLStringList aosMetadataDup(aosMetadata);
1222 52 : for (const auto &[pszKey, pszValue] :
1223 66 : cpl::IterateNameValue(aosMetadataDup))
1224 : {
1225 26 : const hid_t hAttrID = H5Aopen_name(poDS->dataset_id, pszKey);
1226 26 : const hid_t hAttrSpace = H5Aget_space(hAttrID);
1227 47 : if (H5Sget_simple_extent_ndims(hAttrSpace) == 1 &&
1228 21 : H5Sget_simple_extent_npoints(hAttrSpace) == nBands)
1229 : {
1230 : CPLStringList aosTokens(
1231 40 : CSLTokenizeString2(pszValue, " ", 0));
1232 20 : if (aosTokens.size() == nBands)
1233 : {
1234 40 : std::string osAttrName(pszKey);
1235 30 : if (osAttrName.size() > strlen("_coefficients") &&
1236 30 : osAttrName.substr(osAttrName.size() -
1237 10 : strlen("_coefficients")) ==
1238 : "_coefficients")
1239 : {
1240 5 : osAttrName.pop_back();
1241 : }
1242 25 : else if (osAttrName.size() > strlen("_wavelengths") &&
1243 25 : osAttrName.substr(osAttrName.size() -
1244 10 : strlen("_wavelengths")) ==
1245 : "_wavelengths")
1246 : {
1247 5 : osAttrName.pop_back();
1248 : }
1249 15 : else if (osAttrName.size() > strlen("_list") &&
1250 15 : osAttrName.substr(osAttrName.size() -
1251 5 : strlen("_list")) == "_list")
1252 : {
1253 5 : osAttrName.resize(osAttrName.size() -
1254 : strlen("_list"));
1255 : }
1256 20 : oMapBandSpecificMetadata[osAttrName] =
1257 40 : std::move(aosTokens);
1258 20 : aosMetadata.SetNameValue(pszKey, nullptr);
1259 : }
1260 : }
1261 26 : H5Sclose(hAttrSpace);
1262 26 : H5Aclose(hAttrID);
1263 : }
1264 : }
1265 : }
1266 :
1267 59 : poDS->m_nBlockXSize = poDS->GetRasterXSize();
1268 59 : poDS->m_nBlockYSize = 1;
1269 59 : poDS->m_nBandChunkSize = 1;
1270 :
1271 : // Check for chunksize and set it as the blocksize (optimizes read).
1272 59 : const hid_t listid = H5Dget_create_plist(poDS->dataset_id);
1273 59 : if (listid > 0)
1274 : {
1275 59 : if (H5Pget_layout(listid) == H5D_CHUNKED)
1276 : {
1277 12 : hsize_t panChunkDims[3] = {0, 0, 0};
1278 12 : const int nDimSize = H5Pget_chunk(listid, 3, panChunkDims);
1279 12 : CPL_IGNORE_RET_VAL(nDimSize);
1280 12 : CPLAssert(nDimSize == poDS->ndims);
1281 12 : poDS->m_nBlockXSize =
1282 12 : static_cast<int>(panChunkDims[poDS->GetXIndex()]);
1283 12 : if (poDS->GetYIndex() >= 0)
1284 12 : poDS->m_nBlockYSize =
1285 12 : static_cast<int>(panChunkDims[poDS->GetYIndex()]);
1286 12 : if (nBands > 1)
1287 : {
1288 3 : poDS->m_nBandChunkSize =
1289 3 : static_cast<int>(panChunkDims[poDS->m_nOtherDimIndex]);
1290 :
1291 3 : poDS->SetMetadataItem("BAND_CHUNK_SIZE",
1292 : CPLSPrintf("%d", poDS->m_nBandChunkSize),
1293 3 : "IMAGE_STRUCTURE");
1294 : }
1295 : }
1296 :
1297 59 : const int nFilters = H5Pget_nfilters(listid);
1298 65 : 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("COMPRESSION", "DEFLATE",
1308 5 : "IMAGE_STRUCTURE");
1309 : }
1310 1 : else if (eFilter == H5Z_FILTER_SZIP)
1311 : {
1312 0 : poDS->SetMetadataItem("COMPRESSION", "SZIP", "IMAGE_STRUCTURE");
1313 : }
1314 : }
1315 :
1316 59 : H5Pclose(listid);
1317 : }
1318 :
1319 185 : for (int i = 0; i < nBands; i++)
1320 : {
1321 : HDF5ImageRasterBand *const poBand =
1322 126 : new HDF5ImageRasterBand(poDS, i + 1, eGDALDataType);
1323 :
1324 126 : poDS->SetBand(i + 1, poBand);
1325 :
1326 126 : if (poDS->poH5Objects->nType == H5G_DATASET)
1327 : {
1328 126 : poBand->SetMetadata(aosMetadata.List());
1329 166 : for (const auto &oIter : oMapBandSpecificMetadata)
1330 : {
1331 40 : poBand->SetMetadataItem(oIter.first.c_str(), oIter.second[i]);
1332 : }
1333 : }
1334 : }
1335 :
1336 59 : if (!poDS->GetMetadata("GEOLOCATION"))
1337 52 : poDS->CreateProjections();
1338 :
1339 : // Setup/check for pam .aux.xml.
1340 59 : poDS->TryLoadXML();
1341 :
1342 : // If the PAM .aux.xml file contains the serialized GEOLOCATION metadata
1343 : // domain, make sure to patch the X_DATASET and Y_DATASET values, when
1344 : // they point to the current file, to use the filename with which we have
1345 : // opened the dataset. Helps in scenarios where the .aux.xml file has
1346 : // been saved with a relative filename for example.
1347 : // Or scenarios like https://github.com/OSGeo/gdal/issues/12824 mixing
1348 : // native Windows and WSL use.
1349 : CSLConstList papszGeoLocationAfterLoadXML =
1350 59 : poDS->GetMetadata("GEOLOCATION");
1351 59 : if (papszGeoLocationAfterLoadXML)
1352 : {
1353 21 : for (const char *pszItem : {"X_DATASET", "Y_DATASET"})
1354 : {
1355 : const char *pszSubdataset =
1356 14 : CSLFetchNameValue(papszGeoLocationAfterLoadXML, pszItem);
1357 28 : if (pszSubdataset && STARTS_WITH(pszSubdataset, "HDF5:\"") &&
1358 14 : strstr(pszSubdataset, CPLGetFilename(osFilename.c_str())))
1359 : {
1360 14 : auto hSubDSInfo = GDALGetSubdatasetInfo(pszSubdataset);
1361 14 : if (hSubDSInfo)
1362 : {
1363 : char *pszOriPath =
1364 14 : GDALSubdatasetInfoGetPathComponent(hSubDSInfo);
1365 14 : if (EQUAL(CPLGetFilename(pszOriPath),
1366 : CPLGetFilename(osFilename.c_str())))
1367 : {
1368 14 : char *pszNewVal = GDALSubdatasetInfoModifyPathComponent(
1369 : hSubDSInfo, osFilename.c_str());
1370 14 : poDS->SetMetadataItem(pszItem, pszNewVal,
1371 14 : "GEOLOCATION");
1372 14 : CPLFree(pszNewVal);
1373 : }
1374 14 : CPLFree(pszOriPath);
1375 14 : GDALDestroySubdatasetInfo(hSubDSInfo);
1376 : }
1377 : }
1378 : }
1379 : }
1380 :
1381 59 : poDS->SetPamFlags(poDS->GetPamFlags() & ~GPF_DIRTY);
1382 :
1383 : // Setup overviews.
1384 59 : poDS->oOvManager.Initialize(poDS, ":::VIRTUAL:::");
1385 :
1386 59 : return poDS;
1387 : }
1388 :
1389 : /************************************************************************/
1390 : /* HDF5ImageDatasetDriverUnload() */
1391 : /************************************************************************/
1392 :
1393 323 : static void HDF5ImageDatasetDriverUnload(GDALDriver *)
1394 : {
1395 323 : HDF5UnloadFileDriver();
1396 323 : }
1397 :
1398 : /************************************************************************/
1399 : /* GDALRegister_HDF5Image() */
1400 : /************************************************************************/
1401 355 : void GDALRegister_HDF5Image()
1402 :
1403 : {
1404 355 : if (!GDAL_CHECK_VERSION("HDF5Image driver"))
1405 0 : return;
1406 :
1407 355 : if (GDALGetDriverByName(HDF5_IMAGE_DRIVER_NAME) != nullptr)
1408 0 : return;
1409 :
1410 355 : GDALDriver *poDriver = new GDALDriver();
1411 :
1412 355 : HDF5ImageDriverSetCommonMetadata(poDriver);
1413 :
1414 355 : poDriver->pfnOpen = HDF5ImageDataset::Open;
1415 355 : poDriver->pfnUnloadDriver = HDF5ImageDatasetDriverUnload;
1416 :
1417 355 : GetGDALDriverManager()->RegisterDriver(poDriver);
1418 : }
1419 :
1420 : /************************************************************************/
1421 : /* CreateODIMH5Projection() */
1422 : /************************************************************************/
1423 :
1424 : // Reference:
1425 : // http://www.knmi.nl/opera/opera3/OPERA_2008_03_WP2.1b_ODIM_H5_v2.1.pdf
1426 : //
1427 : // 4.3.2 where for geographically referenced image Groups
1428 : // We don't use the where_xscale and where_yscale parameters, but recompute them
1429 : // from the lower-left and upper-right coordinates. There's some difference.
1430 : // As all those parameters are linked together, I'm not sure which one should be
1431 : // considered as the reference.
1432 :
1433 0 : CPLErr HDF5ImageDataset::CreateODIMH5Projection()
1434 : {
1435 0 : const char *const pszProj4String = GetMetadataItem("where_projdef");
1436 0 : const char *const pszLL_lon = GetMetadataItem("where_LL_lon");
1437 0 : const char *const pszLL_lat = GetMetadataItem("where_LL_lat");
1438 0 : const char *const pszUR_lon = GetMetadataItem("where_UR_lon");
1439 0 : const char *const pszUR_lat = GetMetadataItem("where_UR_lat");
1440 0 : if (pszProj4String == nullptr || pszLL_lon == nullptr ||
1441 0 : pszLL_lat == nullptr || pszUR_lon == nullptr || pszUR_lat == nullptr)
1442 0 : return CE_Failure;
1443 :
1444 0 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1445 0 : if (m_oSRS.importFromProj4(pszProj4String) != OGRERR_NONE)
1446 0 : return CE_Failure;
1447 :
1448 0 : OGRSpatialReference oSRSWGS84;
1449 0 : oSRSWGS84.SetWellKnownGeogCS("WGS84");
1450 0 : oSRSWGS84.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1451 :
1452 : OGRCoordinateTransformation *poCT =
1453 0 : OGRCreateCoordinateTransformation(&oSRSWGS84, &m_oSRS);
1454 0 : if (poCT == nullptr)
1455 0 : return CE_Failure;
1456 :
1457 : // Reproject corners from long,lat WGS84 to the target SRS.
1458 0 : double dfLLX = CPLAtof(pszLL_lon);
1459 0 : double dfLLY = CPLAtof(pszLL_lat);
1460 0 : double dfURX = CPLAtof(pszUR_lon);
1461 0 : double dfURY = CPLAtof(pszUR_lat);
1462 0 : if (!poCT->Transform(1, &dfLLX, &dfLLY) ||
1463 0 : !poCT->Transform(1, &dfURX, &dfURY))
1464 : {
1465 0 : delete poCT;
1466 0 : return CE_Failure;
1467 : }
1468 0 : delete poCT;
1469 :
1470 : // Compute the geotransform now.
1471 0 : const double dfPixelX = (dfURX - dfLLX) / nRasterXSize;
1472 0 : const double dfPixelY = (dfURY - dfLLY) / nRasterYSize;
1473 :
1474 0 : bHasGeoTransform = true;
1475 0 : m_gt[0] = dfLLX;
1476 0 : m_gt[1] = dfPixelX;
1477 0 : m_gt[2] = 0;
1478 0 : m_gt[3] = dfURY;
1479 0 : m_gt[4] = 0;
1480 0 : m_gt[5] = -dfPixelY;
1481 :
1482 0 : return CE_None;
1483 : }
1484 :
1485 : /************************************************************************/
1486 : /* CreateProjections() */
1487 : /************************************************************************/
1488 52 : CPLErr HDF5ImageDataset::CreateProjections()
1489 : {
1490 52 : switch (iSubdatasetType)
1491 : {
1492 2 : case CSK_PRODUCT:
1493 : {
1494 2 : int productType = PROD_UNKNOWN;
1495 :
1496 2 : if (GetMetadataItem("Product_Type") != nullptr)
1497 : {
1498 : // Get the format's level.
1499 : const char *osMissionLevel =
1500 2 : HDF5Dataset::GetMetadataItem("Product_Type");
1501 :
1502 2 : if (STARTS_WITH_CI(osMissionLevel, "RAW"))
1503 0 : productType = PROD_CSK_L0;
1504 :
1505 2 : if (STARTS_WITH_CI(osMissionLevel, "SSC"))
1506 0 : productType = PROD_CSK_L1A;
1507 :
1508 2 : if (STARTS_WITH_CI(osMissionLevel, "DGM"))
1509 1 : productType = PROD_CSK_L1B;
1510 :
1511 2 : if (STARTS_WITH_CI(osMissionLevel, "GEC"))
1512 1 : productType = PROD_CSK_L1C;
1513 :
1514 2 : if (STARTS_WITH_CI(osMissionLevel, "GTC"))
1515 0 : productType = PROD_CSK_L1D;
1516 : }
1517 :
1518 2 : CaptureCSKGeoTransform(productType);
1519 2 : CaptureCSKGeolocation(productType);
1520 2 : CaptureCSKGCPs(productType);
1521 :
1522 2 : break;
1523 : }
1524 50 : case UNKNOWN_PRODUCT:
1525 : {
1526 50 : constexpr int NBGCPLAT = 100;
1527 50 : constexpr int NBGCPLON = 30;
1528 :
1529 50 : const int nDeltaLat = nRasterYSize / NBGCPLAT;
1530 50 : const int nDeltaLon = nRasterXSize / NBGCPLON;
1531 :
1532 50 : if (nDeltaLat == 0 || nDeltaLon == 0)
1533 42 : return CE_None;
1534 :
1535 : // Create HDF5 Data Hierarchy in a link list.
1536 8 : poH5Objects = HDF5FindDatasetObjects(poH5RootGroup, "Latitude");
1537 8 : if (!poH5Objects)
1538 : {
1539 8 : if (GetMetadataItem("where_projdef") != nullptr)
1540 0 : return CreateODIMH5Projection();
1541 8 : return CE_None;
1542 : }
1543 :
1544 : // The Latitude and Longitude arrays must have a rank of 2 to
1545 : // retrieve GCPs.
1546 0 : if (poH5Objects->nRank != 2 ||
1547 0 : poH5Objects->paDims[0] != static_cast<size_t>(nRasterYSize) ||
1548 0 : poH5Objects->paDims[1] != static_cast<size_t>(nRasterXSize))
1549 : {
1550 0 : return CE_None;
1551 : }
1552 :
1553 : // Retrieve HDF5 data information.
1554 : const hid_t LatitudeDatasetID =
1555 0 : H5Dopen(m_hHDF5, poH5Objects->pszPath);
1556 : // LatitudeDataspaceID = H5Dget_space(dataset_id);
1557 :
1558 0 : poH5Objects = HDF5FindDatasetObjects(poH5RootGroup, "Longitude");
1559 : // GCPs.
1560 0 : if (poH5Objects == nullptr || poH5Objects->nRank != 2 ||
1561 0 : poH5Objects->paDims[0] != static_cast<size_t>(nRasterYSize) ||
1562 0 : poH5Objects->paDims[1] != static_cast<size_t>(nRasterXSize))
1563 : {
1564 0 : if (LatitudeDatasetID > 0)
1565 0 : H5Dclose(LatitudeDatasetID);
1566 0 : return CE_None;
1567 : }
1568 :
1569 : const hid_t LongitudeDatasetID =
1570 0 : H5Dopen(m_hHDF5, poH5Objects->pszPath);
1571 : // LongitudeDataspaceID = H5Dget_space(dataset_id);
1572 :
1573 0 : if (LatitudeDatasetID > 0 && LongitudeDatasetID > 0)
1574 : {
1575 : float *const Latitude =
1576 0 : static_cast<float *>(VSI_MALLOC3_VERBOSE(
1577 : nRasterYSize, nRasterXSize, sizeof(float)));
1578 : float *const Longitude =
1579 0 : static_cast<float *>(VSI_MALLOC3_VERBOSE(
1580 : nRasterYSize, nRasterXSize, sizeof(float)));
1581 0 : if (!Latitude || !Longitude)
1582 : {
1583 0 : CPLFree(Latitude);
1584 0 : CPLFree(Longitude);
1585 0 : H5Dclose(LatitudeDatasetID);
1586 0 : H5Dclose(LongitudeDatasetID);
1587 0 : return CE_Failure;
1588 : }
1589 0 : memset(Latitude, 0,
1590 0 : static_cast<size_t>(nRasterXSize) * nRasterYSize *
1591 : sizeof(float));
1592 0 : memset(Longitude, 0,
1593 0 : static_cast<size_t>(nRasterXSize) * nRasterYSize *
1594 : sizeof(float));
1595 :
1596 : // netCDF convention for nodata
1597 0 : double dfLatNoData = 0;
1598 0 : bool bHasLatNoData = GH5_FetchAttribute(
1599 : LatitudeDatasetID, "_FillValue", dfLatNoData);
1600 :
1601 0 : double dfLongNoData = 0;
1602 0 : bool bHasLongNoData = GH5_FetchAttribute(
1603 : LongitudeDatasetID, "_FillValue", dfLongNoData);
1604 :
1605 0 : H5Dread(LatitudeDatasetID, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
1606 : H5P_DEFAULT, Latitude);
1607 :
1608 0 : H5Dread(LongitudeDatasetID, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
1609 : H5P_DEFAULT, Longitude);
1610 :
1611 0 : m_oSRS.Clear();
1612 0 : m_oGCPSRS.SetWellKnownGeogCS("WGS84");
1613 :
1614 0 : const int nYLimit =
1615 0 : (static_cast<int>(nRasterYSize) / nDeltaLat) * nDeltaLat;
1616 0 : const int nXLimit =
1617 0 : (static_cast<int>(nRasterXSize) / nDeltaLon) * nDeltaLon;
1618 :
1619 : // The original code in
1620 : // https://trac.osgeo.org/gdal/changeset/8066 always add +180 to
1621 : // the longitudes, but without justification I suspect this
1622 : // might be due to handling products crossing the antimeridian.
1623 : // Trying to do it just when needed through a heuristics.
1624 0 : bool bHasLonNearMinus180 = false;
1625 0 : bool bHasLonNearPlus180 = false;
1626 0 : bool bHasLonNearZero = false;
1627 0 : for (int j = 0; j < nYLimit; j += nDeltaLat)
1628 : {
1629 0 : for (int i = 0; i < nXLimit; i += nDeltaLon)
1630 : {
1631 0 : const int iGCP = j * nRasterXSize + i;
1632 0 : if ((bHasLatNoData && static_cast<float>(dfLatNoData) ==
1633 0 : Latitude[iGCP]) ||
1634 0 : (bHasLongNoData &&
1635 0 : static_cast<float>(dfLongNoData) ==
1636 0 : Longitude[iGCP]))
1637 0 : continue;
1638 0 : if (Longitude[iGCP] > 170 && Longitude[iGCP] <= 180)
1639 0 : bHasLonNearPlus180 = true;
1640 0 : if (Longitude[iGCP] < -170 && Longitude[iGCP] >= -180)
1641 0 : bHasLonNearMinus180 = true;
1642 0 : if (fabs(Longitude[iGCP]) < 90)
1643 0 : bHasLonNearZero = true;
1644 : }
1645 : }
1646 :
1647 : // Fill the GCPs list.
1648 : const char *pszShiftGCP =
1649 0 : CPLGetConfigOption("HDF5_SHIFT_GCPX_BY_180", nullptr);
1650 : const bool bAdd180 =
1651 0 : (bHasLonNearPlus180 && bHasLonNearMinus180 &&
1652 0 : !bHasLonNearZero && pszShiftGCP == nullptr) ||
1653 0 : (pszShiftGCP != nullptr && CPLTestBool(pszShiftGCP));
1654 :
1655 0 : for (int j = 0; j < nYLimit; j += nDeltaLat)
1656 : {
1657 0 : for (int i = 0; i < nXLimit; i += nDeltaLon)
1658 : {
1659 0 : const int iGCP = j * nRasterXSize + i;
1660 0 : if ((bHasLatNoData && static_cast<float>(dfLatNoData) ==
1661 0 : Latitude[iGCP]) ||
1662 0 : (bHasLongNoData &&
1663 0 : static_cast<float>(dfLongNoData) ==
1664 0 : Longitude[iGCP]))
1665 0 : continue;
1666 0 : double dfGCPX = static_cast<double>(Longitude[iGCP]);
1667 0 : if (bAdd180)
1668 0 : dfGCPX += 180.0;
1669 0 : const double dfGCPY =
1670 0 : static_cast<double>(Latitude[iGCP]);
1671 :
1672 0 : m_aoGCPs.emplace_back("", "", i + 0.5, j + 0.5, dfGCPX,
1673 0 : dfGCPY);
1674 : }
1675 : }
1676 :
1677 0 : CPLFree(Latitude);
1678 0 : CPLFree(Longitude);
1679 : }
1680 :
1681 0 : if (LatitudeDatasetID > 0)
1682 0 : H5Dclose(LatitudeDatasetID);
1683 0 : if (LongitudeDatasetID > 0)
1684 0 : H5Dclose(LongitudeDatasetID);
1685 :
1686 0 : break;
1687 : }
1688 : }
1689 :
1690 2 : return CE_None;
1691 : }
1692 :
1693 : /************************************************************************/
1694 : /* GetMetadataItem() */
1695 : /************************************************************************/
1696 :
1697 123 : const char *HDF5ImageDataset::GetMetadataItem(const char *pszName,
1698 : const char *pszDomain)
1699 : {
1700 123 : if (pszDomain && EQUAL(pszDomain, "__DEBUG__") &&
1701 66 : EQUAL(pszName, "WholeBandChunkOptim"))
1702 : {
1703 66 : switch (m_eWholeBandChunkOptim)
1704 : {
1705 4 : case WBC_DETECTION_IN_PROGRESS:
1706 4 : return "DETECTION_IN_PROGRESS";
1707 3 : case WBC_DISABLED:
1708 3 : return "DISABLED";
1709 59 : case WBC_ENABLED:
1710 59 : return "ENABLED";
1711 : }
1712 : }
1713 57 : return GDALPamDataset::GetMetadataItem(pszName, pszDomain);
1714 : }
1715 :
1716 : /************************************************************************/
1717 : /* GetSpatialRef() */
1718 : /************************************************************************/
1719 :
1720 12 : const OGRSpatialReference *HDF5ImageDataset::GetSpatialRef() const
1721 : {
1722 12 : if (!m_oSRS.IsEmpty())
1723 11 : return &m_oSRS;
1724 1 : return GDALPamDataset::GetSpatialRef();
1725 : }
1726 :
1727 : /************************************************************************/
1728 : /* GetGCPCount() */
1729 : /************************************************************************/
1730 :
1731 3 : int HDF5ImageDataset::GetGCPCount()
1732 :
1733 : {
1734 3 : if (!m_aoGCPs.empty())
1735 1 : return static_cast<int>(m_aoGCPs.size());
1736 :
1737 2 : return GDALPamDataset::GetGCPCount();
1738 : }
1739 :
1740 : /************************************************************************/
1741 : /* GetGCPSpatialRef() */
1742 : /************************************************************************/
1743 :
1744 1 : const OGRSpatialReference *HDF5ImageDataset::GetGCPSpatialRef() const
1745 :
1746 : {
1747 1 : if (!m_aoGCPs.empty() && !m_oGCPSRS.IsEmpty())
1748 1 : return &m_oGCPSRS;
1749 :
1750 0 : return GDALPamDataset::GetGCPSpatialRef();
1751 : }
1752 :
1753 : /************************************************************************/
1754 : /* GetGCPs() */
1755 : /************************************************************************/
1756 :
1757 2 : const GDAL_GCP *HDF5ImageDataset::GetGCPs()
1758 : {
1759 2 : if (!m_aoGCPs.empty())
1760 1 : return gdal::GCP::c_ptr(m_aoGCPs);
1761 :
1762 1 : return GDALPamDataset::GetGCPs();
1763 : }
1764 :
1765 : /************************************************************************/
1766 : /* GetGeoTransform() */
1767 : /************************************************************************/
1768 :
1769 8 : CPLErr HDF5ImageDataset::GetGeoTransform(GDALGeoTransform >) const
1770 : {
1771 8 : if (bHasGeoTransform)
1772 : {
1773 7 : gt = m_gt;
1774 7 : return CE_None;
1775 : }
1776 :
1777 1 : return GDALPamDataset::GetGeoTransform(gt);
1778 : }
1779 :
1780 : /************************************************************************/
1781 : /* IdentifyProductType() */
1782 : /************************************************************************/
1783 :
1784 : /**
1785 : * Identify if the subdataset has a known product format
1786 : * It stores a product identifier in iSubdatasetType,
1787 : * UNKNOWN_PRODUCT, if it isn't a recognizable format.
1788 : */
1789 59 : void HDF5ImageDataset::IdentifyProductType()
1790 : {
1791 59 : iSubdatasetType = UNKNOWN_PRODUCT;
1792 :
1793 : // COSMO-SKYMED
1794 :
1795 : // Get the Mission Id as a char *, because the field may not exist.
1796 59 : const char *const pszMissionId = HDF5Dataset::GetMetadataItem("Mission_ID");
1797 :
1798 : // If there is a Mission_ID field.
1799 59 : if (pszMissionId != nullptr && strstr(GetDescription(), "QLK") == nullptr)
1800 : {
1801 : // Check if the mission type is CSK, KMPS or CSG.
1802 : // KMPS: Komsat-5 is Korean mission with a SAR instrument.
1803 : // CSG: Cosmo Skymed 2nd Generation
1804 2 : if (EQUAL(pszMissionId, "CSK") || EQUAL(pszMissionId, "KMPS") ||
1805 0 : EQUAL(pszMissionId, "CSG"))
1806 : {
1807 2 : iSubdatasetType = CSK_PRODUCT;
1808 :
1809 2 : if (GetMetadataItem("Product_Type") != nullptr)
1810 : {
1811 : // Get the format's level.
1812 : const char *osMissionLevel =
1813 2 : HDF5Dataset::GetMetadataItem("Product_Type");
1814 :
1815 2 : if (STARTS_WITH_CI(osMissionLevel, "RAW"))
1816 0 : iCSKProductType = PROD_CSK_L0;
1817 :
1818 2 : if (STARTS_WITH_CI(osMissionLevel, "SCS"))
1819 0 : iCSKProductType = PROD_CSK_L1A;
1820 :
1821 2 : if (STARTS_WITH_CI(osMissionLevel, "DGM"))
1822 1 : iCSKProductType = PROD_CSK_L1B;
1823 :
1824 2 : if (STARTS_WITH_CI(osMissionLevel, "GEC"))
1825 1 : iCSKProductType = PROD_CSK_L1C;
1826 :
1827 2 : if (STARTS_WITH_CI(osMissionLevel, "GTC"))
1828 0 : iCSKProductType = PROD_CSK_L1D;
1829 : }
1830 : }
1831 : }
1832 59 : }
1833 :
1834 : /************************************************************************/
1835 : /* CaptureCSKGeolocation() */
1836 : /************************************************************************/
1837 :
1838 : /**
1839 : * Captures Geolocation information from a COSMO-SKYMED
1840 : * file.
1841 : * The geoid will always be WGS84
1842 : * The projection type may be UTM or UPS, depending on the
1843 : * latitude from the center of the image.
1844 : * @param iProductType type of CSK subproduct, see HDF5CSKProduct
1845 : */
1846 2 : void HDF5ImageDataset::CaptureCSKGeolocation(int iProductType)
1847 : {
1848 : // Set the ellipsoid to WGS84.
1849 2 : m_oSRS.SetWellKnownGeogCS("WGS84");
1850 :
1851 2 : if (iProductType == PROD_CSK_L1C || iProductType == PROD_CSK_L1D)
1852 : {
1853 1 : double *dfProjFalseEastNorth = nullptr;
1854 1 : double *dfProjScaleFactor = nullptr;
1855 1 : double *dfCenterCoord = nullptr;
1856 :
1857 : // Check if all the metadata attributes are present.
1858 1 : if (HDF5ReadDoubleAttr("Map Projection False East-North",
1859 1 : &dfProjFalseEastNorth) == CE_Failure ||
1860 1 : HDF5ReadDoubleAttr("Map Projection Scale Factor",
1861 1 : &dfProjScaleFactor) == CE_Failure ||
1862 1 : HDF5ReadDoubleAttr("Map Projection Centre", &dfCenterCoord) ==
1863 2 : CE_Failure ||
1864 1 : GetMetadataItem("Projection_ID") == nullptr)
1865 : {
1866 0 : m_oSRS.Clear();
1867 0 : m_oGCPSRS.Clear();
1868 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1869 : "The CSK hdf5 file geolocation information is "
1870 : "malformed");
1871 : }
1872 : else
1873 : {
1874 : // Fetch projection Type.
1875 2 : CPLString osProjectionID = GetMetadataItem("Projection_ID");
1876 :
1877 : // If the projection is UTM.
1878 1 : if (EQUAL(osProjectionID, "UTM"))
1879 : {
1880 : // @TODO: use SetUTM
1881 1 : m_oSRS.SetProjCS(SRS_PT_TRANSVERSE_MERCATOR);
1882 1 : m_oSRS.SetTM(dfCenterCoord[0], dfCenterCoord[1],
1883 : dfProjScaleFactor[0], dfProjFalseEastNorth[0],
1884 1 : dfProjFalseEastNorth[1]);
1885 : }
1886 : else
1887 : {
1888 : // TODO: No UPS projected files to test!
1889 : // If the projection is UPS.
1890 0 : if (EQUAL(osProjectionID, "UPS"))
1891 : {
1892 0 : m_oSRS.SetProjCS(SRS_PT_POLAR_STEREOGRAPHIC);
1893 0 : m_oSRS.SetPS(dfCenterCoord[0], dfCenterCoord[1],
1894 : dfProjScaleFactor[0], dfProjFalseEastNorth[0],
1895 0 : dfProjFalseEastNorth[1]);
1896 : }
1897 : }
1898 :
1899 1 : CPLFree(dfCenterCoord);
1900 1 : CPLFree(dfProjScaleFactor);
1901 1 : CPLFree(dfProjFalseEastNorth);
1902 1 : }
1903 : }
1904 : else
1905 : {
1906 1 : m_oGCPSRS = m_oSRS;
1907 : }
1908 2 : }
1909 :
1910 : /************************************************************************/
1911 : /* CaptureCSKGeoTransform() */
1912 : /************************************************************************/
1913 :
1914 : /**
1915 : * Get Geotransform information for COSMO-SKYMED files
1916 : * In case of success it stores the transformation
1917 : * in m_gt. In case of failure it doesn't
1918 : * modify m_gt
1919 : * @param iProductType type of CSK subproduct, see HDF5CSKProduct
1920 : */
1921 2 : void HDF5ImageDataset::CaptureCSKGeoTransform(int iProductType)
1922 : {
1923 2 : const char *pszSubdatasetName = GetSubdatasetName();
1924 :
1925 2 : bHasGeoTransform = false;
1926 : // If the product level is not L1C or L1D then
1927 : // it doesn't have a valid projection.
1928 2 : if (iProductType == PROD_CSK_L1C || iProductType == PROD_CSK_L1D)
1929 : {
1930 : // If there is a subdataset.
1931 1 : if (pszSubdatasetName != nullptr)
1932 : {
1933 2 : CPLString osULPath = pszSubdatasetName;
1934 1 : osULPath += "/Top Left East-North";
1935 :
1936 2 : CPLString osLineSpacingPath = pszSubdatasetName;
1937 1 : osLineSpacingPath += "/Line Spacing";
1938 :
1939 2 : CPLString osColumnSpacingPath = pszSubdatasetName;
1940 1 : osColumnSpacingPath += "/Column Spacing";
1941 :
1942 1 : double *pdOutUL = nullptr;
1943 1 : double *pdLineSpacing = nullptr;
1944 1 : double *pdColumnSpacing = nullptr;
1945 :
1946 : // If it could find the attributes on the metadata.
1947 1 : if (HDF5ReadDoubleAttr(osULPath.c_str(), &pdOutUL) == CE_Failure ||
1948 1 : HDF5ReadDoubleAttr(osLineSpacingPath.c_str(), &pdLineSpacing) ==
1949 2 : CE_Failure ||
1950 1 : HDF5ReadDoubleAttr(osColumnSpacingPath.c_str(),
1951 : &pdColumnSpacing) == CE_Failure)
1952 : {
1953 0 : bHasGeoTransform = false;
1954 : }
1955 : else
1956 : {
1957 : // geotransform[1] : width of pixel
1958 : // geotransform[4] : rotational coefficient, zero for north up
1959 : // images.
1960 : // geotransform[2] : rotational coefficient, zero for north up
1961 : // images.
1962 : // geotransform[5] : height of pixel (but negative)
1963 :
1964 1 : m_gt[0] = pdOutUL[0];
1965 1 : m_gt[1] = pdLineSpacing[0];
1966 1 : m_gt[2] = 0;
1967 1 : m_gt[3] = pdOutUL[1];
1968 1 : m_gt[4] = 0;
1969 1 : m_gt[5] = -pdColumnSpacing[0];
1970 :
1971 1 : CPLFree(pdOutUL);
1972 1 : CPLFree(pdLineSpacing);
1973 1 : CPLFree(pdColumnSpacing);
1974 :
1975 1 : bHasGeoTransform = true;
1976 : }
1977 : }
1978 : }
1979 2 : }
1980 :
1981 : /************************************************************************/
1982 : /* CaptureCSKGCPs() */
1983 : /************************************************************************/
1984 :
1985 : /**
1986 : * Retrieves and stores the GCPs from a COSMO-SKYMED dataset
1987 : * It only retrieves the GCPs for L0, L1A and L1B products
1988 : * for L1C and L1D products, geotransform is provided.
1989 : * The GCPs provided will be the Image's corners.
1990 : * @param iProductType type of CSK product @see HDF5CSKProductEnum
1991 : */
1992 2 : void HDF5ImageDataset::CaptureCSKGCPs(int iProductType)
1993 : {
1994 : // Only retrieve GCPs for L0,L1A and L1B products.
1995 2 : if (iProductType == PROD_CSK_L0 || iProductType == PROD_CSK_L1A ||
1996 : iProductType == PROD_CSK_L1B)
1997 : {
1998 10 : CPLString osCornerName[4];
1999 1 : double pdCornerPixel[4] = {0.0, 0.0, 0.0, 0.0};
2000 1 : double pdCornerLine[4] = {0.0, 0.0, 0.0, 0.0};
2001 :
2002 1 : const char *const pszSubdatasetName = GetSubdatasetName();
2003 :
2004 : // Load the subdataset name first.
2005 5 : for (int i = 0; i < 4; i++)
2006 4 : osCornerName[i] = pszSubdatasetName;
2007 :
2008 : // Load the attribute name, and raster coordinates for
2009 : // all the corners.
2010 1 : osCornerName[0] += "/Top Left Geodetic Coordinates";
2011 1 : pdCornerPixel[0] = 0;
2012 1 : pdCornerLine[0] = 0;
2013 :
2014 1 : osCornerName[1] += "/Top Right Geodetic Coordinates";
2015 1 : pdCornerPixel[1] = GetRasterXSize();
2016 1 : pdCornerLine[1] = 0;
2017 :
2018 1 : osCornerName[2] += "/Bottom Left Geodetic Coordinates";
2019 1 : pdCornerPixel[2] = 0;
2020 1 : pdCornerLine[2] = GetRasterYSize();
2021 :
2022 1 : osCornerName[3] += "/Bottom Right Geodetic Coordinates";
2023 1 : pdCornerPixel[3] = GetRasterXSize();
2024 1 : pdCornerLine[3] = GetRasterYSize();
2025 :
2026 : // For all the image's corners.
2027 5 : for (int i = 0; i < 4; i++)
2028 : {
2029 4 : double *pdCornerCoordinates = nullptr;
2030 :
2031 : // Retrieve the attributes.
2032 4 : if (HDF5ReadDoubleAttr(osCornerName[i].c_str(),
2033 4 : &pdCornerCoordinates) == CE_Failure)
2034 : {
2035 0 : CPLError(CE_Failure, CPLE_OpenFailed,
2036 : "Error retrieving CSK GCPs");
2037 0 : m_aoGCPs.clear();
2038 0 : break;
2039 : }
2040 :
2041 0 : m_aoGCPs.emplace_back(osCornerName[i].c_str(), "", pdCornerPixel[i],
2042 4 : pdCornerLine[i],
2043 4 : /* X = */ pdCornerCoordinates[1],
2044 : /* Y = */ pdCornerCoordinates[0],
2045 4 : /* Z = */ pdCornerCoordinates[2]);
2046 :
2047 : // Free the returned coordinates.
2048 4 : CPLFree(pdCornerCoordinates);
2049 : }
2050 : }
2051 2 : }
|
