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