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