1 : /****************************************************************************
2 : * $Id: gs7bgdataset.cpp 17664 2009-09-21 21:16:45Z rouault $
3 : *
4 : * Project: GDAL
5 : * Purpose: Implements the Golden Software Surfer 7 Binary Grid Format.
6 : * Author: Adam Guernsey, adam@ctech.com
7 : * (Based almost entirely on gsbgdataset.cpp by Kevin Locke)
8 : *
9 : ****************************************************************************
10 : * Copyright (c) 2007, Adam Guernsey <adam@ctech.com>
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 <float.h>
32 : #include <limits.h>
33 : #include <math.h>
34 : #include <assert.h>
35 :
36 : #include "gdal_pam.h"
37 :
38 : #ifndef DBL_MAX
39 : # ifdef __DBL_MAX__
40 : # define DBL_MAX __DBL_MAX__
41 : # else
42 : # define DBL_MAX 1.7976931348623157E+308
43 : # endif /* __DBL_MAX__ */
44 : #endif /* DBL_MAX */
45 :
46 : #ifndef FLT_MAX
47 : # ifdef __FLT_MAX__
48 : # define FLT_MAX __FLT_MAX__
49 : # else
50 : # define FLT_MAX 3.40282347E+38F
51 : # endif /* __FLT_MAX__ */
52 : #endif /* FLT_MAX */
53 :
54 : #ifndef INT_MAX
55 : # define INT_MAX 2147483647
56 : #endif /* INT_MAX */
57 :
58 : #ifndef SHRT_MAX
59 : # define SHRT_MAX 32767
60 : #endif /* SHRT_MAX */
61 :
62 : CPL_CVSID("$Id: gs7bgdataset.cpp 17664 2009-09-21 21:16:45Z rouault $");
63 :
64 : CPL_C_START
65 : void GDALRegister_GS7BG(void);
66 : CPL_C_END
67 :
68 : /************************************************************************/
69 : /* ==================================================================== */
70 : /* GS7BGDataset */
71 : /* ==================================================================== */
72 : /************************************************************************/
73 :
74 : class GS7BGRasterBand;
75 :
76 : class GS7BGDataset : public GDALPamDataset
77 1 : {
78 : friend class GS7BGRasterBand;
79 :
80 : static double dfNoData_Value;
81 : static size_t nData_Position;
82 :
83 : FILE *fp;
84 :
85 : public:
86 : ~GS7BGDataset();
87 :
88 : static GDALDataset *Open( GDALOpenInfo * );
89 :
90 : CPLErr GetGeoTransform( double *padfGeoTransform );
91 : };
92 :
93 : /* NOTE: This is not mentioned in the spec, but Surfer 8 uses this value */
94 : /* 0x7effffee (Little Endian: eeffff7e) */
95 : double GS7BGDataset::dfNoData_Value = 1.701410009187828e+38f;
96 :
97 : size_t GS7BGDataset::nData_Position = 0;
98 :
99 : const long nHEADER_TAG = 0x42525344;
100 : const long nGRID_TAG = 0x44495247;
101 : const long nDATA_TAG = 0x41544144;
102 : const long nFAULT_TAG = 0x49544c46;
103 :
104 : /************************************************************************/
105 : /* ==================================================================== */
106 : /* GS7BGRasterBand */
107 : /* ==================================================================== */
108 : /************************************************************************/
109 :
110 : class GS7BGRasterBand : public GDALPamRasterBand
111 2 : {
112 : friend class GS7BGDataset;
113 :
114 : double dfMinX;
115 : double dfMaxX;
116 : double dfMinY;
117 : double dfMaxY;
118 : double dfMinZ;
119 : double dfMaxZ;
120 :
121 : public:
122 :
123 : GS7BGRasterBand( GS7BGDataset *, int );
124 :
125 : CPLErr IReadBlock( int, int, void * );
126 :
127 : double GetNoDataValue( int *pbSuccess = NULL );
128 : };
129 :
130 : /************************************************************************/
131 : /* GS7BGRasterBand() */
132 : /************************************************************************/
133 :
134 1 : GS7BGRasterBand::GS7BGRasterBand( GS7BGDataset *poDS, int nBand )
135 : {
136 1 : this->poDS = poDS;
137 1 : nBand = nBand;
138 :
139 1 : eDataType = GDT_Float64;
140 :
141 1 : nBlockXSize = poDS->GetRasterXSize();
142 1 : nBlockYSize = 1;
143 1 : }
144 :
145 : /************************************************************************/
146 : /* IReadBlock() */
147 : /************************************************************************/
148 :
149 20 : CPLErr GS7BGRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
150 : void * pImage )
151 :
152 : {
153 20 : if( nBlockYOff < 0 || nBlockYOff > nRasterYSize - 1 || nBlockXOff != 0 )
154 0 : return CE_Failure;
155 :
156 20 : GS7BGDataset *poGDS = (GS7BGDataset *) ( poDS );
157 :
158 20 : if( VSIFSeekL( poGDS->fp,
159 : ( GS7BGDataset::nData_Position +
160 : sizeof(double) * nRasterXSize * (nRasterYSize - nBlockYOff - 1) ),
161 : SEEK_SET ) != 0 )
162 : {
163 : CPLError( CE_Failure, CPLE_FileIO,
164 0 : "Unable to seek to beginning of grid row.\n" );
165 0 : return CE_Failure;
166 : }
167 :
168 20 : if( VSIFReadL( pImage, sizeof(double), nBlockXSize,
169 : poGDS->fp ) != static_cast<unsigned>(nBlockXSize) )
170 : {
171 : CPLError( CE_Failure, CPLE_FileIO,
172 0 : "Unable to read block from grid file.\n" );
173 0 : return CE_Failure;
174 : }
175 :
176 : double *pfImage;
177 20 : pfImage = (double *)pImage;
178 20 : for( int iPixel=0; iPixel<nBlockXSize; iPixel++ )
179 : CPL_LSBPTR64( pfImage + iPixel );
180 :
181 20 : return CE_None;
182 : }
183 :
184 : /************************************************************************/
185 : /* GetNoDataValue() */
186 : /************************************************************************/
187 :
188 0 : double GS7BGRasterBand::GetNoDataValue( int * pbSuccess )
189 : {
190 0 : if( pbSuccess )
191 0 : *pbSuccess = TRUE;
192 :
193 0 : return GS7BGDataset::dfNoData_Value;
194 : }
195 :
196 : /************************************************************************/
197 : /* ==================================================================== */
198 : /* GS7BGDataset */
199 : /* ==================================================================== */
200 : /************************************************************************/
201 :
202 2 : GS7BGDataset::~GS7BGDataset()
203 :
204 : {
205 1 : FlushCache();
206 1 : if( fp != NULL )
207 1 : VSIFCloseL( fp );
208 2 : }
209 :
210 : /************************************************************************/
211 : /* Open() */
212 : /************************************************************************/
213 :
214 8785 : GDALDataset *GS7BGDataset::Open( GDALOpenInfo * poOpenInfo )
215 :
216 : {
217 : /* Check for signature */
218 8785 : if( poOpenInfo->nHeaderBytes < 4
219 : || !EQUALN((const char *) poOpenInfo->pabyHeader,"DSRB",4) )
220 : {
221 8784 : return NULL;
222 : }
223 :
224 : /* -------------------------------------------------------------------- */
225 : /* Confirm the requested access is supported. */
226 : /* -------------------------------------------------------------------- */
227 1 : if( poOpenInfo->eAccess == GA_Update )
228 : {
229 : CPLError( CE_Failure, CPLE_NotSupported,
230 : "The GS7BG driver does not support update access to existing"
231 0 : " datasets.\n" );
232 0 : return NULL;
233 : }
234 : /* ------------------------------------------------------------------- */
235 : /* Create a corresponding GDALDataset. */
236 : /* ------------------------------------------------------------------- */
237 1 : GS7BGDataset *poDS = new GS7BGDataset();
238 :
239 : /* ------------------------------------------------------------------- */
240 : /* Open file with large file API. */
241 : /* ------------------------------------------------------------------- */
242 1 : poDS->eAccess = poOpenInfo->eAccess;
243 1 : if( poOpenInfo->eAccess == GA_ReadOnly )
244 1 : poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
245 : else
246 0 : poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "r+b" );
247 :
248 1 : if( poDS->fp == NULL )
249 : {
250 0 : delete poDS;
251 : CPLError( CE_Failure, CPLE_OpenFailed,
252 : "VSIFOpenL(%s) failed unexpectedly.",
253 0 : poOpenInfo->pszFilename );
254 0 : return NULL;
255 : }
256 :
257 : /* ------------------------------------------------------------------- */
258 : /* Read the header. The Header section must be the first section */
259 : /* in the file. */
260 : /* ------------------------------------------------------------------- */
261 1 : if( VSIFSeekL( poDS->fp, 0, SEEK_SET ) != 0 )
262 : {
263 0 : delete poDS;
264 : CPLError( CE_Failure, CPLE_FileIO,
265 0 : "Unable to seek to start of grid file header.\n" );
266 0 : return NULL;
267 : }
268 :
269 : GInt32 nTag;
270 : GInt32 nSize;
271 : GInt32 nVersion;
272 :
273 1 : if( VSIFReadL( (void *)&nTag, sizeof(GInt32), 1, poDS->fp ) != 1 )
274 : {
275 0 : delete poDS;
276 0 : CPLError( CE_Failure, CPLE_FileIO, "Unable to read Tag.\n" );
277 0 : return NULL;
278 : }
279 :
280 : CPL_LSBPTR32( &nTag );
281 :
282 1 : if(nTag != nHEADER_TAG)
283 : {
284 0 : delete poDS;
285 0 : CPLError( CE_Failure, CPLE_FileIO, "Header tag not found.\n" );
286 0 : return NULL;
287 : }
288 :
289 1 : if( VSIFReadL( (void *)&nSize, sizeof(GInt32), 1, poDS->fp ) != 1 )
290 : {
291 0 : delete poDS;
292 : CPLError( CE_Failure, CPLE_FileIO,
293 0 : "Unable to read file section size.\n" );
294 0 : return NULL;
295 : }
296 :
297 : CPL_LSBPTR32( &nSize );
298 :
299 1 : if( VSIFReadL( (void *)&nVersion, sizeof(GInt32), 1, poDS->fp ) != 1 )
300 : {
301 0 : delete poDS;
302 : CPLError( CE_Failure, CPLE_FileIO,
303 0 : "Unable to read file version.\n" );
304 0 : return NULL;
305 : }
306 :
307 : CPL_LSBPTR32( &nVersion );
308 :
309 1 : if(nVersion != 1 && nVersion != 2)
310 : {
311 0 : delete poDS;
312 : CPLError( CE_Failure, CPLE_FileIO,
313 0 : "Incorrect file version (%d).", nVersion );
314 0 : return NULL;
315 : }
316 :
317 : // advance until the grid tag is found
318 3 : while(nTag != nGRID_TAG)
319 : {
320 1 : if( VSIFReadL( (void *)&nTag, sizeof(GInt32), 1, poDS->fp ) != 1 )
321 : {
322 0 : delete poDS;
323 0 : CPLError( CE_Failure, CPLE_FileIO, "Unable to read Tag.\n" );
324 0 : return NULL;
325 : }
326 :
327 : CPL_LSBPTR32( &nTag );
328 :
329 1 : if( VSIFReadL( (void *)&nSize, sizeof(GInt32), 1, poDS->fp ) != 1 )
330 : {
331 0 : delete poDS;
332 : CPLError( CE_Failure, CPLE_FileIO,
333 0 : "Unable to read file section size.\n" );
334 0 : return NULL;
335 : }
336 :
337 : CPL_LSBPTR32( &nSize );
338 :
339 1 : if(nTag != nGRID_TAG)
340 : {
341 0 : if( VSIFSeekL( poDS->fp, nSize, SEEK_SET ) != 0 )
342 : {
343 0 : delete poDS;
344 : CPLError( CE_Failure, CPLE_FileIO,
345 0 : "Unable to seek to end of file section.\n" );
346 0 : return NULL;
347 : }
348 : }
349 : }
350 :
351 : /* --------------------------------------------------------------------*/
352 : /* Read the grid. */
353 : /* --------------------------------------------------------------------*/
354 : /* Parse number of Y axis grid rows */
355 : GInt32 nRows;
356 1 : if( VSIFReadL( (void *)&nRows, sizeof(GInt32), 1, poDS->fp ) != 1 )
357 : {
358 0 : delete poDS;
359 : CPLError( CE_Failure, CPLE_FileIO,
360 0 : "Unable to read raster Y size.\n" );
361 0 : return NULL;
362 : }
363 : CPL_LSBPTR32( &nRows );
364 1 : poDS->nRasterYSize = nRows;
365 :
366 : /* Parse number of X axis grid columns */
367 : GInt32 nCols;
368 1 : if( VSIFReadL( (void *)&nCols, sizeof(GInt32), 1, poDS->fp ) != 1 )
369 : {
370 0 : delete poDS;
371 : CPLError( CE_Failure, CPLE_FileIO,
372 0 : "Unable to read raster X size.\n" );
373 0 : return NULL;
374 : }
375 : CPL_LSBPTR32( &nCols );
376 1 : poDS->nRasterXSize = nCols;
377 :
378 1 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
379 : {
380 0 : delete poDS;
381 0 : return NULL;
382 : }
383 :
384 : /* --------------------------------------------------------------------*/
385 : /* Create band information objects. */
386 : /* --------------------------------------------------------------------*/
387 1 : GS7BGRasterBand *poBand = new GS7BGRasterBand( poDS, 1 );
388 :
389 : // find the min X Value of the grid
390 : double dfTemp;
391 1 : if( VSIFReadL( (void *)&dfTemp, sizeof(double), 1, poDS->fp ) != 1 )
392 : {
393 0 : delete poDS;
394 : CPLError( CE_Failure, CPLE_FileIO,
395 0 : "Unable to read minimum X value.\n" );
396 0 : return NULL;
397 : }
398 : CPL_LSBPTR64( &dfTemp );
399 1 : poBand->dfMinX = dfTemp;
400 :
401 : // find the min Y value of the grid
402 1 : if( VSIFReadL( (void *)&dfTemp, sizeof(double), 1, poDS->fp ) != 1 )
403 : {
404 0 : delete poDS;
405 : CPLError( CE_Failure, CPLE_FileIO,
406 0 : "Unable to read minimum X value.\n" );
407 0 : return NULL;
408 : }
409 : CPL_LSBPTR64( &dfTemp );
410 1 : poBand->dfMinY = dfTemp;
411 :
412 : // find the spacing between adjacent nodes in the X direction
413 : // (between columns)
414 1 : if( VSIFReadL( (void *)&dfTemp, sizeof(double), 1, poDS->fp ) != 1 )
415 : {
416 0 : delete poDS;
417 : CPLError( CE_Failure, CPLE_FileIO,
418 0 : "Unable to read spacing in X value.\n" );
419 0 : return NULL;
420 : }
421 : CPL_LSBPTR64( &dfTemp );
422 1 : poBand->dfMaxX = poBand->dfMinX + (dfTemp * (nCols - 1));
423 :
424 : // find the spacing between adjacent nodes in the Y direction
425 : // (between rows)
426 1 : if( VSIFReadL( (void *)&dfTemp, sizeof(double), 1, poDS->fp ) != 1 )
427 : {
428 0 : delete poDS;
429 : CPLError( CE_Failure, CPLE_FileIO,
430 0 : "Unable to read spacing in Y value.\n" );
431 0 : return NULL;
432 : }
433 : CPL_LSBPTR64( &dfTemp );
434 1 : poBand->dfMaxY = poBand->dfMinY + (dfTemp * (nRows - 1));
435 :
436 : // set the z min
437 1 : if( VSIFReadL( (void *)&dfTemp, sizeof(double), 1, poDS->fp ) != 1 )
438 : {
439 0 : delete poDS;
440 : CPLError( CE_Failure, CPLE_FileIO,
441 0 : "Unable to read Z min value.\n" );
442 0 : return NULL;
443 : }
444 : CPL_LSBPTR64( &dfTemp );
445 1 : poBand->dfMinZ = dfTemp;
446 :
447 : // set the z max
448 1 : if( VSIFReadL( (void *)&dfTemp, sizeof(double), 1, poDS->fp ) != 1 )
449 : {
450 0 : delete poDS;
451 : CPLError( CE_Failure, CPLE_FileIO,
452 0 : "Unable to read Z max value.\n" );
453 0 : return NULL;
454 : }
455 : CPL_LSBPTR64( &dfTemp );
456 1 : poBand->dfMaxZ = dfTemp;
457 1 : poDS->SetBand( 1, poBand );
458 :
459 : // read and ignore the rotation value
460 : //(This is not used in the current version).
461 1 : if( VSIFReadL( (void *)&dfTemp, sizeof(double), 1, poDS->fp ) != 1 )
462 : {
463 0 : delete poDS;
464 : CPLError( CE_Failure, CPLE_FileIO,
465 0 : "Unable to read rotation value.\n" );
466 0 : return NULL;
467 : }
468 :
469 : // read and set the cell blank value
470 1 : if( VSIFReadL( (void *)&dfTemp, sizeof(double), 1, poDS->fp ) != 1 )
471 : {
472 0 : delete poDS;
473 : CPLError( CE_Failure, CPLE_FileIO,
474 0 : "Unable to Blank value.\n" );
475 0 : return NULL;
476 : }
477 : CPL_LSBPTR64( &dfTemp );
478 1 : poDS->dfNoData_Value = dfTemp;
479 :
480 : /* --------------------------------------------------------------------*/
481 : /* Set the current offset of the grid data. */
482 : /* --------------------------------------------------------------------*/
483 1 : if( VSIFReadL( (void *)&nTag, sizeof(GInt32), 1, poDS->fp ) != 1 )
484 : {
485 0 : delete poDS;
486 0 : CPLError( CE_Failure, CPLE_FileIO, "Unable to read Tag.\n" );
487 0 : return NULL;
488 : }
489 :
490 : CPL_LSBPTR32( &nTag );
491 1 : if(nTag != nDATA_TAG)
492 : {
493 0 : delete poDS;
494 0 : CPLError( CE_Failure, CPLE_FileIO, "Data tag not found.\n" );
495 0 : return NULL;
496 : }
497 :
498 1 : if( VSIFReadL( (void *)&nSize, sizeof(GInt32), 1, poDS->fp ) != 1 )
499 : {
500 0 : delete poDS;
501 : CPLError( CE_Failure, CPLE_FileIO,
502 0 : "Unable to data section size.\n" );
503 0 : return NULL;
504 : }
505 :
506 1 : poDS->nData_Position = VSIFTellL(poDS->fp);
507 :
508 :
509 : /* --------------------------------------------------------------------*/
510 : /* Initialize any PAM information. */
511 : /* --------------------------------------------------------------------*/
512 1 : poDS->SetDescription( poOpenInfo->pszFilename );
513 1 : poDS->TryLoadXML();
514 :
515 1 : return poDS;
516 : }
517 :
518 : /************************************************************************/
519 : /* GetGeoTransform() */
520 : /************************************************************************/
521 :
522 1 : CPLErr GS7BGDataset::GetGeoTransform( double *padfGeoTransform )
523 : {
524 1 : if( padfGeoTransform == NULL )
525 0 : return CE_Failure;
526 :
527 1 : GS7BGRasterBand *poGRB = (GS7BGRasterBand *)GetRasterBand( 1 );
528 :
529 1 : if( poGRB == NULL )
530 : {
531 0 : padfGeoTransform[0] = 0;
532 0 : padfGeoTransform[1] = 1;
533 0 : padfGeoTransform[2] = 0;
534 0 : padfGeoTransform[3] = 0;
535 0 : padfGeoTransform[4] = 0;
536 0 : padfGeoTransform[5] = 1;
537 0 : return CE_Failure;
538 : }
539 :
540 : /* check if we have a PAM GeoTransform stored */
541 1 : CPLPushErrorHandler( CPLQuietErrorHandler );
542 1 : CPLErr eErr = GDALPamDataset::GetGeoTransform( padfGeoTransform );
543 1 : CPLPopErrorHandler();
544 :
545 1 : if( eErr == CE_None )
546 0 : return CE_None;
547 :
548 : /* calculate pixel size first */
549 1 : padfGeoTransform[1] = (poGRB->dfMaxX - poGRB->dfMinX)/(nRasterXSize - 1);
550 1 : padfGeoTransform[5] = (poGRB->dfMinY - poGRB->dfMaxY)/(nRasterYSize - 1);
551 :
552 : /* then calculate image origin */
553 1 : padfGeoTransform[0] = poGRB->dfMinX - padfGeoTransform[1] / 2;
554 1 : padfGeoTransform[3] = poGRB->dfMaxY - padfGeoTransform[5] / 2;
555 :
556 : /* tilt/rotation does not supported by the GS grids */
557 1 : padfGeoTransform[4] = 0.0;
558 1 : padfGeoTransform[2] = 0.0;
559 :
560 1 : return CE_None;
561 : }
562 :
563 : /************************************************************************/
564 : /* GDALRegister_GS7BG() */
565 : /************************************************************************/
566 338 : void GDALRegister_GS7BG()
567 :
568 : {
569 : GDALDriver *poDriver;
570 :
571 338 : if( GDALGetDriverByName( "GS7BG" ) == NULL )
572 : {
573 336 : poDriver = new GDALDriver();
574 :
575 336 : poDriver->SetDescription( "GS7BG" );
576 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
577 336 : "Golden Software 7 Binary Grid (.grd)" );
578 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
579 336 : "frmt_various.html#GS7BG" );
580 336 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "grd" );
581 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
582 336 : "Byte Int16 UInt16 Float32 Float64" );
583 :
584 336 : poDriver->pfnOpen = GS7BGDataset::Open;
585 :
586 336 : GetGDALDriverManager()->RegisterDriver( poDriver );
587 : }
588 338 : }
589 :
|
