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 180 : Hdf5ProductType GetSubdatasetType() const
134 : {
135 180 : return iSubdatasetType;
136 : }
137 :
138 6 : HDF5CSKProductEnum GetCSKProductType() const
139 : {
140 6 : return iCSKProductType;
141 : }
142 :
143 180 : int IsComplexCSKL1A() const
144 : {
145 186 : return GetSubdatasetType() == CSK_PRODUCT &&
146 186 : GetCSKProductType() == PROD_CSK_L1A && ndims == 3;
147 : }
148 :
149 568 : int GetYIndex() const
150 : {
151 568 : return m_nYIndex;
152 : }
153 :
154 701 : int GetXIndex() const
155 : {
156 701 : 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 61 : 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 61 : bHasGeoTransform(false)
205 : {
206 61 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
207 61 : m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
208 61 : }
209 :
210 : /************************************************************************/
211 : /* ~HDF5ImageDataset() */
212 : /************************************************************************/
213 122 : HDF5ImageDataset::~HDF5ImageDataset()
214 : {
215 : HDF5_GLOBAL_LOCK();
216 :
217 61 : FlushCache(true);
218 :
219 61 : if (dataset_id > 0)
220 61 : H5Dclose(dataset_id);
221 61 : if (dataspace_id > 0)
222 61 : H5Sclose(dataspace_id);
223 61 : if (native > 0)
224 60 : H5Tclose(native);
225 :
226 61 : CPLFree(dims);
227 61 : CPLFree(maxdims);
228 122 : }
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 254 : HDF5ImageRasterBand::~HDF5ImageRasterBand()
269 : {
270 254 : }
271 :
272 : /************************************************************************/
273 : /* HDF5ImageRasterBand() */
274 : /************************************************************************/
275 127 : HDF5ImageRasterBand::HDF5ImageRasterBand(HDF5ImageDataset *poDSIn, int nBandIn,
276 127 : GDALDataType eType)
277 : {
278 127 : poDS = poDSIn;
279 127 : nBand = nBandIn;
280 127 : eDataType = eType;
281 127 : nBlockXSize = poDSIn->m_nBlockXSize;
282 127 : nBlockYSize = poDSIn->m_nBlockYSize;
283 :
284 : // netCDF convention for nodata
285 127 : bNoDataSet =
286 127 : GH5_FetchAttribute(poDSIn->dataset_id, "_FillValue", dfNoDataValue);
287 127 : if (!bNoDataSet)
288 118 : dfNoDataValue = -9999.0;
289 :
290 : // netCDF conventions for scale and offset
291 127 : m_bHasOffset =
292 127 : GH5_FetchAttribute(poDSIn->dataset_id, "add_offset", m_dfOffset);
293 127 : if (!m_bHasOffset)
294 126 : m_dfOffset = 0.0;
295 127 : m_bHasScale =
296 127 : GH5_FetchAttribute(poDSIn->dataset_id, "scale_factor", m_dfScale);
297 127 : if (!m_bHasScale)
298 126 : m_dfScale = 1.0;
299 127 : }
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 61 : GDALDataset *HDF5ImageDataset::Open(GDALOpenInfo *poOpenInfo)
963 : {
964 61 : 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 61 : if (poOpenInfo->eAccess == GA_Update)
971 : {
972 0 : ReportUpdateNotSupportedByDriver("HDF5");
973 0 : return nullptr;
974 : }
975 :
976 61 : HDF5ImageDataset *poDS = new HDF5ImageDataset();
977 :
978 : // Create a corresponding GDALDataset.
979 : char **papszName =
980 61 : CSLTokenizeString2(poOpenInfo->pszFilename, ":",
981 : CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES);
982 :
983 61 : if (!(CSLCount(papszName) == 3 || CSLCount(papszName) == 4))
984 : {
985 0 : CSLDestroy(papszName);
986 0 : delete poDS;
987 0 : return nullptr;
988 : }
989 :
990 61 : poDS->SetDescription(poOpenInfo->pszFilename);
991 :
992 : // Check for drive name in windows HDF5:"D:\...
993 122 : CPLString osSubdatasetName;
994 :
995 122 : CPLString osFilename(papszName[1]);
996 :
997 61 : if ((strlen(papszName[1]) == 1 && papszName[3] != nullptr) ||
998 61 : (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 61 : osSubdatasetName = papszName[2];
1007 : }
1008 :
1009 61 : poDS->SetSubdatasetName(osSubdatasetName);
1010 :
1011 61 : CSLDestroy(papszName);
1012 61 : papszName = nullptr;
1013 :
1014 61 : poDS->SetPhysicalFilename(osFilename);
1015 :
1016 : // Try opening the dataset.
1017 61 : poDS->m_hHDF5 = GDAL_HDF5Open(osFilename);
1018 61 : if (poDS->m_hHDF5 < 0)
1019 : {
1020 0 : delete poDS;
1021 0 : return nullptr;
1022 : }
1023 :
1024 61 : poDS->hGroupID = H5Gopen(poDS->m_hHDF5, "/");
1025 61 : if (poDS->hGroupID < 0)
1026 : {
1027 0 : delete poDS;
1028 0 : return nullptr;
1029 : }
1030 :
1031 : // This is an HDF5 file.
1032 61 : poDS->ReadGlobalAttributes(FALSE);
1033 :
1034 : // Create HDF5 Data Hierarchy in a link list.
1035 61 : poDS->poH5Objects = poDS->HDF5FindDatasetObjectsbyPath(poDS->poH5RootGroup,
1036 : osSubdatasetName);
1037 :
1038 61 : if (poDS->poH5Objects == nullptr)
1039 : {
1040 0 : delete poDS;
1041 0 : return nullptr;
1042 : }
1043 :
1044 : // Retrieve HDF5 data information.
1045 61 : poDS->dataset_id = H5Dopen(poDS->m_hHDF5, poDS->poH5Objects->pszPath);
1046 61 : poDS->dataspace_id = H5Dget_space(poDS->dataset_id);
1047 61 : poDS->ndims = H5Sget_simple_extent_ndims(poDS->dataspace_id);
1048 61 : if (poDS->ndims <= 0)
1049 : {
1050 0 : delete poDS;
1051 0 : return nullptr;
1052 : }
1053 61 : poDS->dims =
1054 61 : static_cast<hsize_t *>(CPLCalloc(poDS->ndims, sizeof(hsize_t)));
1055 61 : poDS->maxdims =
1056 61 : static_cast<hsize_t *>(CPLCalloc(poDS->ndims, sizeof(hsize_t)));
1057 61 : poDS->dimensions = H5Sget_simple_extent_dims(poDS->dataspace_id, poDS->dims,
1058 : poDS->maxdims);
1059 195 : for (int i = 0; i < poDS->dimensions; ++i)
1060 : {
1061 270 : if (poDS->dims[i] >
1062 135 : 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 60 : auto datatype = H5Dget_type(poDS->dataset_id);
1072 60 : poDS->native = H5Tget_native_type(datatype, H5T_DIR_ASCEND);
1073 60 : H5Tclose(datatype);
1074 :
1075 60 : const auto eGDALDataType = poDS->GetDataType(poDS->native);
1076 60 : 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 60 : poDS->SetMetadata(poDS->m_aosMetadata.List());
1094 :
1095 : // Check if the hdf5 is a well known product type
1096 60 : poDS->IdentifyProductType();
1097 :
1098 60 : poDS->m_nYIndex = poDS->IsComplexCSKL1A() ? 0 : poDS->ndims - 2;
1099 60 : poDS->m_nXIndex = poDS->IsComplexCSKL1A() ? 1 : poDS->ndims - 1;
1100 :
1101 60 : if (poDS->IsComplexCSKL1A())
1102 : {
1103 0 : poDS->m_nOtherDimIndex = 2;
1104 : }
1105 60 : else if (poDS->ndims == 3)
1106 : {
1107 15 : poDS->m_nOtherDimIndex = 0;
1108 : }
1109 :
1110 60 : if (HDF5EOSParser::HasHDFEOS(poDS->hGroupID))
1111 : {
1112 40 : HDF5EOSParser oHDFEOSParser;
1113 20 : if (oHDFEOSParser.Parse(poDS->hGroupID))
1114 : {
1115 20 : CPLDebug("HDF5", "Successfully parsed HDFEOS metadata");
1116 40 : HDF5EOSParser::GridDataFieldMetadata oGridDataFieldMetadata;
1117 40 : HDF5EOSParser::SwathFieldMetadata oSwathFieldMetadata;
1118 20 : if (oHDFEOSParser.GetDataModel() ==
1119 6 : HDF5EOSParser::DataModel::GRID &&
1120 6 : oHDFEOSParser.GetGridDataFieldMetadata(
1121 26 : 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 14 : else if (oHDFEOSParser.GetDataModel() ==
1146 14 : HDF5EOSParser::DataModel::SWATH &&
1147 14 : oHDFEOSParser.GetSwathFieldMetadata(
1148 14 : osSubdatasetName.c_str(), oSwathFieldMetadata) &&
1149 : static_cast<int>(
1150 14 : oSwathFieldMetadata.aoDimensions.size()) ==
1151 14 : poDS->ndims &&
1152 42 : oSwathFieldMetadata.iXDim >= 0 &&
1153 14 : oSwathFieldMetadata.iYDim >= 0)
1154 : {
1155 14 : poDS->m_nXIndex = oSwathFieldMetadata.iXDim;
1156 14 : poDS->m_nYIndex = oSwathFieldMetadata.iYDim;
1157 14 : poDS->m_nOtherDimIndex = oSwathFieldMetadata.iOtherDim;
1158 14 : if (!oSwathFieldMetadata.osLongitudeSubdataset.empty())
1159 : {
1160 28 : CPLStringList aosGeolocation;
1161 :
1162 : // Arbitrary
1163 14 : aosGeolocation.AddNameValue("SRS", SRS_WKT_WGS84_LAT_LONG);
1164 : aosGeolocation.AddNameValue(
1165 28 : "X_DATASET", ("HDF5:\"" + osFilename + "\":" +
1166 : oSwathFieldMetadata.osLongitudeSubdataset)
1167 14 : .c_str());
1168 14 : aosGeolocation.AddNameValue("X_BAND", "1");
1169 : aosGeolocation.AddNameValue(
1170 28 : "Y_DATASET", ("HDF5:\"" + osFilename + "\":" +
1171 : oSwathFieldMetadata.osLatitudeSubdataset)
1172 14 : .c_str());
1173 14 : aosGeolocation.AddNameValue("Y_BAND", "1");
1174 : aosGeolocation.AddNameValue(
1175 : "PIXEL_OFFSET",
1176 14 : CPLSPrintf("%d", oSwathFieldMetadata.nPixelOffset));
1177 : aosGeolocation.AddNameValue(
1178 : "PIXEL_STEP",
1179 14 : CPLSPrintf("%d", oSwathFieldMetadata.nPixelStep));
1180 : aosGeolocation.AddNameValue(
1181 : "LINE_OFFSET",
1182 14 : CPLSPrintf("%d", oSwathFieldMetadata.nLineOffset));
1183 : aosGeolocation.AddNameValue(
1184 : "LINE_STEP",
1185 14 : CPLSPrintf("%d", oSwathFieldMetadata.nLineStep));
1186 : // Not totally sure about that
1187 : aosGeolocation.AddNameValue("GEOREFERENCING_CONVENTION",
1188 14 : "PIXEL_CENTER");
1189 14 : poDS->SetMetadata(aosGeolocation.List(), "GEOLOCATION");
1190 : }
1191 : }
1192 : }
1193 : }
1194 :
1195 60 : poDS->nRasterYSize =
1196 60 : poDS->GetYIndex() < 0
1197 60 : ? 1
1198 59 : : static_cast<int>(poDS->dims[poDS->GetYIndex()]); // nRows
1199 60 : poDS->nRasterXSize =
1200 60 : static_cast<int>(poDS->dims[poDS->GetXIndex()]); // nCols
1201 60 : int nBands = 1;
1202 60 : if (poDS->m_nOtherDimIndex >= 0)
1203 : {
1204 15 : nBands = static_cast<int>(poDS->dims[poDS->m_nOtherDimIndex]);
1205 : }
1206 :
1207 120 : CPLStringList aosMetadata;
1208 60 : std::map<std::string, CPLStringList> oMapBandSpecificMetadata;
1209 60 : if (poDS->poH5Objects->nType == H5G_DATASET)
1210 : {
1211 60 : HDF5Dataset::CreateMetadata(poDS->m_hHDF5, poDS->poH5Objects,
1212 : H5G_DATASET, false, aosMetadata);
1213 60 : if (nBands > 1 && poDS->nRasterXSize != nBands &&
1214 14 : poDS->nRasterYSize != nBands)
1215 : {
1216 : // Heuristics to detect non-scalar attributes, that are intended
1217 : // to be attached to a specific band.
1218 28 : const CPLStringList aosMetadataDup(aosMetadata);
1219 52 : for (const auto &[pszKey, pszValue] :
1220 66 : cpl::IterateNameValue(aosMetadataDup))
1221 : {
1222 26 : const hid_t hAttrID = H5Aopen_name(poDS->dataset_id, pszKey);
1223 26 : const hid_t hAttrSpace = H5Aget_space(hAttrID);
1224 47 : if (H5Sget_simple_extent_ndims(hAttrSpace) == 1 &&
1225 21 : H5Sget_simple_extent_npoints(hAttrSpace) == nBands)
1226 : {
1227 : CPLStringList aosTokens(
1228 40 : CSLTokenizeString2(pszValue, " ", 0));
1229 20 : if (aosTokens.size() == nBands)
1230 : {
1231 40 : std::string osAttrName(pszKey);
1232 30 : if (osAttrName.size() > strlen("_coefficients") &&
1233 30 : osAttrName.substr(osAttrName.size() -
1234 10 : strlen("_coefficients")) ==
1235 : "_coefficients")
1236 : {
1237 5 : osAttrName.pop_back();
1238 : }
1239 25 : else if (osAttrName.size() > strlen("_wavelengths") &&
1240 25 : osAttrName.substr(osAttrName.size() -
1241 10 : strlen("_wavelengths")) ==
1242 : "_wavelengths")
1243 : {
1244 5 : osAttrName.pop_back();
1245 : }
1246 15 : else if (osAttrName.size() > strlen("_list") &&
1247 15 : osAttrName.substr(osAttrName.size() -
1248 5 : strlen("_list")) == "_list")
1249 : {
1250 5 : osAttrName.resize(osAttrName.size() -
1251 : strlen("_list"));
1252 : }
1253 20 : oMapBandSpecificMetadata[osAttrName] =
1254 40 : std::move(aosTokens);
1255 20 : aosMetadata.SetNameValue(pszKey, nullptr);
1256 : }
1257 : }
1258 26 : H5Sclose(hAttrSpace);
1259 26 : H5Aclose(hAttrID);
1260 : }
1261 : }
1262 : }
1263 :
1264 60 : poDS->m_nBlockXSize = poDS->GetRasterXSize();
1265 60 : poDS->m_nBlockYSize = 1;
1266 60 : poDS->m_nBandChunkSize = 1;
1267 :
1268 : // Check for chunksize and set it as the blocksize (optimizes read).
1269 60 : const hid_t listid = H5Dget_create_plist(poDS->dataset_id);
1270 60 : if (listid > 0)
1271 : {
1272 60 : if (H5Pget_layout(listid) == H5D_CHUNKED)
1273 : {
1274 12 : hsize_t panChunkDims[3] = {0, 0, 0};
1275 12 : const int nDimSize = H5Pget_chunk(listid, 3, panChunkDims);
1276 12 : CPL_IGNORE_RET_VAL(nDimSize);
1277 12 : CPLAssert(nDimSize == poDS->ndims);
1278 12 : poDS->m_nBlockXSize =
1279 12 : static_cast<int>(panChunkDims[poDS->GetXIndex()]);
1280 12 : if (poDS->GetYIndex() >= 0)
1281 12 : poDS->m_nBlockYSize =
1282 12 : static_cast<int>(panChunkDims[poDS->GetYIndex()]);
1283 12 : if (nBands > 1)
1284 : {
1285 3 : poDS->m_nBandChunkSize =
1286 3 : static_cast<int>(panChunkDims[poDS->m_nOtherDimIndex]);
1287 :
1288 3 : if (poDS->m_nBandChunkSize > 1)
1289 3 : poDS->SetMetadataItem("INTERLEAVE", "PIXEL",
1290 3 : "IMAGE_STRUCTURE");
1291 :
1292 3 : poDS->SetMetadataItem("BAND_CHUNK_SIZE",
1293 : CPLSPrintf("%d", poDS->m_nBandChunkSize),
1294 3 : "IMAGE_STRUCTURE");
1295 : }
1296 : }
1297 :
1298 60 : const int nFilters = H5Pget_nfilters(listid);
1299 66 : for (int i = 0; i < nFilters; ++i)
1300 : {
1301 6 : unsigned int flags = 0;
1302 6 : size_t cd_nelmts = 0;
1303 6 : char szName[64 + 1] = {0};
1304 6 : const auto eFilter = H5Pget_filter(listid, i, &flags, &cd_nelmts,
1305 : nullptr, 64, szName);
1306 6 : if (eFilter == H5Z_FILTER_DEFLATE)
1307 : {
1308 5 : poDS->SetMetadataItem("COMPRESSION", "DEFLATE",
1309 5 : "IMAGE_STRUCTURE");
1310 : }
1311 1 : else if (eFilter == H5Z_FILTER_SZIP)
1312 : {
1313 0 : poDS->SetMetadataItem("COMPRESSION", "SZIP", "IMAGE_STRUCTURE");
1314 : }
1315 : }
1316 :
1317 60 : H5Pclose(listid);
1318 : }
1319 :
1320 187 : for (int i = 0; i < nBands; i++)
1321 : {
1322 : HDF5ImageRasterBand *const poBand =
1323 127 : new HDF5ImageRasterBand(poDS, i + 1, eGDALDataType);
1324 :
1325 127 : poDS->SetBand(i + 1, poBand);
1326 :
1327 127 : if (poDS->poH5Objects->nType == H5G_DATASET)
1328 : {
1329 127 : poBand->SetMetadata(aosMetadata.List());
1330 167 : for (const auto &oIter : oMapBandSpecificMetadata)
1331 : {
1332 40 : poBand->SetMetadataItem(oIter.first.c_str(), oIter.second[i]);
1333 : }
1334 : }
1335 : }
1336 :
1337 60 : if (!poDS->GetMetadata("GEOLOCATION"))
1338 46 : poDS->CreateProjections();
1339 :
1340 : // Setup/check for pam .aux.xml.
1341 60 : poDS->TryLoadXML();
1342 :
1343 : // If the PAM .aux.xml file contains the serialized GEOLOCATION metadata
1344 : // domain, make sure to patch the X_DATASET and Y_DATASET values, when
1345 : // they point to the current file, to use the filename with which we have
1346 : // opened the dataset. Helps in scenarios where the .aux.xml file has
1347 : // been saved with a relative filename for example.
1348 : // Or scenarios like https://github.com/OSGeo/gdal/issues/12824 mixing
1349 : // native Windows and WSL use.
1350 : CSLConstList papszGeoLocationAfterLoadXML =
1351 60 : poDS->GetMetadata("GEOLOCATION");
1352 60 : if (papszGeoLocationAfterLoadXML)
1353 : {
1354 42 : for (const char *pszItem : {"X_DATASET", "Y_DATASET"})
1355 : {
1356 : const char *pszSubdataset =
1357 28 : CSLFetchNameValue(papszGeoLocationAfterLoadXML, pszItem);
1358 56 : if (pszSubdataset && STARTS_WITH(pszSubdataset, "HDF5:\"") &&
1359 28 : strstr(pszSubdataset, CPLGetFilename(osFilename.c_str())))
1360 : {
1361 28 : auto hSubDSInfo = GDALGetSubdatasetInfo(pszSubdataset);
1362 28 : if (hSubDSInfo)
1363 : {
1364 : char *pszOriPath =
1365 28 : GDALSubdatasetInfoGetPathComponent(hSubDSInfo);
1366 28 : if (EQUAL(CPLGetFilename(pszOriPath),
1367 : CPLGetFilename(osFilename.c_str())))
1368 : {
1369 28 : char *pszNewVal = GDALSubdatasetInfoModifyPathComponent(
1370 : hSubDSInfo, osFilename.c_str());
1371 28 : poDS->SetMetadataItem(pszItem, pszNewVal,
1372 28 : "GEOLOCATION");
1373 28 : CPLFree(pszNewVal);
1374 : }
1375 28 : CPLFree(pszOriPath);
1376 28 : GDALDestroySubdatasetInfo(hSubDSInfo);
1377 : }
1378 : }
1379 : }
1380 : }
1381 :
1382 60 : poDS->SetPamFlags(poDS->GetPamFlags() & ~GPF_DIRTY);
1383 :
1384 : // Setup overviews.
1385 60 : poDS->oOvManager.Initialize(poDS, ":::VIRTUAL:::");
1386 :
1387 60 : return poDS;
1388 : }
1389 :
1390 : /************************************************************************/
1391 : /* HDF5ImageDatasetDriverUnload() */
1392 : /************************************************************************/
1393 :
1394 9 : static void HDF5ImageDatasetDriverUnload(GDALDriver *)
1395 : {
1396 9 : HDF5UnloadFileDriver();
1397 9 : }
1398 :
1399 : /************************************************************************/
1400 : /* GDALRegister_HDF5Image() */
1401 : /************************************************************************/
1402 14 : void GDALRegister_HDF5Image()
1403 :
1404 : {
1405 14 : if (!GDAL_CHECK_VERSION("HDF5Image driver"))
1406 0 : return;
1407 :
1408 14 : if (GDALGetDriverByName(HDF5_IMAGE_DRIVER_NAME) != nullptr)
1409 0 : return;
1410 :
1411 14 : GDALDriver *poDriver = new GDALDriver();
1412 :
1413 14 : HDF5ImageDriverSetCommonMetadata(poDriver);
1414 :
1415 14 : poDriver->pfnOpen = HDF5ImageDataset::Open;
1416 14 : poDriver->pfnUnloadDriver = HDF5ImageDatasetDriverUnload;
1417 :
1418 14 : GetGDALDriverManager()->RegisterDriver(poDriver);
1419 : }
1420 :
1421 : /************************************************************************/
1422 : /* CreateODIMH5Projection() */
1423 : /************************************************************************/
1424 :
1425 : // Reference:
1426 : // http://www.knmi.nl/opera/opera3/OPERA_2008_03_WP2.1b_ODIM_H5_v2.1.pdf
1427 : //
1428 : // 4.3.2 where for geographically referenced image Groups
1429 : // We don't use the where_xscale and where_yscale parameters, but recompute them
1430 : // from the lower-left and upper-right coordinates. There's some difference.
1431 : // As all those parameters are linked together, I'm not sure which one should be
1432 : // considered as the reference.
1433 :
1434 0 : CPLErr HDF5ImageDataset::CreateODIMH5Projection()
1435 : {
1436 0 : const char *const pszProj4String = GetMetadataItem("where_projdef");
1437 0 : const char *const pszLL_lon = GetMetadataItem("where_LL_lon");
1438 0 : const char *const pszLL_lat = GetMetadataItem("where_LL_lat");
1439 0 : const char *const pszUR_lon = GetMetadataItem("where_UR_lon");
1440 0 : const char *const pszUR_lat = GetMetadataItem("where_UR_lat");
1441 0 : if (pszProj4String == nullptr || pszLL_lon == nullptr ||
1442 0 : pszLL_lat == nullptr || pszUR_lon == nullptr || pszUR_lat == nullptr)
1443 0 : return CE_Failure;
1444 :
1445 0 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1446 0 : if (m_oSRS.importFromProj4(pszProj4String) != OGRERR_NONE)
1447 0 : return CE_Failure;
1448 :
1449 0 : OGRSpatialReference oSRSWGS84;
1450 0 : oSRSWGS84.SetWellKnownGeogCS("WGS84");
1451 0 : oSRSWGS84.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1452 :
1453 : OGRCoordinateTransformation *poCT =
1454 0 : OGRCreateCoordinateTransformation(&oSRSWGS84, &m_oSRS);
1455 0 : if (poCT == nullptr)
1456 0 : return CE_Failure;
1457 :
1458 : // Reproject corners from long,lat WGS84 to the target SRS.
1459 0 : double dfLLX = CPLAtof(pszLL_lon);
1460 0 : double dfLLY = CPLAtof(pszLL_lat);
1461 0 : double dfURX = CPLAtof(pszUR_lon);
1462 0 : double dfURY = CPLAtof(pszUR_lat);
1463 0 : if (!poCT->Transform(1, &dfLLX, &dfLLY) ||
1464 0 : !poCT->Transform(1, &dfURX, &dfURY))
1465 : {
1466 0 : delete poCT;
1467 0 : return CE_Failure;
1468 : }
1469 0 : delete poCT;
1470 :
1471 : // Compute the geotransform now.
1472 0 : const double dfPixelX = (dfURX - dfLLX) / nRasterXSize;
1473 0 : const double dfPixelY = (dfURY - dfLLY) / nRasterYSize;
1474 :
1475 0 : bHasGeoTransform = true;
1476 0 : m_gt[0] = dfLLX;
1477 0 : m_gt[1] = dfPixelX;
1478 0 : m_gt[2] = 0;
1479 0 : m_gt[3] = dfURY;
1480 0 : m_gt[4] = 0;
1481 0 : m_gt[5] = -dfPixelY;
1482 :
1483 0 : return CE_None;
1484 : }
1485 :
1486 : /************************************************************************/
1487 : /* CreateProjections() */
1488 : /************************************************************************/
1489 46 : CPLErr HDF5ImageDataset::CreateProjections()
1490 : {
1491 46 : switch (iSubdatasetType)
1492 : {
1493 2 : case CSK_PRODUCT:
1494 : {
1495 2 : int productType = PROD_UNKNOWN;
1496 :
1497 2 : if (GetMetadataItem("Product_Type") != nullptr)
1498 : {
1499 : // Get the format's level.
1500 : const char *osMissionLevel =
1501 2 : HDF5Dataset::GetMetadataItem("Product_Type");
1502 :
1503 2 : if (STARTS_WITH_CI(osMissionLevel, "RAW"))
1504 0 : productType = PROD_CSK_L0;
1505 :
1506 2 : if (STARTS_WITH_CI(osMissionLevel, "SSC"))
1507 0 : productType = PROD_CSK_L1A;
1508 :
1509 2 : if (STARTS_WITH_CI(osMissionLevel, "DGM"))
1510 1 : productType = PROD_CSK_L1B;
1511 :
1512 2 : if (STARTS_WITH_CI(osMissionLevel, "GEC"))
1513 1 : productType = PROD_CSK_L1C;
1514 :
1515 2 : if (STARTS_WITH_CI(osMissionLevel, "GTC"))
1516 0 : productType = PROD_CSK_L1D;
1517 : }
1518 :
1519 2 : CaptureCSKGeoTransform(productType);
1520 2 : CaptureCSKGeolocation(productType);
1521 2 : CaptureCSKGCPs(productType);
1522 :
1523 2 : break;
1524 : }
1525 44 : case UNKNOWN_PRODUCT:
1526 : {
1527 44 : constexpr int NBGCPLAT = 100;
1528 44 : constexpr int NBGCPLON = 30;
1529 :
1530 44 : const int nDeltaLat = nRasterYSize / NBGCPLAT;
1531 44 : const int nDeltaLon = nRasterXSize / NBGCPLON;
1532 :
1533 44 : if (nDeltaLat == 0 || nDeltaLon == 0)
1534 36 : return CE_None;
1535 :
1536 : // Create HDF5 Data Hierarchy in a link list.
1537 8 : poH5Objects = HDF5FindDatasetObjects(poH5RootGroup, "Latitude");
1538 8 : if (!poH5Objects)
1539 : {
1540 8 : if (GetMetadataItem("where_projdef") != nullptr)
1541 0 : return CreateODIMH5Projection();
1542 8 : return CE_None;
1543 : }
1544 :
1545 : // The Latitude and Longitude arrays must have a rank of 2 to
1546 : // retrieve GCPs.
1547 0 : if (poH5Objects->nRank != 2 ||
1548 0 : poH5Objects->paDims[0] != static_cast<size_t>(nRasterYSize) ||
1549 0 : poH5Objects->paDims[1] != static_cast<size_t>(nRasterXSize))
1550 : {
1551 0 : return CE_None;
1552 : }
1553 :
1554 : // Retrieve HDF5 data information.
1555 : const hid_t LatitudeDatasetID =
1556 0 : H5Dopen(m_hHDF5, poH5Objects->pszPath);
1557 : // LatitudeDataspaceID = H5Dget_space(dataset_id);
1558 :
1559 0 : poH5Objects = HDF5FindDatasetObjects(poH5RootGroup, "Longitude");
1560 : // GCPs.
1561 0 : if (poH5Objects == nullptr || poH5Objects->nRank != 2 ||
1562 0 : poH5Objects->paDims[0] != static_cast<size_t>(nRasterYSize) ||
1563 0 : poH5Objects->paDims[1] != static_cast<size_t>(nRasterXSize))
1564 : {
1565 0 : if (LatitudeDatasetID > 0)
1566 0 : H5Dclose(LatitudeDatasetID);
1567 0 : return CE_None;
1568 : }
1569 :
1570 : const hid_t LongitudeDatasetID =
1571 0 : H5Dopen(m_hHDF5, poH5Objects->pszPath);
1572 : // LongitudeDataspaceID = H5Dget_space(dataset_id);
1573 :
1574 0 : if (LatitudeDatasetID > 0 && LongitudeDatasetID > 0)
1575 : {
1576 : float *const Latitude =
1577 0 : static_cast<float *>(VSI_MALLOC3_VERBOSE(
1578 : nRasterYSize, nRasterXSize, sizeof(float)));
1579 : float *const Longitude =
1580 0 : static_cast<float *>(VSI_MALLOC3_VERBOSE(
1581 : nRasterYSize, nRasterXSize, sizeof(float)));
1582 0 : if (!Latitude || !Longitude)
1583 : {
1584 0 : CPLFree(Latitude);
1585 0 : CPLFree(Longitude);
1586 0 : H5Dclose(LatitudeDatasetID);
1587 0 : H5Dclose(LongitudeDatasetID);
1588 0 : return CE_Failure;
1589 : }
1590 0 : memset(Latitude, 0,
1591 0 : static_cast<size_t>(nRasterXSize) * nRasterYSize *
1592 : sizeof(float));
1593 0 : memset(Longitude, 0,
1594 0 : static_cast<size_t>(nRasterXSize) * nRasterYSize *
1595 : sizeof(float));
1596 :
1597 : // netCDF convention for nodata
1598 0 : double dfLatNoData = 0;
1599 0 : bool bHasLatNoData = GH5_FetchAttribute(
1600 : LatitudeDatasetID, "_FillValue", dfLatNoData);
1601 :
1602 0 : double dfLongNoData = 0;
1603 0 : bool bHasLongNoData = GH5_FetchAttribute(
1604 : LongitudeDatasetID, "_FillValue", dfLongNoData);
1605 :
1606 0 : H5Dread(LatitudeDatasetID, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
1607 : H5P_DEFAULT, Latitude);
1608 :
1609 0 : H5Dread(LongitudeDatasetID, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
1610 : H5P_DEFAULT, Longitude);
1611 :
1612 0 : m_oSRS.Clear();
1613 0 : m_oGCPSRS.SetWellKnownGeogCS("WGS84");
1614 :
1615 0 : const int nYLimit =
1616 0 : (static_cast<int>(nRasterYSize) / nDeltaLat) * nDeltaLat;
1617 0 : const int nXLimit =
1618 0 : (static_cast<int>(nRasterXSize) / nDeltaLon) * nDeltaLon;
1619 :
1620 : // The original code in
1621 : // https://trac.osgeo.org/gdal/changeset/8066 always add +180 to
1622 : // the longitudes, but without justification I suspect this
1623 : // might be due to handling products crossing the antimeridian.
1624 : // Trying to do it just when needed through a heuristics.
1625 0 : bool bHasLonNearMinus180 = false;
1626 0 : bool bHasLonNearPlus180 = false;
1627 0 : bool bHasLonNearZero = false;
1628 0 : for (int j = 0; j < nYLimit; j += nDeltaLat)
1629 : {
1630 0 : for (int i = 0; i < nXLimit; i += nDeltaLon)
1631 : {
1632 0 : const int iGCP = j * nRasterXSize + i;
1633 0 : if ((bHasLatNoData && static_cast<float>(dfLatNoData) ==
1634 0 : Latitude[iGCP]) ||
1635 0 : (bHasLongNoData &&
1636 0 : static_cast<float>(dfLongNoData) ==
1637 0 : Longitude[iGCP]))
1638 0 : continue;
1639 0 : if (Longitude[iGCP] > 170 && Longitude[iGCP] <= 180)
1640 0 : bHasLonNearPlus180 = true;
1641 0 : if (Longitude[iGCP] < -170 && Longitude[iGCP] >= -180)
1642 0 : bHasLonNearMinus180 = true;
1643 0 : if (fabs(Longitude[iGCP]) < 90)
1644 0 : bHasLonNearZero = true;
1645 : }
1646 : }
1647 :
1648 : // Fill the GCPs list.
1649 : const char *pszShiftGCP =
1650 0 : CPLGetConfigOption("HDF5_SHIFT_GCPX_BY_180", nullptr);
1651 : const bool bAdd180 =
1652 0 : (bHasLonNearPlus180 && bHasLonNearMinus180 &&
1653 0 : !bHasLonNearZero && pszShiftGCP == nullptr) ||
1654 0 : (pszShiftGCP != nullptr && CPLTestBool(pszShiftGCP));
1655 :
1656 0 : for (int j = 0; j < nYLimit; j += nDeltaLat)
1657 : {
1658 0 : for (int i = 0; i < nXLimit; i += nDeltaLon)
1659 : {
1660 0 : const int iGCP = j * nRasterXSize + i;
1661 0 : if ((bHasLatNoData && static_cast<float>(dfLatNoData) ==
1662 0 : Latitude[iGCP]) ||
1663 0 : (bHasLongNoData &&
1664 0 : static_cast<float>(dfLongNoData) ==
1665 0 : Longitude[iGCP]))
1666 0 : continue;
1667 0 : double dfGCPX = static_cast<double>(Longitude[iGCP]);
1668 0 : if (bAdd180)
1669 0 : dfGCPX += 180.0;
1670 0 : const double dfGCPY =
1671 0 : static_cast<double>(Latitude[iGCP]);
1672 :
1673 0 : m_aoGCPs.emplace_back("", "", i + 0.5, j + 0.5, dfGCPX,
1674 0 : dfGCPY);
1675 : }
1676 : }
1677 :
1678 0 : CPLFree(Latitude);
1679 0 : CPLFree(Longitude);
1680 : }
1681 :
1682 0 : if (LatitudeDatasetID > 0)
1683 0 : H5Dclose(LatitudeDatasetID);
1684 0 : if (LongitudeDatasetID > 0)
1685 0 : H5Dclose(LongitudeDatasetID);
1686 :
1687 0 : break;
1688 : }
1689 : }
1690 :
1691 2 : return CE_None;
1692 : }
1693 :
1694 : /************************************************************************/
1695 : /* GetMetadataItem() */
1696 : /************************************************************************/
1697 :
1698 123 : const char *HDF5ImageDataset::GetMetadataItem(const char *pszName,
1699 : const char *pszDomain)
1700 : {
1701 123 : if (pszDomain && EQUAL(pszDomain, "__DEBUG__") &&
1702 66 : EQUAL(pszName, "WholeBandChunkOptim"))
1703 : {
1704 66 : switch (m_eWholeBandChunkOptim)
1705 : {
1706 4 : case WBC_DETECTION_IN_PROGRESS:
1707 4 : return "DETECTION_IN_PROGRESS";
1708 3 : case WBC_DISABLED:
1709 3 : return "DISABLED";
1710 59 : case WBC_ENABLED:
1711 59 : return "ENABLED";
1712 : }
1713 : }
1714 57 : return GDALPamDataset::GetMetadataItem(pszName, pszDomain);
1715 : }
1716 :
1717 : /************************************************************************/
1718 : /* GetSpatialRef() */
1719 : /************************************************************************/
1720 :
1721 12 : const OGRSpatialReference *HDF5ImageDataset::GetSpatialRef() const
1722 : {
1723 12 : if (!m_oSRS.IsEmpty())
1724 11 : return &m_oSRS;
1725 1 : return GDALPamDataset::GetSpatialRef();
1726 : }
1727 :
1728 : /************************************************************************/
1729 : /* GetGCPCount() */
1730 : /************************************************************************/
1731 :
1732 4 : int HDF5ImageDataset::GetGCPCount()
1733 :
1734 : {
1735 4 : if (!m_aoGCPs.empty())
1736 1 : return static_cast<int>(m_aoGCPs.size());
1737 :
1738 3 : return GDALPamDataset::GetGCPCount();
1739 : }
1740 :
1741 : /************************************************************************/
1742 : /* GetGCPSpatialRef() */
1743 : /************************************************************************/
1744 :
1745 1 : const OGRSpatialReference *HDF5ImageDataset::GetGCPSpatialRef() const
1746 :
1747 : {
1748 1 : if (!m_aoGCPs.empty() && !m_oGCPSRS.IsEmpty())
1749 1 : return &m_oGCPSRS;
1750 :
1751 0 : return GDALPamDataset::GetGCPSpatialRef();
1752 : }
1753 :
1754 : /************************************************************************/
1755 : /* GetGCPs() */
1756 : /************************************************************************/
1757 :
1758 3 : const GDAL_GCP *HDF5ImageDataset::GetGCPs()
1759 : {
1760 3 : if (!m_aoGCPs.empty())
1761 1 : return gdal::GCP::c_ptr(m_aoGCPs);
1762 :
1763 2 : return GDALPamDataset::GetGCPs();
1764 : }
1765 :
1766 : /************************************************************************/
1767 : /* GetGeoTransform() */
1768 : /************************************************************************/
1769 :
1770 8 : CPLErr HDF5ImageDataset::GetGeoTransform(GDALGeoTransform >) const
1771 : {
1772 8 : if (bHasGeoTransform)
1773 : {
1774 7 : gt = m_gt;
1775 7 : return CE_None;
1776 : }
1777 :
1778 1 : return GDALPamDataset::GetGeoTransform(gt);
1779 : }
1780 :
1781 : /************************************************************************/
1782 : /* IdentifyProductType() */
1783 : /************************************************************************/
1784 :
1785 : /**
1786 : * Identify if the subdataset has a known product format
1787 : * It stores a product identifier in iSubdatasetType,
1788 : * UNKNOWN_PRODUCT, if it isn't a recognizable format.
1789 : */
1790 60 : void HDF5ImageDataset::IdentifyProductType()
1791 : {
1792 60 : iSubdatasetType = UNKNOWN_PRODUCT;
1793 :
1794 : // COSMO-SKYMED
1795 :
1796 : // Get the Mission Id as a char *, because the field may not exist.
1797 60 : const char *const pszMissionId = HDF5Dataset::GetMetadataItem("Mission_ID");
1798 :
1799 : // If there is a Mission_ID field.
1800 60 : if (pszMissionId != nullptr && strstr(GetDescription(), "QLK") == nullptr)
1801 : {
1802 : // Check if the mission type is CSK, KMPS or CSG.
1803 : // KMPS: Komsat-5 is Korean mission with a SAR instrument.
1804 : // CSG: Cosmo Skymed 2nd Generation
1805 2 : if (EQUAL(pszMissionId, "CSK") || EQUAL(pszMissionId, "KMPS") ||
1806 0 : EQUAL(pszMissionId, "CSG"))
1807 : {
1808 2 : iSubdatasetType = CSK_PRODUCT;
1809 :
1810 2 : if (GetMetadataItem("Product_Type") != nullptr)
1811 : {
1812 : // Get the format's level.
1813 : const char *osMissionLevel =
1814 2 : HDF5Dataset::GetMetadataItem("Product_Type");
1815 :
1816 2 : if (STARTS_WITH_CI(osMissionLevel, "RAW"))
1817 0 : iCSKProductType = PROD_CSK_L0;
1818 :
1819 2 : if (STARTS_WITH_CI(osMissionLevel, "SCS"))
1820 0 : iCSKProductType = PROD_CSK_L1A;
1821 :
1822 2 : if (STARTS_WITH_CI(osMissionLevel, "DGM"))
1823 1 : iCSKProductType = PROD_CSK_L1B;
1824 :
1825 2 : if (STARTS_WITH_CI(osMissionLevel, "GEC"))
1826 1 : iCSKProductType = PROD_CSK_L1C;
1827 :
1828 2 : if (STARTS_WITH_CI(osMissionLevel, "GTC"))
1829 0 : iCSKProductType = PROD_CSK_L1D;
1830 : }
1831 : }
1832 : }
1833 60 : }
1834 :
1835 : /************************************************************************/
1836 : /* CaptureCSKGeolocation() */
1837 : /************************************************************************/
1838 :
1839 : /**
1840 : * Captures Geolocation information from a COSMO-SKYMED
1841 : * file.
1842 : * The geoid will always be WGS84
1843 : * The projection type may be UTM or UPS, depending on the
1844 : * latitude from the center of the image.
1845 : * @param iProductType type of CSK subproduct, see HDF5CSKProduct
1846 : */
1847 2 : void HDF5ImageDataset::CaptureCSKGeolocation(int iProductType)
1848 : {
1849 : // Set the ellipsoid to WGS84.
1850 2 : m_oSRS.SetWellKnownGeogCS("WGS84");
1851 :
1852 2 : if (iProductType == PROD_CSK_L1C || iProductType == PROD_CSK_L1D)
1853 : {
1854 1 : double *dfProjFalseEastNorth = nullptr;
1855 1 : double *dfProjScaleFactor = nullptr;
1856 1 : double *dfCenterCoord = nullptr;
1857 :
1858 : // Check if all the metadata attributes are present.
1859 1 : if (HDF5ReadDoubleAttr("Map Projection False East-North",
1860 1 : &dfProjFalseEastNorth) == CE_Failure ||
1861 1 : HDF5ReadDoubleAttr("Map Projection Scale Factor",
1862 1 : &dfProjScaleFactor) == CE_Failure ||
1863 1 : HDF5ReadDoubleAttr("Map Projection Centre", &dfCenterCoord) ==
1864 2 : CE_Failure ||
1865 1 : GetMetadataItem("Projection_ID") == nullptr)
1866 : {
1867 0 : m_oSRS.Clear();
1868 0 : m_oGCPSRS.Clear();
1869 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1870 : "The CSK hdf5 file geolocation information is "
1871 : "malformed");
1872 : }
1873 : else
1874 : {
1875 : // Fetch projection Type.
1876 2 : CPLString osProjectionID = GetMetadataItem("Projection_ID");
1877 :
1878 : // If the projection is UTM.
1879 1 : if (EQUAL(osProjectionID, "UTM"))
1880 : {
1881 : // @TODO: use SetUTM
1882 1 : m_oSRS.SetProjCS(SRS_PT_TRANSVERSE_MERCATOR);
1883 1 : m_oSRS.SetTM(dfCenterCoord[0], dfCenterCoord[1],
1884 : dfProjScaleFactor[0], dfProjFalseEastNorth[0],
1885 1 : dfProjFalseEastNorth[1]);
1886 : }
1887 : else
1888 : {
1889 : // TODO: No UPS projected files to test!
1890 : // If the projection is UPS.
1891 0 : if (EQUAL(osProjectionID, "UPS"))
1892 : {
1893 0 : m_oSRS.SetProjCS(SRS_PT_POLAR_STEREOGRAPHIC);
1894 0 : m_oSRS.SetPS(dfCenterCoord[0], dfCenterCoord[1],
1895 : dfProjScaleFactor[0], dfProjFalseEastNorth[0],
1896 0 : dfProjFalseEastNorth[1]);
1897 : }
1898 : }
1899 :
1900 1 : CPLFree(dfCenterCoord);
1901 1 : CPLFree(dfProjScaleFactor);
1902 1 : CPLFree(dfProjFalseEastNorth);
1903 1 : }
1904 : }
1905 : else
1906 : {
1907 1 : m_oGCPSRS = m_oSRS;
1908 : }
1909 2 : }
1910 :
1911 : /************************************************************************/
1912 : /* CaptureCSKGeoTransform() */
1913 : /************************************************************************/
1914 :
1915 : /**
1916 : * Get Geotransform information for COSMO-SKYMED files
1917 : * In case of success it stores the transformation
1918 : * in m_gt. In case of failure it doesn't
1919 : * modify m_gt
1920 : * @param iProductType type of CSK subproduct, see HDF5CSKProduct
1921 : */
1922 2 : void HDF5ImageDataset::CaptureCSKGeoTransform(int iProductType)
1923 : {
1924 2 : const char *pszSubdatasetName = GetSubdatasetName();
1925 :
1926 2 : bHasGeoTransform = false;
1927 : // If the product level is not L1C or L1D then
1928 : // it doesn't have a valid projection.
1929 2 : if (iProductType == PROD_CSK_L1C || iProductType == PROD_CSK_L1D)
1930 : {
1931 : // If there is a subdataset.
1932 1 : if (pszSubdatasetName != nullptr)
1933 : {
1934 2 : CPLString osULPath = pszSubdatasetName;
1935 1 : osULPath += "/Top Left East-North";
1936 :
1937 2 : CPLString osLineSpacingPath = pszSubdatasetName;
1938 1 : osLineSpacingPath += "/Line Spacing";
1939 :
1940 2 : CPLString osColumnSpacingPath = pszSubdatasetName;
1941 1 : osColumnSpacingPath += "/Column Spacing";
1942 :
1943 1 : double *pdOutUL = nullptr;
1944 1 : double *pdLineSpacing = nullptr;
1945 1 : double *pdColumnSpacing = nullptr;
1946 :
1947 : // If it could find the attributes on the metadata.
1948 1 : if (HDF5ReadDoubleAttr(osULPath.c_str(), &pdOutUL) == CE_Failure ||
1949 1 : HDF5ReadDoubleAttr(osLineSpacingPath.c_str(), &pdLineSpacing) ==
1950 2 : CE_Failure ||
1951 1 : HDF5ReadDoubleAttr(osColumnSpacingPath.c_str(),
1952 : &pdColumnSpacing) == CE_Failure)
1953 : {
1954 0 : bHasGeoTransform = false;
1955 : }
1956 : else
1957 : {
1958 : // geotransform[1] : width of pixel
1959 : // geotransform[4] : rotational coefficient, zero for north up
1960 : // images.
1961 : // geotransform[2] : rotational coefficient, zero for north up
1962 : // images.
1963 : // geotransform[5] : height of pixel (but negative)
1964 :
1965 1 : m_gt[0] = pdOutUL[0];
1966 1 : m_gt[1] = pdLineSpacing[0];
1967 1 : m_gt[2] = 0;
1968 1 : m_gt[3] = pdOutUL[1];
1969 1 : m_gt[4] = 0;
1970 1 : m_gt[5] = -pdColumnSpacing[0];
1971 :
1972 1 : CPLFree(pdOutUL);
1973 1 : CPLFree(pdLineSpacing);
1974 1 : CPLFree(pdColumnSpacing);
1975 :
1976 1 : bHasGeoTransform = true;
1977 : }
1978 : }
1979 : }
1980 2 : }
1981 :
1982 : /************************************************************************/
1983 : /* CaptureCSKGCPs() */
1984 : /************************************************************************/
1985 :
1986 : /**
1987 : * Retrieves and stores the GCPs from a COSMO-SKYMED dataset
1988 : * It only retrieves the GCPs for L0, L1A and L1B products
1989 : * for L1C and L1D products, geotransform is provided.
1990 : * The GCPs provided will be the Image's corners.
1991 : * @param iProductType type of CSK product @see HDF5CSKProductEnum
1992 : */
1993 2 : void HDF5ImageDataset::CaptureCSKGCPs(int iProductType)
1994 : {
1995 : // Only retrieve GCPs for L0,L1A and L1B products.
1996 2 : if (iProductType == PROD_CSK_L0 || iProductType == PROD_CSK_L1A ||
1997 : iProductType == PROD_CSK_L1B)
1998 : {
1999 10 : CPLString osCornerName[4];
2000 1 : double pdCornerPixel[4] = {0.0, 0.0, 0.0, 0.0};
2001 1 : double pdCornerLine[4] = {0.0, 0.0, 0.0, 0.0};
2002 :
2003 1 : const char *const pszSubdatasetName = GetSubdatasetName();
2004 :
2005 : // Load the subdataset name first.
2006 5 : for (int i = 0; i < 4; i++)
2007 4 : osCornerName[i] = pszSubdatasetName;
2008 :
2009 : // Load the attribute name, and raster coordinates for
2010 : // all the corners.
2011 1 : osCornerName[0] += "/Top Left Geodetic Coordinates";
2012 1 : pdCornerPixel[0] = 0;
2013 1 : pdCornerLine[0] = 0;
2014 :
2015 1 : osCornerName[1] += "/Top Right Geodetic Coordinates";
2016 1 : pdCornerPixel[1] = GetRasterXSize();
2017 1 : pdCornerLine[1] = 0;
2018 :
2019 1 : osCornerName[2] += "/Bottom Left Geodetic Coordinates";
2020 1 : pdCornerPixel[2] = 0;
2021 1 : pdCornerLine[2] = GetRasterYSize();
2022 :
2023 1 : osCornerName[3] += "/Bottom Right Geodetic Coordinates";
2024 1 : pdCornerPixel[3] = GetRasterXSize();
2025 1 : pdCornerLine[3] = GetRasterYSize();
2026 :
2027 : // For all the image's corners.
2028 5 : for (int i = 0; i < 4; i++)
2029 : {
2030 4 : double *pdCornerCoordinates = nullptr;
2031 :
2032 : // Retrieve the attributes.
2033 4 : if (HDF5ReadDoubleAttr(osCornerName[i].c_str(),
2034 4 : &pdCornerCoordinates) == CE_Failure)
2035 : {
2036 0 : CPLError(CE_Failure, CPLE_OpenFailed,
2037 : "Error retrieving CSK GCPs");
2038 0 : m_aoGCPs.clear();
2039 0 : break;
2040 : }
2041 :
2042 0 : m_aoGCPs.emplace_back(osCornerName[i].c_str(), "", pdCornerPixel[i],
2043 4 : pdCornerLine[i],
2044 4 : /* X = */ pdCornerCoordinates[1],
2045 : /* Y = */ pdCornerCoordinates[0],
2046 4 : /* Z = */ pdCornerCoordinates[2]);
2047 :
2048 : // Free the returned coordinates.
2049 4 : CPLFree(pdCornerCoordinates);
2050 : }
2051 : }
2052 2 : }
|
