LCOV - code coverage report
Current view: top level - frmts/hdf4 - hdf4imagedataset.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 508 1704 29.8 %
Date: 2025-09-10 17:48:50 Functions: 21 36 58.3 %

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

Generated by: LCOV version 1.14