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