Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Polarimetric Workstation
4 : * Purpose: Convair PolGASP data (.img/.hdr format).
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2004, Frank Warmerdam <warmerdam@pobox.com>
9 : * Copyright (c) 2009, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_string.h"
15 : #include "gdal_frmts.h"
16 : #include "ogr_spatialref.h"
17 : #include "rawdataset.h"
18 :
19 : #include <cinttypes>
20 : #include <vector>
21 :
22 : /************************************************************************/
23 : /* ==================================================================== */
24 : /* CPGDataset */
25 : /* ==================================================================== */
26 : /************************************************************************/
27 :
28 : class SIRC_QSLCRasterBand;
29 : class CPG_STOKESRasterBand;
30 :
31 : class CPGDataset final : public RawDataset
32 : {
33 : friend class SIRC_QSLCRasterBand;
34 : friend class CPG_STOKESRasterBand;
35 :
36 : static constexpr int NUMBER_OF_BANDS = 4;
37 : std::vector<VSILFILE *> afpImage{NUMBER_OF_BANDS};
38 : std::vector<CPLString> aosImageFilenames{};
39 :
40 : int nGCPCount;
41 : GDAL_GCP *pasGCPList;
42 : OGRSpatialReference m_oGCPSRS{};
43 :
44 : GDALGeoTransform m_gt{};
45 : OGRSpatialReference m_oSRS{};
46 :
47 : int nLoadedStokesLine;
48 : float *padfStokesMatrix;
49 :
50 : Interleave eInterleave = Interleave::BSQ;
51 : static int AdjustFilename(char **, const char *, const char *);
52 : static int FindType1(const char *pszWorkname);
53 : static int FindType2(const char *pszWorkname);
54 : #ifdef notdef
55 : static int FindType3(const char *pszWorkname);
56 : #endif
57 : static GDALDataset *InitializeType1Or2Dataset(const char *pszWorkname);
58 : #ifdef notdef
59 : static GDALDataset *InitializeType3Dataset(const char *pszWorkname);
60 : #endif
61 : CPLErr LoadStokesLine(int iLine, int bNativeOrder);
62 :
63 : CPL_DISALLOW_COPY_ASSIGN(CPGDataset)
64 :
65 : CPLErr Close() override;
66 :
67 : public:
68 : CPGDataset();
69 : ~CPGDataset() override;
70 :
71 : int GetGCPCount() override;
72 :
73 0 : const OGRSpatialReference *GetGCPSpatialRef() const override
74 : {
75 0 : return m_oGCPSRS.IsEmpty() ? nullptr : &m_oGCPSRS;
76 : }
77 :
78 : const GDAL_GCP *GetGCPs() override;
79 :
80 0 : const OGRSpatialReference *GetSpatialRef() const override
81 : {
82 0 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
83 : }
84 :
85 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
86 :
87 : char **GetFileList() override;
88 :
89 : static GDALDataset *Open(GDALOpenInfo *);
90 : };
91 :
92 : /************************************************************************/
93 : /* CPGDataset() */
94 : /************************************************************************/
95 :
96 2 : CPGDataset::CPGDataset()
97 : : nGCPCount(0), pasGCPList(nullptr), nLoadedStokesLine(-1),
98 2 : padfStokesMatrix(nullptr)
99 : {
100 2 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
101 2 : m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
102 2 : }
103 :
104 : /************************************************************************/
105 : /* ~CPGDataset() */
106 : /************************************************************************/
107 :
108 4 : CPGDataset::~CPGDataset()
109 :
110 : {
111 2 : CPGDataset::Close();
112 4 : }
113 :
114 : /************************************************************************/
115 : /* Close() */
116 : /************************************************************************/
117 :
118 4 : CPLErr CPGDataset::Close()
119 : {
120 4 : CPLErr eErr = CE_None;
121 4 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
122 : {
123 2 : if (CPGDataset::FlushCache(true) != CE_None)
124 0 : eErr = CE_Failure;
125 :
126 10 : for (auto poFile : afpImage)
127 : {
128 8 : if (poFile != nullptr)
129 2 : VSIFCloseL(poFile);
130 : }
131 :
132 2 : if (nGCPCount > 0)
133 : {
134 2 : GDALDeinitGCPs(nGCPCount, pasGCPList);
135 2 : CPLFree(pasGCPList);
136 : }
137 :
138 2 : CPLFree(padfStokesMatrix);
139 :
140 2 : if (GDALPamDataset::Close() != CE_None)
141 0 : eErr = CE_Failure;
142 : }
143 4 : return eErr;
144 : }
145 :
146 : /************************************************************************/
147 : /* GetFileList() */
148 : /************************************************************************/
149 :
150 1 : char **CPGDataset::GetFileList()
151 :
152 : {
153 1 : char **papszFileList = RawDataset::GetFileList();
154 2 : for (size_t i = 0; i < aosImageFilenames.size(); ++i)
155 1 : papszFileList = CSLAddString(papszFileList, aosImageFilenames[i]);
156 1 : return papszFileList;
157 : }
158 :
159 : /************************************************************************/
160 : /* ==================================================================== */
161 : /* SIRC_QSLCPRasterBand */
162 : /* ==================================================================== */
163 : /************************************************************************/
164 :
165 : class SIRC_QSLCRasterBand final : public GDALRasterBand
166 : {
167 : friend class CPGDataset;
168 :
169 : public:
170 : SIRC_QSLCRasterBand(CPGDataset *, int, GDALDataType);
171 :
172 16 : ~SIRC_QSLCRasterBand() override
173 8 : {
174 16 : }
175 :
176 : CPLErr IReadBlock(int, int, void *) override;
177 : };
178 :
179 : constexpr int M11 = 0;
180 : // constexpr int M12 = 1;
181 : constexpr int M13 = 2;
182 : constexpr int M14 = 3;
183 : // constexpr int M21 = 4;
184 : constexpr int M22 = 5;
185 : constexpr int M23 = 6;
186 : constexpr int M24 = 7;
187 : constexpr int M31 = 8;
188 : constexpr int M32 = 9;
189 : constexpr int M33 = 10;
190 : constexpr int M34 = 11;
191 : constexpr int M41 = 12;
192 : constexpr int M42 = 13;
193 : constexpr int M43 = 14;
194 : constexpr int M44 = 15;
195 :
196 : /************************************************************************/
197 : /* ==================================================================== */
198 : /* CPG_STOKESRasterBand */
199 : /* ==================================================================== */
200 : /************************************************************************/
201 :
202 : class CPG_STOKESRasterBand final : public GDALRasterBand
203 : {
204 : friend class CPGDataset;
205 :
206 : int bNativeOrder;
207 :
208 : public:
209 : CPG_STOKESRasterBand(GDALDataset *poDS, GDALDataType eType,
210 : int bNativeOrder);
211 :
212 0 : ~CPG_STOKESRasterBand() override
213 0 : {
214 0 : }
215 :
216 : CPLErr IReadBlock(int, int, void *) override;
217 : };
218 :
219 : /************************************************************************/
220 : /* AdjustFilename() */
221 : /* */
222 : /* Try to find the file with the request polarization and */
223 : /* extension and update the passed filename accordingly. */
224 : /* */
225 : /* Return TRUE if file found otherwise FALSE. */
226 : /************************************************************************/
227 :
228 8 : int CPGDataset::AdjustFilename(char **pszFilename, const char *pszPolarization,
229 : const char *pszExtension)
230 :
231 : {
232 : // TODO: Eventually we should handle upper/lower case.
233 8 : if (EQUAL(pszPolarization, "stokes"))
234 : {
235 : const std::string osNewName =
236 0 : CPLResetExtensionSafe(*pszFilename, pszExtension);
237 0 : CPLFree(*pszFilename);
238 0 : *pszFilename = CPLStrdup(osNewName.c_str());
239 : }
240 8 : else if (strlen(pszPolarization) == 2)
241 : {
242 2 : char *subptr = strstr(*pszFilename, "hh");
243 2 : if (subptr == nullptr)
244 2 : subptr = strstr(*pszFilename, "hv");
245 2 : if (subptr == nullptr)
246 2 : subptr = strstr(*pszFilename, "vv");
247 2 : if (subptr == nullptr)
248 2 : subptr = strstr(*pszFilename, "vh");
249 2 : if (subptr == nullptr)
250 2 : return FALSE;
251 :
252 0 : strncpy(subptr, pszPolarization, 2);
253 : const std::string osNewName =
254 0 : CPLResetExtensionSafe(*pszFilename, pszExtension);
255 0 : CPLFree(*pszFilename);
256 0 : *pszFilename = CPLStrdup(osNewName.c_str());
257 : }
258 : else
259 : {
260 : const std::string osNewName =
261 6 : CPLResetExtensionSafe(*pszFilename, pszExtension);
262 6 : CPLFree(*pszFilename);
263 6 : *pszFilename = CPLStrdup(osNewName.c_str());
264 : }
265 : VSIStatBufL sStatBuf;
266 6 : return VSIStatL(*pszFilename, &sStatBuf) == 0;
267 : }
268 :
269 : /************************************************************************/
270 : /* Search for the various types of Convair filesets */
271 : /* Return TRUE for a match, FALSE for no match */
272 : /************************************************************************/
273 33520 : int CPGDataset::FindType1(const char *pszFilename)
274 : {
275 33520 : const int nNameLen = static_cast<int>(strlen(pszFilename));
276 :
277 33520 : if ((strstr(pszFilename, "sso") == nullptr) &&
278 33513 : (strstr(pszFilename, "polgasp") == nullptr))
279 33513 : return FALSE;
280 :
281 7 : if ((strlen(pszFilename) < 5) ||
282 0 : (!EQUAL(pszFilename + nNameLen - 4, ".hdr") &&
283 0 : !EQUAL(pszFilename + nNameLen - 4, ".img")))
284 0 : return FALSE;
285 :
286 : /* Expect all bands and headers to be present */
287 7 : char *pszTemp = CPLStrdup(pszFilename);
288 :
289 0 : const bool bNotFound = !AdjustFilename(&pszTemp, "hh", "img") ||
290 0 : !AdjustFilename(&pszTemp, "hh", "hdr") ||
291 0 : !AdjustFilename(&pszTemp, "hv", "img") ||
292 0 : !AdjustFilename(&pszTemp, "hv", "hdr") ||
293 0 : !AdjustFilename(&pszTemp, "vh", "img") ||
294 0 : !AdjustFilename(&pszTemp, "vh", "hdr") ||
295 0 : !AdjustFilename(&pszTemp, "vv", "img") ||
296 0 : !AdjustFilename(&pszTemp, "vv", "hdr");
297 :
298 0 : CPLFree(pszTemp);
299 :
300 0 : return !bNotFound;
301 : }
302 :
303 33520 : int CPGDataset::FindType2(const char *pszFilename)
304 : {
305 33520 : const int nNameLen = static_cast<int>(strlen(pszFilename));
306 :
307 33520 : if ((strlen(pszFilename) < 9) ||
308 32812 : (!EQUAL(pszFilename + nNameLen - 8, "SIRC.hdr") &&
309 32812 : !EQUAL(pszFilename + nNameLen - 8, "SIRC.img")))
310 33513 : return FALSE;
311 :
312 7 : char *pszTemp = CPLStrdup(pszFilename);
313 4 : const bool bNotFound = !AdjustFilename(&pszTemp, "", "img") ||
314 2 : !AdjustFilename(&pszTemp, "", "hdr");
315 2 : CPLFree(pszTemp);
316 :
317 2 : return !bNotFound;
318 : }
319 :
320 : #ifdef notdef
321 : int CPGDataset::FindType3(const char *pszFilename)
322 : {
323 : const int nNameLen = static_cast<int>(strlen(pszFilename));
324 :
325 : if ((strstr(pszFilename, "sso") == NULL) &&
326 : (strstr(pszFilename, "polgasp") == NULL))
327 : return FALSE;
328 :
329 : if ((strlen(pszFilename) < 9) ||
330 : (!EQUAL(pszFilename + nNameLen - 4, ".img") &&
331 : !EQUAL(pszFilename + nNameLen - 8, ".img_def")))
332 : return FALSE;
333 :
334 : char *pszTemp = CPLStrdup(pszFilename);
335 : const bool bNotFound = !AdjustFilename(&pszTemp, "stokes", "img") ||
336 : !AdjustFilename(&pszTemp, "stokes", "img_def");
337 : CPLFree(pszTemp);
338 :
339 : return !bNotFound;
340 : }
341 : #endif
342 :
343 : /************************************************************************/
344 : /* LoadStokesLine() */
345 : /************************************************************************/
346 :
347 0 : CPLErr CPGDataset::LoadStokesLine(int iLine, int bNativeOrder)
348 :
349 : {
350 0 : if (iLine == nLoadedStokesLine)
351 0 : return CE_None;
352 :
353 0 : constexpr int nDataSize = static_cast<int>(sizeof(float));
354 :
355 : /* -------------------------------------------------------------------- */
356 : /* allocate working buffers if we don't have them already. */
357 : /* -------------------------------------------------------------------- */
358 0 : if (padfStokesMatrix == nullptr)
359 : {
360 0 : padfStokesMatrix = reinterpret_cast<float *>(
361 0 : CPLMalloc(sizeof(float) * nRasterXSize * 16));
362 : }
363 :
364 : /* -------------------------------------------------------------------- */
365 : /* Load all the pixel data associated with this scanline. */
366 : /* Retains same interleaving as original dataset. */
367 : /* -------------------------------------------------------------------- */
368 0 : if (eInterleave == Interleave::BIP)
369 : {
370 0 : const int offset = nRasterXSize * iLine * nDataSize * 16;
371 0 : const int nBytesToRead = nDataSize * nRasterXSize * 16;
372 0 : if ((VSIFSeekL(afpImage[0], offset, SEEK_SET) != 0) ||
373 : static_cast<int>(
374 0 : VSIFReadL(reinterpret_cast<GByte *>(padfStokesMatrix), 1,
375 0 : nBytesToRead, afpImage[0])) != nBytesToRead)
376 : {
377 0 : CPLError(CE_Failure, CPLE_FileIO,
378 : "Error reading %d bytes of Stokes Convair at offset %d.\n"
379 : "Reading file %s failed.",
380 0 : nBytesToRead, offset, GetDescription());
381 0 : CPLFree(padfStokesMatrix);
382 0 : padfStokesMatrix = nullptr;
383 0 : nLoadedStokesLine = -1;
384 0 : return CE_Failure;
385 : }
386 : }
387 0 : else if (eInterleave == Interleave::BIL)
388 : {
389 0 : for (int band_index = 0; band_index < 16; band_index++)
390 : {
391 0 : const int offset =
392 0 : nDataSize * (nRasterXSize * iLine + nRasterXSize * band_index);
393 0 : const int nBytesToRead = nDataSize * nRasterXSize;
394 0 : if ((VSIFSeekL(afpImage[0], offset, SEEK_SET) != 0) ||
395 : static_cast<int>(
396 0 : VSIFReadL(reinterpret_cast<GByte *>(
397 0 : padfStokesMatrix + nBytesToRead * band_index),
398 0 : 1, nBytesToRead, afpImage[0])) != nBytesToRead)
399 : {
400 0 : CPLError(
401 : CE_Failure, CPLE_FileIO,
402 : "Error reading %d bytes of Stokes Convair at offset %d.\n"
403 : "Reading file %s failed.",
404 0 : nBytesToRead, offset, GetDescription());
405 0 : CPLFree(padfStokesMatrix);
406 0 : padfStokesMatrix = nullptr;
407 0 : nLoadedStokesLine = -1;
408 0 : return CE_Failure;
409 : }
410 : }
411 : }
412 : else
413 : {
414 0 : for (int band_index = 0; band_index < 16; band_index++)
415 : {
416 0 : const int offset =
417 0 : nDataSize * (nRasterXSize * iLine +
418 0 : nRasterXSize * nRasterYSize * band_index);
419 0 : const int nBytesToRead = nDataSize * nRasterXSize;
420 0 : if ((VSIFSeekL(afpImage[0], offset, SEEK_SET) != 0) ||
421 : static_cast<int>(
422 0 : VSIFReadL(reinterpret_cast<GByte *>(
423 0 : padfStokesMatrix + nBytesToRead * band_index),
424 0 : 1, nBytesToRead, afpImage[0])) != nBytesToRead)
425 : {
426 0 : CPLError(
427 : CE_Failure, CPLE_FileIO,
428 : "Error reading %d bytes of Stokes Convair at offset %d.\n"
429 : "Reading file %s failed.",
430 0 : nBytesToRead, offset, GetDescription());
431 0 : CPLFree(padfStokesMatrix);
432 0 : padfStokesMatrix = nullptr;
433 0 : nLoadedStokesLine = -1;
434 0 : return CE_Failure;
435 : }
436 : }
437 : }
438 :
439 0 : if (!bNativeOrder)
440 0 : GDALSwapWords(padfStokesMatrix, nDataSize, nRasterXSize * 16,
441 : nDataSize);
442 :
443 0 : nLoadedStokesLine = iLine;
444 :
445 0 : return CE_None;
446 : }
447 :
448 : /************************************************************************/
449 : /* Parse header information and initialize dataset for the */
450 : /* appropriate Convair dataset style. */
451 : /* Returns dataset if successful; NULL if there was a problem. */
452 : /************************************************************************/
453 :
454 2 : GDALDataset *CPGDataset::InitializeType1Or2Dataset(const char *pszFilename)
455 : {
456 :
457 : /* -------------------------------------------------------------------- */
458 : /* Read the .hdr file (the hh one for the .sso and polgasp cases) */
459 : /* and parse it. */
460 : /* -------------------------------------------------------------------- */
461 2 : int nLines = 0;
462 2 : int nSamples = 0;
463 2 : int nError = 0;
464 :
465 : /* Parameters required for pseudo-geocoding. GCPs map */
466 : /* slant range to ground range at 16 points. */
467 2 : int iGeoParamsFound = 0;
468 2 : int itransposed = 0;
469 2 : double dfaltitude = 0.0;
470 2 : double dfnear_srd = 0.0;
471 2 : double dfsample_size = 0.0;
472 2 : double dfsample_size_az = 0.0;
473 :
474 : /* Parameters in geogratis geocoded images */
475 2 : int iUTMParamsFound = 0;
476 2 : int iUTMZone = 0;
477 : // int iCorner = 0;
478 2 : double dfnorth = 0.0;
479 2 : double dfeast = 0.0;
480 :
481 4 : std::string osWorkName;
482 : {
483 2 : char *pszWorkname = CPLStrdup(pszFilename);
484 2 : AdjustFilename(&pszWorkname, "hh", "hdr");
485 2 : osWorkName = pszWorkname;
486 2 : CPLFree(pszWorkname);
487 : }
488 :
489 2 : char **papszHdrLines = CSLLoad(osWorkName.c_str());
490 :
491 16 : for (int iLine = 0; papszHdrLines && papszHdrLines[iLine] != nullptr;
492 : iLine++)
493 : {
494 14 : char **papszTokens = CSLTokenizeString(papszHdrLines[iLine]);
495 :
496 : /* Note: some cv580 file seem to have comments with #, hence the >=
497 : * instead of = for token checking, and the equalN for the corner.
498 : */
499 :
500 14 : if (CSLCount(papszTokens) < 2)
501 : {
502 : /* ignore */;
503 : }
504 14 : else if ((CSLCount(papszTokens) >= 3) &&
505 14 : EQUAL(papszTokens[0], "reference") &&
506 0 : EQUAL(papszTokens[1], "north"))
507 : {
508 0 : dfnorth = CPLAtof(papszTokens[2]);
509 0 : iUTMParamsFound++;
510 : }
511 14 : else if ((CSLCount(papszTokens) >= 3) &&
512 14 : EQUAL(papszTokens[0], "reference") &&
513 0 : EQUAL(papszTokens[1], "east"))
514 : {
515 0 : dfeast = CPLAtof(papszTokens[2]);
516 0 : iUTMParamsFound++;
517 : }
518 14 : else if ((CSLCount(papszTokens) >= 5) &&
519 0 : EQUAL(papszTokens[0], "reference") &&
520 0 : EQUAL(papszTokens[1], "projection") &&
521 14 : EQUAL(papszTokens[2], "UTM") && EQUAL(papszTokens[3], "zone"))
522 : {
523 0 : iUTMZone = atoi(papszTokens[4]);
524 0 : iUTMParamsFound++;
525 : }
526 14 : else if ((CSLCount(papszTokens) >= 3) &&
527 0 : EQUAL(papszTokens[0], "reference") &&
528 14 : EQUAL(papszTokens[1], "corner") &&
529 0 : STARTS_WITH_CI(papszTokens[2], "Upper_Left"))
530 : {
531 : /* iCorner = 0; */
532 0 : iUTMParamsFound++;
533 : }
534 14 : else if (EQUAL(papszTokens[0], "number_lines"))
535 2 : nLines = atoi(papszTokens[1]);
536 :
537 12 : else if (EQUAL(papszTokens[0], "number_samples"))
538 2 : nSamples = atoi(papszTokens[1]);
539 :
540 10 : else if ((EQUAL(papszTokens[0], "header_offset") &&
541 0 : atoi(papszTokens[1]) != 0) ||
542 10 : (EQUAL(papszTokens[0], "number_channels") &&
543 0 : (atoi(papszTokens[1]) != 1) &&
544 0 : (atoi(papszTokens[1]) != 10)) ||
545 10 : (EQUAL(papszTokens[0], "datatype") &&
546 0 : atoi(papszTokens[1]) != 1) ||
547 10 : (EQUAL(papszTokens[0], "number_format") &&
548 0 : !EQUAL(papszTokens[1], "float32") &&
549 0 : !EQUAL(papszTokens[1], "int8")))
550 : {
551 0 : CPLError(CE_Failure, CPLE_AppDefined,
552 : "Keyword %s has value %s which does not match CPG driver "
553 : "expectation.",
554 0 : papszTokens[0], papszTokens[1]);
555 0 : nError = 1;
556 : }
557 10 : else if (EQUAL(papszTokens[0], "altitude"))
558 : {
559 2 : dfaltitude = CPLAtof(papszTokens[1]);
560 2 : iGeoParamsFound++;
561 : }
562 8 : else if (EQUAL(papszTokens[0], "near_srd"))
563 : {
564 2 : dfnear_srd = CPLAtof(papszTokens[1]);
565 2 : iGeoParamsFound++;
566 : }
567 :
568 6 : else if (EQUAL(papszTokens[0], "sample_size"))
569 : {
570 2 : dfsample_size = CPLAtof(papszTokens[1]);
571 2 : iGeoParamsFound++;
572 2 : iUTMParamsFound++;
573 : }
574 4 : else if (EQUAL(papszTokens[0], "sample_size_az"))
575 : {
576 2 : dfsample_size_az = CPLAtof(papszTokens[1]);
577 2 : iGeoParamsFound++;
578 2 : iUTMParamsFound++;
579 : }
580 2 : else if (EQUAL(papszTokens[0], "transposed"))
581 : {
582 2 : itransposed = atoi(papszTokens[1]);
583 2 : iGeoParamsFound++;
584 2 : iUTMParamsFound++;
585 : }
586 :
587 14 : CSLDestroy(papszTokens);
588 : }
589 2 : CSLDestroy(papszHdrLines);
590 : /* -------------------------------------------------------------------- */
591 : /* Check for successful completion. */
592 : /* -------------------------------------------------------------------- */
593 2 : if (nError)
594 : {
595 0 : return nullptr;
596 : }
597 :
598 2 : if (nLines <= 0 || nSamples <= 0)
599 : {
600 0 : CPLError(
601 : CE_Failure, CPLE_AppDefined,
602 : "Did not find valid number_lines or number_samples keywords in %s.",
603 : osWorkName.c_str());
604 0 : return nullptr;
605 : }
606 :
607 : /* -------------------------------------------------------------------- */
608 : /* Initialize dataset. */
609 : /* -------------------------------------------------------------------- */
610 4 : auto poDS = std::make_unique<CPGDataset>();
611 :
612 2 : poDS->nRasterXSize = nSamples;
613 2 : poDS->nRasterYSize = nLines;
614 :
615 : /* -------------------------------------------------------------------- */
616 : /* Open the four bands. */
617 : /* -------------------------------------------------------------------- */
618 : static const char *const apszPolarizations[NUMBER_OF_BANDS] = {"hh", "hv",
619 : "vv", "vh"};
620 :
621 2 : const int nNameLen = static_cast<int>(osWorkName.size());
622 :
623 2 : if (EQUAL(osWorkName.c_str() + nNameLen - 7, "IRC.hdr") ||
624 0 : EQUAL(osWorkName.c_str() + nNameLen - 7, "IRC.img"))
625 : {
626 : {
627 2 : char *pszWorkname = CPLStrdup(osWorkName.c_str());
628 2 : AdjustFilename(&pszWorkname, "", "img");
629 2 : osWorkName = pszWorkname;
630 2 : CPLFree(pszWorkname);
631 : }
632 :
633 2 : poDS->afpImage[0] = VSIFOpenL(osWorkName.c_str(), "rb");
634 2 : if (poDS->afpImage[0] == nullptr)
635 : {
636 0 : CPLError(CE_Failure, CPLE_OpenFailed,
637 : "Failed to open .img file: %s", osWorkName.c_str());
638 0 : return nullptr;
639 : }
640 2 : poDS->aosImageFilenames.push_back(osWorkName.c_str());
641 10 : for (int iBand = 0; iBand < NUMBER_OF_BANDS; iBand++)
642 : {
643 : SIRC_QSLCRasterBand *poBand =
644 8 : new SIRC_QSLCRasterBand(poDS.get(), iBand + 1, GDT_CFloat32);
645 8 : poDS->SetBand(iBand + 1, poBand);
646 8 : poBand->SetMetadataItem("POLARIMETRIC_INTERP",
647 8 : apszPolarizations[iBand]);
648 : }
649 : }
650 : else
651 : {
652 0 : constexpr int SIZEOF_CFLOAT32 = 2 * static_cast<int>(sizeof(float));
653 0 : if (nSamples > INT_MAX / SIZEOF_CFLOAT32)
654 : {
655 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too large nBlockXSize");
656 0 : return nullptr;
657 : }
658 :
659 0 : CPLAssert(poDS->afpImage.size() == NUMBER_OF_BANDS);
660 0 : for (int iBand = 0; iBand < NUMBER_OF_BANDS; iBand++)
661 : {
662 : {
663 0 : char *pszWorkname = CPLStrdup(osWorkName.c_str());
664 0 : AdjustFilename(&pszWorkname, apszPolarizations[iBand], "img");
665 0 : osWorkName = pszWorkname;
666 0 : CPLFree(pszWorkname);
667 : }
668 :
669 0 : poDS->afpImage[iBand] = VSIFOpenL(osWorkName.c_str(), "rb");
670 0 : if (poDS->afpImage[iBand] == nullptr)
671 : {
672 0 : CPLError(CE_Failure, CPLE_OpenFailed,
673 : "Failed to open .img file: %s", osWorkName.c_str());
674 0 : return nullptr;
675 : }
676 0 : poDS->aosImageFilenames.push_back(osWorkName.c_str());
677 :
678 : auto poBand = RawRasterBand::Create(
679 0 : poDS.get(), iBand + 1, poDS->afpImage[iBand], 0, 8,
680 : SIZEOF_CFLOAT32 * nSamples, GDT_CFloat32,
681 : RawRasterBand::ByteOrder::ORDER_BIG_ENDIAN,
682 0 : RawRasterBand::OwnFP::NO);
683 0 : if (!poBand)
684 0 : return nullptr;
685 0 : poBand->SetMetadataItem("POLARIMETRIC_INTERP",
686 0 : apszPolarizations[iBand]);
687 0 : poDS->SetBand(iBand + 1, std::move(poBand));
688 : }
689 : }
690 :
691 : /* Set an appropriate matrix representation metadata item for the set */
692 2 : poDS->SetMetadataItem("MATRIX_REPRESENTATION", "SCATTERING");
693 :
694 : /* -------------------------------------------------------------------------
695 : */
696 : /* Add georeferencing or pseudo-geocoding, if enough information found. */
697 : /* -------------------------------------------------------------------------
698 : */
699 2 : if (iUTMParamsFound == 7)
700 : {
701 0 : poDS->m_gt[1] = 0.0;
702 0 : poDS->m_gt[2] = 0.0;
703 0 : poDS->m_gt[4] = 0.0;
704 0 : poDS->m_gt[5] = 0.0;
705 :
706 : double dfnorth_center;
707 0 : if (itransposed == 1)
708 : {
709 0 : CPLError(CE_Warning, CPLE_AppDefined,
710 : "Did not have a convair SIRC-style test dataset\n"
711 : "with transposed=1 for testing. Georeferencing may be "
712 : "wrong.\n");
713 0 : dfnorth_center = dfnorth - nSamples * dfsample_size / 2.0;
714 0 : poDS->m_gt[0] = dfeast;
715 0 : poDS->m_gt[2] = dfsample_size_az;
716 0 : poDS->m_gt[3] = dfnorth;
717 0 : poDS->m_gt[4] = -1 * dfsample_size;
718 : }
719 : else
720 : {
721 0 : dfnorth_center = dfnorth - nLines * dfsample_size / 2.0;
722 0 : poDS->m_gt[0] = dfeast;
723 0 : poDS->m_gt[1] = dfsample_size_az;
724 0 : poDS->m_gt[3] = dfnorth;
725 0 : poDS->m_gt[5] = -1 * dfsample_size;
726 : }
727 :
728 0 : if (dfnorth_center < 0)
729 0 : poDS->m_oSRS.SetUTM(iUTMZone, 0);
730 : else
731 0 : poDS->m_oSRS.SetUTM(iUTMZone, 1);
732 :
733 : /* Assuming WGS84 */
734 0 : poDS->m_oSRS.SetWellKnownGeogCS("WGS84");
735 : }
736 2 : else if (iGeoParamsFound == 5)
737 : {
738 : double dfgcpLine, dfgcpPixel, dfgcpX, dfgcpY, dftemp;
739 :
740 2 : poDS->nGCPCount = 16;
741 4 : poDS->pasGCPList = reinterpret_cast<GDAL_GCP *>(
742 2 : CPLCalloc(sizeof(GDAL_GCP), poDS->nGCPCount));
743 2 : GDALInitGCPs(poDS->nGCPCount, poDS->pasGCPList);
744 :
745 34 : for (int ngcp = 0; ngcp < 16; ngcp++)
746 : {
747 : char szID[32];
748 :
749 32 : snprintf(szID, sizeof(szID), "%d", ngcp + 1);
750 32 : if (itransposed == 1)
751 : {
752 0 : if (ngcp < 4)
753 0 : dfgcpPixel = 0.0;
754 0 : else if (ngcp < 8)
755 0 : dfgcpPixel = nSamples / 3.0;
756 0 : else if (ngcp < 12)
757 0 : dfgcpPixel = 2.0 * nSamples / 3.0;
758 : else
759 0 : dfgcpPixel = nSamples;
760 :
761 0 : dfgcpLine = nLines * (ngcp % 4) / 3.0;
762 :
763 0 : dftemp = dfnear_srd + (dfsample_size * dfgcpLine);
764 : /* -1 so that 0,0 maps to largest Y */
765 0 : dfgcpY = -1 * sqrt(dftemp * dftemp - dfaltitude * dfaltitude);
766 0 : dfgcpX = dfgcpPixel * dfsample_size_az;
767 : }
768 : else
769 : {
770 32 : if (ngcp < 4)
771 8 : dfgcpLine = 0.0;
772 24 : else if (ngcp < 8)
773 8 : dfgcpLine = nLines / 3.0;
774 16 : else if (ngcp < 12)
775 8 : dfgcpLine = 2.0 * nLines / 3.0;
776 : else
777 8 : dfgcpLine = nLines;
778 :
779 32 : dfgcpPixel = nSamples * ((ngcp % 4) / 3.0);
780 :
781 32 : dftemp = dfnear_srd + (dfsample_size * dfgcpPixel);
782 32 : dfgcpX = sqrt(dftemp * dftemp - dfaltitude * dfaltitude);
783 32 : dfgcpY = (nLines - dfgcpLine) * dfsample_size_az;
784 : }
785 32 : poDS->pasGCPList[ngcp].dfGCPX = dfgcpX;
786 32 : poDS->pasGCPList[ngcp].dfGCPY = dfgcpY;
787 32 : poDS->pasGCPList[ngcp].dfGCPZ = 0.0;
788 :
789 32 : poDS->pasGCPList[ngcp].dfGCPPixel = dfgcpPixel;
790 32 : poDS->pasGCPList[ngcp].dfGCPLine = dfgcpLine;
791 :
792 32 : CPLFree(poDS->pasGCPList[ngcp].pszId);
793 32 : poDS->pasGCPList[ngcp].pszId = CPLStrdup(szID);
794 : }
795 :
796 2 : poDS->m_oGCPSRS.importFromWkt(
797 : "LOCAL_CS[\"Ground range view / unreferenced meters\","
798 : "UNIT[\"Meter\",1.0]]");
799 : }
800 :
801 2 : return poDS.release();
802 : }
803 :
804 : #ifdef notdef
805 : GDALDataset *CPGDataset::InitializeType3Dataset(const char *pszFilename)
806 : {
807 : int iBytesPerPixel = 0;
808 : Interleave::eInterleave = Interleave::BSQ;
809 : bool bInterleaveSpecified = false;
810 : int nLines = 0;
811 : int nSamples = 0;
812 : int nBands = 0;
813 : int nError = 0;
814 :
815 : /* Parameters in geogratis geocoded images */
816 : int iUTMParamsFound = 0;
817 : int iUTMZone = 0;
818 : double dfnorth = 0.0;
819 : double dfeast = 0.0;
820 : double dfOffsetX = 0.0;
821 : double dfOffsetY = 0.0;
822 : double dfxsize = 0.0;
823 : double dfysize = 0.0;
824 :
825 : char *pszWorkname = CPLStrdup(pszFilename);
826 : AdjustFilename(&pszWorkname, "stokes", "img_def");
827 : char **papszHdrLines = CSLLoad(pszWorkname);
828 :
829 : for (int iLine = 0; papszHdrLines && papszHdrLines[iLine] != NULL; iLine++)
830 : {
831 : char **papszTokens =
832 : CSLTokenizeString2(papszHdrLines[iLine], " \t",
833 : CSLT_HONOURSTRINGS & CSLT_ALLOWEMPTYTOKENS);
834 :
835 : /* Note: some cv580 file seem to have comments with #, hence the >=
836 : * instead of = for token checking, and the equalN for the corner.
837 : */
838 :
839 : if ((CSLCount(papszTokens) >= 3) && EQUAL(papszTokens[0], "data") &&
840 : EQUAL(papszTokens[1], "organization:"))
841 : {
842 :
843 : if (STARTS_WITH_CI(papszTokens[2], "BSQ"))
844 : {
845 : bInterleaveSpecified = true;
846 : eInterleave = Interleave::BSQ;
847 : }
848 : else if (STARTS_WITH_CI(papszTokens[2], "BIL"))
849 : {
850 : bInterleaveSpecified = true;
851 : eInterleave = Interleave::BIL;
852 : }
853 : else if (STARTS_WITH_CI(papszTokens[2], "BIP"))
854 : {
855 : bInterleaveSpecified = true;
856 : eInterleave = Interleave::BIP;
857 : }
858 : else
859 : {
860 : CPLError(
861 : CE_Failure, CPLE_AppDefined,
862 : "The interleaving type of the file (%s) is not supported.",
863 : papszTokens[2]);
864 : nError = 1;
865 : }
866 : }
867 : else if ((CSLCount(papszTokens) >= 3) &&
868 : EQUAL(papszTokens[0], "data") &&
869 : EQUAL(papszTokens[1], "state:"))
870 : {
871 :
872 : if (!STARTS_WITH_CI(papszTokens[2], "RAW") &&
873 : !STARTS_WITH_CI(papszTokens[2], "GEO"))
874 : {
875 : CPLError(CE_Failure, CPLE_AppDefined,
876 : "The data state of the file (%s) is not "
877 : "supported.\n. Only RAW and GEO are currently "
878 : "recognized.",
879 : papszTokens[2]);
880 : nError = 1;
881 : }
882 : }
883 : else if ((CSLCount(papszTokens) >= 4) &&
884 : EQUAL(papszTokens[0], "data") &&
885 : EQUAL(papszTokens[1], "origin") &&
886 : EQUAL(papszTokens[2], "point:"))
887 : {
888 : if (!STARTS_WITH_CI(papszTokens[3], "Upper_Left"))
889 : {
890 : CPLError(CE_Failure, CPLE_AppDefined,
891 : "Unexpected value (%s) for data origin point- expect "
892 : "Upper_Left.",
893 : papszTokens[3]);
894 : nError = 1;
895 : }
896 : iUTMParamsFound++;
897 : }
898 : else if ((CSLCount(papszTokens) >= 5) && EQUAL(papszTokens[0], "map") &&
899 : EQUAL(papszTokens[1], "projection:") &&
900 : EQUAL(papszTokens[2], "UTM") && EQUAL(papszTokens[3], "zone"))
901 : {
902 : iUTMZone = atoi(papszTokens[4]);
903 : iUTMParamsFound++;
904 : }
905 : else if ((CSLCount(papszTokens) >= 4) &&
906 : EQUAL(papszTokens[0], "project") &&
907 : EQUAL(papszTokens[1], "origin:"))
908 : {
909 : dfeast = CPLAtof(papszTokens[2]);
910 : dfnorth = CPLAtof(papszTokens[3]);
911 : iUTMParamsFound += 2;
912 : }
913 : else if ((CSLCount(papszTokens) >= 4) &&
914 : EQUAL(papszTokens[0], "file") &&
915 : EQUAL(papszTokens[1], "start:"))
916 : {
917 : dfOffsetX = CPLAtof(papszTokens[2]);
918 : dfOffsetY = CPLAtof(papszTokens[3]);
919 : iUTMParamsFound += 2;
920 : }
921 : else if ((CSLCount(papszTokens) >= 6) &&
922 : EQUAL(papszTokens[0], "pixel") &&
923 : EQUAL(papszTokens[1], "size") && EQUAL(papszTokens[2], "on") &&
924 : EQUAL(papszTokens[3], "ground:"))
925 : {
926 : dfxsize = CPLAtof(papszTokens[4]);
927 : dfysize = CPLAtof(papszTokens[5]);
928 : iUTMParamsFound += 2;
929 : }
930 : else if ((CSLCount(papszTokens) >= 4) &&
931 : EQUAL(papszTokens[0], "number") &&
932 : EQUAL(papszTokens[1], "of") &&
933 : EQUAL(papszTokens[2], "pixels:"))
934 : {
935 : nSamples = atoi(papszTokens[3]);
936 : }
937 : else if ((CSLCount(papszTokens) >= 4) &&
938 : EQUAL(papszTokens[0], "number") &&
939 : EQUAL(papszTokens[1], "of") && EQUAL(papszTokens[2], "lines:"))
940 : {
941 : nLines = atoi(papszTokens[3]);
942 : }
943 : else if ((CSLCount(papszTokens) >= 4) &&
944 : EQUAL(papszTokens[0], "number") &&
945 : EQUAL(papszTokens[1], "of") && EQUAL(papszTokens[2], "bands:"))
946 : {
947 : nBands = atoi(papszTokens[3]);
948 : if (nBands != 16)
949 : {
950 : CPLError(CE_Failure, CPLE_AppDefined,
951 : "Number of bands has a value %s which does not match "
952 : "CPG driver\nexpectation (expect a value of 16).",
953 : papszTokens[3]);
954 : nError = 1;
955 : }
956 : }
957 : else if ((CSLCount(papszTokens) >= 4) &&
958 : EQUAL(papszTokens[0], "bytes") &&
959 : EQUAL(papszTokens[1], "per") &&
960 : EQUAL(papszTokens[2], "pixel:"))
961 : {
962 : iBytesPerPixel = atoi(papszTokens[3]);
963 : if (iBytesPerPixel != 4)
964 : {
965 : CPLError(CE_Failure, CPLE_AppDefined,
966 : "Bytes per pixel has a value %s which does not match "
967 : "CPG driver\nexpectation (expect a value of 4).",
968 : papszTokens[1]);
969 : nError = 1;
970 : }
971 : }
972 : CSLDestroy(papszTokens);
973 : }
974 :
975 : CSLDestroy(papszHdrLines);
976 :
977 : /* -------------------------------------------------------------------- */
978 : /* Check for successful completion. */
979 : /* -------------------------------------------------------------------- */
980 : if (nError)
981 : {
982 : CPLFree(pszWorkname);
983 : return NULL;
984 : }
985 :
986 : if (!GDALCheckDatasetDimensions(nSamples, nLines) ||
987 : !GDALCheckBandCount(nBands, FALSE) || !bInterleaveSpecified ||
988 : iBytesPerPixel == 0)
989 : {
990 : CPLError(CE_Failure, CPLE_AppDefined,
991 : "%s is missing a required parameter (number of pixels, "
992 : "number of lines,\number of bands, bytes per pixel, or "
993 : "data organization).",
994 : pszWorkname);
995 : CPLFree(pszWorkname);
996 : return NULL;
997 : }
998 :
999 : /* -------------------------------------------------------------------- */
1000 : /* Initialize dataset. */
1001 : /* -------------------------------------------------------------------- */
1002 : CPGDataset *poDS = new CPGDataset();
1003 :
1004 : poDS->nRasterXSize = nSamples;
1005 : poDS->nRasterYSize = nLines;
1006 : poDS->eInterleave = eInterleave;
1007 :
1008 : /* -------------------------------------------------------------------- */
1009 : /* Open the 16 bands. */
1010 : /* -------------------------------------------------------------------- */
1011 :
1012 : AdjustFilename(&pszWorkname, "stokes", "img");
1013 : poDS->afpImage[0] = VSIFOpenL(pszWorkname, "rb");
1014 : if (poDS->afpImage[0] == NULL)
1015 : {
1016 : CPLError(CE_Failure, CPLE_OpenFailed, "Failed to open .img file: %s",
1017 : pszWorkname);
1018 : CPLFree(pszWorkname);
1019 : delete poDS;
1020 : return NULL;
1021 : }
1022 : aosImageFilenames.push_back(pszWorkname);
1023 : for (int iBand = 0; iBand < 16; iBand++)
1024 : {
1025 : CPG_STOKESRasterBand *poBand =
1026 : new CPG_STOKESRasterBand(poDS, GDT_CFloat32, !CPL_IS_LSB);
1027 : poDS->SetBand(iBand + 1, poBand);
1028 : }
1029 :
1030 : /* -------------------------------------------------------------------- */
1031 : /* Set appropriate MATRIX_REPRESENTATION. */
1032 : /* -------------------------------------------------------------------- */
1033 : if (poDS->GetRasterCount() == 6)
1034 : {
1035 : poDS->SetMetadataItem("MATRIX_REPRESENTATION", "COVARIANCE");
1036 : }
1037 :
1038 : /* -------------------------------------------------------------------------
1039 : */
1040 : /* Add georeferencing, if enough information found. */
1041 : /* -------------------------------------------------------------------------
1042 : */
1043 : if (iUTMParamsFound == 8)
1044 : {
1045 : double dfnorth_center = dfnorth - nLines * dfysize / 2.0;
1046 : poDS->m_gt[0] = dfeast + dfOffsetX;
1047 : poDS->m_gt[1] = dfxsize;
1048 : poDS->m_gt[2] = 0.0;
1049 : poDS->m_gt[3] = dfnorth + dfOffsetY;
1050 : poDS->m_gt[4] = 0.0;
1051 : poDS->m_gt[5] = -1 * dfysize;
1052 :
1053 : OGRSpatialReference oUTM;
1054 : if (dfnorth_center < 0)
1055 : oUTM.SetUTM(iUTMZone, 0);
1056 : else
1057 : oUTM.SetUTM(iUTMZone, 1);
1058 :
1059 : /* Assuming WGS84 */
1060 : oUTM.SetWellKnownGeogCS("WGS84");
1061 : CPLFree(poDS->pszProjection);
1062 : poDS->pszProjection = NULL;
1063 : oUTM.exportToWkt(&(poDS->pszProjection));
1064 : }
1065 :
1066 : return poDS;
1067 : }
1068 : #endif
1069 :
1070 : /************************************************************************/
1071 : /* Open() */
1072 : /************************************************************************/
1073 :
1074 33518 : GDALDataset *CPGDataset::Open(GDALOpenInfo *poOpenInfo)
1075 :
1076 : {
1077 : /* -------------------------------------------------------------------- */
1078 : /* Is this a PolGASP fileset? We expect fileset to follow */
1079 : /* one of these patterns: */
1080 : /* 1) <stuff>hh<stuff2>.img, <stuff>hh<stuff2>.hdr, */
1081 : /* <stuff>hv<stuff2>.img, <stuff>hv<stuff2>.hdr, */
1082 : /* <stuff>vh<stuff2>.img, <stuff>vh<stuff2>.hdr, */
1083 : /* <stuff>vv<stuff2>.img, <stuff>vv<stuff2>.hdr, */
1084 : /* where <stuff> or <stuff2> should contain the */
1085 : /* substring "sso" or "polgasp" */
1086 : /* 2) <stuff>SIRC.hdr and <stuff>SIRC.img */
1087 : /* 3) <stuff>.img and <stuff>.img_def */
1088 : /* where <stuff> should contain the */
1089 : /* substring "sso" or "polgasp" */
1090 : /* -------------------------------------------------------------------- */
1091 33518 : int CPGType = 0;
1092 33518 : if (FindType1(poOpenInfo->pszFilename))
1093 0 : CPGType = 1;
1094 33514 : else if (FindType2(poOpenInfo->pszFilename))
1095 2 : CPGType = 2;
1096 :
1097 : /* Stokes matrix convair data: not quite working yet- something
1098 : * is wrong in the interpretation of the matrix elements in terms
1099 : * of hh, hv, vv, vh. Data will load if the next two lines are
1100 : * uncommented, but values will be incorrect. Expect C11 = hh*conj(hh),
1101 : * C12 = hh*conj(hv), etc. Used geogratis data in both scattering
1102 : * matrix and stokes format for comparison.
1103 : */
1104 : // else if ( FindType3( poOpenInfo->pszFilename ))
1105 : // CPGType = 3;
1106 :
1107 : /* Set working name back to original */
1108 :
1109 33516 : if (CPGType == 0)
1110 : {
1111 33510 : int nNameLen = static_cast<int>(strlen(poOpenInfo->pszFilename));
1112 33510 : if ((nNameLen > 8) &&
1113 32809 : ((strstr(poOpenInfo->pszFilename, "sso") != nullptr) ||
1114 32809 : (strstr(poOpenInfo->pszFilename, "polgasp") != nullptr)) &&
1115 0 : (EQUAL(poOpenInfo->pszFilename + nNameLen - 4, "img") ||
1116 0 : EQUAL(poOpenInfo->pszFilename + nNameLen - 4, "hdr") ||
1117 0 : EQUAL(poOpenInfo->pszFilename + nNameLen - 7, "img_def")))
1118 : {
1119 0 : CPLError(
1120 : CE_Failure, CPLE_OpenFailed,
1121 : "Apparent attempt to open Convair PolGASP data failed as\n"
1122 : "one or more of the required files is missing (eight files\n"
1123 : "are expected for scattering matrix format, two for Stokes).");
1124 : }
1125 33510 : else if ((nNameLen > 8) &&
1126 32812 : (strstr(poOpenInfo->pszFilename, "SIRC") != nullptr) &&
1127 0 : (EQUAL(poOpenInfo->pszFilename + nNameLen - 4, "img") ||
1128 0 : EQUAL(poOpenInfo->pszFilename + nNameLen - 4, "hdr")))
1129 : {
1130 0 : CPLError(
1131 : CE_Failure, CPLE_OpenFailed,
1132 : "Apparent attempt to open SIRC Convair PolGASP data failed \n"
1133 : "as one of the expected files is missing (hdr or img)!");
1134 : }
1135 33515 : return nullptr;
1136 : }
1137 :
1138 : /* -------------------------------------------------------------------- */
1139 : /* Confirm the requested access is supported. */
1140 : /* -------------------------------------------------------------------- */
1141 6 : if (poOpenInfo->eAccess == GA_Update)
1142 : {
1143 0 : ReportUpdateNotSupportedByDriver("CPG");
1144 0 : return nullptr;
1145 : }
1146 :
1147 : /* Read the header info and create the dataset */
1148 : #ifdef notdef
1149 : if (CPGType < 3)
1150 : #endif
1151 : CPGDataset *poDS = reinterpret_cast<CPGDataset *>(
1152 6 : InitializeType1Or2Dataset(poOpenInfo->pszFilename));
1153 : #ifdef notdef
1154 : else
1155 : poDS = reinterpret_cast<CPGDataset *>(
1156 : InitializeType3Dataset(poOpenInfo->pszFilename));
1157 : #endif
1158 2 : if (poDS == nullptr)
1159 0 : return nullptr;
1160 :
1161 : /* -------------------------------------------------------------------- */
1162 : /* Check for overviews. */
1163 : /* -------------------------------------------------------------------- */
1164 : // Need to think about this.
1165 : // poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
1166 :
1167 : /* -------------------------------------------------------------------- */
1168 : /* Initialize any PAM information. */
1169 : /* -------------------------------------------------------------------- */
1170 2 : poDS->SetDescription(poOpenInfo->pszFilename);
1171 2 : poDS->TryLoadXML();
1172 :
1173 2 : return poDS;
1174 : }
1175 :
1176 : /************************************************************************/
1177 : /* GetGCPCount() */
1178 : /************************************************************************/
1179 :
1180 0 : int CPGDataset::GetGCPCount()
1181 :
1182 : {
1183 0 : return nGCPCount;
1184 : }
1185 :
1186 : /************************************************************************/
1187 : /* GetGCPs() */
1188 : /************************************************************************/
1189 :
1190 0 : const GDAL_GCP *CPGDataset::GetGCPs()
1191 :
1192 : {
1193 0 : return pasGCPList;
1194 : }
1195 :
1196 : /************************************************************************/
1197 : /* GetGeoTransform() */
1198 : /************************************************************************/
1199 :
1200 0 : CPLErr CPGDataset::GetGeoTransform(GDALGeoTransform >) const
1201 :
1202 : {
1203 0 : gt = m_gt;
1204 0 : return CE_None;
1205 : }
1206 :
1207 : /************************************************************************/
1208 : /* SIRC_QSLCRasterBand() */
1209 : /************************************************************************/
1210 :
1211 8 : SIRC_QSLCRasterBand::SIRC_QSLCRasterBand(CPGDataset *poGDSIn, int nBandIn,
1212 8 : GDALDataType eType)
1213 :
1214 : {
1215 8 : poDS = poGDSIn;
1216 8 : nBand = nBandIn;
1217 :
1218 8 : eDataType = eType;
1219 :
1220 8 : nBlockXSize = poGDSIn->nRasterXSize;
1221 8 : nBlockYSize = 1;
1222 :
1223 8 : if (nBand == 1)
1224 2 : SetMetadataItem("POLARIMETRIC_INTERP", "HH");
1225 6 : else if (nBand == 2)
1226 2 : SetMetadataItem("POLARIMETRIC_INTERP", "HV");
1227 4 : else if (nBand == 3)
1228 2 : SetMetadataItem("POLARIMETRIC_INTERP", "VH");
1229 2 : else if (nBand == 4)
1230 2 : SetMetadataItem("POLARIMETRIC_INTERP", "VV");
1231 8 : }
1232 :
1233 : /************************************************************************/
1234 : /* IReadBlock() */
1235 : /************************************************************************/
1236 :
1237 : /* From: http://southport.jpl.nasa.gov/software/dcomp/dcomp.html
1238 :
1239 : ysca = sqrt{ [ (Byte(2) / 254 ) + 1.5] 2Byte(1) }
1240 : Re(SHH) = byte(3) ysca/127
1241 : Im(SHH) = byte(4) ysca/127
1242 : Re(SHV) = byte(5) ysca/127
1243 : Im(SHV) = byte(6) ysca/127
1244 : Re(SVH) = byte(7) ysca/127
1245 : Im(SVH) = byte(8) ysca/127
1246 : Re(SVV) = byte(9) ysca/127
1247 : Im(SVV) = byte(10) ysca/127
1248 : */
1249 :
1250 1 : CPLErr SIRC_QSLCRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff,
1251 : int nBlockYOff, void *pImage)
1252 : {
1253 1 : const int nBytesPerSample = 10;
1254 1 : CPGDataset *poGDS = cpl::down_cast<CPGDataset *>(poDS);
1255 1 : const vsi_l_offset offset =
1256 1 : static_cast<vsi_l_offset>(nBlockXSize) * nBlockYOff * nBytesPerSample;
1257 :
1258 : /* -------------------------------------------------------------------- */
1259 : /* Load all the pixel data associated with this scanline. */
1260 : /* -------------------------------------------------------------------- */
1261 1 : if (nBytesPerSample > INT_MAX / nBlockXSize)
1262 : {
1263 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too large nBlockXSize");
1264 0 : return CE_Failure;
1265 : }
1266 1 : const int nBytesToRead = nBytesPerSample * nBlockXSize;
1267 :
1268 : signed char *pabyRecord =
1269 1 : static_cast<signed char *>(VSI_MALLOC_VERBOSE(nBytesToRead));
1270 1 : if (!pabyRecord)
1271 0 : return CE_Failure;
1272 :
1273 2 : if (VSIFSeekL(poGDS->afpImage[0], offset, SEEK_SET) != 0 ||
1274 1 : static_cast<int>(VSIFReadL(pabyRecord, 1, nBytesToRead,
1275 2 : poGDS->afpImage[0])) != nBytesToRead)
1276 : {
1277 0 : CPLError(CE_Failure, CPLE_FileIO,
1278 : "Error reading %d bytes of SIRC Convair at offset %" PRIu64
1279 : ".\n"
1280 : "Reading file %s failed.",
1281 : nBytesToRead, static_cast<uint64_t>(offset),
1282 0 : poGDS->GetDescription());
1283 0 : CPLFree(pabyRecord);
1284 0 : return CE_Failure;
1285 : }
1286 :
1287 : /* -------------------------------------------------------------------- */
1288 : /* Initialize our power table if this is our first time through. */
1289 : /* -------------------------------------------------------------------- */
1290 : static float afPowTable[256];
1291 1 : [[maybe_unused]] static bool bPowTableInitialized = []()
1292 : {
1293 257 : for (int i = 0; i < 256; i++)
1294 : {
1295 256 : afPowTable[i] = static_cast<float>(pow(2.0, i - 128));
1296 : }
1297 1 : return true;
1298 1 : }();
1299 :
1300 : /* -------------------------------------------------------------------- */
1301 : /* Copy the desired band out based on the size of the type, and */
1302 : /* the interleaving mode. */
1303 : /* -------------------------------------------------------------------- */
1304 2 : for (int iX = 0; iX < nBlockXSize; iX++)
1305 : {
1306 : /* A ones based alias */
1307 1 : const signed char *pabyIn = pabyRecord + iX * nBytesPerSample - 1;
1308 :
1309 : /* coverity[tainted_data] */
1310 : const float fScale = static_cast<float>(
1311 1 : sqrt((static_cast<double>(pabyIn[2]) / 254 + 1.5) *
1312 1 : afPowTable[pabyIn[1] + 128]) /
1313 1 : 127.0);
1314 :
1315 1 : float *pafImage = static_cast<float *>(pImage);
1316 :
1317 1 : if (nBand == 1)
1318 : {
1319 1 : const float fReSHH = static_cast<float>(pabyIn[3] * fScale);
1320 1 : const float fImSHH = static_cast<float>(pabyIn[4] * fScale);
1321 :
1322 1 : pafImage[iX * 2] = fReSHH;
1323 1 : pafImage[iX * 2 + 1] = fImSHH;
1324 : }
1325 0 : else if (nBand == 2)
1326 : {
1327 0 : const float fReSHV = static_cast<float>(pabyIn[5] * fScale);
1328 0 : const float fImSHV = static_cast<float>(pabyIn[6] * fScale);
1329 :
1330 0 : pafImage[iX * 2] = fReSHV;
1331 0 : pafImage[iX * 2 + 1] = fImSHV;
1332 : }
1333 0 : else if (nBand == 3)
1334 : {
1335 0 : const float fReSVH = static_cast<float>(pabyIn[7] * fScale);
1336 0 : const float fImSVH = static_cast<float>(pabyIn[8] * fScale);
1337 :
1338 0 : pafImage[iX * 2] = fReSVH;
1339 0 : pafImage[iX * 2 + 1] = fImSVH;
1340 : }
1341 0 : else if (nBand == 4)
1342 : {
1343 0 : const float fReSVV = static_cast<float>(pabyIn[9] * fScale);
1344 0 : const float fImSVV = static_cast<float>(pabyIn[10] * fScale);
1345 :
1346 0 : pafImage[iX * 2] = fReSVV;
1347 0 : pafImage[iX * 2 + 1] = fImSVV;
1348 : }
1349 : }
1350 :
1351 1 : CPLFree(pabyRecord);
1352 :
1353 1 : return CE_None;
1354 : }
1355 :
1356 : /************************************************************************/
1357 : /* CPG_STOKESRasterBand() */
1358 : /************************************************************************/
1359 :
1360 0 : CPG_STOKESRasterBand::CPG_STOKESRasterBand(GDALDataset *poDSIn,
1361 : GDALDataType eType,
1362 0 : int bNativeOrderIn)
1363 0 : : bNativeOrder(bNativeOrderIn)
1364 : {
1365 : static const char *const apszPolarizations[16] = {
1366 : "Covariance_11", "Covariance_12", "Covariance_13", "Covariance_14",
1367 : "Covariance_21", "Covariance_22", "Covariance_23", "Covariance_24",
1368 : "Covariance_31", "Covariance_32", "Covariance_33", "Covariance_34",
1369 : "Covariance_41", "Covariance_42", "Covariance_43", "Covariance_44"};
1370 :
1371 0 : poDS = poDSIn;
1372 0 : eDataType = eType;
1373 :
1374 0 : nBlockXSize = poDSIn->GetRasterXSize();
1375 0 : nBlockYSize = 1;
1376 :
1377 0 : SetMetadataItem("POLARIMETRIC_INTERP", apszPolarizations[nBand - 1]);
1378 0 : SetDescription(apszPolarizations[nBand - 1]);
1379 0 : }
1380 :
1381 : /************************************************************************/
1382 : /* IReadBlock() */
1383 : /************************************************************************/
1384 :
1385 : /* Convert from Stokes to Covariance representation */
1386 :
1387 0 : CPLErr CPG_STOKESRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff,
1388 : int nBlockYOff, void *pImage)
1389 :
1390 : {
1391 0 : CPLAssert(nBlockXOff == 0);
1392 :
1393 : int m11, /* m12, */ m13, m14, /* m21, */ m22, m23, m24, step;
1394 : int m31, m32, m33, m34, m41, m42, m43, m44;
1395 0 : CPGDataset *poGDS = reinterpret_cast<CPGDataset *>(poDS);
1396 :
1397 0 : CPLErr eErr = poGDS->LoadStokesLine(nBlockYOff, bNativeOrder);
1398 0 : if (eErr != CE_None)
1399 0 : return eErr;
1400 :
1401 0 : float *M = poGDS->padfStokesMatrix;
1402 0 : float *pafLine = reinterpret_cast<float *>(pImage);
1403 :
1404 0 : if (poGDS->eInterleave == RawDataset::Interleave::BIP)
1405 : {
1406 0 : step = 16;
1407 0 : m11 = M11;
1408 : // m12 = M12;
1409 0 : m13 = M13;
1410 0 : m14 = M14;
1411 : // m21 = M21;
1412 0 : m22 = M22;
1413 0 : m23 = M23;
1414 0 : m24 = M24;
1415 0 : m31 = M31;
1416 0 : m32 = M32;
1417 0 : m33 = M33;
1418 0 : m34 = M34;
1419 0 : m41 = M41;
1420 0 : m42 = M42;
1421 0 : m43 = M43;
1422 0 : m44 = M44;
1423 : }
1424 : else
1425 : {
1426 0 : step = 1;
1427 0 : m11 = 0;
1428 : // m12=nRasterXSize;
1429 0 : m13 = nRasterXSize * 2;
1430 0 : m14 = nRasterXSize * 3;
1431 : // m21=nRasterXSize*4;
1432 0 : m22 = nRasterXSize * 5;
1433 0 : m23 = nRasterXSize * 6;
1434 0 : m24 = nRasterXSize * 7;
1435 0 : m31 = nRasterXSize * 8;
1436 0 : m32 = nRasterXSize * 9;
1437 0 : m33 = nRasterXSize * 10;
1438 0 : m34 = nRasterXSize * 11;
1439 0 : m41 = nRasterXSize * 12;
1440 0 : m42 = nRasterXSize * 13;
1441 0 : m43 = nRasterXSize * 14;
1442 0 : m44 = nRasterXSize * 15;
1443 : }
1444 0 : if (nBand == 1) /* C11 */
1445 : {
1446 0 : for (int iPixel = 0; iPixel < nRasterXSize; iPixel++)
1447 : {
1448 0 : pafLine[iPixel * 2 + 0] = M[m11] - M[m22] - M[m33] + M[m44];
1449 0 : pafLine[iPixel * 2 + 1] = 0.0;
1450 0 : m11 += step;
1451 0 : m22 += step;
1452 0 : m33 += step;
1453 0 : m44 += step;
1454 : }
1455 : }
1456 0 : else if (nBand == 2) /* C12 */
1457 : {
1458 0 : for (int iPixel = 0; iPixel < nRasterXSize; iPixel++)
1459 : {
1460 0 : pafLine[iPixel * 2 + 0] = M[m13] - M[m23];
1461 0 : pafLine[iPixel * 2 + 1] = M[m14] - M[m24];
1462 0 : m13 += step;
1463 0 : m23 += step;
1464 0 : m14 += step;
1465 0 : m24 += step;
1466 : }
1467 : }
1468 0 : else if (nBand == 3) /* C13 */
1469 : {
1470 0 : for (int iPixel = 0; iPixel < nRasterXSize; iPixel++)
1471 : {
1472 0 : pafLine[iPixel * 2 + 0] = M[m33] - M[m44];
1473 0 : pafLine[iPixel * 2 + 1] = M[m43] + M[m34];
1474 0 : m33 += step;
1475 0 : m44 += step;
1476 0 : m43 += step;
1477 0 : m34 += step;
1478 : }
1479 : }
1480 0 : else if (nBand == 4) /* C14 */
1481 : {
1482 0 : for (int iPixel = 0; iPixel < nRasterXSize; iPixel++)
1483 : {
1484 0 : pafLine[iPixel * 2 + 0] = M[m31] - M[m32];
1485 0 : pafLine[iPixel * 2 + 1] = M[m41] - M[m42];
1486 0 : m31 += step;
1487 0 : m32 += step;
1488 0 : m41 += step;
1489 0 : m42 += step;
1490 : }
1491 : }
1492 0 : else if (nBand == 5) /* C21 */
1493 : {
1494 0 : for (int iPixel = 0; iPixel < nRasterXSize; iPixel++)
1495 : {
1496 0 : pafLine[iPixel * 2 + 0] = M[m13] - M[m23];
1497 0 : pafLine[iPixel * 2 + 1] = M[m24] - M[m14];
1498 0 : m13 += step;
1499 0 : m23 += step;
1500 0 : m14 += step;
1501 0 : m24 += step;
1502 : }
1503 : }
1504 0 : else if (nBand == 6) /* C22 */
1505 : {
1506 0 : for (int iPixel = 0; iPixel < nRasterXSize; iPixel++)
1507 : {
1508 0 : pafLine[iPixel * 2 + 0] = M[m11] + M[m22] - M[m33] - M[m44];
1509 0 : pafLine[iPixel * 2 + 1] = 0.0;
1510 0 : m11 += step;
1511 0 : m22 += step;
1512 0 : m33 += step;
1513 0 : m44 += step;
1514 : }
1515 : }
1516 0 : else if (nBand == 7) /* C23 */
1517 : {
1518 0 : for (int iPixel = 0; iPixel < nRasterXSize; iPixel++)
1519 : {
1520 0 : pafLine[iPixel * 2 + 0] = M[m31] + M[m32];
1521 0 : pafLine[iPixel * 2 + 1] = M[m41] + M[m42];
1522 0 : m31 += step;
1523 0 : m32 += step;
1524 0 : m41 += step;
1525 0 : m42 += step;
1526 : }
1527 : }
1528 0 : else if (nBand == 8) /* C24 */
1529 : {
1530 0 : for (int iPixel = 0; iPixel < nRasterXSize; iPixel++)
1531 : {
1532 0 : pafLine[iPixel * 2 + 0] = M[m33] + M[m44];
1533 0 : pafLine[iPixel * 2 + 1] = M[m43] - M[m34];
1534 0 : m33 += step;
1535 0 : m44 += step;
1536 0 : m43 += step;
1537 0 : m34 += step;
1538 : }
1539 : }
1540 0 : else if (nBand == 9) /* C31 */
1541 : {
1542 0 : for (int iPixel = 0; iPixel < nRasterXSize; iPixel++)
1543 : {
1544 0 : pafLine[iPixel * 2 + 0] = M[m33] - M[m44];
1545 0 : pafLine[iPixel * 2 + 1] = -1 * M[m43] - M[m34];
1546 0 : m33 += step;
1547 0 : m44 += step;
1548 0 : m43 += step;
1549 0 : m34 += step;
1550 : }
1551 : }
1552 0 : else if (nBand == 10) /* C32 */
1553 : {
1554 0 : for (int iPixel = 0; iPixel < nRasterXSize; iPixel++)
1555 : {
1556 0 : pafLine[iPixel * 2 + 0] = M[m31] + M[m32];
1557 0 : pafLine[iPixel * 2 + 1] = -1 * M[m41] - M[m42];
1558 0 : m31 += step;
1559 0 : m32 += step;
1560 0 : m41 += step;
1561 0 : m42 += step;
1562 : }
1563 : }
1564 0 : else if (nBand == 11) /* C33 */
1565 : {
1566 0 : for (int iPixel = 0; iPixel < nRasterXSize; iPixel++)
1567 : {
1568 0 : pafLine[iPixel * 2 + 0] = M[m11] + M[m22] + M[m33] + M[m44];
1569 0 : pafLine[iPixel * 2 + 1] = 0.0;
1570 0 : m11 += step;
1571 0 : m22 += step;
1572 0 : m33 += step;
1573 0 : m44 += step;
1574 : }
1575 : }
1576 0 : else if (nBand == 12) /* C34 */
1577 : {
1578 0 : for (int iPixel = 0; iPixel < nRasterXSize; iPixel++)
1579 : {
1580 0 : pafLine[iPixel * 2 + 0] = M[m13] - M[m23];
1581 0 : pafLine[iPixel * 2 + 1] = -1 * M[m14] - M[m24];
1582 0 : m13 += step;
1583 0 : m23 += step;
1584 0 : m14 += step;
1585 0 : m24 += step;
1586 : }
1587 : }
1588 0 : else if (nBand == 13) /* C41 */
1589 : {
1590 0 : for (int iPixel = 0; iPixel < nRasterXSize; iPixel++)
1591 : {
1592 0 : pafLine[iPixel * 2 + 0] = M[m31] - M[m32];
1593 0 : pafLine[iPixel * 2 + 1] = M[m42] - M[m41];
1594 0 : m31 += step;
1595 0 : m32 += step;
1596 0 : m41 += step;
1597 0 : m42 += step;
1598 : }
1599 : }
1600 0 : else if (nBand == 14) /* C42 */
1601 : {
1602 0 : for (int iPixel = 0; iPixel < nRasterXSize; iPixel++)
1603 : {
1604 0 : pafLine[iPixel * 2 + 0] = M[m33] + M[m44];
1605 0 : pafLine[iPixel * 2 + 1] = M[m34] - M[m43];
1606 0 : m33 += step;
1607 0 : m44 += step;
1608 0 : m43 += step;
1609 0 : m34 += step;
1610 : }
1611 : }
1612 0 : else if (nBand == 15) /* C43 */
1613 : {
1614 0 : for (int iPixel = 0; iPixel < nRasterXSize; iPixel++)
1615 : {
1616 0 : pafLine[iPixel * 2 + 0] = M[m13] - M[m23];
1617 0 : pafLine[iPixel * 2 + 1] = M[m14] + M[m24];
1618 0 : m13 += step;
1619 0 : m23 += step;
1620 0 : m14 += step;
1621 0 : m24 += step;
1622 : }
1623 : }
1624 : else /* C44 */
1625 : {
1626 0 : for (int iPixel = 0; iPixel < nRasterXSize; iPixel++)
1627 : {
1628 0 : pafLine[iPixel * 2 + 0] = M[m11] - M[m22] + M[m33] - M[m44];
1629 0 : pafLine[iPixel * 2 + 1] = 0.0;
1630 0 : m11 += step;
1631 0 : m22 += step;
1632 0 : m33 += step;
1633 0 : m44 += step;
1634 : }
1635 : }
1636 :
1637 0 : return CE_None;
1638 : }
1639 :
1640 : /************************************************************************/
1641 : /* GDALRegister_CPG() */
1642 : /************************************************************************/
1643 :
1644 1935 : void GDALRegister_CPG()
1645 :
1646 : {
1647 1935 : if (GDALGetDriverByName("CPG") != nullptr)
1648 282 : return;
1649 :
1650 1653 : GDALDriver *poDriver = new GDALDriver();
1651 :
1652 1653 : poDriver->SetDescription("CPG");
1653 1653 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1654 1653 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Convair PolGASP");
1655 1653 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1656 :
1657 1653 : poDriver->pfnOpen = CPGDataset::Open;
1658 :
1659 1653 : GetGDALDriverManager()->RegisterDriver(poDriver);
1660 : }
|
