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