LCOV - code coverage report
Current view: top level - frmts/grib - gribdataset.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 1222 1431 85.4 %
Date: 2025-09-10 17:48:50 Functions: 60 64 93.8 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GRIB Driver
       4             :  * Purpose:  GDALDataset driver for GRIB translator for read support
       5             :  * Author:   Bas Retsios, retsios@itc.nl
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2007, ITC
       9             :  * Copyright (c) 2008-2017, Even Rouault <even dot rouault at spatialys dot com>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ******************************************************************************
      13             :  *
      14             :  */
      15             : 
      16             : #include "cpl_port.h"
      17             : #include "gribdataset.h"
      18             : #include "gribdrivercore.h"
      19             : 
      20             : #include <cassert>
      21             : #include <cerrno>
      22             : #include <cmath>
      23             : #include <cstddef>
      24             : #include <cstdio>
      25             : #include <cstdlib>
      26             : #include <cstring>
      27             : 
      28             : #include <algorithm>
      29             : #include <mutex>
      30             : #include <set>
      31             : #include <string>
      32             : #include <vector>
      33             : 
      34             : #include "cpl_conv.h"
      35             : #include "cpl_error.h"
      36             : #include "cpl_multiproc.h"
      37             : #include "cpl_string.h"
      38             : #include "cpl_vsi.h"
      39             : #include "cpl_time.h"
      40             : #include "degrib/degrib/degrib2.h"
      41             : #include "degrib/degrib/inventory.h"
      42             : #include "degrib/degrib/meta.h"
      43             : #include "degrib/degrib/metaname.h"
      44             : #include "degrib/degrib/myerror.h"
      45             : #include "degrib/degrib/type.h"
      46             : CPL_C_START
      47             : #include "degrib/g2clib/grib2.h"
      48             : #include "degrib/g2clib/pdstemplates.h"
      49             : CPL_C_END
      50             : #include "gdal.h"
      51             : #include "gdal_frmts.h"
      52             : #include "gdal_pam.h"
      53             : #include "gdal_priv.h"
      54             : #include "ogr_spatialref.h"
      55             : #include "memdataset.h"
      56             : 
      57             : static CPLMutex *hGRIBMutex = nullptr;
      58             : 
      59             : /************************************************************************/
      60             : /*                         ConvertUnitInText()                          */
      61             : /************************************************************************/
      62             : 
      63         922 : static CPLString ConvertUnitInText(bool bMetricUnits, const char *pszTxt)
      64             : {
      65         922 :     if (pszTxt == nullptr)
      66           0 :         return CPLString();
      67         922 :     if (!bMetricUnits)
      68           4 :         return pszTxt;
      69             : 
      70        1836 :     CPLString osRes(pszTxt);
      71         918 :     size_t iPos = osRes.find("[K]");
      72         918 :     if (iPos != std::string::npos)
      73          96 :         osRes = osRes.substr(0, iPos) + "[C]" + osRes.substr(iPos + 3);
      74         918 :     return osRes;
      75             : }
      76             : 
      77             : /************************************************************************/
      78             : /*                         Lon360to180()                               */
      79             : /************************************************************************/
      80             : 
      81         442 : static inline double Lon360to180(double lon)
      82             : {
      83         442 :     if (lon == 180)
      84           0 :         return 180;
      85         442 :     return fmod(lon + 180, 360) - 180;
      86             : }
      87             : 
      88             : /************************************************************************/
      89             : /*                           GRIBRasterBand()                            */
      90             : /************************************************************************/
      91             : 
      92         444 : GRIBRasterBand::GRIBRasterBand(GRIBDataset *poDSIn, int nBandIn,
      93         444 :                                inventoryType *psInv)
      94         444 :     : start(psInv->start), subgNum(psInv->subgNum),
      95         888 :       longFstLevel(CPLStrdup(psInv->longFstLevel)), m_Grib_Data(nullptr),
      96         444 :       m_Grib_MetaData(nullptr), nGribDataXSize(poDSIn->nRasterXSize),
      97         444 :       nGribDataYSize(poDSIn->nRasterYSize), m_nGribVersion(psInv->GribVersion),
      98         444 :       m_bHasLookedForNoData(false), m_dfNoData(0.0), m_bHasNoData(false)
      99             : 
     100             : {
     101         444 :     poDS = poDSIn;
     102         444 :     nBand = nBandIn;
     103             : 
     104             :     // Let user do -ot Float32 if needed for saving space, GRIB contains
     105             :     // Float64 (though not fully utilized most of the time).
     106         444 :     eDataType = GDT_Float64;
     107             : 
     108         444 :     nBlockXSize = poDSIn->nRasterXSize;
     109         444 :     nBlockYSize = 1;
     110             : 
     111         444 :     if (psInv->unitName != nullptr && psInv->comment != nullptr &&
     112         412 :         psInv->element != nullptr)
     113             :     {
     114         412 :         bLoadedMetadata = true;
     115             :         const char *pszGribNormalizeUnits =
     116         412 :             CPLGetConfigOption("GRIB_NORMALIZE_UNITS", "YES");
     117         412 :         bool bMetricUnits = CPLTestBool(pszGribNormalizeUnits);
     118             : 
     119         412 :         SetMetadataItem("GRIB_UNIT",
     120         824 :                         ConvertUnitInText(bMetricUnits, psInv->unitName));
     121         412 :         SetMetadataItem("GRIB_COMMENT",
     122         824 :                         ConvertUnitInText(bMetricUnits, psInv->comment));
     123         412 :         SetMetadataItem("GRIB_ELEMENT", psInv->element);
     124         412 :         SetMetadataItem("GRIB_SHORT_NAME", psInv->shortFstLevel);
     125         412 :         SetMetadataItem("GRIB_REF_TIME",
     126         824 :                         CPLString().Printf("%.0f", psInv->refTime));
     127         412 :         SetMetadataItem("GRIB_VALID_TIME",
     128         824 :                         CPLString().Printf("%.0f", psInv->validTime));
     129         412 :         SetMetadataItem("GRIB_FORECAST_SECONDS",
     130         824 :                         CPLString().Printf("%.0f", psInv->foreSec));
     131             :     }
     132         444 : }
     133             : 
     134             : /************************************************************************/
     135             : /*                           FindMetaData()                             */
     136             : /************************************************************************/
     137             : 
     138         775 : void GRIBRasterBand::FindMetaData()
     139             : {
     140         775 :     if (bLoadedMetadata)
     141         749 :         return;
     142          26 :     if (m_Grib_MetaData == nullptr)
     143             :     {
     144             :         grib_MetaData *metaData;
     145           5 :         GRIBDataset *poGDS = static_cast<GRIBDataset *>(poDS);
     146           5 :         GRIBRasterBand::ReadGribData(poGDS->fp, start, subgNum, nullptr,
     147             :                                      &metaData);
     148           5 :         if (metaData == nullptr)
     149           0 :             return;
     150           5 :         m_Grib_MetaData = metaData;
     151             :     }
     152          26 :     bLoadedMetadata = true;
     153          26 :     m_nGribVersion = m_Grib_MetaData->GribVersion;
     154             : 
     155             :     const char *pszGribNormalizeUnits =
     156          26 :         CPLGetConfigOption("GRIB_NORMALIZE_UNITS", "YES");
     157          26 :     bool bMetricUnits = CPLTestBool(pszGribNormalizeUnits);
     158             : 
     159          26 :     GDALRasterBand::SetMetadataItem(
     160             :         "GRIB_UNIT",
     161          52 :         ConvertUnitInText(bMetricUnits, m_Grib_MetaData->unitName));
     162          26 :     GDALRasterBand::SetMetadataItem(
     163             :         "GRIB_COMMENT",
     164          52 :         ConvertUnitInText(bMetricUnits, m_Grib_MetaData->comment));
     165             : 
     166          26 :     GDALRasterBand::SetMetadataItem("GRIB_ELEMENT", m_Grib_MetaData->element);
     167          26 :     GDALRasterBand::SetMetadataItem("GRIB_SHORT_NAME",
     168          26 :                                     m_Grib_MetaData->shortFstLevel);
     169             : 
     170          26 :     if (m_nGribVersion == 2)
     171             :     {
     172          14 :         GDALRasterBand::SetMetadataItem(
     173             :             "GRIB_REF_TIME",
     174          28 :             CPLString().Printf("%.0f", m_Grib_MetaData->pds2.refTime));
     175          14 :         GDALRasterBand::SetMetadataItem(
     176             :             "GRIB_VALID_TIME",
     177          28 :             CPLString().Printf("%.0f", m_Grib_MetaData->pds2.sect4.validTime));
     178             :     }
     179          12 :     else if (m_nGribVersion == 1)
     180             :     {
     181          12 :         GDALRasterBand::SetMetadataItem(
     182             :             "GRIB_REF_TIME",
     183          24 :             CPLString().Printf("%.0f", m_Grib_MetaData->pds1.refTime));
     184          12 :         GDALRasterBand::SetMetadataItem(
     185             :             "GRIB_VALID_TIME",
     186          24 :             CPLString().Printf("%.0f", m_Grib_MetaData->pds1.validTime));
     187             :     }
     188             : 
     189          26 :     GDALRasterBand::SetMetadataItem(
     190             :         "GRIB_FORECAST_SECONDS",
     191          52 :         CPLString().Printf("%d", m_Grib_MetaData->deltTime));
     192             : }
     193             : 
     194             : /************************************************************************/
     195             : /*                          FindTrueStart()                             */
     196             : /*                                                                      */
     197             : /*      Scan after the official start of the message to find its        */
     198             : /*      true starting offset.                                           */
     199             : /************************************************************************/
     200         866 : vsi_l_offset GRIBRasterBand::FindTrueStart(VSILFILE *fp, vsi_l_offset start)
     201             : {
     202             :     // GRIB messages can be preceded by "garbage". GRIB2Inventory()
     203             :     // does not return the offset to the real start of the message
     204             :     char szHeader[1024 + 1];
     205         866 :     VSIFSeekL(fp, start, SEEK_SET);
     206             :     const int nRead =
     207         866 :         static_cast<int>(VSIFReadL(szHeader, 1, sizeof(szHeader) - 1, fp));
     208         866 :     szHeader[nRead] = 0;
     209             :     // Find the real offset of the fist message
     210         866 :     int nOffsetFirstMessage = 0;
     211        2548 :     for (int j = 0; j + 3 < nRead; j++)
     212             :     {
     213        2548 :         if (STARTS_WITH_CI(szHeader + j, "GRIB")
     214             : #ifdef ENABLE_TDLP
     215             :             || STARTS_WITH_CI(szHeader + j, "TDLP")
     216             : #endif
     217             :         )
     218             :         {
     219         866 :             nOffsetFirstMessage = j;
     220         866 :             break;
     221             :         }
     222             :     }
     223         866 :     return start + nOffsetFirstMessage;
     224             : }
     225             : 
     226             : /************************************************************************/
     227             : /*                      FindPDSTemplateGRIB2()                          */
     228             : /*                                                                      */
     229             : /*      Scan the file for the PDS template info and represent it as     */
     230             : /*      metadata.                                                       */
     231             : /************************************************************************/
     232             : 
     233         886 : void GRIBRasterBand::FindPDSTemplateGRIB2()
     234             : 
     235             : {
     236         886 :     CPLAssert(m_nGribVersion == 2);
     237             : 
     238         886 :     if (bLoadedPDS)
     239         558 :         return;
     240         328 :     bLoadedPDS = true;
     241             : 
     242         328 :     GRIBDataset *poGDS = static_cast<GRIBDataset *>(poDS);
     243         328 :     start = FindTrueStart(poGDS->fp, start);
     244             : 
     245             :     // Read section 0
     246             :     GByte abySection0[16];
     247         656 :     if (VSIFSeekL(poGDS->fp, start, SEEK_SET) != 0 ||
     248         328 :         VSIFReadL(abySection0, 16, 1, poGDS->fp) != 1)
     249             :     {
     250           0 :         CPLDebug("GRIB", "Cannot read leading bytes of section 0");
     251           0 :         return;
     252             :     }
     253         328 :     GByte nDiscipline = abySection0[7 - 1];
     254         328 :     CPLString osDiscipline;
     255         328 :     osDiscipline = CPLString().Printf("%d", nDiscipline);
     256             :     static const char *const table00[] = {"Meteorological",
     257             :                                           "Hydrological",
     258             :                                           "Land Surface",
     259             :                                           "Space products",
     260             :                                           "Space products",
     261             :                                           "Reserved",
     262             :                                           "Reserved",
     263             :                                           "Reserved",
     264             :                                           "Reserved",
     265             :                                           "Reserved",
     266             :                                           "Oceanographic Products"};
     267         328 :     m_nDisciplineCode = nDiscipline;
     268         328 :     if (nDiscipline < CPL_ARRAYSIZE(table00))
     269             :     {
     270         327 :         m_osDisciplineName = table00[nDiscipline];
     271         654 :         osDiscipline += CPLString("(") +
     272        1308 :                         CPLString(table00[nDiscipline]).replaceAll(' ', '_') +
     273         327 :                         ")";
     274             :     }
     275             : 
     276         328 :     GDALRasterBand::SetMetadataItem("GRIB_DISCIPLINE", osDiscipline.c_str());
     277             : 
     278         328 :     GByte abyHead[5] = {0};
     279             : 
     280         328 :     if (VSIFReadL(abyHead, 5, 1, poGDS->fp) != 1)
     281             :     {
     282           0 :         CPLDebug("GRIB", "Cannot read 5 leading bytes past section 0");
     283           0 :         return;
     284             :     }
     285             : 
     286         328 :     GUInt32 nSectSize = 0;
     287         328 :     if (abyHead[4] == 1)
     288             :     {
     289         328 :         memcpy(&nSectSize, abyHead, 4);
     290         328 :         CPL_MSBPTR32(&nSectSize);
     291         328 :         if (nSectSize >= 21 && nSectSize <= 100000 /* arbitrary upper limit */)
     292             :         {
     293         328 :             GByte *pabyBody = static_cast<GByte *>(CPLMalloc(nSectSize));
     294         328 :             memcpy(pabyBody, abyHead, 5);
     295         328 :             VSIFReadL(pabyBody + 5, 1, nSectSize - 5, poGDS->fp);
     296             : 
     297         656 :             CPLString osIDS;
     298         328 :             unsigned short nCenter = static_cast<unsigned short>(
     299         328 :                 pabyBody[6 - 1] * 256 + pabyBody[7 - 1]);
     300         328 :             if (nCenter != GRIB2MISSING_u1 && nCenter != GRIB2MISSING_u2)
     301             :             {
     302          93 :                 osIDS += "CENTER=";
     303          93 :                 m_nCenter = nCenter;
     304          93 :                 osIDS += CPLSPrintf("%d", nCenter);
     305          93 :                 const char *pszCenter = centerLookup(nCenter);
     306          93 :                 if (pszCenter)
     307             :                 {
     308          93 :                     m_osCenterName = pszCenter;
     309          93 :                     osIDS += CPLString("(") + pszCenter + ")";
     310             :                 }
     311             :             }
     312             : 
     313         328 :             unsigned short nSubCenter = static_cast<unsigned short>(
     314         328 :                 pabyBody[8 - 1] * 256 + pabyBody[9 - 1]);
     315         328 :             if (nSubCenter != GRIB2MISSING_u2)
     316             :             {
     317         119 :                 if (!osIDS.empty())
     318          78 :                     osIDS += " ";
     319         119 :                 osIDS += "SUBCENTER=";
     320         119 :                 osIDS += CPLSPrintf("%d", nSubCenter);
     321         119 :                 m_nSubCenter = nSubCenter;
     322         119 :                 const char *pszSubCenter = subCenterLookup(nCenter, nSubCenter);
     323         119 :                 if (pszSubCenter)
     324             :                 {
     325           3 :                     m_osSubCenterName = pszSubCenter;
     326           3 :                     osIDS += CPLString("(") + pszSubCenter + ")";
     327             :                 }
     328             :             }
     329             : 
     330         328 :             if (!osIDS.empty())
     331         134 :                 osIDS += " ";
     332         328 :             osIDS += "MASTER_TABLE=";
     333         328 :             osIDS += CPLSPrintf("%d", pabyBody[10 - 1]);
     334         328 :             osIDS += " ";
     335         328 :             osIDS += "LOCAL_TABLE=";
     336         328 :             osIDS += CPLSPrintf("%d", pabyBody[11 - 1]);
     337         328 :             osIDS += " ";
     338         328 :             osIDS += "SIGNF_REF_TIME=";
     339         328 :             unsigned nSignRefTime = pabyBody[12 - 1];
     340         328 :             osIDS += CPLSPrintf("%d", nSignRefTime);
     341             :             static const char *const table12[] = {
     342             :                 "Analysis", "Start of Forecast", "Verifying time of forecast",
     343             :                 "Observation time"};
     344         328 :             if (nSignRefTime < CPL_ARRAYSIZE(table12))
     345             :             {
     346         328 :                 m_osSignRefTimeName = table12[nSignRefTime];
     347         656 :                 osIDS += CPLString("(") +
     348        1312 :                          CPLString(table12[nSignRefTime]).replaceAll(' ', '_') +
     349         328 :                          ")";
     350             :             }
     351         328 :             osIDS += " ";
     352         328 :             osIDS += "REF_TIME=";
     353             :             m_osRefTime =
     354             :                 CPLSPrintf("%04d-%02d-%02dT%02d:%02d:%02dZ",
     355         328 :                            pabyBody[13 - 1] * 256 + pabyBody[14 - 1],
     356         328 :                            pabyBody[15 - 1], pabyBody[16 - 1], pabyBody[17 - 1],
     357         328 :                            pabyBody[18 - 1], pabyBody[19 - 1]);
     358         328 :             osIDS += m_osRefTime;
     359         328 :             osIDS += " ";
     360         328 :             osIDS += "PROD_STATUS=";
     361         328 :             unsigned nProdStatus = pabyBody[20 - 1];
     362         328 :             osIDS += CPLSPrintf("%d", nProdStatus);
     363             :             static const char *const table13[] = {
     364             :                 "Operational",     "Operational test",
     365             :                 "Research",        "Re-analysis",
     366             :                 "TIGGE",           "TIGGE test",
     367             :                 "S2S operational", "S2S test",
     368             :                 "UERRA",           "UERRA test"};
     369         328 :             if (nProdStatus < CPL_ARRAYSIZE(table13))
     370             :             {
     371         133 :                 m_osProductionStatus = table13[nProdStatus];
     372         266 :                 osIDS += CPLString("(") +
     373         532 :                          CPLString(table13[nProdStatus]).replaceAll(' ', '_') +
     374         133 :                          ")";
     375             :             }
     376         328 :             osIDS += " ";
     377         328 :             osIDS += "TYPE=";
     378         328 :             unsigned nType = pabyBody[21 - 1];
     379         328 :             osIDS += CPLSPrintf("%d", nType);
     380             :             static const char *const table14[] = {
     381             :                 "Analysis",
     382             :                 "Forecast",
     383             :                 "Analysis and forecast",
     384             :                 "Control forecast",
     385             :                 "Perturbed forecast",
     386             :                 "Control and perturbed forecast",
     387             :                 "Processed satellite observations",
     388             :                 "Processed radar observations",
     389             :                 "Event Probability"};
     390         328 :             if (nType < CPL_ARRAYSIZE(table14))
     391             :             {
     392         133 :                 m_osType = table14[nType];
     393         266 :                 osIDS += CPLString("(") +
     394         399 :                          CPLString(table14[nType]).replaceAll(' ', '_') + ")";
     395             :             }
     396             : 
     397         328 :             GDALRasterBand::SetMetadataItem("GRIB_IDS", osIDS);
     398             : 
     399         328 :             CPLFree(pabyBody);
     400             :         }
     401             : 
     402         328 :         if (VSIFReadL(abyHead, 5, 1, poGDS->fp) != 1)
     403             :         {
     404           0 :             CPLDebug("GRIB", "Cannot read 5 leading bytes past section 1");
     405           0 :             return;
     406             :         }
     407             :     }
     408             : 
     409         328 :     if (subgNum > 0)
     410             :     {
     411             :         // If we are a subgrid, then iterate over all preceding subgrids
     412          14 :         for (int iSubMessage = 0; iSubMessage < subgNum;)
     413             :         {
     414          12 :             memcpy(&nSectSize, abyHead, 4);
     415          12 :             CPL_MSBPTR32(&nSectSize);
     416          12 :             if (nSectSize < 5)
     417             :             {
     418           0 :                 CPLDebug("GRIB", "Invalid section size for iSubMessage = %d",
     419             :                          iSubMessage);
     420           0 :                 return;
     421             :             }
     422          12 :             if (VSIFSeekL(poGDS->fp, nSectSize - 5, SEEK_CUR) != 0)
     423             :             {
     424           0 :                 CPLDebug("GRIB",
     425             :                          "Cannot read past section for iSubMessage = %d",
     426             :                          iSubMessage);
     427           0 :                 return;
     428             :             }
     429          12 :             if (abyHead[4] < 2 || abyHead[4] > 7)
     430             :             {
     431           0 :                 CPLDebug("GRIB", "Invalid section number for iSubMessage = %d",
     432             :                          iSubMessage);
     433           0 :                 return;
     434             :             }
     435          12 :             if (abyHead[4] == 7)
     436             :             {
     437           2 :                 ++iSubMessage;
     438             :             }
     439          12 :             if (VSIFReadL(abyHead, 5, 1, poGDS->fp) != 1)
     440             :             {
     441           0 :                 CPLDebug("GRIB",
     442             :                          "Cannot read 5 leading bytes for iSubMessage = %d",
     443             :                          iSubMessage);
     444           0 :                 return;
     445             :             }
     446             :         }
     447             :     }
     448             : 
     449             :     // Skip to section 4
     450         957 :     while (abyHead[4] != 4)
     451             :     {
     452         629 :         memcpy(&nSectSize, abyHead, 4);
     453         629 :         CPL_MSBPTR32(&nSectSize);
     454             : 
     455         629 :         const int nCurSection = abyHead[4];
     456         629 :         if (nSectSize < 5)
     457             :         {
     458           0 :             CPLDebug("GRIB", "Invalid section size for section %d",
     459             :                      nCurSection);
     460           0 :             return;
     461             :         }
     462        1258 :         if (VSIFSeekL(poGDS->fp, nSectSize - 5, SEEK_CUR) != 0 ||
     463         629 :             VSIFReadL(abyHead, 5, 1, poGDS->fp) != 1)
     464             :         {
     465           0 :             CPLDebug("GRIB", "Cannot read section %d", nCurSection);
     466           0 :             return;
     467             :         }
     468             :     }
     469             : 
     470             :     // Collect section 4 octet information.  We read the file
     471             :     // ourselves since the GRIB API does not appear to preserve all
     472             :     // this for us.
     473             :     // if( abyHead[4] == 4 )
     474             :     {
     475         328 :         memcpy(&nSectSize, abyHead, 4);
     476         328 :         CPL_MSBPTR32(&nSectSize);
     477         328 :         if (nSectSize >= 9 && nSectSize <= 100000 /* arbitrary upper limit */)
     478             :         {
     479         328 :             GByte *pabyBody = static_cast<GByte *>(CPLMalloc(nSectSize));
     480         328 :             memcpy(pabyBody, abyHead, 5);
     481         328 :             if (VSIFReadL(pabyBody + 5, 1, nSectSize - 5, poGDS->fp) !=
     482         328 :                 nSectSize - 5)
     483             :             {
     484           0 :                 CPLDebug("GRIB", "Cannot read section 4");
     485           0 :                 CPLFree(pabyBody);
     486           0 :                 return;
     487             :             }
     488             : 
     489         328 :             GUInt16 nCoordCount = 0;
     490         328 :             memcpy(&nCoordCount, pabyBody + 6 - 1, 2);
     491         328 :             CPL_MSBPTR16(&nCoordCount);
     492             : 
     493         328 :             GUInt16 nPDTN = 0;
     494         328 :             memcpy(&nPDTN, pabyBody + 8 - 1, 2);
     495         328 :             CPL_MSBPTR16(&nPDTN);
     496             : 
     497         328 :             GDALRasterBand::SetMetadataItem("GRIB_PDS_PDTN",
     498         656 :                                             CPLString().Printf("%d", nPDTN));
     499         328 :             m_nPDTN = nPDTN;
     500             : 
     501         656 :             CPLString osOctet;
     502         328 :             const int nTemplateFoundByteCount =
     503         328 :                 nSectSize - 9U >= nCoordCount * 4U
     504         328 :                     ? static_cast<int>(nSectSize - 9 - nCoordCount * 4)
     505             :                     : 0;
     506        9274 :             for (int i = 0; i < nTemplateFoundByteCount; i++)
     507             :             {
     508        8946 :                 char szByte[10] = {'\0'};
     509             : 
     510        8946 :                 if (i == 0)
     511         328 :                     snprintf(szByte, sizeof(szByte), "%d", pabyBody[i + 9]);
     512             :                 else
     513        8618 :                     snprintf(szByte, sizeof(szByte), " %d", pabyBody[i + 9]);
     514        8946 :                 osOctet += szByte;
     515             :             }
     516             : 
     517         328 :             GDALRasterBand::SetMetadataItem("GRIB_PDS_TEMPLATE_NUMBERS",
     518             :                                             osOctet);
     519             : 
     520         328 :             g2int iofst = 0;
     521         328 :             g2int pdsnum = 0;
     522         328 :             g2int *pdstempl = nullptr;
     523         328 :             g2int mappdslen = 0;
     524         328 :             g2float *coordlist = nullptr;
     525         328 :             g2int numcoord = 0;
     526         328 :             if (getpdsindex(nPDTN) < 0)
     527             :             {
     528           3 :                 CPLError(CE_Warning, CPLE_NotSupported,
     529             :                          "Template 4.%d is not recognized currently", nPDTN);
     530             :             }
     531         325 :             else if (g2_unpack4(pabyBody, nSectSize, &iofst, &pdsnum, &pdstempl,
     532         325 :                                 &mappdslen, &coordlist, &numcoord) == 0)
     533             :             {
     534         325 :                 gtemplate *mappds = extpdstemplate(pdsnum, pdstempl);
     535         325 :                 if (mappds)
     536             :                 {
     537         325 :                     int nTemplateByteCount = 0;
     538        5592 :                     for (int i = 0; i < mappds->maplen; i++)
     539        5267 :                         nTemplateByteCount += abs(mappds->map[i]);
     540         390 :                     for (int i = 0; i < mappds->extlen; i++)
     541          65 :                         nTemplateByteCount += abs(mappds->ext[i]);
     542         325 :                     if (nTemplateByteCount == nTemplateFoundByteCount)
     543             :                     {
     544         646 :                         CPLString osValues;
     545        5623 :                         for (g2int i = 0; i < mappds->maplen + mappds->extlen;
     546             :                              i++)
     547             :                         {
     548        5300 :                             if (i > 0)
     549        4977 :                                 osValues += " ";
     550        5300 :                             const int nEltSize =
     551        5300 :                                 (i < mappds->maplen)
     552        5300 :                                     ? mappds->map[i]
     553          65 :                                     : mappds->ext[i - mappds->maplen];
     554        5300 :                             if (nEltSize == 4)
     555             :                             {
     556          88 :                                 m_anPDSTemplateAssembledValues.push_back(
     557          88 :                                     static_cast<GUInt32>(pdstempl[i]));
     558             :                                 osValues += CPLSPrintf(
     559          88 :                                     "%u", static_cast<GUInt32>(pdstempl[i]));
     560             :                             }
     561             :                             else
     562             :                             {
     563        5212 :                                 m_anPDSTemplateAssembledValues.push_back(
     564        5212 :                                     pdstempl[i]);
     565        5212 :                                 osValues += CPLSPrintf("%d", pdstempl[i]);
     566             :                             }
     567             :                         }
     568         323 :                         GDALRasterBand::SetMetadataItem(
     569             :                             "GRIB_PDS_TEMPLATE_ASSEMBLED_VALUES", osValues);
     570             :                     }
     571             :                     else
     572             :                     {
     573           2 :                         CPLDebug(
     574             :                             "GRIB",
     575             :                             "Cannot expose GRIB_PDS_TEMPLATE_ASSEMBLED_VALUES "
     576             :                             "as we would expect %d bytes from the "
     577             :                             "tables, but %d are available",
     578             :                             nTemplateByteCount, nTemplateFoundByteCount);
     579             :                     }
     580             : 
     581         325 :                     free(mappds->ext);
     582         325 :                     free(mappds);
     583             :                 }
     584             :             }
     585         328 :             free(pdstempl);
     586         328 :             free(coordlist);
     587             : 
     588         328 :             CPLFree(pabyBody);
     589             : 
     590         656 :             FindNoDataGrib2(false);
     591             :         }
     592             :         else
     593             :         {
     594           0 :             CPLDebug("GRIB", "Invalid section size for section %d", 4);
     595             :         }
     596             :     }
     597             : }
     598             : 
     599             : /************************************************************************/
     600             : /*                        FindNoDataGrib2()                             */
     601             : /************************************************************************/
     602             : 
     603         328 : void GRIBRasterBand::FindNoDataGrib2(bool bSeekToStart)
     604             : {
     605             :     // There is no easy way in the degrib API to retrieve the nodata value
     606             :     // without decompressing the data point section (which is slow), so
     607             :     // retrieve nodata value by parsing section 5 (Data Representation Section)
     608             :     // We also check section 6 to see if there is a bitmap
     609         328 :     GRIBDataset *poGDS = static_cast<GRIBDataset *>(poDS);
     610         328 :     CPLAssert(m_nGribVersion == 2);
     611             : 
     612         328 :     if (m_bHasLookedForNoData)
     613           0 :         return;
     614         328 :     m_bHasLookedForNoData = true;
     615             : 
     616         328 :     if (bSeekToStart)
     617             :     {
     618             :         // Skip over section 0
     619           0 :         VSIFSeekL(poGDS->fp, start + 16, SEEK_SET);
     620             :     }
     621             : 
     622         328 :     GByte abyHead[5] = {0};
     623         328 :     VSIFReadL(abyHead, 5, 1, poGDS->fp);
     624             : 
     625             :     // Skip to section 5
     626         328 :     GUInt32 nSectSize = 0;
     627         328 :     while (abyHead[4] != 5)
     628             :     {
     629           0 :         memcpy(&nSectSize, abyHead, 4);
     630           0 :         CPL_MSBPTR32(&nSectSize);
     631             : 
     632           0 :         if (nSectSize < 5 ||
     633           0 :             VSIFSeekL(poGDS->fp, nSectSize - 5, SEEK_CUR) != 0 ||
     634           0 :             VSIFReadL(abyHead, 5, 1, poGDS->fp) != 1)
     635           0 :             break;
     636             :     }
     637             : 
     638             :     // See http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_sect5.shtml
     639         328 :     if (abyHead[4] == 5)
     640             :     {
     641         328 :         memcpy(&nSectSize, abyHead, 4);
     642         328 :         CPL_MSBPTR32(&nSectSize);
     643         328 :         if (nSectSize >= 11 && nSectSize <= 100000 /* arbitrary upper limit */)
     644             :         {
     645         328 :             GByte *pabyBody = static_cast<GByte *>(CPLMalloc(nSectSize));
     646         328 :             memcpy(pabyBody, abyHead, 5);
     647         328 :             VSIFReadL(pabyBody + 5, 1, nSectSize - 5, poGDS->fp);
     648             : 
     649         328 :             GUInt16 nDRTN = 0;
     650         328 :             memcpy(&nDRTN, pabyBody + 10 - 1, 2);
     651         328 :             CPL_MSBPTR16(&nDRTN);
     652             : 
     653         328 :             GDALRasterBand::SetMetadataItem("DRS_DRTN", CPLSPrintf("%d", nDRTN),
     654             :                                             "GRIB");
     655         328 :             if ((nDRTN == GS5_SIMPLE || nDRTN == GS5_CMPLX ||
     656         182 :                  nDRTN == GS5_CMPLXSEC || nDRTN == GS5_JPEG2000 ||
     657         126 :                  nDRTN == GS5_PNG) &&
     658         232 :                 nSectSize >= 20)
     659             :             {
     660             :                 float fRef;
     661         232 :                 memcpy(&fRef, pabyBody + 12 - 1, 4);
     662         232 :                 CPL_MSBPTR32(&fRef);
     663         232 :                 GDALRasterBand::SetMetadataItem(
     664             :                     "DRS_REF_VALUE", CPLSPrintf("%.10f", fRef), "GRIB");
     665             : 
     666             :                 GUInt16 nBinaryScaleFactorUnsigned;
     667         232 :                 memcpy(&nBinaryScaleFactorUnsigned, pabyBody + 16 - 1, 2);
     668         232 :                 CPL_MSBPTR16(&nBinaryScaleFactorUnsigned);
     669         232 :                 const int nBSF =
     670             :                     (nBinaryScaleFactorUnsigned & 0x8000)
     671          25 :                         ? -static_cast<int>(nBinaryScaleFactorUnsigned & 0x7FFF)
     672         232 :                         : static_cast<int>(nBinaryScaleFactorUnsigned);
     673         232 :                 GDALRasterBand::SetMetadataItem("DRS_BINARY_SCALE_FACTOR",
     674             :                                                 CPLSPrintf("%d", nBSF), "GRIB");
     675             : 
     676             :                 GUInt16 nDecimalScaleFactorUnsigned;
     677         232 :                 memcpy(&nDecimalScaleFactorUnsigned, pabyBody + 18 - 1, 2);
     678         232 :                 CPL_MSBPTR16(&nDecimalScaleFactorUnsigned);
     679         232 :                 const int nDSF =
     680             :                     (nDecimalScaleFactorUnsigned & 0x8000)
     681          19 :                         ? -static_cast<int>(nDecimalScaleFactorUnsigned &
     682             :                                             0x7FFF)
     683         232 :                         : static_cast<int>(nDecimalScaleFactorUnsigned);
     684         232 :                 GDALRasterBand::SetMetadataItem("DRS_DECIMAL_SCALE_FACTOR",
     685             :                                                 CPLSPrintf("%d", nDSF), "GRIB");
     686             : 
     687         232 :                 const int nBits = pabyBody[20 - 1];
     688         232 :                 GDALRasterBand::SetMetadataItem(
     689             :                     "DRS_NBITS", CPLSPrintf("%d", nBits), "GRIB");
     690             :             }
     691             : 
     692             :             // 2 = Grid Point Data - Complex Packing
     693             :             // 3 = Grid Point Data - Complex Packing and Spatial Differencing
     694         328 :             if ((nDRTN == GS5_CMPLX || nDRTN == GS5_CMPLXSEC) &&
     695          54 :                 nSectSize >= 31)
     696             :             {
     697          54 :                 const int nMiss = pabyBody[23 - 1];
     698          54 :                 if (nMiss == 1 || nMiss == 2)
     699             :                 {
     700          35 :                     const int original_field_type = pabyBody[21 - 1];
     701          35 :                     if (original_field_type == 0)  // Floating Point
     702             :                     {
     703             :                         float fTemp;
     704          27 :                         memcpy(&fTemp, &pabyBody[24 - 1], 4);
     705          27 :                         CPL_MSBPTR32(&fTemp);
     706          27 :                         m_dfNoData = fTemp;
     707          27 :                         m_bHasNoData = true;
     708          27 :                         if (nMiss == 2)
     709             :                         {
     710           0 :                             memcpy(&fTemp, &pabyBody[28 - 1], 4);
     711           0 :                             CPL_MSBPTR32(&fTemp);
     712           0 :                             CPLDebug("GRIB",
     713             :                                      "Secondary missing value also set for "
     714             :                                      "band %d : %f",
     715             :                                      nBand, fTemp);
     716             :                         }
     717             :                     }
     718           8 :                     else if (original_field_type == 1)  // Integer
     719             :                     {
     720             :                         int iTemp;
     721           8 :                         memcpy(&iTemp, &pabyBody[24 - 1], 4);
     722           8 :                         CPL_MSBPTR32(&iTemp);
     723           8 :                         m_dfNoData = iTemp;
     724           8 :                         m_bHasNoData = true;
     725           8 :                         if (nMiss == 2)
     726             :                         {
     727           0 :                             memcpy(&iTemp, &pabyBody[28 - 1], 4);
     728           0 :                             CPL_MSBPTR32(&iTemp);
     729           0 :                             CPLDebug("GRIB",
     730             :                                      "Secondary missing value also set for "
     731             :                                      "band %d : %d",
     732             :                                      nBand, iTemp);
     733             :                         }
     734             :                     }
     735             :                     else
     736             :                     {
     737             :                         // FIXME What to do? Blindly convert to float?
     738           0 :                         CPLDebug("GRIB",
     739             :                                  "Complex Packing - Type of Original Field "
     740             :                                  "Values for band %d:  %u",
     741             :                                  nBand, original_field_type);
     742             :                     }
     743             :                 }
     744             :             }
     745             : 
     746         328 :             if (nDRTN == GS5_CMPLXSEC && nSectSize >= 48)
     747             :             {
     748          21 :                 const int nOrder = pabyBody[48 - 1];
     749          21 :                 GDALRasterBand::SetMetadataItem(
     750             :                     "DRS_SPATIAL_DIFFERENCING_ORDER", CPLSPrintf("%d", nOrder),
     751             :                     "GRIB");
     752             :             }
     753             : 
     754         328 :             CPLFree(pabyBody);
     755             :         }
     756           0 :         else if (nSectSize > 5)
     757             :         {
     758           0 :             VSIFSeekL(poGDS->fp, nSectSize - 5, SEEK_CUR);
     759             :         }
     760             :     }
     761             : 
     762         328 :     if (!m_bHasNoData)
     763             :     {
     764             :         // Check bitmap section
     765         293 :         GByte abySection6[6] = {0};
     766         293 :         VSIFReadL(abySection6, 6, 1, poGDS->fp);
     767             :         // Is there a bitmap ?
     768         293 :         if (abySection6[4] == 6 && abySection6[5] == 0)
     769             :         {
     770           6 :             m_dfNoData = 9999.0;  // Same value as in metaparse.cpp:ParseGrid()
     771           6 :             m_bHasNoData = true;
     772             :         }
     773             :     }
     774             : }
     775             : 
     776             : /************************************************************************/
     777             : /*                         GetDescription()                             */
     778             : /************************************************************************/
     779             : 
     780          29 : const char *GRIBRasterBand::GetDescription() const
     781             : {
     782          29 :     if (longFstLevel == nullptr)
     783           0 :         return GDALPamRasterBand::GetDescription();
     784             : 
     785          29 :     return longFstLevel;
     786             : }
     787             : 
     788             : /************************************************************************/
     789             : /*                             LoadData()                               */
     790             : /************************************************************************/
     791             : 
     792        9152 : CPLErr GRIBRasterBand::LoadData()
     793             : 
     794             : {
     795        9152 :     if (!m_Grib_Data)
     796             :     {
     797         184 :         GRIBDataset *poGDS = static_cast<GRIBDataset *>(poDS);
     798             : 
     799         184 :         if (poGDS->bCacheOnlyOneBand)
     800             :         {
     801             :             // In "one-band-at-a-time" strategy, if the last recently used
     802             :             // band is not that one, uncache it. We could use a smarter strategy
     803             :             // based on a LRU, but that's a bit overkill for now.
     804           0 :             poGDS->poLastUsedBand->UncacheData();
     805           0 :             poGDS->nCachedBytes = 0;
     806             :         }
     807             :         else
     808             :         {
     809             :             // Once we have cached more than nCachedBytesThreshold bytes, we
     810             :             // will switch to "one-band-at-a-time" strategy, instead of caching
     811             :             // all bands that have been accessed.
     812         184 :             if (poGDS->nCachedBytes > poGDS->nCachedBytesThreshold)
     813             :             {
     814             :                 GUIntBig nMinCacheSize =
     815           0 :                     1 + static_cast<GUIntBig>(poGDS->nRasterXSize) *
     816           0 :                             poGDS->nRasterYSize * poGDS->nBands *
     817           0 :                             GDALGetDataTypeSizeBytes(eDataType) / 1024 / 1024;
     818           0 :                 CPLDebug("GRIB",
     819             :                          "Maximum band cache size reached for this dataset. "
     820             :                          "Caching only one band at a time from now, which can "
     821             :                          "negatively affect performance. Consider "
     822             :                          "increasing GRIB_CACHEMAX to a higher value (in MB), "
     823             :                          "at least " CPL_FRMT_GUIB " in that instance",
     824             :                          nMinCacheSize);
     825           0 :                 for (int i = 0; i < poGDS->nBands; i++)
     826             :                 {
     827             :                     reinterpret_cast<GRIBRasterBand *>(
     828           0 :                         poGDS->GetRasterBand(i + 1))
     829           0 :                         ->UncacheData();
     830             :                 }
     831           0 :                 poGDS->nCachedBytes = 0;
     832           0 :                 poGDS->bCacheOnlyOneBand = TRUE;
     833             :             }
     834             :         }
     835             : 
     836             :         // we don't seem to have any way to detect errors in this!
     837         184 :         if (m_Grib_MetaData != nullptr)
     838             :         {
     839         137 :             MetaFree(m_Grib_MetaData);
     840         137 :             delete m_Grib_MetaData;
     841         137 :             m_Grib_MetaData = nullptr;
     842             :         }
     843         184 :         ReadGribData(poGDS->fp, start, subgNum, &m_Grib_Data, &m_Grib_MetaData);
     844         184 :         if (!m_Grib_Data)
     845             :         {
     846           0 :             CPLError(CE_Failure, CPLE_AppDefined, "Out of memory.");
     847           0 :             if (m_Grib_MetaData != nullptr)
     848             :             {
     849           0 :                 MetaFree(m_Grib_MetaData);
     850           0 :                 delete m_Grib_MetaData;
     851           0 :                 m_Grib_MetaData = nullptr;
     852             :             }
     853           0 :             return CE_Failure;
     854             :         }
     855             : 
     856             :         // Check the band matches the dataset as a whole, size wise. (#3246)
     857         184 :         nGribDataXSize = m_Grib_MetaData->gds.Nx;
     858         184 :         nGribDataYSize = m_Grib_MetaData->gds.Ny;
     859         184 :         if (nGribDataXSize <= 0 || nGribDataYSize <= 0)
     860             :         {
     861           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     862             :                      "Band %d of GRIB dataset is %dx%d.", nBand, nGribDataXSize,
     863             :                      nGribDataYSize);
     864           0 :             MetaFree(m_Grib_MetaData);
     865           0 :             delete m_Grib_MetaData;
     866           0 :             m_Grib_MetaData = nullptr;
     867           0 :             return CE_Failure;
     868             :         }
     869             : 
     870         184 :         poGDS->nCachedBytes += static_cast<GIntBig>(nGribDataXSize) *
     871         184 :                                nGribDataYSize * sizeof(double);
     872         184 :         poGDS->poLastUsedBand = this;
     873             : 
     874         184 :         if (nGribDataXSize != nRasterXSize || nGribDataYSize != nRasterYSize)
     875             :         {
     876          24 :             CPLError(CE_Warning, CPLE_AppDefined,
     877             :                      "Band %d of GRIB dataset is %dx%d, while the first band "
     878             :                      "and dataset is %dx%d.  Georeferencing of band %d may "
     879             :                      "be incorrect, and data access may be incomplete.",
     880             :                      nBand, nGribDataXSize, nGribDataYSize, nRasterXSize,
     881             :                      nRasterYSize, nBand);
     882             :         }
     883             :     }
     884             : 
     885        9152 :     return CE_None;
     886             : }
     887             : 
     888             : /************************************************************************/
     889             : /*                       IsGdalinfoInteractive()                        */
     890             : /************************************************************************/
     891             : 
     892             : #ifdef BUILD_APPS
     893           0 : static bool IsGdalinfoInteractive()
     894             : {
     895           0 :     static const bool bIsGdalinfoInteractive = []()
     896             :     {
     897           0 :         if (CPLIsInteractive(stdout))
     898             :         {
     899           0 :             std::string osPath;
     900           0 :             osPath.resize(1024);
     901           0 :             if (CPLGetExecPath(&osPath[0], static_cast<int>(osPath.size())))
     902             :             {
     903           0 :                 osPath = CPLGetBasenameSafe(osPath.c_str());
     904             :             }
     905           0 :             return osPath == "gdalinfo";
     906             :         }
     907           0 :         return false;
     908           0 :     }();
     909           0 :     return bIsGdalinfoInteractive;
     910             : }
     911             : #endif
     912             : 
     913             : /************************************************************************/
     914             : /*                             GetMetaData()                            */
     915             : /************************************************************************/
     916          87 : char **GRIBRasterBand::GetMetadata(const char *pszDomain)
     917             : {
     918          87 :     FindMetaData();
     919         138 :     if ((pszDomain == nullptr || pszDomain[0] == 0) && m_nGribVersion == 2 &&
     920          51 :         CPLTestBool(CPLGetConfigOption("GRIB_PDS_ALL_BANDS", "ON")))
     921             :     {
     922             : #ifdef BUILD_APPS
     923             :         // Detect slow execution of e.g.
     924             :         // "gdalinfo /vsis3/noaa-hrrr-bdp-pds/hrrr.20220804/conus/hrrr.t00z.wrfsfcf01.grib2"
     925          51 :         GRIBDataset *poGDS = static_cast<GRIBDataset *>(poDS);
     926          11 :         if (poGDS->m_bSideCarIdxUsed && !poGDS->m_bWarnedGdalinfoNomd &&
     927          11 :             poGDS->GetRasterCount() > 10 &&
     928          62 :             !VSIIsLocal(poGDS->GetDescription()) && IsGdalinfoInteractive())
     929             :         {
     930           0 :             if (poGDS->m_nFirstMetadataQueriedTimeStamp)
     931             :             {
     932           0 :                 if (time(nullptr) - poGDS->m_nFirstMetadataQueriedTimeStamp > 2)
     933             :                 {
     934           0 :                     poGDS->m_bWarnedGdalinfoNomd = true;
     935             : 
     936           0 :                     CPLError(
     937             :                         CE_Warning, CPLE_AppDefined,
     938             :                         "If metadata does not matter, faster result could be "
     939             :                         "obtained by adding the -nomd switch to gdalinfo");
     940             :                 }
     941             :             }
     942             :             else
     943             :             {
     944           0 :                 poGDS->m_nFirstMetadataQueriedTimeStamp = time(nullptr);
     945             :             }
     946             :         }
     947             : #endif
     948             : 
     949          51 :         FindPDSTemplateGRIB2();
     950             :     }
     951          87 :     return GDALPamRasterBand::GetMetadata(pszDomain);
     952             : }
     953             : 
     954             : /************************************************************************/
     955             : /*                             GetMetaDataItem()                        */
     956             : /************************************************************************/
     957         693 : const char *GRIBRasterBand::GetMetadataItem(const char *pszName,
     958             :                                             const char *pszDomain)
     959             : {
     960         693 :     if (!((!pszDomain || pszDomain[0] == 0) &&
     961         681 :           (EQUAL(pszName, "STATISTICS_MINIMUM") ||
     962         677 :            EQUAL(pszName, "STATISTICS_MAXIMUM"))))
     963             :     {
     964         688 :         FindMetaData();
     965        1221 :         if (m_nGribVersion == 2 &&
     966         533 :             CPLTestBool(CPLGetConfigOption("GRIB_PDS_ALL_BANDS", "ON")))
     967             :         {
     968         532 :             FindPDSTemplateGRIB2();
     969             :         }
     970             :     }
     971         693 :     return GDALPamRasterBand::GetMetadataItem(pszName, pszDomain);
     972             : }
     973             : 
     974             : /************************************************************************/
     975             : /*                             IReadBlock()                             */
     976             : /************************************************************************/
     977             : 
     978        9152 : CPLErr GRIBRasterBand::IReadBlock(int /* nBlockXOff */, int nBlockYOff,
     979             :                                   void *pImage)
     980             : 
     981             : {
     982        9152 :     CPLErr eErr = LoadData();
     983        9152 :     if (eErr != CE_None)
     984           0 :         return eErr;
     985             : 
     986        9152 :     GRIBDataset *poGDS = static_cast<GRIBDataset *>(poDS);
     987             : 
     988             :     // The image as read is always upside down to our normal
     989             :     // orientation so we need to effectively flip it at this
     990             :     // point.  We also need to deal with bands that are a different
     991             :     // size than the dataset as a whole.
     992             : 
     993        9152 :     if (nGribDataXSize == nRasterXSize && nGribDataYSize == nRasterYSize &&
     994        8660 :         poGDS->nSplitAndSwapColumn == 0)
     995             :     {
     996             :         // Simple 1:1 case.
     997        7001 :         memcpy(pImage,
     998        7001 :                m_Grib_Data + static_cast<size_t>(nRasterXSize) *
     999        7001 :                                  (nRasterYSize - nBlockYOff - 1),
    1000        7001 :                nRasterXSize * sizeof(double));
    1001             : 
    1002        7001 :         return CE_None;
    1003             :     }
    1004             : 
    1005        2151 :     memset(pImage, 0, sizeof(double) * nRasterXSize);
    1006             : 
    1007        2151 :     if (nBlockYOff >= nGribDataYSize)  // Off image?
    1008          57 :         return CE_None;
    1009             : 
    1010        2094 :     int nSplitAndSwapColumn = poGDS->nSplitAndSwapColumn;
    1011        2094 :     if (nRasterXSize != nGribDataXSize)
    1012         435 :         nSplitAndSwapColumn = 0;
    1013             : 
    1014        2094 :     const int nCopyWords = std::min(nRasterXSize, nGribDataXSize);
    1015             : 
    1016        2094 :     memcpy(pImage,
    1017        2094 :            m_Grib_Data +
    1018        2094 :                static_cast<size_t>(nGribDataXSize) *
    1019        2094 :                    (nGribDataYSize - nBlockYOff - 1) +
    1020             :                nSplitAndSwapColumn,
    1021        2094 :            (nCopyWords - nSplitAndSwapColumn) * sizeof(double));
    1022             : 
    1023        2094 :     if (nSplitAndSwapColumn > 0)
    1024        1659 :         memcpy(reinterpret_cast<void *>(reinterpret_cast<double *>(pImage) +
    1025        1659 :                                         nCopyWords - nSplitAndSwapColumn),
    1026        1659 :                m_Grib_Data + static_cast<size_t>(nGribDataXSize) *
    1027        1659 :                                  (nGribDataYSize - nBlockYOff - 1),
    1028        1659 :                nSplitAndSwapColumn * sizeof(double));
    1029             : 
    1030        2094 :     return CE_None;
    1031             : }
    1032             : 
    1033             : /************************************************************************/
    1034             : /*                           GetNoDataValue()                           */
    1035             : /************************************************************************/
    1036             : 
    1037         125 : double GRIBRasterBand::GetNoDataValue(int *pbSuccess)
    1038             : {
    1039         125 :     if (m_bHasLookedForNoData)
    1040             :     {
    1041         105 :         if (pbSuccess)
    1042         105 :             *pbSuccess = m_bHasNoData;
    1043         105 :         return m_dfNoData;
    1044             :     }
    1045             : 
    1046          20 :     m_bHasLookedForNoData = true;
    1047          20 :     if (m_Grib_MetaData == nullptr)
    1048             :     {
    1049          18 :         GRIBDataset *poGDS = static_cast<GRIBDataset *>(poDS);
    1050             : 
    1051             : #ifdef BUILD_APPS
    1052             :         // Detect slow execution of e.g.
    1053             :         // "gdalinfo /vsis3/noaa-hrrr-bdp-pds/hrrr.20220804/conus/hrrr.t00z.wrfsfcf01.grib2"
    1054             : 
    1055           0 :         if (poGDS->m_bSideCarIdxUsed && !poGDS->m_bWarnedGdalinfoNonodata &&
    1056           0 :             poGDS->GetRasterCount() > 10 &&
    1057          18 :             !VSIIsLocal(poGDS->GetDescription()) && IsGdalinfoInteractive())
    1058             :         {
    1059           0 :             if (poGDS->m_nFirstNodataQueriedTimeStamp)
    1060             :             {
    1061           0 :                 if (time(nullptr) - poGDS->m_nFirstNodataQueriedTimeStamp > 2)
    1062             :                 {
    1063           0 :                     poGDS->m_bWarnedGdalinfoNonodata = true;
    1064             : 
    1065           0 :                     CPLError(CE_Warning, CPLE_AppDefined,
    1066             :                              "If nodata value does not matter, faster result "
    1067             :                              "could be obtained by adding the -nonodata switch "
    1068             :                              "to gdalinfo");
    1069             :                 }
    1070             :             }
    1071             :             else
    1072             :             {
    1073           0 :                 poGDS->m_nFirstNodataQueriedTimeStamp = time(nullptr);
    1074             :             }
    1075             :         }
    1076             : #endif
    1077             : 
    1078          18 :         ReadGribData(poGDS->fp, start, subgNum, nullptr, &m_Grib_MetaData);
    1079          18 :         if (m_Grib_MetaData == nullptr)
    1080             :         {
    1081           0 :             m_bHasNoData = false;
    1082           0 :             m_dfNoData = 0;
    1083           0 :             if (pbSuccess)
    1084           0 :                 *pbSuccess = m_bHasNoData;
    1085           0 :             return m_dfNoData;
    1086             :         }
    1087             :     }
    1088             : 
    1089          20 :     if (m_Grib_MetaData->gridAttrib.f_miss == 0)
    1090             :     {
    1091           4 :         m_bHasNoData = false;
    1092           4 :         m_dfNoData = 0;
    1093           4 :         if (pbSuccess)
    1094           4 :             *pbSuccess = m_bHasNoData;
    1095           4 :         return m_dfNoData;
    1096             :     }
    1097             : 
    1098          16 :     if (m_Grib_MetaData->gridAttrib.f_miss == 2)
    1099             :     {
    1100             :         // What TODO?
    1101           0 :         CPLDebug("GRIB", "Secondary missing value also set for band %d : %f",
    1102           0 :                  nBand, m_Grib_MetaData->gridAttrib.missSec);
    1103             :     }
    1104             : 
    1105          16 :     m_bHasNoData = true;
    1106          16 :     m_dfNoData = m_Grib_MetaData->gridAttrib.missPri;
    1107          16 :     if (pbSuccess)
    1108          16 :         *pbSuccess = m_bHasNoData;
    1109          16 :     return m_dfNoData;
    1110             : }
    1111             : 
    1112             : /************************************************************************/
    1113             : /*                            ReadGribData()                            */
    1114             : /************************************************************************/
    1115             : 
    1116         538 : void GRIBRasterBand::ReadGribData(VSILFILE *fp, vsi_l_offset start, int subgNum,
    1117             :                                   double **data, grib_MetaData **metaData)
    1118             : {
    1119             :     // Initialization, for calling the ReadGrib2Record function.
    1120         538 :     sInt4 f_endMsg = 1;  // 1 if we read the last grid in a GRIB message, or we
    1121             :                          // haven't read any messages.
    1122             :     // int subgNum = 0; // The subgrid in the message that we are interested in.
    1123         538 :     sChar f_unit = 2;       // None = 0, English = 1, Metric = 2
    1124         538 :     double majEarth = 0.0;  // -radEarth if < 6000 ignore, otherwise use this
    1125             :                             // to override the radEarth in the GRIB1 or GRIB2
    1126             :                             // message.  Needed because NCEP uses 6371.2 but
    1127             :                             // GRIB1 could only state 6367.47.
    1128         538 :     double minEarth = 0.0;  // -minEarth if < 6000 ignore, otherwise use this
    1129             :                             // to override the minEarth in the GRIB1 or GRIB2
    1130             :                             // message.
    1131         538 :     sChar f_SimpleVer = 4;  // Which version of the simple NDFD Weather table
    1132             :                             // to use. (1 is 6/2003) (2 is 1/2004) (3 is
    1133             :                             // 2/2004) (4 is 11/2004) (default 4)
    1134             :     LatLon lwlf;            // Lower left corner (cookie slicing) -lwlf
    1135             :     LatLon uprt;            // Upper right corner (cookie slicing) -uprt
    1136             :     IS_dataType is;  // Un-parsed meta data for this GRIB2 message. As well
    1137             :                      // as some memory used by the unpacker.
    1138             : 
    1139         538 :     lwlf.lat = -100;  // lat == -100 instructs the GRIB decoder that we don't
    1140             :                       // want a subgrid
    1141             : 
    1142         538 :     IS_Init(&is);
    1143             : 
    1144             :     const char *pszGribNormalizeUnits =
    1145         538 :         CPLGetConfigOption("GRIB_NORMALIZE_UNITS", "YES");
    1146         538 :     if (!CPLTestBool(pszGribNormalizeUnits))
    1147           2 :         f_unit = 0;  // Do not normalize units to metric.
    1148             : 
    1149         538 :     start = FindTrueStart(fp, start);
    1150             :     // Read GRIB message from file position "start".
    1151         538 :     VSIFSeekL(fp, start, SEEK_SET);
    1152         538 :     uInt4 grib_DataLen = 0;  // Size of Grib_Data.
    1153         538 :     *metaData = new grib_MetaData();
    1154         538 :     MetaInit(*metaData);
    1155         538 :     const int simpWWA = 0;  // seem to be unused in degrib
    1156         538 :     ReadGrib2Record(fp, f_unit, data, &grib_DataLen, *metaData, &is, subgNum,
    1157             :                     majEarth, minEarth, f_SimpleVer, simpWWA, &f_endMsg, &lwlf,
    1158             :                     &uprt);
    1159             : 
    1160             :     // No intention to show errors, just swallow it and free the memory.
    1161         538 :     char *errMsg = errSprintf(nullptr);
    1162         538 :     if (errMsg != nullptr)
    1163         341 :         CPLDebug("GRIB", "%s", errMsg);
    1164         538 :     free(errMsg);
    1165         538 :     IS_Free(&is);
    1166         538 : }
    1167             : 
    1168             : /************************************************************************/
    1169             : /*                            UncacheData()                             */
    1170             : /************************************************************************/
    1171             : 
    1172         444 : void GRIBRasterBand::UncacheData()
    1173             : {
    1174         444 :     if (m_Grib_Data)
    1175         184 :         free(m_Grib_Data);
    1176         444 :     m_Grib_Data = nullptr;
    1177         444 :     if (m_Grib_MetaData)
    1178             :     {
    1179         373 :         MetaFree(m_Grib_MetaData);
    1180         373 :         delete m_Grib_MetaData;
    1181             :     }
    1182         444 :     m_Grib_MetaData = nullptr;
    1183         444 : }
    1184             : 
    1185             : /************************************************************************/
    1186             : /*                           ~GRIBRasterBand()                          */
    1187             : /************************************************************************/
    1188             : 
    1189         865 : GRIBRasterBand::~GRIBRasterBand()
    1190             : {
    1191         444 :     if (longFstLevel != nullptr)
    1192         444 :         CPLFree(longFstLevel);
    1193         444 :     UncacheData();
    1194         865 : }
    1195             : 
    1196             : gdal::grib::InventoryWrapper::~InventoryWrapper() = default;
    1197             : 
    1198             : /************************************************************************/
    1199             : /*                           InventoryWrapperGrib                       */
    1200             : /************************************************************************/
    1201             : class InventoryWrapperGrib final : public gdal::grib::InventoryWrapper
    1202             : {
    1203             :   public:
    1204         310 :     explicit InventoryWrapperGrib(VSILFILE *fp) : gdal::grib::InventoryWrapper()
    1205             :     {
    1206         310 :         result_ = GRIB2Inventory(fp, &inv_, &inv_len_, 0 /* all messages */,
    1207             :                                  &num_messages_);
    1208         310 :     }
    1209             : 
    1210             :     ~InventoryWrapperGrib() override;
    1211             : };
    1212             : 
    1213         620 : InventoryWrapperGrib::~InventoryWrapperGrib()
    1214             : {
    1215         310 :     if (inv_ == nullptr)
    1216           0 :         return;
    1217         732 :     for (uInt4 i = 0; i < inv_len_; i++)
    1218             :     {
    1219         422 :         GRIB2InventoryFree(inv_ + i);
    1220             :     }
    1221         310 :     free(inv_);
    1222         620 : }
    1223             : 
    1224             : /************************************************************************/
    1225             : /*                           InventoryWrapperSidecar                    */
    1226             : /************************************************************************/
    1227             : 
    1228             : class InventoryWrapperSidecar final : public gdal::grib::InventoryWrapper
    1229             : {
    1230             :   public:
    1231           6 :     explicit InventoryWrapperSidecar(VSILFILE *fp, uint64_t nStartOffset,
    1232             :                                      int64_t nSize)
    1233           6 :         : gdal::grib::InventoryWrapper()
    1234             :     {
    1235           6 :         result_ = -1;
    1236           6 :         VSIFSeekL(fp, 0, SEEK_END);
    1237           6 :         size_t length = static_cast<size_t>(VSIFTellL(fp));
    1238           6 :         if (length > 4 * 1024 * 1024)
    1239           0 :             return;
    1240           6 :         std::string osSidecar;
    1241           6 :         osSidecar.resize(length);
    1242           6 :         VSIFSeekL(fp, 0, SEEK_SET);
    1243           6 :         if (VSIFReadL(&osSidecar[0], length, 1, fp) != 1)
    1244           0 :             return;
    1245             : 
    1246             :         const CPLStringList aosMsgs(
    1247             :             CSLTokenizeString2(osSidecar.c_str(), "\n",
    1248           6 :                                CSLT_PRESERVEQUOTES | CSLT_STRIPLEADSPACES));
    1249           6 :         inv_ = static_cast<inventoryType *>(
    1250           6 :             CPLCalloc(aosMsgs.size(), sizeof(inventoryType)));
    1251             : 
    1252          42 :         for (const char *pszMsg : aosMsgs)
    1253             :         {
    1254             :             // We are parsing
    1255             :             // "msgNum[.subgNum]:start:dontcare:name1:name2:name3" For NOMADS:
    1256             :             // "msgNum[.subgNum]:start:reftime:var:level:time"
    1257             :             const CPLStringList aosTokens(CSLTokenizeString2(
    1258          38 :                 pszMsg, ":", CSLT_PRESERVEQUOTES | CSLT_ALLOWEMPTYTOKENS));
    1259          38 :             CPLStringList aosNum;
    1260             : 
    1261          38 :             if (aosTokens.size() < 6)
    1262           0 :                 goto err_sidecar;
    1263             : 
    1264          38 :             aosNum = CPLStringList(CSLTokenizeString2(aosTokens[0], ".", 0));
    1265          38 :             if (aosNum.size() < 1)
    1266           0 :                 goto err_sidecar;
    1267             : 
    1268             :             // FindMetaData will retrieve the correct version number
    1269             :             char *endptr;
    1270          38 :             strtol(aosNum[0], &endptr, 10);
    1271          38 :             if (*endptr != 0)
    1272           0 :                 goto err_sidecar;
    1273             : 
    1274          38 :             if (aosNum.size() < 2)
    1275          36 :                 inv_[inv_len_].subgNum = 0;
    1276             :             else
    1277             :             {
    1278           2 :                 auto subgNum = strtol(aosNum[1], &endptr, 10);
    1279           2 :                 if (*endptr != 0)
    1280           0 :                     goto err_sidecar;
    1281           2 :                 if (subgNum <= 0 || subgNum > 65536)
    1282           0 :                     goto err_sidecar;
    1283             :                 // .idx file use a 1-based indexing, whereas DEGRIB uses a
    1284             :                 // 0-based one
    1285           2 :                 subgNum--;
    1286           2 :                 inv_[inv_len_].subgNum = static_cast<unsigned short>(subgNum);
    1287             :             }
    1288             : 
    1289          38 :             inv_[inv_len_].start = strtoll(aosTokens[1], &endptr, 10);
    1290          38 :             if (*endptr != 0)
    1291           0 :                 goto err_sidecar;
    1292             : 
    1293          38 :             if (inv_[inv_len_].start < nStartOffset)
    1294           4 :                 continue;
    1295          34 :             if (nSize > 0 && inv_[inv_len_].start >= nStartOffset + nSize)
    1296           2 :                 break;
    1297             : 
    1298          32 :             inv_[inv_len_].start -= nStartOffset;
    1299             : 
    1300          32 :             inv_[inv_len_].unitName = nullptr;
    1301          32 :             inv_[inv_len_].comment = nullptr;
    1302          32 :             inv_[inv_len_].element = nullptr;
    1303          32 :             inv_[inv_len_].shortFstLevel = nullptr;
    1304             :             // This is going into the description field ->
    1305             :             // the only one available before loading the metadata
    1306          32 :             inv_[inv_len_].longFstLevel = VSIStrdup(CPLSPrintf(
    1307             :                 "%s:%s:%s", aosTokens[3], aosTokens[4], aosTokens[5]));
    1308          32 :             ++inv_len_;
    1309             : 
    1310          32 :             continue;
    1311             : 
    1312           0 :         err_sidecar:
    1313           0 :             CPLDebug("GRIB",
    1314             :                      "Failed parsing sidecar entry '%s', "
    1315             :                      "falling back to constructing an inventory",
    1316             :                      pszMsg);
    1317           0 :             return;
    1318             :         }
    1319             : 
    1320           6 :         result_ = inv_len_;
    1321             :     }
    1322             : 
    1323             :     ~InventoryWrapperSidecar() override;
    1324             : };
    1325             : 
    1326          12 : InventoryWrapperSidecar::~InventoryWrapperSidecar()
    1327             : {
    1328           6 :     if (inv_ == nullptr)
    1329           0 :         return;
    1330             : 
    1331          38 :     for (unsigned i = 0; i < inv_len_; i++)
    1332          32 :         VSIFree(inv_[i].longFstLevel);
    1333             : 
    1334           6 :     VSIFree(inv_);
    1335          12 : }
    1336             : 
    1337             : /************************************************************************/
    1338             : /* ==================================================================== */
    1339             : /*                              GRIBDataset                             */
    1340             : /* ==================================================================== */
    1341             : /************************************************************************/
    1342             : 
    1343         316 : GRIBDataset::GRIBDataset()
    1344             :     : fp(nullptr), nCachedBytes(0),
    1345             :       // Switch caching strategy once 100 MB threshold is reached.
    1346             :       // Why 100 MB? --> Why not.
    1347         632 :       nCachedBytesThreshold(static_cast<GIntBig>(atoi(
    1348             :                                 CPLGetConfigOption("GRIB_CACHEMAX", "100"))) *
    1349         316 :                             1024 * 1024),
    1350         316 :       bCacheOnlyOneBand(FALSE), nSplitAndSwapColumn(0), poLastUsedBand(nullptr)
    1351             : {
    1352         316 : }
    1353             : 
    1354             : /************************************************************************/
    1355             : /*                            ~GRIBDataset()                             */
    1356             : /************************************************************************/
    1357             : 
    1358         632 : GRIBDataset::~GRIBDataset()
    1359             : 
    1360             : {
    1361         316 :     FlushCache(true);
    1362         316 :     if (fp != nullptr)
    1363         312 :         VSIFCloseL(fp);
    1364         632 : }
    1365             : 
    1366             : /************************************************************************/
    1367             : /*                          GetGeoTransform()                           */
    1368             : /************************************************************************/
    1369             : 
    1370         135 : CPLErr GRIBDataset::GetGeoTransform(GDALGeoTransform &gt) const
    1371             : 
    1372             : {
    1373         135 :     gt = m_gt;
    1374         135 :     return CE_None;
    1375             : }
    1376             : 
    1377             : /************************************************************************/
    1378             : /*                                Inventory()                           */
    1379             : /************************************************************************/
    1380             : 
    1381             : std::unique_ptr<gdal::grib::InventoryWrapper>
    1382         312 : GRIBDataset::Inventory(GDALOpenInfo *poOpenInfo)
    1383             : {
    1384         312 :     std::unique_ptr<gdal::grib::InventoryWrapper> pInventories;
    1385             : 
    1386         312 :     VSIFSeekL(fp, 0, SEEK_SET);
    1387         624 :     std::string osSideCarFilename(poOpenInfo->pszFilename);
    1388         312 :     uint64_t nStartOffset = 0;
    1389         312 :     int64_t nSize = -1;
    1390         312 :     if (STARTS_WITH(poOpenInfo->pszFilename, "/vsisubfile/"))
    1391             :     {
    1392           5 :         const char *pszPtr = poOpenInfo->pszFilename + strlen("/vsisubfile/");
    1393           5 :         const char *pszComma = strchr(pszPtr, ',');
    1394           5 :         if (pszComma)
    1395             :         {
    1396             :             const CPLStringList aosTokens(CSLTokenizeString2(
    1397          15 :                 std::string(pszPtr, pszComma - pszPtr).c_str(), "_", 0));
    1398           5 :             if (aosTokens.size() == 2)
    1399             :             {
    1400           5 :                 nStartOffset = std::strtoull(aosTokens[0], nullptr, 10);
    1401           5 :                 nSize = std::strtoll(aosTokens[1], nullptr, 10);
    1402           5 :                 osSideCarFilename = pszComma + 1;
    1403             :             }
    1404             :         }
    1405             :     }
    1406         312 :     osSideCarFilename += ".idx";
    1407         312 :     VSILFILE *fpSideCar = nullptr;
    1408         312 :     if (CPLTestBool(CSLFetchNameValueDef(poOpenInfo->papszOpenOptions,
    1409         619 :                                          "USE_IDX", "YES")) &&
    1410         307 :         ((fpSideCar = VSIFOpenL(osSideCarFilename.c_str(), "rb")) != nullptr))
    1411             :     {
    1412           6 :         CPLDebug("GRIB", "Reading inventories from sidecar file %s",
    1413             :                  osSideCarFilename.c_str());
    1414             :         // Contains an GRIB2 message inventory of the file.
    1415          12 :         pInventories = std::make_unique<InventoryWrapperSidecar>(
    1416           6 :             fpSideCar, nStartOffset, nSize);
    1417           6 :         if (pInventories->result() <= 0 || pInventories->length() == 0)
    1418           0 :             pInventories = nullptr;
    1419           6 :         VSIFCloseL(fpSideCar);
    1420             : #ifdef BUILD_APPS
    1421           6 :         m_bSideCarIdxUsed = true;
    1422             : #endif
    1423             :     }
    1424             :     else
    1425         306 :         CPLDebug("GRIB", "Failed opening sidecar %s",
    1426             :                  osSideCarFilename.c_str());
    1427             : 
    1428         312 :     if (pInventories == nullptr)
    1429             :     {
    1430         306 :         CPLDebug("GRIB", "Reading inventories from GRIB file %s",
    1431             :                  poOpenInfo->pszFilename);
    1432             :         // Contains an GRIB2 message inventory of the file.
    1433         306 :         pInventories = std::make_unique<InventoryWrapperGrib>(fp);
    1434             :     }
    1435             : 
    1436         624 :     return pInventories;
    1437             : }
    1438             : 
    1439             : /************************************************************************/
    1440             : /*                                Open()                                */
    1441             : /************************************************************************/
    1442             : 
    1443         318 : GDALDataset *GRIBDataset::Open(GDALOpenInfo *poOpenInfo)
    1444             : 
    1445             : {
    1446             : #ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
    1447             :     // During fuzzing, do not use Identify to reject crazy content.
    1448         318 :     if (!GRIBDriverIdentify(poOpenInfo))
    1449           2 :         return nullptr;
    1450             : #endif
    1451         316 :     if (poOpenInfo->fpL == nullptr)
    1452           0 :         return nullptr;
    1453             : 
    1454             :     // A fast "probe" on the header that is partially read in memory.
    1455         316 :     char *buff = nullptr;
    1456         316 :     uInt4 buffLen = 0;
    1457         316 :     sInt4 sect0[SECT0LEN_WORD] = {0};
    1458         316 :     uInt4 gribLen = 0;
    1459         316 :     int version = 0;
    1460             : 
    1461             :     // grib is not thread safe, make sure not to cause problems
    1462             :     // for other thread safe formats
    1463         632 :     CPLMutexHolderD(&hGRIBMutex);
    1464             : 
    1465         632 :     VSILFILE *memfp = VSIFileFromMemBuffer(nullptr, poOpenInfo->pabyHeader,
    1466         316 :                                            poOpenInfo->nHeaderBytes, FALSE);
    1467         632 :     if (memfp == nullptr ||
    1468         316 :         ReadSECT0(memfp, &buff, &buffLen, -1, sect0, &gribLen, &version) < 0)
    1469             :     {
    1470           0 :         if (memfp != nullptr)
    1471             :         {
    1472           0 :             VSIFCloseL(memfp);
    1473             :         }
    1474           0 :         free(buff);
    1475           0 :         char *errMsg = errSprintf(nullptr);
    1476           0 :         if (errMsg != nullptr && strstr(errMsg, "Ran out of file") == nullptr)
    1477           0 :             CPLDebug("GRIB", "%s", errMsg);
    1478           0 :         free(errMsg);
    1479           0 :         return nullptr;
    1480             :     }
    1481         316 :     VSIFCloseL(memfp);
    1482         316 :     free(buff);
    1483             : 
    1484             :     // Confirm the requested access is supported.
    1485         316 :     if (poOpenInfo->eAccess == GA_Update)
    1486             :     {
    1487           0 :         ReportUpdateNotSupportedByDriver("GRIB");
    1488           0 :         return nullptr;
    1489             :     }
    1490             : 
    1491         316 :     if (poOpenInfo->nOpenFlags & GDAL_OF_MULTIDIM_RASTER)
    1492             :     {
    1493           4 :         return OpenMultiDim(poOpenInfo);
    1494             :     }
    1495             : 
    1496             :     // Create a corresponding GDALDataset.
    1497         312 :     GRIBDataset *poDS = new GRIBDataset();
    1498             : 
    1499         312 :     poDS->fp = poOpenInfo->fpL;
    1500         312 :     poOpenInfo->fpL = nullptr;
    1501             : 
    1502             :     // Make an inventory of the GRIB file.
    1503             :     // The inventory does not contain all the information needed for
    1504             :     // creating the RasterBands (especially the x and y size), therefore
    1505             :     // the first GRIB band is also read for some additional metadata.
    1506             :     // The band-data that is read is stored into the first RasterBand,
    1507             :     // simply so that the same portion of the file is not read twice.
    1508             : 
    1509         624 :     auto pInventories = poDS->Inventory(poOpenInfo);
    1510         312 :     if (pInventories->result() <= 0)
    1511             :     {
    1512           9 :         char *errMsg = errSprintf(nullptr);
    1513           9 :         if (errMsg != nullptr)
    1514           8 :             CPLDebug("GRIB", "%s", errMsg);
    1515           9 :         free(errMsg);
    1516             : 
    1517           9 :         CPLError(CE_Failure, CPLE_OpenFailed,
    1518             :                  "%s is a grib file, "
    1519             :                  "but no raster dataset was successfully identified.",
    1520             :                  poOpenInfo->pszFilename);
    1521             :         // Release hGRIBMutex otherwise we'll deadlock with GDALDataset own
    1522             :         // hGRIBMutex.
    1523           9 :         CPLReleaseMutex(hGRIBMutex);
    1524           9 :         delete poDS;
    1525           9 :         CPLAcquireMutex(hGRIBMutex, 1000.0);
    1526           9 :         return nullptr;
    1527             :     }
    1528             : 
    1529             :     // Create band objects.
    1530         303 :     const uInt4 nCount = std::min(pInventories->length(), 65536U);
    1531         724 :     for (uInt4 i = 0; i < nCount; ++i)
    1532             :     {
    1533         421 :         inventoryType *psInv = pInventories->get(i);
    1534         421 :         GRIBRasterBand *gribBand = nullptr;
    1535         421 :         const uInt4 bandNr = i + 1;
    1536         421 :         assert(bandNr <= 65536);
    1537             : 
    1538         421 :         if (bandNr == 1)
    1539             :         {
    1540             :             // Important: set DataSet extents before creating first RasterBand
    1541             :             // in it.
    1542         303 :             grib_MetaData *metaData = nullptr;
    1543         303 :             GRIBRasterBand::ReadGribData(poDS->fp, 0, psInv->subgNum, nullptr,
    1544             :                                          &metaData);
    1545         303 :             if (metaData == nullptr || metaData->gds.Nx < 1 ||
    1546         303 :                 metaData->gds.Ny < 1)
    1547             :             {
    1548           0 :                 CPLError(CE_Failure, CPLE_OpenFailed,
    1549             :                          "%s is a grib file, "
    1550             :                          "but no raster dataset was successfully identified.",
    1551             :                          poOpenInfo->pszFilename);
    1552             :                 // Release hGRIBMutex otherwise we'll deadlock with GDALDataset
    1553             :                 // own hGRIBMutex.
    1554           0 :                 CPLReleaseMutex(hGRIBMutex);
    1555           0 :                 delete poDS;
    1556           0 :                 CPLAcquireMutex(hGRIBMutex, 1000.0);
    1557           0 :                 if (metaData != nullptr)
    1558             :                 {
    1559           0 :                     MetaFree(metaData);
    1560           0 :                     delete metaData;
    1561             :                 }
    1562           0 :                 return nullptr;
    1563             :             }
    1564         303 :             psInv->GribVersion = metaData->GribVersion;
    1565             : 
    1566             :             // Set the DataSet's x,y size, georeference and projection from
    1567             :             // the first GRIB band.
    1568         303 :             poDS->SetGribMetaData(metaData);
    1569         303 :             gribBand = new GRIBRasterBand(poDS, bandNr, psInv);
    1570             : 
    1571         303 :             if (psInv->GribVersion == 2)
    1572         296 :                 gribBand->FindPDSTemplateGRIB2();
    1573             : 
    1574         303 :             gribBand->m_Grib_MetaData = metaData;
    1575             :         }
    1576             :         else
    1577             :         {
    1578         118 :             gribBand = new GRIBRasterBand(poDS, bandNr, psInv);
    1579             :         }
    1580         421 :         poDS->SetBand(bandNr, gribBand);
    1581             :     }
    1582             : 
    1583             :     // Initialize any PAM information.
    1584         303 :     poDS->SetDescription(poOpenInfo->pszFilename);
    1585             : 
    1586             :     // Release hGRIBMutex otherwise we'll deadlock with GDALDataset own
    1587             :     // hGRIBMutex.
    1588         303 :     CPLReleaseMutex(hGRIBMutex);
    1589         303 :     poDS->TryLoadXML(poOpenInfo->GetSiblingFiles());
    1590             : 
    1591             :     // Check for external overviews.
    1592         303 :     poDS->oOvManager.Initialize(poDS, poOpenInfo);
    1593         303 :     CPLAcquireMutex(hGRIBMutex, 1000.0);
    1594             : 
    1595         303 :     return poDS;
    1596             : }
    1597             : 
    1598             : /************************************************************************/
    1599             : /*                          GRIBSharedResource                          */
    1600             : /************************************************************************/
    1601             : 
    1602             : struct GRIBSharedResource
    1603             : {
    1604             :     VSILFILE *m_fp = nullptr;
    1605             :     vsi_l_offset m_nOffsetCurData = static_cast<vsi_l_offset>(-1);
    1606             :     std::vector<double> m_adfCurData{};
    1607             :     std::string m_osFilename;
    1608             :     std::shared_ptr<GDALPamMultiDim> m_poPAM{};
    1609             : 
    1610             :     GRIBSharedResource(const std::string &osFilename, VSILFILE *fp);
    1611             :     ~GRIBSharedResource();
    1612             : 
    1613             :     const std::vector<double> &LoadData(vsi_l_offset nOffset, int subgNum);
    1614             : 
    1615          23 :     const std::shared_ptr<GDALPamMultiDim> &GetPAM()
    1616             :     {
    1617          23 :         return m_poPAM;
    1618             :     }
    1619             : };
    1620             : 
    1621           4 : GRIBSharedResource::GRIBSharedResource(const std::string &osFilename,
    1622           4 :                                        VSILFILE *fp)
    1623             :     : m_fp(fp), m_osFilename(osFilename),
    1624           4 :       m_poPAM(std::make_shared<GDALPamMultiDim>(osFilename))
    1625             : {
    1626           4 : }
    1627             : 
    1628           4 : GRIBSharedResource::~GRIBSharedResource()
    1629             : {
    1630           4 :     if (m_fp)
    1631           4 :         VSIFCloseL(m_fp);
    1632           4 : }
    1633             : 
    1634             : /************************************************************************/
    1635             : /*                                GRIBGroup                             */
    1636             : /************************************************************************/
    1637             : 
    1638             : class GRIBArray;
    1639             : 
    1640             : class GRIBGroup final : public GDALGroup
    1641             : {
    1642             :     friend class GRIBArray;
    1643             :     std::shared_ptr<GRIBSharedResource> m_poShared{};
    1644             :     std::vector<std::shared_ptr<GDALMDArray>> m_poArrays{};
    1645             :     std::vector<std::shared_ptr<GDALDimension>> m_dims{};
    1646             :     std::map<std::string, std::shared_ptr<GDALDimension>> m_oMapDims{};
    1647             :     int m_nHorizDimCounter = 0;
    1648             :     std::shared_ptr<GDALGroup> m_memRootGroup{};
    1649             : 
    1650             :   public:
    1651           4 :     explicit GRIBGroup(const std::shared_ptr<GRIBSharedResource> &poShared)
    1652           4 :         : GDALGroup(std::string(), "/"), m_poShared(poShared)
    1653             :     {
    1654             :         std::unique_ptr<GDALDataset> poTmpDS(
    1655           4 :             MEMDataset::CreateMultiDimensional("", nullptr, nullptr));
    1656           4 :         m_memRootGroup = poTmpDS->GetRootGroup();
    1657           4 :     }
    1658             : 
    1659          36 :     void AddArray(const std::shared_ptr<GDALMDArray> &array)
    1660             :     {
    1661          36 :         m_poArrays.emplace_back(array);
    1662          36 :     }
    1663             : 
    1664             :     std::vector<std::string>
    1665             :     GetMDArrayNames(CSLConstList papszOptions) const override;
    1666             :     std::shared_ptr<GDALMDArray>
    1667             :     OpenMDArray(const std::string &osName,
    1668             :                 CSLConstList papszOptions) const override;
    1669             : 
    1670             :     std::vector<std::shared_ptr<GDALDimension>>
    1671           3 :     GetDimensions(CSLConstList) const override
    1672             :     {
    1673           3 :         return m_dims;
    1674             :     }
    1675             : };
    1676             : 
    1677             : /************************************************************************/
    1678             : /*                                GRIBArray                             */
    1679             : /************************************************************************/
    1680             : 
    1681             : class GRIBArray final : public GDALPamMDArray
    1682             : {
    1683             :     std::shared_ptr<GRIBSharedResource> m_poShared;
    1684             :     std::vector<std::shared_ptr<GDALDimension>> m_dims{};
    1685             :     GDALExtendedDataType m_dt = GDALExtendedDataType::Create(GDT_Float64);
    1686             :     std::shared_ptr<OGRSpatialReference> m_poSRS{};
    1687             :     std::vector<vsi_l_offset> m_anOffsets{};
    1688             :     std::vector<int> m_anSubgNums{};
    1689             :     std::vector<double> m_adfTimes{};
    1690             :     std::vector<std::shared_ptr<GDALAttribute>> m_attributes{};
    1691             :     std::string m_osUnit{};
    1692             :     std::vector<GByte> m_abyNoData{};
    1693             : 
    1694             :     GRIBArray(const std::string &osName,
    1695             :               const std::shared_ptr<GRIBSharedResource> &poShared);
    1696             : 
    1697             :   protected:
    1698             :     bool IRead(const GUInt64 *arrayStartIdx, const size_t *count,
    1699             :                const GInt64 *arrayStep, const GPtrDiff_t *bufferStride,
    1700             :                const GDALExtendedDataType &bufferDataType,
    1701             :                void *pDstBuffer) const override;
    1702             : 
    1703             :   public:
    1704             :     static std::shared_ptr<GRIBArray>
    1705          23 :     Create(const std::string &osName,
    1706             :            const std::shared_ptr<GRIBSharedResource> &poShared)
    1707             :     {
    1708          23 :         auto ar(std::shared_ptr<GRIBArray>(new GRIBArray(osName, poShared)));
    1709          23 :         ar->SetSelf(ar);
    1710          23 :         return ar;
    1711             :     }
    1712             : 
    1713             :     void Init(GRIBGroup *poGroup, GRIBDataset *poDS, GRIBRasterBand *poBand,
    1714             :               inventoryType *psInv);
    1715             :     void ExtendTimeDim(vsi_l_offset nOffset, int subgNum, double dfValidTime);
    1716             :     void Finalize(GRIBGroup *poGroup, inventoryType *psInv);
    1717             : 
    1718           0 :     bool IsWritable() const override
    1719             :     {
    1720           0 :         return false;
    1721             :     }
    1722             : 
    1723           4 :     const std::string &GetFilename() const override
    1724             :     {
    1725           4 :         return m_poShared->m_osFilename;
    1726             :     }
    1727             : 
    1728             :     const std::vector<std::shared_ptr<GDALDimension>> &
    1729          30 :     GetDimensions() const override
    1730             :     {
    1731          30 :         return m_dims;
    1732             :     }
    1733             : 
    1734          11 :     const GDALExtendedDataType &GetDataType() const override
    1735             :     {
    1736          11 :         return m_dt;
    1737             :     }
    1738             : 
    1739           1 :     std::shared_ptr<OGRSpatialReference> GetSpatialRef() const override
    1740             :     {
    1741           1 :         return m_poSRS;
    1742             :     }
    1743             : 
    1744             :     std::vector<std::shared_ptr<GDALAttribute>>
    1745          14 :     GetAttributes(CSLConstList) const override
    1746             :     {
    1747          14 :         return m_attributes;
    1748             :     }
    1749             : 
    1750           1 :     const std::string &GetUnit() const override
    1751             :     {
    1752           1 :         return m_osUnit;
    1753             :     }
    1754             : 
    1755           1 :     const void *GetRawNoDataValue() const override
    1756             :     {
    1757           1 :         return m_abyNoData.empty() ? nullptr : m_abyNoData.data();
    1758             :     }
    1759             : };
    1760             : 
    1761             : /************************************************************************/
    1762             : /*                          GetMDArrayNames()                           */
    1763             : /************************************************************************/
    1764             : 
    1765           3 : std::vector<std::string> GRIBGroup::GetMDArrayNames(CSLConstList) const
    1766             : {
    1767           3 :     std::vector<std::string> ret;
    1768          31 :     for (const auto &array : m_poArrays)
    1769             :     {
    1770          28 :         ret.push_back(array->GetName());
    1771             :     }
    1772           3 :     return ret;
    1773             : }
    1774             : 
    1775             : /************************************************************************/
    1776             : /*                            OpenMDArray()                             */
    1777             : /************************************************************************/
    1778             : 
    1779           3 : std::shared_ptr<GDALMDArray> GRIBGroup::OpenMDArray(const std::string &osName,
    1780             :                                                     CSLConstList) const
    1781             : {
    1782          12 :     for (const auto &array : m_poArrays)
    1783             :     {
    1784          11 :         if (array->GetName() == osName)
    1785           2 :             return array;
    1786             :     }
    1787           1 :     return nullptr;
    1788             : }
    1789             : 
    1790             : /************************************************************************/
    1791             : /*                             GRIBArray()                              */
    1792             : /************************************************************************/
    1793             : 
    1794          23 : GRIBArray::GRIBArray(const std::string &osName,
    1795          23 :                      const std::shared_ptr<GRIBSharedResource> &poShared)
    1796             :     : GDALAbstractMDArray("/", osName),
    1797          23 :       GDALPamMDArray("/", osName, poShared->GetPAM()), m_poShared(poShared)
    1798             : {
    1799          23 : }
    1800             : 
    1801             : /************************************************************************/
    1802             : /*                                Init()                                */
    1803             : /************************************************************************/
    1804             : 
    1805          23 : void GRIBArray::Init(GRIBGroup *poGroup, GRIBDataset *poDS,
    1806             :                      GRIBRasterBand *poBand, inventoryType *psInv)
    1807             : {
    1808          23 :     std::shared_ptr<GDALDimension> poDimX;
    1809          23 :     std::shared_ptr<GDALDimension> poDimY;
    1810             : 
    1811          23 :     GDALGeoTransform gt;
    1812          23 :     poDS->GetGeoTransform(gt);
    1813             : 
    1814          40 :     for (int i = 1; i <= poGroup->m_nHorizDimCounter; i++)
    1815             :     {
    1816          34 :         std::string osXLookup("X");
    1817          34 :         std::string osYLookup("Y");
    1818          34 :         if (i > 1)
    1819             :         {
    1820          15 :             osXLookup += CPLSPrintf("%d", i);
    1821          15 :             osYLookup += CPLSPrintf("%d", i);
    1822             :         }
    1823          34 :         auto oIterX = poGroup->m_oMapDims.find(osXLookup);
    1824          34 :         auto oIterY = poGroup->m_oMapDims.find(osYLookup);
    1825          34 :         CPLAssert(oIterX != poGroup->m_oMapDims.end());
    1826          34 :         CPLAssert(oIterY != poGroup->m_oMapDims.end());
    1827          34 :         if (oIterX->second->GetSize() ==
    1828          51 :                 static_cast<size_t>(poDS->GetRasterXSize()) &&
    1829          17 :             oIterY->second->GetSize() ==
    1830          17 :                 static_cast<size_t>(poDS->GetRasterYSize()))
    1831             :         {
    1832          17 :             bool bOK = true;
    1833          17 :             auto poVar = oIterX->second->GetIndexingVariable();
    1834          17 :             constexpr double EPSILON = 1e-10;
    1835          17 :             if (poVar)
    1836             :             {
    1837          17 :                 GUInt64 nStart = 0;
    1838          17 :                 size_t nCount = 1;
    1839          17 :                 double dfVal = 0;
    1840          17 :                 poVar->Read(&nStart, &nCount, nullptr, nullptr, m_dt, &dfVal);
    1841          17 :                 if (std::fabs(dfVal - (gt[0] + 0.5 * gt[1])) >
    1842          17 :                     EPSILON * std::fabs(dfVal))
    1843             :                 {
    1844           0 :                     bOK = false;
    1845             :                 }
    1846             :             }
    1847          17 :             if (bOK)
    1848             :             {
    1849          17 :                 poVar = oIterY->second->GetIndexingVariable();
    1850          17 :                 if (poVar)
    1851             :                 {
    1852          17 :                     GUInt64 nStart = 0;
    1853          17 :                     size_t nCount = 1;
    1854          17 :                     double dfVal = 0;
    1855          17 :                     poVar->Read(&nStart, &nCount, nullptr, nullptr, m_dt,
    1856             :                                 &dfVal);
    1857          17 :                     if (std::fabs(dfVal - (gt[3] + poDS->nRasterYSize * gt[5] -
    1858          17 :                                            0.5 * gt[5])) >
    1859          17 :                         EPSILON * std::fabs(dfVal))
    1860             :                     {
    1861           0 :                         bOK = false;
    1862             :                     }
    1863             :                 }
    1864             :             }
    1865          17 :             if (bOK)
    1866             :             {
    1867          17 :                 poDimX = oIterX->second;
    1868          17 :                 poDimY = oIterY->second;
    1869          17 :                 break;
    1870             :             }
    1871             :         }
    1872             :     }
    1873             : 
    1874          23 :     if (!poDimX || !poDimY)
    1875             :     {
    1876           6 :         poGroup->m_nHorizDimCounter++;
    1877             :         {
    1878          12 :             std::string osName("Y");
    1879           6 :             if (poGroup->m_nHorizDimCounter >= 2)
    1880             :             {
    1881           2 :                 osName = CPLSPrintf("Y%d", poGroup->m_nHorizDimCounter);
    1882             :             }
    1883             : 
    1884          12 :             poDimY = std::make_shared<GDALDimensionWeakIndexingVar>(
    1885           6 :                 poGroup->GetFullName(), osName, GDAL_DIM_TYPE_HORIZONTAL_Y,
    1886          18 :                 std::string(), poDS->GetRasterYSize());
    1887           6 :             poGroup->m_oMapDims[osName] = poDimY;
    1888           6 :             poGroup->m_dims.emplace_back(poDimY);
    1889             : 
    1890             :             auto var = GDALMDArrayRegularlySpaced::Create(
    1891             :                 "/", poDimY->GetName(), poDimY,
    1892          12 :                 gt[3] + poDS->GetRasterYSize() * gt[5], -gt[5], 0.5);
    1893           6 :             poDimY->SetIndexingVariable(var);
    1894           6 :             poGroup->AddArray(var);
    1895             :         }
    1896             : 
    1897             :         {
    1898          12 :             std::string osName("X");
    1899           6 :             if (poGroup->m_nHorizDimCounter >= 2)
    1900             :             {
    1901           2 :                 osName = CPLSPrintf("X%d", poGroup->m_nHorizDimCounter);
    1902             :             }
    1903             : 
    1904          12 :             poDimX = std::make_shared<GDALDimensionWeakIndexingVar>(
    1905           6 :                 poGroup->GetFullName(), osName, GDAL_DIM_TYPE_HORIZONTAL_X,
    1906          18 :                 std::string(), poDS->GetRasterXSize());
    1907           6 :             poGroup->m_oMapDims[osName] = poDimX;
    1908           6 :             poGroup->m_dims.emplace_back(poDimX);
    1909             : 
    1910             :             auto var = GDALMDArrayRegularlySpaced::Create(
    1911          12 :                 "/", poDimX->GetName(), poDimX, gt[0], gt[1], 0.5);
    1912           6 :             poDimX->SetIndexingVariable(var);
    1913           6 :             poGroup->AddArray(var);
    1914             :         }
    1915             :     }
    1916             : 
    1917          23 :     m_dims.emplace_back(poDimY);
    1918          23 :     m_dims.emplace_back(poDimX);
    1919          23 :     if (poDS->m_poSRS.get())
    1920             :     {
    1921          23 :         m_poSRS.reset(poDS->m_poSRS->Clone());
    1922          23 :         if (poDS->m_poSRS->GetDataAxisToSRSAxisMapping() ==
    1923          46 :             std::vector<int>{2, 1})
    1924          22 :             m_poSRS->SetDataAxisToSRSAxisMapping({1, 2});
    1925             :         else
    1926           1 :             m_poSRS->SetDataAxisToSRSAxisMapping({2, 1});
    1927             :     }
    1928             : 
    1929             :     const char *pszGribNormalizeUnits =
    1930          23 :         CPLGetConfigOption("GRIB_NORMALIZE_UNITS", "YES");
    1931          23 :     const bool bMetricUnits = CPLTestBool(pszGribNormalizeUnits);
    1932             : 
    1933          46 :     m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
    1934          46 :         GetFullName(), "name", psInv->element));
    1935          46 :     m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
    1936          23 :         GetFullName(), "long_name",
    1937          69 :         ConvertUnitInText(bMetricUnits, psInv->comment)));
    1938          23 :     m_osUnit = ConvertUnitInText(bMetricUnits, psInv->unitName);
    1939          23 :     if (!m_osUnit.empty() && m_osUnit[0] == '[' && m_osUnit.back() == ']')
    1940             :     {
    1941          23 :         m_osUnit = m_osUnit.substr(1, m_osUnit.size() - 2);
    1942             :     }
    1943          46 :     m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
    1944          46 :         GetFullName(), "first_level", psInv->shortFstLevel));
    1945             : 
    1946          23 :     if (poBand->m_nDisciplineCode >= 0)
    1947             :     {
    1948          14 :         m_attributes.emplace_back(std::make_shared<GDALAttributeNumeric>(
    1949          14 :             GetFullName(), "discipline_code", poBand->m_nDisciplineCode));
    1950             :     }
    1951          23 :     if (!poBand->m_osDisciplineName.empty())
    1952             :     {
    1953          14 :         m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
    1954          14 :             GetFullName(), "discipline_name", poBand->m_osDisciplineName));
    1955             :     }
    1956          23 :     if (poBand->m_nCenter >= 0)
    1957             :     {
    1958          14 :         m_attributes.emplace_back(std::make_shared<GDALAttributeNumeric>(
    1959          14 :             GetFullName(), "center_code", poBand->m_nCenter));
    1960             :     }
    1961          23 :     if (!poBand->m_osCenterName.empty())
    1962             :     {
    1963          14 :         m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
    1964          14 :             GetFullName(), "center_name", poBand->m_osCenterName));
    1965             :     }
    1966          23 :     if (poBand->m_nSubCenter >= 0)
    1967             :     {
    1968          12 :         m_attributes.emplace_back(std::make_shared<GDALAttributeNumeric>(
    1969          12 :             GetFullName(), "subcenter_code", poBand->m_nSubCenter));
    1970             :     }
    1971          23 :     if (!poBand->m_osSubCenterName.empty())
    1972             :     {
    1973           0 :         m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
    1974           0 :             GetFullName(), "subcenter_name", poBand->m_osSubCenterName));
    1975             :     }
    1976          23 :     if (!poBand->m_osSignRefTimeName.empty())
    1977             :     {
    1978          14 :         m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
    1979           7 :             GetFullName(), "signification_of_ref_time",
    1980          14 :             poBand->m_osSignRefTimeName));
    1981             :     }
    1982          23 :     if (!poBand->m_osRefTime.empty())
    1983             :     {
    1984          14 :         m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
    1985          14 :             GetFullName(), "reference_time_iso8601", poBand->m_osRefTime));
    1986             :     }
    1987          23 :     if (!poBand->m_osProductionStatus.empty())
    1988             :     {
    1989          14 :         m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
    1990          14 :             GetFullName(), "production_status", poBand->m_osProductionStatus));
    1991             :     }
    1992          23 :     if (!poBand->m_osType.empty())
    1993             :     {
    1994          14 :         m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
    1995          14 :             GetFullName(), "type", poBand->m_osType));
    1996             :     }
    1997          23 :     if (poBand->m_nPDTN >= 0)
    1998             :     {
    1999          14 :         m_attributes.emplace_back(std::make_shared<GDALAttributeNumeric>(
    2000           7 :             GetFullName(), "product_definition_template_number",
    2001          14 :             poBand->m_nPDTN));
    2002             :     }
    2003          23 :     if (!poBand->m_anPDSTemplateAssembledValues.empty())
    2004             :     {
    2005          14 :         m_attributes.emplace_back(std::make_shared<GDALAttributeNumeric>(
    2006           7 :             GetFullName(), "product_definition_numbers",
    2007          14 :             poBand->m_anPDSTemplateAssembledValues));
    2008             :     }
    2009             : 
    2010          23 :     int bHasNoData = FALSE;
    2011          23 :     double dfNoData = poBand->GetNoDataValue(&bHasNoData);
    2012          23 :     if (bHasNoData)
    2013             :     {
    2014          14 :         m_abyNoData.resize(sizeof(double));
    2015          14 :         memcpy(&m_abyNoData[0], &dfNoData, sizeof(double));
    2016             :     }
    2017          23 : }
    2018             : 
    2019             : /************************************************************************/
    2020             : /*                         ExtendTimeDim()                              */
    2021             : /************************************************************************/
    2022             : 
    2023          24 : void GRIBArray::ExtendTimeDim(vsi_l_offset nOffset, int subgNum,
    2024             :                               double dfValidTime)
    2025             : {
    2026          24 :     m_anOffsets.push_back(nOffset);
    2027          24 :     m_anSubgNums.push_back(subgNum);
    2028          24 :     m_adfTimes.push_back(dfValidTime);
    2029          24 : }
    2030             : 
    2031             : /************************************************************************/
    2032             : /*                           Finalize()                                 */
    2033             : /************************************************************************/
    2034             : 
    2035          23 : void GRIBArray::Finalize(GRIBGroup *poGroup, inventoryType *psInv)
    2036             : {
    2037          23 :     CPLAssert(!m_adfTimes.empty());
    2038          23 :     CPLAssert(m_anOffsets.size() == m_adfTimes.size());
    2039             : 
    2040          23 :     if (m_adfTimes.size() == 1)
    2041             :     {
    2042          44 :         m_attributes.emplace_back(std::make_shared<GDALAttributeNumeric>(
    2043          44 :             GetFullName(), "forecast_time", psInv->foreSec));
    2044          44 :         m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
    2045          44 :             GetFullName(), "forecast_time_unit", "sec"));
    2046          44 :         m_attributes.emplace_back(std::make_shared<GDALAttributeNumeric>(
    2047          44 :             GetFullName(), "reference_time", psInv->refTime));
    2048          44 :         m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
    2049          44 :             GetFullName(), "reference_time_unit", "sec UTC"));
    2050          44 :         m_attributes.emplace_back(std::make_shared<GDALAttributeNumeric>(
    2051          44 :             GetFullName(), "validity_time", m_adfTimes[0]));
    2052          44 :         m_attributes.emplace_back(std::make_shared<GDALAttributeString>(
    2053          44 :             GetFullName(), "validity_time_unit", "sec UTC"));
    2054          22 :         return;
    2055             :     }
    2056             : 
    2057           1 :     std::shared_ptr<GDALDimension> poDimTime;
    2058             : 
    2059           3 :     for (const auto &poDim : poGroup->m_dims)
    2060             :     {
    2061           2 :         if (STARTS_WITH(poDim->GetName().c_str(), "TIME") &&
    2062           0 :             poDim->GetSize() == m_adfTimes.size())
    2063             :         {
    2064           0 :             auto poVar = poDim->GetIndexingVariable();
    2065           0 :             if (poVar)
    2066             :             {
    2067           0 :                 GUInt64 nStart = 0;
    2068           0 :                 size_t nCount = 1;
    2069           0 :                 double dfStartTime = 0;
    2070           0 :                 poVar->Read(&nStart, &nCount, nullptr, nullptr, m_dt,
    2071             :                             &dfStartTime);
    2072           0 :                 if (dfStartTime == m_adfTimes[0])
    2073             :                 {
    2074           0 :                     poDimTime = poDim;
    2075           0 :                     break;
    2076             :                 }
    2077             :             }
    2078             :         }
    2079             :     }
    2080             : 
    2081           1 :     if (!poDimTime)
    2082             :     {
    2083           2 :         std::string osName("TIME");
    2084           1 :         int counter = 2;
    2085           1 :         while (poGroup->m_oMapDims.find(osName) != poGroup->m_oMapDims.end())
    2086             :         {
    2087           0 :             osName = CPLSPrintf("TIME%d", counter);
    2088           0 :             counter++;
    2089             :         }
    2090             : 
    2091           2 :         poDimTime = std::make_shared<GDALDimensionWeakIndexingVar>(
    2092           1 :             poGroup->GetFullName(), osName, GDAL_DIM_TYPE_TEMPORAL,
    2093           3 :             std::string(), m_adfTimes.size());
    2094           1 :         poGroup->m_oMapDims[osName] = poDimTime;
    2095           1 :         poGroup->m_dims.emplace_back(poDimTime);
    2096             : 
    2097           1 :         auto var = poGroup->m_memRootGroup->CreateMDArray(
    2098             :             poDimTime->GetName(),
    2099           4 :             std::vector<std::shared_ptr<GDALDimension>>{poDimTime},
    2100           3 :             GDALExtendedDataType::Create(GDT_Float64), nullptr);
    2101           1 :         poDimTime->SetIndexingVariable(var);
    2102           1 :         poGroup->AddArray(var);
    2103             : 
    2104           1 :         GUInt64 nStart = 0;
    2105           1 :         size_t nCount = m_adfTimes.size();
    2106           1 :         var->SetUnit("sec UTC");
    2107           1 :         const GUInt64 anStart[] = {nStart};
    2108           1 :         const size_t anCount[] = {nCount};
    2109           2 :         var->Write(anStart, anCount, nullptr, nullptr, var->GetDataType(),
    2110           1 :                    &m_adfTimes[0]);
    2111           1 :         auto attr = var->CreateAttribute("long_name", {},
    2112           3 :                                          GDALExtendedDataType::CreateString());
    2113           1 :         attr->Write("validity_time");
    2114             :     }
    2115             : 
    2116             : #if defined(__GNUC__)
    2117             : #pragma GCC diagnostic push
    2118             : #pragma GCC diagnostic ignored "-Wnull-dereference"
    2119             : #endif
    2120           1 :     m_dims.insert(m_dims.begin(), poDimTime);
    2121             : #if defined(__GNUC__)
    2122             : #pragma GCC diagnostic pop
    2123             : #endif
    2124           1 :     if (m_poSRS)
    2125             :     {
    2126           2 :         auto mapping = m_poSRS->GetDataAxisToSRSAxisMapping();
    2127           3 :         for (auto &v : mapping)
    2128           2 :             v += 1;
    2129           1 :         m_poSRS->SetDataAxisToSRSAxisMapping(mapping);
    2130             :     }
    2131             : }
    2132             : 
    2133             : /************************************************************************/
    2134             : /*                              LoadData()                              */
    2135             : /************************************************************************/
    2136             : 
    2137           6 : const std::vector<double> &GRIBSharedResource::LoadData(vsi_l_offset nOffset,
    2138             :                                                         int subgNum)
    2139             : {
    2140           6 :     if (m_nOffsetCurData == nOffset)
    2141           1 :         return m_adfCurData;
    2142             : 
    2143           5 :     grib_MetaData *metadata = nullptr;
    2144           5 :     double *data = nullptr;
    2145           5 :     GRIBRasterBand::ReadGribData(m_fp, nOffset, subgNum, &data, &metadata);
    2146           5 :     if (data == nullptr || metadata == nullptr)
    2147             :     {
    2148           0 :         if (metadata != nullptr)
    2149             :         {
    2150           0 :             MetaFree(metadata);
    2151           0 :             delete metadata;
    2152             :         }
    2153           0 :         free(data);
    2154           0 :         m_adfCurData.clear();
    2155           0 :         return m_adfCurData;
    2156             :     }
    2157           5 :     const int nx = metadata->gds.Nx;
    2158           5 :     const int ny = metadata->gds.Ny;
    2159           5 :     MetaFree(metadata);
    2160           5 :     delete metadata;
    2161           5 :     if (nx <= 0 || ny <= 0)
    2162             :     {
    2163           0 :         free(data);
    2164           0 :         m_adfCurData.clear();
    2165           0 :         return m_adfCurData;
    2166             :     }
    2167           5 :     const size_t nPointCount = static_cast<size_t>(nx) * ny;
    2168           5 :     const size_t nByteCount = nPointCount * sizeof(double);
    2169             : #ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
    2170             :     if (nByteCount > static_cast<size_t>(INT_MAX))
    2171             :     {
    2172             :         CPLError(CE_Failure, CPLE_OutOfMemory,
    2173             :                  "Too large memory allocation attempt");
    2174             :         free(data);
    2175             :         m_adfCurData.clear();
    2176             :         return m_adfCurData;
    2177             :     }
    2178             : #endif
    2179             :     try
    2180             :     {
    2181           5 :         m_adfCurData.resize(nPointCount);
    2182             :     }
    2183           0 :     catch (const std::exception &e)
    2184             :     {
    2185           0 :         CPLError(CE_Failure, CPLE_OutOfMemory, "%s", e.what());
    2186           0 :         free(data);
    2187           0 :         m_adfCurData.clear();
    2188           0 :         return m_adfCurData;
    2189             :     }
    2190           5 :     m_nOffsetCurData = nOffset;
    2191           5 :     memcpy(&m_adfCurData[0], data, nByteCount);
    2192           5 :     free(data);
    2193           5 :     return m_adfCurData;
    2194             : }
    2195             : 
    2196             : /************************************************************************/
    2197             : /*                             IRead()                                  */
    2198             : /************************************************************************/
    2199             : 
    2200           4 : bool GRIBArray::IRead(const GUInt64 *arrayStartIdx, const size_t *count,
    2201             :                       const GInt64 *arrayStep, const GPtrDiff_t *bufferStride,
    2202             :                       const GDALExtendedDataType &bufferDataType,
    2203             :                       void *pDstBuffer) const
    2204             : {
    2205           4 :     const size_t nBufferDTSize(bufferDataType.GetSize());
    2206           4 :     if (m_dims.size() == 2)
    2207             :     {
    2208           2 :         auto &vals = m_poShared->LoadData(m_anOffsets[0], m_anSubgNums[0]);
    2209           2 :         constexpr int Y_IDX = 0;
    2210           2 :         constexpr int X_IDX = 1;
    2211           4 :         if (vals.empty() ||
    2212           2 :             vals.size() != m_dims[Y_IDX]->GetSize() * m_dims[X_IDX]->GetSize())
    2213           0 :             return false;
    2214           2 :         const size_t nWidth(static_cast<size_t>(m_dims[X_IDX]->GetSize()));
    2215           2 :         const bool bDirectCopy = m_dt == bufferDataType &&
    2216           3 :                                  arrayStep[X_IDX] == 1 &&
    2217           1 :                                  bufferStride[X_IDX] == 1;
    2218         150 :         for (size_t j = 0; j < count[Y_IDX]; j++)
    2219             :         {
    2220         148 :             const size_t y = static_cast<size_t>(arrayStartIdx[Y_IDX] +
    2221         148 :                                                  j * arrayStep[Y_IDX]);
    2222         148 :             GByte *pabyDstPtr = static_cast<GByte *>(pDstBuffer) +
    2223         148 :                                 j * bufferStride[Y_IDX] * nBufferDTSize;
    2224         148 :             const size_t x = static_cast<size_t>(arrayStartIdx[X_IDX]);
    2225         148 :             const double *srcPtr = &vals[y * nWidth + x];
    2226         148 :             if (bDirectCopy)
    2227             :             {
    2228          74 :                 memcpy(pabyDstPtr, srcPtr, count[X_IDX] * sizeof(double));
    2229             :             }
    2230             :             else
    2231             :             {
    2232          74 :                 const auto dstPtrInc = bufferStride[X_IDX] * nBufferDTSize;
    2233        4958 :                 for (size_t i = 0; i < count[X_IDX]; i++)
    2234             :                 {
    2235        4884 :                     GDALExtendedDataType::CopyValue(srcPtr, m_dt, pabyDstPtr,
    2236             :                                                     bufferDataType);
    2237        4884 :                     srcPtr += static_cast<std::ptrdiff_t>(arrayStep[X_IDX]);
    2238        4884 :                     pabyDstPtr += dstPtrInc;
    2239             :                 }
    2240             :             }
    2241             :         }
    2242           2 :         return true;
    2243             :     }
    2244             : 
    2245           2 :     constexpr int T_IDX = 0;
    2246           2 :     constexpr int Y_IDX = 1;
    2247           2 :     constexpr int X_IDX = 2;
    2248           2 :     const size_t nWidth(static_cast<size_t>(m_dims[X_IDX]->GetSize()));
    2249           3 :     const bool bDirectCopy = m_dt == bufferDataType && arrayStep[X_IDX] == 1 &&
    2250           1 :                              bufferStride[X_IDX] == 1;
    2251           6 :     for (size_t k = 0; k < count[T_IDX]; k++)
    2252             :     {
    2253           4 :         const size_t tIdx =
    2254           4 :             static_cast<size_t>(arrayStartIdx[T_IDX] + k * arrayStep[T_IDX]);
    2255           4 :         CPLAssert(tIdx < m_anOffsets.size());
    2256             :         auto &vals =
    2257           4 :             m_poShared->LoadData(m_anOffsets[tIdx], m_anSubgNums[tIdx]);
    2258           8 :         if (vals.empty() ||
    2259           4 :             vals.size() != m_dims[Y_IDX]->GetSize() * m_dims[X_IDX]->GetSize())
    2260           0 :             return false;
    2261         520 :         for (size_t j = 0; j < count[Y_IDX]; j++)
    2262             :         {
    2263         516 :             const size_t y = static_cast<size_t>(arrayStartIdx[Y_IDX] +
    2264         516 :                                                  j * arrayStep[Y_IDX]);
    2265         516 :             GByte *pabyDstPtr =
    2266             :                 static_cast<GByte *>(pDstBuffer) +
    2267         516 :                 (k * bufferStride[T_IDX] + j * bufferStride[Y_IDX]) *
    2268             :                     nBufferDTSize;
    2269         516 :             const size_t x = static_cast<size_t>(arrayStartIdx[X_IDX]);
    2270         516 :             const double *srcPtr = &vals[y * nWidth + x];
    2271         516 :             if (bDirectCopy)
    2272             :             {
    2273         258 :                 memcpy(pabyDstPtr, srcPtr, count[X_IDX] * sizeof(double));
    2274             :             }
    2275             :             else
    2276             :             {
    2277         258 :                 const auto dstPtrInc = bufferStride[X_IDX] * nBufferDTSize;
    2278       45924 :                 for (size_t i = 0; i < count[X_IDX]; i++)
    2279             :                 {
    2280       45666 :                     GDALExtendedDataType::CopyValue(srcPtr, m_dt, pabyDstPtr,
    2281             :                                                     bufferDataType);
    2282       45666 :                     srcPtr += static_cast<std::ptrdiff_t>(arrayStep[X_IDX]);
    2283       45666 :                     pabyDstPtr += dstPtrInc;
    2284             :                 }
    2285             :             }
    2286             :         }
    2287             :     }
    2288             : 
    2289           2 :     return true;
    2290             : }
    2291             : 
    2292             : /************************************************************************/
    2293             : /*                          OpenMultiDim()                              */
    2294             : /************************************************************************/
    2295             : 
    2296           4 : GDALDataset *GRIBDataset::OpenMultiDim(GDALOpenInfo *poOpenInfo)
    2297             : 
    2298             : {
    2299             :     auto poShared = std::make_shared<GRIBSharedResource>(
    2300           8 :         poOpenInfo->pszFilename, poOpenInfo->fpL);
    2301           8 :     auto poRootGroup = std::make_shared<GRIBGroup>(poShared);
    2302           4 :     poOpenInfo->fpL = nullptr;
    2303             : 
    2304           4 :     VSIFSeekL(poShared->m_fp, 0, SEEK_SET);
    2305             : 
    2306             :     // Contains an GRIB2 message inventory of the file.
    2307             :     // We can't use the potential .idx file
    2308           8 :     auto pInventories = std::make_unique<InventoryWrapperGrib>(poShared->m_fp);
    2309             : 
    2310           4 :     if (pInventories->result() <= 0)
    2311             :     {
    2312           0 :         char *errMsg = errSprintf(nullptr);
    2313           0 :         if (errMsg != nullptr)
    2314           0 :             CPLDebug("GRIB", "%s", errMsg);
    2315           0 :         free(errMsg);
    2316             : 
    2317           0 :         CPLError(CE_Failure, CPLE_OpenFailed,
    2318             :                  "%s is a grib file, "
    2319             :                  "but no raster dataset was successfully identified.",
    2320             :                  poOpenInfo->pszFilename);
    2321           0 :         return nullptr;
    2322             :     }
    2323             : 
    2324             :     // Create band objects.
    2325           4 :     GRIBDataset *poDS = new GRIBDataset();
    2326           4 :     poDS->fp = poShared->m_fp;
    2327             : 
    2328           4 :     std::shared_ptr<GRIBArray> poArray;
    2329           8 :     std::set<std::string> oSetArrayNames;
    2330           8 :     std::string osElement;
    2331           8 :     std::string osShortFstLevel;
    2332           4 :     double dfRefTime = 0;
    2333           4 :     double dfForecastTime = 0;
    2334          28 :     for (uInt4 i = 0; i < pInventories->length(); ++i)
    2335             :     {
    2336          24 :         inventoryType *psInv = pInventories->get(i);
    2337          24 :         uInt4 bandNr = i + 1;
    2338          24 :         assert(bandNr <= 65536);
    2339             : 
    2340             :         // GRIB messages can be preceded by "garbage". GRIB2Inventory()
    2341             :         // does not return the offset to the real start of the message
    2342             :         GByte abyHeader[1024 + 1];
    2343          24 :         VSIFSeekL(poShared->m_fp, psInv->start, SEEK_SET);
    2344             :         size_t nRead =
    2345          24 :             VSIFReadL(abyHeader, 1, sizeof(abyHeader) - 1, poShared->m_fp);
    2346          24 :         abyHeader[nRead] = 0;
    2347             :         // Find the real offset of the fist message
    2348          24 :         const char *pasHeader = reinterpret_cast<char *>(abyHeader);
    2349          24 :         int nOffsetFirstMessage = 0;
    2350         144 :         for (int j = 0; j < poOpenInfo->nHeaderBytes - 3; j++)
    2351             :         {
    2352         144 :             if (STARTS_WITH_CI(pasHeader + j, "GRIB")
    2353             : #ifdef ENABLE_TDLP
    2354             :                 || STARTS_WITH_CI(pasHeader + j, "TDLP")
    2355             : #endif
    2356             :             )
    2357             :             {
    2358          24 :                 nOffsetFirstMessage = j;
    2359          24 :                 break;
    2360             :             }
    2361             :         }
    2362          24 :         psInv->start += nOffsetFirstMessage;
    2363             : 
    2364          46 :         if (poArray == nullptr || osElement != psInv->element ||
    2365          46 :             osShortFstLevel != psInv->shortFstLevel ||
    2366           1 :             !((dfRefTime == psInv->refTime &&
    2367           1 :                dfForecastTime != psInv->foreSec) ||
    2368           0 :               (dfRefTime != psInv->refTime &&
    2369           0 :                dfForecastTime == psInv->foreSec)))
    2370             :         {
    2371          23 :             if (poArray)
    2372             :             {
    2373          19 :                 poArray->Finalize(poRootGroup.get(), pInventories->get(i - 1));
    2374          19 :                 poRootGroup->AddArray(poArray);
    2375             :             }
    2376             : 
    2377          23 :             grib_MetaData *metaData = nullptr;
    2378          23 :             GRIBRasterBand::ReadGribData(poShared->m_fp, psInv->start,
    2379          23 :                                          psInv->subgNum, nullptr, &metaData);
    2380          23 :             if (metaData == nullptr || metaData->gds.Nx < 1 ||
    2381          23 :                 metaData->gds.Ny < 1)
    2382             :             {
    2383           0 :                 CPLError(CE_Failure, CPLE_OpenFailed,
    2384             :                          "%s is a grib file, "
    2385             :                          "but no raster dataset was successfully identified.",
    2386             :                          poOpenInfo->pszFilename);
    2387             : 
    2388           0 :                 if (metaData != nullptr)
    2389             :                 {
    2390           0 :                     MetaFree(metaData);
    2391           0 :                     delete metaData;
    2392             :                 }
    2393           0 :                 poDS->fp = nullptr;
    2394           0 :                 CPLReleaseMutex(hGRIBMutex);
    2395           0 :                 delete poDS;
    2396           0 :                 CPLAcquireMutex(hGRIBMutex, 1000.0);
    2397           0 :                 return nullptr;
    2398             :             }
    2399          23 :             psInv->GribVersion = metaData->GribVersion;
    2400             : 
    2401             :             // Set the DataSet's x,y size, georeference and projection from
    2402             :             // the first GRIB band.
    2403          23 :             poDS->SetGribMetaData(metaData);
    2404             : 
    2405          46 :             GRIBRasterBand gribBand(poDS, bandNr, psInv);
    2406          23 :             if (psInv->GribVersion == 2)
    2407           7 :                 gribBand.FindPDSTemplateGRIB2();
    2408          23 :             osElement = psInv->element;
    2409          23 :             osShortFstLevel = psInv->shortFstLevel;
    2410          23 :             dfRefTime = psInv->refTime;
    2411          23 :             dfForecastTime = psInv->foreSec;
    2412             : 
    2413          46 :             std::string osRadix(osElement + '_' + osShortFstLevel);
    2414          46 :             std::string osName(osRadix);
    2415          23 :             int nCounter = 2;
    2416          23 :             while (oSetArrayNames.find(osName) != oSetArrayNames.end())
    2417             :             {
    2418           0 :                 osName = osRadix + CPLSPrintf("_%d", nCounter);
    2419           0 :                 nCounter++;
    2420             :             }
    2421          23 :             oSetArrayNames.insert(osName);
    2422          23 :             poArray = GRIBArray::Create(osName, poShared);
    2423          23 :             poArray->Init(poRootGroup.get(), poDS, &gribBand, psInv);
    2424          23 :             poArray->ExtendTimeDim(psInv->start, psInv->subgNum,
    2425             :                                    psInv->validTime);
    2426             : 
    2427          23 :             MetaFree(metaData);
    2428          23 :             delete metaData;
    2429             :         }
    2430             :         else
    2431             :         {
    2432           1 :             poArray->ExtendTimeDim(psInv->start, psInv->subgNum,
    2433             :                                    psInv->validTime);
    2434             :         }
    2435             :     }
    2436             : 
    2437           4 :     if (poArray)
    2438             :     {
    2439           8 :         poArray->Finalize(poRootGroup.get(),
    2440           4 :                           pInventories->get(pInventories->length() - 1));
    2441           4 :         poRootGroup->AddArray(poArray);
    2442             :     }
    2443             : 
    2444           4 :     poDS->fp = nullptr;
    2445           4 :     poDS->m_poRootGroup = poRootGroup;
    2446             : 
    2447           4 :     poDS->SetDescription(poOpenInfo->pszFilename);
    2448             : 
    2449             :     // Release hGRIBMutex otherwise we'll deadlock with GDALDataset own
    2450             :     // hGRIBMutex.
    2451           4 :     CPLReleaseMutex(hGRIBMutex);
    2452           4 :     poDS->TryLoadXML();
    2453           4 :     CPLAcquireMutex(hGRIBMutex, 1000.0);
    2454             : 
    2455           4 :     return poDS;
    2456             : }
    2457             : 
    2458             : /************************************************************************/
    2459             : /*                            SetMetadata()                             */
    2460             : /************************************************************************/
    2461             : 
    2462         326 : void GRIBDataset::SetGribMetaData(grib_MetaData *meta)
    2463             : {
    2464         326 :     nRasterXSize = meta->gds.Nx;
    2465         326 :     nRasterYSize = meta->gds.Ny;
    2466             : 
    2467             :     // Image projection.
    2468         326 :     OGRSpatialReference oSRS;
    2469         326 :     oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    2470             : 
    2471         326 :     switch (meta->gds.projType)
    2472             :     {
    2473         134 :         case GS3_LATLON:
    2474             :         case GS3_GAUSSIAN_LATLON:
    2475             :             // No projection, only latlon system (geographic).
    2476         134 :             break;
    2477           2 :         case GS3_ROTATED_LATLON:
    2478             :             // will apply pole rotation afterwards
    2479           2 :             break;
    2480          26 :         case GS3_MERCATOR:
    2481          26 :             if (meta->gds.orientLon == 0.0)
    2482             :             {
    2483          26 :                 if (meta->gds.meshLat == 0.0)
    2484           5 :                     oSRS.SetMercator(0.0, 0.0, 1.0, 0.0, 0.0);
    2485             :                 else
    2486          21 :                     oSRS.SetMercator2SP(meta->gds.meshLat, 0.0, 0.0, 0.0, 0.0);
    2487             :             }
    2488             :             else
    2489             :             {
    2490           0 :                 CPLError(CE_Warning, CPLE_NotSupported,
    2491             :                          "Orientation of the grid != 0 not supported");
    2492           0 :                 return;
    2493             :             }
    2494          26 :             break;
    2495         135 :         case GS3_TRANSVERSE_MERCATOR:
    2496         135 :             oSRS.SetTM(meta->gds.latitude_of_origin,
    2497             :                        Lon360to180(meta->gds.central_meridian),
    2498         135 :                        std::abs(meta->gds.scaleLat1 - 0.9996) < 1e-8
    2499             :                            ? 0.9996
    2500             :                            : meta->gds.scaleLat1,
    2501             :                        meta->gds.x0, meta->gds.y0);
    2502         135 :             break;
    2503           9 :         case GS3_POLAR:
    2504           9 :             oSRS.SetPS(meta->gds.meshLat, meta->gds.orientLon, 1.0, 0.0, 0.0);
    2505           9 :             break;
    2506           9 :         case GS3_LAMBERT:
    2507           9 :             oSRS.SetLCC(meta->gds.scaleLat1, meta->gds.scaleLat2,
    2508             :                         meta->gds.meshLat, Lon360to180(meta->gds.orientLon),
    2509             :                         0.0, 0.0);
    2510           9 :             break;
    2511           5 :         case GS3_ALBERS_EQUAL_AREA:
    2512           5 :             oSRS.SetACEA(meta->gds.scaleLat1, meta->gds.scaleLat2,
    2513             :                          meta->gds.meshLat, Lon360to180(meta->gds.orientLon),
    2514             :                          0.0, 0.0);
    2515           5 :             break;
    2516             : 
    2517           0 :         case GS3_ORTHOGRAPHIC:
    2518             : 
    2519             :             // oSRS.SetOrthographic( 0.0, meta->gds.orientLon,
    2520             :             //                       meta->gds.lon2, meta->gds.lat2);
    2521             : 
    2522             :             // oSRS.SetGEOS( meta->gds.orientLon, meta->gds.stretchFactor,
    2523             :             //               meta->gds.lon2, meta->gds.lat2);
    2524             : 
    2525             :             // TODO: Hardcoded for now. How to parse the meta->gds section?
    2526           0 :             oSRS.SetGEOS(0, 35785831, 0, 0);
    2527           0 :             break;
    2528           6 :         case GS3_LAMBERT_AZIMUTHAL:
    2529           6 :             oSRS.SetLAEA(meta->gds.meshLat, Lon360to180(meta->gds.orientLon),
    2530             :                          0.0, 0.0);
    2531           6 :             break;
    2532             : 
    2533           0 :         case GS3_EQUATOR_EQUIDIST:
    2534           0 :             break;
    2535           0 :         case GS3_AZIMUTH_RANGE:
    2536           0 :             break;
    2537             :     }
    2538             : 
    2539         326 :     if (oSRS.IsProjected())
    2540             :     {
    2541         190 :         oSRS.SetLinearUnits("Metre", 1.0);
    2542             :     }
    2543             : 
    2544         326 :     const bool bHaveEarthModel =
    2545         326 :         meta->gds.majEarth > 0.0 && meta->gds.minEarth > 0.0;
    2546             :     // In meters.
    2547             :     const double a = bHaveEarthModel
    2548         326 :                          ? meta->gds.majEarth * 1.0e3
    2549           2 :                          : CPLAtof(CPLGetConfigOption("GRIB_DEFAULT_SEMI_MAJOR",
    2550         326 :                                                       "6377563.396"));
    2551             :     const double b =
    2552             :         bHaveEarthModel
    2553         328 :             ? meta->gds.minEarth * 1.0e3
    2554           2 :             : (meta->gds.f_sphere
    2555           2 :                    ? a
    2556           0 :                    : CPLAtof(CPLGetConfigOption("GRIB_DEFAULT_SEMI_MINOR",
    2557         326 :                                                 "6356256.910")));
    2558         326 :     if (meta->gds.majEarth == 0 || meta->gds.minEarth == 0)
    2559             :     {
    2560           0 :         CPLDebug("GRIB", "No earth model. Assuming a=%f and b=%f", a, b);
    2561             :     }
    2562         326 :     else if (meta->gds.majEarth < 0 || meta->gds.minEarth < 0)
    2563             :     {
    2564             :         const char *pszUseDefaultSpheroid =
    2565           2 :             CPLGetConfigOption("GRIB_USE_DEFAULT_SPHEROID", nullptr);
    2566           2 :         if (!pszUseDefaultSpheroid)
    2567             :         {
    2568           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    2569             :                      "The GRIB file contains invalid values for the spheroid. "
    2570             :                      "You may set the GRIB_USE_DEFAULT_SPHEROID configuration "
    2571             :                      "option to YES to use a default spheroid with "
    2572             :                      "a=%f and b=%f",
    2573             :                      a, b);
    2574           1 :             return;
    2575             :         }
    2576           1 :         else if (!CPLTestBool(pszUseDefaultSpheroid))
    2577             :         {
    2578           0 :             return;
    2579             :         }
    2580           1 :         CPLDebug("GRIB", "Invalid earth model. Assuming a=%f and b=%f", a, b);
    2581             :     }
    2582             : 
    2583         325 :     if (meta->gds.f_sphere || (a == b))
    2584             :     {
    2585          88 :         oSRS.SetGeogCS("Coordinate System imported from GRIB file", nullptr,
    2586             :                        "Sphere", a, 0.0);
    2587             :     }
    2588             :     else
    2589             :     {
    2590         237 :         const double fInv = a / (a - b);
    2591         339 :         if (std::abs(a - 6378137.0) < 0.01 &&
    2592         102 :             std::abs(fInv - 298.257223563) < 1e-9)  // WGS84
    2593             :         {
    2594         100 :             if (meta->gds.projType == GS3_LATLON)
    2595          68 :                 oSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
    2596             :             else
    2597             :             {
    2598          32 :                 oSRS.SetGeogCS("Coordinate System imported from GRIB file",
    2599             :                                "WGS_1984", "WGS 84", 6378137., 298.257223563);
    2600             :             }
    2601             :         }
    2602         139 :         else if (std::abs(a - 6378137.0) < 0.01 &&
    2603           2 :                  std::abs(fInv - 298.257222101) < 1e-9)  // GRS80
    2604             :         {
    2605           2 :             oSRS.SetGeogCS("Coordinate System imported from GRIB file", nullptr,
    2606             :                            "GRS80", 6378137., 298.257222101);
    2607             :         }
    2608             :         else
    2609             :         {
    2610         135 :             oSRS.SetGeogCS("Coordinate System imported from GRIB file", nullptr,
    2611             :                            "Spheroid imported from GRIB file", a, fInv);
    2612             :         }
    2613             :     }
    2614             : 
    2615         325 :     if (meta->gds.projType == GS3_ROTATED_LATLON)
    2616             :     {
    2617           2 :         oSRS.SetDerivedGeogCRSWithPoleRotationGRIBConvention(
    2618             :             oSRS.GetName(), meta->gds.southLat, Lon360to180(meta->gds.southLon),
    2619             :             meta->gds.angleRotate);
    2620             :     }
    2621             : 
    2622         325 :     OGRSpatialReference oLL;  // Construct the "geographic" part of oSRS.
    2623         325 :     oLL.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    2624         325 :     oLL.CopyGeogCSFrom(&oSRS);
    2625             : 
    2626         325 :     double rMinX = 0.0;
    2627         325 :     double rMaxY = 0.0;
    2628         325 :     double rPixelSizeX = 0.0;
    2629         325 :     double rPixelSizeY = 0.0;
    2630         325 :     bool bError = false;
    2631         325 :     if (meta->gds.projType == GS3_ORTHOGRAPHIC)
    2632             :     {
    2633             :         // This is what should work, but it doesn't. Dx seems to have an
    2634             :         // inverse relation with pixel size.
    2635             :         // rMinX = -meta->gds.Dx * (meta->gds.Nx / 2);
    2636             :         // rMaxY = meta->gds.Dy * (meta->gds.Ny / 2);
    2637             :         // Hardcoded for now, assumption: GEOS projection, full disc (like MSG).
    2638           0 :         const double geosExtentInMeters = 11137496.552;
    2639           0 :         rMinX = -(geosExtentInMeters / 2);
    2640           0 :         rMaxY = geosExtentInMeters / 2;
    2641           0 :         rPixelSizeX = geosExtentInMeters / meta->gds.Nx;
    2642           0 :         rPixelSizeY = geosExtentInMeters / meta->gds.Ny;
    2643             :     }
    2644         325 :     else if (meta->gds.projType == GS3_TRANSVERSE_MERCATOR)
    2645             :     {
    2646         134 :         rMinX = meta->gds.x1;
    2647         134 :         rMaxY = meta->gds.y2;
    2648         134 :         rPixelSizeX = meta->gds.Dx;
    2649         134 :         rPixelSizeY = meta->gds.Dy;
    2650             :     }
    2651         191 :     else if (oSRS.IsProjected() && meta->gds.projType != GS3_ROTATED_LATLON)
    2652             :     {
    2653             :         // Longitude in degrees, to be transformed to meters (or degrees in
    2654             :         // case of latlon).
    2655          55 :         rMinX = Lon360to180(meta->gds.lon1);
    2656             :         // Latitude in degrees, to be transformed to meters.
    2657          55 :         double dfGridOriY = meta->gds.lat1;
    2658             : 
    2659          55 :         if (m_poSRS == nullptr || m_poLL == nullptr ||
    2660          55 :             !m_poSRS->IsSame(&oSRS) || !m_poLL->IsSame(&oLL))
    2661             :         {
    2662         110 :             m_poCT = std::unique_ptr<OGRCoordinateTransformation>(
    2663          55 :                 OGRCreateCoordinateTransformation(&(oLL), &(oSRS)));
    2664             :         }
    2665             : 
    2666             :         // Transform it to meters.
    2667          55 :         if ((m_poCT != nullptr) && m_poCT->Transform(1, &rMinX, &dfGridOriY))
    2668             :         {
    2669          55 :             if (meta->gds.scan == GRIB2BIT_2)  // Y is minY, GDAL wants maxY.
    2670             :             {
    2671          55 :                 const char *pszConfigOpt = CPLGetConfigOption(
    2672             :                     "GRIB_LATITUDE_OF_FIRST_GRID_POINT_IS_SOUTHERN_MOST",
    2673             :                     nullptr);
    2674             :                 bool bLatOfFirstPointIsSouthernMost =
    2675          55 :                     !pszConfigOpt || CPLTestBool(pszConfigOpt);
    2676             : 
    2677             :                 // Hack for a file called MANAL_2023030103.grb2 that
    2678             :                 // uses LCC and has Latitude of false origin = 30
    2679             :                 // Longitude of false origin = 140
    2680             :                 // Latitude of 1st standard parallel = 60
    2681             :                 // Latitude of 2nd standard parallel = 30
    2682             :                 // but whose (meta->gds.lon1, meta->gds.lat1) qualifies the
    2683             :                 // northern-most point of the grid and not the bottom-most one
    2684             :                 // as it should given the scan == GRIB2BIT_2
    2685          55 :                 if (!pszConfigOpt && meta->gds.projType == GS3_LAMBERT &&
    2686           9 :                     std::fabs(meta->gds.scaleLat1 - 60) <= 1e-8 &&
    2687           1 :                     std::fabs(meta->gds.scaleLat2 - 30) <= 1e-8 &&
    2688         111 :                     std::fabs(meta->gds.meshLat - 30) <= 1e-8 &&
    2689           1 :                     std::fabs(Lon360to180(meta->gds.orientLon) - 140) <= 1e-8)
    2690             :                 {
    2691           1 :                     double dfXCenterProj = Lon360to180(meta->gds.orientLon);
    2692           1 :                     double dfYCenterProj = meta->gds.meshLat;
    2693           1 :                     if (m_poCT->Transform(1, &dfXCenterProj, &dfYCenterProj))
    2694             :                     {
    2695           1 :                         double dfXCenterGridNominal =
    2696           1 :                             rMinX + nRasterXSize * meta->gds.Dx / 2;
    2697           1 :                         double dfYCenterGridNominal =
    2698           1 :                             dfGridOriY + nRasterYSize * meta->gds.Dy / 2;
    2699           1 :                         double dfXCenterGridBuggy = dfXCenterGridNominal;
    2700           1 :                         double dfYCenterGridBuggy =
    2701           1 :                             dfGridOriY - nRasterYSize * meta->gds.Dy / 2;
    2702           5 :                         const auto SQR = [](double x) { return x * x; };
    2703           1 :                         if (SQR(dfXCenterGridBuggy - dfXCenterProj) +
    2704           1 :                                 SQR(dfYCenterGridBuggy - dfYCenterProj) <
    2705           1 :                             SQR(10) *
    2706           2 :                                 (SQR(dfXCenterGridNominal - dfXCenterProj) +
    2707           1 :                                  SQR(dfYCenterGridNominal - dfYCenterProj)))
    2708             :                         {
    2709           1 :                             CPLError(
    2710             :                                 CE_Warning, CPLE_AppDefined,
    2711             :                                 "Likely buggy grid registration for GRIB2 "
    2712             :                                 "product: heuristics shows that the "
    2713             :                                 "latitudeOfFirstGridPoint is likely to qualify "
    2714             :                                 "the latitude of the northern-most grid point "
    2715             :                                 "instead of the southern-most grid point as "
    2716             :                                 "expected. Please report to data producer. "
    2717             :                                 "This heuristics can be disabled by setting "
    2718             :                                 "the "
    2719             :                                 "GRIB_LATITUDE_OF_FIRST_GRID_POINT_IS_SOUTHERN_"
    2720             :                                 "MOST configuration option to YES.");
    2721           1 :                             bLatOfFirstPointIsSouthernMost = false;
    2722             :                         }
    2723             :                     }
    2724             :                 }
    2725          55 :                 if (bLatOfFirstPointIsSouthernMost)
    2726             :                 {
    2727             :                     // -1 because we GDAL needs the coordinates of the centre of
    2728             :                     // the pixel.
    2729          54 :                     rMaxY = dfGridOriY + (meta->gds.Ny - 1) * meta->gds.Dy;
    2730             :                 }
    2731             :                 else
    2732             :                 {
    2733           1 :                     rMaxY = dfGridOriY;
    2734             :                 }
    2735             :             }
    2736             :             else
    2737             :             {
    2738           0 :                 rMaxY = dfGridOriY;
    2739             :             }
    2740          55 :             rPixelSizeX = meta->gds.Dx;
    2741          55 :             rPixelSizeY = meta->gds.Dy;
    2742             :         }
    2743             :         else
    2744             :         {
    2745           0 :             rMinX = 0.0;
    2746             :             // rMaxY = 0.0;
    2747             : 
    2748           0 :             rPixelSizeX = 1.0;
    2749           0 :             rPixelSizeY = -1.0;
    2750             : 
    2751           0 :             bError = true;
    2752           0 :             CPLError(CE_Warning, CPLE_AppDefined,
    2753             :                      "Unable to perform coordinate transformations, so the "
    2754             :                      "correct projected geotransform could not be deduced "
    2755             :                      "from the lat/long control points.  "
    2756             :                      "Defaulting to ungeoreferenced.");
    2757             :         }
    2758             :     }
    2759             :     else
    2760             :     {
    2761             :         // Longitude in degrees, to be transformed to meters (or degrees in
    2762             :         // case of latlon).
    2763         136 :         rMinX = meta->gds.lon1;
    2764             :         // Latitude in degrees, to be transformed to meters.
    2765         136 :         rMaxY = meta->gds.lat1;
    2766             : 
    2767         136 :         double rMinY = meta->gds.lat2;
    2768         136 :         double rMaxX = meta->gds.lon2;
    2769         136 :         if (meta->gds.lat2 > rMaxY)
    2770             :         {
    2771         105 :             rMaxY = meta->gds.lat2;
    2772         105 :             rMinY = meta->gds.lat1;
    2773             :         }
    2774             : 
    2775         136 :         if (meta->gds.Nx == 1)
    2776          16 :             rPixelSizeX = meta->gds.Dx;
    2777         120 :         else if (meta->gds.lon1 > meta->gds.lon2)
    2778          13 :             rPixelSizeX = (360.0 - (meta->gds.lon1 - meta->gds.lon2)) /
    2779          13 :                           (meta->gds.Nx - 1);
    2780             :         else
    2781         107 :             rPixelSizeX =
    2782         107 :                 (meta->gds.lon2 - meta->gds.lon1) / (meta->gds.Nx - 1);
    2783             : 
    2784         136 :         if (meta->gds.Ny == 1)
    2785          17 :             rPixelSizeY = meta->gds.Dy;
    2786             :         else
    2787         119 :             rPixelSizeY = (rMaxY - rMinY) / (meta->gds.Ny - 1);
    2788             : 
    2789             :         // Do some sanity checks for cases that can't be handled by the above
    2790             :         // pixel size corrections. GRIB1 has a minimum precision of 0.001
    2791             :         // for latitudes and longitudes, so we'll allow a bit higher than that.
    2792         136 :         if (rPixelSizeX < 0 || fabs(rPixelSizeX - meta->gds.Dx) > 0.002)
    2793           0 :             rPixelSizeX = meta->gds.Dx;
    2794             : 
    2795         136 :         if (rPixelSizeY < 0 || fabs(rPixelSizeY - meta->gds.Dy) > 0.002)
    2796           0 :             rPixelSizeY = meta->gds.Dy;
    2797             : 
    2798             :         // GRIB2 files have longitudes in the [0-360] range
    2799             :         // Shift them to the traditional [-180,180] longitude range
    2800             :         // See https://github.com/OSGeo/gdal/issues/4524
    2801         194 :         if ((rMinX + rPixelSizeX >= 180 || rMaxX - rPixelSizeX >= 180) &&
    2802          58 :             CPLTestBool(
    2803             :                 CPLGetConfigOption("GRIB_ADJUST_LONGITUDE_RANGE", "YES")))
    2804             :         {
    2805          56 :             if (rPixelSizeX * nRasterXSize > 360 + rPixelSizeX / 4)
    2806           0 :                 CPLDebug("GRIB", "Cannot properly handle GRIB2 files with "
    2807             :                                  "overlaps and 0-360 longitudes");
    2808          56 :             else if (rMinX == 180)
    2809             :             {
    2810             :                 // Case of https://github.com/OSGeo/gdal/issues/10655
    2811           2 :                 CPLDebug("GRIB", "Shifting longitudes from %lf:%lf to %lf:%lf",
    2812             :                          rMinX, rMaxX, -180.0, Lon360to180(rMaxX));
    2813           2 :                 rMinX = -180;
    2814             :             }
    2815          54 :             else if (fabs(360 - rPixelSizeX * nRasterXSize) < rPixelSizeX / 4 &&
    2816          23 :                      rMinX <= 180 && meta->gds.projType == GS3_LATLON)
    2817             :             {
    2818             :                 // Find the first row number east of the antimeridian
    2819           8 :                 const int nSplitAndSwapColumnCandidate =
    2820           8 :                     static_cast<int>(ceil((180 - rMinX) / rPixelSizeX));
    2821           8 :                 if (nSplitAndSwapColumnCandidate < nRasterXSize)
    2822             :                 {
    2823           8 :                     nSplitAndSwapColumn = nSplitAndSwapColumnCandidate;
    2824           8 :                     CPLDebug("GRIB",
    2825             :                              "Rewrapping around the antimeridian at column %d",
    2826             :                              nSplitAndSwapColumn);
    2827           8 :                     rMinX = -180;
    2828           8 :                 }
    2829             :             }
    2830          46 :             else if (Lon360to180(rMinX) > Lon360to180(rMaxX))
    2831             :             {
    2832           2 :                 CPLDebug("GRIB", "GRIB with 0-360 longitudes spanning across "
    2833             :                                  "the antimeridian");
    2834           2 :                 rMinX = Lon360to180(rMinX);
    2835             :             }
    2836             :             else
    2837             :             {
    2838          44 :                 CPLDebug("GRIB", "Shifting longitudes from %lf:%lf to %lf:%lf",
    2839             :                          rMinX, rMaxX, Lon360to180(rMinX), Lon360to180(rMaxX));
    2840          44 :                 rMinX = Lon360to180(rMinX);
    2841             :             }
    2842             :         }
    2843             :     }
    2844             : 
    2845             :     // https://gdal.org/user/raster_data_model.html :
    2846             :     //   we need the top left corner of the top left pixel.
    2847             :     //   At the moment we have the center of the pixel.
    2848         325 :     rMinX -= rPixelSizeX / 2;
    2849         325 :     rMaxY += rPixelSizeY / 2;
    2850             : 
    2851         325 :     m_gt[0] = rMinX;
    2852         325 :     m_gt[3] = rMaxY;
    2853         325 :     m_gt[1] = rPixelSizeX;
    2854         325 :     m_gt[5] = -rPixelSizeY;
    2855             : 
    2856         325 :     if (bError)
    2857           0 :         m_poSRS.reset();
    2858             :     else
    2859         325 :         m_poSRS.reset(oSRS.Clone());
    2860         325 :     m_poLL = std::unique_ptr<OGRSpatialReference>(oLL.Clone());
    2861             : }
    2862             : 
    2863             : /************************************************************************/
    2864             : /*                       GDALDeregister_GRIB()                          */
    2865             : /************************************************************************/
    2866             : 
    2867        1121 : static void GDALDeregister_GRIB(GDALDriver *)
    2868             : {
    2869        1121 :     if (hGRIBMutex != nullptr)
    2870             :     {
    2871           1 :         MetanameCleanup();
    2872           1 :         CPLDestroyMutex(hGRIBMutex);
    2873           1 :         hGRIBMutex = nullptr;
    2874             :     }
    2875        1121 : }
    2876             : 
    2877             : /************************************************************************/
    2878             : /*                          GDALGRIBDriver                              */
    2879             : /************************************************************************/
    2880             : 
    2881             : class GDALGRIBDriver final : public GDALDriver
    2882             : {
    2883             :     std::recursive_mutex m_oMutex{};
    2884             :     bool m_bHasFullInitMetadata = false;
    2885             :     void InitializeMetadata();
    2886             : 
    2887             :   public:
    2888        1741 :     GDALGRIBDriver() = default;
    2889             : 
    2890             :     char **GetMetadata(const char *pszDomain = "") override;
    2891             :     const char *GetMetadataItem(const char *pszName,
    2892             :                                 const char *pszDomain) override;
    2893             : };
    2894             : 
    2895             : /************************************************************************/
    2896             : /*                         InitializeMetadata()                         */
    2897             : /************************************************************************/
    2898             : 
    2899         884 : void GDALGRIBDriver::InitializeMetadata()
    2900             : {
    2901             :     {
    2902             :         // Defer until necessary the setting of the CreationOptionList
    2903             :         // to let a chance to JPEG2000 drivers to have been loaded.
    2904         884 :         if (!m_bHasFullInitMetadata)
    2905             :         {
    2906         210 :             m_bHasFullInitMetadata = true;
    2907             : 
    2908         420 :             std::vector<CPLString> aosJ2KDrivers;
    2909        1050 :             for (size_t i = 0; i < CPL_ARRAYSIZE(apszJ2KDrivers); i++)
    2910             :             {
    2911         840 :                 if (GDALGetDriverByName(apszJ2KDrivers[i]) != nullptr)
    2912             :                 {
    2913         420 :                     aosJ2KDrivers.push_back(apszJ2KDrivers[i]);
    2914             :                 }
    2915             :             }
    2916             : 
    2917             :             CPLString osCreationOptionList(
    2918             :                 "<CreationOptionList>"
    2919             :                 "   <Option name='DATA_ENCODING' type='string-select' "
    2920             :                 "default='AUTO' "
    2921             :                 "description='How data is encoded internally'>"
    2922             :                 "       <Value>AUTO</Value>"
    2923             :                 "       <Value>SIMPLE_PACKING</Value>"
    2924             :                 "       <Value>COMPLEX_PACKING</Value>"
    2925         420 :                 "       <Value>IEEE_FLOATING_POINT</Value>");
    2926         210 :             if (GDALGetDriverByName("PNG") != nullptr)
    2927         210 :                 osCreationOptionList += "       <Value>PNG</Value>";
    2928         210 :             if (!aosJ2KDrivers.empty())
    2929         210 :                 osCreationOptionList += "       <Value>JPEG2000</Value>";
    2930             :             osCreationOptionList +=
    2931             :                 "   </Option>"
    2932             :                 "   <Option name='NBITS' type='int' default='0' "
    2933             :                 "description='Number of bits per value'/>"
    2934             :                 "   <Option name='DECIMAL_SCALE_FACTOR' type='int' default='0' "
    2935             :                 "description='Value such that raw values are multiplied by "
    2936             :                 "10^DECIMAL_SCALE_FACTOR before integer encoding'/>"
    2937             :                 "   <Option name='SPATIAL_DIFFERENCING_ORDER' type='int' "
    2938             :                 "default='0' "
    2939         210 :                 "description='Order of spatial difference' min='0' max='2'/>";
    2940         210 :             if (!aosJ2KDrivers.empty())
    2941             :             {
    2942             :                 osCreationOptionList +=
    2943             :                     "   <Option name='COMPRESSION_RATIO' type='int' "
    2944             :                     "default='1' min='1' max='100' "
    2945             :                     "description='N:1 target compression ratio for JPEG2000'/>"
    2946             :                     "   <Option name='JPEG2000_DRIVER' type='string-select' "
    2947         210 :                     "description='Explicitly select a JPEG2000 driver'>";
    2948         630 :                 for (size_t i = 0; i < aosJ2KDrivers.size(); i++)
    2949             :                 {
    2950             :                     osCreationOptionList +=
    2951         420 :                         "       <Value>" + aosJ2KDrivers[i] + "</Value>";
    2952             :                 }
    2953         210 :                 osCreationOptionList += "   </Option>";
    2954             :             }
    2955             :             osCreationOptionList +=
    2956             :                 "   <Option name='DISCIPLINE' type='int' "
    2957             :                 "description='Discipline of the processed data'/>"
    2958             :                 "   <Option name='IDS' type='string' "
    2959             :                 "description='String equivalent to the GRIB_IDS metadata "
    2960             :                 "item'/>"
    2961             :                 "   <Option name='IDS_CENTER' type='int' "
    2962             :                 "description='Originating/generating center'/>"
    2963             :                 "   <Option name='IDS_SUBCENTER' type='int' "
    2964             :                 "description='Originating/generating subcenter'/>"
    2965             :                 "   <Option name='IDS_MASTER_TABLE' type='int' "
    2966             :                 "description='GRIB master tables version number'/>"
    2967             :                 "   <Option name='IDS_SIGNF_REF_TIME' type='int' "
    2968             :                 "description='Significance of Reference Time'/>"
    2969             :                 "   <Option name='IDS_REF_TIME' type='string' "
    2970             :                 "description='Reference time as YYYY-MM-DDTHH:MM:SSZ'/>"
    2971             :                 "   <Option name='IDS_PROD_STATUS' type='int' "
    2972             :                 "description='Production Status of Processed data'/>"
    2973             :                 "   <Option name='IDS_TYPE' type='int' "
    2974             :                 "description='Type of processed data'/>"
    2975             :                 "   <Option name='PDS_PDTN' type='int' "
    2976             :                 "description='Product Definition Template Number'/>"
    2977             :                 "   <Option name='PDS_TEMPLATE_NUMBERS' type='string' "
    2978             :                 "description='Product definition template raw numbers'/>"
    2979             :                 "   <Option name='PDS_TEMPLATE_ASSEMBLED_VALUES' type='string' "
    2980             :                 "description='Product definition template assembled values'/>"
    2981             :                 "   <Option name='INPUT_UNIT' type='string' "
    2982             :                 "description='Unit of input values. Only for temperatures. C "
    2983             :                 "or K'/>"
    2984             :                 "   <Option name='BAND_*' type='string' "
    2985             :                 "description='Override options at band level'/>"
    2986         210 :                 "</CreationOptionList>";
    2987             : 
    2988         210 :             SetMetadataItem(GDAL_DMD_CREATIONOPTIONLIST, osCreationOptionList);
    2989             :         }
    2990             :     }
    2991         884 : }
    2992             : 
    2993             : /************************************************************************/
    2994             : /*                            GetMetadata()                             */
    2995             : /************************************************************************/
    2996             : 
    2997         520 : char **GDALGRIBDriver::GetMetadata(const char *pszDomain)
    2998             : {
    2999        1040 :     std::lock_guard oLock(m_oMutex);
    3000         520 :     if (pszDomain == nullptr || EQUAL(pszDomain, ""))
    3001             :     {
    3002         520 :         InitializeMetadata();
    3003             :     }
    3004        1040 :     return GDALDriver::GetMetadata(pszDomain);
    3005             : }
    3006             : 
    3007             : /************************************************************************/
    3008             : /*                          GetMetadataItem()                           */
    3009             : /************************************************************************/
    3010             : 
    3011       60842 : const char *GDALGRIBDriver::GetMetadataItem(const char *pszName,
    3012             :                                             const char *pszDomain)
    3013             : {
    3014      121684 :     std::lock_guard oLock(m_oMutex);
    3015       60842 :     if (pszDomain == nullptr || EQUAL(pszDomain, ""))
    3016             :     {
    3017             :         // Defer until necessary the setting of the CreationOptionList
    3018             :         // to let a chance to JPEG2000 drivers to have been loaded.
    3019       60842 :         if (EQUAL(pszName, GDAL_DMD_CREATIONOPTIONLIST))
    3020         364 :             InitializeMetadata();
    3021             :     }
    3022      121684 :     return GDALDriver::GetMetadataItem(pszName, pszDomain);
    3023             : }
    3024             : 
    3025             : /************************************************************************/
    3026             : /*                         GDALRegister_GRIB()                          */
    3027             : /************************************************************************/
    3028             : 
    3029        2024 : void GDALRegister_GRIB()
    3030             : 
    3031             : {
    3032        2024 :     if (GDALGetDriverByName(DRIVER_NAME) != nullptr)
    3033         283 :         return;
    3034             : 
    3035        1741 :     GDALDriver *poDriver = new GDALGRIBDriver();
    3036        1741 :     GRIBDriverSetCommonMetadata(poDriver);
    3037             : 
    3038        1741 :     poDriver->pfnOpen = GRIBDataset::Open;
    3039        1741 :     poDriver->pfnCreateCopy = GRIBDataset::CreateCopy;
    3040        1741 :     poDriver->pfnUnloadDriver = GDALDeregister_GRIB;
    3041             : 
    3042             : #ifdef USE_AEC
    3043        1741 :     poDriver->SetMetadataItem("HAVE_AEC", "YES");
    3044             : #endif
    3045             : 
    3046        1741 :     GetGDALDriverManager()->RegisterDriver(poDriver);
    3047             : }

Generated by: LCOV version 1.14