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