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