Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Sentinel SAFE products
4 : * Purpose: Sentinel Products (manifest.safe) driver
5 : * Author: Delfim Rego, delfimrego@gmail.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2015, Delfim Rego <delfimrego@gmail.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "safedataset.h"
14 :
15 : #include "cpl_time.h"
16 :
17 : #ifdef USE_OMP
18 : #include <omp.h>
19 : #endif
20 :
21 : #ifdef USE_OMP
22 : static int GetNumThreadsToUse()
23 : {
24 : unsigned int nCores = std::thread::hardware_concurrency();
25 : return (nCores / 2 > 1) ? nCores / 2 : 1;
26 : }
27 : #endif
28 :
29 : /************************************************************************/
30 : /* SAFERasterBand */
31 : /************************************************************************/
32 :
33 14 : SAFERasterBand::SAFERasterBand(SAFEDataset *poDSIn, GDALDataType eDataTypeIn,
34 : const CPLString &osSwath,
35 : const CPLString &osPolarization,
36 14 : std::unique_ptr<GDALDataset> &&poBandFileIn)
37 14 : : poBandFile(std::move(poBandFileIn))
38 : {
39 14 : poDS = poDSIn;
40 14 : GDALRasterBand *poSrcBand = poBandFile->GetRasterBand(1);
41 14 : poSrcBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
42 14 : eDataType = eDataTypeIn;
43 :
44 14 : if (!osSwath.empty())
45 14 : SetMetadataItem("SWATH", osSwath.c_str());
46 :
47 14 : if (!osPolarization.empty())
48 14 : SetMetadataItem("POLARIZATION", osPolarization.c_str());
49 14 : }
50 :
51 : /************************************************************************/
52 : /* IReadBlock() */
53 : /************************************************************************/
54 :
55 55 : CPLErr SAFERasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage)
56 :
57 : {
58 : /* -------------------------------------------------------------------- */
59 : /* If the last strip is partial, we need to avoid */
60 : /* over-requesting. We also need to initialize the extra part */
61 : /* of the block to zero. */
62 : /* -------------------------------------------------------------------- */
63 : int nRequestYSize;
64 55 : if ((nBlockYOff + 1) * nBlockYSize > nRasterYSize)
65 : {
66 5 : nRequestYSize = nRasterYSize - nBlockYOff * nBlockYSize;
67 5 : memset(pImage, 0,
68 5 : static_cast<size_t>(GDALGetDataTypeSizeBytes(eDataType)) *
69 5 : nBlockXSize * nBlockYSize);
70 : }
71 : else
72 : {
73 50 : nRequestYSize = nBlockYSize;
74 : }
75 :
76 : /*-------------------------------------------------------------------- */
77 : /* If the input imagery is tiled, also need to avoid over- */
78 : /* requesting in the X-direction. */
79 : /* ------------------------------------------------------------------- */
80 : int nRequestXSize;
81 55 : if ((nBlockXOff + 1) * nBlockXSize > nRasterXSize)
82 : {
83 0 : nRequestXSize = nRasterXSize - nBlockXOff * nBlockXSize;
84 0 : memset(pImage, 0,
85 0 : static_cast<size_t>(GDALGetDataTypeSizeBytes(eDataType)) *
86 0 : nBlockXSize * nBlockYSize);
87 : }
88 : else
89 : {
90 55 : nRequestXSize = nBlockXSize;
91 : }
92 55 : if (eDataType == GDT_CInt16 && poBandFile->GetRasterCount() == 2)
93 0 : return poBandFile->RasterIO(
94 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
95 : nRequestXSize, nRequestYSize, pImage, nRequestXSize, nRequestYSize,
96 0 : GDT_Int16, 2, nullptr, 4, nBlockXSize * 4, 2, nullptr);
97 :
98 : /* -------------------------------------------------------------------- */
99 : /* File has one sample marked as sample format void, a 32bits. */
100 : /* -------------------------------------------------------------------- */
101 55 : else if (eDataType == GDT_CInt16 && poBandFile->GetRasterCount() == 1)
102 : {
103 0 : CPLErr eErr = poBandFile->RasterIO(
104 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
105 : nRequestXSize, nRequestYSize, pImage, nRequestXSize, nRequestYSize,
106 0 : GDT_CInt16, 1, nullptr, 4, nBlockXSize * 4, 0, nullptr);
107 :
108 0 : return eErr;
109 : }
110 :
111 : /* -------------------------------------------------------------------- */
112 : /* The 16bit case is straight forward. The underlying file */
113 : /* looks like a 16bit unsigned data too. */
114 : /* -------------------------------------------------------------------- */
115 55 : else if (eDataType == GDT_UInt16)
116 55 : return poBandFile->RasterIO(
117 55 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
118 : nRequestXSize, nRequestYSize, pImage, nRequestXSize, nRequestYSize,
119 55 : GDT_UInt16, 1, nullptr, 2, nBlockXSize * 2, 0, nullptr);
120 :
121 0 : else if (eDataType == GDT_Byte)
122 0 : return poBandFile->RasterIO(
123 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
124 : nRequestXSize, nRequestYSize, pImage, nRequestXSize, nRequestYSize,
125 0 : GDT_Byte, 1, nullptr, 1, nBlockXSize, 0, nullptr);
126 :
127 0 : CPLAssert(false);
128 : return CE_Failure;
129 : }
130 :
131 : /************************************************************************/
132 : /* SAFESLCRasterBand */
133 : /************************************************************************/
134 :
135 3 : SAFESLCRasterBand::SAFESLCRasterBand(
136 : SAFEDataset *poDSIn, GDALDataType eDataTypeIn, const CPLString &osSwath,
137 : const CPLString &osPolarization,
138 3 : std::unique_ptr<GDALDataset> &&poBandFileIn, BandType eBandType)
139 3 : : poBandFile(std::move(poBandFileIn))
140 : {
141 3 : poDS = poDSIn;
142 3 : eDataType = eDataTypeIn;
143 3 : m_eInputDataType = eDataTypeIn;
144 3 : GDALRasterBand *poSrcBand = poBandFile->GetRasterBand(1);
145 3 : poSrcBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
146 3 : m_eBandType = eBandType;
147 :
148 3 : if (!osSwath.empty())
149 3 : SetMetadataItem("SWATH", osSwath.c_str());
150 :
151 3 : if (!osPolarization.empty())
152 3 : SetMetadataItem("POLARIZATION", osPolarization.c_str());
153 :
154 : // For intensity band
155 3 : if (m_eBandType == INTENSITY)
156 2 : eDataType = GDT_Float32;
157 : else
158 : // For complex bands
159 1 : eDataType = GDT_CInt16;
160 3 : }
161 :
162 : /************************************************************************/
163 : /* IReadBlock() */
164 : /************************************************************************/
165 :
166 0 : CPLErr SAFESLCRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
167 : void *pImage)
168 :
169 : {
170 : /* -------------------------------------------------------------------- */
171 : /* If the last strip is partial, we need to avoid */
172 : /* over-requesting. We also need to initialize the extra part */
173 : /* of the block to zero. */
174 : /* -------------------------------------------------------------------- */
175 : int nRequestYSize;
176 0 : if ((nBlockYOff + 1) * nBlockYSize > nRasterYSize)
177 : {
178 0 : nRequestYSize = nRasterYSize - nBlockYOff * nBlockYSize;
179 0 : memset(pImage, 0,
180 0 : static_cast<size_t>(GDALGetDataTypeSizeBytes(eDataType)) *
181 0 : nBlockXSize * nBlockYSize);
182 : }
183 : else
184 : {
185 0 : nRequestYSize = nBlockYSize;
186 : }
187 :
188 : /*-------------------------------------------------------------------- */
189 : /* If the input imagery is tiled, also need to avoid over- */
190 : /* requesting in the X-direction. */
191 : /* ------------------------------------------------------------------- */
192 : int nRequestXSize;
193 0 : if ((nBlockXOff + 1) * nBlockXSize > nRasterXSize)
194 : {
195 0 : nRequestXSize = nRasterXSize - nBlockXOff * nBlockXSize;
196 0 : memset(pImage, 0,
197 0 : static_cast<size_t>(GDALGetDataTypeSizeBytes(eDataType)) *
198 0 : nBlockXSize * nBlockYSize);
199 : }
200 : else
201 0 : nRequestXSize = nBlockXSize;
202 :
203 0 : if (m_eInputDataType == GDT_CInt16 && poBandFile->GetRasterCount() == 2)
204 : {
205 0 : return poBandFile->RasterIO(
206 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
207 : nRequestXSize, nRequestYSize, pImage, nRequestXSize, nRequestYSize,
208 0 : GDT_Int16, 2, nullptr, 4, nBlockXSize * 4, 2, nullptr);
209 : }
210 : // File has one sample marked as sample format void, a 32bits.
211 0 : else if (m_eInputDataType == GDT_CInt16 &&
212 0 : poBandFile->GetRasterCount() == 1)
213 : {
214 0 : if (m_eBandType == COMPLEX)
215 : {
216 0 : CPLErr eErr = poBandFile->RasterIO(
217 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
218 : nRequestXSize, nRequestYSize, pImage, nRequestXSize,
219 0 : nRequestYSize, GDT_CInt16, 1, nullptr, 4, nBlockXSize * 4, 0,
220 : nullptr);
221 0 : if (eErr != CE_None)
222 : {
223 0 : return eErr;
224 : }
225 : }
226 0 : else if (m_eBandType == INTENSITY)
227 : {
228 0 : GInt16 *pnImageTmp = static_cast<GInt16 *>(VSI_MALLOC3_VERBOSE(
229 : 2 * sizeof(int16_t), nBlockXSize, nBlockYSize));
230 0 : if (!pnImageTmp)
231 : {
232 0 : return CE_Failure;
233 : }
234 :
235 0 : CPLErr eErr = poBandFile->RasterIO(
236 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
237 : nRequestXSize, nRequestYSize, pnImageTmp, nRequestXSize,
238 0 : nRequestYSize, GDT_CInt16, 1, nullptr, 4, nBlockXSize * 4, 0,
239 : nullptr);
240 0 : if (eErr != CE_None)
241 : {
242 0 : CPLFree(pnImageTmp);
243 0 : return eErr;
244 : }
245 :
246 0 : float *pfBuffer = static_cast<float *>(pImage);
247 : #ifdef USE_OMP
248 : omp_set_num_threads(GetNumThreadsToUse());
249 : #pragma omp parallel
250 : #endif
251 0 : for (int i = 0; i < nBlockYSize; i++)
252 : {
253 : #ifdef USE_OMP
254 : #pragma omp for nowait
255 : #endif
256 0 : for (int j = 0; j < nBlockXSize; j++)
257 : {
258 0 : int nPixOff = (2 * (i * nBlockXSize)) + (j * 2);
259 0 : int nOutPixOff = (i * nBlockXSize) + j;
260 0 : pfBuffer[nOutPixOff] = static_cast<float>(
261 0 : static_cast<double>(pnImageTmp[nPixOff] *
262 0 : pnImageTmp[nPixOff]) +
263 0 : static_cast<double>(pnImageTmp[nPixOff + 1] *
264 0 : pnImageTmp[nPixOff + 1]));
265 : }
266 : }
267 0 : CPLFree(pnImageTmp);
268 : }
269 0 : return CE_None;
270 : }
271 :
272 0 : CPLAssert(false);
273 : return CE_Failure;
274 : }
275 :
276 : /************************************************************************/
277 : /* SAFECalibRasterBand */
278 : /************************************************************************/
279 :
280 0 : SAFECalibratedRasterBand::SAFECalibratedRasterBand(
281 : SAFEDataset *poDSIn, GDALDataType eDataTypeIn, const CPLString &osSwath,
282 : const CPLString &osPolarization,
283 : std::unique_ptr<GDALDataset> &&poBandDatasetIn,
284 0 : const char *pszCalibrationFilename, CalibrationType eCalibrationType)
285 0 : : poBandDataset(std::move(poBandDatasetIn))
286 : {
287 0 : poDS = poDSIn;
288 0 : GDALRasterBand *poSrcBand = poBandDataset->GetRasterBand(1);
289 0 : poSrcBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
290 0 : eDataType = eDataTypeIn;
291 :
292 0 : if (!osSwath.empty())
293 0 : SetMetadataItem("SWATH", osSwath.c_str());
294 :
295 0 : if (!osPolarization.empty())
296 0 : SetMetadataItem("POLARIZATION", osPolarization.c_str());
297 :
298 0 : m_osCalibrationFilename = pszCalibrationFilename;
299 0 : m_eInputDataType = eDataTypeIn;
300 0 : eDataType = GDT_Float32;
301 0 : m_eCalibrationType = eCalibrationType;
302 0 : }
303 :
304 : /************************************************************************/
305 : /* ReadLUT() */
306 : /************************************************************************/
307 : /* Read the provided LUT in to m_ndTable */
308 : /************************************************************************/
309 0 : bool SAFECalibratedRasterBand::ReadLUT()
310 : {
311 0 : const char *const papszCalibrationNodes[3] = {
312 : "=calibrationVector.sigmaNought", "=calibrationVector.betaNought",
313 : "=calibrationVector.gamma"};
314 0 : CPLString osCalibrationNodeName = papszCalibrationNodes[m_eCalibrationType];
315 0 : const char *pszEndSpace = " ";
316 0 : CPLString osStartTime, osEndTime;
317 0 : CPLXMLNode *poLUT = CPLParseXMLFile(m_osCalibrationFilename);
318 0 : if (!poLUT)
319 0 : return false;
320 :
321 0 : CPLString osParseLUT;
322 0 : CPLString osParsePixelLUT;
323 0 : CPLString osParseAzimuthLUT;
324 0 : CPLString osParseLineNoLUT;
325 0 : for (CPLXMLNode *psNode = poLUT; psNode != nullptr;)
326 : {
327 0 : if (psNode->psNext != nullptr)
328 0 : psNode = psNode->psNext;
329 :
330 0 : if (EQUAL(psNode->pszValue, "calibration"))
331 : {
332 0 : for (psNode = psNode->psChild; psNode != nullptr;)
333 : {
334 0 : if (EQUAL(psNode->pszValue, "adsHeader"))
335 : {
336 : osStartTime =
337 0 : CPLGetXMLValue(psNode, "=adsHeader.startTime", " ");
338 : osEndTime =
339 0 : CPLGetXMLValue(psNode, "=adsHeader.stopTime", " ");
340 : }
341 :
342 0 : if (psNode->psNext != nullptr)
343 0 : psNode = psNode->psNext;
344 :
345 0 : if (EQUAL(psNode->pszValue, "calibrationVectorList"))
346 : {
347 0 : for (psNode = psNode->psChild; psNode != nullptr;
348 0 : psNode = psNode->psNext)
349 : {
350 0 : if (EQUAL(psNode->pszValue, "calibrationVector"))
351 : {
352 : osParseAzimuthLUT += CPLGetXMLValue(
353 0 : psNode, "=calibrationVector.azimuthTime", " ");
354 0 : osParseAzimuthLUT += pszEndSpace;
355 : osParseLineNoLUT += CPLGetXMLValue(
356 0 : psNode, "=calibrationVector.line", " ");
357 0 : osParseLineNoLUT += pszEndSpace;
358 : osParsePixelLUT =
359 0 : CPLGetXMLValue(psNode, "pixel", " ");
360 0 : m_nNumPixels = static_cast<int>(CPLAtof(
361 : CPLGetXMLValue(psNode, "pixel.count", " ")));
362 : osParseLUT += CPLGetXMLValue(
363 0 : psNode, osCalibrationNodeName, " ");
364 0 : osParseLUT += pszEndSpace;
365 : }
366 : }
367 : }
368 : }
369 : }
370 : }
371 0 : CPLDestroyXMLNode(poLUT);
372 :
373 0 : osParsePixelLUT += pszEndSpace;
374 :
375 : CPLStringList oStartTimeList(
376 0 : CSLTokenizeString2(osStartTime, " ", CSLT_HONOURSTRINGS));
377 0 : if (!oStartTimeList.size())
378 0 : return false;
379 0 : m_oStartTimePoint = getTimePoint(oStartTimeList[0]);
380 : CPLStringList oEndTimeList(
381 0 : CSLTokenizeString2(osEndTime, " ", CSLT_HONOURSTRINGS));
382 0 : if (!oEndTimeList.size())
383 0 : return false;
384 0 : m_oStopTimePoint = getTimePoint(oEndTimeList[0]);
385 : m_oAzimuthList.Assign(
386 0 : CSLTokenizeString2(osParseAzimuthLUT, " ", CSLT_HONOURSTRINGS));
387 : CPLStringList oLUTList(
388 0 : CSLTokenizeString2(osParseLUT, " ", CSLT_HONOURSTRINGS));
389 : CPLStringList oPixelList(
390 0 : CSLTokenizeString2(osParsePixelLUT, " ", CSLT_HONOURSTRINGS));
391 : CPLStringList oLineNoList(
392 0 : CSLTokenizeString2(osParseLineNoLUT, " ", CSLT_HONOURSTRINGS));
393 :
394 0 : m_anPixelLUT.resize(m_nNumPixels);
395 0 : for (int i = 0; i < m_nNumPixels; i++)
396 0 : m_anPixelLUT[i] = static_cast<int>(CPLAtof(oPixelList[i]));
397 :
398 0 : int nTableSize = oLUTList.size();
399 0 : m_afTable.resize(nTableSize);
400 0 : for (int i = 0; i < nTableSize; i++)
401 0 : m_afTable[i] = static_cast<float>(CPLAtof(oLUTList[i]));
402 :
403 0 : int nLineListSize = oLineNoList.size();
404 0 : m_anLineLUT.resize(nLineListSize);
405 0 : for (int i = 0; i < nLineListSize; i++)
406 0 : m_anLineLUT[i] = static_cast<int>(CPLAtof(oLineNoList[i]));
407 :
408 0 : return true;
409 : }
410 :
411 : /************************************************************************/
412 : /* IReadBlock() */
413 : /************************************************************************/
414 :
415 0 : CPLErr SAFECalibratedRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
416 : void *pImage)
417 :
418 : {
419 : /* -------------------------------------------------------------------- */
420 : /* If the last strip is partial, we need to avoid */
421 : /* over-requesting. We also need to initialize the extra part */
422 : /* of the block to zero. */
423 : /* -------------------------------------------------------------------- */
424 : int nRequestYSize;
425 0 : if ((nBlockYOff + 1) * nBlockYSize > nRasterYSize)
426 : {
427 0 : nRequestYSize = nRasterYSize - nBlockYOff * nBlockYSize;
428 0 : memset(pImage, 0,
429 0 : static_cast<size_t>(GDALGetDataTypeSizeBytes(eDataType)) *
430 0 : nBlockXSize * nBlockYSize);
431 : }
432 : else
433 : {
434 0 : nRequestYSize = nBlockYSize;
435 : }
436 :
437 : // Check LUT values and fail before reading
438 0 : int nLineCalVecIdx = getCalibrationVectorIndex(nBlockYOff);
439 0 : const char *pszVec0Str = m_oAzimuthList[nLineCalVecIdx];
440 0 : const char *pszVec1Str = m_oAzimuthList[nLineCalVecIdx + 1];
441 0 : if (((m_eInputDataType == GDT_CInt16) || (m_eInputDataType == GDT_Int16)) &&
442 0 : (!pszVec0Str || !pszVec1Str))
443 0 : return CE_Failure;
444 :
445 : /*-------------------------------------------------------------------- */
446 : /* If the input imagery is tiled, also need to avoid over- */
447 : /* requesting in the X-direction. */
448 : /* ------------------------------------------------------------------- */
449 : int nRequestXSize;
450 0 : if ((nBlockXOff + 1) * nBlockXSize > nRasterXSize)
451 : {
452 0 : nRequestXSize = nRasterXSize - nBlockXOff * nBlockXSize;
453 0 : memset(pImage, 0,
454 0 : static_cast<size_t>(GDALGetDataTypeSizeBytes(eDataType)) *
455 0 : nBlockXSize * nBlockYSize);
456 : }
457 : else
458 : {
459 0 : nRequestXSize = nBlockXSize;
460 : }
461 :
462 : TimePoint azTime = getazTime(m_oStartTimePoint, m_oStopTimePoint,
463 0 : nRasterYSize, nBlockYOff);
464 0 : TimePoint oVec0Time = getTimePoint(pszVec0Str);
465 0 : TimePoint oVec1Time = getTimePoint(pszVec1Str);
466 : double dfMuY =
467 0 : getTimeDiff(oVec0Time, azTime) / getTimeDiff(oVec0Time, oVec1Time);
468 :
469 0 : if (m_eInputDataType == GDT_CInt16)
470 : {
471 0 : CPLErr eErr = CE_None;
472 : GInt16 *pnImageTmp = static_cast<GInt16 *>(
473 0 : VSI_MALLOC3_VERBOSE(2 * sizeof(int16_t), nBlockXSize, nBlockYSize));
474 0 : if (!pnImageTmp)
475 0 : return CE_Failure;
476 :
477 0 : if (poBandDataset->GetRasterCount() == 2)
478 : {
479 0 : eErr = poBandDataset->RasterIO(
480 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
481 : nRequestXSize, nRequestYSize, pnImageTmp, nRequestXSize,
482 0 : nRequestYSize, GDT_Int16, 2, nullptr, 4, nBlockXSize * 4, 2,
483 : nullptr);
484 : }
485 : /* --------------------------------------------------------------------
486 : */
487 : /* File has one sample marked as sample format void, a 32bits. */
488 : /* --------------------------------------------------------------------
489 : */
490 0 : else if (poBandDataset->GetRasterCount() == 1)
491 : {
492 0 : eErr = poBandDataset->RasterIO(
493 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
494 : nRequestXSize, nRequestYSize, pnImageTmp, nRequestXSize,
495 0 : nRequestYSize, GDT_CInt16, 1, nullptr, 4, nBlockXSize * 4, 0,
496 : nullptr);
497 : }
498 :
499 : // Interpolation of LUT Value
500 : #ifdef USE_OMP
501 : omp_set_num_threads(GetNumThreadsToUse());
502 : #pragma omp parallel
503 : #endif
504 0 : for (int i = 0; i < nBlockYSize; i++)
505 : {
506 : #ifdef USE_OMP
507 : #pragma omp for nowait
508 : #endif
509 0 : for (int j = 0; j < nBlockXSize; j++)
510 : {
511 0 : int nPixOff = (2 * (i * nBlockXSize)) + (j * 2);
512 0 : int nOutPixOff = (i * nBlockXSize) + j;
513 0 : int nPixelCalvecIdx = getPixelIndex(j);
514 0 : double dfMuX = (double)(j - m_anPixelLUT[nPixelCalvecIdx]) /
515 0 : (double)(m_anPixelLUT[nPixelCalvecIdx + 1] -
516 0 : m_anPixelLUT[nPixelCalvecIdx]);
517 0 : int lutIdx1 = (nLineCalVecIdx * m_nNumPixels) + nPixelCalvecIdx;
518 0 : int lutIdx2 =
519 0 : (nLineCalVecIdx * m_nNumPixels) + (nPixelCalvecIdx + 1);
520 0 : int lutIdx3 =
521 0 : ((nLineCalVecIdx + 1) * m_nNumPixels) + nPixelCalvecIdx;
522 0 : int lutIdx4 = ((nLineCalVecIdx + 1) * m_nNumPixels) +
523 0 : (nPixelCalvecIdx + 1);
524 : double dfLutValue =
525 0 : ((1 - dfMuY) * (((1 - dfMuX) * m_afTable[lutIdx1]) +
526 0 : (dfMuX * m_afTable[lutIdx2]))) +
527 0 : (dfMuY * (((1 - dfMuX) * m_afTable[lutIdx3]) +
528 0 : (dfMuX * m_afTable[lutIdx4])));
529 0 : double dfNum = static_cast<double>(
530 0 : (pnImageTmp[nPixOff] * pnImageTmp[nPixOff]) +
531 0 : (pnImageTmp[nPixOff + 1] * pnImageTmp[nPixOff + 1]));
532 0 : double dfCalibValue = dfNum / (dfLutValue * dfLutValue);
533 0 : ((float *)pImage)[nOutPixOff] = (float)dfCalibValue;
534 : }
535 : }
536 0 : CPLFree(pnImageTmp);
537 0 : return eErr;
538 : }
539 0 : else if (m_eInputDataType == GDT_UInt16)
540 : {
541 0 : CPLErr eErr = CE_None;
542 0 : GUInt16 *pnImageTmp = static_cast<GUInt16 *>(VSI_MALLOC3_VERBOSE(
543 : nBlockXSize, nBlockYSize, GDALGetDataTypeSizeBytes(GDT_UInt16)));
544 0 : if (!pnImageTmp)
545 0 : return CE_Failure;
546 0 : eErr = poBandDataset->RasterIO(GF_Read, nBlockXOff * nBlockXSize,
547 0 : nBlockYOff * nBlockYSize, nRequestXSize,
548 : nRequestYSize, pnImageTmp, nRequestXSize,
549 : nRequestYSize, GDT_UInt16, 1, nullptr, 2,
550 0 : nBlockXSize * 2, 0, nullptr);
551 :
552 : #ifdef USE_OMP
553 : omp_set_num_threads(GetNumThreadsToUse());
554 : #pragma omp parallel
555 : #endif
556 0 : for (int i = 0; i < nBlockYSize; i++)
557 : {
558 : #ifdef USE_OMP
559 : #pragma omp for nowait
560 : #endif
561 0 : for (int j = 0; j < nBlockXSize; j++)
562 : {
563 0 : int nPixOff = (i * nBlockXSize) + j;
564 0 : int nPixelCalvecIdx = getPixelIndex(j);
565 0 : double dfMuX = (double)(j - m_anPixelLUT[nPixelCalvecIdx]) /
566 0 : (double)(m_anPixelLUT[nPixelCalvecIdx + 1] -
567 0 : m_anPixelLUT[nPixelCalvecIdx]);
568 0 : int lutIdx1 = (nLineCalVecIdx * m_nNumPixels) + nPixelCalvecIdx;
569 0 : int lutIdx2 =
570 0 : (nLineCalVecIdx * m_nNumPixels) + (nPixelCalvecIdx + 1);
571 0 : int lutIdx3 =
572 0 : ((nLineCalVecIdx + 1) * m_nNumPixels) + (nPixelCalvecIdx);
573 0 : int lutIdx4 = ((nLineCalVecIdx + 1) * m_nNumPixels) +
574 0 : (nPixelCalvecIdx + 1);
575 : double dfLutValue =
576 0 : ((1 - dfMuY) * (((1 - dfMuX) * m_afTable[lutIdx1]) +
577 0 : (dfMuX * m_afTable[lutIdx2]))) +
578 0 : (dfMuY * (((1 - dfMuX) * m_afTable[lutIdx3]) +
579 0 : (dfMuX * m_afTable[lutIdx4])));
580 0 : double dfCalibValue =
581 0 : (double)(pnImageTmp[nPixOff] * pnImageTmp[nPixOff]) /
582 0 : (dfLutValue * dfLutValue);
583 0 : ((float *)pImage)[nPixOff] = (float)dfCalibValue;
584 : }
585 : }
586 0 : CPLFree(pnImageTmp);
587 0 : return eErr;
588 : }
589 0 : else if (eDataType == GDT_Byte) // Check if this is required.
590 0 : return poBandDataset->RasterIO(
591 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
592 : nRequestXSize, nRequestYSize, pImage, nRequestXSize, nRequestYSize,
593 0 : GDT_Byte, 1, nullptr, 1, nBlockXSize, 0, nullptr);
594 :
595 0 : CPLAssert(false);
596 : return CE_Failure;
597 : }
598 :
599 : /************************************************************************/
600 : /* ==================================================================== */
601 : /* SAFEDataset */
602 : /* ==================================================================== */
603 : /************************************************************************/
604 :
605 : /************************************************************************/
606 : /* SAFEDataset() */
607 : /************************************************************************/
608 :
609 11 : SAFEDataset::SAFEDataset()
610 : {
611 11 : m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
612 11 : }
613 :
614 : /************************************************************************/
615 : /* ~SAFEDataset() */
616 : /************************************************************************/
617 :
618 22 : SAFEDataset::~SAFEDataset()
619 :
620 : {
621 11 : SAFEDataset::FlushCache(true);
622 :
623 11 : if (nGCPCount > 0)
624 : {
625 11 : GDALDeinitGCPs(nGCPCount, pasGCPList);
626 11 : CPLFree(pasGCPList);
627 : }
628 :
629 11 : SAFEDataset::CloseDependentDatasets();
630 :
631 11 : CSLDestroy(papszSubDatasets);
632 11 : CSLDestroy(papszExtraFiles);
633 22 : }
634 :
635 : /************************************************************************/
636 : /* CloseDependentDatasets() */
637 : /************************************************************************/
638 :
639 11 : int SAFEDataset::CloseDependentDatasets()
640 : {
641 11 : int bHasDroppedRef = GDALPamDataset::CloseDependentDatasets();
642 :
643 11 : if (nBands != 0)
644 11 : bHasDroppedRef = TRUE;
645 :
646 28 : for (int iBand = 0; iBand < nBands; iBand++)
647 : {
648 17 : delete papoBands[iBand];
649 : }
650 11 : nBands = 0;
651 :
652 11 : return bHasDroppedRef;
653 : }
654 :
655 : /************************************************************************/
656 : /* GetMetaDataObject() */
657 : /************************************************************************/
658 :
659 : const CPLXMLNode *
660 121 : SAFEDataset::GetMetaDataObject(const CPLXMLNode *psMetaDataObjects,
661 : const char *metadataObjectId)
662 : {
663 : /* -------------------------------------------------------------------- */
664 : /* Look for DataObject Element by ID. */
665 : /* -------------------------------------------------------------------- */
666 1848 : for (const CPLXMLNode *psMDO = psMetaDataObjects->psChild; psMDO != nullptr;
667 1727 : psMDO = psMDO->psNext)
668 : {
669 1848 : if (psMDO->eType != CXT_Element ||
670 1848 : !(EQUAL(psMDO->pszValue, "metadataObject")))
671 : {
672 0 : continue;
673 : }
674 :
675 1848 : const char *pszElementID = CPLGetXMLValue(psMDO, "ID", "");
676 :
677 1848 : if (EQUAL(pszElementID, metadataObjectId))
678 : {
679 121 : return psMDO;
680 : }
681 : }
682 :
683 0 : CPLError(CE_Warning, CPLE_AppDefined, "MetadataObject not found with ID=%s",
684 : metadataObjectId);
685 :
686 0 : return nullptr;
687 : }
688 :
689 : /************************************************************************/
690 : /* GetDataObject() */
691 : /************************************************************************/
692 :
693 88 : const CPLXMLNode *SAFEDataset::GetDataObject(const CPLXMLNode *psDataObjects,
694 : const char *dataObjectId)
695 : {
696 : /* -------------------------------------------------------------------- */
697 : /* Look for DataObject Element by ID. */
698 : /* -------------------------------------------------------------------- */
699 836 : for (const CPLXMLNode *psDO = psDataObjects->psChild; psDO != nullptr;
700 748 : psDO = psDO->psNext)
701 : {
702 836 : if (psDO->eType != CXT_Element ||
703 836 : !(EQUAL(psDO->pszValue, "dataObject")))
704 : {
705 0 : continue;
706 : }
707 :
708 836 : const char *pszElementID = CPLGetXMLValue(psDO, "ID", "");
709 :
710 836 : if (EQUAL(pszElementID, dataObjectId))
711 : {
712 88 : return psDO;
713 : }
714 : }
715 :
716 0 : CPLError(CE_Warning, CPLE_AppDefined, "DataObject not found with ID=%s",
717 : dataObjectId);
718 :
719 0 : return nullptr;
720 : }
721 :
722 : const CPLXMLNode *
723 66 : SAFEDataset::GetDataObject(const CPLXMLNode *psMetaDataObjects,
724 : const CPLXMLNode *psDataObjects,
725 : const char *metadataObjectId)
726 : {
727 : /* -------------------------------------------------------------------- */
728 : /* Look for MetadataObject Element by ID. */
729 : /* -------------------------------------------------------------------- */
730 : const CPLXMLNode *psMDO =
731 66 : SAFEDataset::GetMetaDataObject(psMetaDataObjects, metadataObjectId);
732 :
733 66 : if (psMDO != nullptr)
734 : {
735 : const char *dataObjectId =
736 66 : CPLGetXMLValue(psMDO, "dataObjectPointer.dataObjectID", "");
737 66 : if (*dataObjectId != '\0')
738 : {
739 66 : return SAFEDataset::GetDataObject(psDataObjects, dataObjectId);
740 : }
741 : }
742 :
743 0 : CPLError(CE_Warning, CPLE_AppDefined, "DataObject not found with MetaID=%s",
744 : metadataObjectId);
745 :
746 0 : return nullptr;
747 : }
748 :
749 : /************************************************************************/
750 : /* GetFileList() */
751 : /************************************************************************/
752 :
753 5 : char **SAFEDataset::GetFileList()
754 :
755 : {
756 5 : char **papszFileList = GDALPamDataset::GetFileList();
757 :
758 5 : papszFileList = CSLInsertStrings(papszFileList, -1, papszExtraFiles);
759 :
760 5 : return papszFileList;
761 : }
762 :
763 : /************************************************************************/
764 : /* Identify() */
765 : /************************************************************************/
766 :
767 55600 : int SAFEDataset::Identify(GDALOpenInfo *poOpenInfo)
768 : {
769 : /* Check for the case where we're trying to read the calibrated data: */
770 55600 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL1_CALIB:"))
771 : {
772 12 : return TRUE;
773 : }
774 :
775 : /* Check for directory access when there is a manifest.safe file in the
776 : directory. */
777 55588 : if (poOpenInfo->bIsDirectory)
778 : {
779 : VSIStatBufL sStat;
780 417 : CPLString osMDFilename = CPLFormCIFilename(poOpenInfo->pszFilename,
781 834 : "manifest.safe", nullptr);
782 :
783 417 : if (VSIStatL(osMDFilename, &sStat) == 0 && VSI_ISREG(sStat.st_mode))
784 : {
785 6 : GDALOpenInfo oOpenInfo(osMDFilename, GA_ReadOnly, nullptr);
786 3 : return Identify(&oOpenInfo);
787 : }
788 414 : return FALSE;
789 : }
790 :
791 : /* otherwise, do our normal stuff */
792 55171 : if (!EQUAL(CPLGetFilename(poOpenInfo->pszFilename), "manifest.safe"))
793 55155 : return FALSE;
794 :
795 16 : if (poOpenInfo->nHeaderBytes < 100)
796 0 : return FALSE;
797 :
798 16 : const char *pszHeader =
799 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader);
800 16 : if (!strstr(pszHeader, "<xfdu:XFDU"))
801 0 : return FALSE;
802 :
803 : // This driver doesn't handle Sentinel-2 or RCM (RADARSAT Constellation Mission) data
804 16 : if (strstr(pszHeader, "sentinel-2") ||
805 13 : strstr(pszHeader, "rcm_prod_manifest.xsd"))
806 4 : return FALSE;
807 :
808 12 : return TRUE;
809 : }
810 :
811 : /************************************************************************/
812 : /* Open() */
813 : /************************************************************************/
814 :
815 12 : GDALDataset *SAFEDataset::Open(GDALOpenInfo *poOpenInfo)
816 :
817 : {
818 :
819 : // Is this a SENTINEL-1 manifest.safe definition?
820 12 : if (!SAFEDataset::Identify(poOpenInfo))
821 : {
822 0 : return nullptr;
823 : }
824 :
825 : /* -------------------------------------------------------------------- */
826 : /* Get subdataset information, if relevant */
827 : /* -------------------------------------------------------------------- */
828 24 : CPLString osMDFilename;
829 12 : bool bIsSubDS = false;
830 12 : bool bIsSLC = false;
831 : // Subdataset 1st level selection (ex: for swath selection)
832 24 : CPLString osSelectedSubDS1;
833 : // Subdataset 2nd level selection (ex: for polarization selection)
834 24 : CPLString osSelectedSubDS2;
835 24 : CPLString osSelectedSubDS3;
836 : // 0 for SIGMA , 1 for BETA, 2 for GAMMA and # for UNCALIB dataset
837 12 : SAFECalibratedRasterBand::CalibrationType eCalibrationType =
838 : SAFECalibratedRasterBand::SIGMA_NOUGHT;
839 12 : bool bCalibrated = false;
840 :
841 : // 0 for amplitude, 1 for complex (2 band : I , Q) and 2 for INTENSITY
842 : typedef enum
843 : {
844 : UNKNOWN = -1,
845 : AMPLITUDE,
846 : COMPLEX,
847 : INTENSITY
848 : } RequestDataType;
849 :
850 12 : RequestDataType eRequestType = UNKNOWN;
851 : // Calibration Information selection
852 24 : CPLString osSubdatasetName;
853 :
854 12 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL1_CALIB:"))
855 : {
856 6 : bIsSubDS = true;
857 6 : osMDFilename = poOpenInfo->pszFilename + strlen("SENTINEL1_CALIB:");
858 6 : const char *pszSelectionCalib = strchr(osMDFilename.c_str(), ':');
859 12 : if (pszSelectionCalib == nullptr ||
860 6 : pszSelectionCalib == osMDFilename.c_str())
861 : {
862 0 : CPLError(CE_Failure, CPLE_AppDefined,
863 : "Invalid syntax for SENTINEL1_CALIB:");
864 1 : return nullptr;
865 : }
866 :
867 6 : CPLString osCalibrationValue = osMDFilename;
868 6 : osCalibrationValue.resize(pszSelectionCalib - osMDFilename.c_str());
869 6 : osMDFilename = pszSelectionCalib + strlen(":");
870 6 : if (EQUAL(osCalibrationValue.c_str(), "UNCALIB"))
871 : {
872 3 : bCalibrated = false;
873 : }
874 3 : else if (EQUAL(osCalibrationValue.c_str(), "SIGMA0"))
875 : {
876 3 : bCalibrated = true;
877 3 : eCalibrationType = SAFECalibratedRasterBand::SIGMA_NOUGHT;
878 : }
879 0 : else if (EQUAL(osCalibrationValue.c_str(), "BETA0"))
880 : {
881 0 : bCalibrated = true;
882 0 : eCalibrationType = SAFECalibratedRasterBand::BETA_NOUGHT;
883 : }
884 0 : else if (EQUAL(osCalibrationValue.c_str(), "GAMMA"))
885 : {
886 0 : bCalibrated = true;
887 0 : eCalibrationType = SAFECalibratedRasterBand::GAMMA;
888 : }
889 : else
890 : {
891 0 : CPLError(CE_Failure, CPLE_AppDefined,
892 : "Invalid syntax for SENTINEL1_CALIB:");
893 0 : return nullptr;
894 : }
895 :
896 6 : const auto nSelectionUnitPos = osMDFilename.rfind(':');
897 6 : if (nSelectionUnitPos == std::string::npos || nSelectionUnitPos == 0)
898 : {
899 0 : CPLError(CE_Failure, CPLE_AppDefined,
900 : "Invalid syntax for SENTINEL1_CALIB:");
901 0 : return nullptr;
902 : }
903 :
904 6 : CPLString osUnitValue = osMDFilename.substr(nSelectionUnitPos + 1);
905 6 : osMDFilename.resize(nSelectionUnitPos);
906 6 : if (EQUAL(osUnitValue.c_str(), "AMPLITUDE"))
907 3 : eRequestType = AMPLITUDE;
908 3 : else if (EQUAL(osUnitValue.c_str(), "COMPLEX"))
909 0 : eRequestType = COMPLEX;
910 3 : else if (EQUAL(osUnitValue.c_str(), "INTENSITY"))
911 2 : eRequestType = INTENSITY;
912 : else
913 : {
914 1 : CPLError(CE_Failure, CPLE_AppDefined,
915 : "Invalid syntax for SENTINEL1_CALIB:");
916 1 : return nullptr;
917 : }
918 :
919 5 : const auto nSelection1Pos = osMDFilename.rfind(':');
920 5 : if (nSelection1Pos == std::string::npos || nSelection1Pos == 0)
921 : {
922 0 : CPLError(CE_Failure, CPLE_AppDefined,
923 : "Invalid syntax for SENTINEL1_CALIB:");
924 0 : return nullptr;
925 : }
926 5 : osSelectedSubDS1 = osMDFilename.substr(nSelection1Pos + 1);
927 5 : osMDFilename.resize(nSelection1Pos);
928 :
929 5 : const auto nSelection2Pos = osSelectedSubDS1.find('_');
930 5 : if (nSelection2Pos != std::string::npos && nSelection2Pos != 0)
931 : {
932 4 : osSelectedSubDS2 = osSelectedSubDS1.substr(nSelection2Pos + 1);
933 4 : osSelectedSubDS1.resize(nSelection2Pos);
934 4 : const auto nSelection3Pos = osSelectedSubDS2.find('_');
935 4 : if (nSelection3Pos != std::string::npos && nSelection3Pos != 0)
936 : {
937 2 : osSelectedSubDS3 = osSelectedSubDS2.substr(nSelection3Pos + 1);
938 2 : osSelectedSubDS2.resize(nSelection3Pos);
939 : }
940 : }
941 :
942 : // update directory check:
943 : VSIStatBufL sStat;
944 5 : if (VSIStatL(osMDFilename.c_str(), &sStat) == 0)
945 5 : poOpenInfo->bIsDirectory = VSI_ISDIR(sStat.st_mode);
946 5 : if (!bCalibrated)
947 3 : osSubdatasetName = "UNCALIB";
948 2 : else if (eCalibrationType == SAFECalibratedRasterBand::SIGMA_NOUGHT)
949 2 : osSubdatasetName = "SIGMA0";
950 0 : else if (eCalibrationType == SAFECalibratedRasterBand::BETA_NOUGHT)
951 0 : osSubdatasetName = "BETA0";
952 0 : else if (eCalibrationType == SAFECalibratedRasterBand::GAMMA)
953 0 : osSubdatasetName = "GAMMA";
954 5 : osSubdatasetName += ":";
955 5 : if (!osUnitValue.empty())
956 : {
957 5 : osSubdatasetName += osUnitValue;
958 5 : osSubdatasetName += ":";
959 : }
960 5 : if (!osSelectedSubDS1.empty())
961 : {
962 5 : osSubdatasetName += osSelectedSubDS1;
963 5 : osSubdatasetName += ":";
964 : }
965 5 : if (!osSelectedSubDS2.empty())
966 : {
967 4 : osSubdatasetName += osSelectedSubDS2;
968 4 : osSubdatasetName += ":";
969 : }
970 5 : if (!osSelectedSubDS3.empty())
971 : {
972 2 : osSubdatasetName += osSelectedSubDS3;
973 2 : osSubdatasetName += ":";
974 : }
975 5 : if (!osSubdatasetName.empty())
976 : {
977 5 : if (osSubdatasetName.back() == ':')
978 5 : osSubdatasetName.pop_back();
979 : }
980 : }
981 : else
982 : {
983 6 : osMDFilename = poOpenInfo->pszFilename;
984 : }
985 :
986 11 : if (poOpenInfo->bIsDirectory)
987 : {
988 : osMDFilename =
989 2 : CPLFormCIFilename(osMDFilename.c_str(), "manifest.safe", nullptr);
990 : }
991 :
992 : /* -------------------------------------------------------------------- */
993 : /* Ingest the manifest.safe file. */
994 : /* -------------------------------------------------------------------- */
995 :
996 22 : auto psManifest = CPLXMLTreeCloser(CPLParseXMLFile(osMDFilename));
997 11 : if (psManifest == nullptr)
998 0 : return nullptr;
999 :
1000 22 : CPLString osPath(CPLGetPath(osMDFilename));
1001 :
1002 : /* -------------------------------------------------------------------- */
1003 : /* Confirm the requested access is supported. */
1004 : /* -------------------------------------------------------------------- */
1005 11 : if (poOpenInfo->eAccess == GA_Update)
1006 : {
1007 0 : CPLError(CE_Failure, CPLE_NotSupported,
1008 : "The SAFE driver does not support update access to existing"
1009 : " datasets.\n");
1010 0 : return nullptr;
1011 : }
1012 :
1013 : /* -------------------------------------------------------------------- */
1014 : /* Get contentUnit parent element. */
1015 : /* -------------------------------------------------------------------- */
1016 11 : const CPLXMLNode *psContentUnits = CPLGetXMLNode(
1017 : psManifest.get(), "=xfdu:XFDU.informationPackageMap.xfdu:contentUnit");
1018 11 : if (psContentUnits == nullptr)
1019 : {
1020 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1021 : "Failed to find <xfdu:XFDU><informationPackageMap>"
1022 : "<xfdu:contentUnit> in manifest file.");
1023 0 : return nullptr;
1024 : }
1025 :
1026 : /* -------------------------------------------------------------------- */
1027 : /* Get Metadata Objects element. */
1028 : /* -------------------------------------------------------------------- */
1029 : const CPLXMLNode *psMetaDataObjects =
1030 11 : CPLGetXMLNode(psManifest.get(), "=xfdu:XFDU.metadataSection");
1031 11 : if (psMetaDataObjects == nullptr)
1032 : {
1033 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1034 : "Failed to find <xfdu:XFDU><metadataSection>"
1035 : "in manifest file.");
1036 0 : return nullptr;
1037 : }
1038 :
1039 : /* -------------------------------------------------------------------- */
1040 : /* Get Data Objects element. */
1041 : /* -------------------------------------------------------------------- */
1042 : const CPLXMLNode *psDataObjects =
1043 11 : CPLGetXMLNode(psManifest.get(), "=xfdu:XFDU.dataObjectSection");
1044 11 : if (psDataObjects == nullptr)
1045 : {
1046 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1047 : "Failed to find <xfdu:XFDU><dataObjectSection> in document.");
1048 0 : return nullptr;
1049 : }
1050 :
1051 : /* -------------------------------------------------------------------- */
1052 : /* Create the dataset. */
1053 : /* -------------------------------------------------------------------- */
1054 22 : auto poDS = std::make_unique<SAFEDataset>();
1055 :
1056 11 : poDS->psManifest = std::move(psManifest);
1057 :
1058 : /* -------------------------------------------------------------------- */
1059 : /* Look for "Measurement Data Unit" contentUnit elements. */
1060 : /* -------------------------------------------------------------------- */
1061 : // Map with all measures aggregated by swath
1062 22 : std::map<CPLString, std::set<CPLString>> oMapSwaths2Pols;
1063 22 : std::vector<CPLString> oImageNumberSwPol;
1064 11 : bool isWave = false;
1065 :
1066 11 : for (CPLXMLNode *psContentUnit = psContentUnits->psChild;
1067 374 : psContentUnit != nullptr; psContentUnit = psContentUnit->psNext)
1068 : {
1069 363 : if (psContentUnit->eType != CXT_Element ||
1070 316 : !(EQUAL(psContentUnit->pszValue, "xfdu:contentUnit")))
1071 : {
1072 47 : continue;
1073 : }
1074 :
1075 316 : const char *pszUnitType = CPLGetXMLValue(psContentUnit, "unitType", "");
1076 :
1077 316 : const char *pszAnnotation = nullptr;
1078 316 : const char *pszCalibration = nullptr;
1079 316 : const char *pszMeasurement = nullptr;
1080 :
1081 316 : if (EQUAL(pszUnitType, "Measurement Data Unit"))
1082 : {
1083 : /* Get dmdID and dataObjectID */
1084 30 : const char *pszDmdID = CPLGetXMLValue(psContentUnit, "dmdID", "");
1085 30 : const char *pszDataObjectID = CPLGetXMLValue(
1086 : psContentUnit, "dataObjectPointer.dataObjectID", "");
1087 30 : if (*pszDataObjectID == '\0' || *pszDmdID == '\0')
1088 13 : continue;
1089 :
1090 : const CPLXMLNode *psDataObject =
1091 22 : SAFEDataset::GetDataObject(psDataObjects, pszDataObjectID);
1092 :
1093 22 : const char *pszRepId = CPLGetXMLValue(psDataObject, "repID", "");
1094 22 : if (!EQUAL(pszRepId, "s1Level1MeasurementSchema"))
1095 0 : continue;
1096 :
1097 22 : pszMeasurement = CPLGetXMLValue(psDataObject,
1098 : "byteStream.fileLocation.href", "");
1099 22 : if (*pszMeasurement == '\0')
1100 0 : continue;
1101 :
1102 22 : char **papszTokens = CSLTokenizeString2(pszDmdID, " ",
1103 : CSLT_ALLOWEMPTYTOKENS |
1104 : CSLT_STRIPLEADSPACES |
1105 : CSLT_STRIPENDSPACES);
1106 :
1107 110 : for (int j = 0; j < CSLCount(papszTokens); j++)
1108 : {
1109 88 : const char *pszId = papszTokens[j];
1110 88 : if (*pszId == '\0')
1111 22 : continue;
1112 :
1113 : // Map the metadata ID to the object element
1114 66 : const CPLXMLNode *psDO = SAFEDataset::GetDataObject(
1115 : psMetaDataObjects, psDataObjects, pszId);
1116 66 : if (psDO == nullptr)
1117 0 : continue;
1118 :
1119 : // check object type
1120 66 : pszRepId = CPLGetXMLValue(psDO, "repID", "");
1121 66 : if (EQUAL(pszRepId, "s1Level1ProductSchema"))
1122 : {
1123 : /* Get annotation filename */
1124 22 : pszAnnotation = CPLGetXMLValue(
1125 : psDO, "byteStream.fileLocation.href", "");
1126 22 : if (*pszAnnotation == '\0')
1127 0 : continue;
1128 : }
1129 44 : else if (EQUAL(pszRepId, "s1Level1CalibrationSchema"))
1130 : {
1131 22 : pszCalibration = CPLGetXMLValue(
1132 : psDO, "byteStream.fileLocation.href", "");
1133 22 : if (*pszCalibration == '\0')
1134 0 : continue;
1135 : }
1136 : else
1137 : {
1138 22 : continue;
1139 : }
1140 : }
1141 :
1142 22 : CSLDestroy(papszTokens);
1143 :
1144 22 : if (pszAnnotation == nullptr || pszCalibration == nullptr)
1145 0 : continue;
1146 :
1147 : // open Annotation XML file
1148 : const CPLString osAnnotationFilePath =
1149 22 : CPLFormFilename(osPath, pszAnnotation, nullptr);
1150 : const CPLString osCalibrationFilePath =
1151 22 : CPLFormFilename(osPath, pszCalibration, nullptr);
1152 :
1153 : CPLXMLTreeCloser psAnnotation(
1154 22 : CPLParseXMLFile(osAnnotationFilePath));
1155 22 : if (psAnnotation.get() == nullptr)
1156 0 : continue;
1157 : CPLXMLTreeCloser psCalibration(
1158 22 : CPLParseXMLFile(osCalibrationFilePath));
1159 22 : if (psCalibration.get() == nullptr)
1160 0 : continue;
1161 :
1162 : /* --------------------------------------------------------------------
1163 : */
1164 : /* Get overall image information. */
1165 : /* --------------------------------------------------------------------
1166 : */
1167 : const CPLString osProductType = CPLGetXMLValue(
1168 22 : psAnnotation.get(), "=product.adsHeader.productType", "UNK");
1169 : const CPLString osMissionId = CPLGetXMLValue(
1170 22 : psAnnotation.get(), "=product.adsHeader.missionId", "UNK");
1171 : const CPLString osPolarization = CPLGetXMLValue(
1172 22 : psAnnotation.get(), "=product.adsHeader.polarisation", "UNK");
1173 : const CPLString osMode = CPLGetXMLValue(
1174 22 : psAnnotation.get(), "=product.adsHeader.mode", "UNK");
1175 : const CPLString osSwath = CPLGetXMLValue(
1176 22 : psAnnotation.get(), "=product.adsHeader.swath", "UNK");
1177 : const CPLString osImageNumber = CPLGetXMLValue(
1178 22 : psAnnotation.get(), "=product.adsHeader.imageNumber", "UNK");
1179 :
1180 22 : oMapSwaths2Pols[osSwath].insert(osPolarization);
1181 22 : oImageNumberSwPol.push_back(osImageNumber + " " + osSwath + " " +
1182 : osPolarization);
1183 22 : if (EQUAL(osMode.c_str(), "WV"))
1184 6 : isWave = true;
1185 :
1186 22 : if (EQUAL(osProductType.c_str(), "SLC"))
1187 6 : bIsSLC = true;
1188 :
1189 : // if the dataunit is amplitude or complex and there is calibration
1190 : // applied it's not possible as calibrated datasets are intensity.
1191 22 : if (eRequestType != INTENSITY && bCalibrated)
1192 0 : continue;
1193 :
1194 22 : if (osSelectedSubDS1.empty())
1195 : {
1196 : // If not subdataset was selected,
1197 : // open the first one we can find.
1198 6 : osSelectedSubDS1 = osSwath;
1199 : }
1200 :
1201 22 : if (osSelectedSubDS3.empty() && isWave)
1202 : {
1203 : // If the selected mode is Wave mode (different file structure)
1204 : // open the first vignette in the dataset.
1205 1 : osSelectedSubDS3 = osImageNumber;
1206 : }
1207 22 : if (!EQUAL(osSelectedSubDS1.c_str(), osSwath.c_str()))
1208 : {
1209 : // do not mix swath, otherwise it does not work for SLC products
1210 3 : continue;
1211 : }
1212 :
1213 25 : if (!osSelectedSubDS2.empty() &&
1214 6 : (osSelectedSubDS2.find(osPolarization) == std::string::npos))
1215 : {
1216 : // Add only selected polarizations.
1217 2 : continue;
1218 : }
1219 :
1220 20 : if (!osSelectedSubDS3.empty() &&
1221 3 : !EQUAL(osSelectedSubDS3.c_str(), osImageNumber.c_str()))
1222 : {
1223 : // Add only selected image number (for Wave)
1224 0 : continue;
1225 : }
1226 :
1227 : // Changed the position of this code segment till nullptr
1228 17 : poDS->nRasterXSize = atoi(CPLGetXMLValue(
1229 17 : psAnnotation.get(),
1230 : "=product.imageAnnotation.imageInformation.numberOfSamples",
1231 : "-1"));
1232 17 : poDS->nRasterYSize = atoi(CPLGetXMLValue(
1233 17 : psAnnotation.get(),
1234 : "=product.imageAnnotation.imageInformation.numberOfLines",
1235 : "-1"));
1236 17 : if (poDS->nRasterXSize <= 1 || poDS->nRasterYSize <= 1)
1237 : {
1238 0 : CPLError(
1239 : CE_Failure, CPLE_OpenFailed,
1240 : "Non-sane raster dimensions provided in manifest.safe. "
1241 : "If this is a valid SENTINEL-1 scene, please contact your "
1242 : "data provider for a corrected dataset.");
1243 0 : return nullptr;
1244 : }
1245 :
1246 17 : poDS->SetMetadataItem("PRODUCT_TYPE", osProductType.c_str());
1247 17 : poDS->SetMetadataItem("MISSION_ID", osMissionId.c_str());
1248 17 : poDS->SetMetadataItem("MODE", osMode.c_str());
1249 17 : poDS->SetMetadataItem("SWATH", osSwath.c_str());
1250 :
1251 : /* --------------------------------------------------------------------
1252 : */
1253 : /* Get dataType (so we can recognize complex data), and the */
1254 : /* bitsPerSample. */
1255 : /* --------------------------------------------------------------------
1256 : */
1257 :
1258 17 : const char *pszDataType = CPLGetXMLValue(
1259 17 : psAnnotation.get(),
1260 : "=product.imageAnnotation.imageInformation.outputPixels", "");
1261 :
1262 : GDALDataType eDataType;
1263 17 : if (EQUAL(pszDataType, "16 bit Signed Integer"))
1264 3 : eDataType = GDT_CInt16;
1265 14 : else if (EQUAL(pszDataType, "16 bit Unsigned Integer"))
1266 14 : eDataType = GDT_UInt16;
1267 : else
1268 : {
1269 0 : CPLError(CE_Failure, CPLE_AppDefined,
1270 : "dataType=%s: not a supported configuration.",
1271 : pszDataType);
1272 0 : return nullptr;
1273 : }
1274 :
1275 : /* Extract pixel spacing information */
1276 17 : const char *pszPixelSpacing = CPLGetXMLValue(
1277 17 : psAnnotation.get(),
1278 : "=product.imageAnnotation.imageInformation.rangePixelSpacing",
1279 : "UNK");
1280 17 : poDS->SetMetadataItem("PIXEL_SPACING", pszPixelSpacing);
1281 :
1282 17 : const char *pszLineSpacing = CPLGetXMLValue(
1283 17 : psAnnotation.get(),
1284 : "=product.imageAnnotation.imageInformation.azimuthPixelSpacing",
1285 : "UNK");
1286 17 : poDS->SetMetadataItem("LINE_SPACING", pszLineSpacing);
1287 :
1288 : /* --------------------------------------------------------------------
1289 : */
1290 : /* Form full filename (path of manifest.safe + measurement
1291 : * file). */
1292 : /* --------------------------------------------------------------------
1293 : */
1294 : const std::string osFullFilename(
1295 17 : CPLFormFilename(osPath, pszMeasurement, nullptr));
1296 :
1297 : /* --------------------------------------------------------------------
1298 : */
1299 : /* Try and open the file. */
1300 : /* --------------------------------------------------------------------
1301 : */
1302 17 : CPLTurnFailureIntoWarning(true);
1303 : auto poBandFile = std::unique_ptr<GDALDataset>(
1304 : GDALDataset::Open(osFullFilename.c_str(),
1305 17 : GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR));
1306 17 : CPLTurnFailureIntoWarning(false);
1307 17 : if (poBandFile == nullptr)
1308 : {
1309 : // NOP
1310 : }
1311 17 : else if (poBandFile->GetRasterCount() == 0)
1312 : {
1313 0 : poBandFile.reset();
1314 : }
1315 : else
1316 : {
1317 34 : poDS->papszExtraFiles =
1318 17 : CSLAddString(poDS->papszExtraFiles, osAnnotationFilePath);
1319 34 : poDS->papszExtraFiles =
1320 17 : CSLAddString(poDS->papszExtraFiles, osCalibrationFilePath);
1321 34 : poDS->papszExtraFiles =
1322 17 : CSLAddString(poDS->papszExtraFiles, osFullFilename.c_str());
1323 : /* --------------------------------------------------------------------
1324 : */
1325 : /* Collect Annotation Processing Information */
1326 : /* --------------------------------------------------------------------
1327 : */
1328 17 : CPLXMLNode *psProcessingInfo = CPLGetXMLNode(
1329 : psAnnotation.get(),
1330 : "=product.imageAnnotation.processingInformation");
1331 :
1332 17 : if (psProcessingInfo != nullptr)
1333 : {
1334 28 : OGRSpatialReference oLL, oPrj;
1335 :
1336 : const char *pszEllipsoidName =
1337 14 : CPLGetXMLValue(psProcessingInfo, "ellipsoidName", "");
1338 14 : const double minor_axis = CPLAtof(CPLGetXMLValue(
1339 : psProcessingInfo, "ellipsoidSemiMinorAxis", "0.0"));
1340 14 : const double major_axis = CPLAtof(CPLGetXMLValue(
1341 : psProcessingInfo, "ellipsoidSemiMajorAxis", "0.0"));
1342 :
1343 14 : if (EQUAL(pszEllipsoidName, "") || (minor_axis == 0.0) ||
1344 : (major_axis == 0.0))
1345 : {
1346 0 : CPLError(CE_Warning, CPLE_AppDefined,
1347 : "Warning- incomplete"
1348 : " ellipsoid information. Using wgs-84 "
1349 : "parameters.\n");
1350 0 : oLL.SetWellKnownGeogCS("WGS84");
1351 0 : oPrj.SetWellKnownGeogCS("WGS84");
1352 : }
1353 14 : else if (EQUAL(pszEllipsoidName, "WGS84"))
1354 : {
1355 14 : oLL.SetWellKnownGeogCS("WGS84");
1356 14 : oPrj.SetWellKnownGeogCS("WGS84");
1357 : }
1358 : else
1359 : {
1360 0 : const double inv_flattening =
1361 0 : major_axis / (major_axis - minor_axis);
1362 0 : oLL.SetGeogCS("", "", pszEllipsoidName, major_axis,
1363 : inv_flattening);
1364 0 : oPrj.SetGeogCS("", "", pszEllipsoidName, major_axis,
1365 : inv_flattening);
1366 : }
1367 :
1368 14 : oLL.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1369 14 : poDS->m_oGCPSRS = std::move(oLL);
1370 : }
1371 :
1372 : /* --------------------------------------------------------------------
1373 : */
1374 : /* Collect GCPs. */
1375 : /* --------------------------------------------------------------------
1376 : */
1377 17 : CPLXMLNode *psGeoGrid = CPLGetXMLNode(
1378 : psAnnotation.get(),
1379 : "=product.geolocationGrid.geolocationGridPointList");
1380 :
1381 17 : if (psGeoGrid != nullptr)
1382 : {
1383 17 : if (poDS->nGCPCount > 0)
1384 : {
1385 6 : GDALDeinitGCPs(poDS->nGCPCount, poDS->pasGCPList);
1386 6 : CPLFree(poDS->pasGCPList);
1387 : }
1388 :
1389 : /* count GCPs */
1390 17 : poDS->nGCPCount = 0;
1391 :
1392 17 : for (CPLXMLNode *psNode = psGeoGrid->psChild;
1393 52 : psNode != nullptr; psNode = psNode->psNext)
1394 : {
1395 35 : if (EQUAL(psNode->pszValue, "geolocationGridPoint"))
1396 18 : poDS->nGCPCount++;
1397 : }
1398 :
1399 17 : if (poDS->nGCPCount > 0)
1400 : {
1401 34 : poDS->pasGCPList = static_cast<GDAL_GCP *>(
1402 17 : CPLCalloc(sizeof(GDAL_GCP), poDS->nGCPCount));
1403 :
1404 17 : poDS->nGCPCount = 0;
1405 :
1406 17 : for (CPLXMLNode *psNode = psGeoGrid->psChild;
1407 52 : psNode != nullptr; psNode = psNode->psNext)
1408 : {
1409 : GDAL_GCP *psGCP =
1410 35 : poDS->pasGCPList + poDS->nGCPCount;
1411 :
1412 35 : if (!EQUAL(psNode->pszValue,
1413 : "geolocationGridPoint"))
1414 17 : continue;
1415 :
1416 18 : poDS->nGCPCount++;
1417 :
1418 : char szID[32];
1419 18 : snprintf(szID, sizeof(szID), "%d", poDS->nGCPCount);
1420 18 : psGCP->pszId = CPLStrdup(szID);
1421 18 : psGCP->pszInfo = CPLStrdup("");
1422 18 : psGCP->dfGCPPixel =
1423 18 : CPLAtof(CPLGetXMLValue(psNode, "pixel", "0"));
1424 18 : psGCP->dfGCPLine =
1425 18 : CPLAtof(CPLGetXMLValue(psNode, "line", "0"));
1426 18 : psGCP->dfGCPX = CPLAtof(
1427 : CPLGetXMLValue(psNode, "longitude", ""));
1428 18 : psGCP->dfGCPY =
1429 18 : CPLAtof(CPLGetXMLValue(psNode, "latitude", ""));
1430 18 : psGCP->dfGCPZ =
1431 18 : CPLAtof(CPLGetXMLValue(psNode, "height", ""));
1432 : }
1433 : }
1434 : }
1435 :
1436 : // Create bands
1437 17 : if (EQUAL(osProductType.c_str(), "SLC"))
1438 : {
1439 : // only add bands if no subDS or uncalibrated subdataset
1440 : // with complex data. (Calibrated will always be intensity
1441 : // only)
1442 3 : if (!bCalibrated &&
1443 0 : (eRequestType == UNKNOWN || eRequestType == COMPLEX))
1444 : {
1445 1 : poBandFile->MarkAsShared();
1446 : SAFESLCRasterBand *poBand0 = new SAFESLCRasterBand(
1447 1 : poDS.get(), eDataType, osSwath, osPolarization,
1448 1 : std::move(poBandFile), SAFESLCRasterBand::COMPLEX);
1449 1 : poDS->SetBand(poDS->GetRasterCount() + 1, poBand0);
1450 : }
1451 2 : else if (eRequestType == INTENSITY) // Intensity
1452 : {
1453 : SAFESLCRasterBand *poBand1 = new SAFESLCRasterBand(
1454 2 : poDS.get(), eDataType, osSwath, osPolarization,
1455 2 : std::move(poBandFile),
1456 2 : SAFESLCRasterBand::INTENSITY);
1457 2 : poDS->SetBand(poDS->GetRasterCount() + 1, poBand1);
1458 : }
1459 : }
1460 14 : else if (!bCalibrated &&
1461 4 : (eRequestType == UNKNOWN || eRequestType == AMPLITUDE))
1462 : {
1463 : SAFERasterBand *poBand = new SAFERasterBand(
1464 14 : poDS.get(), eDataType, osSwath, osPolarization,
1465 14 : std::move(poBandFile));
1466 14 : poDS->SetBand(poDS->GetRasterCount() + 1, poBand);
1467 : }
1468 0 : else if (bCalibrated &&
1469 0 : (eRequestType == UNKNOWN || eRequestType == COMPLEX))
1470 : {
1471 : auto poBand = std::make_unique<SAFECalibratedRasterBand>(
1472 0 : poDS.get(), eDataType, osSwath, osPolarization,
1473 0 : std::move(poBandFile), osCalibrationFilePath,
1474 0 : eCalibrationType);
1475 0 : if (!poBand->ReadLUT())
1476 : {
1477 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1478 : "Reading calibration LUT(s) failed: %s.",
1479 : osCalibrationFilePath.c_str());
1480 0 : return nullptr;
1481 : }
1482 0 : poDS->SetBand(poDS->GetRasterCount() + 1,
1483 0 : std::move(poBand));
1484 : }
1485 : }
1486 : }
1487 : }
1488 :
1489 : // loop through all Swath/pols to add subdatasets
1490 11 : if (!bIsSubDS)
1491 : {
1492 : const CPLString aosCalibrationValues[4] = {"SIGMA0", "BETA0", "GAMMA",
1493 36 : "UNCALIB"};
1494 : const CPLString aosDataUnitValues[3] = {"AMPLITUDE", "COMPLEX",
1495 30 : "INTENSITY"};
1496 6 : if (!isWave)
1497 : {
1498 10 : for (const auto &iterSwath : oMapSwaths2Pols)
1499 : {
1500 10 : CPLString osSubDS1 = iterSwath.first;
1501 10 : CPLString osSubDS2;
1502 :
1503 15 : for (const auto &pol : iterSwath.second)
1504 : {
1505 10 : if (!osSubDS2.empty())
1506 5 : osSubDS2 += "+";
1507 10 : osSubDS2 += pol;
1508 : // Create single band or multiband complex SubDataset
1509 10 : int i = 0;
1510 10 : if (bIsSLC)
1511 : {
1512 0 : for (i = 0; i < 3; i++)
1513 : {
1514 0 : CPLString osCalibTemp = aosCalibrationValues[i];
1515 0 : poDS->AddSubDataset(
1516 : CPLSPrintf("SENTINEL1_CALIB:%s:%s:%s_%s:%s",
1517 : osCalibTemp.c_str(),
1518 : osMDFilename.c_str(),
1519 : osSubDS1.c_str(), pol.c_str(),
1520 : aosDataUnitValues[2].c_str()),
1521 : CPLSPrintf("Single band with %s swath and %s "
1522 : "polarization and %s calibration",
1523 : osSubDS1.c_str(), pol.c_str(),
1524 : osCalibTemp.c_str()));
1525 : }
1526 :
1527 0 : CPLString osCalibTemp = aosCalibrationValues[i];
1528 0 : poDS->AddSubDataset(
1529 : CPLSPrintf("SENTINEL1_CALIB:%s:%s:%s_%s:%s",
1530 : osCalibTemp.c_str(),
1531 : osMDFilename.c_str(), osSubDS1.c_str(),
1532 : pol.c_str(),
1533 : aosDataUnitValues[1].c_str()),
1534 : CPLSPrintf("Single band with %s swath and %s "
1535 : "polarization and %s calibration",
1536 : osSubDS1.c_str(), pol.c_str(),
1537 : osCalibTemp.c_str()));
1538 0 : poDS->AddSubDataset(
1539 : CPLSPrintf("SENTINEL1_CALIB:%s:%s:%s_%s:%s",
1540 : osCalibTemp.c_str(),
1541 : osMDFilename.c_str(), osSubDS1.c_str(),
1542 : pol.c_str(),
1543 : aosDataUnitValues[2].c_str()),
1544 : CPLSPrintf("Single band with %s swath and %s "
1545 : "polarization and %s calibration",
1546 : osSubDS1.c_str(), pol.c_str(),
1547 : osCalibTemp.c_str()));
1548 : }
1549 : else
1550 : {
1551 10 : i = 3;
1552 10 : CPLString osCalibTemp = aosCalibrationValues[i];
1553 :
1554 10 : poDS->AddSubDataset(
1555 : CPLSPrintf("SENTINEL1_CALIB:%s:%s:%s_%s:%s",
1556 : osCalibTemp.c_str(),
1557 : osMDFilename.c_str(), osSubDS1.c_str(),
1558 : pol.c_str(),
1559 : aosDataUnitValues[0].c_str()),
1560 : CPLSPrintf("Single band with %s swath and %s "
1561 : "polarization and %s calibration",
1562 : osSubDS1.c_str(), pol.c_str(),
1563 : osCalibTemp.c_str()));
1564 : }
1565 : }
1566 :
1567 5 : if (iterSwath.second.size() > 1)
1568 : {
1569 : // Create single band subdataset with all polarizations
1570 5 : int i = 0;
1571 5 : if (bIsSLC)
1572 : {
1573 0 : for (i = 0; i < 3; i++)
1574 : {
1575 0 : CPLString osCalibTemp = aosCalibrationValues[i];
1576 :
1577 0 : poDS->AddSubDataset(
1578 : CPLSPrintf("SENTINEL1_CALIB:%s:%s:%s:%s",
1579 : osCalibTemp.c_str(),
1580 : osMDFilename.c_str(),
1581 : osSubDS1.c_str(),
1582 : aosDataUnitValues[2].c_str()),
1583 : CPLSPrintf(
1584 : "%s swath with all polarizations (%s) as "
1585 : "bands and %s calibration",
1586 : osSubDS1.c_str(), osSubDS2.c_str(),
1587 : osCalibTemp.c_str()));
1588 : }
1589 :
1590 0 : CPLString osCalibTemp = aosCalibrationValues[i];
1591 0 : poDS->AddSubDataset(
1592 : CPLSPrintf("SENTINEL1_CALIB:%s:%s:%s:%s",
1593 : osCalibTemp.c_str(),
1594 : osMDFilename.c_str(), osSubDS1.c_str(),
1595 : aosDataUnitValues[1].c_str()),
1596 : CPLSPrintf(
1597 : "%s swath with all polarizations (%s) as "
1598 : "bands and %s calibration",
1599 : osSubDS1.c_str(), osSubDS2.c_str(),
1600 : osCalibTemp.c_str()));
1601 0 : poDS->AddSubDataset(
1602 : CPLSPrintf("SENTINEL1_CALIB:%s:%s:%s:%s",
1603 : osCalibTemp.c_str(),
1604 : osMDFilename.c_str(), osSubDS1.c_str(),
1605 : aosDataUnitValues[2].c_str()),
1606 : CPLSPrintf(
1607 : "%s swath with all polarizations (%s) as "
1608 : "bands and %s calibration",
1609 : osSubDS1.c_str(), osSubDS2.c_str(),
1610 : osCalibTemp.c_str()));
1611 : }
1612 : else
1613 : {
1614 5 : i = 3;
1615 5 : CPLString osCalibTemp = aosCalibrationValues[i];
1616 5 : poDS->AddSubDataset(
1617 : CPLSPrintf("SENTINEL1_CALIB:%s:%s:%s:%s",
1618 : osCalibTemp.c_str(),
1619 : osMDFilename.c_str(), osSubDS1.c_str(),
1620 : aosDataUnitValues[0].c_str()),
1621 : CPLSPrintf(
1622 : "%s swath with all polarizations (%s) as "
1623 : "bands and %s calibration",
1624 : osSubDS1.c_str(), osSubDS2.c_str(),
1625 : osCalibTemp.c_str()));
1626 : }
1627 : }
1628 : }
1629 : }
1630 : else
1631 : {
1632 3 : for (const CPLString &osImgSwPol : oImageNumberSwPol)
1633 : {
1634 4 : CPLString osImgSwPolTmp = osImgSwPol;
1635 2 : const char *pszImage = strchr(osImgSwPolTmp.c_str(), ' ');
1636 4 : CPLString osImage, osSwath, osPolarization;
1637 2 : if (pszImage != nullptr)
1638 : {
1639 2 : osImage = osImgSwPolTmp;
1640 2 : osImage.resize(pszImage - osImgSwPolTmp.c_str());
1641 2 : osImgSwPolTmp = pszImage + strlen(" ");
1642 2 : const char *pszSwath = strchr(osImgSwPolTmp.c_str(), ' ');
1643 2 : if (pszSwath != nullptr)
1644 : {
1645 2 : osSwath = osImgSwPolTmp;
1646 2 : osSwath.resize(pszSwath - osImgSwPolTmp.c_str());
1647 2 : osPolarization = pszSwath + strlen(" ");
1648 2 : int i = 0;
1649 :
1650 2 : if (bIsSLC)
1651 : {
1652 8 : for (i = 0; i < 3; i++)
1653 : {
1654 6 : CPLString osCalibTemp = aosCalibrationValues[i];
1655 :
1656 6 : poDS->AddSubDataset(
1657 : CPLSPrintf(
1658 : "SENTINEL1_CALIB:%s:%s:%s_%s_%s:%s",
1659 : osCalibTemp.c_str(),
1660 : osMDFilename.c_str(), osSwath.c_str(),
1661 : osPolarization.c_str(), osImage.c_str(),
1662 : aosDataUnitValues[2].c_str()),
1663 : CPLSPrintf(
1664 : "Single band with %s swath and %s "
1665 : "polarization and %s calibration",
1666 : osSwath.c_str(), osPolarization.c_str(),
1667 : osCalibTemp.c_str()));
1668 : }
1669 :
1670 2 : CPLString osCalibTemp = aosCalibrationValues[i];
1671 :
1672 2 : poDS->AddSubDataset(
1673 : CPLSPrintf(
1674 : "SENTINEL1_CALIB:%s:%s:%s_%s_%s:%s",
1675 : osCalibTemp.c_str(), osMDFilename.c_str(),
1676 : osSwath.c_str(), osPolarization.c_str(),
1677 : osImage.c_str(),
1678 : aosDataUnitValues[1].c_str()),
1679 : CPLSPrintf("Single band with %s swath and %s "
1680 : "polarization and %s calibration",
1681 : osSwath.c_str(),
1682 : osPolarization.c_str(),
1683 : osCalibTemp.c_str()));
1684 :
1685 2 : poDS->AddSubDataset(
1686 : CPLSPrintf(
1687 : "SENTINEL1_CALIB:%s:%s:%s_%s_%s:%s",
1688 : osCalibTemp.c_str(), osMDFilename.c_str(),
1689 : osSwath.c_str(), osPolarization.c_str(),
1690 : osImage.c_str(),
1691 : aosDataUnitValues[2].c_str()),
1692 : CPLSPrintf("Single band with %s swath and %s "
1693 : "polarization and %s calibration",
1694 : osSwath.c_str(),
1695 : osPolarization.c_str(),
1696 : osCalibTemp.c_str()));
1697 : }
1698 : else
1699 : {
1700 0 : i = 3;
1701 0 : CPLString osCalibTemp = aosCalibrationValues[i];
1702 :
1703 0 : poDS->AddSubDataset(
1704 : CPLSPrintf(
1705 : "SENTINEL1_CALIB:%s:%s:%s_%s_%s:%s",
1706 : osCalibTemp.c_str(), osMDFilename.c_str(),
1707 : osSwath.c_str(), osPolarization.c_str(),
1708 : osImage.c_str(),
1709 : aosDataUnitValues[0].c_str()),
1710 : CPLSPrintf("Single band with %s swath and %s "
1711 : "polarization and %s calibration",
1712 : osSwath.c_str(),
1713 : osPolarization.c_str(),
1714 : osCalibTemp.c_str()));
1715 : }
1716 : }
1717 : }
1718 : }
1719 : }
1720 : }
1721 11 : if (poDS->GetRasterCount() == 0)
1722 : {
1723 0 : CPLError(CE_Failure, CPLE_OpenFailed, "Measurement bands not found.");
1724 0 : return nullptr;
1725 : }
1726 :
1727 : /* -------------------------------------------------------------------- */
1728 : /* Collect more metadata elements */
1729 : /* -------------------------------------------------------------------- */
1730 :
1731 : /* -------------------------------------------------------------------- */
1732 : /* Platform information */
1733 : /* -------------------------------------------------------------------- */
1734 : const CPLXMLNode *psPlatformAttrs =
1735 11 : SAFEDataset::GetMetaDataObject(psMetaDataObjects, "platform");
1736 :
1737 11 : if (psPlatformAttrs != nullptr)
1738 : {
1739 : const char *pszItem =
1740 11 : CPLGetXMLValue(psPlatformAttrs,
1741 : "metadataWrap.xmlData.safe:platform"
1742 : ".safe:familyName",
1743 : "");
1744 11 : poDS->SetMetadataItem("SATELLITE_IDENTIFIER", pszItem);
1745 :
1746 : pszItem =
1747 11 : CPLGetXMLValue(psPlatformAttrs,
1748 : "metadataWrap.xmlData.safe:platform"
1749 : ".safe:instrument.safe:familyName.abbreviation",
1750 : "");
1751 11 : poDS->SetMetadataItem("SENSOR_IDENTIFIER", pszItem);
1752 :
1753 11 : pszItem = CPLGetXMLValue(psPlatformAttrs,
1754 : "metadataWrap.xmlData.safe:platform"
1755 : ".safe:instrument.safe:extension"
1756 : ".s1sarl1:instrumentMode.s1sarl1:mode",
1757 : "UNK");
1758 11 : poDS->SetMetadataItem("BEAM_MODE", pszItem);
1759 :
1760 11 : pszItem = CPLGetXMLValue(psPlatformAttrs,
1761 : "metadataWrap.xmlData.safe:platform"
1762 : ".safe:instrument.safe:extension"
1763 : ".s1sarl1:instrumentMode.s1sarl1:swath",
1764 : "UNK");
1765 11 : poDS->SetMetadataItem("BEAM_SWATH", pszItem);
1766 : }
1767 :
1768 : /* -------------------------------------------------------------------- */
1769 : /* Acquisition Period information */
1770 : /* -------------------------------------------------------------------- */
1771 : const CPLXMLNode *psAcquisitionAttrs =
1772 11 : SAFEDataset::GetMetaDataObject(psMetaDataObjects, "acquisitionPeriod");
1773 :
1774 11 : if (psAcquisitionAttrs != nullptr)
1775 : {
1776 : const char *pszItem =
1777 11 : CPLGetXMLValue(psAcquisitionAttrs,
1778 : "metadataWrap.xmlData.safe:acquisitionPeriod"
1779 : ".safe:startTime",
1780 : "UNK");
1781 11 : poDS->SetMetadataItem("ACQUISITION_START_TIME", pszItem);
1782 11 : pszItem = CPLGetXMLValue(psAcquisitionAttrs,
1783 : "metadataWrap.xmlData.safe:acquisitionPeriod"
1784 : ".safe:stopTime",
1785 : "UNK");
1786 11 : poDS->SetMetadataItem("ACQUISITION_STOP_TIME", pszItem);
1787 : }
1788 :
1789 : /* -------------------------------------------------------------------- */
1790 : /* Processing information */
1791 : /* -------------------------------------------------------------------- */
1792 : const CPLXMLNode *psProcessingAttrs =
1793 11 : SAFEDataset::GetMetaDataObject(psMetaDataObjects, "processing");
1794 :
1795 11 : if (psProcessingAttrs != nullptr)
1796 : {
1797 11 : const char *pszItem = CPLGetXMLValue(
1798 : psProcessingAttrs,
1799 : "metadataWrap.xmlData.safe:processing.safe:facility.name", "UNK");
1800 11 : poDS->SetMetadataItem("FACILITY_IDENTIFIER", pszItem);
1801 : }
1802 :
1803 : /* -------------------------------------------------------------------- */
1804 : /* Measurement Orbit Reference information */
1805 : /* -------------------------------------------------------------------- */
1806 11 : const CPLXMLNode *psOrbitAttrs = SAFEDataset::GetMetaDataObject(
1807 : psMetaDataObjects, "measurementOrbitReference");
1808 :
1809 11 : if (psOrbitAttrs != nullptr)
1810 : {
1811 : const char *pszItem =
1812 11 : CPLGetXMLValue(psOrbitAttrs,
1813 : "metadataWrap.xmlData.safe:orbitReference"
1814 : ".safe:orbitNumber",
1815 : "UNK");
1816 11 : poDS->SetMetadataItem("ORBIT_NUMBER", pszItem);
1817 11 : pszItem = CPLGetXMLValue(psOrbitAttrs,
1818 : "metadataWrap.xmlData.safe:orbitReference"
1819 : ".safe:extension.s1:orbitProperties.s1:pass",
1820 : "UNK");
1821 11 : poDS->SetMetadataItem("ORBIT_DIRECTION", pszItem);
1822 : }
1823 :
1824 : /* -------------------------------------------------------------------- */
1825 : /* Footprint */
1826 : /* -------------------------------------------------------------------- */
1827 11 : const CPLXMLNode *psFrameSet = SAFEDataset::GetMetaDataObject(
1828 : psMetaDataObjects, "measurementFrameSet");
1829 :
1830 11 : if (psFrameSet)
1831 : {
1832 11 : const auto psFootPrint = CPLGetXMLNode(psFrameSet, "metadataWrap."
1833 : "xmlData."
1834 : "safe:frameSet."
1835 : "safe:frame."
1836 : "safe:footPrint");
1837 11 : if (psFootPrint)
1838 : {
1839 : const char *pszSRSName =
1840 11 : CPLGetXMLValue(psFootPrint, "srsName", nullptr);
1841 : const char *pszCoordinates =
1842 11 : CPLGetXMLValue(psFootPrint, "gml:coordinates", nullptr);
1843 11 : if (pszSRSName &&
1844 11 : EQUAL(pszSRSName,
1845 11 : "http://www.opengis.net/gml/srs/epsg.xml#4326") &&
1846 : pszCoordinates)
1847 : {
1848 : const CPLStringList aosValues(
1849 22 : CSLTokenizeString2(pszCoordinates, " ,", 0));
1850 11 : if (aosValues.size() == 8)
1851 : {
1852 22 : poDS->SetMetadataItem(
1853 : "FOOTPRINT",
1854 : CPLSPrintf("POLYGON((%s %s,%s %s,%s %s,%s %s, %s %s))",
1855 : aosValues[1], aosValues[0], aosValues[3],
1856 : aosValues[2], aosValues[5], aosValues[4],
1857 : aosValues[7], aosValues[6], aosValues[1],
1858 11 : aosValues[0]));
1859 : }
1860 : }
1861 : }
1862 : }
1863 :
1864 : /* -------------------------------------------------------------------- */
1865 : /* Initialize any PAM information. */
1866 : /* -------------------------------------------------------------------- */
1867 22 : const CPLString osDescription = osMDFilename;
1868 :
1869 : /* -------------------------------------------------------------------- */
1870 : /* Initialize any PAM information. */
1871 : /* -------------------------------------------------------------------- */
1872 11 : poDS->SetDescription(osDescription);
1873 :
1874 11 : poDS->SetPhysicalFilename(osMDFilename);
1875 11 : if (!osSubdatasetName.empty())
1876 : {
1877 5 : poDS->SetDescription(poOpenInfo->pszFilename);
1878 5 : poDS->SetSubdatasetName(osSubdatasetName);
1879 : }
1880 :
1881 11 : poDS->TryLoadXML();
1882 :
1883 : /* -------------------------------------------------------------------- */
1884 : /* Check for overviews. */
1885 : /* -------------------------------------------------------------------- */
1886 11 : poDS->oOvManager.Initialize(poDS.get(), ":::VIRTUAL:::");
1887 :
1888 11 : return poDS.release();
1889 : }
1890 :
1891 : /************************************************************************/
1892 : /* AddSubDataset() */
1893 : /************************************************************************/
1894 25 : void SAFEDataset::AddSubDataset(const CPLString &osName,
1895 : const CPLString &osDesc)
1896 : {
1897 25 : ++m_nSubDSNum;
1898 25 : papszSubDatasets = CSLAddNameValue(
1899 : papszSubDatasets, CPLSPrintf("SUBDATASET_%d_NAME", m_nSubDSNum),
1900 : osName.c_str());
1901 25 : papszSubDatasets = CSLAddNameValue(
1902 : papszSubDatasets, CPLSPrintf("SUBDATASET_%d_DESC", m_nSubDSNum),
1903 : osDesc.c_str());
1904 25 : }
1905 :
1906 : /************************************************************************/
1907 : /* GetGCPCount() */
1908 : /************************************************************************/
1909 :
1910 2 : int SAFEDataset::GetGCPCount()
1911 : {
1912 2 : return nGCPCount;
1913 : }
1914 :
1915 : /************************************************************************/
1916 : /* GetGCPSpatialRef() */
1917 : /************************************************************************/
1918 :
1919 0 : const OGRSpatialReference *SAFEDataset::GetGCPSpatialRef() const
1920 :
1921 : {
1922 0 : return m_oGCPSRS.IsEmpty() ? nullptr : &m_oGCPSRS;
1923 : }
1924 :
1925 : /************************************************************************/
1926 : /* GetGCPs() */
1927 : /************************************************************************/
1928 :
1929 2 : const GDAL_GCP *SAFEDataset::GetGCPs()
1930 :
1931 : {
1932 2 : return pasGCPList;
1933 : }
1934 :
1935 : /************************************************************************/
1936 : /* GetGeoTransform() */
1937 : /************************************************************************/
1938 :
1939 0 : CPLErr SAFEDataset::GetGeoTransform(double *padfTransform)
1940 : {
1941 0 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
1942 :
1943 0 : if (bHaveGeoTransform)
1944 0 : return CE_None;
1945 :
1946 0 : return CE_Failure;
1947 : }
1948 :
1949 : /************************************************************************/
1950 : /* GetMetadataDomainList() */
1951 : /************************************************************************/
1952 :
1953 0 : char **SAFEDataset::GetMetadataDomainList()
1954 : {
1955 0 : return BuildMetadataDomainList(GDALDataset::GetMetadataDomainList(), TRUE,
1956 0 : "SUBDATASETS", nullptr);
1957 : }
1958 :
1959 : /************************************************************************/
1960 : /* GetMetadata() */
1961 : /************************************************************************/
1962 :
1963 3 : char **SAFEDataset::GetMetadata(const char *pszDomain)
1964 : {
1965 3 : if (pszDomain != nullptr && STARTS_WITH_CI(pszDomain, "SUBDATASETS"))
1966 3 : return papszSubDatasets;
1967 :
1968 0 : return GDALDataset::GetMetadata(pszDomain);
1969 : }
1970 :
1971 : /************************************************************************/
1972 : /* GDALRegister_SAFE() */
1973 : /************************************************************************/
1974 :
1975 1595 : void GDALRegister_SAFE()
1976 : {
1977 1595 : if (GDALGetDriverByName("SAFE") != nullptr)
1978 302 : return;
1979 :
1980 1293 : GDALDriver *poDriver = new GDALDriver();
1981 :
1982 1293 : poDriver->SetDescription("SAFE");
1983 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1984 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1985 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Sentinel-1 SAR SAFE Product");
1986 1293 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/safe.html");
1987 :
1988 1293 : poDriver->pfnOpen = SAFEDataset::Open;
1989 1293 : poDriver->pfnIdentify = SAFEDataset::Identify;
1990 :
1991 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
1992 : }
1993 :
1994 : /************************************************************************/
1995 : /* Azimuth time handling functions */
1996 : /************************************************************************/
1997 :
1998 : SAFECalibratedRasterBand::TimePoint
1999 0 : SAFECalibratedRasterBand::getTimePoint(const char *pszTime)
2000 : {
2001 : int nYear, nMonth, nDay, nHour, nMinute, nSeconds;
2002 : long nMicroSeconds;
2003 0 : sscanf(pszTime, "%d-%d-%dT%d:%d:%d.%ld", &nYear, &nMonth, &nDay, &nHour,
2004 : &nMinute, &nSeconds, &nMicroSeconds);
2005 :
2006 : struct tm oTm;
2007 0 : oTm.tm_sec = nSeconds;
2008 0 : oTm.tm_min = nMinute;
2009 0 : oTm.tm_hour = nHour;
2010 0 : oTm.tm_mday = nDay;
2011 0 : oTm.tm_mon = nMonth - 1;
2012 0 : oTm.tm_year = nYear - 1900;
2013 0 : oTm.tm_isdst = -1;
2014 :
2015 0 : std::time_t oTt = static_cast<std::time_t>(CPLYMDHMSToUnixTime(&oTm));
2016 :
2017 0 : TimePoint oTp = std::chrono::system_clock::from_time_t(oTt);
2018 0 : TimePoint oTp1 = oTp + std::chrono::microseconds(nMicroSeconds);
2019 0 : return oTp1;
2020 : }
2021 :
2022 0 : double SAFECalibratedRasterBand::getTimeDiff(TimePoint oT1, TimePoint oT2)
2023 : {
2024 0 : std::chrono::duration<double> oResult = (oT2 - oT1);
2025 0 : return oResult.count();
2026 : }
2027 :
2028 : SAFECalibratedRasterBand::TimePoint
2029 0 : SAFECalibratedRasterBand::getazTime(TimePoint oStart, TimePoint oStop,
2030 : long nNumOfLines, int nOffset)
2031 : {
2032 0 : double dfTemp = getTimeDiff(oStart, oStop);
2033 0 : dfTemp /= (nNumOfLines - 1);
2034 0 : unsigned long nTimeDiffMicro = static_cast<unsigned long>(dfTemp * 1000000);
2035 : TimePoint oResult =
2036 0 : oStart + (nOffset * std::chrono::microseconds(nTimeDiffMicro));
2037 0 : return oResult;
2038 : }
2039 :
2040 : /************************************************************************/
2041 : /* Utility functions used in interpolation */
2042 : /************************************************************************/
2043 :
2044 0 : int SAFECalibratedRasterBand::getCalibrationVectorIndex(int nLineNo)
2045 : {
2046 0 : for (size_t i = 1; i < m_anLineLUT.size(); i++)
2047 : {
2048 0 : if (nLineNo < m_anLineLUT[i])
2049 0 : return static_cast<int>(i - 1);
2050 : }
2051 0 : return 0;
2052 : }
2053 :
2054 0 : int SAFECalibratedRasterBand::getPixelIndex(int nPixelNo)
2055 : {
2056 0 : for (int i = 1; i < m_nNumPixels; i++)
2057 : {
2058 0 : if (nPixelNo < m_anPixelLUT[i])
2059 0 : return i - 1;
2060 : }
2061 0 : return 0;
2062 : }
|
