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