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