Coverage Report

Created: 2025-06-09 08:44

/src/gdal/frmts/bsb/bsbdataset.cpp
Line
Count
Source (jump to first uncovered line)
1
/******************************************************************************
2
 *
3
 * Project:  BSB Reader
4
 * Purpose:  BSBDataset implementation for BSB format.
5
 * Author:   Frank Warmerdam, warmerdam@pobox.com
6
 *
7
 ******************************************************************************
8
 * Copyright (c) 2001, Frank Warmerdam <warmerdam@pobox.com>
9
 * Copyright (c) 2008-2012, Even Rouault <even dot rouault at spatialys.com>
10
 *
11
 * SPDX-License-Identifier: MIT
12
 ****************************************************************************/
13
14
#include "bsb_read.h"
15
#include "cpl_string.h"
16
#include "gdal_frmts.h"
17
#include "gdal_pam.h"
18
#include "ogr_spatialref.h"
19
20
#include <cstdlib>
21
#include <algorithm>
22
23
// Disabled as people may worry about the BSB patent
24
// #define BSB_CREATE
25
26
/************************************************************************/
27
/* ==================================================================== */
28
/*                              BSBDataset                              */
29
/* ==================================================================== */
30
/************************************************************************/
31
32
class BSBRasterBand;
33
34
class BSBDataset final : public GDALPamDataset
35
{
36
    int nGCPCount;
37
    GDAL_GCP *pasGCPList;
38
    OGRSpatialReference m_oGCPSRS{};
39
40
    double adfGeoTransform[6];
41
    int bGeoTransformSet;
42
43
    void ScanForGCPs(bool isNos, const char *pszFilename);
44
    void ScanForGCPsNos(const char *pszFilename);
45
    void ScanForGCPsBSB();
46
47
    void ScanForCutline();
48
49
    static int IdentifyInternal(GDALOpenInfo *, bool &isNosOut);
50
51
  public:
52
    BSBDataset();
53
    ~BSBDataset() override;
54
55
    BSBInfo *psInfo;
56
57
    static GDALDataset *Open(GDALOpenInfo *);
58
    static int Identify(GDALOpenInfo *);
59
60
    int GetGCPCount() override;
61
    const OGRSpatialReference *GetSpatialRef() const override;
62
    const GDAL_GCP *GetGCPs() override;
63
64
    CPLErr GetGeoTransform(double *padfTransform) override;
65
    const OGRSpatialReference *GetGCPSpatialRef() const override;
66
};
67
68
/************************************************************************/
69
/* ==================================================================== */
70
/*                            BSBRasterBand                             */
71
/* ==================================================================== */
72
/************************************************************************/
73
74
class BSBRasterBand final : public GDALPamRasterBand
75
{
76
    GDALColorTable oCT;
77
78
  public:
79
    explicit BSBRasterBand(BSBDataset *);
80
81
    CPLErr IReadBlock(int, int, void *) override;
82
    GDALColorTable *GetColorTable() override;
83
    GDALColorInterp GetColorInterpretation() override;
84
};
85
86
/************************************************************************/
87
/*                           BSBRasterBand()                            */
88
/************************************************************************/
89
90
BSBRasterBand::BSBRasterBand(BSBDataset *poDSIn)
91
92
59
{
93
59
    poDS = poDSIn;
94
59
    nBand = 1;
95
96
59
    eDataType = GDT_Byte;
97
98
59
    nBlockXSize = poDS->GetRasterXSize();
99
59
    nBlockYSize = 1;
100
101
    // Note that the first color table entry is dropped, everything is
102
    // shifted down.
103
3.10k
    for (int i = 0; i < poDSIn->psInfo->nPCTSize - 1; i++)
104
3.04k
    {
105
3.04k
        GDALColorEntry oColor = {poDSIn->psInfo->pabyPCT[i * 3 + 0 + 3],
106
3.04k
                                 poDSIn->psInfo->pabyPCT[i * 3 + 1 + 3],
107
3.04k
                                 poDSIn->psInfo->pabyPCT[i * 3 + 2 + 3], 255};
108
109
3.04k
        oCT.SetColorEntry(i, &oColor);
110
3.04k
    }
111
59
}
112
113
/************************************************************************/
114
/*                             IReadBlock()                             */
115
/************************************************************************/
116
117
CPLErr BSBRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
118
                                 void *pImage)
119
2.90k
{
120
2.90k
    BSBDataset *poGDS = cpl::down_cast<BSBDataset *>(poDS);
121
2.90k
    GByte *pabyScanline = (GByte *)pImage;
122
123
2.90k
    if (BSBReadScanline(poGDS->psInfo, nBlockYOff, pabyScanline))
124
2.85k
    {
125
2.03M
        for (int i = 0; i < nBlockXSize; i++)
126
2.03M
        {
127
            /* The indices start at 1, except in case of some charts */
128
            /* where there are missing values, which are filled to 0 */
129
            /* by BSBReadScanline */
130
2.03M
            if (pabyScanline[i] > 0)
131
1.67M
                pabyScanline[i] -= 1;
132
2.03M
        }
133
134
2.85k
        return CE_None;
135
2.85k
    }
136
137
50
    return CE_Failure;
138
2.90k
}
139
140
/************************************************************************/
141
/*                           GetColorTable()                            */
142
/************************************************************************/
143
144
GDALColorTable *BSBRasterBand::GetColorTable()
145
146
0
{
147
0
    return &oCT;
148
0
}
149
150
/************************************************************************/
151
/*                       GetColorInterpretation()                       */
152
/************************************************************************/
153
154
GDALColorInterp BSBRasterBand::GetColorInterpretation()
155
156
0
{
157
0
    return GCI_PaletteIndex;
158
0
}
159
160
/************************************************************************/
161
/* ==================================================================== */
162
/*                              BSBDataset                              */
163
/* ==================================================================== */
164
/************************************************************************/
165
166
/************************************************************************/
167
/*                           BSBDataset()                               */
168
/************************************************************************/
169
170
BSBDataset::BSBDataset()
171
747
    : nGCPCount(0), pasGCPList(nullptr), bGeoTransformSet(FALSE),
172
747
      psInfo(nullptr)
173
747
{
174
747
    m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
175
747
    m_oGCPSRS.importFromWkt(SRS_WKT_WGS84_LAT_LONG);
176
177
747
    adfGeoTransform[0] = 0.0; /* X Origin (top left corner) */
178
747
    adfGeoTransform[1] = 1.0; /* X Pixel size */
179
747
    adfGeoTransform[2] = 0.0;
180
747
    adfGeoTransform[3] = 0.0; /* Y Origin (top left corner) */
181
747
    adfGeoTransform[4] = 0.0;
182
747
    adfGeoTransform[5] = 1.0; /* Y Pixel Size */
183
747
}
184
185
/************************************************************************/
186
/*                            ~BSBDataset()                             */
187
/************************************************************************/
188
189
BSBDataset::~BSBDataset()
190
191
747
{
192
747
    FlushCache(true);
193
194
747
    GDALDeinitGCPs(nGCPCount, pasGCPList);
195
747
    CPLFree(pasGCPList);
196
197
747
    if (psInfo != nullptr)
198
59
        BSBClose(psInfo);
199
747
}
200
201
/************************************************************************/
202
/*                          GetGeoTransform()                           */
203
/************************************************************************/
204
205
CPLErr BSBDataset::GetGeoTransform(double *padfTransform)
206
207
59
{
208
59
    memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
209
210
59
    if (bGeoTransformSet)
211
6
        return CE_None;
212
213
53
    return CE_Failure;
214
59
}
215
216
/************************************************************************/
217
/*                          GetSpatialRef()                             */
218
/************************************************************************/
219
220
const OGRSpatialReference *BSBDataset::GetSpatialRef() const
221
222
59
{
223
59
    if (bGeoTransformSet)
224
6
        return &m_oGCPSRS;
225
226
53
    return nullptr;
227
59
}
228
229
/************************************************************************/
230
/*                     GDALHeuristicDatelineWrap()                      */
231
/************************************************************************/
232
233
static void GDALHeuristicDatelineWrap(int nPointCount, double *padfX)
234
235
40
{
236
40
    if (nPointCount < 2)
237
0
        return;
238
239
    /* -------------------------------------------------------------------- */
240
    /*      Work out what the longitude range will be centering on the      */
241
    /*      prime meridian (-180 to 180) and centering on the dateline      */
242
    /*      (0 to 360).                                                     */
243
    /* -------------------------------------------------------------------- */
244
    /* Following inits are useless but keep GCC happy */
245
40
    double dfX_PM_Min = 0.0;
246
40
    double dfX_PM_Max = 0.0;
247
40
    double dfX_Dateline_Min = 0.0;
248
40
    double dfX_Dateline_Max = 0.0;
249
250
1.41k
    for (int i = 0; i < nPointCount; i++)
251
1.37k
    {
252
1.37k
        double dfX_PM = padfX[i];
253
1.37k
        if (dfX_PM > 180)
254
19
            dfX_PM -= 360.0;
255
256
1.37k
        double dfX_Dateline = padfX[i];
257
1.37k
        if (dfX_Dateline < 0)
258
1
            dfX_Dateline += 360.0;
259
260
1.37k
        if (i == 0)
261
40
        {
262
40
            dfX_PM_Min = dfX_PM;
263
40
            dfX_PM_Max = dfX_PM;
264
40
            dfX_Dateline_Min = dfX_Dateline;
265
40
            dfX_Dateline_Max = dfX_Dateline;
266
40
        }
267
1.33k
        else
268
1.33k
        {
269
1.33k
            dfX_PM_Min = std::min(dfX_PM_Min, dfX_PM);
270
1.33k
            dfX_PM_Max = std::max(dfX_PM_Max, dfX_PM);
271
1.33k
            dfX_Dateline_Min = std::min(dfX_Dateline_Min, dfX_Dateline);
272
1.33k
            dfX_Dateline_Max = std::max(dfX_Dateline_Max, dfX_Dateline);
273
1.33k
        }
274
1.37k
    }
275
276
    /* -------------------------------------------------------------------- */
277
    /*      Do nothing if the range is always fairly small - no apparent    */
278
    /*      wrapping issues.                                                */
279
    /* -------------------------------------------------------------------- */
280
40
    if ((dfX_PM_Max - dfX_PM_Min) < 270.0 &&
281
40
        (dfX_Dateline_Max - dfX_Dateline_Min) < 270.0)
282
31
        return;
283
284
    /* -------------------------------------------------------------------- */
285
    /*      Do nothing if both approach have a wide range - best not to    */
286
    /*      fiddle if we aren't sure we are improving things.               */
287
    /* -------------------------------------------------------------------- */
288
9
    if ((dfX_PM_Max - dfX_PM_Min) > 270.0 &&
289
9
        (dfX_Dateline_Max - dfX_Dateline_Min) > 270.0)
290
3
        return;
291
292
    /* -------------------------------------------------------------------- */
293
    /*      Pick which way to transform things.                             */
294
    /* -------------------------------------------------------------------- */
295
6
    bool bUsePMWrap;
296
297
6
    if ((dfX_PM_Max - dfX_PM_Min) > 270.0 &&
298
6
        (dfX_Dateline_Max - dfX_Dateline_Min) < 270.0)
299
0
        bUsePMWrap = false;
300
6
    else
301
6
        bUsePMWrap = true;
302
303
    /* -------------------------------------------------------------------- */
304
    /*      Apply rewrapping.                                               */
305
    /* -------------------------------------------------------------------- */
306
64
    for (int i = 0; i < nPointCount; i++)
307
58
    {
308
58
        if (bUsePMWrap)
309
58
        {
310
58
            if (padfX[i] > 180)
311
14
                padfX[i] -= 360.0;
312
58
        }
313
0
        else
314
0
        {
315
0
            if (padfX[i] < 0)
316
0
                padfX[i] += 360.0;
317
0
        }
318
58
    }
319
6
}
320
321
/************************************************************************/
322
/*                   GDALHeuristicDatelineWrapGCPs()                    */
323
/************************************************************************/
324
325
static void GDALHeuristicDatelineWrapGCPs(int nPointCount, GDAL_GCP *pasGCPList)
326
40
{
327
40
    std::vector<double> oadfX;
328
329
40
    oadfX.resize(nPointCount);
330
1.41k
    for (int i = 0; i < nPointCount; i++)
331
1.37k
        oadfX[i] = pasGCPList[i].dfGCPX;
332
333
40
    GDALHeuristicDatelineWrap(nPointCount, &(oadfX[0]));
334
335
1.41k
    for (int i = 0; i < nPointCount; i++)
336
1.37k
        pasGCPList[i].dfGCPX = oadfX[i];
337
40
}
338
339
/************************************************************************/
340
/*                            ScanForGCPs()                             */
341
/************************************************************************/
342
343
void BSBDataset::ScanForGCPs(bool isNos, const char *pszFilename)
344
345
59
{
346
    /* -------------------------------------------------------------------- */
347
    /*      Collect GCPs as appropriate to source.                          */
348
    /* -------------------------------------------------------------------- */
349
59
    nGCPCount = 0;
350
351
59
    if (isNos)
352
0
    {
353
0
        ScanForGCPsNos(pszFilename);
354
0
    }
355
59
    else
356
59
    {
357
59
        ScanForGCPsBSB();
358
59
    }
359
360
    /* -------------------------------------------------------------------- */
361
    /*      Apply heuristics to re-wrap GCPs to maintain continuity        */
362
    /*      over the international dateline.                                */
363
    /* -------------------------------------------------------------------- */
364
59
    if (nGCPCount > 1)
365
40
        GDALHeuristicDatelineWrapGCPs(nGCPCount, pasGCPList);
366
367
    /* -------------------------------------------------------------------- */
368
    /*      Collect coordinate system related parameters from header.       */
369
    /* -------------------------------------------------------------------- */
370
59
    const char *pszKNP = nullptr;
371
59
    const char *pszKNQ = nullptr;
372
373
5.54k
    for (int i = 0; psInfo->papszHeader[i] != nullptr; i++)
374
5.48k
    {
375
5.48k
        if (STARTS_WITH_CI(psInfo->papszHeader[i], "KNP/"))
376
99
        {
377
99
            pszKNP = psInfo->papszHeader[i];
378
99
            SetMetadataItem("BSB_KNP", pszKNP + 4);
379
99
        }
380
5.48k
        if (STARTS_WITH_CI(psInfo->papszHeader[i], "KNQ/"))
381
171
        {
382
171
            pszKNQ = psInfo->papszHeader[i];
383
171
            SetMetadataItem("BSB_KNQ", pszKNQ + 4);
384
171
        }
385
5.48k
    }
386
387
    /* -------------------------------------------------------------------- */
388
    /*      Can we derive a reasonable coordinate system definition for     */
389
    /*      this file?  For now we keep it simple, just handling            */
390
    /*      mercator. In the future we should consider others.              */
391
    /* -------------------------------------------------------------------- */
392
59
    CPLString osUnderlyingSRS;
393
59
    if (pszKNP != nullptr)
394
53
    {
395
53
        const char *pszPR = strstr(pszKNP, "PR=");
396
53
        const char *pszGD = strstr(pszKNP, "GD=");
397
53
        const char *pszGEOGCS = SRS_WKT_WGS84_LAT_LONG;
398
53
        CPLString osPP;
399
400
        // Capture the PP string.
401
53
        const char *pszValue = strstr(pszKNP, "PP=");
402
53
        const char *pszEnd = pszValue ? strstr(pszValue, ",") : nullptr;
403
53
        if (pszValue && pszEnd)
404
48
            osPP.assign(pszValue + 3, pszEnd - pszValue - 3);
405
406
        // Look at the datum
407
53
        if (pszGD == nullptr)
408
1
        {
409
            /* no match. We'll default to EPSG:4326 */
410
1
        }
411
52
        else if (STARTS_WITH_CI(pszGD, "GD=European 1950"))
412
0
        {
413
0
            pszGEOGCS =
414
0
                "GEOGCS[\"ED50\",DATUM[\"European_Datum_1950\",SPHEROID["
415
0
                "\"International "
416
0
                "1924\",6378388,297,AUTHORITY[\"EPSG\",\"7022\"]],TOWGS84[-87,-"
417
0
                "98,-121,0,0,0,0],AUTHORITY[\"EPSG\",\"6230\"]],PRIMEM["
418
0
                "\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\","
419
0
                "0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY["
420
0
                "\"EPSG\",\"4230\"]]";
421
0
        }
422
423
        // Look at the projection
424
53
        if (pszPR == nullptr)
425
1
        {
426
            /* no match */
427
1
        }
428
52
        else if (STARTS_WITH_CI(pszPR, "PR=MERCATOR") && nGCPCount > 0)
429
6
        {
430
            // We somewhat arbitrarily select our first GCPX as our
431
            // central meridian.  This is mostly helpful to ensure
432
            // that regions crossing the dateline will be contiguous
433
            // in mercator.
434
6
            osUnderlyingSRS.Printf(
435
6
                "PROJCS[\"Global "
436
6
                "Mercator\",%s,PROJECTION[\"Mercator_2SP\"],PARAMETER["
437
6
                "\"standard_parallel_1\",0],PARAMETER[\"latitude_of_origin\",0]"
438
6
                ",PARAMETER[\"central_meridian\",%d],PARAMETER[\"false_"
439
6
                "easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]"
440
6
                "]",
441
6
                pszGEOGCS, (int)pasGCPList[0].dfGCPX);
442
6
        }
443
444
46
        else if (STARTS_WITH_CI(pszPR, "PR=TRANSVERSE MERCATOR") &&
445
46
                 !osPP.empty())
446
9
        {
447
448
9
            osUnderlyingSRS.Printf(
449
9
                "PROJCS[\"unnamed\",%s,PROJECTION[\"Transverse_Mercator\"],"
450
9
                "PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_"
451
9
                "meridian\",%s],PARAMETER[\"scale_factor\",1],PARAMETER["
452
9
                "\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT["
453
9
                "\"Meter\",1]]",
454
9
                pszGEOGCS, osPP.c_str());
455
9
        }
456
457
37
        else if (STARTS_WITH_CI(pszPR, "PR=UNIVERSAL TRANSVERSE MERCATOR") &&
458
37
                 !osPP.empty())
459
0
        {
460
            // This is not *really* UTM unless the central meridian
461
            // matches a zone which it does not in some (most?) maps.
462
0
            osUnderlyingSRS.Printf(
463
0
                "PROJCS[\"unnamed\",%s,PROJECTION[\"Transverse_Mercator\"],"
464
0
                "PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_"
465
0
                "meridian\",%s],PARAMETER[\"scale_factor\",0.9996],PARAMETER["
466
0
                "\"false_easting\",500000],PARAMETER[\"false_northing\",0],"
467
0
                "UNIT[\"Meter\",1]]",
468
0
                pszGEOGCS, osPP.c_str());
469
0
        }
470
471
37
        else if (STARTS_WITH_CI(pszPR, "PR=POLYCONIC") && !osPP.empty())
472
0
        {
473
0
            osUnderlyingSRS.Printf(
474
0
                "PROJCS[\"unnamed\",%s,PROJECTION[\"Polyconic\"],PARAMETER["
475
0
                "\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",%s],"
476
0
                "PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0]"
477
0
                ",UNIT[\"Meter\",1]]",
478
0
                pszGEOGCS, osPP.c_str());
479
0
        }
480
481
37
        else if (STARTS_WITH_CI(pszPR, "PR=LAMBERT CONFORMAL CONIC") &&
482
37
                 !osPP.empty() && pszKNQ != nullptr)
483
14
        {
484
14
            CPLString osP2, osP3;
485
486
            // Capture the KNQ/P2 string.
487
14
            pszValue = strstr(pszKNQ, "P2=");
488
14
            if (pszValue)
489
13
                pszEnd = strstr(pszValue, ",");
490
14
            if (pszValue && pszEnd)
491
12
                osP2.assign(pszValue + 3, pszEnd - pszValue - 3);
492
493
            // Capture the KNQ/P3 string.
494
14
            pszValue = strstr(pszKNQ, "P3=");
495
14
            if (pszValue)
496
7
                pszEnd = strstr(pszValue, ",");
497
14
            if (pszValue)
498
7
            {
499
7
                if (pszEnd)
500
7
                    osP3.assign(pszValue + 3, pszEnd - pszValue - 3);
501
0
                else
502
0
                    osP3.assign(pszValue + 3);
503
7
            }
504
505
14
            if (!osP2.empty() && !osP3.empty())
506
7
                osUnderlyingSRS.Printf(
507
7
                    "PROJCS[\"unnamed\",%s,PROJECTION[\"Lambert_Conformal_"
508
7
                    "Conic_2SP\"],PARAMETER[\"standard_parallel_1\",%s],"
509
7
                    "PARAMETER[\"standard_parallel_2\",%s],PARAMETER["
510
7
                    "\"latitude_of_origin\",0.0],PARAMETER[\"central_"
511
7
                    "meridian\",%s],PARAMETER[\"false_easting\",0.0],PARAMETER["
512
7
                    "\"false_northing\",0.0],UNIT[\"Meter\",1]]",
513
7
                    pszGEOGCS, osP2.c_str(), osP3.c_str(), osPP.c_str());
514
14
        }
515
53
    }
516
517
    /* -------------------------------------------------------------------- */
518
    /*      If we got an alternate underlying coordinate system, try        */
519
    /*      converting the GCPs to that coordinate system.                  */
520
    /* -------------------------------------------------------------------- */
521
59
    if (osUnderlyingSRS.length() > 0)
522
22
    {
523
22
        OGRSpatialReference oGeog_SRS, oProjected_SRS;
524
525
22
        oProjected_SRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
526
22
        oProjected_SRS.SetFromUserInput(osUnderlyingSRS);
527
22
        oGeog_SRS.CopyGeogCSFrom(&oProjected_SRS);
528
22
        oGeog_SRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
529
530
22
        OGRCoordinateTransformation *poCT =
531
22
            OGRCreateCoordinateTransformation(&oGeog_SRS, &oProjected_SRS);
532
22
        if (poCT != nullptr)
533
20
        {
534
1.27k
            for (int i = 0; i < nGCPCount; i++)
535
1.25k
            {
536
1.25k
                poCT->Transform(1, &(pasGCPList[i].dfGCPX),
537
1.25k
                                &(pasGCPList[i].dfGCPY),
538
1.25k
                                &(pasGCPList[i].dfGCPZ));
539
1.25k
            }
540
541
20
            m_oGCPSRS.importFromWkt(osUnderlyingSRS.c_str());
542
543
20
            delete poCT;
544
20
        }
545
2
        else
546
2
            CPLErrorReset();
547
22
    }
548
549
    /* -------------------------------------------------------------------- */
550
    /*      Attempt to prepare a geotransform from the GCPs.                */
551
    /* -------------------------------------------------------------------- */
552
59
    if (GDALGCPsToGeoTransform(nGCPCount, pasGCPList, adfGeoTransform, FALSE))
553
6
    {
554
6
        bGeoTransformSet = TRUE;
555
6
    }
556
59
}
557
558
/************************************************************************/
559
/*                           ScanForGCPsNos()                           */
560
/*                                                                      */
561
/*      Nos files have an accompanying .geo file, that contains some    */
562
/*      of the information normally contained in the header section     */
563
/*      with BSB files. we try and open a file with the same name,      */
564
/*      but a .geo extension, and look for lines like...                */
565
/*      PointX=long lat line pixel    (using the same naming system     */
566
/*      as BSB) Point1=-22.0000 64.250000 197 744                       */
567
/************************************************************************/
568
569
void BSBDataset::ScanForGCPsNos(const char *pszFilename)
570
0
{
571
0
    const std::string extension = CPLGetExtensionSafe(pszFilename);
572
573
    // pseudointelligently try and guess whether we want a .geo or a .GEO
574
0
    std::string geofile;
575
0
    if (extension.size() >= 2 && extension[1] == 'O')
576
0
    {
577
0
        geofile = CPLResetExtensionSafe(pszFilename, "GEO");
578
0
    }
579
0
    else
580
0
    {
581
0
        geofile = CPLResetExtensionSafe(pszFilename, "geo");
582
0
    }
583
584
0
    FILE *gfp = VSIFOpen(geofile.c_str(), "r");  // Text files
585
0
    if (gfp == nullptr)
586
0
    {
587
0
        CPLError(CE_Failure, CPLE_OpenFailed,
588
0
                 "Couldn't find a matching .GEO file: %s", geofile.c_str());
589
0
        return;
590
0
    }
591
592
0
    char *thisLine = (char *)CPLMalloc(80);  // FIXME
593
594
    // Count the GCPs (reference points) and seek the file pointer 'gfp' to the
595
    // starting point
596
0
    int fileGCPCount = 0;
597
0
    while (fgets(thisLine, 80, gfp))
598
0
    {
599
0
        if (STARTS_WITH_CI(thisLine, "Point"))
600
0
            fileGCPCount++;
601
0
    }
602
0
    VSIRewind(gfp);
603
604
    // Memory has not been allocated to fileGCPCount yet
605
0
    pasGCPList = (GDAL_GCP *)CPLCalloc(sizeof(GDAL_GCP), fileGCPCount + 1);
606
607
0
    while (fgets(thisLine, 80, gfp))
608
0
    {
609
0
        if (STARTS_WITH_CI(thisLine, "Point"))
610
0
        {
611
            // got a point line, turn it into a gcp
612
0
            char **Tokens =
613
0
                CSLTokenizeStringComplex(thisLine, "= ", FALSE, FALSE);
614
0
            if (CSLCount(Tokens) >= 5)
615
0
            {
616
0
                GDALInitGCPs(1, pasGCPList + nGCPCount);
617
0
                pasGCPList[nGCPCount].dfGCPX = CPLAtof(Tokens[1]);
618
0
                pasGCPList[nGCPCount].dfGCPY = CPLAtof(Tokens[2]);
619
0
                pasGCPList[nGCPCount].dfGCPPixel = CPLAtof(Tokens[4]);
620
0
                pasGCPList[nGCPCount].dfGCPLine = CPLAtof(Tokens[3]);
621
622
0
                CPLFree(pasGCPList[nGCPCount].pszId);
623
0
                char szName[50];
624
0
                snprintf(szName, sizeof(szName), "GCP_%d", nGCPCount + 1);
625
0
                pasGCPList[nGCPCount].pszId = CPLStrdup(szName);
626
627
0
                nGCPCount++;
628
0
            }
629
0
            CSLDestroy(Tokens);
630
0
        }
631
0
    }
632
633
0
    CPLFree(thisLine);
634
0
    VSIFClose(gfp);
635
0
}
636
637
/************************************************************************/
638
/*                            ScanForGCPsBSB()                          */
639
/************************************************************************/
640
641
void BSBDataset::ScanForGCPsBSB()
642
59
{
643
    /* -------------------------------------------------------------------- */
644
    /*      Collect standalone GCPs.  They look like:                       */
645
    /*                                                                      */
646
    /*      REF/1,115,2727,32.346666666667,-60.881666666667                 */
647
    /*      REF/n,pixel,line,lat,long                                       */
648
    /* -------------------------------------------------------------------- */
649
59
    int fileGCPCount = 0;
650
651
    // Count the GCPs (reference points) in psInfo->papszHeader
652
5.54k
    for (int i = 0; psInfo->papszHeader[i] != nullptr; i++)
653
5.48k
        if (STARTS_WITH_CI(psInfo->papszHeader[i], "REF/"))
654
1.40k
            fileGCPCount++;
655
656
    // Memory has not been allocated to fileGCPCount yet
657
59
    pasGCPList = (GDAL_GCP *)CPLCalloc(sizeof(GDAL_GCP), fileGCPCount + 1);
658
659
5.54k
    for (int i = 0; psInfo->papszHeader[i] != nullptr; i++)
660
5.48k
    {
661
662
5.48k
        if (!STARTS_WITH_CI(psInfo->papszHeader[i], "REF/"))
663
4.07k
            continue;
664
665
1.40k
        char **papszTokens = CSLTokenizeStringComplex(
666
1.40k
            psInfo->papszHeader[i] + 4, ",", FALSE, FALSE);
667
668
1.40k
        if (CSLCount(papszTokens) > 4)
669
1.38k
        {
670
1.38k
            GDALInitGCPs(1, pasGCPList + nGCPCount);
671
672
1.38k
            pasGCPList[nGCPCount].dfGCPX = CPLAtof(papszTokens[4]);
673
1.38k
            pasGCPList[nGCPCount].dfGCPY = CPLAtof(papszTokens[3]);
674
1.38k
            pasGCPList[nGCPCount].dfGCPPixel = CPLAtof(papszTokens[1]);
675
1.38k
            pasGCPList[nGCPCount].dfGCPLine = CPLAtof(papszTokens[2]);
676
677
1.38k
            CPLFree(pasGCPList[nGCPCount].pszId);
678
1.38k
            if (CSLCount(papszTokens) > 5)
679
51
            {
680
51
                pasGCPList[nGCPCount].pszId = CPLStrdup(papszTokens[5]);
681
51
            }
682
1.33k
            else
683
1.33k
            {
684
1.33k
                char szName[50];
685
1.33k
                snprintf(szName, sizeof(szName), "GCP_%d", nGCPCount + 1);
686
1.33k
                pasGCPList[nGCPCount].pszId = CPLStrdup(szName);
687
1.33k
            }
688
689
1.38k
            nGCPCount++;
690
1.38k
        }
691
1.40k
        CSLDestroy(papszTokens);
692
1.40k
    }
693
59
}
694
695
/************************************************************************/
696
/*                            ScanForCutline()                          */
697
/************************************************************************/
698
699
void BSBDataset::ScanForCutline()
700
59
{
701
    /* PLY: Border Polygon Record - coordinates of the panel within the
702
     * raster image, given in chart datum lat/long. Any shape polygon.
703
     * They look like:
704
     *      PLY/1,32.346666666667,-60.881666666667
705
     *      PLY/n,lat,long
706
     *
707
     * If found then we return it via a BSB_CUTLINE metadata item as a WKT
708
     * POLYGON.
709
     */
710
711
59
    std::string wkt;
712
5.54k
    for (int i = 0; psInfo->papszHeader[i] != nullptr; i++)
713
5.48k
    {
714
5.48k
        if (!STARTS_WITH_CI(psInfo->papszHeader[i], "PLY/"))
715
5.08k
            continue;
716
717
395
        const CPLStringList aosTokens(
718
395
            CSLTokenizeString2(psInfo->papszHeader[i] + 4, ",", 0));
719
720
395
        if (aosTokens.size() >= 3)
721
317
        {
722
317
            if (wkt.empty())
723
33
                wkt = "POLYGON ((";
724
284
            else
725
284
                wkt += ',';
726
317
            wkt += aosTokens[2];
727
317
            wkt += ' ';
728
317
            wkt += aosTokens[1];
729
317
        }
730
395
    }
731
732
59
    if (!wkt.empty())
733
33
    {
734
33
        wkt += "))";
735
33
        SetMetadataItem("BSB_CUTLINE", wkt.c_str());
736
33
    }
737
59
}
738
739
/************************************************************************/
740
/*                          IdentifyInternal()                          */
741
/************************************************************************/
742
743
int BSBDataset::IdentifyInternal(GDALOpenInfo *poOpenInfo, bool &isNosOut)
744
745
914k
{
746
    /* -------------------------------------------------------------------- */
747
    /*      Check for BSB/ keyword.                                         */
748
    /* -------------------------------------------------------------------- */
749
914k
    isNosOut = false;
750
751
914k
    if (poOpenInfo->nHeaderBytes < 1000)
752
845k
        return FALSE;
753
754
69.4k
    int i = 0;
755
74.9M
    for (; i < poOpenInfo->nHeaderBytes - 4; i++)
756
74.9M
    {
757
74.9M
        if (poOpenInfo->pabyHeader[i + 0] == 'B' &&
758
74.9M
            poOpenInfo->pabyHeader[i + 1] == 'S' &&
759
74.9M
            poOpenInfo->pabyHeader[i + 2] == 'B' &&
760
74.9M
            poOpenInfo->pabyHeader[i + 3] == '/')
761
1.71k
            break;
762
74.9M
        if (poOpenInfo->pabyHeader[i + 0] == 'N' &&
763
74.9M
            poOpenInfo->pabyHeader[i + 1] == 'O' &&
764
74.9M
            poOpenInfo->pabyHeader[i + 2] == 'S' &&
765
74.9M
            poOpenInfo->pabyHeader[i + 3] == '/')
766
0
        {
767
0
            isNosOut = true;
768
0
            break;
769
0
        }
770
74.9M
        if (poOpenInfo->pabyHeader[i + 0] == 'W' &&
771
74.9M
            poOpenInfo->pabyHeader[i + 1] == 'X' &&
772
74.9M
            poOpenInfo->pabyHeader[i + 2] == '\\' &&
773
74.9M
            poOpenInfo->pabyHeader[i + 3] == '8')
774
0
            break;
775
74.9M
    }
776
777
69.4k
    if (i == poOpenInfo->nHeaderBytes - 4)
778
67.6k
        return FALSE;
779
780
    /* Additional test to avoid false positive. See #2881 */
781
1.71k
    const char *pszHeader =
782
1.71k
        reinterpret_cast<const char *>(poOpenInfo->pabyHeader);
783
1.71k
    const char *pszShiftedHeader = pszHeader + i;
784
1.71k
    const char *pszRA = strstr(pszShiftedHeader, "RA=");
785
1.71k
    if (pszRA == nullptr) /* This may be a NO1 file */
786
789
        pszRA = strstr(pszShiftedHeader, "[JF");
787
1.71k
    if (pszRA == nullptr)
788
186
        return FALSE;
789
1.52k
    if (pszRA - pszShiftedHeader > 100 && !strstr(pszHeader, "VER/") &&
790
1.52k
        !strstr(pszHeader, "KNP/") && !strstr(pszHeader, "KNQ/") &&
791
1.52k
        !strstr(pszHeader, "RGB/"))
792
31
        return FALSE;
793
794
1.49k
    return TRUE;
795
1.52k
}
796
797
/************************************************************************/
798
/*                              Identify()                              */
799
/************************************************************************/
800
801
int BSBDataset::Identify(GDALOpenInfo *poOpenInfo)
802
803
913k
{
804
913k
    bool isNos;
805
913k
    return IdentifyInternal(poOpenInfo, isNos);
806
913k
}
807
808
/************************************************************************/
809
/*                                Open()                                */
810
/************************************************************************/
811
812
GDALDataset *BSBDataset::Open(GDALOpenInfo *poOpenInfo)
813
814
747
{
815
747
    bool isNos = false;
816
747
    if (!IdentifyInternal(poOpenInfo, isNos))
817
0
        return nullptr;
818
819
747
    if (poOpenInfo->eAccess == GA_Update)
820
0
    {
821
0
        ReportUpdateNotSupportedByDriver("BSB");
822
0
        return nullptr;
823
0
    }
824
825
    /* -------------------------------------------------------------------- */
826
    /*      Create a corresponding GDALDataset.                             */
827
    /* -------------------------------------------------------------------- */
828
747
    BSBDataset *poDS = new BSBDataset();
829
830
    /* -------------------------------------------------------------------- */
831
    /*      Open the file.                                                  */
832
    /* -------------------------------------------------------------------- */
833
747
    poDS->psInfo = BSBOpen(poOpenInfo->pszFilename);
834
747
    if (poDS->psInfo == nullptr)
835
688
    {
836
688
        delete poDS;
837
688
        return nullptr;
838
688
    }
839
840
59
    poDS->nRasterXSize = poDS->psInfo->nXSize;
841
59
    poDS->nRasterYSize = poDS->psInfo->nYSize;
842
843
    /* -------------------------------------------------------------------- */
844
    /*      Create band information objects.                                */
845
    /* -------------------------------------------------------------------- */
846
59
    poDS->SetBand(1, new BSBRasterBand(poDS));
847
848
59
    poDS->ScanForGCPs(isNos, poOpenInfo->pszFilename);
849
850
    /* -------------------------------------------------------------------- */
851
    /*      Set CUTLINE metadata if a bounding polygon is available         */
852
    /* -------------------------------------------------------------------- */
853
59
    poDS->ScanForCutline();
854
855
    /* -------------------------------------------------------------------- */
856
    /*      Initialize any PAM information.                                 */
857
    /* -------------------------------------------------------------------- */
858
59
    poDS->SetDescription(poOpenInfo->pszFilename);
859
59
    poDS->TryLoadXML();
860
861
59
    poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
862
863
59
    return poDS;
864
747
}
865
866
/************************************************************************/
867
/*                            GetGCPCount()                             */
868
/************************************************************************/
869
870
int BSBDataset::GetGCPCount()
871
872
59
{
873
59
    return nGCPCount;
874
59
}
875
876
/************************************************************************/
877
/*                          GetGCPSpatialRef()                          */
878
/************************************************************************/
879
880
const OGRSpatialReference *BSBDataset::GetGCPSpatialRef() const
881
882
59
{
883
59
    return m_oGCPSRS.IsEmpty() ? nullptr : &m_oGCPSRS;
884
59
}
885
886
/************************************************************************/
887
/*                               GetGCP()                               */
888
/************************************************************************/
889
890
const GDAL_GCP *BSBDataset::GetGCPs()
891
892
59
{
893
59
    return pasGCPList;
894
59
}
895
896
#ifdef BSB_CREATE
897
898
/************************************************************************/
899
/*                             BSBIsSRSOK()                             */
900
/************************************************************************/
901
902
static int BSBIsSRSOK(const char *pszWKT)
903
{
904
    bool bOK = false;
905
    OGRSpatialReference oSRS, oSRS_WGS84, oSRS_NAD83;
906
907
    if (pszWKT != NULL && pszWKT[0] != '\0')
908
    {
909
        oSRS.importFromWkt(pszWKT);
910
911
        oSRS_WGS84.SetWellKnownGeogCS("WGS84");
912
        oSRS_NAD83.SetWellKnownGeogCS("NAD83");
913
        if ((oSRS.IsSameGeogCS(&oSRS_WGS84) ||
914
             oSRS.IsSameGeogCS(&oSRS_NAD83)) &&
915
            oSRS.IsGeographic() && oSRS.GetPrimeMeridian() == 0.0)
916
        {
917
            bOK = true;
918
        }
919
    }
920
921
    if (!bOK)
922
    {
923
        CPLError(CE_Warning, CPLE_NotSupported,
924
                 "BSB only supports WGS84 or NAD83 geographic projections.\n");
925
    }
926
927
    return bOK;
928
}
929
930
/************************************************************************/
931
/*                           BSBCreateCopy()                            */
932
/************************************************************************/
933
934
static GDALDataset *BSBCreateCopy(const char *pszFilename, GDALDataset *poSrcDS,
935
                                  int bStrict, char ** /*papszOptions*/,
936
                                  GDALProgressFunc /*pfnProgress*/,
937
                                  void * /*pProgressData*/)
938
939
{
940
    /* -------------------------------------------------------------------- */
941
    /*      Some some rudimentary checks                                    */
942
    /* -------------------------------------------------------------------- */
943
    const int nBands = poSrcDS->GetRasterCount();
944
    if (nBands != 1)
945
    {
946
        CPLError(CE_Failure, CPLE_NotSupported,
947
                 "BSB driver only supports one band images.\n");
948
949
        return NULL;
950
    }
951
952
    if (poSrcDS->GetRasterBand(1)->GetRasterDataType() != GDT_Byte && bStrict)
953
    {
954
        CPLError(CE_Failure, CPLE_NotSupported,
955
                 "BSB driver doesn't support data type %s. "
956
                 "Only eight bit bands supported.\n",
957
                 GDALGetDataTypeName(
958
                     poSrcDS->GetRasterBand(1)->GetRasterDataType()));
959
960
        return NULL;
961
    }
962
963
    const int nXSize = poSrcDS->GetRasterXSize();
964
    const int nYSize = poSrcDS->GetRasterYSize();
965
966
    /* -------------------------------------------------------------------- */
967
    /*      Open the output file.                                           */
968
    /* -------------------------------------------------------------------- */
969
    BSBInfo *psBSB = BSBCreate(pszFilename, 0, 200, nXSize, nYSize);
970
    if (psBSB == NULL)
971
        return NULL;
972
973
    /* -------------------------------------------------------------------- */
974
    /*      Prepare initial color table.colortable.                         */
975
    /* -------------------------------------------------------------------- */
976
    GDALRasterBand *poBand = poSrcDS->GetRasterBand(1);
977
    unsigned char abyPCT[771];
978
    int nPCTSize;
979
    int anRemap[256];
980
981
    abyPCT[0] = 0;
982
    abyPCT[1] = 0;
983
    abyPCT[2] = 0;
984
985
    if (poBand->GetColorTable() == NULL)
986
    {
987
        /* map greyscale down to 63 grey levels. */
988
        for (int iColor = 0; iColor < 256; iColor++)
989
        {
990
            int nOutValue = (int)(iColor / 4.1) + 1;
991
992
            anRemap[iColor] = nOutValue;
993
            abyPCT[nOutValue * 3 + 0] = (unsigned char)iColor;
994
            abyPCT[nOutValue * 3 + 1] = (unsigned char)iColor;
995
            abyPCT[nOutValue * 3 + 2] = (unsigned char)iColor;
996
        }
997
        nPCTSize = 64;
998
    }
999
    else
1000
    {
1001
        GDALColorTable *poCT = poBand->GetColorTable();
1002
        int nColorTableSize = poCT->GetColorEntryCount();
1003
        if (nColorTableSize > 255)
1004
            nColorTableSize = 255;
1005
1006
        for (int iColor = 0; iColor < nColorTableSize; iColor++)
1007
        {
1008
            GDALColorEntry sEntry;
1009
1010
            poCT->GetColorEntryAsRGB(iColor, &sEntry);
1011
1012
            anRemap[iColor] = iColor + 1;
1013
            abyPCT[(iColor + 1) * 3 + 0] = (unsigned char)sEntry.c1;
1014
            abyPCT[(iColor + 1) * 3 + 1] = (unsigned char)sEntry.c2;
1015
            abyPCT[(iColor + 1) * 3 + 2] = (unsigned char)sEntry.c3;
1016
        }
1017
1018
        nPCTSize = nColorTableSize + 1;
1019
1020
        // Add entries for pixel values which apparently will not occur.
1021
        for (int iColor = nPCTSize; iColor < 256; iColor++)
1022
            anRemap[iColor] = 1;
1023
    }
1024
1025
    /* -------------------------------------------------------------------- */
1026
    /*      Boil out all duplicate entries.                                 */
1027
    /* -------------------------------------------------------------------- */
1028
    for (int i = 1; i < nPCTSize - 1; i++)
1029
    {
1030
        for (int j = i + 1; j < nPCTSize; j++)
1031
        {
1032
            if (abyPCT[i * 3 + 0] == abyPCT[j * 3 + 0] &&
1033
                abyPCT[i * 3 + 1] == abyPCT[j * 3 + 1] &&
1034
                abyPCT[i * 3 + 2] == abyPCT[j * 3 + 2])
1035
            {
1036
                nPCTSize--;
1037
                abyPCT[j * 3 + 0] = abyPCT[nPCTSize * 3 + 0];
1038
                abyPCT[j * 3 + 1] = abyPCT[nPCTSize * 3 + 1];
1039
                abyPCT[j * 3 + 2] = abyPCT[nPCTSize * 3 + 2];
1040
1041
                for (int k = 0; k < 256; k++)
1042
                {
1043
                    // merge matching entries.
1044
                    if (anRemap[k] == j)
1045
                        anRemap[k] = i;
1046
1047
                    // shift the last PCT entry into the new hole.
1048
                    if (anRemap[k] == nPCTSize)
1049
                        anRemap[k] = j;
1050
                }
1051
            }
1052
        }
1053
    }
1054
1055
    /* -------------------------------------------------------------------- */
1056
    /*      Boil out all duplicate entries.                                 */
1057
    /* -------------------------------------------------------------------- */
1058
    if (nPCTSize > 128)
1059
    {
1060
        CPLError(CE_Warning, CPLE_AppDefined,
1061
                 "Having to merge color table entries to reduce %d real\n"
1062
                 "color table entries down to 127 values.",
1063
                 nPCTSize);
1064
    }
1065
1066
    while (nPCTSize > 128)
1067
    {
1068
        int nBestRange = 768;
1069
        int iBestMatch1 = -1;
1070
        int iBestMatch2 = -1;
1071
1072
        // Find the closest pair of color table entries.
1073
1074
        for (int i = 1; i < nPCTSize - 1; i++)
1075
        {
1076
            for (int j = i + 1; j < nPCTSize; j++)
1077
            {
1078
                int nRange = std::abs(abyPCT[i * 3 + 0] - abyPCT[j * 3 + 0]) +
1079
                             std::abs(abyPCT[i * 3 + 1] - abyPCT[j * 3 + 1]) +
1080
                             std::abs(abyPCT[i * 3 + 2] - abyPCT[j * 3 + 2]);
1081
1082
                if (nRange < nBestRange)
1083
                {
1084
                    iBestMatch1 = i;
1085
                    iBestMatch2 = j;
1086
                    nBestRange = nRange;
1087
                }
1088
            }
1089
        }
1090
1091
        // Merge the second entry into the first.
1092
        nPCTSize--;
1093
        abyPCT[iBestMatch2 * 3 + 0] = abyPCT[nPCTSize * 3 + 0];
1094
        abyPCT[iBestMatch2 * 3 + 1] = abyPCT[nPCTSize * 3 + 1];
1095
        abyPCT[iBestMatch2 * 3 + 2] = abyPCT[nPCTSize * 3 + 2];
1096
1097
        for (int i = 0; i < 256; i++)
1098
        {
1099
            // merge matching entries.
1100
            if (anRemap[i] == iBestMatch2)
1101
                anRemap[i] = iBestMatch1;
1102
1103
            // shift the last PCT entry into the new hole.
1104
            if (anRemap[i] == nPCTSize)
1105
                anRemap[i] = iBestMatch2;
1106
        }
1107
    }
1108
1109
    /* -------------------------------------------------------------------- */
1110
    /*      Write the PCT.                                                  */
1111
    /* -------------------------------------------------------------------- */
1112
    if (!BSBWritePCT(psBSB, nPCTSize, abyPCT))
1113
    {
1114
        BSBClose(psBSB);
1115
        return NULL;
1116
    }
1117
1118
    /* -------------------------------------------------------------------- */
1119
    /*      Write the GCPs.                                                 */
1120
    /* -------------------------------------------------------------------- */
1121
    double adfGeoTransform[6];
1122
    int nGCPCount = poSrcDS->GetGCPCount();
1123
    if (nGCPCount)
1124
    {
1125
        const char *pszGCPProjection = poSrcDS->GetGCPProjection();
1126
        if (BSBIsSRSOK(pszGCPProjection))
1127
        {
1128
            const GDAL_GCP *pasGCPList = poSrcDS->GetGCPs();
1129
            for (int i = 0; i < nGCPCount; i++)
1130
            {
1131
                VSIFPrintfL(psBSB->fp, "REF/%d,%f,%f,%f,%f\n", i + 1,
1132
                            pasGCPList[i].dfGCPPixel, pasGCPList[i].dfGCPLine,
1133
                            pasGCPList[i].dfGCPY, pasGCPList[i].dfGCPX);
1134
            }
1135
        }
1136
    }
1137
    else if (poSrcDS->GetGeoTransform(adfGeoTransform) == CE_None)
1138
    {
1139
        const char *pszProjection = poSrcDS->GetProjectionRef();
1140
        if (BSBIsSRSOK(pszProjection))
1141
        {
1142
            VSIFPrintfL(psBSB->fp, "REF/%d,%d,%d,%f,%f\n", 1, 0, 0,
1143
                        adfGeoTransform[3] + 0 * adfGeoTransform[4] +
1144
                            0 * adfGeoTransform[5],
1145
                        adfGeoTransform[0] + 0 * adfGeoTransform[1] +
1146
                            0 * adfGeoTransform[2]);
1147
            VSIFPrintfL(psBSB->fp, "REF/%d,%d,%d,%f,%f\n", 2, nXSize, 0,
1148
                        adfGeoTransform[3] + nXSize * adfGeoTransform[4] +
1149
                            0 * adfGeoTransform[5],
1150
                        adfGeoTransform[0] + nXSize * adfGeoTransform[1] +
1151
                            0 * adfGeoTransform[2]);
1152
            VSIFPrintfL(psBSB->fp, "REF/%d,%d,%d,%f,%f\n", 3, nXSize, nYSize,
1153
                        adfGeoTransform[3] + nXSize * adfGeoTransform[4] +
1154
                            nYSize * adfGeoTransform[5],
1155
                        adfGeoTransform[0] + nXSize * adfGeoTransform[1] +
1156
                            nYSize * adfGeoTransform[2]);
1157
            VSIFPrintfL(psBSB->fp, "REF/%d,%d,%d,%f,%f\n", 4, 0, nYSize,
1158
                        adfGeoTransform[3] + 0 * adfGeoTransform[4] +
1159
                            nYSize * adfGeoTransform[5],
1160
                        adfGeoTransform[0] + 0 * adfGeoTransform[1] +
1161
                            nYSize * adfGeoTransform[2]);
1162
        }
1163
    }
1164
1165
    /* -------------------------------------------------------------------- */
1166
    /*      Loop over image, copying image data.                            */
1167
    /* -------------------------------------------------------------------- */
1168
    CPLErr eErr = CE_None;
1169
1170
    GByte *pabyScanline = (GByte *)CPLMalloc(nXSize);
1171
1172
    for (int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++)
1173
    {
1174
        eErr =
1175
            poBand->RasterIO(GF_Read, 0, iLine, nXSize, 1, pabyScanline, nXSize,
1176
                             1, GDT_Byte, nBands, nBands * nXSize, NULL);
1177
        if (eErr == CE_None)
1178
        {
1179
            for (int i = 0; i < nXSize; i++)
1180
                pabyScanline[i] = (GByte)anRemap[pabyScanline[i]];
1181
1182
            if (!BSBWriteScanline(psBSB, pabyScanline))
1183
                eErr = CE_Failure;
1184
        }
1185
    }
1186
1187
    CPLFree(pabyScanline);
1188
1189
    /* -------------------------------------------------------------------- */
1190
    /*      cleanup                                                         */
1191
    /* -------------------------------------------------------------------- */
1192
    BSBClose(psBSB);
1193
1194
    if (eErr != CE_None)
1195
    {
1196
        VSIUnlink(pszFilename);
1197
        return NULL;
1198
    }
1199
1200
    return (GDALDataset *)GDALOpen(pszFilename, GA_ReadOnly);
1201
}
1202
#endif
1203
1204
/************************************************************************/
1205
/*                        GDALRegister_BSB()                            */
1206
/************************************************************************/
1207
1208
void GDALRegister_BSB()
1209
1210
24
{
1211
24
    if (GDALGetDriverByName("BSB") != nullptr)
1212
0
        return;
1213
1214
24
    GDALDriver *poDriver = new GDALDriver();
1215
1216
24
    poDriver->SetDescription("BSB");
1217
24
    poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1218
24
    poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Maptech BSB Nautical Charts");
1219
24
    poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/bsb.html");
1220
24
    poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1221
24
    poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "kap");
1222
#ifdef BSB_CREATE
1223
    poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Byte");
1224
#endif
1225
24
    poDriver->pfnOpen = BSBDataset::Open;
1226
24
    poDriver->pfnIdentify = BSBDataset::Identify;
1227
#ifdef BSB_CREATE
1228
    poDriver->pfnCreateCopy = BSBCreateCopy;
1229
#endif
1230
1231
24
    GetGDALDriverManager()->RegisterDriver(poDriver);
1232
24
}