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