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