LCOV - code coverage report
Current view: top level - frmts/aaigrid - aaigriddataset.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 690 852 81.0 %
Date: 2024-05-04 12:52:34 Functions: 32 33 97.0 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  Implements Arc/Info ASCII Grid Format.
       5             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2001, Frank Warmerdam (warmerdam@pobox.com)
       9             :  * Copyright (c) 2007-2012, Even Rouault <even dot rouault at spatialys.com>
      10             :  * Copyright (c) 2014, Kyle Shannon <kyle at pobox dot com>
      11             :  *
      12             :  * Permission is hereby granted, free of charge, to any person obtaining a
      13             :  * copy of this software and associated documentation files (the "Software"),
      14             :  * to deal in the Software without restriction, including without limitation
      15             :  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
      16             :  * and/or sell copies of the Software, and to permit persons to whom the
      17             :  * Software is furnished to do so, subject to the following conditions:
      18             :  *
      19             :  * The above copyright notice and this permission notice shall be included
      20             :  * in all copies or substantial portions of the Software.
      21             :  *
      22             :  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
      23             :  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      24             :  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
      25             :  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      26             :  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
      27             :  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
      28             :  * DEALINGS IN THE SOFTWARE.
      29             :  ****************************************************************************/
      30             : 
      31             : // We need cpl_port as first include to avoid VSIStatBufL being not
      32             : // defined on i586-mingw32msvc.
      33             : #include "cpl_port.h"
      34             : #include "aaigriddataset.h"
      35             : #include "gdal_frmts.h"
      36             : 
      37             : #include <cassert>
      38             : #include <cctype>
      39             : #include <climits>
      40             : #include <cmath>
      41             : #include <cstddef>
      42             : #include <cstdio>
      43             : #include <cstdlib>
      44             : #include <cstring>
      45             : #if HAVE_FCNTL_H
      46             : #include <fcntl.h>
      47             : #endif
      48             : 
      49             : #include <algorithm>
      50             : #include <limits>
      51             : #include <string>
      52             : 
      53             : #include "cpl_conv.h"
      54             : #include "cpl_error.h"
      55             : #include "cpl_progress.h"
      56             : #include "cpl_string.h"
      57             : #include "cpl_vsi.h"
      58             : #include "gdal.h"
      59             : #include "gdal_pam.h"
      60             : #include "gdal_priv.h"
      61             : #include "ogr_core.h"
      62             : #include "ogr_spatialref.h"
      63             : 
      64             : namespace
      65             : {
      66             : 
      67        2629 : float DoubleToFloatClamp(double dfValue)
      68             : {
      69        2629 :     if (dfValue <= std::numeric_limits<float>::lowest())
      70           0 :         return std::numeric_limits<float>::lowest();
      71        2629 :     if (dfValue >= std::numeric_limits<float>::max())
      72           0 :         return std::numeric_limits<float>::max();
      73        2629 :     return static_cast<float>(dfValue);
      74             : }
      75             : 
      76             : // Cast to float and back for make sure the NoData value matches
      77             : // that expressed by a float value.  Clamps to the range of a float
      78             : // if the value is too large.  Preserves +/-inf and NaN.
      79             : // TODO(schwehr): This should probably be moved to port as it is likely
      80             : // to be needed for other formats.
      81          13 : double MapNoDataToFloat(double dfNoDataValue)
      82             : {
      83          13 :     if (CPLIsInf(dfNoDataValue) || CPLIsNan(dfNoDataValue))
      84           0 :         return dfNoDataValue;
      85             : 
      86          13 :     if (dfNoDataValue >= std::numeric_limits<float>::max())
      87           0 :         return std::numeric_limits<float>::max();
      88             : 
      89          13 :     if (dfNoDataValue <= -std::numeric_limits<float>::max())
      90           0 :         return -std::numeric_limits<float>::max();
      91             : 
      92          13 :     return static_cast<double>(static_cast<float>(dfNoDataValue));
      93             : }
      94             : 
      95             : }  // namespace
      96             : 
      97             : static CPLString OSR_GDS(char **papszNV, const char *pszField,
      98             :                          const char *pszDefaultValue);
      99             : 
     100             : /************************************************************************/
     101             : /*                           AAIGRasterBand()                           */
     102             : /************************************************************************/
     103             : 
     104         163 : AAIGRasterBand::AAIGRasterBand(AAIGDataset *poDSIn, int nDataStart)
     105         163 :     : panLineOffset(nullptr)
     106             : {
     107         163 :     poDS = poDSIn;
     108             : 
     109         163 :     nBand = 1;
     110         163 :     eDataType = poDSIn->eDataType;
     111             : 
     112         163 :     nBlockXSize = poDSIn->nRasterXSize;
     113         163 :     nBlockYSize = 1;
     114             : 
     115         163 :     panLineOffset = static_cast<GUIntBig *>(
     116         163 :         VSI_CALLOC_VERBOSE(poDSIn->nRasterYSize, sizeof(GUIntBig)));
     117         163 :     if (panLineOffset == nullptr)
     118             :     {
     119           0 :         return;
     120             :     }
     121         163 :     panLineOffset[0] = nDataStart;
     122             : }
     123             : 
     124             : /************************************************************************/
     125             : /*                          ~AAIGRasterBand()                           */
     126             : /************************************************************************/
     127             : 
     128         326 : AAIGRasterBand::~AAIGRasterBand()
     129             : {
     130         163 :     CPLFree(panLineOffset);
     131         326 : }
     132             : 
     133             : /************************************************************************/
     134             : /*                             IReadBlock()                             */
     135             : /************************************************************************/
     136             : 
     137         997 : CPLErr AAIGRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage)
     138             : 
     139             : {
     140         997 :     AAIGDataset *poODS = static_cast<AAIGDataset *>(poDS);
     141             : 
     142         997 :     if (nBlockYOff < 0 || nBlockYOff > poODS->nRasterYSize - 1 ||
     143         997 :         nBlockXOff != 0 || panLineOffset == nullptr || poODS->fp == nullptr)
     144           0 :         return CE_Failure;
     145             : 
     146         997 :     if (panLineOffset[nBlockYOff] == 0)
     147             :     {
     148          22 :         for (int iPrevLine = 1; iPrevLine <= nBlockYOff; iPrevLine++)
     149          18 :             if (panLineOffset[iPrevLine] == 0)
     150          18 :                 IReadBlock(nBlockXOff, iPrevLine - 1, nullptr);
     151             :     }
     152             : 
     153         997 :     if (panLineOffset[nBlockYOff] == 0)
     154           0 :         return CE_Failure;
     155             : 
     156         997 :     if (poODS->Seek(panLineOffset[nBlockYOff]) != 0)
     157             :     {
     158           0 :         ReportError(CE_Failure, CPLE_FileIO,
     159             :                     "Can't seek to offset %lu in input file to read data.",
     160           0 :                     static_cast<long unsigned int>(panLineOffset[nBlockYOff]));
     161           0 :         return CE_Failure;
     162             :     }
     163             : 
     164       33814 :     for (int iPixel = 0; iPixel < poODS->nRasterXSize;)
     165             :     {
     166             :         // Suck up any pre-white space.
     167       32817 :         char chNext = '\0';
     168       11022 :         do
     169             :         {
     170       43839 :             chNext = poODS->Getc();
     171       43839 :         } while (isspace(static_cast<unsigned char>(chNext)));
     172             : 
     173       32817 :         char szToken[500] = {'\0'};
     174       32817 :         int iTokenChar = 0;
     175      133097 :         while (chNext != '\0' && !isspace((unsigned char)chNext))
     176             :         {
     177      100280 :             if (iTokenChar == sizeof(szToken) - 2)
     178             :             {
     179           0 :                 ReportError(CE_Failure, CPLE_FileIO,
     180             :                             "Token too long at scanline %d.", nBlockYOff);
     181           0 :                 return CE_Failure;
     182             :             }
     183             : 
     184      100280 :             szToken[iTokenChar++] = chNext;
     185      100280 :             chNext = poODS->Getc();
     186             :         }
     187             : 
     188       32817 :         if (chNext == '\0' && (iPixel != poODS->nRasterXSize - 1 ||
     189          26 :                                nBlockYOff != poODS->nRasterYSize - 1))
     190             :         {
     191           0 :             ReportError(CE_Failure, CPLE_FileIO,
     192             :                         "File short, can't read line %d.", nBlockYOff);
     193           0 :             return CE_Failure;
     194             :         }
     195             : 
     196       32817 :         szToken[iTokenChar] = '\0';
     197             : 
     198       32817 :         if (pImage != nullptr)
     199             :         {
     200             :             // "null" seems to be specific of D12 software
     201             :             // See https://github.com/OSGeo/gdal/issues/5095
     202       32619 :             if (eDataType == GDT_Float64)
     203             :             {
     204         104 :                 if (strcmp(szToken, "null") == 0)
     205           2 :                     reinterpret_cast<double *>(pImage)[iPixel] =
     206           2 :                         -std::numeric_limits<double>::max();
     207             :                 else
     208         204 :                     reinterpret_cast<double *>(pImage)[iPixel] =
     209         102 :                         CPLAtofM(szToken);
     210             :             }
     211       32515 :             else if (eDataType == GDT_Float32)
     212             :             {
     213        2631 :                 if (strcmp(szToken, "null") == 0)
     214           2 :                     reinterpret_cast<float *>(pImage)[iPixel] =
     215           2 :                         -std::numeric_limits<float>::max();
     216             :                 else
     217        2629 :                     reinterpret_cast<float *>(pImage)[iPixel] =
     218        2629 :                         DoubleToFloatClamp(CPLAtofM(szToken));
     219             :             }
     220             :             else
     221       29884 :                 reinterpret_cast<GInt32 *>(pImage)[iPixel] =
     222       29884 :                     static_cast<GInt32>(atoi(szToken));
     223             :         }
     224             : 
     225       32817 :         iPixel++;
     226             :     }
     227             : 
     228         997 :     if (nBlockYOff < poODS->nRasterYSize - 1)
     229         905 :         panLineOffset[nBlockYOff + 1] = poODS->Tell();
     230             : 
     231         997 :     return CE_None;
     232             : }
     233             : 
     234             : /************************************************************************/
     235             : /*                           GetNoDataValue()                           */
     236             : /************************************************************************/
     237             : 
     238         169 : double AAIGRasterBand::GetNoDataValue(int *pbSuccess)
     239             : 
     240             : {
     241         169 :     AAIGDataset *poODS = static_cast<AAIGDataset *>(poDS);
     242             : 
     243         169 :     if (pbSuccess)
     244         151 :         *pbSuccess = poODS->bNoDataSet;
     245             : 
     246         169 :     return poODS->dfNoDataValue;
     247             : }
     248             : 
     249             : /************************************************************************/
     250             : /*                           SetNoDataValue()                           */
     251             : /************************************************************************/
     252             : 
     253           0 : CPLErr AAIGRasterBand::SetNoDataValue(double dfNoData)
     254             : 
     255             : {
     256           0 :     AAIGDataset *poODS = static_cast<AAIGDataset *>(poDS);
     257             : 
     258           0 :     poODS->bNoDataSet = true;
     259           0 :     poODS->dfNoDataValue = dfNoData;
     260             : 
     261           0 :     return CE_None;
     262             : }
     263             : 
     264             : /************************************************************************/
     265             : /* ==================================================================== */
     266             : /*                            AAIGDataset                               */
     267             : /* ==================================================================== */
     268             : /************************************************************************/
     269             : 
     270             : /************************************************************************/
     271             : /*                            AAIGDataset()                            */
     272             : /************************************************************************/
     273             : 
     274         164 : AAIGDataset::AAIGDataset()
     275             :     : fp(nullptr), papszPrj(nullptr), nBufferOffset(0), nOffsetInBuffer(256),
     276         164 :       eDataType(GDT_Int32), bNoDataSet(false), dfNoDataValue(-9999.0)
     277             : {
     278         164 :     m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     279         164 :     adfGeoTransform[0] = 0.0;
     280         164 :     adfGeoTransform[1] = 1.0;
     281         164 :     adfGeoTransform[2] = 0.0;
     282         164 :     adfGeoTransform[3] = 0.0;
     283         164 :     adfGeoTransform[4] = 0.0;
     284         164 :     adfGeoTransform[5] = 1.0;
     285         164 :     memset(achReadBuf, 0, sizeof(achReadBuf));
     286         164 : }
     287             : 
     288             : /************************************************************************/
     289             : /*                           ~AAIGDataset()                            */
     290             : /************************************************************************/
     291             : 
     292         319 : AAIGDataset::~AAIGDataset()
     293             : 
     294             : {
     295         164 :     FlushCache(true);
     296             : 
     297         164 :     if (fp != nullptr)
     298             :     {
     299         163 :         if (VSIFCloseL(fp) != 0)
     300             :         {
     301           0 :             ReportError(CE_Failure, CPLE_FileIO, "I/O error");
     302             :         }
     303             :     }
     304             : 
     305         164 :     CSLDestroy(papszPrj);
     306         319 : }
     307             : 
     308             : /************************************************************************/
     309             : /*                                Tell()                                */
     310             : /************************************************************************/
     311             : 
     312         905 : GUIntBig AAIGDataset::Tell() const
     313             : {
     314         905 :     return nBufferOffset + nOffsetInBuffer;
     315             : }
     316             : 
     317             : /************************************************************************/
     318             : /*                                Seek()                                */
     319             : /************************************************************************/
     320             : 
     321         997 : int AAIGDataset::Seek(GUIntBig nNewOffset)
     322             : 
     323             : {
     324         997 :     nOffsetInBuffer = sizeof(achReadBuf);
     325         997 :     return VSIFSeekL(fp, nNewOffset, SEEK_SET);
     326             : }
     327             : 
     328             : /************************************************************************/
     329             : /*                                Getc()                                */
     330             : /*                                                                      */
     331             : /*      Read a single character from the input file (efficiently we     */
     332             : /*      hope).                                                          */
     333             : /************************************************************************/
     334             : 
     335      144119 : char AAIGDataset::Getc()
     336             : 
     337             : {
     338      144119 :     if (nOffsetInBuffer < static_cast<int>(sizeof(achReadBuf)))
     339      142922 :         return achReadBuf[nOffsetInBuffer++];
     340             : 
     341        1197 :     nBufferOffset = VSIFTellL(fp);
     342             :     const int nRead =
     343        1197 :         static_cast<int>(VSIFReadL(achReadBuf, 1, sizeof(achReadBuf), fp));
     344       56215 :     for (unsigned int i = nRead; i < sizeof(achReadBuf); i++)
     345       55018 :         achReadBuf[i] = '\0';
     346             : 
     347        1197 :     nOffsetInBuffer = 0;
     348             : 
     349        1197 :     return achReadBuf[nOffsetInBuffer++];
     350             : }
     351             : 
     352             : /************************************************************************/
     353             : /*                            GetFileList()                             */
     354             : /************************************************************************/
     355             : 
     356          32 : char **AAIGDataset::GetFileList()
     357             : 
     358             : {
     359          32 :     char **papszFileList = GDALPamDataset::GetFileList();
     360             : 
     361          32 :     if (papszPrj != nullptr)
     362          12 :         papszFileList = CSLAddString(papszFileList, osPrjFilename);
     363             : 
     364          32 :     return papszFileList;
     365             : }
     366             : 
     367             : /************************************************************************/
     368             : /*                            Identify()                                */
     369             : /************************************************************************/
     370             : 
     371       54139 : int AAIGDataset::Identify(GDALOpenInfo *poOpenInfo)
     372             : 
     373             : {
     374             :     // Does this look like an AI grid file?
     375       54139 :     if (poOpenInfo->nHeaderBytes < 40 ||
     376        7326 :         !(STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "ncols") ||
     377        7015 :           STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "nrows") ||
     378        7015 :           STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "xllcorner") ||
     379        7015 :           STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "yllcorner") ||
     380        7015 :           STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "xllcenter") ||
     381        7015 :           STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "yllcenter") ||
     382        7015 :           STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "dx") ||
     383        7015 :           STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "dy") ||
     384        7015 :           STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "cellsize")))
     385       53828 :         return FALSE;
     386             : 
     387         311 :     return TRUE;
     388             : }
     389             : 
     390             : /************************************************************************/
     391             : /*                            Identify()                                */
     392             : /************************************************************************/
     393             : 
     394       53783 : int GRASSASCIIDataset::Identify(GDALOpenInfo *poOpenInfo)
     395             : 
     396             : {
     397             :     // Does this look like a GRASS ASCII grid file?
     398       53783 :     if (poOpenInfo->nHeaderBytes < 40 ||
     399        7018 :         !(STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "north:") ||
     400        7014 :           STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "south:") ||
     401        7014 :           STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "east:") ||
     402        7014 :           STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "west:") ||
     403        7014 :           STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "rows:") ||
     404        7014 :           STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "cols:")))
     405       53778 :         return FALSE;
     406             : 
     407           5 :     return TRUE;
     408             : }
     409             : 
     410             : /************************************************************************/
     411             : /*                            Identify()                                */
     412             : /************************************************************************/
     413             : 
     414       53786 : int ISGDataset::Identify(GDALOpenInfo *poOpenInfo)
     415             : 
     416             : {
     417             :     // Does this look like a ISG grid file?
     418       53786 :     if (poOpenInfo->nHeaderBytes < 40 ||
     419        7021 :         !strstr((const char *)poOpenInfo->pabyHeader, "model name"))
     420             :     {
     421       53772 :         return FALSE;
     422             :     }
     423          15 :     for (int i = 0; i < 2; ++i)
     424             :     {
     425          15 :         if (strstr((const char *)poOpenInfo->pabyHeader, "lat min") !=
     426          15 :                 nullptr &&
     427          15 :             strstr((const char *)poOpenInfo->pabyHeader, "lat max") !=
     428          15 :                 nullptr &&
     429          15 :             strstr((const char *)poOpenInfo->pabyHeader, "lon min") !=
     430          15 :                 nullptr &&
     431          15 :             strstr((const char *)poOpenInfo->pabyHeader, "lon max") !=
     432          15 :                 nullptr &&
     433          15 :             strstr((const char *)poOpenInfo->pabyHeader, "nrows") != nullptr &&
     434          14 :             strstr((const char *)poOpenInfo->pabyHeader, "ncols") != nullptr)
     435             :         {
     436          14 :             return TRUE;
     437             :         }
     438             :         // Some files like https://isgeoid.polimi.it/Geoid/Europe/Slovenia/public/Slovenia_2016_SLO_VRP2016_Koper_hybrQ_20221122.isg
     439             :         // have initial comment lines, so we may need to ingest more bytes
     440           1 :         if (i == 0)
     441             :         {
     442           1 :             if (poOpenInfo->nHeaderBytes >= 8192)
     443           0 :                 break;
     444           1 :             poOpenInfo->TryToIngest(8192);
     445             :         }
     446             :     }
     447             : 
     448           0 :     return TRUE;
     449             : }
     450             : 
     451             : /************************************************************************/
     452             : /*                                Open()                                */
     453             : /************************************************************************/
     454             : 
     455         155 : GDALDataset *AAIGDataset::Open(GDALOpenInfo *poOpenInfo)
     456             : {
     457             : #ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
     458             :     // During fuzzing, do not use Identify to reject crazy content.
     459         155 :     if (!Identify(poOpenInfo))
     460           0 :         return nullptr;
     461             : #endif
     462             : 
     463         155 :     return CommonOpen(poOpenInfo, FORMAT_AAIG);
     464             : }
     465             : 
     466             : /************************************************************************/
     467             : /*                          ParseHeader()                               */
     468             : /************************************************************************/
     469             : 
     470         155 : int AAIGDataset::ParseHeader(const char *pszHeader, const char *pszDataType)
     471             : {
     472         155 :     char **papszTokens = CSLTokenizeString2(pszHeader, " \n\r\t", 0);
     473         155 :     const int nTokens = CSLCount(papszTokens);
     474             : 
     475         155 :     int i = 0;
     476         155 :     if ((i = CSLFindString(papszTokens, "ncols")) < 0 || i + 1 >= nTokens)
     477             :     {
     478           0 :         CSLDestroy(papszTokens);
     479           0 :         return FALSE;
     480             :     }
     481         155 :     nRasterXSize = atoi(papszTokens[i + 1]);
     482         155 :     if ((i = CSLFindString(papszTokens, "nrows")) < 0 || i + 1 >= nTokens)
     483             :     {
     484           0 :         CSLDestroy(papszTokens);
     485           0 :         return FALSE;
     486             :     }
     487         155 :     nRasterYSize = atoi(papszTokens[i + 1]);
     488             : 
     489         155 :     if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize))
     490             :     {
     491           0 :         CSLDestroy(papszTokens);
     492           0 :         return FALSE;
     493             :     }
     494             : 
     495             :     // TODO(schwehr): Would be good to also factor the file size into the max.
     496             :     // TODO(schwehr): Allow the user to disable this check.
     497             :     // The driver allocates a panLineOffset array based on nRasterYSize
     498         155 :     constexpr int kMaxDimSize = 10000000;  // 1e7 cells.
     499         155 :     if (nRasterXSize > kMaxDimSize || nRasterYSize > kMaxDimSize)
     500             :     {
     501           0 :         CSLDestroy(papszTokens);
     502           0 :         return FALSE;
     503             :     }
     504             : 
     505         155 :     double dfCellDX = 0.0;
     506         155 :     double dfCellDY = 0.0;
     507         155 :     if ((i = CSLFindString(papszTokens, "cellsize")) < 0)
     508             :     {
     509             :         int iDX, iDY;
     510           4 :         if ((iDX = CSLFindString(papszTokens, "dx")) < 0 ||
     511           4 :             (iDY = CSLFindString(papszTokens, "dy")) < 0 ||
     512           8 :             iDX + 1 >= nTokens || iDY + 1 >= nTokens)
     513             :         {
     514           0 :             CSLDestroy(papszTokens);
     515           0 :             return FALSE;
     516             :         }
     517             : 
     518           4 :         dfCellDX = CPLAtofM(papszTokens[iDX + 1]);
     519           4 :         dfCellDY = CPLAtofM(papszTokens[iDY + 1]);
     520             :     }
     521             :     else
     522             :     {
     523         151 :         if (i + 1 >= nTokens)
     524             :         {
     525           0 :             CSLDestroy(papszTokens);
     526           0 :             return FALSE;
     527             :         }
     528         151 :         dfCellDY = CPLAtofM(papszTokens[i + 1]);
     529         151 :         dfCellDX = dfCellDY;
     530             :     }
     531             : 
     532         155 :     int j = 0;
     533         155 :     if ((i = CSLFindString(papszTokens, "xllcorner")) >= 0 &&
     534         310 :         (j = CSLFindString(papszTokens, "yllcorner")) >= 0 && i + 1 < nTokens &&
     535         155 :         j + 1 < nTokens)
     536             :     {
     537         155 :         adfGeoTransform[0] = CPLAtofM(papszTokens[i + 1]);
     538             : 
     539             :         // Small hack to compensate from insufficient precision in cellsize
     540             :         // parameter in datasets of
     541             :         // http://ccafs-climate.org/data/A2a_2020s/hccpr_hadcm3
     542         155 :         if ((nRasterXSize % 360) == 0 &&
     543           0 :             fabs(adfGeoTransform[0] - (-180.0)) < 1e-12 &&
     544           0 :             dfCellDX == dfCellDY &&
     545           0 :             fabs(dfCellDX - (360.0 / nRasterXSize)) < 1e-9)
     546             :         {
     547           0 :             dfCellDY = 360.0 / nRasterXSize;
     548           0 :             dfCellDX = dfCellDY;
     549             :         }
     550             : 
     551         155 :         adfGeoTransform[1] = dfCellDX;
     552         155 :         adfGeoTransform[2] = 0.0;
     553         155 :         adfGeoTransform[3] =
     554         155 :             CPLAtofM(papszTokens[j + 1]) + nRasterYSize * dfCellDY;
     555         155 :         adfGeoTransform[4] = 0.0;
     556         155 :         adfGeoTransform[5] = -dfCellDY;
     557             :     }
     558           0 :     else if ((i = CSLFindString(papszTokens, "xllcenter")) >= 0 &&
     559           0 :              (j = CSLFindString(papszTokens, "yllcenter")) >= 0 &&
     560           0 :              i + 1 < nTokens && j + 1 < nTokens)
     561             :     {
     562           0 :         SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
     563             : 
     564           0 :         adfGeoTransform[0] = CPLAtofM(papszTokens[i + 1]) - 0.5 * dfCellDX;
     565           0 :         adfGeoTransform[1] = dfCellDX;
     566           0 :         adfGeoTransform[2] = 0.0;
     567           0 :         adfGeoTransform[3] = CPLAtofM(papszTokens[j + 1]) - 0.5 * dfCellDY +
     568           0 :                              nRasterYSize * dfCellDY;
     569           0 :         adfGeoTransform[4] = 0.0;
     570           0 :         adfGeoTransform[5] = -dfCellDY;
     571             :     }
     572             :     else
     573             :     {
     574           0 :         adfGeoTransform[0] = 0.0;
     575           0 :         adfGeoTransform[1] = dfCellDX;
     576           0 :         adfGeoTransform[2] = 0.0;
     577           0 :         adfGeoTransform[3] = 0.0;
     578           0 :         adfGeoTransform[4] = 0.0;
     579           0 :         adfGeoTransform[5] = -dfCellDY;
     580             :     }
     581             : 
     582         220 :     if ((i = CSLFindString(papszTokens, "NODATA_value")) >= 0 &&
     583          65 :         i + 1 < nTokens)
     584             :     {
     585          65 :         const char *pszNoData = papszTokens[i + 1];
     586             : 
     587          65 :         bNoDataSet = true;
     588          65 :         if (strcmp(pszNoData, "null") == 0)
     589             :         {
     590             :             // "null" seems to be specific of D12 software
     591             :             // See https://github.com/OSGeo/gdal/issues/5095
     592           2 :             if (pszDataType == nullptr || eDataType == GDT_Float32)
     593             :             {
     594           1 :                 dfNoDataValue = -std::numeric_limits<float>::max();
     595           1 :                 eDataType = GDT_Float32;
     596             :             }
     597             :             else
     598             :             {
     599           1 :                 dfNoDataValue = -std::numeric_limits<double>::max();
     600           1 :                 eDataType = GDT_Float64;
     601             :             }
     602             :         }
     603             :         else
     604             :         {
     605          63 :             dfNoDataValue = CPLAtofM(pszNoData);
     606         121 :             if (pszDataType == nullptr &&
     607          58 :                 (strchr(pszNoData, '.') != nullptr ||
     608          50 :                  strchr(pszNoData, ',') != nullptr ||
     609          50 :                  std::numeric_limits<int>::min() > dfNoDataValue ||
     610          50 :                  dfNoDataValue > std::numeric_limits<int>::max()))
     611             :             {
     612           8 :                 eDataType = GDT_Float32;
     613          16 :                 if (!CPLIsInf(dfNoDataValue) &&
     614           8 :                     (fabs(dfNoDataValue) < std::numeric_limits<float>::min() ||
     615           7 :                      fabs(dfNoDataValue) > std::numeric_limits<float>::max()))
     616             :                 {
     617           1 :                     eDataType = GDT_Float64;
     618             :                 }
     619             :             }
     620          63 :             if (eDataType == GDT_Float32)
     621             :             {
     622           7 :                 dfNoDataValue = MapNoDataToFloat(dfNoDataValue);
     623             :             }
     624             :         }
     625             :     }
     626             : 
     627         155 :     CSLDestroy(papszTokens);
     628             : 
     629         155 :     return TRUE;
     630             : }
     631             : 
     632             : /************************************************************************/
     633             : /*                                Open()                                */
     634             : /************************************************************************/
     635             : 
     636           2 : GDALDataset *GRASSASCIIDataset::Open(GDALOpenInfo *poOpenInfo)
     637             : {
     638             : #ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
     639             :     // During fuzzing, do not use Identify to reject crazy content.
     640           2 :     if (!Identify(poOpenInfo))
     641           0 :         return nullptr;
     642             : #endif
     643             : 
     644           2 :     return CommonOpen(poOpenInfo, FORMAT_GRASSASCII);
     645             : }
     646             : 
     647             : /************************************************************************/
     648             : /*                          ParseHeader()                               */
     649             : /************************************************************************/
     650             : 
     651           2 : int GRASSASCIIDataset::ParseHeader(const char *pszHeader,
     652             :                                    const char *pszDataType)
     653             : {
     654           2 :     char **papszTokens = CSLTokenizeString2(pszHeader, " \n\r\t:", 0);
     655           2 :     const int nTokens = CSLCount(papszTokens);
     656           2 :     int i = 0;
     657           2 :     if ((i = CSLFindString(papszTokens, "cols")) < 0 || i + 1 >= nTokens)
     658             :     {
     659           0 :         CSLDestroy(papszTokens);
     660           0 :         return FALSE;
     661             :     }
     662           2 :     nRasterXSize = atoi(papszTokens[i + 1]);
     663           2 :     if ((i = CSLFindString(papszTokens, "rows")) < 0 || i + 1 >= nTokens)
     664             :     {
     665           0 :         CSLDestroy(papszTokens);
     666           0 :         return FALSE;
     667             :     }
     668           2 :     nRasterYSize = atoi(papszTokens[i + 1]);
     669             : 
     670           2 :     if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize))
     671             :     {
     672           0 :         CSLDestroy(papszTokens);
     673           0 :         return FALSE;
     674             :     }
     675             : 
     676             :     // TODO(schwehr): Would be good to also factor the file size into the max.
     677             :     // TODO(schwehr): Allow the user to disable this check.
     678             :     // The driver allocates a panLineOffset array based on nRasterYSize
     679           2 :     constexpr int kMaxDimSize = 10000000;  // 1e7 cells.
     680           2 :     if (nRasterXSize > kMaxDimSize || nRasterYSize > kMaxDimSize)
     681             :     {
     682           0 :         CSLDestroy(papszTokens);
     683           0 :         return FALSE;
     684             :     }
     685             : 
     686           2 :     const int iNorth = CSLFindString(papszTokens, "north");
     687           2 :     const int iSouth = CSLFindString(papszTokens, "south");
     688           2 :     const int iEast = CSLFindString(papszTokens, "east");
     689           2 :     const int iWest = CSLFindString(papszTokens, "west");
     690             : 
     691           4 :     if (iNorth == -1 || iSouth == -1 || iEast == -1 || iWest == -1 ||
     692           2 :         std::max(std::max(iNorth, iSouth), std::max(iEast, iWest)) + 1 >=
     693             :             nTokens)
     694             :     {
     695           0 :         CSLDestroy(papszTokens);
     696           0 :         return FALSE;
     697             :     }
     698             : 
     699           2 :     const double dfNorth = CPLAtofM(papszTokens[iNorth + 1]);
     700           2 :     const double dfSouth = CPLAtofM(papszTokens[iSouth + 1]);
     701           2 :     const double dfEast = CPLAtofM(papszTokens[iEast + 1]);
     702           2 :     const double dfWest = CPLAtofM(papszTokens[iWest + 1]);
     703           2 :     const double dfPixelXSize = (dfEast - dfWest) / nRasterXSize;
     704           2 :     const double dfPixelYSize = (dfNorth - dfSouth) / nRasterYSize;
     705             : 
     706           2 :     adfGeoTransform[0] = dfWest;
     707           2 :     adfGeoTransform[1] = dfPixelXSize;
     708           2 :     adfGeoTransform[2] = 0.0;
     709           2 :     adfGeoTransform[3] = dfNorth;
     710           2 :     adfGeoTransform[4] = 0.0;
     711           2 :     adfGeoTransform[5] = -dfPixelYSize;
     712             : 
     713           2 :     if ((i = CSLFindString(papszTokens, "null")) >= 0 && i + 1 < nTokens)
     714             :     {
     715           0 :         const char *pszNoData = papszTokens[i + 1];
     716             : 
     717           0 :         bNoDataSet = true;
     718           0 :         dfNoDataValue = CPLAtofM(pszNoData);
     719           0 :         if (pszDataType == nullptr &&
     720           0 :             (strchr(pszNoData, '.') != nullptr ||
     721           0 :              strchr(pszNoData, ',') != nullptr ||
     722           0 :              std::numeric_limits<int>::min() > dfNoDataValue ||
     723           0 :              dfNoDataValue > std::numeric_limits<int>::max()))
     724             :         {
     725           0 :             eDataType = GDT_Float32;
     726             :         }
     727           0 :         if (eDataType == GDT_Float32)
     728             :         {
     729           0 :             dfNoDataValue = MapNoDataToFloat(dfNoDataValue);
     730             :         }
     731             :     }
     732             : 
     733           2 :     if ((i = CSLFindString(papszTokens, "type")) >= 0 && i + 1 < nTokens)
     734             :     {
     735           0 :         const char *pszType = papszTokens[i + 1];
     736           0 :         if (EQUAL(pszType, "int"))
     737           0 :             eDataType = GDT_Int32;
     738           0 :         else if (EQUAL(pszType, "float"))
     739           0 :             eDataType = GDT_Float32;
     740           0 :         else if (EQUAL(pszType, "double"))
     741           0 :             eDataType = GDT_Float64;
     742             :         else
     743             :         {
     744           0 :             ReportError(CE_Warning, CPLE_AppDefined,
     745             :                         "Invalid value for type parameter : %s", pszType);
     746             :         }
     747             :     }
     748             : 
     749           2 :     CSLDestroy(papszTokens);
     750             : 
     751           2 :     return TRUE;
     752             : }
     753             : 
     754             : /************************************************************************/
     755             : /*                                Open()                                */
     756             : /************************************************************************/
     757             : 
     758           7 : GDALDataset *ISGDataset::Open(GDALOpenInfo *poOpenInfo)
     759             : {
     760             : #ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
     761             :     // During fuzzing, do not use Identify to reject crazy content.
     762           7 :     if (!Identify(poOpenInfo))
     763           0 :         return nullptr;
     764             : #endif
     765             : 
     766           7 :     return CommonOpen(poOpenInfo, FORMAT_ISG);
     767             : }
     768             : 
     769             : /************************************************************************/
     770             : /*                          ParseHeader()                               */
     771             : /************************************************************************/
     772             : 
     773           7 : int ISGDataset::ParseHeader(const char *pszHeader, const char *)
     774             : {
     775             :     // See https://www.isgeoid.polimi.it/Geoid/ISG_format_v10_20160121.pdf
     776             :     //     https://www.isgeoid.polimi.it/Geoid/ISG_format_v101_20180915.pdf
     777             :     //     https://www.isgeoid.polimi.it/Geoid/ISG_format_v20_20200625.pdf
     778             : 
     779          14 :     CPLStringList aosLines(CSLTokenizeString2(pszHeader, "\n\r", 0));
     780          14 :     CPLString osLatMin;
     781          14 :     CPLString osLatMax;
     782          14 :     CPLString osLonMin;
     783          14 :     CPLString osLonMax;
     784          14 :     CPLString osDeltaLat;
     785          14 :     CPLString osDeltaLon;
     786          14 :     CPLString osRows;
     787          14 :     CPLString osCols;
     788          14 :     CPLString osNodata;
     789          14 :     std::string osISGFormat;
     790          14 :     std::string osDataFormat;    // ISG 2.0
     791          14 :     std::string osDataOrdering;  // ISG 2.0
     792          14 :     std::string osCoordType;     // ISG 2.0
     793          14 :     std::string osCoordUnits;    // ISG 2.0
     794         151 :     for (int iLine = 0; iLine < aosLines.size(); iLine++)
     795             :     {
     796         288 :         CPLStringList aosTokens(CSLTokenizeString2(aosLines[iLine], ":=", 0));
     797         144 :         if (aosTokens.size() == 2)
     798             :         {
     799         228 :             CPLString osLeft(aosTokens[0]);
     800         114 :             osLeft.Trim();
     801         228 :             const CPLString osRight(CPLString(aosTokens[1]).Trim());
     802         114 :             if (osLeft == "lat min")
     803           7 :                 osLatMin = osRight;
     804         107 :             else if (osLeft == "lat max")
     805           7 :                 osLatMax = osRight;
     806         100 :             else if (osLeft == "lon min")
     807           7 :                 osLonMin = osRight;
     808          93 :             else if (osLeft == "lon max")
     809           7 :                 osLonMax = osRight;
     810          86 :             else if (osLeft == "delta lat")
     811           7 :                 osDeltaLat = osRight;
     812          79 :             else if (osLeft == "delta lon")
     813           7 :                 osDeltaLon = osRight;
     814          72 :             else if (osLeft == "nrows")
     815           7 :                 osRows = osRight;
     816          65 :             else if (osLeft == "ncols")
     817           7 :                 osCols = osRight;
     818          58 :             else if (osLeft == "nodata")
     819           7 :                 osNodata = osRight;
     820          51 :             else if (osLeft == "model name")
     821           7 :                 SetMetadataItem("MODEL_NAME", osRight);
     822          44 :             else if (osLeft == "model type")
     823           7 :                 SetMetadataItem("MODEL_TYPE", osRight);
     824          37 :             else if (osLeft == "units" || osLeft == "data units")
     825           7 :                 osUnits = osRight;
     826          30 :             else if (osLeft == "ISG format")
     827           7 :                 osISGFormat = osRight;
     828          23 :             else if (osLeft == "data format")
     829           1 :                 osDataFormat = osRight;
     830          22 :             else if (osLeft == "data ordering")
     831           1 :                 osDataOrdering = osRight;
     832          21 :             else if (osLeft == "coord type")
     833           1 :                 osCoordType = osRight;
     834          20 :             else if (osLeft == "coord units")
     835           1 :                 osCoordUnits = osRight;
     836             :         }
     837             :     }
     838             :     const double dfVersion =
     839           7 :         osISGFormat.empty() ? 0.0 : CPLAtof(osISGFormat.c_str());
     840          21 :     if (osLatMin.empty() || osLatMax.empty() || osLonMin.empty() ||
     841          21 :         osLonMax.empty() || osDeltaLat.empty() || osDeltaLon.empty() ||
     842          21 :         osRows.empty() || osCols.empty())
     843             :     {
     844           0 :         return FALSE;
     845             :     }
     846           7 :     if (!osDataFormat.empty() && osDataFormat != "grid")
     847             :     {
     848           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     849             :                  "ISG: data format = %s not supported", osDataFormat.c_str());
     850           0 :         return FALSE;
     851             :     }
     852           7 :     if (!osDataOrdering.empty() && osDataOrdering != "N-to-S, W-to-E")
     853             :     {
     854           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     855             :                  "ISG: data ordering = %s not supported",
     856             :                  osDataOrdering.c_str());
     857           0 :         return FALSE;
     858             :     }
     859           7 :     if (!osCoordType.empty() && osCoordType != "geodetic")
     860             :     {
     861           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     862             :                  "ISG: coord type = %s not supported", osCoordType.c_str());
     863           0 :         return FALSE;
     864             :     }
     865           7 :     if (!osCoordUnits.empty() && osCoordUnits != "deg")
     866             :     {
     867           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     868             :                  "ISG: coord units = %s not supported", osCoordUnits.c_str());
     869           0 :         return FALSE;
     870             :     }
     871           7 :     double dfLatMin = CPLAtof(osLatMin);
     872           7 :     double dfLatMax = CPLAtof(osLatMax);
     873           7 :     double dfLonMin = CPLAtof(osLonMin);
     874           7 :     double dfLonMax = CPLAtof(osLonMax);
     875           7 :     double dfDeltaLon = CPLAtof(osDeltaLon);
     876           7 :     double dfDeltaLat = CPLAtof(osDeltaLat);
     877           7 :     if (dfVersion >= 2.0)
     878             :     {
     879           1 :         dfLatMin -= dfDeltaLat / 2.0;
     880           1 :         dfLatMax += dfDeltaLat / 2.0;
     881           1 :         dfLonMin -= dfDeltaLon / 2.0;
     882           1 :         dfLonMax += dfDeltaLon / 2.0;
     883             :     }
     884           7 :     const int nRows = atoi(osRows);
     885           7 :     const int nCols = atoi(osCols);
     886           7 :     if (nRows <= 0 || nCols <= 0 ||
     887           7 :         !(dfDeltaLat > 0 && dfDeltaLon > 0 && dfDeltaLat < 180 &&
     888           7 :           dfDeltaLon < 360))
     889             :     {
     890           0 :         return FALSE;
     891             :     }
     892             : 
     893             :     // Correct rounding errors.
     894             : 
     895          12 :     const auto TryRoundTo = [](double &dfDelta, double dfRoundedDelta,
     896             :                                double &dfMin, double &dfMax, int nVals,
     897             :                                double dfRelTol)
     898             :     {
     899          12 :         double dfMinTry = dfMin;
     900          12 :         double dfMaxTry = dfMax;
     901          12 :         double dfDeltaTry = dfDelta;
     902          12 :         if (dfRoundedDelta != dfDelta &&
     903           4 :             fabs(fabs(dfMin / dfRoundedDelta) -
     904           4 :                  (floor(fabs(dfMin / dfRoundedDelta)) + 0.5)) < dfRelTol &&
     905           4 :             fabs(fabs(dfMax / dfRoundedDelta) -
     906           4 :                  (floor(fabs(dfMax / dfRoundedDelta)) + 0.5)) < dfRelTol)
     907             :         {
     908             :             {
     909           3 :                 double dfVal = (floor(fabs(dfMin / dfRoundedDelta)) + 0.5) *
     910             :                                dfRoundedDelta;
     911           3 :                 dfMinTry = (dfMin < 0) ? -dfVal : dfVal;
     912             :             }
     913             :             {
     914           3 :                 double dfVal = (floor(fabs(dfMax / dfRoundedDelta)) + 0.5) *
     915             :                                dfRoundedDelta;
     916           3 :                 dfMaxTry = (dfMax < 0) ? -dfVal : dfVal;
     917             :             }
     918           3 :             dfDeltaTry = dfRoundedDelta;
     919             :         }
     920           9 :         else if (dfRoundedDelta != dfDelta &&
     921           1 :                  fabs(fabs(dfMin / dfRoundedDelta) -
     922           1 :                       (floor(fabs(dfMin / dfRoundedDelta) + 0.5) + 0.)) <
     923           0 :                      dfRelTol &&
     924           0 :                  fabs(fabs(dfMax / dfRoundedDelta) -
     925           0 :                       (floor(fabs(dfMax / dfRoundedDelta) + 0.5) + 0.)) <
     926             :                      dfRelTol)
     927             :         {
     928             :             {
     929           0 :                 double dfVal =
     930           0 :                     (floor(fabs(dfMin / dfRoundedDelta) + 0.5) + 0.) *
     931             :                     dfRoundedDelta;
     932           0 :                 dfMinTry = (dfMin < 0) ? -dfVal : dfVal;
     933             :             }
     934             :             {
     935           0 :                 double dfVal =
     936           0 :                     (floor(fabs(dfMax / dfRoundedDelta) + 0.5) + 0.) *
     937             :                     dfRoundedDelta;
     938           0 :                 dfMaxTry = (dfMax < 0) ? -dfVal : dfVal;
     939             :             }
     940           0 :             dfDeltaTry = dfRoundedDelta;
     941             :         }
     942          12 :         if (fabs(dfMinTry + dfDeltaTry * nVals - dfMaxTry) <
     943          12 :             dfRelTol * dfDeltaTry)
     944             :         {
     945           8 :             dfMin = dfMinTry;
     946           8 :             dfMax = dfMaxTry;
     947           8 :             dfDelta = dfDeltaTry;
     948           8 :             return true;
     949             :         }
     950           4 :         return false;
     951             :     };
     952             : 
     953             :     const double dfRoundedDeltaLon =
     954           7 :         (osDeltaLon == "0.0167" ||
     955           6 :          (dfDeltaLon < 1 &&
     956           6 :           fabs(1. / dfDeltaLon - floor(1. / dfDeltaLon + 0.5)) < 0.06))
     957          13 :             ? 1. / floor(1. / dfDeltaLon + 0.5)
     958           7 :             : dfDeltaLon;
     959             : 
     960             :     const double dfRoundedDeltaLat =
     961           7 :         (osDeltaLat == "0.0167" ||
     962           3 :          (dfDeltaLat < 1 &&
     963           3 :           fabs(1. / dfDeltaLat - floor(1. / dfDeltaLat + 0.5)) < 0.06))
     964          10 :             ? 1. / floor(1. / dfDeltaLat + 0.5)
     965           7 :             : dfDeltaLat;
     966             : 
     967           7 :     bool bOK = TryRoundTo(dfDeltaLon, dfRoundedDeltaLon, dfLonMin, dfLonMax,
     968          10 :                           nCols, 1e-2) &&
     969           3 :                TryRoundTo(dfDeltaLat, dfRoundedDeltaLat, dfLatMin, dfLatMax,
     970           7 :                           nRows, 1e-2);
     971           7 :     if (!bOK && osDeltaLon == "0.0167" && osDeltaLat == "0.0167")
     972             :     {
     973             :         // For https://www.isgeoid.polimi.it/Geoid/America/Argentina/public/GEOIDEAR16_20160419.isg
     974           1 :         bOK =
     975           2 :             TryRoundTo(dfDeltaLon, 0.016667, dfLonMin, dfLonMax, nCols, 1e-1) &&
     976           1 :             TryRoundTo(dfDeltaLat, 0.016667, dfLatMin, dfLatMax, nRows, 1e-1);
     977             :     }
     978           7 :     if (!bOK)
     979             :     {
     980             :         // 0.005 is what would be needed for the above GEOIDEAR16_20160419.isg
     981             :         // file without the specific fine tuning done.
     982           6 :         if ((fabs((dfLonMax - dfLonMin) / nCols - dfDeltaLon) <
     983           4 :                  0.005 * dfDeltaLon &&
     984           1 :              fabs((dfLatMax - dfLatMin) / nRows - dfDeltaLat) <
     985           5 :                  0.005 * dfDeltaLat) ||
     986           2 :             CPLTestBool(
     987             :                 CPLGetConfigOption("ISG_SKIP_GEOREF_CONSISTENCY_CHECK", "NO")))
     988             :         {
     989           2 :             CPLError(CE_Warning, CPLE_AppDefined,
     990             :                      "Georeference might be slightly approximate due to "
     991             :                      "rounding of coordinates and resolution in file header.");
     992           2 :             dfDeltaLon = (dfLonMax - dfLonMin) / nCols;
     993           2 :             dfDeltaLat = (dfLatMax - dfLatMin) / nRows;
     994             :         }
     995             :         else
     996             :         {
     997           1 :             CPLError(CE_Failure, CPLE_AppDefined,
     998             :                      "Inconsistent extent/resolution/raster dimension, or "
     999             :                      "rounding of coordinates and resolution in file header "
    1000             :                      "higher than accepted. You may skip this consistency "
    1001             :                      "check by setting the ISG_SKIP_GEOREF_CONSISTENCY_CHECK "
    1002             :                      "configuration option to YES.");
    1003           1 :             return false;
    1004             :         }
    1005             :     }
    1006           6 :     nRasterXSize = nCols;
    1007           6 :     nRasterYSize = nRows;
    1008           6 :     adfGeoTransform[0] = dfLonMin;
    1009           6 :     adfGeoTransform[1] = dfDeltaLon;
    1010           6 :     adfGeoTransform[2] = 0.0;
    1011           6 :     adfGeoTransform[3] = dfLatMax;
    1012           6 :     adfGeoTransform[4] = 0.0;
    1013           6 :     adfGeoTransform[5] = -dfDeltaLat;
    1014           6 :     if (!osNodata.empty())
    1015             :     {
    1016           6 :         bNoDataSet = true;
    1017           6 :         dfNoDataValue = MapNoDataToFloat(CPLAtof(osNodata));
    1018             :     }
    1019           6 :     return TRUE;
    1020             : }
    1021             : 
    1022             : /************************************************************************/
    1023             : /*                           CommonOpen()                               */
    1024             : /************************************************************************/
    1025             : 
    1026         164 : GDALDataset *AAIGDataset::CommonOpen(GDALOpenInfo *poOpenInfo,
    1027             :                                      GridFormat eFormat)
    1028             : {
    1029         164 :     if (poOpenInfo->fpL == nullptr)
    1030           0 :         return nullptr;
    1031             : 
    1032             :     // Create a corresponding GDALDataset.
    1033         164 :     AAIGDataset *poDS = nullptr;
    1034             : 
    1035         164 :     if (eFormat == FORMAT_AAIG)
    1036         155 :         poDS = new AAIGDataset();
    1037           9 :     else if (eFormat == FORMAT_GRASSASCII)
    1038           2 :         poDS = new GRASSASCIIDataset();
    1039             :     else
    1040             :     {
    1041           7 :         poDS = new ISGDataset();
    1042           7 :         poDS->eDataType = GDT_Float32;
    1043             :     }
    1044             : 
    1045         173 :     const char *pszDataTypeOption = eFormat == FORMAT_AAIG ? "AAIGRID_DATATYPE"
    1046             :                                     : eFormat == FORMAT_GRASSASCII
    1047           9 :                                         ? "GRASSASCIIGRID_DATATYPE"
    1048             :                                         : nullptr;
    1049             : 
    1050             :     const char *pszDataType =
    1051         164 :         pszDataTypeOption ? CPLGetConfigOption(pszDataTypeOption, nullptr)
    1052         164 :                           : nullptr;
    1053         164 :     if (pszDataType == nullptr)
    1054             :     {
    1055             :         pszDataType =
    1056         162 :             CSLFetchNameValue(poOpenInfo->papszOpenOptions, "DATATYPE");
    1057             :     }
    1058         164 :     if (pszDataType != nullptr)
    1059             :     {
    1060           6 :         poDS->eDataType = GDALGetDataTypeByName(pszDataType);
    1061           6 :         if (!(poDS->eDataType == GDT_Int32 || poDS->eDataType == GDT_Float32 ||
    1062           6 :               poDS->eDataType == GDT_Float64))
    1063             :         {
    1064           0 :             ReportError(poOpenInfo->pszFilename, CE_Warning, CPLE_NotSupported,
    1065             :                         "Unsupported value for %s : %s", pszDataTypeOption,
    1066             :                         pszDataType);
    1067           0 :             poDS->eDataType = GDT_Int32;
    1068           0 :             pszDataType = nullptr;
    1069             :         }
    1070             :     }
    1071             : 
    1072             :     // Parse the header.
    1073         164 :     if (!poDS->ParseHeader((const char *)poOpenInfo->pabyHeader, pszDataType))
    1074             :     {
    1075           1 :         delete poDS;
    1076           1 :         return nullptr;
    1077             :     }
    1078             : 
    1079         163 :     poDS->fp = poOpenInfo->fpL;
    1080         163 :     poOpenInfo->fpL = nullptr;
    1081             : 
    1082             :     // Find the start of real data.
    1083         163 :     int nStartOfData = 0;
    1084             : 
    1085         163 :     if (eFormat == FORMAT_ISG)
    1086             :     {
    1087             :         const char *pszEOH =
    1088           6 :             strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
    1089             :                    "end_of_head");
    1090           6 :         if (pszEOH == nullptr)
    1091             :         {
    1092           0 :             delete poDS;
    1093           0 :             return nullptr;
    1094             :         }
    1095         380 :         for (int i = 0; pszEOH[i]; i++)
    1096             :         {
    1097         380 :             if (pszEOH[i] == '\n' || pszEOH[i] == '\r')
    1098             :             {
    1099           6 :                 nStartOfData =
    1100           6 :                     static_cast<int>(pszEOH - reinterpret_cast<const char *>(
    1101           6 :                                                   poOpenInfo->pabyHeader)) +
    1102             :                     i;
    1103           6 :                 break;
    1104             :             }
    1105             :         }
    1106           6 :         if (nStartOfData == 0)
    1107             :         {
    1108           0 :             delete poDS;
    1109           0 :             return nullptr;
    1110             :         }
    1111           6 :         if (poOpenInfo->pabyHeader[nStartOfData] == '\n' ||
    1112           0 :             poOpenInfo->pabyHeader[nStartOfData] == '\r')
    1113             :         {
    1114           6 :             nStartOfData++;
    1115             :         }
    1116             : 
    1117           6 :         poDS->m_oSRS.importFromWkt(SRS_WKT_WGS84_LAT_LONG);
    1118             :     }
    1119             :     else
    1120             :     {
    1121       18544 :         for (int i = 2; true; i++)
    1122             :         {
    1123       18544 :             if (poOpenInfo->pabyHeader[i] == '\0')
    1124             :             {
    1125           0 :                 ReportError(poOpenInfo->pszFilename, CE_Failure,
    1126             :                             CPLE_AppDefined,
    1127             :                             "Couldn't find data values in ASCII Grid file.");
    1128           0 :                 delete poDS;
    1129           0 :                 return nullptr;
    1130             :             }
    1131             : 
    1132       18544 :             if (poOpenInfo->pabyHeader[i - 1] == '\n' ||
    1133       17688 :                 poOpenInfo->pabyHeader[i - 2] == '\n' ||
    1134       16989 :                 poOpenInfo->pabyHeader[i - 1] == '\r' ||
    1135       16905 :                 poOpenInfo->pabyHeader[i - 2] == '\r')
    1136             :             {
    1137        1639 :                 if ((!isalpha(static_cast<unsigned char>(
    1138        1639 :                          poOpenInfo->pabyHeader[i])) ||
    1139             :                      // null seems to be specific of D12 software
    1140             :                      // See https://github.com/OSGeo/gdal/issues/5095
    1141        1401 :                      (i + 5 < poOpenInfo->nHeaderBytes &&
    1142        1401 :                       memcmp(poOpenInfo->pabyHeader + i, "null ", 5) == 0) ||
    1143        1399 :                      (i + 4 < poOpenInfo->nHeaderBytes &&
    1144        1399 :                       EQUALN(reinterpret_cast<const char *>(
    1145             :                                  poOpenInfo->pabyHeader + i),
    1146         241 :                              "nan ", 4))) &&
    1147         241 :                     poOpenInfo->pabyHeader[i] != '\n' &&
    1148         157 :                     poOpenInfo->pabyHeader[i] != '\r')
    1149             :                 {
    1150         157 :                     nStartOfData = i;
    1151             : 
    1152             :                     // Beginning of real data found.
    1153         157 :                     break;
    1154             :                 }
    1155             :             }
    1156             :         }
    1157             :     }
    1158             : 
    1159             :     // Recognize the type of data.
    1160         163 :     CPLAssert(nullptr != poDS->fp);
    1161             : 
    1162         163 :     if (pszDataType == nullptr && poDS->eDataType != GDT_Float32 &&
    1163         143 :         poDS->eDataType != GDT_Float64)
    1164             :     {
    1165             :         // Allocate 100K chunk + 1 extra byte for NULL character.
    1166         142 :         constexpr size_t nChunkSize = 1024 * 100;
    1167             :         GByte *pabyChunk = static_cast<GByte *>(
    1168         142 :             VSI_CALLOC_VERBOSE(nChunkSize + 1, sizeof(GByte)));
    1169         142 :         if (pabyChunk == nullptr)
    1170             :         {
    1171           0 :             delete poDS;
    1172           0 :             return nullptr;
    1173             :         }
    1174         142 :         pabyChunk[nChunkSize] = '\0';
    1175             : 
    1176         142 :         if (VSIFSeekL(poDS->fp, nStartOfData, SEEK_SET) < 0)
    1177             :         {
    1178           0 :             delete poDS;
    1179           0 :             VSIFree(pabyChunk);
    1180           0 :             return nullptr;
    1181             :         }
    1182             : 
    1183             :         // Scan for dot in subsequent chunks of data.
    1184         284 :         while (!VSIFEofL(poDS->fp))
    1185             :         {
    1186         142 :             const size_t nLen = VSIFReadL(pabyChunk, 1, nChunkSize, poDS->fp);
    1187             : 
    1188      173607 :             for (size_t i = 0; i < nLen; i++)
    1189             :             {
    1190      173511 :                 const GByte ch = pabyChunk[i];
    1191      173511 :                 if (ch == '.' || ch == ',' || ch == 'e' || ch == 'E')
    1192             :                 {
    1193          46 :                     poDS->eDataType = GDT_Float32;
    1194          46 :                     break;
    1195             :                 }
    1196             :             }
    1197             :         }
    1198             : 
    1199             :         // Deallocate chunk.
    1200         142 :         VSIFree(pabyChunk);
    1201             :     }
    1202             : 
    1203             :     // Create band information objects.
    1204         163 :     AAIGRasterBand *band = new AAIGRasterBand(poDS, nStartOfData);
    1205         163 :     poDS->SetBand(1, band);
    1206         163 :     if (band->panLineOffset == nullptr)
    1207             :     {
    1208           0 :         delete poDS;
    1209           0 :         return nullptr;
    1210             :     }
    1211         163 :     if (!poDS->osUnits.empty())
    1212             :     {
    1213           6 :         poDS->GetRasterBand(1)->SetUnitType(poDS->osUnits);
    1214             :     }
    1215             : 
    1216             :     // Try to read projection file.
    1217         163 :     char *const pszDirname = CPLStrdup(CPLGetPath(poOpenInfo->pszFilename));
    1218             :     char *const pszBasename =
    1219         163 :         CPLStrdup(CPLGetBasename(poOpenInfo->pszFilename));
    1220             : 
    1221         163 :     poDS->osPrjFilename = CPLFormFilename(pszDirname, pszBasename, "prj");
    1222         163 :     int nRet = 0;
    1223             :     {
    1224             :         VSIStatBufL sStatBuf;
    1225         163 :         nRet = VSIStatL(poDS->osPrjFilename, &sStatBuf);
    1226             :     }
    1227         163 :     if (nRet != 0 && VSIIsCaseSensitiveFS(poDS->osPrjFilename))
    1228             :     {
    1229         115 :         poDS->osPrjFilename = CPLFormFilename(pszDirname, pszBasename, "PRJ");
    1230             : 
    1231             :         VSIStatBufL sStatBuf;
    1232         115 :         nRet = VSIStatL(poDS->osPrjFilename, &sStatBuf);
    1233             :     }
    1234             : 
    1235         163 :     if (nRet == 0)
    1236             :     {
    1237          50 :         poDS->papszPrj = CSLLoad(poDS->osPrjFilename);
    1238             : 
    1239          50 :         CPLDebug("AAIGrid", "Loaded SRS from %s", poDS->osPrjFilename.c_str());
    1240             : 
    1241         100 :         OGRSpatialReference oSRS;
    1242          50 :         oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    1243          50 :         if (oSRS.importFromESRI(poDS->papszPrj) == OGRERR_NONE)
    1244             :         {
    1245             :             // If geographic values are in seconds, we must transform.
    1246             :             // Is there a code for minutes too?
    1247          61 :             if (oSRS.IsGeographic() &&
    1248          61 :                 EQUAL(OSR_GDS(poDS->papszPrj, "Units", ""), "DS"))
    1249             :             {
    1250           0 :                 poDS->adfGeoTransform[0] /= 3600.0;
    1251           0 :                 poDS->adfGeoTransform[1] /= 3600.0;
    1252           0 :                 poDS->adfGeoTransform[2] /= 3600.0;
    1253           0 :                 poDS->adfGeoTransform[3] /= 3600.0;
    1254           0 :                 poDS->adfGeoTransform[4] /= 3600.0;
    1255           0 :                 poDS->adfGeoTransform[5] /= 3600.0;
    1256             :             }
    1257             : 
    1258          50 :             poDS->m_oSRS = std::move(oSRS);
    1259             :         }
    1260             :     }
    1261             : 
    1262         163 :     CPLFree(pszDirname);
    1263         163 :     CPLFree(pszBasename);
    1264             : 
    1265             :     // Initialize any PAM information.
    1266         163 :     poDS->SetDescription(poOpenInfo->pszFilename);
    1267         163 :     poDS->TryLoadXML();
    1268             : 
    1269             :     // Check for external overviews.
    1270         163 :     poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename,
    1271             :                                 poOpenInfo->GetSiblingFiles());
    1272             : 
    1273         163 :     return poDS;
    1274             : }
    1275             : 
    1276             : /************************************************************************/
    1277             : /*                          GetGeoTransform()                           */
    1278             : /************************************************************************/
    1279             : 
    1280         109 : CPLErr AAIGDataset::GetGeoTransform(double *padfTransform)
    1281             : 
    1282             : {
    1283         109 :     memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
    1284         109 :     return CE_None;
    1285             : }
    1286             : 
    1287             : /************************************************************************/
    1288             : /*                          GetSpatialRef()                             */
    1289             : /************************************************************************/
    1290             : 
    1291          72 : const OGRSpatialReference *AAIGDataset::GetSpatialRef() const
    1292             : {
    1293          72 :     return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
    1294             : }
    1295             : 
    1296             : /************************************************************************/
    1297             : /*                          CreateCopy()                                */
    1298             : /************************************************************************/
    1299             : 
    1300          49 : GDALDataset *AAIGDataset::CreateCopy(const char *pszFilename,
    1301             :                                      GDALDataset *poSrcDS, int /* bStrict */,
    1302             :                                      char **papszOptions,
    1303             :                                      GDALProgressFunc pfnProgress,
    1304             :                                      void *pProgressData)
    1305             : {
    1306          49 :     const int nBands = poSrcDS->GetRasterCount();
    1307          49 :     const int nXSize = poSrcDS->GetRasterXSize();
    1308          49 :     const int nYSize = poSrcDS->GetRasterYSize();
    1309             : 
    1310             :     // Some rudimentary checks.
    1311          49 :     if (nBands != 1)
    1312             :     {
    1313           5 :         ReportError(pszFilename, CE_Failure, CPLE_NotSupported,
    1314             :                     "AAIG driver doesn't support %d bands.  Must be 1 band.",
    1315             :                     nBands);
    1316             : 
    1317           5 :         return nullptr;
    1318             :     }
    1319             : 
    1320          44 :     if (!pfnProgress(0.0, nullptr, pProgressData))
    1321           0 :         return nullptr;
    1322             : 
    1323             :     // Create the dataset.
    1324          44 :     VSILFILE *fpImage = VSIFOpenL(pszFilename, "wt");
    1325          44 :     if (fpImage == nullptr)
    1326             :     {
    1327           3 :         ReportError(pszFilename, CE_Failure, CPLE_OpenFailed,
    1328             :                     "Unable to create file.");
    1329           3 :         return nullptr;
    1330             :     }
    1331             : 
    1332             :     // Write ASCII Grid file header.
    1333          41 :     double adfGeoTransform[6] = {};
    1334          41 :     char szHeader[2000] = {};
    1335             :     const char *pszForceCellsize =
    1336          41 :         CSLFetchNameValue(papszOptions, "FORCE_CELLSIZE");
    1337             : 
    1338          41 :     poSrcDS->GetGeoTransform(adfGeoTransform);
    1339             : 
    1340          41 :     const double dfYLLCorner =
    1341          41 :         adfGeoTransform[5] < 0
    1342          41 :             ? adfGeoTransform[3] + nYSize * adfGeoTransform[5]
    1343             :             : adfGeoTransform[3];
    1344          43 :     if (std::abs(adfGeoTransform[1] + adfGeoTransform[5]) < 0.0000001 ||
    1345          43 :         std::abs(adfGeoTransform[1] - adfGeoTransform[5]) < 0.0000001 ||
    1346           0 :         (pszForceCellsize && CPLTestBool(pszForceCellsize)))
    1347             :     {
    1348          40 :         CPLsnprintf(szHeader, sizeof(szHeader),
    1349             :                     "ncols        %d\n"
    1350             :                     "nrows        %d\n"
    1351             :                     "xllcorner    %.12f\n"
    1352             :                     "yllcorner    %.12f\n"
    1353             :                     "cellsize     %.12f\n",
    1354             :                     nXSize, nYSize, adfGeoTransform[0], dfYLLCorner,
    1355             :                     adfGeoTransform[1]);
    1356             :     }
    1357             :     else
    1358             :     {
    1359           1 :         if (pszForceCellsize == nullptr)
    1360           1 :             ReportError(pszFilename, CE_Warning, CPLE_AppDefined,
    1361             :                         "Producing a Golden Surfer style file with DX and DY "
    1362             :                         "instead of CELLSIZE since the input pixels are "
    1363             :                         "non-square.  Use the FORCE_CELLSIZE=TRUE creation "
    1364             :                         "option to force use of DX for even though this will "
    1365             :                         "be distorted.  Most ASCII Grid readers (ArcGIS "
    1366             :                         "included) do not support the DX and DY parameters.");
    1367           1 :         CPLsnprintf(szHeader, sizeof(szHeader),
    1368             :                     "ncols        %d\n"
    1369             :                     "nrows        %d\n"
    1370             :                     "xllcorner    %.12f\n"
    1371             :                     "yllcorner    %.12f\n"
    1372             :                     "dx           %.12f\n"
    1373             :                     "dy           %.12f\n",
    1374             :                     nXSize, nYSize, adfGeoTransform[0], dfYLLCorner,
    1375           1 :                     adfGeoTransform[1], fabs(adfGeoTransform[5]));
    1376             :     }
    1377             : 
    1378             :     // Builds the format string used for printing float values.
    1379          41 :     char szFormatFloat[32] = {'\0'};
    1380          41 :     strcpy(szFormatFloat, "%.20g");
    1381             :     const char *pszDecimalPrecision =
    1382          41 :         CSLFetchNameValue(papszOptions, "DECIMAL_PRECISION");
    1383             :     const char *pszSignificantDigits =
    1384          41 :         CSLFetchNameValue(papszOptions, "SIGNIFICANT_DIGITS");
    1385          41 :     bool bIgnoreSigDigits = false;
    1386          41 :     if (pszDecimalPrecision && pszSignificantDigits)
    1387             :     {
    1388           0 :         ReportError(pszFilename, CE_Warning, CPLE_AppDefined,
    1389             :                     "Conflicting precision arguments, using DECIMAL_PRECISION");
    1390           0 :         bIgnoreSigDigits = true;
    1391             :     }
    1392             :     int nPrecision;
    1393          41 :     if (pszSignificantDigits && !bIgnoreSigDigits)
    1394             :     {
    1395           2 :         nPrecision = atoi(pszSignificantDigits);
    1396           2 :         if (nPrecision >= 0)
    1397           2 :             snprintf(szFormatFloat, sizeof(szFormatFloat), "%%.%dg",
    1398             :                      nPrecision);
    1399           2 :         CPLDebug("AAIGrid", "Setting precision format: %s", szFormatFloat);
    1400             :     }
    1401          39 :     else if (pszDecimalPrecision)
    1402             :     {
    1403           2 :         nPrecision = atoi(pszDecimalPrecision);
    1404           2 :         if (nPrecision >= 0)
    1405           2 :             snprintf(szFormatFloat, sizeof(szFormatFloat), "%%.%df",
    1406             :                      nPrecision);
    1407           2 :         CPLDebug("AAIGrid", "Setting precision format: %s", szFormatFloat);
    1408             :     }
    1409             : 
    1410             :     // Handle nodata (optionally).
    1411          41 :     GDALRasterBand *poBand = poSrcDS->GetRasterBand(1);
    1412          41 :     const bool bReadAsInt = poBand->GetRasterDataType() == GDT_Byte ||
    1413          24 :                             poBand->GetRasterDataType() == GDT_Int16 ||
    1414          85 :                             poBand->GetRasterDataType() == GDT_UInt16 ||
    1415          20 :                             poBand->GetRasterDataType() == GDT_Int32;
    1416             : 
    1417             :     // Write `nodata' value to header if it is exists in source dataset
    1418          41 :     int bSuccess = FALSE;
    1419          41 :     const double dfNoData = poBand->GetNoDataValue(&bSuccess);
    1420          41 :     if (bSuccess)
    1421             :     {
    1422           2 :         snprintf(szHeader + strlen(szHeader),
    1423           2 :                  sizeof(szHeader) - strlen(szHeader), "%s", "NODATA_value ");
    1424           2 :         if (bReadAsInt)
    1425           0 :             snprintf(szHeader + strlen(szHeader),
    1426           0 :                      sizeof(szHeader) - strlen(szHeader), "%d",
    1427             :                      static_cast<int>(dfNoData));
    1428             :         else
    1429           2 :             CPLsnprintf(szHeader + strlen(szHeader),
    1430           2 :                         sizeof(szHeader) - strlen(szHeader), szFormatFloat,
    1431             :                         dfNoData);
    1432           2 :         snprintf(szHeader + strlen(szHeader),
    1433           2 :                  sizeof(szHeader) - strlen(szHeader), "%s", "\n");
    1434             :     }
    1435             : 
    1436          41 :     if (VSIFWriteL(szHeader, strlen(szHeader), 1, fpImage) != 1)
    1437             :     {
    1438           3 :         CPL_IGNORE_RET_VAL(VSIFCloseL(fpImage));
    1439           3 :         return nullptr;
    1440             :     }
    1441             : 
    1442             :     // Loop over image, copying image data.
    1443             : 
    1444             :     // Write scanlines to output file
    1445             :     int *panScanline = bReadAsInt
    1446          38 :                            ? static_cast<int *>(CPLMalloc(sizeof(int) * nXSize))
    1447          38 :                            : nullptr;
    1448             : 
    1449             :     double *padfScanline =
    1450          38 :         bReadAsInt ? nullptr
    1451          14 :                    : static_cast<double *>(CPLMalloc(sizeof(double) * nXSize));
    1452             : 
    1453          38 :     CPLErr eErr = CE_None;
    1454             : 
    1455          38 :     bool bHasOutputDecimalDot = false;
    1456         548 :     for (int iLine = 0; eErr == CE_None && iLine < nYSize; iLine++)
    1457             :     {
    1458        1020 :         CPLString osBuf;
    1459         510 :         const int iSrcLine =
    1460         510 :             adfGeoTransform[5] < 0 ? iLine : nYSize - 1 - iLine;
    1461         510 :         eErr = poBand->RasterIO(GF_Read, 0, iSrcLine, nXSize, 1,
    1462             :                                 bReadAsInt ? static_cast<void *>(panScanline)
    1463             :                                            : static_cast<void *>(padfScanline),
    1464             :                                 nXSize, 1, bReadAsInt ? GDT_Int32 : GDT_Float64,
    1465             :                                 0, 0, nullptr);
    1466             : 
    1467         510 :         if (bReadAsInt)
    1468             :         {
    1469       14657 :             for (int iPixel = 0; iPixel < nXSize; iPixel++)
    1470             :             {
    1471       14316 :                 snprintf(szHeader, sizeof(szHeader), "%d", panScanline[iPixel]);
    1472       14316 :                 osBuf += szHeader;
    1473       14316 :                 osBuf += ' ';
    1474       14316 :                 if ((iPixel > 0 && (iPixel % 1024) == 0) ||
    1475       14316 :                     iPixel == nXSize - 1)
    1476             :                 {
    1477         696 :                     if (VSIFWriteL(osBuf, static_cast<int>(osBuf.size()), 1,
    1478         696 :                                    fpImage) != 1)
    1479             :                     {
    1480           7 :                         eErr = CE_Failure;
    1481           7 :                         ReportError(pszFilename, CE_Failure, CPLE_AppDefined,
    1482             :                                     "Write failed, disk full?");
    1483           7 :                         break;
    1484             :                     }
    1485         341 :                     osBuf = "";
    1486             :                 }
    1487             :             }
    1488             :         }
    1489             :         else
    1490             :         {
    1491         162 :             assert(padfScanline);
    1492             : 
    1493        2514 :             for (int iPixel = 0; iPixel < nXSize; iPixel++)
    1494             :             {
    1495        2352 :                 CPLsnprintf(szHeader, sizeof(szHeader), szFormatFloat,
    1496        2352 :                             padfScanline[iPixel]);
    1497             : 
    1498             :                 // Make sure that as least one value has a decimal point (#6060)
    1499        2352 :                 if (!bHasOutputDecimalDot)
    1500             :                 {
    1501          14 :                     if (strchr(szHeader, '.') || strchr(szHeader, 'e') ||
    1502          11 :                         strchr(szHeader, 'E'))
    1503             :                     {
    1504           3 :                         bHasOutputDecimalDot = true;
    1505             :                     }
    1506          11 :                     else if (!CPLIsInf(padfScanline[iPixel]) &&
    1507          11 :                              !CPLIsNan(padfScanline[iPixel]))
    1508             :                     {
    1509          11 :                         strcat(szHeader, ".0");
    1510          11 :                         bHasOutputDecimalDot = true;
    1511             :                     }
    1512             :                 }
    1513             : 
    1514        2352 :                 osBuf += szHeader;
    1515        2352 :                 osBuf += ' ';
    1516        2352 :                 if ((iPixel > 0 && (iPixel % 1024) == 0) ||
    1517        2352 :                     iPixel == nXSize - 1)
    1518             :                 {
    1519         324 :                     if (VSIFWriteL(osBuf, static_cast<int>(osBuf.size()), 1,
    1520         324 :                                    fpImage) != 1)
    1521             :                     {
    1522           0 :                         eErr = CE_Failure;
    1523           0 :                         ReportError(pszFilename, CE_Failure, CPLE_AppDefined,
    1524             :                                     "Write failed, disk full?");
    1525           0 :                         break;
    1526             :                     }
    1527         162 :                     osBuf = "";
    1528             :                 }
    1529             :             }
    1530             :         }
    1531         510 :         if (VSIFWriteL("\n", 1, 1, fpImage) != 1)
    1532           1 :             eErr = CE_Failure;
    1533             : 
    1534        1013 :         if (eErr == CE_None &&
    1535         503 :             !pfnProgress((iLine + 1) / static_cast<double>(nYSize), nullptr,
    1536             :                          pProgressData))
    1537             :         {
    1538           0 :             eErr = CE_Failure;
    1539           0 :             ReportError(pszFilename, CE_Failure, CPLE_UserInterrupt,
    1540             :                         "User terminated CreateCopy()");
    1541             :         }
    1542             :     }
    1543             : 
    1544          38 :     CPLFree(panScanline);
    1545          38 :     CPLFree(padfScanline);
    1546          38 :     if (VSIFCloseL(fpImage) != 0)
    1547           0 :         eErr = CE_Failure;
    1548             : 
    1549          38 :     if (eErr != CE_None)
    1550           7 :         return nullptr;
    1551             : 
    1552             :     // Try to write projection file.
    1553          31 :     const char *pszOriginalProjection = poSrcDS->GetProjectionRef();
    1554          31 :     if (!EQUAL(pszOriginalProjection, ""))
    1555             :     {
    1556          21 :         char *pszDirname = CPLStrdup(CPLGetPath(pszFilename));
    1557          21 :         char *pszBasename = CPLStrdup(CPLGetBasename(pszFilename));
    1558             :         char *pszPrjFilename =
    1559          21 :             CPLStrdup(CPLFormFilename(pszDirname, pszBasename, "prj"));
    1560          21 :         VSILFILE *fp = VSIFOpenL(pszPrjFilename, "wt");
    1561          21 :         if (fp != nullptr)
    1562             :         {
    1563          42 :             OGRSpatialReference oSRS;
    1564          21 :             oSRS.importFromWkt(pszOriginalProjection);
    1565          21 :             oSRS.morphToESRI();
    1566          21 :             char *pszESRIProjection = nullptr;
    1567          21 :             oSRS.exportToWkt(&pszESRIProjection);
    1568          21 :             CPL_IGNORE_RET_VAL(VSIFWriteL(pszESRIProjection, 1,
    1569             :                                           strlen(pszESRIProjection), fp));
    1570             : 
    1571          21 :             CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
    1572          21 :             CPLFree(pszESRIProjection);
    1573             :         }
    1574             :         else
    1575             :         {
    1576           0 :             ReportError(pszFilename, CE_Failure, CPLE_FileIO,
    1577             :                         "Unable to create file %s.", pszPrjFilename);
    1578             :         }
    1579          21 :         CPLFree(pszDirname);
    1580          21 :         CPLFree(pszBasename);
    1581          21 :         CPLFree(pszPrjFilename);
    1582             :     }
    1583             : 
    1584             :     // Re-open dataset, and copy any auxiliary pam information.
    1585             : 
    1586             :     // If writing to stdout, we can't reopen it, so return
    1587             :     // a fake dataset to make the caller happy.
    1588          31 :     CPLPushErrorHandler(CPLQuietErrorHandler);
    1589             :     GDALPamDataset *poDS =
    1590          31 :         reinterpret_cast<GDALPamDataset *>(GDALOpen(pszFilename, GA_ReadOnly));
    1591          31 :     CPLPopErrorHandler();
    1592          31 :     if (poDS)
    1593             :     {
    1594          31 :         poDS->CloneInfo(poSrcDS, GCIF_PAM_DEFAULT);
    1595          31 :         return poDS;
    1596             :     }
    1597             : 
    1598           0 :     CPLErrorReset();
    1599             : 
    1600           0 :     AAIGDataset *poAAIG_DS = new AAIGDataset();
    1601           0 :     poAAIG_DS->nRasterXSize = nXSize;
    1602           0 :     poAAIG_DS->nRasterYSize = nYSize;
    1603           0 :     poAAIG_DS->nBands = 1;
    1604           0 :     poAAIG_DS->SetBand(1, new AAIGRasterBand(poAAIG_DS, 1));
    1605           0 :     return poAAIG_DS;
    1606             : }
    1607             : 
    1608             : /************************************************************************/
    1609             : /*                              OSR_GDS()                               */
    1610             : /************************************************************************/
    1611             : 
    1612          11 : static CPLString OSR_GDS(char **papszNV, const char *pszField,
    1613             :                          const char *pszDefaultValue)
    1614             : 
    1615             : {
    1616          11 :     if (papszNV == nullptr || papszNV[0] == nullptr)
    1617           0 :         return pszDefaultValue;
    1618             : 
    1619          11 :     int iLine = 0;  // Used after for.
    1620          22 :     for (; papszNV[iLine] != nullptr &&
    1621          11 :            !EQUALN(papszNV[iLine], pszField, strlen(pszField));
    1622             :          iLine++)
    1623             :     {
    1624             :     }
    1625             : 
    1626          11 :     if (papszNV[iLine] == nullptr)
    1627          11 :         return pszDefaultValue;
    1628             :     else
    1629             :     {
    1630           0 :         char **papszTokens = CSLTokenizeString(papszNV[iLine]);
    1631             : 
    1632           0 :         CPLString osResult;
    1633           0 :         if (CSLCount(papszTokens) > 1)
    1634           0 :             osResult = papszTokens[1];
    1635             :         else
    1636           0 :             osResult = pszDefaultValue;
    1637             : 
    1638           0 :         CSLDestroy(papszTokens);
    1639           0 :         return osResult;
    1640             :     }
    1641             : }
    1642             : 
    1643             : /************************************************************************/
    1644             : /*                        GDALRegister_AAIGrid()                        */
    1645             : /************************************************************************/
    1646             : 
    1647        1520 : void GDALRegister_AAIGrid()
    1648             : 
    1649             : {
    1650        1520 :     if (GDALGetDriverByName("AAIGrid") != nullptr)
    1651         301 :         return;
    1652             : 
    1653        1219 :     GDALDriver *poDriver = new GDALDriver();
    1654             : 
    1655        1219 :     poDriver->SetDescription("AAIGrid");
    1656        1219 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    1657        1219 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Arc/Info ASCII Grid");
    1658        1219 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
    1659        1219 :                               "drivers/raster/aaigrid.html");
    1660        1219 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "asc");
    1661        1219 :     poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
    1662        1219 :                               "Byte UInt16 Int16 Int32 Float32");
    1663             : 
    1664        1219 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
    1665        1219 :     poDriver->SetMetadataItem(
    1666             :         GDAL_DMD_CREATIONOPTIONLIST,
    1667             :         "<CreationOptionList>\n"
    1668             :         "   <Option name='FORCE_CELLSIZE' type='boolean' description='Force "
    1669             :         "use of CELLSIZE, default is FALSE.'/>\n"
    1670             :         "   <Option name='DECIMAL_PRECISION' type='int' description='Number of "
    1671             :         "decimal when writing floating-point numbers(%f).'/>\n"
    1672             :         "   <Option name='SIGNIFICANT_DIGITS' type='int' description='Number "
    1673             :         "of significant digits when writing floating-point numbers(%g).'/>\n"
    1674        1219 :         "</CreationOptionList>\n");
    1675        1219 :     poDriver->SetMetadataItem(GDAL_DMD_OPENOPTIONLIST,
    1676             :                               "<OpenOptionList>\n"
    1677             :                               "   <Option name='DATATYPE' type='string-select' "
    1678             :                               "description='Data type to be used.'>\n"
    1679             :                               "       <Value>Int32</Value>\n"
    1680             :                               "       <Value>Float32</Value>\n"
    1681             :                               "       <Value>Float64</Value>\n"
    1682             :                               "   </Option>\n"
    1683        1219 :                               "</OpenOptionList>\n");
    1684             : 
    1685        1219 :     poDriver->pfnOpen = AAIGDataset::Open;
    1686        1219 :     poDriver->pfnIdentify = AAIGDataset::Identify;
    1687        1219 :     poDriver->pfnCreateCopy = AAIGDataset::CreateCopy;
    1688             : 
    1689        1219 :     GetGDALDriverManager()->RegisterDriver(poDriver);
    1690             : }
    1691             : 
    1692             : /************************************************************************/
    1693             : /*                   GDALRegister_GRASSASCIIGrid()                      */
    1694             : /************************************************************************/
    1695             : 
    1696        1520 : void GDALRegister_GRASSASCIIGrid()
    1697             : 
    1698             : {
    1699        1520 :     if (GDALGetDriverByName("GRASSASCIIGrid") != nullptr)
    1700         301 :         return;
    1701             : 
    1702        1219 :     GDALDriver *poDriver = new GDALDriver();
    1703             : 
    1704        1219 :     poDriver->SetDescription("GRASSASCIIGrid");
    1705        1219 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    1706        1219 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "GRASS ASCII Grid");
    1707        1219 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
    1708        1219 :                               "drivers/raster/grassasciigrid.html");
    1709             : 
    1710        1219 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
    1711             : 
    1712        1219 :     poDriver->pfnOpen = GRASSASCIIDataset::Open;
    1713        1219 :     poDriver->pfnIdentify = GRASSASCIIDataset::Identify;
    1714             : 
    1715        1219 :     GetGDALDriverManager()->RegisterDriver(poDriver);
    1716             : }
    1717             : 
    1718             : /************************************************************************/
    1719             : /*                       GDALRegister_ISG()                             */
    1720             : /************************************************************************/
    1721             : 
    1722        1520 : void GDALRegister_ISG()
    1723             : 
    1724             : {
    1725        1520 :     if (GDALGetDriverByName("ISG") != nullptr)
    1726         301 :         return;
    1727             : 
    1728        1219 :     GDALDriver *poDriver = new GDALDriver();
    1729             : 
    1730        1219 :     poDriver->SetDescription("ISG");
    1731        1219 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    1732        1219 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
    1733        1219 :                               "International Service for the Geoid");
    1734        1219 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/isg.html");
    1735        1219 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "isg");
    1736             : 
    1737        1219 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
    1738             : 
    1739        1219 :     poDriver->pfnOpen = ISGDataset::Open;
    1740        1219 :     poDriver->pfnIdentify = ISGDataset::Identify;
    1741             : 
    1742        1219 :     GetGDALDriverManager()->RegisterDriver(poDriver);
    1743             : }

Generated by: LCOV version 1.14