Coverage Report

Created: 2025-06-09 07:07

/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
0
{
93
0
    poDS = poDSIn;
94
0
    nBand = 1;
95
96
0
    eDataType = GDT_Byte;
97
98
0
    nBlockXSize = poDS->GetRasterXSize();
99
0
    nBlockYSize = 1;
100
101
    // Note that the first color table entry is dropped, everything is
102
    // shifted down.
103
0
    for (int i = 0; i < poDSIn->psInfo->nPCTSize - 1; i++)
104
0
    {
105
0
        GDALColorEntry oColor = {poDSIn->psInfo->pabyPCT[i * 3 + 0 + 3],
106
0
                                 poDSIn->psInfo->pabyPCT[i * 3 + 1 + 3],
107
0
                                 poDSIn->psInfo->pabyPCT[i * 3 + 2 + 3], 255};
108
109
0
        oCT.SetColorEntry(i, &oColor);
110
0
    }
111
0
}
112
113
/************************************************************************/
114
/*                             IReadBlock()                             */
115
/************************************************************************/
116
117
CPLErr BSBRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
118
                                 void *pImage)
119
0
{
120
0
    BSBDataset *poGDS = cpl::down_cast<BSBDataset *>(poDS);
121
0
    GByte *pabyScanline = (GByte *)pImage;
122
123
0
    if (BSBReadScanline(poGDS->psInfo, nBlockYOff, pabyScanline))
124
0
    {
125
0
        for (int i = 0; i < nBlockXSize; i++)
126
0
        {
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
0
            if (pabyScanline[i] > 0)
131
0
                pabyScanline[i] -= 1;
132
0
        }
133
134
0
        return CE_None;
135
0
    }
136
137
0
    return CE_Failure;
138
0
}
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
0
    : nGCPCount(0), pasGCPList(nullptr), bGeoTransformSet(FALSE),
172
0
      psInfo(nullptr)
173
0
{
174
0
    m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
175
0
    m_oGCPSRS.importFromWkt(SRS_WKT_WGS84_LAT_LONG);
176
177
0
    adfGeoTransform[0] = 0.0; /* X Origin (top left corner) */
178
0
    adfGeoTransform[1] = 1.0; /* X Pixel size */
179
0
    adfGeoTransform[2] = 0.0;
180
0
    adfGeoTransform[3] = 0.0; /* Y Origin (top left corner) */
181
0
    adfGeoTransform[4] = 0.0;
182
0
    adfGeoTransform[5] = 1.0; /* Y Pixel Size */
183
0
}
184
185
/************************************************************************/
186
/*                            ~BSBDataset()                             */
187
/************************************************************************/
188
189
BSBDataset::~BSBDataset()
190
191
0
{
192
0
    FlushCache(true);
193
194
0
    GDALDeinitGCPs(nGCPCount, pasGCPList);
195
0
    CPLFree(pasGCPList);
196
197
0
    if (psInfo != nullptr)
198
0
        BSBClose(psInfo);
199
0
}
200
201
/************************************************************************/
202
/*                          GetGeoTransform()                           */
203
/************************************************************************/
204
205
CPLErr BSBDataset::GetGeoTransform(double *padfTransform)
206
207
0
{
208
0
    memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
209
210
0
    if (bGeoTransformSet)
211
0
        return CE_None;
212
213
0
    return CE_Failure;
214
0
}
215
216
/************************************************************************/
217
/*                          GetSpatialRef()                             */
218
/************************************************************************/
219
220
const OGRSpatialReference *BSBDataset::GetSpatialRef() const
221
222
0
{
223
0
    if (bGeoTransformSet)
224
0
        return &m_oGCPSRS;
225
226
0
    return nullptr;
227
0
}
228
229
/************************************************************************/
230
/*                     GDALHeuristicDatelineWrap()                      */
231
/************************************************************************/
232
233
static void GDALHeuristicDatelineWrap(int nPointCount, double *padfX)
234
235
0
{
236
0
    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
0
    double dfX_PM_Min = 0.0;
246
0
    double dfX_PM_Max = 0.0;
247
0
    double dfX_Dateline_Min = 0.0;
248
0
    double dfX_Dateline_Max = 0.0;
249
250
0
    for (int i = 0; i < nPointCount; i++)
251
0
    {
252
0
        double dfX_PM = padfX[i];
253
0
        if (dfX_PM > 180)
254
0
            dfX_PM -= 360.0;
255
256
0
        double dfX_Dateline = padfX[i];
257
0
        if (dfX_Dateline < 0)
258
0
            dfX_Dateline += 360.0;
259
260
0
        if (i == 0)
261
0
        {
262
0
            dfX_PM_Min = dfX_PM;
263
0
            dfX_PM_Max = dfX_PM;
264
0
            dfX_Dateline_Min = dfX_Dateline;
265
0
            dfX_Dateline_Max = dfX_Dateline;
266
0
        }
267
0
        else
268
0
        {
269
0
            dfX_PM_Min = std::min(dfX_PM_Min, dfX_PM);
270
0
            dfX_PM_Max = std::max(dfX_PM_Max, dfX_PM);
271
0
            dfX_Dateline_Min = std::min(dfX_Dateline_Min, dfX_Dateline);
272
0
            dfX_Dateline_Max = std::max(dfX_Dateline_Max, dfX_Dateline);
273
0
        }
274
0
    }
275
276
    /* -------------------------------------------------------------------- */
277
    /*      Do nothing if the range is always fairly small - no apparent    */
278
    /*      wrapping issues.                                                */
279
    /* -------------------------------------------------------------------- */
280
0
    if ((dfX_PM_Max - dfX_PM_Min) < 270.0 &&
281
0
        (dfX_Dateline_Max - dfX_Dateline_Min) < 270.0)
282
0
        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
0
    if ((dfX_PM_Max - dfX_PM_Min) > 270.0 &&
289
0
        (dfX_Dateline_Max - dfX_Dateline_Min) > 270.0)
290
0
        return;
291
292
    /* -------------------------------------------------------------------- */
293
    /*      Pick which way to transform things.                             */
294
    /* -------------------------------------------------------------------- */
295
0
    bool bUsePMWrap;
296
297
0
    if ((dfX_PM_Max - dfX_PM_Min) > 270.0 &&
298
0
        (dfX_Dateline_Max - dfX_Dateline_Min) < 270.0)
299
0
        bUsePMWrap = false;
300
0
    else
301
0
        bUsePMWrap = true;
302
303
    /* -------------------------------------------------------------------- */
304
    /*      Apply rewrapping.                                               */
305
    /* -------------------------------------------------------------------- */
306
0
    for (int i = 0; i < nPointCount; i++)
307
0
    {
308
0
        if (bUsePMWrap)
309
0
        {
310
0
            if (padfX[i] > 180)
311
0
                padfX[i] -= 360.0;
312
0
        }
313
0
        else
314
0
        {
315
0
            if (padfX[i] < 0)
316
0
                padfX[i] += 360.0;
317
0
        }
318
0
    }
319
0
}
320
321
/************************************************************************/
322
/*                   GDALHeuristicDatelineWrapGCPs()                    */
323
/************************************************************************/
324
325
static void GDALHeuristicDatelineWrapGCPs(int nPointCount, GDAL_GCP *pasGCPList)
326
0
{
327
0
    std::vector<double> oadfX;
328
329
0
    oadfX.resize(nPointCount);
330
0
    for (int i = 0; i < nPointCount; i++)
331
0
        oadfX[i] = pasGCPList[i].dfGCPX;
332
333
0
    GDALHeuristicDatelineWrap(nPointCount, &(oadfX[0]));
334
335
0
    for (int i = 0; i < nPointCount; i++)
336
0
        pasGCPList[i].dfGCPX = oadfX[i];
337
0
}
338
339
/************************************************************************/
340
/*                            ScanForGCPs()                             */
341
/************************************************************************/
342
343
void BSBDataset::ScanForGCPs(bool isNos, const char *pszFilename)
344
345
0
{
346
    /* -------------------------------------------------------------------- */
347
    /*      Collect GCPs as appropriate to source.                          */
348
    /* -------------------------------------------------------------------- */
349
0
    nGCPCount = 0;
350
351
0
    if (isNos)
352
0
    {
353
0
        ScanForGCPsNos(pszFilename);
354
0
    }
355
0
    else
356
0
    {
357
0
        ScanForGCPsBSB();
358
0
    }
359
360
    /* -------------------------------------------------------------------- */
361
    /*      Apply heuristics to re-wrap GCPs to maintain continuity        */
362
    /*      over the international dateline.                                */
363
    /* -------------------------------------------------------------------- */
364
0
    if (nGCPCount > 1)
365
0
        GDALHeuristicDatelineWrapGCPs(nGCPCount, pasGCPList);
366
367
    /* -------------------------------------------------------------------- */
368
    /*      Collect coordinate system related parameters from header.       */
369
    /* -------------------------------------------------------------------- */
370
0
    const char *pszKNP = nullptr;
371
0
    const char *pszKNQ = nullptr;
372
373
0
    for (int i = 0; psInfo->papszHeader[i] != nullptr; i++)
374
0
    {
375
0
        if (STARTS_WITH_CI(psInfo->papszHeader[i], "KNP/"))
376
0
        {
377
0
            pszKNP = psInfo->papszHeader[i];
378
0
            SetMetadataItem("BSB_KNP", pszKNP + 4);
379
0
        }
380
0
        if (STARTS_WITH_CI(psInfo->papszHeader[i], "KNQ/"))
381
0
        {
382
0
            pszKNQ = psInfo->papszHeader[i];
383
0
            SetMetadataItem("BSB_KNQ", pszKNQ + 4);
384
0
        }
385
0
    }
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
0
    CPLString osUnderlyingSRS;
393
0
    if (pszKNP != nullptr)
394
0
    {
395
0
        const char *pszPR = strstr(pszKNP, "PR=");
396
0
        const char *pszGD = strstr(pszKNP, "GD=");
397
0
        const char *pszGEOGCS = SRS_WKT_WGS84_LAT_LONG;
398
0
        CPLString osPP;
399
400
        // Capture the PP string.
401
0
        const char *pszValue = strstr(pszKNP, "PP=");
402
0
        const char *pszEnd = pszValue ? strstr(pszValue, ",") : nullptr;
403
0
        if (pszValue && pszEnd)
404
0
            osPP.assign(pszValue + 3, pszEnd - pszValue - 3);
405
406
        // Look at the datum
407
0
        if (pszGD == nullptr)
408
0
        {
409
            /* no match. We'll default to EPSG:4326 */
410
0
        }
411
0
        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
0
        if (pszPR == nullptr)
425
0
        {
426
            /* no match */
427
0
        }
428
0
        else if (STARTS_WITH_CI(pszPR, "PR=MERCATOR") && nGCPCount > 0)
429
0
        {
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
0
            osUnderlyingSRS.Printf(
435
0
                "PROJCS[\"Global "
436
0
                "Mercator\",%s,PROJECTION[\"Mercator_2SP\"],PARAMETER["
437
0
                "\"standard_parallel_1\",0],PARAMETER[\"latitude_of_origin\",0]"
438
0
                ",PARAMETER[\"central_meridian\",%d],PARAMETER[\"false_"
439
0
                "easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]"
440
0
                "]",
441
0
                pszGEOGCS, (int)pasGCPList[0].dfGCPX);
442
0
        }
443
444
0
        else if (STARTS_WITH_CI(pszPR, "PR=TRANSVERSE MERCATOR") &&
445
0
                 !osPP.empty())
446
0
        {
447
448
0
            osUnderlyingSRS.Printf(
449
0
                "PROJCS[\"unnamed\",%s,PROJECTION[\"Transverse_Mercator\"],"
450
0
                "PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_"
451
0
                "meridian\",%s],PARAMETER[\"scale_factor\",1],PARAMETER["
452
0
                "\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT["
453
0
                "\"Meter\",1]]",
454
0
                pszGEOGCS, osPP.c_str());
455
0
        }
456
457
0
        else if (STARTS_WITH_CI(pszPR, "PR=UNIVERSAL TRANSVERSE MERCATOR") &&
458
0
                 !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
0
        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
0
        else if (STARTS_WITH_CI(pszPR, "PR=LAMBERT CONFORMAL CONIC") &&
482
0
                 !osPP.empty() && pszKNQ != nullptr)
483
0
        {
484
0
            CPLString osP2, osP3;
485
486
            // Capture the KNQ/P2 string.
487
0
            pszValue = strstr(pszKNQ, "P2=");
488
0
            if (pszValue)
489
0
                pszEnd = strstr(pszValue, ",");
490
0
            if (pszValue && pszEnd)
491
0
                osP2.assign(pszValue + 3, pszEnd - pszValue - 3);
492
493
            // Capture the KNQ/P3 string.
494
0
            pszValue = strstr(pszKNQ, "P3=");
495
0
            if (pszValue)
496
0
                pszEnd = strstr(pszValue, ",");
497
0
            if (pszValue)
498
0
            {
499
0
                if (pszEnd)
500
0
                    osP3.assign(pszValue + 3, pszEnd - pszValue - 3);
501
0
                else
502
0
                    osP3.assign(pszValue + 3);
503
0
            }
504
505
0
            if (!osP2.empty() && !osP3.empty())
506
0
                osUnderlyingSRS.Printf(
507
0
                    "PROJCS[\"unnamed\",%s,PROJECTION[\"Lambert_Conformal_"
508
0
                    "Conic_2SP\"],PARAMETER[\"standard_parallel_1\",%s],"
509
0
                    "PARAMETER[\"standard_parallel_2\",%s],PARAMETER["
510
0
                    "\"latitude_of_origin\",0.0],PARAMETER[\"central_"
511
0
                    "meridian\",%s],PARAMETER[\"false_easting\",0.0],PARAMETER["
512
0
                    "\"false_northing\",0.0],UNIT[\"Meter\",1]]",
513
0
                    pszGEOGCS, osP2.c_str(), osP3.c_str(), osPP.c_str());
514
0
        }
515
0
    }
516
517
    /* -------------------------------------------------------------------- */
518
    /*      If we got an alternate underlying coordinate system, try        */
519
    /*      converting the GCPs to that coordinate system.                  */
520
    /* -------------------------------------------------------------------- */
521
0
    if (osUnderlyingSRS.length() > 0)
522
0
    {
523
0
        OGRSpatialReference oGeog_SRS, oProjected_SRS;
524
525
0
        oProjected_SRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
526
0
        oProjected_SRS.SetFromUserInput(osUnderlyingSRS);
527
0
        oGeog_SRS.CopyGeogCSFrom(&oProjected_SRS);
528
0
        oGeog_SRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
529
530
0
        OGRCoordinateTransformation *poCT =
531
0
            OGRCreateCoordinateTransformation(&oGeog_SRS, &oProjected_SRS);
532
0
        if (poCT != nullptr)
533
0
        {
534
0
            for (int i = 0; i < nGCPCount; i++)
535
0
            {
536
0
                poCT->Transform(1, &(pasGCPList[i].dfGCPX),
537
0
                                &(pasGCPList[i].dfGCPY),
538
0
                                &(pasGCPList[i].dfGCPZ));
539
0
            }
540
541
0
            m_oGCPSRS.importFromWkt(osUnderlyingSRS.c_str());
542
543
0
            delete poCT;
544
0
        }
545
0
        else
546
0
            CPLErrorReset();
547
0
    }
548
549
    /* -------------------------------------------------------------------- */
550
    /*      Attempt to prepare a geotransform from the GCPs.                */
551
    /* -------------------------------------------------------------------- */
552
0
    if (GDALGCPsToGeoTransform(nGCPCount, pasGCPList, adfGeoTransform, FALSE))
553
0
    {
554
0
        bGeoTransformSet = TRUE;
555
0
    }
556
0
}
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
0
{
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
0
    int fileGCPCount = 0;
650
651
    // Count the GCPs (reference points) in psInfo->papszHeader
652
0
    for (int i = 0; psInfo->papszHeader[i] != nullptr; i++)
653
0
        if (STARTS_WITH_CI(psInfo->papszHeader[i], "REF/"))
654
0
            fileGCPCount++;
655
656
    // Memory has not been allocated to fileGCPCount yet
657
0
    pasGCPList = (GDAL_GCP *)CPLCalloc(sizeof(GDAL_GCP), fileGCPCount + 1);
658
659
0
    for (int i = 0; psInfo->papszHeader[i] != nullptr; i++)
660
0
    {
661
662
0
        if (!STARTS_WITH_CI(psInfo->papszHeader[i], "REF/"))
663
0
            continue;
664
665
0
        char **papszTokens = CSLTokenizeStringComplex(
666
0
            psInfo->papszHeader[i] + 4, ",", FALSE, FALSE);
667
668
0
        if (CSLCount(papszTokens) > 4)
669
0
        {
670
0
            GDALInitGCPs(1, pasGCPList + nGCPCount);
671
672
0
            pasGCPList[nGCPCount].dfGCPX = CPLAtof(papszTokens[4]);
673
0
            pasGCPList[nGCPCount].dfGCPY = CPLAtof(papszTokens[3]);
674
0
            pasGCPList[nGCPCount].dfGCPPixel = CPLAtof(papszTokens[1]);
675
0
            pasGCPList[nGCPCount].dfGCPLine = CPLAtof(papszTokens[2]);
676
677
0
            CPLFree(pasGCPList[nGCPCount].pszId);
678
0
            if (CSLCount(papszTokens) > 5)
679
0
            {
680
0
                pasGCPList[nGCPCount].pszId = CPLStrdup(papszTokens[5]);
681
0
            }
682
0
            else
683
0
            {
684
0
                char szName[50];
685
0
                snprintf(szName, sizeof(szName), "GCP_%d", nGCPCount + 1);
686
0
                pasGCPList[nGCPCount].pszId = CPLStrdup(szName);
687
0
            }
688
689
0
            nGCPCount++;
690
0
        }
691
0
        CSLDestroy(papszTokens);
692
0
    }
693
0
}
694
695
/************************************************************************/
696
/*                            ScanForCutline()                          */
697
/************************************************************************/
698
699
void BSBDataset::ScanForCutline()
700
0
{
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
0
    std::string wkt;
712
0
    for (int i = 0; psInfo->papszHeader[i] != nullptr; i++)
713
0
    {
714
0
        if (!STARTS_WITH_CI(psInfo->papszHeader[i], "PLY/"))
715
0
            continue;
716
717
0
        const CPLStringList aosTokens(
718
0
            CSLTokenizeString2(psInfo->papszHeader[i] + 4, ",", 0));
719
720
0
        if (aosTokens.size() >= 3)
721
0
        {
722
0
            if (wkt.empty())
723
0
                wkt = "POLYGON ((";
724
0
            else
725
0
                wkt += ',';
726
0
            wkt += aosTokens[2];
727
0
            wkt += ' ';
728
0
            wkt += aosTokens[1];
729
0
        }
730
0
    }
731
732
0
    if (!wkt.empty())
733
0
    {
734
0
        wkt += "))";
735
0
        SetMetadataItem("BSB_CUTLINE", wkt.c_str());
736
0
    }
737
0
}
738
739
/************************************************************************/
740
/*                          IdentifyInternal()                          */
741
/************************************************************************/
742
743
int BSBDataset::IdentifyInternal(GDALOpenInfo *poOpenInfo, bool &isNosOut)
744
745
16.3k
{
746
    /* -------------------------------------------------------------------- */
747
    /*      Check for BSB/ keyword.                                         */
748
    /* -------------------------------------------------------------------- */
749
16.3k
    isNosOut = false;
750
751
16.3k
    if (poOpenInfo->nHeaderBytes < 1000)
752
3.90k
        return FALSE;
753
754
12.3k
    int i = 0;
755
14.6M
    for (; i < poOpenInfo->nHeaderBytes - 4; i++)
756
14.5M
    {
757
14.5M
        if (poOpenInfo->pabyHeader[i + 0] == 'B' &&
758
14.5M
            poOpenInfo->pabyHeader[i + 1] == 'S' &&
759
14.5M
            poOpenInfo->pabyHeader[i + 2] == 'B' &&
760
14.5M
            poOpenInfo->pabyHeader[i + 3] == '/')
761
0
            break;
762
14.5M
        if (poOpenInfo->pabyHeader[i + 0] == 'N' &&
763
14.5M
            poOpenInfo->pabyHeader[i + 1] == 'O' &&
764
14.5M
            poOpenInfo->pabyHeader[i + 2] == 'S' &&
765
14.5M
            poOpenInfo->pabyHeader[i + 3] == '/')
766
0
        {
767
0
            isNosOut = true;
768
0
            break;
769
0
        }
770
14.5M
        if (poOpenInfo->pabyHeader[i + 0] == 'W' &&
771
14.5M
            poOpenInfo->pabyHeader[i + 1] == 'X' &&
772
14.5M
            poOpenInfo->pabyHeader[i + 2] == '\\' &&
773
14.5M
            poOpenInfo->pabyHeader[i + 3] == '8')
774
0
            break;
775
14.5M
    }
776
777
12.3k
    if (i == poOpenInfo->nHeaderBytes - 4)
778
12.3k
        return FALSE;
779
780
    /* Additional test to avoid false positive. See #2881 */
781
0
    const char *pszHeader =
782
0
        reinterpret_cast<const char *>(poOpenInfo->pabyHeader);
783
0
    const char *pszShiftedHeader = pszHeader + i;
784
0
    const char *pszRA = strstr(pszShiftedHeader, "RA=");
785
0
    if (pszRA == nullptr) /* This may be a NO1 file */
786
0
        pszRA = strstr(pszShiftedHeader, "[JF");
787
0
    if (pszRA == nullptr)
788
0
        return FALSE;
789
0
    if (pszRA - pszShiftedHeader > 100 && !strstr(pszHeader, "VER/") &&
790
0
        !strstr(pszHeader, "KNP/") && !strstr(pszHeader, "KNQ/") &&
791
0
        !strstr(pszHeader, "RGB/"))
792
0
        return FALSE;
793
794
0
    return TRUE;
795
0
}
796
797
/************************************************************************/
798
/*                              Identify()                              */
799
/************************************************************************/
800
801
int BSBDataset::Identify(GDALOpenInfo *poOpenInfo)
802
803
16.3k
{
804
16.3k
    bool isNos;
805
16.3k
    return IdentifyInternal(poOpenInfo, isNos);
806
16.3k
}
807
808
/************************************************************************/
809
/*                                Open()                                */
810
/************************************************************************/
811
812
GDALDataset *BSBDataset::Open(GDALOpenInfo *poOpenInfo)
813
814
0
{
815
0
    bool isNos = false;
816
0
    if (!IdentifyInternal(poOpenInfo, isNos))
817
0
        return nullptr;
818
819
0
    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
0
    BSBDataset *poDS = new BSBDataset();
829
830
    /* -------------------------------------------------------------------- */
831
    /*      Open the file.                                                  */
832
    /* -------------------------------------------------------------------- */
833
0
    poDS->psInfo = BSBOpen(poOpenInfo->pszFilename);
834
0
    if (poDS->psInfo == nullptr)
835
0
    {
836
0
        delete poDS;
837
0
        return nullptr;
838
0
    }
839
840
0
    poDS->nRasterXSize = poDS->psInfo->nXSize;
841
0
    poDS->nRasterYSize = poDS->psInfo->nYSize;
842
843
    /* -------------------------------------------------------------------- */
844
    /*      Create band information objects.                                */
845
    /* -------------------------------------------------------------------- */
846
0
    poDS->SetBand(1, new BSBRasterBand(poDS));
847
848
0
    poDS->ScanForGCPs(isNos, poOpenInfo->pszFilename);
849
850
    /* -------------------------------------------------------------------- */
851
    /*      Set CUTLINE metadata if a bounding polygon is available         */
852
    /* -------------------------------------------------------------------- */
853
0
    poDS->ScanForCutline();
854
855
    /* -------------------------------------------------------------------- */
856
    /*      Initialize any PAM information.                                 */
857
    /* -------------------------------------------------------------------- */
858
0
    poDS->SetDescription(poOpenInfo->pszFilename);
859
0
    poDS->TryLoadXML();
860
861
0
    poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
862
863
0
    return poDS;
864
0
}
865
866
/************************************************************************/
867
/*                            GetGCPCount()                             */
868
/************************************************************************/
869
870
int BSBDataset::GetGCPCount()
871
872
0
{
873
0
    return nGCPCount;
874
0
}
875
876
/************************************************************************/
877
/*                          GetGCPSpatialRef()                          */
878
/************************************************************************/
879
880
const OGRSpatialReference *BSBDataset::GetGCPSpatialRef() const
881
882
0
{
883
0
    return m_oGCPSRS.IsEmpty() ? nullptr : &m_oGCPSRS;
884
0
}
885
886
/************************************************************************/
887
/*                               GetGCP()                               */
888
/************************************************************************/
889
890
const GDAL_GCP *BSBDataset::GetGCPs()
891
892
0
{
893
0
    return pasGCPList;
894
0
}
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
2
{
1211
2
    if (GDALGetDriverByName("BSB") != nullptr)
1212
0
        return;
1213
1214
2
    GDALDriver *poDriver = new GDALDriver();
1215
1216
2
    poDriver->SetDescription("BSB");
1217
2
    poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1218
2
    poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Maptech BSB Nautical Charts");
1219
2
    poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/bsb.html");
1220
2
    poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1221
2
    poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "kap");
1222
#ifdef BSB_CREATE
1223
    poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Byte");
1224
#endif
1225
2
    poDriver->pfnOpen = BSBDataset::Open;
1226
2
    poDriver->pfnIdentify = BSBDataset::Identify;
1227
#ifdef BSB_CREATE
1228
    poDriver->pfnCreateCopy = BSBCreateCopy;
1229
#endif
1230
1231
2
    GetGDALDriverManager()->RegisterDriver(poDriver);
1232
2
}