LCOV - code coverage report
Current view: top level - frmts/srtmhgt - srtmhgtdataset.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 265 298 88.9 %
Date: 2025-01-18 12:42:00 Functions: 16 16 100.0 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  SRTM HGT Driver
       4             :  * Purpose:  SRTM HGT File Read Support.
       5             :  *           http://dds.cr.usgs.gov/srtm/version2_1/Documentation/SRTM_Topo.pdf
       6             :  *           http://www2.jpl.nasa.gov/srtm/faq.html
       7             :  *           http://dds.cr.usgs.gov/srtm/version2_1
       8             :  * Authors:  Michael Mazzella, Even Rouault
       9             :  *
      10             :  ******************************************************************************
      11             :  * Copyright (c) 2005, Frank Warmerdam <warmerdam@pobox.com>
      12             :  * Copyright (c) 2007-2011, Even Rouault <even dot rouault at spatialys.com>
      13             :  *
      14             :  * SPDX-License-Identifier: MIT
      15             :  ****************************************************************************/
      16             : 
      17             : #include "cpl_port.h"
      18             : #include "cpl_string.h"
      19             : #include "gdal_frmts.h"
      20             : #include "gdal_pam.h"
      21             : #include "ogr_spatialref.h"
      22             : 
      23             : #include <cmath>
      24             : 
      25             : constexpr GInt16 SRTMHG_NODATA_VALUE = -32768;
      26             : 
      27             : /************************************************************************/
      28             : /* ==================================================================== */
      29             : /*                              SRTMHGTDataset                          */
      30             : /* ==================================================================== */
      31             : /************************************************************************/
      32             : 
      33             : class SRTMHGTRasterBand;
      34             : 
      35             : class SRTMHGTDataset final : public GDALPamDataset
      36             : {
      37             :     friend class SRTMHGTRasterBand;
      38             : 
      39             :     VSILFILE *fpImage = nullptr;
      40             :     double adfGeoTransform[6];
      41             :     GByte *pabyBuffer = nullptr;
      42             :     OGRSpatialReference m_oSRS{};
      43             : 
      44             :     static GDALPamDataset *OpenPAM(GDALOpenInfo *);
      45             : 
      46             :   public:
      47             :     SRTMHGTDataset();
      48             :     virtual ~SRTMHGTDataset();
      49             : 
      50             :     const OGRSpatialReference *GetSpatialRef() const override;
      51             :     virtual CPLErr GetGeoTransform(double *) override;
      52             : 
      53             :     static int Identify(GDALOpenInfo *poOpenInfo);
      54             :     static GDALDataset *Open(GDALOpenInfo *);
      55             :     static GDALDataset *CreateCopy(const char *pszFilename,
      56             :                                    GDALDataset *poSrcDS, int bStrict,
      57             :                                    char **papszOptions,
      58             :                                    GDALProgressFunc pfnProgress,
      59             :                                    void *pProgressData);
      60             : };
      61             : 
      62             : /************************************************************************/
      63             : /* ==================================================================== */
      64             : /*                            SRTMHGTRasterBand                             */
      65             : /* ==================================================================== */
      66             : /************************************************************************/
      67             : 
      68             : class SRTMHGTRasterBand final : public GDALPamRasterBand
      69             : {
      70             :     friend class SRTMHGTDataset;
      71             : 
      72             :     int bNoDataSet;
      73             :     double dfNoDataValue;
      74             : 
      75             :   public:
      76             :     SRTMHGTRasterBand(SRTMHGTDataset *, int, GDALDataType);
      77             : 
      78             :     virtual CPLErr IReadBlock(int, int, void *) override;
      79             :     virtual CPLErr IWriteBlock(int nBlockXOff, int nBlockYOff,
      80             :                                void *pImage) override;
      81             : 
      82             :     virtual GDALColorInterp GetColorInterpretation() override;
      83             : 
      84             :     virtual double GetNoDataValue(int *pbSuccess = nullptr) override;
      85             : 
      86             :     virtual const char *GetUnitType() override;
      87             : };
      88             : 
      89             : /************************************************************************/
      90             : /*                         SRTMHGTRasterBand()                          */
      91             : /************************************************************************/
      92             : 
      93          24 : SRTMHGTRasterBand::SRTMHGTRasterBand(SRTMHGTDataset *poDSIn, int nBandIn,
      94          24 :                                      GDALDataType eDT)
      95          24 :     : bNoDataSet(TRUE), dfNoDataValue(SRTMHG_NODATA_VALUE)
      96             : {
      97          24 :     poDS = poDSIn;
      98          24 :     nBand = nBandIn;
      99          24 :     eDataType = eDT;
     100          24 :     nBlockXSize = poDSIn->nRasterXSize;
     101          24 :     nBlockYSize = 1;
     102          24 : }
     103             : 
     104             : /************************************************************************/
     105             : /*                             IReadBlock()                             */
     106             : /************************************************************************/
     107             : 
     108       37214 : CPLErr SRTMHGTRasterBand::IReadBlock(int /*nBlockXOff*/, int nBlockYOff,
     109             :                                      void *pImage)
     110             : {
     111       37214 :     SRTMHGTDataset *poGDS = reinterpret_cast<SRTMHGTDataset *>(poDS);
     112             : 
     113             :     /* -------------------------------------------------------------------- */
     114             :     /*      Load the desired data into the working buffer.                  */
     115             :     /* -------------------------------------------------------------------- */
     116       37214 :     const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
     117       37214 :     VSIFSeekL(poGDS->fpImage,
     118       37214 :               static_cast<size_t>(nBlockYOff) * nBlockXSize * nDTSize,
     119             :               SEEK_SET);
     120       37214 :     VSIFReadL((unsigned char *)pImage, nBlockXSize, nDTSize, poGDS->fpImage);
     121             : #ifdef CPL_LSB
     122       37214 :     GDALSwapWords(pImage, nDTSize, nBlockXSize, nDTSize);
     123             : #endif
     124             : 
     125       37214 :     return CE_None;
     126             : }
     127             : 
     128             : /************************************************************************/
     129             : /*                             IWriteBlock()                            */
     130             : /************************************************************************/
     131             : 
     132        1201 : CPLErr SRTMHGTRasterBand::IWriteBlock(int /*nBlockXOff*/, int nBlockYOff,
     133             :                                       void *pImage)
     134             : {
     135        1201 :     SRTMHGTDataset *poGDS = reinterpret_cast<SRTMHGTDataset *>(poDS);
     136             : 
     137        1201 :     if (poGDS->eAccess != GA_Update)
     138           0 :         return CE_Failure;
     139             : 
     140        1201 :     const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
     141        1201 :     VSIFSeekL(poGDS->fpImage,
     142        1201 :               static_cast<size_t>(nBlockYOff) * nBlockXSize * nDTSize,
     143             :               SEEK_SET);
     144             : 
     145             : #ifdef CPL_LSB
     146        1201 :     if (nDTSize > 1)
     147             :     {
     148        1201 :         memcpy(poGDS->pabyBuffer, pImage,
     149        1201 :                static_cast<size_t>(nBlockXSize) * nDTSize);
     150        1201 :         GDALSwapWords(poGDS->pabyBuffer, nDTSize, nBlockXSize, nDTSize);
     151        1201 :         VSIFWriteL(reinterpret_cast<unsigned char *>(poGDS->pabyBuffer),
     152        1201 :                    nBlockXSize, nDTSize, poGDS->fpImage);
     153             :     }
     154             :     else
     155             : #endif
     156             :     {
     157           0 :         VSIFWriteL(reinterpret_cast<unsigned char *>(pImage), nBlockXSize,
     158             :                    nDTSize, poGDS->fpImage);
     159             :     }
     160             : 
     161        1201 :     return CE_None;
     162             : }
     163             : 
     164             : /************************************************************************/
     165             : /*                           GetNoDataValue()                           */
     166             : /************************************************************************/
     167             : 
     168          41 : double SRTMHGTRasterBand::GetNoDataValue(int *pbSuccess)
     169             : 
     170             : {
     171          41 :     if (eDataType == GDT_Byte)
     172           3 :         return GDALPamRasterBand::GetNoDataValue(pbSuccess);
     173             : 
     174          38 :     if (pbSuccess)
     175          30 :         *pbSuccess = bNoDataSet;
     176             : 
     177          38 :     return dfNoDataValue;
     178             : }
     179             : 
     180             : /************************************************************************/
     181             : /*                             GetUnitType()                            */
     182             : /************************************************************************/
     183             : 
     184          24 : const char *SRTMHGTRasterBand::GetUnitType()
     185             : {
     186          48 :     const std::string osExt = CPLGetExtensionSafe(poDS->GetDescription());
     187          24 :     const char *pszExt = osExt.c_str();
     188          24 :     if (EQUAL(pszExt, "err") || EQUAL(pszExt, "img") || EQUAL(pszExt, "num") ||
     189          24 :         EQUAL(pszExt, "swb"))
     190             :     {
     191           0 :         return "";
     192             :     }
     193          24 :     return "m";
     194             : }
     195             : 
     196             : /************************************************************************/
     197             : /*                       GetColorInterpretation()                       */
     198             : /************************************************************************/
     199             : 
     200          11 : GDALColorInterp SRTMHGTRasterBand::GetColorInterpretation()
     201             : {
     202          11 :     return GCI_Undefined;
     203             : }
     204             : 
     205             : /************************************************************************/
     206             : /* ==================================================================== */
     207             : /*                             SRTMHGTDataset                               */
     208             : /* ==================================================================== */
     209             : /************************************************************************/
     210             : 
     211             : /************************************************************************/
     212             : /*                            SRTMHGTDataset()                              */
     213             : /************************************************************************/
     214             : 
     215          24 : SRTMHGTDataset::SRTMHGTDataset()
     216             : {
     217          24 :     m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     218          24 :     if (CPLTestBool(CPLGetConfigOption("REPORT_COMPD_CS", "NO")))
     219             :     {
     220           0 :         m_oSRS.importFromWkt(
     221             :             "COMPD_CS[\"WGS 84 + EGM96 geoid height\", GEOGCS[\"WGS 84\", "
     222             :             "DATUM[\"WGS_1984\", SPHEROID[\"WGS 84\",6378137,298.257223563, "
     223             :             "AUTHORITY[\"EPSG\",\"7030\"]], AUTHORITY[\"EPSG\",\"6326\"]], "
     224             :             "PRIMEM[\"Greenwich\",0, AUTHORITY[\"EPSG\",\"8901\"]], "
     225             :             "UNIT[\"degree\",0.0174532925199433, "
     226             :             "AUTHORITY[\"EPSG\",\"9122\"]], AUTHORITY[\"EPSG\",\"4326\"]], "
     227             :             "VERT_CS[\"EGM96 geoid height\", VERT_DATUM[\"EGM96 geoid\",2005, "
     228             :             "AUTHORITY[\"EPSG\",\"5171\"]], UNIT[\"metre\",1, "
     229             :             "AUTHORITY[\"EPSG\",\"9001\"]], AXIS[\"Up\",UP], "
     230             :             "AUTHORITY[\"EPSG\",\"5773\"]]]");
     231             :     }
     232             :     else
     233             :     {
     234          24 :         m_oSRS.importFromWkt(SRS_WKT_WGS84_LAT_LONG);
     235             :     }
     236          24 :     adfGeoTransform[0] = 0.0;
     237          24 :     adfGeoTransform[1] = 1.0;
     238          24 :     adfGeoTransform[2] = 0.0;
     239          24 :     adfGeoTransform[3] = 0.0;
     240          24 :     adfGeoTransform[4] = 0.0;
     241          24 :     adfGeoTransform[5] = 1.0;
     242          24 : }
     243             : 
     244             : /************************************************************************/
     245             : /*                           ~SRTMHGTDataset()                            */
     246             : /************************************************************************/
     247             : 
     248          48 : SRTMHGTDataset::~SRTMHGTDataset()
     249             : {
     250          24 :     FlushCache(true);
     251          24 :     if (fpImage != nullptr)
     252          24 :         VSIFCloseL(fpImage);
     253          24 :     CPLFree(pabyBuffer);
     254          48 : }
     255             : 
     256             : /************************************************************************/
     257             : /*                          GetGeoTransform()                           */
     258             : /************************************************************************/
     259             : 
     260          33 : CPLErr SRTMHGTDataset::GetGeoTransform(double *padfTransform)
     261             : {
     262          33 :     memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
     263          33 :     return CE_None;
     264             : }
     265             : 
     266             : /************************************************************************/
     267             : /*                          GetSpatialRef()                             */
     268             : /************************************************************************/
     269             : 
     270          27 : const OGRSpatialReference *SRTMHGTDataset::GetSpatialRef() const
     271             : 
     272             : {
     273          27 :     return &m_oSRS;
     274             : }
     275             : 
     276             : /************************************************************************/
     277             : /*                              Identify()                              */
     278             : /************************************************************************/
     279             : 
     280       55999 : int SRTMHGTDataset::Identify(GDALOpenInfo *poOpenInfo)
     281             : 
     282             : {
     283       55999 :     const char *fileName = CPLGetFilename(poOpenInfo->pszFilename);
     284       55999 :     if (strlen(fileName) < 11 || fileName[7] != '.')
     285       52643 :         return FALSE;
     286        6712 :     CPLString osLCFilename(CPLString(fileName).tolower());
     287        3837 :     if ((osLCFilename[0] != 'n' && osLCFilename[0] != 's') ||
     288         481 :         (osLCFilename[3] != 'e' && osLCFilename[3] != 'w'))
     289        3226 :         return FALSE;
     290         130 :     if (!STARTS_WITH(fileName, "/vsizip/") && osLCFilename.endsWith(".hgt.zip"))
     291             :     {
     292           4 :         CPLString osNewName("/vsizip/");
     293           2 :         osNewName += poOpenInfo->pszFilename;
     294           2 :         osNewName += "/";
     295           2 :         osNewName += CPLString(fileName).substr(0, 7);
     296           2 :         osNewName += ".hgt";
     297           4 :         GDALOpenInfo oOpenInfo(osNewName, GA_ReadOnly);
     298           2 :         return Identify(&oOpenInfo);
     299             :     }
     300             : 
     301         256 :     if (!STARTS_WITH(fileName, "/vsizip/") &&
     302         256 :         osLCFilename.endsWith(".srtmswbd.raw.zip"))
     303             :     {
     304           4 :         CPLString osNewName("/vsizip/");
     305           2 :         osNewName += poOpenInfo->pszFilename;
     306           2 :         osNewName += "/";
     307           2 :         osNewName += CPLString(fileName).substr(0, 7);
     308           2 :         osNewName += ".raw";
     309           4 :         GDALOpenInfo oOpenInfo(osNewName, GA_ReadOnly);
     310           2 :         return Identify(&oOpenInfo);
     311             :     }
     312             : 
     313             :     // .hgts and .err files from
     314             :     // https://e4ftl01.cr.usgs.gov/MEASURES/NASADEM_SHHP.001/2000.02.11/ .img
     315             :     // and .img.num files from
     316             :     // https://e4ftl01.cr.usgs.gov/MEASURES/NASADEM_SIM.001/2000.02.11/
     317         257 :     if (!osLCFilename.endsWith(".hgt") && !osLCFilename.endsWith(".hgts") &&
     318         129 :         !osLCFilename.endsWith(".err") && !osLCFilename.endsWith(".img") &&
     319         129 :         !osLCFilename.endsWith(".num") &&  // .img.num or .num
     320         129 :         !osLCFilename.endsWith(".raw") &&
     321         126 :         !osLCFilename.endsWith(
     322         257 :             ".swb") &&  // https://e4ftl01.cr.usgs.gov/MEASURES/NASADEM_HGT.001/2000.02.11/
     323         126 :         !osLCFilename.endsWith(".hgt.gz"))
     324           0 :         return FALSE;
     325             : 
     326             :     /* -------------------------------------------------------------------- */
     327             :     /*      We check the file size to see if it is                          */
     328             :     /*      SRTM1 (below or above lat 50) or SRTM 3                         */
     329             :     /* -------------------------------------------------------------------- */
     330             :     VSIStatBufL fileStat;
     331             : 
     332         126 :     if (VSIStatL(poOpenInfo->pszFilename, &fileStat) != 0)
     333          87 :         return FALSE;
     334          39 :     if (fileStat.st_size != 1201 * 1201 * 2 &&
     335          20 :         fileStat.st_size != 1801 * 3601 * 2 &&
     336          17 :         fileStat.st_size != 3601 * 3601 &&
     337          12 :         fileStat.st_size != 3601 * 3601 * 2 &&
     338           7 :         fileStat.st_size != 3601 * 3601 * 4 &&  // .hgts
     339           3 :         fileStat.st_size != 7201 * 7201 * 2)
     340           0 :         return FALSE;
     341             : 
     342          39 :     return TRUE;
     343             : }
     344             : 
     345             : /************************************************************************/
     346             : /*                                Open()                                */
     347             : /************************************************************************/
     348             : 
     349          13 : GDALDataset *SRTMHGTDataset::Open(GDALOpenInfo *poOpenInfo)
     350             : {
     351          13 :     return OpenPAM(poOpenInfo);
     352             : }
     353             : 
     354             : /************************************************************************/
     355             : /*                              OpenPAM()                               */
     356             : /************************************************************************/
     357             : 
     358          26 : GDALPamDataset *SRTMHGTDataset::OpenPAM(GDALOpenInfo *poOpenInfo)
     359             : {
     360          26 :     if (!Identify(poOpenInfo))
     361           0 :         return nullptr;
     362             : 
     363          26 :     const char *fileName = CPLGetFilename(poOpenInfo->pszFilename);
     364          52 :     CPLString osLCFilename(CPLString(fileName).tolower());
     365          26 :     if (!STARTS_WITH(fileName, "/vsizip/") && osLCFilename.endsWith(".hgt.zip"))
     366             :     {
     367           2 :         CPLString osFilename("/vsizip/");
     368           1 :         osFilename += poOpenInfo->pszFilename;
     369           1 :         osFilename += "/";
     370           1 :         osFilename += CPLString(fileName).substr(0, 7);
     371           1 :         osFilename += ".hgt";
     372           1 :         GDALOpenInfo oOpenInfo(osFilename, poOpenInfo->eAccess);
     373           1 :         auto poDS = OpenPAM(&oOpenInfo);
     374           1 :         if (poDS != nullptr)
     375             :         {
     376             :             // override description with the main one
     377           1 :             poDS->SetDescription(poOpenInfo->pszFilename);
     378             :         }
     379           1 :         return poDS;
     380             :     }
     381             : 
     382          50 :     if (!STARTS_WITH(fileName, "/vsizip/") &&
     383          50 :         osLCFilename.endsWith(".srtmswbd.raw.zip"))
     384             :     {
     385           2 :         CPLString osFilename("/vsizip/");
     386           1 :         osFilename += poOpenInfo->pszFilename;
     387           1 :         osFilename += "/";
     388           1 :         osFilename += CPLString(fileName).substr(0, 7);
     389           1 :         osFilename += ".raw";
     390           1 :         GDALOpenInfo oOpenInfo(osFilename, poOpenInfo->eAccess);
     391           1 :         auto poDS = OpenPAM(&oOpenInfo);
     392           1 :         if (poDS != nullptr)
     393             :         {
     394             :             // override description with the main one
     395           1 :             poDS->SetDescription(poOpenInfo->pszFilename);
     396             :         }
     397           1 :         return poDS;
     398             :     }
     399             : 
     400             :     char latLonValueString[4];
     401          24 :     memset(latLonValueString, 0, 4);
     402          24 :     strncpy(latLonValueString, &fileName[1], 2);
     403          24 :     int southWestLat = atoi(latLonValueString);
     404          24 :     memset(latLonValueString, 0, 4);
     405             :     // cppcheck-suppress redundantCopy
     406          24 :     strncpy(latLonValueString, &fileName[4], 3);
     407          24 :     int southWestLon = atoi(latLonValueString);
     408             : 
     409          24 :     if (fileName[0] == 'N' || fileName[0] == 'n')
     410             :         /*southWestLat = southWestLat */;
     411           0 :     else if (fileName[0] == 'S' || fileName[0] == 's')
     412           0 :         southWestLat = southWestLat * -1;
     413             :     else
     414           0 :         return nullptr;
     415             : 
     416          24 :     if (fileName[3] == 'E' || fileName[3] == 'e')
     417             :         /*southWestLon = southWestLon */;
     418          11 :     else if (fileName[3] == 'W' || fileName[3] == 'w')
     419          11 :         southWestLon = southWestLon * -1;
     420             :     else
     421           0 :         return nullptr;
     422             : 
     423             :     /* -------------------------------------------------------------------- */
     424             :     /*      Create a corresponding GDALDataset.                             */
     425             :     /* -------------------------------------------------------------------- */
     426          48 :     auto poDS = std::make_unique<SRTMHGTDataset>();
     427             : 
     428          24 :     poDS->fpImage = poOpenInfo->fpL;
     429          24 :     poOpenInfo->fpL = nullptr;
     430             : 
     431             :     VSIStatBufL fileStat;
     432          24 :     if (VSIStatL(poOpenInfo->pszFilename, &fileStat) != 0)
     433             :     {
     434           0 :         return nullptr;
     435             :     }
     436             : 
     437             :     int numPixels_x, numPixels_y;
     438             : 
     439          24 :     GDALDataType eDT = GDT_Int16;
     440          24 :     switch (fileStat.st_size)
     441             :     {
     442          12 :         case 1201 * 1201 * 2:
     443          12 :             numPixels_x = numPixels_y = 1201;
     444          12 :             break;
     445           2 :         case 1801 * 3601 * 2:
     446           2 :             numPixels_x = 1801;
     447           2 :             numPixels_y = 3601;
     448           2 :             break;
     449           2 :         case 3601 * 3601:
     450           2 :             numPixels_x = numPixels_y = 3601;
     451           2 :             eDT = GDT_Byte;
     452           2 :             break;
     453           4 :         case 3601 * 3601 * 2:
     454           4 :             numPixels_x = numPixels_y = 3601;
     455           4 :             break;
     456           2 :         case 3601 * 3601 * 4:  // .hgts
     457           2 :             numPixels_x = numPixels_y = 3601;
     458           2 :             eDT = GDT_Float32;
     459           2 :             break;
     460           2 :         case 7201 * 7201 * 2:
     461           2 :             numPixels_x = numPixels_y = 7201;
     462           2 :             break;
     463           0 :         default:
     464           0 :             numPixels_x = numPixels_y = 0;
     465           0 :             break;
     466             :     }
     467             : 
     468          24 :     poDS->eAccess = poOpenInfo->eAccess;
     469             : #ifdef CPL_LSB
     470          24 :     if (poDS->eAccess == GA_Update && eDT != GDT_Byte)
     471             :     {
     472           1 :         poDS->pabyBuffer =
     473           1 :             static_cast<GByte *>(CPLMalloc(numPixels_x * sizeof(eDT)));
     474             :     }
     475             : #endif
     476             : 
     477             :     /* -------------------------------------------------------------------- */
     478             :     /*      Capture some information from the file that is of interest.     */
     479             :     /* -------------------------------------------------------------------- */
     480          24 :     poDS->nRasterXSize = numPixels_x;
     481          24 :     poDS->nRasterYSize = numPixels_y;
     482          24 :     poDS->nBands = 1;
     483             : 
     484          24 :     poDS->adfGeoTransform[0] = southWestLon - 0.5 / (numPixels_x - 1);
     485          24 :     poDS->adfGeoTransform[1] = 1.0 / (numPixels_x - 1);
     486          24 :     poDS->adfGeoTransform[2] = 0.0;
     487          24 :     poDS->adfGeoTransform[3] = southWestLat + 1 + 0.5 / (numPixels_y - 1);
     488          24 :     poDS->adfGeoTransform[4] = 0.0;
     489          24 :     poDS->adfGeoTransform[5] = -1.0 / (numPixels_y - 1);
     490             : 
     491          24 :     poDS->SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
     492             : 
     493             :     /* -------------------------------------------------------------------- */
     494             :     /*      Create band information object.                                 */
     495             :     /* -------------------------------------------------------------------- */
     496          24 :     SRTMHGTRasterBand *tmpBand = new SRTMHGTRasterBand(poDS.get(), 1, eDT);
     497          24 :     poDS->SetBand(1, tmpBand);
     498             : 
     499             :     /* -------------------------------------------------------------------- */
     500             :     /*      Initialize any PAM information.                                 */
     501             :     /* -------------------------------------------------------------------- */
     502          24 :     poDS->SetDescription(poOpenInfo->pszFilename);
     503          24 :     poDS->TryLoadXML();
     504             : 
     505             :     /* -------------------------------------------------------------------- */
     506             :     /*      Support overviews.                                              */
     507             :     /* -------------------------------------------------------------------- */
     508          24 :     poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
     509             : 
     510          24 :     return poDS.release();
     511             : }
     512             : 
     513             : /************************************************************************/
     514             : /*                              CreateCopy()                            */
     515             : /************************************************************************/
     516             : 
     517          29 : GDALDataset *SRTMHGTDataset::CreateCopy(const char *pszFilename,
     518             :                                         GDALDataset *poSrcDS, int bStrict,
     519             :                                         char ** /* papszOptions*/,
     520             :                                         GDALProgressFunc pfnProgress,
     521             :                                         void *pProgressData)
     522             : {
     523             :     /* -------------------------------------------------------------------- */
     524             :     /*      Some some rudimentary checks                                    */
     525             :     /* -------------------------------------------------------------------- */
     526          29 :     const int nBands = poSrcDS->GetRasterCount();
     527          29 :     if (nBands == 0)
     528             :     {
     529           1 :         CPLError(
     530             :             CE_Failure, CPLE_NotSupported,
     531             :             "SRTMHGT driver does not support source dataset with zero band.\n");
     532           1 :         return nullptr;
     533             :     }
     534          28 :     else if (nBands != 1)
     535             :     {
     536           4 :         CPLError((bStrict) ? CE_Failure : CE_Warning, CPLE_NotSupported,
     537             :                  "SRTMHGT driver only uses the first band of the dataset.\n");
     538           4 :         if (bStrict)
     539           4 :             return nullptr;
     540             :     }
     541             : 
     542             :     /* -------------------------------------------------------------------- */
     543             :     /*      Checks the input SRS                                            */
     544             :     /* -------------------------------------------------------------------- */
     545          48 :     OGRSpatialReference ogrsr_input;
     546          24 :     ogrsr_input.importFromWkt(poSrcDS->GetProjectionRef());
     547             : 
     548          48 :     OGRSpatialReference ogrsr_wgs84;
     549          24 :     ogrsr_wgs84.SetWellKnownGeogCS("WGS84");
     550             : 
     551          24 :     if (ogrsr_input.IsSameGeogCS(&ogrsr_wgs84) == FALSE)
     552             :     {
     553           0 :         CPLError(CE_Warning, CPLE_AppDefined,
     554             :                  "The source projection coordinate system is %s. Only WGS 84 "
     555             :                  "is supported.\nThe SRTMHGT driver will generate a file as "
     556             :                  "if the source was WGS 84 projection coordinate system.",
     557             :                  poSrcDS->GetProjectionRef());
     558             :     }
     559             : 
     560             :     /* -------------------------------------------------------------------- */
     561             :     /*      Work out the LL origin.                                         */
     562             :     /* -------------------------------------------------------------------- */
     563             :     double adfGeoTransform[6];
     564          24 :     if (poSrcDS->GetGeoTransform(adfGeoTransform) != CE_None)
     565             :     {
     566           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     567             :                  "Source image must have a geo transform matrix.");
     568           0 :         return nullptr;
     569             :     }
     570             : 
     571             :     const int nLLOriginLat = static_cast<int>(
     572          48 :         std::floor(adfGeoTransform[3] +
     573          24 :                    poSrcDS->GetRasterYSize() * adfGeoTransform[5] + 0.5));
     574             : 
     575          24 :     int nLLOriginLong = static_cast<int>(std::floor(adfGeoTransform[0] + 0.5));
     576             : 
     577          72 :     if (std::abs(nLLOriginLat -
     578          24 :                  (adfGeoTransform[3] + (poSrcDS->GetRasterYSize() - 0.5) *
     579          35 :                                            adfGeoTransform[5])) > 1e-10 ||
     580          11 :         std::abs(nLLOriginLong -
     581          11 :                  (adfGeoTransform[0] + 0.5 * adfGeoTransform[1])) > 1e-10)
     582             :     {
     583          13 :         CPLError(CE_Warning, CPLE_AppDefined,
     584             :                  "The corner coordinates of the source are not properly "
     585             :                  "aligned on plain latitude/longitude boundaries.");
     586             :     }
     587             : 
     588             :     /* -------------------------------------------------------------------- */
     589             :     /*      Check image dimensions.                                         */
     590             :     /* -------------------------------------------------------------------- */
     591          24 :     const int nXSize = poSrcDS->GetRasterXSize();
     592          24 :     const int nYSize = poSrcDS->GetRasterYSize();
     593             : 
     594          36 :     if (!((nXSize == 1201 && nYSize == 1201) ||
     595          16 :           (nXSize == 1801 && nYSize == 3601) ||
     596           3 :           (nXSize == 3601 && nYSize == 3601) ||
     597           1 :           (nXSize == 7201 && nYSize == 7201)))
     598             :     {
     599          11 :         CPLError(CE_Failure, CPLE_AppDefined,
     600             :                  "Image dimensions should be 1201x1201, 1801x3601, 3601x3601 "
     601             :                  "or 7201x7201.");
     602          11 :         return nullptr;
     603             :     }
     604             : 
     605             :     /* -------------------------------------------------------------------- */
     606             :     /*      Check filename.                                                 */
     607             :     /* -------------------------------------------------------------------- */
     608             :     char expectedFileName[12];
     609             : 
     610          26 :     CPLsnprintf(expectedFileName, sizeof(expectedFileName), "%c%02d%c%03d.HGT",
     611             :                 (nLLOriginLat >= 0) ? 'N' : 'S',
     612          13 :                 (nLLOriginLat >= 0) ? nLLOriginLat : -nLLOriginLat,
     613             :                 (nLLOriginLong >= 0) ? 'E' : 'W',
     614             :                 (nLLOriginLong >= 0) ? nLLOriginLong : -nLLOriginLong);
     615             : 
     616          13 :     if (!EQUAL(expectedFileName, CPLGetFilename(pszFilename)))
     617             :     {
     618           0 :         CPLError(CE_Warning, CPLE_AppDefined, "Expected output filename is %s.",
     619             :                  expectedFileName);
     620             :     }
     621             : 
     622             :     /* -------------------------------------------------------------------- */
     623             :     /*      Write output file.                                              */
     624             :     /* -------------------------------------------------------------------- */
     625          13 :     VSILFILE *fp = VSIFOpenL(pszFilename, "wb");
     626          13 :     if (fp == nullptr)
     627             :     {
     628           2 :         CPLError(CE_Failure, CPLE_FileIO, "Cannot create file %s", pszFilename);
     629           2 :         return nullptr;
     630             :     }
     631             : 
     632             :     GInt16 *panData =
     633          11 :         reinterpret_cast<GInt16 *>(CPLMalloc(sizeof(GInt16) * nXSize));
     634          11 :     GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(1);
     635             : 
     636             :     int bSrcBandHasNoData;
     637          11 :     double srcBandNoData = poSrcBand->GetNoDataValue(&bSrcBandHasNoData);
     638             : 
     639       28822 :     for (int iY = 0; iY < nYSize; iY++)
     640             :     {
     641       28811 :         if (poSrcBand->RasterIO(GF_Read, 0, iY, nXSize, 1,
     642             :                                 reinterpret_cast<void *>(panData), nXSize, 1,
     643       28811 :                                 GDT_Int16, 0, 0, nullptr) != CE_None)
     644             :         {
     645           0 :             VSIFCloseL(fp);
     646           0 :             CPLFree(panData);
     647           0 :             return nullptr;
     648             :         }
     649             : 
     650             :         /* Translate nodata values */
     651       28811 :         if (bSrcBandHasNoData && srcBandNoData != SRTMHG_NODATA_VALUE)
     652             :         {
     653           0 :             for (int iX = 0; iX < nXSize; iX++)
     654             :             {
     655           0 :                 if (panData[iX] == srcBandNoData)
     656           0 :                     panData[iX] = SRTMHG_NODATA_VALUE;
     657             :             }
     658             :         }
     659             : 
     660             : #ifdef CPL_LSB
     661       28811 :         GDALSwapWords(panData, 2, nXSize, 2);
     662             : #endif
     663             : 
     664       28811 :         if (VSIFWriteL(panData, sizeof(GInt16) * nXSize, 1, fp) != 1)
     665             :         {
     666           0 :             CPLError(CE_Failure, CPLE_FileIO,
     667             :                      "Failed to write line %d in SRTMHGT dataset.\n", iY);
     668           0 :             VSIFCloseL(fp);
     669           0 :             CPLFree(panData);
     670           0 :             return nullptr;
     671             :         }
     672             : 
     673       28811 :         if (pfnProgress && !pfnProgress((iY + 1) / static_cast<double>(nYSize),
     674             :                                         nullptr, pProgressData))
     675             :         {
     676           0 :             CPLError(CE_Failure, CPLE_UserInterrupt,
     677             :                      "User terminated CreateCopy()");
     678           0 :             VSIFCloseL(fp);
     679           0 :             CPLFree(panData);
     680           0 :             return nullptr;
     681             :         }
     682             :     }
     683             : 
     684          11 :     CPLFree(panData);
     685          11 :     VSIFCloseL(fp);
     686             : 
     687             :     /* -------------------------------------------------------------------- */
     688             :     /*      Reopen and copy missing information into a PAM file.            */
     689             :     /* -------------------------------------------------------------------- */
     690          11 :     GDALOpenInfo oOpenInfo(pszFilename, GA_ReadOnly);
     691          11 :     auto poDS = OpenPAM(&oOpenInfo);
     692             : 
     693          11 :     if (poDS)
     694          11 :         poDS->CloneInfo(poSrcDS, GCIF_PAM_DEFAULT);
     695             : 
     696          11 :     return poDS;
     697             : }
     698             : 
     699             : /************************************************************************/
     700             : /*                         GDALRegister_SRTMHGT()                       */
     701             : /************************************************************************/
     702        1682 : void GDALRegister_SRTMHGT()
     703             : {
     704        1682 :     if (GDALGetDriverByName("SRTMHGT") != nullptr)
     705         301 :         return;
     706             : 
     707        1381 :     GDALDriver *poDriver = new GDALDriver();
     708             : 
     709        1381 :     poDriver->SetDescription("SRTMHGT");
     710        1381 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     711        1381 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "SRTMHGT File Format");
     712        1381 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "hgt");
     713        1381 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
     714        1381 :                               "drivers/raster/srtmhgt.html");
     715        1381 :     poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Byte Int16 UInt16");
     716        1381 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     717             : 
     718        1381 :     poDriver->pfnIdentify = SRTMHGTDataset::Identify;
     719        1381 :     poDriver->pfnOpen = SRTMHGTDataset::Open;
     720        1381 :     poDriver->pfnCreateCopy = SRTMHGTDataset::CreateCopy;
     721             : 
     722        1381 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     723             : }

Generated by: LCOV version 1.14