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