1 : /******************************************************************************
2 : * $Id: gsbgdataset.cpp 24585 2012-06-16 10:27:04Z rouault $
3 : *
4 : * Project: GDAL
5 : * Purpose: Implements the Golden Software Binary Grid Format.
6 : * Author: Kevin Locke, kwl7@cornell.edu
7 : * (Based largely on aaigriddataset.cpp by Frank Warmerdam)
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 2006, Kevin Locke <kwl7@cornell.edu>
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 "cpl_conv.h"
32 :
33 : #include <float.h>
34 : #include <limits.h>
35 : #include <assert.h>
36 :
37 : #include "gdal_pam.h"
38 :
39 : #ifndef DBL_MAX
40 : # ifdef __DBL_MAX__
41 : # define DBL_MAX __DBL_MAX__
42 : # else
43 : # define DBL_MAX 1.7976931348623157E+308
44 : # endif /* __DBL_MAX__ */
45 : #endif /* DBL_MAX */
46 :
47 : #ifndef FLT_MAX
48 : # ifdef __FLT_MAX__
49 : # define FLT_MAX __FLT_MAX__
50 : # else
51 : # define FLT_MAX 3.40282347E+38F
52 : # endif /* __FLT_MAX__ */
53 : #endif /* FLT_MAX */
54 :
55 : #ifndef INT_MAX
56 : # define INT_MAX 2147483647
57 : #endif /* INT_MAX */
58 :
59 : #ifndef SHRT_MAX
60 : # define SHRT_MAX 32767
61 : #endif /* SHRT_MAX */
62 :
63 : CPL_CVSID("$Id: gsbgdataset.cpp 24585 2012-06-16 10:27:04Z rouault $");
64 :
65 : CPL_C_START
66 : void GDALRegister_GSBG(void);
67 : CPL_C_END
68 :
69 : /************************************************************************/
70 : /* ==================================================================== */
71 : /* GSBGDataset */
72 : /* ==================================================================== */
73 : /************************************************************************/
74 :
75 : class GSBGRasterBand;
76 :
77 : class GSBGDataset : public GDALPamDataset
78 44 : {
79 : friend class GSBGRasterBand;
80 :
81 : static const float fNODATA_VALUE;
82 : static const size_t nHEADER_SIZE;
83 :
84 : static CPLErr WriteHeader( VSILFILE *fp, GInt16 nXSize, GInt16 nYSize,
85 : double dfMinX, double dfMaxX,
86 : double dfMinY, double dfMaxY,
87 : double dfMinZ, double dfMaxZ );
88 :
89 : VSILFILE *fp;
90 :
91 : public:
92 : ~GSBGDataset();
93 :
94 : static int Identify( GDALOpenInfo * );
95 : static GDALDataset *Open( GDALOpenInfo * );
96 : static GDALDataset *Create( const char * pszFilename,
97 : int nXSize, int nYSize, int nBands,
98 : GDALDataType eType,
99 : char **papszParmList );
100 : static GDALDataset *CreateCopy( const char *pszFilename,
101 : GDALDataset *poSrcDS,
102 : int bStrict, char **papszOptions,
103 : GDALProgressFunc pfnProgress,
104 : void *pProgressData );
105 :
106 : CPLErr GetGeoTransform( double *padfGeoTransform );
107 : CPLErr SetGeoTransform( double *padfGeoTransform );
108 : };
109 :
110 : /* NOTE: This is not mentioned in the spec, but Surfer 8 uses this value */
111 : /* 0x7effffee (Little Endian: eeffff7e) */
112 : const float GSBGDataset::fNODATA_VALUE = 1.701410009187828e+38f;
113 :
114 : const size_t GSBGDataset::nHEADER_SIZE = 56;
115 :
116 : /************************************************************************/
117 : /* ==================================================================== */
118 : /* GSBGRasterBand */
119 : /* ==================================================================== */
120 : /************************************************************************/
121 :
122 : class GSBGRasterBand : public GDALPamRasterBand
123 : {
124 : friend class GSBGDataset;
125 :
126 : double dfMinX;
127 : double dfMaxX;
128 : double dfMinY;
129 : double dfMaxY;
130 : double dfMinZ;
131 : double dfMaxZ;
132 :
133 : float *pafRowMinZ;
134 : float *pafRowMaxZ;
135 : int nMinZRow;
136 : int nMaxZRow;
137 :
138 : CPLErr ScanForMinMaxZ();
139 :
140 : public:
141 :
142 : GSBGRasterBand( GSBGDataset *, int );
143 : ~GSBGRasterBand();
144 :
145 : CPLErr IReadBlock( int, int, void * );
146 : CPLErr IWriteBlock( int, int, void * );
147 :
148 : double GetNoDataValue( int *pbSuccess = NULL );
149 : double GetMinimum( int *pbSuccess = NULL );
150 : double GetMaximum( int *pbSuccess = NULL );
151 : };
152 :
153 : /************************************************************************/
154 : /* GSBGRasterBand() */
155 : /************************************************************************/
156 :
157 44 : GSBGRasterBand::GSBGRasterBand( GSBGDataset *poDS, int nBand ) :
158 : pafRowMinZ(NULL),
159 : pafRowMaxZ(NULL),
160 : nMinZRow(-1),
161 44 : nMaxZRow(-1)
162 :
163 : {
164 44 : this->poDS = poDS;
165 44 : this->nBand = nBand;
166 :
167 44 : eDataType = GDT_Float32;
168 :
169 44 : nBlockXSize = poDS->GetRasterXSize();
170 44 : nBlockYSize = 1;
171 44 : }
172 :
173 : /************************************************************************/
174 : /* ~GSBGRasterBand() */
175 : /************************************************************************/
176 :
177 44 : GSBGRasterBand::~GSBGRasterBand( )
178 :
179 : {
180 44 : if( pafRowMinZ != NULL )
181 1 : CPLFree( pafRowMinZ );
182 44 : if( pafRowMaxZ != NULL )
183 1 : CPLFree( pafRowMaxZ );
184 44 : }
185 :
186 : /************************************************************************/
187 : /* ScanForMinMaxZ() */
188 : /************************************************************************/
189 :
190 1 : CPLErr GSBGRasterBand::ScanForMinMaxZ()
191 :
192 : {
193 1 : float *pafRowVals = (float *)VSIMalloc2( nRasterXSize, 4 );
194 :
195 1 : if( pafRowVals == NULL )
196 : {
197 : CPLError( CE_Failure, CPLE_OutOfMemory,
198 0 : "Unable to allocate row buffer to scan grid file.\n" );
199 0 : return CE_Failure;
200 : }
201 :
202 1 : double dfNewMinZ = DBL_MAX;
203 1 : double dfNewMaxZ = -DBL_MAX;
204 1 : int nNewMinZRow = 0;
205 1 : int nNewMaxZRow = 0;
206 :
207 : /* Since we have to scan, lets calc. statistics too */
208 1 : double dfSum = 0.0;
209 1 : double dfSum2 = 0.0;
210 1 : unsigned long nValuesRead = 0;
211 21 : for( int iRow=0; iRow<nRasterYSize; iRow++ )
212 : {
213 20 : CPLErr eErr = IReadBlock( 0, iRow, pafRowVals );
214 20 : if( eErr != CE_None )
215 : {
216 0 : VSIFree( pafRowVals );
217 0 : return CE_Failure;
218 : }
219 :
220 20 : pafRowMinZ[iRow] = FLT_MAX;
221 20 : pafRowMaxZ[iRow] = -FLT_MAX;
222 420 : for( int iCol=0; iCol<nRasterXSize; iCol++ )
223 : {
224 400 : if( pafRowVals[iCol] == GSBGDataset::fNODATA_VALUE )
225 400 : continue;
226 :
227 0 : if( pafRowVals[iCol] < pafRowMinZ[iRow] )
228 0 : pafRowMinZ[iRow] = pafRowVals[iCol];
229 :
230 0 : if( pafRowVals[iCol] > pafRowMinZ[iRow] )
231 0 : pafRowMaxZ[iRow] = pafRowVals[iCol];
232 :
233 0 : dfSum += pafRowVals[iCol];
234 0 : dfSum2 += pafRowVals[iCol] * pafRowVals[iCol];
235 0 : nValuesRead++;
236 : }
237 :
238 20 : if( pafRowMinZ[iRow] < dfNewMinZ )
239 : {
240 1 : dfNewMinZ = pafRowMinZ[iRow];
241 1 : nNewMinZRow = iRow;
242 : }
243 :
244 20 : if( pafRowMaxZ[iRow] > dfNewMaxZ )
245 : {
246 1 : dfNewMaxZ = pafRowMaxZ[iRow];
247 1 : nNewMaxZRow = iRow;
248 : }
249 : }
250 :
251 1 : VSIFree( pafRowVals );
252 :
253 1 : if( nValuesRead == 0 )
254 : {
255 1 : dfMinZ = 0.0;
256 1 : dfMaxZ = 0.0;
257 1 : nMinZRow = 0;
258 1 : nMaxZRow = 0;
259 1 : return CE_None;
260 : }
261 :
262 0 : dfMinZ = dfNewMinZ;
263 0 : dfMaxZ = dfNewMaxZ;
264 0 : nMinZRow = nNewMinZRow;
265 0 : nMaxZRow = nNewMaxZRow;
266 :
267 0 : double dfMean = dfSum / nValuesRead;
268 0 : double dfStdDev = sqrt((dfSum2 / nValuesRead) - (dfMean * dfMean));
269 0 : SetStatistics( dfMinZ, dfMaxZ, dfMean, dfStdDev );
270 :
271 0 : return CE_None;
272 : }
273 :
274 : /************************************************************************/
275 : /* IReadBlock() */
276 : /************************************************************************/
277 :
278 140 : CPLErr GSBGRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
279 : void * pImage )
280 :
281 : {
282 140 : if( nBlockYOff < 0 || nBlockYOff > nRasterYSize - 1 || nBlockXOff != 0 )
283 0 : return CE_Failure;
284 :
285 140 : GSBGDataset *poGDS = dynamic_cast<GSBGDataset *>(poDS);
286 140 : if( VSIFSeekL( poGDS->fp,
287 : GSBGDataset::nHEADER_SIZE +
288 : 4 * nRasterXSize * (nRasterYSize - nBlockYOff - 1),
289 : SEEK_SET ) != 0 )
290 : {
291 : CPLError( CE_Failure, CPLE_FileIO,
292 0 : "Unable to seek to beginning of grid row.\n" );
293 0 : return CE_Failure;
294 : }
295 :
296 140 : if( VSIFReadL( pImage, sizeof(float), nBlockXSize,
297 : poGDS->fp ) != static_cast<unsigned>(nBlockXSize) )
298 : {
299 : CPLError( CE_Failure, CPLE_FileIO,
300 0 : "Unable to read block from grid file.\n" );
301 0 : return CE_Failure;
302 : }
303 :
304 : float *pfImage;
305 140 : pfImage = (float *)pImage;
306 140 : for( int iPixel=0; iPixel<nBlockXSize; iPixel++ )
307 : CPL_LSBPTR32( pfImage+iPixel );
308 :
309 140 : return CE_None;
310 : }
311 :
312 : /************************************************************************/
313 : /* IWriteBlock() */
314 : /************************************************************************/
315 :
316 20 : CPLErr GSBGRasterBand::IWriteBlock( int nBlockXOff, int nBlockYOff,
317 : void *pImage )
318 :
319 : {
320 20 : if( eAccess == GA_ReadOnly )
321 : {
322 : CPLError( CE_Failure, CPLE_NoWriteAccess,
323 0 : "Unable to write block, dataset opened read only.\n" );
324 0 : return CE_Failure;
325 : }
326 :
327 20 : if( nBlockYOff < 0 || nBlockYOff > nRasterYSize - 1 || nBlockXOff != 0 )
328 0 : return CE_Failure;
329 :
330 20 : GSBGDataset *poGDS = dynamic_cast<GSBGDataset *>(poDS);
331 20 : assert( poGDS != NULL );
332 :
333 20 : if( pafRowMinZ == NULL || pafRowMaxZ == NULL
334 : || nMinZRow < 0 || nMaxZRow < 0 )
335 : {
336 1 : pafRowMinZ = (float *)VSIMalloc2( nRasterYSize,sizeof(float) );
337 1 : if( pafRowMinZ == NULL )
338 : {
339 : CPLError( CE_Failure, CPLE_OutOfMemory,
340 0 : "Unable to allocate space for row minimums array.\n" );
341 0 : return CE_Failure;
342 : }
343 :
344 1 : pafRowMaxZ = (float *)VSIMalloc2( nRasterYSize,sizeof(float) );
345 1 : if( pafRowMaxZ == NULL )
346 : {
347 0 : VSIFree( pafRowMinZ );
348 0 : pafRowMinZ = NULL;
349 : CPLError( CE_Failure, CPLE_OutOfMemory,
350 0 : "Unable to allocate space for row maximums array.\n" );
351 0 : return CE_Failure;
352 : }
353 :
354 1 : CPLErr eErr = ScanForMinMaxZ();
355 1 : if( eErr != CE_None )
356 0 : return eErr;
357 : }
358 :
359 20 : if( VSIFSeekL( poGDS->fp,
360 : GSBGDataset::nHEADER_SIZE +
361 : 4 * nRasterXSize * (nRasterYSize - nBlockYOff - 1),
362 : SEEK_SET ) != 0 )
363 : {
364 : CPLError( CE_Failure, CPLE_FileIO,
365 0 : "Unable to seek to beginning of grid row.\n" );
366 0 : return CE_Failure;
367 : }
368 :
369 20 : float *pfImage = (float *)pImage;
370 20 : pafRowMinZ[nBlockYOff] = FLT_MAX;
371 20 : pafRowMaxZ[nBlockYOff] = -FLT_MAX;
372 420 : for( int iPixel=0; iPixel<nBlockXSize; iPixel++ )
373 : {
374 400 : if( pfImage[iPixel] != GSBGDataset::fNODATA_VALUE )
375 : {
376 400 : if( pfImage[iPixel] < pafRowMinZ[nBlockYOff] )
377 78 : pafRowMinZ[nBlockYOff] = pfImage[iPixel];
378 :
379 400 : if( pfImage[iPixel] > pafRowMaxZ[nBlockYOff] )
380 43 : pafRowMaxZ[nBlockYOff] = pfImage[iPixel];
381 : }
382 :
383 : CPL_LSBPTR32( pfImage+iPixel );
384 : }
385 :
386 20 : if( VSIFWriteL( pImage, sizeof(float), nBlockXSize,
387 : poGDS->fp ) != static_cast<unsigned>(nBlockXSize) )
388 : {
389 : CPLError( CE_Failure, CPLE_FileIO,
390 0 : "Unable to write block to grid file.\n" );
391 0 : return CE_Failure;
392 : }
393 :
394 : /* Update min/max Z values as appropriate */
395 20 : bool bHeaderNeedsUpdate = false;
396 20 : if( nMinZRow == nBlockYOff && pafRowMinZ[nBlockYOff] > dfMinZ )
397 : {
398 1 : double dfNewMinZ = DBL_MAX;
399 21 : for( int iRow=0; iRow<nRasterYSize; iRow++ )
400 : {
401 20 : if( pafRowMinZ[iRow] < dfNewMinZ )
402 : {
403 1 : dfNewMinZ = pafRowMinZ[iRow];
404 1 : nMinZRow = iRow;
405 : }
406 : }
407 :
408 1 : if( dfNewMinZ != dfMinZ )
409 : {
410 1 : dfMinZ = dfNewMinZ;
411 1 : bHeaderNeedsUpdate = true;
412 : }
413 : }
414 :
415 20 : if( nMaxZRow == nBlockYOff && pafRowMaxZ[nBlockYOff] < dfMaxZ )
416 : {
417 0 : double dfNewMaxZ = -DBL_MAX;
418 0 : for( int iRow=0; iRow<nRasterYSize; iRow++ )
419 : {
420 0 : if( pafRowMaxZ[iRow] > dfNewMaxZ )
421 : {
422 0 : dfNewMaxZ = pafRowMaxZ[iRow];
423 0 : nMaxZRow = iRow;
424 : }
425 : }
426 :
427 0 : if( dfNewMaxZ != dfMaxZ )
428 : {
429 0 : dfMaxZ = dfNewMaxZ;
430 0 : bHeaderNeedsUpdate = true;
431 : }
432 : }
433 :
434 20 : if( pafRowMinZ[nBlockYOff] < dfMinZ || pafRowMaxZ[nBlockYOff] > dfMaxZ )
435 : {
436 6 : if( pafRowMinZ[nBlockYOff] < dfMinZ )
437 : {
438 3 : dfMinZ = pafRowMinZ[nBlockYOff];
439 3 : nMinZRow = nBlockYOff;
440 : }
441 :
442 6 : if( pafRowMaxZ[nBlockYOff] > dfMaxZ )
443 : {
444 4 : dfMaxZ = pafRowMaxZ[nBlockYOff];
445 4 : nMaxZRow = nBlockYOff;
446 : }
447 :
448 6 : bHeaderNeedsUpdate = true;
449 : }
450 :
451 20 : if( bHeaderNeedsUpdate && dfMaxZ > dfMinZ )
452 : {
453 : CPLErr eErr = poGDS->WriteHeader( poGDS->fp,
454 : (GInt16) nRasterXSize,
455 : (GInt16) nRasterYSize,
456 : dfMinX, dfMaxX,
457 : dfMinY, dfMaxY,
458 6 : dfMinZ, dfMaxZ );
459 6 : return eErr;
460 : }
461 :
462 14 : return CE_None;
463 : }
464 :
465 : /************************************************************************/
466 : /* GetNoDataValue() */
467 : /************************************************************************/
468 :
469 11 : double GSBGRasterBand::GetNoDataValue( int * pbSuccess )
470 : {
471 11 : if( pbSuccess )
472 10 : *pbSuccess = TRUE;
473 :
474 11 : return GSBGDataset::fNODATA_VALUE;
475 : }
476 :
477 : /************************************************************************/
478 : /* GetMinimum() */
479 : /************************************************************************/
480 :
481 0 : double GSBGRasterBand::GetMinimum( int *pbSuccess )
482 : {
483 0 : if( pbSuccess )
484 0 : *pbSuccess = TRUE;
485 :
486 0 : return dfMinZ;
487 : }
488 :
489 : /************************************************************************/
490 : /* GetMaximum() */
491 : /************************************************************************/
492 :
493 0 : double GSBGRasterBand::GetMaximum( int *pbSuccess )
494 : {
495 0 : if( pbSuccess )
496 0 : *pbSuccess = TRUE;
497 :
498 0 : return dfMaxZ;
499 : }
500 :
501 : /************************************************************************/
502 : /* ==================================================================== */
503 : /* GSBGDataset */
504 : /* ==================================================================== */
505 : /************************************************************************/
506 :
507 44 : GSBGDataset::~GSBGDataset()
508 :
509 : {
510 44 : FlushCache();
511 44 : if( fp != NULL )
512 44 : VSIFCloseL( fp );
513 44 : }
514 :
515 : /************************************************************************/
516 : /* Identify() */
517 : /************************************************************************/
518 :
519 12474 : int GSBGDataset::Identify( GDALOpenInfo * poOpenInfo )
520 :
521 : {
522 : /* Check for signature */
523 12474 : if( poOpenInfo->nHeaderBytes < 4
524 : || !EQUALN((const char *) poOpenInfo->pabyHeader,"DSBB",4) )
525 : {
526 12430 : return FALSE;
527 : }
528 :
529 44 : return TRUE;
530 : }
531 :
532 : /************************************************************************/
533 : /* Open() */
534 : /************************************************************************/
535 :
536 2381 : GDALDataset *GSBGDataset::Open( GDALOpenInfo * poOpenInfo )
537 :
538 : {
539 2381 : if( !Identify(poOpenInfo) )
540 : {
541 2337 : return NULL;
542 : }
543 :
544 : /* -------------------------------------------------------------------- */
545 : /* Create a corresponding GDALDataset. */
546 : /* -------------------------------------------------------------------- */
547 44 : GSBGDataset *poDS = new GSBGDataset();
548 :
549 : /* -------------------------------------------------------------------- */
550 : /* Open file with large file API. */
551 : /* -------------------------------------------------------------------- */
552 44 : poDS->eAccess = poOpenInfo->eAccess;
553 44 : if( poOpenInfo->eAccess == GA_ReadOnly )
554 19 : poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
555 : else
556 25 : poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "r+b" );
557 :
558 44 : if( poDS->fp == NULL )
559 : {
560 0 : delete poDS;
561 : CPLError( CE_Failure, CPLE_OpenFailed,
562 : "VSIFOpenL(%s) failed unexpectedly.",
563 0 : poOpenInfo->pszFilename );
564 0 : return NULL;
565 : }
566 :
567 : /* -------------------------------------------------------------------- */
568 : /* Read the header. */
569 : /* -------------------------------------------------------------------- */
570 44 : if( VSIFSeekL( poDS->fp, 4, SEEK_SET ) != 0 )
571 : {
572 0 : delete poDS;
573 : CPLError( CE_Failure, CPLE_FileIO,
574 0 : "Unable to seek to start of grid file header.\n" );
575 0 : return NULL;
576 : }
577 :
578 : /* Parse number of X axis grid rows */
579 : GInt16 nTemp;
580 44 : if( VSIFReadL( (void *)&nTemp, 2, 1, poDS->fp ) != 1 )
581 : {
582 0 : delete poDS;
583 0 : CPLError( CE_Failure, CPLE_FileIO, "Unable to read raster X size.\n" );
584 0 : return NULL;
585 : }
586 44 : poDS->nRasterXSize = CPL_LSBWORD16( nTemp );
587 :
588 44 : if( VSIFReadL( (void *)&nTemp, 2, 1, poDS->fp ) != 1 )
589 : {
590 0 : delete poDS;
591 0 : CPLError( CE_Failure, CPLE_FileIO, "Unable to read raster Y size.\n" );
592 0 : return NULL;
593 : }
594 44 : poDS->nRasterYSize = CPL_LSBWORD16( nTemp );
595 :
596 44 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
597 : {
598 0 : delete poDS;
599 0 : return NULL;
600 : }
601 :
602 : /* -------------------------------------------------------------------- */
603 : /* Create band information objects. */
604 : /* -------------------------------------------------------------------- */
605 44 : GSBGRasterBand *poBand = new GSBGRasterBand( poDS, 1 );
606 :
607 : double dfTemp;
608 44 : if( VSIFReadL( (void *)&dfTemp, 8, 1, poDS->fp ) != 1 )
609 : {
610 0 : delete poDS;
611 : CPLError( CE_Failure, CPLE_FileIO,
612 0 : "Unable to read minimum X value.\n" );
613 0 : return NULL;
614 : }
615 : CPL_LSBPTR64( &dfTemp );
616 44 : poBand->dfMinX = dfTemp;
617 :
618 44 : if( VSIFReadL( (void *)&dfTemp, 8, 1, poDS->fp ) != 1 )
619 : {
620 0 : delete poDS;
621 : CPLError( CE_Failure, CPLE_FileIO,
622 0 : "Unable to read maximum X value.\n" );
623 0 : return NULL;
624 : }
625 : CPL_LSBPTR64( &dfTemp );
626 44 : poBand->dfMaxX = dfTemp;
627 :
628 44 : if( VSIFReadL( (void *)&dfTemp, 8, 1, poDS->fp ) != 1 )
629 : {
630 0 : delete poDS;
631 : CPLError( CE_Failure, CPLE_FileIO,
632 0 : "Unable to read minimum Y value.\n" );
633 0 : return NULL;
634 : }
635 : CPL_LSBPTR64( &dfTemp );
636 44 : poBand->dfMinY = dfTemp;
637 :
638 44 : if( VSIFReadL( (void *)&dfTemp, 8, 1, poDS->fp ) != 1 )
639 : {
640 0 : delete poDS;
641 : CPLError( CE_Failure, CPLE_FileIO,
642 0 : "Unable to read maximum Y value.\n" );
643 0 : return NULL;
644 : }
645 : CPL_LSBPTR64( &dfTemp );
646 44 : poBand->dfMaxY = dfTemp;
647 :
648 44 : if( VSIFReadL( (void *)&dfTemp, 8, 1, poDS->fp ) != 1 )
649 : {
650 0 : delete poDS;
651 : CPLError( CE_Failure, CPLE_FileIO,
652 0 : "Unable to read minimum Z value.\n" );
653 0 : return NULL;
654 : }
655 : CPL_LSBPTR64( &dfTemp );
656 44 : poBand->dfMinZ = dfTemp;
657 :
658 44 : if( VSIFReadL( (void *)&dfTemp, 8, 1, poDS->fp ) != 1 )
659 : {
660 0 : delete poDS;
661 : CPLError( CE_Failure, CPLE_FileIO,
662 0 : "Unable to read maximum Z value.\n" );
663 0 : return NULL;
664 : }
665 : CPL_LSBPTR64( &dfTemp );
666 44 : poBand->dfMaxZ = dfTemp;
667 :
668 44 : poDS->SetBand( 1, poBand );
669 :
670 : /* -------------------------------------------------------------------- */
671 : /* Initialize any PAM information. */
672 : /* -------------------------------------------------------------------- */
673 44 : poDS->SetDescription( poOpenInfo->pszFilename );
674 44 : poDS->TryLoadXML();
675 :
676 : /* -------------------------------------------------------------------- */
677 : /* Check for external overviews. */
678 : /* -------------------------------------------------------------------- */
679 44 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename, poOpenInfo->papszSiblingFiles );
680 :
681 44 : return poDS;
682 : }
683 :
684 : /************************************************************************/
685 : /* GetGeoTransform() */
686 : /************************************************************************/
687 :
688 29 : CPLErr GSBGDataset::GetGeoTransform( double *padfGeoTransform )
689 : {
690 29 : if( padfGeoTransform == NULL )
691 0 : return CE_Failure;
692 :
693 29 : GSBGRasterBand *poGRB = dynamic_cast<GSBGRasterBand *>(GetRasterBand( 1 ));
694 :
695 29 : if( poGRB == NULL )
696 : {
697 0 : padfGeoTransform[0] = 0;
698 0 : padfGeoTransform[1] = 1;
699 0 : padfGeoTransform[2] = 0;
700 0 : padfGeoTransform[3] = 0;
701 0 : padfGeoTransform[4] = 0;
702 0 : padfGeoTransform[5] = 1;
703 0 : return CE_Failure;
704 : }
705 :
706 : /* check if we have a PAM GeoTransform stored */
707 29 : CPLPushErrorHandler( CPLQuietErrorHandler );
708 29 : CPLErr eErr = GDALPamDataset::GetGeoTransform( padfGeoTransform );
709 29 : CPLPopErrorHandler();
710 :
711 29 : if( eErr == CE_None )
712 0 : return CE_None;
713 :
714 : /* calculate pixel size first */
715 29 : padfGeoTransform[1] = (poGRB->dfMaxX - poGRB->dfMinX)/(nRasterXSize - 1);
716 29 : padfGeoTransform[5] = (poGRB->dfMinY - poGRB->dfMaxY)/(nRasterYSize - 1);
717 :
718 : /* then calculate image origin */
719 29 : padfGeoTransform[0] = poGRB->dfMinX - padfGeoTransform[1] / 2;
720 29 : padfGeoTransform[3] = poGRB->dfMaxY - padfGeoTransform[5] / 2;
721 :
722 : /* tilt/rotation does not supported by the GS grids */
723 29 : padfGeoTransform[4] = 0.0;
724 29 : padfGeoTransform[2] = 0.0;
725 :
726 29 : return CE_None;
727 : }
728 :
729 : /************************************************************************/
730 : /* SetGeoTransform() */
731 : /************************************************************************/
732 :
733 12 : CPLErr GSBGDataset::SetGeoTransform( double *padfGeoTransform )
734 : {
735 12 : if( eAccess == GA_ReadOnly )
736 : {
737 : CPLError( CE_Failure, CPLE_NoWriteAccess,
738 0 : "Unable to set GeoTransform, dataset opened read only.\n" );
739 0 : return CE_Failure;
740 : }
741 :
742 12 : GSBGRasterBand *poGRB = dynamic_cast<GSBGRasterBand *>(GetRasterBand( 1 ));
743 :
744 12 : if( poGRB == NULL || padfGeoTransform == NULL)
745 0 : return CE_Failure;
746 :
747 : /* non-zero transform 2 or 4 or negative 1 or 5 not supported natively */
748 12 : CPLErr eErr = CE_None;
749 : /*if( padfGeoTransform[2] != 0.0 || padfGeoTransform[4] != 0.0
750 : || padfGeoTransform[1] < 0.0 || padfGeoTransform[5] < 0.0 )
751 : eErr = GDALPamDataset::SetGeoTransform( padfGeoTransform );
752 :
753 : if( eErr != CE_None )
754 : return eErr;*/
755 :
756 12 : double dfMinX = padfGeoTransform[0] + padfGeoTransform[1] / 2;
757 : double dfMaxX =
758 12 : padfGeoTransform[1] * (nRasterXSize - 0.5) + padfGeoTransform[0];
759 : double dfMinY =
760 12 : padfGeoTransform[5] * (nRasterYSize - 0.5) + padfGeoTransform[3];
761 12 : double dfMaxY = padfGeoTransform[3] + padfGeoTransform[5] / 2;
762 :
763 : eErr = WriteHeader( fp,
764 : (GInt16) poGRB->nRasterXSize,
765 : (GInt16) poGRB->nRasterYSize,
766 : dfMinX, dfMaxX, dfMinY, dfMaxY,
767 12 : poGRB->dfMinZ, poGRB->dfMaxZ );
768 :
769 12 : if( eErr == CE_None )
770 : {
771 12 : poGRB->dfMinX = dfMinX;
772 12 : poGRB->dfMaxX = dfMaxX;
773 12 : poGRB->dfMinY = dfMinY;
774 12 : poGRB->dfMaxY = dfMaxY;
775 : }
776 :
777 12 : return eErr;
778 : }
779 :
780 : /************************************************************************/
781 : /* WriteHeader() */
782 : /************************************************************************/
783 :
784 55 : CPLErr GSBGDataset::WriteHeader( VSILFILE *fp, GInt16 nXSize, GInt16 nYSize,
785 : double dfMinX, double dfMaxX,
786 : double dfMinY, double dfMaxY,
787 : double dfMinZ, double dfMaxZ )
788 :
789 : {
790 55 : if( VSIFSeekL( fp, 0, SEEK_SET ) != 0 )
791 : {
792 : CPLError( CE_Failure, CPLE_FileIO,
793 0 : "Unable to seek to start of grid file.\n" );
794 0 : return CE_Failure;
795 : }
796 :
797 55 : if( VSIFWriteL( (void *)"DSBB", 1, 4, fp ) != 4 )
798 : {
799 : CPLError( CE_Failure, CPLE_FileIO,
800 0 : "Unable to write signature to grid file.\n" );
801 0 : return CE_Failure;
802 : }
803 :
804 55 : GInt16 nTemp = CPL_LSBWORD16(nXSize);
805 55 : if( VSIFWriteL( (void *)&nTemp, 2, 1, fp ) != 1 )
806 : {
807 : CPLError( CE_Failure, CPLE_FileIO,
808 0 : "Unable to write raster X size to grid file.\n" );
809 0 : return CE_Failure;
810 : }
811 :
812 55 : nTemp = CPL_LSBWORD16(nYSize);
813 55 : if( VSIFWriteL( (void *)&nTemp, 2, 1, fp ) != 1 )
814 : {
815 : CPLError( CE_Failure, CPLE_FileIO,
816 0 : "Unable to write raster Y size to grid file.\n" );
817 0 : return CE_Failure;
818 : }
819 :
820 55 : double dfTemp = dfMinX;
821 : CPL_LSBPTR64( &dfTemp );
822 55 : if( VSIFWriteL( (void *)&dfTemp, 8, 1, fp ) != 1 )
823 : {
824 : CPLError( CE_Failure, CPLE_FileIO,
825 0 : "Unable to write minimum X value to grid file.\n" );
826 0 : return CE_Failure;
827 : }
828 :
829 55 : dfTemp = dfMaxX;
830 : CPL_LSBPTR64( &dfTemp );
831 55 : if( VSIFWriteL( (void *)&dfTemp, 8, 1, fp ) != 1 )
832 : {
833 : CPLError( CE_Failure, CPLE_FileIO,
834 0 : "Unable to write maximum X value to grid file.\n" );
835 0 : return CE_Failure;
836 : }
837 :
838 55 : dfTemp = dfMinY;
839 : CPL_LSBPTR64( &dfTemp );
840 55 : if( VSIFWriteL( (void *)&dfTemp, 8, 1, fp ) != 1 )
841 : {
842 : CPLError( CE_Failure, CPLE_FileIO,
843 0 : "Unable to write minimum Y value to grid file.\n" );
844 0 : return CE_Failure;
845 : }
846 :
847 55 : dfTemp = dfMaxY;
848 : CPL_LSBPTR64( &dfTemp );
849 55 : if( VSIFWriteL( (void *)&dfTemp, 8, 1, fp ) != 1 )
850 : {
851 : CPLError( CE_Failure, CPLE_FileIO,
852 0 : "Unable to write maximum Y value to grid file.\n" );
853 0 : return CE_Failure;
854 : }
855 :
856 55 : dfTemp = dfMinZ;
857 : CPL_LSBPTR64( &dfTemp );
858 55 : if( VSIFWriteL( (void *)&dfTemp, 8, 1, fp ) != 1 )
859 : {
860 : CPLError( CE_Failure, CPLE_FileIO,
861 0 : "Unable to write minimum Z value to grid file.\n" );
862 0 : return CE_Failure;
863 : }
864 :
865 55 : dfTemp = dfMaxZ;
866 : CPL_LSBPTR64( &dfTemp );
867 55 : if( VSIFWriteL( (void *)&dfTemp, 8, 1, fp ) != 1 )
868 : {
869 : CPLError( CE_Failure, CPLE_FileIO,
870 0 : "Unable to write maximum Z value to grid file.\n" );
871 0 : return CE_Failure;
872 : }
873 :
874 55 : return CE_None;
875 : }
876 :
877 : /************************************************************************/
878 : /* Create() */
879 : /************************************************************************/
880 :
881 27 : GDALDataset *GSBGDataset::Create( const char * pszFilename,
882 : int nXSize, int nYSize, int nBands,
883 : GDALDataType eType,
884 : char **papszParmList )
885 :
886 : {
887 27 : if( nXSize <= 0 || nYSize <= 0 )
888 : {
889 : CPLError( CE_Failure, CPLE_IllegalArg,
890 : "Unable to create grid, both X and Y size must be "
891 0 : "non-negative.\n" );
892 :
893 0 : return NULL;
894 : }
895 27 : else if( nXSize > SHRT_MAX
896 : || nYSize > SHRT_MAX )
897 : {
898 : CPLError( CE_Failure, CPLE_IllegalArg,
899 : "Unable to create grid, Golden Software Binary Grid format "
900 : "only supports sizes up to %dx%d. %dx%d not supported.\n",
901 0 : SHRT_MAX, SHRT_MAX, nXSize, nYSize );
902 :
903 0 : return NULL;
904 : }
905 :
906 27 : if( eType != GDT_Byte && eType != GDT_Float32 && eType != GDT_UInt16
907 : && eType != GDT_Int16 )
908 : {
909 : CPLError( CE_Failure, CPLE_AppDefined,
910 : "Golden Software ASCII Grid only supports Byte, Int16, "
911 : "Uint16, and Float32 datatypes. Unable to create with "
912 14 : "type %s.\n", GDALGetDataTypeName( eType ) );
913 :
914 14 : return NULL;
915 : }
916 :
917 13 : VSILFILE *fp = VSIFOpenL( pszFilename, "w+b" );
918 :
919 13 : if( fp == NULL )
920 : {
921 : CPLError( CE_Failure, CPLE_OpenFailed,
922 : "Attempt to create file '%s' failed.\n",
923 0 : pszFilename );
924 0 : return NULL;
925 : }
926 :
927 : CPLErr eErr = WriteHeader( fp, (GInt16) nXSize, (GInt16) nYSize,
928 13 : 0.0, nXSize, 0.0, nYSize, 0.0, 0.0 );
929 13 : if( eErr != CE_None )
930 : {
931 0 : VSIFCloseL( fp );
932 0 : return NULL;
933 : }
934 :
935 13 : float fVal = fNODATA_VALUE;
936 : CPL_LSBPTR32( &fVal );
937 1233 : for( int iRow = 0; iRow < nYSize; iRow++ )
938 : {
939 121620 : for( int iCol=0; iCol<nXSize; iCol++ )
940 : {
941 120400 : if( VSIFWriteL( (void *)&fVal, 4, 1, fp ) != 1 )
942 : {
943 0 : VSIFCloseL( fp );
944 : CPLError( CE_Failure, CPLE_FileIO,
945 0 : "Unable to write grid cell. Disk full?\n" );
946 0 : return NULL;
947 : }
948 : }
949 : }
950 :
951 13 : VSIFCloseL( fp );
952 :
953 13 : return (GDALDataset *)GDALOpen( pszFilename, GA_Update );
954 : }
955 :
956 : /************************************************************************/
957 : /* CreateCopy() */
958 : /************************************************************************/
959 :
960 30 : GDALDataset *GSBGDataset::CreateCopy( const char *pszFilename,
961 : GDALDataset *poSrcDS,
962 : int bStrict, char **papszOptions,
963 : GDALProgressFunc pfnProgress,
964 : void *pProgressData )
965 : {
966 30 : if( pfnProgress == NULL )
967 0 : pfnProgress = GDALDummyProgress;
968 :
969 30 : int nBands = poSrcDS->GetRasterCount();
970 30 : if (nBands == 0)
971 : {
972 : CPLError( CE_Failure, CPLE_NotSupported,
973 1 : "GSBG driver does not support source dataset with zero band.\n");
974 1 : return NULL;
975 : }
976 29 : else if (nBands > 1)
977 : {
978 4 : if( bStrict )
979 : {
980 : CPLError( CE_Failure, CPLE_NotSupported,
981 : "Unable to create copy, Golden Software Binary Grid "
982 4 : "format only supports one raster band.\n" );
983 4 : return NULL;
984 : }
985 : else
986 : CPLError( CE_Warning, CPLE_NotSupported,
987 : "Golden Software Binary Grid format only supports one "
988 0 : "raster band, first band will be copied.\n" );
989 : }
990 :
991 25 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( 1 );
992 25 : if( poSrcBand->GetXSize() > SHRT_MAX
993 : || poSrcBand->GetYSize() > SHRT_MAX )
994 : {
995 : CPLError( CE_Failure, CPLE_IllegalArg,
996 : "Unable to create grid, Golden Software Binary Grid format "
997 : "only supports sizes up to %dx%d. %dx%d not supported.\n",
998 : SHRT_MAX, SHRT_MAX,
999 0 : poSrcBand->GetXSize(), poSrcBand->GetYSize() );
1000 :
1001 0 : return NULL;
1002 : }
1003 :
1004 25 : if( !pfnProgress( 0.0, NULL, pProgressData ) )
1005 : {
1006 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated\n" );
1007 0 : return NULL;
1008 : }
1009 :
1010 25 : VSILFILE *fp = VSIFOpenL( pszFilename, "w+b" );
1011 :
1012 25 : if( fp == NULL )
1013 : {
1014 : CPLError( CE_Failure, CPLE_OpenFailed,
1015 : "Attempt to create file '%s' failed.\n",
1016 13 : pszFilename );
1017 13 : return NULL;
1018 : }
1019 :
1020 12 : GInt16 nXSize = (GInt16) poSrcBand->GetXSize();
1021 12 : GInt16 nYSize = (GInt16) poSrcBand->GetYSize();
1022 : double adfGeoTransform[6];
1023 :
1024 12 : poSrcDS->GetGeoTransform( adfGeoTransform );
1025 :
1026 12 : double dfMinX = adfGeoTransform[0] + adfGeoTransform[1] / 2;
1027 12 : double dfMaxX = adfGeoTransform[1] * (nXSize - 0.5) + adfGeoTransform[0];
1028 12 : double dfMinY = adfGeoTransform[5] * (nYSize - 0.5) + adfGeoTransform[3];
1029 12 : double dfMaxY = adfGeoTransform[3] + adfGeoTransform[5] / 2;
1030 : CPLErr eErr = WriteHeader( fp, nXSize, nYSize,
1031 12 : dfMinX, dfMaxX, dfMinY, dfMaxY, 0.0, 0.0 );
1032 :
1033 12 : if( eErr != CE_None )
1034 : {
1035 0 : VSIFCloseL( fp );
1036 0 : return NULL;
1037 : }
1038 :
1039 : /* -------------------------------------------------------------------- */
1040 : /* Copy band data. */
1041 : /* -------------------------------------------------------------------- */
1042 12 : float *pfData = (float *)VSIMalloc2( nXSize, sizeof( float ) );
1043 12 : if( pfData == NULL )
1044 : {
1045 0 : VSIFCloseL( fp );
1046 : CPLError( CE_Failure, CPLE_OutOfMemory,
1047 0 : "Unable to create copy, unable to allocate line buffer.\n" );
1048 0 : return NULL;
1049 : }
1050 :
1051 : int bSrcHasNDValue;
1052 12 : float fSrcNoDataValue = (float) poSrcBand->GetNoDataValue( &bSrcHasNDValue );
1053 12 : double dfMinZ = DBL_MAX;
1054 12 : double dfMaxZ = -DBL_MAX;
1055 142 : for( GInt16 iRow = nYSize - 1; iRow >= 0; iRow-- )
1056 : {
1057 : eErr = poSrcBand->RasterIO( GF_Read, 0, iRow,
1058 : nXSize, 1, pfData,
1059 130 : nXSize, 1, GDT_Float32, 0, 0 );
1060 :
1061 130 : if( eErr != CE_None )
1062 : {
1063 0 : VSIFCloseL( fp );
1064 0 : VSIFree( pfData );
1065 0 : return NULL;
1066 : }
1067 :
1068 1630 : for( int iCol=0; iCol<nXSize; iCol++ )
1069 : {
1070 1500 : if( bSrcHasNDValue && pfData[iCol] == fSrcNoDataValue )
1071 : {
1072 0 : pfData[iCol] = fNODATA_VALUE;
1073 : }
1074 : else
1075 : {
1076 1500 : if( pfData[iCol] > dfMaxZ )
1077 14 : dfMaxZ = pfData[iCol];
1078 :
1079 1500 : if( pfData[iCol] < dfMinZ )
1080 19 : dfMinZ = pfData[iCol];
1081 : }
1082 :
1083 : CPL_LSBPTR32( pfData+iCol );
1084 : }
1085 :
1086 130 : if( VSIFWriteL( (void *)pfData, 4, nXSize,
1087 : fp ) != static_cast<unsigned>(nXSize) )
1088 : {
1089 0 : VSIFCloseL( fp );
1090 0 : VSIFree( pfData );
1091 : CPLError( CE_Failure, CPLE_FileIO,
1092 0 : "Unable to write grid row. Disk full?\n" );
1093 0 : return NULL;
1094 : }
1095 :
1096 130 : if( !pfnProgress( static_cast<double>(nYSize - iRow)/nYSize,
1097 : NULL, pProgressData ) )
1098 : {
1099 0 : VSIFCloseL( fp );
1100 0 : VSIFree( pfData );
1101 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1102 0 : return NULL;
1103 : }
1104 : }
1105 :
1106 12 : VSIFree( pfData );
1107 :
1108 : /* write out the min and max values */
1109 : eErr = WriteHeader( fp, nXSize, nYSize,
1110 12 : dfMinX, dfMaxX, dfMinY, dfMaxY, dfMinZ, dfMaxZ );
1111 :
1112 12 : if( eErr != CE_None )
1113 : {
1114 0 : VSIFCloseL( fp );
1115 0 : return NULL;
1116 : }
1117 :
1118 12 : VSIFCloseL( fp );
1119 :
1120 : GDALPamDataset *poDS = (GDALPamDataset *)GDALOpen( pszFilename,
1121 12 : GA_Update );
1122 12 : if (poDS)
1123 : {
1124 12 : poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
1125 : }
1126 12 : return poDS;
1127 : }
1128 :
1129 : /************************************************************************/
1130 : /* GDALRegister_GSBG() */
1131 : /************************************************************************/
1132 :
1133 582 : void GDALRegister_GSBG()
1134 :
1135 : {
1136 : GDALDriver *poDriver;
1137 :
1138 582 : if( GDALGetDriverByName( "GSBG" ) == NULL )
1139 : {
1140 561 : poDriver = new GDALDriver();
1141 :
1142 561 : poDriver->SetDescription( "GSBG" );
1143 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
1144 561 : "Golden Software Binary Grid (.grd)" );
1145 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
1146 561 : "frmt_various.html#GSBG" );
1147 561 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "grd" );
1148 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
1149 561 : "Byte Int16 UInt16 Float32" );
1150 561 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
1151 :
1152 561 : poDriver->pfnIdentify = GSBGDataset::Identify;
1153 561 : poDriver->pfnOpen = GSBGDataset::Open;
1154 561 : poDriver->pfnCreate = GSBGDataset::Create;
1155 561 : poDriver->pfnCreateCopy = GSBGDataset::CreateCopy;
1156 :
1157 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
1158 : }
1159 582 : }
|
