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