← Back to branch summary

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

« back to all changes in this revision

Viewing changes to frmts/ngsgeoid/ngsgeoiddataset.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/ngsgeoid/ngsgeoiddataset.cpp
 
1
/******************************************************************************
 
2
 * $Id: ngsgeoiddataset.cpp 23334 2011-11-05 22:40:46Z rouault $
 
3
 *
 
4
 * Project:  NGSGEOID driver
 
5
 * Purpose:  GDALDataset driver for NGSGEOID dataset.
 
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 "cpl_vsi_virtual.h"
 
31
#include "cpl_string.h"
 
32
#include "gdal_pam.h"
 
33
#include "ogr_srs_api.h"
 
34
 
 
35
CPL_CVSID("$Id: ngsgeoiddataset.cpp 23334 2011-11-05 22:40:46Z rouault $");
 
36
 
 
37
CPL_C_START
 
38
void    GDALRegister_NGSGEOID(void);
 
39
CPL_C_END
 
40
 
 
41
#define HEADER_SIZE (4 * 8 + 3 * 4)
 
42
 
 
43
/************************************************************************/
 
44
/* ==================================================================== */
 
45
/*                            NGSGEOIDDataset                           */
 
46
/* ==================================================================== */
 
47
/************************************************************************/
 
48
 
 
49
class NGSGEOIDRasterBand;
 
50
 
 
51
class NGSGEOIDDataset : public GDALPamDataset
 
52
{
 
53
    friend class NGSGEOIDRasterBand;
 
54
 
 
55
    VSILFILE   *fp;
 
56
    double      adfGeoTransform[6];
 
57
    int         bIsLittleEndian;
 
58
 
 
59
    static int   GetHeaderInfo( const GByte* pBuffer,
 
60
                                double* padfGeoTransform,
 
61
                                int* pnRows,
 
62
                                int* pnCols,
 
63
                                int* pbIsLittleEndian );
 
64
 
 
65
  public:
 
66
                 NGSGEOIDDataset();
 
67
    virtual     ~NGSGEOIDDataset();
 
68
 
 
69
    virtual CPLErr GetGeoTransform( double * );
 
70
    virtual const char* GetProjectionRef();
 
71
 
 
72
    static GDALDataset *Open( GDALOpenInfo * );
 
73
    static int          Identify( GDALOpenInfo * );
 
74
};
 
75
 
 
76
/************************************************************************/
 
77
/* ==================================================================== */
 
78
/*                          NGSGEOIDRasterBand                          */
 
79
/* ==================================================================== */
 
80
/************************************************************************/
 
81
 
 
82
class NGSGEOIDRasterBand : public GDALPamRasterBand
 
83
{
 
84
    friend class NGSGEOIDDataset;
 
85
 
 
86
  public:
 
87
 
 
88
                NGSGEOIDRasterBand( NGSGEOIDDataset * );
 
89
 
 
90
    virtual CPLErr IReadBlock( int, int, void * );
 
91
 
 
92
    virtual const char* GetUnitType() { return "m"; }
 
93
};
 
94
 
 
95
 
 
96
/************************************************************************/
 
97
/*                        NGSGEOIDRasterBand()                          */
 
98
/************************************************************************/
 
99
 
 
100
NGSGEOIDRasterBand::NGSGEOIDRasterBand( NGSGEOIDDataset *poDS )
 
101
 
 
102
{
 
103
    this->poDS = poDS;
 
104
    this->nBand = 1;
 
105
 
 
106
    eDataType = GDT_Float32;
 
107
 
 
108
    nBlockXSize = poDS->GetRasterXSize();
 
109
    nBlockYSize = 1;
 
110
}
 
111
 
 
112
/************************************************************************/
 
113
/*                             IReadBlock()                             */
 
114
/************************************************************************/
 
115
 
 
116
CPLErr NGSGEOIDRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
 
117
                                       void * pImage )
 
118
 
 
119
{
 
120
    NGSGEOIDDataset *poGDS = (NGSGEOIDDataset *) poDS;
 
121
 
 
122
    /* First values in the file corresponds to the south-most line of the imagery */
 
123
    VSIFSeekL(poGDS->fp,
 
124
              HEADER_SIZE + (nRasterYSize - 1 - nBlockYOff) * nRasterXSize * 4,
 
125
              SEEK_SET);
 
126
 
 
127
    if ((int)VSIFReadL(pImage, 4, nRasterXSize, poGDS->fp) != nRasterXSize)
 
128
        return CE_Failure;
 
129
 
 
130
    if (poGDS->bIsLittleEndian)
 
131
    {
 
132
#ifdef CPL_MSB
 
133
        GDALSwapWords( pImage, 4, nRasterXSize, 4 );
 
134
#endif
 
135
    }
 
136
    else
 
137
    {
 
138
#ifdef CPL_LSB
 
139
        GDALSwapWords( pImage, 4, nRasterXSize, 4 );
 
140
#endif
 
141
    }
 
142
 
 
143
    return CE_None;
 
144
}
 
145
 
 
146
/************************************************************************/
 
147
/*                          ~NGSGEOIDDataset()                          */
 
148
/************************************************************************/
 
149
 
 
150
NGSGEOIDDataset::NGSGEOIDDataset()
 
151
{
 
152
    fp = NULL;
 
153
    adfGeoTransform[0] = 0;
 
154
    adfGeoTransform[1] = 1;
 
155
    adfGeoTransform[2] = 0;
 
156
    adfGeoTransform[3] = 0;
 
157
    adfGeoTransform[4] = 0;
 
158
    adfGeoTransform[5] = 1;
 
159
    bIsLittleEndian = TRUE;
 
160
}
 
161
 
 
162
/************************************************************************/
 
163
/*                           ~NGSGEOIDDataset()                         */
 
164
/************************************************************************/
 
165
 
 
166
NGSGEOIDDataset::~NGSGEOIDDataset()
 
167
 
 
168
{
 
169
    FlushCache();
 
170
    if (fp)
 
171
        VSIFCloseL(fp);
 
172
}
 
173
 
 
174
/************************************************************************/
 
175
/*                            GetHeaderInfo()                           */
 
176
/************************************************************************/
 
177
 
 
178
int NGSGEOIDDataset::GetHeaderInfo( const GByte* pBuffer,
 
179
                                    double* padfGeoTransform,
 
180
                                    int* pnRows,
 
181
                                    int* pnCols,
 
182
                                    int* pbIsLittleEndian )
 
183
{
 
184
    double dfSLAT;
 
185
    double dfWLON;
 
186
    double dfDLAT;
 
187
    double dfDLON;
 
188
    int nNLAT;
 
189
    int nNLON;
 
190
    int nIKIND;
 
191
 
 
192
    /* First check IKIND marker to determine if the file */
 
193
    /* is in little or big-endian order, and if it is a valid */
 
194
    /* NGSGEOID dataset */
 
195
    memcpy(&nIKIND, pBuffer + HEADER_SIZE - 4, 4);
 
196
    CPL_LSBPTR32(&nIKIND);
 
197
    if (nIKIND == 1)
 
198
    {
 
199
        *pbIsLittleEndian = TRUE;
 
200
    }
 
201
    else
 
202
    {
 
203
        memcpy(&nIKIND, pBuffer + HEADER_SIZE - 4, 4);
 
204
        CPL_MSBPTR32(&nIKIND);
 
205
        if (nIKIND == 1)
 
206
        {
 
207
            *pbIsLittleEndian = FALSE;
 
208
        }
 
209
        else
 
210
        {
 
211
            return FALSE;
 
212
        }
 
213
    }
 
214
 
 
215
    memcpy(&dfSLAT, pBuffer, 8);
 
216
    if (*pbIsLittleEndian)
 
217
    {
 
218
        CPL_LSBPTR64(&dfSLAT);
 
219
    }
 
220
    else
 
221
    {
 
222
        CPL_MSBPTR64(&dfSLAT);
 
223
    }
 
224
    pBuffer += 8;
 
225
    memcpy(&dfWLON, pBuffer, 8);
 
226
    if (*pbIsLittleEndian)
 
227
    {
 
228
        CPL_LSBPTR64(&dfWLON);
 
229
    }
 
230
    else
 
231
    {
 
232
        CPL_MSBPTR64(&dfWLON);
 
233
    }
 
234
    pBuffer += 8;
 
235
    memcpy(&dfDLAT, pBuffer, 8);
 
236
    if (*pbIsLittleEndian)
 
237
    {
 
238
        CPL_LSBPTR64(&dfDLAT);
 
239
    }
 
240
    else
 
241
    {
 
242
        CPL_MSBPTR64(&dfDLAT);
 
243
    }
 
244
    pBuffer += 8;
 
245
    memcpy(&dfDLON, pBuffer, 8);
 
246
    if (*pbIsLittleEndian)
 
247
    {
 
248
        CPL_LSBPTR64(&dfDLON);
 
249
    }
 
250
    else
 
251
    {
 
252
        CPL_MSBPTR64(&dfDLON);
 
253
    }
 
254
    pBuffer += 8;
 
255
    memcpy(&nNLAT, pBuffer, 4);
 
256
    if (*pbIsLittleEndian)
 
257
    {
 
258
        CPL_LSBPTR32(&nNLAT);
 
259
    }
 
260
    else
 
261
    {
 
262
        CPL_MSBPTR32(&nNLAT);
 
263
    }
 
264
    pBuffer += 4;
 
265
    memcpy(&nNLON, pBuffer, 4);
 
266
    if (*pbIsLittleEndian)
 
267
    {
 
268
        CPL_LSBPTR32(&nNLON);
 
269
    }
 
270
    else
 
271
    {
 
272
        CPL_MSBPTR32(&nNLON);
 
273
    }
 
274
    pBuffer += 4;
 
275
 
 
276
    /*CPLDebug("NGSGEOID", "SLAT=%f, WLON=%f, DLAT=%f, DLON=%f, NLAT=%d, NLON=%d, IKIND=%d",
 
277
             dfSLAT, dfWLON, dfDLAT, dfDLON, nNLAT, nNLON, nIKIND);*/
 
278
 
 
279
    if (nNLAT <= 0 || nNLON <= 0 || dfDLAT <= 0.0 || dfDLON <= 0.0)
 
280
        return FALSE;
 
281
 
 
282
    /* Grids go over +180 in longitude */
 
283
    if (dfSLAT < -90.0 || dfSLAT + nNLAT * dfDLAT > 90.0 ||
 
284
        dfWLON < -180.0 || dfWLON + nNLON * dfDLON > 360.0)
 
285
        return FALSE;
 
286
 
 
287
    padfGeoTransform[0] = dfWLON - dfDLON / 2;
 
288
    padfGeoTransform[1] = dfDLON;
 
289
    padfGeoTransform[2] = 0.0;
 
290
    padfGeoTransform[3] = dfSLAT + nNLAT * dfDLAT - dfDLAT / 2;
 
291
    padfGeoTransform[4] = 0.0;
 
292
    padfGeoTransform[5] = -dfDLAT;
 
293
 
 
294
    *pnRows = nNLAT;
 
295
    *pnCols = nNLON;
 
296
 
 
297
    return TRUE;
 
298
}
 
299
 
 
300
/************************************************************************/
 
301
/*                             Identify()                               */
 
302
/************************************************************************/
 
303
 
 
304
int NGSGEOIDDataset::Identify( GDALOpenInfo * poOpenInfo )
 
305
{
 
306
    if (poOpenInfo->nHeaderBytes < HEADER_SIZE)
 
307
        return FALSE;
 
308
 
 
309
    double adfGeoTransform[6];
 
310
    int nRows, nCols;
 
311
    int bIsLittleEndian;
 
312
    if ( !GetHeaderInfo( poOpenInfo->pabyHeader,
 
313
                         adfGeoTransform,
 
314
                         &nRows, &nCols, &bIsLittleEndian ) )
 
315
        return FALSE;
 
316
 
 
317
    return TRUE;
 
318
}
 
319
 
 
320
 
 
321
/************************************************************************/
 
322
/*                                Open()                                */
 
323
/************************************************************************/
 
324
 
 
325
GDALDataset *NGSGEOIDDataset::Open( GDALOpenInfo * poOpenInfo )
 
326
 
 
327
{
 
328
    if (!Identify(poOpenInfo))
 
329
        return NULL;
 
330
 
 
331
    if (poOpenInfo->eAccess == GA_Update)
 
332
    {
 
333
        CPLError( CE_Failure, CPLE_NotSupported,
 
334
                  "The NGSGEOID driver does not support update access to existing"
 
335
                  " datasets.\n" );
 
336
        return NULL;
 
337
    }
 
338
 
 
339
    VSILFILE* fp = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
 
340
    if (fp == NULL)
 
341
        return NULL;
 
342
 
 
343
/* -------------------------------------------------------------------- */
 
344
/*      Create a corresponding GDALDataset.                             */
 
345
/* -------------------------------------------------------------------- */
 
346
    NGSGEOIDDataset         *poDS;
 
347
 
 
348
    poDS = new NGSGEOIDDataset();
 
349
    poDS->fp = fp;
 
350
 
 
351
    int nRows, nCols;
 
352
    GetHeaderInfo( poOpenInfo->pabyHeader,
 
353
                   poDS->adfGeoTransform,
 
354
                   &nRows,
 
355
                   &nCols,
 
356
                   &poDS->bIsLittleEndian );
 
357
    poDS->nRasterXSize = nCols;
 
358
    poDS->nRasterYSize = nRows;
 
359
 
 
360
/* -------------------------------------------------------------------- */
 
361
/*      Create band information objects.                                */
 
362
/* -------------------------------------------------------------------- */
 
363
    poDS->nBands = 1;
 
364
    poDS->SetBand( 1, new NGSGEOIDRasterBand( poDS ) );
 
365
 
 
366
/* -------------------------------------------------------------------- */
 
367
/*      Initialize any PAM information.                                 */
 
368
/* -------------------------------------------------------------------- */
 
369
    poDS->SetDescription( poOpenInfo->pszFilename );
 
370
    poDS->TryLoadXML();
 
371
 
 
372
/* -------------------------------------------------------------------- */
 
373
/*      Support overviews.                                              */
 
374
/* -------------------------------------------------------------------- */
 
375
    poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
 
376
    return( poDS );
 
377
}
 
378
 
 
379
/************************************************************************/
 
380
/*                          GetGeoTransform()                           */
 
381
/************************************************************************/
 
382
 
 
383
CPLErr NGSGEOIDDataset::GetGeoTransform( double * padfTransform )
 
384
 
 
385
{
 
386
    memcpy(padfTransform, adfGeoTransform, 6 * sizeof(double));
 
387
 
 
388
    return( CE_None );
 
389
}
 
390
 
 
391
/************************************************************************/
 
392
/*                         GetProjectionRef()                           */
 
393
/************************************************************************/
 
394
 
 
395
const char* NGSGEOIDDataset::GetProjectionRef()
 
396
{
 
397
    return SRS_WKT_WGS84;
 
398
}
 
399
 
 
400
/************************************************************************/
 
401
/*                       GDALRegister_NGSGEOID()                        */
 
402
/************************************************************************/
 
403
 
 
404
void GDALRegister_NGSGEOID()
 
405
 
 
406
{
 
407
    GDALDriver  *poDriver;
 
408
 
 
409
    if( GDALGetDriverByName( "NGSGEOID" ) == NULL )
 
410
    {
 
411
        poDriver = new GDALDriver();
 
412
 
 
413
        poDriver->SetDescription( "NGSGEOID" );
 
414
        poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
 
415
                                   "NOAA NGS Geoid Height Grids" );
 
416
        poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
 
417
                                   "frmt_ngsgeoid.html" );
 
418
        poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "bin" );
 
419
 
 
420
        poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
 
421
 
 
422
        poDriver->pfnOpen = NGSGEOIDDataset::Open;
 
423
        poDriver->pfnIdentify = NGSGEOIDDataset::Identify;
 
424
 
 
425
        GetGDALDriverManager()->RegisterDriver( poDriver );
 
426
    }
 
427
}
 
428