/src/gdal/frmts/raw/snodasdataset.cpp

Source (jump to first uncovered line)
/******************************************************************************
 *
 * Project:  SNODAS driver
 * Purpose:  Implementation of SNODASDataset
 * Author:   Even Rouault, <even dot rouault at spatialys.com>
 *
 ******************************************************************************
 * Copyright (c) 2011, Even Rouault <even dot rouault at spatialys.com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "cpl_string.h"
#include "gdal_frmts.h"
#include "ogr_srs_api.h"
#include "rawdataset.h"

/************************************************************************/
/* ==================================================================== */
/*                            SNODASDataset                             */
/* ==================================================================== */
/************************************************************************/

class SNODASRasterBand;

class SNODASDataset final : public RawDataset
{
    CPLString osDataFilename{};
    bool bGotTransform;
    double adfGeoTransform[6];
    bool bHasNoData;
    double dfNoData;
    bool bHasMin;
    double dfMin;
    int bHasMax;
    double dfMax;
    OGRSpatialReference m_oSRS{};

    friend class SNODASRasterBand;

    CPL_DISALLOW_COPY_ASSIGN(SNODASDataset)

    CPLErr Close() override;

  public:
    SNODASDataset();
    ~SNODASDataset() override;

    CPLErr GetGeoTransform(double *padfTransform) override;

    const OGRSpatialReference *GetSpatialRef() const override
    {
        return &m_oSRS;
    }

    char **GetFileList() override;

    static GDALDataset *Open(GDALOpenInfo *);
    static int Identify(GDALOpenInfo *);
};

/************************************************************************/
/* ==================================================================== */
/*                            SNODASRasterBand                          */
/* ==================================================================== */
/************************************************************************/

class SNODASRasterBand final : public RawRasterBand
{
    CPL_DISALLOW_COPY_ASSIGN(SNODASRasterBand)

  public:
    SNODASRasterBand(VSILFILE *fpRaw, int nXSize, int nYSize);

    ~SNODASRasterBand() override
    {
    }

    double GetNoDataValue(int *pbSuccess = nullptr) override;
    double GetMinimum(int *pbSuccess = nullptr) override;
    double GetMaximum(int *pbSuccess = nullptr) override;
};

/************************************************************************/
/*                         SNODASRasterBand()                           */
/************************************************************************/

SNODASRasterBand::SNODASRasterBand(VSILFILE *fpRawIn, int nXSize, int nYSize)
    : RawRasterBand(fpRawIn, 0, 2, nXSize * 2, GDT_Int16, !CPL_IS_LSB, nXSize,
                    nYSize, RawRasterBand::OwnFP::YES)
{
}

/************************************************************************/
/*                          GetNoDataValue()                            */
/************************************************************************/

double SNODASRasterBand::GetNoDataValue(int *pbSuccess)
{
    SNODASDataset *poGDS = reinterpret_cast<SNODASDataset *>(poDS);
    if (pbSuccess)
        *pbSuccess = poGDS->bHasNoData;

    if (poGDS->bHasNoData)
        return poGDS->dfNoData;

    return RawRasterBand::GetNoDataValue(pbSuccess);
}

/************************************************************************/
/*                            GetMinimum()                              */
/************************************************************************/

double SNODASRasterBand::GetMinimum(int *pbSuccess)
{
    SNODASDataset *poGDS = reinterpret_cast<SNODASDataset *>(poDS);
    if (pbSuccess)
        *pbSuccess = poGDS->bHasMin;

    if (poGDS->bHasMin)
        return poGDS->dfMin;

    return RawRasterBand::GetMinimum(pbSuccess);
}

/************************************************************************/
/*                            GetMaximum()                             */
/************************************************************************/

double SNODASRasterBand::GetMaximum(int *pbSuccess)
{
    SNODASDataset *poGDS = reinterpret_cast<SNODASDataset *>(poDS);
    if (pbSuccess)
        *pbSuccess = poGDS->bHasMax;

    if (poGDS->bHasMax)
        return poGDS->dfMax;

    return RawRasterBand::GetMaximum(pbSuccess);
}

/************************************************************************/
/* ==================================================================== */
/*                            SNODASDataset                             */
/* ==================================================================== */
/************************************************************************/

/************************************************************************/
/*                           SNODASDataset()                            */
/************************************************************************/

SNODASDataset::SNODASDataset()
    : bGotTransform(false), bHasNoData(false), dfNoData(0.0), bHasMin(false),
      dfMin(0.0), bHasMax(false), dfMax(0.0)
{
    m_oSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
    m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);

    adfGeoTransform[0] = 0.0;
    adfGeoTransform[1] = 1.0;
    adfGeoTransform[2] = 0.0;
    adfGeoTransform[3] = 0.0;
    adfGeoTransform[4] = 0.0;
    adfGeoTransform[5] = 1.0;
}

/************************************************************************/
/*                           ~SNODASDataset()                           */
/************************************************************************/

SNODASDataset::~SNODASDataset()

{
    SNODASDataset::Close();
}

/************************************************************************/
/*                              Close()                                 */
/************************************************************************/

CPLErr SNODASDataset::Close()
{
    CPLErr eErr = CE_None;
    if (nOpenFlags != OPEN_FLAGS_CLOSED)
    {
        if (SNODASDataset::FlushCache(true) != CE_None)
            eErr = CE_Failure;

        if (GDALPamDataset::Close() != CE_None)
            eErr = CE_Failure;
    }
    return eErr;
}

/************************************************************************/
/*                          GetGeoTransform()                           */
/************************************************************************/

CPLErr SNODASDataset::GetGeoTransform(double *padfTransform)

{
    if (bGotTransform)
    {
        memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
        return CE_None;
    }

    return GDALPamDataset::GetGeoTransform(padfTransform);
}

/************************************************************************/
/*                            GetFileList()                             */
/************************************************************************/

char **SNODASDataset::GetFileList()

{
    char **papszFileList = GDALPamDataset::GetFileList();

    papszFileList = CSLAddString(papszFileList, osDataFilename);

    return papszFileList;
}

/************************************************************************/
/*                            Identify()                                */
/************************************************************************/

int SNODASDataset::Identify(GDALOpenInfo *poOpenInfo)
{
    if (poOpenInfo->nHeaderBytes == 0)
        return FALSE;

    return STARTS_WITH_CI(reinterpret_cast<char *>(poOpenInfo->pabyHeader),
                          "Format version: NOHRSC GIS/RS raster file v1.1");
}

/************************************************************************/
/*                                Open()                                */
/************************************************************************/

GDALDataset *SNODASDataset::Open(GDALOpenInfo *poOpenInfo)

{
    if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
        return nullptr;

    /* -------------------------------------------------------------------- */
    /*      Confirm the requested access is supported.                      */
    /* -------------------------------------------------------------------- */
    if (poOpenInfo->eAccess == GA_Update)
    {
        ReportUpdateNotSupportedByDriver("SNODAS");
        return nullptr;
    }

    int nRows = -1;
    int nCols = -1;
    CPLString osDataFilename;
    bool bIsInteger = false;
    bool bIs2Bytes = false;
    double dfNoData = 0;
    bool bHasNoData = false;
    double dfMin = 0;
    bool bHasMin = false;
    double dfMax = 0;
    bool bHasMax = false;
    double dfMinX = 0.0;
    double dfMinY = 0.0;
    double dfMaxX = 0.0;
    double dfMaxY = 0.0;
    bool bHasMinX = false;
    bool bHasMinY = false;
    bool bHasMaxX = false;
    bool bHasMaxY = false;
    bool bNotProjected = false;
    bool bIsWGS84 = false;
    CPLString osDataUnits;
    CPLString osDescription;
    int nStartYear = -1;
    int nStartMonth = -1;
    int nStartDay = -1;
    int nStartHour = -1;
    int nStartMinute = -1;
    int nStartSecond = -1;
    int nStopYear = -1;
    int nStopMonth = -1;
    int nStopDay = -1;
    int nStopHour = -1;
    int nStopMinute = -1;
    int nStopSecond = -1;

    const char *pszLine = nullptr;
    while ((pszLine = CPLReadLine2L(poOpenInfo->fpL, 1024, nullptr)) != nullptr)
    {
        char **papszTokens =
            CSLTokenizeStringComplex(pszLine, ":", TRUE, FALSE);
        if (CSLCount(papszTokens) != 2)
        {
            CSLDestroy(papszTokens);
            continue;
        }
        if (papszTokens[1][0] == ' ')
            memmove(papszTokens[1], papszTokens[1] + 1,
                    strlen(papszTokens[1] + 1) + 1);

        if (EQUAL(papszTokens[0], "Data file pathname"))
        {
            osDataFilename = papszTokens[1];
        }
        else if (EQUAL(papszTokens[0], "Description"))
        {
            osDescription = papszTokens[1];
        }
        else if (EQUAL(papszTokens[0], "Data units"))
        {
            osDataUnits = papszTokens[1];
        }

        else if (EQUAL(papszTokens[0], "Start year"))
            nStartYear = atoi(papszTokens[1]);
        else if (EQUAL(papszTokens[0], "Start month"))
            nStartMonth = atoi(papszTokens[1]);
        else if (EQUAL(papszTokens[0], "Start day"))
            nStartDay = atoi(papszTokens[1]);
        else if (EQUAL(papszTokens[0], "Start hour"))
            nStartHour = atoi(papszTokens[1]);
        else if (EQUAL(papszTokens[0], " Start minute"))
            nStartMinute = atoi(papszTokens[1]);
        else if (EQUAL(papszTokens[0], "Start second"))
            nStartSecond = atoi(papszTokens[1]);

        else if (EQUAL(papszTokens[0], "Stop year"))
            nStopYear = atoi(papszTokens[1]);
        else if (EQUAL(papszTokens[0], "Stop month"))
            nStopMonth = atoi(papszTokens[1]);
        else if (EQUAL(papszTokens[0], "Stop day"))
            nStopDay = atoi(papszTokens[1]);
        else if (EQUAL(papszTokens[0], "Stop hour"))
            nStopHour = atoi(papszTokens[1]);
        else if (EQUAL(papszTokens[0], "Stop minute"))
            nStopMinute = atoi(papszTokens[1]);
        else if (EQUAL(papszTokens[0], "Stop second"))
            nStopSecond = atoi(papszTokens[1]);

        else if (EQUAL(papszTokens[0], "Number of columns"))
        {
            nCols = atoi(papszTokens[1]);
        }
        else if (EQUAL(papszTokens[0], "Number of rows"))
        {
            nRows = atoi(papszTokens[1]);
        }
        else if (EQUAL(papszTokens[0], "Data type"))
        {
            bIsInteger = EQUAL(papszTokens[1], "integer");
        }
        else if (EQUAL(papszTokens[0], "Data bytes per pixel"))
        {
            bIs2Bytes = EQUAL(papszTokens[1], "2");
        }
        else if (EQUAL(papszTokens[0], "Projected"))
        {
            bNotProjected = EQUAL(papszTokens[1], "no");
        }
        else if (EQUAL(papszTokens[0], "Horizontal datum"))
        {
            bIsWGS84 = EQUAL(papszTokens[1], "WGS84");
        }
        else if (EQUAL(papszTokens[0], "No data value"))
        {
            bHasNoData = true;
            dfNoData = CPLAtofM(papszTokens[1]);
        }
        else if (EQUAL(papszTokens[0], "Minimum data value"))
        {
            bHasMin = true;
            dfMin = CPLAtofM(papszTokens[1]);
        }
        else if (EQUAL(papszTokens[0], "Maximum data value"))
        {
            bHasMax = true;
            dfMax = CPLAtofM(papszTokens[1]);
        }
        else if (EQUAL(papszTokens[0], "Minimum x-axis coordinate"))
        {
            bHasMinX = true;
            dfMinX = CPLAtofM(papszTokens[1]);
        }
        else if (EQUAL(papszTokens[0], "Minimum y-axis coordinate"))
        {
            bHasMinY = true;
            dfMinY = CPLAtofM(papszTokens[1]);
        }
        else if (EQUAL(papszTokens[0], "Maximum x-axis coordinate"))
        {
            bHasMaxX = true;
            dfMaxX = CPLAtofM(papszTokens[1]);
        }
        else if (EQUAL(papszTokens[0], "Maximum y-axis coordinate"))
        {
            bHasMaxY = true;
            dfMaxY = CPLAtofM(papszTokens[1]);
        }

        CSLDestroy(papszTokens);
    }

    CPL_IGNORE_RET_VAL(VSIFCloseL(poOpenInfo->fpL));
    poOpenInfo->fpL = nullptr;

    /* -------------------------------------------------------------------- */
    /*      Did we get the required keywords?  If not we return with        */
    /*      this never having been considered to be a match. This isn't     */
    /*      an error!                                                       */
    /* -------------------------------------------------------------------- */
    if (nRows == -1 || nCols == -1 || !bIsInteger || !bIs2Bytes)
        return nullptr;

    if (!bNotProjected || !bIsWGS84)
        return nullptr;

    if (osDataFilename.empty())
        return nullptr;

    if (!GDALCheckDatasetDimensions(nCols, nRows))
        return nullptr;

    /* -------------------------------------------------------------------- */
    /*      Open target binary file.                                        */
    /* -------------------------------------------------------------------- */
    const std::string osPath = CPLGetPathSafe(poOpenInfo->pszFilename);
    osDataFilename =
        CPLFormFilenameSafe(osPath.c_str(), osDataFilename, nullptr);

    VSILFILE *fpRaw = VSIFOpenL(osDataFilename, "rb");

    if (fpRaw == nullptr)
        return nullptr;

    /* -------------------------------------------------------------------- */
    /*      Create a corresponding GDALDataset.                             */
    /* -------------------------------------------------------------------- */
    auto poDS = std::make_unique<SNODASDataset>();

    poDS->nRasterXSize = nCols;
    poDS->nRasterYSize = nRows;
    poDS->osDataFilename = std::move(osDataFilename);
    poDS->bHasNoData = bHasNoData;
    poDS->dfNoData = dfNoData;
    poDS->bHasMin = bHasMin;
    poDS->dfMin = dfMin;
    poDS->bHasMax = bHasMax;
    poDS->dfMax = dfMax;
    if (bHasMinX && bHasMinY && bHasMaxX && bHasMaxY)
    {
        poDS->bGotTransform = true;
        poDS->adfGeoTransform[0] = dfMinX;
        poDS->adfGeoTransform[1] = (dfMaxX - dfMinX) / nCols;
        poDS->adfGeoTransform[2] = 0.0;
        poDS->adfGeoTransform[3] = dfMaxY;
        poDS->adfGeoTransform[4] = 0.0;
        poDS->adfGeoTransform[5] = -(dfMaxY - dfMinY) / nRows;
    }

    if (!osDescription.empty())
        poDS->SetMetadataItem("Description", osDescription);
    if (!osDataUnits.empty())
        poDS->SetMetadataItem("Data_Units", osDataUnits);
    if (nStartYear != -1 && nStartMonth != -1 && nStartDay != -1 &&
        nStartHour != -1 && nStartMinute != -1 && nStartSecond != -1)
        poDS->SetMetadataItem(
            "Start_Date",
            CPLSPrintf("%04d/%02d/%02d %02d:%02d:%02d", nStartYear, nStartMonth,
                       nStartDay, nStartHour, nStartMinute, nStartSecond));
    if (nStopYear != -1 && nStopMonth != -1 && nStopDay != -1 &&
        nStopHour != -1 && nStopMinute != -1 && nStopSecond != -1)
        poDS->SetMetadataItem("Stop_Date",
                              CPLSPrintf("%04d/%02d/%02d %02d:%02d:%02d",
                                         nStopYear, nStopMonth, nStopDay,
                                         nStopHour, nStopMinute, nStopSecond));

    /* -------------------------------------------------------------------- */
    /*      Create band information objects.                                */
    /* -------------------------------------------------------------------- */
    auto poBand = std::make_unique<SNODASRasterBand>(fpRaw, nCols, nRows);
    if (!poBand->IsValid())
        return nullptr;
    poDS->SetBand(1, std::move(poBand));

    /* -------------------------------------------------------------------- */
    /*      Initialize any PAM information.                                 */
    /* -------------------------------------------------------------------- */
    poDS->SetDescription(poOpenInfo->pszFilename);
    poDS->TryLoadXML();

    /* -------------------------------------------------------------------- */
    /*      Check for overviews.                                            */
    /* -------------------------------------------------------------------- */
    poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);

    return poDS.release();
}

/************************************************************************/
/*                       GDALRegister_SNODAS()                          */
/************************************************************************/

void GDALRegister_SNODAS()

{
    if (GDALGetDriverByName("SNODAS") != nullptr)
        return;

    GDALDriver *poDriver = new GDALDriver();

    poDriver->SetDescription("SNODAS");
    poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
                              "Snow Data Assimilation System");
    poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/snodas.html");
    poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "hdr");
    poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");

    poDriver->pfnOpen = SNODASDataset::Open;
    poDriver->pfnIdentify = SNODASDataset::Identify;

    GetGDALDriverManager()->RegisterDriver(poDriver);
}

Coverage Report

Created: 2025-06-09 07:43

Line	Count	Source (jump to first uncovered line)
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		const OGRSpatialReference *GetSpatialRef() const override
52	0	{
53	0	return &m_oSRS;
54	0	}
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		~SNODASRasterBand() override
76	0	{
77	0	}
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		SNODASRasterBand::SNODASRasterBand(VSILFILE *fpRawIn, int nXSize, int nYSize)
89	0	: RawRasterBand(fpRawIn, 0, 2, nXSize * 2, GDT_Int16, !CPL_IS_LSB, nXSize,
90	0	nYSize, RawRasterBand::OwnFP::YES)
91	0	{
92	0	}
93
94		/************************************************************************/
95		/* GetNoDataValue() */
96		/************************************************************************/
97
98		double SNODASRasterBand::GetNoDataValue(int *pbSuccess)
99	0	{
100	0	SNODASDataset poGDS = reinterpret_cast<SNODASDataset >(poDS);
101	0	if (pbSuccess)
102	0	*pbSuccess = poGDS->bHasNoData;
103
104	0	if (poGDS->bHasNoData)
105	0	return poGDS->dfNoData;
106
107	0	return RawRasterBand::GetNoDataValue(pbSuccess);
108	0	}
109
110		/************************************************************************/
111		/* GetMinimum() */
112		/************************************************************************/
113
114		double SNODASRasterBand::GetMinimum(int *pbSuccess)
115	0	{
116	0	SNODASDataset poGDS = reinterpret_cast<SNODASDataset >(poDS);
117	0	if (pbSuccess)
118	0	*pbSuccess = poGDS->bHasMin;
119
120	0	if (poGDS->bHasMin)
121	0	return poGDS->dfMin;
122
123	0	return RawRasterBand::GetMinimum(pbSuccess);
124	0	}
125
126		/************************************************************************/
127		/* GetMaximum() */
128		/************************************************************************/
129
130		double SNODASRasterBand::GetMaximum(int *pbSuccess)
131	0	{
132	0	SNODASDataset poGDS = reinterpret_cast<SNODASDataset >(poDS);
133	0	if (pbSuccess)
134	0	*pbSuccess = poGDS->bHasMax;
135
136	0	if (poGDS->bHasMax)
137	0	return poGDS->dfMax;
138
139	0	return RawRasterBand::GetMaximum(pbSuccess);
140	0	}
141
142		/************************************************************************/
143		/* ==================================================================== */
144		/* SNODASDataset */
145		/* ==================================================================== */
146		/************************************************************************/
147
148		/************************************************************************/
149		/* SNODASDataset() */
150		/************************************************************************/
151
152		SNODASDataset::SNODASDataset()
153	0	: bGotTransform(false), bHasNoData(false), dfNoData(0.0), bHasMin(false),
154	0	dfMin(0.0), bHasMax(false), dfMax(0.0)
155	0	{
156	0	m_oSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
157	0	m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
158
159	0	adfGeoTransform[0] = 0.0;
160	0	adfGeoTransform[1] = 1.0;
161	0	adfGeoTransform[2] = 0.0;
162	0	adfGeoTransform[3] = 0.0;
163	0	adfGeoTransform[4] = 0.0;
164	0	adfGeoTransform[5] = 1.0;
165	0	}
166
167		/************************************************************************/
168		/* ~SNODASDataset() */
169		/************************************************************************/
170
171		SNODASDataset::~SNODASDataset()
172
173	0	{
174	0	SNODASDataset::Close();
175	0	}
176
177		/************************************************************************/
178		/* Close() */
179		/************************************************************************/
180
181		CPLErr SNODASDataset::Close()
182	0	{
183	0	CPLErr eErr = CE_None;
184	0	if (nOpenFlags != OPEN_FLAGS_CLOSED)
185	0	{
186	0	if (SNODASDataset::FlushCache(true) != CE_None)
187	0	eErr = CE_Failure;
188
189	0	if (GDALPamDataset::Close() != CE_None)
190	0	eErr = CE_Failure;
191	0	}
192	0	return eErr;
193	0	}
194
195		/************************************************************************/
196		/* GetGeoTransform() */
197		/************************************************************************/
198
199		CPLErr SNODASDataset::GetGeoTransform(double *padfTransform)
200
201	0	{
202	0	if (bGotTransform)
203	0	{
204	0	memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
205	0	return CE_None;
206	0	}
207
208	0	return GDALPamDataset::GetGeoTransform(padfTransform);
209	0	}
210
211		/************************************************************************/
212		/* GetFileList() */
213		/************************************************************************/
214
215		char **SNODASDataset::GetFileList()
216
217	0	{
218	0	char **papszFileList = GDALPamDataset::GetFileList();
219
220	0	papszFileList = CSLAddString(papszFileList, osDataFilename);
221
222	0	return papszFileList;
223	0	}
224
225		/************************************************************************/
226		/* Identify() */
227		/************************************************************************/
228
229		int SNODASDataset::Identify(GDALOpenInfo *poOpenInfo)
230	177k	{
231	177k	if (poOpenInfo->nHeaderBytes == 0)
232	173k	return FALSE;
233
234	3.57k	return STARTS_WITH_CI(reinterpret_cast<char *>(poOpenInfo->pabyHeader),
235	177k	"Format version: NOHRSC GIS/RS raster file v1.1");
236	177k	}
237
238		/************************************************************************/
239		/* Open() */
240		/************************************************************************/
241
242		GDALDataset SNODASDataset::Open(GDALOpenInfo poOpenInfo)
243
244	80	{
245	80	if (!Identify(poOpenInfo) \|\| poOpenInfo->fpL == nullptr)
246	0	return nullptr;
247
248		/* -------------------------------------------------------------------- */
249		/* Confirm the requested access is supported. */
250		/* -------------------------------------------------------------------- */
251	80	if (poOpenInfo->eAccess == GA_Update)
252	0	{
253	0	ReportUpdateNotSupportedByDriver("SNODAS");
254	0	return nullptr;
255	0	}
256
257	80	int nRows = -1;
258	80	int nCols = -1;
259	80	CPLString osDataFilename;
260	80	bool bIsInteger = false;
261	80	bool bIs2Bytes = false;
262	80	double dfNoData = 0;
263	80	bool bHasNoData = false;
264	80	double dfMin = 0;
265	80	bool bHasMin = false;
266	80	double dfMax = 0;
267	80	bool bHasMax = false;
268	80	double dfMinX = 0.0;
269	80	double dfMinY = 0.0;
270	80	double dfMaxX = 0.0;
271	80	double dfMaxY = 0.0;
272	80	bool bHasMinX = false;
273	80	bool bHasMinY = false;
274	80	bool bHasMaxX = false;
275	80	bool bHasMaxY = false;
276	80	bool bNotProjected = false;
277	80	bool bIsWGS84 = false;
278	80	CPLString osDataUnits;
279	80	CPLString osDescription;
280	80	int nStartYear = -1;
281	80	int nStartMonth = -1;
282	80	int nStartDay = -1;
283	80	int nStartHour = -1;
284	80	int nStartMinute = -1;
285	80	int nStartSecond = -1;
286	80	int nStopYear = -1;
287	80	int nStopMonth = -1;
288	80	int nStopDay = -1;
289	80	int nStopHour = -1;
290	80	int nStopMinute = -1;
291	80	int nStopSecond = -1;
292
293	80	const char *pszLine = nullptr;
294	49.6k	while ((pszLine = CPLReadLine2L(poOpenInfo->fpL, 1024, nullptr)) != nullptr)
295	49.5k	{
296	49.5k	char **papszTokens =
297	49.5k	CSLTokenizeStringComplex(pszLine, ":", TRUE, FALSE);
298	49.5k	if (CSLCount(papszTokens) != 2)
299	11.3k	{
300	11.3k	CSLDestroy(papszTokens);
301	11.3k	continue;
302	11.3k	}
303	38.1k	if (papszTokens[1][0] == ' ')
304	36.9k	memmove(papszTokens[1], papszTokens[1] + 1,
305	36.9k	strlen(papszTokens[1] + 1) + 1);
306
307	38.1k	if (EQUAL(papszTokens[0], "Data file pathname"))
308	844	{
309	844	osDataFilename = papszTokens[1];
310	844	}
311	37.3k	else if (EQUAL(papszTokens[0], "Description"))
312	714	{
313	714	osDescription = papszTokens[1];
314	714	}
315	36.6k	else if (EQUAL(papszTokens[0], "Data units"))
316	611	{
317	611	osDataUnits = papszTokens[1];
318	611	}
319
320	35.9k	else if (EQUAL(papszTokens[0], "Start year"))
321	279	nStartYear = atoi(papszTokens[1]);
322	35.7k	else if (EQUAL(papszTokens[0], "Start month"))
323	280	nStartMonth = atoi(papszTokens[1]);
324	35.4k	else if (EQUAL(papszTokens[0], "Start day"))
325	398	nStartDay = atoi(papszTokens[1]);
326	35.0k	else if (EQUAL(papszTokens[0], "Start hour"))
327	338	nStartHour = atoi(papszTokens[1]);
328	34.7k	else if (EQUAL(papszTokens[0], " Start minute"))
329	0	nStartMinute = atoi(papszTokens[1]);
330	34.7k	else if (EQUAL(papszTokens[0], "Start second"))
331	161	nStartSecond = atoi(papszTokens[1]);
332
333	34.5k	else if (EQUAL(papszTokens[0], "Stop year"))
334	103	nStopYear = atoi(papszTokens[1]);
335	34.4k	else if (EQUAL(papszTokens[0], "Stop month"))
336	181	nStopMonth = atoi(papszTokens[1]);
337	34.2k	else if (EQUAL(papszTokens[0], "Stop day"))
338	511	nStopDay = atoi(papszTokens[1]);
339	33.7k	else if (EQUAL(papszTokens[0], "Stop hour"))
340	355	nStopHour = atoi(papszTokens[1]);
341	33.3k	else if (EQUAL(papszTokens[0], "Stop minute"))
342	617	nStopMinute = atoi(papszTokens[1]);
343	32.7k	else if (EQUAL(papszTokens[0], "Stop second"))
344	614	nStopSecond = atoi(papszTokens[1]);
345
346	32.1k	else if (EQUAL(papszTokens[0], "Number of columns"))
347	762	{
348	762	nCols = atoi(papszTokens[1]);
349	762	}
350	31.4k	else if (EQUAL(papszTokens[0], "Number of rows"))
351	441	{
352	441	nRows = atoi(papszTokens[1]);
353	441	}
354	30.9k	else if (EQUAL(papszTokens[0], "Data type"))
355	1.18k	{
356	1.18k	bIsInteger = EQUAL(papszTokens[1], "integer");
357	1.18k	}
358	29.7k	else if (EQUAL(papszTokens[0], "Data bytes per pixel"))
359	935	{
360	935	bIs2Bytes = EQUAL(papszTokens[1], "2");
361	935	}
362	28.8k	else if (EQUAL(papszTokens[0], "Projected"))
363	398	{
364	398	bNotProjected = EQUAL(papszTokens[1], "no");
365	398	}
366	28.4k	else if (EQUAL(papszTokens[0], "Horizontal datum"))
367	393	{
368	393	bIsWGS84 = EQUAL(papszTokens[1], "WGS84");
369	393	}
370	28.0k	else if (EQUAL(papszTokens[0], "No data value"))
371	748	{
372	748	bHasNoData = true;
373	748	dfNoData = CPLAtofM(papszTokens[1]);
374	748	}
375	27.2k	else if (EQUAL(papszTokens[0], "Minimum data value"))
376	716	{
377	716	bHasMin = true;
378	716	dfMin = CPLAtofM(papszTokens[1]);
379	716	}
380	26.5k	else if (EQUAL(papszTokens[0], "Maximum data value"))
381	726	{
382	726	bHasMax = true;
383	726	dfMax = CPLAtofM(papszTokens[1]);
384	726	}
385	25.8k	else if (EQUAL(papszTokens[0], "Minimum x-axis coordinate"))
386	120	{
387	120	bHasMinX = true;
388	120	dfMinX = CPLAtofM(papszTokens[1]);
389	120	}
390	25.7k	else if (EQUAL(papszTokens[0], "Minimum y-axis coordinate"))
391	236	{
392	236	bHasMinY = true;
393	236	dfMinY = CPLAtofM(papszTokens[1]);
394	236	}
395	25.4k	else if (EQUAL(papszTokens[0], "Maximum x-axis coordinate"))
396	134	{
397	134	bHasMaxX = true;
398	134	dfMaxX = CPLAtofM(papszTokens[1]);
399	134	}
400	25.3k	else if (EQUAL(papszTokens[0], "Maximum y-axis coordinate"))
401	123	{
402	123	bHasMaxY = true;
403	123	dfMaxY = CPLAtofM(papszTokens[1]);
404	123	}
405
406	38.1k	CSLDestroy(papszTokens);
407	38.1k	}
408
409	80	CPL_IGNORE_RET_VAL(VSIFCloseL(poOpenInfo->fpL));
410	80	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	80	if (nRows == -1 \|\| nCols == -1 \|\| !bIsInteger \|\| !bIs2Bytes)
418	51	return nullptr;
419
420	29	if (!bNotProjected \|\| !bIsWGS84)
421	14	return nullptr;
422
423	15	if (osDataFilename.empty())
424	0	return nullptr;
425
426	15	if (!GDALCheckDatasetDimensions(nCols, nRows))
427	0	return nullptr;
428
429		/* -------------------------------------------------------------------- */
430		/* Open target binary file. */
431		/* -------------------------------------------------------------------- */
432	15	const std::string osPath = CPLGetPathSafe(poOpenInfo->pszFilename);
433	15	osDataFilename =
434	15	CPLFormFilenameSafe(osPath.c_str(), osDataFilename, nullptr);
435
436	15	VSILFILE *fpRaw = VSIFOpenL(osDataFilename, "rb");
437
438	15	if (fpRaw == nullptr)
439	15	return nullptr;
440
441		/* -------------------------------------------------------------------- */
442		/* Create a corresponding GDALDataset. */
443		/* -------------------------------------------------------------------- */
444	0	auto poDS = std::make_unique<SNODASDataset>();
445
446	0	poDS->nRasterXSize = nCols;
447	0	poDS->nRasterYSize = nRows;
448	0	poDS->osDataFilename = std::move(osDataFilename);
449	0	poDS->bHasNoData = bHasNoData;
450	0	poDS->dfNoData = dfNoData;
451	0	poDS->bHasMin = bHasMin;
452	0	poDS->dfMin = dfMin;
453	0	poDS->bHasMax = bHasMax;
454	0	poDS->dfMax = dfMax;
455	0	if (bHasMinX && bHasMinY && bHasMaxX && bHasMaxY)
456	0	{
457	0	poDS->bGotTransform = true;
458	0	poDS->adfGeoTransform[0] = dfMinX;
459	0	poDS->adfGeoTransform[1] = (dfMaxX - dfMinX) / nCols;
460	0	poDS->adfGeoTransform[2] = 0.0;
461	0	poDS->adfGeoTransform[3] = dfMaxY;
462	0	poDS->adfGeoTransform[4] = 0.0;
463	0	poDS->adfGeoTransform[5] = -(dfMaxY - dfMinY) / nRows;
464	0	}
465
466	0	if (!osDescription.empty())
467	0	poDS->SetMetadataItem("Description", osDescription);
468	0	if (!osDataUnits.empty())
469	0	poDS->SetMetadataItem("Data_Units", osDataUnits);
470	0	if (nStartYear != -1 && nStartMonth != -1 && nStartDay != -1 &&
471	0	nStartHour != -1 && nStartMinute != -1 && nStartSecond != -1)
472	0	poDS->SetMetadataItem(
473	0	"Start_Date",
474	0	CPLSPrintf("%04d/%02d/%02d %02d:%02d:%02d", nStartYear, nStartMonth,
475	0	nStartDay, nStartHour, nStartMinute, nStartSecond));
476	0	if (nStopYear != -1 && nStopMonth != -1 && nStopDay != -1 &&
477	0	nStopHour != -1 && nStopMinute != -1 && nStopSecond != -1)
478	0	poDS->SetMetadataItem("Stop_Date",
479	0	CPLSPrintf("%04d/%02d/%02d %02d:%02d:%02d",
480	0	nStopYear, nStopMonth, nStopDay,
481	0	nStopHour, nStopMinute, nStopSecond));
482
483		/* -------------------------------------------------------------------- */
484		/* Create band information objects. */
485		/* -------------------------------------------------------------------- */
486	0	auto poBand = std::make_unique<SNODASRasterBand>(fpRaw, nCols, nRows);
487	0	if (!poBand->IsValid())
488	0	return nullptr;
489	0	poDS->SetBand(1, std::move(poBand));
490
491		/* -------------------------------------------------------------------- */
492		/* Initialize any PAM information. */
493		/* -------------------------------------------------------------------- */
494	0	poDS->SetDescription(poOpenInfo->pszFilename);
495	0	poDS->TryLoadXML();
496
497		/* -------------------------------------------------------------------- */
498		/* Check for overviews. */
499		/* -------------------------------------------------------------------- */
500	0	poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
501
502	0	return poDS.release();
503	0	}
504
505		/************************************************************************/
506		/* GDALRegister_SNODAS() */
507		/************************************************************************/
508
509		void GDALRegister_SNODAS()
510
511	2	{
512	2	if (GDALGetDriverByName("SNODAS") != nullptr)
513	0	return;
514
515	2	GDALDriver *poDriver = new GDALDriver();
516
517	2	poDriver->SetDescription("SNODAS");
518	2	poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
519	2	poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
520	2	"Snow Data Assimilation System");
521	2	poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/snodas.html");
522	2	poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "hdr");
523	2	poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
524
525	2	poDriver->pfnOpen = SNODASDataset::Open;
526	2	poDriver->pfnIdentify = SNODASDataset::Identify;
527
528	2	GetGDALDriverManager()->RegisterDriver(poDriver);
529	2	}