1 : /******************************************************************************
2 : * $Id: snodasdataset.cpp 21984 2011-03-19 17:22:58Z rouault $
3 : *
4 : * Project: SNODAS driver
5 : * Purpose: Implementation of SNODASDataset
6 : * Author: Even Rouault, <even dot rouault at mines dash paris dot org>
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2011, Even Rouault
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "rawdataset.h"
31 : #include "ogr_srs_api.h"
32 : #include "cpl_string.h"
33 :
34 : CPL_CVSID("$Id: snodasdataset.cpp 21984 2011-03-19 17:22:58Z rouault $");
35 :
36 : // g++ -g -Wall -fPIC frmts/raw/snodasdataset.cpp -shared -o gdal_SNODAS.so -Iport -Igcore -Ifrmts/raw -Iogr -L. -lgdal
37 :
38 : CPL_C_START
39 : void GDALRegister_SNODAS(void);
40 : CPL_C_END
41 :
42 : /************************************************************************/
43 : /* ==================================================================== */
44 : /* SNODASDataset */
45 : /* ==================================================================== */
46 : /************************************************************************/
47 :
48 : class SNODASRasterBand;
49 :
50 : class SNODASDataset : public RawDataset
51 : {
52 : CPLString osDataFilename;
53 : int bGotTransform;
54 : double adfGeoTransform[6];
55 : int bHasNoData;
56 : double dfNoData;
57 : int bHasMin;
58 : double dfMin;
59 : int bHasMax;
60 : double dfMax;
61 :
62 : friend class SNODASRasterBand;
63 :
64 : public:
65 : SNODASDataset();
66 : ~SNODASDataset();
67 :
68 : virtual CPLErr GetGeoTransform( double * padfTransform );
69 : virtual const char *GetProjectionRef(void);
70 :
71 : virtual char **GetFileList();
72 :
73 : static GDALDataset *Open( GDALOpenInfo * );
74 : static int Identify( GDALOpenInfo * );
75 : };
76 :
77 :
78 : /************************************************************************/
79 : /* ==================================================================== */
80 : /* SNODASRasterBand */
81 : /* ==================================================================== */
82 : /************************************************************************/
83 :
84 : class SNODASRasterBand : public RawRasterBand
85 4 : {
86 : public:
87 : SNODASRasterBand(VSILFILE* fpRaw, int nXSize, int nYSize);
88 :
89 : virtual double GetNoDataValue( int *pbSuccess = NULL );
90 : virtual double GetMinimum( int *pbSuccess = NULL );
91 : virtual double GetMaximum(int *pbSuccess = NULL );
92 : };
93 :
94 :
95 : /************************************************************************/
96 : /* SNODASRasterBand() */
97 : /************************************************************************/
98 :
99 4 : SNODASRasterBand::SNODASRasterBand(VSILFILE* fpRaw,
100 : int nXSize, int nYSize) :
101 : RawRasterBand( fpRaw, 0, 2,
102 : nXSize * 2, GDT_Int16,
103 4 : !CPL_IS_LSB, nXSize, nYSize, TRUE, TRUE)
104 : {
105 4 : }
106 :
107 : /************************************************************************/
108 : /* GetNoDataValue() */
109 : /************************************************************************/
110 :
111 2 : double SNODASRasterBand::GetNoDataValue( int *pbSuccess )
112 : {
113 2 : SNODASDataset* poGDS = (SNODASDataset*) poDS;
114 2 : if (pbSuccess)
115 2 : *pbSuccess = poGDS->bHasNoData;
116 2 : if (poGDS->bHasNoData)
117 2 : return poGDS->dfNoData;
118 : else
119 0 : return RawRasterBand::GetNoDataValue(pbSuccess);
120 : }
121 :
122 : /************************************************************************/
123 : /* GetMinimum() */
124 : /************************************************************************/
125 :
126 2 : double SNODASRasterBand::GetMinimum( int *pbSuccess )
127 : {
128 2 : SNODASDataset* poGDS = (SNODASDataset*) poDS;
129 2 : if (pbSuccess)
130 2 : *pbSuccess = poGDS->bHasMin;
131 2 : if (poGDS->bHasMin)
132 2 : return poGDS->dfMin;
133 : else
134 0 : return RawRasterBand::GetMinimum(pbSuccess);
135 : }
136 :
137 : /************************************************************************/
138 : /* GetMaximum() */
139 : /************************************************************************/
140 :
141 2 : double SNODASRasterBand::GetMaximum( int *pbSuccess )
142 : {
143 2 : SNODASDataset* poGDS = (SNODASDataset*) poDS;
144 2 : if (pbSuccess)
145 2 : *pbSuccess = poGDS->bHasMax;
146 2 : if (poGDS->bHasMax)
147 2 : return poGDS->dfMax;
148 : else
149 0 : return RawRasterBand::GetMaximum(pbSuccess);
150 : }
151 :
152 : /************************************************************************/
153 : /* ==================================================================== */
154 : /* SNODASDataset */
155 : /* ==================================================================== */
156 : /************************************************************************/
157 :
158 : /************************************************************************/
159 : /* SNODASDataset() */
160 : /************************************************************************/
161 :
162 4 : SNODASDataset::SNODASDataset()
163 : {
164 4 : bGotTransform = FALSE;
165 4 : adfGeoTransform[0] = 0.0;
166 4 : adfGeoTransform[1] = 1.0;
167 4 : adfGeoTransform[2] = 0.0;
168 4 : adfGeoTransform[3] = 0.0;
169 4 : adfGeoTransform[4] = 0.0;
170 4 : adfGeoTransform[5] = 1.0;
171 4 : bHasNoData = FALSE;
172 4 : dfNoData = 0.0;
173 4 : bHasMin = FALSE;
174 4 : dfMin = 0.0;
175 4 : bHasMax = FALSE;
176 4 : dfMax = 0.0;
177 4 : }
178 :
179 : /************************************************************************/
180 : /* ~SNODASDataset() */
181 : /************************************************************************/
182 :
183 4 : SNODASDataset::~SNODASDataset()
184 :
185 : {
186 4 : FlushCache();
187 4 : }
188 :
189 : /************************************************************************/
190 : /* GetProjectionRef() */
191 : /************************************************************************/
192 :
193 2 : const char *SNODASDataset::GetProjectionRef()
194 :
195 : {
196 2 : return SRS_WKT_WGS84;
197 : }
198 :
199 : /************************************************************************/
200 : /* GetGeoTransform() */
201 : /************************************************************************/
202 :
203 2 : CPLErr SNODASDataset::GetGeoTransform( double * padfTransform )
204 :
205 : {
206 2 : if( bGotTransform )
207 : {
208 2 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
209 2 : return CE_None;
210 : }
211 : else
212 : {
213 0 : return GDALPamDataset::GetGeoTransform( padfTransform );
214 : }
215 : }
216 :
217 :
218 : /************************************************************************/
219 : /* GetFileList() */
220 : /************************************************************************/
221 :
222 2 : char **SNODASDataset::GetFileList()
223 :
224 : {
225 2 : char **papszFileList = GDALPamDataset::GetFileList();
226 :
227 2 : papszFileList = CSLAddString(papszFileList, osDataFilename);
228 :
229 2 : return papszFileList;
230 : }
231 :
232 : /************************************************************************/
233 : /* Identify() */
234 : /************************************************************************/
235 :
236 22424 : int SNODASDataset::Identify( GDALOpenInfo * poOpenInfo )
237 : {
238 22424 : if (poOpenInfo->nHeaderBytes == 0)
239 21220 : return FALSE;
240 :
241 1204 : return EQUALN((const char*)poOpenInfo->pabyHeader,
242 : "Format version: NOHRSC GIS/RS raster file v1.1",
243 : strlen("Format version: NOHRSC GIS/RS raster file v1.1"));
244 : }
245 :
246 : /************************************************************************/
247 : /* Open() */
248 : /************************************************************************/
249 :
250 3514 : GDALDataset *SNODASDataset::Open( GDALOpenInfo * poOpenInfo )
251 :
252 : {
253 3514 : if( !Identify(poOpenInfo) )
254 3510 : return NULL;
255 :
256 : VSILFILE *fp;
257 :
258 4 : fp = VSIFOpenL( poOpenInfo->pszFilename, "r" );
259 :
260 4 : if( fp == NULL )
261 : {
262 0 : return NULL;
263 : }
264 :
265 : const char * pszLine;
266 4 : int nRows = -1, nCols = -1;
267 4 : CPLString osDataFilename;
268 4 : int bIsInteger = FALSE, bIs2Bytes = FALSE;
269 4 : double dfNoData = 0;
270 4 : int bHasNoData = FALSE;
271 4 : double dfMin = 0;
272 4 : int bHasMin = FALSE;
273 4 : double dfMax = 0;
274 4 : int bHasMax = FALSE;
275 4 : double dfMinX = 0.0, dfMinY = 0.0, dfMaxX = 0.0, dfMaxY = 0.0;
276 4 : int bHasMinX = FALSE, bHasMinY = FALSE, bHasMaxX = FALSE, bHasMaxY = FALSE;
277 4 : int bNotProjected = FALSE, bIsWGS84 = FALSE;
278 4 : CPLString osDescription, osDataUnits;
279 4 : int nStartYear = -1, nStartMonth = -1, nStartDay = -1,
280 4 : nStartHour = -1, nStartMinute = -1, nStartSecond = -1;
281 4 : int nStopYear = -1, nStopMonth = -1, nStopDay = -1,
282 4 : nStopHour = -1, nStopMinute = -1, nStopSecond = -1;
283 :
284 436 : while( (pszLine = CPLReadLine2L( fp, 256, NULL )) != NULL )
285 : {
286 428 : char** papszTokens = CSLTokenizeStringComplex( pszLine, ":", TRUE, FALSE );
287 428 : if( CSLCount( papszTokens ) != 2 )
288 : {
289 4 : CSLDestroy( papszTokens );
290 4 : continue;
291 : }
292 424 : if( papszTokens[1][0] == ' ' )
293 424 : memmove(papszTokens[1], papszTokens[1] + 1, strlen(papszTokens[1] + 1) + 1);
294 :
295 424 : if( EQUAL(papszTokens[0],"Data file pathname") )
296 : {
297 4 : osDataFilename = papszTokens[1];
298 : }
299 420 : else if( EQUAL(papszTokens[0],"Description") )
300 : {
301 4 : osDescription = papszTokens[1];
302 : }
303 416 : else if( EQUAL(papszTokens[0],"Data units") )
304 : {
305 4 : osDataUnits= papszTokens[1];
306 : }
307 :
308 412 : else if( EQUAL(papszTokens[0],"Start year") )
309 4 : nStartYear = atoi(papszTokens[1]);
310 408 : else if( EQUAL(papszTokens[0],"Start month") )
311 4 : nStartMonth = atoi(papszTokens[1]);
312 404 : else if( EQUAL(papszTokens[0],"Start day") )
313 4 : nStartDay = atoi(papszTokens[1]);
314 400 : else if( EQUAL(papszTokens[0],"Start hour") )
315 4 : nStartHour = atoi(papszTokens[1]);
316 396 : else if( EQUAL(papszTokens[0],"Start minute") )
317 4 : nStartMinute = atoi(papszTokens[1]);
318 392 : else if( EQUAL(papszTokens[0],"Start second") )
319 4 : nStartSecond = atoi(papszTokens[1]);
320 :
321 388 : else if( EQUAL(papszTokens[0],"Stop year") )
322 4 : nStopYear = atoi(papszTokens[1]);
323 384 : else if( EQUAL(papszTokens[0],"Stop month") )
324 4 : nStopMonth = atoi(papszTokens[1]);
325 380 : else if( EQUAL(papszTokens[0],"Stop day") )
326 4 : nStopDay = atoi(papszTokens[1]);
327 376 : else if( EQUAL(papszTokens[0],"Stop hour") )
328 4 : nStopHour = atoi(papszTokens[1]);
329 372 : else if( EQUAL(papszTokens[0],"Stop minute") )
330 4 : nStopMinute = atoi(papszTokens[1]);
331 368 : else if( EQUAL(papszTokens[0],"Stop second") )
332 4 : nStopSecond = atoi(papszTokens[1]);
333 :
334 364 : else if( EQUAL(papszTokens[0],"Number of columns") )
335 : {
336 4 : nCols = atoi(papszTokens[1]);
337 : }
338 360 : else if( EQUAL(papszTokens[0],"Number of rows") )
339 : {
340 4 : nRows = atoi(papszTokens[1]);
341 : }
342 356 : else if( EQUAL(papszTokens[0],"Data type"))
343 : {
344 4 : bIsInteger = EQUAL(papszTokens[1],"integer");
345 : }
346 352 : else if( EQUAL(papszTokens[0],"Data bytes per pixel"))
347 : {
348 4 : bIs2Bytes = EQUAL(papszTokens[1],"2");
349 : }
350 348 : else if( EQUAL(papszTokens[0],"Projected"))
351 : {
352 4 : bNotProjected = EQUAL(papszTokens[1],"no");
353 : }
354 344 : else if( EQUAL(papszTokens[0],"Horizontal datum"))
355 : {
356 4 : bIsWGS84 = EQUAL(papszTokens[1],"WGS84");
357 : }
358 340 : else if( EQUAL(papszTokens[0],"No data value"))
359 : {
360 4 : bHasNoData = TRUE;
361 4 : dfNoData = CPLAtofM(papszTokens[1]);
362 : }
363 336 : else if( EQUAL(papszTokens[0],"Minimum data value"))
364 : {
365 4 : bHasMin = TRUE;
366 4 : dfMin = CPLAtofM(papszTokens[1]);
367 : }
368 332 : else if( EQUAL(papszTokens[0],"Maximum data value"))
369 : {
370 4 : bHasMax = TRUE;
371 4 : dfMax = CPLAtofM(papszTokens[1]);
372 : }
373 328 : else if( EQUAL(papszTokens[0],"Minimum x-axis coordinate") )
374 : {
375 4 : bHasMinX = TRUE;
376 4 : dfMinX = CPLAtofM(papszTokens[1]);
377 : }
378 324 : else if( EQUAL(papszTokens[0],"Minimum y-axis coordinate") )
379 : {
380 4 : bHasMinY = TRUE;
381 4 : dfMinY = CPLAtofM(papszTokens[1]);
382 : }
383 320 : else if( EQUAL(papszTokens[0],"Maximum x-axis coordinate") )
384 : {
385 4 : bHasMaxX = TRUE;
386 4 : dfMaxX = CPLAtofM(papszTokens[1]);
387 : }
388 316 : else if( EQUAL(papszTokens[0],"Maximum y-axis coordinate") )
389 : {
390 4 : bHasMaxY = TRUE;
391 4 : dfMaxY = CPLAtofM(papszTokens[1]);
392 : }
393 :
394 424 : CSLDestroy( papszTokens );
395 : }
396 :
397 4 : VSIFCloseL( fp );
398 :
399 : /* -------------------------------------------------------------------- */
400 : /* Did we get the required keywords? If not we return with */
401 : /* this never having been considered to be a match. This isn't */
402 : /* an error! */
403 : /* -------------------------------------------------------------------- */
404 4 : if( nRows == -1 || nCols == -1 || !bIsInteger || !bIs2Bytes )
405 0 : return NULL;
406 :
407 4 : if( !bNotProjected || !bIsWGS84 )
408 0 : return NULL;
409 :
410 4 : if( osDataFilename.size() == 0 )
411 0 : return NULL;
412 :
413 4 : if (!GDALCheckDatasetDimensions(nCols, nRows))
414 0 : return NULL;
415 :
416 : /* -------------------------------------------------------------------- */
417 : /* Open target binary file. */
418 : /* -------------------------------------------------------------------- */
419 4 : const char* pszPath = CPLGetPath(poOpenInfo->pszFilename);
420 4 : osDataFilename = CPLFormFilename(pszPath, osDataFilename, NULL);
421 :
422 4 : VSILFILE* fpRaw = VSIFOpenL( osDataFilename, "rb" );
423 :
424 4 : if( fpRaw == NULL )
425 0 : return NULL;
426 :
427 : /* -------------------------------------------------------------------- */
428 : /* Create a corresponding GDALDataset. */
429 : /* -------------------------------------------------------------------- */
430 : SNODASDataset *poDS;
431 :
432 4 : poDS = new SNODASDataset();
433 :
434 4 : poDS->nRasterXSize = nCols;
435 4 : poDS->nRasterYSize = nRows;
436 4 : poDS->osDataFilename = osDataFilename;
437 4 : poDS->bHasNoData = bHasNoData;
438 4 : poDS->dfNoData = dfNoData;
439 4 : poDS->bHasMin = bHasMin;
440 4 : poDS->dfMin = dfMin;
441 4 : poDS->bHasMax = bHasMax;
442 4 : poDS->dfMax = dfMax;
443 4 : if (bHasMinX && bHasMinY && bHasMaxX && bHasMaxY)
444 : {
445 4 : poDS->bGotTransform = TRUE;
446 4 : poDS->adfGeoTransform[0] = dfMinX;
447 4 : poDS->adfGeoTransform[1] = (dfMaxX - dfMinX) / nCols;
448 4 : poDS->adfGeoTransform[2] = 0.0;
449 4 : poDS->adfGeoTransform[3] = dfMaxY;
450 4 : poDS->adfGeoTransform[4] = 0.0;
451 4 : poDS->adfGeoTransform[5] = - (dfMaxY - dfMinY) / nRows;
452 : }
453 :
454 4 : if (osDescription.size())
455 4 : poDS->SetMetadataItem("Description", osDescription);
456 4 : if (osDataUnits.size())
457 4 : poDS->SetMetadataItem("Data_Units", osDataUnits);
458 4 : if (nStartYear != -1 && nStartMonth != -1 && nStartDay != -1 &&
459 : nStartHour != -1 && nStartMinute != -1 && nStartSecond != -1)
460 : poDS->SetMetadataItem("Start_Date",
461 : CPLSPrintf("%04d/%02d/%02d %02d:%02d:%02d",
462 : nStartYear, nStartMonth, nStartDay,
463 4 : nStartHour, nStartMinute, nStartSecond));
464 4 : if (nStopYear != -1 && nStopMonth != -1 && nStopDay != -1 &&
465 : nStopHour != -1 && nStopMinute != -1 && nStopSecond != -1)
466 : poDS->SetMetadataItem("Stop_Date",
467 : CPLSPrintf("%04d/%02d/%02d %02d:%02d:%02d",
468 : nStopYear, nStopMonth, nStopDay,
469 4 : nStopHour, nStopMinute, nStopSecond));
470 :
471 : /* -------------------------------------------------------------------- */
472 : /* Create band information objects. */
473 : /* -------------------------------------------------------------------- */
474 4 : poDS->SetBand( 1, new SNODASRasterBand( fpRaw, nCols, nRows) );
475 :
476 : /* -------------------------------------------------------------------- */
477 : /* Initialize any PAM information. */
478 : /* -------------------------------------------------------------------- */
479 4 : poDS->SetDescription( poOpenInfo->pszFilename );
480 4 : poDS->TryLoadXML();
481 :
482 : /* -------------------------------------------------------------------- */
483 : /* Check for overviews. */
484 : /* -------------------------------------------------------------------- */
485 4 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
486 :
487 4 : return( poDS );
488 : }
489 :
490 : /************************************************************************/
491 : /* GDALRegister_SNODAS() */
492 : /************************************************************************/
493 :
494 1135 : void GDALRegister_SNODAS()
495 :
496 : {
497 : GDALDriver *poDriver;
498 :
499 1135 : if( GDALGetDriverByName( "SNODAS" ) == NULL )
500 : {
501 1093 : poDriver = new GDALDriver();
502 :
503 1093 : poDriver->SetDescription( "SNODAS" );
504 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
505 1093 : "Snow Data Assimilation System" );
506 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
507 1093 : "frmt_various.html#SNODAS" );
508 1093 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "hdr" );
509 :
510 1093 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
511 1093 : poDriver->pfnOpen = SNODASDataset::Open;
512 1093 : poDriver->pfnIdentify = SNODASDataset::Identify;
513 :
514 1093 : GetGDALDriverManager()->RegisterDriver( poDriver );
515 : }
516 1135 : }
|
