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