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