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             :     double adfGeoTransform[6];
      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(double *padfTransform) 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             : 
     159           3 :     adfGeoTransform[0] = 0.0;
     160           3 :     adfGeoTransform[1] = 1.0;
     161           3 :     adfGeoTransform[2] = 0.0;
     162           3 :     adfGeoTransform[3] = 0.0;
     163           3 :     adfGeoTransform[4] = 0.0;
     164           3 :     adfGeoTransform[5] = 1.0;
     165           3 : }
     166             : 
     167             : /************************************************************************/
     168             : /*                           ~SNODASDataset()                           */
     169             : /************************************************************************/
     170             : 
     171           6 : SNODASDataset::~SNODASDataset()
     172             : 
     173             : {
     174           3 :     SNODASDataset::Close();
     175           6 : }
     176             : 
     177             : /************************************************************************/
     178             : /*                              Close()                                 */
     179             : /************************************************************************/
     180             : 
     181           6 : CPLErr SNODASDataset::Close()
     182             : {
     183           6 :     CPLErr eErr = CE_None;
     184           6 :     if (nOpenFlags != OPEN_FLAGS_CLOSED)
     185             :     {
     186           3 :         if (SNODASDataset::FlushCache(true) != CE_None)
     187           0 :             eErr = CE_Failure;
     188             : 
     189           3 :         if (GDALPamDataset::Close() != CE_None)
     190           0 :             eErr = CE_Failure;
     191             :     }
     192           6 :     return eErr;
     193             : }
     194             : 
     195             : /************************************************************************/
     196             : /*                          GetGeoTransform()                           */
     197             : /************************************************************************/
     198             : 
     199           1 : CPLErr SNODASDataset::GetGeoTransform(double *padfTransform)
     200             : 
     201             : {
     202           1 :     if (bGotTransform)
     203             :     {
     204           1 :         memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
     205           1 :         return CE_None;
     206             :     }
     207             : 
     208           0 :     return GDALPamDataset::GetGeoTransform(padfTransform);
     209             : }
     210             : 
     211             : /************************************************************************/
     212             : /*                            GetFileList()                             */
     213             : /************************************************************************/
     214             : 
     215           2 : char **SNODASDataset::GetFileList()
     216             : 
     217             : {
     218           2 :     char **papszFileList = GDALPamDataset::GetFileList();
     219             : 
     220           2 :     papszFileList = CSLAddString(papszFileList, osDataFilename);
     221             : 
     222           2 :     return papszFileList;
     223             : }
     224             : 
     225             : /************************************************************************/
     226             : /*                            Identify()                                */
     227             : /************************************************************************/
     228             : 
     229       51413 : int SNODASDataset::Identify(GDALOpenInfo *poOpenInfo)
     230             : {
     231       51413 :     if (poOpenInfo->nHeaderBytes == 0)
     232       47800 :         return FALSE;
     233             : 
     234        3613 :     return STARTS_WITH_CI(reinterpret_cast<char *>(poOpenInfo->pabyHeader),
     235             :                           "Format version: NOHRSC GIS/RS raster file v1.1");
     236             : }
     237             : 
     238             : /************************************************************************/
     239             : /*                                Open()                                */
     240             : /************************************************************************/
     241             : 
     242           3 : GDALDataset *SNODASDataset::Open(GDALOpenInfo *poOpenInfo)
     243             : 
     244             : {
     245           3 :     if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
     246           0 :         return nullptr;
     247             : 
     248             :     /* -------------------------------------------------------------------- */
     249             :     /*      Confirm the requested access is supported.                      */
     250             :     /* -------------------------------------------------------------------- */
     251           3 :     if (poOpenInfo->eAccess == GA_Update)
     252             :     {
     253           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     254             :                  "The SNODAS driver does not support update access to existing"
     255             :                  " datasets.");
     256           0 :         return nullptr;
     257             :     }
     258             : 
     259           3 :     int nRows = -1;
     260           3 :     int nCols = -1;
     261           6 :     CPLString osDataFilename;
     262           3 :     bool bIsInteger = false;
     263           3 :     bool bIs2Bytes = false;
     264           3 :     double dfNoData = 0;
     265           3 :     bool bHasNoData = false;
     266           3 :     double dfMin = 0;
     267           3 :     bool bHasMin = false;
     268           3 :     double dfMax = 0;
     269           3 :     bool bHasMax = false;
     270           3 :     double dfMinX = 0.0;
     271           3 :     double dfMinY = 0.0;
     272           3 :     double dfMaxX = 0.0;
     273           3 :     double dfMaxY = 0.0;
     274           3 :     bool bHasMinX = false;
     275           3 :     bool bHasMinY = false;
     276           3 :     bool bHasMaxX = false;
     277           3 :     bool bHasMaxY = false;
     278           3 :     bool bNotProjected = false;
     279           3 :     bool bIsWGS84 = false;
     280           6 :     CPLString osDataUnits;
     281           6 :     CPLString osDescription;
     282           3 :     int nStartYear = -1;
     283           3 :     int nStartMonth = -1;
     284           3 :     int nStartDay = -1;
     285           3 :     int nStartHour = -1;
     286           3 :     int nStartMinute = -1;
     287           3 :     int nStartSecond = -1;
     288           3 :     int nStopYear = -1;
     289           3 :     int nStopMonth = -1;
     290           3 :     int nStopDay = -1;
     291           3 :     int nStopHour = -1;
     292           3 :     int nStopMinute = -1;
     293           3 :     int nStopSecond = -1;
     294             : 
     295           3 :     const char *pszLine = nullptr;
     296         324 :     while ((pszLine = CPLReadLine2L(poOpenInfo->fpL, 1024, nullptr)) != nullptr)
     297             :     {
     298             :         char **papszTokens =
     299         321 :             CSLTokenizeStringComplex(pszLine, ":", TRUE, FALSE);
     300         321 :         if (CSLCount(papszTokens) != 2)
     301             :         {
     302           3 :             CSLDestroy(papszTokens);
     303           3 :             continue;
     304             :         }
     305         318 :         if (papszTokens[1][0] == ' ')
     306         318 :             memmove(papszTokens[1], papszTokens[1] + 1,
     307         318 :                     strlen(papszTokens[1] + 1) + 1);
     308             : 
     309         318 :         if (EQUAL(papszTokens[0], "Data file pathname"))
     310             :         {
     311           3 :             osDataFilename = papszTokens[1];
     312             :         }
     313         315 :         else if (EQUAL(papszTokens[0], "Description"))
     314             :         {
     315           3 :             osDescription = papszTokens[1];
     316             :         }
     317         312 :         else if (EQUAL(papszTokens[0], "Data units"))
     318             :         {
     319           3 :             osDataUnits = papszTokens[1];
     320             :         }
     321             : 
     322         309 :         else if (EQUAL(papszTokens[0], "Start year"))
     323           3 :             nStartYear = atoi(papszTokens[1]);
     324         306 :         else if (EQUAL(papszTokens[0], "Start month"))
     325           3 :             nStartMonth = atoi(papszTokens[1]);
     326         303 :         else if (EQUAL(papszTokens[0], "Start day"))
     327           3 :             nStartDay = atoi(papszTokens[1]);
     328         300 :         else if (EQUAL(papszTokens[0], "Start hour"))
     329           3 :             nStartHour = atoi(papszTokens[1]);
     330         297 :         else if (EQUAL(papszTokens[0], " Start minute"))
     331           0 :             nStartMinute = atoi(papszTokens[1]);
     332         297 :         else if (EQUAL(papszTokens[0], "Start second"))
     333           3 :             nStartSecond = atoi(papszTokens[1]);
     334             : 
     335         294 :         else if (EQUAL(papszTokens[0], "Stop year"))
     336           3 :             nStopYear = atoi(papszTokens[1]);
     337         291 :         else if (EQUAL(papszTokens[0], "Stop month"))
     338           3 :             nStopMonth = atoi(papszTokens[1]);
     339         288 :         else if (EQUAL(papszTokens[0], "Stop day"))
     340           3 :             nStopDay = atoi(papszTokens[1]);
     341         285 :         else if (EQUAL(papszTokens[0], "Stop hour"))
     342           3 :             nStopHour = atoi(papszTokens[1]);
     343         282 :         else if (EQUAL(papszTokens[0], "Stop minute"))
     344           3 :             nStopMinute = atoi(papszTokens[1]);
     345         279 :         else if (EQUAL(papszTokens[0], "Stop second"))
     346           3 :             nStopSecond = atoi(papszTokens[1]);
     347             : 
     348         276 :         else if (EQUAL(papszTokens[0], "Number of columns"))
     349             :         {
     350           3 :             nCols = atoi(papszTokens[1]);
     351             :         }
     352         273 :         else if (EQUAL(papszTokens[0], "Number of rows"))
     353             :         {
     354           3 :             nRows = atoi(papszTokens[1]);
     355             :         }
     356         270 :         else if (EQUAL(papszTokens[0], "Data type"))
     357             :         {
     358           3 :             bIsInteger = EQUAL(papszTokens[1], "integer");
     359             :         }
     360         267 :         else if (EQUAL(papszTokens[0], "Data bytes per pixel"))
     361             :         {
     362           3 :             bIs2Bytes = EQUAL(papszTokens[1], "2");
     363             :         }
     364         264 :         else if (EQUAL(papszTokens[0], "Projected"))
     365             :         {
     366           3 :             bNotProjected = EQUAL(papszTokens[1], "no");
     367             :         }
     368         261 :         else if (EQUAL(papszTokens[0], "Horizontal datum"))
     369             :         {
     370           3 :             bIsWGS84 = EQUAL(papszTokens[1], "WGS84");
     371             :         }
     372         258 :         else if (EQUAL(papszTokens[0], "No data value"))
     373             :         {
     374           3 :             bHasNoData = true;
     375           3 :             dfNoData = CPLAtofM(papszTokens[1]);
     376             :         }
     377         255 :         else if (EQUAL(papszTokens[0], "Minimum data value"))
     378             :         {
     379           3 :             bHasMin = true;
     380           3 :             dfMin = CPLAtofM(papszTokens[1]);
     381             :         }
     382         252 :         else if (EQUAL(papszTokens[0], "Maximum data value"))
     383             :         {
     384           3 :             bHasMax = true;
     385           3 :             dfMax = CPLAtofM(papszTokens[1]);
     386             :         }
     387         249 :         else if (EQUAL(papszTokens[0], "Minimum x-axis coordinate"))
     388             :         {
     389           3 :             bHasMinX = true;
     390           3 :             dfMinX = CPLAtofM(papszTokens[1]);
     391             :         }
     392         246 :         else if (EQUAL(papszTokens[0], "Minimum y-axis coordinate"))
     393             :         {
     394           3 :             bHasMinY = true;
     395           3 :             dfMinY = CPLAtofM(papszTokens[1]);
     396             :         }
     397         243 :         else if (EQUAL(papszTokens[0], "Maximum x-axis coordinate"))
     398             :         {
     399           3 :             bHasMaxX = true;
     400           3 :             dfMaxX = CPLAtofM(papszTokens[1]);
     401             :         }
     402         240 :         else if (EQUAL(papszTokens[0], "Maximum y-axis coordinate"))
     403             :         {
     404           3 :             bHasMaxY = true;
     405           3 :             dfMaxY = CPLAtofM(papszTokens[1]);
     406             :         }
     407             : 
     408         318 :         CSLDestroy(papszTokens);
     409             :     }
     410             : 
     411           3 :     CPL_IGNORE_RET_VAL(VSIFCloseL(poOpenInfo->fpL));
     412           3 :     poOpenInfo->fpL = nullptr;
     413             : 
     414             :     /* -------------------------------------------------------------------- */
     415             :     /*      Did we get the required keywords?  If not we return with        */
     416             :     /*      this never having been considered to be a match. This isn't     */
     417             :     /*      an error!                                                       */
     418             :     /* -------------------------------------------------------------------- */
     419           3 :     if (nRows == -1 || nCols == -1 || !bIsInteger || !bIs2Bytes)
     420           0 :         return nullptr;
     421             : 
     422           3 :     if (!bNotProjected || !bIsWGS84)
     423           0 :         return nullptr;
     424             : 
     425           3 :     if (osDataFilename.empty())
     426           0 :         return nullptr;
     427             : 
     428           3 :     if (!GDALCheckDatasetDimensions(nCols, nRows))
     429           0 :         return nullptr;
     430             : 
     431             :     /* -------------------------------------------------------------------- */
     432             :     /*      Open target binary file.                                        */
     433             :     /* -------------------------------------------------------------------- */
     434           3 :     const char *pszPath = CPLGetPath(poOpenInfo->pszFilename);
     435           3 :     osDataFilename = CPLFormFilename(pszPath, osDataFilename, nullptr);
     436             : 
     437           3 :     VSILFILE *fpRaw = VSIFOpenL(osDataFilename, "rb");
     438             : 
     439           3 :     if (fpRaw == nullptr)
     440           0 :         return nullptr;
     441             : 
     442             :     /* -------------------------------------------------------------------- */
     443             :     /*      Create a corresponding GDALDataset.                             */
     444             :     /* -------------------------------------------------------------------- */
     445           6 :     auto poDS = std::make_unique<SNODASDataset>();
     446             : 
     447           3 :     poDS->nRasterXSize = nCols;
     448           3 :     poDS->nRasterYSize = nRows;
     449           3 :     poDS->osDataFilename = std::move(osDataFilename);
     450           3 :     poDS->bHasNoData = bHasNoData;
     451           3 :     poDS->dfNoData = dfNoData;
     452           3 :     poDS->bHasMin = bHasMin;
     453           3 :     poDS->dfMin = dfMin;
     454           3 :     poDS->bHasMax = bHasMax;
     455           3 :     poDS->dfMax = dfMax;
     456           3 :     if (bHasMinX && bHasMinY && bHasMaxX && bHasMaxY)
     457             :     {
     458           3 :         poDS->bGotTransform = true;
     459           3 :         poDS->adfGeoTransform[0] = dfMinX;
     460           3 :         poDS->adfGeoTransform[1] = (dfMaxX - dfMinX) / nCols;
     461           3 :         poDS->adfGeoTransform[2] = 0.0;
     462           3 :         poDS->adfGeoTransform[3] = dfMaxY;
     463           3 :         poDS->adfGeoTransform[4] = 0.0;
     464           3 :         poDS->adfGeoTransform[5] = -(dfMaxY - dfMinY) / nRows;
     465             :     }
     466             : 
     467           3 :     if (!osDescription.empty())
     468           3 :         poDS->SetMetadataItem("Description", osDescription);
     469           3 :     if (!osDataUnits.empty())
     470           3 :         poDS->SetMetadataItem("Data_Units", osDataUnits);
     471           3 :     if (nStartYear != -1 && nStartMonth != -1 && nStartDay != -1 &&
     472           3 :         nStartHour != -1 && nStartMinute != -1 && nStartSecond != -1)
     473           0 :         poDS->SetMetadataItem(
     474             :             "Start_Date",
     475             :             CPLSPrintf("%04d/%02d/%02d %02d:%02d:%02d", nStartYear, nStartMonth,
     476           0 :                        nStartDay, nStartHour, nStartMinute, nStartSecond));
     477           3 :     if (nStopYear != -1 && nStopMonth != -1 && nStopDay != -1 &&
     478           3 :         nStopHour != -1 && nStopMinute != -1 && nStopSecond != -1)
     479           6 :         poDS->SetMetadataItem("Stop_Date",
     480             :                               CPLSPrintf("%04d/%02d/%02d %02d:%02d:%02d",
     481             :                                          nStopYear, nStopMonth, nStopDay,
     482           3 :                                          nStopHour, nStopMinute, nStopSecond));
     483             : 
     484             :     /* -------------------------------------------------------------------- */
     485             :     /*      Create band information objects.                                */
     486             :     /* -------------------------------------------------------------------- */
     487           6 :     auto poBand = std::make_unique<SNODASRasterBand>(fpRaw, nCols, nRows);
     488           3 :     if (!poBand->IsValid())
     489           0 :         return nullptr;
     490           3 :     poDS->SetBand(1, std::move(poBand));
     491             : 
     492             :     /* -------------------------------------------------------------------- */
     493             :     /*      Initialize any PAM information.                                 */
     494             :     /* -------------------------------------------------------------------- */
     495           3 :     poDS->SetDescription(poOpenInfo->pszFilename);
     496           3 :     poDS->TryLoadXML();
     497             : 
     498             :     /* -------------------------------------------------------------------- */
     499             :     /*      Check for overviews.                                            */
     500             :     /* -------------------------------------------------------------------- */
     501           3 :     poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
     502             : 
     503           3 :     return poDS.release();
     504             : }
     505             : 
     506             : /************************************************************************/
     507             : /*                       GDALRegister_SNODAS()                          */
     508             : /************************************************************************/
     509             : 
     510        1595 : void GDALRegister_SNODAS()
     511             : 
     512             : {
     513        1595 :     if (GDALGetDriverByName("SNODAS") != nullptr)
     514         302 :         return;
     515             : 
     516        1293 :     GDALDriver *poDriver = new GDALDriver();
     517             : 
     518        1293 :     poDriver->SetDescription("SNODAS");
     519        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     520        1293 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
     521        1293 :                               "Snow Data Assimilation System");
     522        1293 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/snodas.html");
     523        1293 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "hdr");
     524        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     525             : 
     526        1293 :     poDriver->pfnOpen = SNODASDataset::Open;
     527        1293 :     poDriver->pfnIdentify = SNODASDataset::Identify;
     528             : 
     529        1293 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     530             : }