1 : /******************************************************************************
2 : * $Id: gtxdataset.cpp 24534 2012-06-03 16:10:29Z rouault $
3 : *
4 : * Project: Vertical Datum Transformation
5 : * Purpose: Implementation of NOAA .gtx vertical datum shift file format.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2010, Frank Warmerdam
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "rawdataset.h"
31 : #include "cpl_string.h"
32 : #include "ogr_srs_api.h"
33 :
34 : CPL_CVSID("$Id: gtxdataset.cpp 24534 2012-06-03 16:10:29Z rouault $");
35 :
36 : /**
37 :
38 : NOAA .GTX Vertical Datum Grid Shift Format
39 :
40 : All values are bigendian
41 :
42 : Header
43 : ------
44 :
45 : float64 latitude_of_origin
46 : float64 longitude_of_origin (0-360)
47 : float64 cell size (x?y?)
48 : float64 cell size (x?y?)
49 : int32 length in pixels
50 : int32 width in pixels
51 :
52 : Data
53 : ----
54 :
55 : float32 * width in pixels * length in pixels
56 :
57 : Values are an offset in meters between two vertical datums.
58 :
59 : **/
60 :
61 : /************************************************************************/
62 : /* ==================================================================== */
63 : /* GTXDataset */
64 : /* ==================================================================== */
65 : /************************************************************************/
66 :
67 : class GTXDataset : public RawDataset
68 : {
69 : public:
70 : VSILFILE *fpImage; // image data file.
71 :
72 : double adfGeoTransform[6];
73 :
74 : public:
75 : GTXDataset();
76 : ~GTXDataset();
77 :
78 : virtual CPLErr GetGeoTransform( double * padfTransform );
79 : virtual CPLErr SetGeoTransform( double * padfTransform );
80 : virtual const char *GetProjectionRef();
81 :
82 : static GDALDataset *Open( GDALOpenInfo * );
83 : static int Identify( GDALOpenInfo * );
84 : static GDALDataset *Create( const char * pszFilename,
85 : int nXSize, int nYSize, int nBands,
86 : GDALDataType eType, char ** papszOptions );
87 : };
88 :
89 : /************************************************************************/
90 : /* ==================================================================== */
91 : /* GTXDataset */
92 : /* ==================================================================== */
93 : /************************************************************************/
94 :
95 : /************************************************************************/
96 : /* GTXDataset() */
97 : /************************************************************************/
98 :
99 6 : GTXDataset::GTXDataset()
100 : {
101 6 : fpImage = NULL;
102 6 : }
103 :
104 : /************************************************************************/
105 : /* ~GTXDataset() */
106 : /************************************************************************/
107 :
108 6 : GTXDataset::~GTXDataset()
109 :
110 : {
111 6 : FlushCache();
112 :
113 6 : if( fpImage != NULL )
114 6 : VSIFCloseL( fpImage );
115 6 : }
116 :
117 : /************************************************************************/
118 : /* Identify() */
119 : /************************************************************************/
120 :
121 11951 : int GTXDataset::Identify( GDALOpenInfo *poOpenInfo )
122 :
123 : {
124 11951 : if( poOpenInfo->nHeaderBytes < 40 )
125 11357 : return FALSE;
126 :
127 594 : if( !EQUAL(CPLGetExtension(poOpenInfo->pszFilename),"gtx") )
128 588 : return FALSE;
129 :
130 6 : return TRUE;
131 : }
132 :
133 : /************************************************************************/
134 : /* Open() */
135 : /************************************************************************/
136 :
137 1858 : GDALDataset *GTXDataset::Open( GDALOpenInfo * poOpenInfo )
138 :
139 : {
140 1858 : if( !Identify( poOpenInfo ) )
141 1852 : return NULL;
142 :
143 : /* -------------------------------------------------------------------- */
144 : /* Create a corresponding GDALDataset. */
145 : /* -------------------------------------------------------------------- */
146 : GTXDataset *poDS;
147 :
148 6 : poDS = new GTXDataset();
149 :
150 6 : poDS->eAccess = poOpenInfo->eAccess;
151 :
152 : /* -------------------------------------------------------------------- */
153 : /* Open the file. */
154 : /* -------------------------------------------------------------------- */
155 6 : if( poOpenInfo->eAccess == GA_ReadOnly )
156 3 : poDS->fpImage = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
157 : else
158 3 : poDS->fpImage = VSIFOpenL( poOpenInfo->pszFilename, "rb+" );
159 :
160 6 : if( poDS->fpImage == NULL )
161 : {
162 0 : delete poDS;
163 0 : return NULL;
164 : }
165 :
166 : /* -------------------------------------------------------------------- */
167 : /* Read the header. */
168 : /* -------------------------------------------------------------------- */
169 6 : poDS->adfGeoTransform[2] = 0.0;
170 6 : poDS->adfGeoTransform[4] = 0.0;
171 :
172 6 : VSIFReadL( poDS->adfGeoTransform+3, 8, 1, poDS->fpImage );
173 6 : VSIFReadL( poDS->adfGeoTransform+0, 8, 1, poDS->fpImage );
174 6 : VSIFReadL( poDS->adfGeoTransform+5, 8, 1, poDS->fpImage );
175 6 : VSIFReadL( poDS->adfGeoTransform+1, 8, 1, poDS->fpImage );
176 :
177 6 : VSIFReadL( &(poDS->nRasterYSize), 4, 1, poDS->fpImage );
178 6 : VSIFReadL( &(poDS->nRasterXSize), 4, 1, poDS->fpImage );
179 :
180 6 : CPL_MSBPTR32( &(poDS->nRasterYSize) );
181 6 : CPL_MSBPTR32( &(poDS->nRasterXSize) );
182 :
183 6 : CPL_MSBPTR64( poDS->adfGeoTransform + 0 );
184 6 : CPL_MSBPTR64( poDS->adfGeoTransform + 1 );
185 6 : CPL_MSBPTR64( poDS->adfGeoTransform + 3 );
186 6 : CPL_MSBPTR64( poDS->adfGeoTransform + 5 );
187 :
188 : poDS->adfGeoTransform[3] +=
189 6 : poDS->adfGeoTransform[5] * (poDS->nRasterYSize-1);
190 :
191 6 : poDS->adfGeoTransform[0] -= poDS->adfGeoTransform[1] * 0.5;
192 6 : poDS->adfGeoTransform[3] += poDS->adfGeoTransform[5] * 0.5;
193 :
194 6 : poDS->adfGeoTransform[5] *= -1;
195 :
196 6 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
197 : {
198 0 : delete poDS;
199 0 : return NULL;
200 : }
201 :
202 : /* -------------------------------------------------------------------- */
203 : /* Guess the data type. Since October 1, 2009, it should be */
204 : /* Float32. Before it was double. */
205 : /* -------------------------------------------------------------------- */
206 6 : GDALDataType eDT = GDT_Float32;
207 6 : VSIFSeekL(poDS->fpImage, 0, SEEK_END);
208 6 : vsi_l_offset nSize = VSIFTellL(poDS->fpImage);
209 6 : if( nSize == 40 + 8 * poDS->nRasterXSize * poDS->nRasterYSize )
210 0 : eDT = GDT_Float64;
211 6 : int nDTSize = GDALGetDataTypeSize(eDT) / 8;
212 :
213 : /* -------------------------------------------------------------------- */
214 : /* Create band information object. */
215 : /* -------------------------------------------------------------------- */
216 : poDS->SetBand(
217 : 1, new RawRasterBand( poDS, 1, poDS->fpImage,
218 : (poDS->nRasterYSize-1)*poDS->nRasterXSize*nDTSize + 40,
219 : nDTSize, poDS->nRasterXSize * -nDTSize,
220 : eDT,
221 6 : !CPL_IS_LSB, TRUE, FALSE ) );
222 :
223 : /* -------------------------------------------------------------------- */
224 : /* Initialize any PAM information. */
225 : /* -------------------------------------------------------------------- */
226 6 : poDS->SetDescription( poOpenInfo->pszFilename );
227 6 : poDS->TryLoadXML();
228 :
229 : /* -------------------------------------------------------------------- */
230 : /* Check for overviews. */
231 : /* -------------------------------------------------------------------- */
232 6 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
233 :
234 6 : return( poDS );
235 : }
236 :
237 : /************************************************************************/
238 : /* GetGeoTransform() */
239 : /************************************************************************/
240 :
241 3 : CPLErr GTXDataset::GetGeoTransform( double * padfTransform )
242 :
243 : {
244 3 : memcpy( padfTransform, adfGeoTransform, sizeof(double)*6 );
245 3 : return CE_None;
246 : }
247 :
248 : /************************************************************************/
249 : /* SetGeoTransform() */
250 : /************************************************************************/
251 :
252 3 : CPLErr GTXDataset::SetGeoTransform( double * padfTransform )
253 :
254 : {
255 3 : if( padfTransform[2] != 0.0 || padfTransform[4] != 0.0 )
256 : {
257 : CPLError( CE_Failure, CPLE_AppDefined,
258 0 : "Attempt to write skewed or rotated geotransform to gtx." );
259 0 : return CE_Failure;
260 : }
261 :
262 3 : memcpy( adfGeoTransform, padfTransform, sizeof(double)*6 );
263 :
264 : unsigned char header[32];
265 : double dfXOrigin, dfYOrigin, dfWidth, dfHeight;
266 :
267 3 : dfXOrigin = adfGeoTransform[0] + 0.5 * adfGeoTransform[1];
268 3 : dfYOrigin = adfGeoTransform[3] + (nRasterYSize-0.5) * adfGeoTransform[5];
269 3 : dfWidth = adfGeoTransform[1];
270 3 : dfHeight = - adfGeoTransform[5];
271 :
272 :
273 3 : memcpy( header + 0, &dfYOrigin, 8 );
274 3 : CPL_MSBPTR64( header + 0 );
275 :
276 3 : memcpy( header + 8, &dfXOrigin, 8 );
277 3 : CPL_MSBPTR64( header + 8 );
278 :
279 3 : memcpy( header + 16, &dfHeight, 8 );
280 3 : CPL_MSBPTR64( header + 16 );
281 :
282 3 : memcpy( header + 24, &dfWidth, 8 );
283 3 : CPL_MSBPTR64( header + 24 );
284 :
285 3 : if( VSIFSeekL( fpImage, SEEK_SET, 0 ) != 0
286 : || VSIFWriteL( header, 32, 1, fpImage ) != 1 )
287 : {
288 : CPLError( CE_Failure, CPLE_AppDefined,
289 0 : "Attempt to write geotrasform header to gtx failed." );
290 0 : return CE_Failure;
291 : }
292 :
293 3 : return CE_None;
294 : }
295 :
296 : /************************************************************************/
297 : /* GetProjectionRef() */
298 : /************************************************************************/
299 :
300 1 : const char *GTXDataset::GetProjectionRef()
301 :
302 : {
303 1 : return SRS_WKT_WGS84;
304 : }
305 :
306 : /************************************************************************/
307 : /* Create() */
308 : /************************************************************************/
309 :
310 44 : GDALDataset *GTXDataset::Create( const char * pszFilename,
311 : int nXSize, int nYSize, int nBands,
312 : GDALDataType eType,
313 : char ** papszOptions )
314 :
315 : {
316 44 : if( eType != GDT_Float32 )
317 : {
318 : CPLError( CE_Failure, CPLE_AppDefined,
319 : "Attempt to create gtx file with unsupported data type '%s'.",
320 38 : GDALGetDataTypeName( eType ) );
321 38 : return NULL;
322 : }
323 :
324 6 : if( !EQUAL(CPLGetExtension(pszFilename),"gtx") )
325 : {
326 : CPLError( CE_Failure, CPLE_AppDefined,
327 0 : "Attempt to create gtx file with extension other than gtx." );
328 0 : return NULL;
329 : }
330 :
331 : /* -------------------------------------------------------------------- */
332 : /* Try to create the file. */
333 : /* -------------------------------------------------------------------- */
334 : VSILFILE *fp;
335 :
336 6 : fp = VSIFOpenL( pszFilename, "wb" );
337 :
338 6 : if( fp == NULL )
339 : {
340 : CPLError( CE_Failure, CPLE_OpenFailed,
341 : "Attempt to create file `%s' failed.\n",
342 3 : pszFilename );
343 3 : return NULL;
344 : }
345 :
346 : /* -------------------------------------------------------------------- */
347 : /* Write out the header with stub georeferencing. */
348 : /* -------------------------------------------------------------------- */
349 : unsigned char header[40];
350 3 : double dfXOrigin=0, dfYOrigin=0, dfXSize=0.01, dfYSize=0.01;
351 3 : GInt32 nXSize32 = nXSize, nYSize32 = nYSize;
352 :
353 3 : memcpy( header + 0, &dfYOrigin, 8 );
354 3 : CPL_MSBPTR64( header + 0 );
355 :
356 3 : memcpy( header + 8, &dfXOrigin, 8 );
357 3 : CPL_MSBPTR64( header + 8 );
358 :
359 3 : memcpy( header + 16, &dfYSize, 8 );
360 3 : CPL_MSBPTR64( header + 16 );
361 :
362 3 : memcpy( header + 24, &dfXSize, 8 );
363 3 : CPL_MSBPTR64( header + 24 );
364 :
365 3 : memcpy( header + 32, &nYSize32, 4 );
366 3 : CPL_MSBPTR32( header + 32 );
367 3 : memcpy( header + 36, &nXSize32, 4 );
368 3 : CPL_MSBPTR32( header + 36 );
369 :
370 3 : VSIFWriteL( header, 40, 1, fp );
371 3 : VSIFCloseL( fp );
372 :
373 3 : return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
374 : }
375 :
376 :
377 : /************************************************************************/
378 : /* GDALRegister_GTX() */
379 : /************************************************************************/
380 :
381 582 : void GDALRegister_GTX()
382 :
383 : {
384 : GDALDriver *poDriver;
385 :
386 582 : if( GDALGetDriverByName( "GTX" ) == NULL )
387 : {
388 561 : poDriver = new GDALDriver();
389 :
390 561 : poDriver->SetDescription( "GTX" );
391 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
392 561 : "NOAA Vertical Datum .GTX" );
393 561 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "gtx" );
394 561 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
395 : // poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
396 : // "frmt_various.html#GTX" );
397 :
398 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
399 561 : "Float32" );
400 :
401 561 : poDriver->pfnOpen = GTXDataset::Open;
402 561 : poDriver->pfnIdentify = GTXDataset::Identify;
403 561 : poDriver->pfnCreate = GTXDataset::Create;
404 :
405 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
406 : }
407 582 : }
408 :
|
