LCOV - code coverage report
Current view: top level - frmts/grib - gribcreatecopy.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 1137 1339 84.9 %
Date: 2025-07-09 17:50:03 Functions: 42 43 97.7 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GRIB Driver
       4             :  * Purpose:  GDALDataset driver for GRIB translator for write support
       5             :  * Author:   Even Rouault <even dot rouault at spatialys dot com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2017, Even Rouault <even dot rouault at spatialys dot com>
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ******************************************************************************
      12             :  *
      13             :  */
      14             : 
      15             : /* Support for GRIB2 write capabilities has been funded by Meteorological */
      16             : /* Service of Canada */
      17             : 
      18             : #include "cpl_port.h"
      19             : #include "gribdataset.h"
      20             : #include "gdal_priv_templates.hpp"
      21             : #include "memdataset.h"
      22             : 
      23             : #include <algorithm>
      24             : #include <cmath>
      25             : #include <limits>
      26             : 
      27             : #include "degrib/degrib/meta.h"
      28             : CPL_C_START
      29             : #include "degrib/g2clib/grib2.h"
      30             : CPL_C_END
      31             : 
      32             : /************************************************************************/
      33             : /*                         Lon180to360()                                */
      34             : /************************************************************************/
      35             : 
      36          96 : static inline double Lon180to360(double lon)
      37             : {
      38          96 :     if (lon == 180)
      39           0 :         return 180;
      40          96 :     return fmod(fmod(lon, 360) + 360, 360);
      41             : }
      42             : 
      43             : /************************************************************************/
      44             : /*                             WriteByte()                              */
      45             : /************************************************************************/
      46             : 
      47        7707 : static bool WriteByte(VSILFILE *fp, int nVal)
      48             : {
      49        7707 :     GByte byVal = static_cast<GByte>(nVal);
      50        7707 :     return VSIFWriteL(&byVal, 1, sizeof(byVal), fp) == sizeof(byVal);
      51             : }
      52             : 
      53             : /************************************************************************/
      54             : /*                            WriteSByte()                              */
      55             : /************************************************************************/
      56             : 
      57          16 : static bool WriteSByte(VSILFILE *fp, int nVal)
      58             : {
      59          16 :     signed char sVal = static_cast<signed char>(nVal);
      60          16 :     if (sVal == std::numeric_limits<signed char>::min())
      61           1 :         sVal = std::numeric_limits<signed char>::min() + 1;
      62          16 :     GByte nUnsignedVal = (sVal < 0) ? static_cast<GByte>(-sVal) | 0x80U
      63             :                                     : static_cast<GByte>(sVal);
      64          16 :     return VSIFWriteL(&nUnsignedVal, 1, sizeof(nUnsignedVal), fp) ==
      65          16 :            sizeof(nUnsignedVal);
      66             : }
      67             : 
      68             : /************************************************************************/
      69             : /*                            WriteUInt16()                             */
      70             : /************************************************************************/
      71             : 
      72        1262 : static bool WriteUInt16(VSILFILE *fp, int nVal)
      73             : {
      74        1262 :     GUInt16 usVal = static_cast<GUInt16>(nVal);
      75        1262 :     CPL_MSBPTR16(&usVal);
      76        1262 :     return VSIFWriteL(&usVal, 1, sizeof(usVal), fp) == sizeof(usVal);
      77             : }
      78             : 
      79             : /************************************************************************/
      80             : /*                             WriteInt16()                             */
      81             : /************************************************************************/
      82             : 
      83         231 : static bool WriteInt16(VSILFILE *fp, int nVal)
      84             : {
      85         231 :     GInt16 sVal = static_cast<GInt16>(nVal);
      86         231 :     if (sVal == std::numeric_limits<GInt16>::min())
      87           1 :         sVal = std::numeric_limits<GInt16>::min() + 1;
      88         231 :     GUInt16 nUnsignedVal = (sVal < 0) ? static_cast<GUInt16>(-sVal) | 0x8000U
      89             :                                       : static_cast<GUInt16>(sVal);
      90         231 :     CPL_MSBPTR16(&nUnsignedVal);
      91         231 :     return VSIFWriteL(&nUnsignedVal, 1, sizeof(nUnsignedVal), fp) ==
      92         231 :            sizeof(nUnsignedVal);
      93             : }
      94             : 
      95             : /************************************************************************/
      96             : /*                             WriteUInt32()                            */
      97             : /************************************************************************/
      98             : 
      99        3765 : static bool WriteUInt32(VSILFILE *fp, GUInt32 nVal)
     100             : {
     101        3765 :     CPL_MSBPTR32(&nVal);
     102        3765 :     return VSIFWriteL(&nVal, 1, sizeof(nVal), fp) == sizeof(nVal);
     103             : }
     104             : 
     105             : /************************************************************************/
     106             : /*                             WriteInt32()                             */
     107             : /************************************************************************/
     108             : 
     109        1331 : static bool WriteInt32(VSILFILE *fp, GInt32 nVal)
     110             : {
     111        1331 :     if (nVal == std::numeric_limits<GInt32>::min())
     112           0 :         nVal = std::numeric_limits<GInt32>::min() + 1;
     113        1331 :     GUInt32 nUnsignedVal = (nVal < 0)
     114        1331 :                                ? static_cast<GUInt32>(-nVal) | 0x80000000U
     115             :                                : static_cast<GUInt32>(nVal);
     116        1331 :     CPL_MSBPTR32(&nUnsignedVal);
     117        1331 :     return VSIFWriteL(&nUnsignedVal, 1, sizeof(nUnsignedVal), fp) ==
     118        1331 :            sizeof(nUnsignedVal);
     119             : }
     120             : 
     121             : /************************************************************************/
     122             : /*                            WriteFloat32()                            */
     123             : /************************************************************************/
     124             : 
     125         200 : static bool WriteFloat32(VSILFILE *fp, float fVal)
     126             : {
     127         200 :     CPL_MSBPTR32(&fVal);
     128         200 :     return VSIFWriteL(&fVal, 1, sizeof(fVal), fp) == sizeof(fVal);
     129             : }
     130             : 
     131             : /************************************************************************/
     132             : /*                         PatchSectionSize()                           */
     133             : /************************************************************************/
     134             : 
     135         348 : static void PatchSectionSize(VSILFILE *fp, vsi_l_offset nStartSection)
     136             : {
     137         348 :     vsi_l_offset nCurOffset = VSIFTellL(fp);
     138         348 :     VSIFSeekL(fp, nStartSection, SEEK_SET);
     139         348 :     GUInt32 nSect3Size = static_cast<GUInt32>(nCurOffset - nStartSection);
     140         348 :     WriteUInt32(fp, nSect3Size);
     141         348 :     VSIFSeekL(fp, nCurOffset, SEEK_SET);
     142         348 : }
     143             : 
     144             : /************************************************************************/
     145             : /*                         GRIB2Section3Writer                          */
     146             : /************************************************************************/
     147             : 
     148             : class GRIB2Section3Writer
     149             : {
     150             :     VSILFILE *fp;
     151             :     GDALDataset *poSrcDS;
     152             :     OGRSpatialReference oSRS;
     153             :     const char *pszProjection;
     154             :     double dfLLX, dfLLY, dfURX, dfURY;
     155             :     GDALGeoTransform m_gt{};
     156             :     int nSplitAndSwapColumn = 0;
     157             : 
     158             :     bool WriteScaled(double dfVal, double dfUnit);
     159             :     bool TransformToGeo(double &dfX, double &dfY);
     160             :     bool WriteEllipsoidAndRasterSize();
     161             : 
     162             :     bool WriteGeographic();
     163             :     bool WriteRotatedLatLon(double dfLatSouthernPole, double dfLonSouthernPole,
     164             :                             double dfAxisRotation);
     165             :     bool WriteMercator1SP();
     166             :     bool WriteMercator2SP(OGRSpatialReference *poSRS = nullptr);
     167             :     bool WriteTransverseMercator();
     168             :     bool WritePolarStereographic();
     169             :     bool WriteLCC1SP();
     170             :     bool WriteLCC2SPOrAEA(OGRSpatialReference *poSRS = nullptr);
     171             :     bool WriteLAEA();
     172             : 
     173             :   public:
     174             :     GRIB2Section3Writer(VSILFILE *fpIn, GDALDataset *poSrcDSIn);
     175             : 
     176         162 :     inline int SplitAndSwap() const
     177             :     {
     178         162 :         return nSplitAndSwapColumn;
     179             :     }
     180             : 
     181             :     bool Write();
     182             : };
     183             : 
     184             : /************************************************************************/
     185             : /*                        GRIB2Section3Writer()                         */
     186             : /************************************************************************/
     187             : 
     188         162 : GRIB2Section3Writer::GRIB2Section3Writer(VSILFILE *fpIn, GDALDataset *poSrcDSIn)
     189         162 :     : fp(fpIn), poSrcDS(poSrcDSIn)
     190             : {
     191         162 :     oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     192         162 :     oSRS.importFromWkt(poSrcDS->GetProjectionRef());
     193         162 :     pszProjection = oSRS.GetAttrValue("PROJECTION");
     194             : 
     195         162 :     poSrcDS->GetGeoTransform(m_gt);
     196             : 
     197         162 :     dfLLX = m_gt[0] + m_gt[1] / 2;
     198         162 :     dfLLY = m_gt[3] + m_gt[5] / 2 + (poSrcDS->GetRasterYSize() - 1) * m_gt[5];
     199         162 :     dfURX = m_gt[0] + m_gt[1] / 2 + (poSrcDS->GetRasterXSize() - 1) * m_gt[1];
     200         162 :     dfURY = m_gt[3] + m_gt[5] / 2;
     201         162 :     if (dfURY < dfLLY)
     202             :     {
     203           0 :         double dfTemp = dfURY;
     204           0 :         dfURY = dfLLY;
     205           0 :         dfLLY = dfTemp;
     206             :     }
     207         162 : }
     208             : 
     209             : /************************************************************************/
     210             : /*                     WriteEllipsoidAndRasterSize()                    */
     211             : /************************************************************************/
     212             : 
     213         162 : bool GRIB2Section3Writer::WriteEllipsoidAndRasterSize()
     214             : {
     215         162 :     const double dfSemiMajor = oSRS.GetSemiMajor();
     216         162 :     const double dfSemiMinor = oSRS.GetSemiMinor();
     217         162 :     const double dfInvFlattening = oSRS.GetInvFlattening();
     218         231 :     if (std::abs(dfSemiMajor - 6378137.0) < 0.01 &&
     219          69 :         std::abs(dfInvFlattening - 298.257223563) < 1e-9)  // WGS84
     220             :     {
     221          68 :         WriteByte(fp, 5);  // WGS84
     222          68 :         WriteByte(fp, GRIB2MISSING_u1);
     223          68 :         WriteUInt32(fp, GRIB2MISSING_u4);
     224          68 :         WriteByte(fp, GRIB2MISSING_u1);
     225          68 :         WriteUInt32(fp, GRIB2MISSING_u4);
     226          68 :         WriteByte(fp, GRIB2MISSING_u1);
     227          68 :         WriteUInt32(fp, GRIB2MISSING_u4);
     228             :     }
     229          95 :     else if (std::abs(dfSemiMajor - 6378137.0) < 0.01 &&
     230           1 :              std::abs(dfInvFlattening - 298.257222101) < 1e-9)  // GRS80
     231             :     {
     232           1 :         WriteByte(fp, 4);  // GRS80
     233           1 :         WriteByte(fp, GRIB2MISSING_u1);
     234           1 :         WriteUInt32(fp, GRIB2MISSING_u4);
     235           1 :         WriteByte(fp, GRIB2MISSING_u1);
     236           1 :         WriteUInt32(fp, GRIB2MISSING_u4);
     237           1 :         WriteByte(fp, GRIB2MISSING_u1);
     238           1 :         WriteUInt32(fp, GRIB2MISSING_u4);
     239             :     }
     240          93 :     else if (dfInvFlattening == 0)
     241             :     {
     242             :         // Earth assumed spherical with radius specified (in m)
     243             :         // by data producer
     244          14 :         WriteByte(fp, 1);
     245          14 :         WriteByte(fp, 2);  // scale = * 100
     246          14 :         WriteUInt32(fp, static_cast<GUInt32>(dfSemiMajor * 100 + 0.5));
     247          14 :         WriteByte(fp, GRIB2MISSING_u1);
     248          14 :         WriteUInt32(fp, GRIB2MISSING_u4);
     249          14 :         WriteByte(fp, GRIB2MISSING_u1);
     250          14 :         WriteUInt32(fp, GRIB2MISSING_u4);
     251             :     }
     252             :     else
     253             :     {
     254             :         // Earth assumed oblate spheroid with major and minor axes
     255             :         // specified (in m) by data producer
     256          79 :         WriteByte(fp, 7);
     257          79 :         WriteByte(fp, GRIB2MISSING_u1);
     258          79 :         WriteUInt32(fp, GRIB2MISSING_u4);
     259          79 :         WriteByte(fp, 2);  // scale = * 100
     260          79 :         WriteUInt32(fp, static_cast<GUInt32>(dfSemiMajor * 100 + 0.5));
     261          79 :         WriteByte(fp, 2);  // scale = * 100
     262          79 :         WriteUInt32(fp, static_cast<GUInt32>(dfSemiMinor * 100 + 0.5));
     263             :     }
     264         162 :     WriteUInt32(fp, poSrcDS->GetRasterXSize());
     265         162 :     WriteUInt32(fp, poSrcDS->GetRasterYSize());
     266             : 
     267         162 :     return true;
     268             : }
     269             : 
     270             : /************************************************************************/
     271             : /*                            WriteScaled()                             */
     272             : /************************************************************************/
     273             : 
     274        1297 : bool GRIB2Section3Writer::WriteScaled(double dfVal, double dfUnit)
     275             : {
     276        1297 :     return WriteInt32(fp, static_cast<GInt32>(floor(dfVal / dfUnit + 0.5)));
     277             : }
     278             : 
     279             : /************************************************************************/
     280             : /*                          WriteGeographic()                           */
     281             : /************************************************************************/
     282             : 
     283          71 : bool GRIB2Section3Writer::WriteGeographic()
     284             : {
     285          71 :     WriteUInt16(fp, GS3_LATLON);  // Grid template number
     286             : 
     287          71 :     WriteEllipsoidAndRasterSize();
     288             : 
     289          77 :     if (dfLLX < 0 &&
     290           6 :         CPLTestBool(CPLGetConfigOption("GRIB_ADJUST_LONGITUDE_RANGE", "YES")))
     291             :     {
     292           6 :         CPLDebug("GRIB", "Source longitude range is %lf to %lf", dfLLX, dfURX);
     293           6 :         double dfOrigLLX = dfLLX;
     294           6 :         dfLLX = Lon180to360(dfLLX);
     295           6 :         dfURX = Lon180to360(dfURX);
     296             : 
     297           6 :         if (dfLLX > dfURX)
     298             :         {
     299           4 :             if (fabs(360 - poSrcDS->GetRasterXSize() * m_gt[1]) < m_gt[1] / 4)
     300             :             {
     301             :                 // Find the first row number east of the prime meridian
     302           4 :                 nSplitAndSwapColumn =
     303           4 :                     static_cast<int>(ceil((0 - dfOrigLLX) / m_gt[1]));
     304           4 :                 CPLDebug("GRIB",
     305             :                          "Rewrapping around the prime meridian at column %d",
     306             :                          nSplitAndSwapColumn);
     307           4 :                 dfLLX = 0;
     308           4 :                 dfURX = 360 - m_gt[1];
     309             :             }
     310             :             else
     311             :             {
     312           0 :                 CPLDebug("GRIB", "Writing a GRIB with 0-360 longitudes "
     313             :                                  "crossing the prime meridian");
     314             :             }
     315             :         }
     316           6 :         CPLDebug("GRIB", "Target longitudes range is %lf %lf", dfLLX, dfURX);
     317             :     }
     318             : 
     319          71 :     WriteUInt32(fp, 0);  // Basic angle. 0 equivalent of 1
     320             :     // Subdivisions of basic angle used. ~0 equivalent of 10^6
     321          71 :     WriteUInt32(fp, GRIB2MISSING_u4);
     322          71 :     const double dfAngUnit = 1e-6;
     323          71 :     WriteScaled(dfLLY, dfAngUnit);
     324          71 :     WriteScaled(dfLLX, dfAngUnit);
     325          71 :     WriteByte(fp, GRIB2BIT_3 | GRIB2BIT_4);  // Resolution and component flags
     326          71 :     WriteScaled(dfURY, dfAngUnit);
     327          71 :     WriteScaled(dfURX, dfAngUnit);
     328          71 :     WriteScaled(m_gt[1], dfAngUnit);
     329          71 :     WriteScaled(fabs(m_gt[5]), dfAngUnit);
     330          71 :     WriteByte(fp, GRIB2BIT_2);  // Scanning mode: bottom-to-top
     331             : 
     332          71 :     return true;
     333             : }
     334             : 
     335             : /************************************************************************/
     336             : /*                         WriteRotatedLatLon()                         */
     337             : /************************************************************************/
     338             : 
     339           0 : bool GRIB2Section3Writer::WriteRotatedLatLon(double dfLatSouthernPole,
     340             :                                              double dfLonSouthernPole,
     341             :                                              double dfAxisRotation)
     342             : {
     343           0 :     WriteUInt16(fp, GS3_ROTATED_LATLON);  // Grid template number
     344             : 
     345           0 :     WriteEllipsoidAndRasterSize();
     346             : 
     347           0 :     if (dfLLX < 0 &&
     348           0 :         CPLTestBool(CPLGetConfigOption("GRIB_ADJUST_LONGITUDE_RANGE", "YES")))
     349             :     {
     350           0 :         CPLDebug("GRIB", "Source longitude range is %lf to %lf", dfLLX, dfURX);
     351           0 :         double dfOrigLLX = dfLLX;
     352           0 :         dfLLX = Lon180to360(dfLLX);
     353           0 :         dfURX = Lon180to360(dfURX);
     354             : 
     355           0 :         if (dfLLX > dfURX)
     356             :         {
     357           0 :             if (fabs(360 - poSrcDS->GetRasterXSize() * m_gt[1]) < m_gt[1] / 4)
     358             :             {
     359             :                 // Find the first row number east of the prime meridian
     360           0 :                 nSplitAndSwapColumn =
     361           0 :                     static_cast<int>(ceil((0 - dfOrigLLX) / m_gt[1]));
     362           0 :                 CPLDebug("GRIB",
     363             :                          "Rewrapping around the prime meridian at column %d",
     364             :                          nSplitAndSwapColumn);
     365           0 :                 dfLLX = 0;
     366           0 :                 dfURX = 360 - m_gt[1];
     367             :             }
     368             :             else
     369             :             {
     370           0 :                 CPLDebug("GRIB", "Writing a GRIB with 0-360 longitudes "
     371             :                                  "crossing the prime meridian");
     372             :             }
     373             :         }
     374           0 :         CPLDebug("GRIB", "Target longitudes range is %lf %lf", dfLLX, dfURX);
     375             :     }
     376             : 
     377           0 :     WriteUInt32(fp, 0);  // Basic angle. 0 equivalent of 1
     378             :     // Subdivisions of basic angle used. ~0 equivalent of 10^6
     379           0 :     WriteUInt32(fp, GRIB2MISSING_u4);
     380           0 :     const double dfAngUnit = 1e-6;
     381           0 :     WriteScaled(dfLLY, dfAngUnit);
     382           0 :     WriteScaled(dfLLX, dfAngUnit);
     383           0 :     WriteByte(fp, GRIB2BIT_3 | GRIB2BIT_4);  // Resolution and component flags
     384           0 :     WriteScaled(dfURY, dfAngUnit);
     385           0 :     WriteScaled(dfURX, dfAngUnit);
     386           0 :     WriteScaled(m_gt[1], dfAngUnit);
     387           0 :     WriteScaled(fabs(m_gt[5]), dfAngUnit);
     388           0 :     WriteByte(fp, GRIB2BIT_2);  // Scanning mode: bottom-to-top
     389           0 :     WriteScaled(dfLatSouthernPole, dfAngUnit);
     390           0 :     WriteScaled(Lon180to360(dfLonSouthernPole), dfAngUnit);
     391           0 :     WriteScaled(dfAxisRotation, dfAngUnit);
     392             : 
     393           0 :     return true;
     394             : }
     395             : 
     396             : /************************************************************************/
     397             : /*                           TransformToGeo()                           */
     398             : /************************************************************************/
     399             : 
     400          21 : bool GRIB2Section3Writer::TransformToGeo(double &dfX, double &dfY)
     401             : {
     402          42 :     OGRSpatialReference oLL;  // Construct the "geographic" part of oSRS.
     403          21 :     oLL.CopyGeogCSFrom(&oSRS);
     404          21 :     oLL.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     405             :     OGRCoordinateTransformation *poTransformSRSToLL =
     406          21 :         OGRCreateCoordinateTransformation(&(oSRS), &(oLL));
     407          42 :     if (poTransformSRSToLL == nullptr ||
     408          21 :         !poTransformSRSToLL->Transform(1, &dfX, &dfY))
     409             :     {
     410           0 :         delete poTransformSRSToLL;
     411           0 :         return false;
     412             :     }
     413          21 :     delete poTransformSRSToLL;
     414          21 :     if (dfX < 0.0)
     415          20 :         dfX += 360.0;
     416          21 :     return true;
     417             : }
     418             : 
     419             : /************************************************************************/
     420             : /*                           WriteMercator1SP()                         */
     421             : /************************************************************************/
     422             : 
     423           2 : bool GRIB2Section3Writer::WriteMercator1SP()
     424             : {
     425           2 :     if (oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0) != 0.0)
     426             :     {
     427           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     428             :                  "Mercator_1SP with central_meridian != 0 not supported");
     429           0 :         return false;
     430             :     }
     431           2 :     if (oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0) != 0.0)
     432             :     {
     433           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     434             :                  "Mercator_1SP with latitude_of_origin != 0 not supported");
     435           0 :         return false;
     436             :     }
     437             : 
     438             :     OGRSpatialReference *poMerc2SP =
     439           2 :         oSRS.convertToOtherProjection(SRS_PT_MERCATOR_2SP);
     440           2 :     if (poMerc2SP == nullptr)
     441             :     {
     442           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     443             :                  "Cannot get Mercator_2SP formulation");
     444           0 :         return false;
     445             :     }
     446             : 
     447           2 :     bool bRet = WriteMercator2SP(poMerc2SP);
     448           2 :     delete poMerc2SP;
     449           2 :     return bRet;
     450             : }
     451             : 
     452             : /************************************************************************/
     453             : /*                           WriteMercator2SP()                         */
     454             : /************************************************************************/
     455             : 
     456           7 : bool GRIB2Section3Writer::WriteMercator2SP(OGRSpatialReference *poSRS)
     457             : {
     458           7 :     if (poSRS == nullptr)
     459           5 :         poSRS = &oSRS;
     460             : 
     461           7 :     if (poSRS->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0) != 0.0)
     462             :     {
     463           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     464             :                  "Mercator_2SP with central_meridian != 0 not supported");
     465           0 :         return false;
     466             :     }
     467           7 :     if (poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0) != 0.0)
     468             :     {
     469           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     470             :                  "Mercator_2SP with latitude_of_origin != 0 not supported");
     471           0 :         return false;
     472             :     }
     473             : 
     474           7 :     WriteUInt16(fp, GS3_MERCATOR);  // Grid template number
     475             : 
     476           7 :     WriteEllipsoidAndRasterSize();
     477             : 
     478           7 :     if (!TransformToGeo(dfLLX, dfLLY) || !TransformToGeo(dfURX, dfURY))
     479           0 :         return false;
     480             : 
     481           7 :     const double dfAngUnit = 1e-6;
     482           7 :     WriteScaled(dfLLY, dfAngUnit);
     483           7 :     WriteScaled(dfLLX, dfAngUnit);
     484           7 :     WriteByte(fp, GRIB2BIT_3 | GRIB2BIT_4);  // Resolution and component flags
     485           7 :     WriteScaled(poSRS->GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0),
     486             :                 dfAngUnit);
     487           7 :     WriteScaled(dfURY, dfAngUnit);
     488           7 :     WriteScaled(dfURX, dfAngUnit);
     489           7 :     WriteByte(fp, GRIB2BIT_2);  // Scanning mode: bottom-to-top
     490           7 :     WriteInt32(fp, 0);          // angle of the grid
     491           7 :     const double dfLinearUnit = 1e-3;
     492           7 :     WriteScaled(m_gt[1], dfLinearUnit);
     493           7 :     WriteScaled(fabs(m_gt[5]), dfLinearUnit);
     494             : 
     495           7 :     return true;
     496             : }
     497             : 
     498             : /************************************************************************/
     499             : /*                      WriteTransverseMercator()                       */
     500             : /************************************************************************/
     501             : 
     502          78 : bool GRIB2Section3Writer::WriteTransverseMercator()
     503             : {
     504          78 :     WriteUInt16(fp, GS3_TRANSVERSE_MERCATOR);  // Grid template number
     505          78 :     WriteEllipsoidAndRasterSize();
     506             : 
     507          78 :     const double dfAngUnit = 1e-6;
     508          78 :     WriteScaled(oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0),
     509             :                 dfAngUnit);
     510          78 :     WriteScaled(Lon180to360(oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0)),
     511             :                 dfAngUnit);
     512          78 :     WriteByte(fp, GRIB2BIT_3 | GRIB2BIT_4);  // Resolution and component flags
     513             :     float fScale =
     514          78 :         static_cast<float>(oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 0.0));
     515          78 :     WriteFloat32(fp, fScale);
     516          78 :     const double dfLinearUnit = 1e-2;
     517          78 :     WriteScaled(oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0), dfLinearUnit);
     518          78 :     WriteScaled(oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0), dfLinearUnit);
     519          78 :     WriteByte(fp, GRIB2BIT_2);  // Scanning mode: bottom-to-top
     520          78 :     WriteScaled(m_gt[1], dfLinearUnit);
     521          78 :     WriteScaled(fabs(m_gt[5]), dfLinearUnit);
     522          78 :     WriteScaled(dfLLX, dfLinearUnit);
     523          78 :     WriteScaled(dfLLY, dfLinearUnit);
     524          78 :     WriteScaled(dfURX, dfLinearUnit);
     525          78 :     WriteScaled(dfURY, dfLinearUnit);
     526             : 
     527          78 :     return true;
     528             : }
     529             : 
     530             : /************************************************************************/
     531             : /*                       WritePolarStereographic()                       */
     532             : /************************************************************************/
     533             : 
     534           2 : bool GRIB2Section3Writer::WritePolarStereographic()
     535             : {
     536           2 :     WriteUInt16(fp, GS3_POLAR);  // Grid template number
     537           2 :     WriteEllipsoidAndRasterSize();
     538             : 
     539           2 :     if (!TransformToGeo(dfLLX, dfLLY))
     540           0 :         return false;
     541             : 
     542           2 :     const double dfAngUnit = 1e-6;
     543           2 :     WriteScaled(dfLLY, dfAngUnit);
     544           2 :     WriteScaled(dfLLX, dfAngUnit);
     545           2 :     WriteByte(fp, GRIB2BIT_3 | GRIB2BIT_4);  // Resolution and component flags
     546             :     const double dfLatOrigin =
     547           2 :         oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
     548           2 :     WriteScaled(dfLatOrigin, dfAngUnit);
     549           2 :     WriteScaled(Lon180to360(oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0)),
     550             :                 dfAngUnit);
     551           2 :     const double dfLinearUnit = 1e-3;
     552           2 :     WriteScaled(m_gt[1], dfLinearUnit);
     553           2 :     WriteScaled(fabs(m_gt[5]), dfLinearUnit);
     554             :     // Projection center flag: BIT1=0 North Pole, BIT1=1 South Pole
     555           2 :     WriteByte(fp, (dfLatOrigin < 0) ? GRIB2BIT_1 : 0);
     556           2 :     WriteByte(fp, GRIB2BIT_2);  // Scanning mode: bottom-to-top
     557             : 
     558           2 :     return true;
     559             : }
     560             : 
     561             : /************************************************************************/
     562             : /*                            WriteLCC1SP()                             */
     563             : /************************************************************************/
     564             : 
     565           1 : bool GRIB2Section3Writer::WriteLCC1SP()
     566             : {
     567             :     OGRSpatialReference *poLCC2SP =
     568           1 :         oSRS.convertToOtherProjection(SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP);
     569           1 :     if (poLCC2SP == nullptr)
     570             :     {
     571           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     572             :                  "Cannot get Lambert_Conformal_Conic_2SP formulation");
     573           0 :         return false;
     574             :     }
     575             : 
     576           1 :     bool bRet = WriteLCC2SPOrAEA(poLCC2SP);
     577           1 :     delete poLCC2SP;
     578           1 :     return bRet;
     579             : }
     580             : 
     581             : /************************************************************************/
     582             : /*                            WriteLCC2SPOrAEA()                        */
     583             : /************************************************************************/
     584             : 
     585           3 : bool GRIB2Section3Writer::WriteLCC2SPOrAEA(OGRSpatialReference *poSRS)
     586             : {
     587           3 :     if (poSRS == nullptr)
     588           2 :         poSRS = &oSRS;
     589           3 :     if (EQUAL(poSRS->GetAttrValue("PROJECTION"),
     590             :               SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP))
     591           2 :         WriteUInt16(fp, GS3_LAMBERT);  // Grid template number
     592             :     else
     593           1 :         WriteUInt16(fp, GS3_ALBERS_EQUAL_AREA);  // Grid template number
     594             : 
     595           3 :     WriteEllipsoidAndRasterSize();
     596             : 
     597           3 :     if (!TransformToGeo(dfLLX, dfLLY))
     598           0 :         return false;
     599             : 
     600           3 :     const double dfAngUnit = 1e-6;
     601           3 :     WriteScaled(dfLLY, dfAngUnit);
     602           3 :     WriteScaled(dfLLX, dfAngUnit);
     603             :     // Resolution and component flags. "not applicable" ==> 0 ?
     604           3 :     WriteByte(fp, 0);
     605           3 :     WriteScaled(poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0),
     606             :                 dfAngUnit);
     607           3 :     WriteScaled(Lon180to360(oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0)),
     608             :                 dfAngUnit);
     609           3 :     const double dfLinearUnit = 1e-3;
     610           3 :     WriteScaled(m_gt[1], dfLinearUnit);
     611           3 :     WriteScaled(fabs(m_gt[5]), dfLinearUnit);
     612           3 :     WriteByte(fp, 0);           // Projection centre flag
     613           3 :     WriteByte(fp, GRIB2BIT_2);  // Scanning mode: bottom-to-top
     614           3 :     WriteScaled(poSRS->GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0),
     615             :                 dfAngUnit);
     616           3 :     WriteScaled(poSRS->GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2, 0.0),
     617             :                 dfAngUnit);
     618             :     // Latitude of the southern pole of projection
     619           3 :     WriteUInt32(fp, GRIB2MISSING_u4);
     620             :     // Longitude of the southern pole of projection
     621           3 :     WriteUInt32(fp, GRIB2MISSING_u4);
     622           3 :     return true;
     623             : }
     624             : 
     625             : /************************************************************************/
     626             : /*                              WriteLAEA()                             */
     627             : /************************************************************************/
     628             : 
     629           1 : bool GRIB2Section3Writer::WriteLAEA()
     630             : {
     631           1 :     WriteUInt16(fp, GS3_LAMBERT_AZIMUTHAL);  // Grid template number
     632             : 
     633           1 :     WriteEllipsoidAndRasterSize();
     634             : 
     635           1 :     if (!TransformToGeo(dfLLX, dfLLY) || !TransformToGeo(dfURX, dfURY))
     636           0 :         return false;
     637             : 
     638             :     const bool bNormalizeLongitude =
     639           1 :         CPLTestBool(CPLGetConfigOption("GRIB_ADJUST_LONGITUDE_RANGE", "YES"));
     640             : 
     641           1 :     const double dfAngUnit = 1e-6;
     642           1 :     WriteScaled(dfLLY, dfAngUnit);
     643           1 :     if (!bNormalizeLongitude && dfLLX > 360)
     644           0 :         dfLLX -= 360;
     645           1 :     WriteScaled(dfLLX, dfAngUnit);
     646           1 :     WriteScaled(oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_CENTER, 0.0),
     647             :                 dfAngUnit);
     648             :     const double dfLonCenter =
     649           1 :         oSRS.GetNormProjParm(SRS_PP_LONGITUDE_OF_CENTER, 0.0);
     650           1 :     WriteScaled(bNormalizeLongitude ? Lon180to360(dfLonCenter) : dfLonCenter,
     651             :                 dfAngUnit);
     652           1 :     WriteByte(fp, GRIB2BIT_3 | GRIB2BIT_4);  // Resolution and component flags
     653           1 :     const double dfLinearUnit = 1e-3;
     654           1 :     WriteScaled(m_gt[1], dfLinearUnit);
     655           1 :     WriteScaled(fabs(m_gt[5]), dfLinearUnit);
     656           1 :     WriteByte(fp, GRIB2BIT_2);  // Scanning mode: bottom-to-top
     657           1 :     return true;
     658             : }
     659             : 
     660             : /************************************************************************/
     661             : /*                                Write()                               */
     662             : /************************************************************************/
     663             : 
     664         162 : bool GRIB2Section3Writer::Write()
     665             : {
     666             :     // Section 3: Grid Definition Section
     667         162 :     vsi_l_offset nStartSection = VSIFTellL(fp);
     668             : 
     669         162 :     WriteUInt32(fp, GRIB2MISSING_u4);  // section size
     670             : 
     671         162 :     WriteByte(fp, 3);  // section number
     672             : 
     673             :     // Source of grid definition = Specified in Code Table 3.1
     674         162 :     WriteByte(fp, 0);
     675             : 
     676             :     const GUInt32 nDataPoints =
     677         162 :         static_cast<GUInt32>(poSrcDS->GetRasterXSize()) *
     678         162 :         poSrcDS->GetRasterYSize();
     679         162 :     WriteUInt32(fp, nDataPoints);
     680             : 
     681             :     // Number of octets for optional list of numbers defining number of points
     682         162 :     WriteByte(fp, 0);
     683             : 
     684             :     // Interpretation of list of numbers defining number of points =
     685             :     // No appended list
     686         162 :     WriteByte(fp, 0);
     687             : 
     688         162 :     bool bRet = false;
     689         162 :     if (oSRS.IsGeographic())
     690             :     {
     691          71 :         if (oSRS.IsDerivedGeographic())
     692             :         {
     693             :             const OGR_SRSNode *poConversion =
     694           0 :                 oSRS.GetAttrNode("DERIVINGCONVERSION");
     695           0 :             const char *pszMethod = oSRS.GetAttrValue("METHOD");
     696           0 :             if (!pszMethod)
     697           0 :                 pszMethod = "unknown";
     698             : 
     699           0 :             std::map<std::string, double> oValMap;
     700           0 :             if (poConversion)
     701             :             {
     702           0 :                 for (int iChild = 0; iChild < poConversion->GetChildCount();
     703             :                      iChild++)
     704             :                 {
     705           0 :                     const OGR_SRSNode *poNode = poConversion->GetChild(iChild);
     706           0 :                     if (!EQUAL(poNode->GetValue(), "PARAMETER") ||
     707           0 :                         poNode->GetChildCount() <= 2)
     708           0 :                         continue;
     709           0 :                     const char *pszParamStr = poNode->GetChild(0)->GetValue();
     710           0 :                     const char *pszParamVal = poNode->GetChild(1)->GetValue();
     711           0 :                     oValMap[pszParamStr] = CPLAtof(pszParamVal);
     712             :                 }
     713             :             }
     714             : 
     715           0 :             if (poConversion && EQUAL(pszMethod, "PROJ ob_tran o_proj=longlat"))
     716             :             {
     717           0 :                 const double dfLon0 = oValMap["lon_0"];
     718           0 :                 const double dfLonp = oValMap["o_lon_p"];
     719           0 :                 const double dfLatp = oValMap["o_lat_p"];
     720             : 
     721           0 :                 const double dfLatSouthernPole = -dfLatp;
     722           0 :                 const double dfLonSouthernPole = dfLon0;
     723           0 :                 const double dfAxisRotation = -dfLonp;
     724           0 :                 bRet = WriteRotatedLatLon(dfLatSouthernPole, dfLonSouthernPole,
     725           0 :                                           dfAxisRotation);
     726             :             }
     727           0 :             else if (poConversion &&
     728           0 :                      EQUAL(pszMethod, "Pole rotation (netCDF CF convention)"))
     729             :             {
     730             :                 const double dfGridNorthPoleLat =
     731           0 :                     oValMap["Grid north pole latitude (netCDF CF convention)"];
     732             :                 const double dfGridNorthPoleLong =
     733           0 :                     oValMap["Grid north pole longitude (netCDF CF convention)"];
     734             :                 const double dfNorthPoleGridLong =
     735           0 :                     oValMap["North pole grid longitude (netCDF CF convention)"];
     736             : 
     737           0 :                 const double dfLon0 = 180.0 + dfGridNorthPoleLong;
     738           0 :                 const double dfLonp = dfNorthPoleGridLong;
     739           0 :                 const double dfLatp = dfGridNorthPoleLat;
     740             : 
     741           0 :                 const double dfLatSouthernPole = -dfLatp;
     742           0 :                 const double dfLonSouthernPole = dfLon0;
     743           0 :                 const double dfAxisRotation = -dfLonp;
     744           0 :                 bRet = WriteRotatedLatLon(dfLatSouthernPole, dfLonSouthernPole,
     745           0 :                                           dfAxisRotation);
     746             :             }
     747           0 :             else if (poConversion &&
     748           0 :                      EQUAL(pszMethod, "Pole rotation (GRIB convention)"))
     749             :             {
     750             :                 const double dfLatSouthernPole =
     751           0 :                     oValMap["Latitude of the southern pole (GRIB convention)"];
     752             :                 const double dfLonSouthernPole =
     753           0 :                     oValMap["Longitude of the southern pole (GRIB convention)"];
     754             :                 const double dfAxisRotation =
     755           0 :                     oValMap["Axis rotation (GRIB convention)"];
     756             : 
     757           0 :                 bRet = WriteRotatedLatLon(dfLatSouthernPole, dfLonSouthernPole,
     758           0 :                                           dfAxisRotation);
     759             :             }
     760             :             else
     761             :             {
     762           0 :                 CPLError(CE_Failure, CPLE_NotSupported,
     763             :                          "Unsupported method for DerivedGeographicCRS: %s",
     764             :                          pszMethod);
     765           0 :                 return false;
     766             :             }
     767             :         }
     768             :         else
     769             :         {
     770          71 :             bRet = WriteGeographic();
     771             :         }
     772             :     }
     773          91 :     else if (pszProjection && EQUAL(pszProjection, SRS_PT_MERCATOR_1SP))
     774             :     {
     775           2 :         bRet = WriteMercator1SP();
     776             :     }
     777          89 :     else if (pszProjection && EQUAL(pszProjection, SRS_PT_MERCATOR_2SP))
     778             :     {
     779           5 :         bRet = WriteMercator2SP();
     780             :     }
     781          84 :     else if (pszProjection && EQUAL(pszProjection, SRS_PT_TRANSVERSE_MERCATOR))
     782             :     {
     783          78 :         bRet = WriteTransverseMercator();
     784             :     }
     785           6 :     else if (pszProjection && EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC))
     786             :     {
     787           2 :         bRet = WritePolarStereographic();
     788             :     }
     789           4 :     else if (pszProjection != nullptr &&
     790           4 :              EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP))
     791             :     {
     792           1 :         bRet = WriteLCC1SP();
     793             :     }
     794           3 :     else if (pszProjection != nullptr &&
     795           3 :              (EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP) ||
     796           2 :               EQUAL(pszProjection, SRS_PT_ALBERS_CONIC_EQUAL_AREA)))
     797             :     {
     798           2 :         bRet = WriteLCC2SPOrAEA();
     799             :     }
     800           1 :     else if (pszProjection &&
     801           1 :              EQUAL(pszProjection, SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA))
     802             :     {
     803           1 :         bRet = WriteLAEA();
     804             :     }
     805             : 
     806         162 :     PatchSectionSize(fp, nStartSection);
     807             : 
     808         162 :     return bRet;
     809             : }
     810             : 
     811             : /************************************************************************/
     812             : /*                         GetBandOption()                              */
     813             : /************************************************************************/
     814             : 
     815        3999 : static const char *GetBandOption(char **papszOptions, GDALDataset *poSrcDS,
     816             :                                  int nBand, const char *pszKey,
     817             :                                  const char *pszDefault)
     818             : {
     819        3999 :     const char *pszVal = CSLFetchNameValue(
     820             :         papszOptions, CPLSPrintf("BAND_%d_%s", nBand, pszKey));
     821        3999 :     if (pszVal == nullptr)
     822             :     {
     823        3998 :         pszVal = CSLFetchNameValue(papszOptions, pszKey);
     824             :     }
     825        3999 :     if (pszVal == nullptr && poSrcDS != nullptr)
     826             :     {
     827        3134 :         pszVal = poSrcDS->GetRasterBand(nBand)->GetMetadataItem(
     828        3134 :             (CPLString("GRIB_") + pszKey).c_str());
     829             :     }
     830        3999 :     if (pszVal == nullptr)
     831             :     {
     832        3613 :         pszVal = pszDefault;
     833             :     }
     834        3999 :     return pszVal;
     835             : }
     836             : 
     837             : /************************************************************************/
     838             : /*                        GRIB2Section567Writer                         */
     839             : /************************************************************************/
     840             : 
     841             : class GRIB2Section567Writer
     842             : {
     843             :     VSILFILE *m_fp;
     844             :     GDALDataset *m_poSrcDS;
     845             :     int m_nBand;
     846             :     int m_nXSize;
     847             :     int m_nYSize;
     848             :     GUInt32 m_nDataPoints;
     849             :     GDALDataType m_eDT;
     850             :     GDALGeoTransform m_gt{};
     851             :     int m_nDecimalScaleFactor;
     852             :     double m_dfDecimalScale;
     853             :     float m_fMin;
     854             :     float m_fMax;
     855             :     double m_dfMinScaled;
     856             :     int m_nBits;
     857             :     bool m_bUseZeroBits;
     858             :     float m_fValOffset;
     859             :     int m_bHasNoData;
     860             :     double m_dfNoData;
     861             :     int m_nSplitAndSwap;
     862             : 
     863             :     float *GetFloatData();
     864             :     bool WriteSimplePacking();
     865             :     bool WriteComplexPacking(int nSpatialDifferencingOrder);
     866             :     bool WriteIEEE(GDALProgressFunc pfnProgress, void *pProgressData);
     867             :     bool WritePNG();
     868             :     bool WriteJPEG2000(char **papszOptions);
     869             : 
     870             :   public:
     871             :     GRIB2Section567Writer(VSILFILE *fp, GDALDataset *poSrcDS, int nBand,
     872             :                           int nSplitAndSwap);
     873             : 
     874             :     bool Write(float fValOffset, char **papszOptions,
     875             :                GDALProgressFunc pfnProgress, void *pProgressData);
     876             :     void WriteComplexPackingNoData();
     877             : };
     878             : 
     879             : /************************************************************************/
     880             : /*                      GRIB2Section567Writer()                         */
     881             : /************************************************************************/
     882             : 
     883         156 : GRIB2Section567Writer::GRIB2Section567Writer(VSILFILE *fp, GDALDataset *poSrcDS,
     884         156 :                                              int nBand, int nSplitAndSwap)
     885             :     : m_fp(fp), m_poSrcDS(poSrcDS), m_nBand(nBand),
     886         468 :       m_nXSize(poSrcDS->GetRasterXSize()), m_nYSize(poSrcDS->GetRasterYSize()),
     887         156 :       m_nDataPoints(static_cast<GUInt32>(m_nXSize) * m_nYSize),
     888         156 :       m_eDT(m_poSrcDS->GetRasterBand(m_nBand)->GetRasterDataType()),
     889             :       m_nDecimalScaleFactor(0), m_dfDecimalScale(1.0), m_fMin(0.0f),
     890             :       m_fMax(0.0f), m_dfMinScaled(0.0), m_nBits(0), m_bUseZeroBits(false),
     891             :       m_fValOffset(0.0), m_bHasNoData(false), m_dfNoData(0.0),
     892         312 :       m_nSplitAndSwap(nSplitAndSwap)
     893             : {
     894         156 :     m_poSrcDS->GetGeoTransform(m_gt);
     895         156 :     m_dfNoData = m_poSrcDS->GetRasterBand(nBand)->GetNoDataValue(&m_bHasNoData);
     896         156 : }
     897             : 
     898             : /************************************************************************/
     899             : /*                          GetFloatData()                              */
     900             : /************************************************************************/
     901             : 
     902         118 : float *GRIB2Section567Writer::GetFloatData()
     903             : {
     904             :     float *pafData =
     905         118 :         static_cast<float *>(VSI_MALLOC2_VERBOSE(m_nDataPoints, sizeof(float)));
     906         118 :     if (pafData == nullptr)
     907             :     {
     908           0 :         return nullptr;
     909             :     }
     910         118 :     CPLErr eErr = m_poSrcDS->GetRasterBand(m_nBand)->RasterIO(
     911         118 :         GF_Read, m_nSplitAndSwap, 0, m_nXSize - m_nSplitAndSwap, m_nYSize,
     912         118 :         pafData + (m_gt[5] < 0 ? (m_nYSize - 1) * m_nXSize : 0),
     913         118 :         m_nXSize - m_nSplitAndSwap, m_nYSize, GDT_Float32, sizeof(float),
     914         118 :         m_gt[5] < 0 ? -static_cast<GSpacing>(m_nXSize * sizeof(float))
     915           0 :                     : static_cast<GSpacing>(m_nXSize * sizeof(float)),
     916             :         nullptr);
     917         118 :     if (eErr != CE_None)
     918             :     {
     919           0 :         VSIFree(pafData);
     920           0 :         return nullptr;
     921             :     }
     922         118 :     if (m_nSplitAndSwap > 0)
     923             :     {
     924           1 :         eErr = m_poSrcDS->GetRasterBand(m_nBand)->RasterIO(
     925             :             GF_Read, 0, 0, m_nSplitAndSwap, m_nYSize,
     926           1 :             pafData + (m_gt[5] < 0 ? (m_nYSize - 1) * m_nXSize : 0) +
     927           1 :                 (m_nXSize - m_nSplitAndSwap),
     928             :             m_nSplitAndSwap, m_nYSize, GDT_Float32, sizeof(float),
     929           1 :             m_gt[5] < 0 ? -static_cast<GSpacing>(m_nXSize * sizeof(float))
     930           0 :                         : static_cast<GSpacing>(m_nXSize * sizeof(float)),
     931             :             nullptr);
     932           1 :         if (eErr != CE_None)
     933             :         {
     934           0 :             VSIFree(pafData);
     935           0 :             return nullptr;
     936             :         }
     937             :     }
     938             : 
     939         118 :     m_fMin = std::numeric_limits<float>::max();
     940         118 :     m_fMax = -std::numeric_limits<float>::max();
     941         118 :     bool bHasNoDataValuePoint = false;
     942         118 :     bool bHasDataValuePoint = false;
     943     1681620 :     for (GUInt32 i = 0; i < m_nDataPoints; i++)
     944             :     {
     945     1681500 :         if (m_bHasNoData && pafData[i] == static_cast<float>(m_dfNoData))
     946             :         {
     947      680030 :             if (!bHasNoDataValuePoint)
     948           8 :                 bHasNoDataValuePoint = true;
     949      680030 :             continue;
     950             :         }
     951     1001470 :         if (!std::isfinite(pafData[i]))
     952             :         {
     953           0 :             CPLError(CE_Failure, CPLE_NotSupported,
     954             :                      "Non-finite values not supported for "
     955             :                      "this data encoding");
     956           0 :             VSIFree(pafData);
     957           0 :             return nullptr;
     958             :         }
     959     1001470 :         if (!bHasDataValuePoint)
     960         116 :             bHasDataValuePoint = true;
     961     1001470 :         pafData[i] += m_fValOffset;
     962     1001470 :         if (pafData[i] < m_fMin)
     963         623 :             m_fMin = pafData[i];
     964     1001470 :         if (pafData[i] > m_fMax)
     965         290 :             m_fMax = pafData[i];
     966             :     }
     967         118 :     if (m_fMin > m_fMax)
     968             :     {
     969           2 :         m_fMin = m_fMax = static_cast<float>(m_dfNoData);
     970             :     }
     971             : 
     972             :     // We check that the actual range of values got from the above RasterIO
     973             :     // request does not go over the expected range of the datatype, as we
     974             :     // later assume that for computing nMaxBitsPerElt.
     975             :     // This shouldn't happen for well-behaved drivers, but this can still
     976             :     // happen in practice, if some drivers don't completely fill buffers etc.
     977         177 :     if (m_fMax > m_fMin && GDALDataTypeIsInteger(m_eDT) &&
     978          59 :         ceil(log(m_fMax - m_fMin) / log(2.0)) > GDALGetDataTypeSizeBits(m_eDT))
     979             :     {
     980           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     981             :                  "Garbage values found when requesting input dataset");
     982           0 :         VSIFree(pafData);
     983           0 :         return nullptr;
     984             :     }
     985             : 
     986         118 :     m_dfMinScaled =
     987         118 :         m_dfDecimalScale == 1.0 ? m_fMin : floor(m_fMin * m_dfDecimalScale);
     988         236 :     if (!(m_dfMinScaled >= -std::numeric_limits<float>::max() &&
     989         118 :           m_dfMinScaled < std::numeric_limits<float>::max()))
     990             :     {
     991           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     992             :                  "Scaled min value not representable on IEEE754 "
     993             :                  "single precision float");
     994           0 :         VSIFree(pafData);
     995           0 :         return nullptr;
     996             :     }
     997             : 
     998         118 :     const double dfScaledMaxDiff = (m_fMax - m_fMin) * m_dfDecimalScale;
     999         118 :     if (GDALDataTypeIsFloating(m_eDT) && m_nBits == 0 && dfScaledMaxDiff > 0 &&
    1000             :         dfScaledMaxDiff <= 256)
    1001             :     {
    1002           8 :         m_nBits = 8;
    1003             :     }
    1004             : 
    1005         118 :     m_bUseZeroBits =
    1006         194 :         (m_fMin == m_fMax && !(bHasDataValuePoint && bHasNoDataValuePoint)) ||
    1007          76 :         (!GDALDataTypeIsFloating(m_eDT) && dfScaledMaxDiff < 1.0);
    1008             : 
    1009         118 :     return pafData;
    1010             : }
    1011             : 
    1012             : /************************************************************************/
    1013             : /*                        WriteSimplePacking()                          */
    1014             : /************************************************************************/
    1015             : 
    1016             : // See http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_temp5-0.shtml
    1017          66 : bool GRIB2Section567Writer::WriteSimplePacking()
    1018             : {
    1019          66 :     float *pafData = GetFloatData();
    1020          66 :     if (pafData == nullptr)
    1021           0 :         return false;
    1022             : 
    1023             :     const int nBitCorrectionForDec =
    1024          66 :         static_cast<int>(ceil(m_nDecimalScaleFactor * log(10.0) / log(2.0)));
    1025             :     const int nMaxBitsPerElt = std::max(
    1026          66 :         1, std::min(31, (m_nBits > 0) ? m_nBits
    1027          56 :                                       : GDALGetDataTypeSizeBits(m_eDT) +
    1028          66 :                                             nBitCorrectionForDec));
    1029          66 :     if (nMaxBitsPerElt > 0 &&
    1030          66 :         m_nDataPoints > static_cast<GUInt32>(INT_MAX) / nMaxBitsPerElt)
    1031             :     {
    1032           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    1033             :                  "Int overflow while computing maximum number of bits");
    1034           0 :         VSIFree(pafData);
    1035           0 :         return false;
    1036             :     }
    1037             : 
    1038          66 :     const int nMaxSize = (m_nDataPoints * nMaxBitsPerElt + 7) / 8;
    1039          66 :     void *pabyData = VSI_MALLOC_VERBOSE(nMaxSize);
    1040          66 :     if (pabyData == nullptr)
    1041             :     {
    1042           0 :         VSIFree(pafData);
    1043           0 :         VSIFree(pabyData);
    1044           0 :         return false;
    1045             :     }
    1046             : 
    1047             :     // Indices expected by simpack()
    1048             :     enum
    1049             :     {
    1050             :         TMPL5_R_IDX = 0,      // Reference value (R)
    1051             :         TMPL5_E_IDX = 1,      // Binary scale factor (E)
    1052             :         TMPL5_D_IDX = 2,      // Decimal scale factor (D)
    1053             :         TMPL5_NBITS_IDX = 3,  // Number of bits used for each packed value
    1054             :         TMPL5_TYPE_IDX = 4    // type of original data
    1055             :     };
    1056             : 
    1057          66 :     g2int idrstmpl[TMPL5_TYPE_IDX + 1] = {0};
    1058          66 :     idrstmpl[TMPL5_R_IDX] = 0;  // reference value, to be filled by simpack
    1059          66 :     idrstmpl[TMPL5_E_IDX] = 0;  // binary scale factor, to be filled by simpack
    1060          66 :     idrstmpl[TMPL5_D_IDX] = m_nDecimalScaleFactor;
    1061             :     // to be filled by simpack if set to 0
    1062          66 :     idrstmpl[TMPL5_NBITS_IDX] = m_nBits;
    1063             :     // to be filled by simpack (and we will ignore it)
    1064          66 :     idrstmpl[TMPL5_TYPE_IDX] = 0;
    1065          66 :     g2int nLengthPacked = 0;
    1066          66 :     simpack(pafData, m_nDataPoints, idrstmpl,
    1067             :             static_cast<unsigned char *>(pabyData), &nLengthPacked);
    1068          66 :     CPLAssert(nLengthPacked <= nMaxSize);
    1069          66 :     if (nLengthPacked < 0)
    1070             :     {
    1071           0 :         CPLError(CE_Failure, CPLE_AppDefined, "Error while packing");
    1072           0 :         VSIFree(pafData);
    1073           0 :         VSIFree(pabyData);
    1074           0 :         return false;
    1075             :     }
    1076             : 
    1077             :     // Section 5: Data Representation Section
    1078          66 :     WriteUInt32(m_fp, 21);  // section size
    1079          66 :     WriteByte(m_fp, 5);     // section number
    1080          66 :     WriteUInt32(m_fp, m_nDataPoints);
    1081          66 :     WriteUInt16(m_fp, GS5_SIMPLE);
    1082             :     float fRefValue;
    1083          66 :     memcpy(&fRefValue, &idrstmpl[TMPL5_R_IDX], 4);
    1084          66 :     WriteFloat32(m_fp, fRefValue);
    1085          66 :     WriteInt16(m_fp, idrstmpl[TMPL5_E_IDX]);
    1086          66 :     WriteInt16(m_fp, idrstmpl[TMPL5_D_IDX]);
    1087          66 :     WriteByte(m_fp, idrstmpl[TMPL5_NBITS_IDX]);
    1088             :     // Type of original data: 0=Floating, 1=Integer
    1089          66 :     WriteByte(m_fp, GDALDataTypeIsFloating(m_eDT) ? 0 : 1);
    1090             : 
    1091             :     // Section 6: Bitmap section
    1092             : #ifdef DEBUG
    1093          66 :     if (CPLTestBool(CPLGetConfigOption("GRIB_WRITE_BITMAP_TEST", "NO")))
    1094             :     {
    1095             :         // Just for the purpose of generating a test product !
    1096             :         static int counter = 0;
    1097           0 :         counter++;
    1098           0 :         if (counter == 1)
    1099             :         {
    1100           0 :             WriteUInt32(m_fp, 6 + ((m_nDataPoints + 7) / 8));  // section size
    1101           0 :             WriteByte(m_fp, 6);                                // section number
    1102           0 :             WriteByte(m_fp, 0);                                // bitmap
    1103           0 :             for (GUInt32 i = 0; i < (m_nDataPoints + 7) / 8; i++)
    1104           0 :                 WriteByte(m_fp, 255);
    1105             :         }
    1106             :         else
    1107             :         {
    1108           0 :             WriteUInt32(m_fp, 6);  // section size
    1109           0 :             WriteByte(m_fp, 6);    // section number
    1110           0 :             WriteByte(m_fp, 254);  // reuse previous bitmap
    1111             :         }
    1112             :     }
    1113             :     else
    1114             : #endif
    1115             :     {
    1116          66 :         WriteUInt32(m_fp, 6);              // section size
    1117          66 :         WriteByte(m_fp, 6);                // section number
    1118          66 :         WriteByte(m_fp, GRIB2MISSING_u1);  // no bitmap
    1119             :     }
    1120             : 
    1121             :     // Section 7: Data Section
    1122          66 :     WriteUInt32(m_fp, 5 + nLengthPacked);  // section size
    1123          66 :     WriteByte(m_fp, 7);                    // section number
    1124          66 :     if (static_cast<int>(VSIFWriteL(pabyData, 1, nLengthPacked, m_fp)) !=
    1125             :         nLengthPacked)
    1126             :     {
    1127           0 :         VSIFree(pafData);
    1128           0 :         VSIFree(pabyData);
    1129           0 :         return false;
    1130             :     }
    1131             : 
    1132          66 :     VSIFree(pafData);
    1133          66 :     VSIFree(pabyData);
    1134             : 
    1135          66 :     return true;
    1136             : }
    1137             : 
    1138             : /************************************************************************/
    1139             : /*                      WriteComplexPackingNoData()                     */
    1140             : /************************************************************************/
    1141             : 
    1142          20 : void GRIB2Section567Writer::WriteComplexPackingNoData()
    1143             : {
    1144          20 :     if (!m_bHasNoData)
    1145             :     {
    1146           9 :         WriteUInt32(m_fp, GRIB2MISSING_u4);
    1147             :     }
    1148          11 :     else if (GDALDataTypeIsFloating(m_eDT))
    1149             :     {
    1150           7 :         WriteFloat32(m_fp, static_cast<float>(m_dfNoData));
    1151             :     }
    1152             :     else
    1153             :     {
    1154           4 :         if (GDALIsValueInRange<int>(m_dfNoData))
    1155             :         {
    1156           4 :             WriteInt32(m_fp, static_cast<int>(m_dfNoData));
    1157             :         }
    1158             :         else
    1159             :         {
    1160           0 :             WriteUInt32(m_fp, GRIB2MISSING_u4);
    1161             :         }
    1162             :     }
    1163          20 : }
    1164             : 
    1165             : /************************************************************************/
    1166             : /*                       WriteComplexPacking()                          */
    1167             : /************************************************************************/
    1168             : 
    1169             : // See http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_temp5-2.shtml
    1170             : // and http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_temp5-3.shtml
    1171             : 
    1172          22 : bool GRIB2Section567Writer::WriteComplexPacking(int nSpatialDifferencingOrder)
    1173             : {
    1174          22 :     if (nSpatialDifferencingOrder < 0 || nSpatialDifferencingOrder > 2)
    1175             :     {
    1176           1 :         CPLError(CE_Failure, CPLE_AppDefined,
    1177             :                  "Unsupported value for SPATIAL_DIFFERENCING_ORDER");
    1178           1 :         return false;
    1179             :     }
    1180             : 
    1181          21 :     float *pafData = GetFloatData();
    1182          21 :     if (pafData == nullptr)
    1183           0 :         return false;
    1184             : 
    1185          21 :     const float fNoData = static_cast<float>(m_dfNoData);
    1186          21 :     if (m_bUseZeroBits)
    1187             :     {
    1188             :         // Case where all values are at nodata or a single value
    1189           8 :         VSIFree(pafData);
    1190             : 
    1191             :         // Section 5: Data Representation Section
    1192           8 :         WriteUInt32(m_fp, 47);  // section size
    1193           8 :         WriteByte(m_fp, 5);     // section number
    1194           8 :         WriteUInt32(m_fp, m_nDataPoints);
    1195           8 :         WriteUInt16(m_fp, GS5_CMPLX);
    1196           8 :         WriteFloat32(m_fp, m_fMin);  // ref value = nodata or single data
    1197           8 :         WriteInt16(m_fp, 0);         // binary scale factor
    1198           8 :         WriteInt16(m_fp, 0);         // decimal scale factor
    1199           8 :         WriteByte(m_fp, 0);          // number of bits
    1200             :         // Type of original data: 0=Floating, 1=Integer
    1201           8 :         WriteByte(m_fp, GDALDataTypeIsFloating(m_eDT) ? 0 : 1);
    1202           8 :         WriteByte(m_fp, 0);
    1203           8 :         WriteByte(m_fp, m_bHasNoData ? 1 : 0);  // 1 missing value
    1204           8 :         WriteComplexPackingNoData();
    1205           8 :         WriteUInt32(m_fp, GRIB2MISSING_u4);
    1206           8 :         WriteUInt32(m_fp, 0);
    1207           8 :         WriteByte(m_fp, 0);
    1208           8 :         WriteByte(m_fp, 0);
    1209           8 :         WriteUInt32(m_fp, 0);
    1210           8 :         WriteByte(m_fp, 0);
    1211           8 :         WriteUInt32(m_fp, 0);
    1212           8 :         WriteByte(m_fp, 0);
    1213             : 
    1214             :         // Section 6: Bitmap section
    1215           8 :         WriteUInt32(m_fp, 6);              // section size
    1216           8 :         WriteByte(m_fp, 6);                // section number
    1217           8 :         WriteByte(m_fp, GRIB2MISSING_u1);  // no bitmap
    1218             : 
    1219             :         // Section 7: Data Section
    1220           8 :         WriteUInt32(m_fp, 5);  // section size
    1221           8 :         WriteByte(m_fp, 7);    // section number
    1222             : 
    1223           8 :         return true;
    1224             :     }
    1225             : 
    1226             :     const int nBitCorrectionForDec =
    1227          13 :         static_cast<int>(ceil(m_nDecimalScaleFactor * log(10.0) / log(2.0)));
    1228             :     const int nMaxBitsPerElt = std::max(
    1229          13 :         1, std::min(31, (m_nBits > 0) ? m_nBits
    1230           6 :                                       : GDALGetDataTypeSizeBits(m_eDT) +
    1231          13 :                                             nBitCorrectionForDec));
    1232          13 :     if (nMaxBitsPerElt > 0 &&
    1233          13 :         m_nDataPoints > static_cast<GUInt32>(INT_MAX) / nMaxBitsPerElt)
    1234             :     {
    1235           0 :         CPLError(CE_Failure, CPLE_AppDefined,
    1236             :                  "Int overflow while computing maximum number of bits");
    1237           0 :         VSIFree(pafData);
    1238           0 :         return false;
    1239             :     }
    1240             : 
    1241             :     // No idea what is the exact maximum bound... Take the value of simple
    1242             :     // packing and multiply by 2, plus some constant...
    1243          13 :     const int nMaxSize = 10000 + 2 * ((m_nDataPoints * nMaxBitsPerElt + 7) / 8);
    1244          13 :     void *pabyData = VSI_MALLOC_VERBOSE(nMaxSize);
    1245          13 :     if (pabyData == nullptr)
    1246             :     {
    1247           0 :         VSIFree(pafData);
    1248           0 :         VSIFree(pabyData);
    1249           0 :         return false;
    1250             :     }
    1251             : 
    1252          13 :     const double dfScaledMaxDiff =
    1253          13 :         (m_fMax == m_fMin) ? 1 : (m_fMax - m_fMin) * m_dfDecimalScale;
    1254          13 :     if (m_nBits == 0)
    1255             :     {
    1256           6 :         double dfTemp = log(ceil(dfScaledMaxDiff)) / log(2.0);
    1257           6 :         m_nBits = std::max(1, std::min(31, static_cast<int>(ceil(dfTemp))));
    1258             :     }
    1259          13 :     const int nMaxNum = (m_nBits == 31) ? INT_MAX : ((1 << m_nBits) - 1);
    1260          13 :     double dfTemp = log(nMaxNum / dfScaledMaxDiff) / log(2.0);
    1261          13 :     int nBinaryScaleFactor = static_cast<GInt16>(ceil(-dfTemp));
    1262             : 
    1263             :     // Indices expected by cmplxpack()
    1264             :     enum
    1265             :     {
    1266             :         TMPL5_R_IDX = 0,                        // reference value
    1267             :         TMPL5_E_IDX = 1,                        // binary scale factor
    1268             :         TMPL5_D_IDX = 2,                        // decimal scale factor
    1269             :         TMPL5_NBITS_IDX = 3,                    // number of bits
    1270             :         TMPL5_TYPE_IDX = 4,                     // type of original data
    1271             :         TMPL5_GROUP_SPLITTING_IDX = 5,          // Group splitting method used
    1272             :         TMPL5_MISSING_VALUE_MGNT_IDX = 6,       // Missing value management used
    1273             :         TMPL5_PRIMARY_MISSING_VALUE_IDX = 7,    // Primary missing value
    1274             :         TMPL5_SECONDARY_MISSING_VALUE_IDX = 8,  // Secondary missing value
    1275             :         TMPL5_NG_IDX = 9,                 // number of groups of data values
    1276             :         TMPL5_REF_GROUP_WIDTHS_IDX = 10,  // Reference for group widths
    1277             :         // Number of bits used for the group widths
    1278             :         TMPL5_NBITS_GROUP_WIDTHS_IDX = 11,
    1279             :         TMPL5_REF_GROUP_LENGTHS_IDX = 12,  // Reference for group lengths
    1280             :         // Length increment for the group lengths
    1281             :         TMPL5_LENGTH_INCR_GROUP_LENGTHS_IDX = 13,
    1282             :         TMPL5_TRUE_LENGTH_LAST_GROUP_IDX = 14,  // True length of last group
    1283             :         // Number of bits used for the scaled group lengths
    1284             :         TMPL5_NBITS_SCALED_GROUP_LENGTHS_IDX = 15,
    1285             :         // Order of spatial differencing
    1286             :         TMPL5_ORDER_SPATIAL_DIFFERENCE_IDX = 16,
    1287             :         // Number of octets required in the data section to specify extra
    1288             :         // descriptors needed for spatial differencing
    1289             :         TMPL5_NB_OCTETS_EXTRA_DESCR_IDX = 17
    1290             :     };
    1291             : 
    1292          13 :     g2int idrstmpl[TMPL5_NB_OCTETS_EXTRA_DESCR_IDX + 1] = {0};
    1293          13 :     idrstmpl[TMPL5_E_IDX] = nBinaryScaleFactor;
    1294          13 :     idrstmpl[TMPL5_D_IDX] = m_nDecimalScaleFactor;
    1295          13 :     idrstmpl[TMPL5_MISSING_VALUE_MGNT_IDX] = m_bHasNoData ? 1 : 0;
    1296          13 :     idrstmpl[TMPL5_ORDER_SPATIAL_DIFFERENCE_IDX] = nSpatialDifferencingOrder;
    1297          13 :     if (m_bHasNoData)
    1298             :     {
    1299           4 :         memcpy(&idrstmpl[TMPL5_PRIMARY_MISSING_VALUE_IDX], &fNoData, 4);
    1300             :     }
    1301          13 :     g2int nLengthPacked = 0;
    1302          13 :     const int nTemplateNumber =
    1303          13 :         (nSpatialDifferencingOrder > 0) ? GS5_CMPLXSEC : GS5_CMPLX;
    1304          13 :     cmplxpack(pafData, m_nDataPoints, nTemplateNumber, idrstmpl,
    1305             :               static_cast<unsigned char *>(pabyData), &nLengthPacked);
    1306          13 :     CPLAssert(nLengthPacked <= nMaxSize);
    1307          13 :     if (nLengthPacked < 0)
    1308             :     {
    1309           1 :         CPLError(CE_Failure, CPLE_AppDefined, "Error while packing");
    1310           1 :         VSIFree(pafData);
    1311           1 :         VSIFree(pabyData);
    1312           1 :         return false;
    1313             :     }
    1314             : 
    1315             :     // Section 5: Data Representation Section
    1316          12 :     WriteUInt32(m_fp, nTemplateNumber == GS5_CMPLX ? 47 : 49);  // section size
    1317          12 :     WriteByte(m_fp, 5);  // section number
    1318          12 :     WriteUInt32(m_fp, m_nDataPoints);
    1319          12 :     WriteUInt16(m_fp, nTemplateNumber);
    1320             :     float fRefValue;
    1321          12 :     memcpy(&fRefValue, &idrstmpl[TMPL5_R_IDX], 4);
    1322          12 :     WriteFloat32(m_fp, fRefValue);
    1323          12 :     WriteInt16(m_fp, idrstmpl[TMPL5_E_IDX]);
    1324          12 :     WriteInt16(m_fp, idrstmpl[TMPL5_D_IDX]);
    1325          12 :     WriteByte(m_fp, idrstmpl[TMPL5_NBITS_IDX]);
    1326             :     // Type of original data: 0=Floating, 1=Integer
    1327          12 :     WriteByte(m_fp, GDALDataTypeIsFloating(m_eDT) ? 0 : 1);
    1328          12 :     WriteByte(m_fp, idrstmpl[TMPL5_GROUP_SPLITTING_IDX]);
    1329          12 :     WriteByte(m_fp, idrstmpl[TMPL5_MISSING_VALUE_MGNT_IDX]);
    1330          12 :     WriteComplexPackingNoData();
    1331          12 :     WriteUInt32(m_fp, GRIB2MISSING_u4);
    1332          12 :     WriteUInt32(m_fp, idrstmpl[TMPL5_NG_IDX]);
    1333          12 :     WriteByte(m_fp, idrstmpl[TMPL5_REF_GROUP_WIDTHS_IDX]);
    1334          12 :     WriteByte(m_fp, idrstmpl[TMPL5_NBITS_GROUP_WIDTHS_IDX]);
    1335          12 :     WriteUInt32(m_fp, idrstmpl[TMPL5_REF_GROUP_LENGTHS_IDX]);
    1336          12 :     WriteByte(m_fp, idrstmpl[TMPL5_LENGTH_INCR_GROUP_LENGTHS_IDX]);
    1337          12 :     WriteUInt32(m_fp, idrstmpl[TMPL5_TRUE_LENGTH_LAST_GROUP_IDX]);
    1338          12 :     WriteByte(m_fp, idrstmpl[TMPL5_NBITS_SCALED_GROUP_LENGTHS_IDX]);
    1339          12 :     if (nTemplateNumber == GS5_CMPLXSEC)
    1340             :     {
    1341           3 :         WriteByte(m_fp, nSpatialDifferencingOrder);
    1342           3 :         WriteByte(m_fp, idrstmpl[TMPL5_NB_OCTETS_EXTRA_DESCR_IDX]);
    1343             :     }
    1344             : 
    1345             :     // Section 6: Bitmap section
    1346          12 :     WriteUInt32(m_fp, 6);              // section size
    1347          12 :     WriteByte(m_fp, 6);                // section number
    1348          12 :     WriteByte(m_fp, GRIB2MISSING_u1);  // no bitmap
    1349             : 
    1350             :     // Section 7: Data Section
    1351          12 :     WriteUInt32(m_fp, 5 + nLengthPacked);  // section size
    1352          12 :     WriteByte(m_fp, 7);                    // section number
    1353          12 :     if (static_cast<int>(VSIFWriteL(pabyData, 1, nLengthPacked, m_fp)) !=
    1354             :         nLengthPacked)
    1355             :     {
    1356           0 :         VSIFree(pafData);
    1357           0 :         VSIFree(pabyData);
    1358           0 :         return false;
    1359             :     }
    1360             : 
    1361          12 :     VSIFree(pafData);
    1362          12 :     VSIFree(pabyData);
    1363             : 
    1364          12 :     return true;
    1365             : }
    1366             : 
    1367             : /************************************************************************/
    1368             : /*                             WriteIEEE()                              */
    1369             : /************************************************************************/
    1370             : 
    1371             : // See http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_temp5-4.shtml
    1372          30 : bool GRIB2Section567Writer::WriteIEEE(GDALProgressFunc pfnProgress,
    1373             :                                       void *pProgressData)
    1374             : {
    1375             :     GDALDataType eReqDT;
    1376          30 :     if (GDALGetDataTypeSizeBytes(m_eDT) <= 2 || m_eDT == GDT_Float32)
    1377          11 :         eReqDT = GDT_Float32;
    1378             :     else
    1379          19 :         eReqDT = GDT_Float64;
    1380             : 
    1381             :     // Section 5: Data Representation Section
    1382          30 :     WriteUInt32(m_fp, 12);  // section size
    1383          30 :     WriteByte(m_fp, 5);     // section number
    1384          30 :     WriteUInt32(m_fp, m_nDataPoints);
    1385          30 :     WriteUInt16(m_fp, GS5_IEEE);
    1386          30 :     WriteByte(m_fp, (eReqDT == GDT_Float32) ? 1 : 2);  // Precision
    1387             : 
    1388             :     // Section 6: Bitmap section
    1389          30 :     WriteUInt32(m_fp, 6);              // section size
    1390          30 :     WriteByte(m_fp, 6);                // section number
    1391          30 :     WriteByte(m_fp, GRIB2MISSING_u1);  // no bitmap
    1392             : 
    1393             :     // Section 7: Data Section
    1394             :     const size_t nBufferSize =
    1395          30 :         static_cast<size_t>(m_nXSize) * GDALGetDataTypeSizeBytes(eReqDT);
    1396             :     // section size
    1397          30 :     WriteUInt32(m_fp, static_cast<GUInt32>(5 + nBufferSize * m_nYSize));
    1398          30 :     WriteByte(m_fp, 7);  // section number
    1399          30 :     void *pData = CPLMalloc(nBufferSize);
    1400          30 :     const int nDenominator = std::max(1, m_poSrcDS->GetRasterCount());
    1401          60 :     void *pScaledProgressData = GDALCreateScaledProgress(
    1402          30 :         static_cast<double>(m_nBand - 1) / nDenominator,
    1403          30 :         static_cast<double>(m_nBand) / nDenominator, pfnProgress,
    1404             :         pProgressData);
    1405         486 :     for (int i = 0; i < m_nYSize; i++)
    1406             :     {
    1407         456 :         int iSrcLine = m_gt[5] < 0 ? m_nYSize - 1 - i : i;
    1408         456 :         CPLErr eErr = m_poSrcDS->GetRasterBand(m_nBand)->RasterIO(
    1409         456 :             GF_Read, m_nSplitAndSwap, iSrcLine, m_nXSize - m_nSplitAndSwap, 1,
    1410         456 :             pData, m_nXSize - m_nSplitAndSwap, 1, eReqDT, 0, 0, nullptr);
    1411         456 :         if (eErr != CE_None)
    1412             :         {
    1413           0 :             CPLFree(pData);
    1414           0 :             GDALDestroyScaledProgress(pScaledProgressData);
    1415           0 :             return false;
    1416             :         }
    1417         456 :         if (m_nSplitAndSwap > 0)
    1418             :         {
    1419          54 :             eErr = m_poSrcDS->GetRasterBand(m_nBand)->RasterIO(
    1420             :                 GF_Read, 0, iSrcLine, m_nSplitAndSwap, 1,
    1421             :                 reinterpret_cast<void *>(reinterpret_cast<GByte *>(pData) +
    1422         108 :                                          (m_nXSize - m_nSplitAndSwap) *
    1423          54 :                                              GDALGetDataTypeSizeBytes(eReqDT)),
    1424             :                 m_nSplitAndSwap, 1, eReqDT, 0, 0, nullptr);
    1425          54 :             if (eErr != CE_None)
    1426             :             {
    1427           0 :                 CPLFree(pData);
    1428           0 :                 GDALDestroyScaledProgress(pScaledProgressData);
    1429           0 :                 return false;
    1430             :             }
    1431             :         }
    1432         456 :         if (m_fValOffset != 0.0)
    1433             :         {
    1434           6 :             if (eReqDT == GDT_Float32)
    1435             :             {
    1436          12 :                 for (int j = 0; j < m_nXSize; j++)
    1437             :                 {
    1438           8 :                     static_cast<float *>(pData)[j] += m_fValOffset;
    1439             :                 }
    1440             :             }
    1441             :             else
    1442             :             {
    1443           6 :                 for (int j = 0; j < m_nXSize; j++)
    1444             :                 {
    1445           4 :                     static_cast<double *>(pData)[j] += m_fValOffset;
    1446             :                 }
    1447             :             }
    1448             :         }
    1449             : #ifdef CPL_LSB
    1450         456 :         GDALSwapWords(pData, GDALGetDataTypeSizeBytes(eReqDT), m_nXSize,
    1451             :                       GDALGetDataTypeSizeBytes(eReqDT));
    1452             : #endif
    1453         456 :         if (VSIFWriteL(pData, 1, nBufferSize, m_fp) != nBufferSize)
    1454             :         {
    1455           0 :             CPLFree(pData);
    1456           0 :             GDALDestroyScaledProgress(pScaledProgressData);
    1457           0 :             return false;
    1458             :         }
    1459         456 :         if (!GDALScaledProgress(static_cast<double>(i + 1) / m_nYSize, nullptr,
    1460             :                                 pScaledProgressData))
    1461             :         {
    1462           0 :             CPLFree(pData);
    1463           0 :             GDALDestroyScaledProgress(pScaledProgressData);
    1464           0 :             return false;
    1465             :         }
    1466             :     }
    1467          30 :     GDALDestroyScaledProgress(pScaledProgressData);
    1468          30 :     CPLFree(pData);
    1469             : 
    1470          30 :     return true;
    1471             : }
    1472             : 
    1473             : /************************************************************************/
    1474             : /*                        WrapArrayAsMemDataset()                       */
    1475             : /************************************************************************/
    1476             : 
    1477          28 : static GDALDataset *WrapArrayAsMemDataset(int nXSize, int nYSize,
    1478             :                                           GDALDataType eReducedDT, void *pData)
    1479             : {
    1480          28 :     CPLAssert(eReducedDT == GDT_Byte || eReducedDT == GDT_UInt16);
    1481             :     auto poMEMDS =
    1482          28 :         MEMDataset::Create("", nXSize, nYSize, 0, eReducedDT, nullptr);
    1483          28 :     GByte *pabyData = static_cast<GByte *>(pData);
    1484             : #ifdef CPL_MSB
    1485             :     if (eReducedDT == GDT_Byte)
    1486             :         pabyData++;
    1487             : #endif
    1488             :     auto hBand =
    1489          28 :         MEMCreateRasterBandEx(poMEMDS, 1, pabyData, eReducedDT, 2, 0, false);
    1490          28 :     poMEMDS->AddMEMBand(hBand);
    1491          28 :     return poMEMDS;
    1492             : }
    1493             : 
    1494             : /************************************************************************/
    1495             : /*                      GetRoundedToUpperPowerOfTwo()                   */
    1496             : /************************************************************************/
    1497             : 
    1498          13 : static int GetRoundedToUpperPowerOfTwo(int nBits)
    1499             : {
    1500          13 :     if (nBits == 3)
    1501           1 :         nBits = 4;
    1502          12 :     else if (nBits > 4 && nBits < 8)
    1503           2 :         nBits = 8;
    1504          10 :     else if (nBits > 8 && nBits < 15)
    1505           1 :         nBits = 16;
    1506          13 :     return nBits;
    1507             : }
    1508             : 
    1509             : /************************************************************************/
    1510             : /*                             GetScaledData()                          */
    1511             : /************************************************************************/
    1512             : 
    1513          28 : static GUInt16 *GetScaledData(GUInt32 nDataPoints, const float *pafData,
    1514             :                               float fMin, float fMax, double dfDecimalScale,
    1515             :                               double dfMinScaled,
    1516             :                               bool bOnlyPowerOfTwoDepthAllowed, int &nBits,
    1517             :                               GInt16 &nBinaryScaleFactor)
    1518             : {
    1519          28 :     bool bDone = false;
    1520          28 :     nBinaryScaleFactor = 0;
    1521             :     GUInt16 *panData = static_cast<GUInt16 *>(
    1522          28 :         VSI_MALLOC2_VERBOSE(nDataPoints, sizeof(GUInt16)));
    1523          28 :     if (panData == nullptr)
    1524             :     {
    1525           0 :         return nullptr;
    1526             :     }
    1527             : 
    1528          28 :     const double dfScaledMaxDiff = (fMax - fMin) * dfDecimalScale;
    1529          28 :     if (nBits == 0)
    1530             :     {
    1531          12 :         nBits = (g2int)ceil(log(ceil(dfScaledMaxDiff)) / log(2.0));
    1532          12 :         if (nBits > 16)
    1533             :         {
    1534           2 :             CPLError(CE_Warning, CPLE_AppDefined,
    1535             :                      "More than 16 bits of integer precision would be "
    1536             :                      "required. Dropping precision to fit on 16 bits");
    1537           2 :             nBits = 16;
    1538             :         }
    1539             :         else
    1540             :         {
    1541          10 :             bDone = true;
    1542      527498 :             for (GUInt32 i = 0; i < nDataPoints; i++)
    1543             :             {
    1544      527488 :                 panData[i] = static_cast<GUInt16>(
    1545      527488 :                     0.5 + (pafData[i] * dfDecimalScale - dfMinScaled));
    1546             :             }
    1547             :         }
    1548             :     }
    1549             : 
    1550          28 :     if (bOnlyPowerOfTwoDepthAllowed)
    1551          13 :         nBits = GetRoundedToUpperPowerOfTwo(nBits);
    1552             : 
    1553          28 :     if (!bDone && nBits != 0)
    1554             :     {
    1555          18 :         if (nBits > 16)
    1556             :         {
    1557           0 :             CPLError(CE_Warning, CPLE_AppDefined,
    1558             :                      "Maximum bit depth supported is 16. Using that");
    1559           0 :             nBits = 16;
    1560             :         }
    1561          18 :         const int nMaxNum = (1 << nBits) - 1;
    1562          18 :         double dfTemp = log(nMaxNum / dfScaledMaxDiff) / log(2.0);
    1563          18 :         nBinaryScaleFactor = static_cast<GInt16>(ceil(-dfTemp));
    1564          18 :         double dfBinaryScale = pow(2.0, -1.0 * nBinaryScaleFactor);
    1565        5634 :         for (GUInt32 i = 0; i < nDataPoints; i++)
    1566             :         {
    1567        5616 :             panData[i] = static_cast<GUInt16>(
    1568        5616 :                 0.5 +
    1569        5616 :                 (pafData[i] * dfDecimalScale - dfMinScaled) * dfBinaryScale);
    1570             :         }
    1571             :     }
    1572             : 
    1573          28 :     return panData;
    1574             : }
    1575             : 
    1576             : /************************************************************************/
    1577             : /*                              WritePNG()                              */
    1578             : /************************************************************************/
    1579             : 
    1580             : // See http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_temp5-41.shtml
    1581          14 : bool GRIB2Section567Writer::WritePNG()
    1582             : {
    1583          14 :     float *pafData = GetFloatData();
    1584          14 :     if (pafData == nullptr)
    1585           0 :         return false;
    1586             : 
    1587          14 :     if (m_bUseZeroBits)
    1588             :     {
    1589             :         // Section 5: Data Representation Section
    1590           1 :         WriteUInt32(m_fp, 21);  // section size
    1591           1 :         WriteByte(m_fp, 5);     // section number
    1592           1 :         WriteUInt32(m_fp, m_nDataPoints);
    1593           1 :         WriteUInt16(m_fp, GS5_PNG);
    1594           1 :         WriteFloat32(
    1595             :             m_fp,
    1596           1 :             static_cast<float>(m_dfMinScaled / m_dfDecimalScale));  // ref value
    1597           1 :         WriteInt16(m_fp, 0);  // Binary scale factor (E)
    1598           1 :         WriteInt16(m_fp, 0);  // Decimal scale factor (D)
    1599           1 :         WriteByte(m_fp, 0);   // Number of bits
    1600             :         // Type of original data: 0=Floating, 1=Integer
    1601           1 :         WriteByte(m_fp, GDALDataTypeIsFloating(m_eDT) ? 0 : 1);
    1602             : 
    1603             :         // Section 6: Bitmap section
    1604           1 :         WriteUInt32(m_fp, 6);              // section size
    1605           1 :         WriteByte(m_fp, 6);                // section number
    1606           1 :         WriteByte(m_fp, GRIB2MISSING_u1);  // no bitmap
    1607             : 
    1608             :         // Section 7: Data Section
    1609           1 :         WriteUInt32(m_fp, 5);  // section size
    1610           1 :         WriteByte(m_fp, 7);    // section number
    1611             : 
    1612           1 :         CPLFree(pafData);
    1613             : 
    1614           1 :         return true;
    1615             :     }
    1616             : 
    1617             :     GDALDriver *poPNGDriver =
    1618          13 :         reinterpret_cast<GDALDriver *>(GDALGetDriverByName("PNG"));
    1619          13 :     if (poPNGDriver == nullptr)
    1620             :     {
    1621           0 :         CPLError(CE_Failure, CPLE_AppDefined, "Cannot find PNG driver");
    1622           0 :         return false;
    1623             :     }
    1624             : 
    1625          13 :     GInt16 nBinaryScaleFactor = 0;
    1626             :     GUInt16 *panData =
    1627          26 :         GetScaledData(m_nDataPoints, pafData, m_fMin, m_fMax, m_dfDecimalScale,
    1628          13 :                       m_dfMinScaled, true, m_nBits, nBinaryScaleFactor);
    1629          13 :     if (panData == nullptr)
    1630             :     {
    1631           0 :         VSIFree(pafData);
    1632           0 :         return false;
    1633             :     }
    1634             : 
    1635          13 :     CPLFree(pafData);
    1636             : 
    1637          26 :     CPLStringList aosPNGOptions;
    1638          13 :     aosPNGOptions.SetNameValue("NBITS", CPLSPrintf("%d", m_nBits));
    1639             : 
    1640          13 :     const GDALDataType eReducedDT = (m_nBits <= 8) ? GDT_Byte : GDT_UInt16;
    1641             :     GDALDataset *poMEMDS =
    1642          13 :         WrapArrayAsMemDataset(m_nXSize, m_nYSize, eReducedDT, panData);
    1643             : 
    1644          26 :     const CPLString osTmpFile(VSIMemGenerateHiddenFilename("grib_driver.png"));
    1645          13 :     GDALDataset *poPNGDS = poPNGDriver->CreateCopy(
    1646          13 :         osTmpFile, poMEMDS, FALSE, aosPNGOptions.List(), nullptr, nullptr);
    1647          13 :     if (poPNGDS == nullptr)
    1648             :     {
    1649           0 :         CPLError(CE_Failure, CPLE_AppDefined, "PNG compression failed");
    1650           0 :         VSIUnlink(osTmpFile);
    1651           0 :         delete poMEMDS;
    1652           0 :         CPLFree(panData);
    1653           0 :         return false;
    1654             :     }
    1655          13 :     delete poPNGDS;
    1656          13 :     delete poMEMDS;
    1657          13 :     CPLFree(panData);
    1658             : 
    1659             :     // Section 5: Data Representation Section
    1660          13 :     WriteUInt32(m_fp, 21);  // section size
    1661          13 :     WriteByte(m_fp, 5);     // section number
    1662          13 :     WriteUInt32(m_fp, m_nDataPoints);
    1663          13 :     WriteUInt16(m_fp, GS5_PNG);
    1664          13 :     WriteFloat32(m_fp, static_cast<float>(m_dfMinScaled));
    1665          13 :     WriteInt16(m_fp, nBinaryScaleFactor);     // Binary scale factor (E)
    1666          13 :     WriteInt16(m_fp, m_nDecimalScaleFactor);  // Decimal scale factor (D)
    1667          13 :     WriteByte(m_fp, m_nBits);                 // Number of bits
    1668             :     // Type of original data: 0=Floating, 1=Integer
    1669          13 :     WriteByte(m_fp, GDALDataTypeIsFloating(m_eDT) ? 0 : 1);
    1670             : 
    1671             :     // Section 6: Bitmap section
    1672          13 :     WriteUInt32(m_fp, 6);              // section size
    1673          13 :     WriteByte(m_fp, 6);                // section number
    1674          13 :     WriteByte(m_fp, GRIB2MISSING_u1);  // no bitmap
    1675             : 
    1676             :     // Section 7: Data Section
    1677          13 :     vsi_l_offset nDataLength = 0;
    1678          13 :     GByte *pabyData = VSIGetMemFileBuffer(osTmpFile, &nDataLength, FALSE);
    1679          13 :     WriteUInt32(m_fp, static_cast<GUInt32>(5 + nDataLength));  // section size
    1680          13 :     WriteByte(m_fp, 7);                                        // section number
    1681          13 :     const size_t nDataLengthSize = static_cast<size_t>(nDataLength);
    1682             :     const bool bOK =
    1683          13 :         VSIFWriteL(pabyData, 1, nDataLengthSize, m_fp) == nDataLengthSize;
    1684             : 
    1685          13 :     VSIUnlink(osTmpFile);
    1686          13 :     VSIUnlink((osTmpFile + ".aux.xml").c_str());
    1687             : 
    1688          13 :     return bOK;
    1689             : }
    1690             : 
    1691             : /************************************************************************/
    1692             : /*                             WriteJPEG2000()                          */
    1693             : /************************************************************************/
    1694             : 
    1695             : // See http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_temp5-40.shtml
    1696          17 : bool GRIB2Section567Writer::WriteJPEG2000(char **papszOptions)
    1697             : {
    1698          17 :     float *pafData = GetFloatData();
    1699          17 :     if (pafData == nullptr)
    1700           0 :         return false;
    1701             : 
    1702          17 :     if (m_bUseZeroBits)
    1703             :     {
    1704             :         // Section 5: Data Representation Section
    1705           1 :         WriteUInt32(m_fp, 23);  // section size
    1706           1 :         WriteByte(m_fp, 5);     // section number
    1707           1 :         WriteUInt32(m_fp, m_nDataPoints);
    1708           1 :         WriteUInt16(m_fp, GS5_JPEG2000);
    1709           1 :         WriteFloat32(
    1710             :             m_fp,
    1711           1 :             static_cast<float>(m_dfMinScaled / m_dfDecimalScale));  // ref val
    1712           1 :         WriteInt16(m_fp, 0);  // Binary scale factor (E)
    1713           1 :         WriteInt16(m_fp, 0);  // Decimal scale factor (D)
    1714           1 :         WriteByte(m_fp, 0);   // Number of bits
    1715             :         // Type of original data: 0=Floating, 1=Integer
    1716           1 :         WriteByte(m_fp, GDALDataTypeIsFloating(m_eDT) ? 0 : 1);
    1717           1 :         WriteByte(m_fp, 0);                // compression type: lossless
    1718           1 :         WriteByte(m_fp, GRIB2MISSING_u1);  // compression ratio
    1719             : 
    1720             :         // Section 6: Bitmap section
    1721           1 :         WriteUInt32(m_fp, 6);              // section size
    1722           1 :         WriteByte(m_fp, 6);                // section number
    1723           1 :         WriteByte(m_fp, GRIB2MISSING_u1);  // no bitmap
    1724             : 
    1725             :         // Section 7: Data Section
    1726           1 :         WriteUInt32(m_fp, 5);  // section size
    1727           1 :         WriteByte(m_fp, 7);    // section number
    1728             : 
    1729           1 :         CPLFree(pafData);
    1730             : 
    1731           1 :         return true;
    1732             :     }
    1733             : 
    1734          16 :     GDALDriver *poJ2KDriver = nullptr;
    1735          16 :     const char *pszJ2KDriver = GetBandOption(papszOptions, nullptr, m_nBand,
    1736             :                                              "JPEG2000_DRIVER", nullptr);
    1737          16 :     if (pszJ2KDriver)
    1738             :     {
    1739             :         poJ2KDriver =
    1740           6 :             reinterpret_cast<GDALDriver *>(GDALGetDriverByName(pszJ2KDriver));
    1741             :     }
    1742             :     else
    1743             :     {
    1744          20 :         for (size_t i = 0; i < CPL_ARRAYSIZE(apszJ2KDrivers); i++)
    1745             :         {
    1746             :             poJ2KDriver = reinterpret_cast<GDALDriver *>(
    1747          20 :                 GDALGetDriverByName(apszJ2KDrivers[i]));
    1748          20 :             if (poJ2KDriver)
    1749             :             {
    1750          10 :                 CPLDebug("GRIB", "Using %s", poJ2KDriver->GetDescription());
    1751          10 :                 break;
    1752             :             }
    1753             :         }
    1754             :     }
    1755          16 :     if (poJ2KDriver == nullptr)
    1756             :     {
    1757           1 :         CPLError(CE_Failure, CPLE_AppDefined, "Cannot find JPEG2000 driver");
    1758           1 :         VSIFree(pafData);
    1759           1 :         return false;
    1760             :     }
    1761             : 
    1762          15 :     GInt16 nBinaryScaleFactor = 0;
    1763             :     GUInt16 *panData =
    1764          30 :         GetScaledData(m_nDataPoints, pafData, m_fMin, m_fMax, m_dfDecimalScale,
    1765          15 :                       m_dfMinScaled, false, m_nBits, nBinaryScaleFactor);
    1766          15 :     if (panData == nullptr)
    1767             :     {
    1768           0 :         VSIFree(pafData);
    1769           0 :         return false;
    1770             :     }
    1771             : 
    1772          15 :     CPLFree(pafData);
    1773             : 
    1774          30 :     CPLStringList aosJ2KOptions;
    1775          15 :     int nCompressionRatio = atoi(GetBandOption(papszOptions, nullptr, m_nBand,
    1776             :                                                "COMPRESSION_RATIO", "1"));
    1777          15 :     if (m_nDataPoints < 10000 && nCompressionRatio > 1)
    1778             :     {
    1779             :         // Lossy compression with too few pixels is really lossy due to how
    1780             :         // codec work
    1781           1 :         CPLDebug("GRIB", "Forcing JPEG2000 lossless mode given "
    1782             :                          "the low number of pixels");
    1783           1 :         nCompressionRatio = 1;
    1784             :     }
    1785          15 :     const bool bLossLess = nCompressionRatio <= 1;
    1786          15 :     if (EQUAL(poJ2KDriver->GetDescription(), "JP2KAK"))
    1787             :     {
    1788           0 :         if (bLossLess)
    1789             :         {
    1790           0 :             aosJ2KOptions.SetNameValue("QUALITY", "100");
    1791             :         }
    1792             :         else
    1793             :         {
    1794             :             aosJ2KOptions.SetNameValue(
    1795             :                 "QUALITY",
    1796           0 :                 CPLSPrintf("%d", std::max(1, 100 / nCompressionRatio)));
    1797             :         }
    1798             :     }
    1799          15 :     else if (EQUAL(poJ2KDriver->GetDescription(), "JP2OPENJPEG"))
    1800             :     {
    1801          12 :         if (bLossLess)
    1802             :         {
    1803          11 :             aosJ2KOptions.SetNameValue("QUALITY", "100");
    1804          11 :             aosJ2KOptions.SetNameValue("REVERSIBLE", "YES");
    1805             :         }
    1806             :         else
    1807             :         {
    1808             :             aosJ2KOptions.SetNameValue(
    1809           1 :                 "QUALITY", CPLSPrintf("%f", 100.0 / nCompressionRatio));
    1810             :         }
    1811             :     }
    1812           3 :     else if (EQUAL(poJ2KDriver->GetDescription(), "JP2ECW"))
    1813             :     {
    1814           2 :         if (bLossLess)
    1815             :         {
    1816           1 :             aosJ2KOptions.SetNameValue("TARGET", "0");
    1817             :         }
    1818             :         else
    1819             :         {
    1820             :             aosJ2KOptions.SetNameValue(
    1821           1 :                 "TARGET", CPLSPrintf("%f", 100.0 - 100.0 / nCompressionRatio));
    1822             :         }
    1823             :     }
    1824          15 :     aosJ2KOptions.SetNameValue("NBITS", CPLSPrintf("%d", m_nBits));
    1825             : 
    1826          15 :     const GDALDataType eReducedDT = (m_nBits <= 8) ? GDT_Byte : GDT_UInt16;
    1827             :     GDALDataset *poMEMDS =
    1828          15 :         WrapArrayAsMemDataset(m_nXSize, m_nYSize, eReducedDT, panData);
    1829             : 
    1830          30 :     const CPLString osTmpFile(VSIMemGenerateHiddenFilename("grib_driver.j2k"));
    1831          15 :     GDALDataset *poJ2KDS = poJ2KDriver->CreateCopy(
    1832          15 :         osTmpFile, poMEMDS, FALSE, aosJ2KOptions.List(), nullptr, nullptr);
    1833          15 :     if (poJ2KDS == nullptr)
    1834             :     {
    1835           1 :         CPLError(CE_Failure, CPLE_AppDefined, "JPEG2000 compression failed");
    1836           1 :         VSIUnlink(osTmpFile);
    1837           1 :         delete poMEMDS;
    1838           1 :         CPLFree(panData);
    1839           1 :         return false;
    1840             :     }
    1841          14 :     delete poJ2KDS;
    1842          14 :     delete poMEMDS;
    1843          14 :     CPLFree(panData);
    1844             : 
    1845             :     // Section 5: Data Representation Section
    1846          14 :     WriteUInt32(m_fp, 23);  // section size
    1847          14 :     WriteByte(m_fp, 5);     // section number
    1848          14 :     WriteUInt32(m_fp, m_nDataPoints);
    1849          14 :     WriteUInt16(m_fp, GS5_JPEG2000);
    1850          14 :     WriteFloat32(m_fp, static_cast<float>(m_dfMinScaled));
    1851          14 :     WriteInt16(m_fp, nBinaryScaleFactor);     // Binary scale factor (E)
    1852          14 :     WriteInt16(m_fp, m_nDecimalScaleFactor);  // Decimal scale factor (D)
    1853          14 :     WriteByte(m_fp, m_nBits);                 // Number of bits
    1854             :     // Type of original data: 0=Floating, 1=Integer
    1855          14 :     WriteByte(m_fp, GDALDataTypeIsFloating(m_eDT) ? 0 : 1);
    1856             :     // compression type: lossless(0) or lossy(1)
    1857          14 :     WriteByte(m_fp, bLossLess ? 0 : 1);
    1858          14 :     WriteByte(m_fp, bLossLess ? GRIB2MISSING_u1
    1859             :                               : nCompressionRatio);  // compression ratio
    1860             : 
    1861             :     // Section 6: Bitmap section
    1862          14 :     WriteUInt32(m_fp, 6);              // section size
    1863          14 :     WriteByte(m_fp, 6);                // section number
    1864          14 :     WriteByte(m_fp, GRIB2MISSING_u1);  // no bitmap
    1865             : 
    1866             :     // Section 7: Data Section
    1867          14 :     vsi_l_offset nDataLength = 0;
    1868          14 :     GByte *pabyData = VSIGetMemFileBuffer(osTmpFile, &nDataLength, FALSE);
    1869          14 :     WriteUInt32(m_fp, static_cast<GUInt32>(5 + nDataLength));  // section size
    1870          14 :     WriteByte(m_fp, 7);                                        // section number
    1871          14 :     const size_t nDataLengthSize = static_cast<size_t>(nDataLength);
    1872             :     const bool bOK =
    1873          14 :         VSIFWriteL(pabyData, 1, nDataLengthSize, m_fp) == nDataLengthSize;
    1874             : 
    1875          14 :     VSIUnlink(osTmpFile);
    1876          14 :     VSIUnlink((osTmpFile + ".aux.xml").c_str());
    1877             : 
    1878          14 :     return bOK;
    1879             : }
    1880             : 
    1881             : /************************************************************************/
    1882             : /*                               Write()                                */
    1883             : /************************************************************************/
    1884             : 
    1885         156 : bool GRIB2Section567Writer::Write(float fValOffset, char **papszOptions,
    1886             :                                   GDALProgressFunc pfnProgress,
    1887             :                                   void *pProgressData)
    1888             : {
    1889         156 :     m_fValOffset = fValOffset;
    1890             : 
    1891             :     typedef enum
    1892             :     {
    1893             :         SIMPLE_PACKING,
    1894             :         COMPLEX_PACKING,
    1895             :         IEEE_FLOATING_POINT,
    1896             :         PNG,
    1897             :         JPEG2000
    1898             :     } GRIBDataEncoding;
    1899             : 
    1900         156 :     if (m_eDT != GDT_Byte && m_eDT != GDT_UInt16 && m_eDT != GDT_Int16 &&
    1901          56 :         m_eDT != GDT_UInt32 && m_eDT != GDT_Int32 && m_eDT != GDT_Float32 &&
    1902          30 :         m_eDT != GDT_Float64)
    1903             :     {
    1904           5 :         CPLError(CE_Failure, CPLE_NotSupported, "Unsupported data type: %s",
    1905             :                  GDALGetDataTypeName(m_eDT));
    1906           5 :         return false;
    1907             :     }
    1908             :     const char *pszDataEncoding =
    1909         151 :         GetBandOption(papszOptions, nullptr, m_nBand, "DATA_ENCODING", "AUTO");
    1910         151 :     GRIBDataEncoding eDataEncoding(SIMPLE_PACKING);
    1911         151 :     const char *pszJ2KDriver = GetBandOption(papszOptions, nullptr, m_nBand,
    1912             :                                              "JPEG2000_DRIVER", nullptr);
    1913         151 :     const char *pszSpatialDifferencingOrder = GetBandOption(
    1914             :         papszOptions, nullptr, m_nBand, "SPATIAL_DIFFERENCING_ORDER", nullptr);
    1915         151 :     if (pszJ2KDriver && pszSpatialDifferencingOrder)
    1916             :     {
    1917           1 :         CPLError(CE_Failure, CPLE_NotSupported,
    1918             :                  "JPEG2000_DRIVER and SPATIAL_DIFFERENCING_ORDER are not "
    1919             :                  "compatible");
    1920           1 :         return false;
    1921             :     }
    1922             : 
    1923         150 :     if (m_bHasNoData && !EQUAL(pszDataEncoding, "COMPLEX_PACKING") &&
    1924             :         pszSpatialDifferencingOrder == nullptr)
    1925             :     {
    1926             :         double *padfVals = static_cast<double *>(
    1927           6 :             VSI_MALLOC2_VERBOSE(m_nXSize, sizeof(double)));
    1928           6 :         if (padfVals == nullptr)
    1929           0 :             return false;
    1930           6 :         bool bFoundNoData = false;
    1931           7 :         for (int j = 0; j < m_nYSize; j++)
    1932             :         {
    1933           6 :             CPLErr eErr = m_poSrcDS->GetRasterBand(m_nBand)->RasterIO(
    1934             :                 GF_Read, 0, j, m_nXSize, 1, padfVals, m_nXSize, 1, GDT_Float64,
    1935             :                 0, 0, nullptr);
    1936           6 :             if (eErr != CE_None)
    1937             :             {
    1938           0 :                 VSIFree(padfVals);
    1939           0 :                 return false;
    1940             :             }
    1941           7 :             for (int i = 0; i < m_nXSize; i++)
    1942             :             {
    1943           6 :                 if (padfVals[i] == m_dfNoData)
    1944             :                 {
    1945           5 :                     bFoundNoData = true;
    1946           5 :                     break;
    1947             :                 }
    1948             :             }
    1949           6 :             if (bFoundNoData)
    1950           5 :                 break;
    1951             :         }
    1952           6 :         VSIFree(padfVals);
    1953             : 
    1954           6 :         if (!bFoundNoData)
    1955             :         {
    1956           1 :             m_bHasNoData = false;
    1957             :         }
    1958             :     }
    1959             : 
    1960         150 :     if (EQUAL(pszDataEncoding, "AUTO"))
    1961             :     {
    1962          81 :         if (m_bHasNoData || pszSpatialDifferencingOrder != nullptr)
    1963             :         {
    1964           5 :             eDataEncoding = COMPLEX_PACKING;
    1965           5 :             CPLDebug("GRIB", "Using COMPLEX_PACKING");
    1966             :         }
    1967          76 :         else if (pszJ2KDriver != nullptr)
    1968             :         {
    1969           6 :             eDataEncoding = JPEG2000;
    1970           6 :             CPLDebug("GRIB", "Using JPEG2000");
    1971             :         }
    1972          70 :         else if (m_eDT == GDT_Float32 || m_eDT == GDT_Float64)
    1973             :         {
    1974          20 :             eDataEncoding = IEEE_FLOATING_POINT;
    1975          20 :             CPLDebug("GRIB", "Using IEEE_FLOATING_POINT");
    1976             :         }
    1977             :         else
    1978             :         {
    1979          50 :             CPLDebug("GRIB", "Using SIMPLE_PACKING");
    1980             :         }
    1981             :     }
    1982          69 :     else if (EQUAL(pszDataEncoding, "SIMPLE_PACKING"))
    1983             :     {
    1984          16 :         eDataEncoding = SIMPLE_PACKING;
    1985             :     }
    1986          53 :     else if (EQUAL(pszDataEncoding, "COMPLEX_PACKING"))
    1987             :     {
    1988          17 :         eDataEncoding = COMPLEX_PACKING;
    1989             :     }
    1990          36 :     else if (EQUAL(pszDataEncoding, "IEEE_FLOATING_POINT"))
    1991             :     {
    1992          10 :         eDataEncoding = IEEE_FLOATING_POINT;
    1993             :     }
    1994          26 :     else if (EQUAL(pszDataEncoding, "PNG"))
    1995             :     {
    1996          14 :         eDataEncoding = PNG;
    1997             :     }
    1998          12 :     else if (EQUAL(pszDataEncoding, "JPEG2000"))
    1999             :     {
    2000          11 :         eDataEncoding = JPEG2000;
    2001             :     }
    2002             :     else
    2003             :     {
    2004           1 :         CPLError(CE_Failure, CPLE_NotSupported, "Unsupported DATA_ENCODING=%s",
    2005             :                  pszDataEncoding);
    2006           1 :         return false;
    2007             :     }
    2008             : 
    2009             :     const char *pszBits =
    2010         149 :         GetBandOption(papszOptions, nullptr, m_nBand, "NBITS", nullptr);
    2011         149 :     if (pszBits == nullptr && eDataEncoding != IEEE_FLOATING_POINT)
    2012             :     {
    2013         200 :         pszBits = m_poSrcDS->GetRasterBand(m_nBand)->GetMetadataItem(
    2014         100 :             "DRS_NBITS", "GRIB");
    2015             :     }
    2016          49 :     else if (pszBits != nullptr && eDataEncoding == IEEE_FLOATING_POINT)
    2017             :     {
    2018           1 :         CPLError(CE_Warning, CPLE_NotSupported,
    2019             :                  "NBITS ignored for DATA_ENCODING = IEEE_FLOATING_POINT");
    2020             :     }
    2021         149 :     if (pszBits == nullptr)
    2022             :     {
    2023         123 :         pszBits = "0";
    2024             :     }
    2025         149 :     m_nBits = std::max(0, atoi(pszBits));
    2026         149 :     if (m_nBits > 31)
    2027             :     {
    2028           1 :         CPLError(CE_Warning, CPLE_NotSupported, "NBITS clamped to 31");
    2029           1 :         m_nBits = 31;
    2030             :     }
    2031             : 
    2032         149 :     const char *pszDecimalScaleFactor = GetBandOption(
    2033             :         papszOptions, nullptr, m_nBand, "DECIMAL_SCALE_FACTOR", nullptr);
    2034         149 :     if (pszDecimalScaleFactor != nullptr)
    2035             :     {
    2036          10 :         m_nDecimalScaleFactor = atoi(pszDecimalScaleFactor);
    2037          10 :         if (m_nDecimalScaleFactor != 0 && eDataEncoding == IEEE_FLOATING_POINT)
    2038             :         {
    2039           1 :             CPLError(CE_Warning, CPLE_NotSupported,
    2040             :                      "DECIMAL_SCALE_FACTOR ignored for "
    2041             :                      "DATA_ENCODING = IEEE_FLOATING_POINT");
    2042             :         }
    2043           9 :         else if (m_nDecimalScaleFactor > 0 && !GDALDataTypeIsFloating(m_eDT))
    2044             :         {
    2045           1 :             CPLError(CE_Warning, CPLE_AppDefined,
    2046             :                      "DECIMAL_SCALE_FACTOR > 0 makes no sense for integer "
    2047             :                      "data types. Ignored");
    2048           1 :             m_nDecimalScaleFactor = 0;
    2049             :         }
    2050             :     }
    2051         139 :     else if (eDataEncoding != IEEE_FLOATING_POINT)
    2052             :     {
    2053             :         pszDecimalScaleFactor =
    2054         110 :             m_poSrcDS->GetRasterBand(m_nBand)->GetMetadataItem(
    2055         110 :                 "DRS_DECIMAL_SCALE_FACTOR", "GRIB");
    2056         110 :         if (pszDecimalScaleFactor != nullptr)
    2057             :         {
    2058           6 :             m_nDecimalScaleFactor = atoi(pszDecimalScaleFactor);
    2059             :         }
    2060             :     }
    2061         149 :     m_dfDecimalScale = pow(10.0, static_cast<double>(m_nDecimalScaleFactor));
    2062             : 
    2063         149 :     if (pszJ2KDriver != nullptr && eDataEncoding != JPEG2000)
    2064             :     {
    2065           2 :         CPLError(CE_Warning, CPLE_AppDefined,
    2066             :                  "JPEG2000_DRIVER option ignored for "
    2067             :                  "non-JPEG2000 DATA_ENCODING");
    2068             :     }
    2069         149 :     if (pszSpatialDifferencingOrder && eDataEncoding != COMPLEX_PACKING)
    2070             :     {
    2071           1 :         CPLError(CE_Warning, CPLE_AppDefined,
    2072             :                  "SPATIAL_DIFFERENCING_ORDER option ignored for "
    2073             :                  "non-COMPLEX_PACKING DATA_ENCODING");
    2074             :     }
    2075         149 :     if (m_bHasNoData && eDataEncoding != COMPLEX_PACKING)
    2076             :     {
    2077           2 :         CPLError(CE_Warning, CPLE_AppDefined,
    2078             :                  "non-COMPLEX_PACKING DATA_ENCODING cannot preserve nodata");
    2079             :     }
    2080             : 
    2081         149 :     if (eDataEncoding == SIMPLE_PACKING)
    2082             :     {
    2083          66 :         return WriteSimplePacking();
    2084             :     }
    2085          83 :     else if (eDataEncoding == COMPLEX_PACKING)
    2086             :     {
    2087          22 :         const int nSpatialDifferencingOrder =
    2088          22 :             pszSpatialDifferencingOrder ? atoi(pszSpatialDifferencingOrder) : 0;
    2089          22 :         return WriteComplexPacking(nSpatialDifferencingOrder);
    2090             :     }
    2091          61 :     else if (eDataEncoding == IEEE_FLOATING_POINT)
    2092             :     {
    2093          30 :         return WriteIEEE(pfnProgress, pProgressData);
    2094             :     }
    2095          31 :     else if (eDataEncoding == PNG)
    2096             :     {
    2097          14 :         return WritePNG();
    2098             :     }
    2099             :     else /* if( eDataEncoding == JPEG2000 ) */
    2100             :     {
    2101          17 :         return WriteJPEG2000(papszOptions);
    2102             :     }
    2103             : }
    2104             : 
    2105             : /************************************************************************/
    2106             : /*                           GetIDSOption()                             */
    2107             : /************************************************************************/
    2108             : 
    2109        1134 : static const char *GetIDSOption(char **papszOptions, GDALDataset *poSrcDS,
    2110             :                                 int nBand, const char *pszKey,
    2111             :                                 const char *pszDefault)
    2112             : {
    2113             :     const char *pszValue =
    2114        1134 :         GetBandOption(papszOptions, nullptr, nBand,
    2115        2268 :                       (CPLString("IDS_") + pszKey).c_str(), nullptr);
    2116        1134 :     if (pszValue == nullptr)
    2117             :     {
    2118             :         const char *pszIDS =
    2119        1133 :             GetBandOption(papszOptions, poSrcDS, nBand, "IDS", nullptr);
    2120        1133 :         if (pszIDS != nullptr)
    2121             :         {
    2122         153 :             char **papszTokens = CSLTokenizeString2(pszIDS, " ", 0);
    2123         153 :             pszValue = CSLFetchNameValue(papszTokens, pszKey);
    2124         153 :             if (pszValue)
    2125         142 :                 pszValue = CPLSPrintf("%s", pszValue);
    2126         153 :             CSLDestroy(papszTokens);
    2127             :         }
    2128             :     }
    2129        1134 :     if (pszValue == nullptr)
    2130         991 :         pszValue = pszDefault;
    2131        1134 :     return pszValue;
    2132             : }
    2133             : 
    2134             : /************************************************************************/
    2135             : /*                           WriteSection1()                            */
    2136             : /************************************************************************/
    2137             : 
    2138         162 : static void WriteSection1(VSILFILE *fp, GDALDataset *poSrcDS, int nBand,
    2139             :                           char **papszOptions)
    2140             : {
    2141             :     // Section 1: Identification Section
    2142         162 :     WriteUInt32(fp, 21);  // section size
    2143         162 :     WriteByte(fp, 1);     // section number
    2144             : 
    2145             :     GUInt16 nCenter = static_cast<GUInt16>(
    2146         162 :         atoi(GetIDSOption(papszOptions, poSrcDS, nBand, "CENTER",
    2147         162 :                           CPLSPrintf("%d", GRIB2MISSING_u1))));
    2148         162 :     WriteUInt16(fp, nCenter);
    2149             : 
    2150             :     GUInt16 nSubCenter = static_cast<GUInt16>(
    2151         162 :         atoi(GetIDSOption(papszOptions, poSrcDS, nBand, "SUBCENTER",
    2152         162 :                           CPLSPrintf("%d", GRIB2MISSING_u2))));
    2153         162 :     WriteUInt16(fp, nSubCenter);
    2154             : 
    2155             :     GByte nMasterTable = static_cast<GByte>(
    2156         162 :         atoi(GetIDSOption(papszOptions, poSrcDS, nBand, "MASTER_TABLE", "2")));
    2157         162 :     WriteByte(fp, nMasterTable);
    2158             : 
    2159         162 :     WriteByte(fp, 0);  // local table
    2160             : 
    2161         162 :     GByte nSignfRefTime = static_cast<GByte>(atoi(
    2162         162 :         GetIDSOption(papszOptions, poSrcDS, nBand, "SIGNF_REF_TIME", "0")));
    2163         162 :     WriteByte(fp, nSignfRefTime);  // Significance of reference time
    2164             : 
    2165             :     const char *pszRefTime =
    2166         162 :         GetIDSOption(papszOptions, poSrcDS, nBand, "REF_TIME", "");
    2167         162 :     int nYear = 1970, nMonth = 1, nDay = 1, nHour = 0, nMinute = 0, nSecond = 0;
    2168         162 :     sscanf(pszRefTime, "%04d-%02d-%02dT%02d:%02d:%02dZ", &nYear, &nMonth, &nDay,
    2169             :            &nHour, &nMinute, &nSecond);
    2170         162 :     WriteUInt16(fp, nYear);
    2171         162 :     WriteByte(fp, nMonth);
    2172         162 :     WriteByte(fp, nDay);
    2173         162 :     WriteByte(fp, nHour);
    2174         162 :     WriteByte(fp, nMinute);
    2175         162 :     WriteByte(fp, nSecond);
    2176             : 
    2177             :     GByte nProdStatus = static_cast<GByte>(
    2178         162 :         atoi(GetIDSOption(papszOptions, poSrcDS, nBand, "PROD_STATUS",
    2179         162 :                           CPLSPrintf("%d", GRIB2MISSING_u1))));
    2180         162 :     WriteByte(fp, nProdStatus);
    2181             : 
    2182             :     GByte nType = static_cast<GByte>(
    2183         162 :         atoi(GetIDSOption(papszOptions, poSrcDS, nBand, "TYPE",
    2184         162 :                           CPLSPrintf("%d", GRIB2MISSING_u1))));
    2185         162 :     WriteByte(fp, nType);
    2186         162 : }
    2187             : 
    2188             : /************************************************************************/
    2189             : /*                        WriteAssembledPDS()                           */
    2190             : /************************************************************************/
    2191             : 
    2192          12 : static void WriteAssembledPDS(VSILFILE *fp, const gtemplate *mappds,
    2193             :                               bool bWriteExt, char **papszTokens,
    2194             :                               std::vector<int> &anVals)
    2195             : {
    2196          12 :     const int iStart = bWriteExt ? mappds->maplen : 0;
    2197          12 :     const int iEnd =
    2198          12 :         bWriteExt ? mappds->maplen + mappds->extlen : mappds->maplen;
    2199         173 :     for (int i = iStart; i < iEnd; i++)
    2200             :     {
    2201         161 :         const int nVal = atoi(papszTokens[i]);
    2202         161 :         anVals.push_back(nVal);
    2203         161 :         const int nEltSize =
    2204         161 :             bWriteExt ? mappds->ext[i - mappds->maplen] : mappds->map[i];
    2205         161 :         if (nEltSize == 1)
    2206             :         {
    2207          90 :             if (nVal < 0 || nVal > 255)
    2208             :             {
    2209           2 :                 CPLError(CE_Warning, CPLE_AppDefined,
    2210             :                          "Value %d of index %d in PDS should be in [0,255] "
    2211             :                          "range",
    2212             :                          nVal, i);
    2213             :             }
    2214          90 :             WriteByte(fp, nVal);
    2215             :         }
    2216          71 :         else if (nEltSize == 2)
    2217             :         {
    2218          25 :             if (nVal < 0 || nVal > 65535)
    2219             :             {
    2220           2 :                 CPLError(CE_Warning, CPLE_AppDefined,
    2221             :                          "Value %d of index %d in PDS should be in [0,65535] "
    2222             :                          "range",
    2223             :                          nVal, i);
    2224             :             }
    2225          25 :             WriteUInt16(fp, nVal);
    2226             :         }
    2227          46 :         else if (nEltSize == 4)
    2228             :         {
    2229           6 :             GIntBig nBigVal = CPLAtoGIntBig(papszTokens[i]);
    2230           6 :             anVals.back() = static_cast<int>(nBigVal);
    2231           6 :             if (nBigVal < 0 || nBigVal > static_cast<GIntBig>(UINT_MAX))
    2232             :             {
    2233           0 :                 CPLError(CE_Warning, CPLE_AppDefined,
    2234             :                          "Value " CPL_FRMT_GIB " of index %d in PDS should be "
    2235             :                          "in [0,%d] range",
    2236             :                          nBigVal, i, INT_MAX);
    2237             :             }
    2238           6 :             WriteUInt32(fp, static_cast<GUInt32>(nBigVal));
    2239             :         }
    2240          40 :         else if (nEltSize == -1)
    2241             :         {
    2242          16 :             if (nVal < -128 || nVal > 127)
    2243             :             {
    2244           2 :                 CPLError(CE_Warning, CPLE_AppDefined,
    2245             :                          "Value %d of index %d in PDS should be in [-128,127] "
    2246             :                          "range",
    2247             :                          nVal, i);
    2248             :             }
    2249          16 :             WriteSByte(fp, nVal);
    2250             :         }
    2251          24 :         else if (nEltSize == -2)
    2252             :         {
    2253           1 :             if (nVal < -32768 || nVal > 32767)
    2254             :             {
    2255           1 :                 CPLError(CE_Warning, CPLE_AppDefined,
    2256             :                          "Value %d of index %d in PDS should be in "
    2257             :                          "[-32768,32767] range",
    2258             :                          nVal, i);
    2259             :             }
    2260           1 :             WriteInt16(fp, nVal);
    2261             :         }
    2262          23 :         else if (nEltSize == -4)
    2263             :         {
    2264          23 :             GIntBig nBigVal = CPLAtoGIntBig(papszTokens[i]);
    2265          23 :             if (nBigVal < INT_MIN || nBigVal > INT_MAX)
    2266             :             {
    2267           1 :                 CPLError(CE_Warning, CPLE_AppDefined,
    2268             :                          "Value " CPL_FRMT_GIB " of index %d in PDS should be "
    2269             :                          "in [%d,%d] range",
    2270             :                          nBigVal, i, INT_MIN, INT_MAX);
    2271             :             }
    2272          23 :             WriteInt32(fp, atoi(papszTokens[i]));
    2273             :         }
    2274             :         else
    2275             :         {
    2276           0 :             CPLAssert(false);
    2277             :         }
    2278             :     }
    2279          12 : }
    2280             : 
    2281             : /************************************************************************/
    2282             : /*                         ComputeValOffset()                           */
    2283             : /************************************************************************/
    2284             : 
    2285          38 : static float ComputeValOffset(int nTokens, char **papszTokens,
    2286             :                               const char *pszInputUnit)
    2287             : {
    2288          38 :     float fValOffset = 0.0f;
    2289             : 
    2290             :     // Parameter category 0 = Temperature
    2291          38 :     if (nTokens >= 2 && atoi(papszTokens[0]) == 0)
    2292             :     {
    2293             :         // Cf
    2294             :         // https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table4-2-0-0.shtml
    2295             :         // PARAMETERS FOR DISCIPLINE 0 CATEGORY 0
    2296          11 :         int nParamNumber = atoi(papszTokens[1]);
    2297          11 :         if ((nParamNumber >= 0 && nParamNumber <= 18 && nParamNumber != 8 &&
    2298          11 :              nParamNumber != 10 && nParamNumber != 11 && nParamNumber != 16) ||
    2299           1 :             nParamNumber == 21 || nParamNumber == 27)
    2300             :         {
    2301          10 :             if (pszInputUnit == nullptr || EQUAL(pszInputUnit, "C") ||
    2302           6 :                 EQUAL(pszInputUnit, "[C]"))
    2303             :             {
    2304           8 :                 fValOffset = 273.15f;
    2305           8 :                 CPLDebug("GRIB",
    2306             :                          "Applying a %f offset to convert from "
    2307             :                          "Celsius to Kelvin",
    2308             :                          fValOffset);
    2309             :             }
    2310             :         }
    2311             :     }
    2312             : 
    2313          38 :     return fValOffset;
    2314             : }
    2315             : 
    2316             : /************************************************************************/
    2317             : /*                           WriteSection4()                            */
    2318             : /************************************************************************/
    2319             : 
    2320         162 : static bool WriteSection4(VSILFILE *fp, GDALDataset *poSrcDS, int nBand,
    2321             :                           char **papszOptions, float &fValOffset)
    2322             : {
    2323             :     // Section 4: Product Definition Section
    2324         162 :     vsi_l_offset nStartSection4 = VSIFTellL(fp);
    2325         162 :     WriteUInt32(fp, GRIB2MISSING_u4);  // section size
    2326         162 :     WriteByte(fp, 4);                  // section number
    2327         162 :     WriteUInt16(fp, 0);  // Number of coordinate values after template
    2328             : 
    2329             :     // 0 = Analysis or forecast at a horizontal level or in a horizontal
    2330             :     // layer at a point in time
    2331             :     int nPDTN =
    2332         162 :         atoi(GetBandOption(papszOptions, poSrcDS, nBand, "PDS_PDTN", "0"));
    2333         162 :     const char *pszPDSTemplateNumbers = GetBandOption(
    2334             :         papszOptions, nullptr, nBand, "PDS_TEMPLATE_NUMBERS", nullptr);
    2335         162 :     const char *pszPDSTemplateAssembledValues = GetBandOption(
    2336             :         papszOptions, nullptr, nBand, "PDS_TEMPLATE_ASSEMBLED_VALUES", nullptr);
    2337         162 :     if (pszPDSTemplateNumbers == nullptr &&
    2338             :         pszPDSTemplateAssembledValues == nullptr)
    2339             :     {
    2340         140 :         pszPDSTemplateNumbers = GetBandOption(papszOptions, poSrcDS, nBand,
    2341             :                                               "PDS_TEMPLATE_NUMBERS", nullptr);
    2342             :     }
    2343         324 :     CPLString osInputUnit;
    2344             :     const char *pszInputUnit =
    2345         162 :         GetBandOption(papszOptions, nullptr, nBand, "INPUT_UNIT", nullptr);
    2346         162 :     if (pszInputUnit == nullptr)
    2347             :     {
    2348             :         const char *pszGribUnit =
    2349         160 :             poSrcDS->GetRasterBand(nBand)->GetMetadataItem("GRIB_UNIT");
    2350         160 :         if (pszGribUnit != nullptr)
    2351             :         {
    2352          21 :             osInputUnit = pszGribUnit;
    2353          21 :             pszInputUnit = osInputUnit.c_str();
    2354             :         }
    2355             :     }
    2356         162 :     WriteUInt16(fp, nPDTN);  // PDTN
    2357         162 :     if (nPDTN == 0 && pszPDSTemplateNumbers == nullptr &&
    2358             :         pszPDSTemplateAssembledValues == nullptr)
    2359             :     {
    2360             :         // See http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_temp4-0.shtml
    2361         120 :         WriteByte(fp, GRIB2MISSING_u1);  // Parameter category = Missing
    2362         120 :         WriteByte(fp, GRIB2MISSING_u1);  // Parameter number = Missing
    2363         120 :         WriteByte(fp, GRIB2MISSING_u1);  // Type of generating process = Missing
    2364         120 :         WriteByte(fp, 0);  // Background generating process identifier
    2365             :         // Analysis or forecast generating process identified
    2366         120 :         WriteByte(fp, GRIB2MISSING_u1);
    2367         120 :         WriteUInt16(fp, 0);  // Hours
    2368         120 :         WriteByte(fp, 0);    // Minutes
    2369         120 :         WriteByte(fp, 0);    // Indicator of unit of time range: 0=Minute
    2370         120 :         WriteUInt32(fp, 0);  // Forecast time in units
    2371         120 :         WriteByte(fp, 0);    // Type of first fixed surface
    2372         120 :         WriteByte(fp, 0);    // Scale factor of first fixed surface
    2373         120 :         WriteUInt32(fp, 0);  // Type of second fixed surface
    2374         120 :         WriteByte(fp, GRIB2MISSING_u1);  // Type of second fixed surface
    2375         120 :         WriteByte(fp, GRIB2MISSING_u1);  // Scale factor of second fixed surface
    2376             :         // Scaled value of second fixed surface
    2377         120 :         WriteUInt32(fp, GRIB2MISSING_u4);
    2378             :     }
    2379          42 :     else if (pszPDSTemplateNumbers == nullptr &&
    2380             :              pszPDSTemplateAssembledValues == nullptr)
    2381             :     {
    2382           1 :         CPLError(CE_Failure, CPLE_NotSupported,
    2383             :                  "PDS_PDTN != 0 specified but both PDS_TEMPLATE_NUMBERS and "
    2384             :                  "PDS_TEMPLATE_ASSEMBLED_VALUES missing");
    2385           1 :         return false;
    2386             :     }
    2387          41 :     else if (pszPDSTemplateNumbers != nullptr &&
    2388             :              pszPDSTemplateAssembledValues != nullptr)
    2389             :     {
    2390           1 :         CPLError(CE_Failure, CPLE_NotSupported,
    2391             :                  "PDS_TEMPLATE_NUMBERS and "
    2392             :                  "PDS_TEMPLATE_ASSEMBLED_VALUES are exclusive");
    2393           1 :         return false;
    2394             :     }
    2395          40 :     else if (pszPDSTemplateNumbers != nullptr)
    2396             :     {
    2397          30 :         char **papszTokens = CSLTokenizeString2(pszPDSTemplateNumbers, " ", 0);
    2398          30 :         const int nTokens = CSLCount(papszTokens);
    2399             : 
    2400          30 :         fValOffset = ComputeValOffset(nTokens, papszTokens, pszInputUnit);
    2401             : 
    2402         952 :         for (int i = 0; papszTokens[i] != nullptr; i++)
    2403             :         {
    2404         922 :             int nVal = atoi(papszTokens[i]);
    2405         922 :             if (nVal < 0 || nVal > 255)
    2406             :             {
    2407           2 :                 CPLError(CE_Warning, CPLE_AppDefined,
    2408             :                          "Value %d of index %d in PDS should be in [0,255] "
    2409             :                          "range",
    2410             :                          nVal, i);
    2411             :             }
    2412         922 :             WriteByte(fp, nVal);
    2413             :         }
    2414          30 :         CSLDestroy(papszTokens);
    2415             : 
    2416             :         // Read back section
    2417          30 :         PatchSectionSize(fp, nStartSection4);
    2418             : 
    2419          30 :         vsi_l_offset nCurOffset = VSIFTellL(fp);
    2420          30 :         VSIFSeekL(fp, nStartSection4, SEEK_SET);
    2421          30 :         size_t nSizeSect4 = static_cast<size_t>(nCurOffset - nStartSection4);
    2422          30 :         GByte *pabySect4 = static_cast<GByte *>(CPLMalloc(nSizeSect4));
    2423          30 :         VSIFReadL(pabySect4, 1, nSizeSect4, fp);
    2424          30 :         VSIFSeekL(fp, nCurOffset, SEEK_SET);
    2425             : 
    2426             :         // Check consistency with template definition
    2427          30 :         g2int iofst = 0;
    2428          30 :         g2int pdsnum = 0;
    2429          30 :         g2int *pdstempl = nullptr;
    2430          30 :         g2int mappdslen = 0;
    2431          30 :         g2float *coordlist = nullptr;
    2432          30 :         g2int numcoord = 0;
    2433             :         int ret =
    2434          30 :             g2_unpack4(pabySect4, static_cast<g2int>(nSizeSect4), &iofst,
    2435             :                        &pdsnum, &pdstempl, &mappdslen, &coordlist, &numcoord);
    2436          30 :         CPLFree(pabySect4);
    2437          30 :         if (ret == 0)
    2438             :         {
    2439          29 :             gtemplate *mappds = extpdstemplate(pdsnum, pdstempl);
    2440          29 :             free(pdstempl);
    2441          29 :             free(coordlist);
    2442          29 :             if (mappds)
    2443             :             {
    2444          29 :                 int nTemplateByteCount = 0;
    2445         564 :                 for (int i = 0; i < mappds->maplen; i++)
    2446         535 :                     nTemplateByteCount += abs(mappds->map[i]);
    2447          39 :                 for (int i = 0; i < mappds->extlen; i++)
    2448          10 :                     nTemplateByteCount += abs(mappds->ext[i]);
    2449          29 :                 if (nTokens < nTemplateByteCount)
    2450             :                 {
    2451           1 :                     CPLError(CE_Failure, CPLE_AppDefined,
    2452             :                              "PDS_PDTN = %d (with provided elements) requires "
    2453             :                              "%d bytes in PDS_TEMPLATE_NUMBERS. "
    2454             :                              "Only %d provided",
    2455             :                              nPDTN, nTemplateByteCount, nTokens);
    2456           1 :                     free(mappds->ext);
    2457           1 :                     free(mappds);
    2458           1 :                     return false;
    2459             :                 }
    2460          28 :                 else if (nTokens > nTemplateByteCount)
    2461             :                 {
    2462           1 :                     CPLError(CE_Warning, CPLE_AppDefined,
    2463             :                              "PDS_PDTN = %d (with provided elements) requires "
    2464             :                              "%d bytes in PDS_TEMPLATE_NUMBERS. "
    2465             :                              "But %d provided. Extra bytes will be ignored "
    2466             :                              "by readers",
    2467             :                              nPDTN, nTemplateByteCount, nTokens);
    2468             :                 }
    2469             : 
    2470          28 :                 free(mappds->ext);
    2471          28 :                 free(mappds);
    2472             :             }
    2473             :         }
    2474             :         else
    2475             :         {
    2476           1 :             free(pdstempl);
    2477           1 :             free(coordlist);
    2478           1 :             CPLError(CE_Warning, CPLE_AppDefined,
    2479             :                      "PDS_PDTN = %d is unknown. Product will not be "
    2480             :                      "correctly read by this driver (but potentially valid "
    2481             :                      "for other readers)",
    2482             :                      nPDTN);
    2483             :         }
    2484             :     }
    2485             :     else
    2486             :     {
    2487          10 :         gtemplate *mappds = getpdstemplate(nPDTN);
    2488          10 :         if (mappds == nullptr)
    2489             :         {
    2490           1 :             CPLError(CE_Failure, CPLE_NotSupported,
    2491             :                      "PDS_PDTN = %d is unknown, so it is not possible to use "
    2492             :                      "PDS_TEMPLATE_ASSEMBLED_VALUES. Use PDS_TEMPLATE_NUMBERS "
    2493             :                      "instead",
    2494             :                      nPDTN);
    2495           3 :             return false;
    2496             :         }
    2497             :         char **papszTokens =
    2498           9 :             CSLTokenizeString2(pszPDSTemplateAssembledValues, " ", 0);
    2499           9 :         const int nTokens = CSLCount(papszTokens);
    2500           9 :         if (nTokens < mappds->maplen)
    2501             :         {
    2502           1 :             CPLError(CE_Failure, CPLE_AppDefined,
    2503             :                      "PDS_PDTN = %d requires at least %d elements in "
    2504             :                      "PDS_TEMPLATE_ASSEMBLED_VALUES. Only %d provided",
    2505             :                      nPDTN, mappds->maplen, nTokens);
    2506           1 :             free(mappds);
    2507           1 :             CSLDestroy(papszTokens);
    2508           1 :             return false;
    2509             :         }
    2510             : 
    2511           8 :         fValOffset = ComputeValOffset(nTokens, papszTokens, pszInputUnit);
    2512             : 
    2513           8 :         std::vector<int> anVals;
    2514           8 :         WriteAssembledPDS(fp, mappds, false, papszTokens, anVals);
    2515             : 
    2516           8 :         if (mappds->needext && !anVals.empty())
    2517             :         {
    2518           5 :             free(mappds);
    2519           5 :             mappds = extpdstemplate(nPDTN, &anVals[0]);
    2520           5 :             if (mappds == nullptr)
    2521             :             {
    2522           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
    2523             :                          "Could not get extended template definition");
    2524           0 :                 CSLDestroy(papszTokens);
    2525           0 :                 return false;
    2526             :             }
    2527           5 :             if (nTokens < mappds->maplen + mappds->extlen)
    2528             :             {
    2529           1 :                 CPLError(CE_Failure, CPLE_AppDefined,
    2530             :                          "PDS_PDTN = %d (with provided elements) requires "
    2531             :                          "%d elements in PDS_TEMPLATE_ASSEMBLED_VALUES. "
    2532             :                          "Only %d provided",
    2533           1 :                          nPDTN, mappds->maplen + mappds->extlen, nTokens);
    2534           1 :                 free(mappds->ext);
    2535           1 :                 free(mappds);
    2536           1 :                 CSLDestroy(papszTokens);
    2537           1 :                 return false;
    2538             :             }
    2539           4 :             else if (nTokens > mappds->maplen + mappds->extlen)
    2540             :             {
    2541           1 :                 CPLError(CE_Warning, CPLE_AppDefined,
    2542             :                          "PDS_PDTN = %d (with provided elements) requires"
    2543             :                          "%d elements in PDS_TEMPLATE_ASSEMBLED_VALUES. "
    2544             :                          "But %d provided. Extra elements will be ignored",
    2545           1 :                          nPDTN, mappds->maplen + mappds->extlen, nTokens);
    2546             :             }
    2547             : 
    2548           4 :             WriteAssembledPDS(fp, mappds, true, papszTokens, anVals);
    2549             :         }
    2550             : 
    2551           7 :         free(mappds->ext);
    2552           7 :         free(mappds);
    2553           7 :         CSLDestroy(papszTokens);
    2554             :     }
    2555         156 :     PatchSectionSize(fp, nStartSection4);
    2556         156 :     return true;
    2557             : }
    2558             : 
    2559             : /************************************************************************/
    2560             : /*                             CreateCopy()                             */
    2561             : /************************************************************************/
    2562             : 
    2563         159 : GDALDataset *GRIBDataset::CreateCopy(const char *pszFilename,
    2564             :                                      GDALDataset *poSrcDS, int /* bStrict */,
    2565             :                                      char **papszOptions,
    2566             :                                      GDALProgressFunc pfnProgress,
    2567             :                                      void *pProgressData)
    2568             : 
    2569             : {
    2570         318 :     if (poSrcDS->GetRasterYSize() == 0 ||
    2571         159 :         poSrcDS->GetRasterXSize() > INT_MAX / poSrcDS->GetRasterXSize())
    2572             :     {
    2573           1 :         CPLError(CE_Failure, CPLE_NotSupported,
    2574             :                  "Cannot create GRIB2 rasters with more than 2 billion pixels");
    2575           1 :         return nullptr;
    2576             :     }
    2577             : 
    2578         158 :     GDALGeoTransform gt;
    2579         158 :     if (poSrcDS->GetGeoTransform(gt) != CE_None)
    2580             :     {
    2581           1 :         CPLError(CE_Failure, CPLE_NotSupported,
    2582             :                  "Source dataset must have a geotransform");
    2583           1 :         return nullptr;
    2584             :     }
    2585         157 :     if (gt[2] != 0.0 || gt[4] != 0.0)
    2586             :     {
    2587           1 :         CPLError(CE_Failure, CPLE_NotSupported,
    2588             :                  "Geotransform with rotation terms not supported");
    2589           1 :         return nullptr;
    2590             :     }
    2591             : 
    2592         312 :     OGRSpatialReference oSRS;
    2593         156 :     oSRS.importFromWkt(poSrcDS->GetProjectionRef());
    2594         156 :     if (oSRS.IsProjected())
    2595             :     {
    2596          90 :         const char *pszProjection = oSRS.GetAttrValue("PROJECTION");
    2597          90 :         if (pszProjection == nullptr ||
    2598          90 :             !(EQUAL(pszProjection, SRS_PT_TRANSVERSE_MERCATOR) ||
    2599          11 :               EQUAL(pszProjection, SRS_PT_MERCATOR_1SP) ||
    2600           9 :               EQUAL(pszProjection, SRS_PT_MERCATOR_2SP) ||
    2601           6 :               EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC) ||
    2602           4 :               EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP) ||
    2603           3 :               EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP) ||
    2604           2 :               EQUAL(pszProjection, SRS_PT_ALBERS_CONIC_EQUAL_AREA) ||
    2605           1 :               EQUAL(pszProjection, SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA)))
    2606             :         {
    2607           0 :             CPLError(CE_Failure, CPLE_NotSupported,
    2608             :                      "Unsupported projection: %s",
    2609             :                      pszProjection ? pszProjection : "");
    2610           0 :             return nullptr;
    2611             :         }
    2612             :     }
    2613          66 :     else if (!oSRS.IsGeographic())
    2614             :     {
    2615           1 :         CPLError(CE_Failure, CPLE_NotSupported,
    2616             :                  "Unsupported or missing spatial reference system");
    2617           1 :         return nullptr;
    2618             :     }
    2619             : 
    2620         155 :     const bool bAppendSubdataset = CPLTestBool(
    2621             :         CSLFetchNameValueDef(papszOptions, "APPEND_SUBDATASET", "NO"));
    2622         155 :     VSILFILE *fp = VSIFOpenL(pszFilename, bAppendSubdataset ? "rb+" : "wb+");
    2623         155 :     if (fp == nullptr)
    2624             :     {
    2625           4 :         CPLError(CE_Failure, CPLE_FileIO, "Cannot create %s", pszFilename);
    2626           4 :         return nullptr;
    2627             :     }
    2628         151 :     VSIFSeekL(fp, 0, SEEK_END);
    2629             : 
    2630         151 :     vsi_l_offset nStartOffset = 0;
    2631         151 :     vsi_l_offset nTotalSizeOffset = 0;
    2632         151 :     int nSplitAndSwapColumn = 0;
    2633             :     // Note: WRITE_SUBGRIDS=YES should not be used blindly currently, as it
    2634             :     // does not check that the content of the DISCIPLINE and IDS are the same.
    2635             :     // A smarter behavior would be to break into separate messages if needed
    2636             :     const bool bWriteSubGrids =
    2637         151 :         CPLTestBool(CSLFetchNameValueDef(papszOptions, "WRITE_SUBGRIDS", "NO"));
    2638         296 :     for (int nBand = 1; nBand <= poSrcDS->GetRasterCount(); nBand++)
    2639             :     {
    2640         162 :         if (nBand == 1 || !bWriteSubGrids)
    2641             :         {
    2642             :             // Section 0: Indicator section
    2643         162 :             nStartOffset = VSIFTellL(fp);
    2644         162 :             VSIFWriteL("GRIB", 4, 1, fp);
    2645         162 :             WriteByte(fp, 0);  // reserved
    2646         162 :             WriteByte(fp, 0);  // reserved
    2647             :             int nDiscipline =
    2648         162 :                 atoi(GetBandOption(papszOptions, poSrcDS, nBand, "DISCIPLINE",
    2649             :                                    "0"));  // 0 = Meteorological
    2650         162 :             WriteByte(fp, nDiscipline);    // discipline
    2651         162 :             WriteByte(fp, 2);              // GRIB edition number
    2652         162 :             nTotalSizeOffset = VSIFTellL(fp);
    2653         162 :             WriteUInt32(fp, GRIB2MISSING_u4);  // dummy file size (high 32 bits)
    2654         162 :             WriteUInt32(fp, GRIB2MISSING_u4);  // dummy file size (low 32 bits)
    2655             : 
    2656             :             // Section 1: Identification Section
    2657         162 :             WriteSection1(fp, poSrcDS, nBand, papszOptions);
    2658             : 
    2659             :             // Section 2: Local use section
    2660         162 :             WriteUInt32(fp, 5);  // section size
    2661         162 :             WriteByte(fp, 2);    // section number
    2662             : 
    2663             :             // Section 3: Grid Definition Section
    2664         162 :             GRIB2Section3Writer oSection3(fp, poSrcDS);
    2665         162 :             if (!oSection3.Write())
    2666             :             {
    2667           0 :                 VSIFCloseL(fp);
    2668           0 :                 return nullptr;
    2669             :             }
    2670         162 :             nSplitAndSwapColumn = oSection3.SplitAndSwap();
    2671             :         }
    2672             : 
    2673             :         // Section 4: Product Definition Section
    2674         162 :         float fValOffset = 0.0f;
    2675         162 :         if (!WriteSection4(fp, poSrcDS, nBand, papszOptions, fValOffset))
    2676             :         {
    2677           6 :             VSIFCloseL(fp);
    2678           6 :             return nullptr;
    2679             :         }
    2680             : 
    2681             :         // Section 5, 6 and 7
    2682         156 :         if (!GRIB2Section567Writer(fp, poSrcDS, nBand, nSplitAndSwapColumn)
    2683         156 :                  .Write(fValOffset, papszOptions, pfnProgress, pProgressData))
    2684             :         {
    2685          11 :             VSIFCloseL(fp);
    2686          11 :             return nullptr;
    2687             :         }
    2688             : 
    2689         145 :         if (nBand == poSrcDS->GetRasterCount() || !bWriteSubGrids)
    2690             :         {
    2691             :             // Section 8: End section
    2692         145 :             VSIFWriteL("7777", 4, 1, fp);
    2693             : 
    2694             :             // Patch total message size at end of section 0
    2695         145 :             vsi_l_offset nCurOffset = VSIFTellL(fp);
    2696         145 :             if (nCurOffset - nStartOffset > INT_MAX)
    2697             :             {
    2698           0 :                 CPLError(CE_Failure, CPLE_NotSupported,
    2699             :                          "GRIB message larger than 2 GB");
    2700           0 :                 VSIFCloseL(fp);
    2701           0 :                 return nullptr;
    2702             :             }
    2703         145 :             GUInt32 nTotalSize =
    2704         145 :                 static_cast<GUInt32>(nCurOffset - nStartOffset);
    2705         145 :             VSIFSeekL(fp, nTotalSizeOffset, SEEK_SET);
    2706         145 :             WriteUInt32(fp, 0);           // file size (high 32 bits)
    2707         145 :             WriteUInt32(fp, nTotalSize);  // file size (low 32 bits)
    2708             : 
    2709         145 :             VSIFSeekL(fp, nCurOffset, SEEK_SET);
    2710             :         }
    2711             : 
    2712         290 :         if (pfnProgress &&
    2713         145 :             !pfnProgress(static_cast<double>(nBand) / poSrcDS->GetRasterCount(),
    2714             :                          nullptr, pProgressData))
    2715             :         {
    2716           0 :             VSIFCloseL(fp);
    2717           0 :             return nullptr;
    2718             :         }
    2719             :     }
    2720             : 
    2721         134 :     VSIFCloseL(fp);
    2722             : 
    2723         268 :     GDALOpenInfo oOpenInfo(pszFilename, GA_ReadOnly);
    2724         134 :     return Open(&oOpenInfo);
    2725             : }

Generated by: LCOV version 1.14