Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  USGS DEM Driver
       4             :  * Purpose:  All reader for USGS DEM Reader
       5             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       6             :  *
       7             :  * Portions of this module derived from the VTP USGS DEM driver by Ben
       8             :  * Discoe, see http://www.vterrain.org
       9             :  *
      10             :  ******************************************************************************
      11             :  * Copyright (c) 2001, Frank Warmerdam <warmerdam@pobox.com>
      12             :  * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
      13             :  *
      14             :  * Permission is hereby granted, free of charge, to any person obtaining a
      15             :  * copy of this software and associated documentation files (the "Software"),
      16             :  * to deal in the Software without restriction, including without limitation
      17             :  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
      18             :  * and/or sell copies of the Software, and to permit persons to whom the
      19             :  * Software is furnished to do so, subject to the following conditions:
      20             :  *
      21             :  * The above copyright notice and this permission notice shall be included
      22             :  * in all copies or substantial portions of the Software.
      23             :  *
      24             :  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
      25             :  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      26             :  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
      27             :  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      28             :  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
      29             :  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
      30             :  * DEALINGS IN THE SOFTWARE.
      31             :  ****************************************************************************/
      32             : 
      33             : #include "gdal_frmts.h"
      34             : #include "gdal_pam.h"
      35             : #include "ogr_spatialref.h"
      36             : 
      37             : #include <algorithm>
      38             : 
      39             : typedef struct
      40             : {
      41             :     double x;
      42             :     double y;
      43             : } DPoint2;
      44             : 
      45             : constexpr int USGSDEM_NODATA = -32767;
      46             : 
      47             : GDALDataset *USGSDEMCreateCopy(const char *, GDALDataset *, int, char **,
      48             :                                GDALProgressFunc pfnProgress,
      49             :                                void *pProgressData);
      50             : 
      51             : /************************************************************************/
      52             : /*                              ReadInt()                               */
      53             : /************************************************************************/
      54             : 
      55         409 : static int ReadInt(VSILFILE *fp)
      56             : {
      57             :     char c;
      58         409 :     int nRead = 0;
      59             :     char szBuffer[12];
      60         409 :     bool bInProlog = true;
      61             : 
      62             :     while (true)
      63             :     {
      64        3833 :         if (VSIFReadL(&c, 1, 1, fp) != 1)
      65             :         {
      66           4 :             return 0;
      67             :         }
      68        3829 :         if (bInProlog)
      69             :         {
      70        3321 :             if (!isspace(static_cast<unsigned char>(c)))
      71             :             {
      72         405 :                 bInProlog = false;
      73             :             }
      74             :         }
      75        3829 :         if (!bInProlog)
      76             :         {
      77         913 :             if (c != '-' && c != '+' && !(c >= '0' && c <= '9'))
      78             :             {
      79         405 :                 CPL_IGNORE_RET_VAL(VSIFSeekL(fp, VSIFTellL(fp) - 1, SEEK_SET));
      80         405 :                 break;
      81             :             }
      82         508 :             if (nRead < 11)
      83         508 :                 szBuffer[nRead] = c;
      84         508 :             nRead++;
      85             :         }
      86             :     }
      87         405 :     szBuffer[std::min(nRead, 11)] = 0;
      88         405 :     return atoi(szBuffer);
      89             : }
      90             : 
      91             : typedef struct
      92             : {
      93             :     VSILFILE *fp;
      94             :     int max_size;
      95             :     char *buffer;
      96             :     int buffer_size;
      97             :     int cur_index;
      98             : } Buffer;
      99             : 
     100             : /************************************************************************/
     101             : /*                       USGSDEMRefillBuffer()                          */
     102             : /************************************************************************/
     103             : 
     104          53 : static void USGSDEMRefillBuffer(Buffer *psBuffer)
     105             : {
     106          53 :     memmove(psBuffer->buffer, psBuffer->buffer + psBuffer->cur_index,
     107          53 :             psBuffer->buffer_size - psBuffer->cur_index);
     108             : 
     109          53 :     psBuffer->buffer_size -= psBuffer->cur_index;
     110          53 :     psBuffer->buffer_size += static_cast<int>(
     111         106 :         VSIFReadL(psBuffer->buffer + psBuffer->buffer_size, 1,
     112          53 :                   psBuffer->max_size - psBuffer->buffer_size, psBuffer->fp));
     113          53 :     psBuffer->cur_index = 0;
     114          53 : }
     115             : 
     116             : /************************************************************************/
     117             : /*                      USGSDEMGetCurrentFilePos()                      */
     118             : /************************************************************************/
     119             : 
     120         258 : static vsi_l_offset USGSDEMGetCurrentFilePos(const Buffer *psBuffer)
     121             : {
     122         258 :     return VSIFTellL(psBuffer->fp) - psBuffer->buffer_size +
     123         258 :            psBuffer->cur_index;
     124             : }
     125             : 
     126             : /************************************************************************/
     127             : /*                      USGSDEMSetCurrentFilePos()                      */
     128             : /************************************************************************/
     129             : 
     130         258 : static void USGSDEMSetCurrentFilePos(Buffer *psBuffer, vsi_l_offset nNewPos)
     131             : {
     132         258 :     vsi_l_offset nCurPosFP = VSIFTellL(psBuffer->fp);
     133         258 :     if (nNewPos >= nCurPosFP - psBuffer->buffer_size && nNewPos < nCurPosFP)
     134             :     {
     135         242 :         psBuffer->cur_index =
     136         242 :             static_cast<int>(nNewPos - (nCurPosFP - psBuffer->buffer_size));
     137             :     }
     138             :     else
     139             :     {
     140          16 :         CPL_IGNORE_RET_VAL(VSIFSeekL(psBuffer->fp, nNewPos, SEEK_SET));
     141          16 :         psBuffer->buffer_size = 0;
     142          16 :         psBuffer->cur_index = 0;
     143             :     }
     144         258 : }
     145             : 
     146             : /************************************************************************/
     147             : /*               USGSDEMReadIntFromBuffer()                             */
     148             : /************************************************************************/
     149             : 
     150     1080060 : static int USGSDEMReadIntFromBuffer(Buffer *psBuffer, int *pbSuccess = nullptr)
     151             : {
     152             :     char c;
     153             : 
     154             :     while (true)
     155             :     {
     156     1080060 :         if (psBuffer->cur_index >= psBuffer->buffer_size)
     157             :         {
     158          49 :             USGSDEMRefillBuffer(psBuffer);
     159          49 :             if (psBuffer->cur_index >= psBuffer->buffer_size)
     160             :             {
     161           0 :                 if (pbSuccess)
     162           0 :                     *pbSuccess = FALSE;
     163           0 :                 return 0;
     164             :             }
     165             :         }
     166             : 
     167     1080060 :         c = psBuffer->buffer[psBuffer->cur_index];
     168     1080060 :         psBuffer->cur_index++;
     169     1080060 :         if (!isspace(static_cast<unsigned char>(c)))
     170       46830 :             break;
     171             :     }
     172             : 
     173       46830 :     GIntBig nVal = 0;
     174       46830 :     int nSign = 1;
     175       46830 :     if (c == '-')
     176        6374 :         nSign = -1;
     177       40456 :     else if (c == '+')
     178           0 :         nSign = 1;
     179       40456 :     else if (c >= '0' && c <= '9')
     180       40456 :         nVal = c - '0';
     181             :     else
     182             :     {
     183           0 :         if (pbSuccess)
     184           0 :             *pbSuccess = FALSE;
     185           0 :         return 0;
     186             :     }
     187             : 
     188             :     while (true)
     189             :     {
     190      136899 :         if (psBuffer->cur_index >= psBuffer->buffer_size)
     191             :         {
     192           0 :             USGSDEMRefillBuffer(psBuffer);
     193           0 :             if (psBuffer->cur_index >= psBuffer->buffer_size)
     194             :             {
     195           0 :                 if (pbSuccess)
     196           0 :                     *pbSuccess = TRUE;
     197           0 :                 return static_cast<int>(nSign * nVal);
     198             :             }
     199             :         }
     200             : 
     201      136899 :         c = psBuffer->buffer[psBuffer->cur_index];
     202      136899 :         if (c >= '0' && c <= '9')
     203             :         {
     204       90069 :             psBuffer->cur_index++;
     205       90069 :             if (nVal * nSign < INT_MAX && nVal * nSign > INT_MIN)
     206             :             {
     207       90069 :                 nVal = nVal * 10 + (c - '0');
     208       90069 :                 if (nVal * nSign > INT_MAX)
     209             :                 {
     210           0 :                     nVal = INT_MAX;
     211           0 :                     nSign = 1;
     212             :                 }
     213       90069 :                 else if (nVal * nSign < INT_MIN)
     214             :                 {
     215           0 :                     nVal = INT_MIN;
     216           0 :                     nSign = 1;
     217             :                 }
     218             :             }
     219             :         }
     220             :         else
     221             :         {
     222       46830 :             if (pbSuccess)
     223       46830 :                 *pbSuccess = TRUE;
     224       46830 :             return static_cast<int>(nSign * nVal);
     225             :         }
     226             :     }
     227             : }
     228             : 
     229             : /************************************************************************/
     230             : /*                USGSDEMReadDoubleFromBuffer()                         */
     231             : /************************************************************************/
     232             : 
     233        6424 : static double USGSDEMReadDoubleFromBuffer(Buffer *psBuffer, int nCharCount,
     234             :                                           int *pbSuccess = nullptr)
     235             : 
     236             : {
     237        6424 :     if (psBuffer->cur_index + nCharCount > psBuffer->buffer_size)
     238             :     {
     239           4 :         USGSDEMRefillBuffer(psBuffer);
     240           4 :         if (psBuffer->cur_index + nCharCount > psBuffer->buffer_size)
     241             :         {
     242           1 :             if (pbSuccess)
     243           1 :                 *pbSuccess = FALSE;
     244           1 :             return 0;
     245             :         }
     246             :     }
     247             : 
     248        6423 :     char *szPtr = psBuffer->buffer + psBuffer->cur_index;
     249        6423 :     char backupC = szPtr[nCharCount];
     250        6423 :     szPtr[nCharCount] = 0;
     251      160575 :     for (int i = 0; i < nCharCount; i++)
     252             :     {
     253      154152 :         if (szPtr[i] == 'D')
     254        6404 :             szPtr[i] = 'E';
     255             :     }
     256             : 
     257        6423 :     double dfVal = CPLAtof(szPtr);
     258        6423 :     szPtr[nCharCount] = backupC;
     259        6423 :     psBuffer->cur_index += nCharCount;
     260             : 
     261        6423 :     if (pbSuccess)
     262        6423 :         *pbSuccess = TRUE;
     263        6423 :     return dfVal;
     264             : }
     265             : 
     266             : /************************************************************************/
     267             : /*                              DConvert()                              */
     268             : /************************************************************************/
     269             : 
     270         497 : static double DConvert(VSILFILE *fp, int nCharCount)
     271             : 
     272             : {
     273             :     char szBuffer[100];
     274             : 
     275         497 :     CPL_IGNORE_RET_VAL(VSIFReadL(szBuffer, nCharCount, 1, fp));
     276         497 :     szBuffer[nCharCount] = '\0';
     277             : 
     278       12869 :     for (int i = 0; i < nCharCount; i++)
     279             :     {
     280       12372 :         if (szBuffer[i] == 'D')
     281         449 :             szBuffer[i] = 'E';
     282             :     }
     283             : 
     284         994 :     return CPLAtof(szBuffer);
     285             : }
     286             : 
     287             : /************************************************************************/
     288             : /* ==================================================================== */
     289             : /*                              USGSDEMDataset                          */
     290             : /* ==================================================================== */
     291             : /************************************************************************/
     292             : 
     293             : class USGSDEMRasterBand;
     294             : 
     295             : class USGSDEMDataset final : public GDALPamDataset
     296             : {
     297             :     friend class USGSDEMRasterBand;
     298             : 
     299             :     int nDataStartOffset;
     300             :     GDALDataType eNaturalDataFormat;
     301             : 
     302             :     double adfGeoTransform[6];
     303             :     OGRSpatialReference m_oSRS{};
     304             : 
     305             :     double fVRes;
     306             : 
     307             :     const char *pszUnits;
     308             : 
     309             :     int LoadFromFile(VSILFILE *);
     310             : 
     311             :     VSILFILE *fp;
     312             : 
     313             :   public:
     314             :     USGSDEMDataset();
     315             :     ~USGSDEMDataset();
     316             : 
     317             :     static int Identify(GDALOpenInfo *);
     318             :     static GDALDataset *Open(GDALOpenInfo *);
     319             :     CPLErr GetGeoTransform(double *padfTransform) override;
     320             :     const OGRSpatialReference *GetSpatialRef() const override;
     321             : };
     322             : 
     323             : /************************************************************************/
     324             : /* ==================================================================== */
     325             : /*                            USGSDEMRasterBand                         */
     326             : /* ==================================================================== */
     327             : /************************************************************************/
     328             : 
     329             : class USGSDEMRasterBand final : public GDALPamRasterBand
     330             : {
     331             :     friend class USGSDEMDataset;
     332             : 
     333             :   public:
     334             :     explicit USGSDEMRasterBand(USGSDEMDataset *);
     335             : 
     336             :     virtual const char *GetUnitType() override;
     337             :     virtual double GetNoDataValue(int *pbSuccess = nullptr) override;
     338             :     virtual CPLErr IReadBlock(int, int, void *) override;
     339             : };
     340             : 
     341             : /************************************************************************/
     342             : /*                           USGSDEMRasterBand()                            */
     343             : /************************************************************************/
     344             : 
     345          37 : USGSDEMRasterBand::USGSDEMRasterBand(USGSDEMDataset *poDSIn)
     346             : 
     347             : {
     348          37 :     this->poDS = poDSIn;
     349          37 :     this->nBand = 1;
     350             : 
     351          37 :     eDataType = poDSIn->eNaturalDataFormat;
     352             : 
     353          37 :     nBlockXSize = poDSIn->GetRasterXSize();
     354          37 :     nBlockYSize = poDSIn->GetRasterYSize();
     355          37 : }
     356             : 
     357             : /************************************************************************/
     358             : /*                             IReadBlock()                             */
     359             : /************************************************************************/
     360             : 
     361          14 : CPLErr USGSDEMRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff,
     362             :                                      CPL_UNUSED int nBlockYOff, void *pImage)
     363             : 
     364             : {
     365             :     /* int bad = FALSE; */
     366          14 :     USGSDEMDataset *poGDS = reinterpret_cast<USGSDEMDataset *>(poDS);
     367             : 
     368             :     /* -------------------------------------------------------------------- */
     369             :     /*      Initialize image buffer to nodata value.                        */
     370             :     /* -------------------------------------------------------------------- */
     371          14 :     GDALCopyWords(&USGSDEM_NODATA, GDT_Int32, 0, pImage, GetRasterDataType(),
     372             :                   GDALGetDataTypeSizeBytes(GetRasterDataType()),
     373          14 :                   GetXSize() * GetYSize());
     374             : 
     375             :     /* -------------------------------------------------------------------- */
     376             :     /*      Seek to data.                                                   */
     377             :     /* -------------------------------------------------------------------- */
     378          14 :     CPL_IGNORE_RET_VAL(VSIFSeekL(poGDS->fp, poGDS->nDataStartOffset, 0));
     379             : 
     380          14 :     double dfYMin = poGDS->adfGeoTransform[3] +
     381          14 :                     (GetYSize() - 0.5) * poGDS->adfGeoTransform[5];
     382             : 
     383             :     /* -------------------------------------------------------------------- */
     384             :     /*      Read all the profiles into the image buffer.                    */
     385             :     /* -------------------------------------------------------------------- */
     386             : 
     387             :     Buffer sBuffer;
     388          14 :     sBuffer.max_size = 32768;
     389          14 :     sBuffer.buffer = reinterpret_cast<char *>(CPLMalloc(sBuffer.max_size + 1));
     390          14 :     sBuffer.fp = poGDS->fp;
     391          14 :     sBuffer.buffer_size = 0;
     392          14 :     sBuffer.cur_index = 0;
     393             : 
     394        1298 :     for (int i = 0; i < GetXSize(); i++)
     395             :     {
     396             :         int bSuccess;
     397        1285 :         const int nRowNumber = USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
     398        1285 :         if (nRowNumber != 1)
     399           1 :             CPLDebug("USGSDEM", "i = %d, nRowNumber = %d", i, nRowNumber);
     400        1285 :         if (bSuccess)
     401             :         {
     402             :             const int nColNumber =
     403        1285 :                 USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
     404        1285 :             if (nColNumber != i + 1)
     405             :             {
     406           5 :                 CPLDebug("USGSDEM", "i = %d, nColNumber = %d", i, nColNumber);
     407             :             }
     408             :         }
     409             :         const int nCPoints =
     410        1285 :             (bSuccess) ? USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess) : 0;
     411             : #ifdef DEBUG_VERBOSE
     412             :         CPLDebug("USGSDEM", "i = %d, nCPoints = %d", i, nCPoints);
     413             : #endif
     414             : 
     415        1285 :         if (bSuccess)
     416             :         {
     417             :             const int nNumberOfCols =
     418        1285 :                 USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
     419        1285 :             if (nNumberOfCols != 1)
     420             :             {
     421           0 :                 CPLDebug("USGSDEM", "i = %d, nNumberOfCols = %d", i,
     422             :                          nNumberOfCols);
     423             :             }
     424             :         }
     425             : 
     426             :         // x-start
     427        1285 :         if (bSuccess)
     428        1285 :             /* dxStart = */ USGSDEMReadDoubleFromBuffer(&sBuffer, 24,
     429             :                                                         &bSuccess);
     430             : 
     431             :         double dyStart =
     432        1285 :             (bSuccess) ? USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess)
     433        1285 :                        : 0;
     434             :         const double dfElevOffset =
     435        1285 :             (bSuccess) ? USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess)
     436        1285 :                        : 0;
     437             : 
     438             :         // min z value
     439        1285 :         if (bSuccess)
     440        1285 :             /* djunk = */ USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess);
     441             : 
     442             :         // max z value
     443        1285 :         if (bSuccess)
     444        1284 :             /* djunk = */ USGSDEMReadDoubleFromBuffer(&sBuffer, 24, &bSuccess);
     445        1285 :         if (!bSuccess)
     446             :         {
     447           1 :             CPLFree(sBuffer.buffer);
     448           1 :             return CE_Failure;
     449             :         }
     450             : 
     451        1284 :         if (poGDS->m_oSRS.IsGeographic())
     452         246 :             dyStart = dyStart / 3600.0;
     453             : 
     454        1284 :         double dygap = (dfYMin - dyStart) / poGDS->adfGeoTransform[5] + 0.5;
     455        1284 :         if (dygap <= INT_MIN || dygap >= INT_MAX || !CPLIsFinite(dygap))
     456             :         {
     457           0 :             CPLFree(sBuffer.buffer);
     458           0 :             return CE_Failure;
     459             :         }
     460        1284 :         int lygap = static_cast<int>(dygap);
     461        1284 :         if (nCPoints <= 0)
     462           0 :             continue;
     463        1284 :         if (lygap > INT_MAX - nCPoints)
     464           0 :             lygap = INT_MAX - nCPoints;
     465        1284 :         if (lygap < 0 && GetYSize() > INT_MAX + lygap)
     466             :         {
     467           0 :             CPLFree(sBuffer.buffer);
     468           0 :             return CE_Failure;
     469             :         }
     470             : 
     471       42974 :         for (int j = lygap; j < (nCPoints + lygap); j++)
     472             :         {
     473       41690 :             const int iY = GetYSize() - j - 1;
     474             : 
     475       41690 :             const int nElev = USGSDEMReadIntFromBuffer(&sBuffer, &bSuccess);
     476             : #ifdef DEBUG_VERBOSE
     477             :             CPLDebug("USGSDEM", "  j - lygap = %d, nElev = %d", j - lygap,
     478             :                      nElev);
     479             : #endif
     480             : 
     481       41690 :             if (!bSuccess)
     482             :             {
     483           0 :                 CPLFree(sBuffer.buffer);
     484           0 :                 return CE_Failure;
     485             :             }
     486             : 
     487       41690 :             if (iY < 0 || iY >= GetYSize())
     488             :             {
     489             :                 /* bad = TRUE; */
     490             :             }
     491       41690 :             else if (nElev == USGSDEM_NODATA)
     492             :                 /* leave in output buffer as nodata */;
     493             :             else
     494             :             {
     495       36298 :                 const float fComputedElev =
     496       36298 :                     static_cast<float>(nElev * poGDS->fVRes + dfElevOffset);
     497             : 
     498       36298 :                 if (GetRasterDataType() == GDT_Int16)
     499             :                 {
     500       72474 :                     GUInt16 nVal = (fComputedElev < -32768) ? -32768
     501             :                                    : (fComputedElev > 32767)
     502             :                                        ? 32767
     503       36237 :                                        : static_cast<GInt16>(fComputedElev);
     504       36237 :                     reinterpret_cast<GInt16 *>(pImage)[i + iY * GetXSize()] =
     505       36237 :                         nVal;
     506             :                 }
     507             :                 else
     508             :                 {
     509          61 :                     reinterpret_cast<float *>(pImage)[i + iY * GetXSize()] =
     510             :                         fComputedElev;
     511             :                 }
     512             :             }
     513             :         }
     514             : 
     515        1284 :         if (poGDS->nDataStartOffset == 1024)
     516             :         {
     517             :             // Seek to the next 1024 byte boundary.
     518             :             // Some files have 'junk' profile values after the valid/declared
     519             :             // ones
     520         258 :             vsi_l_offset nCurPos = USGSDEMGetCurrentFilePos(&sBuffer);
     521         258 :             vsi_l_offset nNewPos = (nCurPos + 1023) / 1024 * 1024;
     522         258 :             if (nNewPos > nCurPos)
     523             :             {
     524         258 :                 USGSDEMSetCurrentFilePos(&sBuffer, nNewPos);
     525             :             }
     526             :         }
     527             :     }
     528          13 :     CPLFree(sBuffer.buffer);
     529             : 
     530          13 :     return CE_None;
     531             : }
     532             : 
     533             : /************************************************************************/
     534             : /*                           GetNoDataValue()                           */
     535             : /************************************************************************/
     536             : 
     537          16 : double USGSDEMRasterBand::GetNoDataValue(int *pbSuccess)
     538             : 
     539             : {
     540          16 :     if (pbSuccess != nullptr)
     541          14 :         *pbSuccess = TRUE;
     542             : 
     543          16 :     return USGSDEM_NODATA;
     544             : }
     545             : 
     546             : /************************************************************************/
     547             : /*                            GetUnitType()                             */
     548             : /************************************************************************/
     549           9 : const char *USGSDEMRasterBand::GetUnitType()
     550             : {
     551           9 :     USGSDEMDataset *poGDS = reinterpret_cast<USGSDEMDataset *>(poDS);
     552             : 
     553           9 :     return poGDS->pszUnits;
     554             : }
     555             : 
     556             : /************************************************************************/
     557             : /* ==================================================================== */
     558             : /*                              USGSDEMDataset                          */
     559             : /* ==================================================================== */
     560             : /************************************************************************/
     561             : 
     562             : /************************************************************************/
     563             : /*                           USGSDEMDataset()                           */
     564             : /************************************************************************/
     565             : 
     566          37 : USGSDEMDataset::USGSDEMDataset()
     567             :     : nDataStartOffset(0), eNaturalDataFormat(GDT_Unknown), fVRes(0.0),
     568          37 :       pszUnits(nullptr), fp(nullptr)
     569             : {
     570          37 :     m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     571          37 :     memset(adfGeoTransform, 0, sizeof(adfGeoTransform));
     572          37 : }
     573             : 
     574             : /************************************************************************/
     575             : /*                            ~USGSDEMDataset()                         */
     576             : /************************************************************************/
     577             : 
     578          74 : USGSDEMDataset::~USGSDEMDataset()
     579             : 
     580             : {
     581          37 :     FlushCache(true);
     582             : 
     583          37 :     if (fp != nullptr)
     584          37 :         CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
     585          74 : }
     586             : 
     587             : /************************************************************************/
     588             : /*                            LoadFromFile()                            */
     589             : /*                                                                      */
     590             : /*      If the data from DEM is in meters, then values are stored as    */
     591             : /*      shorts. If DEM data is in feet, then height data will be        */
     592             : /*      stored in float, to preserve the precision of the original      */
     593             : /*      data. returns true if the file was successfully opened and      */
     594             : /*      read.                                                           */
     595             : /************************************************************************/
     596             : 
     597          37 : int USGSDEMDataset::LoadFromFile(VSILFILE *InDem)
     598             : {
     599             :     // check for version of DEM format
     600          37 :     CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 864, 0));
     601             : 
     602             :     // Read DEM into matrix
     603          37 :     const int nRow = ReadInt(InDem);
     604          37 :     const int nColumn = ReadInt(InDem);
     605             :     const bool bNewFormat =
     606          37 :         VSIFTellL(InDem) >= 1024 || nRow != 1 || nColumn != 1;
     607          37 :     if (bNewFormat)
     608             :     {
     609          37 :         CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 1024, 0));  // New Format
     610          37 :         int i = ReadInt(InDem);
     611          37 :         int j = ReadInt(InDem);
     612          37 :         if (i != 1 || (j != 1 && j != 0))  // File OK?
     613             :         {
     614           4 :             CPL_IGNORE_RET_VAL(
     615           4 :                 VSIFSeekL(InDem, 893, 0));  // Undocumented Format (39109h1.dem)
     616           4 :             i = ReadInt(InDem);
     617           4 :             j = ReadInt(InDem);
     618           4 :             if (i != 1 || j != 1)  // File OK?
     619             :             {
     620           2 :                 CPL_IGNORE_RET_VAL(VSIFSeekL(
     621             :                     InDem, 918, 0));  // Latest iteration of the A record, such
     622             :                                       // as in fema06-140cm_2995441b.dem
     623           2 :                 i = ReadInt(InDem);
     624           2 :                 j = ReadInt(InDem);
     625           2 :                 if (i != 1 || j != 1)  // File OK?
     626             :                 {
     627           0 :                     CPLError(CE_Failure, CPLE_AppDefined,
     628             :                              "Does not appear to be a USGS DEM file.");
     629           0 :                     return FALSE;
     630             :                 }
     631             :                 else
     632           2 :                     nDataStartOffset = 918;
     633             :             }
     634             :             else
     635           2 :                 nDataStartOffset = 893;
     636             :         }
     637             :         else
     638             :         {
     639          33 :             nDataStartOffset = 1024;
     640             : 
     641             :             // Some files use 1025 byte records ending with a newline character.
     642             :             // See https://github.com/OSGeo/gdal/issues/5007
     643          33 :             CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 1024, 0));
     644             :             char c;
     645          66 :             if (VSIFReadL(&c, 1, 1, InDem) == 1 && c == '\n' &&
     646           2 :                 VSIFSeekL(InDem, 1024 + 1024 + 1, 0) == 0 &&
     647          66 :                 VSIFReadL(&c, 1, 1, InDem) == 1 && c == '\n')
     648             :             {
     649           2 :                 nDataStartOffset = 1025;
     650             :             }
     651             :         }
     652             :     }
     653             :     else
     654           0 :         nDataStartOffset = 864;
     655             : 
     656          37 :     CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 156, 0));
     657          37 :     const int nCoordSystem = ReadInt(InDem);
     658          37 :     const int iUTMZone = ReadInt(InDem);
     659             : 
     660          37 :     CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 528, 0));
     661          37 :     const int nGUnit = ReadInt(InDem);
     662          37 :     const int nVUnit = ReadInt(InDem);
     663             : 
     664             :     // Vertical Units in meters
     665          37 :     if (nVUnit == 1)
     666           0 :         pszUnits = "ft";
     667             :     else
     668          37 :         pszUnits = "m";
     669             : 
     670          37 :     CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 816, 0));
     671          37 :     const double dxdelta = DConvert(InDem, 12);
     672          37 :     const double dydelta = DConvert(InDem, 12);
     673          37 :     if (dydelta == 0)
     674           0 :         return FALSE;
     675          37 :     fVRes = DConvert(InDem, 12);
     676             : 
     677             :     /* -------------------------------------------------------------------- */
     678             :     /*      Should we treat this as floating point, or GInt16.              */
     679             :     /* -------------------------------------------------------------------- */
     680          37 :     if (nVUnit == 1 || fVRes < 1.0)
     681           4 :         eNaturalDataFormat = GDT_Float32;
     682             :     else
     683          33 :         eNaturalDataFormat = GDT_Int16;
     684             : 
     685             :     /* -------------------------------------------------------------------- */
     686             :     /*      Read four corner coordinates.                                   */
     687             :     /* -------------------------------------------------------------------- */
     688          37 :     CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 546, 0));
     689             :     DPoint2 corners[4];  // SW, NW, NE, SE
     690         185 :     for (int i = 0; i < 4; i++)
     691             :     {
     692         148 :         corners[i].x = DConvert(InDem, 24);
     693         148 :         corners[i].y = DConvert(InDem, 24);
     694             :     }
     695             : 
     696             :     // find absolute extents of raw vales
     697             :     DPoint2 extent_min, extent_max;
     698          37 :     extent_min.x = std::min(corners[0].x, corners[1].x);
     699          37 :     extent_max.x = std::max(corners[2].x, corners[3].x);
     700          37 :     extent_min.y = std::min(corners[0].y, corners[3].y);
     701          37 :     extent_max.y = std::max(corners[1].y, corners[2].y);
     702             : 
     703          37 :     /* dElevMin = */ DConvert(InDem, 48);
     704          37 :     /* dElevMax = */ DConvert(InDem, 48);
     705             : 
     706          37 :     CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 858, 0));
     707          37 :     const int nProfiles = ReadInt(InDem);
     708             : 
     709             :     /* -------------------------------------------------------------------- */
     710             :     /*      Collect the spatial reference system.                           */
     711             :     /* -------------------------------------------------------------------- */
     712          74 :     OGRSpatialReference sr;
     713          37 :     sr.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     714          37 :     bool bNAD83 = true;
     715             : 
     716             :     // OLD format header ends at byte 864
     717          37 :     if (bNewFormat)
     718             :     {
     719             :         // year of data compilation
     720          37 :         CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 876, 0));
     721             :         char szDateBuffer[5];
     722          37 :         CPL_IGNORE_RET_VAL(VSIFReadL(szDateBuffer, 4, 1, InDem));
     723             :         /* szDateBuffer[4] = 0; */
     724             : 
     725             :         // Horizontal datum
     726             :         // 1=North American Datum 1927 (NAD 27)
     727             :         // 2=World Geodetic System 1972 (WGS 72)
     728             :         // 3=WGS 84
     729             :         // 4=NAD 83
     730             :         // 5=Old Hawaii Datum
     731             :         // 6=Puerto Rico Datum
     732          37 :         CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, 890, 0));
     733             : 
     734             :         char szHorzDatum[3];
     735          37 :         CPL_IGNORE_RET_VAL(VSIFReadL(szHorzDatum, 1, 2, InDem));
     736          37 :         szHorzDatum[2] = '\0';
     737          37 :         const int datum = atoi(szHorzDatum);
     738          37 :         switch (datum)
     739             :         {
     740           4 :             case 1:
     741           4 :                 sr.SetWellKnownGeogCS("NAD27");
     742           4 :                 bNAD83 = false;
     743           4 :                 break;
     744             : 
     745           6 :             case 2:
     746           6 :                 sr.SetWellKnownGeogCS("WGS72");
     747           6 :                 break;
     748             : 
     749          14 :             case 3:
     750          14 :                 sr.SetWellKnownGeogCS("WGS84");
     751          14 :                 break;
     752             : 
     753           3 :             case 4:
     754           3 :                 sr.SetWellKnownGeogCS("NAD83");
     755           3 :                 break;
     756             : 
     757           0 :             case -9:
     758           0 :                 break;
     759             : 
     760          10 :             default:
     761          10 :                 sr.SetWellKnownGeogCS("NAD27");
     762          10 :                 break;
     763             :         }
     764             :     }
     765             :     else
     766             :     {
     767           0 :         sr.SetWellKnownGeogCS("NAD27");
     768           0 :         bNAD83 = false;
     769             :     }
     770             : 
     771          37 :     if (nCoordSystem == 1)  // UTM
     772             :     {
     773          16 :         if (iUTMZone >= -60 && iUTMZone <= 60)
     774             :         {
     775          16 :             sr.SetUTM(abs(iUTMZone), iUTMZone >= 0);
     776          16 :             if (nGUnit == 1)
     777             :             {
     778           0 :                 sr.SetLinearUnitsAndUpdateParameters(
     779             :                     SRS_UL_US_FOOT, CPLAtof(SRS_UL_US_FOOT_CONV));
     780             :                 char szUTMName[128];
     781           0 :                 snprintf(szUTMName, sizeof(szUTMName),
     782             :                          "UTM Zone %d, Northern Hemisphere, us-ft", iUTMZone);
     783           0 :                 sr.SetNode("PROJCS", szUTMName);
     784             :             }
     785             :         }
     786             :     }
     787          21 :     else if (nCoordSystem == 2)  // state plane
     788             :     {
     789           0 :         if (nGUnit == 1)
     790           0 :             sr.SetStatePlane(iUTMZone, bNAD83, "Foot",
     791             :                              CPLAtof(SRS_UL_US_FOOT_CONV));
     792             :         else
     793           0 :             sr.SetStatePlane(iUTMZone, bNAD83);
     794             :     }
     795             : 
     796          37 :     m_oSRS = std::move(sr);
     797             : 
     798             :     /* -------------------------------------------------------------------- */
     799             :     /*      For UTM we use the extents (really the UTM coordinates of       */
     800             :     /*      the lat/long corners of the quad) to determine the size in      */
     801             :     /*      pixels and lines, but we have to make the anchors be modulus    */
     802             :     /*      the pixel size which what really gets used.                     */
     803             :     /* -------------------------------------------------------------------- */
     804          37 :     if (nCoordSystem == 1          // UTM
     805          21 :         || nCoordSystem == 2       // State Plane
     806          21 :         || nCoordSystem == -9999)  // unknown
     807             :     {
     808             :         // expand extents modulus the pixel size.
     809          16 :         extent_min.y = floor(extent_min.y / dydelta) * dydelta;
     810          16 :         extent_max.y = ceil(extent_max.y / dydelta) * dydelta;
     811             : 
     812             :         // Forcibly compute X extents based on first profile and pixelsize.
     813          16 :         CPL_IGNORE_RET_VAL(VSIFSeekL(InDem, nDataStartOffset, 0));
     814          16 :         /* njunk = */ ReadInt(InDem);
     815          16 :         /* njunk = */ ReadInt(InDem);
     816          16 :         /* njunk = */ ReadInt(InDem);
     817          16 :         /* njunk = */ ReadInt(InDem);
     818          16 :         const double dxStart = DConvert(InDem, 24);
     819             : 
     820          16 :         nRasterYSize =
     821          16 :             static_cast<int>((extent_max.y - extent_min.y) / dydelta + 1.5);
     822          16 :         nRasterXSize = nProfiles;
     823             : 
     824          16 :         adfGeoTransform[0] = dxStart - dxdelta / 2.0;
     825          16 :         adfGeoTransform[1] = dxdelta;
     826          16 :         adfGeoTransform[2] = 0.0;
     827          16 :         adfGeoTransform[3] = extent_max.y + dydelta / 2.0;
     828          16 :         adfGeoTransform[4] = 0.0;
     829          16 :         adfGeoTransform[5] = -dydelta;
     830             :     }
     831             :     /* -------------------------------------------------------------------- */
     832             :     /*      Geographic -- use corners directly.                             */
     833             :     /* -------------------------------------------------------------------- */
     834             :     else
     835             :     {
     836          21 :         nRasterYSize =
     837          21 :             static_cast<int>((extent_max.y - extent_min.y) / dydelta + 1.5);
     838          21 :         nRasterXSize = nProfiles;
     839             : 
     840             :         // Translate extents from arc-seconds to decimal degrees.
     841          21 :         adfGeoTransform[0] = (extent_min.x - dxdelta / 2.0) / 3600.0;
     842          21 :         adfGeoTransform[1] = dxdelta / 3600.0;
     843          21 :         adfGeoTransform[2] = 0.0;
     844          21 :         adfGeoTransform[3] = (extent_max.y + dydelta / 2.0) / 3600.0;
     845          21 :         adfGeoTransform[4] = 0.0;
     846          21 :         adfGeoTransform[5] = (-dydelta) / 3600.0;
     847             :     }
     848             : 
     849             :     // IReadBlock() not ready for more than INT_MAX pixels, and that
     850             :     // would behave badly
     851          74 :     if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize) ||
     852          37 :         nRasterXSize > INT_MAX / nRasterYSize)
     853             :     {
     854           0 :         return FALSE;
     855             :     }
     856             : 
     857          37 :     return TRUE;
     858             : }
     859             : 
     860             : /************************************************************************/
     861             : /*                          GetGeoTransform()                           */
     862             : /************************************************************************/
     863             : 
     864          33 : CPLErr USGSDEMDataset::GetGeoTransform(double *padfTransform)
     865             : 
     866             : {
     867          33 :     memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
     868          33 :     return CE_None;
     869             : }
     870             : 
     871             : /************************************************************************/
     872             : /*                          GetSpatialRef()                             */
     873             : /************************************************************************/
     874             : 
     875          32 : const OGRSpatialReference *USGSDEMDataset::GetSpatialRef() const
     876             : 
     877             : {
     878          32 :     return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
     879             : }
     880             : 
     881             : /************************************************************************/
     882             : /*                              Identify()                              */
     883             : /************************************************************************/
     884             : 
     885       48646 : int USGSDEMDataset::Identify(GDALOpenInfo *poOpenInfo)
     886             : 
     887             : {
     888       48646 :     if (poOpenInfo->nHeaderBytes < 200)
     889       45828 :         return FALSE;
     890             : 
     891        2818 :     if (!STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, "     0") &&
     892        2776 :         !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, "     1") &&
     893        2744 :         !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, "     2") &&
     894        2744 :         !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, "     3") &&
     895        2744 :         !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 156, " -9999"))
     896        2744 :         return FALSE;
     897             : 
     898          74 :     if (!STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 150, "     1") &&
     899           6 :         !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 150, "     4"))
     900           0 :         return FALSE;
     901             : 
     902          74 :     return TRUE;
     903             : }
     904             : 
     905             : /************************************************************************/
     906             : /*                                Open()                                */
     907             : /************************************************************************/
     908             : 
     909          37 : GDALDataset *USGSDEMDataset::Open(GDALOpenInfo *poOpenInfo)
     910             : 
     911             : {
     912          37 :     if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
     913           0 :         return nullptr;
     914             : 
     915             :     /* -------------------------------------------------------------------- */
     916             :     /*      Create a corresponding GDALDataset.                             */
     917             :     /* -------------------------------------------------------------------- */
     918          37 :     USGSDEMDataset *poDS = new USGSDEMDataset();
     919             : 
     920          37 :     poDS->fp = poOpenInfo->fpL;
     921          37 :     poOpenInfo->fpL = nullptr;
     922             : 
     923             :     /* -------------------------------------------------------------------- */
     924             :     /*      Read the file.                                                  */
     925             :     /* -------------------------------------------------------------------- */
     926          37 :     if (!poDS->LoadFromFile(poDS->fp))
     927             :     {
     928           0 :         delete poDS;
     929           0 :         return nullptr;
     930             :     }
     931             : 
     932             :     /* -------------------------------------------------------------------- */
     933             :     /*      Confirm the requested access is supported.                      */
     934             :     /* -------------------------------------------------------------------- */
     935          37 :     if (poOpenInfo->eAccess == GA_Update)
     936             :     {
     937           0 :         delete poDS;
     938           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     939             :                  "The USGSDEM driver does not support update access to existing"
     940             :                  " datasets.\n");
     941           0 :         return nullptr;
     942             :     }
     943             : 
     944             :     /* -------------------------------------------------------------------- */
     945             :     /*      Create band information objects.                                */
     946             :     /* -------------------------------------------------------------------- */
     947          37 :     poDS->SetBand(1, new USGSDEMRasterBand(poDS));
     948             : 
     949          37 :     poDS->SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
     950             : 
     951             :     /* -------------------------------------------------------------------- */
     952             :     /*      Initialize any PAM information.                                 */
     953             :     /* -------------------------------------------------------------------- */
     954          37 :     poDS->SetDescription(poOpenInfo->pszFilename);
     955          37 :     poDS->TryLoadXML();
     956             : 
     957             :     /* -------------------------------------------------------------------- */
     958             :     /*      Open overviews.                                                 */
     959             :     /* -------------------------------------------------------------------- */
     960          37 :     poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
     961             : 
     962          37 :     return poDS;
     963             : }
     964             : 
     965             : /************************************************************************/
     966             : /*                        GDALRegister_USGSDEM()                        */
     967             : /************************************************************************/
     968             : 
     969        1520 : void GDALRegister_USGSDEM()
     970             : 
     971             : {
     972        1520 :     if (GDALGetDriverByName("USGSDEM") != nullptr)
     973         301 :         return;
     974             : 
     975        1219 :     GDALDriver *poDriver = new GDALDriver();
     976             : 
     977        1219 :     poDriver->SetDescription("USGSDEM");
     978        1219 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     979        1219 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "dem");
     980        1219 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
     981        1219 :                               "USGS Optional ASCII DEM (and CDED)");
     982        1219 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
     983        1219 :                               "drivers/raster/usgsdem.html");
     984        1219 :     poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Int16");
     985        1219 :     poDriver->SetMetadataItem(
     986             :         GDAL_DMD_CREATIONOPTIONLIST,
     987             :         "<CreationOptionList>"
     988             :         "   <Option name='PRODUCT' type='string-select' description='Specific "
     989             :         "Product Type'>"
     990             :         "       <Value>DEFAULT</Value>"
     991             :         "       <Value>CDED50K</Value>"
     992             :         "   </Option>"
     993             :         "   <Option name='TOPLEFT' type='string' description='Top left product "
     994             :         "corner (i.e. 117d15w,52d30n'/>"
     995             :         "   <Option name='RESAMPLE' type='string-select' "
     996             :         "description='Resampling kernel to use if resampled.'>"
     997             :         "       <Value>Nearest</Value>"
     998             :         "       <Value>Bilinear</Value>"
     999             :         "       <Value>Cubic</Value>"
    1000             :         "       <Value>CubicSpline</Value>"
    1001             :         "   </Option>"
    1002             :         "   <Option name='TEMPLATE' type='string' description='File to default "
    1003             :         "metadata from.'/>"
    1004             :         "   <Option name='DEMLevelCode' type='int' description='DEM Level (1, "
    1005             :         "2 or 3 if set)'/>"
    1006             :         "   <Option name='DataSpecVersion' type='int' description='Data and "
    1007             :         "Specification version/revision (eg. 1020)'/>"
    1008             :         "   <Option name='PRODUCER' type='string' description='Producer Agency "
    1009             :         "(up to 60 characters)'/>"
    1010             :         "   <Option name='OriginCode' type='string' description='Origin code "
    1011             :         "(up to 4 characters, YT for Yukon)'/>"
    1012             :         "   <Option name='ProcessCode' type='string' description='Processing "
    1013             :         "Code (8=ANUDEM, 9=FME, A=TopoGrid)'/>"
    1014             :         "   <Option name='ZRESOLUTION' type='float' description='Scaling "
    1015             :         "factor for elevation values'/>"
    1016             :         "   <Option name='NTS' type='string' description='NTS Mapsheet name, "
    1017             :         "used to derive TOPLEFT.'/>"
    1018             :         "   <Option name='INTERNALNAME' type='string' description='Dataset "
    1019             :         "name written into file header.'/>"
    1020        1219 :         "</CreationOptionList>");
    1021        1219 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
    1022             : 
    1023        1219 :     poDriver->pfnOpen = USGSDEMDataset::Open;
    1024        1219 :     poDriver->pfnCreateCopy = USGSDEMCreateCopy;
    1025        1219 :     poDriver->pfnIdentify = USGSDEMDataset::Identify;
    1026             : 
    1027        1219 :     GetGDALDriverManager()->RegisterDriver(poDriver);
    1028             : }