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       52417 : int SNODASDataset::Identify(GDALOpenInfo *poOpenInfo)
     230             : {
     231       52417 :     if (poOpenInfo->nHeaderBytes == 0)
     232       48726 :         return FALSE;
     233             : 
     234        3691 :     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 :         ReportUpdateNotSupportedByDriver("SNODAS");
     254           0 :         return nullptr;
     255             :     }
     256             : 
     257           3 :     int nRows = -1;
     258           3 :     int nCols = -1;
     259           6 :     CPLString osDataFilename;
     260           3 :     bool bIsInteger = false;
     261           3 :     bool bIs2Bytes = false;
     262           3 :     double dfNoData = 0;
     263           3 :     bool bHasNoData = false;
     264           3 :     double dfMin = 0;
     265           3 :     bool bHasMin = false;
     266           3 :     double dfMax = 0;
     267           3 :     bool bHasMax = false;
     268           3 :     double dfMinX = 0.0;
     269           3 :     double dfMinY = 0.0;
     270           3 :     double dfMaxX = 0.0;
     271           3 :     double dfMaxY = 0.0;
     272           3 :     bool bHasMinX = false;
     273           3 :     bool bHasMinY = false;
     274           3 :     bool bHasMaxX = false;
     275           3 :     bool bHasMaxY = false;
     276           3 :     bool bNotProjected = false;
     277           3 :     bool bIsWGS84 = false;
     278           6 :     CPLString osDataUnits;
     279           6 :     CPLString osDescription;
     280           3 :     int nStartYear = -1;
     281           3 :     int nStartMonth = -1;
     282           3 :     int nStartDay = -1;
     283           3 :     int nStartHour = -1;
     284           3 :     int nStartMinute = -1;
     285           3 :     int nStartSecond = -1;
     286           3 :     int nStopYear = -1;
     287           3 :     int nStopMonth = -1;
     288           3 :     int nStopDay = -1;
     289           3 :     int nStopHour = -1;
     290           3 :     int nStopMinute = -1;
     291           3 :     int nStopSecond = -1;
     292             : 
     293           3 :     const char *pszLine = nullptr;
     294         324 :     while ((pszLine = CPLReadLine2L(poOpenInfo->fpL, 1024, nullptr)) != nullptr)
     295             :     {
     296             :         char **papszTokens =
     297         321 :             CSLTokenizeStringComplex(pszLine, ":", TRUE, FALSE);
     298         321 :         if (CSLCount(papszTokens) != 2)
     299             :         {
     300           3 :             CSLDestroy(papszTokens);
     301           3 :             continue;
     302             :         }
     303         318 :         if (papszTokens[1][0] == ' ')
     304         318 :             memmove(papszTokens[1], papszTokens[1] + 1,
     305         318 :                     strlen(papszTokens[1] + 1) + 1);
     306             : 
     307         318 :         if (EQUAL(papszTokens[0], "Data file pathname"))
     308             :         {
     309           3 :             osDataFilename = papszTokens[1];
     310             :         }
     311         315 :         else if (EQUAL(papszTokens[0], "Description"))
     312             :         {
     313           3 :             osDescription = papszTokens[1];
     314             :         }
     315         312 :         else if (EQUAL(papszTokens[0], "Data units"))
     316             :         {
     317           3 :             osDataUnits = papszTokens[1];
     318             :         }
     319             : 
     320         309 :         else if (EQUAL(papszTokens[0], "Start year"))
     321           3 :             nStartYear = atoi(papszTokens[1]);
     322         306 :         else if (EQUAL(papszTokens[0], "Start month"))
     323           3 :             nStartMonth = atoi(papszTokens[1]);
     324         303 :         else if (EQUAL(papszTokens[0], "Start day"))
     325           3 :             nStartDay = atoi(papszTokens[1]);
     326         300 :         else if (EQUAL(papszTokens[0], "Start hour"))
     327           3 :             nStartHour = atoi(papszTokens[1]);
     328         297 :         else if (EQUAL(papszTokens[0], " Start minute"))
     329           0 :             nStartMinute = atoi(papszTokens[1]);
     330         297 :         else if (EQUAL(papszTokens[0], "Start second"))
     331           3 :             nStartSecond = atoi(papszTokens[1]);
     332             : 
     333         294 :         else if (EQUAL(papszTokens[0], "Stop year"))
     334           3 :             nStopYear = atoi(papszTokens[1]);
     335         291 :         else if (EQUAL(papszTokens[0], "Stop month"))
     336           3 :             nStopMonth = atoi(papszTokens[1]);
     337         288 :         else if (EQUAL(papszTokens[0], "Stop day"))
     338           3 :             nStopDay = atoi(papszTokens[1]);
     339         285 :         else if (EQUAL(papszTokens[0], "Stop hour"))
     340           3 :             nStopHour = atoi(papszTokens[1]);
     341         282 :         else if (EQUAL(papszTokens[0], "Stop minute"))
     342           3 :             nStopMinute = atoi(papszTokens[1]);
     343         279 :         else if (EQUAL(papszTokens[0], "Stop second"))
     344           3 :             nStopSecond = atoi(papszTokens[1]);
     345             : 
     346         276 :         else if (EQUAL(papszTokens[0], "Number of columns"))
     347             :         {
     348           3 :             nCols = atoi(papszTokens[1]);
     349             :         }
     350         273 :         else if (EQUAL(papszTokens[0], "Number of rows"))
     351             :         {
     352           3 :             nRows = atoi(papszTokens[1]);
     353             :         }
     354         270 :         else if (EQUAL(papszTokens[0], "Data type"))
     355             :         {
     356           3 :             bIsInteger = EQUAL(papszTokens[1], "integer");
     357             :         }
     358         267 :         else if (EQUAL(papszTokens[0], "Data bytes per pixel"))
     359             :         {
     360           3 :             bIs2Bytes = EQUAL(papszTokens[1], "2");
     361             :         }
     362         264 :         else if (EQUAL(papszTokens[0], "Projected"))
     363             :         {
     364           3 :             bNotProjected = EQUAL(papszTokens[1], "no");
     365             :         }
     366         261 :         else if (EQUAL(papszTokens[0], "Horizontal datum"))
     367             :         {
     368           3 :             bIsWGS84 = EQUAL(papszTokens[1], "WGS84");
     369             :         }
     370         258 :         else if (EQUAL(papszTokens[0], "No data value"))
     371             :         {
     372           3 :             bHasNoData = true;
     373           3 :             dfNoData = CPLAtofM(papszTokens[1]);
     374             :         }
     375         255 :         else if (EQUAL(papszTokens[0], "Minimum data value"))
     376             :         {
     377           3 :             bHasMin = true;
     378           3 :             dfMin = CPLAtofM(papszTokens[1]);
     379             :         }
     380         252 :         else if (EQUAL(papszTokens[0], "Maximum data value"))
     381             :         {
     382           3 :             bHasMax = true;
     383           3 :             dfMax = CPLAtofM(papszTokens[1]);
     384             :         }
     385         249 :         else if (EQUAL(papszTokens[0], "Minimum x-axis coordinate"))
     386             :         {
     387           3 :             bHasMinX = true;
     388           3 :             dfMinX = CPLAtofM(papszTokens[1]);
     389             :         }
     390         246 :         else if (EQUAL(papszTokens[0], "Minimum y-axis coordinate"))
     391             :         {
     392           3 :             bHasMinY = true;
     393           3 :             dfMinY = CPLAtofM(papszTokens[1]);
     394             :         }
     395         243 :         else if (EQUAL(papszTokens[0], "Maximum x-axis coordinate"))
     396             :         {
     397           3 :             bHasMaxX = true;
     398           3 :             dfMaxX = CPLAtofM(papszTokens[1]);
     399             :         }
     400         240 :         else if (EQUAL(papszTokens[0], "Maximum y-axis coordinate"))
     401             :         {
     402           3 :             bHasMaxY = true;
     403           3 :             dfMaxY = CPLAtofM(papszTokens[1]);
     404             :         }
     405             : 
     406         318 :         CSLDestroy(papszTokens);
     407             :     }
     408             : 
     409           3 :     CPL_IGNORE_RET_VAL(VSIFCloseL(poOpenInfo->fpL));
     410           3 :     poOpenInfo->fpL = nullptr;
     411             : 
     412             :     /* -------------------------------------------------------------------- */
     413             :     /*      Did we get the required keywords?  If not we return with        */
     414             :     /*      this never having been considered to be a match. This isn't     */
     415             :     /*      an error!                                                       */
     416             :     /* -------------------------------------------------------------------- */
     417           3 :     if (nRows == -1 || nCols == -1 || !bIsInteger || !bIs2Bytes)
     418           0 :         return nullptr;
     419             : 
     420           3 :     if (!bNotProjected || !bIsWGS84)
     421           0 :         return nullptr;
     422             : 
     423           3 :     if (osDataFilename.empty())
     424           0 :         return nullptr;
     425             : 
     426           3 :     if (!GDALCheckDatasetDimensions(nCols, nRows))
     427           0 :         return nullptr;
     428             : 
     429             :     /* -------------------------------------------------------------------- */
     430             :     /*      Open target binary file.                                        */
     431             :     /* -------------------------------------------------------------------- */
     432           6 :     const std::string osPath = CPLGetPathSafe(poOpenInfo->pszFilename);
     433             :     osDataFilename =
     434           3 :         CPLFormFilenameSafe(osPath.c_str(), osDataFilename, nullptr);
     435             : 
     436           3 :     VSILFILE *fpRaw = VSIFOpenL(osDataFilename, "rb");
     437             : 
     438           3 :     if (fpRaw == nullptr)
     439           0 :         return nullptr;
     440             : 
     441             :     /* -------------------------------------------------------------------- */
     442             :     /*      Create a corresponding GDALDataset.                             */
     443             :     /* -------------------------------------------------------------------- */
     444           6 :     auto poDS = std::make_unique<SNODASDataset>();
     445             : 
     446           3 :     poDS->nRasterXSize = nCols;
     447           3 :     poDS->nRasterYSize = nRows;
     448           3 :     poDS->osDataFilename = std::move(osDataFilename);
     449           3 :     poDS->bHasNoData = bHasNoData;
     450           3 :     poDS->dfNoData = dfNoData;
     451           3 :     poDS->bHasMin = bHasMin;
     452           3 :     poDS->dfMin = dfMin;
     453           3 :     poDS->bHasMax = bHasMax;
     454           3 :     poDS->dfMax = dfMax;
     455           3 :     if (bHasMinX && bHasMinY && bHasMaxX && bHasMaxY)
     456             :     {
     457           3 :         poDS->bGotTransform = true;
     458           3 :         poDS->adfGeoTransform[0] = dfMinX;
     459           3 :         poDS->adfGeoTransform[1] = (dfMaxX - dfMinX) / nCols;
     460           3 :         poDS->adfGeoTransform[2] = 0.0;
     461           3 :         poDS->adfGeoTransform[3] = dfMaxY;
     462           3 :         poDS->adfGeoTransform[4] = 0.0;
     463           3 :         poDS->adfGeoTransform[5] = -(dfMaxY - dfMinY) / nRows;
     464             :     }
     465             : 
     466           3 :     if (!osDescription.empty())
     467           3 :         poDS->SetMetadataItem("Description", osDescription);
     468           3 :     if (!osDataUnits.empty())
     469           3 :         poDS->SetMetadataItem("Data_Units", osDataUnits);
     470           3 :     if (nStartYear != -1 && nStartMonth != -1 && nStartDay != -1 &&
     471           3 :         nStartHour != -1 && nStartMinute != -1 && nStartSecond != -1)
     472           0 :         poDS->SetMetadataItem(
     473             :             "Start_Date",
     474             :             CPLSPrintf("%04d/%02d/%02d %02d:%02d:%02d", nStartYear, nStartMonth,
     475           0 :                        nStartDay, nStartHour, nStartMinute, nStartSecond));
     476           3 :     if (nStopYear != -1 && nStopMonth != -1 && nStopDay != -1 &&
     477           3 :         nStopHour != -1 && nStopMinute != -1 && nStopSecond != -1)
     478           6 :         poDS->SetMetadataItem("Stop_Date",
     479             :                               CPLSPrintf("%04d/%02d/%02d %02d:%02d:%02d",
     480             :                                          nStopYear, nStopMonth, nStopDay,
     481           3 :                                          nStopHour, nStopMinute, nStopSecond));
     482             : 
     483             :     /* -------------------------------------------------------------------- */
     484             :     /*      Create band information objects.                                */
     485             :     /* -------------------------------------------------------------------- */
     486           6 :     auto poBand = std::make_unique<SNODASRasterBand>(fpRaw, nCols, nRows);
     487           3 :     if (!poBand->IsValid())
     488           0 :         return nullptr;
     489           3 :     poDS->SetBand(1, std::move(poBand));
     490             : 
     491             :     /* -------------------------------------------------------------------- */
     492             :     /*      Initialize any PAM information.                                 */
     493             :     /* -------------------------------------------------------------------- */
     494           3 :     poDS->SetDescription(poOpenInfo->pszFilename);
     495           3 :     poDS->TryLoadXML();
     496             : 
     497             :     /* -------------------------------------------------------------------- */
     498             :     /*      Check for overviews.                                            */
     499             :     /* -------------------------------------------------------------------- */
     500           3 :     poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
     501             : 
     502           3 :     return poDS.release();
     503             : }
     504             : 
     505             : /************************************************************************/
     506             : /*                       GDALRegister_SNODAS()                          */
     507             : /************************************************************************/
     508             : 
     509        1686 : void GDALRegister_SNODAS()
     510             : 
     511             : {
     512        1686 :     if (GDALGetDriverByName("SNODAS") != nullptr)
     513         302 :         return;
     514             : 
     515        1384 :     GDALDriver *poDriver = new GDALDriver();
     516             : 
     517        1384 :     poDriver->SetDescription("SNODAS");
     518        1384 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     519        1384 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
     520        1384 :                               "Snow Data Assimilation System");
     521        1384 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/snodas.html");
     522        1384 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "hdr");
     523        1384 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     524             : 
     525        1384 :     poDriver->pfnOpen = SNODASDataset::Open;
     526        1384 :     poDriver->pfnIdentify = SNODASDataset::Identify;
     527             : 
     528        1384 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     529             : }