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

Source
/******************************************************************************
 *
 * 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 "gdal_priv.h"
#include "ogr_srs_api.h"
#include "rawdataset.h"

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

class SNODASRasterBand;

class SNODASDataset final : public RawDataset
{
    CPLString osDataFilename{};
    bool bGotTransform;
    GDALGeoTransform m_gt{};
    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(GDALProgressFunc = nullptr, void * = nullptr) override;

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

    CPLErr GetGeoTransform(GDALGeoTransform &gt) const 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 = cpl::down_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 = cpl::down_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 = cpl::down_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);
}

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

SNODASDataset::~SNODASDataset()

{
    SNODASDataset::Close();
}

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

CPLErr SNODASDataset::Close(GDALProgressFunc, void *)
{
    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(GDALGeoTransform &gt) const

{
    if (bGotTransform)
    {
        gt = m_gt;
        return CE_None;
    }

    return GDALPamDataset::GetGeoTransform(gt);
}

/************************************************************************/
/*                            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 (CPLHasPathTraversal(osDataFilename.c_str()))
    {
        CPLError(CE_Failure, CPLE_AppDefined, "Path traversal detected in %s",
                 osDataFilename.c_str());
        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->m_gt.xorig = dfMinX;
        poDS->m_gt.xscale = (dfMaxX - dfMinX) / nCols;
        poDS->m_gt.xrot = 0.0;
        poDS->m_gt.yorig = dfMaxY;
        poDS->m_gt.yrot = 0.0;
        poDS->m_gt.yscale = -(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: 2026-06-30 08:33

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