LCOV - code coverage report
Current view: top level - frmts/xyz - xyzdataset.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 716 889 80.5 %
Date: 2025-07-11 10:11:13 Functions: 14 14 100.0 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  XYZ driver
       4             :  * Purpose:  GDALDataset driver for XYZ dataset.
       5             :  * Author:   Even Rouault, <even dot rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2010-2013, Even Rouault <even dot rouault at spatialys.com>
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "cpl_string.h"
      14             : #include "cpl_vsi_virtual.h"
      15             : #include "gdal_frmts.h"
      16             : #include "gdal_pam.h"
      17             : 
      18             : #include <algorithm>
      19             : #include <cmath>
      20             : #include <mutex>
      21             : #include <vector>
      22             : 
      23             : constexpr double RELATIVE_ERROR = 1e-3;
      24             : 
      25             : class XYZDataset;
      26             : 
      27             : // Global cache when we must ingest all grid points
      28             : static std::mutex gMutex;
      29             : static XYZDataset *gpoActiveDS = nullptr;
      30             : static std::vector<short> gasValues;
      31             : static std::vector<float> gafValues;
      32             : 
      33             : /************************************************************************/
      34             : /* ==================================================================== */
      35             : /*                              XYZDataset                              */
      36             : /* ==================================================================== */
      37             : /************************************************************************/
      38             : 
      39             : class XYZRasterBand;
      40             : 
      41             : class XYZDataset final : public GDALPamDataset
      42             : {
      43             :     friend class XYZRasterBand;
      44             : 
      45             :     VSILFILE *fp;
      46             :     int bHasHeaderLine;
      47             :     int nCommentLineCount;
      48             :     char chDecimalSep;
      49             :     int nXIndex;
      50             :     int nYIndex;
      51             :     int nZIndex;
      52             :     int nMinTokens;
      53             :     GIntBig nLineNum;     /* any line */
      54             :     GIntBig nDataLineNum; /* line with values (header line and empty lines
      55             :                              ignored) */
      56             :     GDALGeoTransform m_gt{};
      57             :     int bSameNumberOfValuesPerLine;
      58             :     double dfMinZ;
      59             :     double dfMaxZ;
      60             :     bool bEOF;
      61             :     bool bIngestAll = false;
      62             : 
      63             :     static int IdentifyEx(GDALOpenInfo *, int &, int &nCommentLineCount,
      64             :                           int &nXIndex, int &nYIndex, int &nZIndex);
      65             : 
      66             :     CPL_DISALLOW_COPY_ASSIGN(XYZDataset)
      67             : 
      68             :   public:
      69             :     XYZDataset();
      70             :     virtual ~XYZDataset();
      71             : 
      72             :     virtual CPLErr GetGeoTransform(GDALGeoTransform &gt) const override;
      73             : 
      74             :     static GDALDataset *Open(GDALOpenInfo *);
      75             :     static int Identify(GDALOpenInfo *);
      76             :     static GDALDataset *CreateCopy(const char *pszFilename,
      77             :                                    GDALDataset *poSrcDS, int bStrict,
      78             :                                    char **papszOptions,
      79             :                                    GDALProgressFunc pfnProgress,
      80             :                                    void *pProgressData);
      81             : };
      82             : 
      83             : /************************************************************************/
      84             : /* ==================================================================== */
      85             : /*                            XYZRasterBand                             */
      86             : /* ==================================================================== */
      87             : /************************************************************************/
      88             : 
      89             : class XYZRasterBand final : public GDALPamRasterBand
      90             : {
      91             :     friend class XYZDataset;
      92             : 
      93             :     int nLastYOff;
      94             : 
      95             :   public:
      96             :     XYZRasterBand(XYZDataset *, int, GDALDataType);
      97             : 
      98             :     virtual CPLErr IReadBlock(int, int, void *) override;
      99             :     virtual double GetMinimum(int *pbSuccess = nullptr) override;
     100             :     virtual double GetMaximum(int *pbSuccess = nullptr) override;
     101             :     virtual double GetNoDataValue(int *pbSuccess = nullptr) override;
     102             : };
     103             : 
     104             : /************************************************************************/
     105             : /*                           XYZRasterBand()                            */
     106             : /************************************************************************/
     107             : 
     108          38 : XYZRasterBand::XYZRasterBand(XYZDataset *poDSIn, int nBandIn, GDALDataType eDT)
     109          38 :     : nLastYOff(-1)
     110             : {
     111          38 :     poDS = poDSIn;
     112          38 :     nBand = nBandIn;
     113             : 
     114          38 :     eDataType = eDT;
     115             : 
     116          38 :     nBlockXSize = poDSIn->GetRasterXSize();
     117          38 :     nBlockYSize = 1;
     118          38 : }
     119             : 
     120             : /************************************************************************/
     121             : /*                             IReadBlock()                             */
     122             : /************************************************************************/
     123             : 
     124         318 : CPLErr XYZRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
     125             :                                  void *pImage)
     126             : {
     127         318 :     XYZDataset *poGDS = reinterpret_cast<XYZDataset *>(poDS);
     128             : 
     129         318 :     if (poGDS->fp == nullptr)
     130           0 :         return CE_Failure;
     131             : 
     132         318 :     if (poGDS->bIngestAll)
     133             :     {
     134          12 :         CPLAssert(eDataType == GDT_Int16 || eDataType == GDT_Float32);
     135             : 
     136          24 :         std::lock_guard<std::mutex> guard(gMutex);
     137             : 
     138          12 :         if (gpoActiveDS != poGDS || (gasValues.empty() && gafValues.empty()))
     139             :         {
     140           4 :             gpoActiveDS = poGDS;
     141             : 
     142           4 :             const int nGridSize = nRasterXSize * nRasterYSize;
     143             :             try
     144             :             {
     145           4 :                 if (eDataType == GDT_Int16)
     146           2 :                     gasValues.resize(nGridSize);
     147             :                 else
     148           2 :                     gafValues.resize(nGridSize);
     149             :             }
     150           0 :             catch (const std::exception &)
     151             :             {
     152           0 :                 CPLError(CE_Failure, CPLE_OutOfMemory, "Cannot allocate grid");
     153           0 :                 return CE_Failure;
     154             :             }
     155             : 
     156           4 :             poGDS->nDataLineNum = 0;
     157           4 :             poGDS->nLineNum = 0;
     158           4 :             poGDS->bEOF = false;
     159           4 :             VSIFSeekL(poGDS->fp, 0, SEEK_SET);
     160             : 
     161           4 :             for (int i = 0; i < poGDS->nCommentLineCount; i++)
     162             :             {
     163           0 :                 if (CPLReadLine2L(poGDS->fp, 100, nullptr) == nullptr)
     164             :                 {
     165           0 :                     poGDS->bEOF = true;
     166           0 :                     return CE_Failure;
     167             :                 }
     168           0 :                 poGDS->nLineNum++;
     169             :             }
     170             : 
     171           4 :             if (poGDS->bHasHeaderLine)
     172             :             {
     173           4 :                 const char *pszLine = CPLReadLine2L(poGDS->fp, 100, nullptr);
     174           4 :                 if (pszLine == nullptr)
     175             :                 {
     176           0 :                     poGDS->bEOF = true;
     177           0 :                     return CE_Failure;
     178             :                 }
     179           4 :                 poGDS->nLineNum++;
     180             :             }
     181             : 
     182          34 :             for (int i = 0; i < nGridSize; i++)
     183             :             {
     184          30 :                 const char *pszLine = CPLReadLine2L(poGDS->fp, 100, nullptr);
     185          30 :                 if (pszLine == nullptr)
     186             :                 {
     187           0 :                     poGDS->bEOF = true;
     188           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
     189             :                              "Cannot read line " CPL_FRMT_GIB,
     190           0 :                              poGDS->nLineNum + 1);
     191           0 :                     return CE_Failure;
     192             :                 }
     193          30 :                 poGDS->nLineNum++;
     194             : 
     195          30 :                 const char *pszPtr = pszLine;
     196             :                 char ch;
     197          30 :                 int nCol = 0;
     198          30 :                 bool bLastWasSep = true;
     199          30 :                 double dfX = 0.0;
     200          30 :                 double dfY = 0.0;
     201          30 :                 double dfZ = 0.0;
     202          30 :                 int nUsefulColsFound = 0;
     203         326 :                 while ((ch = *pszPtr) != '\0')
     204             :                 {
     205         296 :                     if (ch == ' ')
     206             :                     {
     207          60 :                         if (!bLastWasSep)
     208          60 :                             nCol++;
     209          60 :                         bLastWasSep = true;
     210             :                     }
     211         236 :                     else if ((ch == ',' && poGDS->chDecimalSep != ',') ||
     212         236 :                              ch == '\t' || ch == ';')
     213             :                     {
     214           0 :                         nCol++;
     215           0 :                         bLastWasSep = true;
     216             :                     }
     217             :                     else
     218             :                     {
     219         236 :                         if (bLastWasSep)
     220             :                         {
     221          90 :                             if (nCol == poGDS->nXIndex)
     222             :                             {
     223          30 :                                 nUsefulColsFound++;
     224          30 :                                 dfX = CPLAtofDelim(pszPtr, poGDS->chDecimalSep);
     225             :                             }
     226          60 :                             else if (nCol == poGDS->nYIndex)
     227             :                             {
     228          30 :                                 nUsefulColsFound++;
     229          30 :                                 dfY = CPLAtofDelim(pszPtr, poGDS->chDecimalSep);
     230             :                             }
     231          30 :                             else if (nCol == poGDS->nZIndex)
     232             :                             {
     233          30 :                                 nUsefulColsFound++;
     234          30 :                                 dfZ = CPLAtofDelim(pszPtr, poGDS->chDecimalSep);
     235             :                             }
     236             :                         }
     237         236 :                         bLastWasSep = false;
     238             :                     }
     239         296 :                     pszPtr++;
     240             :                 }
     241             : 
     242             :                 /* Skip empty line */
     243          30 :                 if (nCol == 0 && bLastWasSep)
     244           0 :                     continue;
     245             : 
     246          30 :                 if (nUsefulColsFound != 3)
     247             :                 {
     248           0 :                     CPLError(
     249             :                         CE_Failure, CPLE_AppDefined,
     250             :                         "Unexpected number of values at line " CPL_FRMT_GIB,
     251             :                         poGDS->nLineNum);
     252           0 :                     return CE_Failure;
     253             :                 }
     254             : 
     255          30 :                 poGDS->nDataLineNum++;
     256             : 
     257             :                 const int nX = static_cast<int>(
     258          30 :                     (dfX - 0.5 * poGDS->m_gt[1] - poGDS->m_gt[0]) /
     259          30 :                         poGDS->m_gt[1] +
     260          30 :                     0.5);
     261             :                 const int nY = static_cast<int>(
     262          30 :                     (dfY - 0.5 * poGDS->m_gt[5] - poGDS->m_gt[3]) /
     263          30 :                         poGDS->m_gt[5] +
     264          30 :                     0.5);
     265          30 :                 if (nX < 0 || nX >= nRasterXSize)
     266             :                 {
     267           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
     268             :                              "Unexpected X value at line " CPL_FRMT_GIB,
     269             :                              poGDS->nLineNum);
     270           0 :                     return CE_Failure;
     271             :                 }
     272          30 :                 if (nY < 0 || nY >= nRasterYSize)
     273             :                 {
     274           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
     275             :                              "Unexpected Y value at line " CPL_FRMT_GIB,
     276             :                              poGDS->nLineNum);
     277           0 :                     return CE_Failure;
     278             :                 }
     279          30 :                 const int nIdx = nX + nY * nRasterXSize;
     280          30 :                 if (eDataType == GDT_Int16)
     281          15 :                     gasValues[nIdx] = static_cast<short>(0.5 + dfZ);
     282             :                 else
     283          15 :                     gafValues[nIdx] = static_cast<float>(dfZ);
     284             :             }
     285             :         }
     286             : 
     287          12 :         if (eDataType == GDT_Int16)
     288           6 :             memcpy(pImage,
     289           6 :                    &gasValues[static_cast<size_t>(nBlockYOff) * nBlockXSize],
     290           6 :                    sizeof(short) * nBlockXSize);
     291             :         else
     292           6 :             memcpy(pImage,
     293           6 :                    &gafValues[static_cast<size_t>(nBlockYOff) * nBlockXSize],
     294           6 :                    sizeof(float) * nBlockXSize);
     295          12 :         return CE_None;
     296             :     }
     297             : 
     298         306 :     if (pImage)
     299             :     {
     300         306 :         int bSuccess = FALSE;
     301         306 :         double dfNoDataValue = GetNoDataValue(&bSuccess);
     302         306 :         if (!bSuccess)
     303         284 :             dfNoDataValue = 0.0;
     304         306 :         GDALCopyWords(&dfNoDataValue, GDT_Float64, 0, pImage, eDataType,
     305             :                       GDALGetDataTypeSizeBytes(eDataType), nRasterXSize);
     306             :     }
     307             : 
     308             :     // Only valid if bSameNumberOfValuesPerLine.
     309         306 :     const GIntBig nLineInFile = static_cast<GIntBig>(nBlockYOff) * nBlockXSize;
     310         306 :     if ((poGDS->bSameNumberOfValuesPerLine &&
     311         284 :          poGDS->nDataLineNum > nLineInFile) ||
     312         289 :         (!poGDS->bSameNumberOfValuesPerLine &&
     313          22 :          (nLastYOff == -1 || nBlockYOff == 0)))
     314             :     {
     315          22 :         poGDS->nDataLineNum = 0;
     316          22 :         poGDS->nLineNum = 0;
     317          22 :         poGDS->bEOF = false;
     318          22 :         VSIFSeekL(poGDS->fp, 0, SEEK_SET);
     319             : 
     320          22 :         for (int i = 0; i < poGDS->nCommentLineCount; i++)
     321             :         {
     322           0 :             if (CPLReadLine2L(poGDS->fp, 100, nullptr) == nullptr)
     323             :             {
     324           0 :                 poGDS->bEOF = true;
     325           0 :                 return CE_Failure;
     326             :             }
     327           0 :             poGDS->nLineNum++;
     328             :         }
     329             : 
     330          22 :         if (poGDS->bHasHeaderLine)
     331             :         {
     332          14 :             const char *pszLine = CPLReadLine2L(poGDS->fp, 100, nullptr);
     333          14 :             if (pszLine == nullptr)
     334             :             {
     335           0 :                 poGDS->bEOF = true;
     336           0 :                 return CE_Failure;
     337             :             }
     338          14 :             poGDS->nLineNum++;
     339             :         }
     340             :     }
     341             : 
     342         306 :     if (!poGDS->bSameNumberOfValuesPerLine)
     343             :     {
     344          22 :         if (nBlockYOff < nLastYOff)
     345             :         {
     346           0 :             nLastYOff = -1;
     347           0 :             for (int iY = 0; iY < nBlockYOff; iY++)
     348             :             {
     349           0 :                 if (IReadBlock(0, iY, nullptr) != CE_None)
     350           0 :                     return CE_Failure;
     351             :             }
     352             :         }
     353             :         else
     354             :         {
     355          22 :             if (poGDS->bEOF)
     356             :             {
     357           0 :                 return CE_Failure;
     358             :             }
     359          22 :             for (int iY = nLastYOff + 1; iY < nBlockYOff; iY++)
     360             :             {
     361           0 :                 if (IReadBlock(0, iY, nullptr) != CE_None)
     362           0 :                     return CE_Failure;
     363             :             }
     364             :         }
     365             :     }
     366             :     else
     367             :     {
     368         284 :         if (poGDS->bEOF)
     369             :         {
     370           0 :             return CE_Failure;
     371             :         }
     372         310 :         while (poGDS->nDataLineNum < nLineInFile)
     373             :         {
     374          26 :             const char *pszLine = CPLReadLine2L(poGDS->fp, 100, nullptr);
     375          26 :             if (pszLine == nullptr)
     376             :             {
     377           0 :                 poGDS->bEOF = true;
     378           0 :                 return CE_Failure;
     379             :             }
     380          26 :             poGDS->nLineNum++;
     381             : 
     382          26 :             const char *pszPtr = pszLine;
     383             :             char ch;
     384          26 :             int nCol = 0;
     385          26 :             bool bLastWasSep = true;
     386         146 :             while ((ch = *pszPtr) != '\0')
     387             :             {
     388         120 :                 if (ch == ' ')
     389             :                 {
     390          40 :                     if (!bLastWasSep)
     391          40 :                         nCol++;
     392          40 :                     bLastWasSep = true;
     393             :                 }
     394          80 :                 else if ((ch == ',' && poGDS->chDecimalSep != ',') ||
     395          80 :                          ch == '\t' || ch == ';')
     396             :                 {
     397           0 :                     nCol++;
     398           0 :                     bLastWasSep = true;
     399             :                 }
     400             :                 else
     401             :                 {
     402          80 :                     bLastWasSep = false;
     403             :                 }
     404         120 :                 pszPtr++;
     405             :             }
     406             : 
     407             :             /* Skip empty line */
     408          26 :             if (nCol == 0 && bLastWasSep)
     409           6 :                 continue;
     410             : 
     411          20 :             poGDS->nDataLineNum++;
     412             :         }
     413             :     }
     414             : 
     415             :     const double dfExpectedY =
     416         306 :         poGDS->m_gt[3] + (0.5 + nBlockYOff) * poGDS->m_gt[5];
     417             : 
     418         306 :     int idx = -1;
     419             :     while (true)
     420             :     {
     421             :         int nCol;
     422             :         bool bLastWasSep;
     423           7 :         do
     424             :         {
     425       41702 :             const vsi_l_offset nOffsetBefore = VSIFTellL(poGDS->fp);
     426       41702 :             const char *pszLine = CPLReadLine2L(poGDS->fp, 100, nullptr);
     427       41702 :             if (pszLine == nullptr)
     428             :             {
     429           1 :                 poGDS->bEOF = true;
     430           1 :                 if (poGDS->bSameNumberOfValuesPerLine)
     431             :                 {
     432           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
     433             :                              "Cannot read line " CPL_FRMT_GIB,
     434           0 :                              poGDS->nLineNum + 1);
     435           0 :                     return CE_Failure;
     436             :                 }
     437             :                 else
     438             :                 {
     439           1 :                     nLastYOff = nBlockYOff;
     440             :                 }
     441           1 :                 return CE_None;
     442             :             }
     443       41701 :             poGDS->nLineNum++;
     444             : 
     445       41701 :             const char *pszPtr = pszLine;
     446             :             char ch;
     447       41701 :             nCol = 0;
     448       41701 :             bLastWasSep = true;
     449       41701 :             double dfX = 0.0;
     450       41701 :             double dfY = 0.0;
     451       41701 :             double dfZ = 0.0;
     452       41701 :             int nUsefulColsFound = 0;
     453     1353710 :             while ((ch = *pszPtr) != '\0')
     454             :             {
     455     1312010 :                 if (ch == ' ')
     456             :                 {
     457        2706 :                     if (!bLastWasSep)
     458        2586 :                         nCol++;
     459        2706 :                     bLastWasSep = true;
     460             :                 }
     461     1309300 :                 else if ((ch == ',' && poGDS->chDecimalSep != ',') ||
     462     1228500 :                          ch == '\t' || ch == ';')
     463             :                 {
     464       80802 :                     nCol++;
     465       80802 :                     bLastWasSep = true;
     466             :                 }
     467             :                 else
     468             :                 {
     469     1228500 :                     if (bLastWasSep)
     470             :                     {
     471      125082 :                         if (nCol == poGDS->nXIndex)
     472             :                         {
     473       41694 :                             nUsefulColsFound++;
     474       41694 :                             if (!poGDS->bSameNumberOfValuesPerLine)
     475          36 :                                 dfX = CPLAtofDelim(pszPtr, poGDS->chDecimalSep);
     476             :                         }
     477       83388 :                         else if (nCol == poGDS->nYIndex)
     478             :                         {
     479       41694 :                             nUsefulColsFound++;
     480       41694 :                             if (!poGDS->bSameNumberOfValuesPerLine)
     481          36 :                                 dfY = CPLAtofDelim(pszPtr, poGDS->chDecimalSep);
     482             :                         }
     483       41694 :                         else if (nCol == poGDS->nZIndex)
     484             :                         {
     485       41694 :                             nUsefulColsFound++;
     486       41694 :                             dfZ = CPLAtofDelim(pszPtr, poGDS->chDecimalSep);
     487             :                         }
     488             :                     }
     489     1228500 :                     bLastWasSep = false;
     490             :                 }
     491     1312010 :                 pszPtr++;
     492             :             }
     493       41701 :             nCol++;
     494             : 
     495       41701 :             if (nUsefulColsFound == 3)
     496             :             {
     497       41694 :                 if (poGDS->bSameNumberOfValuesPerLine)
     498             :                 {
     499       41658 :                     idx++;
     500             :                 }
     501             :                 else
     502             :                 {
     503          36 :                     if (fabs((dfY - dfExpectedY) / poGDS->m_gt[5]) >
     504             :                         RELATIVE_ERROR)
     505             :                     {
     506           5 :                         if (idx < 0)
     507             :                         {
     508             :                             const double dfYDeltaOrigin =
     509           3 :                                 dfY + 0.5 * poGDS->m_gt[5] - poGDS->m_gt[3];
     510           6 :                             if (!(fabs(dfYDeltaOrigin) > fabs(poGDS->m_gt[5]) &&
     511           6 :                                   fabs(std::round(dfYDeltaOrigin /
     512           3 :                                                   poGDS->m_gt[5]) -
     513           3 :                                        (dfYDeltaOrigin / poGDS->m_gt[5])) <=
     514             :                                       RELATIVE_ERROR))
     515             :                             {
     516           0 :                                 CPLError(CE_Failure, CPLE_AppDefined,
     517             :                                          "At line " CPL_FRMT_GIB
     518             :                                          ", found Y=%f instead of %f "
     519             :                                          "for nBlockYOff = %d",
     520             :                                          poGDS->nLineNum, dfY, dfExpectedY,
     521             :                                          nBlockYOff);
     522           0 :                                 return CE_Failure;
     523             :                             }
     524             :                         }
     525           5 :                         VSIFSeekL(poGDS->fp, nOffsetBefore, SEEK_SET);
     526           5 :                         nLastYOff = nBlockYOff;
     527           5 :                         poGDS->nLineNum--;
     528           5 :                         return CE_None;
     529             :                     }
     530             : 
     531          31 :                     idx = static_cast<int>(
     532          31 :                         (dfX - 0.5 * poGDS->m_gt[1] - poGDS->m_gt[0]) /
     533          31 :                             poGDS->m_gt[1] +
     534             :                         0.5);
     535             :                 }
     536       41689 :                 CPLAssert(idx >= 0 && idx < nRasterXSize);
     537             : 
     538       41689 :                 if (pImage)
     539             :                 {
     540       41689 :                     if (eDataType == GDT_Float32)
     541             :                     {
     542       40416 :                         reinterpret_cast<float *>(pImage)[idx] =
     543       40416 :                             static_cast<float>(dfZ);
     544             :                     }
     545        1273 :                     else if (eDataType == GDT_Int32)
     546             :                     {
     547         400 :                         reinterpret_cast<GInt32 *>(pImage)[idx] =
     548         400 :                             static_cast<GInt32>(dfZ);
     549             :                     }
     550         873 :                     else if (eDataType == GDT_Int16)
     551             :                     {
     552          13 :                         reinterpret_cast<GInt16 *>(pImage)[idx] =
     553          13 :                             static_cast<GInt16>(dfZ);
     554             :                     }
     555             :                     else
     556             :                     {
     557         860 :                         reinterpret_cast<GByte *>(pImage)[idx] =
     558         860 :                             static_cast<GByte>(dfZ);
     559             :                     }
     560             :                 }
     561             :             }
     562             :             /* Skip empty line */
     563       41696 :         } while (nCol == 1 && bLastWasSep);
     564             : 
     565       41689 :         poGDS->nDataLineNum++;
     566       41689 :         if (nCol < poGDS->nMinTokens)
     567           0 :             return CE_Failure;
     568             : 
     569       41689 :         if (idx + 1 == nRasterXSize)
     570         300 :             break;
     571       41389 :     }
     572             : 
     573         300 :     if (poGDS->bSameNumberOfValuesPerLine)
     574             :     {
     575         284 :         if (poGDS->nDataLineNum !=
     576         284 :             static_cast<GIntBig>(nBlockYOff + 1) * nBlockXSize)
     577             :         {
     578           0 :             CPLError(CE_Failure, CPLE_AssertionFailed,
     579             :                      "The file does not have the same number of values per "
     580             :                      "line as initially thought. It must be somehow corrupted");
     581           0 :             return CE_Failure;
     582             :         }
     583             :     }
     584             : 
     585         300 :     nLastYOff = nBlockYOff;
     586             : 
     587         300 :     return CE_None;
     588             : }
     589             : 
     590             : /************************************************************************/
     591             : /*                            GetMinimum()                              */
     592             : /************************************************************************/
     593             : 
     594           1 : double XYZRasterBand::GetMinimum(int *pbSuccess)
     595             : {
     596           1 :     XYZDataset *poGDS = reinterpret_cast<XYZDataset *>(poDS);
     597           1 :     if (pbSuccess)
     598           1 :         *pbSuccess = TRUE;
     599           1 :     return poGDS->dfMinZ;
     600             : }
     601             : 
     602             : /************************************************************************/
     603             : /*                            GetMaximum()                              */
     604             : /************************************************************************/
     605             : 
     606           1 : double XYZRasterBand::GetMaximum(int *pbSuccess)
     607             : {
     608           1 :     XYZDataset *poGDS = reinterpret_cast<XYZDataset *>(poDS);
     609           1 :     if (pbSuccess)
     610           1 :         *pbSuccess = TRUE;
     611           1 :     return poGDS->dfMaxZ;
     612             : }
     613             : 
     614             : /************************************************************************/
     615             : /*                          GetNoDataValue()                            */
     616             : /************************************************************************/
     617             : 
     618         318 : double XYZRasterBand::GetNoDataValue(int *pbSuccess)
     619             : {
     620         318 :     XYZDataset *poGDS = reinterpret_cast<XYZDataset *>(poDS);
     621         318 :     if (!poGDS->bSameNumberOfValuesPerLine && poGDS->dfMinZ > -32768 &&
     622          23 :         eDataType != GDT_Byte)
     623             :     {
     624           0 :         if (pbSuccess)
     625           0 :             *pbSuccess = TRUE;
     626           0 :         return (poGDS->dfMinZ > 0) ? 0 : -32768;
     627             :     }
     628         318 :     else if (!poGDS->bSameNumberOfValuesPerLine && poGDS->dfMinZ > 0 &&
     629          23 :              eDataType == GDT_Byte)
     630             :     {
     631          23 :         if (pbSuccess)
     632          23 :             *pbSuccess = TRUE;
     633          23 :         return 0;
     634             :     }
     635             : 
     636         295 :     return GDALPamRasterBand::GetNoDataValue(pbSuccess);
     637             : }
     638             : 
     639             : /************************************************************************/
     640             : /*                            ~XYZDataset()                            */
     641             : /************************************************************************/
     642             : 
     643          38 : XYZDataset::XYZDataset()
     644             :     : fp(nullptr), bHasHeaderLine(FALSE), nCommentLineCount(0),
     645             :       chDecimalSep('.'), nXIndex(-1), nYIndex(-1), nZIndex(-1), nMinTokens(0),
     646             :       nLineNum(0), nDataLineNum(GINTBIG_MAX), bSameNumberOfValuesPerLine(TRUE),
     647          38 :       dfMinZ(0), dfMaxZ(0), bEOF(false)
     648             : {
     649          38 : }
     650             : 
     651             : /************************************************************************/
     652             : /*                            ~XYZDataset()                            */
     653             : /************************************************************************/
     654             : 
     655          76 : XYZDataset::~XYZDataset()
     656             : 
     657             : {
     658          38 :     FlushCache(true);
     659          38 :     if (fp)
     660          38 :         VSIFCloseL(fp);
     661             : 
     662             :     {
     663          76 :         std::lock_guard<std::mutex> guard(gMutex);
     664          38 :         if (gpoActiveDS == this)
     665             :         {
     666           4 :             gpoActiveDS = nullptr;
     667           4 :             gasValues.clear();
     668           4 :             gafValues.clear();
     669             :         }
     670             :     }
     671          76 : }
     672             : 
     673             : /************************************************************************/
     674             : /*                             Identify()                               */
     675             : /************************************************************************/
     676             : 
     677       56940 : int XYZDataset::Identify(GDALOpenInfo *poOpenInfo)
     678             : {
     679             :     int bHasHeaderLine, nCommentLineCount;
     680             :     int nXIndex;
     681             :     int nYIndex;
     682             :     int nZIndex;
     683       56940 :     return IdentifyEx(poOpenInfo, bHasHeaderLine, nCommentLineCount, nXIndex,
     684      113874 :                       nYIndex, nZIndex);
     685             : }
     686             : 
     687             : /************************************************************************/
     688             : /*                            IdentifyEx()                              */
     689             : /************************************************************************/
     690             : 
     691       56965 : int XYZDataset::IdentifyEx(GDALOpenInfo *poOpenInfo, int &bHasHeaderLine,
     692             :                            int &nCommentLineCount, int &nXIndex, int &nYIndex,
     693             :                            int &nZIndex)
     694             : 
     695             : {
     696       56965 :     bHasHeaderLine = FALSE;
     697       56965 :     nCommentLineCount = 0;
     698             : 
     699      113927 :     CPLString osFilename(poOpenInfo->pszFilename);
     700       56965 :     if (EQUAL(CPLGetExtensionSafe(osFilename).c_str(), "GRA") &&
     701           0 :         !poOpenInfo->IsSingleAllowedDriver("XYZ"))
     702             :     {
     703             :         // IGNFHeightASCIIGRID .GRA
     704           0 :         return FALSE;
     705             :     }
     706             : 
     707       56961 :     std::unique_ptr<GDALOpenInfo> poOpenInfoToDelete;  // keep in this scope
     708             :     /*  GZipped .xyz files are common, so automagically open them */
     709             :     /*  if the /vsigzip/ has not been explicitly passed */
     710       56963 :     if (strlen(poOpenInfo->pszFilename) > 6 &&
     711       55720 :         EQUAL(poOpenInfo->pszFilename + strlen(poOpenInfo->pszFilename) - 6,
     712           0 :               "xyz.gz") &&
     713           0 :         !STARTS_WITH_CI(poOpenInfo->pszFilename, "/vsigzip/"))
     714             :     {
     715           0 :         osFilename = "/vsigzip/";
     716           0 :         osFilename += poOpenInfo->pszFilename;
     717           0 :         poOpenInfoToDelete = std::make_unique<GDALOpenInfo>(
     718           0 :             osFilename.c_str(), GA_ReadOnly, poOpenInfo->GetSiblingFiles());
     719           0 :         poOpenInfo = poOpenInfoToDelete.get();
     720             :     }
     721             : 
     722       56961 :     if (poOpenInfo->nHeaderBytes == 0)
     723             :     {
     724       53717 :         return FALSE;
     725             :     }
     726             : 
     727             :     /* -------------------------------------------------------------------- */
     728             :     /*      Check that it looks roughly as an XYZ dataset                   */
     729             :     /* -------------------------------------------------------------------- */
     730        3244 :     const char *pszData =
     731             :         reinterpret_cast<const char *>(poOpenInfo->pabyHeader);
     732             : 
     733        3245 :     if (poOpenInfo->nHeaderBytes >= 4 && STARTS_WITH(pszData, "DSAA") &&
     734           1 :         !poOpenInfo->IsSingleAllowedDriver("XYZ"))
     735             :     {
     736             :         // Do not match GSAG datasets
     737           1 :         return FALSE;
     738             :     }
     739             : 
     740             :     /* Skip comments line at the beginning such as in */
     741             :     /* http://pubs.usgs.gov/of/2003/ofr-03-230/DATA/NSLCU.XYZ */
     742        3243 :     int i = 0;
     743        3243 :     if (pszData[i] == '/')
     744             :     {
     745           0 :         nCommentLineCount++;
     746             : 
     747           0 :         i++;
     748           0 :         for (; i < poOpenInfo->nHeaderBytes; i++)
     749             :         {
     750           0 :             const char ch = pszData[i];
     751           0 :             if (ch == 13 || ch == 10)
     752             :             {
     753           0 :                 if (ch == 13 && pszData[i + 1] == 10)
     754           0 :                     i++;
     755           0 :                 if (pszData[i + 1] == '/')
     756             :                 {
     757           0 :                     nCommentLineCount++;
     758           0 :                     i++;
     759             :                 }
     760             :                 else
     761           0 :                     break;
     762             :             }
     763             :         }
     764             :     }
     765             : 
     766        3243 :     int iStartLine = i;
     767       28521 :     for (; i < poOpenInfo->nHeaderBytes; i++)
     768             :     {
     769       28469 :         const char ch = pszData[i];
     770       28469 :         if (ch == 13 || ch == 10)
     771             :         {
     772             :             break;
     773             :         }
     774       28298 :         else if (ch == ' ' || ch == ',' || ch == '\t' || ch == ';')
     775             :             ;
     776       17307 :         else if ((ch >= '0' && ch <= '9') || ch == '.' || ch == '+' ||
     777       11703 :                  ch == '-' || ch == 'e' || ch == 'E')
     778             :             ;
     779       11029 :         else if (ch == '"' || (ch >= 'a' && ch <= 'z') ||
     780        3993 :                  (ch >= 'A' && ch <= 'Z'))
     781        8009 :             bHasHeaderLine = TRUE;
     782             :         else
     783             :         {
     784        3020 :             return FALSE;
     785             :         }
     786             :     }
     787             : 
     788         223 :     nXIndex = -1;
     789         223 :     nYIndex = -1;
     790         223 :     nZIndex = -1;
     791         446 :     const char *pszColumnOrder = CSLFetchNameValueDef(
     792         223 :         poOpenInfo->papszOpenOptions, "COLUMN_ORDER", "AUTO");
     793         223 :     if (EQUAL(pszColumnOrder, "XYZ"))
     794             :     {
     795           2 :         nXIndex = 0;
     796           2 :         nYIndex = 1;
     797           2 :         nZIndex = 2;
     798           2 :         return TRUE;
     799             :     }
     800         221 :     else if (EQUAL(pszColumnOrder, "YXZ"))
     801             :     {
     802           4 :         nXIndex = 1;
     803           4 :         nYIndex = 0;
     804           4 :         nZIndex = 2;
     805           4 :         return TRUE;
     806             :     }
     807         217 :     else if (!EQUAL(pszColumnOrder, "AUTO"))
     808             :     {
     809           1 :         CPLError(CE_Failure, CPLE_IllegalArg,
     810             :                  "Option COLUMN_ORDER can only be XYZ, YXZ and AUTO."
     811             :                  "%s is not valid",
     812             :                  pszColumnOrder);
     813           1 :         return FALSE;
     814             :     }
     815             : 
     816         216 :     if (bHasHeaderLine)
     817             :     {
     818         118 :         CPLString osHeaderLine;
     819         118 :         osHeaderLine.assign(pszData + iStartLine, i - iStartLine);
     820             :         char **papszTokens =
     821         118 :             CSLTokenizeString2(osHeaderLine, " ,\t;", CSLT_HONOURSTRINGS);
     822         118 :         int nTokens = CSLCount(papszTokens);
     823         859 :         for (int iToken = 0; iToken < nTokens; iToken++)
     824             :         {
     825         741 :             const char *pszToken = papszTokens[iToken];
     826         741 :             if (EQUAL(pszToken, "x") || STARTS_WITH_CI(pszToken, "lon") ||
     827         711 :                 STARTS_WITH_CI(pszToken, "east"))
     828          30 :                 nXIndex = iToken;
     829         711 :             else if (EQUAL(pszToken, "y") || STARTS_WITH_CI(pszToken, "lat") ||
     830         683 :                      STARTS_WITH_CI(pszToken, "north"))
     831          28 :                 nYIndex = iToken;
     832         683 :             else if (EQUAL(pszToken, "z") || STARTS_WITH_CI(pszToken, "alt") ||
     833         655 :                      EQUAL(pszToken, "height"))
     834          28 :                 nZIndex = iToken;
     835             :         }
     836         118 :         CSLDestroy(papszTokens);
     837         118 :         if (nXIndex >= 0 && nYIndex >= 0 && nZIndex >= 0)
     838             :         {
     839          28 :             return TRUE;
     840             :         }
     841             :     }
     842             : 
     843         188 :     bool bHasFoundNewLine = false;
     844         188 :     bool bPrevWasSep = true;
     845         188 :     int nCols = 0;
     846         188 :     int nMaxCols = 0;
     847       16795 :     for (; i < poOpenInfo->nHeaderBytes; i++)
     848             :     {
     849       16704 :         char ch = pszData[i];
     850       16704 :         if (ch == 13 || ch == 10)
     851             :         {
     852        1145 :             bHasFoundNewLine = true;
     853        1145 :             if (!bPrevWasSep)
     854             :             {
     855         989 :                 nCols++;
     856         989 :                 if (nCols > nMaxCols)
     857          46 :                     nMaxCols = nCols;
     858             :             }
     859        1145 :             bPrevWasSep = true;
     860        1145 :             nCols = 0;
     861             :         }
     862       15559 :         else if (ch == ' ' || ch == ',' || ch == '\t' || ch == ';')
     863             :         {
     864        2568 :             if (!bPrevWasSep)
     865             :             {
     866        2507 :                 nCols++;
     867        2507 :                 if (nCols > nMaxCols)
     868         122 :                     nMaxCols = nCols;
     869             :             }
     870        2568 :             bPrevWasSep = true;
     871             :         }
     872       12991 :         else if ((ch >= '0' && ch <= '9') || ch == '.' || ch == '+' ||
     873          97 :                  ch == '-' || ch == 'e' || ch == 'E')
     874             :         {
     875       12894 :             bPrevWasSep = false;
     876             :         }
     877             :         else
     878             :         {
     879          97 :             return FALSE;
     880             :         }
     881             :     }
     882             : 
     883          91 :     return bHasFoundNewLine && nMaxCols >= 3;
     884             : }
     885             : 
     886             : /************************************************************************/
     887             : /*                                Open()                                */
     888             : /************************************************************************/
     889             : 
     890          25 : GDALDataset *XYZDataset::Open(GDALOpenInfo *poOpenInfo)
     891             : 
     892             : {
     893             :     int bHasHeaderLine;
     894          25 :     int nCommentLineCount = 0;
     895             : 
     896          25 :     int nXIndex = -1;
     897          25 :     int nYIndex = -1;
     898          25 :     int nZIndex = -1;
     899          25 :     if (!IdentifyEx(poOpenInfo, bHasHeaderLine, nCommentLineCount, nXIndex,
     900             :                     nYIndex, nZIndex))
     901           0 :         return nullptr;
     902             : 
     903          50 :     CPLString osFilename(poOpenInfo->pszFilename);
     904             : 
     905             :     /*  GZipped .xyz files are common, so automagically open them */
     906             :     /*  if the /vsigzip/ has not been explicitly passed */
     907          25 :     if (strlen(poOpenInfo->pszFilename) > 6 &&
     908          25 :         EQUAL(poOpenInfo->pszFilename + strlen(poOpenInfo->pszFilename) - 6,
     909           0 :               "xyz.gz") &&
     910           0 :         !STARTS_WITH_CI(poOpenInfo->pszFilename, "/vsigzip/"))
     911             :     {
     912           0 :         osFilename = "/vsigzip/";
     913           0 :         osFilename += poOpenInfo->pszFilename;
     914             :     }
     915             : 
     916             :     /* -------------------------------------------------------------------- */
     917             :     /*      Find dataset characteristics                                    */
     918             :     /* -------------------------------------------------------------------- */
     919          25 :     VSILFILE *fp = VSIFOpenL(osFilename.c_str(), "rb");
     920          25 :     if (fp == nullptr)
     921           0 :         return nullptr;
     922             : 
     923             :     /* For better performance of CPLReadLine2L() we create a buffered reader */
     924             :     /* (except for /vsigzip/ since it has one internally) */
     925          25 :     if (!STARTS_WITH_CI(poOpenInfo->pszFilename, "/vsigzip/"))
     926          25 :         fp = VSICreateBufferedReaderHandle(fp);
     927             : 
     928          25 :     int nMinTokens = 0;
     929             : 
     930          25 :     for (int i = 0; i < nCommentLineCount; i++)
     931             :     {
     932           0 :         if (CPLReadLine2L(fp, 100, nullptr) == nullptr)
     933             :         {
     934           0 :             VSIFCloseL(fp);
     935           0 :             return nullptr;
     936             :         }
     937             :     }
     938             : 
     939             :     /* -------------------------------------------------------------------- */
     940             :     /*      Parse header line                                               */
     941             :     /* -------------------------------------------------------------------- */
     942          25 :     if (bHasHeaderLine)
     943             :     {
     944          17 :         const char *pszLine = CPLReadLine2L(fp, 100, nullptr);
     945          17 :         if (pszLine == nullptr)
     946             :         {
     947           0 :             VSIFCloseL(fp);
     948           0 :             return nullptr;
     949             :         }
     950          17 :         if (nXIndex < 0 || nYIndex < 0 || nZIndex < 0)
     951             :         {
     952           1 :             CPLError(CE_Warning, CPLE_AppDefined,
     953             :                      "Could not find one of the X, Y or Z column names in "
     954             :                      "header line. Defaulting to the first 3 columns");
     955           1 :             nXIndex = 0;
     956           1 :             nYIndex = 1;
     957           1 :             nZIndex = 2;
     958             :         }
     959             :     }
     960           8 :     else if (nXIndex < 0 || nYIndex < 0 || nZIndex < 0)
     961             :     {
     962           7 :         nXIndex = 0;
     963           7 :         nYIndex = 1;
     964           7 :         nZIndex = 2;
     965             :     }
     966          25 :     nMinTokens = 1 + std::max(std::max(nXIndex, nYIndex), nZIndex);
     967             : 
     968             :     /* -------------------------------------------------------------------- */
     969             :     /*      Parse data lines                                                */
     970             :     /* -------------------------------------------------------------------- */
     971             : 
     972          25 :     GIntBig nLineNum = 0;
     973          25 :     GIntBig nDataLineNum = 0;
     974          25 :     double dfX = 0.0;
     975          25 :     double dfY = 0.0;
     976          25 :     double dfZ = 0.0;
     977          25 :     double dfMinX = 0.0;
     978          25 :     double dfMinY = 0.0;
     979          25 :     double dfMaxX = 0.0;
     980          25 :     double dfMaxY = 0.0;
     981          25 :     double dfMinZ = 0.0;
     982          25 :     double dfMaxZ = 0.0;
     983          25 :     double dfLastX = 0.0;
     984          25 :     double dfLastY = 0.0;
     985          50 :     std::vector<double> adfStepX;
     986          50 :     std::vector<double> adfStepY;
     987          25 :     GDALDataType eDT = GDT_Byte;
     988          25 :     bool bSameNumberOfValuesPerLine = true;
     989          25 :     char chDecimalSep = '\0';
     990          25 :     int nStepYSign = 0;
     991          25 :     bool bColOrganization = false;
     992             : 
     993             :     const char *pszLine;
     994          25 :     GIntBig nCountStepX = 0;
     995          25 :     GIntBig nCountStepY = 0;
     996       41785 :     while ((pszLine = CPLReadLine2L(fp, 100, nullptr)) != nullptr)
     997             :     {
     998       41760 :         nLineNum++;
     999             : 
    1000       41760 :         const char *pszPtr = pszLine;
    1001             :         char ch;
    1002       41760 :         int nCol = 0;
    1003       41760 :         bool bLastWasSep = true;
    1004       41760 :         if (chDecimalSep == '\0')
    1005             :         {
    1006         872 :             int nCountComma = 0;
    1007         872 :             int nCountFieldSep = 0;
    1008       15582 :             while ((ch = *pszPtr) != '\0')
    1009             :             {
    1010       14722 :                 if (ch == '.')
    1011             :                 {
    1012          12 :                     chDecimalSep = '.';
    1013          12 :                     break;
    1014             :                 }
    1015       14710 :                 else if (ch == ',')
    1016             :                 {
    1017           3 :                     nCountComma++;
    1018           3 :                     bLastWasSep = false;
    1019             :                 }
    1020       14707 :                 else if (ch == ' ')
    1021             :                 {
    1022        1701 :                     if (!bLastWasSep)
    1023        1698 :                         nCountFieldSep++;
    1024        1701 :                     bLastWasSep = true;
    1025             :                 }
    1026       13006 :                 else if (ch == '\t' || ch == ';')
    1027             :                 {
    1028           4 :                     nCountFieldSep++;
    1029           4 :                     bLastWasSep = true;
    1030             :                 }
    1031             :                 else
    1032       13002 :                     bLastWasSep = false;
    1033       14710 :                 pszPtr++;
    1034             :             }
    1035         872 :             if (chDecimalSep == '\0')
    1036             :             {
    1037             :                 /* 1,2,3 */
    1038         860 :                 if (nCountComma >= 2 && nCountFieldSep == 0)
    1039           1 :                     chDecimalSep = '.';
    1040             :                 /* 23,5;33;45 */
    1041         859 :                 else if (nCountComma > 0 && nCountFieldSep > 0)
    1042           1 :                     chDecimalSep = ',';
    1043             :             }
    1044         872 :             pszPtr = pszLine;
    1045         872 :             bLastWasSep = true;
    1046             :         }
    1047             : 
    1048       41760 :         char chLocalDecimalSep = chDecimalSep ? chDecimalSep : '.';
    1049       41760 :         int nUsefulColsFound = 0;
    1050     1354170 :         while ((ch = *pszPtr) != '\0')
    1051             :         {
    1052     1312410 :             if (ch == ' ')
    1053             :             {
    1054        2784 :                 if (!bLastWasSep)
    1055        2664 :                     nCol++;
    1056        2784 :                 bLastWasSep = true;
    1057             :             }
    1058     1309620 :             else if ((ch == ',' && chLocalDecimalSep != ',') || ch == '\t' ||
    1059             :                      ch == ';')
    1060             :             {
    1061       80826 :                 nCol++;
    1062       80826 :                 bLastWasSep = true;
    1063             :             }
    1064             :             else
    1065             :             {
    1066     1228800 :                 if (bLastWasSep)
    1067             :                 {
    1068      125235 :                     if (nCol == nXIndex)
    1069             :                     {
    1070       41745 :                         nUsefulColsFound++;
    1071       41745 :                         dfX = CPLAtofDelim(pszPtr, chLocalDecimalSep);
    1072             :                     }
    1073       83490 :                     else if (nCol == nYIndex)
    1074             :                     {
    1075       41745 :                         nUsefulColsFound++;
    1076       41745 :                         dfY = CPLAtofDelim(pszPtr, chLocalDecimalSep);
    1077             :                     }
    1078       41745 :                     else if (nCol == nZIndex)
    1079             :                     {
    1080       41745 :                         nUsefulColsFound++;
    1081       41745 :                         dfZ = CPLAtofDelim(pszPtr, chLocalDecimalSep);
    1082       41745 :                         if (nDataLineNum == 0)
    1083             :                         {
    1084          25 :                             dfMinZ = dfZ;
    1085          25 :                             dfMaxZ = dfZ;
    1086             :                         }
    1087       41720 :                         else if (dfZ < dfMinZ)
    1088             :                         {
    1089          22 :                             dfMinZ = dfZ;
    1090             :                         }
    1091       41698 :                         else if (dfZ > dfMaxZ)
    1092             :                         {
    1093         306 :                             dfMaxZ = dfZ;
    1094             :                         }
    1095             : 
    1096       41745 :                         if (dfZ < INT_MIN || dfZ > INT_MAX)
    1097             :                         {
    1098           0 :                             eDT = GDT_Float32;
    1099             :                         }
    1100             :                         else
    1101             :                         {
    1102       41745 :                             int nZ = static_cast<int>(dfZ);
    1103       41745 :                             if (static_cast<double>(nZ) != dfZ)
    1104             :                             {
    1105       28327 :                                 eDT = GDT_Float32;
    1106             :                             }
    1107       13418 :                             else if ((eDT == GDT_Byte || eDT == GDT_Int16)
    1108             :                                      // cppcheck-suppress
    1109             :                                      // knownConditionTrueFalse
    1110        5588 :                                      && (nZ < 0 || nZ > 255))
    1111             :                             {
    1112          10 :                                 if (nZ < -32768 || nZ > 32767)
    1113           0 :                                     eDT = GDT_Int32;
    1114             :                                 else
    1115          10 :                                     eDT = GDT_Int16;
    1116             :                             }
    1117             :                         }
    1118             :                     }
    1119             :                 }
    1120     1228800 :                 bLastWasSep = false;
    1121             :             }
    1122     1312410 :             pszPtr++;
    1123             :         }
    1124             :         /* skip empty lines */
    1125       41760 :         if (bLastWasSep && nCol == 0)
    1126             :         {
    1127          15 :             continue;
    1128             :         }
    1129       41745 :         nDataLineNum++;
    1130       41745 :         nCol++;
    1131       41745 :         if (nCol < nMinTokens)
    1132             :         {
    1133           0 :             CPLError(CE_Failure, CPLE_AppDefined,
    1134             :                      "At line " CPL_FRMT_GIB
    1135             :                      ", found %d tokens. Expected %d at least",
    1136             :                      nLineNum, nCol, nMinTokens);
    1137           0 :             VSIFCloseL(fp);
    1138           0 :             return nullptr;
    1139             :         }
    1140       41745 :         if (nUsefulColsFound != 3)
    1141             :         {
    1142           0 :             CPLError(CE_Failure, CPLE_AppDefined,
    1143             :                      "At line " CPL_FRMT_GIB
    1144             :                      ", did not find X, Y and/or Z values",
    1145             :                      nLineNum);
    1146           0 :             VSIFCloseL(fp);
    1147           0 :             return nullptr;
    1148             :         }
    1149             : 
    1150       41745 :         if (nDataLineNum == 1)
    1151             :         {
    1152          25 :             dfMinX = dfX;
    1153          25 :             dfMaxX = dfX;
    1154          25 :             dfMinY = dfY;
    1155          25 :             dfMaxY = dfY;
    1156             :         }
    1157       41720 :         else if (nDataLineNum == 2 && dfX == dfLastX)
    1158             :         {
    1159             :             // Detect datasets organized by columns
    1160           7 :             if (dfY == dfLastY)
    1161             :             {
    1162           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
    1163             :                          "Ungridded dataset: At line " CPL_FRMT_GIB
    1164             :                          ", Failed to detect grid layout.",
    1165             :                          nLineNum);
    1166           0 :                 VSIFCloseL(fp);
    1167           0 :                 return nullptr;
    1168             :             }
    1169             : 
    1170           7 :             bColOrganization = true;
    1171           7 :             const double dfStepY = dfY - dfLastY;
    1172           7 :             adfStepY.push_back(fabs(dfStepY));
    1173           7 :             nStepYSign = dfStepY > 0 ? 1 : -1;
    1174             :         }
    1175       41713 :         else if (bColOrganization)
    1176             :         {
    1177          29 :             const double dfStepX = dfX - dfLastX;
    1178          29 :             if (dfStepX == 0)
    1179             :             {
    1180          19 :                 const double dfStepY = dfY - dfLastY;
    1181          19 :                 const double dfExpectedStepY = adfStepY.back() * nStepYSign;
    1182          19 :                 if (fabs((dfStepY - dfExpectedStepY) / dfExpectedStepY) >
    1183             :                     RELATIVE_ERROR)
    1184             :                 {
    1185           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
    1186             :                              "Ungridded dataset: At line " CPL_FRMT_GIB
    1187             :                              ", Y spacing was %f. Expected %f",
    1188             :                              nLineNum, dfStepY, dfExpectedStepY);
    1189           0 :                     VSIFCloseL(fp);
    1190           0 :                     return nullptr;
    1191             :                 }
    1192             :             }
    1193          10 :             else if (dfStepX > 0)
    1194             :             {
    1195           8 :                 if (adfStepX.empty())
    1196             :                 {
    1197           5 :                     adfStepX.push_back(dfStepX);
    1198             :                 }
    1199           3 :                 else if (fabs((dfStepX - adfStepX.back()) / adfStepX.back()) >
    1200             :                          RELATIVE_ERROR)
    1201             :                 {
    1202           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
    1203             :                              "Ungridded dataset: At line " CPL_FRMT_GIB
    1204             :                              ", X spacing was %f. Expected %f",
    1205           0 :                              nLineNum, dfStepX, adfStepX.back());
    1206           0 :                     VSIFCloseL(fp);
    1207           0 :                     return nullptr;
    1208             :                 }
    1209             :             }
    1210           2 :             else if (nDataLineNum == 3)
    1211             :             {
    1212           1 :                 const double dfStepY = dfY - dfLastY;
    1213           1 :                 const double dfLastSignedStepY = nStepYSign * adfStepY.back();
    1214           1 :                 if (dfStepY * dfLastSignedStepY > 0 &&
    1215           1 :                     (fabs(dfStepY - dfLastSignedStepY) <=
    1216           1 :                          RELATIVE_ERROR * fabs(dfLastSignedStepY)
    1217             : #ifdef multiple_of_step_y_not_yet_supported
    1218             :                      ||
    1219             :                      (fabs(dfStepY) > fabs(dfLastSignedStepY) &&
    1220             :                       fabs(std::round(dfStepY / dfLastSignedStepY) -
    1221             :                            (dfStepY / dfLastSignedStepY)) <= RELATIVE_ERROR) ||
    1222             :                      (fabs(dfLastSignedStepY) > fabs(dfStepY) &&
    1223             :                       fabs(std::round(dfLastSignedStepY / dfStepY) -
    1224             :                            (dfLastSignedStepY / dfStepY)) <= RELATIVE_ERROR)
    1225             : #endif
    1226             :                          ))
    1227             :                 {
    1228             :                     // Assume it is a file starting with something like:
    1229             :                     // 371999.50 5806917.50 41.21
    1230             :                     // 371999.50 5806918.50 51.99
    1231             :                     // 371998.50 5806919.50 53.50
    1232             :                     // 371999.50 5806919.50 53.68
    1233           1 :                     adfStepX.push_back(dfLastX - dfX);
    1234           1 :                     bColOrganization = false;
    1235             :                 }
    1236             :                 else
    1237             :                 {
    1238           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
    1239             :                              "Ungridded dataset: At line " CPL_FRMT_GIB
    1240             :                              ", X spacing was %f. Expected >0 value",
    1241             :                              nLineNum, dfX - dfLastX);
    1242           0 :                     VSIFCloseL(fp);
    1243           0 :                     return nullptr;
    1244             :                 }
    1245             :             }
    1246           1 :             else if (!adfStepX.empty() &&
    1247           0 :                      fabs(std::round(-dfStepX / adfStepX[0]) -
    1248           0 :                           (-dfStepX / adfStepX[0])) <= RELATIVE_ERROR)
    1249             :             {
    1250           0 :                 bColOrganization = false;
    1251             :             }
    1252           1 :             else if (adfStepX.empty())
    1253             :             {
    1254           1 :                 adfStepX.push_back(fabs(dfStepX));
    1255           1 :                 bColOrganization = false;
    1256             :             }
    1257             :             else
    1258             :             {
    1259           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
    1260             :                          "Ungridded dataset: At line " CPL_FRMT_GIB
    1261             :                          ", X spacing was %f. Expected a multiple of %f",
    1262           0 :                          nLineNum, dfStepX, adfStepX[0]);
    1263           0 :                 VSIFCloseL(fp);
    1264           0 :                 return nullptr;
    1265             :             }
    1266             :         }
    1267             :         else
    1268             :         {
    1269       41684 :             double dfStepY = dfY - dfLastY;
    1270       41684 :             if (dfStepY == 0.0)
    1271             :             {
    1272       41398 :                 const double dfStepX = dfX - dfLastX;
    1273       41398 :                 if (dfStepX <= 0)
    1274             :                 {
    1275           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
    1276             :                              "Ungridded dataset: At line " CPL_FRMT_GIB
    1277             :                              ", X spacing was %f. Expected >0 value",
    1278             :                              nLineNum, dfStepX);
    1279           0 :                     VSIFCloseL(fp);
    1280           0 :                     return nullptr;
    1281             :                 }
    1282       41398 :                 if (std::find(adfStepX.begin(), adfStepX.end(), dfStepX) ==
    1283       82796 :                     adfStepX.end())
    1284             :                 {
    1285          40 :                     bool bAddNewValue = true;
    1286          40 :                     std::vector<double>::iterator oIter = adfStepX.begin();
    1287          40 :                     std::vector<double> adfStepXNew;
    1288          41 :                     while (oIter != adfStepX.end())
    1289             :                     {
    1290          22 :                         if (fabs((dfStepX - *oIter) / dfStepX) < RELATIVE_ERROR)
    1291             :                         {
    1292          20 :                             double dfNewVal = *oIter;
    1293          20 :                             if (nCountStepX > 0)
    1294             :                             {
    1295             :                                 // Update mean step
    1296             :                                 /* n * mean(n) = (n-1) * mean(n-1) + val(n)
    1297             :                                 mean(n) = mean(n-1) + (val(n) - mean(n-1)) / n
    1298             :                               */
    1299          20 :                                 nCountStepX++;
    1300          20 :                                 dfNewVal += (dfStepX - *oIter) / nCountStepX;
    1301             :                             }
    1302             : 
    1303          20 :                             adfStepXNew.push_back(dfNewVal);
    1304          20 :                             bAddNewValue = false;
    1305          20 :                             break;
    1306             :                         }
    1307           3 :                         else if (dfStepX < *oIter &&
    1308           1 :                                  fabs(*oIter -
    1309           1 :                                       static_cast<int>(*oIter / dfStepX + 0.5) *
    1310           1 :                                           dfStepX) /
    1311           1 :                                          dfStepX <
    1312             :                                      RELATIVE_ERROR)
    1313             :                         {
    1314           1 :                             nCountStepX = -1;  // disable update of mean
    1315           1 :                             ++oIter;
    1316             :                         }
    1317           2 :                         else if (dfStepX > *oIter &&
    1318           2 :                                  fabs(dfStepX -
    1319           1 :                                       static_cast<int>(dfStepX / *oIter + 0.5) *
    1320           1 :                                           (*oIter)) /
    1321           1 :                                          dfStepX <
    1322             :                                      RELATIVE_ERROR)
    1323             :                         {
    1324           1 :                             nCountStepX = -1;  // disable update of mean
    1325           1 :                             bAddNewValue = false;
    1326           1 :                             adfStepXNew.push_back(*oIter);
    1327           1 :                             break;
    1328             :                         }
    1329             :                         else
    1330             :                         {
    1331           0 :                             adfStepXNew.push_back(*oIter);
    1332           0 :                             ++oIter;
    1333             :                         }
    1334             :                     }
    1335          40 :                     adfStepX = std::move(adfStepXNew);
    1336          40 :                     if (bAddNewValue)
    1337             :                     {
    1338          19 :                         CPLDebug("XYZ", "New stepX=%.15f", dfStepX);
    1339          19 :                         adfStepX.push_back(dfStepX);
    1340          19 :                         if (adfStepX.size() == 1 && nCountStepX == 0)
    1341             :                         {
    1342          18 :                             nCountStepX++;
    1343             :                         }
    1344           1 :                         else if (adfStepX.size() == 2)
    1345             :                         {
    1346           0 :                             nCountStepX = -1;  // disable update of mean
    1347             :                         }
    1348           1 :                         else if (adfStepX.size() == 10)
    1349             :                         {
    1350           0 :                             CPLError(
    1351             :                                 CE_Failure, CPLE_AppDefined,
    1352             :                                 "Ungridded dataset: too many stepX values");
    1353           0 :                             VSIFCloseL(fp);
    1354           0 :                             return nullptr;
    1355             :                         }
    1356             :                     }
    1357             :                 }
    1358             :             }
    1359             :             else
    1360             :             {
    1361         286 :                 int bNewStepYSign = (dfStepY < 0.0) ? -1 : 1;
    1362         286 :                 if (nStepYSign == 0)
    1363          18 :                     nStepYSign = bNewStepYSign;
    1364         268 :                 else if (nStepYSign != bNewStepYSign)
    1365             :                 {
    1366           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
    1367             :                              "Ungridded dataset: At line " CPL_FRMT_GIB
    1368             :                              ", change of Y direction",
    1369             :                              nLineNum);
    1370           0 :                     VSIFCloseL(fp);
    1371           0 :                     return nullptr;
    1372             :                 }
    1373         286 :                 if (bNewStepYSign < 0)
    1374         263 :                     dfStepY = -dfStepY;
    1375         286 :                 nCountStepY++;
    1376         286 :                 if (adfStepY.empty())
    1377             :                 {
    1378          18 :                     adfStepY.push_back(dfStepY);
    1379             :                 }
    1380         268 :                 else if (fabs((adfStepY[0] - dfStepY) / dfStepY) >
    1381             :                          RELATIVE_ERROR)
    1382             :                 {
    1383           4 :                     if (dfStepY > adfStepY[0] &&
    1384           2 :                         fabs(std::round(dfStepY / adfStepY[0]) -
    1385           2 :                              (dfStepY / adfStepY[0])) <= RELATIVE_ERROR)
    1386             :                     {
    1387             :                         // The new step is a multiple of the previous one,
    1388             :                         // which means we have a missing line: OK
    1389             :                     }
    1390             :                     else
    1391             :                     {
    1392           0 :                         CPLDebug("XYZ", "New stepY=%.15f prev stepY=%.15f",
    1393           0 :                                  dfStepY, adfStepY[0]);
    1394           0 :                         CPLError(CE_Failure, CPLE_AppDefined,
    1395             :                                  "Ungridded dataset: At line " CPL_FRMT_GIB
    1396             :                                  ", too many stepY values",
    1397             :                                  nLineNum);
    1398           0 :                         VSIFCloseL(fp);
    1399           0 :                         return nullptr;
    1400             :                     }
    1401             :                 }
    1402             :                 else
    1403             :                 {
    1404             :                     // Update mean step
    1405         266 :                     adfStepY[0] += (dfStepY - adfStepY[0]) / nCountStepY;
    1406             :                 }
    1407             :             }
    1408             :         }
    1409             : 
    1410       41745 :         if (dfX < dfMinX)
    1411           4 :             dfMinX = dfX;
    1412       41745 :         if (dfX > dfMaxX)
    1413         285 :             dfMaxX = dfX;
    1414       41745 :         if (dfY < dfMinY)
    1415         267 :             dfMinY = dfY;
    1416       41745 :         if (dfY > dfMaxY)
    1417          33 :             dfMaxY = dfY;
    1418             : 
    1419       41745 :         dfLastX = dfX;
    1420       41745 :         dfLastY = dfY;
    1421             :     }
    1422             : 
    1423          25 :     if (adfStepX.size() != 1 || adfStepX[0] == 0)
    1424             :     {
    1425           0 :         CPLError(CE_Failure, CPLE_AppDefined, "Couldn't determine X spacing");
    1426           0 :         VSIFCloseL(fp);
    1427           0 :         return nullptr;
    1428             :     }
    1429             : 
    1430          25 :     if (adfStepY.size() != 1 || adfStepY[0] == 0)
    1431             :     {
    1432           0 :         CPLError(CE_Failure, CPLE_AppDefined, "Couldn't determine Y spacing");
    1433           0 :         VSIFCloseL(fp);
    1434           0 :         return nullptr;
    1435             :     }
    1436             : 
    1437             :     // Decide for a north-up organization
    1438          25 :     if (bColOrganization)
    1439           5 :         nStepYSign = -1;
    1440             : 
    1441          25 :     const double dfXSize = 1 + ((dfMaxX - dfMinX) / adfStepX[0] + 0.5);
    1442          25 :     const double dfYSize = 1 + ((dfMaxY - dfMinY) / adfStepY[0] + 0.5);
    1443             :     // Test written such as to detect NaN values
    1444          25 :     if (!(dfXSize > 0 && dfXSize < INT_MAX) ||
    1445          25 :         !(dfYSize > 0 && dfYSize < INT_MAX))
    1446             :     {
    1447           0 :         CPLError(CE_Failure, CPLE_AppDefined, "Invalid dimensions");
    1448           0 :         VSIFCloseL(fp);
    1449           0 :         return nullptr;
    1450             :     }
    1451          25 :     const int nXSize = static_cast<int>(dfXSize);
    1452          25 :     const int nYSize = static_cast<int>(dfYSize);
    1453          25 :     const double dfStepX = (dfMaxX - dfMinX) / (nXSize - 1);
    1454          25 :     const double dfStepY = (dfMaxY - dfMinY) / (nYSize - 1) * nStepYSign;
    1455             : 
    1456             : #ifdef DEBUG_VERBOSE
    1457             :     CPLDebug("XYZ", "minx=%f maxx=%f stepx=%f", dfMinX, dfMaxX, dfStepX);
    1458             :     CPLDebug("XYZ", "miny=%f maxy=%f stepy=%f", dfMinY, dfMaxY, dfStepY);
    1459             : #endif
    1460             : 
    1461          25 :     if (nDataLineNum != static_cast<GIntBig>(nXSize) * nYSize)
    1462             :     {
    1463           5 :         if (bColOrganization)
    1464             :         {
    1465           0 :             CPLError(CE_Failure, CPLE_NotSupported,
    1466             :                      "The XYZ driver does not support datasets organized by "
    1467             :                      "columns with missing values");
    1468           0 :             VSIFCloseL(fp);
    1469           0 :             return nullptr;
    1470             :         }
    1471           5 :         bSameNumberOfValuesPerLine = false;
    1472             :     }
    1473          20 :     else if (bColOrganization && nDataLineNum > 100 * 1000 * 1000)
    1474             :     {
    1475           0 :         CPLError(CE_Failure, CPLE_NotSupported,
    1476             :                  "The XYZ driver cannot load datasets organized by "
    1477             :                  "columns with more than 100 million points");
    1478           0 :         VSIFCloseL(fp);
    1479           0 :         return nullptr;
    1480             :     }
    1481             : 
    1482          25 :     const bool bIngestAll = bColOrganization;
    1483          25 :     if (bIngestAll)
    1484             :     {
    1485           5 :         if (eDT == GDT_Int32)
    1486           0 :             eDT = GDT_Float32;
    1487           5 :         else if (eDT == GDT_Byte)
    1488           1 :             eDT = GDT_Int16;
    1489           5 :         CPLAssert(eDT == GDT_Int16 || eDT == GDT_Float32);
    1490             :     }
    1491             : 
    1492          25 :     if (poOpenInfo->eAccess == GA_Update)
    1493             :     {
    1494           0 :         ReportUpdateNotSupportedByDriver("XYZ");
    1495           0 :         VSIFCloseL(fp);
    1496           0 :         return nullptr;
    1497             :     }
    1498             : 
    1499             :     /* -------------------------------------------------------------------- */
    1500             :     /*      Create a corresponding GDALDataset.                             */
    1501             :     /* -------------------------------------------------------------------- */
    1502          25 :     XYZDataset *poDS = new XYZDataset();
    1503          25 :     poDS->fp = fp;
    1504          25 :     poDS->bHasHeaderLine = bHasHeaderLine;
    1505          25 :     poDS->nCommentLineCount = nCommentLineCount;
    1506          25 :     poDS->chDecimalSep = chDecimalSep ? chDecimalSep : '.';
    1507          25 :     poDS->nXIndex = nXIndex;
    1508          25 :     poDS->nYIndex = nYIndex;
    1509          25 :     poDS->nZIndex = nZIndex;
    1510          25 :     poDS->nMinTokens = nMinTokens;
    1511          25 :     poDS->nRasterXSize = nXSize;
    1512          25 :     poDS->nRasterYSize = nYSize;
    1513          25 :     poDS->m_gt[0] = dfMinX - dfStepX / 2;
    1514          25 :     poDS->m_gt[1] = dfStepX;
    1515          25 :     poDS->m_gt[3] = (dfStepY < 0) ? dfMaxY - dfStepY / 2 : dfMinY - dfStepY / 2;
    1516          25 :     poDS->m_gt[5] = dfStepY;
    1517          25 :     poDS->bSameNumberOfValuesPerLine = bSameNumberOfValuesPerLine;
    1518          25 :     poDS->dfMinZ = dfMinZ;
    1519          25 :     poDS->dfMaxZ = dfMaxZ;
    1520          25 :     poDS->bIngestAll = bIngestAll;
    1521             : #ifdef DEBUG_VERBOSE
    1522             :     CPLDebug("XYZ", "bSameNumberOfValuesPerLine = %d",
    1523             :              bSameNumberOfValuesPerLine);
    1524             : #endif
    1525             : 
    1526          25 :     if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
    1527             :     {
    1528           0 :         delete poDS;
    1529           0 :         return nullptr;
    1530             :     }
    1531             : 
    1532             :     /* -------------------------------------------------------------------- */
    1533             :     /*      Create band information objects.                                */
    1534             :     /* -------------------------------------------------------------------- */
    1535          25 :     poDS->nBands = 1;
    1536          50 :     for (int i = 0; i < poDS->nBands; i++)
    1537          25 :         poDS->SetBand(i + 1, new XYZRasterBand(poDS, i + 1, eDT));
    1538             : 
    1539             :     /* -------------------------------------------------------------------- */
    1540             :     /*      Initialize any PAM information.                                 */
    1541             :     /* -------------------------------------------------------------------- */
    1542          25 :     poDS->SetDescription(poOpenInfo->pszFilename);
    1543          25 :     poDS->TryLoadXML();
    1544             : 
    1545             :     /* -------------------------------------------------------------------- */
    1546             :     /*      Support overviews.                                              */
    1547             :     /* -------------------------------------------------------------------- */
    1548          25 :     poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
    1549          25 :     return poDS;
    1550             : }
    1551             : 
    1552             : /************************************************************************/
    1553             : /*                             CreateCopy()                             */
    1554             : /************************************************************************/
    1555             : 
    1556          31 : GDALDataset *XYZDataset::CreateCopy(const char *pszFilename,
    1557             :                                     GDALDataset *poSrcDS, int bStrict,
    1558             :                                     char **papszOptions,
    1559             :                                     GDALProgressFunc pfnProgress,
    1560             :                                     void *pProgressData)
    1561             : {
    1562             :     /* -------------------------------------------------------------------- */
    1563             :     /*      Some some rudimentary checks                                    */
    1564             :     /* -------------------------------------------------------------------- */
    1565          31 :     int nBands = poSrcDS->GetRasterCount();
    1566          31 :     if (nBands == 0)
    1567             :     {
    1568           1 :         CPLError(
    1569             :             CE_Failure, CPLE_NotSupported,
    1570             :             "XYZ driver does not support source dataset with zero band.\n");
    1571           1 :         return nullptr;
    1572             :     }
    1573             : 
    1574          30 :     if (nBands != 1)
    1575             :     {
    1576           4 :         CPLError((bStrict) ? CE_Failure : CE_Warning, CPLE_NotSupported,
    1577             :                  "XYZ driver only uses the first band of the dataset.\n");
    1578           4 :         if (bStrict)
    1579           4 :             return nullptr;
    1580             :     }
    1581             : 
    1582          26 :     if (pfnProgress && !pfnProgress(0.0, nullptr, pProgressData))
    1583           0 :         return nullptr;
    1584             : 
    1585             :     /* -------------------------------------------------------------------- */
    1586             :     /*      Get source dataset info                                         */
    1587             :     /* -------------------------------------------------------------------- */
    1588             : 
    1589          26 :     int nXSize = poSrcDS->GetRasterXSize();
    1590          26 :     int nYSize = poSrcDS->GetRasterYSize();
    1591          26 :     GDALGeoTransform gt;
    1592          26 :     poSrcDS->GetGeoTransform(gt);
    1593          26 :     if (gt[2] != 0 || gt[4] != 0)
    1594             :     {
    1595           0 :         CPLError(CE_Failure, CPLE_NotSupported,
    1596             :                  "XYZ driver does not support CreateCopy() from skewed or "
    1597             :                  "rotated dataset.\n");
    1598           0 :         return nullptr;
    1599             :     }
    1600             : 
    1601          26 :     const GDALDataType eSrcDT = poSrcDS->GetRasterBand(1)->GetRasterDataType();
    1602             :     GDALDataType eReqDT;
    1603          26 :     if (eSrcDT == GDT_Byte || eSrcDT == GDT_Int16 || eSrcDT == GDT_UInt16 ||
    1604             :         eSrcDT == GDT_Int32)
    1605          18 :         eReqDT = GDT_Int32;
    1606             :     else
    1607           8 :         eReqDT = GDT_Float32;
    1608             : 
    1609             :     /* -------------------------------------------------------------------- */
    1610             :     /*      Create target file                                              */
    1611             :     /* -------------------------------------------------------------------- */
    1612             : 
    1613          26 :     VSILFILE *fp = VSIFOpenL(pszFilename, "wb");
    1614          26 :     if (fp == nullptr)
    1615             :     {
    1616           3 :         CPLError(CE_Failure, CPLE_AppDefined, "Cannot create %s", pszFilename);
    1617           3 :         return nullptr;
    1618             :     }
    1619             : 
    1620             :     /* -------------------------------------------------------------------- */
    1621             :     /*      Read creation options                                           */
    1622             :     /* -------------------------------------------------------------------- */
    1623          23 :     const char *pszColSep = CSLFetchNameValue(papszOptions, "COLUMN_SEPARATOR");
    1624          23 :     if (pszColSep == nullptr)
    1625          22 :         pszColSep = " ";
    1626           1 :     else if (EQUAL(pszColSep, "COMMA"))
    1627           0 :         pszColSep = ",";
    1628           1 :     else if (EQUAL(pszColSep, "SPACE"))
    1629           0 :         pszColSep = " ";
    1630           1 :     else if (EQUAL(pszColSep, "SEMICOLON"))
    1631           0 :         pszColSep = ";";
    1632           1 :     else if (EQUAL(pszColSep, "\\t") || EQUAL(pszColSep, "TAB"))
    1633           0 :         pszColSep = "\t";
    1634             : #ifdef DEBUG_VERBOSE
    1635             :     else
    1636             :         CPLDebug("XYZ", "Using raw column separator: '%s' ", pszColSep);
    1637             : #endif
    1638             : 
    1639             :     const char *pszAddHeaderLine =
    1640          23 :         CSLFetchNameValue(papszOptions, "ADD_HEADER_LINE");
    1641          23 :     if (pszAddHeaderLine != nullptr && CPLTestBool(pszAddHeaderLine))
    1642             :     {
    1643           1 :         VSIFPrintfL(fp, "X%sY%sZ\n", pszColSep, pszColSep);
    1644             :     }
    1645             : 
    1646             :     /* -------------------------------------------------------------------- */
    1647             :     /*      Copy imagery                                                    */
    1648             :     /* -------------------------------------------------------------------- */
    1649          23 :     char szFormat[50] = {'\0'};
    1650          23 :     if (eReqDT == GDT_Int32)
    1651          15 :         strcpy(szFormat, "%.17g%c%.17g%c%d\n");
    1652             :     else
    1653           8 :         strcpy(szFormat, "%.17g%c%.17g%c%.17g\n");
    1654             :     const char *pszDecimalPrecision =
    1655          23 :         CSLFetchNameValue(papszOptions, "DECIMAL_PRECISION");
    1656             :     const char *pszSignificantDigits =
    1657          23 :         CSLFetchNameValue(papszOptions, "SIGNIFICANT_DIGITS");
    1658          23 :     bool bIgnoreSigDigits = false;
    1659          23 :     if (pszDecimalPrecision && pszSignificantDigits)
    1660             :     {
    1661           0 :         CPLError(CE_Warning, CPLE_AppDefined,
    1662             :                  "Conflicting precision arguments, using DECIMAL_PRECISION");
    1663           0 :         bIgnoreSigDigits = true;
    1664             :     }
    1665             :     int nPrecision;
    1666          23 :     if (pszSignificantDigits && !bIgnoreSigDigits)
    1667             :     {
    1668           0 :         nPrecision = atoi(pszSignificantDigits);
    1669           0 :         if (nPrecision >= 0)
    1670             :         {
    1671           0 :             if (eReqDT == GDT_Int32)
    1672           0 :                 snprintf(szFormat, sizeof(szFormat), "%%.%dg%%c%%.%dg%%c%%d\n",
    1673             :                          nPrecision, nPrecision);
    1674             :             else
    1675           0 :                 snprintf(szFormat, sizeof(szFormat),
    1676             :                          "%%.%dg%%c%%.%dg%%c%%.%dg\n", nPrecision, nPrecision,
    1677             :                          nPrecision);
    1678             :         }
    1679           0 :         CPLDebug("XYZ", "Setting precision format: %s", szFormat);
    1680             :     }
    1681          23 :     else if (pszDecimalPrecision)
    1682             :     {
    1683           0 :         nPrecision = atoi(pszDecimalPrecision);
    1684           0 :         if (nPrecision >= 0)
    1685             :         {
    1686           0 :             if (eReqDT == GDT_Int32)
    1687           0 :                 snprintf(szFormat, sizeof(szFormat), "%%.%df%%c%%.%df%%c%%d\n",
    1688             :                          nPrecision, nPrecision);
    1689             :             else
    1690           0 :                 snprintf(szFormat, sizeof(szFormat),
    1691             :                          "%%.%df%%c%%.%df%%c%%.%df\n", nPrecision, nPrecision,
    1692             :                          nPrecision);
    1693             :         }
    1694           0 :         CPLDebug("XYZ", "Setting precision format: %s", szFormat);
    1695             :     }
    1696             :     void *pLineBuffer =
    1697          23 :         reinterpret_cast<void *>(CPLMalloc(nXSize * sizeof(int)));
    1698          23 :     CPLErr eErr = CE_None;
    1699         405 :     for (int j = 0; j < nYSize && eErr == CE_None; j++)
    1700             :     {
    1701         382 :         eErr = poSrcDS->GetRasterBand(1)->RasterIO(GF_Read, 0, j, nXSize, 1,
    1702             :                                                    pLineBuffer, nXSize, 1,
    1703             :                                                    eReqDT, 0, 0, nullptr);
    1704         382 :         if (eErr != CE_None)
    1705           0 :             break;
    1706         382 :         const double dfY = gt[3] + (j + 0.5) * gt[5];
    1707         382 :         CPLString osBuf;
    1708       42747 :         for (int i = 0; i < nXSize; i++)
    1709             :         {
    1710       42375 :             const double dfX = gt[0] + (i + 0.5) * gt[1];
    1711             :             char szBuf[256];
    1712       42375 :             if (eReqDT == GDT_Int32)
    1713        1274 :                 CPLsnprintf(szBuf, sizeof(szBuf), szFormat, dfX, pszColSep[0],
    1714        1274 :                             dfY, pszColSep[0],
    1715        1274 :                             reinterpret_cast<int *>(pLineBuffer)[i]);
    1716             :             else
    1717       41101 :                 CPLsnprintf(szBuf, sizeof(szBuf), szFormat, dfX, pszColSep[0],
    1718       41101 :                             dfY, pszColSep[0],
    1719       41101 :                             reinterpret_cast<float *>(pLineBuffer)[i]);
    1720       42375 :             osBuf += szBuf;
    1721       42375 :             if ((i & 1023) == 0 || i == nXSize - 1)
    1722             :             {
    1723         760 :                 if (VSIFWriteL(osBuf, static_cast<int>(osBuf.size()), 1, fp) !=
    1724             :                     1)
    1725             :                 {
    1726          10 :                     eErr = CE_Failure;
    1727          10 :                     CPLError(CE_Failure, CPLE_AppDefined,
    1728             :                              "Write failed, disk full?\n");
    1729          10 :                     break;
    1730             :                 }
    1731         750 :                 osBuf = "";
    1732             :             }
    1733             :         }
    1734         764 :         if (pfnProgress &&
    1735         382 :             !pfnProgress((j + 1) * 1.0 / nYSize, nullptr, pProgressData))
    1736             :         {
    1737           0 :             eErr = CE_Failure;
    1738           0 :             break;
    1739             :         }
    1740             :     }
    1741          23 :     CPLFree(pLineBuffer);
    1742          23 :     VSIFCloseL(fp);
    1743             : 
    1744          23 :     if (eErr != CE_None)
    1745          10 :         return nullptr;
    1746             : 
    1747             :     /* -------------------------------------------------------------------- */
    1748             :     /*      We don't want to call GDALOpen() since it will be expensive,    */
    1749             :     /*      so we "hand prepare" an XYZ dataset referencing our file.       */
    1750             :     /* -------------------------------------------------------------------- */
    1751          13 :     XYZDataset *poXYZ_DS = new XYZDataset();
    1752          13 :     poXYZ_DS->nRasterXSize = nXSize;
    1753          13 :     poXYZ_DS->nRasterYSize = nYSize;
    1754          13 :     poXYZ_DS->nBands = 1;
    1755          13 :     poXYZ_DS->SetBand(1, new XYZRasterBand(poXYZ_DS, 1, eReqDT));
    1756             :     /* If writing to stdout, we can't reopen it --> silence warning */
    1757          13 :     CPLPushErrorHandler(CPLQuietErrorHandler);
    1758          13 :     poXYZ_DS->fp = VSIFOpenL(pszFilename, "rb");
    1759          13 :     CPLPopErrorHandler();
    1760          13 :     poXYZ_DS->m_gt = gt;
    1761          13 :     poXYZ_DS->nXIndex = 0;
    1762          13 :     poXYZ_DS->nYIndex = 1;
    1763          13 :     poXYZ_DS->nZIndex = 2;
    1764          13 :     if (pszAddHeaderLine)
    1765             :     {
    1766           1 :         poXYZ_DS->nDataLineNum = 1;
    1767           1 :         poXYZ_DS->bHasHeaderLine = TRUE;
    1768             :     }
    1769             : 
    1770          13 :     return poXYZ_DS;
    1771             : }
    1772             : 
    1773             : /************************************************************************/
    1774             : /*                          GetGeoTransform()                           */
    1775             : /************************************************************************/
    1776             : 
    1777          19 : CPLErr XYZDataset::GetGeoTransform(GDALGeoTransform &gt) const
    1778             : 
    1779             : {
    1780          19 :     gt = m_gt;
    1781             : 
    1782          19 :     return CE_None;
    1783             : }
    1784             : 
    1785             : /************************************************************************/
    1786             : /*                         GDALRegister_XYZ()                           */
    1787             : /************************************************************************/
    1788             : 
    1789        1935 : void GDALRegister_XYZ()
    1790             : 
    1791             : {
    1792        1935 :     if (GDALGetDriverByName("XYZ") != nullptr)
    1793         282 :         return;
    1794             : 
    1795        1653 :     GDALDriver *poDriver = new GDALDriver();
    1796             : 
    1797        1653 :     poDriver->SetDescription("XYZ");
    1798        1653 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    1799        1653 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "ASCII Gridded XYZ");
    1800        1653 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/xyz.html");
    1801        1653 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "xyz");
    1802        1653 :     poDriver->SetMetadataItem(
    1803             :         GDAL_DMD_CREATIONOPTIONLIST,
    1804             :         "<CreationOptionList>"
    1805             :         "   <Option name='COLUMN_SEPARATOR' type='string' default=' ' "
    1806             :         "description='Separator between fields.'/>"
    1807             :         "   <Option name='ADD_HEADER_LINE' type='boolean' default='false' "
    1808             :         "description='Add an header line with column names.'/>"
    1809             :         "   <Option name='SIGNIFICANT_DIGITS' type='int' description='Number "
    1810             :         "of significant digits when writing floating-point numbers (%g format; "
    1811             :         "default with 18).'/>\n"
    1812             :         "   <Option name='DECIMAL_PRECISION' type='int' description='Number of "
    1813             :         "decimal places when writing floating-point numbers (%f format).'/>\n"
    1814        1653 :         "</CreationOptionList>");
    1815        1653 :     poDriver->SetMetadataItem(
    1816             :         GDAL_DMD_OPENOPTIONLIST,
    1817             :         "<OpenOptionList>"
    1818             :         "   <Option name='COLUMN_ORDER' type='string-select' default='AUTO' "
    1819             :         "description='Specifies the order of the columns. It overrides the "
    1820             :         "header.'>"
    1821             :         "       <Value>AUTO</Value>"
    1822             :         "       <Value>XYZ</Value>"
    1823             :         "       <Value>YXZ</Value>"
    1824             :         "   </Option>"
    1825        1653 :         "</OpenOptionList>");
    1826             : 
    1827        1653 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
    1828             : 
    1829        1653 :     poDriver->pfnOpen = XYZDataset::Open;
    1830        1653 :     poDriver->pfnIdentify = XYZDataset::Identify;
    1831        1653 :     poDriver->pfnCreateCopy = XYZDataset::CreateCopy;
    1832             : 
    1833        1653 :     GetGDALDriverManager()->RegisterDriver(poDriver);
    1834             : }

Generated by: LCOV version 1.14