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