← Back to branch summary

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

« back to all changes in this revision

Viewing changes to frmts/raw/ntv2dataset.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/ntv2dataset.cpp
 
1
/******************************************************************************
 
2
 * $Id: ntv2dataset.cpp 21682 2011-02-11 21:25:05Z warmerdam $
 
3
 *
 
4
 * Project:  Horizontal Datum Formats
 
5
 * Purpose:  Implementation of NTv2 datum shift format used in Canada, France, 
 
6
 *           Australia and elsewhere.
 
7
 * Author:   Frank Warmerdam, warmerdam@pobox.com
 
8
 * Financial Support: i-cubed (http://www.i-cubed.com)
 
9
 *
 
10
 ******************************************************************************
 
11
 * Copyright (c) 2010, Frank Warmerdam
 
12
 *
 
13
 * Permission is hereby granted, free of charge, to any person obtaining a
 
14
 * copy of this software and associated documentation files (the "Software"),
 
15
 * to deal in the Software without restriction, including without limitation
 
16
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 
17
 * and/or sell copies of the Software, and to permit persons to whom the
 
18
 * Software is furnished to do so, subject to the following conditions:
 
19
 *
 
20
 * The above copyright notice and this permission notice shall be included
 
21
 * in all copies or substantial portions of the Software.
 
22
 *
 
23
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 
24
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 
25
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 
26
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 
27
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 
28
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 
29
 * DEALINGS IN THE SOFTWARE.
 
30
 ****************************************************************************/
 
31
 
 
32
#include "rawdataset.h"
 
33
#include "cpl_string.h"
 
34
#include "ogr_srs_api.h"
 
35
 
 
36
CPL_CVSID("$Id: ntv2dataset.cpp 21682 2011-02-11 21:25:05Z warmerdam $");
 
37
 
 
38
/** 
 
39
 * The header for the file, and each grid consists of 11 16byte records.
 
40
 * The first half is an ASCII label, and the second half is the value
 
41
 * often in a little endian int or float. 
 
42
 *
 
43
 * Example:
 
44
 
 
45
00000000  4e 55 4d 5f 4f 52 45 43  0b 00 00 00 00 00 00 00  |NUM_OREC........|
 
46
00000010  4e 55 4d 5f 53 52 45 43  0b 00 00 00 00 00 00 00  |NUM_SREC........|
 
47
00000020  4e 55 4d 5f 46 49 4c 45  01 00 00 00 00 00 00 00  |NUM_FILE........|
 
48
00000030  47 53 5f 54 59 50 45 20  53 45 43 4f 4e 44 53 20  |GS_TYPE SECONDS |
 
49
00000040  56 45 52 53 49 4f 4e 20  49 47 4e 30 37 5f 30 31  |VERSION IGN07_01|
 
50
00000050  53 59 53 54 45 4d 5f 46  4e 54 46 20 20 20 20 20  |SYSTEM_FNTF     |
 
51
00000060  53 59 53 54 45 4d 5f 54  52 47 46 39 33 20 20 20  |SYSTEM_TRGF93   |
 
52
00000070  4d 41 4a 4f 52 5f 46 20  cd cc cc 4c c2 54 58 41  |MAJOR_F ...L.TXA|
 
53
00000080  4d 49 4e 4f 52 5f 46 20  00 00 00 c0 88 3f 58 41  |MINOR_F .....?XA|
 
54
00000090  4d 41 4a 4f 52 5f 54 20  00 00 00 40 a6 54 58 41  |MAJOR_T ...@.TXA|
 
55
000000a0  4d 49 4e 4f 52 5f 54 20  27 e0 1a 14 c4 3f 58 41  |MINOR_T '....?XA|
 
56
000000b0  53 55 42 5f 4e 41 4d 45  46 52 41 4e 43 45 20 20  |SUB_NAMEFRANCE  |
 
57
000000c0  50 41 52 45 4e 54 20 20  4e 4f 4e 45 20 20 20 20  |PARENT  NONE    |
 
58
000000d0  43 52 45 41 54 45 44 20  33 31 2f 31 30 2f 30 37  |CREATED 31/10/07|
 
59
000000e0  55 50 44 41 54 45 44 20  20 20 20 20 20 20 20 20  |UPDATED         |
 
60
000000f0  53 5f 4c 41 54 20 20 20  00 00 00 00 80 04 02 41  |S_LAT   .......A|
 
61
00000100  4e 5f 4c 41 54 20 20 20  00 00 00 00 00 da 06 41  |N_LAT   .......A|
 
62
00000110  45 5f 4c 4f 4e 47 20 20  00 00 00 00 00 94 e1 c0  |E_LONG  ........|
 
63
00000120  57 5f 4c 4f 4e 47 20 20  00 00 00 00 00 56 d3 40  |W_LONG  .....V.@|
 
64
00000130  4c 41 54 5f 49 4e 43 20  00 00 00 00 00 80 76 40  |LAT_INC ......v@|
 
65
00000140  4c 4f 4e 47 5f 49 4e 43  00 00 00 00 00 80 76 40  |LONG_INC......v@|
 
66
00000150  47 53 5f 43 4f 55 4e 54  a4 43 00 00 00 00 00 00  |GS_COUNT.C......|
 
67
00000160  94 f7 c1 3e 70 ee a3 3f  2a c7 84 3d ff 42 af 3d  |...>p..?*..=.B.=|
 
68
 
 
69
the actual grid data is a raster with 4 float32 bands (lat offset, long
 
70
offset, lat error, long error).  The offset values are in arc seconds.
 
71
The grid is flipped in the x and y axis from our usual GDAL orientation.
 
72
That is, the first pixel is the south east corner with scanlines going
 
73
east to west, and rows from south to north.  As a GDAL dataset we represent
 
74
these both in the more conventional orientation.
 
75
 */
 
76
 
 
77
/************************************************************************/
 
78
/* ==================================================================== */
 
79
/*                              NTv2Dataset                             */
 
80
/* ==================================================================== */
 
81
/************************************************************************/
 
82
 
 
83
class NTv2Dataset : public RawDataset
 
84
{
 
85
  public:
 
86
    VSILFILE    *fpImage;       // image data file.
 
87
 
 
88
    int         nRecordLength;
 
89
    vsi_l_offset nGridOffset;
 
90
    
 
91
    double      adfGeoTransform[6];
 
92
 
 
93
    void        CaptureMetadataItem( char *pszItem );
 
94
 
 
95
    int         OpenGrid( char *pachGridHeader, vsi_l_offset nDataStart );
 
96
 
 
97
  public:
 
98
                NTv2Dataset();
 
99
                ~NTv2Dataset();
 
100
    
 
101
    virtual CPLErr SetGeoTransform( double * padfTransform );
 
102
    virtual CPLErr GetGeoTransform( double * padfTransform );
 
103
    virtual const char *GetProjectionRef();
 
104
    virtual void   FlushCache(void);
 
105
 
 
106
    static GDALDataset *Open( GDALOpenInfo * );
 
107
    static int          Identify( GDALOpenInfo * );
 
108
    static GDALDataset *Create( const char * pszFilename,
 
109
                                int nXSize, int nYSize, int nBands,
 
110
                                GDALDataType eType, char ** papszOptions );
 
111
};
 
112
 
 
113
/************************************************************************/
 
114
/* ==================================================================== */
 
115
/*                              NTv2Dataset                             */
 
116
/* ==================================================================== */
 
117
/************************************************************************/
 
118
 
 
119
/************************************************************************/
 
120
/*                             NTv2Dataset()                          */
 
121
/************************************************************************/
 
122
 
 
123
NTv2Dataset::NTv2Dataset()
 
124
{
 
125
    fpImage = NULL;
 
126
}
 
127
 
 
128
/************************************************************************/
 
129
/*                            ~NTv2Dataset()                          */
 
130
/************************************************************************/
 
131
 
 
132
NTv2Dataset::~NTv2Dataset()
 
133
 
 
134
{
 
135
    FlushCache();
 
136
 
 
137
    if( fpImage != NULL )
 
138
        VSIFCloseL( fpImage );
 
139
}
 
140
 
 
141
/************************************************************************/
 
142
/*                             FlushCache()                             */
 
143
/************************************************************************/
 
144
 
 
145
void NTv2Dataset::FlushCache()
 
146
 
 
147
{
 
148
/* -------------------------------------------------------------------- */
 
149
/*      Nothing to do in readonly mode, or if nothing seems to have     */
 
150
/*      changed metadata wise.                                          */
 
151
/* -------------------------------------------------------------------- */
 
152
    if( eAccess != GA_Update || !(GetPamFlags() & GPF_DIRTY) )
 
153
    {
 
154
        RawDataset::FlushCache();
 
155
        return;
 
156
    }
 
157
 
 
158
/* -------------------------------------------------------------------- */
 
159
/*      Load grid and file headers.                                     */
 
160
/* -------------------------------------------------------------------- */
 
161
    char achFileHeader[11*16];
 
162
    char achGridHeader[11*16];
 
163
 
 
164
    VSIFSeekL( fpImage, 0, SEEK_SET );
 
165
    VSIFReadL( achFileHeader, 11, 16, fpImage );
 
166
 
 
167
    VSIFSeekL( fpImage, nGridOffset, SEEK_SET );
 
168
    VSIFReadL( achGridHeader, 11, 16, fpImage );
 
169
 
 
170
/* -------------------------------------------------------------------- */
 
171
/*      Update the grid, and file headers with any available            */
 
172
/*      metadata.  If all available metadata is recognised then mark    */
 
173
/*      things "clean" from a PAM point of view.                        */
 
174
/* -------------------------------------------------------------------- */
 
175
    char **papszMD = GetMetadata();
 
176
    int i;
 
177
    int bSomeLeftOver = FALSE;
 
178
 
 
179
    for( i = 0; papszMD != NULL && papszMD[i] != NULL; i++ )
 
180
    {
 
181
        char *pszKey = NULL;
 
182
        const char *pszValue = CPLParseNameValue( papszMD[i], &pszKey );
 
183
 
 
184
        if( EQUAL(pszKey,"GS_TYPE") )
 
185
        {
 
186
            memcpy( achFileHeader + 3*16+8, "        ", 8 );
 
187
            memcpy( achFileHeader + 3*16+8, pszValue, MIN(8,strlen(pszValue)) );
 
188
        }
 
189
        else if( EQUAL(pszKey,"VERSION") )
 
190
        {
 
191
            memcpy( achFileHeader + 4*16+8, "        ", 8 );
 
192
            memcpy( achFileHeader + 4*16+8, pszValue, MIN(8,strlen(pszValue)) );
 
193
        }
 
194
        else if( EQUAL(pszKey,"SYSTEM_F") )
 
195
        {
 
196
            memcpy( achFileHeader + 5*16+8, "        ", 8 );
 
197
            memcpy( achFileHeader + 5*16+8, pszValue, MIN(8,strlen(pszValue)) );
 
198
        }
 
199
        else if( EQUAL(pszKey,"SYSTEM_T") )
 
200
        {
 
201
            memcpy( achFileHeader + 6*16+8, "        ", 8 );
 
202
            memcpy( achFileHeader + 6*16+8, pszValue, MIN(8,strlen(pszValue)) );
 
203
        }
 
204
        else if( EQUAL(pszKey,"MAJOR_F") )
 
205
        {
 
206
            double dfValue = atof(pszValue);
 
207
            CPL_LSBPTR64( &dfValue );
 
208
            memcpy( achFileHeader + 7*16+8, &dfValue, 8 );
 
209
        }
 
210
        else if( EQUAL(pszKey,"MINOR_F") )
 
211
        {
 
212
            double dfValue = atof(pszValue);
 
213
            CPL_LSBPTR64( &dfValue );
 
214
            memcpy( achFileHeader + 8*16+8, &dfValue, 8 );
 
215
        }
 
216
        else if( EQUAL(pszKey,"MAJOR_T") )
 
217
        {
 
218
            double dfValue = atof(pszValue);
 
219
            CPL_LSBPTR64( &dfValue );
 
220
            memcpy( achFileHeader + 9*16+8, &dfValue, 8 );
 
221
        }
 
222
        else if( EQUAL(pszKey,"MINOR_T") )
 
223
        {
 
224
            double dfValue = atof(pszValue);
 
225
            CPL_LSBPTR64( &dfValue );
 
226
            memcpy( achFileHeader + 10*16+8, &dfValue, 8 );
 
227
        }
 
228
        else if( EQUAL(pszKey,"SUB_NAME") )
 
229
        {
 
230
            memcpy( achGridHeader + 0*16+8, "        ", 8 );
 
231
            memcpy( achGridHeader + 0*16+8, pszValue, MIN(8,strlen(pszValue)) );
 
232
        }
 
233
        else if( EQUAL(pszKey,"PARENT") )
 
234
        {
 
235
            memcpy( achGridHeader + 1*16+8, "        ", 8 );
 
236
            memcpy( achGridHeader + 1*16+8, pszValue, MIN(8,strlen(pszValue)) );
 
237
        }
 
238
        else if( EQUAL(pszKey,"CREATED") )
 
239
        {
 
240
            memcpy( achGridHeader + 2*16+8, "        ", 8 );
 
241
            memcpy( achGridHeader + 2*16+8, pszValue, MIN(8,strlen(pszValue)) );
 
242
        }
 
243
        else if( EQUAL(pszKey,"UPDATED") )
 
244
        {
 
245
            memcpy( achGridHeader + 3*16+8, "        ", 8 );
 
246
            memcpy( achGridHeader + 3*16+8, pszValue, MIN(8,strlen(pszValue)) );
 
247
        }
 
248
        else
 
249
        {
 
250
            bSomeLeftOver = TRUE;
 
251
        }
 
252
        
 
253
        CPLFree( pszKey );
 
254
    }
 
255
 
 
256
/* -------------------------------------------------------------------- */
 
257
/*      Load grid and file headers.                                     */
 
258
/* -------------------------------------------------------------------- */
 
259
    VSIFSeekL( fpImage, 0, SEEK_SET );
 
260
    VSIFWriteL( achFileHeader, 11, 16, fpImage );
 
261
 
 
262
    VSIFSeekL( fpImage, nGridOffset, SEEK_SET );
 
263
    VSIFWriteL( achGridHeader, 11, 16, fpImage );
 
264
 
 
265
/* -------------------------------------------------------------------- */
 
266
/*      Clear flags if we got everything, then let pam and below do     */
 
267
/*      their flushing.                                                 */
 
268
/* -------------------------------------------------------------------- */
 
269
    if( !bSomeLeftOver )
 
270
        SetPamFlags( GetPamFlags() & (~GPF_DIRTY) );
 
271
 
 
272
    RawDataset::FlushCache();
 
273
}
 
274
 
 
275
/************************************************************************/
 
276
/*                              Identify()                              */
 
277
/************************************************************************/
 
278
 
 
279
int NTv2Dataset::Identify( GDALOpenInfo *poOpenInfo )
 
280
 
 
281
{
 
282
    if( EQUALN(poOpenInfo->pszFilename,"NTv2:",5) )
 
283
        return TRUE;
 
284
    
 
285
    if( poOpenInfo->nHeaderBytes < 64 )
 
286
        return FALSE;
 
287
 
 
288
    if( !EQUALN((const char *)poOpenInfo->pabyHeader + 0, "NUM_OREC", 8 ) )
 
289
        return FALSE;
 
290
 
 
291
    if( !EQUALN((const char *)poOpenInfo->pabyHeader +16, "NUM_SREC", 8 ) )
 
292
        return FALSE;
 
293
 
 
294
    return TRUE;
 
295
}
 
296
 
 
297
/************************************************************************/
 
298
/*                                Open()                                */
 
299
/************************************************************************/
 
300
 
 
301
GDALDataset *NTv2Dataset::Open( GDALOpenInfo * poOpenInfo )
 
302
 
 
303
{
 
304
    if( !Identify( poOpenInfo ) )
 
305
        return NULL;
 
306
        
 
307
/* -------------------------------------------------------------------- */
 
308
/*      Are we targetting a particular grid?                            */
 
309
/* -------------------------------------------------------------------- */
 
310
    CPLString osFilename;
 
311
    int iTargetGrid = -1;
 
312
 
 
313
    if( EQUALN(poOpenInfo->pszFilename,"NTv2:",5) )
 
314
    {
 
315
        const char *pszRest = poOpenInfo->pszFilename+5;
 
316
        
 
317
        iTargetGrid = atoi(pszRest);
 
318
        while( *pszRest != '\0' && *pszRest != ':' )
 
319
            pszRest++;
 
320
     
 
321
        if( *pszRest == ':' )
 
322
            pszRest++;
 
323
        
 
324
        osFilename = pszRest;
 
325
    }
 
326
    else
 
327
        osFilename = poOpenInfo->pszFilename;
 
328
    
 
329
/* -------------------------------------------------------------------- */
 
330
/*      Create a corresponding GDALDataset.                             */
 
331
/* -------------------------------------------------------------------- */
 
332
    NTv2Dataset         *poDS;
 
333
 
 
334
    poDS = new NTv2Dataset();
 
335
    poDS->eAccess = poOpenInfo->eAccess;
 
336
 
 
337
/* -------------------------------------------------------------------- */
 
338
/*      Open the file.                                                  */
 
339
/* -------------------------------------------------------------------- */
 
340
    if( poOpenInfo->eAccess == GA_ReadOnly )
 
341
        poDS->fpImage = VSIFOpenL( osFilename, "rb" );
 
342
    else
 
343
        poDS->fpImage = VSIFOpenL( osFilename, "rb+" );
 
344
 
 
345
    if( poDS->fpImage == NULL )
 
346
    {
 
347
        delete poDS;
 
348
        return NULL;
 
349
    }
 
350
 
 
351
/* -------------------------------------------------------------------- */
 
352
/*      Read the file header.                                           */
 
353
/* -------------------------------------------------------------------- */
 
354
    char  achHeader[11*16];
 
355
    GInt32 nSubFileCount;
 
356
    double dfValue;
 
357
    CPLString osFValue;
 
358
 
 
359
    VSIFSeekL( poDS->fpImage, 0, SEEK_SET );
 
360
    VSIFReadL( achHeader, 11, 16, poDS->fpImage );
 
361
 
 
362
    CPL_LSBPTR32( achHeader + 2*16 + 8 );
 
363
    memcpy( &nSubFileCount, achHeader + 2*16 + 8, 4 );
 
364
    if (nSubFileCount <= 0 || nSubFileCount >= 1024)
 
365
    {
 
366
        CPLError(CE_Failure, CPLE_AppDefined,
 
367
                  "Invalid value for NUM_FILE : %d", nSubFileCount);
 
368
        delete poDS;
 
369
        return NULL;
 
370
    }
 
371
 
 
372
    poDS->CaptureMetadataItem( achHeader + 3*16 );
 
373
    poDS->CaptureMetadataItem( achHeader + 4*16 );
 
374
    poDS->CaptureMetadataItem( achHeader + 5*16 );
 
375
    poDS->CaptureMetadataItem( achHeader + 6*16 );
 
376
 
 
377
    memcpy( &dfValue, achHeader + 7*16 + 8, 8 );
 
378
    CPL_LSBPTR64( &dfValue );
 
379
    osFValue.Printf( "%.15g", dfValue );
 
380
    poDS->SetMetadataItem( "MAJOR_F", osFValue );
 
381
    
 
382
    memcpy( &dfValue, achHeader + 8*16 + 8, 8 );
 
383
    CPL_LSBPTR64( &dfValue );
 
384
    osFValue.Printf( "%.15g", dfValue );
 
385
    poDS->SetMetadataItem( "MINOR_F", osFValue );
 
386
    
 
387
    memcpy( &dfValue, achHeader + 9*16 + 8, 8 );
 
388
    CPL_LSBPTR64( &dfValue );
 
389
    osFValue.Printf( "%.15g", dfValue );
 
390
    poDS->SetMetadataItem( "MAJOR_T", osFValue );
 
391
    
 
392
    memcpy( &dfValue, achHeader + 10*16 + 8, 8 );
 
393
    CPL_LSBPTR64( &dfValue );
 
394
    osFValue.Printf( "%.15g", dfValue );
 
395
    poDS->SetMetadataItem( "MINOR_T", osFValue );
 
396
 
 
397
/* ==================================================================== */
 
398
/*      Loop over grids.                                                */
 
399
/* ==================================================================== */
 
400
    int iGrid;
 
401
    vsi_l_offset nGridOffset = sizeof(achHeader);
 
402
 
 
403
    for( iGrid = 0; iGrid < nSubFileCount; iGrid++ )
 
404
    {
 
405
        CPLString  osSubName;
 
406
        int i;
 
407
        GUInt32 nGSCount;
 
408
 
 
409
        VSIFSeekL( poDS->fpImage, nGridOffset, SEEK_SET );
 
410
        if (VSIFReadL( achHeader, 11, 16, poDS->fpImage ) != 16)
 
411
        {
 
412
            CPLError(CE_Failure, CPLE_AppDefined,
 
413
                     "Cannot read header for subfile %d", iGrid);
 
414
            delete poDS;
 
415
            return NULL;
 
416
        }
 
417
 
 
418
        for( i = 4; i <= 9; i++ )
 
419
            CPL_LSBPTR64( achHeader + i*16 + 8 );
 
420
        
 
421
        CPL_LSBPTR32( achHeader + 10*16 + 8 );
 
422
        
 
423
        memcpy( &nGSCount, achHeader + 10*16 + 8, 4 );
 
424
 
 
425
        osSubName.assign( achHeader + 8, 8 );
 
426
        osSubName.Trim();
 
427
 
 
428
        // If this is our target grid, open it as a dataset.
 
429
        if( iTargetGrid == iGrid || (iTargetGrid == -1 && iGrid == 0) )
 
430
        {
 
431
            if( !poDS->OpenGrid( achHeader, nGridOffset ) )
 
432
            {
 
433
                delete poDS;
 
434
                return NULL;
 
435
            }
 
436
        }
 
437
 
 
438
        // If we are opening the file as a whole, list subdatasets.
 
439
        if( iTargetGrid == -1 )
 
440
        {
 
441
            CPLString osKey, osValue;
 
442
 
 
443
            osKey.Printf( "SUBDATASET_%d_NAME", iGrid );
 
444
            osValue.Printf( "NTv2:%d:%s", iGrid, osFilename.c_str() );
 
445
            poDS->SetMetadataItem( osKey, osValue, "SUBDATASETS" );
 
446
 
 
447
            osKey.Printf( "SUBDATASET_%d_DESC", iGrid );
 
448
            osValue.Printf( "%s", osSubName.c_str() );
 
449
            poDS->SetMetadataItem( osKey, osValue, "SUBDATASETS" );
 
450
        }
 
451
 
 
452
        nGridOffset += (11+(vsi_l_offset)nGSCount) * 16;
 
453
    }
 
454
 
 
455
/* -------------------------------------------------------------------- */
 
456
/*      Initialize any PAM information.                                 */
 
457
/* -------------------------------------------------------------------- */
 
458
    poDS->SetDescription( poOpenInfo->pszFilename );
 
459
    poDS->TryLoadXML();
 
460
 
 
461
/* -------------------------------------------------------------------- */
 
462
/*      Check for overviews.                                            */
 
463
/* -------------------------------------------------------------------- */
 
464
    poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
 
465
 
 
466
    return( poDS );
 
467
}
 
468
 
 
469
/************************************************************************/
 
470
/*                              OpenGrid()                              */
 
471
/*                                                                      */
 
472
/*      Note that the caller will already have byte swapped needed      */
 
473
/*      portions of the header.                                         */
 
474
/************************************************************************/
 
475
 
 
476
int NTv2Dataset::OpenGrid( char *pachHeader, vsi_l_offset nGridOffset )
 
477
 
 
478
{
 
479
    this->nGridOffset = nGridOffset;
 
480
 
 
481
/* -------------------------------------------------------------------- */
 
482
/*      Read the grid header.                                           */
 
483
/* -------------------------------------------------------------------- */
 
484
    double s_lat, n_lat, e_long, w_long, lat_inc, long_inc;
 
485
 
 
486
    CaptureMetadataItem( pachHeader + 0*16 );
 
487
    CaptureMetadataItem( pachHeader + 1*16 );
 
488
    CaptureMetadataItem( pachHeader + 2*16 );
 
489
    CaptureMetadataItem( pachHeader + 3*16 );
 
490
 
 
491
    memcpy( &s_lat,  pachHeader + 4*16 + 8, 8 );
 
492
    memcpy( &n_lat,  pachHeader + 5*16 + 8, 8 );
 
493
    memcpy( &e_long, pachHeader + 6*16 + 8, 8 );
 
494
    memcpy( &w_long, pachHeader + 7*16 + 8, 8 );
 
495
    memcpy( &lat_inc, pachHeader + 8*16 + 8, 8 );
 
496
    memcpy( &long_inc, pachHeader + 9*16 + 8, 8 );
 
497
 
 
498
    e_long *= -1;
 
499
    w_long *= -1;
 
500
 
 
501
    nRasterXSize = (int) floor((e_long - w_long) / long_inc + 1.5);
 
502
    nRasterYSize = (int) floor((n_lat - s_lat) / lat_inc + 1.5);
 
503
 
 
504
    if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize))
 
505
        return FALSE;
 
506
 
 
507
/* -------------------------------------------------------------------- */
 
508
/*      Create band information object.                                 */
 
509
/*                                                                      */
 
510
/*      We use unusual offsets to remap from bottom to top, to top      */
 
511
/*      to bottom orientation, and also to remap east to west, to       */
 
512
/*      west to east.                                                   */
 
513
/* -------------------------------------------------------------------- */
 
514
    int iBand;
 
515
    
 
516
    for( iBand = 0; iBand < 4; iBand++ )
 
517
    {
 
518
        RawRasterBand *poBand = 
 
519
            new RawRasterBand( this, iBand+1, fpImage, 
 
520
                               nGridOffset + 4*iBand + 11*16
 
521
                               + (nRasterXSize-1) * 16
 
522
                               + (nRasterYSize-1) * 16 * nRasterXSize,
 
523
                               -16, -16 * nRasterXSize,
 
524
                               GDT_Float32, CPL_IS_LSB, TRUE, FALSE );
 
525
        SetBand( iBand+1, poBand );
 
526
    }
 
527
    
 
528
    GetRasterBand(1)->SetDescription( "Latitude Offset" );
 
529
    GetRasterBand(2)->SetDescription( "Longitude Offset" );
 
530
    GetRasterBand(3)->SetDescription( "Latitude Error" );
 
531
    GetRasterBand(4)->SetDescription( "Longitude Error" );
 
532
    
 
533
/* -------------------------------------------------------------------- */
 
534
/*      Setup georeferencing.                                           */
 
535
/* -------------------------------------------------------------------- */
 
536
    adfGeoTransform[0] = (w_long - long_inc*0.5) / 3600.0;
 
537
    adfGeoTransform[1] = long_inc / 3600.0;
 
538
    adfGeoTransform[2] = 0.0;
 
539
    adfGeoTransform[3] = (n_lat + lat_inc*0.5) / 3600.0;
 
540
    adfGeoTransform[4] = 0.0;
 
541
    adfGeoTransform[5] = (-1 * lat_inc) / 3600.0;
 
542
 
 
543
    return TRUE;
 
544
}
 
545
 
 
546
/************************************************************************/
 
547
/*                        CaptureMetadataItem()                         */
 
548
/************************************************************************/
 
549
 
 
550
void NTv2Dataset::CaptureMetadataItem( char *pszItem )
 
551
 
 
552
{
 
553
    CPLString osKey, osValue;
 
554
 
 
555
    osKey.assign( pszItem, 8 );
 
556
    osValue.assign( pszItem+8, 8 );
 
557
 
 
558
    SetMetadataItem( osKey.Trim(), osValue.Trim() );
 
559
}
 
560
 
 
561
/************************************************************************/
 
562
/*                          GetGeoTransform()                           */
 
563
/************************************************************************/
 
564
 
 
565
CPLErr NTv2Dataset::GetGeoTransform( double * padfTransform )
 
566
 
 
567
{
 
568
    memcpy( padfTransform, adfGeoTransform, sizeof(double)*6 );
 
569
    return CE_None;
 
570
}
 
571
 
 
572
/************************************************************************/
 
573
/*                          SetGeoTransform()                           */
 
574
/************************************************************************/
 
575
 
 
576
CPLErr NTv2Dataset::SetGeoTransform( double * padfTransform )
 
577
 
 
578
{
 
579
    if( eAccess == GA_ReadOnly )
 
580
    {
 
581
        CPLError( CE_Failure, CPLE_NoWriteAccess,
 
582
                  "Unable to update geotransform on readonly file." ); 
 
583
        return CE_Failure;
 
584
    }
 
585
 
 
586
    if( padfTransform[2] != 0.0 || padfTransform[4] != 0.0 )
 
587
    {
 
588
        CPLError( CE_Failure, CPLE_AppDefined,
 
589
                  "Rotated and sheared geotransforms not supported for NTv2."); 
 
590
        return CE_Failure;
 
591
    }
 
592
 
 
593
    memcpy( adfGeoTransform, padfTransform, sizeof(double)*6 );
 
594
 
 
595
/* -------------------------------------------------------------------- */
 
596
/*      Update grid header.                                             */
 
597
/* -------------------------------------------------------------------- */
 
598
    double dfValue;
 
599
    char   achHeader[11*16];
 
600
 
 
601
    // read grid header
 
602
    VSIFSeekL( fpImage, nGridOffset, SEEK_SET );
 
603
    VSIFReadL( achHeader, 11, 16, fpImage );
 
604
 
 
605
    // S_LAT
 
606
    dfValue = 3600 * (adfGeoTransform[3] + (nRasterYSize-0.5) * adfGeoTransform[5]);
 
607
    CPL_LSBPTR64( &dfValue );
 
608
    memcpy( achHeader +  4*16 + 8, &dfValue, 8 );
 
609
 
 
610
    // N_LAT
 
611
    dfValue = 3600 * (adfGeoTransform[3] + 0.5 * adfGeoTransform[5]);
 
612
    CPL_LSBPTR64( &dfValue );
 
613
    memcpy( achHeader +  5*16 + 8, &dfValue, 8 );
 
614
 
 
615
    // E_LONG
 
616
    dfValue = -3600 * (adfGeoTransform[0] + (nRasterXSize-0.5)*adfGeoTransform[1]);
 
617
    CPL_LSBPTR64( &dfValue );
 
618
    memcpy( achHeader +  6*16 + 8, &dfValue, 8 );
 
619
 
 
620
    // W_LONG
 
621
    dfValue = -3600 * (adfGeoTransform[0] + 0.5 * adfGeoTransform[1]);
 
622
    CPL_LSBPTR64( &dfValue );
 
623
    memcpy( achHeader +  7*16 + 8, &dfValue, 8 );
 
624
 
 
625
    // LAT_INC
 
626
    dfValue = -3600 * adfGeoTransform[5];
 
627
    CPL_LSBPTR64( &dfValue );
 
628
    memcpy( achHeader +  8*16 + 8, &dfValue, 8 );
 
629
    
 
630
    // LONG_INC
 
631
    dfValue = 3600 * adfGeoTransform[1];
 
632
    CPL_LSBPTR64( &dfValue );
 
633
    memcpy( achHeader +  9*16 + 8, &dfValue, 8 );
 
634
    
 
635
    // write grid header.
 
636
    VSIFSeekL( fpImage, nGridOffset, SEEK_SET );
 
637
    VSIFWriteL( achHeader, 11, 16, fpImage );
 
638
 
 
639
    return CE_None;
 
640
}
 
641
 
 
642
 
 
643
/************************************************************************/
 
644
/*                          GetProjectionRef()                          */
 
645
/************************************************************************/
 
646
 
 
647
const char *NTv2Dataset::GetProjectionRef()
 
648
 
 
649
{
 
650
    return SRS_WKT_WGS84;
 
651
}
 
652
 
 
653
/************************************************************************/
 
654
/*                               Create()                               */
 
655
/************************************************************************/
 
656
 
 
657
GDALDataset *NTv2Dataset::Create( const char * pszFilename,
 
658
                                  int nXSize, int nYSize, int nBands,
 
659
                                  GDALDataType eType,
 
660
                                  char ** papszOptions )
 
661
 
 
662
{
 
663
    if( eType != GDT_Float32 )
 
664
    {
 
665
        CPLError(CE_Failure, CPLE_AppDefined,
 
666
                 "Attempt to create NTv2 file with unsupported data type '%s'.",
 
667
                 GDALGetDataTypeName( eType ) );
 
668
        return NULL;
 
669
    }
 
670
 
 
671
/* -------------------------------------------------------------------- */
 
672
/*      Are we extending an existing file?                              */
 
673
/* -------------------------------------------------------------------- */
 
674
    VSILFILE    *fp;
 
675
    GUInt32   nNumFile = 1;
 
676
 
 
677
    int bAppend = CSLFetchBoolean(papszOptions,"APPEND_SUBDATASET",FALSE);
 
678
    
 
679
/* -------------------------------------------------------------------- */
 
680
/*      Try to open or create file.                                     */
 
681
/* -------------------------------------------------------------------- */
 
682
    if( bAppend )
 
683
        fp = VSIFOpenL( pszFilename, "rb+" );
 
684
    else
 
685
        fp = VSIFOpenL( pszFilename, "wb" );
 
686
    
 
687
    if( fp == NULL )
 
688
    {
 
689
        CPLError( CE_Failure, CPLE_OpenFailed,
 
690
                  "Attempt to open/create file `%s' failed.\n",
 
691
                  pszFilename );
 
692
        return NULL;
 
693
    }
 
694
 
 
695
/* -------------------------------------------------------------------- */
 
696
/*      Create a file level header if we are creating new.              */
 
697
/* -------------------------------------------------------------------- */
 
698
    char achHeader[11*16];
 
699
    const char *pszValue;
 
700
 
 
701
    if( !bAppend )
 
702
    {
 
703
        memset( achHeader, 0, sizeof(achHeader) );
 
704
        
 
705
        memcpy( achHeader +  0*16, "NUM_OREC", 8 );
 
706
        achHeader[ 0*16 + 8] = 0xb;
 
707
 
 
708
        memcpy( achHeader +  1*16, "NUM_SREC", 8 );
 
709
        achHeader[ 1*16 + 8] = 0xb;
 
710
 
 
711
        memcpy( achHeader +  2*16, "NUM_FILE", 8 );
 
712
        achHeader[ 2*16 + 8] = 0x1;
 
713
 
 
714
        memcpy( achHeader +  3*16, "GS_TYPE         ", 16 );
 
715
        pszValue = CSLFetchNameValueDef( papszOptions, "GS_TYPE", "SECONDS");
 
716
        memcpy( achHeader +  3*16+8, pszValue, MIN(16,strlen(pszValue)) );
 
717
 
 
718
        memcpy( achHeader +  4*16, "VERSION         ", 16 );
 
719
        pszValue = CSLFetchNameValueDef( papszOptions, "VERSION", "" );
 
720
        memcpy( achHeader +  4*16+8, pszValue, MIN(16,strlen(pszValue)) );
 
721
 
 
722
        memcpy( achHeader +  5*16, "SYSTEM_F        ", 16 );
 
723
        pszValue = CSLFetchNameValueDef( papszOptions, "SYSTEM_F", "" );
 
724
        memcpy( achHeader +  5*16+8, pszValue, MIN(16,strlen(pszValue)) );
 
725
 
 
726
        memcpy( achHeader +  6*16, "SYSTEM_T        ", 16 );
 
727
        pszValue = CSLFetchNameValueDef( papszOptions, "SYSTEM_T", "" );
 
728
        memcpy( achHeader +  6*16+8, pszValue, MIN(16,strlen(pszValue)) );
 
729
 
 
730
        memcpy( achHeader +  7*16, "MAJOR_F ", 8);
 
731
        memcpy( achHeader +  8*16, "MINOR_F ", 8 );
 
732
        memcpy( achHeader +  9*16, "MAJOR_T ", 8 );
 
733
        memcpy( achHeader + 10*16, "MINOR_T ", 8 );
 
734
 
 
735
        VSIFWriteL( achHeader, 1, sizeof(achHeader), fp );
 
736
    }
 
737
 
 
738
/* -------------------------------------------------------------------- */
 
739
/*      Otherwise update the header with an increased subfile count,    */
 
740
/*      and advanced to the last record of the file.                    */
 
741
/* -------------------------------------------------------------------- */
 
742
    else
 
743
    {
 
744
        VSIFSeekL( fp, 2*16 + 8, SEEK_SET );
 
745
        VSIFReadL( &nNumFile, 1, 4, fp );
 
746
        CPL_LSBPTR32( &nNumFile );
 
747
 
 
748
        nNumFile++;
 
749
        
 
750
        CPL_LSBPTR32( &nNumFile );
 
751
        VSIFSeekL( fp, 2*16 + 8, SEEK_SET );
 
752
        VSIFWriteL( &nNumFile, 1, 4, fp );
 
753
 
 
754
        vsi_l_offset nEnd;
 
755
 
 
756
        VSIFSeekL( fp, 0, SEEK_END );
 
757
        nEnd = VSIFTellL( fp );
 
758
        VSIFSeekL( fp, nEnd-16, SEEK_SET );
 
759
    }
 
760
 
 
761
/* -------------------------------------------------------------------- */
 
762
/*      Write the grid header.                                          */
 
763
/* -------------------------------------------------------------------- */
 
764
    memset( achHeader, 0, sizeof(achHeader) );
 
765
 
 
766
    memcpy( achHeader +  0*16, "SUB_NAME        ", 16 );
 
767
    pszValue = CSLFetchNameValueDef( papszOptions, "SUB_NAME", "" );
 
768
    memcpy( achHeader +  0*16+8, pszValue, MIN(16,strlen(pszValue)) );
 
769
    
 
770
    memcpy( achHeader +  1*16, "PARENT          ", 16 );
 
771
    pszValue = CSLFetchNameValueDef( papszOptions, "PARENT", "NONE" );
 
772
    memcpy( achHeader +  1*16+8, pszValue, MIN(16,strlen(pszValue)) );
 
773
    
 
774
    memcpy( achHeader +  2*16, "CREATED         ", 16 );
 
775
    pszValue = CSLFetchNameValueDef( papszOptions, "CREATED", "" );
 
776
    memcpy( achHeader +  2*16+8, pszValue, MIN(16,strlen(pszValue)) );
 
777
    
 
778
    memcpy( achHeader +  3*16, "UPDATED         ", 16 );
 
779
    pszValue = CSLFetchNameValueDef( papszOptions, "UPDATED", "" );
 
780
    memcpy( achHeader +  3*16+8, pszValue, MIN(16,strlen(pszValue)) );
 
781
 
 
782
    double dfValue;
 
783
 
 
784
    memcpy( achHeader +  4*16, "S_LAT   ", 8 );
 
785
    dfValue = 0;
 
786
    CPL_LSBPTR64( &dfValue );
 
787
    memcpy( achHeader +  4*16 + 8, &dfValue, 8 );
 
788
 
 
789
    memcpy( achHeader +  5*16, "N_LAT   ", 8 );
 
790
    dfValue = nYSize-1;
 
791
    CPL_LSBPTR64( &dfValue );
 
792
    memcpy( achHeader +  5*16 + 8, &dfValue, 8 );
 
793
 
 
794
    memcpy( achHeader +  6*16, "E_LONG  ", 8 );
 
795
    dfValue = -1*(nXSize-1);
 
796
    CPL_LSBPTR64( &dfValue );
 
797
    memcpy( achHeader +  6*16 + 8, &dfValue, 8 );
 
798
 
 
799
    memcpy( achHeader +  7*16, "W_LONG  ", 8 );
 
800
    dfValue = 0;
 
801
    CPL_LSBPTR64( &dfValue );
 
802
    memcpy( achHeader +  7*16 + 8, &dfValue, 8 );
 
803
 
 
804
    memcpy( achHeader +  8*16, "LAT_INC ", 8 );
 
805
    dfValue = 1;
 
806
    CPL_LSBPTR64( &dfValue );
 
807
    memcpy( achHeader +  8*16 + 8, &dfValue, 8 );
 
808
    
 
809
    memcpy( achHeader +  9*16, "LONG_INC", 8 );
 
810
    memcpy( achHeader +  9*16 + 8, &dfValue, 8 );
 
811
    
 
812
    memcpy( achHeader + 10*16, "GS_COUNT", 8 );
 
813
    GUInt32 nGSCount = nXSize * nYSize;
 
814
    CPL_LSBPTR32( &nGSCount );
 
815
    memcpy( achHeader + 10*16+8, &nGSCount, 4 );
 
816
    
 
817
    VSIFWriteL( achHeader, 1, sizeof(achHeader), fp );
 
818
 
 
819
/* -------------------------------------------------------------------- */
 
820
/*      Write zeroed grid data.                                         */
 
821
/* -------------------------------------------------------------------- */
 
822
    int i;
 
823
 
 
824
    memset( achHeader, 0, 16 );
 
825
 
 
826
    // Use -1 (0x000080bf) as the default error value.
 
827
    memset( achHeader + 10, 0x80, 1 );
 
828
    memset( achHeader + 11, 0xbf, 1 );
 
829
    memset( achHeader + 14, 0x80, 1 );
 
830
    memset( achHeader + 15, 0xbf, 1 );
 
831
 
 
832
    for( i = 0; i < nXSize * nYSize; i++ )
 
833
        VSIFWriteL( achHeader, 1, 16, fp );
 
834
    
 
835
/* -------------------------------------------------------------------- */
 
836
/*      Write the end record.                                           */
 
837
/* -------------------------------------------------------------------- */
 
838
    memset( achHeader, 0, 16 );
 
839
    memcpy( achHeader, "END     ", 8 );
 
840
    VSIFWriteL( achHeader, 1, 16, fp );
 
841
    VSIFCloseL( fp );
 
842
 
 
843
    if( nNumFile == 1 )
 
844
        return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
 
845
    else
 
846
    {
 
847
        CPLString osSubDSName;
 
848
        osSubDSName.Printf( "NTv2:%d:%s", nNumFile-1, pszFilename );
 
849
        return (GDALDataset *) GDALOpen( osSubDSName, GA_Update );
 
850
    }
 
851
}
 
852
 
 
853
/************************************************************************/
 
854
/*                         GDALRegister_NTv2()                          */
 
855
/************************************************************************/
 
856
 
 
857
void GDALRegister_NTv2()
 
858
 
 
859
{
 
860
    GDALDriver  *poDriver;
 
861
 
 
862
    if( GDALGetDriverByName( "NTv2" ) == NULL )
 
863
    {
 
864
        poDriver = new GDALDriver();
 
865
        
 
866
        poDriver->SetDescription( "NTv2" );
 
867
        poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, 
 
868
                                   "NTv2 Datum Grid Shift" );
 
869
        poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "gsb" );
 
870
        poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
 
871
 
 
872
        poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES, 
 
873
                                   "Float32" );
 
874
 
 
875
        poDriver->pfnOpen = NTv2Dataset::Open;
 
876
        poDriver->pfnIdentify = NTv2Dataset::Identify;
 
877
        poDriver->pfnCreate = NTv2Dataset::Create;
 
878
 
 
879
        GetGDALDriverManager()->RegisterDriver( poDriver );
 
880
    }
 
881
}
 
882