1 : /******************************************************************************
2 : * $Id: bsbdataset.cpp 18181 2009-12-05 01:07:47Z warmerdam $
3 : *
4 : * Project: BSB Reader
5 : * Purpose: BSBDataset implementation for BSB format.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2001, Frank Warmerdam <warmerdam@pobox.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 "gdal_pam.h"
31 : #include "bsb_read.h"
32 : #include "cpl_string.h"
33 : #include "ogr_spatialref.h"
34 :
35 : CPL_CVSID("$Id: bsbdataset.cpp 18181 2009-12-05 01:07:47Z warmerdam $");
36 :
37 : CPL_C_START
38 : void GDALRegister_BSB(void);
39 : CPL_C_END
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 : 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 : public:
65 : BSBDataset();
66 : ~BSBDataset();
67 :
68 : BSBInfo *psInfo;
69 :
70 : static GDALDataset *Open( GDALOpenInfo * );
71 :
72 : virtual int GetGCPCount();
73 : virtual const char *GetGCPProjection();
74 : virtual const GDAL_GCP *GetGCPs();
75 :
76 : CPLErr GetGeoTransform( double * padfTransform );
77 : const char *GetProjectionRef();
78 : };
79 :
80 : /************************************************************************/
81 : /* ==================================================================== */
82 : /* BSBRasterBand */
83 : /* ==================================================================== */
84 : /************************************************************************/
85 :
86 : class BSBRasterBand : public GDALPamRasterBand
87 10 : {
88 : GDALColorTable oCT;
89 :
90 : public:
91 : BSBRasterBand( BSBDataset * );
92 :
93 : virtual CPLErr IReadBlock( int, int, void * );
94 : virtual GDALColorTable *GetColorTable();
95 : virtual GDALColorInterp GetColorInterpretation();
96 : };
97 :
98 :
99 : /************************************************************************/
100 : /* BSBRasterBand() */
101 : /************************************************************************/
102 :
103 5 : BSBRasterBand::BSBRasterBand( BSBDataset *poDS )
104 :
105 : {
106 5 : this->poDS = poDS;
107 5 : this->nBand = 1;
108 :
109 5 : eDataType = GDT_Byte;
110 :
111 5 : nBlockXSize = poDS->GetRasterXSize();
112 5 : nBlockYSize = 1;
113 :
114 : // Note that the first color table entry is dropped, everything is
115 : // shifted down.
116 640 : for( int i = 0; i < poDS->psInfo->nPCTSize-1; i++ )
117 : {
118 : GDALColorEntry oColor;
119 :
120 635 : oColor.c1 = poDS->psInfo->pabyPCT[i*3+0+3];
121 635 : oColor.c2 = poDS->psInfo->pabyPCT[i*3+1+3];
122 635 : oColor.c3 = poDS->psInfo->pabyPCT[i*3+2+3];
123 635 : oColor.c4 = 255;
124 :
125 635 : oCT.SetColorEntry( i, &oColor );
126 : }
127 5 : }
128 :
129 : /************************************************************************/
130 : /* IReadBlock() */
131 : /************************************************************************/
132 :
133 250 : CPLErr BSBRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
134 : void * pImage )
135 :
136 : {
137 250 : BSBDataset *poGDS = (BSBDataset *) poDS;
138 250 : GByte *pabyScanline = (GByte*) pImage;
139 :
140 250 : if( BSBReadScanline( poGDS->psInfo, nBlockYOff, pabyScanline ) )
141 : {
142 12648 : for( int i = 0; i < nBlockXSize; i++ )
143 12400 : pabyScanline[i] -= 1;
144 :
145 248 : return CE_None;
146 : }
147 : else
148 2 : return CE_Failure;
149 : }
150 :
151 : /************************************************************************/
152 : /* GetColorTable() */
153 : /************************************************************************/
154 :
155 0 : GDALColorTable *BSBRasterBand::GetColorTable()
156 :
157 : {
158 0 : return &oCT;
159 : }
160 :
161 : /************************************************************************/
162 : /* GetColorInterpretation() */
163 : /************************************************************************/
164 :
165 0 : GDALColorInterp BSBRasterBand::GetColorInterpretation()
166 :
167 : {
168 0 : return GCI_PaletteIndex;
169 : }
170 :
171 : /************************************************************************/
172 : /* ==================================================================== */
173 : /* BSBDataset */
174 : /* ==================================================================== */
175 : /************************************************************************/
176 :
177 : /************************************************************************/
178 : /* BSBDataset() */
179 : /************************************************************************/
180 :
181 5 : BSBDataset::BSBDataset()
182 :
183 : {
184 5 : psInfo = NULL;
185 :
186 5 : bGeoTransformSet = FALSE;
187 :
188 5 : nGCPCount = 0;
189 5 : pasGCPList = NULL;
190 5 : osGCPProjection = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",7030]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",6326]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",8901]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",9108]],AUTHORITY[\"EPSG\",4326]]";
191 :
192 5 : adfGeoTransform[0] = 0.0; /* X Origin (top left corner) */
193 5 : adfGeoTransform[1] = 1.0; /* X Pixel size */
194 5 : adfGeoTransform[2] = 0.0;
195 5 : adfGeoTransform[3] = 0.0; /* Y Origin (top left corner) */
196 5 : adfGeoTransform[4] = 0.0;
197 5 : adfGeoTransform[5] = 1.0; /* Y Pixel Size */
198 :
199 5 : }
200 :
201 : /************************************************************************/
202 : /* ~BSBDataset() */
203 : /************************************************************************/
204 :
205 10 : BSBDataset::~BSBDataset()
206 :
207 : {
208 5 : FlushCache();
209 :
210 5 : GDALDeinitGCPs( nGCPCount, pasGCPList );
211 5 : CPLFree( pasGCPList );
212 :
213 5 : if( psInfo != NULL )
214 5 : BSBClose( psInfo );
215 10 : }
216 :
217 :
218 : /************************************************************************/
219 : /* GetGeoTransform() */
220 : /************************************************************************/
221 :
222 0 : CPLErr BSBDataset::GetGeoTransform( double * padfTransform )
223 :
224 : {
225 0 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
226 :
227 0 : if( bGeoTransformSet )
228 0 : return CE_None;
229 : else
230 0 : return CE_Failure;
231 : }
232 :
233 : /************************************************************************/
234 : /* GetProjectionRef() */
235 : /************************************************************************/
236 :
237 0 : const char *BSBDataset::GetProjectionRef()
238 :
239 : {
240 0 : if( bGeoTransformSet )
241 0 : return osGCPProjection;
242 : else
243 0 : return "";
244 : }
245 :
246 : /************************************************************************/
247 : /* GDALHeuristicDatelineWrap() */
248 : /************************************************************************/
249 :
250 : static void
251 0 : GDALHeuristicDatelineWrap( int nPointCount, double *padfX )
252 :
253 : {
254 : int i;
255 : /* Following inits are useless but keep GCC happy */
256 0 : double dfX_PM_Min = 0, dfX_PM_Max = 0, dfX_Dateline_Min = 0, dfX_Dateline_Max = 0;
257 : int bUsePMWrap;
258 :
259 0 : if( nPointCount < 2 )
260 0 : return;
261 :
262 : /* -------------------------------------------------------------------- */
263 : /* Work out what the longitude range will be centering on the */
264 : /* prime meridian (-180 to 180) and centering on the dateline */
265 : /* (0 to 360). */
266 : /* -------------------------------------------------------------------- */
267 0 : for( i = 0; i < nPointCount; i++ )
268 : {
269 : double dfX_PM, dfX_Dateline;
270 :
271 0 : dfX_PM = padfX[i];
272 0 : if( dfX_PM > 180 )
273 0 : dfX_PM -= 360.0;
274 :
275 0 : dfX_Dateline = padfX[i];
276 0 : if( dfX_Dateline < 0 )
277 0 : dfX_Dateline += 360.0;
278 :
279 0 : if( i == 0 )
280 : {
281 0 : dfX_PM_Min = dfX_PM_Max = dfX_PM;
282 0 : dfX_Dateline_Min = dfX_Dateline_Max = dfX_Dateline;
283 : }
284 : else
285 : {
286 0 : dfX_PM_Min = MIN(dfX_PM_Min,dfX_PM);
287 0 : dfX_PM_Max = MAX(dfX_PM_Max,dfX_PM);
288 0 : dfX_Dateline_Min = MIN(dfX_Dateline_Min,dfX_Dateline);
289 0 : dfX_Dateline_Max = MAX(dfX_Dateline_Max,dfX_Dateline);
290 : }
291 : }
292 :
293 : /* -------------------------------------------------------------------- */
294 : /* Do nothing if the range is always fairly small - no apparent */
295 : /* wrapping issues. */
296 : /* -------------------------------------------------------------------- */
297 0 : if( (dfX_PM_Max - dfX_PM_Min) < 270.0
298 : && (dfX_Dateline_Max - dfX_Dateline_Min) < 270.0 )
299 0 : return;
300 :
301 : /* -------------------------------------------------------------------- */
302 : /* Do nothing if both appproach have a wide range - best not to */
303 : /* fiddle if we aren't sure we are improving things. */
304 : /* -------------------------------------------------------------------- */
305 0 : if( (dfX_PM_Max - dfX_PM_Min) > 270.0
306 : && (dfX_Dateline_Max - dfX_Dateline_Min) > 270.0 )
307 0 : return;
308 :
309 : /* -------------------------------------------------------------------- */
310 : /* Pick which way to transform things. */
311 : /* -------------------------------------------------------------------- */
312 0 : if( (dfX_PM_Max - dfX_PM_Min) > 270.0
313 : && (dfX_Dateline_Max - dfX_Dateline_Min) < 270.0 )
314 0 : bUsePMWrap = FALSE;
315 : else
316 0 : bUsePMWrap = TRUE;
317 :
318 :
319 : /* -------------------------------------------------------------------- */
320 : /* Apply rewrapping. */
321 : /* -------------------------------------------------------------------- */
322 0 : for( i = 0; i < nPointCount; i++ )
323 : {
324 0 : if( bUsePMWrap )
325 : {
326 0 : if( padfX[i] > 180 )
327 0 : padfX[i] -= 360.0;
328 : }
329 : else
330 : {
331 0 : if( padfX[i] < 0 )
332 0 : padfX[i] += 360.0;
333 : }
334 : }
335 : }
336 :
337 : /************************************************************************/
338 : /* GDALHeuristicDatelineWrapGCPs() */
339 : /************************************************************************/
340 :
341 : static void
342 0 : GDALHeuristicDatelineWrapGCPs( int nPointCount, GDAL_GCP *pasGCPList )
343 : {
344 0 : std::vector<double> oadfX;
345 : int i;
346 :
347 0 : oadfX.resize( nPointCount );
348 0 : for( i = 0; i < nPointCount; i++ )
349 0 : oadfX[i] = pasGCPList[i].dfGCPX;
350 :
351 0 : GDALHeuristicDatelineWrap( nPointCount, &(oadfX[0]) );
352 :
353 0 : for( i = 0; i < nPointCount; i++ )
354 0 : pasGCPList[i].dfGCPX = oadfX[i];
355 0 : }
356 :
357 : /************************************************************************/
358 : /* ScanForGCPs() */
359 : /************************************************************************/
360 :
361 5 : void BSBDataset::ScanForGCPs( bool isNos, const char *pszFilename )
362 :
363 : {
364 : /* -------------------------------------------------------------------- */
365 : /* Collect GCPs as appropriate to source. */
366 : /* -------------------------------------------------------------------- */
367 5 : nGCPCount = 0;
368 :
369 5 : if ( isNos )
370 : {
371 0 : ScanForGCPsNos(pszFilename);
372 : } else {
373 5 : ScanForGCPsBSB();
374 : }
375 :
376 : /* -------------------------------------------------------------------- */
377 : /* Apply heuristics to re-wrap GCPs to maintain continguity */
378 : /* over the international dateline. */
379 : /* -------------------------------------------------------------------- */
380 5 : if( nGCPCount > 1 )
381 0 : GDALHeuristicDatelineWrapGCPs( nGCPCount, pasGCPList );
382 :
383 : /* -------------------------------------------------------------------- */
384 : /* Can we derive a reasonable coordinate system definition for */
385 : /* this file? For now we keep it simple, just handling */
386 : /* mercator. In the future we should consider others. */
387 : /* -------------------------------------------------------------------- */
388 5 : CPLString osUnderlyingSRS;
389 : int i;
390 :
391 15 : for( i = 0; psInfo->papszHeader[i] != NULL; i++ )
392 : {
393 15 : if( EQUALN(psInfo->papszHeader[i],"KNP/",4) )
394 : {
395 5 : const char *pszPR = strstr(psInfo->papszHeader[i],"PR=");
396 :
397 : // Capture whole line as metadata so some apps can do more
398 : // specific processing.
399 5 : SetMetadataItem( "BSB_KNP", psInfo->papszHeader[i] + 4 );
400 :
401 5 : if( pszPR == NULL )
402 : {
403 : /* no match */
404 : }
405 5 : else if( EQUALN(pszPR,"PR=MERCATOR", 11) )
406 : {
407 : // We somewhat arbitrarily select our first GCPX as our
408 : // central meridian. This is mostly helpful to ensure
409 : // that regions crossing the dateline will be contiguous
410 : // in mercator.
411 5 : osUnderlyingSRS.Printf( "PROJCS[\"Global Mercator\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.01745329251994328]],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]]", (int) pasGCPList[0].dfGCPX );
412 : }
413 :
414 5 : break;
415 : }
416 : }
417 :
418 : /* -------------------------------------------------------------------- */
419 : /* If we got an alternate underlying coordinate system, try */
420 : /* converting the GCPs to that coordinate system. */
421 : /* -------------------------------------------------------------------- */
422 5 : if( osUnderlyingSRS.length() > 0 )
423 : {
424 5 : OGRSpatialReference oWGS84_SRS, oProjected_SRS;
425 : OGRCoordinateTransformation *poCT;
426 :
427 5 : oWGS84_SRS.SetWellKnownGeogCS( "WGS84" );
428 5 : oProjected_SRS.SetFromUserInput( osUnderlyingSRS );
429 :
430 : poCT = OGRCreateCoordinateTransformation( &oWGS84_SRS,
431 5 : &oProjected_SRS );
432 5 : if( poCT != NULL )
433 : {
434 5 : for( i = 0; i < nGCPCount; i++ )
435 : {
436 : poCT->Transform( 1,
437 0 : &(pasGCPList[i].dfGCPX),
438 0 : &(pasGCPList[i].dfGCPY),
439 0 : &(pasGCPList[i].dfGCPZ) );
440 : }
441 :
442 5 : osGCPProjection = osUnderlyingSRS;
443 :
444 5 : delete poCT;
445 : }
446 : else
447 0 : CPLErrorReset();
448 : }
449 :
450 : /* -------------------------------------------------------------------- */
451 : /* Attempt to prepare a geotransform from the GCPs. */
452 : /* -------------------------------------------------------------------- */
453 5 : if( GDALGCPsToGeoTransform( nGCPCount, pasGCPList, adfGeoTransform,
454 : FALSE ) )
455 : {
456 0 : bGeoTransformSet = TRUE;
457 5 : }
458 5 : }
459 :
460 : /************************************************************************/
461 : /* ScanForGCPsNos() */
462 : /* */
463 : /* Nos files have an accompanying .geo file, that contains some */
464 : /* of the information normally contained in the header section */
465 : /* with BSB files. we try and open a file with the same name, */
466 : /* but a .geo extension, and look for lines like... */
467 : /* PointX=long lat line pixel (using the same naming system */
468 : /* as BSB) Point1=-22.0000 64.250000 197 744 */
469 : /************************************************************************/
470 :
471 0 : void BSBDataset::ScanForGCPsNos( const char *pszFilename )
472 : {
473 : char **Tokens;
474 : const char *geofile;
475 : const char *extension;
476 0 : int fileGCPCount=0;
477 :
478 0 : extension = CPLGetExtension(pszFilename);
479 :
480 : // pseudointelligently try and guess whether we want a .geo or a .GEO
481 0 : if (extension[1] == 'O')
482 : {
483 0 : geofile = CPLResetExtension( pszFilename, "GEO");
484 : } else {
485 0 : geofile = CPLResetExtension( pszFilename, "geo");
486 : }
487 :
488 0 : FILE *gfp = VSIFOpen( geofile, "r" ); // Text files
489 0 : if( gfp == NULL )
490 : {
491 : CPLError( CE_Failure, CPLE_OpenFailed,
492 0 : "Couldn't find a matching .GEO file: %s", geofile );
493 0 : return;
494 : }
495 :
496 0 : char *thisLine = (char *) CPLMalloc( 80 ); // FIXME
497 :
498 : // Count the GCPs (reference points) and seek the file pointer 'gfp' to the starting point
499 0 : while (fgets(thisLine, 80, gfp))
500 : {
501 0 : if( EQUALN(thisLine, "Point", 5) )
502 0 : fileGCPCount++;
503 : }
504 0 : VSIRewind( gfp );
505 :
506 : // Memory has not been allocated to fileGCPCount yet
507 0 : pasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),fileGCPCount+1);
508 :
509 0 : while (fgets(thisLine, 80, gfp))
510 : {
511 0 : if( EQUALN(thisLine, "Point", 5) )
512 : {
513 : // got a point line, turn it into a gcp
514 0 : Tokens = CSLTokenizeStringComplex(thisLine, "= ", FALSE, FALSE);
515 0 : if (CSLCount(Tokens) >= 5)
516 : {
517 0 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
518 0 : pasGCPList[nGCPCount].dfGCPX = atof(Tokens[1]);
519 0 : pasGCPList[nGCPCount].dfGCPY = atof(Tokens[2]);
520 0 : pasGCPList[nGCPCount].dfGCPPixel = atof(Tokens[4]);
521 0 : pasGCPList[nGCPCount].dfGCPLine = atof(Tokens[3]);
522 :
523 0 : CPLFree( pasGCPList[nGCPCount].pszId );
524 : char szName[50];
525 0 : sprintf( szName, "GCP_%d", nGCPCount+1 );
526 0 : pasGCPList[nGCPCount].pszId = CPLStrdup( szName );
527 :
528 0 : nGCPCount++;
529 : }
530 0 : CSLDestroy(Tokens);
531 : }
532 : }
533 :
534 0 : CPLFree(thisLine);
535 0 : VSIFClose(gfp);
536 : }
537 :
538 :
539 : /************************************************************************/
540 : /* ScanForGCPsBSB() */
541 : /************************************************************************/
542 :
543 5 : void BSBDataset::ScanForGCPsBSB()
544 : {
545 : /* -------------------------------------------------------------------- */
546 : /* Collect standalone GCPs. They look like: */
547 : /* */
548 : /* REF/1,115,2727,32.346666666667,-60.881666666667 */
549 : /* REF/n,pixel,line,lat,long */
550 : /* -------------------------------------------------------------------- */
551 5 : int fileGCPCount=0;
552 : int i;
553 :
554 : // Count the GCPs (reference points) in psInfo->papszHeader
555 655 : for( i = 0; psInfo->papszHeader[i] != NULL; i++ )
556 650 : if( EQUALN(psInfo->papszHeader[i],"REF/",4) )
557 0 : fileGCPCount++;
558 :
559 : // Memory has not been allocated to fileGCPCount yet
560 5 : pasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),fileGCPCount+1);
561 :
562 655 : for( i = 0; psInfo->papszHeader[i] != NULL; i++ )
563 : {
564 : char **papszTokens;
565 : char szName[50];
566 :
567 650 : if( !EQUALN(psInfo->papszHeader[i],"REF/",4) )
568 650 : continue;
569 :
570 : papszTokens =
571 0 : CSLTokenizeStringComplex( psInfo->papszHeader[i]+4, ",",
572 0 : FALSE, FALSE );
573 :
574 0 : if( CSLCount(papszTokens) > 4 )
575 : {
576 0 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
577 :
578 0 : pasGCPList[nGCPCount].dfGCPX = atof(papszTokens[4]);
579 0 : pasGCPList[nGCPCount].dfGCPY = atof(papszTokens[3]);
580 0 : pasGCPList[nGCPCount].dfGCPPixel = atof(papszTokens[1]);
581 0 : pasGCPList[nGCPCount].dfGCPLine = atof(papszTokens[2]);
582 :
583 0 : CPLFree( pasGCPList[nGCPCount].pszId );
584 0 : if( CSLCount(papszTokens) > 5 )
585 : {
586 0 : pasGCPList[nGCPCount].pszId = CPLStrdup(papszTokens[5]);
587 : }
588 : else
589 : {
590 0 : sprintf( szName, "GCP_%d", nGCPCount+1 );
591 0 : pasGCPList[nGCPCount].pszId = CPLStrdup( szName );
592 : }
593 :
594 0 : nGCPCount++;
595 : }
596 0 : CSLDestroy( papszTokens );
597 : }
598 5 : }
599 :
600 : /************************************************************************/
601 : /* Open() */
602 : /************************************************************************/
603 :
604 9629 : GDALDataset *BSBDataset::Open( GDALOpenInfo * poOpenInfo )
605 :
606 : {
607 : /* -------------------------------------------------------------------- */
608 : /* Check for BSB/ keyword. */
609 : /* -------------------------------------------------------------------- */
610 : int i;
611 9629 : bool isNos = false;
612 :
613 9629 : if( poOpenInfo->nHeaderBytes < 1000 )
614 8888 : return NULL;
615 :
616 751596 : for( i = 0; i < poOpenInfo->nHeaderBytes - 4; i++ )
617 : {
618 754407 : if( poOpenInfo->pabyHeader[i+0] == 'B'
619 3532 : && poOpenInfo->pabyHeader[i+1] == 'S'
620 10 : && poOpenInfo->pabyHeader[i+2] == 'B'
621 5 : && poOpenInfo->pabyHeader[i+3] == '/' )
622 5 : break;
623 752879 : if( poOpenInfo->pabyHeader[i+0] == 'N'
624 1954 : && poOpenInfo->pabyHeader[i+1] == 'O'
625 63 : && poOpenInfo->pabyHeader[i+2] == 'S'
626 7 : && poOpenInfo->pabyHeader[i+3] == '/' )
627 : {
628 0 : isNos = true;
629 0 : break;
630 : }
631 751575 : if( poOpenInfo->pabyHeader[i+0] == 'W'
632 704 : && poOpenInfo->pabyHeader[i+1] == 'X'
633 16 : && poOpenInfo->pabyHeader[i+2] == '\\'
634 0 : && poOpenInfo->pabyHeader[i+3] == '8' )
635 0 : break;
636 : }
637 :
638 741 : if( i == poOpenInfo->nHeaderBytes - 4 )
639 736 : return NULL;
640 :
641 : /* Additional test to avoid false positive. See #2881 */
642 5 : const char* pszRA = strstr((const char*)poOpenInfo->pabyHeader + i, "RA=");
643 5 : if (pszRA == NULL) /* This may be a NO1 file */
644 0 : pszRA = strstr((const char*)poOpenInfo->pabyHeader + i, "[JF");
645 5 : if (pszRA == NULL || pszRA - ((const char*)poOpenInfo->pabyHeader + i) > 100 )
646 0 : return NULL;
647 :
648 5 : if( poOpenInfo->eAccess == GA_Update )
649 : {
650 : CPLError( CE_Failure, CPLE_NotSupported,
651 : "The BSB driver does not support update access to existing"
652 0 : " datasets.\n" );
653 0 : return NULL;
654 : }
655 :
656 : /* -------------------------------------------------------------------- */
657 : /* Create a corresponding GDALDataset. */
658 : /* -------------------------------------------------------------------- */
659 : BSBDataset *poDS;
660 :
661 5 : poDS = new BSBDataset();
662 :
663 : /* -------------------------------------------------------------------- */
664 : /* Open the file. */
665 : /* -------------------------------------------------------------------- */
666 5 : poDS->psInfo = BSBOpen( poOpenInfo->pszFilename );
667 5 : if( poDS->psInfo == NULL )
668 : {
669 0 : delete poDS;
670 0 : return NULL;
671 : }
672 :
673 5 : poDS->nRasterXSize = poDS->psInfo->nXSize;
674 5 : poDS->nRasterYSize = poDS->psInfo->nYSize;
675 :
676 : /* -------------------------------------------------------------------- */
677 : /* Create band information objects. */
678 : /* -------------------------------------------------------------------- */
679 5 : poDS->SetBand( 1, new BSBRasterBand( poDS ));
680 :
681 5 : poDS->ScanForGCPs( isNos, poOpenInfo->pszFilename );
682 :
683 : /* -------------------------------------------------------------------- */
684 : /* Initialize any PAM information. */
685 : /* -------------------------------------------------------------------- */
686 5 : poDS->SetDescription( poOpenInfo->pszFilename );
687 5 : poDS->TryLoadXML();
688 :
689 5 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
690 :
691 5 : return( poDS );
692 : }
693 :
694 : /************************************************************************/
695 : /* GetGCPCount() */
696 : /************************************************************************/
697 :
698 0 : int BSBDataset::GetGCPCount()
699 :
700 : {
701 0 : return nGCPCount;
702 : }
703 :
704 : /************************************************************************/
705 : /* GetGCPProjection() */
706 : /************************************************************************/
707 :
708 0 : const char *BSBDataset::GetGCPProjection()
709 :
710 : {
711 0 : return osGCPProjection;
712 : }
713 :
714 : /************************************************************************/
715 : /* GetGCP() */
716 : /************************************************************************/
717 :
718 0 : const GDAL_GCP *BSBDataset::GetGCPs()
719 :
720 : {
721 0 : return pasGCPList;
722 : }
723 :
724 : #ifdef BSB_CREATE
725 :
726 : /************************************************************************/
727 : /* BSBIsSRSOK() */
728 : /************************************************************************/
729 :
730 : static int BSBIsSRSOK(const char *pszWKT)
731 : {
732 : int bOK = FALSE;
733 : OGRSpatialReference oSRS, oSRS_WGS84, oSRS_NAD83;
734 :
735 : if( pszWKT != NULL && pszWKT[0] != '\0' )
736 : {
737 : char* pszTmpWKT = (char*)pszWKT;
738 : oSRS.importFromWkt( &pszTmpWKT );
739 :
740 : oSRS_WGS84.SetWellKnownGeogCS( "WGS84" );
741 : oSRS_NAD83.SetWellKnownGeogCS( "NAD83" );
742 : if ( (oSRS.IsSameGeogCS(&oSRS_WGS84) || oSRS.IsSameGeogCS(&oSRS_NAD83)) &&
743 : oSRS.IsGeographic() && oSRS.GetPrimeMeridian() == 0.0 )
744 : {
745 : bOK = TRUE;
746 : }
747 : }
748 :
749 : if (!bOK)
750 : {
751 : CPLError(CE_Warning, CPLE_NotSupported,
752 : "BSB only supports WGS84 or NAD83 geographic projections.\n");
753 : }
754 :
755 : return bOK;
756 : }
757 :
758 : /************************************************************************/
759 : /* BSBCreateCopy() */
760 : /************************************************************************/
761 :
762 : static GDALDataset *
763 : BSBCreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
764 : int bStrict, char ** papszOptions,
765 : GDALProgressFunc pfnProgress, void * pProgressData )
766 :
767 : {
768 : int nBands = poSrcDS->GetRasterCount();
769 : int nXSize = poSrcDS->GetRasterXSize();
770 : int nYSize = poSrcDS->GetRasterYSize();
771 :
772 : /* -------------------------------------------------------------------- */
773 : /* Some some rudimentary checks */
774 : /* -------------------------------------------------------------------- */
775 : if( nBands != 1 )
776 : {
777 : CPLError( CE_Failure, CPLE_NotSupported,
778 : "BSB driver only supports one band images.\n" );
779 :
780 : return NULL;
781 : }
782 :
783 : if( poSrcDS->GetRasterBand(1)->GetRasterDataType() != GDT_Byte
784 : && bStrict )
785 : {
786 : CPLError( CE_Failure, CPLE_NotSupported,
787 : "BSB driver doesn't support data type %s. "
788 : "Only eight bit bands supported.\n",
789 : GDALGetDataTypeName(
790 : poSrcDS->GetRasterBand(1)->GetRasterDataType()) );
791 :
792 : return NULL;
793 : }
794 :
795 : /* -------------------------------------------------------------------- */
796 : /* Open the output file. */
797 : /* -------------------------------------------------------------------- */
798 : BSBInfo *psBSB;
799 :
800 : psBSB = BSBCreate( pszFilename, 0, 200, nXSize, nYSize );
801 : if( psBSB == NULL )
802 : return NULL;
803 :
804 : /* -------------------------------------------------------------------- */
805 : /* Prepare initial color table.colortable. */
806 : /* -------------------------------------------------------------------- */
807 : GDALRasterBand *poBand = poSrcDS->GetRasterBand(1);
808 : int iColor;
809 : unsigned char abyPCT[771];
810 : int nPCTSize;
811 : int anRemap[256];
812 :
813 : abyPCT[0] = 0;
814 : abyPCT[1] = 0;
815 : abyPCT[2] = 0;
816 :
817 : if( poBand->GetColorTable() == NULL )
818 : {
819 : /* map greyscale down to 63 grey levels. */
820 : for( iColor = 0; iColor < 256; iColor++ )
821 : {
822 : int nOutValue = (int) (iColor / 4.1) + 1;
823 :
824 : anRemap[iColor] = nOutValue;
825 : abyPCT[nOutValue*3 + 0] = (unsigned char) iColor;
826 : abyPCT[nOutValue*3 + 1] = (unsigned char) iColor;
827 : abyPCT[nOutValue*3 + 2] = (unsigned char) iColor;
828 : }
829 : nPCTSize = 64;
830 : }
831 : else
832 : {
833 : GDALColorTable *poCT = poBand->GetColorTable();
834 : int nColorTableSize = poCT->GetColorEntryCount();
835 : if (nColorTableSize > 255)
836 : nColorTableSize = 255;
837 :
838 : for( iColor = 0; iColor < nColorTableSize; iColor++ )
839 : {
840 : GDALColorEntry sEntry;
841 :
842 : poCT->GetColorEntryAsRGB( iColor, &sEntry );
843 :
844 : anRemap[iColor] = iColor + 1;
845 : abyPCT[(iColor+1)*3 + 0] = (unsigned char) sEntry.c1;
846 : abyPCT[(iColor+1)*3 + 1] = (unsigned char) sEntry.c2;
847 : abyPCT[(iColor+1)*3 + 2] = (unsigned char) sEntry.c3;
848 : }
849 :
850 : nPCTSize = nColorTableSize + 1;
851 :
852 : // Add entries for pixel values which apparently will not occur.
853 : for( iColor = nPCTSize; iColor < 256; iColor++ )
854 : anRemap[iColor] = 1;
855 : }
856 :
857 : /* -------------------------------------------------------------------- */
858 : /* Boil out all duplicate entries. */
859 : /* -------------------------------------------------------------------- */
860 : int i;
861 :
862 : for( i = 1; i < nPCTSize-1; i++ )
863 : {
864 : int j;
865 :
866 : for( j = i+1; j < nPCTSize; j++ )
867 : {
868 : if( abyPCT[i*3+0] == abyPCT[j*3+0]
869 : && abyPCT[i*3+1] == abyPCT[j*3+1]
870 : && abyPCT[i*3+2] == abyPCT[j*3+2] )
871 : {
872 : int k;
873 :
874 : nPCTSize--;
875 : abyPCT[j*3+0] = abyPCT[nPCTSize*3+0];
876 : abyPCT[j*3+1] = abyPCT[nPCTSize*3+1];
877 : abyPCT[j*3+2] = abyPCT[nPCTSize*3+2];
878 :
879 : for( k = 0; k < 256; k++ )
880 : {
881 : // merge matching entries.
882 : if( anRemap[k] == j )
883 : anRemap[k] = i;
884 :
885 : // shift the last PCT entry into the new hole.
886 : if( anRemap[k] == nPCTSize )
887 : anRemap[k] = j;
888 : }
889 : }
890 : }
891 : }
892 :
893 : /* -------------------------------------------------------------------- */
894 : /* Boil out all duplicate entries. */
895 : /* -------------------------------------------------------------------- */
896 : if( nPCTSize > 128 )
897 : {
898 : CPLError( CE_Warning, CPLE_AppDefined,
899 : "Having to merge color table entries to reduce %d real\n"
900 : "color table entries down to 127 values.",
901 : nPCTSize );
902 : }
903 :
904 : while( nPCTSize > 128 )
905 : {
906 : int nBestRange = 768;
907 : int iBestMatch1=-1, iBestMatch2=-1;
908 :
909 : // Find the closest pair of color table entries.
910 :
911 : for( i = 1; i < nPCTSize-1; i++ )
912 : {
913 : int j;
914 :
915 : for( j = i+1; j < nPCTSize; j++ )
916 : {
917 : int nRange = ABS(abyPCT[i*3+0] - abyPCT[j*3+0])
918 : + ABS(abyPCT[i*3+1] - abyPCT[j*3+1])
919 : + ABS(abyPCT[i*3+2] - abyPCT[j*3+2]);
920 :
921 : if( nRange < nBestRange )
922 : {
923 : iBestMatch1 = i;
924 : iBestMatch2 = j;
925 : nBestRange = nRange;
926 : }
927 : }
928 : }
929 :
930 : // Merge the second entry into the first.
931 : nPCTSize--;
932 : abyPCT[iBestMatch2*3+0] = abyPCT[nPCTSize*3+0];
933 : abyPCT[iBestMatch2*3+1] = abyPCT[nPCTSize*3+1];
934 : abyPCT[iBestMatch2*3+2] = abyPCT[nPCTSize*3+2];
935 :
936 : for( i = 0; i < 256; i++ )
937 : {
938 : // merge matching entries.
939 : if( anRemap[i] == iBestMatch2 )
940 : anRemap[i] = iBestMatch1;
941 :
942 : // shift the last PCT entry into the new hole.
943 : if( anRemap[i] == nPCTSize )
944 : anRemap[i] = iBestMatch2;
945 : }
946 : }
947 :
948 : /* -------------------------------------------------------------------- */
949 : /* Write the PCT. */
950 : /* -------------------------------------------------------------------- */
951 : if( !BSBWritePCT( psBSB, nPCTSize, abyPCT ) )
952 : {
953 : BSBClose( psBSB );
954 : return NULL;
955 : }
956 :
957 : /* -------------------------------------------------------------------- */
958 : /* Write the GCPs. */
959 : /* -------------------------------------------------------------------- */
960 : double adfGeoTransform[6];
961 : int nGCPCount = poSrcDS->GetGCPCount();
962 : if (nGCPCount)
963 : {
964 : const char* pszGCPProjection = poSrcDS->GetGCPProjection();
965 : if ( BSBIsSRSOK(pszGCPProjection) )
966 : {
967 : const GDAL_GCP * pasGCPList = poSrcDS->GetGCPs();
968 : for( i = 0; i < nGCPCount; i++ )
969 : {
970 : VSIFPrintfL( psBSB->fp,
971 : "REF/%d,%f,%f,%f,%f\n",
972 : i+1,
973 : pasGCPList[i].dfGCPPixel, pasGCPList[i].dfGCPLine,
974 : pasGCPList[i].dfGCPY, pasGCPList[i].dfGCPX);
975 : }
976 : }
977 : }
978 : else if (poSrcDS->GetGeoTransform(adfGeoTransform) == CE_None)
979 : {
980 : const char* pszProjection = poSrcDS->GetProjectionRef();
981 : if ( BSBIsSRSOK(pszProjection) )
982 : {
983 : VSIFPrintfL( psBSB->fp,
984 : "REF/%d,%d,%d,%f,%f\n",
985 : 1,
986 : 0, 0,
987 : adfGeoTransform[3] + 0 * adfGeoTransform[4] + 0 * adfGeoTransform[5],
988 : adfGeoTransform[0] + 0 * adfGeoTransform[1] + 0 * adfGeoTransform[2]);
989 : VSIFPrintfL( psBSB->fp,
990 : "REF/%d,%d,%d,%f,%f\n",
991 : 2,
992 : nXSize, 0,
993 : adfGeoTransform[3] + nXSize * adfGeoTransform[4] + 0 * adfGeoTransform[5],
994 : adfGeoTransform[0] + nXSize * adfGeoTransform[1] + 0 * adfGeoTransform[2]);
995 : VSIFPrintfL( psBSB->fp,
996 : "REF/%d,%d,%d,%f,%f\n",
997 : 3,
998 : nXSize, nYSize,
999 : adfGeoTransform[3] + nXSize * adfGeoTransform[4] + nYSize * adfGeoTransform[5],
1000 : adfGeoTransform[0] + nXSize * adfGeoTransform[1] + nYSize * adfGeoTransform[2]);
1001 : VSIFPrintfL( psBSB->fp,
1002 : "REF/%d,%d,%d,%f,%f\n",
1003 : 4,
1004 : 0, nYSize,
1005 : adfGeoTransform[3] + 0 * adfGeoTransform[4] + nYSize * adfGeoTransform[5],
1006 : adfGeoTransform[0] + 0 * adfGeoTransform[1] + nYSize * adfGeoTransform[2]);
1007 : }
1008 : }
1009 :
1010 : /* -------------------------------------------------------------------- */
1011 : /* Loop over image, copying image data. */
1012 : /* -------------------------------------------------------------------- */
1013 : GByte *pabyScanline;
1014 : CPLErr eErr = CE_None;
1015 :
1016 : pabyScanline = (GByte *) CPLMalloc( nXSize );
1017 :
1018 : for( int iLine = 0; iLine < nYSize && eErr == CE_None; iLine++ )
1019 : {
1020 : eErr = poBand->RasterIO( GF_Read, 0, iLine, nXSize, 1,
1021 : pabyScanline, nXSize, 1, GDT_Byte,
1022 : nBands, nBands * nXSize );
1023 : if( eErr == CE_None )
1024 : {
1025 : for( i = 0; i < nXSize; i++ )
1026 : pabyScanline[i] = (GByte) anRemap[pabyScanline[i]];
1027 :
1028 : if( !BSBWriteScanline( psBSB, pabyScanline ) )
1029 : eErr = CE_Failure;
1030 : }
1031 : }
1032 :
1033 : CPLFree( pabyScanline );
1034 :
1035 : /* -------------------------------------------------------------------- */
1036 : /* cleanup */
1037 : /* -------------------------------------------------------------------- */
1038 : BSBClose( psBSB );
1039 :
1040 : if( eErr != CE_None )
1041 : {
1042 : VSIUnlink( pszFilename );
1043 : return NULL;
1044 : }
1045 : else
1046 : return (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );
1047 : }
1048 : #endif
1049 :
1050 : /************************************************************************/
1051 : /* GDALRegister_BSB() */
1052 : /************************************************************************/
1053 :
1054 338 : void GDALRegister_BSB()
1055 :
1056 : {
1057 : GDALDriver *poDriver;
1058 :
1059 338 : if( GDALGetDriverByName( "BSB" ) == NULL )
1060 : {
1061 336 : poDriver = new GDALDriver();
1062 :
1063 336 : poDriver->SetDescription( "BSB" );
1064 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
1065 336 : "Maptech BSB Nautical Charts" );
1066 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
1067 336 : "frmt_various.html#BSB" );
1068 : #ifdef BSB_CREATE
1069 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES, "Byte" );
1070 : #endif
1071 336 : poDriver->pfnOpen = BSBDataset::Open;
1072 : #ifdef BSB_CREATE
1073 : poDriver->pfnCreateCopy = BSBCreateCopy;
1074 : #endif
1075 :
1076 336 : GetGDALDriverManager()->RegisterDriver( poDriver );
1077 : }
1078 338 : }
|
