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