Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Hierarchical Data Format Release 4 (HDF4)
4 : * Purpose: Read subdatasets of HDF4 file.
5 : * This driver initially based on code supplied by Markus Neteler
6 : * Author: Andrey Kiselev, dron@ak4719.spb.edu
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2002, Andrey Kiselev <dron@ak4719.spb.edu>
10 : * Copyright (c) 2009-2013, Even Rouault <even dot rouault at spatialys.com>
11 : *
12 : * SPDX-License-Identifier: MIT
13 : ****************************************************************************/
14 :
15 : #if defined(_WIN32) && !defined(NOMINMAX)
16 : // min/max are defined here on Windows, so block them.
17 : // TODO: Move this to someplace more appropriate.
18 : #define NOMINMAX
19 : #endif
20 :
21 : #include <string.h>
22 : #include <math.h>
23 :
24 : #include "cpl_multiproc.h"
25 : #include "cpl_string.h"
26 : #include "gdal_frmts.h"
27 : #include "gdal_priv.h"
28 : #include "ogr_spatialref.h"
29 :
30 : #include "hdf.h"
31 : #include "mfhdf.h"
32 :
33 : #include "HdfEosDef.h"
34 :
35 : #include "hdf4compat.h"
36 : #include "hdf4dataset.h"
37 :
38 : #include "nasakeywordhandler.h"
39 :
40 : #include "hdf4drivercore.h"
41 :
42 : #include <algorithm>
43 :
44 : constexpr int HDF4_SDS_MAXNAMELEN = 65;
45 :
46 : extern const char *const pszGDALSignature;
47 :
48 : // Signature to recognize files written by GDAL.
49 : const char *const pszGDALSignature =
50 : "Created with GDAL (http://www.remotesensing.org/gdal/)";
51 :
52 : extern CPLMutex *hHDF4Mutex;
53 :
54 : constexpr int N_BUF_SIZE = 8192;
55 :
56 : /************************************************************************/
57 : /* ==================================================================== */
58 : /* List of HDF-EOS Swath product types. */
59 : /* ==================================================================== */
60 : /************************************************************************/
61 :
62 : enum HDF4EOSProduct
63 : {
64 : PROD_UNKNOWN,
65 : PROD_ASTER_L1A,
66 : PROD_ASTER_L1B,
67 : PROD_ASTER_L2,
68 : PROD_ASTER_L3,
69 : PROD_AST14DEM,
70 : PROD_MODIS_L1B,
71 : PROD_MODIS_L2
72 : };
73 :
74 : /************************************************************************/
75 : /* ==================================================================== */
76 : /* HDF4ImageDataset */
77 : /* ==================================================================== */
78 : /************************************************************************/
79 :
80 : constexpr int N_COLOR_ENTRIES = 256;
81 :
82 : class HDF4ImageDataset final : public HDF4Dataset
83 : {
84 : friend class HDF4ImageRasterBand;
85 :
86 : char *pszFilename;
87 : int32 hHDF4;
88 : int32 iGR;
89 : int32 iPal;
90 : int32 iDataset;
91 : int32 iRank;
92 : int32 iNumType;
93 : int32 nAttrs;
94 : int32 iInterlaceMode;
95 : int32 iPalInterlaceMode;
96 : int32 iPalDataType;
97 : int32 nComps;
98 : int32 nPalEntries;
99 : int32 aiDimSizes[H4_MAX_VAR_DIMS];
100 : int iXDim;
101 : int iYDim;
102 : int iBandDim;
103 : int i4Dim;
104 : int nBandCount;
105 : char **papszLocalMetadata{};
106 : uint8 aiPaletteData[N_COLOR_ENTRIES][3]; // XXX: Static array for now
107 : char szName[HDF4_SDS_MAXNAMELEN];
108 : char *pszSubdatasetName;
109 : char *pszFieldName;
110 :
111 : GDALColorTable *poColorTable;
112 :
113 : OGRSpatialReference m_oSRS{};
114 : OGRSpatialReference m_oGCPSRS{};
115 : bool bHasGeoTransform;
116 : GDALGeoTransform m_gt{};
117 : std::vector<gdal::GCP> m_aoGCPs{};
118 :
119 : HDF4DatasetType iDatasetType;
120 :
121 : int32 iSDS;
122 :
123 : int nBlockPreferredXSize;
124 : int nBlockPreferredYSize;
125 : bool bReadTile;
126 :
127 : void ToGeoref(double *, double *);
128 : void GetImageDimensions(char *);
129 : void GetSwatAttrs(int32);
130 : void GetGridAttrs(int32 hGD);
131 : void CaptureNRLGeoTransform(void);
132 : void CaptureL1GMTLInfo(void);
133 : void CaptureCoastwatchGCTPInfo(void);
134 : void ProcessModisSDSGeolocation(void);
135 : int ProcessSwathGeolocation(int32, char **);
136 :
137 : static long USGSMnemonicToCode(const char *);
138 : static void ReadCoordinates(const char *, double *, double *);
139 :
140 : CPL_DISALLOW_COPY_ASSIGN(HDF4ImageDataset)
141 :
142 : public:
143 : HDF4ImageDataset();
144 : ~HDF4ImageDataset() override;
145 :
146 : static GDALDataset *Open(GDALOpenInfo *);
147 : static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
148 : int nBandsIn, GDALDataType eType,
149 : char **papszParamList);
150 : CPLErr FlushCache(bool bAtClosing) override;
151 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
152 : CPLErr SetGeoTransform(const GDALGeoTransform >) override;
153 : const OGRSpatialReference *GetSpatialRef() const override;
154 : CPLErr SetSpatialRef(const OGRSpatialReference *poSRS) override;
155 : int GetGCPCount() override;
156 : const OGRSpatialReference *GetGCPSpatialRef() const override;
157 : const GDAL_GCP *GetGCPs() override;
158 : };
159 :
160 : /************************************************************************/
161 : /* ==================================================================== */
162 : /* HDF4ImageRasterBand */
163 : /* ==================================================================== */
164 : /************************************************************************/
165 :
166 : class HDF4ImageRasterBand final : public GDALPamRasterBand
167 : {
168 : friend class HDF4ImageDataset;
169 :
170 : bool bNoDataSet{false};
171 : double dfNoDataValue{-9999.0};
172 :
173 : bool bHaveScale{false};
174 : bool bHaveOffset{false};
175 : double dfScale{1.0};
176 : double dfOffset{0.0};
177 :
178 : CPLString osUnitType{};
179 :
180 : CPL_DISALLOW_COPY_ASSIGN(HDF4ImageRasterBand)
181 :
182 : public:
183 : HDF4ImageRasterBand(HDF4ImageDataset *, int, GDALDataType);
184 :
185 : CPLErr IReadBlock(int, int, void *) override;
186 : CPLErr IWriteBlock(int, int, void *) override;
187 : GDALColorInterp GetColorInterpretation() override;
188 : GDALColorTable *GetColorTable() override;
189 : double GetNoDataValue(int *) override;
190 : CPLErr SetNoDataValue(double) override;
191 : double GetOffset(int *pbSuccess) override;
192 : double GetScale(int *pbSuccess) override;
193 : const char *GetUnitType() override;
194 : };
195 :
196 : /************************************************************************/
197 : /* HDF4ImageRasterBand() */
198 : /************************************************************************/
199 :
200 544 : HDF4ImageRasterBand::HDF4ImageRasterBand(HDF4ImageDataset *poDSIn, int nBandIn,
201 544 : GDALDataType eType)
202 : {
203 544 : poDS = poDSIn;
204 544 : nBand = nBandIn;
205 544 : eDataType = eType;
206 :
207 544 : nBlockXSize = poDSIn->GetRasterXSize();
208 :
209 : // Aim for a block of about 1000000 pixels. Chunking up substantially
210 : // improves performance in some situations. For now we only chunk up for
211 : // SDS and EOS based datasets since other variations haven't been
212 : // tested. #2208
213 544 : if (poDSIn->iDatasetType == HDF4_SDS || poDSIn->iDatasetType == HDF4_EOS)
214 : {
215 : const int nChunkSize =
216 541 : atoi(CPLGetConfigOption("HDF4_BLOCK_PIXELS", "1000000"));
217 :
218 541 : nBlockYSize = nChunkSize / poDSIn->GetRasterXSize();
219 541 : nBlockYSize =
220 541 : std::max(1, std::min(nBlockYSize, poDSIn->GetRasterYSize()));
221 : }
222 : else
223 : {
224 3 : nBlockYSize = 1;
225 : }
226 :
227 : // HDF4_EOS:EOS_GRID case. We ensure that the block size matches
228 : // the raster width, as the IReadBlock() code can only handle multiple
229 : // blocks per row.
230 544 : if (poDSIn->nBlockPreferredXSize == nBlockXSize &&
231 0 : poDSIn->nBlockPreferredYSize > 0)
232 : {
233 0 : if (poDSIn->nBlockPreferredYSize == 1)
234 : {
235 : // Avoid defaulting to tile reading when the preferred height is 1
236 : // as it leads to very poor performance with:
237 : // ftp://e4ftl01u.ecs.nasa.gov/MODIS_Composites/MOLT/MOD13Q1.005/2006.06.10/MOD13Q1.A2006161.h21v13.005.2008234103220.hd
238 0 : poDSIn->bReadTile = false;
239 : }
240 : else
241 : {
242 0 : nBlockYSize = poDSIn->nBlockPreferredYSize;
243 : }
244 : }
245 :
246 : /* -------------------------------------------------------------------- */
247 : /* We need to avoid using the tile based api if we aren't */
248 : /* matching the tile size. (#4672) */
249 : /* -------------------------------------------------------------------- */
250 544 : if (nBlockXSize != poDSIn->nBlockPreferredXSize ||
251 0 : nBlockYSize != poDSIn->nBlockPreferredYSize)
252 : {
253 544 : poDSIn->bReadTile = false;
254 : }
255 544 : }
256 :
257 : /************************************************************************/
258 : /* IReadBlock() */
259 : /************************************************************************/
260 :
261 62 : CPLErr HDF4ImageRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
262 : void *pImage)
263 : {
264 62 : CPLAssert(nBlockXOff == 0);
265 62 : HDF4ImageDataset *poGDS = cpl::down_cast<HDF4ImageDataset *>(poDS);
266 :
267 124 : CPLMutexHolderD(&hHDF4Mutex);
268 :
269 62 : if (poGDS->eAccess == GA_Update)
270 : {
271 0 : memset(pImage, 0,
272 0 : nBlockXSize * nBlockYSize * GDALGetDataTypeSizeBytes(eDataType));
273 0 : return CE_None;
274 : }
275 :
276 : /* -------------------------------------------------------------------- */
277 : /* Work out some block oriented details. */
278 : /* -------------------------------------------------------------------- */
279 62 : const int nYOff = nBlockYOff * nBlockYSize;
280 : const int nYSize =
281 62 : std::min(nYOff + nBlockYSize, poDS->GetRasterYSize()) - nYOff;
282 :
283 : /* -------------------------------------------------------------------- */
284 : /* HDF files with external data files, such as some landsat */
285 : /* products (eg. data/hdf/L1G) need to be told what directory */
286 : /* to look in to find the external files. Normally this is the */
287 : /* directory holding the hdf file. */
288 : /* -------------------------------------------------------------------- */
289 62 : HXsetdir(CPLGetPathSafe(poGDS->pszFilename).c_str());
290 :
291 : /* -------------------------------------------------------------------- */
292 : /* Handle different configurations. */
293 : /* -------------------------------------------------------------------- */
294 62 : CPLErr eErr = CE_None;
295 62 : int32 aiStart[H4_MAX_NC_DIMS] = {};
296 62 : int32 aiEdges[H4_MAX_NC_DIMS] = {};
297 :
298 62 : switch (poGDS->iDatasetType)
299 : {
300 52 : case HDF4_SDS:
301 : {
302 : // We avoid doing SDselect() / SDendaccess() for each block access
303 : // as this is very slow when zlib compression is used.
304 :
305 52 : if (poGDS->iSDS == FAIL)
306 52 : poGDS->iSDS = SDselect(poGDS->hSD, poGDS->iDataset);
307 :
308 : /* HDF rank:
309 : A rank 2 dataset is an image read in scan-line order (2D).
310 : A rank 3 dataset is a series of images which are read in
311 : an image at a time to form a volume.
312 : A rank 4 dataset may be thought of as a series of volumes.
313 :
314 : The "aiStart" array specifies the multi-dimensional index of the
315 : starting corner of the hyperslab to read. The values are zero
316 : based.
317 :
318 : The "edge" array specifies the number of values to read along
319 : each dimension of the hyperslab.
320 :
321 : The "iStride" array allows for sub-sampling along each
322 : dimension. If a iStride value is specified for a dimension,
323 : that many values will be skipped over when reading along that
324 : dimension. Specifying iStride = NULL in the C interface or
325 : iStride = 1 in either interface specifies contiguous reading
326 : of data. If the iStride values are set to 0, SDreaddata
327 : returns FAIL (or -1). No matter what iStride value is
328 : provided, data is always placed contiguously in buffer.
329 : */
330 52 : switch (poGDS->iRank)
331 : {
332 0 : case 4: // 4Dim: volume-time
333 : // FIXME: needs sample file. Does not work currently.
334 0 : aiStart[3] = 0; // range: 0--aiDimSizes[3]-1
335 0 : aiEdges[3] = 1;
336 0 : aiStart[2] = 0; // range: 0--aiDimSizes[2]-1
337 0 : aiEdges[2] = 1;
338 0 : aiStart[1] = nYOff;
339 0 : aiEdges[1] = nYSize;
340 0 : aiStart[0] = nBlockXOff;
341 0 : aiEdges[0] = nBlockXSize;
342 0 : break;
343 25 : case 3: // 3Dim: volume
344 25 : aiStart[poGDS->iBandDim] = nBand - 1;
345 25 : aiEdges[poGDS->iBandDim] = 1;
346 :
347 25 : aiStart[poGDS->iYDim] = nYOff;
348 25 : aiEdges[poGDS->iYDim] = nYSize;
349 :
350 25 : aiStart[poGDS->iXDim] = nBlockXOff;
351 25 : aiEdges[poGDS->iXDim] = nBlockXSize;
352 25 : break;
353 27 : case 2: // 2Dim: rows/cols
354 27 : aiStart[poGDS->iYDim] = nYOff;
355 27 : aiEdges[poGDS->iYDim] = nYSize;
356 :
357 27 : aiStart[poGDS->iXDim] = nBlockXOff;
358 27 : aiEdges[poGDS->iXDim] = nBlockXSize;
359 27 : break;
360 0 : case 1: // 1Dim:
361 0 : aiStart[poGDS->iXDim] = nBlockXOff;
362 0 : aiEdges[poGDS->iXDim] = nBlockXSize;
363 0 : break;
364 : }
365 :
366 : // Read HDF SDS array
367 52 : if (SDreaddata(poGDS->iSDS, aiStart, nullptr, aiEdges, pImage) < 0)
368 : {
369 0 : CPLError(CE_Failure, CPLE_AppDefined,
370 : "SDreaddata() failed for block.");
371 0 : eErr = CE_Failure;
372 : }
373 :
374 : // SDendaccess( l_iSDS );
375 : }
376 52 : break;
377 :
378 10 : case HDF4_GR:
379 : {
380 : const int nDataTypeSize =
381 10 : GDALGetDataTypeSizeBytes(poGDS->GetDataType(poGDS->iNumType));
382 20 : GByte *pbBuffer = static_cast<GByte *>(CPLMalloc(
383 10 : nBlockXSize * nBlockYSize * poGDS->iRank * nDataTypeSize));
384 :
385 10 : aiStart[poGDS->iYDim] = nYOff;
386 10 : aiEdges[poGDS->iYDim] = nYSize;
387 :
388 10 : aiStart[poGDS->iXDim] = nBlockXOff;
389 10 : aiEdges[poGDS->iXDim] = nBlockXSize;
390 :
391 10 : if (GRreadimage(poGDS->iGR, aiStart, nullptr, aiEdges, pbBuffer) <
392 : 0)
393 : {
394 0 : CPLError(CE_Failure, CPLE_AppDefined,
395 : "GRreaddata() failed for block.");
396 0 : eErr = CE_Failure;
397 : }
398 : else
399 : {
400 10 : for (int i = 0, j = (nBand - 1) * nDataTypeSize;
401 110 : i < nBlockXSize * nDataTypeSize;
402 100 : i += nDataTypeSize, j += poGDS->nBands * nDataTypeSize)
403 100 : memcpy(reinterpret_cast<GByte *>(pImage) + i, pbBuffer + j,
404 : nDataTypeSize);
405 : }
406 :
407 10 : CPLFree(pbBuffer);
408 : }
409 10 : break;
410 :
411 0 : case HDF4_EOS:
412 : {
413 0 : switch (poGDS->iSubdatasetType)
414 : {
415 0 : case H4ST_EOS_GRID:
416 : {
417 : const int32 hGD =
418 0 : GDattach(poGDS->hHDF4, poGDS->pszSubdatasetName);
419 0 : switch (poGDS->iRank)
420 : {
421 0 : case 4: // 4Dim: volume
422 0 : aiStart[poGDS->i4Dim] =
423 0 : (nBand - 1) /
424 0 : poGDS->aiDimSizes[poGDS->iBandDim];
425 0 : aiEdges[poGDS->i4Dim] = 1;
426 :
427 0 : aiStart[poGDS->iBandDim] =
428 0 : (nBand - 1) %
429 0 : poGDS->aiDimSizes[poGDS->iBandDim];
430 0 : aiEdges[poGDS->iBandDim] = 1;
431 :
432 0 : aiStart[poGDS->iYDim] = nYOff;
433 0 : aiEdges[poGDS->iYDim] = nYSize;
434 :
435 0 : aiStart[poGDS->iXDim] = nBlockXOff;
436 0 : aiEdges[poGDS->iXDim] = nBlockXSize;
437 0 : break;
438 0 : case 3: // 3Dim: volume
439 0 : aiStart[poGDS->iBandDim] = nBand - 1;
440 0 : aiEdges[poGDS->iBandDim] = 1;
441 :
442 0 : aiStart[poGDS->iYDim] = nYOff;
443 0 : aiEdges[poGDS->iYDim] = nYSize;
444 :
445 0 : aiStart[poGDS->iXDim] = nBlockXOff;
446 0 : aiEdges[poGDS->iXDim] = nBlockXSize;
447 0 : break;
448 0 : case 2: // 2Dim: rows/cols
449 0 : aiStart[poGDS->iYDim] = nYOff;
450 0 : aiEdges[poGDS->iYDim] = nYSize;
451 :
452 0 : aiStart[poGDS->iXDim] = nBlockXOff;
453 0 : aiEdges[poGDS->iXDim] = nBlockXSize;
454 0 : break;
455 : }
456 :
457 : /* Ensure that we don't overlap the bottom or right edges */
458 : /* of the dataset in order to use the GDreadtile() API */
459 0 : if (poGDS->bReadTile &&
460 0 : (nBlockXOff + 1) * nBlockXSize <= nRasterXSize &&
461 0 : (nBlockYOff + 1) * nBlockYSize <= nRasterYSize)
462 : {
463 0 : int32 tilecoords[] = {nBlockYOff, nBlockXOff};
464 0 : if (GDreadtile(hGD, poGDS->pszFieldName, tilecoords,
465 0 : pImage) != 0)
466 : {
467 0 : CPLError(CE_Failure, CPLE_AppDefined,
468 : "GDreadtile() failed for block.");
469 0 : eErr = CE_Failure;
470 0 : }
471 : }
472 0 : else if (GDreadfield(hGD, poGDS->pszFieldName, aiStart,
473 0 : nullptr, aiEdges, pImage) < 0)
474 : {
475 0 : CPLError(CE_Failure, CPLE_AppDefined,
476 : "GDreadfield() failed for block.");
477 0 : eErr = CE_Failure;
478 : }
479 0 : GDdetach(hGD);
480 : }
481 0 : break;
482 :
483 0 : case H4ST_EOS_SWATH:
484 : case H4ST_EOS_SWATH_GEOL:
485 : {
486 : const int32 hSW =
487 0 : SWattach(poGDS->hHDF4, poGDS->pszSubdatasetName);
488 0 : switch (poGDS->iRank)
489 : {
490 0 : case 3: // 3Dim: volume
491 0 : aiStart[poGDS->iBandDim] = nBand - 1;
492 0 : aiEdges[poGDS->iBandDim] = 1;
493 :
494 0 : aiStart[poGDS->iYDim] = nYOff;
495 0 : aiEdges[poGDS->iYDim] = nYSize;
496 :
497 0 : aiStart[poGDS->iXDim] = nBlockXOff;
498 0 : aiEdges[poGDS->iXDim] = nBlockXSize;
499 0 : break;
500 0 : case 2: // 2Dim: rows/cols
501 0 : aiStart[poGDS->iYDim] = nYOff;
502 0 : aiEdges[poGDS->iYDim] = nYSize;
503 :
504 0 : aiStart[poGDS->iXDim] = nBlockXOff;
505 0 : aiEdges[poGDS->iXDim] = nBlockXSize;
506 0 : break;
507 : }
508 0 : if (SWreadfield(hSW, poGDS->pszFieldName, aiStart, nullptr,
509 0 : aiEdges, pImage) < 0)
510 : {
511 0 : CPLError(CE_Failure, CPLE_AppDefined,
512 : "SWreadfield() failed for block.");
513 0 : eErr = CE_Failure;
514 : }
515 0 : SWdetach(hSW);
516 : }
517 0 : break;
518 0 : default:
519 0 : break;
520 : }
521 : }
522 0 : break;
523 :
524 0 : default:
525 0 : eErr = CE_Failure;
526 0 : break;
527 : }
528 :
529 62 : return eErr;
530 : }
531 :
532 : /************************************************************************/
533 : /* IWriteBlock() */
534 : /************************************************************************/
535 :
536 57 : CPLErr HDF4ImageRasterBand::IWriteBlock(int nBlockXOff, int nBlockYOff,
537 : void *pImage)
538 : {
539 57 : CPLAssert(nBlockXOff == 0);
540 57 : CPLAssert(nBlockYOff >= 0);
541 57 : CPLAssert(pImage != nullptr);
542 :
543 57 : HDF4ImageDataset *poGDS = cpl::down_cast<HDF4ImageDataset *>(poDS);
544 57 : CPLAssert(poGDS != nullptr);
545 :
546 57 : int32 aiStart[H4_MAX_NC_DIMS] = {};
547 57 : int32 aiEdges[H4_MAX_NC_DIMS] = {};
548 57 : CPLErr eErr = CE_None;
549 :
550 57 : CPLMutexHolderD(&hHDF4Mutex);
551 :
552 : /* -------------------------------------------------------------------- */
553 : /* Work out some block oriented details. */
554 : /* -------------------------------------------------------------------- */
555 57 : const int nYOff = nBlockYOff * nBlockYSize;
556 : const int nYSize =
557 57 : std::min(nYOff + nBlockYSize, poDS->GetRasterYSize()) - nYOff;
558 :
559 : /* -------------------------------------------------------------------- */
560 : /* Process based on rank. */
561 : /* -------------------------------------------------------------------- */
562 57 : switch (poGDS->iRank)
563 : {
564 39 : case 3:
565 : {
566 39 : const int32 l_iSDS = SDselect(poGDS->hSD, poGDS->iDataset);
567 :
568 39 : aiStart[poGDS->iBandDim] = nBand - 1;
569 39 : aiEdges[poGDS->iBandDim] = 1;
570 :
571 39 : aiStart[poGDS->iYDim] = nYOff;
572 39 : aiEdges[poGDS->iYDim] = nYSize;
573 :
574 39 : aiStart[poGDS->iXDim] = nBlockXOff;
575 39 : aiEdges[poGDS->iXDim] = nBlockXSize;
576 :
577 39 : if ((SDwritedata(l_iSDS, aiStart, nullptr, aiEdges,
578 39 : static_cast<VOIDP>(pImage))) < 0)
579 0 : eErr = CE_Failure;
580 :
581 39 : SDendaccess(l_iSDS);
582 : }
583 39 : break;
584 :
585 18 : case 2:
586 : {
587 18 : const int32 l_iSDS = SDselect(poGDS->hSD, nBand - 1);
588 18 : aiStart[poGDS->iYDim] = nYOff;
589 18 : aiEdges[poGDS->iYDim] = nYSize;
590 :
591 18 : aiStart[poGDS->iXDim] = nBlockXOff;
592 18 : aiEdges[poGDS->iXDim] = nBlockXSize;
593 :
594 18 : if ((SDwritedata(l_iSDS, aiStart, nullptr, aiEdges,
595 18 : static_cast<VOIDP>(pImage))) < 0)
596 0 : eErr = CE_Failure;
597 :
598 18 : SDendaccess(l_iSDS);
599 : }
600 18 : break;
601 :
602 0 : default:
603 0 : eErr = CE_Failure;
604 0 : break;
605 : }
606 :
607 114 : return eErr;
608 : }
609 :
610 : /************************************************************************/
611 : /* GetColorTable() */
612 : /************************************************************************/
613 :
614 2 : GDALColorTable *HDF4ImageRasterBand::GetColorTable()
615 : {
616 2 : HDF4ImageDataset *poGDS = cpl::down_cast<HDF4ImageDataset *>(poDS);
617 :
618 2 : return poGDS->poColorTable;
619 : }
620 :
621 : /************************************************************************/
622 : /* GetColorInterpretation() */
623 : /************************************************************************/
624 :
625 18 : GDALColorInterp HDF4ImageRasterBand::GetColorInterpretation()
626 : {
627 18 : HDF4ImageDataset *poGDS = cpl::down_cast<HDF4ImageDataset *>(poDS);
628 :
629 18 : if (poGDS->iDatasetType == HDF4_SDS)
630 : {
631 18 : return GCI_GrayIndex;
632 : }
633 0 : else if (poGDS->iDatasetType == HDF4_GR)
634 : {
635 0 : if (poGDS->poColorTable != nullptr)
636 : {
637 0 : return GCI_PaletteIndex;
638 : }
639 0 : else if (poGDS->nBands != 1)
640 : {
641 0 : if (nBand == 1)
642 0 : return GCI_RedBand;
643 0 : else if (nBand == 2)
644 0 : return GCI_GreenBand;
645 0 : else if (nBand == 3)
646 0 : return GCI_BlueBand;
647 0 : else if (nBand == 4)
648 0 : return GCI_AlphaBand;
649 : else
650 0 : return GCI_Undefined;
651 : }
652 : else
653 : {
654 0 : return GCI_GrayIndex;
655 : }
656 : }
657 :
658 0 : return GCI_GrayIndex;
659 : }
660 :
661 : /************************************************************************/
662 : /* GetNoDataValue() */
663 : /************************************************************************/
664 :
665 162 : double HDF4ImageRasterBand::GetNoDataValue(int *pbSuccess)
666 :
667 : {
668 162 : if (pbSuccess)
669 162 : *pbSuccess = bNoDataSet;
670 :
671 162 : return dfNoDataValue;
672 : }
673 :
674 : /************************************************************************/
675 : /* SetNoDataValue() */
676 : /************************************************************************/
677 :
678 54 : CPLErr HDF4ImageRasterBand::SetNoDataValue(double dfNoData)
679 :
680 : {
681 54 : bNoDataSet = true;
682 54 : dfNoDataValue = dfNoData;
683 :
684 54 : return CE_None;
685 : }
686 :
687 : /************************************************************************/
688 : /* GetUnitType() */
689 : /************************************************************************/
690 :
691 0 : const char *HDF4ImageRasterBand::GetUnitType()
692 :
693 : {
694 0 : if (!osUnitType.empty())
695 0 : return osUnitType;
696 :
697 0 : return GDALRasterBand::GetUnitType();
698 : }
699 :
700 : /************************************************************************/
701 : /* GetOffset() */
702 : /************************************************************************/
703 :
704 0 : double HDF4ImageRasterBand::GetOffset(int *pbSuccess)
705 :
706 : {
707 0 : if (bHaveOffset)
708 : {
709 0 : if (pbSuccess != nullptr)
710 0 : *pbSuccess = TRUE;
711 0 : return dfOffset;
712 : }
713 :
714 0 : return GDALRasterBand::GetOffset(pbSuccess);
715 : }
716 :
717 : /************************************************************************/
718 : /* GetScale() */
719 : /************************************************************************/
720 :
721 0 : double HDF4ImageRasterBand::GetScale(int *pbSuccess)
722 :
723 : {
724 0 : if (bHaveScale)
725 : {
726 0 : if (pbSuccess != nullptr)
727 0 : *pbSuccess = TRUE;
728 0 : return dfScale;
729 : }
730 :
731 0 : return GDALRasterBand::GetScale(pbSuccess);
732 : }
733 :
734 : /************************************************************************/
735 : /* ==================================================================== */
736 : /* HDF4ImageDataset */
737 : /* ==================================================================== */
738 : /************************************************************************/
739 :
740 : /************************************************************************/
741 : /* HDF4ImageDataset() */
742 : /************************************************************************/
743 :
744 473 : HDF4ImageDataset::HDF4ImageDataset()
745 : : pszFilename(nullptr), hHDF4(0), iGR(0), iPal(0), iDataset(0), iRank(0),
746 : iNumType(0), nAttrs(0), iInterlaceMode(0), iPalInterlaceMode(0),
747 : iPalDataType(0), nComps(0), nPalEntries(0), iXDim(0), iYDim(0),
748 : iBandDim(-1), i4Dim(0), nBandCount(0), pszSubdatasetName(nullptr),
749 : pszFieldName(nullptr), poColorTable(nullptr), bHasGeoTransform(false),
750 : iDatasetType(HDF4_UNKNOWN), iSDS(FAIL), nBlockPreferredXSize(-1),
751 473 : nBlockPreferredYSize(-1), bReadTile(false)
752 : {
753 473 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
754 473 : m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
755 473 : memset(aiDimSizes, 0, sizeof(aiDimSizes));
756 473 : memset(aiPaletteData, 0, sizeof(aiPaletteData));
757 473 : memset(szName, 0, sizeof(szName));
758 473 : }
759 :
760 : /************************************************************************/
761 : /* ~HDF4ImageDataset() */
762 : /************************************************************************/
763 :
764 946 : HDF4ImageDataset::~HDF4ImageDataset()
765 : {
766 946 : CPLMutexHolderD(&hHDF4Mutex);
767 :
768 473 : HDF4ImageDataset::FlushCache(true);
769 :
770 473 : CPLFree(pszFilename);
771 473 : if (iSDS != FAIL)
772 52 : SDendaccess(iSDS);
773 473 : if (hSD > 0)
774 471 : SDend(hSD);
775 473 : hSD = 0;
776 473 : if (iGR > 0)
777 2 : GRendaccess(iGR);
778 473 : if (hGR > 0)
779 2 : GRend(hGR);
780 473 : hGR = 0;
781 473 : CPLFree(pszSubdatasetName);
782 473 : CPLFree(pszFieldName);
783 473 : if (papszLocalMetadata)
784 299 : CSLDestroy(papszLocalMetadata);
785 473 : if (poColorTable != nullptr)
786 1 : delete poColorTable;
787 :
788 473 : if (hHDF4 > 0)
789 : {
790 301 : switch (iDatasetType)
791 : {
792 0 : case HDF4_EOS:
793 0 : switch (iSubdatasetType)
794 : {
795 0 : case H4ST_EOS_SWATH:
796 : case H4ST_EOS_SWATH_GEOL:
797 0 : SWclose(hHDF4);
798 0 : break;
799 0 : case H4ST_EOS_GRID:
800 0 : GDclose(hHDF4);
801 0 : break;
802 0 : default:
803 0 : break;
804 : }
805 0 : break;
806 301 : case HDF4_SDS:
807 : case HDF4_GR:
808 301 : hHDF4 = Hclose(hHDF4);
809 301 : break;
810 0 : default:
811 0 : break;
812 : }
813 : }
814 946 : }
815 :
816 : /************************************************************************/
817 : /* GetGeoTransform() */
818 : /************************************************************************/
819 :
820 45 : CPLErr HDF4ImageDataset::GetGeoTransform(GDALGeoTransform >) const
821 : {
822 45 : gt = m_gt;
823 :
824 45 : if (!bHasGeoTransform)
825 0 : return CE_Failure;
826 :
827 45 : return CE_None;
828 : }
829 :
830 : /************************************************************************/
831 : /* SetGeoTransform() */
832 : /************************************************************************/
833 :
834 96 : CPLErr HDF4ImageDataset::SetGeoTransform(const GDALGeoTransform >)
835 : {
836 96 : bHasGeoTransform = true;
837 96 : m_gt = gt;
838 :
839 96 : return CE_None;
840 : }
841 :
842 : /************************************************************************/
843 : /* GetSpatialRef() */
844 : /************************************************************************/
845 :
846 18 : const OGRSpatialReference *HDF4ImageDataset::GetSpatialRef() const
847 :
848 : {
849 18 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
850 : }
851 :
852 : /************************************************************************/
853 : /* SetSpatialRef() */
854 : /************************************************************************/
855 :
856 78 : CPLErr HDF4ImageDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
857 :
858 : {
859 78 : m_oSRS.Clear();
860 78 : if (poSRS)
861 78 : m_oSRS = *poSRS;
862 :
863 78 : return CE_None;
864 : }
865 :
866 : /************************************************************************/
867 : /* GetGCPCount() */
868 : /************************************************************************/
869 :
870 0 : int HDF4ImageDataset::GetGCPCount()
871 :
872 : {
873 0 : return static_cast<int>(m_aoGCPs.size());
874 : }
875 :
876 : /************************************************************************/
877 : /* GetGCPSpatialRef() */
878 : /************************************************************************/
879 :
880 0 : const OGRSpatialReference *HDF4ImageDataset::GetGCPSpatialRef() const
881 :
882 : {
883 0 : return m_oSRS.IsEmpty() || m_aoGCPs.empty() ? nullptr : &m_oGCPSRS;
884 : }
885 :
886 : /************************************************************************/
887 : /* GetGCPs() */
888 : /************************************************************************/
889 :
890 0 : const GDAL_GCP *HDF4ImageDataset::GetGCPs()
891 : {
892 0 : return gdal::GCP::c_ptr(m_aoGCPs);
893 : }
894 :
895 : /************************************************************************/
896 : /* FlushCache() */
897 : /************************************************************************/
898 :
899 491 : CPLErr HDF4ImageDataset::FlushCache(bool bAtClosing)
900 :
901 : {
902 982 : CPLMutexHolderD(&hHDF4Mutex);
903 :
904 491 : CPLErr eErr = GDALDataset::FlushCache(bAtClosing);
905 :
906 491 : if (eAccess == GA_ReadOnly)
907 323 : return eErr;
908 :
909 : // Write out transformation matrix.
910 : const char *pszValue =
911 168 : CPLSPrintf("%f, %f, %f, %f, %f, %f", m_gt[0], m_gt[1], m_gt[2], m_gt[3],
912 168 : m_gt[4], m_gt[5]);
913 336 : if ((SDsetattr(hSD, "TransformationMatrix", DFNT_CHAR8,
914 168 : static_cast<int>(strlen(pszValue)) + 1, pszValue)) < 0)
915 : {
916 0 : eErr = CE_Failure;
917 0 : CPLDebug("HDF4Image",
918 : "Cannot write transformation matrix to output file");
919 : }
920 :
921 : // Write out projection
922 168 : if (!m_oSRS.IsEmpty())
923 : {
924 78 : char *pszWKT = nullptr;
925 78 : m_oSRS.exportToWkt(&pszWKT);
926 78 : if (pszWKT)
927 : {
928 156 : if ((SDsetattr(hSD, "Projection", DFNT_CHAR8,
929 78 : static_cast<int>(strlen(pszWKT)) + 1, pszWKT)) < 0)
930 : {
931 0 : eErr = CE_Failure;
932 0 : CPLDebug("HDF4Image",
933 : "Cannot write projection information to output file");
934 : }
935 78 : CPLFree(pszWKT);
936 : }
937 : }
938 :
939 : // Store all metadata from source dataset as HDF attributes
940 168 : if (GetMetadata())
941 : {
942 51 : char **papszMeta = GetMetadata();
943 :
944 102 : while (*papszMeta)
945 : {
946 51 : char *pszName = nullptr;
947 51 : pszValue = CPLParseNameValue(*papszMeta++, &pszName);
948 87 : if (pszName != nullptr &&
949 36 : (SDsetattr(hSD, pszName, DFNT_CHAR8,
950 36 : static_cast<int>(strlen(pszValue)) + 1, pszValue)) <
951 : 0)
952 : {
953 0 : eErr = CE_Failure;
954 0 : CPLDebug("HDF4Image",
955 : "Cannot write metadata information to output file");
956 : }
957 :
958 51 : CPLFree(pszName);
959 : }
960 : }
961 :
962 : // Write out NoData values
963 378 : for (int iBand = 1; iBand <= nBands; iBand++)
964 : {
965 : HDF4ImageRasterBand *poBand =
966 210 : reinterpret_cast<HDF4ImageRasterBand *>(GetRasterBand(iBand));
967 :
968 210 : if (poBand->bNoDataSet)
969 : {
970 18 : char *pszName = CPLStrdup(CPLSPrintf("NoDataValue%d", iBand));
971 18 : pszValue = CPLSPrintf("%f", poBand->dfNoDataValue);
972 36 : if ((SDsetattr(hSD, pszName, DFNT_CHAR8,
973 18 : static_cast<int>(strlen(pszValue)) + 1, pszValue)) <
974 : 0)
975 : {
976 0 : eErr = CE_Failure;
977 0 : CPLDebug("HDF4Image",
978 : "Cannot write NoData value for band %d "
979 : "to output file",
980 : iBand);
981 : }
982 :
983 18 : CPLFree(pszName);
984 : }
985 : }
986 :
987 : // Write out band descriptions
988 378 : for (int iBand = 1; iBand <= nBands; iBand++)
989 : {
990 : HDF4ImageRasterBand *poBand =
991 210 : reinterpret_cast<HDF4ImageRasterBand *>(GetRasterBand(iBand));
992 :
993 210 : char *pszName = CPLStrdup(CPLSPrintf("BandDesc%d", iBand));
994 210 : pszValue = poBand->GetDescription();
995 210 : if (pszValue != nullptr && !EQUAL(pszValue, ""))
996 : {
997 36 : if (SDsetattr(hSD, pszName, DFNT_CHAR8,
998 18 : static_cast<int>(strlen(pszValue)) + 1, pszValue) < 0)
999 : {
1000 0 : eErr = CE_Failure;
1001 0 : CPLDebug("HDF4Image",
1002 : "Cannot write band's %d description to output file",
1003 : iBand);
1004 : }
1005 : }
1006 :
1007 210 : CPLFree(pszName);
1008 : }
1009 168 : return eErr;
1010 : }
1011 :
1012 : /************************************************************************/
1013 : /* USGSMnemonicToCode() */
1014 : /************************************************************************/
1015 :
1016 0 : long HDF4ImageDataset::USGSMnemonicToCode(const char *pszMnemonic)
1017 : {
1018 0 : if (EQUAL(pszMnemonic, "UTM"))
1019 0 : return 1L;
1020 0 : else if (EQUAL(pszMnemonic, "LAMCC"))
1021 0 : return 4L;
1022 0 : else if (EQUAL(pszMnemonic, "PS"))
1023 0 : return 6L;
1024 0 : else if (EQUAL(pszMnemonic, "PC"))
1025 0 : return 7L;
1026 0 : else if (EQUAL(pszMnemonic, "TM"))
1027 0 : return 9L;
1028 0 : else if (EQUAL(pszMnemonic, "EQRECT"))
1029 0 : return 17L;
1030 0 : else if (EQUAL(pszMnemonic, "OM"))
1031 0 : return 20L;
1032 0 : else if (EQUAL(pszMnemonic, "SOM"))
1033 0 : return 22L;
1034 : else
1035 0 : return 1L; // UTM by default
1036 : }
1037 :
1038 : /************************************************************************/
1039 : /* ToGeoref() */
1040 : /************************************************************************/
1041 :
1042 0 : void HDF4ImageDataset::ToGeoref(double *pdfGeoX, double *pdfGeoY)
1043 : {
1044 0 : OGRSpatialReference *poLatLong = m_oSRS.CloneGeogCS();
1045 0 : OGRCoordinateTransformation *poTransform = nullptr;
1046 0 : if (poLatLong)
1047 : {
1048 0 : poLatLong->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1049 0 : poTransform = OGRCreateCoordinateTransformation(poLatLong, &m_oSRS);
1050 : }
1051 :
1052 0 : if (poTransform != nullptr)
1053 0 : poTransform->Transform(1, pdfGeoX, pdfGeoY, nullptr);
1054 :
1055 0 : if (poTransform != nullptr)
1056 0 : delete poTransform;
1057 :
1058 0 : if (poLatLong != nullptr)
1059 0 : delete poLatLong;
1060 0 : }
1061 :
1062 : /************************************************************************/
1063 : /* ReadCoordinates() */
1064 : /************************************************************************/
1065 :
1066 0 : void HDF4ImageDataset::ReadCoordinates(const char *pszString,
1067 : double *pdfCenterY, double *pdfCenterX)
1068 : {
1069 0 : char **papszStrList = CSLTokenizeString2(pszString, ", ", 0);
1070 0 : *pdfCenterY = CPLAtof(papszStrList[0]); /* lat */
1071 0 : *pdfCenterX = CPLAtof(papszStrList[1]); /* lon */
1072 0 : CSLDestroy(papszStrList);
1073 0 : }
1074 :
1075 : /************************************************************************/
1076 : /* CaptureL1GMTLInfo() */
1077 : /************************************************************************/
1078 :
1079 : /* FILE L71002025_02520010722_M_MTL.L1G
1080 :
1081 : GROUP = L1_METADATA_FILE
1082 : ...
1083 : GROUP = PRODUCT_METADATA
1084 : PRODUCT_TYPE = "L1G"
1085 : PROCESSING_SOFTWARE = "IAS_5.1"
1086 : EPHEMERIS_TYPE = "DEFINITIVE"
1087 : SPACECRAFT_ID = "Landsat7"
1088 : SENSOR_ID = "ETM+"
1089 : ACQUISITION_DATE = 2001-07-22
1090 : WRS_PATH = 002
1091 : STARTING_ROW = 025
1092 : ENDING_ROW = 025
1093 : BAND_COMBINATION = "12345--7-"
1094 : PRODUCT_UL_CORNER_LAT = 51.2704805
1095 : PRODUCT_UL_CORNER_LON = -53.8914311
1096 : PRODUCT_UR_CORNER_LAT = 50.8458100
1097 : PRODUCT_UR_CORNER_LON = -50.9869091
1098 : PRODUCT_LL_CORNER_LAT = 49.6960897
1099 : PRODUCT_LL_CORNER_LON = -54.4047933
1100 : PRODUCT_LR_CORNER_LAT = 49.2841436
1101 : PRODUCT_LR_CORNER_LON = -51.5900428
1102 : PRODUCT_UL_CORNER_MAPX = 298309.894
1103 : PRODUCT_UL_CORNER_MAPY = 5683875.631
1104 : PRODUCT_UR_CORNER_MAPX = 500921.624
1105 : PRODUCT_UR_CORNER_MAPY = 5632678.683
1106 : PRODUCT_LL_CORNER_MAPX = 254477.193
1107 : PRODUCT_LL_CORNER_MAPY = 5510407.880
1108 : PRODUCT_LR_CORNER_MAPX = 457088.923
1109 : PRODUCT_LR_CORNER_MAPY = 5459210.932
1110 : PRODUCT_SAMPLES_REF = 6967
1111 : PRODUCT_LINES_REF = 5965
1112 : BAND1_FILE_NAME = "L71002025_02520010722_B10.L1G"
1113 : BAND2_FILE_NAME = "L71002025_02520010722_B20.L1G"
1114 : BAND3_FILE_NAME = "L71002025_02520010722_B30.L1G"
1115 : BAND4_FILE_NAME = "L71002025_02520010722_B40.L1G"
1116 : BAND5_FILE_NAME = "L71002025_02520010722_B50.L1G"
1117 : BAND7_FILE_NAME = "L72002025_02520010722_B70.L1G"
1118 : METADATA_L1_FILE_NAME = "L71002025_02520010722_MTL.L1G"
1119 : CPF_FILE_NAME = "L7CPF20010701_20010930_06"
1120 : HDF_DIR_FILE_NAME = "L71002025_02520010722_HDF.L1G"
1121 : END_GROUP = PRODUCT_METADATA
1122 : ...
1123 : GROUP = PROJECTION_PARAMETERS
1124 : REFERENCE_DATUM = "NAD83"
1125 : REFERENCE_ELLIPSOID = "GRS80"
1126 : GRID_CELL_SIZE_PAN = 15.000000
1127 : GRID_CELL_SIZE_THM = 60.000000
1128 : GRID_CELL_SIZE_REF = 30.000000
1129 : ORIENTATION = "NOM"
1130 : RESAMPLING_OPTION = "CC"
1131 : MAP_PROJECTION = "UTM"
1132 : END_GROUP = PROJECTION_PARAMETERS
1133 : GROUP = UTM_PARAMETERS
1134 : ZONE_NUMBER = 22
1135 : END_GROUP = UTM_PARAMETERS
1136 : END_GROUP = L1_METADATA_FILE
1137 : END
1138 : */
1139 :
1140 299 : void HDF4ImageDataset::CaptureL1GMTLInfo()
1141 :
1142 : {
1143 : /* -------------------------------------------------------------------- */
1144 : /* Does the physical file look like it matches our expected */
1145 : /* name pattern? */
1146 : /* -------------------------------------------------------------------- */
1147 299 : if (strlen(pszFilename) < 8 ||
1148 299 : !EQUAL(pszFilename + strlen(pszFilename) - 8, "_HDF.L1G"))
1149 299 : return;
1150 :
1151 : /* -------------------------------------------------------------------- */
1152 : /* Construct the name of the corresponding MTL file. We should */
1153 : /* likely be able to extract that from the HDF itself but I'm */
1154 : /* not sure where to find it. */
1155 : /* -------------------------------------------------------------------- */
1156 0 : CPLString osMTLFilename = pszFilename;
1157 0 : osMTLFilename.resize(osMTLFilename.length() - 8);
1158 0 : osMTLFilename += "_MTL.L1G";
1159 :
1160 : /* -------------------------------------------------------------------- */
1161 : /* Ingest the MTL using the NASAKeywordHandler written for the */
1162 : /* PDS driver. */
1163 : /* -------------------------------------------------------------------- */
1164 0 : VSILFILE *fp = VSIFOpenL(osMTLFilename, "r");
1165 0 : if (fp == nullptr)
1166 0 : return;
1167 :
1168 0 : NASAKeywordHandler oMTL;
1169 0 : if (!oMTL.Ingest(fp, 0))
1170 : {
1171 0 : VSIFCloseL(fp);
1172 0 : return;
1173 : }
1174 :
1175 0 : VSIFCloseL(fp);
1176 :
1177 : /* -------------------------------------------------------------------- */
1178 : /* Note: Different variation of MTL files use different group names. */
1179 : /* Check for LPGS_METADATA_FILE and L1_METADATA_FILE. */
1180 : /* -------------------------------------------------------------------- */
1181 0 : CPLString osPrefix;
1182 :
1183 0 : if (oMTL.GetKeyword("LPGS_METADATA_FILE.PRODUCT_METADATA"
1184 : ".PRODUCT_UL_CORNER_LON",
1185 0 : nullptr))
1186 0 : osPrefix = "LPGS_METADATA_FILE.PRODUCT_METADATA.PRODUCT_";
1187 0 : else if (oMTL.GetKeyword("L1_METADATA_FILE.PRODUCT_METADATA"
1188 : ".PRODUCT_UL_CORNER_LON",
1189 0 : nullptr))
1190 0 : osPrefix = "L1_METADATA_FILE.PRODUCT_METADATA.PRODUCT_";
1191 : else
1192 0 : return;
1193 :
1194 : const double dfULX =
1195 0 : CPLAtof(oMTL.GetKeyword((osPrefix + "UL_CORNER_LON").c_str(), "0"));
1196 : const double dfULY =
1197 0 : CPLAtof(oMTL.GetKeyword((osPrefix + "UL_CORNER_LAT").c_str(), "0"));
1198 : const double dfLRX =
1199 0 : CPLAtof(oMTL.GetKeyword((osPrefix + "LR_CORNER_LON").c_str(), "0"));
1200 : const double dfLRY =
1201 0 : CPLAtof(oMTL.GetKeyword((osPrefix + "LR_CORNER_LAT").c_str(), "0"));
1202 : const double dfLLX =
1203 0 : CPLAtof(oMTL.GetKeyword((osPrefix + "LL_CORNER_LON").c_str(), "0"));
1204 : const double dfLLY =
1205 0 : CPLAtof(oMTL.GetKeyword((osPrefix + "LL_CORNER_LAT").c_str(), "0"));
1206 : const double dfURX =
1207 0 : CPLAtof(oMTL.GetKeyword((osPrefix + "UR_CORNER_LON").c_str(), "0"));
1208 : const double dfURY =
1209 0 : CPLAtof(oMTL.GetKeyword((osPrefix + "UR_CORNER_LAT").c_str(), "0"));
1210 :
1211 0 : m_oGCPSRS.importFromWkt(
1212 : "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,"
1213 : "298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],TOWGS84[0,0,0,0,0,0,0],"
1214 : "AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,"
1215 : "AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,"
1216 : "AUTHORITY[\"EPSG\",\"9108\"]],AXIS[\"Lat\",NORTH],AXIS[\"Long\",EAST]"
1217 : ",AUTHORITY[\"EPSG\",\"4326\"]]");
1218 :
1219 0 : m_aoGCPs.emplace_back(nullptr, nullptr, 0.0, 0.0, dfULX, dfULY);
1220 0 : m_aoGCPs.emplace_back(nullptr, nullptr, GetRasterXSize(), 0.0, dfURX,
1221 0 : dfURY);
1222 0 : m_aoGCPs.emplace_back(nullptr, nullptr, 0.0, GetRasterYSize(), dfLLX,
1223 0 : dfLLY);
1224 0 : m_aoGCPs.emplace_back(nullptr, nullptr, GetRasterXSize(), GetRasterYSize(),
1225 0 : dfLRX, dfLRY);
1226 : }
1227 :
1228 : /************************************************************************/
1229 : /* CaptureNRLGeoTransform() */
1230 : /* */
1231 : /* Capture geotransform and coordinate system from NRL (Naval */
1232 : /* Research Laboratory, Stennis Space Center) metadata. */
1233 : /************************************************************************/
1234 :
1235 : /* Example metadata:
1236 : Metadata:
1237 : createTime=Fri Oct 1 18:00:07 2004
1238 : createSoftware=APS v2.8.4
1239 : createPlatform=i686-pc-linux-gnu
1240 : createAgency=Naval Research Laboratory, Stennis Space Center
1241 : sensor=MODIS
1242 : sensorPlatform=TERRA-AM
1243 : sensorAgency=NASA
1244 : sensorType=whiskbroom
1245 : sensorSpectrum=Visible/Thermal
1246 : sensorNumberOfBands=36
1247 : sensorBandUnits=nano meters
1248 : sensorBands=645, 858.5, 469, 555, 1240, 1640, 2130, 412.5, 443, 488, 531, 551,
1249 : 667, 678, 748, 869.5, 905, 936, 940, 3750, 3959, 3959, 4050, 4465.5, 4515.5, 13
1250 : 75, 6715, 7325, 8550, 9730, 11130, 12020, 13335, 13635, 13935, 14235
1251 : sensorBandWidths=50, 35, 20, 20, 20, 24, 50, 15, 10, 10, 10, 10, 10, 10, 10, 1
1252 : 5, 30, 10, 50, 90, 60, 60, 60, 65, 67, 30, 360, 300, 300, 300, 500, 500, 300, 30
1253 : 0, 300, 300
1254 : sensorNominalAltitudeInKM=705
1255 : sensorScanWidthInKM=2330
1256 : sensorResolutionInKM=1
1257 : sensorPlatformType=Polar-orbiting Satellite
1258 : timeStartYear=2004
1259 : timeStartDay=275
1260 : timeStartTime=56400000
1261 : timeStart=Fri Oct 1 15:40:00 2004
1262 : timeDayNight=Day
1263 : timeEndYear=2004
1264 : timeEndDay=275
1265 : timeEndTime=56700000
1266 : timeEnd=Fri Oct 1 15:45:00 2004
1267 : inputMasks=HIGLINT,CLDICE,LAND,ATMFAIL
1268 : inputMasksInt=523
1269 : processedVersion=1.2
1270 : file=MODAM2004275.L3_Mosaic_NOAA_GMX
1271 : fileTitle=NRL Level-3 Mosaic
1272 : fileVersion=3.0
1273 : fileClassification=UNCLASSIFIED
1274 : fileStatus=EXPERIMENTAL
1275 : navType=mapped
1276 : mapProjectionSystem=NRL(USGS)
1277 : mapProjection=Gomex
1278 : mapUpperLeft=31, -99
1279 : mapUpperRight=31, -79
1280 : mapLowerLeft=14.9844128048645, -99
1281 : mapLowerRight=14.9844128048645, -79
1282 : inputFiles=MODAM2004275154000.L3_NOAA_GMX
1283 : ...
1284 : */
1285 :
1286 0 : void HDF4ImageDataset::CaptureNRLGeoTransform()
1287 :
1288 : {
1289 : /* -------------------------------------------------------------------- */
1290 : /* Collect the four corners. */
1291 : /* -------------------------------------------------------------------- */
1292 0 : double adfXY[8] = {};
1293 : static const char *const apszItems[] = {"mapUpperLeft", "mapUpperRight",
1294 : "mapLowerLeft", "mapLowerRight"};
1295 0 : bool bLLPossible = true;
1296 :
1297 0 : for (int iCorner = 0; iCorner < 4; iCorner++)
1298 : {
1299 : const char *pszCornerLoc =
1300 0 : CSLFetchNameValue(papszGlobalMetadata, apszItems[iCorner]);
1301 :
1302 0 : if (pszCornerLoc == nullptr)
1303 0 : return;
1304 :
1305 : char **papszTokens =
1306 0 : CSLTokenizeStringComplex(pszCornerLoc, ",", FALSE, FALSE);
1307 0 : if (CSLCount(papszTokens) != 2)
1308 : {
1309 0 : CSLDestroy(papszTokens);
1310 0 : return;
1311 : }
1312 :
1313 0 : adfXY[iCorner * 2 + 0] = CPLAtof(papszTokens[1]);
1314 0 : adfXY[iCorner * 2 + 1] = CPLAtof(papszTokens[0]);
1315 :
1316 0 : if (adfXY[iCorner * 2 + 0] < -360 || adfXY[iCorner * 2 + 0] > 360 ||
1317 0 : adfXY[iCorner * 2 + 1] < -90 || adfXY[iCorner * 2 + 1] > 90)
1318 0 : bLLPossible = false;
1319 :
1320 0 : CSLDestroy(papszTokens);
1321 : }
1322 :
1323 : /* -------------------------------------------------------------------- */
1324 : /* Does this look like nice clean "northup" lat/long data? */
1325 : /* -------------------------------------------------------------------- */
1326 0 : if (adfXY[0 * 2 + 0] == adfXY[2 * 2 + 0] &&
1327 0 : adfXY[0 * 2 + 1] == adfXY[1 * 2 + 1] && bLLPossible)
1328 : {
1329 0 : bHasGeoTransform = true;
1330 0 : m_gt[0] = adfXY[0 * 2 + 0];
1331 0 : m_gt[1] = (adfXY[1 * 2 + 0] - adfXY[0 * 2 + 0]) / nRasterXSize;
1332 0 : m_gt[2] = 0.0;
1333 0 : m_gt[3] = adfXY[0 * 2 + 1];
1334 0 : m_gt[4] = 0.0;
1335 0 : m_gt[5] = (adfXY[2 * 2 + 1] - adfXY[0 * 2 + 1]) / nRasterYSize;
1336 :
1337 0 : m_oSRS.SetWellKnownGeogCS("WGS84");
1338 : }
1339 :
1340 : /* -------------------------------------------------------------------- */
1341 : /* Can we find the USGS Projection Parameters? */
1342 : /* -------------------------------------------------------------------- */
1343 0 : bool bGotGCTPProjection = false;
1344 0 : int l_iSDSIndex = FAIL;
1345 0 : int l_iSDS = FAIL;
1346 : const char *mapProjection =
1347 0 : CSLFetchNameValue(papszGlobalMetadata, "mapProjection");
1348 :
1349 0 : if (mapProjection)
1350 0 : l_iSDSIndex = SDnametoindex(hSD, mapProjection);
1351 :
1352 0 : if (l_iSDSIndex != FAIL)
1353 0 : l_iSDS = SDselect(hSD, l_iSDSIndex);
1354 :
1355 0 : if (l_iSDS != FAIL)
1356 : {
1357 0 : char l_szName[HDF4_SDS_MAXNAMELEN] = {};
1358 0 : int32 l_iRank = 0;
1359 0 : int32 l_iNumType = 0;
1360 0 : int32 l_nAttrs = 0;
1361 0 : int32 l_aiDimSizes[H4_MAX_VAR_DIMS] = {};
1362 :
1363 0 : double adfGCTP[29] = {};
1364 0 : int32 aiStart[H4_MAX_NC_DIMS] = {};
1365 0 : int32 aiEdges[H4_MAX_NC_DIMS] = {};
1366 :
1367 0 : aiStart[0] = 0;
1368 0 : aiEdges[0] = 29;
1369 :
1370 0 : if (SDgetinfo(l_iSDS, l_szName, &l_iRank, l_aiDimSizes, &l_iNumType,
1371 0 : &l_nAttrs) == 0 &&
1372 0 : l_iNumType == DFNT_FLOAT64 && l_iRank == 1 &&
1373 0 : l_aiDimSizes[0] >= 29 &&
1374 0 : SDreaddata(l_iSDS, aiStart, nullptr, aiEdges, adfGCTP) == 0 &&
1375 0 : m_oSRS.importFromUSGS(static_cast<long>(adfGCTP[1]),
1376 0 : static_cast<long>(adfGCTP[2]), adfGCTP + 4,
1377 0 : static_cast<long>(adfGCTP[3])) == OGRERR_NONE)
1378 : {
1379 0 : CPLDebug("HDF4Image",
1380 : "GCTP Params = %g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,"
1381 : "%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g",
1382 : adfGCTP[0], adfGCTP[1], adfGCTP[2], adfGCTP[3], adfGCTP[4],
1383 : adfGCTP[5], adfGCTP[6], adfGCTP[7], adfGCTP[8], adfGCTP[9],
1384 : adfGCTP[10], adfGCTP[11], adfGCTP[12], adfGCTP[13],
1385 : adfGCTP[14], adfGCTP[15], adfGCTP[16], adfGCTP[17],
1386 : adfGCTP[18], adfGCTP[19], adfGCTP[20], adfGCTP[21],
1387 : adfGCTP[22], adfGCTP[23], adfGCTP[24], adfGCTP[25],
1388 : adfGCTP[26], adfGCTP[27], adfGCTP[28]);
1389 :
1390 0 : bGotGCTPProjection = true;
1391 : }
1392 :
1393 0 : SDendaccess(l_iSDS);
1394 : }
1395 :
1396 : /* -------------------------------------------------------------------- */
1397 : /* If we derived a GCTP based projection, then we need to */
1398 : /* transform the lat/long corners into this projection and use */
1399 : /* them to establish the geotransform. */
1400 : /* -------------------------------------------------------------------- */
1401 0 : if (bLLPossible && bGotGCTPProjection)
1402 : {
1403 0 : OGRSpatialReference oWGS84;
1404 :
1405 0 : oWGS84.SetWellKnownGeogCS("WGS84");
1406 0 : oWGS84.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1407 :
1408 : OGRCoordinateTransformation *poCT =
1409 0 : OGRCreateCoordinateTransformation(&oWGS84, &m_oSRS);
1410 :
1411 0 : double dfULX = adfXY[0 * 2 + 0];
1412 0 : double dfULY = adfXY[0 * 2 + 1];
1413 :
1414 0 : double dfLRX = adfXY[3 * 2 + 0];
1415 0 : double dfLRY = adfXY[3 * 2 + 1];
1416 :
1417 0 : if (poCT->Transform(1, &dfULX, &dfULY) &&
1418 0 : poCT->Transform(1, &dfLRX, &dfLRY))
1419 : {
1420 0 : bHasGeoTransform = true;
1421 0 : m_gt[0] = dfULX;
1422 0 : m_gt[1] = (dfLRX - dfULX) / nRasterXSize;
1423 0 : m_gt[2] = 0.0;
1424 0 : m_gt[3] = dfULY;
1425 0 : m_gt[4] = 0.0;
1426 0 : m_gt[5] = (dfLRY - dfULY) / nRasterYSize;
1427 : }
1428 :
1429 0 : delete poCT;
1430 : }
1431 : }
1432 :
1433 : /************************************************************************/
1434 : /* CaptureCoastwatchGCTPInfo() */
1435 : /************************************************************************/
1436 :
1437 : /* Example Metadata from:
1438 :
1439 : http://coastwatch.noaa.gov/interface/most_recent.php?sensor=MODIS&product=chlorNASA
1440 :
1441 : Definitions at:
1442 : http://coastwatch.noaa.gov/cw_form_hdf.html
1443 :
1444 : Metadata:
1445 : satellite=Aqua
1446 : sensor=MODIS
1447 : origin=USDOC/NOAA/NESDIS CoastWatch
1448 : history=PGE01:4.1.12;PGE02:4.3.1.12;SeaDAS Version ?.?, MSl12 4.0.2,
1449 : Linux 2.4.21-27.0.1.EL cwregister GulfOfMexicoSinusoidal.hdf
1450 : MODSCW.P2005023.1835.swath09.hdf MODSCW.P2005023.1835.GM16.mapped09.hdf
1451 : cwgraphics MODSCW.P2005023.1835.GM16.closest.hdf
1452 : cwmath --template chlor_a --expr
1453 : chlor_a=select(and(l2_flags,514)!=0,nan,chlor_a)
1454 : /data/aps/browse/lvl3/seadas/coastwatch/hdf/MODSCW_P2005023_1835_GM16_closest.hdf
1455 : /data/aps/browse/lvl3/seadas/coastwatch/maskhdf/MODSCW_P2005023_1835_GM16_closest_chlora.hdf
1456 : cwmath --template latitude --expr latitude=latitude
1457 : /data/aps/browse/lvl3/seadas/coastwatch/hdf/MODSCW_P2005023_1835_GM16_closest.hdf
1458 : /data/aps/browse/lvl3/seadas/coastwatch/maskhdf/MODSCW_P2005023_1835_GM16_closest_chlora.hdf
1459 : cwmath --template longitude --expr longitude=longitude
1460 : /data/aps/browse/lvl3/seadas/coastwatch/hdf/MODSCW_P2005023_1835_GM16_closest.hdf
1461 : /data/aps/browse/lvl3/seadas/coastwatch/maskhdf/MODSCW_P2005023_1835_GM16_closest_chlora.hdf
1462 : cwhdf_version=3.2
1463 : pass_type=day
1464 : pass_date=12806
1465 : start_time=66906
1466 : temporal_extent=298
1467 : projection_type=mapped
1468 : projection=Sinusoidal
1469 : gctp_sys=16
1470 : gctp_zone=62
1471 : gctp_parm=6378137, 0, 0, 0, -89000000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
1472 : gctp_datum=12
1473 : et_affine=0, -1008.74836097881, 1008.74836097881, 0, -953126.102425113,
1474 : 3447041.10282512 rows=1540 cols=2000 polygon_latitude=31, 31, 31, 31,
1475 : 31, 27.5095879249529, 24.0191758499058, 20.5287637748587, 17.0383516998116, 17.0383516998116,
1476 : 17.0383516998116, 17.0383516998116, 17.0383516998116, 20.5287637748587, 24.0191758499058,
1477 : 27.5095879249529, 31 polygon_longitude=-99, -93.7108573344442,
1478 : -88.4217146688883, -83.1325720033325, -77.8434293377767, -78.217853417453,
1479 : -78.5303805448579, -78.7884829057512, -78.9979508907244, -83.7397542896832,
1480 : -88.481557688642, -93.2233610876007, -97.9651644865595, -98.1529175079091,
1481 : -98.3842631146439, -98.664391423662, -99 orbit_type=ascending
1482 : raster_type=RasterPixelIsArea
1483 : swath_sync_lines=1
1484 :
1485 : */
1486 :
1487 0 : void HDF4ImageDataset::CaptureCoastwatchGCTPInfo()
1488 :
1489 : {
1490 0 : if (CSLFetchNameValue(papszGlobalMetadata, "gctp_sys") == nullptr ||
1491 0 : CSLFetchNameValue(papszGlobalMetadata, "gctp_zone") == nullptr ||
1492 0 : CSLFetchNameValue(papszGlobalMetadata, "gctp_parm") == nullptr ||
1493 0 : CSLFetchNameValue(papszGlobalMetadata, "gctp_datum") == nullptr ||
1494 0 : CSLFetchNameValue(papszGlobalMetadata, "et_affine") == nullptr)
1495 0 : return;
1496 :
1497 : /* -------------------------------------------------------------------- */
1498 : /* Grab USGS/GCTP Parameters. */
1499 : /* -------------------------------------------------------------------- */
1500 0 : const int nSys = atoi(CSLFetchNameValue(papszGlobalMetadata, "gctp_sys"));
1501 0 : const int nZone = atoi(CSLFetchNameValue(papszGlobalMetadata, "gctp_zone"));
1502 : const int nDatum =
1503 0 : atoi(CSLFetchNameValue(papszGlobalMetadata, "gctp_datum"));
1504 :
1505 0 : char **papszTokens = CSLTokenizeStringComplex(
1506 0 : CSLFetchNameValue(papszGlobalMetadata, "gctp_parm"), ",", FALSE, FALSE);
1507 0 : if (CSLCount(papszTokens) < 15)
1508 : {
1509 0 : CSLDestroy(papszTokens);
1510 0 : return;
1511 : }
1512 :
1513 : double adfParams[15];
1514 0 : for (int iParam = 0; iParam < 15; iParam++)
1515 0 : adfParams[iParam] = CPLAtof(papszTokens[iParam]);
1516 0 : CSLDestroy(papszTokens);
1517 :
1518 : /* -------------------------------------------------------------------- */
1519 : /* Convert into an SRS. */
1520 : /* -------------------------------------------------------------------- */
1521 :
1522 0 : if (m_oSRS.importFromUSGS(nSys, nZone, adfParams, nDatum) != OGRERR_NONE)
1523 0 : return;
1524 :
1525 : /* -------------------------------------------------------------------- */
1526 : /* Capture the affine transform info. */
1527 : /* -------------------------------------------------------------------- */
1528 :
1529 0 : papszTokens = CSLTokenizeStringComplex(
1530 0 : CSLFetchNameValue(papszGlobalMetadata, "et_affine"), ",", FALSE, FALSE);
1531 0 : if (CSLCount(papszTokens) != 6)
1532 : {
1533 0 : CSLDestroy(papszTokens);
1534 0 : return;
1535 : }
1536 :
1537 : // We don't seem to have proper ef_affine docs so I don't
1538 : // know which of these two coefficients goes where.
1539 0 : if (CPLAtof(papszTokens[0]) != 0.0 || CPLAtof(papszTokens[3]) != 0.0)
1540 : {
1541 0 : CSLDestroy(papszTokens);
1542 0 : return;
1543 : }
1544 :
1545 0 : bHasGeoTransform = true;
1546 0 : m_gt[0] = CPLAtof(papszTokens[4]);
1547 0 : m_gt[1] = CPLAtof(papszTokens[2]);
1548 0 : m_gt[2] = 0.0;
1549 0 : m_gt[3] = CPLAtof(papszTokens[5]);
1550 0 : m_gt[4] = 0.0;
1551 0 : m_gt[5] = CPLAtof(papszTokens[1]);
1552 :
1553 : // Middle of pixel adjustment.
1554 0 : m_gt[0] -= m_gt[1] * 0.5;
1555 0 : m_gt[3] -= m_gt[5] * 0.5;
1556 :
1557 0 : CSLDestroy(papszTokens);
1558 : }
1559 :
1560 : /************************************************************************/
1561 : /* GetImageDimensions() */
1562 : /************************************************************************/
1563 :
1564 0 : void HDF4ImageDataset::GetImageDimensions(char *pszDimList)
1565 : {
1566 : char **papszDimList =
1567 0 : CSLTokenizeString2(pszDimList, ",", CSLT_HONOURSTRINGS);
1568 0 : const int nDimCount = CSLCount(papszDimList);
1569 :
1570 : // TODO: check whether nDimCount is > 1 and do something if it isn't.
1571 :
1572 : // Search for the "Band" word in the name of dimension
1573 : // or take the first one as a number of bands
1574 0 : if (iRank == 2)
1575 : {
1576 0 : nBandCount = 1;
1577 : }
1578 : else
1579 : {
1580 0 : for (int i = 0; i < nDimCount; i++)
1581 : {
1582 0 : if (strstr(papszDimList[i], "band"))
1583 : {
1584 0 : iBandDim = i;
1585 0 : nBandCount = aiDimSizes[i];
1586 : // Handle 4D datasets
1587 0 : if (iRank > 3 && i < nDimCount - 1)
1588 : {
1589 : // FIXME: is there a better way to search for
1590 : // the 4th dimension?
1591 0 : i4Dim = i + 1;
1592 0 : nBandCount *= aiDimSizes[i4Dim];
1593 : }
1594 0 : break;
1595 : }
1596 : }
1597 : }
1598 :
1599 : // Search for the starting "X" and "Y" in the names or take
1600 : // the last two dimensions as X and Y sizes.
1601 0 : iXDim = nDimCount - 1;
1602 0 : iYDim = nDimCount - 2;
1603 :
1604 0 : for (int i = 0; i < nDimCount; i++)
1605 : {
1606 0 : if (STARTS_WITH_CI(papszDimList[i], "X") && iBandDim != i)
1607 0 : iXDim = i;
1608 0 : else if (STARTS_WITH_CI(papszDimList[i], "Y") && iBandDim != i)
1609 0 : iYDim = i;
1610 : }
1611 :
1612 : // If didn't get a band dimension yet, but have an extra
1613 : // dimension, use it as the band dimension.
1614 0 : if (iRank > 2 && iBandDim == -1)
1615 : {
1616 0 : if (iXDim != 0 && iYDim != 0)
1617 0 : iBandDim = 0;
1618 0 : else if (iXDim != 1 && iYDim != 1)
1619 0 : iBandDim = 1;
1620 0 : else if (iXDim != 2 && iYDim != 2)
1621 0 : iBandDim = 2;
1622 :
1623 0 : nBandCount = aiDimSizes[iBandDim];
1624 : }
1625 :
1626 0 : CSLDestroy(papszDimList);
1627 0 : }
1628 :
1629 : /************************************************************************/
1630 : /* GetSwatAttrs() */
1631 : /************************************************************************/
1632 :
1633 0 : void HDF4ImageDataset::GetSwatAttrs(int32 hSW)
1634 : {
1635 : /* -------------------------------------------------------------------- */
1636 : /* At the start we will fetch the global HDF attributes. */
1637 : /* -------------------------------------------------------------------- */
1638 0 : int32 hDummy = 0;
1639 :
1640 0 : EHidinfo(hHDF4, &hDummy, &hSD);
1641 0 : ReadGlobalAttributes(hSD);
1642 0 : papszLocalMetadata = CSLDuplicate(papszGlobalMetadata);
1643 :
1644 : /* -------------------------------------------------------------------- */
1645 : /* Fetch the esoteric HDF-EOS attributes then. */
1646 : /* -------------------------------------------------------------------- */
1647 0 : int32 nStrBufSize = 0;
1648 :
1649 0 : if (SWinqattrs(hSW, nullptr, &nStrBufSize) > 0 && nStrBufSize > 0)
1650 : {
1651 0 : char *pszAttrList = static_cast<char *>(CPLMalloc(nStrBufSize + 1));
1652 0 : SWinqattrs(hSW, pszAttrList, &nStrBufSize);
1653 :
1654 : #ifdef DEBUG
1655 0 : CPLDebug("HDF4Image", "List of attributes in swath \"%s\": %s",
1656 : pszFieldName, pszAttrList);
1657 : #endif
1658 :
1659 : char **papszAttributes =
1660 0 : CSLTokenizeString2(pszAttrList, ",", CSLT_HONOURSTRINGS);
1661 0 : const int l_nAttrs = CSLCount(papszAttributes);
1662 0 : for (int i = 0; i < l_nAttrs; i++)
1663 : {
1664 0 : int32 l_iNumType = 0;
1665 0 : int32 nValues = 0;
1666 :
1667 0 : if (SWattrinfo(hSW, papszAttributes[i], &l_iNumType, &nValues) < 0)
1668 0 : continue;
1669 0 : const int nDataTypeSize = GetDataTypeSize(l_iNumType);
1670 0 : if (nDataTypeSize == 0)
1671 0 : continue;
1672 0 : CPLAssert((nValues % nDataTypeSize) == 0);
1673 :
1674 0 : void *pData = CPLMalloc(nValues + 1);
1675 0 : SWreadattr(hSW, papszAttributes[i], pData);
1676 :
1677 0 : if (l_iNumType == DFNT_CHAR8 || l_iNumType == DFNT_UCHAR8)
1678 : {
1679 0 : reinterpret_cast<char *>(pData)[nValues] = '\0';
1680 0 : papszLocalMetadata = CSLAddNameValue(
1681 0 : papszLocalMetadata, papszAttributes[i],
1682 : const_cast<const char *>(reinterpret_cast<char *>(pData)));
1683 : }
1684 : else
1685 : {
1686 0 : char *pszTemp = SPrintArray(GetDataType(l_iNumType), pData,
1687 : nValues / nDataTypeSize, ", ");
1688 0 : papszLocalMetadata = CSLAddNameValue(
1689 0 : papszLocalMetadata, papszAttributes[i], pszTemp);
1690 0 : CPLFree(pszTemp);
1691 : }
1692 :
1693 0 : CPLFree(pData);
1694 : }
1695 :
1696 0 : CSLDestroy(papszAttributes);
1697 0 : CPLFree(pszAttrList);
1698 : }
1699 :
1700 : /* -------------------------------------------------------------------- */
1701 : /* After fetching HDF-EOS specific stuff we will read the generic */
1702 : /* HDF attributes and append them to the list of metadata. */
1703 : /* -------------------------------------------------------------------- */
1704 0 : int32 l_iSDS = 0;
1705 0 : if (SWsdid(hSW, pszFieldName, &l_iSDS) != -1)
1706 : {
1707 0 : int32 l_iRank = 0;
1708 0 : int32 l_iNumType = 0;
1709 0 : int32 l_nAttrs = 0;
1710 0 : char l_szName[HDF4_SDS_MAXNAMELEN] = {};
1711 0 : int32 l_aiDimSizes[H4_MAX_VAR_DIMS] = {};
1712 :
1713 0 : if (SDgetinfo(l_iSDS, l_szName, &l_iRank, l_aiDimSizes, &l_iNumType,
1714 0 : &l_nAttrs) == 0)
1715 : {
1716 0 : for (int32 iAttribute = 0; iAttribute < l_nAttrs; iAttribute++)
1717 : {
1718 0 : char szAttrName[H4_MAX_NC_NAME] = {};
1719 0 : int32 nValues = 0;
1720 0 : SDattrinfo(l_iSDS, iAttribute, szAttrName, &l_iNumType,
1721 : &nValues);
1722 0 : papszLocalMetadata = TranslateHDF4Attributes(
1723 : l_iSDS, iAttribute, szAttrName, l_iNumType, nValues,
1724 : papszLocalMetadata);
1725 : }
1726 : }
1727 : }
1728 :
1729 : /* -------------------------------------------------------------------- */
1730 : /* Finally make the whole list visible. */
1731 : /* -------------------------------------------------------------------- */
1732 0 : SetMetadata(papszLocalMetadata);
1733 0 : }
1734 :
1735 : /************************************************************************/
1736 : /* GetGridAttrs() */
1737 : /************************************************************************/
1738 :
1739 0 : void HDF4ImageDataset::GetGridAttrs(int32 hGD)
1740 : {
1741 : /* -------------------------------------------------------------------- */
1742 : /* At the start we will fetch the global HDF attributes. */
1743 : /* -------------------------------------------------------------------- */
1744 0 : int32 hDummy = 0;
1745 :
1746 0 : EHidinfo(hHDF4, &hDummy, &hSD);
1747 0 : ReadGlobalAttributes(hSD);
1748 0 : papszLocalMetadata = CSLDuplicate(papszGlobalMetadata);
1749 :
1750 : /* -------------------------------------------------------------------- */
1751 : /* Fetch the esoteric HDF-EOS attributes then. */
1752 : /* -------------------------------------------------------------------- */
1753 0 : int32 nStrBufSize = 0;
1754 :
1755 0 : if (GDinqattrs(hGD, nullptr, &nStrBufSize) > 0 && nStrBufSize > 0)
1756 : {
1757 0 : char *pszAttrList = static_cast<char *>(CPLMalloc(nStrBufSize + 1));
1758 0 : GDinqattrs(hGD, pszAttrList, &nStrBufSize);
1759 :
1760 : #ifdef DEBUG
1761 0 : CPLDebug("HDF4Image", "List of attributes in grid %s: %s", pszFieldName,
1762 : pszAttrList);
1763 : #endif
1764 :
1765 : char **papszAttributes =
1766 0 : CSLTokenizeString2(pszAttrList, ",", CSLT_HONOURSTRINGS);
1767 0 : const int l_nAttrs = CSLCount(papszAttributes);
1768 0 : for (int i = 0; i < l_nAttrs; i++)
1769 : {
1770 0 : int32 l_iNumType = 0;
1771 0 : int32 nValues = 0;
1772 :
1773 0 : GDattrinfo(hGD, papszAttributes[i], &l_iNumType, &nValues);
1774 0 : const int nDataTypeSize = GetDataTypeSize(l_iNumType);
1775 0 : if (nDataTypeSize == 0)
1776 0 : continue;
1777 0 : CPLAssert((nValues % nDataTypeSize) == 0);
1778 :
1779 0 : void *pData = CPLMalloc(nValues + 1);
1780 0 : GDreadattr(hGD, papszAttributes[i], pData);
1781 :
1782 0 : if (l_iNumType == DFNT_CHAR8 || l_iNumType == DFNT_UCHAR8)
1783 : {
1784 0 : static_cast<char *>(pData)[nValues] = '\0';
1785 0 : papszLocalMetadata =
1786 0 : CSLAddNameValue(papszLocalMetadata, papszAttributes[i],
1787 : static_cast<const char *>(pData));
1788 : }
1789 : else
1790 : {
1791 0 : char *pszTemp = SPrintArray(GetDataType(l_iNumType), pData,
1792 : nValues / nDataTypeSize, ", ");
1793 0 : papszLocalMetadata = CSLAddNameValue(
1794 0 : papszLocalMetadata, papszAttributes[i], pszTemp);
1795 0 : CPLFree(pszTemp);
1796 : }
1797 :
1798 0 : CPLFree(pData);
1799 : }
1800 :
1801 0 : CSLDestroy(papszAttributes);
1802 0 : CPLFree(pszAttrList);
1803 : }
1804 :
1805 : /* -------------------------------------------------------------------- */
1806 : /* After fetching HDF-EOS specific stuff we will read the generic */
1807 : /* HDF attributes and append them to the list of metadata. */
1808 : /* -------------------------------------------------------------------- */
1809 0 : int32 l_iSDS = 0;
1810 0 : if (GDsdid(hGD, pszFieldName, &l_iSDS) != -1)
1811 : {
1812 0 : int32 l_iRank = 0;
1813 0 : int32 l_iNumType = 0;
1814 0 : int32 l_nAttrs = 0;
1815 0 : int32 nValues = 0;
1816 0 : char l_szName[HDF4_SDS_MAXNAMELEN] = {};
1817 0 : int32 l_aiDimSizes[H4_MAX_VAR_DIMS] = {};
1818 :
1819 0 : if (SDgetinfo(l_iSDS, l_szName, &l_iRank, l_aiDimSizes, &l_iNumType,
1820 0 : &l_nAttrs) == 0)
1821 : {
1822 0 : for (int32 iAttribute = 0; iAttribute < l_nAttrs; iAttribute++)
1823 : {
1824 0 : char szAttrName[H4_MAX_NC_NAME] = {};
1825 0 : SDattrinfo(l_iSDS, iAttribute, szAttrName, &l_iNumType,
1826 : &nValues);
1827 0 : papszLocalMetadata = TranslateHDF4Attributes(
1828 : l_iSDS, iAttribute, szAttrName, l_iNumType, nValues,
1829 : papszLocalMetadata);
1830 : }
1831 : }
1832 : }
1833 :
1834 : /* -------------------------------------------------------------------- */
1835 : /* Finally make the whole list visible. */
1836 : /* -------------------------------------------------------------------- */
1837 0 : SetMetadata(papszLocalMetadata);
1838 0 : }
1839 :
1840 : /************************************************************************/
1841 : /* ProcessModisSDSGeolocation() */
1842 : /* */
1843 : /* Recognise latitude and longitude geolocation arrays in */
1844 : /* simple SDS datasets like: */
1845 : /* */
1846 : /* download.osgeo.org/gdal/data/hdf4/A2006005182000.L2_LAC_SST.x.hdf */
1847 : /* */
1848 : /* As reported in ticket #1895. */
1849 : /************************************************************************/
1850 :
1851 299 : void HDF4ImageDataset::ProcessModisSDSGeolocation(void)
1852 :
1853 : {
1854 : // No point in assigning geolocation to the geolocation SDSes themselves.
1855 299 : if (EQUAL(szName, "longitude") || EQUAL(szName, "latitude"))
1856 299 : return;
1857 :
1858 299 : if (nRasterYSize == 1)
1859 0 : return;
1860 :
1861 : /* -------------------------------------------------------------------- */
1862 : /* Scan for latitude and longitude sections. */
1863 : /* -------------------------------------------------------------------- */
1864 299 : int32 nDatasets = 0;
1865 299 : int32 nAttributes = 0;
1866 :
1867 299 : if (SDfileinfo(hSD, &nDatasets, &nAttributes) != 0)
1868 0 : return;
1869 :
1870 299 : int nLongitudeWidth = 0;
1871 299 : int nLongitudeHeight = 0;
1872 299 : int nLatitudeWidth = 0;
1873 299 : int nLatitudeHeight = 0;
1874 299 : int iXIndex = -1;
1875 299 : int iYIndex = -1;
1876 598 : for (int iDSIndex = 0; iDSIndex < nDatasets; iDSIndex++)
1877 : {
1878 299 : int32 l_iRank = 0;
1879 299 : int32 l_iNumType = 0;
1880 299 : int32 l_nAttrs = 0;
1881 299 : char l_szName[HDF4_SDS_MAXNAMELEN] = {};
1882 299 : int32 l_aiDimSizes[H4_MAX_VAR_DIMS] = {};
1883 :
1884 299 : const int32 l_iSDS = SDselect(hSD, iDSIndex);
1885 :
1886 299 : if (SDgetinfo(l_iSDS, l_szName, &l_iRank, l_aiDimSizes, &l_iNumType,
1887 299 : &l_nAttrs) == 0)
1888 : {
1889 299 : if (EQUAL(l_szName, "latitude"))
1890 : {
1891 0 : iYIndex = iDSIndex;
1892 0 : if (l_iRank == 2)
1893 : {
1894 0 : nLatitudeWidth = l_aiDimSizes[1];
1895 0 : nLatitudeHeight = l_aiDimSizes[0];
1896 : }
1897 : }
1898 :
1899 299 : if (EQUAL(l_szName, "longitude"))
1900 : {
1901 0 : iXIndex = iDSIndex;
1902 0 : if (l_iRank == 2)
1903 : {
1904 0 : nLongitudeWidth = l_aiDimSizes[1];
1905 0 : nLongitudeHeight = l_aiDimSizes[0];
1906 : }
1907 : }
1908 : }
1909 :
1910 299 : SDendaccess(l_iSDS);
1911 : }
1912 :
1913 299 : if (iXIndex == -1 || iYIndex == -1)
1914 299 : return;
1915 :
1916 0 : int nPixelOffset = 0;
1917 0 : int nLineOffset = 0;
1918 0 : int nPixelStep = 1;
1919 0 : int nLineStep = 1;
1920 0 : if (nLongitudeWidth != nLatitudeWidth ||
1921 : nLongitudeHeight != nLatitudeHeight)
1922 : {
1923 0 : CPLDebug("HDF4", "Longitude and latitude subdatasets don't have same "
1924 : "dimensions...");
1925 : }
1926 0 : else if (nLongitudeWidth > 0 && nLongitudeHeight > 0)
1927 : {
1928 0 : nPixelStep =
1929 0 : static_cast<int>(0.5 + 1.0 * nRasterXSize / nLongitudeWidth);
1930 0 : nLineStep =
1931 0 : static_cast<int>(0.5 + 1.0 * nRasterYSize / nLongitudeHeight);
1932 0 : nPixelOffset = (nPixelStep - 1) / 2;
1933 0 : nLineOffset = (nLineStep - 1) / 2;
1934 : }
1935 :
1936 : /* -------------------------------------------------------------------- */
1937 : /* We found geolocation information. Record it as metadata. */
1938 : /* -------------------------------------------------------------------- */
1939 :
1940 0 : SetMetadataItem("SRS", SRS_WKT_WGS84_LAT_LONG, "GEOLOCATION");
1941 :
1942 0 : CPLString osWrk;
1943 0 : osWrk.Printf("HDF4_SDS:UNKNOWN:\"%s\":%d", pszFilename, iXIndex);
1944 0 : SetMetadataItem("X_DATASET", osWrk, "GEOLOCATION");
1945 0 : SetMetadataItem("X_BAND", "1", "GEOLOCATION");
1946 :
1947 0 : osWrk.Printf("HDF4_SDS:UNKNOWN:\"%s\":%d", pszFilename, iYIndex);
1948 0 : SetMetadataItem("Y_DATASET", osWrk, "GEOLOCATION");
1949 0 : SetMetadataItem("Y_BAND", "1", "GEOLOCATION");
1950 :
1951 0 : SetMetadataItem("PIXEL_OFFSET", CPLSPrintf("%d", nPixelOffset),
1952 0 : "GEOLOCATION");
1953 0 : SetMetadataItem("PIXEL_STEP", CPLSPrintf("%d", nPixelStep), "GEOLOCATION");
1954 :
1955 0 : SetMetadataItem("LINE_OFFSET", CPLSPrintf("%d", nLineOffset),
1956 0 : "GEOLOCATION");
1957 0 : SetMetadataItem("LINE_STEP", CPLSPrintf("%d", nLineStep), "GEOLOCATION");
1958 : }
1959 :
1960 : /************************************************************************/
1961 : /* ProcessSwathGeolocation() */
1962 : /* */
1963 : /* Handle the swath geolocation data for a swath. Attach */
1964 : /* geolocation metadata corresponding to it (if there is no */
1965 : /* lattice), and also attach it as GCPs. This is only invoked */
1966 : /* for EOS_SWATH, not EOS_SWATH_GEOL datasets. */
1967 : /************************************************************************/
1968 :
1969 0 : int HDF4ImageDataset::ProcessSwathGeolocation(int32 hSW, char **papszDimList)
1970 : {
1971 :
1972 : /* -------------------------------------------------------------------- */
1973 : /* Determine a product name. */
1974 : /* -------------------------------------------------------------------- */
1975 0 : const char *pszProduct = CSLFetchNameValue(papszLocalMetadata, "SHORTNAME");
1976 :
1977 0 : HDF4EOSProduct eProduct = PROD_UNKNOWN;
1978 0 : if (pszProduct)
1979 : {
1980 0 : if (STARTS_WITH_CI(pszProduct, "ASTL1A"))
1981 0 : eProduct = PROD_ASTER_L1A;
1982 0 : else if (STARTS_WITH_CI(pszProduct, "ASTL1B"))
1983 0 : eProduct = PROD_ASTER_L1B;
1984 0 : else if (STARTS_WITH_CI(pszProduct, "AST_04") ||
1985 0 : STARTS_WITH_CI(pszProduct, "AST_05") ||
1986 0 : STARTS_WITH_CI(pszProduct, "AST_06") ||
1987 0 : STARTS_WITH_CI(pszProduct, "AST_07") ||
1988 0 : STARTS_WITH_CI(pszProduct, "AST_08") ||
1989 0 : STARTS_WITH_CI(pszProduct, "AST_09") ||
1990 0 : STARTS_WITH_CI(pszProduct, "AST13") ||
1991 0 : STARTS_WITH_CI(pszProduct, "AST3"))
1992 0 : eProduct = PROD_ASTER_L2;
1993 0 : else if (STARTS_WITH_CI(pszProduct, "AST14"))
1994 0 : eProduct = PROD_ASTER_L3;
1995 0 : else if (STARTS_WITH_CI(pszProduct, "MOD02") ||
1996 0 : STARTS_WITH_CI(pszProduct, "MYD02"))
1997 0 : eProduct = PROD_MODIS_L1B;
1998 0 : else if (STARTS_WITH_CI(pszProduct, "MOD07_L2"))
1999 0 : eProduct = PROD_MODIS_L2;
2000 : }
2001 :
2002 : /* -------------------------------------------------------------------- */
2003 : /* Read names of geolocation fields and corresponding */
2004 : /* geolocation maps. */
2005 : /* -------------------------------------------------------------------- */
2006 0 : int32 nStrBufSize = 0;
2007 0 : const int32 nDataFields = SWnentries(hSW, HDFE_NENTGFLD, &nStrBufSize);
2008 0 : if (nDataFields < 0 || nDataFields > 1024 * 1024)
2009 0 : return FALSE;
2010 0 : char *pszGeoList = static_cast<char *>(CPLMalloc(nStrBufSize + 1));
2011 : int32 *paiRank =
2012 0 : static_cast<int32 *>(CPLMalloc(nDataFields * sizeof(int32)));
2013 : int32 *paiNumType =
2014 0 : static_cast<int32 *>(CPLMalloc(nDataFields * sizeof(int32)));
2015 :
2016 0 : if (nDataFields != SWinqgeofields(hSW, pszGeoList, paiRank, paiNumType))
2017 : {
2018 0 : CPLDebug("HDF4Image",
2019 : "Can't get the list of geolocation fields in swath \"%s\"",
2020 : pszSubdatasetName);
2021 : }
2022 :
2023 : #ifdef DEBUG
2024 : else
2025 : {
2026 0 : CPLDebug("HDF4Image",
2027 : "Number of geolocation fields in swath \"%s\": %ld",
2028 : pszSubdatasetName, static_cast<long>(nDataFields));
2029 0 : CPLDebug("HDF4Image", "List of geolocation fields in swath \"%s\": %s",
2030 : pszSubdatasetName, pszGeoList);
2031 0 : char *pszTmp = SPrintArray(GDT_UInt32, paiRank, nDataFields, ",");
2032 0 : CPLDebug("HDF4Image", "Geolocation fields ranks: %s", pszTmp);
2033 0 : CPLFree(pszTmp);
2034 : }
2035 : #endif
2036 :
2037 0 : CPLFree(paiRank);
2038 0 : CPLFree(paiNumType);
2039 :
2040 : /* -------------------------------------------------------------------- */
2041 : /* Read geolocation data. */
2042 : /* -------------------------------------------------------------------- */
2043 0 : char szXGeo[N_BUF_SIZE] = "";
2044 0 : char szYGeo[N_BUF_SIZE] = "";
2045 0 : char szPixel[N_BUF_SIZE] = "";
2046 0 : char szLine[N_BUF_SIZE] = "";
2047 0 : int32 *paiOffset = nullptr;
2048 0 : int32 *paiIncrement = nullptr;
2049 :
2050 0 : int32 nDimMaps = SWnentries(hSW, HDFE_NENTMAP, &nStrBufSize);
2051 0 : if (nDimMaps <= 0)
2052 : {
2053 :
2054 : #ifdef DEBUG
2055 0 : CPLDebug("HDF4Image", "No geolocation maps in swath \"%s\"",
2056 : pszSubdatasetName);
2057 0 : CPLDebug("HDF4Image",
2058 : "Suppose one-to-one mapping. X field is \"%s\", "
2059 : "Y field is \"%s\"",
2060 0 : papszDimList[iXDim], papszDimList[iYDim]);
2061 : #endif
2062 :
2063 0 : snprintf(szPixel, sizeof(szPixel), "%s", papszDimList[iXDim]);
2064 :
2065 0 : snprintf(szLine, sizeof(szLine), "%s", papszDimList[iYDim]);
2066 :
2067 0 : snprintf(szXGeo, sizeof(szXGeo), "%s", papszDimList[iXDim]);
2068 :
2069 0 : snprintf(szYGeo, sizeof(szYGeo), "%s", papszDimList[iYDim]);
2070 :
2071 0 : paiOffset = static_cast<int32 *>(CPLCalloc(2, sizeof(int32)));
2072 0 : paiIncrement = static_cast<int32 *>(CPLCalloc(2, sizeof(int32)));
2073 0 : paiOffset[0] = 0;
2074 0 : paiOffset[1] = 0;
2075 0 : paiIncrement[0] = 1;
2076 0 : paiIncrement[1] = 1;
2077 : }
2078 : else
2079 : {
2080 0 : char *pszDimMaps = static_cast<char *>(CPLMalloc(nStrBufSize + 1));
2081 0 : paiOffset = static_cast<int32 *>(CPLCalloc(nDimMaps, sizeof(int32)));
2082 0 : paiIncrement = static_cast<int32 *>(CPLCalloc(nDimMaps, sizeof(int32)));
2083 :
2084 0 : *pszDimMaps = '\0';
2085 0 : if (nDimMaps != SWinqmaps(hSW, pszDimMaps, paiOffset, paiIncrement))
2086 : {
2087 0 : CPLDebug("HDF4Image",
2088 : "Can't get the list of geolocation maps in swath \"%s\"",
2089 : pszSubdatasetName);
2090 : }
2091 :
2092 : #ifdef DEBUG
2093 : else
2094 : {
2095 :
2096 0 : CPLDebug("HDF4Image",
2097 : "List of geolocation maps in swath \"%s\": %s",
2098 : pszSubdatasetName, pszDimMaps);
2099 0 : char *pszTmp = SPrintArray(GDT_Int32, paiOffset, nDimMaps, ",");
2100 0 : CPLDebug("HDF4Image", "Geolocation map offsets: %s", pszTmp);
2101 0 : CPLFree(pszTmp);
2102 0 : pszTmp = SPrintArray(GDT_Int32, paiIncrement, nDimMaps, ",");
2103 0 : CPLDebug("HDF4Image", "Geolocation map increments: %s", pszTmp);
2104 0 : CPLFree(pszTmp);
2105 : }
2106 : #endif
2107 :
2108 : char **papszDimMap =
2109 0 : CSLTokenizeString2(pszDimMaps, ",", CSLT_HONOURSTRINGS);
2110 0 : const int nDimMapCount = CSLCount(papszDimMap);
2111 :
2112 0 : for (int i = 0; i < nDimMapCount; i++)
2113 : {
2114 0 : if (strstr(papszDimMap[i], papszDimList[iXDim]))
2115 : {
2116 0 : snprintf(szPixel, sizeof(szPixel), "%s", papszDimList[iXDim]);
2117 :
2118 0 : snprintf(szXGeo, sizeof(szXGeo), "%s", papszDimMap[i]);
2119 :
2120 0 : char *pszTemp = strchr(szXGeo, '/');
2121 0 : if (pszTemp)
2122 0 : *pszTemp = '\0';
2123 : }
2124 0 : else if (strstr(papszDimMap[i], papszDimList[iYDim]))
2125 : {
2126 0 : snprintf(szLine, sizeof(szLine), "%s", papszDimList[iYDim]);
2127 :
2128 0 : snprintf(szYGeo, sizeof(szYGeo), "%s", papszDimMap[i]);
2129 :
2130 0 : char *pszTemp = strchr(szYGeo, '/');
2131 0 : if (pszTemp)
2132 0 : *pszTemp = '\0';
2133 : }
2134 : }
2135 :
2136 0 : CSLDestroy(papszDimMap);
2137 0 : CPLFree(pszDimMaps);
2138 : }
2139 :
2140 0 : if (*szXGeo == 0 || *szYGeo == 0)
2141 : {
2142 0 : CPLFree(paiOffset);
2143 0 : CPLFree(paiIncrement);
2144 0 : CPLFree(pszGeoList);
2145 0 : return FALSE;
2146 : }
2147 :
2148 : /* -------------------------------------------------------------------- */
2149 : /* Read geolocation fields. */
2150 : /* -------------------------------------------------------------------- */
2151 0 : char szGeoDimList[N_BUF_SIZE] = "";
2152 : char **papszGeolocations =
2153 0 : CSLTokenizeString2(pszGeoList, ",", CSLT_HONOURSTRINGS);
2154 0 : const int nGeolocationsCount = CSLCount(papszGeolocations);
2155 0 : int32 l_aiDimSizes[H4_MAX_VAR_DIMS] = {};
2156 :
2157 0 : int32 iWrkNumType = 0;
2158 0 : void *pLat = nullptr;
2159 0 : void *pLong = nullptr;
2160 :
2161 0 : int32 l_iRank = 0;
2162 0 : int32 nLatCount = 0;
2163 0 : int32 nLongCount = 0;
2164 0 : int32 nXPoints = 0;
2165 0 : int32 nYPoints = 0;
2166 0 : int iDataSize = 0;
2167 :
2168 0 : int iPixelDim = -1;
2169 0 : int iLineDim = -1;
2170 0 : int iLongDim = -1;
2171 0 : int iLatDim = -1;
2172 0 : double dfLatNoData = 0;
2173 0 : bool bLatNoDataSet = false;
2174 0 : double dfLongNoData = 0;
2175 0 : bool bLongNoDataSet = false;
2176 0 : for (int i = 0; i < nGeolocationsCount; i++)
2177 : {
2178 : // Skip "SceneLineNumber" table if present in the list of geolocation
2179 : // fields. It is not needed to fetch geocoding data.
2180 0 : if (EQUAL(papszGeolocations[i], "SceneLineNumber"))
2181 0 : continue;
2182 :
2183 0 : if (SWfieldinfo(hSW, papszGeolocations[i], &l_iRank, l_aiDimSizes,
2184 0 : &iWrkNumType, szGeoDimList) < 0)
2185 : {
2186 :
2187 0 : CPLDebug("HDF4Image",
2188 : "Can't read attributes of geolocation field \"%s\"",
2189 0 : papszGeolocations[i]);
2190 0 : CSLDestroy(papszGeolocations);
2191 0 : CPLFree(paiOffset);
2192 0 : CPLFree(paiIncrement);
2193 0 : CPLFree(pszGeoList);
2194 0 : return FALSE;
2195 : }
2196 :
2197 0 : CPLDebug("HDF4Image",
2198 : "List of dimensions in geolocation field \"%s\": %s",
2199 0 : papszGeolocations[i], szGeoDimList);
2200 :
2201 : char **papszGeoDimList =
2202 0 : CSLTokenizeString2(szGeoDimList, ",", CSLT_HONOURSTRINGS);
2203 :
2204 0 : const int iXGeo = CSLFindString(papszGeoDimList, szXGeo);
2205 0 : const int iYGeo = CSLFindString(papszGeoDimList, szYGeo);
2206 0 : if (CSLCount(papszGeoDimList) > H4_MAX_VAR_DIMS || iXGeo < 0 ||
2207 : iYGeo < 0)
2208 : {
2209 0 : CSLDestroy(papszGeoDimList);
2210 0 : CSLDestroy(papszGeolocations);
2211 0 : CPLFree(paiOffset);
2212 0 : CPLFree(paiIncrement);
2213 0 : CPLFree(pszGeoList);
2214 0 : return FALSE;
2215 : }
2216 :
2217 0 : nXPoints = l_aiDimSizes[iXGeo];
2218 0 : nYPoints = l_aiDimSizes[iYGeo];
2219 :
2220 0 : if (EQUAL(szPixel, papszDimList[iXDim]))
2221 : {
2222 0 : iPixelDim = 1;
2223 0 : iLineDim = 0;
2224 : }
2225 : else
2226 : {
2227 0 : iPixelDim = 0;
2228 0 : iLineDim = 1;
2229 : }
2230 :
2231 0 : iDataSize = GetDataTypeSize(iWrkNumType);
2232 0 : if (strstr(papszGeolocations[i], "Latitude"))
2233 : {
2234 0 : iLatDim = i;
2235 0 : nLatCount = nXPoints * nYPoints;
2236 0 : pLat = CPLMalloc(nLatCount * iDataSize);
2237 0 : if (SWreadfield(hSW, papszGeolocations[i], nullptr, nullptr,
2238 0 : nullptr, static_cast<VOIDP>(pLat)) < 0)
2239 : {
2240 0 : CPLDebug("HDF4Image", "Can't read geolocation field %s",
2241 0 : papszGeolocations[i]);
2242 0 : CPLFree(pLat);
2243 0 : pLat = nullptr;
2244 : }
2245 :
2246 0 : void *pNoDataValue = CPLMalloc(iDataSize);
2247 0 : if (SWgetfillvalue(hSW, papszGeolocations[i], pNoDataValue) != -1)
2248 : {
2249 0 : dfLatNoData = AnyTypeToDouble(iWrkNumType, pNoDataValue);
2250 0 : bLatNoDataSet = true;
2251 : }
2252 0 : CPLFree(pNoDataValue);
2253 : }
2254 0 : else if (strstr(papszGeolocations[i], "Longitude"))
2255 : {
2256 0 : iLongDim = i;
2257 0 : nLongCount = nXPoints * nYPoints;
2258 0 : pLong = CPLMalloc(nLongCount * iDataSize);
2259 0 : if (SWreadfield(hSW, papszGeolocations[i], nullptr, nullptr,
2260 0 : nullptr, static_cast<VOIDP>(pLong)) < 0)
2261 : {
2262 0 : CPLDebug("HDF4Image", "Can't read geolocation field %s",
2263 0 : papszGeolocations[i]);
2264 0 : CPLFree(pLong);
2265 0 : pLong = nullptr;
2266 : }
2267 :
2268 0 : void *pNoDataValue = CPLMalloc(iDataSize);
2269 0 : if (SWgetfillvalue(hSW, papszGeolocations[i], pNoDataValue) != -1)
2270 : {
2271 0 : dfLongNoData = AnyTypeToDouble(iWrkNumType, pNoDataValue);
2272 0 : bLongNoDataSet = true;
2273 : }
2274 0 : CPLFree(pNoDataValue);
2275 : }
2276 :
2277 0 : CSLDestroy(papszGeoDimList);
2278 : }
2279 :
2280 : /* -------------------------------------------------------------------- */
2281 : /* Do we have a lattice table? */
2282 : /* -------------------------------------------------------------------- */
2283 0 : void *pLatticeX = nullptr;
2284 0 : void *pLatticeY = nullptr;
2285 0 : int32 iLatticeType = 0;
2286 0 : int32 iLatticeDataSize = 0;
2287 0 : char pszLatticePoint[] = "LatticePoint";
2288 0 : if (SWfieldinfo(hSW, pszLatticePoint, &l_iRank, l_aiDimSizes, &iLatticeType,
2289 0 : szGeoDimList) == 0 &&
2290 0 : l_iRank == 3 && nXPoints == l_aiDimSizes[1] &&
2291 0 : nYPoints == l_aiDimSizes[0] && l_aiDimSizes[2] == 2)
2292 : {
2293 0 : iLatticeDataSize = GetDataTypeSize(iLatticeType);
2294 :
2295 0 : int32 iStart[H4_MAX_NC_DIMS] = {};
2296 0 : int32 iEdges[H4_MAX_NC_DIMS] = {};
2297 0 : iStart[1] = 0;
2298 0 : iEdges[1] = nXPoints;
2299 :
2300 0 : iStart[0] = 0;
2301 0 : iEdges[0] = nYPoints;
2302 :
2303 0 : iStart[2] = 0;
2304 0 : iEdges[2] = 1;
2305 :
2306 0 : pLatticeX = CPLMalloc(nLatCount * iLatticeDataSize);
2307 0 : if (SWreadfield(hSW, pszLatticePoint, iStart, nullptr, iEdges,
2308 0 : static_cast<VOIDP>(pLatticeX)) < 0)
2309 : {
2310 0 : CPLDebug("HDF4Image", "Can't read lattice field");
2311 0 : CPLFree(pLatticeX);
2312 0 : pLatticeX = nullptr;
2313 : }
2314 :
2315 0 : iStart[2] = 1;
2316 0 : iEdges[2] = 1;
2317 :
2318 0 : pLatticeY = CPLMalloc(nLatCount * iLatticeDataSize);
2319 0 : if (SWreadfield(hSW, pszLatticePoint, iStart, nullptr, iEdges,
2320 0 : static_cast<VOIDP>(pLatticeY)) < 0)
2321 : {
2322 0 : CPLDebug("HDF4Image", "Can't read lattice field");
2323 0 : CPLFree(pLatticeY);
2324 0 : pLatticeY = nullptr;
2325 : }
2326 : }
2327 :
2328 : /* -------------------------------------------------------------------- */
2329 : /* Determine whether to use no, partial or full GCPs. */
2330 : /* -------------------------------------------------------------------- */
2331 0 : const char *pszGEOL_AS_GCPS = CPLGetConfigOption("GEOL_AS_GCPS", "PARTIAL");
2332 0 : int iGCPStepX = 0;
2333 0 : int iGCPStepY = 0;
2334 :
2335 0 : if (EQUAL(pszGEOL_AS_GCPS, "NONE"))
2336 : {
2337 : // Leave as is: iGCPStepX = iGCPStepY = 0;
2338 : }
2339 0 : else if (EQUAL(pszGEOL_AS_GCPS, "FULL"))
2340 : {
2341 0 : iGCPStepX = 1;
2342 0 : iGCPStepY = 1;
2343 : }
2344 : else
2345 : {
2346 : // Aim for 10x10 grid or so.
2347 0 : iGCPStepX = std::max(static_cast<int32>(1), ((nXPoints - 1) / 11));
2348 0 : iGCPStepY = std::max(static_cast<int32>(1), ((nYPoints - 1) / 11));
2349 : }
2350 :
2351 : /* -------------------------------------------------------------------- */
2352 : /* Fetch projection information for various datasets. */
2353 : /* -------------------------------------------------------------------- */
2354 0 : if (nLatCount && nLongCount && nLatCount == nLongCount && pLat && pLong)
2355 : {
2356 : // ASTER Level 1A
2357 0 : if (eProduct == PROD_ASTER_L1A)
2358 : {
2359 0 : m_oGCPSRS.importFromWkt(
2360 : "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\","
2361 : "6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],"
2362 : "TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6326\"]],"
2363 : "PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],"
2364 : "UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\","
2365 : "\"9108\"]],AXIS[\"Lat\",NORTH],AXIS[\"Long\",EAST],"
2366 : "AUTHORITY[\"EPSG\",\"4326\"]]");
2367 : }
2368 :
2369 : // ASTER Level 1B, Level 2
2370 0 : else if (eProduct == PROD_ASTER_L1B || eProduct == PROD_ASTER_L2)
2371 : {
2372 : // Construct the metadata keys.
2373 : // A band number is taken from the field name.
2374 0 : const char *pszBand = strpbrk(pszFieldName, "0123456789");
2375 :
2376 0 : if (!pszBand)
2377 0 : pszBand = "";
2378 :
2379 0 : char *pszProjLine = CPLStrdup(CPLSPrintf("MPMETHOD%s", pszBand));
2380 : char *pszParamsLine =
2381 0 : CPLStrdup(CPLSPrintf("PROJECTIONPARAMETERS%s", pszBand));
2382 0 : char *pszZoneLine = CPLStrdup(CPLSPrintf("UTMZONECODE%s", pszBand));
2383 : char *pszEllipsoidLine =
2384 0 : CPLStrdup(CPLSPrintf("ELLIPSOIDANDDATUM%s", pszBand));
2385 :
2386 : // Fetch projection related values from the
2387 : // metadata.
2388 : const char *pszProj =
2389 0 : CSLFetchNameValue(papszLocalMetadata, pszProjLine);
2390 : const char *pszParams =
2391 0 : CSLFetchNameValue(papszLocalMetadata, pszParamsLine);
2392 : const char *pszZone =
2393 0 : CSLFetchNameValue(papszLocalMetadata, pszZoneLine);
2394 :
2395 : #ifdef DEBUG
2396 : const char *pszEllipsoid =
2397 0 : CSLFetchNameValue(papszLocalMetadata, pszEllipsoidLine);
2398 :
2399 0 : CPLDebug("HDF4Image",
2400 : "Projection %s=%s, parameters %s=%s, "
2401 : "zone %s=%s",
2402 : pszProjLine, pszProj, pszParamsLine, pszParams,
2403 : pszZoneLine, pszZone);
2404 0 : CPLDebug("HDF4Image", "Ellipsoid %s=%s", pszEllipsoidLine,
2405 : pszEllipsoid);
2406 : #endif
2407 :
2408 : // Transform all mnemonic codes in the values.
2409 : // Projection is UTM by default
2410 0 : const long iProjSys = pszProj ? USGSMnemonicToCode(pszProj) : 1L;
2411 0 : const long iZone = (pszZone && iProjSys == 1L) ? atoi(pszZone) : 0L;
2412 : #if 0 // Not needed without the WGS84 check.
2413 : char **papszEllipsoid = (pszEllipsoid) ?
2414 : CSLTokenizeString2( pszEllipsoid, ",",
2415 : CSLT_HONOURSTRINGS ) : NULL;
2416 : #endif
2417 :
2418 0 : const long iEllipsoid = 8L; // WGS84 by default
2419 : #if 0 // This block is redundant.
2420 : if( papszEllipsoid && CSLCount(papszEllipsoid) > 0 )
2421 : {
2422 : if (EQUAL( papszEllipsoid[0], "WGS84"))
2423 : iEllipsoid = 8L;
2424 : }
2425 : #endif
2426 : char **papszParams =
2427 : pszParams
2428 0 : ? CSLTokenizeString2(pszParams, ",", CSLT_HONOURSTRINGS)
2429 0 : : nullptr;
2430 0 : std::vector<double> adfProjParams(15);
2431 0 : if (papszParams)
2432 : {
2433 0 : for (int i = 0; i < 15 && papszParams[i] != nullptr; i++)
2434 0 : adfProjParams[i] = CPLAtof(papszParams[i]);
2435 : }
2436 :
2437 : // Create projection definition
2438 0 : m_oSRS.importFromUSGS(iProjSys, iZone, adfProjParams.data(),
2439 : iEllipsoid);
2440 0 : m_oSRS.SetLinearUnits(SRS_UL_METER, 1.0);
2441 :
2442 0 : CSLDestroy(papszParams);
2443 0 : CPLFree(pszEllipsoidLine);
2444 0 : CPLFree(pszZoneLine);
2445 0 : CPLFree(pszParamsLine);
2446 0 : CPLFree(pszProjLine);
2447 : }
2448 :
2449 : // ASTER Level 3 (DEM)
2450 0 : else if (eProduct == PROD_ASTER_L3)
2451 : {
2452 0 : double dfCenterX = 0.0;
2453 0 : double dfCenterY = 0.0;
2454 :
2455 0 : ReadCoordinates(
2456 0 : CSLFetchNameValue(papszGlobalMetadata, "SCENECENTER"),
2457 : &dfCenterY, &dfCenterX);
2458 :
2459 : // Calculate UTM zone from scene center coordinates
2460 0 : const int iZone = 30 + static_cast<int>((dfCenterX + 6.0) / 6.0);
2461 :
2462 : // Create projection definition
2463 0 : if (dfCenterY > 0)
2464 0 : m_oSRS.SetUTM(iZone, TRUE);
2465 : else
2466 0 : m_oSRS.SetUTM(-iZone, FALSE);
2467 0 : m_oSRS.SetWellKnownGeogCS("WGS84");
2468 0 : m_oSRS.SetLinearUnits(SRS_UL_METER, 1.0);
2469 : }
2470 :
2471 : // MODIS L1B
2472 0 : else if (eProduct == PROD_MODIS_L1B || eProduct == PROD_MODIS_L2)
2473 : {
2474 0 : m_oGCPSRS.importFromWkt(SRS_WKT_WGS84_LAT_LONG);
2475 : }
2476 :
2477 : /* --------------------------------------------------------------------
2478 : */
2479 : /* Fill the GCPs list. */
2480 : /* --------------------------------------------------------------------
2481 : */
2482 0 : if (iGCPStepX > 0)
2483 : {
2484 0 : for (int i = 0; i < nYPoints; i += iGCPStepY)
2485 : {
2486 0 : for (int j = 0; j < nXPoints; j += iGCPStepX)
2487 : {
2488 0 : const int iGeoOff = i * nXPoints + j;
2489 :
2490 0 : double dfGCPX = AnyTypeToDouble(
2491 : iWrkNumType, reinterpret_cast<void *>(
2492 : reinterpret_cast<char *>(pLong) +
2493 0 : iGeoOff * iDataSize));
2494 0 : double dfGCPY = AnyTypeToDouble(
2495 : iWrkNumType, reinterpret_cast<void *>(
2496 : reinterpret_cast<char *>(pLat) +
2497 0 : iGeoOff * iDataSize));
2498 :
2499 0 : if (bLongNoDataSet && dfGCPX == dfLongNoData)
2500 0 : continue;
2501 0 : if (bLatNoDataSet && dfGCPY == dfLatNoData)
2502 0 : continue;
2503 :
2504 : // GCPs in Level 1A/1B dataset are in geocentric
2505 : // coordinates. Convert them in geodetic (we
2506 : // will convert latitudes only, longitudes
2507 : // do not need to be converted, because
2508 : // they are the same).
2509 : // This calculation valid for WGS84 datum only.
2510 0 : if (eProduct == PROD_ASTER_L1A ||
2511 : eProduct == PROD_ASTER_L1B)
2512 : {
2513 0 : dfGCPY = atan(tan(dfGCPY * PI / 180) / 0.99330562) *
2514 0 : 180 / PI;
2515 : }
2516 :
2517 0 : ToGeoref(&dfGCPX, &dfGCPY);
2518 :
2519 0 : double dfGCPPixel = 0.0;
2520 0 : double dfGCPLine = 0.0;
2521 0 : if (pLatticeX && pLatticeY)
2522 : {
2523 0 : dfGCPPixel =
2524 0 : AnyTypeToDouble(
2525 : iLatticeType,
2526 : reinterpret_cast<void *>(
2527 : reinterpret_cast<char *>(pLatticeX) +
2528 0 : iGeoOff * iLatticeDataSize)) +
2529 : 0.5;
2530 0 : dfGCPLine =
2531 0 : AnyTypeToDouble(
2532 : iLatticeType,
2533 : reinterpret_cast<void *>(
2534 : reinterpret_cast<char *>(pLatticeY) +
2535 0 : iGeoOff * iLatticeDataSize)) +
2536 : 0.5;
2537 : }
2538 0 : else if (paiOffset && paiIncrement)
2539 : {
2540 0 : dfGCPPixel = paiOffset[iPixelDim] +
2541 0 : j * paiIncrement[iPixelDim] + 0.5;
2542 0 : dfGCPLine = paiOffset[iLineDim] +
2543 0 : i * paiIncrement[iLineDim] + 0.5;
2544 : }
2545 :
2546 : m_aoGCPs.emplace_back("", "", dfGCPPixel, dfGCPLine, dfGCPX,
2547 0 : dfGCPY);
2548 : }
2549 : }
2550 : }
2551 :
2552 : /* --------------------------------------------------------------------
2553 : */
2554 : /* Establish geolocation metadata, but only if there is no */
2555 : /* lattice. The lattice destroys the regularity of the grid. */
2556 : /* --------------------------------------------------------------------
2557 : */
2558 0 : if (pLatticeX == nullptr && iLatDim != -1 && iLongDim != -1 &&
2559 0 : iPixelDim != -1 && iLineDim != -1)
2560 : {
2561 0 : char *pszWKT = nullptr;
2562 0 : m_oGCPSRS.exportToWkt(&pszWKT);
2563 0 : if (pszWKT)
2564 0 : SetMetadataItem("SRS", pszWKT, "GEOLOCATION");
2565 0 : CPLFree(pszWKT);
2566 :
2567 0 : CPLString osWrk;
2568 : osWrk.Printf("HDF4_EOS:EOS_SWATH_GEOL:\"%s\":%s:%s", pszFilename,
2569 0 : pszSubdatasetName, papszGeolocations[iLongDim]);
2570 0 : SetMetadataItem("X_DATASET", osWrk, "GEOLOCATION");
2571 0 : SetMetadataItem("X_BAND", "1", "GEOLOCATION");
2572 :
2573 : osWrk.Printf("HDF4_EOS:EOS_SWATH_GEOL:\"%s\":%s:%s", pszFilename,
2574 0 : pszSubdatasetName, papszGeolocations[iLatDim]);
2575 0 : SetMetadataItem("Y_DATASET", osWrk, "GEOLOCATION");
2576 0 : SetMetadataItem("Y_BAND", "1", "GEOLOCATION");
2577 :
2578 0 : if (paiOffset && paiIncrement)
2579 : {
2580 0 : osWrk.Printf("%ld", static_cast<long>(paiOffset[iPixelDim]));
2581 0 : SetMetadataItem("PIXEL_OFFSET", osWrk, "GEOLOCATION");
2582 0 : osWrk.Printf("%ld", static_cast<long>(paiIncrement[iPixelDim]));
2583 0 : SetMetadataItem("PIXEL_STEP", osWrk, "GEOLOCATION");
2584 :
2585 0 : osWrk.Printf("%ld", static_cast<long>(paiOffset[iLineDim]));
2586 0 : SetMetadataItem("LINE_OFFSET", osWrk, "GEOLOCATION");
2587 0 : osWrk.Printf("%ld", static_cast<long>(paiIncrement[iLineDim]));
2588 0 : SetMetadataItem("LINE_STEP", osWrk, "GEOLOCATION");
2589 : }
2590 : }
2591 :
2592 : /* --------------------------------------------------------------------
2593 : */
2594 : /* Cleanup */
2595 : /* --------------------------------------------------------------------
2596 : */
2597 0 : CPLFree(pLatticeX);
2598 0 : CPLFree(pLatticeY);
2599 0 : CPLFree(pLat);
2600 0 : CPLFree(pLong);
2601 :
2602 0 : if (iGCPStepX == 0)
2603 : {
2604 0 : m_oGCPSRS.Clear();
2605 : }
2606 : }
2607 :
2608 : /* -------------------------------------------------------------------- */
2609 : /* Cleanup */
2610 : /* -------------------------------------------------------------------- */
2611 0 : CPLFree(paiOffset);
2612 0 : CPLFree(paiIncrement);
2613 0 : CPLFree(pszGeoList);
2614 0 : CSLDestroy(papszGeolocations);
2615 :
2616 0 : return TRUE;
2617 : }
2618 :
2619 : /************************************************************************/
2620 : /* Open() */
2621 : /************************************************************************/
2622 :
2623 301 : GDALDataset *HDF4ImageDataset::Open(GDALOpenInfo *poOpenInfo)
2624 : {
2625 301 : if (!HDF4ImageDatasetIdentify(poOpenInfo))
2626 0 : return nullptr;
2627 :
2628 : /* -------------------------------------------------------------------- */
2629 : /* Create a corresponding GDALDataset. */
2630 : /* -------------------------------------------------------------------- */
2631 301 : if (poOpenInfo->fpL != nullptr)
2632 : {
2633 0 : VSIFCloseL(poOpenInfo->fpL);
2634 0 : poOpenInfo->fpL = nullptr;
2635 : }
2636 :
2637 301 : HDF4ImageDataset *poDS = new HDF4ImageDataset();
2638 602 : CPLMutexHolderD(&hHDF4Mutex);
2639 :
2640 602 : char **papszSubdatasetName = CSLTokenizeString2(
2641 301 : poOpenInfo->pszFilename, ":",
2642 : CSLT_HONOURSTRINGS | CSLT_PRESERVEQUOTES | CSLT_PRESERVEESCAPES);
2643 301 : if (CSLCount(papszSubdatasetName) != 4 &&
2644 301 : CSLCount(papszSubdatasetName) != 5 &&
2645 0 : CSLCount(papszSubdatasetName) != 6)
2646 : {
2647 0 : CSLDestroy(papszSubdatasetName);
2648 : // Release mutex otherwise we deadlock with GDALDataset own mutex.
2649 0 : CPLReleaseMutex(hHDF4Mutex);
2650 0 : delete poDS;
2651 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
2652 0 : return nullptr;
2653 : }
2654 :
2655 : {
2656 : // Un-quote filename
2657 301 : size_t nLenPart2 = strlen(papszSubdatasetName[2]);
2658 301 : if (papszSubdatasetName[2][0] == '"' &&
2659 301 : papszSubdatasetName[2][nLenPart2 - 1] == '"')
2660 : {
2661 301 : papszSubdatasetName[2][nLenPart2 - 1] = 0;
2662 301 : memmove(papszSubdatasetName[2], papszSubdatasetName[2] + 1,
2663 : nLenPart2 - 1);
2664 : }
2665 : }
2666 :
2667 : /* -------------------------------------------------------------------- */
2668 : /* Check for drive name in windows HDF4_xx:TYPE:D:\... */
2669 : /* -------------------------------------------------------------------- */
2670 301 : if (strlen(papszSubdatasetName[2]) == 1)
2671 : {
2672 0 : const size_t nLen = 2 + strlen(papszSubdatasetName[3]) + 1;
2673 0 : char *pszFilename = static_cast<char *>(CPLMalloc(nLen));
2674 0 : snprintf(pszFilename, nLen, "%s:%s", papszSubdatasetName[2],
2675 0 : papszSubdatasetName[3]);
2676 0 : CPLFree(papszSubdatasetName[2]);
2677 0 : CPLFree(papszSubdatasetName[3]);
2678 0 : papszSubdatasetName[2] = pszFilename;
2679 :
2680 : /* Move following arguments one rank upper */
2681 0 : papszSubdatasetName[3] = papszSubdatasetName[4];
2682 0 : if (papszSubdatasetName[4] != nullptr)
2683 : {
2684 0 : papszSubdatasetName[4] = papszSubdatasetName[5];
2685 0 : papszSubdatasetName[5] = nullptr;
2686 : }
2687 : }
2688 :
2689 602 : for (int i = 3; papszSubdatasetName[i] != nullptr; i++)
2690 : {
2691 : // Un-quote and unescape components after filename
2692 301 : size_t nLenPart = strlen(papszSubdatasetName[i]);
2693 301 : if (papszSubdatasetName[i][0] == '"' &&
2694 0 : papszSubdatasetName[i][nLenPart - 1] == '"')
2695 : {
2696 0 : CPLString osStr(papszSubdatasetName[i]);
2697 0 : osStr.replaceAll("\\\\", '\\');
2698 0 : osStr.replaceAll("\\\"", '"');
2699 0 : if (osStr[0] == '"' && osStr.back() == '"')
2700 : {
2701 0 : osStr = osStr.substr(1, osStr.size() - 2);
2702 0 : CPLFree(papszSubdatasetName[i]);
2703 0 : papszSubdatasetName[i] = CPLStrdup(osStr);
2704 : }
2705 : }
2706 : }
2707 :
2708 301 : poDS->pszFilename = CPLStrdup(papszSubdatasetName[2]);
2709 :
2710 301 : if (EQUAL(papszSubdatasetName[0], "HDF4_SDS"))
2711 299 : poDS->iDatasetType = HDF4_SDS;
2712 2 : else if (EQUAL(papszSubdatasetName[0], "HDF4_GR"))
2713 2 : poDS->iDatasetType = HDF4_GR;
2714 0 : else if (EQUAL(papszSubdatasetName[0], "HDF4_EOS"))
2715 0 : poDS->iDatasetType = HDF4_EOS;
2716 : else
2717 0 : poDS->iDatasetType = HDF4_UNKNOWN;
2718 :
2719 301 : if (EQUAL(papszSubdatasetName[1], "GDAL_HDF4"))
2720 299 : poDS->iSubdatasetType = H4ST_GDAL;
2721 2 : else if (EQUAL(papszSubdatasetName[1], "EOS_GRID"))
2722 0 : poDS->iSubdatasetType = H4ST_EOS_GRID;
2723 2 : else if (EQUAL(papszSubdatasetName[1], "EOS_SWATH"))
2724 0 : poDS->iSubdatasetType = H4ST_EOS_SWATH;
2725 2 : else if (EQUAL(papszSubdatasetName[1], "EOS_SWATH_GEOL"))
2726 0 : poDS->iSubdatasetType = H4ST_EOS_SWATH_GEOL;
2727 2 : else if (EQUAL(papszSubdatasetName[1], "SEAWIFS_L3"))
2728 0 : poDS->iSubdatasetType = H4ST_SEAWIFS_L3;
2729 2 : else if (EQUAL(papszSubdatasetName[1], "HYPERION_L1"))
2730 0 : poDS->iSubdatasetType = H4ST_HYPERION_L1;
2731 : else
2732 2 : poDS->iSubdatasetType = H4ST_UNKNOWN;
2733 :
2734 : /* -------------------------------------------------------------------- */
2735 : /* Is our file still here? */
2736 : /* -------------------------------------------------------------------- */
2737 301 : if (!Hishdf(poDS->pszFilename))
2738 : {
2739 0 : CSLDestroy(papszSubdatasetName);
2740 : // Release mutex otherwise we deadlock with GDALDataset own mutex.
2741 0 : CPLReleaseMutex(hHDF4Mutex);
2742 0 : delete poDS;
2743 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
2744 0 : return nullptr;
2745 : }
2746 :
2747 : /* -------------------------------------------------------------------- */
2748 : /* Collect the remain (post filename) components to treat as */
2749 : /* the subdataset name. */
2750 : /* -------------------------------------------------------------------- */
2751 602 : CPLString osSubdatasetName = papszSubdatasetName[3];
2752 301 : if (papszSubdatasetName[4] != nullptr)
2753 : {
2754 0 : osSubdatasetName += ":";
2755 0 : osSubdatasetName += papszSubdatasetName[4];
2756 : }
2757 :
2758 : /* -------------------------------------------------------------------- */
2759 : /* Try opening the dataset. */
2760 : /* -------------------------------------------------------------------- */
2761 301 : double dfNoData = 0.0;
2762 301 : double dfScale = 1.0;
2763 301 : double dfOffset = 0.0;
2764 301 : bool bNoDataSet = false;
2765 301 : bool bHaveScale = false;
2766 301 : bool bHaveOffset = false;
2767 301 : const char *pszUnits = nullptr;
2768 301 : const char *pszDescription = nullptr;
2769 :
2770 : /* -------------------------------------------------------------------- */
2771 : /* Select SDS or GR to read from. */
2772 : /* -------------------------------------------------------------------- */
2773 301 : if (poDS->iDatasetType == HDF4_EOS)
2774 : {
2775 0 : if (papszSubdatasetName[4] == nullptr)
2776 : {
2777 : // Release mutex. Otherwise it will deadlock with GDALDataset's own
2778 : // mutex.
2779 0 : CPLReleaseMutex(hHDF4Mutex);
2780 0 : delete poDS;
2781 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
2782 0 : return nullptr;
2783 : }
2784 0 : poDS->pszSubdatasetName = CPLStrdup(papszSubdatasetName[3]);
2785 0 : poDS->pszFieldName = CPLStrdup(papszSubdatasetName[4]);
2786 : }
2787 : else
2788 : {
2789 301 : CPLAssert(papszSubdatasetName[3]);
2790 301 : poDS->iDataset = atoi(papszSubdatasetName[3]);
2791 : }
2792 301 : CSLDestroy(papszSubdatasetName);
2793 :
2794 301 : switch (poDS->iDatasetType)
2795 : {
2796 0 : case HDF4_EOS:
2797 : {
2798 0 : switch (poDS->iSubdatasetType)
2799 : {
2800 :
2801 : /* --------------------------------------------------------------------
2802 : */
2803 : /* HDF-EOS Swath. */
2804 : /* --------------------------------------------------------------------
2805 : */
2806 0 : case H4ST_EOS_SWATH:
2807 : case H4ST_EOS_SWATH_GEOL:
2808 : {
2809 0 : if (poOpenInfo->eAccess == GA_ReadOnly)
2810 0 : poDS->hHDF4 = SWopen(poDS->pszFilename, DFACC_READ);
2811 : else
2812 0 : poDS->hHDF4 = SWopen(poDS->pszFilename, DFACC_WRITE);
2813 :
2814 0 : if (poDS->hHDF4 <= 0)
2815 : {
2816 0 : CPLDebug("HDF4Image",
2817 : "Can't open file \"%s\" for swath reading",
2818 : poDS->pszFilename);
2819 : // Release mutex otherwise we deadlock with GDALDataset
2820 : // own mutex.
2821 0 : CPLReleaseMutex(hHDF4Mutex);
2822 0 : delete poDS;
2823 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
2824 0 : return nullptr;
2825 : }
2826 :
2827 : const int32 hSW =
2828 0 : SWattach(poDS->hHDF4, poDS->pszSubdatasetName);
2829 0 : if (hSW < 0)
2830 : {
2831 0 : CPLDebug("HDF4Image", "Can't attach to subdataset %s",
2832 : poDS->pszSubdatasetName);
2833 : // Release mutex otherwise we deadlock with GDALDataset
2834 : // own mutex.
2835 0 : CPLReleaseMutex(hHDF4Mutex);
2836 0 : delete poDS;
2837 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
2838 0 : return nullptr;
2839 : }
2840 :
2841 : /* --------------------------------------------------------------------
2842 : */
2843 : /* Decode the dimension map. */
2844 : /* --------------------------------------------------------------------
2845 : */
2846 0 : int32 nStrBufSize = 0;
2847 :
2848 0 : if (SWnentries(hSW, HDFE_NENTDIM, &nStrBufSize) < 0 ||
2849 0 : nStrBufSize <= 0)
2850 : {
2851 0 : CPLDebug("HDF4Image",
2852 : "Can't read a number of dimension maps.");
2853 : // Release mutex otherwise we deadlock with GDALDataset
2854 : // own mutex.
2855 0 : CPLReleaseMutex(hHDF4Mutex);
2856 0 : delete poDS;
2857 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
2858 0 : return nullptr;
2859 : }
2860 :
2861 : char *pszDimList =
2862 0 : static_cast<char *>(CPLMalloc(nStrBufSize + 1));
2863 0 : if (SWfieldinfo(hSW, poDS->pszFieldName, &poDS->iRank,
2864 0 : poDS->aiDimSizes, &poDS->iNumType,
2865 0 : pszDimList) < 0)
2866 : {
2867 0 : CPLDebug("HDF4Image", "Can't read dimension maps.");
2868 0 : CPLFree(pszDimList);
2869 : // Release mutex otherwise we deadlock with GDALDataset
2870 : // own mutex.
2871 0 : CPLReleaseMutex(hHDF4Mutex);
2872 0 : delete poDS;
2873 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
2874 0 : return nullptr;
2875 : }
2876 0 : pszDimList[nStrBufSize] = '\0';
2877 :
2878 : #ifdef DEBUG
2879 0 : CPLDebug("HDF4Image",
2880 : "List of dimensions in swath \"%s\": %s",
2881 : poDS->pszFieldName, pszDimList);
2882 : #endif
2883 :
2884 0 : poDS->GetImageDimensions(pszDimList);
2885 :
2886 : #ifdef DEBUG
2887 0 : CPLDebug("HDF4Image",
2888 : "X dimension is %d, Y dimension is %d",
2889 : poDS->iXDim, poDS->iYDim);
2890 : #endif
2891 :
2892 : /* --------------------------------------------------------------------
2893 : */
2894 : /* Fetch metadata. */
2895 : /* --------------------------------------------------------------------
2896 : */
2897 : // Not H4ST_EOS_SWATH_GEOL.
2898 0 : if (poDS->iSubdatasetType == H4ST_EOS_SWATH)
2899 0 : poDS->GetSwatAttrs(hSW);
2900 :
2901 : /* --------------------------------------------------------------------
2902 : */
2903 : /* Fetch NODATA value. */
2904 : /* --------------------------------------------------------------------
2905 : */
2906 : void *pNoDataValue =
2907 0 : CPLMalloc(poDS->GetDataTypeSize(poDS->iNumType));
2908 0 : if (SWgetfillvalue(hSW, poDS->pszFieldName, pNoDataValue) !=
2909 : -1)
2910 : {
2911 : dfNoData =
2912 0 : poDS->AnyTypeToDouble(poDS->iNumType, pNoDataValue);
2913 0 : bNoDataSet = true;
2914 : }
2915 : else
2916 : {
2917 0 : const char *pszNoData = CSLFetchNameValue(
2918 0 : poDS->papszLocalMetadata, "_FillValue");
2919 0 : if (pszNoData)
2920 : {
2921 0 : dfNoData = CPLAtof(pszNoData);
2922 0 : bNoDataSet = true;
2923 : }
2924 : }
2925 0 : CPLFree(pNoDataValue);
2926 :
2927 : /* --------------------------------------------------------------------
2928 : */
2929 : /* Handle Geolocation processing. */
2930 : /* --------------------------------------------------------------------
2931 : */
2932 : // Not H4ST_SWATH_GEOL.
2933 0 : if (poDS->iSubdatasetType == H4ST_EOS_SWATH)
2934 : {
2935 0 : char **papszDimList = CSLTokenizeString2(
2936 : pszDimList, ",", CSLT_HONOURSTRINGS);
2937 0 : if (!poDS->ProcessSwathGeolocation(hSW, papszDimList))
2938 : {
2939 0 : CPLDebug(
2940 : "HDF4Image",
2941 : "No geolocation available for this swath.");
2942 : }
2943 0 : CSLDestroy(papszDimList);
2944 : }
2945 :
2946 : /* --------------------------------------------------------------------
2947 : */
2948 : /* Cleanup. */
2949 : /* --------------------------------------------------------------------
2950 : */
2951 0 : CPLFree(pszDimList);
2952 0 : SWdetach(hSW);
2953 : }
2954 0 : break;
2955 :
2956 : /* --------------------------------------------------------------------
2957 : */
2958 : /* HDF-EOS Grid. */
2959 : /* --------------------------------------------------------------------
2960 : */
2961 0 : case H4ST_EOS_GRID:
2962 : {
2963 0 : if (poOpenInfo->eAccess == GA_ReadOnly)
2964 0 : poDS->hHDF4 = GDopen(poDS->pszFilename, DFACC_READ);
2965 : else
2966 0 : poDS->hHDF4 = GDopen(poDS->pszFilename, DFACC_WRITE);
2967 :
2968 0 : if (poDS->hHDF4 <= 0)
2969 : {
2970 0 : CPLDebug("HDF4Image",
2971 : "Can't open file \"%s\" for grid reading",
2972 : poDS->pszFilename);
2973 : // Release mutex otherwise we deadlock with GDALDataset
2974 : // own mutex.
2975 0 : CPLReleaseMutex(hHDF4Mutex);
2976 0 : delete poDS;
2977 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
2978 0 : return nullptr;
2979 : }
2980 :
2981 : const int32 hGD =
2982 0 : GDattach(poDS->hHDF4, poDS->pszSubdatasetName);
2983 :
2984 : /* --------------------------------------------------------------------
2985 : */
2986 : /* Decode the dimension map. */
2987 : /* --------------------------------------------------------------------
2988 : */
2989 0 : char szDimList[N_BUF_SIZE] = {};
2990 0 : GDfieldinfo(hGD, poDS->pszFieldName, &poDS->iRank,
2991 0 : poDS->aiDimSizes, &poDS->iNumType, szDimList);
2992 : #ifdef DEBUG
2993 0 : CPLDebug("HDF4Image", "List of dimensions in grid %s: %s",
2994 : poDS->pszFieldName, szDimList);
2995 : #endif
2996 0 : poDS->GetImageDimensions(szDimList);
2997 :
2998 0 : int32 tilecode = 0;
2999 0 : int32 tilerank = 0;
3000 0 : if (GDtileinfo(hGD, poDS->pszFieldName, &tilecode,
3001 0 : &tilerank, nullptr) == 0)
3002 : {
3003 0 : if (tilecode == HDFE_TILE)
3004 : {
3005 : int32 *tiledims = reinterpret_cast<int32 *>(
3006 0 : CPLCalloc(tilerank, sizeof(int32)));
3007 0 : GDtileinfo(hGD, poDS->pszFieldName, &tilecode,
3008 : &tilerank, tiledims);
3009 0 : if ((tilerank == 2) && (poDS->iRank == tilerank))
3010 : {
3011 0 : poDS->nBlockPreferredXSize = tiledims[1];
3012 0 : poDS->nBlockPreferredYSize = tiledims[0];
3013 0 : poDS->bReadTile = true;
3014 : #ifdef DEBUG
3015 0 : CPLDebug("HDF4_EOS:EOS_GRID:",
3016 : "tilerank in grid %s: %d",
3017 : poDS->pszFieldName,
3018 : static_cast<int>(tilerank));
3019 0 : CPLDebug("HDF4_EOS:EOS_GRID:",
3020 : "tiledimens in grid %s: %d,%d",
3021 : poDS->pszFieldName,
3022 : static_cast<int>(tiledims[0]),
3023 0 : static_cast<int>(tiledims[1]));
3024 : #endif
3025 : }
3026 : #ifdef DEBUG
3027 : else
3028 : {
3029 0 : CPLDebug(
3030 : "HDF4_EOS:EOS_GRID:",
3031 : "tilerank in grid %s: %d not supported",
3032 : poDS->pszFieldName,
3033 : static_cast<int>(tilerank));
3034 : }
3035 : #endif
3036 0 : CPLFree(tiledims);
3037 : }
3038 : else
3039 : {
3040 : #ifdef DEBUG
3041 0 : CPLDebug("HDF4_EOS:EOS_GRID:",
3042 : "tilecode == HDFE_NOTILE in grid %s: %d",
3043 : poDS->pszFieldName,
3044 0 : static_cast<int>(poDS->iRank));
3045 : #endif
3046 : }
3047 : }
3048 :
3049 : #ifdef DEBUG
3050 : else
3051 : {
3052 0 : CPLDebug("HDF4_EOS:EOS_GRID:", "ERROR GDtileinfo %s",
3053 : poDS->pszFieldName);
3054 : }
3055 : #endif
3056 :
3057 : /* --------------------------------------------------------------------
3058 : */
3059 : /* Fetch projection information */
3060 : /* --------------------------------------------------------------------
3061 : */
3062 0 : int32 iProjCode = 0;
3063 0 : int32 iZoneCode = 0;
3064 0 : int32 iSphereCode = 0;
3065 : double adfProjParams[15];
3066 :
3067 0 : if (GDprojinfo(hGD, &iProjCode, &iZoneCode, &iSphereCode,
3068 0 : adfProjParams) >= 0)
3069 : {
3070 : #ifdef DEBUG
3071 0 : CPLDebug("HDF4Image",
3072 : "Grid projection: "
3073 : "projection code: %ld, zone code %ld, "
3074 : "sphere code %ld",
3075 : static_cast<long>(iProjCode),
3076 : static_cast<long>(iZoneCode),
3077 : static_cast<long>(iSphereCode));
3078 : #endif
3079 0 : poDS->m_oSRS.importFromUSGS(iProjCode, iZoneCode,
3080 : adfProjParams, iSphereCode,
3081 : USGS_ANGLE_RADIANS);
3082 : }
3083 :
3084 : /* --------------------------------------------------------------------
3085 : */
3086 : /* Fetch geotransformation matrix */
3087 : /* --------------------------------------------------------------------
3088 : */
3089 0 : int32 nXSize = 0;
3090 0 : int32 nYSize = 0;
3091 : double adfUpLeft[2];
3092 : double adfLowRight[2];
3093 0 : if (GDgridinfo(hGD, &nXSize, &nYSize, adfUpLeft,
3094 0 : adfLowRight) >= 0)
3095 : {
3096 : #ifdef DEBUG
3097 0 : CPLDebug("HDF4Image",
3098 : "Grid geolocation: "
3099 : "top left X %f, top left Y %f, "
3100 : "low right X %f, low right Y %f, "
3101 : "cols %ld, rows %ld",
3102 : adfUpLeft[0], adfUpLeft[1], adfLowRight[0],
3103 : adfLowRight[1], static_cast<long>(nXSize),
3104 : static_cast<long>(nYSize));
3105 : #endif
3106 0 : if (iProjCode)
3107 : {
3108 : // For projected systems coordinates are in meters.
3109 0 : poDS->m_gt[1] =
3110 0 : (adfLowRight[0] - adfUpLeft[0]) / nXSize;
3111 0 : poDS->m_gt[5] =
3112 0 : (adfLowRight[1] - adfUpLeft[1]) / nYSize;
3113 0 : poDS->m_gt[0] = adfUpLeft[0];
3114 0 : poDS->m_gt[3] = adfUpLeft[1];
3115 : }
3116 : else
3117 : {
3118 : // Handle angular geographic coordinates here.
3119 0 : poDS->m_gt[1] = (CPLPackedDMSToDec(adfLowRight[0]) -
3120 0 : CPLPackedDMSToDec(adfUpLeft[0])) /
3121 : nXSize;
3122 0 : poDS->m_gt[5] = (CPLPackedDMSToDec(adfLowRight[1]) -
3123 0 : CPLPackedDMSToDec(adfUpLeft[1])) /
3124 : nYSize;
3125 0 : poDS->m_gt[0] = CPLPackedDMSToDec(adfUpLeft[0]);
3126 0 : poDS->m_gt[3] = CPLPackedDMSToDec(adfUpLeft[1]);
3127 : }
3128 0 : poDS->m_gt[2] = 0.0;
3129 0 : poDS->m_gt[4] = 0.0;
3130 0 : poDS->bHasGeoTransform = true;
3131 : }
3132 :
3133 : /* --------------------------------------------------------------------
3134 : */
3135 : /* Fetch metadata. */
3136 : /* --------------------------------------------------------------------
3137 : */
3138 0 : poDS->GetGridAttrs(hGD);
3139 :
3140 0 : GDdetach(hGD);
3141 :
3142 : /* --------------------------------------------------------------------
3143 : */
3144 : /* Fetch NODATA value. */
3145 : /* --------------------------------------------------------------------
3146 : */
3147 : void *pNoDataValue =
3148 0 : CPLMalloc(poDS->GetDataTypeSize(poDS->iNumType));
3149 0 : if (GDgetfillvalue(hGD, poDS->pszFieldName, pNoDataValue) !=
3150 : -1)
3151 : {
3152 : dfNoData =
3153 0 : poDS->AnyTypeToDouble(poDS->iNumType, pNoDataValue);
3154 0 : bNoDataSet = true;
3155 : }
3156 : else
3157 : {
3158 0 : const char *pszNoData = CSLFetchNameValue(
3159 0 : poDS->papszLocalMetadata, "_FillValue");
3160 0 : if (pszNoData)
3161 : {
3162 0 : dfNoData = CPLAtof(pszNoData);
3163 0 : bNoDataSet = true;
3164 : }
3165 : }
3166 0 : CPLFree(pNoDataValue);
3167 : }
3168 0 : break;
3169 :
3170 0 : default:
3171 0 : break;
3172 : }
3173 :
3174 : /* --------------------------------------------------------------------
3175 : */
3176 : /* Fetch unit type, scale, offset and description */
3177 : /* Should be similar among various HDF-EOS kinds. */
3178 : /* --------------------------------------------------------------------
3179 : */
3180 : {
3181 : const char *pszTmp =
3182 0 : CSLFetchNameValue(poDS->papszLocalMetadata, "scale_factor");
3183 0 : if (pszTmp)
3184 : {
3185 0 : dfScale = CPLAtof(pszTmp);
3186 : // Some producers (i.e. lndcsm from LEDAPS) emit
3187 : // files with scale_factor=0 which is crazy to carry
3188 : // through.
3189 0 : if (dfScale == 0.0)
3190 0 : dfScale = 1.0;
3191 :
3192 0 : bHaveScale = (dfScale != 0.0);
3193 : }
3194 :
3195 : pszTmp =
3196 0 : CSLFetchNameValue(poDS->papszLocalMetadata, "add_offset");
3197 0 : if (pszTmp)
3198 : {
3199 0 : dfOffset = CPLAtof(pszTmp);
3200 0 : bHaveOffset = true;
3201 : }
3202 :
3203 0 : pszUnits = CSLFetchNameValue(poDS->papszLocalMetadata, "units");
3204 : pszDescription =
3205 0 : CSLFetchNameValue(poDS->papszLocalMetadata, "long_name");
3206 : }
3207 : }
3208 0 : break;
3209 :
3210 : /* --------------------------------------------------------------------
3211 : */
3212 : /* 'Plain' HDF scientific datasets. */
3213 : /* --------------------------------------------------------------------
3214 : */
3215 299 : case HDF4_SDS:
3216 : {
3217 : #ifdef HDF4_HAS_MAXOPENFILES
3218 : // Attempt to increase maximum number of opened HDF files
3219 299 : intn nCurrMax = 0;
3220 299 : intn nSysLimit = 0;
3221 :
3222 598 : if (SDget_maxopenfiles(&nCurrMax, &nSysLimit) >= 0 &&
3223 299 : nCurrMax < nSysLimit)
3224 : {
3225 0 : SDreset_maxopenfiles(nSysLimit);
3226 : }
3227 : #endif // HDF4_HAS_MAXOPENFILES
3228 :
3229 299 : if (poOpenInfo->eAccess == GA_ReadOnly)
3230 299 : poDS->hHDF4 = Hopen(poDS->pszFilename, DFACC_READ, 0);
3231 : else
3232 0 : poDS->hHDF4 = Hopen(poDS->pszFilename, DFACC_WRITE, 0);
3233 :
3234 299 : if (poDS->hHDF4 <= 0)
3235 : {
3236 : // Release mutex otherwise we deadlock with GDALDataset own
3237 : // mutex.
3238 0 : CPLReleaseMutex(hHDF4Mutex);
3239 0 : delete poDS;
3240 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3241 0 : return nullptr;
3242 : }
3243 :
3244 299 : poDS->hSD = SDstart(poDS->pszFilename, DFACC_READ);
3245 299 : if (poDS->hSD == -1)
3246 : {
3247 : // Release mutex otherwise we deadlock with GDALDataset own
3248 : // mutex.
3249 0 : CPLReleaseMutex(hHDF4Mutex);
3250 0 : delete poDS;
3251 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3252 0 : return nullptr;
3253 : }
3254 :
3255 299 : if (poDS->ReadGlobalAttributes(poDS->hSD) != CE_None)
3256 : {
3257 : // Release mutex otherwise we deadlock with GDALDataset own
3258 : // mutex.
3259 0 : CPLReleaseMutex(hHDF4Mutex);
3260 0 : delete poDS;
3261 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3262 0 : return nullptr;
3263 : }
3264 :
3265 299 : int32 nDatasets = 0;
3266 299 : int32 l_nAttrs = 0;
3267 299 : if (SDfileinfo(poDS->hSD, &nDatasets, &l_nAttrs) != 0)
3268 : {
3269 : // Release mutex otherwise we deadlock with GDALDataset own
3270 : // mutex.
3271 0 : CPLReleaseMutex(hHDF4Mutex);
3272 0 : delete poDS;
3273 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3274 0 : return nullptr;
3275 : }
3276 :
3277 299 : if (poDS->iDataset < 0 || poDS->iDataset >= nDatasets)
3278 : {
3279 0 : CPLError(CE_Failure, CPLE_AppDefined,
3280 : "Subdataset index should be between 0 and %ld",
3281 0 : static_cast<long int>(nDatasets) - 1);
3282 : // Release mutex otherwise we deadlock with GDALDataset own
3283 : // mutex.
3284 0 : CPLReleaseMutex(hHDF4Mutex);
3285 0 : delete poDS;
3286 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3287 0 : return nullptr;
3288 : }
3289 :
3290 299 : memset(poDS->aiDimSizes, 0, sizeof(int32) * H4_MAX_VAR_DIMS);
3291 299 : const int32 iSDS = SDselect(poDS->hSD, poDS->iDataset);
3292 299 : SDgetinfo(iSDS, poDS->szName, &poDS->iRank, poDS->aiDimSizes,
3293 : &poDS->iNumType, &poDS->nAttrs);
3294 :
3295 : // We will duplicate global metadata for every subdataset.
3296 299 : poDS->papszLocalMetadata = CSLDuplicate(poDS->papszGlobalMetadata);
3297 :
3298 299 : for (int32 iAttribute = 0; iAttribute < poDS->nAttrs; iAttribute++)
3299 : {
3300 0 : char szAttrName[H4_MAX_NC_NAME] = {};
3301 0 : int32 iAttrNumType = 0;
3302 0 : int32 nValues = 0;
3303 0 : SDattrinfo(iSDS, iAttribute, szAttrName, &iAttrNumType,
3304 : &nValues);
3305 0 : poDS->papszLocalMetadata = poDS->TranslateHDF4Attributes(
3306 : iSDS, iAttribute, szAttrName, iAttrNumType, nValues,
3307 : poDS->papszLocalMetadata);
3308 : }
3309 299 : poDS->SetMetadata(poDS->papszLocalMetadata, "");
3310 299 : SDendaccess(iSDS);
3311 :
3312 : #ifdef DEBUG
3313 299 : CPLDebug("HDF4Image",
3314 : "aiDimSizes[0]=%ld, aiDimSizes[1]=%ld, "
3315 : "aiDimSizes[2]=%ld, aiDimSizes[3]=%ld",
3316 299 : static_cast<long>(poDS->aiDimSizes[0]),
3317 299 : static_cast<long>(poDS->aiDimSizes[1]),
3318 299 : static_cast<long>(poDS->aiDimSizes[2]),
3319 299 : static_cast<long>(poDS->aiDimSizes[3]));
3320 : #endif
3321 299 : switch (poDS->iRank)
3322 : {
3323 0 : case 1:
3324 0 : poDS->nBandCount = 1;
3325 0 : poDS->iXDim = 0;
3326 0 : poDS->iYDim = -1;
3327 0 : break;
3328 135 : case 2:
3329 135 : poDS->nBandCount = 1;
3330 135 : poDS->iXDim = 1;
3331 135 : poDS->iYDim = 0;
3332 135 : break;
3333 164 : case 3:
3334 : // TODO: We should probably remove the following test as
3335 : // there are valid datasets where the height is lower than
3336 : // the band number. For example:
3337 : // http://www.iapmw.unibe.ch/research/projects/FriOWL/data/otd/LISOTD_HRAC_V2.2.hdf
3338 : // which is a 720x360 x 365 bands.
3339 : // Use a HACK for now.
3340 164 : if (poDS->aiDimSizes[0] < poDS->aiDimSizes[2] &&
3341 0 : !(poDS->aiDimSizes[0] == 360 &&
3342 0 : poDS->aiDimSizes[1] == 720 &&
3343 0 : poDS->aiDimSizes[2] == 365))
3344 : {
3345 0 : poDS->iBandDim = 0;
3346 0 : poDS->iXDim = 2;
3347 0 : poDS->iYDim = 1;
3348 : }
3349 : else
3350 : {
3351 164 : if (poDS->aiDimSizes[1] <= poDS->aiDimSizes[0] &&
3352 164 : poDS->aiDimSizes[1] <= poDS->aiDimSizes[2])
3353 : {
3354 0 : poDS->iBandDim = 1;
3355 0 : poDS->iXDim = 2;
3356 0 : poDS->iYDim = 0;
3357 : }
3358 : else
3359 : {
3360 164 : poDS->iBandDim = 2;
3361 164 : poDS->iXDim = 1;
3362 164 : poDS->iYDim = 0;
3363 : }
3364 : }
3365 164 : poDS->nBandCount = poDS->aiDimSizes[poDS->iBandDim];
3366 164 : break;
3367 0 : case 4: // FIXME
3368 0 : poDS->nBandCount =
3369 0 : poDS->aiDimSizes[2] * poDS->aiDimSizes[3];
3370 0 : break;
3371 0 : default:
3372 0 : break;
3373 : }
3374 :
3375 : // We preset this because CaptureNRLGeoTransform needs it.
3376 299 : poDS->nRasterXSize = poDS->aiDimSizes[poDS->iXDim];
3377 299 : if (poDS->iYDim >= 0)
3378 299 : poDS->nRasterYSize = poDS->aiDimSizes[poDS->iYDim];
3379 : else
3380 0 : poDS->nRasterYSize = 1;
3381 :
3382 : // Special case projection info for NRL generated files.
3383 598 : const char *pszMapProjectionSystem = CSLFetchNameValue(
3384 299 : poDS->papszGlobalMetadata, "mapProjectionSystem");
3385 299 : if (pszMapProjectionSystem != nullptr &&
3386 0 : EQUAL(pszMapProjectionSystem, "NRL(USGS)"))
3387 : {
3388 0 : poDS->CaptureNRLGeoTransform();
3389 : }
3390 :
3391 : // Special case for coastwatch hdf files.
3392 299 : if (CSLFetchNameValue(poDS->papszGlobalMetadata, "gctp_sys") !=
3393 : nullptr)
3394 0 : poDS->CaptureCoastwatchGCTPInfo();
3395 :
3396 : // Special case for MODIS geolocation
3397 299 : poDS->ProcessModisSDSGeolocation();
3398 :
3399 : // Special case for NASA/CCRS Landsat in HDF
3400 299 : poDS->CaptureL1GMTLInfo();
3401 : }
3402 299 : break;
3403 :
3404 : /* --------------------------------------------------------------------
3405 : */
3406 : /* 'Plain' HDF rasters. */
3407 : /* --------------------------------------------------------------------
3408 : */
3409 2 : case HDF4_GR:
3410 : {
3411 : // Attempt to increase maximum number of opened HDF files.
3412 : #ifdef HDF4_HAS_MAXOPENFILES
3413 2 : intn nCurrMax = 0;
3414 2 : intn nSysLimit = 0;
3415 :
3416 4 : if (SDget_maxopenfiles(&nCurrMax, &nSysLimit) >= 0 &&
3417 2 : nCurrMax < nSysLimit)
3418 : {
3419 0 : SDreset_maxopenfiles(nSysLimit);
3420 : }
3421 : #endif // HDF4_HAS_MAXOPENFILES
3422 :
3423 2 : if (poOpenInfo->eAccess == GA_ReadOnly)
3424 2 : poDS->hHDF4 = Hopen(poDS->pszFilename, DFACC_READ, 0);
3425 : else
3426 0 : poDS->hHDF4 = Hopen(poDS->pszFilename, DFACC_WRITE, 0);
3427 :
3428 2 : if (poDS->hHDF4 <= 0)
3429 : {
3430 : // Release mutex otherwise we deadlock with GDALDataset own
3431 : // mutex.
3432 0 : CPLReleaseMutex(hHDF4Mutex);
3433 0 : delete poDS;
3434 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3435 0 : return nullptr;
3436 : }
3437 :
3438 2 : poDS->hGR = GRstart(poDS->hHDF4);
3439 2 : if (poDS->hGR == -1)
3440 : {
3441 : // Release mutex otherwise we deadlock with GDALDataset own
3442 : // mutex.
3443 0 : CPLReleaseMutex(hHDF4Mutex);
3444 0 : delete poDS;
3445 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3446 0 : return nullptr;
3447 : }
3448 :
3449 2 : poDS->iGR = GRselect(poDS->hGR, poDS->iDataset);
3450 4 : if (GRgetiminfo(poDS->iGR, poDS->szName, &poDS->iRank,
3451 : &poDS->iNumType, &poDS->iInterlaceMode,
3452 2 : poDS->aiDimSizes, &poDS->nAttrs) != 0)
3453 : {
3454 : // Release mutex otherwise we deadlock with GDALDataset own
3455 : // mutex.
3456 0 : CPLReleaseMutex(hHDF4Mutex);
3457 0 : delete poDS;
3458 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3459 0 : return nullptr;
3460 : }
3461 :
3462 : // We will duplicate global metadata for every subdataset.
3463 2 : poDS->papszLocalMetadata = CSLDuplicate(poDS->papszGlobalMetadata);
3464 :
3465 2 : poDS->SetMetadata(poDS->papszLocalMetadata, "");
3466 : // Read colour table
3467 :
3468 2 : poDS->iPal = GRgetlutid(poDS->iGR, 0);
3469 2 : if (poDS->iPal != -1)
3470 : {
3471 2 : GRgetlutinfo(poDS->iPal, &poDS->nComps, &poDS->iPalDataType,
3472 : &poDS->iPalInterlaceMode, &poDS->nPalEntries);
3473 1 : if (poDS->nPalEntries && poDS->nComps == 3 &&
3474 1 : GDALGetDataTypeSizeBytes(GetDataType(poDS->iPalDataType)) ==
3475 3 : 1 &&
3476 1 : poDS->nPalEntries <= 256)
3477 : {
3478 1 : GRreadlut(poDS->iPal, poDS->aiPaletteData);
3479 1 : poDS->poColorTable = new GDALColorTable();
3480 : GDALColorEntry oEntry;
3481 257 : for (int i = 0;
3482 257 : i < std::min(static_cast<int>(poDS->nPalEntries),
3483 257 : N_COLOR_ENTRIES);
3484 : i++)
3485 : {
3486 256 : oEntry.c1 = poDS->aiPaletteData[i][0];
3487 256 : oEntry.c2 = poDS->aiPaletteData[i][1];
3488 256 : oEntry.c3 = poDS->aiPaletteData[i][2];
3489 256 : oEntry.c4 = 255;
3490 :
3491 256 : poDS->poColorTable->SetColorEntry(i, &oEntry);
3492 : }
3493 : }
3494 : }
3495 :
3496 2 : poDS->iXDim = 0;
3497 2 : poDS->iYDim = 1;
3498 2 : poDS->nBandCount = poDS->iRank;
3499 2 : break;
3500 : }
3501 0 : default:
3502 : // Release mutex otherwise we'll deadlock with GDALDataset own mutex
3503 0 : CPLReleaseMutex(hHDF4Mutex);
3504 0 : delete poDS;
3505 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3506 0 : return nullptr;
3507 : }
3508 :
3509 301 : poDS->nRasterXSize = poDS->aiDimSizes[poDS->iXDim];
3510 301 : if (poDS->iYDim >= 0)
3511 301 : poDS->nRasterYSize = poDS->aiDimSizes[poDS->iYDim];
3512 : else
3513 0 : poDS->nRasterYSize = 1;
3514 :
3515 301 : if (poDS->iSubdatasetType == H4ST_HYPERION_L1)
3516 : {
3517 : // XXX: Hyperion SDSs has Height x Bands x Width dimensions scheme
3518 0 : if (poDS->iRank > 2)
3519 : {
3520 0 : poDS->nBandCount = poDS->aiDimSizes[1];
3521 0 : poDS->nRasterXSize = poDS->aiDimSizes[2];
3522 0 : poDS->nRasterYSize = poDS->aiDimSizes[0];
3523 : }
3524 : else
3525 : {
3526 0 : poDS->nBandCount = poDS->aiDimSizes[0];
3527 0 : poDS->nRasterXSize = poDS->aiDimSizes[1];
3528 0 : poDS->nRasterYSize = 1;
3529 : }
3530 : }
3531 :
3532 : /* -------------------------------------------------------------------- */
3533 : /* Create band information objects. */
3534 : /* -------------------------------------------------------------------- */
3535 635 : for (int i = 1; i <= poDS->nBandCount; i++)
3536 : {
3537 : HDF4ImageRasterBand *poBand =
3538 334 : new HDF4ImageRasterBand(poDS, i, poDS->GetDataType(poDS->iNumType));
3539 334 : poDS->SetBand(i, poBand);
3540 :
3541 334 : if (bNoDataSet)
3542 0 : poBand->SetNoDataValue(dfNoData);
3543 334 : if (bHaveScale)
3544 : {
3545 0 : poBand->bHaveScale = true;
3546 0 : poBand->dfScale = dfScale;
3547 : }
3548 334 : if (bHaveOffset)
3549 : {
3550 0 : poBand->bHaveOffset = true;
3551 0 : poBand->dfOffset = dfOffset;
3552 : }
3553 334 : if (pszUnits)
3554 0 : poBand->osUnitType = pszUnits;
3555 334 : if (pszDescription)
3556 0 : poBand->SetDescription(pszDescription);
3557 : }
3558 :
3559 : /* -------------------------------------------------------------------- */
3560 : /* Now we will handle particular types of HDF products. Every */
3561 : /* HDF product has its own structure. */
3562 : /* -------------------------------------------------------------------- */
3563 :
3564 301 : switch (poDS->iSubdatasetType)
3565 : {
3566 : /* --------------------------------------------------------------------
3567 : */
3568 : /* HDF, created by GDAL. */
3569 : /* --------------------------------------------------------------------
3570 : */
3571 299 : case H4ST_GDAL:
3572 : {
3573 299 : CPLDebug("HDF4Image", "Input dataset interpreted as GDAL_HDF4");
3574 :
3575 : const char *pszValue =
3576 299 : CSLFetchNameValue(poDS->papszGlobalMetadata, "Projection");
3577 299 : if (pszValue != nullptr)
3578 : {
3579 119 : poDS->m_oSRS.importFromWkt(pszValue);
3580 : }
3581 299 : if ((pszValue = CSLFetchNameValue(poDS->papszGlobalMetadata,
3582 299 : "TransformationMatrix")) !=
3583 : nullptr)
3584 : {
3585 299 : int i = 0;
3586 299 : char *pszString = const_cast<char *>(pszValue);
3587 2093 : while (*pszValue && i < 6)
3588 : {
3589 1794 : poDS->m_gt[i++] = CPLStrtod(pszString, &pszString);
3590 1794 : pszString++;
3591 : }
3592 299 : poDS->bHasGeoTransform = true;
3593 : }
3594 630 : for (int i = 1; i <= poDS->nBands; i++)
3595 : {
3596 331 : if ((pszValue = CSLFetchNameValue(
3597 331 : poDS->papszGlobalMetadata,
3598 331 : CPLSPrintf("BandDesc%d", i))) != nullptr)
3599 36 : poDS->GetRasterBand(i)->SetDescription(pszValue);
3600 : }
3601 630 : for (int i = 1; i <= poDS->nBands; i++)
3602 : {
3603 331 : if ((pszValue = CSLFetchNameValue(
3604 331 : poDS->papszGlobalMetadata,
3605 331 : CPLSPrintf("NoDataValue%d", i))) != nullptr)
3606 36 : poDS->GetRasterBand(i)->SetNoDataValue(CPLAtof(pszValue));
3607 : }
3608 : }
3609 299 : break;
3610 :
3611 : /* --------------------------------------------------------------------
3612 : */
3613 : /* SeaWiFS Level 3 Standard Mapped Image Products. */
3614 : /* Organized similar to MODIS Level 3 products. */
3615 : /* --------------------------------------------------------------------
3616 : */
3617 0 : case H4ST_SEAWIFS_L3:
3618 : {
3619 0 : CPLDebug("HDF4Image", "Input dataset interpreted as SEAWIFS_L3");
3620 :
3621 : // Read band description
3622 0 : for (int i = 1; i <= poDS->nBands; i++)
3623 : {
3624 0 : poDS->GetRasterBand(i)->SetDescription(
3625 0 : CSLFetchNameValue(poDS->papszGlobalMetadata, "Parameter"));
3626 : }
3627 :
3628 : // Read coordinate system and geotransform matrix.
3629 0 : poDS->m_oSRS.SetWellKnownGeogCS("WGS84");
3630 :
3631 0 : if (EQUAL(CSLFetchNameValue(poDS->papszGlobalMetadata,
3632 : "Map Projection"),
3633 : "Equidistant Cylindrical"))
3634 : {
3635 0 : poDS->m_oSRS.SetEquirectangular(0.0, 0.0, 0.0, 0.0);
3636 0 : poDS->m_oSRS.SetLinearUnits(SRS_UL_METER, 1);
3637 : }
3638 :
3639 0 : double dfULX = CPLAtof(CSLFetchNameValue(poDS->papszGlobalMetadata,
3640 0 : "Westernmost Longitude"));
3641 0 : double dfULY = CPLAtof(CSLFetchNameValue(poDS->papszGlobalMetadata,
3642 0 : "Northernmost Latitude"));
3643 0 : double dfLRX = CPLAtof(CSLFetchNameValue(poDS->papszGlobalMetadata,
3644 0 : "Easternmost Longitude"));
3645 0 : double dfLRY = CPLAtof(CSLFetchNameValue(poDS->papszGlobalMetadata,
3646 0 : "Southernmost Latitude"));
3647 0 : poDS->ToGeoref(&dfULX, &dfULY);
3648 0 : poDS->ToGeoref(&dfLRX, &dfLRY);
3649 0 : poDS->m_gt[0] = dfULX;
3650 0 : poDS->m_gt[3] = dfULY;
3651 0 : poDS->m_gt[1] = (dfLRX - dfULX) / poDS->nRasterXSize;
3652 0 : poDS->m_gt[5] = (dfULY - dfLRY) / poDS->nRasterYSize;
3653 0 : if (dfULY > 0) // Northern hemisphere.
3654 0 : poDS->m_gt[5] = -poDS->m_gt[5];
3655 0 : poDS->m_gt[2] = 0.0;
3656 0 : poDS->m_gt[4] = 0.0;
3657 0 : poDS->bHasGeoTransform = true;
3658 : }
3659 0 : break;
3660 :
3661 : /* --------------------------------------------------------------------
3662 : */
3663 : /* Generic SDS */
3664 : /* --------------------------------------------------------------------
3665 : */
3666 2 : case H4ST_UNKNOWN:
3667 : {
3668 : // This is a coastwatch convention.
3669 2 : if (CSLFetchNameValue(poDS->papszLocalMetadata, "missing_value"))
3670 : {
3671 0 : for (int i = 1; i <= poDS->nBands; i++)
3672 : {
3673 0 : poDS->GetRasterBand(i)->SetNoDataValue(
3674 0 : CPLAtof(CSLFetchNameValue(poDS->papszLocalMetadata,
3675 0 : "missing_value")));
3676 : }
3677 : }
3678 :
3679 : // Coastwatch offset and scale.
3680 2 : if (CSLFetchNameValue(poDS->papszLocalMetadata, "scale_factor") &&
3681 0 : CSLFetchNameValue(poDS->papszLocalMetadata, "add_offset"))
3682 : {
3683 0 : for (int i = 1; i <= poDS->nBands; i++)
3684 : {
3685 : HDF4ImageRasterBand *poBand =
3686 : reinterpret_cast<HDF4ImageRasterBand *>(
3687 0 : poDS->GetRasterBand(i));
3688 :
3689 0 : poBand->bHaveScale = true;
3690 0 : poBand->bHaveOffset = true;
3691 0 : poBand->dfScale = CPLAtof(CSLFetchNameValue(
3692 0 : poDS->papszLocalMetadata, "scale_factor"));
3693 : // See #4891 regarding offset interpretation.
3694 : // poBand->dfOffset = -1 * poBand->dfScale *
3695 : // CPLAtof( CSLFetchNameValue(
3696 : // poDS->papszLocalMetadata,
3697 : // "add_offset" ) );
3698 0 : poBand->dfOffset = CPLAtof(CSLFetchNameValue(
3699 0 : poDS->papszLocalMetadata, "add_offset"));
3700 : }
3701 : }
3702 :
3703 : // This is a modis level3 convention (data from ACT)
3704 : // e.g. data/hdf/act/modis/MODAM2004280160000.L3_NOAA_GMX
3705 :
3706 2 : if (CSLFetchNameValue(poDS->papszLocalMetadata, "scalingSlope") &&
3707 0 : CSLFetchNameValue(poDS->papszLocalMetadata, "scalingIntercept"))
3708 : {
3709 0 : CPLString osUnits;
3710 :
3711 0 : if (CSLFetchNameValue(poDS->papszLocalMetadata, "productUnits"))
3712 : {
3713 0 : osUnits = CSLFetchNameValue(poDS->papszLocalMetadata,
3714 0 : "productUnits");
3715 : }
3716 :
3717 0 : for (int i = 1; i <= poDS->nBands; i++)
3718 : {
3719 : HDF4ImageRasterBand *poBand =
3720 : reinterpret_cast<HDF4ImageRasterBand *>(
3721 0 : poDS->GetRasterBand(i));
3722 :
3723 0 : poBand->bHaveScale = true;
3724 0 : poBand->bHaveOffset = true;
3725 0 : poBand->dfScale = CPLAtof(CSLFetchNameValue(
3726 0 : poDS->papszLocalMetadata, "scalingSlope"));
3727 0 : poBand->dfOffset = CPLAtof(CSLFetchNameValue(
3728 0 : poDS->papszLocalMetadata, "scalingIntercept"));
3729 :
3730 0 : poBand->osUnitType = osUnits;
3731 : }
3732 : }
3733 : }
3734 2 : break;
3735 :
3736 : /* --------------------------------------------------------------------
3737 : */
3738 : /* Hyperion Level 1. */
3739 : /* --------------------------------------------------------------------
3740 : */
3741 0 : case H4ST_HYPERION_L1:
3742 : {
3743 0 : CPLDebug("HDF4Image", "Input dataset interpreted as HYPERION_L1");
3744 : }
3745 0 : break;
3746 :
3747 0 : default:
3748 :
3749 : #ifdef DEBUG_VERBOSE
3750 : CPLError(CE_Debug, CPLE_AppDefined, "Unknown subdata type %d",
3751 : poDS->iSubdatasetType);
3752 : #endif
3753 0 : break;
3754 : }
3755 :
3756 : /* -------------------------------------------------------------------- */
3757 : /* Setup PAM info for this subdataset */
3758 : /* -------------------------------------------------------------------- */
3759 301 : poDS->SetPhysicalFilename(poDS->pszFilename);
3760 301 : poDS->SetSubdatasetName(osSubdatasetName);
3761 :
3762 : // Release mutex otherwise we'll deadlock with GDALDataset own mutex.
3763 301 : CPLReleaseMutex(hHDF4Mutex);
3764 301 : poDS->TryLoadXML();
3765 :
3766 301 : poDS->oOvManager.Initialize(poDS, ":::VIRTUAL:::");
3767 301 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3768 :
3769 301 : return poDS;
3770 : }
3771 :
3772 : /************************************************************************/
3773 : /* Create() */
3774 : /************************************************************************/
3775 :
3776 176 : GDALDataset *HDF4ImageDataset::Create(const char *pszFilename, int nXSize,
3777 : int nYSize, int nBandsIn,
3778 : GDALDataType eType, char **papszOptions)
3779 :
3780 : {
3781 : /* -------------------------------------------------------------------- */
3782 : /* Create the dataset. */
3783 : /* -------------------------------------------------------------------- */
3784 176 : if (nBandsIn == 0)
3785 : {
3786 1 : CPLError(CE_Failure, CPLE_NotSupported,
3787 : "Unable to export files with zero bands.");
3788 1 : return nullptr;
3789 : }
3790 :
3791 : // Try now to create the file to avoid memory leaks if it is
3792 : // the SDK that fails to do it.
3793 175 : VSILFILE *fpVSIL = VSIFOpenL(pszFilename, "wb");
3794 175 : if (fpVSIL == nullptr)
3795 : {
3796 3 : CPLError(CE_Failure, CPLE_OpenFailed, "Failed to create %s.",
3797 : pszFilename);
3798 3 : return nullptr;
3799 : }
3800 172 : VSIFCloseL(fpVSIL);
3801 172 : VSIUnlink(pszFilename);
3802 :
3803 172 : HDF4ImageDataset *poDS = new HDF4ImageDataset();
3804 :
3805 344 : CPLMutexHolderD(&hHDF4Mutex);
3806 :
3807 : /* -------------------------------------------------------------------- */
3808 : /* Choose rank for the created dataset. */
3809 : /* -------------------------------------------------------------------- */
3810 172 : poDS->iRank = 3;
3811 298 : if (CSLFetchNameValue(papszOptions, "RANK") != nullptr &&
3812 126 : EQUAL(CSLFetchNameValue(papszOptions, "RANK"), "2"))
3813 63 : poDS->iRank = 2;
3814 :
3815 172 : poDS->hSD = SDstart(pszFilename, DFACC_CREATE);
3816 172 : if (poDS->hSD == -1)
3817 : {
3818 0 : CPLError(CE_Failure, CPLE_OpenFailed, "Can't create HDF4 file %s",
3819 : pszFilename);
3820 : // Release mutex otherwise we'll deadlock with GDALDataset own mutex.
3821 0 : CPLReleaseMutex(hHDF4Mutex);
3822 0 : delete poDS;
3823 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3824 0 : return nullptr;
3825 : }
3826 172 : poDS->iXDim = 1;
3827 172 : poDS->iYDim = 0;
3828 172 : poDS->iBandDim = 2;
3829 :
3830 172 : int32 aiDimSizes[H4_MAX_VAR_DIMS] = {};
3831 172 : aiDimSizes[poDS->iXDim] = nXSize;
3832 172 : aiDimSizes[poDS->iYDim] = nYSize;
3833 172 : aiDimSizes[poDS->iBandDim] = nBandsIn;
3834 :
3835 172 : const auto GetHDFType = [](GDALDataType eTypeIn)
3836 : {
3837 172 : switch (eTypeIn)
3838 : {
3839 17 : case GDT_Float64:
3840 17 : return DFNT_FLOAT64;
3841 17 : case GDT_Float32:
3842 17 : return DFNT_FLOAT32;
3843 2 : case GDT_UInt64:
3844 : // SDCreate() doesn't like it
3845 2 : return DFNT_UINT64;
3846 17 : case GDT_UInt32:
3847 17 : return DFNT_UINT32;
3848 17 : case GDT_UInt16:
3849 17 : return DFNT_UINT16;
3850 2 : case GDT_Int64:
3851 : // SDCreate() doesn't like it
3852 2 : return DFNT_INT64;
3853 17 : case GDT_Int32:
3854 17 : return DFNT_INT32;
3855 17 : case GDT_Int16:
3856 17 : return DFNT_INT16;
3857 38 : case GDT_Byte:
3858 38 : return DFNT_UINT8;
3859 16 : case GDT_Int8:
3860 16 : return DFNT_INT8;
3861 12 : default:
3862 12 : CPLError(CE_Warning, CPLE_NotSupported,
3863 : "Datatype %s not supported. Defauting to Byte",
3864 : GDALGetDataTypeName(eTypeIn));
3865 12 : return DFNT_UINT8;
3866 : }
3867 : };
3868 :
3869 172 : const char *pszSDSName = nullptr;
3870 172 : int32 iSDS = -1;
3871 :
3872 172 : if (poDS->iRank == 2)
3873 : {
3874 126 : for (int iBand = 0; iBand < nBandsIn; iBand++)
3875 : {
3876 63 : pszSDSName = CPLSPrintf("Band%d", iBand);
3877 63 : iSDS = SDcreate(poDS->hSD, pszSDSName, GetHDFType(eType),
3878 : poDS->iRank, aiDimSizes);
3879 63 : if (iSDS < 0)
3880 0 : break;
3881 63 : SDendaccess(iSDS);
3882 : }
3883 : }
3884 109 : else if (poDS->iRank == 3)
3885 : {
3886 109 : pszSDSName = "3-dimensional Scientific Dataset";
3887 109 : poDS->iDataset = 0;
3888 109 : iSDS = SDcreate(poDS->hSD, pszSDSName, GetHDFType(eType), poDS->iRank,
3889 : aiDimSizes);
3890 : }
3891 : else
3892 : {
3893 : // Should never happen.
3894 : // Release mutex otherwise we'll deadlock with GDALDataset own mutex.
3895 0 : CPLReleaseMutex(hHDF4Mutex);
3896 0 : delete poDS;
3897 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3898 0 : return nullptr;
3899 : }
3900 :
3901 172 : if (iSDS < 0)
3902 : {
3903 4 : CPLError(CE_Failure, CPLE_AppDefined,
3904 : "Can't create SDS with rank %ld for file %s",
3905 4 : static_cast<long>(poDS->iRank), pszFilename);
3906 : // Release mutex otherwise we'll deadlock with GDALDataset own mutex.
3907 4 : CPLReleaseMutex(hHDF4Mutex);
3908 4 : delete poDS;
3909 4 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3910 4 : return nullptr;
3911 : }
3912 :
3913 168 : poDS->nRasterXSize = nXSize;
3914 168 : poDS->nRasterYSize = nYSize;
3915 168 : poDS->eAccess = GA_Update;
3916 168 : poDS->iDatasetType = HDF4_SDS;
3917 168 : poDS->iSubdatasetType = H4ST_GDAL;
3918 168 : poDS->nBands = nBandsIn;
3919 :
3920 : /* -------------------------------------------------------------------- */
3921 : /* Create band information objects. */
3922 : /* -------------------------------------------------------------------- */
3923 378 : for (int iBand = 1; iBand <= nBandsIn; iBand++)
3924 210 : poDS->SetBand(iBand, new HDF4ImageRasterBand(poDS, iBand, eType));
3925 :
3926 168 : SDsetattr(poDS->hSD, "Signature", DFNT_CHAR8,
3927 : static_cast<int>(strlen(pszGDALSignature)) + 1, pszGDALSignature);
3928 :
3929 168 : return GDALDataset::FromHandle(poDS);
3930 : }
3931 :
3932 : /************************************************************************/
3933 : /* GDALRegister_HDF4Image() */
3934 : /************************************************************************/
3935 :
3936 354 : void GDALRegister_HDF4Image()
3937 :
3938 : {
3939 354 : if (GDALGetDriverByName(HDF4_IMAGE_DRIVER_NAME) != nullptr)
3940 0 : return;
3941 :
3942 354 : GDALDriver *poDriver = new GDALDriver();
3943 354 : HDF4ImageDriverSetCommonMetadata(poDriver);
3944 :
3945 : #ifdef HDF4_HAS_MAXOPENFILES
3946 354 : poDriver->SetMetadataItem("HDF4_HAS_MAXOPENFILES", "YES");
3947 : #endif
3948 354 : poDriver->pfnOpen = HDF4ImageDataset::Open;
3949 354 : poDriver->pfnCreate = HDF4ImageDataset::Create;
3950 :
3951 354 : GetGDALDriverManager()->RegisterDriver(poDriver);
3952 : }
|
