← Back to branch summary

~ubuntu-branches/debian/sid/gdal/sid

« back to all changes in this revision

Viewing changes to frmts/raw/snodasdataset.cpp

  • Committer: Package Import Robot
  • Author(s): Francesco Paolo Lovergine
  • Date: 2012-05-07 15:04:42 UTC
  • mfrom: (5.5.16 experimental)
  • Revision ID: package-import@ubuntu.com-20120507150442-2eks97loeh6rq005
Tags: 1.9.0-1
* Ready for sid, starting transition.
* All symfiles updated to latest builds.
* Added dh_numpy call in debian/rules to depend on numpy ABI.
* Policy bumped to 3.9.3, no changes required.

Show diffs side-by-side

added added

removed removed

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