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