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