LCOV - code coverage report
Current view: top level - frmts/hdf4 - hdf4imagedataset.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 511 1707 29.9 %
Date: 2025-08-01 13:10:27 Functions: 23 38 60.5 %

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

Generated by: LCOV version 1.14