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