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 2 : {
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 2 : SNODASRasterBand::SNODASRasterBand(VSILFILE* fpRaw,
100 : int nXSize, int nYSize) :
101 : RawRasterBand( fpRaw, 0, 2,
102 : nXSize * 2, GDT_Int16,
103 2 : !CPL_IS_LSB, nXSize, nYSize, TRUE, TRUE)
104 : {
105 2 : }
106 :
107 : /************************************************************************/
108 : /* GetNoDataValue() */
109 : /************************************************************************/
110 :
111 1 : double SNODASRasterBand::GetNoDataValue( int *pbSuccess )
112 : {
113 1 : SNODASDataset* poGDS = (SNODASDataset*) poDS;
114 1 : if (pbSuccess)
115 1 : *pbSuccess = poGDS->bHasNoData;
116 1 : if (poGDS->bHasNoData)
117 1 : return poGDS->dfNoData;
118 : else
119 0 : return RawRasterBand::GetNoDataValue(pbSuccess);
120 : }
121 :
122 : /************************************************************************/
123 : /* GetMinimum() */
124 : /************************************************************************/
125 :
126 1 : double SNODASRasterBand::GetMinimum( int *pbSuccess )
127 : {
128 1 : SNODASDataset* poGDS = (SNODASDataset*) poDS;
129 1 : if (pbSuccess)
130 1 : *pbSuccess = poGDS->bHasMin;
131 1 : if (poGDS->bHasMin)
132 1 : return poGDS->dfMin;
133 : else
134 0 : return RawRasterBand::GetMinimum(pbSuccess);
135 : }
136 :
137 : /************************************************************************/
138 : /* GetMaximum() */
139 : /************************************************************************/
140 :
141 1 : double SNODASRasterBand::GetMaximum( int *pbSuccess )
142 : {
143 1 : SNODASDataset* poGDS = (SNODASDataset*) poDS;
144 1 : if (pbSuccess)
145 1 : *pbSuccess = poGDS->bHasMax;
146 1 : if (poGDS->bHasMax)
147 1 : 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 2 : SNODASDataset::SNODASDataset()
163 : {
164 2 : bGotTransform = FALSE;
165 2 : adfGeoTransform[0] = 0.0;
166 2 : adfGeoTransform[1] = 1.0;
167 2 : adfGeoTransform[2] = 0.0;
168 2 : adfGeoTransform[3] = 0.0;
169 2 : adfGeoTransform[4] = 0.0;
170 2 : adfGeoTransform[5] = 1.0;
171 2 : bHasNoData = FALSE;
172 2 : dfNoData = 0.0;
173 2 : bHasMin = FALSE;
174 2 : dfMin = 0.0;
175 2 : bHasMax = FALSE;
176 2 : dfMax = 0.0;
177 2 : }
178 :
179 : /************************************************************************/
180 : /* ~SNODASDataset() */
181 : /************************************************************************/
182 :
183 2 : SNODASDataset::~SNODASDataset()
184 :
185 : {
186 2 : FlushCache();
187 2 : }
188 :
189 : /************************************************************************/
190 : /* GetProjectionRef() */
191 : /************************************************************************/
192 :
193 1 : const char *SNODASDataset::GetProjectionRef()
194 :
195 : {
196 1 : return SRS_WKT_WGS84;
197 : }
198 :
199 : /************************************************************************/
200 : /* GetGeoTransform() */
201 : /************************************************************************/
202 :
203 1 : CPLErr SNODASDataset::GetGeoTransform( double * padfTransform )
204 :
205 : {
206 1 : if( bGotTransform )
207 : {
208 1 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
209 1 : return CE_None;
210 : }
211 : else
212 : {
213 0 : return GDALPamDataset::GetGeoTransform( padfTransform );
214 : }
215 : }
216 :
217 :
218 : /************************************************************************/
219 : /* GetFileList() */
220 : /************************************************************************/
221 :
222 1 : char **SNODASDataset::GetFileList()
223 :
224 : {
225 1 : char **papszFileList = GDALPamDataset::GetFileList();
226 :
227 1 : papszFileList = CSLAddString(papszFileList, osDataFilename);
228 :
229 1 : return papszFileList;
230 : }
231 :
232 : /************************************************************************/
233 : /* Identify() */
234 : /************************************************************************/
235 :
236 11922 : int SNODASDataset::Identify( GDALOpenInfo * poOpenInfo )
237 : {
238 11922 : if (poOpenInfo->nHeaderBytes == 0)
239 11278 : return FALSE;
240 :
241 644 : 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 1829 : GDALDataset *SNODASDataset::Open( GDALOpenInfo * poOpenInfo )
251 :
252 : {
253 1829 : if( !Identify(poOpenInfo) )
254 1827 : return NULL;
255 :
256 : VSILFILE *fp;
257 :
258 2 : fp = VSIFOpenL( poOpenInfo->pszFilename, "r" );
259 :
260 2 : if( fp == NULL )
261 : {
262 0 : return NULL;
263 : }
264 :
265 : const char * pszLine;
266 2 : int nRows = -1, nCols = -1;
267 2 : CPLString osDataFilename;
268 2 : int bIsInteger = FALSE, bIs2Bytes = FALSE;
269 2 : double dfNoData = 0;
270 2 : int bHasNoData = FALSE;
271 2 : double dfMin = 0;
272 2 : int bHasMin = FALSE;
273 2 : double dfMax = 0;
274 2 : int bHasMax = FALSE;
275 2 : double dfMinX = 0.0, dfMinY = 0.0, dfMaxX = 0.0, dfMaxY = 0.0;
276 2 : int bHasMinX = FALSE, bHasMinY = FALSE, bHasMaxX = FALSE, bHasMaxY = FALSE;
277 2 : int bNotProjected = FALSE, bIsWGS84 = FALSE;
278 2 : CPLString osDescription, osDataUnits;
279 2 : int nStartYear = -1, nStartMonth = -1, nStartDay = -1,
280 2 : nStartHour = -1, nStartMinute = -1, nStartSecond = -1;
281 2 : int nStopYear = -1, nStopMonth = -1, nStopDay = -1,
282 2 : nStopHour = -1, nStopMinute = -1, nStopSecond = -1;
283 :
284 218 : while( (pszLine = CPLReadLine2L( fp, 256, NULL )) != NULL )
285 : {
286 214 : char** papszTokens = CSLTokenizeStringComplex( pszLine, ":", TRUE, FALSE );
287 214 : if( CSLCount( papszTokens ) != 2 )
288 : {
289 2 : CSLDestroy( papszTokens );
290 2 : continue;
291 : }
292 212 : if( papszTokens[1][0] == ' ' )
293 212 : memmove(papszTokens[1], papszTokens[1] + 1, strlen(papszTokens[1] + 1) + 1);
294 :
295 212 : if( EQUAL(papszTokens[0],"Data file pathname") )
296 : {
297 2 : osDataFilename = papszTokens[1];
298 : }
299 210 : else if( EQUAL(papszTokens[0],"Description") )
300 : {
301 2 : osDescription = papszTokens[1];
302 : }
303 208 : else if( EQUAL(papszTokens[0],"Data units") )
304 : {
305 2 : osDataUnits= papszTokens[1];
306 : }
307 :
308 206 : else if( EQUAL(papszTokens[0],"Start year") )
309 2 : nStartYear = atoi(papszTokens[1]);
310 204 : else if( EQUAL(papszTokens[0],"Start month") )
311 2 : nStartMonth = atoi(papszTokens[1]);
312 202 : else if( EQUAL(papszTokens[0],"Start day") )
313 2 : nStartDay = atoi(papszTokens[1]);
314 200 : else if( EQUAL(papszTokens[0],"Start hour") )
315 2 : nStartHour = atoi(papszTokens[1]);
316 198 : else if( EQUAL(papszTokens[0],"Start minute") )
317 2 : nStartMinute = atoi(papszTokens[1]);
318 196 : else if( EQUAL(papszTokens[0],"Start second") )
319 2 : nStartSecond = atoi(papszTokens[1]);
320 :
321 194 : else if( EQUAL(papszTokens[0],"Stop year") )
322 2 : nStopYear = atoi(papszTokens[1]);
323 192 : else if( EQUAL(papszTokens[0],"Stop month") )
324 2 : nStopMonth = atoi(papszTokens[1]);
325 190 : else if( EQUAL(papszTokens[0],"Stop day") )
326 2 : nStopDay = atoi(papszTokens[1]);
327 188 : else if( EQUAL(papszTokens[0],"Stop hour") )
328 2 : nStopHour = atoi(papszTokens[1]);
329 186 : else if( EQUAL(papszTokens[0],"Stop minute") )
330 2 : nStopMinute = atoi(papszTokens[1]);
331 184 : else if( EQUAL(papszTokens[0],"Stop second") )
332 2 : nStopSecond = atoi(papszTokens[1]);
333 :
334 182 : else if( EQUAL(papszTokens[0],"Number of columns") )
335 : {
336 2 : nCols = atoi(papszTokens[1]);
337 : }
338 180 : else if( EQUAL(papszTokens[0],"Number of rows") )
339 : {
340 2 : nRows = atoi(papszTokens[1]);
341 : }
342 178 : else if( EQUAL(papszTokens[0],"Data type"))
343 : {
344 2 : bIsInteger = EQUAL(papszTokens[1],"integer");
345 : }
346 176 : else if( EQUAL(papszTokens[0],"Data bytes per pixel"))
347 : {
348 2 : bIs2Bytes = EQUAL(papszTokens[1],"2");
349 : }
350 174 : else if( EQUAL(papszTokens[0],"Projected"))
351 : {
352 2 : bNotProjected = EQUAL(papszTokens[1],"no");
353 : }
354 172 : else if( EQUAL(papszTokens[0],"Horizontal datum"))
355 : {
356 2 : bIsWGS84 = EQUAL(papszTokens[1],"WGS84");
357 : }
358 170 : else if( EQUAL(papszTokens[0],"No data value"))
359 : {
360 2 : bHasNoData = TRUE;
361 2 : dfNoData = CPLAtofM(papszTokens[1]);
362 : }
363 168 : else if( EQUAL(papszTokens[0],"Minimum data value"))
364 : {
365 2 : bHasMin = TRUE;
366 2 : dfMin = CPLAtofM(papszTokens[1]);
367 : }
368 166 : else if( EQUAL(papszTokens[0],"Maximum data value"))
369 : {
370 2 : bHasMax = TRUE;
371 2 : dfMax = CPLAtofM(papszTokens[1]);
372 : }
373 164 : else if( EQUAL(papszTokens[0],"Minimum x-axis coordinate") )
374 : {
375 2 : bHasMinX = TRUE;
376 2 : dfMinX = CPLAtofM(papszTokens[1]);
377 : }
378 162 : else if( EQUAL(papszTokens[0],"Minimum y-axis coordinate") )
379 : {
380 2 : bHasMinY = TRUE;
381 2 : dfMinY = CPLAtofM(papszTokens[1]);
382 : }
383 160 : else if( EQUAL(papszTokens[0],"Maximum x-axis coordinate") )
384 : {
385 2 : bHasMaxX = TRUE;
386 2 : dfMaxX = CPLAtofM(papszTokens[1]);
387 : }
388 158 : else if( EQUAL(papszTokens[0],"Maximum y-axis coordinate") )
389 : {
390 2 : bHasMaxY = TRUE;
391 2 : dfMaxY = CPLAtofM(papszTokens[1]);
392 : }
393 :
394 212 : CSLDestroy( papszTokens );
395 : }
396 :
397 2 : 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 2 : if( nRows == -1 || nCols == -1 || !bIsInteger || !bIs2Bytes )
405 0 : return NULL;
406 :
407 2 : if( !bNotProjected || !bIsWGS84 )
408 0 : return NULL;
409 :
410 2 : if( osDataFilename.size() == 0 )
411 0 : return NULL;
412 :
413 2 : if (!GDALCheckDatasetDimensions(nCols, nRows))
414 0 : return NULL;
415 :
416 : /* -------------------------------------------------------------------- */
417 : /* Open target binary file. */
418 : /* -------------------------------------------------------------------- */
419 2 : const char* pszPath = CPLGetPath(poOpenInfo->pszFilename);
420 2 : osDataFilename = CPLFormFilename(pszPath, osDataFilename, NULL);
421 :
422 2 : VSILFILE* fpRaw = VSIFOpenL( osDataFilename, "rb" );
423 :
424 2 : if( fpRaw == NULL )
425 0 : return NULL;
426 :
427 : /* -------------------------------------------------------------------- */
428 : /* Create a corresponding GDALDataset. */
429 : /* -------------------------------------------------------------------- */
430 : SNODASDataset *poDS;
431 :
432 2 : poDS = new SNODASDataset();
433 :
434 2 : poDS->nRasterXSize = nCols;
435 2 : poDS->nRasterYSize = nRows;
436 2 : poDS->osDataFilename = osDataFilename;
437 2 : poDS->bHasNoData = bHasNoData;
438 2 : poDS->dfNoData = dfNoData;
439 2 : poDS->bHasMin = bHasMin;
440 2 : poDS->dfMin = dfMin;
441 2 : poDS->bHasMax = bHasMax;
442 2 : poDS->dfMax = dfMax;
443 2 : if (bHasMinX && bHasMinY && bHasMaxX && bHasMaxY)
444 : {
445 2 : poDS->bGotTransform = TRUE;
446 2 : poDS->adfGeoTransform[0] = dfMinX;
447 2 : poDS->adfGeoTransform[1] = (dfMaxX - dfMinX) / nCols;
448 2 : poDS->adfGeoTransform[2] = 0.0;
449 2 : poDS->adfGeoTransform[3] = dfMaxY;
450 2 : poDS->adfGeoTransform[4] = 0.0;
451 2 : poDS->adfGeoTransform[5] = - (dfMaxY - dfMinY) / nRows;
452 : }
453 :
454 2 : if (osDescription.size())
455 2 : poDS->SetMetadataItem("Description", osDescription);
456 2 : if (osDataUnits.size())
457 2 : poDS->SetMetadataItem("Data_Units", osDataUnits);
458 2 : 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 2 : nStartHour, nStartMinute, nStartSecond));
464 2 : 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 2 : nStopHour, nStopMinute, nStopSecond));
470 :
471 : /* -------------------------------------------------------------------- */
472 : /* Create band information objects. */
473 : /* -------------------------------------------------------------------- */
474 2 : poDS->SetBand( 1, new SNODASRasterBand( fpRaw, nCols, nRows) );
475 :
476 : /* -------------------------------------------------------------------- */
477 : /* Initialize any PAM information. */
478 : /* -------------------------------------------------------------------- */
479 2 : poDS->SetDescription( poOpenInfo->pszFilename );
480 2 : poDS->TryLoadXML();
481 :
482 : /* -------------------------------------------------------------------- */
483 : /* Check for overviews. */
484 : /* -------------------------------------------------------------------- */
485 2 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
486 :
487 2 : return( poDS );
488 : }
489 :
490 : /************************************************************************/
491 : /* GDALRegister_SNODAS() */
492 : /************************************************************************/
493 :
494 582 : void GDALRegister_SNODAS()
495 :
496 : {
497 : GDALDriver *poDriver;
498 :
499 582 : if( GDALGetDriverByName( "SNODAS" ) == NULL )
500 : {
501 561 : poDriver = new GDALDriver();
502 :
503 561 : poDriver->SetDescription( "SNODAS" );
504 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
505 561 : "Snow Data Assimilation System" );
506 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
507 561 : "frmt_various.html#SNODAS" );
508 561 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "hdr" );
509 :
510 561 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
511 561 : poDriver->pfnOpen = SNODASDataset::Open;
512 561 : poDriver->pfnIdentify = SNODASDataset::Identify;
513 :
514 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
515 : }
516 582 : }
|
