1 : /******************************************************************************
2 : * $Id: landataset.cpp 17117 2009-05-25 19:26:01Z warmerdam $
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: landataset.cpp 17117 2009-05-25 19:26:01Z warmerdam $");
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 : FILE *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 4 : GTXDataset::GTXDataset()
100 : {
101 4 : fpImage = NULL;
102 4 : }
103 :
104 : /************************************************************************/
105 : /* ~GTXDataset() */
106 : /************************************************************************/
107 :
108 4 : GTXDataset::~GTXDataset()
109 :
110 : {
111 4 : FlushCache();
112 :
113 4 : if( fpImage != NULL )
114 4 : VSIFCloseL( fpImage );
115 4 : }
116 :
117 : /************************************************************************/
118 : /* Identify() */
119 : /************************************************************************/
120 :
121 10010 : int GTXDataset::Identify( GDALOpenInfo *poOpenInfo )
122 :
123 : {
124 10010 : if( poOpenInfo->nHeaderBytes < 40 )
125 9513 : return FALSE;
126 :
127 497 : if( !EQUAL(CPLGetExtension(poOpenInfo->pszFilename),"gtx") )
128 493 : return FALSE;
129 :
130 4 : return TRUE;
131 : }
132 :
133 : /************************************************************************/
134 : /* Open() */
135 : /************************************************************************/
136 :
137 1765 : GDALDataset *GTXDataset::Open( GDALOpenInfo * poOpenInfo )
138 :
139 : {
140 1765 : if( !Identify( poOpenInfo ) )
141 1761 : return NULL;
142 :
143 : /* -------------------------------------------------------------------- */
144 : /* Create a corresponding GDALDataset. */
145 : /* -------------------------------------------------------------------- */
146 : GTXDataset *poDS;
147 :
148 4 : poDS = new GTXDataset();
149 :
150 4 : poDS->eAccess = poOpenInfo->eAccess;
151 :
152 : /* -------------------------------------------------------------------- */
153 : /* Open the file. */
154 : /* -------------------------------------------------------------------- */
155 4 : if( poOpenInfo->eAccess == GA_ReadOnly )
156 2 : poDS->fpImage = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
157 : else
158 2 : poDS->fpImage = VSIFOpenL( poOpenInfo->pszFilename, "rb+" );
159 :
160 4 : if( poDS->fpImage == NULL )
161 : {
162 0 : delete poDS;
163 0 : return NULL;
164 : }
165 :
166 : /* -------------------------------------------------------------------- */
167 : /* Read the header. */
168 : /* -------------------------------------------------------------------- */
169 4 : poDS->adfGeoTransform[2] = 0.0;
170 4 : poDS->adfGeoTransform[4] = 0.0;
171 :
172 4 : VSIFReadL( poDS->adfGeoTransform+3, 8, 1, poDS->fpImage );
173 4 : VSIFReadL( poDS->adfGeoTransform+0, 8, 1, poDS->fpImage );
174 4 : VSIFReadL( poDS->adfGeoTransform+5, 8, 1, poDS->fpImage );
175 4 : VSIFReadL( poDS->adfGeoTransform+1, 8, 1, poDS->fpImage );
176 :
177 4 : VSIFReadL( &(poDS->nRasterYSize), 4, 1, poDS->fpImage );
178 4 : VSIFReadL( &(poDS->nRasterXSize), 4, 1, poDS->fpImage );
179 :
180 4 : CPL_MSBPTR32( &(poDS->nRasterYSize) );
181 4 : CPL_MSBPTR32( &(poDS->nRasterXSize) );
182 :
183 4 : CPL_MSBPTR64( poDS->adfGeoTransform + 0 );
184 4 : CPL_MSBPTR64( poDS->adfGeoTransform + 1 );
185 4 : CPL_MSBPTR64( poDS->adfGeoTransform + 3 );
186 4 : CPL_MSBPTR64( poDS->adfGeoTransform + 5 );
187 :
188 : poDS->adfGeoTransform[3] +=
189 4 : poDS->adfGeoTransform[5] * (poDS->nRasterYSize-1);
190 :
191 4 : poDS->adfGeoTransform[0] -= poDS->adfGeoTransform[1] * 0.5;
192 4 : poDS->adfGeoTransform[3] += poDS->adfGeoTransform[5] * 0.5;
193 :
194 4 : poDS->adfGeoTransform[5] *= -1;
195 :
196 4 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
197 : {
198 0 : delete poDS;
199 0 : return NULL;
200 : }
201 :
202 : /* -------------------------------------------------------------------- */
203 : /* Create band information object. */
204 : /* -------------------------------------------------------------------- */
205 : poDS->SetBand(
206 : 1, new RawRasterBand( poDS, 1, poDS->fpImage,
207 : (poDS->nRasterYSize-1)*poDS->nRasterXSize*4 + 40,
208 : 4, poDS->nRasterXSize * -4,
209 : GDT_Float32,
210 4 : !CPL_IS_LSB, TRUE, FALSE ) );
211 :
212 : /* -------------------------------------------------------------------- */
213 : /* Initialize any PAM information. */
214 : /* -------------------------------------------------------------------- */
215 4 : poDS->SetDescription( poOpenInfo->pszFilename );
216 4 : poDS->TryLoadXML();
217 :
218 : /* -------------------------------------------------------------------- */
219 : /* Check for overviews. */
220 : /* -------------------------------------------------------------------- */
221 4 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
222 :
223 4 : return( poDS );
224 : }
225 :
226 : /************************************************************************/
227 : /* GetGeoTransform() */
228 : /************************************************************************/
229 :
230 0 : CPLErr GTXDataset::GetGeoTransform( double * padfTransform )
231 :
232 : {
233 0 : memcpy( padfTransform, adfGeoTransform, sizeof(double)*6 );
234 0 : return CE_None;
235 : }
236 :
237 : /************************************************************************/
238 : /* SetGeoTransform() */
239 : /************************************************************************/
240 :
241 2 : CPLErr GTXDataset::SetGeoTransform( double * padfTransform )
242 :
243 : {
244 2 : if( padfTransform[2] != 0.0 || padfTransform[4] != 0.0 )
245 : {
246 : CPLError( CE_Failure, CPLE_AppDefined,
247 0 : "Attempt to write skewed or rotated geotransform to gtx." );
248 0 : return CE_Failure;
249 : }
250 :
251 2 : memcpy( adfGeoTransform, padfTransform, sizeof(double)*6 );
252 :
253 : unsigned char header[32];
254 : double dfXOrigin, dfYOrigin, dfWidth, dfHeight;
255 :
256 2 : dfXOrigin = adfGeoTransform[0] + 0.5 * adfGeoTransform[1];
257 2 : dfYOrigin = adfGeoTransform[3] + (nRasterYSize-0.5) * adfGeoTransform[5];
258 2 : dfWidth = adfGeoTransform[1];
259 2 : dfHeight = - adfGeoTransform[5];
260 :
261 :
262 2 : memcpy( header + 0, &dfYOrigin, 8 );
263 2 : CPL_MSBPTR64( header + 0 );
264 :
265 2 : memcpy( header + 8, &dfXOrigin, 8 );
266 2 : CPL_MSBPTR64( header + 8 );
267 :
268 2 : memcpy( header + 16, &dfHeight, 8 );
269 2 : CPL_MSBPTR64( header + 16 );
270 :
271 2 : memcpy( header + 24, &dfWidth, 8 );
272 2 : CPL_MSBPTR64( header + 24 );
273 :
274 2 : if( VSIFSeekL( fpImage, SEEK_SET, 0 ) != 0
275 : || VSIFWriteL( header, 32, 1, fpImage ) != 1 )
276 : {
277 : CPLError( CE_Failure, CPLE_AppDefined,
278 0 : "Attempt to write geotrasform header to gtx failed." );
279 0 : return CE_Failure;
280 : }
281 :
282 2 : return CE_None;
283 : }
284 :
285 : /************************************************************************/
286 : /* GetProjectionRef() */
287 : /************************************************************************/
288 :
289 0 : const char *GTXDataset::GetProjectionRef()
290 :
291 : {
292 0 : return SRS_WKT_WGS84;
293 : }
294 :
295 : /************************************************************************/
296 : /* Create() */
297 : /************************************************************************/
298 :
299 : GDALDataset *GTXDataset::Create( const char * pszFilename,
300 : int nXSize, int nYSize, int nBands,
301 : GDALDataType eType,
302 41 : char ** papszOptions )
303 :
304 : {
305 41 : if( eType != GDT_Float32 )
306 : {
307 : CPLError( CE_Failure, CPLE_AppDefined,
308 : "Attempt to create gtx file with unsupported data type '%s'.",
309 38 : GDALGetDataTypeName( eType ) );
310 38 : return NULL;
311 : }
312 :
313 3 : if( !EQUAL(CPLGetExtension(pszFilename),"gtx") )
314 : {
315 : CPLError( CE_Failure, CPLE_AppDefined,
316 1 : "Attempt to create gtx file with extension other than gtx." );
317 1 : return NULL;
318 : }
319 :
320 : /* -------------------------------------------------------------------- */
321 : /* Try to create the file. */
322 : /* -------------------------------------------------------------------- */
323 : FILE *fp;
324 :
325 2 : fp = VSIFOpenL( pszFilename, "wb" );
326 :
327 2 : if( fp == NULL )
328 : {
329 : CPLError( CE_Failure, CPLE_OpenFailed,
330 : "Attempt to create file `%s' failed.\n",
331 0 : pszFilename );
332 0 : return NULL;
333 : }
334 :
335 : /* -------------------------------------------------------------------- */
336 : /* Write out the header with stub georeferencing. */
337 : /* -------------------------------------------------------------------- */
338 : unsigned char header[40];
339 2 : double dfXOrigin=0, dfYOrigin=0, dfXSize=0.01, dfYSize=0.01;
340 2 : GInt32 nXSize32 = nXSize, nYSize32 = nYSize;
341 :
342 2 : memcpy( header + 0, &dfYOrigin, 8 );
343 2 : CPL_MSBPTR64( header + 0 );
344 :
345 2 : memcpy( header + 8, &dfXOrigin, 8 );
346 2 : CPL_MSBPTR64( header + 8 );
347 :
348 2 : memcpy( header + 16, &dfYSize, 8 );
349 2 : CPL_MSBPTR64( header + 16 );
350 :
351 2 : memcpy( header + 24, &dfXSize, 8 );
352 2 : CPL_MSBPTR64( header + 24 );
353 :
354 2 : memcpy( header + 32, &nYSize32, 4 );
355 2 : CPL_MSBPTR32( header + 32 );
356 2 : memcpy( header + 36, &nXSize32, 4 );
357 2 : CPL_MSBPTR32( header + 36 );
358 :
359 2 : VSIFWriteL( header, 40, 1, fp );
360 2 : VSIFCloseL( fp );
361 :
362 2 : return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
363 : }
364 :
365 :
366 : /************************************************************************/
367 : /* GDALRegister_GTX() */
368 : /************************************************************************/
369 :
370 409 : void GDALRegister_GTX()
371 :
372 : {
373 : GDALDriver *poDriver;
374 :
375 409 : if( GDALGetDriverByName( "GTX" ) == NULL )
376 : {
377 392 : poDriver = new GDALDriver();
378 :
379 392 : poDriver->SetDescription( "GTX" );
380 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
381 392 : "NOAA Vertical Datum .GTX" );
382 392 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "gtx" );
383 392 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
384 : // poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
385 : // "frmt_various.html#GTX" );
386 :
387 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
388 392 : "Float32" );
389 :
390 392 : poDriver->pfnOpen = GTXDataset::Open;
391 392 : poDriver->pfnIdentify = GTXDataset::Identify;
392 392 : poDriver->pfnCreate = GTXDataset::Create;
393 :
394 392 : GetGDALDriverManager()->RegisterDriver( poDriver );
395 : }
396 409 : }
397 :
|
