LCOV - code coverage report
Current view: top level - frmts/raw - snodasdataset.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 221 239 92.5 %
Date: 2025-07-09 17:50:03 Functions: 16 16 100.0 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  SNODAS driver
       4             :  * Purpose:  Implementation of SNODASDataset
       5             :  * Author:   Even Rouault, <even dot rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2011, Even Rouault <even dot rouault at spatialys.com>
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "cpl_string.h"
      14             : #include "gdal_frmts.h"
      15             : #include "ogr_srs_api.h"
      16             : #include "rawdataset.h"
      17             : 
      18             : /************************************************************************/
      19             : /* ==================================================================== */
      20             : /*                            SNODASDataset                             */
      21             : /* ==================================================================== */
      22             : /************************************************************************/
      23             : 
      24             : class SNODASRasterBand;
      25             : 
      26             : class SNODASDataset final : public RawDataset
      27             : {
      28             :     CPLString osDataFilename{};
      29             :     bool bGotTransform;
      30             :     GDALGeoTransform m_gt{};
      31             :     bool bHasNoData;
      32             :     double dfNoData;
      33             :     bool bHasMin;
      34             :     double dfMin;
      35             :     int bHasMax;
      36             :     double dfMax;
      37             :     OGRSpatialReference m_oSRS{};
      38             : 
      39             :     friend class SNODASRasterBand;
      40             : 
      41             :     CPL_DISALLOW_COPY_ASSIGN(SNODASDataset)
      42             : 
      43             :     CPLErr Close() override;
      44             : 
      45             :   public:
      46             :     SNODASDataset();
      47             :     ~SNODASDataset() override;
      48             : 
      49             :     CPLErr GetGeoTransform(GDALGeoTransform &gt) const override;
      50             : 
      51           1 :     const OGRSpatialReference *GetSpatialRef() const override
      52             :     {
      53           1 :         return &m_oSRS;
      54             :     }
      55             : 
      56             :     char **GetFileList() override;
      57             : 
      58             :     static GDALDataset *Open(GDALOpenInfo *);
      59             :     static int Identify(GDALOpenInfo *);
      60             : };
      61             : 
      62             : /************************************************************************/
      63             : /* ==================================================================== */
      64             : /*                            SNODASRasterBand                          */
      65             : /* ==================================================================== */
      66             : /************************************************************************/
      67             : 
      68             : class SNODASRasterBand final : public RawRasterBand
      69             : {
      70             :     CPL_DISALLOW_COPY_ASSIGN(SNODASRasterBand)
      71             : 
      72             :   public:
      73             :     SNODASRasterBand(VSILFILE *fpRaw, int nXSize, int nYSize);
      74             : 
      75           6 :     ~SNODASRasterBand() override
      76           3 :     {
      77           6 :     }
      78             : 
      79             :     double GetNoDataValue(int *pbSuccess = nullptr) override;
      80             :     double GetMinimum(int *pbSuccess = nullptr) override;
      81             :     double GetMaximum(int *pbSuccess = nullptr) override;
      82             : };
      83             : 
      84             : /************************************************************************/
      85             : /*                         SNODASRasterBand()                           */
      86             : /************************************************************************/
      87             : 
      88           3 : SNODASRasterBand::SNODASRasterBand(VSILFILE *fpRawIn, int nXSize, int nYSize)
      89             :     : RawRasterBand(fpRawIn, 0, 2, nXSize * 2, GDT_Int16, !CPL_IS_LSB, nXSize,
      90           3 :                     nYSize, RawRasterBand::OwnFP::YES)
      91             : {
      92           3 : }
      93             : 
      94             : /************************************************************************/
      95             : /*                          GetNoDataValue()                            */
      96             : /************************************************************************/
      97             : 
      98           1 : double SNODASRasterBand::GetNoDataValue(int *pbSuccess)
      99             : {
     100           1 :     SNODASDataset *poGDS = reinterpret_cast<SNODASDataset *>(poDS);
     101           1 :     if (pbSuccess)
     102           1 :         *pbSuccess = poGDS->bHasNoData;
     103             : 
     104           1 :     if (poGDS->bHasNoData)
     105           1 :         return poGDS->dfNoData;
     106             : 
     107           0 :     return RawRasterBand::GetNoDataValue(pbSuccess);
     108             : }
     109             : 
     110             : /************************************************************************/
     111             : /*                            GetMinimum()                              */
     112             : /************************************************************************/
     113             : 
     114           1 : double SNODASRasterBand::GetMinimum(int *pbSuccess)
     115             : {
     116           1 :     SNODASDataset *poGDS = reinterpret_cast<SNODASDataset *>(poDS);
     117           1 :     if (pbSuccess)
     118           1 :         *pbSuccess = poGDS->bHasMin;
     119             : 
     120           1 :     if (poGDS->bHasMin)
     121           1 :         return poGDS->dfMin;
     122             : 
     123           0 :     return RawRasterBand::GetMinimum(pbSuccess);
     124             : }
     125             : 
     126             : /************************************************************************/
     127             : /*                            GetMaximum()                             */
     128             : /************************************************************************/
     129             : 
     130           1 : double SNODASRasterBand::GetMaximum(int *pbSuccess)
     131             : {
     132           1 :     SNODASDataset *poGDS = reinterpret_cast<SNODASDataset *>(poDS);
     133           1 :     if (pbSuccess)
     134           1 :         *pbSuccess = poGDS->bHasMax;
     135             : 
     136           1 :     if (poGDS->bHasMax)
     137           1 :         return poGDS->dfMax;
     138             : 
     139           0 :     return RawRasterBand::GetMaximum(pbSuccess);
     140             : }
     141             : 
     142             : /************************************************************************/
     143             : /* ==================================================================== */
     144             : /*                            SNODASDataset                             */
     145             : /* ==================================================================== */
     146             : /************************************************************************/
     147             : 
     148             : /************************************************************************/
     149             : /*                           SNODASDataset()                            */
     150             : /************************************************************************/
     151             : 
     152           3 : SNODASDataset::SNODASDataset()
     153             :     : bGotTransform(false), bHasNoData(false), dfNoData(0.0), bHasMin(false),
     154           3 :       dfMin(0.0), bHasMax(false), dfMax(0.0)
     155             : {
     156           3 :     m_oSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
     157           3 :     m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     158           3 : }
     159             : 
     160             : /************************************************************************/
     161             : /*                           ~SNODASDataset()                           */
     162             : /************************************************************************/
     163             : 
     164           6 : SNODASDataset::~SNODASDataset()
     165             : 
     166             : {
     167           3 :     SNODASDataset::Close();
     168           6 : }
     169             : 
     170             : /************************************************************************/
     171             : /*                              Close()                                 */
     172             : /************************************************************************/
     173             : 
     174           6 : CPLErr SNODASDataset::Close()
     175             : {
     176           6 :     CPLErr eErr = CE_None;
     177           6 :     if (nOpenFlags != OPEN_FLAGS_CLOSED)
     178             :     {
     179           3 :         if (SNODASDataset::FlushCache(true) != CE_None)
     180           0 :             eErr = CE_Failure;
     181             : 
     182           3 :         if (GDALPamDataset::Close() != CE_None)
     183           0 :             eErr = CE_Failure;
     184             :     }
     185           6 :     return eErr;
     186             : }
     187             : 
     188             : /************************************************************************/
     189             : /*                          GetGeoTransform()                           */
     190             : /************************************************************************/
     191             : 
     192           1 : CPLErr SNODASDataset::GetGeoTransform(GDALGeoTransform &gt) const
     193             : 
     194             : {
     195           1 :     if (bGotTransform)
     196             :     {
     197           1 :         gt = m_gt;
     198           1 :         return CE_None;
     199             :     }
     200             : 
     201           0 :     return GDALPamDataset::GetGeoTransform(gt);
     202             : }
     203             : 
     204             : /************************************************************************/
     205             : /*                            GetFileList()                             */
     206             : /************************************************************************/
     207             : 
     208           2 : char **SNODASDataset::GetFileList()
     209             : 
     210             : {
     211           2 :     char **papszFileList = GDALPamDataset::GetFileList();
     212             : 
     213           2 :     papszFileList = CSLAddString(papszFileList, osDataFilename);
     214             : 
     215           2 :     return papszFileList;
     216             : }
     217             : 
     218             : /************************************************************************/
     219             : /*                            Identify()                                */
     220             : /************************************************************************/
     221             : 
     222       57586 : int SNODASDataset::Identify(GDALOpenInfo *poOpenInfo)
     223             : {
     224       57586 :     if (poOpenInfo->nHeaderBytes == 0)
     225       53798 :         return FALSE;
     226             : 
     227        3788 :     return STARTS_WITH_CI(reinterpret_cast<char *>(poOpenInfo->pabyHeader),
     228             :                           "Format version: NOHRSC GIS/RS raster file v1.1");
     229             : }
     230             : 
     231             : /************************************************************************/
     232             : /*                                Open()                                */
     233             : /************************************************************************/
     234             : 
     235           3 : GDALDataset *SNODASDataset::Open(GDALOpenInfo *poOpenInfo)
     236             : 
     237             : {
     238           3 :     if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
     239           0 :         return nullptr;
     240             : 
     241             :     /* -------------------------------------------------------------------- */
     242             :     /*      Confirm the requested access is supported.                      */
     243             :     /* -------------------------------------------------------------------- */
     244           3 :     if (poOpenInfo->eAccess == GA_Update)
     245             :     {
     246           0 :         ReportUpdateNotSupportedByDriver("SNODAS");
     247           0 :         return nullptr;
     248             :     }
     249             : 
     250           3 :     int nRows = -1;
     251           3 :     int nCols = -1;
     252           6 :     CPLString osDataFilename;
     253           3 :     bool bIsInteger = false;
     254           3 :     bool bIs2Bytes = false;
     255           3 :     double dfNoData = 0;
     256           3 :     bool bHasNoData = false;
     257           3 :     double dfMin = 0;
     258           3 :     bool bHasMin = false;
     259           3 :     double dfMax = 0;
     260           3 :     bool bHasMax = false;
     261           3 :     double dfMinX = 0.0;
     262           3 :     double dfMinY = 0.0;
     263           3 :     double dfMaxX = 0.0;
     264           3 :     double dfMaxY = 0.0;
     265           3 :     bool bHasMinX = false;
     266           3 :     bool bHasMinY = false;
     267           3 :     bool bHasMaxX = false;
     268           3 :     bool bHasMaxY = false;
     269           3 :     bool bNotProjected = false;
     270           3 :     bool bIsWGS84 = false;
     271           6 :     CPLString osDataUnits;
     272           6 :     CPLString osDescription;
     273           3 :     int nStartYear = -1;
     274           3 :     int nStartMonth = -1;
     275           3 :     int nStartDay = -1;
     276           3 :     int nStartHour = -1;
     277           3 :     int nStartMinute = -1;
     278           3 :     int nStartSecond = -1;
     279           3 :     int nStopYear = -1;
     280           3 :     int nStopMonth = -1;
     281           3 :     int nStopDay = -1;
     282           3 :     int nStopHour = -1;
     283           3 :     int nStopMinute = -1;
     284           3 :     int nStopSecond = -1;
     285             : 
     286           3 :     const char *pszLine = nullptr;
     287         324 :     while ((pszLine = CPLReadLine2L(poOpenInfo->fpL, 1024, nullptr)) != nullptr)
     288             :     {
     289             :         char **papszTokens =
     290         321 :             CSLTokenizeStringComplex(pszLine, ":", TRUE, FALSE);
     291         321 :         if (CSLCount(papszTokens) != 2)
     292             :         {
     293           3 :             CSLDestroy(papszTokens);
     294           3 :             continue;
     295             :         }
     296         318 :         if (papszTokens[1][0] == ' ')
     297         318 :             memmove(papszTokens[1], papszTokens[1] + 1,
     298         318 :                     strlen(papszTokens[1] + 1) + 1);
     299             : 
     300         318 :         if (EQUAL(papszTokens[0], "Data file pathname"))
     301             :         {
     302           3 :             osDataFilename = papszTokens[1];
     303             :         }
     304         315 :         else if (EQUAL(papszTokens[0], "Description"))
     305             :         {
     306           3 :             osDescription = papszTokens[1];
     307             :         }
     308         312 :         else if (EQUAL(papszTokens[0], "Data units"))
     309             :         {
     310           3 :             osDataUnits = papszTokens[1];
     311             :         }
     312             : 
     313         309 :         else if (EQUAL(papszTokens[0], "Start year"))
     314           3 :             nStartYear = atoi(papszTokens[1]);
     315         306 :         else if (EQUAL(papszTokens[0], "Start month"))
     316           3 :             nStartMonth = atoi(papszTokens[1]);
     317         303 :         else if (EQUAL(papszTokens[0], "Start day"))
     318           3 :             nStartDay = atoi(papszTokens[1]);
     319         300 :         else if (EQUAL(papszTokens[0], "Start hour"))
     320           3 :             nStartHour = atoi(papszTokens[1]);
     321         297 :         else if (EQUAL(papszTokens[0], " Start minute"))
     322           0 :             nStartMinute = atoi(papszTokens[1]);
     323         297 :         else if (EQUAL(papszTokens[0], "Start second"))
     324           3 :             nStartSecond = atoi(papszTokens[1]);
     325             : 
     326         294 :         else if (EQUAL(papszTokens[0], "Stop year"))
     327           3 :             nStopYear = atoi(papszTokens[1]);
     328         291 :         else if (EQUAL(papszTokens[0], "Stop month"))
     329           3 :             nStopMonth = atoi(papszTokens[1]);
     330         288 :         else if (EQUAL(papszTokens[0], "Stop day"))
     331           3 :             nStopDay = atoi(papszTokens[1]);
     332         285 :         else if (EQUAL(papszTokens[0], "Stop hour"))
     333           3 :             nStopHour = atoi(papszTokens[1]);
     334         282 :         else if (EQUAL(papszTokens[0], "Stop minute"))
     335           3 :             nStopMinute = atoi(papszTokens[1]);
     336         279 :         else if (EQUAL(papszTokens[0], "Stop second"))
     337           3 :             nStopSecond = atoi(papszTokens[1]);
     338             : 
     339         276 :         else if (EQUAL(papszTokens[0], "Number of columns"))
     340             :         {
     341           3 :             nCols = atoi(papszTokens[1]);
     342             :         }
     343         273 :         else if (EQUAL(papszTokens[0], "Number of rows"))
     344             :         {
     345           3 :             nRows = atoi(papszTokens[1]);
     346             :         }
     347         270 :         else if (EQUAL(papszTokens[0], "Data type"))
     348             :         {
     349           3 :             bIsInteger = EQUAL(papszTokens[1], "integer");
     350             :         }
     351         267 :         else if (EQUAL(papszTokens[0], "Data bytes per pixel"))
     352             :         {
     353           3 :             bIs2Bytes = EQUAL(papszTokens[1], "2");
     354             :         }
     355         264 :         else if (EQUAL(papszTokens[0], "Projected"))
     356             :         {
     357           3 :             bNotProjected = EQUAL(papszTokens[1], "no");
     358             :         }
     359         261 :         else if (EQUAL(papszTokens[0], "Horizontal datum"))
     360             :         {
     361           3 :             bIsWGS84 = EQUAL(papszTokens[1], "WGS84");
     362             :         }
     363         258 :         else if (EQUAL(papszTokens[0], "No data value"))
     364             :         {
     365           3 :             bHasNoData = true;
     366           3 :             dfNoData = CPLAtofM(papszTokens[1]);
     367             :         }
     368         255 :         else if (EQUAL(papszTokens[0], "Minimum data value"))
     369             :         {
     370           3 :             bHasMin = true;
     371           3 :             dfMin = CPLAtofM(papszTokens[1]);
     372             :         }
     373         252 :         else if (EQUAL(papszTokens[0], "Maximum data value"))
     374             :         {
     375           3 :             bHasMax = true;
     376           3 :             dfMax = CPLAtofM(papszTokens[1]);
     377             :         }
     378         249 :         else if (EQUAL(papszTokens[0], "Minimum x-axis coordinate"))
     379             :         {
     380           3 :             bHasMinX = true;
     381           3 :             dfMinX = CPLAtofM(papszTokens[1]);
     382             :         }
     383         246 :         else if (EQUAL(papszTokens[0], "Minimum y-axis coordinate"))
     384             :         {
     385           3 :             bHasMinY = true;
     386           3 :             dfMinY = CPLAtofM(papszTokens[1]);
     387             :         }
     388         243 :         else if (EQUAL(papszTokens[0], "Maximum x-axis coordinate"))
     389             :         {
     390           3 :             bHasMaxX = true;
     391           3 :             dfMaxX = CPLAtofM(papszTokens[1]);
     392             :         }
     393         240 :         else if (EQUAL(papszTokens[0], "Maximum y-axis coordinate"))
     394             :         {
     395           3 :             bHasMaxY = true;
     396           3 :             dfMaxY = CPLAtofM(papszTokens[1]);
     397             :         }
     398             : 
     399         318 :         CSLDestroy(papszTokens);
     400             :     }
     401             : 
     402           3 :     CPL_IGNORE_RET_VAL(VSIFCloseL(poOpenInfo->fpL));
     403           3 :     poOpenInfo->fpL = nullptr;
     404             : 
     405             :     /* -------------------------------------------------------------------- */
     406             :     /*      Did we get the required keywords?  If not we return with        */
     407             :     /*      this never having been considered to be a match. This isn't     */
     408             :     /*      an error!                                                       */
     409             :     /* -------------------------------------------------------------------- */
     410           3 :     if (nRows == -1 || nCols == -1 || !bIsInteger || !bIs2Bytes)
     411           0 :         return nullptr;
     412             : 
     413           3 :     if (!bNotProjected || !bIsWGS84)
     414           0 :         return nullptr;
     415             : 
     416           3 :     if (osDataFilename.empty())
     417           0 :         return nullptr;
     418             : 
     419           3 :     if (!GDALCheckDatasetDimensions(nCols, nRows))
     420           0 :         return nullptr;
     421             : 
     422             :     /* -------------------------------------------------------------------- */
     423             :     /*      Open target binary file.                                        */
     424             :     /* -------------------------------------------------------------------- */
     425           6 :     const std::string osPath = CPLGetPathSafe(poOpenInfo->pszFilename);
     426             :     osDataFilename =
     427           3 :         CPLFormFilenameSafe(osPath.c_str(), osDataFilename, nullptr);
     428             : 
     429           3 :     VSILFILE *fpRaw = VSIFOpenL(osDataFilename, "rb");
     430             : 
     431           3 :     if (fpRaw == nullptr)
     432           0 :         return nullptr;
     433             : 
     434             :     /* -------------------------------------------------------------------- */
     435             :     /*      Create a corresponding GDALDataset.                             */
     436             :     /* -------------------------------------------------------------------- */
     437           6 :     auto poDS = std::make_unique<SNODASDataset>();
     438             : 
     439           3 :     poDS->nRasterXSize = nCols;
     440           3 :     poDS->nRasterYSize = nRows;
     441           3 :     poDS->osDataFilename = std::move(osDataFilename);
     442           3 :     poDS->bHasNoData = bHasNoData;
     443           3 :     poDS->dfNoData = dfNoData;
     444           3 :     poDS->bHasMin = bHasMin;
     445           3 :     poDS->dfMin = dfMin;
     446           3 :     poDS->bHasMax = bHasMax;
     447           3 :     poDS->dfMax = dfMax;
     448           3 :     if (bHasMinX && bHasMinY && bHasMaxX && bHasMaxY)
     449             :     {
     450           3 :         poDS->bGotTransform = true;
     451           3 :         poDS->m_gt[0] = dfMinX;
     452           3 :         poDS->m_gt[1] = (dfMaxX - dfMinX) / nCols;
     453           3 :         poDS->m_gt[2] = 0.0;
     454           3 :         poDS->m_gt[3] = dfMaxY;
     455           3 :         poDS->m_gt[4] = 0.0;
     456           3 :         poDS->m_gt[5] = -(dfMaxY - dfMinY) / nRows;
     457             :     }
     458             : 
     459           3 :     if (!osDescription.empty())
     460           3 :         poDS->SetMetadataItem("Description", osDescription);
     461           3 :     if (!osDataUnits.empty())
     462           3 :         poDS->SetMetadataItem("Data_Units", osDataUnits);
     463           3 :     if (nStartYear != -1 && nStartMonth != -1 && nStartDay != -1 &&
     464           3 :         nStartHour != -1 && nStartMinute != -1 && nStartSecond != -1)
     465           0 :         poDS->SetMetadataItem(
     466             :             "Start_Date",
     467             :             CPLSPrintf("%04d/%02d/%02d %02d:%02d:%02d", nStartYear, nStartMonth,
     468           0 :                        nStartDay, nStartHour, nStartMinute, nStartSecond));
     469           3 :     if (nStopYear != -1 && nStopMonth != -1 && nStopDay != -1 &&
     470           3 :         nStopHour != -1 && nStopMinute != -1 && nStopSecond != -1)
     471           6 :         poDS->SetMetadataItem("Stop_Date",
     472             :                               CPLSPrintf("%04d/%02d/%02d %02d:%02d:%02d",
     473             :                                          nStopYear, nStopMonth, nStopDay,
     474           3 :                                          nStopHour, nStopMinute, nStopSecond));
     475             : 
     476             :     /* -------------------------------------------------------------------- */
     477             :     /*      Create band information objects.                                */
     478             :     /* -------------------------------------------------------------------- */
     479           6 :     auto poBand = std::make_unique<SNODASRasterBand>(fpRaw, nCols, nRows);
     480           3 :     if (!poBand->IsValid())
     481           0 :         return nullptr;
     482           3 :     poDS->SetBand(1, std::move(poBand));
     483             : 
     484             :     /* -------------------------------------------------------------------- */
     485             :     /*      Initialize any PAM information.                                 */
     486             :     /* -------------------------------------------------------------------- */
     487           3 :     poDS->SetDescription(poOpenInfo->pszFilename);
     488           3 :     poDS->TryLoadXML();
     489             : 
     490             :     /* -------------------------------------------------------------------- */
     491             :     /*      Check for overviews.                                            */
     492             :     /* -------------------------------------------------------------------- */
     493           3 :     poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
     494             : 
     495           3 :     return poDS.release();
     496             : }
     497             : 
     498             : /************************************************************************/
     499             : /*                       GDALRegister_SNODAS()                          */
     500             : /************************************************************************/
     501             : 
     502        1935 : void GDALRegister_SNODAS()
     503             : 
     504             : {
     505        1935 :     if (GDALGetDriverByName("SNODAS") != nullptr)
     506         282 :         return;
     507             : 
     508        1653 :     GDALDriver *poDriver = new GDALDriver();
     509             : 
     510        1653 :     poDriver->SetDescription("SNODAS");
     511        1653 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     512        1653 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
     513        1653 :                               "Snow Data Assimilation System");
     514        1653 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/snodas.html");
     515        1653 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "hdr");
     516        1653 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     517             : 
     518        1653 :     poDriver->pfnOpen = SNODASDataset::Open;
     519        1653 :     poDriver->pfnIdentify = SNODASDataset::Identify;
     520             : 
     521        1653 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     522             : }

Generated by: LCOV version 1.14