1 : /******************************************************************************
2 : * $Id$
3 : *
4 : * Project: Horizontal Datum Formats
5 : * Purpose: Implementation of the CTable2 format, a PROJ.4 specific format
6 : * that is more compact (due to a lack of unused error band) than NTv2
7 : * Author: Frank Warmerdam, warmerdam@pobox.com
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 2012, Frank Warmerdam
11 : *
12 : * Permission is hereby granted, free of charge, to any person obtaining a
13 : * copy of this software and associated documentation files (the "Software"),
14 : * to deal in the Software without restriction, including without limitation
15 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : * and/or sell copies of the Software, and to permit persons to whom the
17 : * Software is furnished to do so, subject to the following conditions:
18 : *
19 : * The above copyright notice and this permission notice shall be included
20 : * in all copies or substantial portions of the Software.
21 : *
22 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 : * DEALINGS IN THE SOFTWARE.
29 : ****************************************************************************/
30 :
31 : #include "rawdataset.h"
32 : #include "cpl_string.h"
33 : #include "ogr_srs_api.h"
34 :
35 : CPL_CVSID("$Id$");
36 :
37 : #ifndef M_PI
38 : #define M_PI 3.14159265358979323846
39 : #endif
40 :
41 : /************************************************************************/
42 : /* ==================================================================== */
43 : /* CTable2Dataset */
44 : /* ==================================================================== */
45 : /************************************************************************/
46 :
47 : class CTable2Dataset : public RawDataset
48 : {
49 : public:
50 : VSILFILE *fpImage; // image data file.
51 :
52 : double adfGeoTransform[6];
53 :
54 : public:
55 : CTable2Dataset();
56 : ~CTable2Dataset();
57 :
58 : virtual CPLErr SetGeoTransform( double * padfTransform );
59 : virtual CPLErr GetGeoTransform( double * padfTransform );
60 : virtual const char *GetProjectionRef();
61 : virtual void FlushCache(void);
62 :
63 : static GDALDataset *Open( GDALOpenInfo * );
64 : static int Identify( GDALOpenInfo * );
65 : static GDALDataset *Create( const char * pszFilename,
66 : int nXSize, int nYSize, int nBands,
67 : GDALDataType eType, char ** papszOptions );
68 : };
69 :
70 : /************************************************************************/
71 : /* ==================================================================== */
72 : /* CTable2Dataset */
73 : /* ==================================================================== */
74 : /************************************************************************/
75 :
76 : /************************************************************************/
77 : /* CTable2Dataset() */
78 : /************************************************************************/
79 :
80 6 : CTable2Dataset::CTable2Dataset()
81 : {
82 6 : fpImage = NULL;
83 6 : }
84 :
85 : /************************************************************************/
86 : /* ~CTable2Dataset() */
87 : /************************************************************************/
88 :
89 6 : CTable2Dataset::~CTable2Dataset()
90 :
91 : {
92 6 : FlushCache();
93 :
94 6 : if( fpImage != NULL )
95 6 : VSIFCloseL( fpImage );
96 6 : }
97 :
98 : /************************************************************************/
99 : /* FlushCache() */
100 : /************************************************************************/
101 :
102 6 : void CTable2Dataset::FlushCache()
103 :
104 : {
105 6 : RawDataset::FlushCache();
106 6 : }
107 :
108 : /************************************************************************/
109 : /* Identify() */
110 : /************************************************************************/
111 :
112 11929 : int CTable2Dataset::Identify( GDALOpenInfo *poOpenInfo )
113 :
114 : {
115 11929 : if( poOpenInfo->nHeaderBytes < 64 )
116 11360 : return FALSE;
117 :
118 569 : if( !EQUALN((const char *)poOpenInfo->pabyHeader + 0, "CTABLE V2", 9 ) )
119 563 : return FALSE;
120 :
121 6 : return TRUE;
122 : }
123 :
124 : /************************************************************************/
125 : /* Open() */
126 : /************************************************************************/
127 :
128 1836 : GDALDataset *CTable2Dataset::Open( GDALOpenInfo * poOpenInfo )
129 :
130 : {
131 1836 : if( !Identify( poOpenInfo ) )
132 1830 : return NULL;
133 :
134 : /* -------------------------------------------------------------------- */
135 : /* Create a corresponding GDALDataset. */
136 : /* -------------------------------------------------------------------- */
137 : CTable2Dataset *poDS;
138 :
139 6 : poDS = new CTable2Dataset();
140 6 : poDS->eAccess = poOpenInfo->eAccess;
141 :
142 : /* -------------------------------------------------------------------- */
143 : /* Open the file. */
144 : /* -------------------------------------------------------------------- */
145 6 : CPLString osFilename;
146 12 : osFilename = poOpenInfo->pszFilename;
147 :
148 6 : if( poOpenInfo->eAccess == GA_ReadOnly )
149 3 : poDS->fpImage = VSIFOpenL( osFilename, "rb" );
150 : else
151 3 : poDS->fpImage = VSIFOpenL( osFilename, "rb+" );
152 :
153 6 : if( poDS->fpImage == NULL )
154 : {
155 0 : delete poDS;
156 0 : return NULL;
157 : }
158 :
159 : /* -------------------------------------------------------------------- */
160 : /* Read the file header. */
161 : /* -------------------------------------------------------------------- */
162 : char achHeader[160];
163 6 : CPLString osDescription;
164 :
165 6 : VSIFSeekL( poDS->fpImage, 0, SEEK_SET );
166 6 : VSIFReadL( achHeader, 1, 160, poDS->fpImage );
167 :
168 6 : achHeader[16+79] = '\0';
169 6 : osDescription = (const char *) achHeader+16;
170 6 : osDescription.Trim();
171 6 : poDS->SetMetadataItem( "DESCRIPTION", osDescription );
172 :
173 : /* -------------------------------------------------------------------- */
174 : /* Convert from LSB to local machine byte order. */
175 : /* -------------------------------------------------------------------- */
176 : CPL_LSBPTR64( achHeader + 96 );
177 : CPL_LSBPTR64( achHeader + 104 );
178 : CPL_LSBPTR64( achHeader + 112 );
179 : CPL_LSBPTR64( achHeader + 120 );
180 : CPL_LSBPTR32( achHeader + 128 );
181 : CPL_LSBPTR32( achHeader + 132 );
182 :
183 : /* -------------------------------------------------------------------- */
184 : /* Extract size, and geotransform. */
185 : /* -------------------------------------------------------------------- */
186 : GUInt32 nRasterXSize, nRasterYSize;
187 :
188 6 : memcpy( &nRasterXSize, achHeader + 128, 4 );
189 6 : memcpy( &nRasterYSize, achHeader + 132, 4 );
190 6 : poDS->nRasterXSize = nRasterXSize;
191 6 : poDS->nRasterYSize = nRasterYSize;
192 :
193 : int i;
194 : double adfValues[4];
195 6 : memcpy( adfValues, achHeader + 96, sizeof(double)*4 );
196 :
197 30 : for( i = 0; i < 4; i++ )
198 24 : adfValues[i] *= 180/M_PI; // Radians to degrees.
199 :
200 :
201 6 : poDS->adfGeoTransform[0] = adfValues[0] - adfValues[2]*0.5;
202 6 : poDS->adfGeoTransform[1] = adfValues[2];
203 6 : poDS->adfGeoTransform[2] = 0.0;
204 6 : poDS->adfGeoTransform[3] = adfValues[1] + adfValues[3]*(nRasterYSize-0.5);
205 6 : poDS->adfGeoTransform[4] = 0.0;
206 6 : poDS->adfGeoTransform[5] = -adfValues[3];
207 :
208 : /* -------------------------------------------------------------------- */
209 : /* Setup the bands. */
210 : /* -------------------------------------------------------------------- */
211 : RawRasterBand *poBand =
212 : new RawRasterBand( poDS, 1, poDS->fpImage,
213 : 160 + 4 + nRasterXSize * (nRasterYSize-1) * 2 * 4,
214 : 8, -8 * nRasterXSize,
215 6 : GDT_Float32, CPL_IS_LSB, TRUE, FALSE );
216 6 : poBand->SetDescription( "Latitude Offset (radians)" );
217 6 : poDS->SetBand( 1, poBand );
218 :
219 : poBand =
220 : new RawRasterBand( poDS, 2, poDS->fpImage,
221 : 160 + nRasterXSize * (nRasterYSize-1) * 2 * 4,
222 : 8, -8 * nRasterXSize,
223 6 : GDT_Float32, CPL_IS_LSB, TRUE, FALSE );
224 6 : poBand->SetDescription( "Longitude Offset (radians)" );
225 6 : poDS->SetBand( 2, poBand );
226 :
227 : /* -------------------------------------------------------------------- */
228 : /* Initialize any PAM information. */
229 : /* -------------------------------------------------------------------- */
230 6 : poDS->SetDescription( poOpenInfo->pszFilename );
231 6 : poDS->TryLoadXML();
232 :
233 : /* -------------------------------------------------------------------- */
234 : /* Check for overviews. */
235 : /* -------------------------------------------------------------------- */
236 6 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
237 :
238 6 : return( poDS );
239 : }
240 :
241 : /************************************************************************/
242 : /* GetGeoTransform() */
243 : /************************************************************************/
244 :
245 2 : CPLErr CTable2Dataset::GetGeoTransform( double * padfTransform )
246 :
247 : {
248 2 : memcpy( padfTransform, adfGeoTransform, sizeof(double)*6 );
249 2 : return CE_None;
250 : }
251 :
252 : /************************************************************************/
253 : /* SetGeoTransform() */
254 : /************************************************************************/
255 :
256 3 : CPLErr CTable2Dataset::SetGeoTransform( double * padfTransform )
257 :
258 : {
259 3 : if( eAccess == GA_ReadOnly )
260 : {
261 : CPLError( CE_Failure, CPLE_NoWriteAccess,
262 0 : "Unable to update geotransform on readonly file." );
263 0 : return CE_Failure;
264 : }
265 :
266 3 : if( padfTransform[2] != 0.0 || padfTransform[4] != 0.0 )
267 : {
268 : CPLError( CE_Failure, CPLE_AppDefined,
269 0 : "Rotated and sheared geotransforms not supported for CTable2.");
270 0 : return CE_Failure;
271 : }
272 :
273 3 : memcpy( adfGeoTransform, padfTransform, sizeof(double)*6 );
274 :
275 : /* -------------------------------------------------------------------- */
276 : /* Update grid header. */
277 : /* -------------------------------------------------------------------- */
278 : double dfValue;
279 : char achHeader[160];
280 3 : double dfDegToRad = M_PI / 180.0;
281 :
282 : // read grid header
283 3 : VSIFSeekL( fpImage, 0, SEEK_SET );
284 3 : VSIFReadL( achHeader, 1, sizeof(achHeader), fpImage );
285 :
286 : // lower left origin (longitude, center of pixel, radians)
287 3 : dfValue = (adfGeoTransform[0] + adfGeoTransform[1]*0.5) * dfDegToRad;
288 : CPL_LSBPTR64( &dfValue );
289 3 : memcpy( achHeader + 96, &dfValue, 8 );
290 :
291 : // lower left origin (latitude, center of pixel, radians)
292 3 : dfValue = (adfGeoTransform[3] + adfGeoTransform[5] * (nRasterYSize-0.5)) * dfDegToRad;
293 : CPL_LSBPTR64( &dfValue );
294 3 : memcpy( achHeader + 104, &dfValue, 8 );
295 :
296 : // pixel width (radians)
297 3 : dfValue = adfGeoTransform[1] * dfDegToRad;
298 : CPL_LSBPTR64( &dfValue );
299 3 : memcpy( achHeader + 112, &dfValue, 8 );
300 :
301 : // pixel height (radians)
302 3 : dfValue = adfGeoTransform[5] * -1 * dfDegToRad;
303 : CPL_LSBPTR64( &dfValue );
304 3 : memcpy( achHeader + 120, &dfValue, 8 );
305 :
306 : // write grid header.
307 3 : VSIFSeekL( fpImage, 0, SEEK_SET );
308 3 : VSIFWriteL( achHeader, 11, 16, fpImage );
309 :
310 3 : return CE_None;
311 : }
312 :
313 :
314 : /************************************************************************/
315 : /* GetProjectionRef() */
316 : /************************************************************************/
317 :
318 0 : const char *CTable2Dataset::GetProjectionRef()
319 :
320 : {
321 0 : return SRS_WKT_WGS84;
322 : }
323 :
324 : /************************************************************************/
325 : /* Create() */
326 : /************************************************************************/
327 :
328 43 : GDALDataset *CTable2Dataset::Create( const char * pszFilename,
329 : int nXSize, int nYSize, int nBands,
330 : GDALDataType eType,
331 : char ** papszOptions )
332 :
333 : {
334 43 : if( eType != GDT_Float32 )
335 : {
336 : CPLError(CE_Failure, CPLE_AppDefined,
337 : "Attempt to create CTable2 file with unsupported data type '%s'.",
338 40 : GDALGetDataTypeName( eType ) );
339 40 : return NULL;
340 : }
341 :
342 : /* -------------------------------------------------------------------- */
343 : /* Try to open or create file. */
344 : /* -------------------------------------------------------------------- */
345 : VSILFILE *fp;
346 :
347 3 : fp = VSIFOpenL( pszFilename, "wb" );
348 :
349 3 : if( fp == NULL )
350 : {
351 : CPLError( CE_Failure, CPLE_OpenFailed,
352 : "Attempt to create file `%s' failed.\n",
353 0 : pszFilename );
354 0 : return NULL;
355 : }
356 :
357 : /* -------------------------------------------------------------------- */
358 : /* Create a file header, with a defaulted georeferencing. */
359 : /* -------------------------------------------------------------------- */
360 : char achHeader[160];
361 : int nValue32;
362 : double dfValue;
363 :
364 3 : memset( achHeader, 0, sizeof(achHeader));
365 :
366 3 : memcpy( achHeader+0, "CTABLE V2.0 ", 16 );
367 :
368 3 : if( CSLFetchNameValue( papszOptions, "DESCRIPTION" ) != NULL )
369 : strncpy( achHeader + 16,
370 : CSLFetchNameValue( papszOptions, "DESCRIPTION" ),
371 0 : 80 );
372 :
373 : // lower left origin (longitude, center of pixel, radians)
374 3 : dfValue = 0;
375 : CPL_LSBPTR64( &dfValue );
376 3 : memcpy( achHeader + 96, &dfValue, 8 );
377 :
378 : // lower left origin (latitude, center of pixel, radians)
379 3 : dfValue = 0;
380 : CPL_LSBPTR64( &dfValue );
381 3 : memcpy( achHeader + 104, &dfValue, 8 );
382 :
383 : // pixel width (radians)
384 3 : dfValue = 0.01 * M_PI / 180.0;
385 : CPL_LSBPTR64( &dfValue );
386 3 : memcpy( achHeader + 112, &dfValue, 8 );
387 :
388 : // pixel height (radians)
389 3 : dfValue = 0.01 * M_PI / 180.0;
390 : CPL_LSBPTR64( &dfValue );
391 3 : memcpy( achHeader + 120, &dfValue, 8 );
392 :
393 : // raster width in pixels
394 3 : nValue32 = nXSize;
395 : CPL_LSBPTR32( &nValue32 );
396 3 : memcpy( achHeader + 128, &nValue32, 4 );
397 :
398 : // raster width in pixels
399 3 : nValue32 = nYSize;
400 : CPL_LSBPTR32( &nValue32 );
401 3 : memcpy( achHeader + 132, &nValue32, 4 );
402 :
403 3 : VSIFWriteL( achHeader, 1, sizeof(achHeader), fp );
404 :
405 : /* -------------------------------------------------------------------- */
406 : /* Write zeroed grid data. */
407 : /* -------------------------------------------------------------------- */
408 3 : float *pafLine = (float *) CPLCalloc(sizeof(float)*2,nXSize);
409 : int i;
410 :
411 213 : for( i = 0; i < nYSize; i++ )
412 : {
413 210 : if( (int)VSIFWriteL( pafLine, sizeof(float)*2, nXSize, fp ) != nXSize )
414 : {
415 : CPLError( CE_Failure, CPLE_FileIO,
416 : "Write failed at line %d, perhaps the disk is full?",
417 0 : i );
418 0 : return NULL;
419 : }
420 : }
421 :
422 : /* -------------------------------------------------------------------- */
423 : /* Cleanup and return. */
424 : /* -------------------------------------------------------------------- */
425 3 : CPLFree( pafLine );
426 :
427 3 : VSIFCloseL( fp );
428 :
429 3 : return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
430 : }
431 :
432 : /************************************************************************/
433 : /* GDALRegister_CTable2() */
434 : /************************************************************************/
435 :
436 582 : void GDALRegister_CTable2()
437 :
438 : {
439 : GDALDriver *poDriver;
440 :
441 582 : if( GDALGetDriverByName( "CTable2" ) == NULL )
442 : {
443 561 : poDriver = new GDALDriver();
444 :
445 561 : poDriver->SetDescription( "CTable2" );
446 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
447 561 : "CTable2 Datum Grid Shift" );
448 561 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
449 :
450 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
451 561 : "Float32" );
452 :
453 561 : poDriver->pfnOpen = CTable2Dataset::Open;
454 561 : poDriver->pfnIdentify = CTable2Dataset::Identify;
455 561 : poDriver->pfnCreate = CTable2Dataset::Create;
456 :
457 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
458 : }
459 582 : }
460 :
|
