1 : /******************************************************************************
2 : * $Id: ntv2dataset.cpp 25283 2012-12-03 21:01:15Z rouault $
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 25283 2012-12-03 21:01:15Z rouault $");
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 15 : NTv2Dataset::NTv2Dataset()
124 : {
125 15 : fpImage = NULL;
126 15 : }
127 :
128 : /************************************************************************/
129 : /* ~NTv2Dataset() */
130 : /************************************************************************/
131 :
132 15 : NTv2Dataset::~NTv2Dataset()
133 :
134 : {
135 15 : FlushCache();
136 :
137 15 : if( fpImage != NULL )
138 15 : VSIFCloseL( fpImage );
139 15 : }
140 :
141 : /************************************************************************/
142 : /* FlushCache() */
143 : /************************************************************************/
144 :
145 15 : 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 15 : if( eAccess != GA_Update || !(GetPamFlags() & GPF_DIRTY) )
153 : {
154 11 : RawDataset::FlushCache();
155 11 : return;
156 : }
157 :
158 : /* -------------------------------------------------------------------- */
159 : /* Load grid and file headers. */
160 : /* -------------------------------------------------------------------- */
161 : char achFileHeader[11*16];
162 : char achGridHeader[11*16];
163 :
164 4 : VSIFSeekL( fpImage, 0, SEEK_SET );
165 4 : VSIFReadL( achFileHeader, 11, 16, fpImage );
166 :
167 4 : VSIFSeekL( fpImage, nGridOffset, SEEK_SET );
168 4 : 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 4 : char **papszMD = GetMetadata();
176 : int i;
177 4 : int bSomeLeftOver = FALSE;
178 :
179 41 : for( i = 0; papszMD != NULL && papszMD[i] != NULL; i++ )
180 : {
181 37 : char *pszKey = NULL;
182 37 : const char *pszValue = CPLParseNameValue( papszMD[i], &pszKey );
183 37 : if( pszKey == NULL )
184 0 : continue;
185 :
186 37 : if( EQUAL(pszKey,"GS_TYPE") )
187 : {
188 3 : memcpy( achFileHeader + 3*16+8, " ", 8 );
189 3 : memcpy( achFileHeader + 3*16+8, pszValue, MIN(8,strlen(pszValue)) );
190 : }
191 34 : else if( EQUAL(pszKey,"VERSION") )
192 : {
193 3 : memcpy( achFileHeader + 4*16+8, " ", 8 );
194 3 : memcpy( achFileHeader + 4*16+8, pszValue, MIN(8,strlen(pszValue)) );
195 : }
196 31 : else if( EQUAL(pszKey,"SYSTEM_F") )
197 : {
198 3 : memcpy( achFileHeader + 5*16+8, " ", 8 );
199 3 : memcpy( achFileHeader + 5*16+8, pszValue, MIN(8,strlen(pszValue)) );
200 : }
201 28 : else if( EQUAL(pszKey,"SYSTEM_T") )
202 : {
203 3 : memcpy( achFileHeader + 6*16+8, " ", 8 );
204 3 : memcpy( achFileHeader + 6*16+8, pszValue, MIN(8,strlen(pszValue)) );
205 : }
206 25 : else if( EQUAL(pszKey,"MAJOR_F") )
207 : {
208 3 : double dfValue = atof(pszValue);
209 : CPL_LSBPTR64( &dfValue );
210 3 : memcpy( achFileHeader + 7*16+8, &dfValue, 8 );
211 : }
212 22 : else if( EQUAL(pszKey,"MINOR_F") )
213 : {
214 3 : double dfValue = atof(pszValue);
215 : CPL_LSBPTR64( &dfValue );
216 3 : memcpy( achFileHeader + 8*16+8, &dfValue, 8 );
217 : }
218 19 : else if( EQUAL(pszKey,"MAJOR_T") )
219 : {
220 3 : double dfValue = atof(pszValue);
221 : CPL_LSBPTR64( &dfValue );
222 3 : memcpy( achFileHeader + 9*16+8, &dfValue, 8 );
223 : }
224 16 : else if( EQUAL(pszKey,"MINOR_T") )
225 : {
226 3 : double dfValue = atof(pszValue);
227 : CPL_LSBPTR64( &dfValue );
228 3 : memcpy( achFileHeader + 10*16+8, &dfValue, 8 );
229 : }
230 13 : else if( EQUAL(pszKey,"SUB_NAME") )
231 : {
232 3 : memcpy( achGridHeader + 0*16+8, " ", 8 );
233 3 : memcpy( achGridHeader + 0*16+8, pszValue, MIN(8,strlen(pszValue)) );
234 : }
235 10 : else if( EQUAL(pszKey,"PARENT") )
236 : {
237 3 : memcpy( achGridHeader + 1*16+8, " ", 8 );
238 3 : memcpy( achGridHeader + 1*16+8, pszValue, MIN(8,strlen(pszValue)) );
239 : }
240 7 : else if( EQUAL(pszKey,"CREATED") )
241 : {
242 3 : memcpy( achGridHeader + 2*16+8, " ", 8 );
243 3 : memcpy( achGridHeader + 2*16+8, pszValue, MIN(8,strlen(pszValue)) );
244 : }
245 4 : else if( EQUAL(pszKey,"UPDATED") )
246 : {
247 3 : memcpy( achGridHeader + 3*16+8, " ", 8 );
248 3 : memcpy( achGridHeader + 3*16+8, pszValue, MIN(8,strlen(pszValue)) );
249 : }
250 : else
251 : {
252 1 : bSomeLeftOver = TRUE;
253 : }
254 :
255 37 : CPLFree( pszKey );
256 : }
257 :
258 : /* -------------------------------------------------------------------- */
259 : /* Load grid and file headers. */
260 : /* -------------------------------------------------------------------- */
261 4 : VSIFSeekL( fpImage, 0, SEEK_SET );
262 4 : VSIFWriteL( achFileHeader, 11, 16, fpImage );
263 :
264 4 : VSIFSeekL( fpImage, nGridOffset, SEEK_SET );
265 4 : VSIFWriteL( achGridHeader, 11, 16, fpImage );
266 :
267 : /* -------------------------------------------------------------------- */
268 : /* Clear flags if we got everything, then let pam and below do */
269 : /* their flushing. */
270 : /* -------------------------------------------------------------------- */
271 4 : if( !bSomeLeftOver )
272 3 : SetPamFlags( GetPamFlags() & (~GPF_DIRTY) );
273 :
274 4 : RawDataset::FlushCache();
275 : }
276 :
277 : /************************************************************************/
278 : /* Identify() */
279 : /************************************************************************/
280 :
281 11944 : int NTv2Dataset::Identify( GDALOpenInfo *poOpenInfo )
282 :
283 : {
284 11944 : if( EQUALN(poOpenInfo->pszFilename,"NTv2:",5) )
285 0 : return TRUE;
286 :
287 11944 : if( poOpenInfo->nHeaderBytes < 64 )
288 11360 : return FALSE;
289 :
290 584 : if( !EQUALN((const char *)poOpenInfo->pabyHeader + 0, "NUM_OREC", 8 ) )
291 569 : return FALSE;
292 :
293 15 : if( !EQUALN((const char *)poOpenInfo->pabyHeader +16, "NUM_SREC", 8 ) )
294 0 : return FALSE;
295 :
296 15 : return TRUE;
297 : }
298 :
299 : /************************************************************************/
300 : /* Open() */
301 : /************************************************************************/
302 :
303 1851 : GDALDataset *NTv2Dataset::Open( GDALOpenInfo * poOpenInfo )
304 :
305 : {
306 1851 : if( !Identify( poOpenInfo ) )
307 1836 : return NULL;
308 :
309 : /* -------------------------------------------------------------------- */
310 : /* Are we targetting a particular grid? */
311 : /* -------------------------------------------------------------------- */
312 15 : CPLString osFilename;
313 15 : int iTargetGrid = -1;
314 :
315 15 : if( EQUALN(poOpenInfo->pszFilename,"NTv2:",5) )
316 : {
317 0 : const char *pszRest = poOpenInfo->pszFilename+5;
318 :
319 0 : iTargetGrid = atoi(pszRest);
320 0 : while( *pszRest != '\0' && *pszRest != ':' )
321 0 : pszRest++;
322 :
323 0 : if( *pszRest == ':' )
324 0 : pszRest++;
325 :
326 0 : osFilename = pszRest;
327 : }
328 : else
329 15 : osFilename = poOpenInfo->pszFilename;
330 :
331 : /* -------------------------------------------------------------------- */
332 : /* Create a corresponding GDALDataset. */
333 : /* -------------------------------------------------------------------- */
334 : NTv2Dataset *poDS;
335 :
336 15 : poDS = new NTv2Dataset();
337 15 : poDS->eAccess = poOpenInfo->eAccess;
338 :
339 : /* -------------------------------------------------------------------- */
340 : /* Open the file. */
341 : /* -------------------------------------------------------------------- */
342 15 : if( poOpenInfo->eAccess == GA_ReadOnly )
343 10 : poDS->fpImage = VSIFOpenL( osFilename, "rb" );
344 : else
345 5 : poDS->fpImage = VSIFOpenL( osFilename, "rb+" );
346 :
347 15 : if( poDS->fpImage == NULL )
348 : {
349 0 : delete poDS;
350 0 : return NULL;
351 : }
352 :
353 : /* -------------------------------------------------------------------- */
354 : /* Read the file header. */
355 : /* -------------------------------------------------------------------- */
356 : char achHeader[11*16];
357 : GInt32 nSubFileCount;
358 : double dfValue;
359 15 : CPLString osFValue;
360 :
361 15 : VSIFSeekL( poDS->fpImage, 0, SEEK_SET );
362 15 : VSIFReadL( achHeader, 11, 16, poDS->fpImage );
363 :
364 : CPL_LSBPTR32( achHeader + 2*16 + 8 );
365 15 : memcpy( &nSubFileCount, achHeader + 2*16 + 8, 4 );
366 15 : if (nSubFileCount <= 0 || nSubFileCount >= 1024)
367 : {
368 : CPLError(CE_Failure, CPLE_AppDefined,
369 0 : "Invalid value for NUM_FILE : %d", nSubFileCount);
370 0 : delete poDS;
371 0 : return NULL;
372 : }
373 :
374 15 : poDS->CaptureMetadataItem( achHeader + 3*16 );
375 15 : poDS->CaptureMetadataItem( achHeader + 4*16 );
376 15 : poDS->CaptureMetadataItem( achHeader + 5*16 );
377 15 : poDS->CaptureMetadataItem( achHeader + 6*16 );
378 :
379 15 : memcpy( &dfValue, achHeader + 7*16 + 8, 8 );
380 : CPL_LSBPTR64( &dfValue );
381 15 : osFValue.Printf( "%.15g", dfValue );
382 15 : poDS->SetMetadataItem( "MAJOR_F", osFValue );
383 :
384 15 : memcpy( &dfValue, achHeader + 8*16 + 8, 8 );
385 : CPL_LSBPTR64( &dfValue );
386 15 : osFValue.Printf( "%.15g", dfValue );
387 15 : poDS->SetMetadataItem( "MINOR_F", osFValue );
388 :
389 15 : memcpy( &dfValue, achHeader + 9*16 + 8, 8 );
390 : CPL_LSBPTR64( &dfValue );
391 15 : osFValue.Printf( "%.15g", dfValue );
392 15 : poDS->SetMetadataItem( "MAJOR_T", osFValue );
393 :
394 15 : memcpy( &dfValue, achHeader + 10*16 + 8, 8 );
395 : CPL_LSBPTR64( &dfValue );
396 15 : osFValue.Printf( "%.15g", dfValue );
397 15 : poDS->SetMetadataItem( "MINOR_T", osFValue );
398 :
399 : /* ==================================================================== */
400 : /* Loop over grids. */
401 : /* ==================================================================== */
402 : int iGrid;
403 15 : vsi_l_offset nGridOffset = sizeof(achHeader);
404 :
405 15 : for( iGrid = 0; iGrid < nSubFileCount; iGrid++ )
406 : {
407 15 : CPLString osSubName;
408 : int i;
409 : GUInt32 nGSCount;
410 :
411 15 : VSIFSeekL( poDS->fpImage, nGridOffset, SEEK_SET );
412 15 : if (VSIFReadL( achHeader, 11, 16, poDS->fpImage ) != 16)
413 : {
414 : CPLError(CE_Failure, CPLE_AppDefined,
415 0 : "Cannot read header for subfile %d", iGrid);
416 0 : delete poDS;
417 0 : return NULL;
418 : }
419 :
420 15 : for( i = 4; i <= 9; i++ )
421 : CPL_LSBPTR64( achHeader + i*16 + 8 );
422 :
423 : CPL_LSBPTR32( achHeader + 10*16 + 8 );
424 :
425 15 : memcpy( &nGSCount, achHeader + 10*16 + 8, 4 );
426 :
427 15 : osSubName.assign( achHeader + 8, 8 );
428 15 : osSubName.Trim();
429 :
430 : // If this is our target grid, open it as a dataset.
431 15 : if( iTargetGrid == iGrid || (iTargetGrid == -1 && iGrid == 0) )
432 : {
433 15 : if( !poDS->OpenGrid( achHeader, nGridOffset ) )
434 : {
435 0 : delete poDS;
436 0 : return NULL;
437 : }
438 : }
439 :
440 : // If we are opening the file as a whole, list subdatasets.
441 15 : if( iTargetGrid == -1 )
442 : {
443 15 : CPLString osKey, osValue;
444 :
445 15 : osKey.Printf( "SUBDATASET_%d_NAME", iGrid );
446 15 : osValue.Printf( "NTv2:%d:%s", iGrid, osFilename.c_str() );
447 15 : poDS->SetMetadataItem( osKey, osValue, "SUBDATASETS" );
448 :
449 15 : osKey.Printf( "SUBDATASET_%d_DESC", iGrid );
450 15 : osValue.Printf( "%s", osSubName.c_str() );
451 15 : poDS->SetMetadataItem( osKey, osValue, "SUBDATASETS" );
452 : }
453 :
454 15 : nGridOffset += (11+(vsi_l_offset)nGSCount) * 16;
455 : }
456 :
457 : /* -------------------------------------------------------------------- */
458 : /* Initialize any PAM information. */
459 : /* -------------------------------------------------------------------- */
460 15 : poDS->SetDescription( poOpenInfo->pszFilename );
461 15 : poDS->TryLoadXML();
462 :
463 : /* -------------------------------------------------------------------- */
464 : /* Check for overviews. */
465 : /* -------------------------------------------------------------------- */
466 15 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
467 :
468 15 : return( poDS );
469 : }
470 :
471 : /************************************************************************/
472 : /* OpenGrid() */
473 : /* */
474 : /* Note that the caller will already have byte swapped needed */
475 : /* portions of the header. */
476 : /************************************************************************/
477 :
478 15 : int NTv2Dataset::OpenGrid( char *pachHeader, vsi_l_offset nGridOffset )
479 :
480 : {
481 15 : this->nGridOffset = nGridOffset;
482 :
483 : /* -------------------------------------------------------------------- */
484 : /* Read the grid header. */
485 : /* -------------------------------------------------------------------- */
486 : double s_lat, n_lat, e_long, w_long, lat_inc, long_inc;
487 :
488 15 : CaptureMetadataItem( pachHeader + 0*16 );
489 15 : CaptureMetadataItem( pachHeader + 1*16 );
490 15 : CaptureMetadataItem( pachHeader + 2*16 );
491 15 : CaptureMetadataItem( pachHeader + 3*16 );
492 :
493 15 : memcpy( &s_lat, pachHeader + 4*16 + 8, 8 );
494 15 : memcpy( &n_lat, pachHeader + 5*16 + 8, 8 );
495 15 : memcpy( &e_long, pachHeader + 6*16 + 8, 8 );
496 15 : memcpy( &w_long, pachHeader + 7*16 + 8, 8 );
497 15 : memcpy( &lat_inc, pachHeader + 8*16 + 8, 8 );
498 15 : memcpy( &long_inc, pachHeader + 9*16 + 8, 8 );
499 :
500 15 : e_long *= -1;
501 15 : w_long *= -1;
502 :
503 15 : nRasterXSize = (int) floor((e_long - w_long) / long_inc + 1.5);
504 15 : nRasterYSize = (int) floor((n_lat - s_lat) / lat_inc + 1.5);
505 :
506 15 : if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize))
507 0 : return FALSE;
508 :
509 : /* -------------------------------------------------------------------- */
510 : /* Create band information object. */
511 : /* */
512 : /* We use unusual offsets to remap from bottom to top, to top */
513 : /* to bottom orientation, and also to remap east to west, to */
514 : /* west to east. */
515 : /* -------------------------------------------------------------------- */
516 : int iBand;
517 :
518 150 : for( iBand = 0; iBand < 4; iBand++ )
519 : {
520 : RawRasterBand *poBand =
521 : new RawRasterBand( this, iBand+1, fpImage,
522 : nGridOffset + 4*iBand + 11*16
523 : + (nRasterXSize-1) * 16
524 : + (nRasterYSize-1) * 16 * nRasterXSize,
525 : -16, -16 * nRasterXSize,
526 60 : GDT_Float32, CPL_IS_LSB, TRUE, FALSE );
527 60 : SetBand( iBand+1, poBand );
528 : }
529 :
530 15 : GetRasterBand(1)->SetDescription( "Latitude Offset (arc seconds)" );
531 15 : GetRasterBand(2)->SetDescription( "Longitude Offset (arc seconds)" );
532 15 : GetRasterBand(3)->SetDescription( "Latitude Error" );
533 15 : GetRasterBand(4)->SetDescription( "Longitude Error" );
534 :
535 : /* -------------------------------------------------------------------- */
536 : /* Setup georeferencing. */
537 : /* -------------------------------------------------------------------- */
538 15 : adfGeoTransform[0] = (w_long - long_inc*0.5) / 3600.0;
539 15 : adfGeoTransform[1] = long_inc / 3600.0;
540 15 : adfGeoTransform[2] = 0.0;
541 15 : adfGeoTransform[3] = (n_lat + lat_inc*0.5) / 3600.0;
542 15 : adfGeoTransform[4] = 0.0;
543 15 : adfGeoTransform[5] = (-1 * lat_inc) / 3600.0;
544 :
545 15 : return TRUE;
546 : }
547 :
548 : /************************************************************************/
549 : /* CaptureMetadataItem() */
550 : /************************************************************************/
551 :
552 120 : void NTv2Dataset::CaptureMetadataItem( char *pszItem )
553 :
554 : {
555 120 : CPLString osKey, osValue;
556 :
557 120 : osKey.assign( pszItem, 8 );
558 120 : osValue.assign( pszItem+8, 8 );
559 :
560 120 : SetMetadataItem( osKey.Trim(), osValue.Trim() );
561 120 : }
562 :
563 : /************************************************************************/
564 : /* GetGeoTransform() */
565 : /************************************************************************/
566 :
567 5 : CPLErr NTv2Dataset::GetGeoTransform( double * padfTransform )
568 :
569 : {
570 5 : memcpy( padfTransform, adfGeoTransform, sizeof(double)*6 );
571 5 : return CE_None;
572 : }
573 :
574 : /************************************************************************/
575 : /* SetGeoTransform() */
576 : /************************************************************************/
577 :
578 4 : CPLErr NTv2Dataset::SetGeoTransform( double * padfTransform )
579 :
580 : {
581 4 : if( eAccess == GA_ReadOnly )
582 : {
583 : CPLError( CE_Failure, CPLE_NoWriteAccess,
584 0 : "Unable to update geotransform on readonly file." );
585 0 : return CE_Failure;
586 : }
587 :
588 4 : if( padfTransform[2] != 0.0 || padfTransform[4] != 0.0 )
589 : {
590 : CPLError( CE_Failure, CPLE_AppDefined,
591 0 : "Rotated and sheared geotransforms not supported for NTv2.");
592 0 : return CE_Failure;
593 : }
594 :
595 4 : memcpy( adfGeoTransform, padfTransform, sizeof(double)*6 );
596 :
597 : /* -------------------------------------------------------------------- */
598 : /* Update grid header. */
599 : /* -------------------------------------------------------------------- */
600 : double dfValue;
601 : char achHeader[11*16];
602 :
603 : // read grid header
604 4 : VSIFSeekL( fpImage, nGridOffset, SEEK_SET );
605 4 : VSIFReadL( achHeader, 11, 16, fpImage );
606 :
607 : // S_LAT
608 4 : dfValue = 3600 * (adfGeoTransform[3] + (nRasterYSize-0.5) * adfGeoTransform[5]);
609 : CPL_LSBPTR64( &dfValue );
610 4 : memcpy( achHeader + 4*16 + 8, &dfValue, 8 );
611 :
612 : // N_LAT
613 4 : dfValue = 3600 * (adfGeoTransform[3] + 0.5 * adfGeoTransform[5]);
614 : CPL_LSBPTR64( &dfValue );
615 4 : memcpy( achHeader + 5*16 + 8, &dfValue, 8 );
616 :
617 : // E_LONG
618 4 : dfValue = -3600 * (adfGeoTransform[0] + (nRasterXSize-0.5)*adfGeoTransform[1]);
619 : CPL_LSBPTR64( &dfValue );
620 4 : memcpy( achHeader + 6*16 + 8, &dfValue, 8 );
621 :
622 : // W_LONG
623 4 : dfValue = -3600 * (adfGeoTransform[0] + 0.5 * adfGeoTransform[1]);
624 : CPL_LSBPTR64( &dfValue );
625 4 : memcpy( achHeader + 7*16 + 8, &dfValue, 8 );
626 :
627 : // LAT_INC
628 4 : dfValue = -3600 * adfGeoTransform[5];
629 : CPL_LSBPTR64( &dfValue );
630 4 : memcpy( achHeader + 8*16 + 8, &dfValue, 8 );
631 :
632 : // LONG_INC
633 4 : dfValue = 3600 * adfGeoTransform[1];
634 : CPL_LSBPTR64( &dfValue );
635 4 : memcpy( achHeader + 9*16 + 8, &dfValue, 8 );
636 :
637 : // write grid header.
638 4 : VSIFSeekL( fpImage, nGridOffset, SEEK_SET );
639 4 : VSIFWriteL( achHeader, 11, 16, fpImage );
640 :
641 4 : return CE_None;
642 : }
643 :
644 :
645 : /************************************************************************/
646 : /* GetProjectionRef() */
647 : /************************************************************************/
648 :
649 5 : const char *NTv2Dataset::GetProjectionRef()
650 :
651 : {
652 5 : return SRS_WKT_WGS84;
653 : }
654 :
655 : /************************************************************************/
656 : /* Create() */
657 : /************************************************************************/
658 :
659 45 : GDALDataset *NTv2Dataset::Create( const char * pszFilename,
660 : int nXSize, int nYSize, int nBands,
661 : GDALDataType eType,
662 : char ** papszOptions )
663 :
664 : {
665 45 : if( eType != GDT_Float32 )
666 : {
667 : CPLError(CE_Failure, CPLE_AppDefined,
668 : "Attempt to create NTv2 file with unsupported data type '%s'.",
669 38 : GDALGetDataTypeName( eType ) );
670 38 : return NULL;
671 : }
672 :
673 : /* -------------------------------------------------------------------- */
674 : /* Are we extending an existing file? */
675 : /* -------------------------------------------------------------------- */
676 : VSILFILE *fp;
677 7 : GUInt32 nNumFile = 1;
678 :
679 7 : int bAppend = CSLFetchBoolean(papszOptions,"APPEND_SUBDATASET",FALSE);
680 :
681 : /* -------------------------------------------------------------------- */
682 : /* Try to open or create file. */
683 : /* -------------------------------------------------------------------- */
684 7 : if( bAppend )
685 0 : fp = VSIFOpenL( pszFilename, "rb+" );
686 : else
687 7 : fp = VSIFOpenL( pszFilename, "wb" );
688 :
689 7 : if( fp == NULL )
690 : {
691 : CPLError( CE_Failure, CPLE_OpenFailed,
692 : "Attempt to open/create file `%s' failed.\n",
693 2 : pszFilename );
694 2 : return NULL;
695 : }
696 :
697 : /* -------------------------------------------------------------------- */
698 : /* Create a file level header if we are creating new. */
699 : /* -------------------------------------------------------------------- */
700 : char achHeader[11*16];
701 : const char *pszValue;
702 :
703 5 : if( !bAppend )
704 : {
705 5 : memset( achHeader, 0, sizeof(achHeader) );
706 :
707 5 : memcpy( achHeader + 0*16, "NUM_OREC", 8 );
708 5 : achHeader[ 0*16 + 8] = 0xb;
709 :
710 5 : memcpy( achHeader + 1*16, "NUM_SREC", 8 );
711 5 : achHeader[ 1*16 + 8] = 0xb;
712 :
713 5 : memcpy( achHeader + 2*16, "NUM_FILE", 8 );
714 5 : achHeader[ 2*16 + 8] = 0x1;
715 :
716 5 : memcpy( achHeader + 3*16, "GS_TYPE ", 16 );
717 5 : pszValue = CSLFetchNameValueDef( papszOptions, "GS_TYPE", "SECONDS");
718 5 : memcpy( achHeader + 3*16+8, pszValue, MIN(16,strlen(pszValue)) );
719 :
720 5 : memcpy( achHeader + 4*16, "VERSION ", 16 );
721 5 : pszValue = CSLFetchNameValueDef( papszOptions, "VERSION", "" );
722 5 : memcpy( achHeader + 4*16+8, pszValue, MIN(16,strlen(pszValue)) );
723 :
724 5 : memcpy( achHeader + 5*16, "SYSTEM_F ", 16 );
725 5 : pszValue = CSLFetchNameValueDef( papszOptions, "SYSTEM_F", "" );
726 5 : memcpy( achHeader + 5*16+8, pszValue, MIN(16,strlen(pszValue)) );
727 :
728 5 : memcpy( achHeader + 6*16, "SYSTEM_T ", 16 );
729 5 : pszValue = CSLFetchNameValueDef( papszOptions, "SYSTEM_T", "" );
730 5 : memcpy( achHeader + 6*16+8, pszValue, MIN(16,strlen(pszValue)) );
731 :
732 5 : memcpy( achHeader + 7*16, "MAJOR_F ", 8);
733 5 : memcpy( achHeader + 8*16, "MINOR_F ", 8 );
734 5 : memcpy( achHeader + 9*16, "MAJOR_T ", 8 );
735 5 : memcpy( achHeader + 10*16, "MINOR_T ", 8 );
736 :
737 5 : VSIFWriteL( achHeader, 1, sizeof(achHeader), fp );
738 : }
739 :
740 : /* -------------------------------------------------------------------- */
741 : /* Otherwise update the header with an increased subfile count, */
742 : /* and advanced to the last record of the file. */
743 : /* -------------------------------------------------------------------- */
744 : else
745 : {
746 0 : VSIFSeekL( fp, 2*16 + 8, SEEK_SET );
747 0 : VSIFReadL( &nNumFile, 1, 4, fp );
748 : CPL_LSBPTR32( &nNumFile );
749 :
750 0 : nNumFile++;
751 :
752 : CPL_LSBPTR32( &nNumFile );
753 0 : VSIFSeekL( fp, 2*16 + 8, SEEK_SET );
754 0 : VSIFWriteL( &nNumFile, 1, 4, fp );
755 :
756 : vsi_l_offset nEnd;
757 :
758 0 : VSIFSeekL( fp, 0, SEEK_END );
759 0 : nEnd = VSIFTellL( fp );
760 0 : VSIFSeekL( fp, nEnd-16, SEEK_SET );
761 : }
762 :
763 : /* -------------------------------------------------------------------- */
764 : /* Write the grid header. */
765 : /* -------------------------------------------------------------------- */
766 5 : memset( achHeader, 0, sizeof(achHeader) );
767 :
768 5 : memcpy( achHeader + 0*16, "SUB_NAME ", 16 );
769 5 : pszValue = CSLFetchNameValueDef( papszOptions, "SUB_NAME", "" );
770 5 : memcpy( achHeader + 0*16+8, pszValue, MIN(16,strlen(pszValue)) );
771 :
772 5 : memcpy( achHeader + 1*16, "PARENT ", 16 );
773 5 : pszValue = CSLFetchNameValueDef( papszOptions, "PARENT", "NONE" );
774 5 : memcpy( achHeader + 1*16+8, pszValue, MIN(16,strlen(pszValue)) );
775 :
776 5 : memcpy( achHeader + 2*16, "CREATED ", 16 );
777 5 : pszValue = CSLFetchNameValueDef( papszOptions, "CREATED", "" );
778 5 : memcpy( achHeader + 2*16+8, pszValue, MIN(16,strlen(pszValue)) );
779 :
780 5 : memcpy( achHeader + 3*16, "UPDATED ", 16 );
781 5 : pszValue = CSLFetchNameValueDef( papszOptions, "UPDATED", "" );
782 5 : memcpy( achHeader + 3*16+8, pszValue, MIN(16,strlen(pszValue)) );
783 :
784 : double dfValue;
785 :
786 5 : memcpy( achHeader + 4*16, "S_LAT ", 8 );
787 5 : dfValue = 0;
788 : CPL_LSBPTR64( &dfValue );
789 5 : memcpy( achHeader + 4*16 + 8, &dfValue, 8 );
790 :
791 5 : memcpy( achHeader + 5*16, "N_LAT ", 8 );
792 5 : dfValue = nYSize-1;
793 : CPL_LSBPTR64( &dfValue );
794 5 : memcpy( achHeader + 5*16 + 8, &dfValue, 8 );
795 :
796 5 : memcpy( achHeader + 6*16, "E_LONG ", 8 );
797 5 : dfValue = -1*(nXSize-1);
798 : CPL_LSBPTR64( &dfValue );
799 5 : memcpy( achHeader + 6*16 + 8, &dfValue, 8 );
800 :
801 5 : memcpy( achHeader + 7*16, "W_LONG ", 8 );
802 5 : dfValue = 0;
803 : CPL_LSBPTR64( &dfValue );
804 5 : memcpy( achHeader + 7*16 + 8, &dfValue, 8 );
805 :
806 5 : memcpy( achHeader + 8*16, "LAT_INC ", 8 );
807 5 : dfValue = 1;
808 : CPL_LSBPTR64( &dfValue );
809 5 : memcpy( achHeader + 8*16 + 8, &dfValue, 8 );
810 :
811 5 : memcpy( achHeader + 9*16, "LONG_INC", 8 );
812 5 : memcpy( achHeader + 9*16 + 8, &dfValue, 8 );
813 :
814 5 : memcpy( achHeader + 10*16, "GS_COUNT", 8 );
815 5 : GUInt32 nGSCount = nXSize * nYSize;
816 : CPL_LSBPTR32( &nGSCount );
817 5 : memcpy( achHeader + 10*16+8, &nGSCount, 4 );
818 :
819 5 : VSIFWriteL( achHeader, 1, sizeof(achHeader), fp );
820 :
821 : /* -------------------------------------------------------------------- */
822 : /* Write zeroed grid data. */
823 : /* -------------------------------------------------------------------- */
824 : int i;
825 :
826 5 : memset( achHeader, 0, 16 );
827 :
828 : // Use -1 (0x000080bf) as the default error value.
829 5 : memset( achHeader + 10, 0x80, 1 );
830 5 : memset( achHeader + 11, 0xbf, 1 );
831 5 : memset( achHeader + 14, 0x80, 1 );
832 5 : memset( achHeader + 15, 0xbf, 1 );
833 :
834 59867 : for( i = 0; i < nXSize * nYSize; i++ )
835 59862 : VSIFWriteL( achHeader, 1, 16, fp );
836 :
837 : /* -------------------------------------------------------------------- */
838 : /* Write the end record. */
839 : /* -------------------------------------------------------------------- */
840 5 : memset( achHeader, 0, 16 );
841 5 : memcpy( achHeader, "END ", 8 );
842 5 : VSIFWriteL( achHeader, 1, 16, fp );
843 5 : VSIFCloseL( fp );
844 :
845 5 : if( nNumFile == 1 )
846 5 : return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
847 : else
848 : {
849 0 : CPLString osSubDSName;
850 0 : osSubDSName.Printf( "NTv2:%d:%s", nNumFile-1, pszFilename );
851 0 : return (GDALDataset *) GDALOpen( osSubDSName, GA_Update );
852 : }
853 : }
854 :
855 : /************************************************************************/
856 : /* GDALRegister_NTv2() */
857 : /************************************************************************/
858 :
859 582 : void GDALRegister_NTv2()
860 :
861 : {
862 : GDALDriver *poDriver;
863 :
864 582 : if( GDALGetDriverByName( "NTv2" ) == NULL )
865 : {
866 561 : poDriver = new GDALDriver();
867 :
868 561 : poDriver->SetDescription( "NTv2" );
869 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
870 561 : "NTv2 Datum Grid Shift" );
871 561 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "gsb" );
872 561 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
873 :
874 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
875 561 : "Float32" );
876 :
877 561 : poDriver->pfnOpen = NTv2Dataset::Open;
878 561 : poDriver->pfnIdentify = NTv2Dataset::Identify;
879 561 : poDriver->pfnCreate = NTv2Dataset::Create;
880 :
881 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
882 : }
883 582 : }
884 :
|
