Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: DRDC Ottawa GEOINT
4 : * Purpose: Radarsat Constellation Mission - XML Products (product.xml) driver
5 : * Author: Roberto Caron, MDA
6 : * on behalf of DRDC Ottawa
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2020, DRDC Ottawa
10 : *
11 : * Based on the RS2 Dataset Class
12 : *
13 : * SPDX-License-Identifier: MIT
14 : ****************************************************************************/
15 :
16 : #include <time.h>
17 : #include <stdio.h>
18 : #include <sstream>
19 :
20 : #include "cpl_minixml.h"
21 : #include "gdal_frmts.h"
22 : #include "gdal_pam.h"
23 : #include "ogr_spatialref.h"
24 : #include "rcmdataset.h"
25 : #include "rcmdrivercore.h"
26 :
27 : #include <limits>
28 :
29 : constexpr int max_space_for_string = 32;
30 :
31 : /* RCM has a special folder that contains all LUT, Incidence Angle and Noise
32 : * Level files */
33 : constexpr const char *CALIBRATION_FOLDER = "calibration";
34 :
35 : /*** Function to format calibration for unique identification for Layer Name
36 : * ***/
37 : /*
38 : * RCM_CALIB : { SIGMA0 | GAMMA0 | BETA0 | UNCALIB } : product.xml full path
39 : */
40 56 : inline CPLString FormatCalibration(const char *pszCalibName,
41 : const char *pszFilename)
42 : {
43 56 : CPLString ptr;
44 :
45 : // Always begin by the layer calibration name
46 56 : ptr.append(szLayerCalibration);
47 :
48 : // A separator is needed before concat calibration name
49 56 : ptr += chLayerSeparator;
50 : // Add calibration name
51 56 : ptr.append(pszCalibName);
52 :
53 : // A separator is needed before concat full filename name
54 56 : ptr += chLayerSeparator;
55 : // Add full filename name
56 56 : ptr.append(pszFilename);
57 :
58 : /* return calibration format */
59 56 : return ptr;
60 : }
61 :
62 : /*** Function to test for valid LUT files ***/
63 72 : static bool IsValidXMLFile(const char *pszPath)
64 : {
65 : /* Return true for valid xml file, false otherwise */
66 72 : CPLXMLTreeCloser psLut(CPLParseXMLFile(pszPath));
67 :
68 72 : if (psLut.get() == nullptr)
69 : {
70 0 : CPLError(CE_Failure, CPLE_OpenFailed,
71 : "ERROR: Failed to open the LUT file %s", pszPath);
72 : }
73 :
74 144 : return psLut.get() != nullptr;
75 : }
76 :
77 14 : static double *InterpolateValues(CSLConstList papszList, int tableSize,
78 : int stepSize, int numberOfValues,
79 : int pixelFirstLutValue)
80 : {
81 : /* Allocate the right LUT size according to the product range pixel */
82 : double *table =
83 14 : static_cast<double *>(VSI_CALLOC_VERBOSE(sizeof(double), tableSize));
84 14 : if (!table)
85 0 : return nullptr;
86 :
87 14 : if (stepSize <= 0)
88 : {
89 : /* When negative, the range of pixel is calculated from the opposite
90 : * starting from the end of gains array */
91 : /* Just step the range with positive value */
92 14 : const int positiveStepSize = abs(stepSize);
93 :
94 14 : int k = 0;
95 :
96 14 : if (positiveStepSize == 1)
97 : {
98 : // Be fast and just copy the values because all gain values
99 : // represent all image wide pixel
100 : /* Start at the end position and store in the opposite */
101 0 : for (int i = pixelFirstLutValue; i >= 0; i--)
102 : {
103 0 : const double value = CPLAtof(papszList[i]);
104 0 : table[k++] = value;
105 : }
106 : }
107 : else
108 : {
109 :
110 : /* Interpolation between 2 numbers */
111 422 : for (int i = numberOfValues - 1; i >= 0; i--)
112 : {
113 : // We will consider the same value to cover the case that we
114 : // will hit the last pixel
115 408 : double valueFrom = CPLAtof(papszList[i]);
116 408 : double valueTo = valueFrom;
117 :
118 408 : if (i > 0)
119 : {
120 : // We have room to pick the previous number to interpolate
121 : // with
122 394 : valueTo = CPLAtof(papszList[i - 1]);
123 : }
124 :
125 : // If the valueFrom minus ValueTo equal 0, it means to finish
126 : // off with the same number until the end of the table size
127 408 : double interp = (valueTo - valueFrom) / positiveStepSize;
128 :
129 : // Always begin with the value FROM found
130 408 : table[k++] = valueFrom;
131 :
132 : // Then add interpolation, don't forget. The stepSize is
133 : // actually counting our valueFrom number thus we add
134 : // interpolation until the last step - 1
135 116341 : for (int j = 0; j < positiveStepSize - 1; j++)
136 : {
137 115933 : valueFrom += interp;
138 115933 : table[k++] = valueFrom;
139 : }
140 : }
141 : }
142 : }
143 : else
144 : {
145 : /* When positive, the range of pixel is calculated from the beginning of
146 : * gains array */
147 0 : if (stepSize == 1)
148 : {
149 : // Be fast and just copy the values because all gain values
150 : // represent all image wide pixel
151 0 : for (int i = 0; i < numberOfValues; i++)
152 : {
153 0 : const double value = CPLAtof(papszList[i]);
154 0 : table[i] = value;
155 : }
156 : }
157 : else
158 : {
159 : /* Interpolation between 2 numbers */
160 0 : int k = 0;
161 0 : for (int i = 0; i < numberOfValues; i++)
162 : {
163 : // We will consider the same value to cover the case that we
164 : // will hit the last pixel
165 0 : double valueFrom = CPLAtof(papszList[i]);
166 0 : double valueTo = valueFrom;
167 :
168 0 : if (i < (numberOfValues)-1)
169 : {
170 : // We have room to pick the next number to interpolate with
171 0 : valueTo = CPLAtof(papszList[i + 1]);
172 : }
173 :
174 : // If the valueFrom minus ValueTo equal 0, it means to finish
175 : // off with the same number until the end of the table size
176 0 : double interp = (valueTo - valueFrom) / stepSize;
177 :
178 : // Always begin with the value FROM found
179 0 : table[k++] = valueFrom;
180 :
181 : // Then add interpolation, don't forget. The stepSize is
182 : // actually counting our valueFrom number thus we add
183 : // interpolation until the last step - 1
184 0 : for (int j = 0; j < stepSize - 1; j++)
185 : {
186 0 : valueFrom += interp;
187 0 : table[k++] = valueFrom;
188 : }
189 : }
190 : }
191 : }
192 :
193 14 : return table;
194 : }
195 :
196 : /*** check that the referenced dataset for each band has the
197 : correct data type and returns whether a 2 band I+Q dataset should
198 : be mapped onto a single complex band.
199 : Returns BANDERROR for error, STRAIGHT for 1:1 mapping, TWOBANDCOMPLEX for 2
200 : bands -> 1 complex band
201 : */
202 : typedef enum
203 : {
204 : BANDERROR,
205 : STRAIGHT,
206 : TWOBANDCOMPLEX
207 : } BandMappingRCM;
208 :
209 16 : static BandMappingRCM checkBandFileMappingRCM(GDALDataType dataType,
210 : GDALDataset *poBandFile,
211 : bool isNITF)
212 : {
213 :
214 16 : GDALRasterBand *band1 = poBandFile->GetRasterBand(1);
215 16 : GDALDataType bandfileType = band1->GetRasterDataType();
216 : // if there is one band and it has the same datatype, the band file gets
217 : // passed straight through
218 16 : if ((poBandFile->GetRasterCount() == 1 ||
219 16 : poBandFile->GetRasterCount() == 4) &&
220 : dataType == bandfileType)
221 16 : return STRAIGHT;
222 :
223 : // if the band file has 2 bands, they should represent I+Q
224 : // and be a compatible data type
225 0 : if (poBandFile->GetRasterCount() == 2 && GDALDataTypeIsComplex(dataType))
226 : {
227 0 : GDALRasterBand *band2 = poBandFile->GetRasterBand(2);
228 :
229 0 : if (bandfileType != band2->GetRasterDataType())
230 0 : return BANDERROR; // both bands must be same datatype
231 :
232 : // check compatible types - there are 4 complex types in GDAL
233 0 : if ((dataType == GDT_CInt16 && bandfileType == GDT_Int16) ||
234 0 : (dataType == GDT_CInt32 && bandfileType == GDT_Int32) ||
235 0 : (dataType == GDT_CFloat32 && bandfileType == GDT_Float32) ||
236 0 : (dataType == GDT_CFloat64 && bandfileType == GDT_Float64))
237 0 : return TWOBANDCOMPLEX;
238 :
239 0 : if ((dataType == GDT_CInt16 && bandfileType == GDT_CInt16) ||
240 0 : (dataType == GDT_CInt32 && bandfileType == GDT_CInt32) ||
241 0 : (dataType == GDT_CFloat32 && bandfileType == GDT_CFloat32) ||
242 0 : (dataType == GDT_CFloat64 && bandfileType == GDT_CFloat64))
243 0 : return TWOBANDCOMPLEX;
244 : }
245 :
246 0 : if (isNITF)
247 : {
248 0 : return STRAIGHT;
249 : }
250 :
251 0 : return BANDERROR; // don't accept any other combinations
252 : }
253 :
254 : /************************************************************************/
255 : /* RCMRasterBand */
256 : /************************************************************************/
257 :
258 10 : RCMRasterBand::RCMRasterBand(RCMDataset *poDSIn, int nBandIn,
259 : GDALDataType eDataTypeIn, const char *pszPole,
260 : GDALDataset *poBandFileIn, bool bTwoBandComplex,
261 10 : bool isOneFilePerPolIn, bool isNITFIn)
262 : : poBandFile(poBandFileIn), poRCMDataset(poDSIn),
263 : twoBandComplex(bTwoBandComplex), isOneFilePerPol(isOneFilePerPolIn),
264 10 : isNITF(isNITFIn)
265 : {
266 10 : poDS = poDSIn;
267 10 : this->nBand = nBandIn;
268 10 : eDataType = eDataTypeIn;
269 :
270 : /*Check image type, whether there is one file per polarization or
271 : *one file containing all polarizations*/
272 10 : if (this->isOneFilePerPol)
273 : {
274 10 : poBand = poBandFile->GetRasterBand(1);
275 : }
276 : else
277 : {
278 0 : poBand = poBandFile->GetRasterBand(this->nBand);
279 : }
280 :
281 10 : poBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
282 :
283 10 : if (pszPole != nullptr && strlen(pszPole) != 0)
284 : {
285 10 : SetMetadataItem("POLARIMETRIC_INTERP", pszPole);
286 : }
287 10 : }
288 :
289 : /************************************************************************/
290 : /* RCMRasterBand() */
291 : /************************************************************************/
292 :
293 20 : RCMRasterBand::~RCMRasterBand()
294 :
295 : {
296 10 : if (poBandFile != nullptr)
297 10 : GDALClose(reinterpret_cast<GDALRasterBandH>(poBandFile));
298 20 : }
299 :
300 : /************************************************************************/
301 : /* IReadBlock() */
302 : /************************************************************************/
303 :
304 13 : CPLErr RCMRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage)
305 :
306 : {
307 13 : int nRequestXSize = 0;
308 13 : int nRequestYSize = 0;
309 13 : GetActualBlockSize(nBlockXOff, nBlockYOff, &nRequestXSize, &nRequestYSize);
310 :
311 : // Zero initial partial right-most and bottom-most blocks
312 13 : if (nRequestXSize < nBlockXSize || nRequestYSize < nBlockYSize)
313 : {
314 1 : memset(pImage, 0,
315 1 : static_cast<size_t>(GDALGetDataTypeSizeBytes(eDataType)) *
316 1 : nBlockXSize * nBlockYSize);
317 : }
318 :
319 13 : int dataTypeSize = GDALGetDataTypeSizeBytes(eDataType);
320 : GDALDataType bandFileType =
321 13 : poBandFile->GetRasterBand(1)->GetRasterDataType();
322 13 : int bandFileSize = GDALGetDataTypeSizeBytes(bandFileType);
323 :
324 : // case: 2 bands representing I+Q -> one complex band
325 13 : if (twoBandComplex && !this->isNITF)
326 : {
327 : // int bandFileSize = GDALGetDataTypeSizeBytes(bandFileType);
328 : // this data type is the complex version of the band file
329 : // Roberto: don't check that for the moment: CPLAssert(dataTypeSize ==
330 : // bandFileSize * 2);
331 :
332 : return
333 : // I and Q from each band are pixel-interleaved into this complex
334 : // band
335 0 : poBandFile->RasterIO(
336 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
337 : nRequestXSize, nRequestYSize, pImage, nRequestXSize,
338 : nRequestYSize, bandFileType, 2, nullptr, dataTypeSize,
339 0 : static_cast<GSpacing>(dataTypeSize) * nBlockXSize, bandFileSize,
340 0 : nullptr);
341 : }
342 13 : else if (twoBandComplex && this->isNITF)
343 : {
344 0 : return poBand->RasterIO(
345 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
346 : nRequestXSize, nRequestYSize, pImage, nRequestXSize, nRequestYSize,
347 0 : eDataType, 0, static_cast<GSpacing>(dataTypeSize) * nBlockXSize,
348 0 : nullptr);
349 : }
350 :
351 13 : if (poRCMDataset->IsComplexData())
352 : {
353 : // this data type is the complex version of the band file
354 : // Roberto: don't check that for the moment: CPLAssert(dataTypeSize ==
355 : // bandFileSize * 2);
356 : return
357 : // I and Q from each band are pixel-interleaved into this complex
358 : // band
359 0 : poBandFile->RasterIO(
360 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
361 : nRequestXSize, nRequestYSize, pImage, nRequestXSize,
362 : nRequestYSize, bandFileType, 2, nullptr, dataTypeSize,
363 0 : static_cast<GSpacing>(dataTypeSize) * nBlockXSize, bandFileSize,
364 0 : nullptr);
365 : }
366 :
367 : // case: band file == this band
368 : // NOTE: if the underlying band is opened with the NITF driver, it may
369 : // combine 2 band I+Q -> complex band
370 13 : else if (poBandFile->GetRasterBand(1)->GetRasterDataType() == eDataType)
371 : {
372 26 : return poBand->RasterIO(
373 13 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
374 : nRequestXSize, nRequestYSize, pImage, nRequestXSize, nRequestYSize,
375 13 : eDataType, 0, static_cast<GSpacing>(dataTypeSize) * nBlockXSize,
376 13 : nullptr);
377 : }
378 : else
379 : {
380 0 : CPLAssert(FALSE);
381 : return CE_Failure;
382 : }
383 : }
384 :
385 : /************************************************************************/
386 : /* ReadLUT() */
387 : /************************************************************************/
388 : /* Read the provided LUT in to m_ndTable */
389 : /* 1. The gains list spans the range extent covered by all */
390 : /* beams(if applicable). */
391 : /* 2. The mapping between the entry of gains */
392 : /* list and the range sample index is : the range sample */
393 : /* index = gains entry index * stepSize + pixelFirstLutValue, */
394 : /* where the gains entry index starts with ‘0’.For ScanSAR SLC, */
395 : /* the range sample index refers to the index on the COPG */
396 : /************************************************************************/
397 6 : void RCMCalibRasterBand::ReadLUT()
398 : {
399 :
400 : char bandNumber[12];
401 6 : snprintf(bandNumber, sizeof(bandNumber), "%d", poDS->GetRasterCount() + 1);
402 :
403 6 : CPLXMLTreeCloser psLUT(CPLParseXMLFile(m_pszLUTFile));
404 6 : if (!psLUT)
405 0 : return;
406 :
407 6 : this->m_nfOffset =
408 6 : CPLAtof(CPLGetXMLValue(psLUT.get(), "=lut.offset", "0.0"));
409 :
410 6 : this->pixelFirstLutValue =
411 6 : atoi(CPLGetXMLValue(psLUT.get(), "=lut.pixelFirstLutValue", "0"));
412 :
413 6 : this->stepSize = atoi(CPLGetXMLValue(psLUT.get(), "=lut.stepSize", "0"));
414 :
415 6 : this->numberOfValues =
416 6 : atoi(CPLGetXMLValue(psLUT.get(), "=lut.numberOfValues", "0"));
417 :
418 6 : if (this->numberOfValues <= 0)
419 : {
420 0 : CPLError(
421 : CE_Failure, CPLE_NotSupported,
422 : "ERROR: The RCM driver does not support the LUT Number Of Values "
423 : "equal or lower than zero.");
424 0 : return;
425 : }
426 :
427 : const CPLStringList aosLUTList(
428 6 : CSLTokenizeString2(CPLGetXMLValue(psLUT.get(), "=lut.gains", ""), " ",
429 6 : CSLT_HONOURSTRINGS));
430 :
431 6 : if (this->stepSize <= 0)
432 : {
433 6 : if (this->pixelFirstLutValue <= 0)
434 : {
435 0 : CPLError(
436 : CE_Failure, CPLE_NotSupported,
437 : "ERROR: The RCM driver does not support LUT Pixel First Lut "
438 : "Value equal or lower than zero when the product is "
439 : "descending.");
440 0 : return;
441 : }
442 : }
443 :
444 : /* Get the Pixel Per range */
445 6 : if (this->stepSize == 0 || this->stepSize == INT_MIN ||
446 6 : this->numberOfValues == INT_MIN ||
447 6 : abs(this->stepSize) > INT_MAX / abs(this->numberOfValues))
448 : {
449 0 : CPLError(CE_Failure, CPLE_AppDefined,
450 : "Bad values of stepSize / numberOfValues");
451 0 : return;
452 : }
453 :
454 6 : this->m_nTableSize = abs(this->stepSize) * abs(this->numberOfValues);
455 :
456 6 : if (this->m_nTableSize < this->m_poBandDataset->GetRasterXSize())
457 : {
458 0 : CPLError(
459 : CE_Failure, CPLE_NotSupported,
460 : "ERROR: The RCM driver does not support range of LUT gain values "
461 : "lower than the full image pixel range.");
462 0 : return;
463 : }
464 :
465 : // Avoid excessive memory allocation
466 6 : if (this->m_nTableSize > 1000 * 1000)
467 : {
468 0 : CPLError(CE_Failure, CPLE_NotSupported, "Too many elements in LUT: %d",
469 : this->m_nTableSize);
470 0 : return;
471 : }
472 :
473 : /* Allocate the right LUT size according to the product range pixel */
474 6 : this->m_nfTable =
475 6 : InterpolateValues(aosLUTList.List(), this->m_nTableSize, this->stepSize,
476 : this->numberOfValues, this->pixelFirstLutValue);
477 6 : if (!this->m_nfTable)
478 0 : return;
479 :
480 : // 32 max + space
481 : char *lut_gains = static_cast<char *>(
482 6 : VSI_CALLOC_VERBOSE(this->m_nTableSize, max_space_for_string));
483 6 : if (!lut_gains)
484 0 : return;
485 :
486 107496 : for (int i = 0; i < this->m_nTableSize; i++)
487 : {
488 : char lut[max_space_for_string];
489 : // 6.123004711900930e+04 %e Scientific annotation
490 107490 : snprintf(lut, sizeof(lut), "%e ", this->m_nfTable[i]);
491 107490 : strcat(lut_gains, lut);
492 : }
493 :
494 6 : poDS->SetMetadataItem(CPLString("LUT_GAINS_").append(bandNumber).c_str(),
495 6 : lut_gains);
496 : // Can free this because the function SetMetadataItem takes a copy
497 6 : CPLFree(lut_gains);
498 :
499 : char snum[256];
500 6 : if (this->m_eCalib == eCalibration::Sigma0)
501 : {
502 2 : poDS->SetMetadataItem(CPLString("LUT_TYPE_").append(bandNumber).c_str(),
503 2 : "SIGMA0");
504 : }
505 4 : else if (this->m_eCalib == eCalibration::Beta0)
506 : {
507 2 : poDS->SetMetadataItem(CPLString("LUT_TYPE_").append(bandNumber).c_str(),
508 2 : "BETA0");
509 : }
510 2 : else if (this->m_eCalib == eCalibration::Gamma)
511 : {
512 2 : poDS->SetMetadataItem(CPLString("LUT_TYPE_").append(bandNumber).c_str(),
513 2 : "GAMMA");
514 : }
515 6 : snprintf(snum, sizeof(snum), "%d", this->m_nTableSize);
516 6 : poDS->SetMetadataItem(CPLString("LUT_SIZE_").append(bandNumber).c_str(),
517 6 : snum);
518 6 : snprintf(snum, sizeof(snum), "%f", this->m_nfOffset);
519 6 : poDS->SetMetadataItem(CPLString("LUT_OFFSET_").append(bandNumber).c_str(),
520 6 : snum);
521 : }
522 :
523 : /************************************************************************/
524 : /* ReadNoiseLevels() */
525 : /************************************************************************/
526 : /* Read the provided LUT in to m_nfTableNoiseLevels */
527 : /* 1. The gains list spans the range extent covered by all */
528 : /* beams(if applicable). */
529 : /* 2. The mapping between the entry of gains */
530 : /* list and the range sample index is : the range sample */
531 : /* index = gains entry index * stepSize + pixelFirstLutValue, */
532 : /* where the gains entry index starts with ‘0’.For ScanSAR SLC, */
533 : /* the range sample index refers to the index on the COPG */
534 : /************************************************************************/
535 6 : void RCMCalibRasterBand::ReadNoiseLevels()
536 : {
537 :
538 6 : this->m_nfTableNoiseLevels = nullptr;
539 :
540 6 : if (this->m_pszNoiseLevelsFile == nullptr)
541 : {
542 0 : return;
543 : }
544 :
545 : char bandNumber[12];
546 6 : snprintf(bandNumber, sizeof(bandNumber), "%d", poDS->GetRasterCount() + 1);
547 :
548 6 : CPLXMLTreeCloser psNoiseLevels(CPLParseXMLFile(this->m_pszNoiseLevelsFile));
549 6 : if (!psNoiseLevels)
550 0 : return;
551 :
552 : // Load Beta Nought, Sigma Nought, Gamma noise levels
553 : // Loop through all nodes with spaces
554 : CPLXMLNode *psReferenceNoiseLevelNode =
555 6 : CPLGetXMLNode(psNoiseLevels.get(), "=noiseLevels");
556 6 : if (!psReferenceNoiseLevelNode)
557 0 : return;
558 :
559 : CPLXMLNode *psNodeInc;
560 234 : for (psNodeInc = psReferenceNoiseLevelNode->psChild; psNodeInc != nullptr;
561 228 : psNodeInc = psNodeInc->psNext)
562 : {
563 228 : if (EQUAL(psNodeInc->pszValue, "referenceNoiseLevel"))
564 : {
565 : CPLXMLNode *psCalibType =
566 18 : CPLGetXMLNode(psNodeInc, "sarCalibrationType");
567 : CPLXMLNode *psPixelFirstNoiseValue =
568 18 : CPLGetXMLNode(psNodeInc, "pixelFirstNoiseValue");
569 18 : CPLXMLNode *psStepSize = CPLGetXMLNode(psNodeInc, "stepSize");
570 : CPLXMLNode *psNumberOfValues =
571 18 : CPLGetXMLNode(psNodeInc, "numberOfValues");
572 : CPLXMLNode *psNoiseLevelValues =
573 18 : CPLGetXMLNode(psNodeInc, "noiseLevelValues");
574 :
575 18 : if (psCalibType != nullptr && psPixelFirstNoiseValue != nullptr &&
576 18 : psStepSize != nullptr && psNumberOfValues != nullptr &&
577 : psNoiseLevelValues != nullptr)
578 : {
579 18 : const char *calibType = CPLGetXMLValue(psCalibType, "", "");
580 18 : this->pixelFirstLutValueNoiseLevels =
581 18 : atoi(CPLGetXMLValue(psPixelFirstNoiseValue, "", "0"));
582 18 : this->stepSizeNoiseLevels =
583 18 : atoi(CPLGetXMLValue(psStepSize, "", "0"));
584 18 : this->numberOfValuesNoiseLevels =
585 18 : atoi(CPLGetXMLValue(psNumberOfValues, "", "0"));
586 : const char *noiseLevelValues =
587 18 : CPLGetXMLValue(psNoiseLevelValues, "", "");
588 18 : if (this->stepSizeNoiseLevels > 0 &&
589 0 : this->numberOfValuesNoiseLevels != INT_MIN &&
590 0 : abs(this->numberOfValuesNoiseLevels) <
591 0 : INT_MAX / this->stepSizeNoiseLevels)
592 : {
593 0 : char **papszNoiseLevelList = CSLTokenizeString2(
594 : noiseLevelValues, " ", CSLT_HONOURSTRINGS);
595 : /* Get the Pixel Per range */
596 0 : this->m_nTableNoiseLevelsSize =
597 0 : abs(this->stepSizeNoiseLevels) *
598 0 : abs(this->numberOfValuesNoiseLevels);
599 :
600 0 : if ((EQUAL(calibType, "Beta Nought") &&
601 0 : this->m_eCalib == Beta0) ||
602 0 : (EQUAL(calibType, "Sigma Nought") &&
603 0 : this->m_eCalib == Sigma0) ||
604 0 : (EQUAL(calibType, "Gamma") && this->m_eCalib == Gamma))
605 : {
606 : /* Allocate the right Noise Levels size according to the
607 : * product range pixel */
608 0 : this->m_nfTableNoiseLevels = InterpolateValues(
609 : papszNoiseLevelList, this->m_nTableNoiseLevelsSize,
610 : this->stepSizeNoiseLevels,
611 : this->numberOfValuesNoiseLevels,
612 : this->pixelFirstLutValueNoiseLevels);
613 : }
614 :
615 0 : CSLDestroy(papszNoiseLevelList);
616 : }
617 :
618 18 : if (this->m_nfTableNoiseLevels != nullptr)
619 : {
620 0 : break; // We are done
621 : }
622 : }
623 : }
624 : }
625 : }
626 :
627 : /************************************************************************/
628 : /* RCMCalibRasterBand() */
629 : /************************************************************************/
630 :
631 6 : RCMCalibRasterBand::RCMCalibRasterBand(
632 : RCMDataset *poDataset, const char *pszPolarization, GDALDataType eType,
633 : GDALDataset *poBandDataset, eCalibration eCalib, const char *pszLUT,
634 6 : const char *pszNoiseLevels, GDALDataType eOriginalType)
635 : : m_eCalib(eCalib), m_poBandDataset(poBandDataset),
636 12 : m_eOriginalType(eOriginalType), m_pszLUTFile(VSIStrdup(pszLUT)),
637 6 : m_pszNoiseLevelsFile(VSIStrdup(pszNoiseLevels))
638 : {
639 6 : this->poDS = poDataset;
640 :
641 6 : if (pszPolarization != nullptr && strlen(pszPolarization) != 0)
642 : {
643 6 : SetMetadataItem("POLARIMETRIC_INTERP", pszPolarization);
644 : }
645 :
646 6 : if ((eType == GDT_CInt16) || (eType == GDT_CFloat32))
647 0 : this->eDataType = GDT_CFloat32;
648 : else
649 6 : this->eDataType = GDT_Float32;
650 :
651 6 : GDALRasterBand *poRasterBand = poBandDataset->GetRasterBand(1);
652 6 : poRasterBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
653 :
654 6 : ReadLUT();
655 6 : ReadNoiseLevels();
656 6 : }
657 :
658 : /************************************************************************/
659 : /* ~RCMCalibRasterBand() */
660 : /************************************************************************/
661 :
662 12 : RCMCalibRasterBand::~RCMCalibRasterBand()
663 : {
664 6 : CPLFree(m_nfTable);
665 6 : CPLFree(m_nfTableNoiseLevels);
666 6 : CPLFree(m_pszLUTFile);
667 6 : CPLFree(m_pszNoiseLevelsFile);
668 6 : GDALClose(m_poBandDataset);
669 12 : }
670 :
671 : /************************************************************************/
672 : /* IReadBlock() */
673 : /************************************************************************/
674 :
675 0 : CPLErr RCMCalibRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
676 : void *pImage)
677 : {
678 : CPLErr eErr;
679 0 : int nRequestXSize = 0;
680 0 : int nRequestYSize = 0;
681 0 : GetActualBlockSize(nBlockXOff, nBlockYOff, &nRequestXSize, &nRequestYSize);
682 :
683 : // Zero initial partial right-most and bottom-most blocks
684 0 : if (nRequestXSize < nBlockXSize || nRequestYSize < nBlockYSize)
685 : {
686 0 : memset(pImage, 0,
687 0 : static_cast<size_t>(GDALGetDataTypeSizeBytes(eDataType)) *
688 0 : nBlockXSize * nBlockYSize);
689 : }
690 :
691 0 : if (this->m_eOriginalType == GDT_CInt16)
692 : {
693 : /* read in complex values */
694 : GInt16 *panImageTmp = static_cast<GInt16 *>(
695 0 : VSI_MALLOC3_VERBOSE(nBlockXSize, nBlockYSize,
696 : GDALGetDataTypeSizeBytes(m_eOriginalType)));
697 0 : if (!panImageTmp)
698 0 : return CE_Failure;
699 :
700 0 : if (m_poBandDataset->GetRasterCount() == 2)
701 : {
702 0 : eErr = m_poBandDataset->RasterIO(
703 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
704 : nRequestXSize, nRequestYSize, panImageTmp, nRequestXSize,
705 : nRequestYSize, this->m_eOriginalType, 2, nullptr, 4,
706 0 : nBlockXSize * 4, 4, nullptr);
707 :
708 : /*
709 : eErr = m_poBandDataset->RasterIO(GF_Read,
710 : nBlockXOff * nBlockXSize,
711 : nBlockYOff * nBlockYSize,
712 : nRequestXSize, nRequestYSize,
713 : panImageTmp, nRequestXSize, nRequestYSize,
714 : GDT_Int32,
715 : 2, nullptr, 4, nBlockXSize * 4, 2, nullptr);
716 : */
717 : }
718 : else
719 : {
720 0 : eErr = m_poBandDataset->RasterIO(
721 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
722 : nRequestXSize, nRequestYSize, panImageTmp, nRequestXSize,
723 : nRequestYSize, this->m_eOriginalType, 1, nullptr, 4,
724 0 : nBlockXSize * 4, 0, nullptr);
725 :
726 : /*
727 : eErr =
728 : m_poBandDataset->RasterIO(GF_Read,
729 : nBlockXOff * nBlockXSize,
730 : nBlockYOff * nBlockYSize,
731 : nRequestXSize, nRequestYSize,
732 : panImageTmp, nRequestXSize, nRequestYSize,
733 : GDT_UInt32,
734 : 1, nullptr, 4, nBlockXSize * 4, 0, nullptr);
735 : */
736 :
737 : #ifdef CPL_LSB
738 : /* First, undo the 32bit swap. */
739 0 : GDALSwapWords(pImage, 4, nBlockXSize * nBlockYSize, 4);
740 :
741 : /* Then apply 16 bit swap. */
742 0 : GDALSwapWords(pImage, 2, nBlockXSize * nBlockYSize * 2, 2);
743 : #endif
744 : }
745 :
746 : /* calibrate the complex values */
747 0 : for (int i = 0; i < nRequestYSize; i++)
748 : {
749 0 : for (int j = 0; j < nRequestXSize; j++)
750 : {
751 : /* calculate pixel offset in memory*/
752 0 : const int nPixOff = 2 * (i * nBlockXSize + j);
753 0 : const int nTruePixOff = (i * nBlockXSize) + j;
754 :
755 : // Formula for Complex Q+J
756 0 : const float real = static_cast<float>(panImageTmp[nPixOff]);
757 0 : const float img = static_cast<float>(panImageTmp[nPixOff + 1]);
758 0 : const float digitalValue = (real * real) + (img * img);
759 0 : const float lutValue =
760 0 : static_cast<float>(m_nfTable[nBlockXOff * nBlockXSize + j]);
761 0 : const float calibValue = digitalValue / (lutValue * lutValue);
762 :
763 0 : reinterpret_cast<float *>(pImage)[nTruePixOff] = calibValue;
764 : }
765 : }
766 :
767 0 : CPLFree(panImageTmp);
768 : }
769 :
770 : // If the underlying file is NITF CFloat32
771 0 : else if (this->m_eOriginalType == GDT_CFloat32 ||
772 0 : this->m_eOriginalType == GDT_CFloat64)
773 : {
774 : /* read in complex values */
775 : const int dataTypeSize =
776 0 : GDALGetDataTypeSizeBytes(this->m_eOriginalType);
777 0 : const GDALDataType bandFileType = this->m_eOriginalType;
778 0 : const int bandFileDataTypeSize = GDALGetDataTypeSizeBytes(bandFileType);
779 :
780 : /* read the original image complex values in a temporary image space */
781 0 : float *pafImageTmp = static_cast<float *>(VSI_MALLOC3_VERBOSE(
782 : nBlockXSize, nBlockYSize, 2 * bandFileDataTypeSize));
783 0 : if (!pafImageTmp)
784 0 : return CE_Failure;
785 :
786 : eErr =
787 : // I and Q from each band are pixel-interleaved into this complex
788 : // band
789 0 : m_poBandDataset->RasterIO(
790 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
791 : nRequestXSize, nRequestYSize, pafImageTmp, nRequestXSize,
792 : nRequestYSize, bandFileType, 2, nullptr, dataTypeSize,
793 0 : static_cast<GSpacing>(dataTypeSize) * nBlockXSize,
794 : bandFileDataTypeSize, nullptr);
795 :
796 : /* calibrate the complex values */
797 0 : for (int i = 0; i < nRequestYSize; i++)
798 : {
799 0 : for (int j = 0; j < nRequestXSize; j++)
800 : {
801 : /* calculate pixel offset in memory*/
802 0 : const int nPixOff = 2 * (i * nBlockXSize + j);
803 0 : const int nTruePixOff = (i * nBlockXSize) + j;
804 :
805 : // Formula for Complex Q+J
806 0 : const float real = pafImageTmp[nPixOff];
807 0 : const float img = pafImageTmp[nPixOff + 1];
808 0 : const float digitalValue = (real * real) + (img * img);
809 0 : const float lutValue =
810 0 : static_cast<float>(m_nfTable[nBlockXOff * nBlockXSize + j]);
811 0 : const float calibValue = digitalValue / (lutValue * lutValue);
812 :
813 0 : reinterpret_cast<float *>(pImage)[nTruePixOff] = calibValue;
814 : }
815 : }
816 :
817 0 : CPLFree(pafImageTmp);
818 : }
819 :
820 0 : else if (this->m_eOriginalType == GDT_Float32)
821 : {
822 : /* A Float32 is actual 4 bytes */
823 0 : eErr = m_poBandDataset->RasterIO(
824 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
825 : nRequestXSize, nRequestYSize, pImage, nRequestXSize, nRequestYSize,
826 0 : GDT_Float32, 1, nullptr, 4, nBlockXSize * 4, 0, nullptr);
827 :
828 : /* iterate over detected values */
829 0 : for (int i = 0; i < nRequestYSize; i++)
830 : {
831 0 : for (int j = 0; j < nRequestXSize; j++)
832 : {
833 0 : int nPixOff = (i * nBlockXSize) + j;
834 :
835 : /* For detected products, in order to convert the digital number
836 : of a given range sample to a calibrated value, the digital value
837 : is first squared, then the offset(B) is added and the result is
838 : divided by the gains value(A) corresponding to the range sample.
839 : RCM-SP-53-0419 Issue 2/5: January 2, 2018 Page 7-56 */
840 0 : const float digitalValue =
841 0 : reinterpret_cast<float *>(pImage)[nPixOff];
842 0 : const float A =
843 0 : static_cast<float>(m_nfTable[nBlockXOff * nBlockXSize + j]);
844 0 : reinterpret_cast<float *>(pImage)[nPixOff] =
845 0 : ((digitalValue * digitalValue) +
846 0 : static_cast<float>(this->m_nfOffset)) /
847 : A;
848 : }
849 : }
850 : }
851 :
852 0 : else if (this->m_eOriginalType == GDT_Float64)
853 : {
854 : /* A Float64 is actual 8 bytes */
855 0 : eErr = m_poBandDataset->RasterIO(
856 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
857 : nRequestXSize, nRequestYSize, pImage, nRequestXSize, nRequestYSize,
858 0 : GDT_Float64, 1, nullptr, 8, nBlockXSize * 8, 0, nullptr);
859 :
860 : /* iterate over detected values */
861 0 : for (int i = 0; i < nRequestYSize; i++)
862 : {
863 0 : for (int j = 0; j < nRequestXSize; j++)
864 : {
865 0 : int nPixOff = (i * nBlockXSize) + j;
866 :
867 : /* For detected products, in order to convert the digital number
868 : of a given range sample to a calibrated value, the digital value
869 : is first squared, then the offset(B) is added and the result is
870 : divided by the gains value(A) corresponding to the range sample.
871 : RCM-SP-53-0419 Issue 2/5: January 2, 2018 Page 7-56 */
872 0 : const float digitalValue =
873 0 : reinterpret_cast<float *>(pImage)[nPixOff];
874 0 : const float A =
875 0 : static_cast<float>(m_nfTable[nBlockXOff * nBlockXSize + j]);
876 0 : reinterpret_cast<float *>(pImage)[nPixOff] =
877 0 : ((digitalValue * digitalValue) +
878 0 : static_cast<float>(this->m_nfOffset)) /
879 : A;
880 : }
881 : }
882 : }
883 :
884 0 : else if (this->m_eOriginalType == GDT_UInt16)
885 : {
886 : /* read in detected values */
887 0 : GUInt16 *panImageTmp = static_cast<GUInt16 *>(VSI_MALLOC3_VERBOSE(
888 : nBlockXSize, nBlockYSize, GDALGetDataTypeSizeBytes(GDT_UInt16)));
889 0 : if (!panImageTmp)
890 0 : return CE_Failure;
891 0 : eErr = m_poBandDataset->RasterIO(
892 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
893 : nRequestXSize, nRequestYSize, panImageTmp, nRequestXSize,
894 0 : nRequestYSize, GDT_UInt16, 1, nullptr, 2, nBlockXSize * 2, 0,
895 : nullptr);
896 :
897 : /* iterate over detected values */
898 0 : for (int i = 0; i < nRequestYSize; i++)
899 : {
900 0 : for (int j = 0; j < nRequestXSize; j++)
901 : {
902 0 : const int nPixOff = (i * nBlockXSize) + j;
903 :
904 0 : const float digitalValue =
905 0 : static_cast<float>(panImageTmp[nPixOff]);
906 0 : const float A =
907 0 : static_cast<float>(m_nfTable[nBlockXOff * nBlockXSize + j]);
908 0 : reinterpret_cast<float *>(pImage)[nPixOff] =
909 0 : ((digitalValue * digitalValue) +
910 0 : static_cast<float>(this->m_nfOffset)) /
911 : A;
912 : }
913 : }
914 0 : CPLFree(panImageTmp);
915 : } /* Ticket #2104: Support for ScanSAR products */
916 :
917 0 : else if (this->m_eOriginalType == GDT_Byte)
918 : {
919 : GByte *pabyImageTmp =
920 0 : static_cast<GByte *>(VSI_MALLOC2_VERBOSE(nBlockXSize, nBlockYSize));
921 0 : if (!pabyImageTmp)
922 0 : return CE_Failure;
923 0 : eErr = m_poBandDataset->RasterIO(
924 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
925 : nRequestXSize, nRequestYSize, pabyImageTmp, nRequestXSize,
926 0 : nRequestYSize, GDT_Byte, 1, nullptr, 1, nBlockXSize, 0, nullptr);
927 :
928 : /* iterate over detected values */
929 0 : for (int i = 0; i < nRequestYSize; i++)
930 : {
931 0 : for (int j = 0; j < nRequestXSize; j++)
932 : {
933 0 : const int nPixOff = (i * nBlockXSize) + j;
934 :
935 0 : const float digitalValue =
936 0 : static_cast<float>(pabyImageTmp[nPixOff]);
937 0 : const float A =
938 0 : static_cast<float>(m_nfTable[nBlockXOff * nBlockXSize + j]);
939 0 : reinterpret_cast<float *>(pImage)[nPixOff] =
940 0 : ((digitalValue * digitalValue) +
941 0 : static_cast<float>(this->m_nfOffset)) /
942 : A;
943 : }
944 : }
945 0 : CPLFree(pabyImageTmp);
946 : }
947 : else
948 : {
949 0 : CPLAssert(FALSE);
950 : return CE_Failure;
951 : }
952 0 : return eErr;
953 : }
954 :
955 : /************************************************************************/
956 : /* ==================================================================== */
957 : /* RCMDataset */
958 : /* ==================================================================== */
959 : /************************************************************************/
960 :
961 : /************************************************************************/
962 : /* RCMDataset() */
963 : /************************************************************************/
964 :
965 8 : RCMDataset::RCMDataset()
966 : {
967 8 : adfGeoTransform[0] = 0.0;
968 8 : adfGeoTransform[1] = 1.0;
969 8 : adfGeoTransform[2] = 0.0;
970 8 : adfGeoTransform[3] = 0.0;
971 8 : adfGeoTransform[4] = 0.0;
972 8 : adfGeoTransform[5] = 1.0;
973 8 : }
974 :
975 : /************************************************************************/
976 : /* ~RCMDataset() */
977 : /************************************************************************/
978 :
979 16 : RCMDataset::~RCMDataset()
980 :
981 : {
982 8 : RCMDataset::FlushCache(true);
983 :
984 8 : if (nGCPCount > 0)
985 : {
986 8 : GDALDeinitGCPs(nGCPCount, pasGCPList);
987 8 : CPLFree(pasGCPList);
988 : }
989 :
990 8 : RCMDataset::CloseDependentDatasets();
991 :
992 8 : if (papszSubDatasets != nullptr)
993 4 : CSLDestroy(papszSubDatasets);
994 :
995 8 : if (papszExtraFiles != nullptr)
996 8 : CSLDestroy(papszExtraFiles);
997 :
998 8 : if (m_nfIncidenceAngleTable != nullptr)
999 8 : CPLFree(m_nfIncidenceAngleTable);
1000 16 : }
1001 :
1002 : /************************************************************************/
1003 : /* CloseDependentDatasets() */
1004 : /************************************************************************/
1005 :
1006 8 : int RCMDataset::CloseDependentDatasets()
1007 : {
1008 8 : int bHasDroppedRef = GDALPamDataset::CloseDependentDatasets();
1009 :
1010 8 : if (nBands != 0)
1011 8 : bHasDroppedRef = TRUE;
1012 :
1013 24 : for (int iBand = 0; iBand < nBands; iBand++)
1014 : {
1015 16 : delete papoBands[iBand];
1016 : }
1017 8 : nBands = 0;
1018 :
1019 8 : return bHasDroppedRef;
1020 : }
1021 :
1022 : /************************************************************************/
1023 : /* GetFileList() */
1024 : /************************************************************************/
1025 :
1026 0 : char **RCMDataset::GetFileList()
1027 :
1028 : {
1029 0 : char **papszFileList = GDALPamDataset::GetFileList();
1030 :
1031 0 : papszFileList = CSLInsertStrings(papszFileList, -1, papszExtraFiles);
1032 :
1033 0 : return papszFileList;
1034 : }
1035 :
1036 : /************************************************************************/
1037 : /* Open() */
1038 : /************************************************************************/
1039 :
1040 10 : GDALDataset *RCMDataset::Open(GDALOpenInfo *poOpenInfo)
1041 :
1042 : {
1043 : /* -------------------------------------------------------------------- */
1044 : /* Is this a RCM Product.xml definition? */
1045 : /* -------------------------------------------------------------------- */
1046 10 : if (!RCMDatasetIdentify(poOpenInfo))
1047 : {
1048 0 : return nullptr;
1049 : }
1050 :
1051 : /* -------------------------------------------------------------------- */
1052 : /* Get subdataset information, if relevant */
1053 : /* -------------------------------------------------------------------- */
1054 10 : const char *pszFilename = poOpenInfo->pszFilename;
1055 10 : eCalibration eCalib = None;
1056 :
1057 10 : if (STARTS_WITH_CI(pszFilename, szLayerCalibration) &&
1058 6 : pszFilename[strlen(szLayerCalibration)] == chLayerSeparator)
1059 : {
1060 : // The calibration name and filename begins after the hard coded layer
1061 : // name
1062 6 : pszFilename += strlen(szLayerCalibration) + 1;
1063 :
1064 6 : if (STARTS_WITH_CI(pszFilename, szBETA0))
1065 : {
1066 1 : eCalib = Beta0;
1067 : }
1068 5 : else if (STARTS_WITH_CI(pszFilename, szSIGMA0))
1069 : {
1070 1 : eCalib = Sigma0;
1071 : }
1072 4 : else if (STARTS_WITH_CI(pszFilename, szGAMMA) ||
1073 3 : STARTS_WITH_CI(pszFilename, "GAMMA0")) // Cover both situation
1074 : {
1075 1 : eCalib = Gamma;
1076 : }
1077 3 : else if (STARTS_WITH_CI(pszFilename, szUNCALIB))
1078 : {
1079 2 : eCalib = Uncalib;
1080 : }
1081 : else
1082 : {
1083 1 : CPLError(CE_Failure, CPLE_NotSupported,
1084 : "Unsupported calibration type");
1085 1 : return nullptr;
1086 : }
1087 :
1088 : /* advance the pointer to the actual filename */
1089 35 : while (*pszFilename != '\0' && *pszFilename != ':')
1090 30 : pszFilename++;
1091 :
1092 5 : if (*pszFilename == ':')
1093 5 : pszFilename++;
1094 :
1095 : // need to redo the directory check:
1096 : // the GDALOpenInfo check would have failed because of the calibration
1097 : // string on the filename
1098 : VSIStatBufL sStat;
1099 5 : if (VSIStatL(pszFilename, &sStat) == 0)
1100 4 : poOpenInfo->bIsDirectory = VSI_ISDIR(sStat.st_mode);
1101 : }
1102 :
1103 18 : CPLString osMDFilename;
1104 9 : if (poOpenInfo->bIsDirectory)
1105 : {
1106 : /* Check for directory access when there is a product.xml file in the
1107 : directory. */
1108 : osMDFilename =
1109 2 : CPLFormCIFilenameSafe(pszFilename, "product.xml", nullptr);
1110 :
1111 : VSIStatBufL sStat;
1112 2 : if (VSIStatL(osMDFilename, &sStat) != 0)
1113 : {
1114 : /* If not, check for directory extra 'metadata' access when there is
1115 : a product.xml file in the directory. */
1116 1 : osMDFilename = CPLFormCIFilenameSafe(pszFilename,
1117 2 : GetMetadataProduct(), nullptr);
1118 : }
1119 : }
1120 : else
1121 7 : osMDFilename = pszFilename;
1122 :
1123 : /* -------------------------------------------------------------------- */
1124 : /* Ingest the Product.xml file. */
1125 : /* -------------------------------------------------------------------- */
1126 18 : CPLXMLTreeCloser psProduct(CPLParseXMLFile(osMDFilename));
1127 9 : if (!psProduct)
1128 1 : return nullptr;
1129 :
1130 : /* -------------------------------------------------------------------- */
1131 : /* Confirm the requested access is supported. */
1132 : /* -------------------------------------------------------------------- */
1133 8 : if (poOpenInfo->eAccess == GA_Update)
1134 : {
1135 0 : CPLError(CE_Failure, CPLE_NotSupported,
1136 : "ERROR: The RCM driver does not support update "
1137 : "access to existing dataset.");
1138 0 : return nullptr;
1139 : }
1140 :
1141 : CPLXMLNode *psSceneAttributes =
1142 8 : CPLGetXMLNode(psProduct.get(), "=product.sceneAttributes");
1143 8 : if (psSceneAttributes == nullptr)
1144 : {
1145 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1146 : "ERROR: Failed to find <sceneAttributes> in document.");
1147 0 : return nullptr;
1148 : }
1149 :
1150 8 : CPLXMLNode *psImageAttributes = CPLGetXMLNode(
1151 : psProduct.get(), "=product.sceneAttributes.imageAttributes");
1152 8 : if (psImageAttributes == nullptr)
1153 : {
1154 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1155 : "ERROR: Failed to find <sceneAttributes.imageAttributes> in "
1156 : "document.");
1157 0 : return nullptr;
1158 : }
1159 :
1160 : int numberOfEntries =
1161 8 : atoi(CPLGetXMLValue(psSceneAttributes, "numberOfEntries", "0"));
1162 8 : if (numberOfEntries != 1)
1163 : {
1164 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1165 : "ERROR: Only RCM with Complex Single-beam is supported.");
1166 0 : return nullptr;
1167 : }
1168 :
1169 : CPLXMLNode *psImageReferenceAttributes =
1170 8 : CPLGetXMLNode(psProduct.get(), "=product.imageReferenceAttributes");
1171 8 : if (psImageReferenceAttributes == nullptr)
1172 : {
1173 0 : CPLError(
1174 : CE_Failure, CPLE_OpenFailed,
1175 : "ERROR: Failed to find <imageReferenceAttributes> in document.");
1176 0 : return nullptr;
1177 : }
1178 :
1179 : CPLXMLNode *psImageGenerationParameters =
1180 8 : CPLGetXMLNode(psProduct.get(), "=product.imageGenerationParameters");
1181 8 : if (psImageGenerationParameters == nullptr)
1182 : {
1183 0 : CPLError(
1184 : CE_Failure, CPLE_OpenFailed,
1185 : "ERROR: Failed to find <imageGenerationParameters> in document.");
1186 0 : return nullptr;
1187 : }
1188 :
1189 : /* -------------------------------------------------------------------- */
1190 : /* Create the dataset. */
1191 : /* -------------------------------------------------------------------- */
1192 16 : auto poDS = std::make_unique<RCMDataset>();
1193 :
1194 8 : poDS->psProduct = std::move(psProduct);
1195 :
1196 : /* -------------------------------------------------------------------- */
1197 : /* Get overall image information. */
1198 : /* -------------------------------------------------------------------- */
1199 8 : poDS->nRasterXSize = atoi(CPLGetXMLValue(
1200 : psSceneAttributes, "imageAttributes.samplesPerLine", "-1"));
1201 8 : poDS->nRasterYSize = atoi(
1202 : CPLGetXMLValue(psSceneAttributes, "imageAttributes.numLines", "-1"));
1203 8 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
1204 : {
1205 0 : return nullptr;
1206 : }
1207 :
1208 : /* -------------------------------------------------------------------- */
1209 : /* Check product type, as to determine if there are LUTs for */
1210 : /* calibration purposes. */
1211 : /* -------------------------------------------------------------------- */
1212 :
1213 : const char *pszItem =
1214 8 : CPLGetXMLValue(psImageGenerationParameters,
1215 : "generalProcessingInformation.productType", "UNK");
1216 8 : poDS->SetMetadataItem("PRODUCT_TYPE", pszItem);
1217 8 : const char *pszProductType = pszItem;
1218 :
1219 : pszItem =
1220 8 : CPLGetXMLValue(poDS->psProduct.get(), "=product.productId", "UNK");
1221 8 : poDS->SetMetadataItem("PRODUCT_ID", pszItem);
1222 :
1223 8 : pszItem = CPLGetXMLValue(
1224 8 : poDS->psProduct.get(),
1225 : "=product.securityAttributes.securityClassification", "UNK");
1226 8 : poDS->SetMetadataItem("SECURITY_CLASSIFICATION", pszItem);
1227 :
1228 : pszItem =
1229 8 : CPLGetXMLValue(poDS->psProduct.get(),
1230 : "=product.sourceAttributes.polarizationDataMode", "UNK");
1231 8 : poDS->SetMetadataItem("POLARIZATION_DATA_MODE", pszItem);
1232 :
1233 8 : pszItem = CPLGetXMLValue(psImageGenerationParameters,
1234 : "generalProcessingInformation.processingFacility",
1235 : "UNK");
1236 8 : poDS->SetMetadataItem("PROCESSING_FACILITY", pszItem);
1237 :
1238 : pszItem =
1239 8 : CPLGetXMLValue(psImageGenerationParameters,
1240 : "generalProcessingInformation.processingTime", "UNK");
1241 8 : poDS->SetMetadataItem("PROCESSING_TIME", pszItem);
1242 :
1243 8 : pszItem = CPLGetXMLValue(psImageGenerationParameters,
1244 : "sarProcessingInformation.satelliteHeight", "UNK");
1245 8 : poDS->SetMetadataItem("SATELLITE_HEIGHT", pszItem);
1246 :
1247 8 : pszItem = CPLGetXMLValue(
1248 : psImageGenerationParameters,
1249 : "sarProcessingInformation.zeroDopplerTimeFirstLine", "UNK");
1250 8 : poDS->SetMetadataItem("FIRST_LINE_TIME", pszItem);
1251 :
1252 8 : pszItem = CPLGetXMLValue(psImageGenerationParameters,
1253 : "sarProcessingInformation.zeroDopplerTimeLastLine",
1254 : "UNK");
1255 8 : poDS->SetMetadataItem("LAST_LINE_TIME", pszItem);
1256 :
1257 8 : pszItem = CPLGetXMLValue(psImageGenerationParameters,
1258 : "sarProcessingInformation.lutApplied", "");
1259 8 : poDS->SetMetadataItem("LUT_APPLIED", pszItem);
1260 :
1261 : /*---------------------------------------------------------------------
1262 : If true, a polarization dependent application LUT has been applied
1263 : for each polarization channel. Otherwise the same application LUT
1264 : has been applied for all polarization channels. Not applicable to
1265 : lookupTable = "Unity*" or if dataType = "Floating-Point".
1266 : --------------------------------------------------------------------- */
1267 8 : pszItem = CPLGetXMLValue(psImageGenerationParameters,
1268 : "sarProcessingInformation.perPolarizationScaling",
1269 : "false");
1270 8 : poDS->SetMetadataItem("PER_POLARIZATION_SCALING", pszItem);
1271 8 : if (EQUAL(pszItem, "true") || EQUAL(pszItem, "TRUE"))
1272 : {
1273 8 : poDS->bPerPolarizationScaling = TRUE;
1274 : }
1275 :
1276 : /* the following cases can be assumed to have no LUTs, as per
1277 : * RN-RP-51-2713, but also common sense
1278 : * SLC represents a SLant range georeferenced Complex product
1279 : * (i.e., equivalent to a Single-Look Complex product for RADARSAT-1 or
1280 : * RADARSAT-2). • GRD or GRC represent GRound range georeferenced Detected
1281 : * or Complex products (GRD is equivalent to an SGX, SCN or SCW product for
1282 : * RADARSAT1 or RADARSAT-2). • GCD or GCC represent GeoCoded Detected or
1283 : * Complex products (GCD is equivalent to an SSG or SPG product for
1284 : * RADARSAT-1 or RADARSAT-2).
1285 : */
1286 8 : bool bCanCalib = false;
1287 8 : if (!(STARTS_WITH_CI(pszProductType, "UNK") ||
1288 8 : STARTS_WITH_CI(pszProductType, "GCD") ||
1289 8 : STARTS_WITH_CI(pszProductType, "GCC")))
1290 : {
1291 8 : bCanCalib = true;
1292 : }
1293 :
1294 : /* -------------------------------------------------------------------- */
1295 : /* Get dataType (so we can recognise complex data), and the */
1296 : /* bitsPerSample. */
1297 : /* -------------------------------------------------------------------- */
1298 8 : const char *pszSampleDataType = CPLGetXMLValue(
1299 : psImageReferenceAttributes, "rasterAttributes.sampleType", "");
1300 8 : poDS->SetMetadataItem("SAMPLE_TYPE", pszSampleDataType);
1301 :
1302 : /* Either Integer (16 bits) or Floating-Point (32 bits) */
1303 8 : const char *pszDataType = CPLGetXMLValue(psImageReferenceAttributes,
1304 : "rasterAttributes.dataType", "");
1305 8 : poDS->SetMetadataItem("DATA_TYPE", pszDataType);
1306 :
1307 : /* 2 entries for complex data 1 entry for magnitude detected data */
1308 8 : const char *pszBitsPerSample = CPLGetXMLValue(
1309 : psImageReferenceAttributes, "rasterAttributes.bitsPerSample", "");
1310 8 : const int nBitsPerSample = atoi(pszBitsPerSample);
1311 8 : poDS->SetMetadataItem("BITS_PER_SAMPLE", pszBitsPerSample);
1312 :
1313 : const char *pszSampledPixelSpacingTime =
1314 8 : CPLGetXMLValue(psImageReferenceAttributes,
1315 : "rasterAttributes.sampledPixelSpacingTime", "UNK");
1316 8 : poDS->SetMetadataItem("SAMPLED_PIXEL_SPACING_TIME",
1317 8 : pszSampledPixelSpacingTime);
1318 :
1319 : const char *pszSampledLineSpacingTime =
1320 8 : CPLGetXMLValue(psImageReferenceAttributes,
1321 : "rasterAttributes.sampledLineSpacingTime", "UNK");
1322 8 : poDS->SetMetadataItem("SAMPLED_LINE_SPACING_TIME",
1323 8 : pszSampledLineSpacingTime);
1324 :
1325 : GDALDataType eDataType;
1326 :
1327 8 : if (EQUAL(pszSampleDataType, "Mixed")) /* RCM MLC has Mixed sampleType */
1328 : {
1329 1 : poDS->isComplexData = false; /* RCM MLC is detected, non-complex */
1330 1 : if (nBitsPerSample == 32)
1331 : {
1332 1 : eDataType = GDT_Float32; /* 32 bits, check read block */
1333 1 : poDS->magnitudeBits = 32;
1334 : }
1335 : else
1336 : {
1337 0 : eDataType = GDT_UInt16; /* 16 bits, check read block */
1338 0 : poDS->magnitudeBits = 16;
1339 : }
1340 : }
1341 7 : else if (EQUAL(pszSampleDataType, "Complex"))
1342 : {
1343 0 : poDS->isComplexData = true;
1344 : /* Usually this is the same bits for both */
1345 0 : poDS->realBitsComplexData = nBitsPerSample;
1346 0 : poDS->imaginaryBitsComplexData = nBitsPerSample;
1347 :
1348 0 : if (nBitsPerSample == 32)
1349 : {
1350 0 : eDataType = GDT_CFloat32; /* 32 bits, check read block */
1351 : }
1352 : else
1353 : {
1354 0 : eDataType = GDT_CInt16; /* 16 bits, check read block */
1355 : }
1356 : }
1357 7 : else if (nBitsPerSample == 32 &&
1358 0 : EQUAL(pszSampleDataType, "Magnitude Detected"))
1359 : {
1360 : /* Actually, I don't need to test the 'dataType'=' Floating-Point', we
1361 : * know that a 32 bits per sample */
1362 0 : eDataType = GDT_Float32;
1363 0 : poDS->isComplexData = false;
1364 0 : poDS->magnitudeBits = 32;
1365 : }
1366 7 : else if (nBitsPerSample == 16 &&
1367 7 : EQUAL(pszSampleDataType, "Magnitude Detected"))
1368 : {
1369 : /* Actually, I don't need to test the 'dataType'=' Integer', we know
1370 : * that a 16 bits per sample */
1371 7 : eDataType = GDT_UInt16;
1372 7 : poDS->isComplexData = false;
1373 7 : poDS->magnitudeBits = 16;
1374 : }
1375 : else
1376 : {
1377 0 : CPLError(CE_Failure, CPLE_AppDefined,
1378 : "ERROR: dataType=%s and bitsPerSample=%d are not a supported "
1379 : "configuration.",
1380 : pszDataType, nBitsPerSample);
1381 0 : return nullptr;
1382 : }
1383 :
1384 : /* Indicates whether pixel number (i.e., range) increases or decreases with
1385 : range time. For GCD and GCC products, this applies to intermediate ground
1386 : range image data prior to geocoding. */
1387 : const char *pszPixelTimeOrdering =
1388 8 : CPLGetXMLValue(psImageReferenceAttributes,
1389 : "rasterAttributes.pixelTimeOrdering", "UNK");
1390 8 : poDS->SetMetadataItem("PIXEL_TIME_ORDERING", pszPixelTimeOrdering);
1391 :
1392 : /* Indicates whether line numbers (i.e., azimuth) increase or decrease with
1393 : azimuth time. For GCD and GCC products, this applies to intermediate ground
1394 : range image data prior to geocoding. */
1395 8 : const char *pszLineTimeOrdering = CPLGetXMLValue(
1396 : psImageReferenceAttributes, "rasterAttributes.lineTimeOrdering", "UNK");
1397 8 : poDS->SetMetadataItem("LINE_TIME_ORDERING", pszLineTimeOrdering);
1398 :
1399 : /* while we're at it, extract the pixel spacing information */
1400 : const char *pszPixelSpacing =
1401 8 : CPLGetXMLValue(psImageReferenceAttributes,
1402 : "rasterAttributes.sampledPixelSpacing", "UNK");
1403 8 : poDS->SetMetadataItem("PIXEL_SPACING", pszPixelSpacing);
1404 :
1405 : const char *pszLineSpacing =
1406 8 : CPLGetXMLValue(psImageReferenceAttributes,
1407 : "rasterAttributes.sampledLineSpacing", "UNK");
1408 8 : poDS->SetMetadataItem("LINE_SPACING", pszLineSpacing);
1409 :
1410 : /* -------------------------------------------------------------------- */
1411 : /* Open each of the data files as a complex band. */
1412 : /* -------------------------------------------------------------------- */
1413 8 : char *pszBeta0LUT = nullptr;
1414 8 : char *pszGammaLUT = nullptr;
1415 8 : char *pszSigma0LUT = nullptr;
1416 : // Same file for all calibrations except the calibration is plit inside the
1417 : // XML
1418 16 : std::string osNoiseLevelsValues;
1419 :
1420 16 : const CPLString osPath = CPLGetPathSafe(osMDFilename);
1421 :
1422 : /* Get a list of all polarizations */
1423 : CPLXMLNode *psSourceAttrs =
1424 8 : CPLGetXMLNode(poDS->psProduct.get(), "=product.sourceAttributes");
1425 8 : if (psSourceAttrs == nullptr)
1426 : {
1427 0 : CPLError(
1428 : CE_Failure, CPLE_OpenFailed,
1429 : "ERROR: RCM source attributes is missing. Please contact your data "
1430 : "provider for a corrected dataset.");
1431 0 : return nullptr;
1432 : }
1433 :
1434 8 : CPLXMLNode *psRadarParameters = CPLGetXMLNode(
1435 8 : poDS->psProduct.get(), "=product.sourceAttributes.radarParameters");
1436 8 : if (psRadarParameters == nullptr)
1437 : {
1438 0 : CPLError(
1439 : CE_Failure, CPLE_OpenFailed,
1440 : "ERROR: RCM radar parameters is missing. Please contact your data "
1441 : "provider for a corrected dataset.");
1442 0 : return nullptr;
1443 : }
1444 :
1445 : const char *pszPolarizations =
1446 8 : CPLGetXMLValue(psRadarParameters, "polarizations", "");
1447 8 : if (pszPolarizations == nullptr || strlen(pszPolarizations) == 0)
1448 : {
1449 0 : CPLError(
1450 : CE_Failure, CPLE_OpenFailed,
1451 : "ERROR: RCM polarizations list is missing. Please contact your "
1452 : "data provider for a corrected dataset.");
1453 0 : return nullptr;
1454 : }
1455 8 : poDS->SetMetadataItem("POLARIZATIONS", pszPolarizations);
1456 :
1457 : const char *psAcquisitionType =
1458 8 : CPLGetXMLValue(psRadarParameters, "acquisitionType", "UNK");
1459 8 : poDS->SetMetadataItem("ACQUISITION_TYPE", psAcquisitionType);
1460 :
1461 8 : const char *psBeams = CPLGetXMLValue(psRadarParameters, "beams", "UNK");
1462 8 : poDS->SetMetadataItem("BEAMS", psBeams);
1463 :
1464 : const CPLStringList aosPolarizationsGrids(
1465 16 : CSLTokenizeString2(pszPolarizations, " ", 0));
1466 16 : CPLStringList imageBandList;
1467 16 : CPLStringList imageBandFileList;
1468 8 : const int nPolarizationsGridCount = aosPolarizationsGrids.size();
1469 :
1470 : /* File names for full resolution IPDFs. For GeoTIFF format, one entry per
1471 : pole; For NITF 2.1 format, only one entry. */
1472 8 : bool bIsNITF = false;
1473 8 : const char *pszNITF_Filename = nullptr;
1474 8 : int imageBandFileCount = 0;
1475 8 : int imageBandCount = 0;
1476 :
1477 : /* Split the polarization string*/
1478 24 : auto iss = std::istringstream((CPLString(pszPolarizations)).c_str());
1479 16 : auto pol = std::string{};
1480 : /* Count number of polarizations*/
1481 24 : while (iss >> pol)
1482 16 : imageBandCount++;
1483 :
1484 8 : CPLXMLNode *psNode = psImageAttributes->psChild;
1485 131 : for (; psNode != nullptr; psNode = psNode->psNext)
1486 : {
1487 : /* Find the tif or ntf filename */
1488 123 : if (psNode->eType != CXT_Element || !(EQUAL(psNode->pszValue, "ipdf")))
1489 107 : continue;
1490 :
1491 : /* --------------------------------------------------------------------
1492 : */
1493 : /* Fetch ipdf image. Could be either tif or ntf. */
1494 : /* Replace / by \\ */
1495 : /* --------------------------------------------------------------------
1496 : */
1497 17 : const char *pszBasedFilename = CPLGetXMLValue(psNode, "", "");
1498 17 : if (*pszBasedFilename == '\0')
1499 0 : continue;
1500 :
1501 : /* Count number of image names within ipdf tag*/
1502 17 : imageBandFileCount++;
1503 :
1504 17 : CPLString pszUpperBasedFilename(CPLString(pszBasedFilename).toupper());
1505 :
1506 : const bool bEndsWithNTF =
1507 34 : strlen(pszUpperBasedFilename.c_str()) > 4 &&
1508 17 : EQUAL(pszUpperBasedFilename.c_str() +
1509 : strlen(pszUpperBasedFilename.c_str()) - 4,
1510 : ".NTF");
1511 :
1512 17 : if (bEndsWithNTF)
1513 : {
1514 : /* Found it! It would not exist one more */
1515 0 : bIsNITF = true;
1516 0 : pszNITF_Filename = pszBasedFilename;
1517 0 : break;
1518 : }
1519 : else
1520 : {
1521 : /* Keep adding polarizations filename according to the pole */
1522 17 : const char *pszPole = CPLGetXMLValue(psNode, "pole", "");
1523 17 : if (*pszPole == '\0')
1524 : {
1525 : /* Guard against case when pole is a null string, skip it */
1526 0 : imageBandFileCount--;
1527 0 : continue;
1528 : }
1529 :
1530 17 : if (EQUAL(pszPole, "XC"))
1531 : {
1532 : /* Skip RCM MLC's 3rd band file ##XC.tif */
1533 1 : imageBandFileCount--;
1534 1 : continue;
1535 : }
1536 :
1537 16 : imageBandList.AddString(CPLString(pszPole).toupper());
1538 16 : imageBandFileList.AddString(pszBasedFilename);
1539 : }
1540 : }
1541 :
1542 : /* -------------------------------------------------------------------- */
1543 : /* Incidence Angle in a sub-folder */
1544 : /* called 'calibration' from the 'metadata' folder */
1545 : /* -------------------------------------------------------------------- */
1546 :
1547 8 : const char *pszIncidenceAngleFileName = CPLGetXMLValue(
1548 : psImageReferenceAttributes, "incidenceAngleFileName", "");
1549 :
1550 8 : if (pszIncidenceAngleFileName != nullptr)
1551 : {
1552 16 : CPLString osIncidenceAngleFilePath = CPLFormFilenameSafe(
1553 8 : CPLFormFilenameSafe(osPath, CALIBRATION_FOLDER, nullptr).c_str(),
1554 16 : pszIncidenceAngleFileName, nullptr);
1555 :
1556 : /* Check if the file exist */
1557 8 : if (IsValidXMLFile(osIncidenceAngleFilePath))
1558 : {
1559 : CPLXMLTreeCloser psIncidenceAngle(
1560 16 : CPLParseXMLFile(osIncidenceAngleFilePath));
1561 :
1562 8 : int pixelFirstLutValue = atoi(
1563 8 : CPLGetXMLValue(psIncidenceAngle.get(),
1564 : "=incidenceAngles.pixelFirstAnglesValue", "0"));
1565 :
1566 8 : const int stepSize = atoi(CPLGetXMLValue(
1567 8 : psIncidenceAngle.get(), "=incidenceAngles.stepSize", "0"));
1568 : const int numberOfValues =
1569 8 : atoi(CPLGetXMLValue(psIncidenceAngle.get(),
1570 : "=incidenceAngles.numberOfValues", "0"));
1571 :
1572 8 : if (!(stepSize == 0 || stepSize == INT_MIN ||
1573 : numberOfValues == INT_MIN ||
1574 8 : abs(numberOfValues) > INT_MAX / abs(stepSize)))
1575 : {
1576 : /* Get the Pixel Per range */
1577 8 : const int tableSize = abs(stepSize) * abs(numberOfValues);
1578 :
1579 16 : CPLString angles;
1580 : // Loop through all nodes with spaces
1581 : CPLXMLNode *psNextNode =
1582 8 : CPLGetXMLNode(psIncidenceAngle.get(), "=incidenceAngles");
1583 :
1584 : CPLXMLNode *psNodeInc;
1585 458 : for (psNodeInc = psNextNode->psChild; psNodeInc != nullptr;
1586 450 : psNodeInc = psNodeInc->psNext)
1587 : {
1588 450 : if (EQUAL(psNodeInc->pszValue, "angles"))
1589 : {
1590 402 : if (angles.length() > 0)
1591 : {
1592 394 : angles.append(" "); /* separator */
1593 : }
1594 : const char *valAngle =
1595 402 : CPLGetXMLValue(psNodeInc, "", "");
1596 402 : angles.append(valAngle);
1597 : }
1598 : }
1599 :
1600 : char **papszAngleList =
1601 8 : CSLTokenizeString2(angles, " ", CSLT_HONOURSTRINGS);
1602 :
1603 : /* Allocate the right LUT size according to the product range pixel
1604 : */
1605 8 : poDS->m_IncidenceAngleTableSize = tableSize;
1606 16 : poDS->m_nfIncidenceAngleTable =
1607 8 : InterpolateValues(papszAngleList, tableSize, stepSize,
1608 : numberOfValues, pixelFirstLutValue);
1609 :
1610 8 : CSLDestroy(papszAngleList);
1611 : }
1612 : }
1613 : }
1614 :
1615 24 : for (int iPoleInx = 0; iPoleInx < nPolarizationsGridCount; iPoleInx++)
1616 : {
1617 : // Search for a specific band name
1618 : const CPLString pszPole =
1619 16 : CPLString(aosPolarizationsGrids[iPoleInx]).toupper();
1620 :
1621 : // Look if the NoiseLevel file xml exist for the
1622 16 : CPLXMLNode *psRefNode = psImageReferenceAttributes->psChild;
1623 232 : for (; psRefNode != nullptr; psRefNode = psRefNode->psNext)
1624 : {
1625 216 : if (EQUAL(psRefNode->pszValue, "noiseLevelFileName") && bCanCalib)
1626 : {
1627 : /* Determine which incidence angle correction this is */
1628 : const char *pszPoleToMatch =
1629 32 : CPLGetXMLValue(psRefNode, "pole", "");
1630 : const char *pszNoiseLevelFile =
1631 32 : CPLGetXMLValue(psRefNode, "", "");
1632 :
1633 32 : if (*pszPoleToMatch == '\0')
1634 16 : continue;
1635 :
1636 32 : if (EQUAL(pszPoleToMatch, "XC"))
1637 : /* Skip noise for RCM MLC's 3rd band XC */
1638 0 : continue;
1639 :
1640 32 : if (EQUAL(pszNoiseLevelFile, ""))
1641 0 : continue;
1642 :
1643 : /* --------------------------------------------------------------------
1644 : */
1645 : /* With RCM, LUT file is different per polarizarion (pole)
1646 : */
1647 : /* The following code make sure to loop through all
1648 : * possible */
1649 : /* 'noiseLevelFileName' and match the <ipdf 'pole' name */
1650 : /* --------------------------------------------------------------------
1651 : */
1652 32 : if (pszPole.compare(pszPoleToMatch) != 0)
1653 : {
1654 16 : continue;
1655 : }
1656 :
1657 : /* --------------------------------------------------------------------
1658 : */
1659 : /* LUT calibration is unique per pole in a sub-folder */
1660 : /* called 'calibration' from the 'metadata' folder */
1661 : /* --------------------------------------------------------------------
1662 : */
1663 :
1664 32 : const CPLString oNoiseLevelPath = CPLFormFilenameSafe(
1665 16 : CPLFormFilenameSafe(osPath, CALIBRATION_FOLDER, nullptr)
1666 : .c_str(),
1667 32 : pszNoiseLevelFile, nullptr);
1668 16 : if (IsValidXMLFile(oNoiseLevelPath))
1669 : {
1670 16 : osNoiseLevelsValues = oNoiseLevelPath;
1671 : }
1672 : }
1673 : }
1674 :
1675 : // Search again with different file
1676 16 : psRefNode = psImageReferenceAttributes->psChild;
1677 232 : for (; psRefNode != nullptr; psRefNode = psRefNode->psNext)
1678 : {
1679 216 : if (EQUAL(psRefNode->pszValue, "lookupTableFileName") && bCanCalib)
1680 : {
1681 : /* Determine which incidence angle correction this is */
1682 : const char *pszLUTType =
1683 102 : CPLGetXMLValue(psRefNode, "sarCalibrationType", "");
1684 : const char *pszPoleToMatch =
1685 102 : CPLGetXMLValue(psRefNode, "pole", "");
1686 102 : const char *pszLUTFile = CPLGetXMLValue(psRefNode, "", "");
1687 :
1688 102 : if (*pszPoleToMatch == '\0')
1689 54 : continue;
1690 :
1691 102 : if (EQUAL(pszPoleToMatch, "XC"))
1692 : /* Skip Calib for RCM MLC's 3rd band XC */
1693 6 : continue;
1694 :
1695 96 : if (*pszLUTType == '\0')
1696 0 : continue;
1697 :
1698 96 : if (EQUAL(pszLUTType, ""))
1699 0 : continue;
1700 :
1701 : /* --------------------------------------------------------------------
1702 : */
1703 : /* With RCM, LUT file is different per polarizarion (pole)
1704 : */
1705 : /* The following code make sure to loop through all
1706 : * possible */
1707 : /* 'lookupTableFileName' and match the <ipdf 'pole' name */
1708 : /* --------------------------------------------------------------------
1709 : */
1710 96 : if (pszPole.compare(pszPoleToMatch) != 0)
1711 : {
1712 48 : continue;
1713 : }
1714 :
1715 : /* --------------------------------------------------------------------
1716 : */
1717 : /* LUT calibration is unique per pole in a sub-folder */
1718 : /* called 'calibration' from the 'metadata' folder */
1719 : /* --------------------------------------------------------------------
1720 : */
1721 96 : const CPLString osLUTFilePath = CPLFormFilenameSafe(
1722 48 : CPLFormFilenameSafe(osPath, CALIBRATION_FOLDER, nullptr)
1723 : .c_str(),
1724 96 : pszLUTFile, nullptr);
1725 :
1726 64 : if (EQUAL(pszLUTType, "Beta Nought") &&
1727 16 : IsValidXMLFile(osLUTFilePath))
1728 : {
1729 32 : poDS->papszExtraFiles =
1730 16 : CSLAddString(poDS->papszExtraFiles, osLUTFilePath);
1731 :
1732 : CPLString pszBuf(
1733 16 : FormatCalibration(szBETA0, osMDFilename.c_str()));
1734 16 : CPLFree(pszBeta0LUT);
1735 16 : pszBeta0LUT = VSIStrdup(osLUTFilePath);
1736 :
1737 : const char *oldLut =
1738 16 : poDS->GetMetadataItem("BETA_NOUGHT_LUT");
1739 16 : if (oldLut == nullptr)
1740 : {
1741 8 : poDS->SetMetadataItem("BETA_NOUGHT_LUT", osLUTFilePath);
1742 : }
1743 : else
1744 : {
1745 : /* Keep adding LUT file for all bands, should be planty
1746 : * of space */
1747 : char *ptrConcatLut =
1748 8 : static_cast<char *>(CPLMalloc(2048));
1749 8 : ptrConcatLut[0] =
1750 : '\0'; /* Just initialize the first byte */
1751 8 : strcat(ptrConcatLut, oldLut);
1752 8 : strcat(ptrConcatLut, ",");
1753 8 : strcat(ptrConcatLut, osLUTFilePath);
1754 8 : poDS->SetMetadataItem("BETA_NOUGHT_LUT", ptrConcatLut);
1755 8 : CPLFree(ptrConcatLut);
1756 : }
1757 :
1758 32 : poDS->papszSubDatasets = CSLSetNameValue(
1759 16 : poDS->papszSubDatasets, "SUBDATASET_3_NAME", pszBuf);
1760 16 : poDS->papszSubDatasets = CSLSetNameValue(
1761 16 : poDS->papszSubDatasets, "SUBDATASET_3_DESC",
1762 : "Beta Nought calibrated");
1763 : }
1764 48 : else if (EQUAL(pszLUTType, "Sigma Nought") &&
1765 16 : IsValidXMLFile(osLUTFilePath))
1766 : {
1767 32 : poDS->papszExtraFiles =
1768 16 : CSLAddString(poDS->papszExtraFiles, osLUTFilePath);
1769 :
1770 : CPLString pszBuf(
1771 16 : FormatCalibration(szSIGMA0, osMDFilename.c_str()));
1772 16 : CPLFree(pszSigma0LUT);
1773 16 : pszSigma0LUT = VSIStrdup(osLUTFilePath);
1774 :
1775 : const char *oldLut =
1776 16 : poDS->GetMetadataItem("SIGMA_NOUGHT_LUT");
1777 16 : if (oldLut == nullptr)
1778 : {
1779 16 : poDS->SetMetadataItem("SIGMA_NOUGHT_LUT",
1780 8 : osLUTFilePath);
1781 : }
1782 : else
1783 : {
1784 : /* Keep adding LUT file for all bands, should be planty
1785 : * of space */
1786 : char *ptrConcatLut =
1787 8 : static_cast<char *>(CPLMalloc(2048));
1788 8 : ptrConcatLut[0] =
1789 : '\0'; /* Just initialize the first byte */
1790 8 : strcat(ptrConcatLut, oldLut);
1791 8 : strcat(ptrConcatLut, ",");
1792 8 : strcat(ptrConcatLut, osLUTFilePath);
1793 8 : poDS->SetMetadataItem("SIGMA_NOUGHT_LUT", ptrConcatLut);
1794 8 : CPLFree(ptrConcatLut);
1795 : }
1796 :
1797 32 : poDS->papszSubDatasets = CSLSetNameValue(
1798 16 : poDS->papszSubDatasets, "SUBDATASET_2_NAME", pszBuf);
1799 16 : poDS->papszSubDatasets = CSLSetNameValue(
1800 16 : poDS->papszSubDatasets, "SUBDATASET_2_DESC",
1801 : "Sigma Nought calibrated");
1802 : }
1803 32 : else if (EQUAL(pszLUTType, "Gamma") &&
1804 16 : IsValidXMLFile(osLUTFilePath))
1805 : {
1806 32 : poDS->papszExtraFiles =
1807 16 : CSLAddString(poDS->papszExtraFiles, osLUTFilePath);
1808 :
1809 : CPLString pszBuf(
1810 16 : FormatCalibration(szGAMMA, osMDFilename.c_str()));
1811 16 : CPLFree(pszGammaLUT);
1812 16 : pszGammaLUT = VSIStrdup(osLUTFilePath);
1813 :
1814 16 : const char *oldLut = poDS->GetMetadataItem("GAMMA_LUT");
1815 16 : if (oldLut == nullptr)
1816 : {
1817 8 : poDS->SetMetadataItem("GAMMA_LUT", osLUTFilePath);
1818 : }
1819 : else
1820 : {
1821 : /* Keep adding LUT file for all bands, should be planty
1822 : * of space */
1823 : char *ptrConcatLut =
1824 8 : static_cast<char *>(CPLMalloc(2048));
1825 8 : ptrConcatLut[0] =
1826 : '\0'; /* Just initialize the first byte */
1827 8 : strcat(ptrConcatLut, oldLut);
1828 8 : strcat(ptrConcatLut, ",");
1829 8 : strcat(ptrConcatLut, osLUTFilePath);
1830 8 : poDS->SetMetadataItem("GAMMA_LUT", ptrConcatLut);
1831 8 : CPLFree(ptrConcatLut);
1832 : }
1833 :
1834 32 : poDS->papszSubDatasets = CSLSetNameValue(
1835 16 : poDS->papszSubDatasets, "SUBDATASET_4_NAME", pszBuf);
1836 16 : poDS->papszSubDatasets = CSLSetNameValue(
1837 16 : poDS->papszSubDatasets, "SUBDATASET_4_DESC",
1838 : "Gamma calibrated");
1839 : }
1840 : }
1841 : }
1842 :
1843 : /* --------------------------------------------------------------------
1844 : */
1845 : /* Fetch ipdf image. Could be either tif or ntf. */
1846 : /* Replace / by \\ */
1847 : /* --------------------------------------------------------------------
1848 : */
1849 : const char *pszBasedFilename;
1850 16 : if (bIsNITF)
1851 : {
1852 0 : pszBasedFilename = pszNITF_Filename;
1853 : }
1854 : else
1855 : {
1856 16 : const int bandPositionIndex = imageBandList.FindString(pszPole);
1857 16 : if (bandPositionIndex < 0)
1858 : {
1859 0 : CPLFree(imageBandList);
1860 0 : CPLFree(imageBandFileList);
1861 :
1862 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1863 : "ERROR: RCM cannot find the polarization %s. Please "
1864 : "contact your data provider for a corrected dataset",
1865 : pszPole.c_str());
1866 0 : return nullptr;
1867 : }
1868 :
1869 16 : pszBasedFilename = imageBandFileList[bandPositionIndex];
1870 : }
1871 :
1872 : /* --------------------------------------------------------------------
1873 : */
1874 : /* Form full filename (path of product.xml + basename). */
1875 : /* --------------------------------------------------------------------
1876 : */
1877 16 : char *pszFullname = CPLStrdup(
1878 32 : CPLFormFilenameSafe(osPath, pszBasedFilename, nullptr).c_str());
1879 :
1880 : /* --------------------------------------------------------------------
1881 : */
1882 : /* Try and open the file. */
1883 : /* --------------------------------------------------------------------
1884 : */
1885 : auto poBandFile = std::unique_ptr<GDALDataset>(GDALDataset::Open(
1886 16 : pszFullname, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR));
1887 16 : if (poBandFile == nullptr)
1888 : {
1889 0 : CPLFree(pszFullname);
1890 0 : continue;
1891 : }
1892 16 : if (poBandFile->GetRasterCount() == 0)
1893 : {
1894 0 : CPLFree(pszFullname);
1895 0 : continue;
1896 : }
1897 :
1898 32 : poDS->papszExtraFiles =
1899 16 : CSLAddString(poDS->papszExtraFiles, pszFullname);
1900 :
1901 : /* Some CFloat32 NITF files have nBitsPerSample incorrectly reported */
1902 : /* as 16, and get misinterpreted as CInt16. Check the underlying NITF
1903 : */
1904 : /* and override if this is the case. */
1905 16 : if (poBandFile->GetRasterBand(1)->GetRasterDataType() == GDT_CFloat32)
1906 0 : eDataType = GDT_CFloat32;
1907 :
1908 : BandMappingRCM b =
1909 16 : checkBandFileMappingRCM(eDataType, poBandFile.get(), bIsNITF);
1910 16 : if (b == BANDERROR)
1911 : {
1912 0 : CPLFree(pszFullname);
1913 0 : CPLError(CE_Failure, CPLE_AppDefined,
1914 : "The underlying band files do not have an appropriate "
1915 : "data type.");
1916 0 : return nullptr;
1917 : }
1918 16 : bool twoBandComplex = b == TWOBANDCOMPLEX;
1919 16 : bool isOneFilePerPol = (imageBandCount == imageBandFileCount);
1920 :
1921 : /* --------------------------------------------------------------------
1922 : */
1923 : /* Create the band. */
1924 : /* --------------------------------------------------------------------
1925 : */
1926 16 : int bandNum = poDS->GetRasterCount() + 1;
1927 16 : if (eCalib == None || eCalib == Uncalib)
1928 : {
1929 : RCMRasterBand *poBand = new RCMRasterBand(
1930 20 : poDS.get(), bandNum, eDataType, pszPole, poBandFile.release(),
1931 20 : twoBandComplex, isOneFilePerPol, bIsNITF);
1932 :
1933 10 : poDS->SetBand(poDS->GetRasterCount() + 1, poBand);
1934 : }
1935 : else
1936 : {
1937 6 : const char *pszLUT = nullptr;
1938 6 : switch (eCalib)
1939 : {
1940 2 : case Sigma0:
1941 2 : pszLUT = pszSigma0LUT;
1942 2 : break;
1943 2 : case Beta0:
1944 2 : pszLUT = pszBeta0LUT;
1945 2 : break;
1946 2 : case Gamma:
1947 2 : pszLUT = pszGammaLUT;
1948 2 : break;
1949 0 : default:
1950 : /* we should bomb gracefully... */
1951 0 : pszLUT = pszSigma0LUT;
1952 : }
1953 6 : if (!pszLUT)
1954 : {
1955 0 : CPLFree(pszFullname);
1956 0 : CPLError(CE_Failure, CPLE_AppDefined, "LUT missing.");
1957 0 : return nullptr;
1958 : }
1959 :
1960 : // The variable 'osNoiseLevelsValues' is always the same for a ban
1961 : // name except the XML contains different calibration name
1962 6 : if (poDS->isComplexData)
1963 : {
1964 : // If Complex, always 32 bits
1965 : RCMCalibRasterBand *poBand = new RCMCalibRasterBand(
1966 0 : poDS.get(), pszPole, GDT_Float32, poBandFile.release(),
1967 0 : eCalib, pszLUT, osNoiseLevelsValues.c_str(), eDataType);
1968 0 : poDS->SetBand(poDS->GetRasterCount() + 1, poBand);
1969 : }
1970 : else
1971 : {
1972 : // Whatever the datatype was previoulsy set
1973 : RCMCalibRasterBand *poBand = new RCMCalibRasterBand(
1974 12 : poDS.get(), pszPole, eDataType, poBandFile.release(),
1975 12 : eCalib, pszLUT, osNoiseLevelsValues.c_str(), eDataType);
1976 6 : poDS->SetBand(poDS->GetRasterCount() + 1, poBand);
1977 : }
1978 : }
1979 :
1980 16 : CPLFree(pszFullname);
1981 : }
1982 :
1983 8 : if (poDS->papszSubDatasets != nullptr && eCalib == None)
1984 : {
1985 4 : CPLString pszBuf = FormatCalibration(szUNCALIB, osMDFilename.c_str());
1986 4 : poDS->papszSubDatasets = CSLSetNameValue(poDS->papszSubDatasets,
1987 : "SUBDATASET_1_NAME", pszBuf);
1988 4 : poDS->papszSubDatasets =
1989 4 : CSLSetNameValue(poDS->papszSubDatasets, "SUBDATASET_1_DESC",
1990 : "Uncalibrated digital numbers");
1991 : }
1992 4 : else if (poDS->papszSubDatasets != nullptr)
1993 : {
1994 4 : CSLDestroy(poDS->papszSubDatasets);
1995 4 : poDS->papszSubDatasets = nullptr;
1996 : }
1997 :
1998 : /* -------------------------------------------------------------------- */
1999 : /* Set the appropriate MATRIX_REPRESENTATION. */
2000 : /* -------------------------------------------------------------------- */
2001 :
2002 8 : if (poDS->GetRasterCount() == 4 &&
2003 0 : (eDataType == GDT_CInt16 || eDataType == GDT_CFloat32))
2004 : {
2005 0 : poDS->SetMetadataItem("MATRIX_REPRESENTATION", "SCATTERING");
2006 : }
2007 :
2008 : /* -------------------------------------------------------------------- */
2009 : /* Collect a few useful metadata items. */
2010 : /* -------------------------------------------------------------------- */
2011 :
2012 8 : pszItem = CPLGetXMLValue(psSourceAttrs, "satellite", "");
2013 8 : poDS->SetMetadataItem("SATELLITE_IDENTIFIER", pszItem);
2014 :
2015 8 : pszItem = CPLGetXMLValue(psSourceAttrs, "sensor", "");
2016 8 : poDS->SetMetadataItem("SENSOR_IDENTIFIER", pszItem);
2017 :
2018 : /* Get beam mode mnemonic */
2019 8 : pszItem = CPLGetXMLValue(psSourceAttrs, "beamMode", "UNK");
2020 8 : poDS->SetMetadataItem("BEAM_MODE", pszItem);
2021 :
2022 8 : pszItem = CPLGetXMLValue(psSourceAttrs, "beamModeMnemonic", "UNK");
2023 8 : poDS->SetMetadataItem("BEAM_MODE_MNEMONIC", pszItem);
2024 :
2025 8 : pszItem = CPLGetXMLValue(psSourceAttrs, "beamModeDefinitionId", "UNK");
2026 8 : poDS->SetMetadataItem("BEAM_MODE_DEFINITION_ID", pszItem);
2027 :
2028 8 : pszItem = CPLGetXMLValue(psSourceAttrs, "rawDataStartTime", "UNK");
2029 8 : poDS->SetMetadataItem("ACQUISITION_START_TIME", pszItem);
2030 :
2031 8 : pszItem = CPLGetXMLValue(psSourceAttrs, "inputDatasetFacilityId", "UNK");
2032 8 : poDS->SetMetadataItem("FACILITY_IDENTIFIER", pszItem);
2033 :
2034 8 : pszItem = CPLGetXMLValue(psSourceAttrs,
2035 : "orbitAndAttitude.orbitInformation.passDirection",
2036 : "UNK");
2037 8 : poDS->SetMetadataItem("ORBIT_DIRECTION", pszItem);
2038 8 : pszItem = CPLGetXMLValue(
2039 : psSourceAttrs, "orbitAndAttitude.orbitInformation.orbitDataSource",
2040 : "UNK");
2041 8 : poDS->SetMetadataItem("ORBIT_DATA_SOURCE", pszItem);
2042 8 : pszItem = CPLGetXMLValue(
2043 : psSourceAttrs, "orbitAndAttitude.orbitInformation.orbitDataFileName",
2044 : "UNK");
2045 8 : poDS->SetMetadataItem("ORBIT_DATA_FILE", pszItem);
2046 :
2047 : /* Get incidence angle information. DONE */
2048 8 : pszItem = CPLGetXMLValue(psSceneAttributes, "imageAttributes.incAngNearRng",
2049 : "UNK");
2050 8 : poDS->SetMetadataItem("NEAR_RANGE_INCIDENCE_ANGLE", pszItem);
2051 :
2052 8 : pszItem = CPLGetXMLValue(psSceneAttributes, "imageAttributes.incAngFarRng",
2053 : "UNK");
2054 8 : poDS->SetMetadataItem("FAR_RANGE_INCIDENCE_ANGLE", pszItem);
2055 :
2056 8 : pszItem = CPLGetXMLValue(psSceneAttributes,
2057 : "imageAttributes.slantRangeNearEdge", "UNK");
2058 8 : poDS->SetMetadataItem("SLANT_RANGE_NEAR_EDGE", pszItem);
2059 :
2060 8 : pszItem = CPLGetXMLValue(psSceneAttributes,
2061 : "imageAttributes.slantRangeFarEdge", "UNK");
2062 8 : poDS->SetMetadataItem("SLANT_RANGE_FAR_EDGE", pszItem);
2063 :
2064 : /*--------------------------------------------------------------------- */
2065 : /* Collect Map projection/Geotransform information, if present.DONE */
2066 : /* In RCM, there is no file that indicates */
2067 : /* -------------------------------------------------------------------- */
2068 8 : CPLXMLNode *psMapProjection = CPLGetXMLNode(
2069 : psImageReferenceAttributes, "geographicInformation.mapProjection");
2070 :
2071 8 : if (psMapProjection != nullptr)
2072 : {
2073 : CPLXMLNode *psPos =
2074 0 : CPLGetXMLNode(psMapProjection, "positioningInformation");
2075 :
2076 : pszItem =
2077 0 : CPLGetXMLValue(psMapProjection, "mapProjectionDescriptor", "UNK");
2078 0 : poDS->SetMetadataItem("MAP_PROJECTION_DESCRIPTOR", pszItem);
2079 :
2080 : pszItem =
2081 0 : CPLGetXMLValue(psMapProjection, "mapProjectionOrientation", "UNK");
2082 0 : poDS->SetMetadataItem("MAP_PROJECTION_ORIENTATION", pszItem);
2083 :
2084 0 : pszItem = CPLGetXMLValue(psMapProjection, "resamplingKernel", "UNK");
2085 0 : poDS->SetMetadataItem("RESAMPLING_KERNEL", pszItem);
2086 :
2087 0 : pszItem = CPLGetXMLValue(psMapProjection, "satelliteHeading", "UNK");
2088 0 : poDS->SetMetadataItem("SATELLITE_HEADING", pszItem);
2089 :
2090 0 : if (psPos != nullptr)
2091 : {
2092 0 : const double tl_x = CPLStrtod(
2093 : CPLGetXMLValue(psPos, "upperLeftCorner.mapCoordinate.easting",
2094 : "0.0"),
2095 : nullptr);
2096 0 : const double tl_y = CPLStrtod(
2097 : CPLGetXMLValue(psPos, "upperLeftCorner.mapCoordinate.northing",
2098 : "0.0"),
2099 : nullptr);
2100 0 : const double bl_x = CPLStrtod(
2101 : CPLGetXMLValue(psPos, "lowerLeftCorner.mapCoordinate.easting",
2102 : "0.0"),
2103 : nullptr);
2104 0 : const double bl_y = CPLStrtod(
2105 : CPLGetXMLValue(psPos, "lowerLeftCorner.mapCoordinate.northing",
2106 : "0.0"),
2107 : nullptr);
2108 0 : const double tr_x = CPLStrtod(
2109 : CPLGetXMLValue(psPos, "upperRightCorner.mapCoordinate.easting",
2110 : "0.0"),
2111 : nullptr);
2112 0 : const double tr_y = CPLStrtod(
2113 : CPLGetXMLValue(psPos, "upperRightCorner.mapCoordinate.northing",
2114 : "0.0"),
2115 : nullptr);
2116 0 : poDS->adfGeoTransform[1] = (tr_x - tl_x) / (poDS->nRasterXSize - 1);
2117 0 : poDS->adfGeoTransform[4] = (tr_y - tl_y) / (poDS->nRasterXSize - 1);
2118 0 : poDS->adfGeoTransform[2] = (bl_x - tl_x) / (poDS->nRasterYSize - 1);
2119 0 : poDS->adfGeoTransform[5] = (bl_y - tl_y) / (poDS->nRasterYSize - 1);
2120 0 : poDS->adfGeoTransform[0] = (tl_x - 0.5 * poDS->adfGeoTransform[1] -
2121 0 : 0.5 * poDS->adfGeoTransform[2]);
2122 0 : poDS->adfGeoTransform[3] = (tl_y - 0.5 * poDS->adfGeoTransform[4] -
2123 0 : 0.5 * poDS->adfGeoTransform[5]);
2124 :
2125 : /* Use bottom right pixel to test geotransform */
2126 0 : const double br_x = CPLStrtod(
2127 : CPLGetXMLValue(psPos, "lowerRightCorner.mapCoordinate.easting",
2128 : "0.0"),
2129 : nullptr);
2130 0 : const double br_y = CPLStrtod(
2131 : CPLGetXMLValue(psPos, "lowerRightCorner.mapCoordinate.northing",
2132 : "0.0"),
2133 : nullptr);
2134 : const double testx =
2135 0 : poDS->adfGeoTransform[0] +
2136 0 : poDS->adfGeoTransform[1] * (poDS->nRasterXSize - 0.5) +
2137 0 : poDS->adfGeoTransform[2] * (poDS->nRasterYSize - 0.5);
2138 : const double testy =
2139 0 : poDS->adfGeoTransform[3] +
2140 0 : poDS->adfGeoTransform[4] * (poDS->nRasterXSize - 0.5) +
2141 0 : poDS->adfGeoTransform[5] * (poDS->nRasterYSize - 0.5);
2142 :
2143 : /* Give 1/4 pixel numerical error leeway in calculating location
2144 : based on affine transform */
2145 0 : if ((fabs(testx - br_x) >
2146 0 : fabs(0.25 *
2147 0 : (poDS->adfGeoTransform[1] + poDS->adfGeoTransform[2]))) ||
2148 0 : (fabs(testy - br_y) > fabs(0.25 * (poDS->adfGeoTransform[4] +
2149 0 : poDS->adfGeoTransform[5]))))
2150 : {
2151 0 : CPLError(CE_Warning, CPLE_AppDefined,
2152 : "WARNING: Unexpected error in calculating affine "
2153 : "transform: corner coordinates inconsistent.");
2154 : }
2155 : else
2156 : {
2157 0 : poDS->bHaveGeoTransform = TRUE;
2158 : }
2159 : }
2160 : }
2161 :
2162 : /* -------------------------------------------------------------------- */
2163 : /* Collect Projection String Information.DONE */
2164 : /* -------------------------------------------------------------------- */
2165 : CPLXMLNode *psEllipsoid =
2166 8 : CPLGetXMLNode(psImageReferenceAttributes,
2167 : "geographicInformation.ellipsoidParameters");
2168 :
2169 8 : if (psEllipsoid != nullptr)
2170 : {
2171 16 : OGRSpatialReference oLL, oPrj;
2172 :
2173 : const char *pszGeodeticTerrainHeight =
2174 8 : CPLGetXMLValue(psEllipsoid, "geodeticTerrainHeight", "UNK");
2175 8 : poDS->SetMetadataItem("GEODETIC_TERRAIN_HEIGHT",
2176 8 : pszGeodeticTerrainHeight);
2177 :
2178 : const char *pszEllipsoidName =
2179 8 : CPLGetXMLValue(psEllipsoid, "ellipsoidName", "");
2180 : double minor_axis =
2181 8 : CPLAtof(CPLGetXMLValue(psEllipsoid, "semiMinorAxis", "0.0"));
2182 : double major_axis =
2183 8 : CPLAtof(CPLGetXMLValue(psEllipsoid, "semiMajorAxis", "0.0"));
2184 :
2185 8 : if (EQUAL(pszEllipsoidName, "") || (minor_axis == 0.0) ||
2186 : (major_axis == 0.0))
2187 : {
2188 0 : oLL.SetWellKnownGeogCS("WGS84");
2189 0 : oPrj.SetWellKnownGeogCS("WGS84");
2190 :
2191 0 : CPLError(CE_Warning, CPLE_AppDefined,
2192 : "WARNING: Incomplete ellipsoid "
2193 : "information. Using wgs-84 parameters.");
2194 : }
2195 8 : else if (EQUAL(pszEllipsoidName, "WGS84") ||
2196 8 : EQUAL(pszEllipsoidName, "WGS 1984"))
2197 : {
2198 8 : oLL.SetWellKnownGeogCS("WGS84");
2199 8 : oPrj.SetWellKnownGeogCS("WGS84");
2200 : }
2201 : else
2202 : {
2203 0 : const double inv_flattening =
2204 0 : major_axis / (major_axis - minor_axis);
2205 0 : oLL.SetGeogCS("", "", pszEllipsoidName, major_axis, inv_flattening);
2206 0 : oPrj.SetGeogCS("", "", pszEllipsoidName, major_axis,
2207 : inv_flattening);
2208 : }
2209 :
2210 8 : if (psMapProjection != nullptr)
2211 : {
2212 : const char *pszProj =
2213 0 : CPLGetXMLValue(psMapProjection, "mapProjectionDescriptor", "");
2214 0 : bool bUseProjInfo = false;
2215 :
2216 : CPLXMLNode *psUtmParams =
2217 0 : CPLGetXMLNode(psMapProjection, "utmProjectionParameters");
2218 :
2219 : CPLXMLNode *psNspParams =
2220 0 : CPLGetXMLNode(psMapProjection, "nspProjectionParameters");
2221 :
2222 0 : if ((psUtmParams != nullptr) && poDS->bHaveGeoTransform)
2223 : {
2224 : /* double origEasting, origNorthing; */
2225 0 : bool bNorth = true;
2226 :
2227 : const int utmZone =
2228 0 : atoi(CPLGetXMLValue(psUtmParams, "utmZone", ""));
2229 : const char *pszHemisphere =
2230 0 : CPLGetXMLValue(psUtmParams, "hemisphere", "");
2231 : #if 0
2232 : origEasting = CPLStrtod(CPLGetXMLValue(
2233 : psUtmParams, "mapOriginFalseEasting", "0.0"), nullptr);
2234 : origNorthing = CPLStrtod(CPLGetXMLValue(
2235 : psUtmParams, "mapOriginFalseNorthing", "0.0"), nullptr);
2236 : #endif
2237 0 : if (STARTS_WITH_CI(pszHemisphere, "southern"))
2238 0 : bNorth = FALSE;
2239 :
2240 0 : if (STARTS_WITH_CI(pszProj, "UTM"))
2241 : {
2242 0 : oPrj.SetUTM(utmZone, bNorth);
2243 0 : bUseProjInfo = true;
2244 : }
2245 : }
2246 0 : else if ((psNspParams != nullptr) && poDS->bHaveGeoTransform)
2247 : {
2248 0 : const double origEasting = CPLStrtod(
2249 : CPLGetXMLValue(psNspParams, "mapOriginFalseEasting", "0.0"),
2250 : nullptr);
2251 : const double origNorthing =
2252 0 : CPLStrtod(CPLGetXMLValue(psNspParams,
2253 : "mapOriginFalseNorthing", "0.0"),
2254 : nullptr);
2255 0 : const double copLong = CPLStrtod(
2256 : CPLGetXMLValue(psNspParams, "centerOfProjectionLongitude",
2257 : "0.0"),
2258 : nullptr);
2259 0 : const double copLat = CPLStrtod(
2260 : CPLGetXMLValue(psNspParams, "centerOfProjectionLatitude",
2261 : "0.0"),
2262 : nullptr);
2263 0 : const double sP1 = CPLStrtod(
2264 : CPLGetXMLValue(psNspParams, "standardParallels1", "0.0"),
2265 : nullptr);
2266 0 : const double sP2 = CPLStrtod(
2267 : CPLGetXMLValue(psNspParams, "standardParallels2", "0.0"),
2268 : nullptr);
2269 :
2270 0 : if (STARTS_WITH_CI(pszProj, "ARC"))
2271 : {
2272 : /* Albers Conical Equal Area */
2273 0 : oPrj.SetACEA(sP1, sP2, copLat, copLong, origEasting,
2274 : origNorthing);
2275 0 : bUseProjInfo = true;
2276 : }
2277 0 : else if (STARTS_WITH_CI(pszProj, "LCC"))
2278 : {
2279 : /* Lambert Conformal Conic */
2280 0 : oPrj.SetLCC(sP1, sP2, copLat, copLong, origEasting,
2281 : origNorthing);
2282 0 : bUseProjInfo = true;
2283 : }
2284 0 : else if (STARTS_WITH_CI(pszProj, "STPL"))
2285 : {
2286 : /* State Plate
2287 : ASSUMPTIONS: "zone" in product.xml matches USGS
2288 : definition as expected by ogr spatial reference; NAD83
2289 : zones (versus NAD27) are assumed. */
2290 :
2291 : const int nSPZone =
2292 0 : atoi(CPLGetXMLValue(psNspParams, "zone", "1"));
2293 :
2294 0 : oPrj.SetStatePlane(nSPZone, TRUE, nullptr, 0.0);
2295 0 : bUseProjInfo = true;
2296 : }
2297 : }
2298 :
2299 0 : if (bUseProjInfo)
2300 : {
2301 0 : poDS->m_oSRS = std::move(oPrj);
2302 : }
2303 : else
2304 : {
2305 0 : CPLError(CE_Warning, CPLE_AppDefined,
2306 : "WARNING: Unable to interpret projection information; "
2307 : "check mapProjection info in product.xml!");
2308 : }
2309 : }
2310 :
2311 8 : poDS->m_oGCPSRS = std::move(oLL);
2312 : }
2313 :
2314 : /* -------------------------------------------------------------------- */
2315 : /* Collect GCPs.DONE */
2316 : /* -------------------------------------------------------------------- */
2317 8 : CPLXMLNode *psGeoGrid = CPLGetXMLNode(
2318 : psImageReferenceAttributes, "geographicInformation.geolocationGrid");
2319 :
2320 8 : if (psGeoGrid != nullptr)
2321 : {
2322 : /* count GCPs */
2323 8 : poDS->nGCPCount = 0;
2324 :
2325 136 : for (psNode = psGeoGrid->psChild; psNode != nullptr;
2326 128 : psNode = psNode->psNext)
2327 : {
2328 128 : if (EQUAL(psNode->pszValue, "imageTiePoint"))
2329 128 : poDS->nGCPCount++;
2330 : }
2331 :
2332 16 : poDS->pasGCPList = reinterpret_cast<GDAL_GCP *>(
2333 8 : CPLCalloc(sizeof(GDAL_GCP), poDS->nGCPCount));
2334 :
2335 8 : poDS->nGCPCount = 0;
2336 :
2337 136 : for (psNode = psGeoGrid->psChild; psNode != nullptr;
2338 128 : psNode = psNode->psNext)
2339 : {
2340 128 : GDAL_GCP *psGCP = poDS->pasGCPList + poDS->nGCPCount;
2341 :
2342 128 : if (!EQUAL(psNode->pszValue, "imageTiePoint"))
2343 0 : continue;
2344 :
2345 128 : poDS->nGCPCount++;
2346 :
2347 : char szID[32];
2348 128 : snprintf(szID, sizeof(szID), "%d", poDS->nGCPCount);
2349 128 : psGCP->pszId = CPLStrdup(szID);
2350 128 : psGCP->pszInfo = CPLStrdup("");
2351 128 : psGCP->dfGCPPixel =
2352 128 : CPLAtof(CPLGetXMLValue(psNode, "imageCoordinate.pixel", "0"));
2353 128 : psGCP->dfGCPLine =
2354 128 : CPLAtof(CPLGetXMLValue(psNode, "imageCoordinate.line", "0"));
2355 128 : psGCP->dfGCPX = CPLAtof(
2356 : CPLGetXMLValue(psNode, "geodeticCoordinate.longitude", ""));
2357 128 : psGCP->dfGCPY = CPLAtof(
2358 : CPLGetXMLValue(psNode, "geodeticCoordinate.latitude", ""));
2359 128 : psGCP->dfGCPZ = CPLAtof(
2360 : CPLGetXMLValue(psNode, "geodeticCoordinate.height", ""));
2361 : }
2362 : }
2363 :
2364 8 : if (pszBeta0LUT)
2365 8 : CPLFree(pszBeta0LUT);
2366 8 : if (pszSigma0LUT)
2367 8 : CPLFree(pszSigma0LUT);
2368 8 : if (pszGammaLUT)
2369 8 : CPLFree(pszGammaLUT);
2370 :
2371 : /* -------------------------------------------------------------------- */
2372 : /* Collect RPC.DONE */
2373 : /* -------------------------------------------------------------------- */
2374 8 : CPLXMLNode *psRationalFunctions = CPLGetXMLNode(
2375 : psImageReferenceAttributes, "geographicInformation.rationalFunctions");
2376 8 : if (psRationalFunctions != nullptr)
2377 : {
2378 8 : char **papszRPC = nullptr;
2379 : static const char *const apszXMLToGDALMapping[] = {
2380 : "biasError",
2381 : "ERR_BIAS",
2382 : "randomError",
2383 : "ERR_RAND",
2384 : //"lineFitQuality", "????",
2385 : //"pixelFitQuality", "????",
2386 : "lineOffset",
2387 : "LINE_OFF",
2388 : "pixelOffset",
2389 : "SAMP_OFF",
2390 : "latitudeOffset",
2391 : "LAT_OFF",
2392 : "longitudeOffset",
2393 : "LONG_OFF",
2394 : "heightOffset",
2395 : "HEIGHT_OFF",
2396 : "lineScale",
2397 : "LINE_SCALE",
2398 : "pixelScale",
2399 : "SAMP_SCALE",
2400 : "latitudeScale",
2401 : "LAT_SCALE",
2402 : "longitudeScale",
2403 : "LONG_SCALE",
2404 : "heightScale",
2405 : "HEIGHT_SCALE",
2406 : "lineNumeratorCoefficients",
2407 : "LINE_NUM_COEFF",
2408 : "lineDenominatorCoefficients",
2409 : "LINE_DEN_COEFF",
2410 : "pixelNumeratorCoefficients",
2411 : "SAMP_NUM_COEFF",
2412 : "pixelDenominatorCoefficients",
2413 : "SAMP_DEN_COEFF",
2414 : };
2415 136 : for (size_t i = 0; i < CPL_ARRAYSIZE(apszXMLToGDALMapping); i += 2)
2416 : {
2417 256 : const char *pszValue = CPLGetXMLValue(
2418 128 : psRationalFunctions, apszXMLToGDALMapping[i], nullptr);
2419 128 : if (pszValue)
2420 128 : papszRPC = CSLSetNameValue(
2421 128 : papszRPC, apszXMLToGDALMapping[i + 1], pszValue);
2422 : }
2423 8 : poDS->GDALDataset::SetMetadata(papszRPC, "RPC");
2424 8 : CSLDestroy(papszRPC);
2425 : }
2426 :
2427 : /* -------------------------------------------------------------------- */
2428 : /* Initialize any PAM information. */
2429 : /* -------------------------------------------------------------------- */
2430 16 : CPLString osDescription;
2431 16 : CPLString osSubdatasetName;
2432 8 : bool useSubdatasets = true;
2433 :
2434 8 : switch (eCalib)
2435 : {
2436 1 : case Sigma0:
2437 : {
2438 1 : osSubdatasetName = szSIGMA0;
2439 1 : osDescription = FormatCalibration(szSIGMA0, osMDFilename.c_str());
2440 : }
2441 1 : break;
2442 1 : case Beta0:
2443 : {
2444 1 : osSubdatasetName = szBETA0;
2445 1 : osDescription = FormatCalibration(szBETA0, osMDFilename.c_str());
2446 : }
2447 1 : break;
2448 1 : case Gamma:
2449 : {
2450 1 : osSubdatasetName = szGAMMA;
2451 1 : osDescription = FormatCalibration(szGAMMA, osMDFilename.c_str());
2452 : }
2453 1 : break;
2454 1 : case Uncalib:
2455 : {
2456 1 : osSubdatasetName = szUNCALIB;
2457 1 : osDescription = FormatCalibration(szUNCALIB, osMDFilename.c_str());
2458 : }
2459 1 : break;
2460 4 : default:
2461 4 : osSubdatasetName = szUNCALIB;
2462 4 : osDescription = osMDFilename;
2463 4 : useSubdatasets = false;
2464 : }
2465 :
2466 8 : CPL_IGNORE_RET_VAL(osSubdatasetName);
2467 :
2468 8 : if (eCalib != None)
2469 4 : poDS->papszExtraFiles =
2470 4 : CSLAddString(poDS->papszExtraFiles, osMDFilename);
2471 :
2472 : /* -------------------------------------------------------------------- */
2473 : /* Initialize any PAM information. */
2474 : /* -------------------------------------------------------------------- */
2475 8 : poDS->SetDescription(osDescription);
2476 :
2477 8 : poDS->SetPhysicalFilename(osMDFilename);
2478 8 : poDS->SetSubdatasetName(osDescription);
2479 :
2480 8 : poDS->TryLoadXML();
2481 :
2482 : /* -------------------------------------------------------------------- */
2483 : /* Check for overviews. */
2484 : /* -------------------------------------------------------------------- */
2485 8 : if (useSubdatasets)
2486 4 : poDS->oOvManager.Initialize(poDS.get(), ":::VIRTUAL:::");
2487 : else
2488 4 : poDS->oOvManager.Initialize(poDS.get(), osMDFilename);
2489 :
2490 8 : return poDS.release();
2491 : }
2492 :
2493 : /************************************************************************/
2494 : /* GetGCPCount() */
2495 : /************************************************************************/
2496 :
2497 6 : int RCMDataset::GetGCPCount()
2498 :
2499 : {
2500 6 : return nGCPCount;
2501 : }
2502 :
2503 : /************************************************************************/
2504 : /* GetGCPSpatialRef() */
2505 : /************************************************************************/
2506 :
2507 1 : const OGRSpatialReference *RCMDataset::GetGCPSpatialRef() const
2508 :
2509 : {
2510 1 : return m_oGCPSRS.IsEmpty() || nGCPCount == 0 ? nullptr : &m_oGCPSRS;
2511 : }
2512 :
2513 : /************************************************************************/
2514 : /* GetGCPs() */
2515 : /************************************************************************/
2516 :
2517 5 : const GDAL_GCP *RCMDataset::GetGCPs()
2518 :
2519 : {
2520 5 : return pasGCPList;
2521 : }
2522 :
2523 : /************************************************************************/
2524 : /* GetSpatialRef() */
2525 : /************************************************************************/
2526 :
2527 0 : const OGRSpatialReference *RCMDataset::GetSpatialRef() const
2528 :
2529 : {
2530 0 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
2531 : }
2532 :
2533 : /************************************************************************/
2534 : /* GetGeoTransform() */
2535 : /************************************************************************/
2536 :
2537 0 : CPLErr RCMDataset::GetGeoTransform(double *padfTransform)
2538 :
2539 : {
2540 0 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
2541 :
2542 0 : if (bHaveGeoTransform)
2543 : {
2544 0 : return CE_None;
2545 : }
2546 :
2547 0 : return CE_Failure;
2548 : }
2549 :
2550 : /************************************************************************/
2551 : /* GetMetadataDomainList() */
2552 : /************************************************************************/
2553 :
2554 0 : char **RCMDataset::GetMetadataDomainList()
2555 : {
2556 0 : return BuildMetadataDomainList(GDALDataset::GetMetadataDomainList(), TRUE,
2557 0 : "SUBDATASETS", nullptr);
2558 : }
2559 :
2560 : /************************************************************************/
2561 : /* GetMetadata() */
2562 : /************************************************************************/
2563 :
2564 2 : char **RCMDataset::GetMetadata(const char *pszDomain)
2565 :
2566 : {
2567 2 : if (pszDomain != nullptr && STARTS_WITH_CI(pszDomain, "SUBDATASETS") &&
2568 0 : papszSubDatasets != nullptr)
2569 0 : return papszSubDatasets;
2570 :
2571 2 : return GDALDataset::GetMetadata(pszDomain);
2572 : }
2573 :
2574 : /************************************************************************/
2575 : /* GDALRegister_RCM() */
2576 : /************************************************************************/
2577 :
2578 1682 : void GDALRegister_RCM()
2579 :
2580 : {
2581 1682 : if (GDALGetDriverByName("RCM") != nullptr)
2582 : {
2583 301 : return;
2584 : }
2585 :
2586 1381 : GDALDriver *poDriver = new GDALDriver();
2587 1381 : RCMDriverSetCommonMetadata(poDriver);
2588 :
2589 1381 : poDriver->pfnOpen = RCMDataset::Open;
2590 :
2591 1381 : GetGDALDriverManager()->RegisterDriver(poDriver);
2592 : }
|
