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 61919 : int SAFEDataset::Identify(GDALOpenInfo *poOpenInfo)
768 : {
769 : /* Check for the case where we're trying to read the calibrated data: */
770 61919 : 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 61907 : if (poOpenInfo->bIsDirectory)
778 : {
779 : VSIStatBufL sStat;
780 0 : CPLString osMDFilename = CPLFormCIFilenameSafe(
781 1232 : poOpenInfo->pszFilename, "manifest.safe", nullptr);
782 :
783 616 : 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 613 : return FALSE;
789 : }
790 :
791 : /* otherwise, do our normal stuff */
792 61291 : if (!EQUAL(CPLGetFilename(poOpenInfo->pszFilename), "manifest.safe"))
793 61278 : return FALSE;
794 :
795 1 : if (poOpenInfo->nHeaderBytes < 100)
796 0 : return FALSE;
797 :
798 1 : const char *pszHeader =
799 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader);
800 1 : 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 1 : if (strstr(pszHeader, "sentinel-2") ||
805 13 : strstr(pszHeader, "rcm_prod_manifest.xsd"))
806 0 : 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 4 : osMDFilename = CPLFormCIFilenameSafe(osMDFilename.c_str(),
989 2 : "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(CPLGetPathSafe(osMDFilename));
1001 :
1002 : /* -------------------------------------------------------------------- */
1003 : /* Confirm the requested access is supported. */
1004 : /* -------------------------------------------------------------------- */
1005 11 : if (poOpenInfo->eAccess == GA_Update)
1006 : {
1007 0 : ReportUpdateNotSupportedByDriver("SAFE");
1008 0 : return nullptr;
1009 : }
1010 :
1011 : /* -------------------------------------------------------------------- */
1012 : /* Get contentUnit parent element. */
1013 : /* -------------------------------------------------------------------- */
1014 11 : const CPLXMLNode *psContentUnits = CPLGetXMLNode(
1015 : psManifest.get(), "=xfdu:XFDU.informationPackageMap.xfdu:contentUnit");
1016 11 : if (psContentUnits == nullptr)
1017 : {
1018 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1019 : "Failed to find <xfdu:XFDU><informationPackageMap>"
1020 : "<xfdu:contentUnit> in manifest file.");
1021 0 : return nullptr;
1022 : }
1023 :
1024 : /* -------------------------------------------------------------------- */
1025 : /* Get Metadata Objects element. */
1026 : /* -------------------------------------------------------------------- */
1027 : const CPLXMLNode *psMetaDataObjects =
1028 11 : CPLGetXMLNode(psManifest.get(), "=xfdu:XFDU.metadataSection");
1029 11 : if (psMetaDataObjects == nullptr)
1030 : {
1031 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1032 : "Failed to find <xfdu:XFDU><metadataSection>"
1033 : "in manifest file.");
1034 0 : return nullptr;
1035 : }
1036 :
1037 : /* -------------------------------------------------------------------- */
1038 : /* Get Data Objects element. */
1039 : /* -------------------------------------------------------------------- */
1040 : const CPLXMLNode *psDataObjects =
1041 11 : CPLGetXMLNode(psManifest.get(), "=xfdu:XFDU.dataObjectSection");
1042 11 : if (psDataObjects == nullptr)
1043 : {
1044 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1045 : "Failed to find <xfdu:XFDU><dataObjectSection> in document.");
1046 0 : return nullptr;
1047 : }
1048 :
1049 : /* -------------------------------------------------------------------- */
1050 : /* Create the dataset. */
1051 : /* -------------------------------------------------------------------- */
1052 22 : auto poDS = std::make_unique<SAFEDataset>();
1053 :
1054 11 : poDS->psManifest = std::move(psManifest);
1055 :
1056 : /* -------------------------------------------------------------------- */
1057 : /* Look for "Measurement Data Unit" contentUnit elements. */
1058 : /* -------------------------------------------------------------------- */
1059 : // Map with all measures aggregated by swath
1060 22 : std::map<CPLString, std::set<CPLString>> oMapSwaths2Pols;
1061 22 : std::vector<CPLString> oImageNumberSwPol;
1062 11 : bool isWave = false;
1063 :
1064 11 : for (CPLXMLNode *psContentUnit = psContentUnits->psChild;
1065 374 : psContentUnit != nullptr; psContentUnit = psContentUnit->psNext)
1066 : {
1067 363 : if (psContentUnit->eType != CXT_Element ||
1068 316 : !(EQUAL(psContentUnit->pszValue, "xfdu:contentUnit")))
1069 : {
1070 47 : continue;
1071 : }
1072 :
1073 316 : const char *pszUnitType = CPLGetXMLValue(psContentUnit, "unitType", "");
1074 :
1075 316 : const char *pszAnnotation = nullptr;
1076 316 : const char *pszCalibration = nullptr;
1077 316 : const char *pszMeasurement = nullptr;
1078 :
1079 316 : if (EQUAL(pszUnitType, "Measurement Data Unit"))
1080 : {
1081 : /* Get dmdID and dataObjectID */
1082 30 : const char *pszDmdID = CPLGetXMLValue(psContentUnit, "dmdID", "");
1083 30 : const char *pszDataObjectID = CPLGetXMLValue(
1084 : psContentUnit, "dataObjectPointer.dataObjectID", "");
1085 30 : if (*pszDataObjectID == '\0' || *pszDmdID == '\0')
1086 13 : continue;
1087 :
1088 : const CPLXMLNode *psDataObject =
1089 22 : SAFEDataset::GetDataObject(psDataObjects, pszDataObjectID);
1090 :
1091 22 : const char *pszRepId = CPLGetXMLValue(psDataObject, "repID", "");
1092 22 : if (!EQUAL(pszRepId, "s1Level1MeasurementSchema"))
1093 0 : continue;
1094 :
1095 22 : pszMeasurement = CPLGetXMLValue(psDataObject,
1096 : "byteStream.fileLocation.href", "");
1097 22 : if (*pszMeasurement == '\0')
1098 0 : continue;
1099 :
1100 22 : if (CPLHasPathTraversal(pszMeasurement))
1101 : {
1102 0 : CPLError(CE_Failure, CPLE_AppDefined,
1103 : "Path traversal detected in %s", pszMeasurement);
1104 0 : return nullptr;
1105 : }
1106 :
1107 : const CPLStringList aosTokens(CSLTokenizeString2(
1108 : pszDmdID, " ",
1109 : CSLT_ALLOWEMPTYTOKENS | CSLT_STRIPLEADSPACES |
1110 22 : CSLT_STRIPENDSPACES));
1111 :
1112 110 : for (const char *pszId : aosTokens)
1113 : {
1114 88 : if (*pszId == '\0')
1115 22 : continue;
1116 :
1117 : // Map the metadata ID to the object element
1118 66 : const CPLXMLNode *psDO = SAFEDataset::GetDataObject(
1119 : psMetaDataObjects, psDataObjects, pszId);
1120 66 : if (psDO == nullptr)
1121 0 : continue;
1122 :
1123 : // check object type
1124 66 : pszRepId = CPLGetXMLValue(psDO, "repID", "");
1125 66 : if (EQUAL(pszRepId, "s1Level1ProductSchema"))
1126 : {
1127 : /* Get annotation filename */
1128 22 : pszAnnotation = CPLGetXMLValue(
1129 : psDO, "byteStream.fileLocation.href", "");
1130 22 : if (*pszAnnotation == '\0')
1131 0 : continue;
1132 : }
1133 44 : else if (EQUAL(pszRepId, "s1Level1CalibrationSchema"))
1134 : {
1135 22 : pszCalibration = CPLGetXMLValue(
1136 : psDO, "byteStream.fileLocation.href", "");
1137 22 : if (*pszCalibration == '\0')
1138 0 : continue;
1139 : }
1140 : else
1141 : {
1142 22 : continue;
1143 : }
1144 : }
1145 :
1146 22 : if (pszAnnotation == nullptr || pszCalibration == nullptr)
1147 0 : continue;
1148 :
1149 22 : if (CPLHasPathTraversal(pszAnnotation))
1150 : {
1151 0 : CPLError(CE_Failure, CPLE_AppDefined,
1152 : "Path traversal detected in %s", pszAnnotation);
1153 0 : return nullptr;
1154 : }
1155 :
1156 22 : if (CPLHasPathTraversal(pszCalibration))
1157 : {
1158 0 : CPLError(CE_Failure, CPLE_AppDefined,
1159 : "Path traversal detected in %s", pszCalibration);
1160 0 : return nullptr;
1161 : }
1162 :
1163 : // open Annotation XML file
1164 : const CPLString osAnnotationFilePath =
1165 22 : CPLFormFilenameSafe(osPath, pszAnnotation, nullptr);
1166 : const CPLString osCalibrationFilePath =
1167 22 : CPLFormFilenameSafe(osPath, pszCalibration, nullptr);
1168 :
1169 : CPLXMLTreeCloser psAnnotation(
1170 22 : CPLParseXMLFile(osAnnotationFilePath));
1171 22 : if (psAnnotation.get() == nullptr)
1172 0 : continue;
1173 : CPLXMLTreeCloser psCalibration(
1174 22 : CPLParseXMLFile(osCalibrationFilePath));
1175 22 : if (psCalibration.get() == nullptr)
1176 0 : continue;
1177 :
1178 : /* --------------------------------------------------------------------
1179 : */
1180 : /* Get overall image information. */
1181 : /* --------------------------------------------------------------------
1182 : */
1183 : const CPLString osProductType = CPLGetXMLValue(
1184 22 : psAnnotation.get(), "=product.adsHeader.productType", "UNK");
1185 : const CPLString osMissionId = CPLGetXMLValue(
1186 22 : psAnnotation.get(), "=product.adsHeader.missionId", "UNK");
1187 : const CPLString osPolarization = CPLGetXMLValue(
1188 22 : psAnnotation.get(), "=product.adsHeader.polarisation", "UNK");
1189 : const CPLString osMode = CPLGetXMLValue(
1190 22 : psAnnotation.get(), "=product.adsHeader.mode", "UNK");
1191 : const CPLString osSwath = CPLGetXMLValue(
1192 22 : psAnnotation.get(), "=product.adsHeader.swath", "UNK");
1193 : const CPLString osImageNumber = CPLGetXMLValue(
1194 22 : psAnnotation.get(), "=product.adsHeader.imageNumber", "UNK");
1195 :
1196 22 : oMapSwaths2Pols[osSwath].insert(osPolarization);
1197 22 : oImageNumberSwPol.push_back(osImageNumber + " " + osSwath + " " +
1198 : osPolarization);
1199 22 : if (EQUAL(osMode.c_str(), "WV"))
1200 6 : isWave = true;
1201 :
1202 22 : if (EQUAL(osProductType.c_str(), "SLC"))
1203 6 : bIsSLC = true;
1204 :
1205 : // if the dataunit is amplitude or complex and there is calibration
1206 : // applied it's not possible as calibrated datasets are intensity.
1207 22 : if (eRequestType != INTENSITY && bCalibrated)
1208 0 : continue;
1209 :
1210 22 : if (osSelectedSubDS1.empty())
1211 : {
1212 : // If not subdataset was selected,
1213 : // open the first one we can find.
1214 6 : osSelectedSubDS1 = osSwath;
1215 : }
1216 :
1217 22 : if (osSelectedSubDS3.empty() && isWave)
1218 : {
1219 : // If the selected mode is Wave mode (different file structure)
1220 : // open the first vignette in the dataset.
1221 1 : osSelectedSubDS3 = osImageNumber;
1222 : }
1223 22 : if (!EQUAL(osSelectedSubDS1.c_str(), osSwath.c_str()))
1224 : {
1225 : // do not mix swath, otherwise it does not work for SLC products
1226 3 : continue;
1227 : }
1228 :
1229 25 : if (!osSelectedSubDS2.empty() &&
1230 6 : (osSelectedSubDS2.find(osPolarization) == std::string::npos))
1231 : {
1232 : // Add only selected polarizations.
1233 2 : continue;
1234 : }
1235 :
1236 20 : if (!osSelectedSubDS3.empty() &&
1237 3 : !EQUAL(osSelectedSubDS3.c_str(), osImageNumber.c_str()))
1238 : {
1239 : // Add only selected image number (for Wave)
1240 0 : continue;
1241 : }
1242 :
1243 : // Changed the position of this code segment till nullptr
1244 17 : poDS->nRasterXSize = atoi(CPLGetXMLValue(
1245 17 : psAnnotation.get(),
1246 : "=product.imageAnnotation.imageInformation.numberOfSamples",
1247 : "-1"));
1248 17 : poDS->nRasterYSize = atoi(CPLGetXMLValue(
1249 17 : psAnnotation.get(),
1250 : "=product.imageAnnotation.imageInformation.numberOfLines",
1251 : "-1"));
1252 17 : if (poDS->nRasterXSize <= 1 || poDS->nRasterYSize <= 1)
1253 : {
1254 0 : CPLError(
1255 : CE_Failure, CPLE_OpenFailed,
1256 : "Non-sane raster dimensions provided in manifest.safe. "
1257 : "If this is a valid SENTINEL-1 scene, please contact your "
1258 : "data provider for a corrected dataset.");
1259 0 : return nullptr;
1260 : }
1261 :
1262 17 : poDS->SetMetadataItem("PRODUCT_TYPE", osProductType.c_str());
1263 17 : poDS->SetMetadataItem("MISSION_ID", osMissionId.c_str());
1264 17 : poDS->SetMetadataItem("MODE", osMode.c_str());
1265 17 : poDS->SetMetadataItem("SWATH", osSwath.c_str());
1266 :
1267 : /* --------------------------------------------------------------------
1268 : */
1269 : /* Get dataType (so we can recognize complex data), and the */
1270 : /* bitsPerSample. */
1271 : /* --------------------------------------------------------------------
1272 : */
1273 :
1274 17 : const char *pszDataType = CPLGetXMLValue(
1275 17 : psAnnotation.get(),
1276 : "=product.imageAnnotation.imageInformation.outputPixels", "");
1277 :
1278 : GDALDataType eDataType;
1279 17 : if (EQUAL(pszDataType, "16 bit Signed Integer"))
1280 3 : eDataType = GDT_CInt16;
1281 14 : else if (EQUAL(pszDataType, "16 bit Unsigned Integer"))
1282 14 : eDataType = GDT_UInt16;
1283 : else
1284 : {
1285 0 : CPLError(CE_Failure, CPLE_AppDefined,
1286 : "dataType=%s: not a supported configuration.",
1287 : pszDataType);
1288 0 : return nullptr;
1289 : }
1290 :
1291 : /* Extract pixel spacing information */
1292 17 : const char *pszPixelSpacing = CPLGetXMLValue(
1293 17 : psAnnotation.get(),
1294 : "=product.imageAnnotation.imageInformation.rangePixelSpacing",
1295 : "UNK");
1296 17 : poDS->SetMetadataItem("PIXEL_SPACING", pszPixelSpacing);
1297 :
1298 17 : const char *pszLineSpacing = CPLGetXMLValue(
1299 17 : psAnnotation.get(),
1300 : "=product.imageAnnotation.imageInformation.azimuthPixelSpacing",
1301 : "UNK");
1302 17 : poDS->SetMetadataItem("LINE_SPACING", pszLineSpacing);
1303 :
1304 : /* --------------------------------------------------------------------
1305 : */
1306 : /* Form full filename (path of manifest.safe + measurement
1307 : * file). */
1308 : /* --------------------------------------------------------------------
1309 : */
1310 : const std::string osFullFilename(
1311 17 : CPLFormFilenameSafe(osPath, pszMeasurement, nullptr));
1312 :
1313 : /* --------------------------------------------------------------------
1314 : */
1315 : /* Try and open the file. */
1316 : /* --------------------------------------------------------------------
1317 : */
1318 0 : std::unique_ptr<GDALDataset> poBandFile;
1319 : {
1320 17 : CPLTurnFailureIntoWarningBackuper oErrorsToWarnings{};
1321 34 : poBandFile = std::unique_ptr<GDALDataset>(
1322 : GDALDataset::Open(osFullFilename.c_str(),
1323 17 : GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR));
1324 : }
1325 :
1326 17 : if (poBandFile == nullptr)
1327 : {
1328 : // NOP
1329 : }
1330 17 : else if (poBandFile->GetRasterCount() == 0)
1331 : {
1332 0 : poBandFile.reset();
1333 : }
1334 : else
1335 : {
1336 34 : poDS->papszExtraFiles =
1337 17 : CSLAddString(poDS->papszExtraFiles, osAnnotationFilePath);
1338 34 : poDS->papszExtraFiles =
1339 17 : CSLAddString(poDS->papszExtraFiles, osCalibrationFilePath);
1340 34 : poDS->papszExtraFiles =
1341 17 : CSLAddString(poDS->papszExtraFiles, osFullFilename.c_str());
1342 : /* --------------------------------------------------------------------
1343 : */
1344 : /* Collect Annotation Processing Information */
1345 : /* --------------------------------------------------------------------
1346 : */
1347 17 : CPLXMLNode *psProcessingInfo = CPLGetXMLNode(
1348 : psAnnotation.get(),
1349 : "=product.imageAnnotation.processingInformation");
1350 :
1351 17 : if (psProcessingInfo != nullptr)
1352 : {
1353 28 : OGRSpatialReference oLL, oPrj;
1354 :
1355 : const char *pszEllipsoidName =
1356 14 : CPLGetXMLValue(psProcessingInfo, "ellipsoidName", "");
1357 14 : const double minor_axis = CPLAtof(CPLGetXMLValue(
1358 : psProcessingInfo, "ellipsoidSemiMinorAxis", "0.0"));
1359 14 : const double major_axis = CPLAtof(CPLGetXMLValue(
1360 : psProcessingInfo, "ellipsoidSemiMajorAxis", "0.0"));
1361 :
1362 14 : if (EQUAL(pszEllipsoidName, "") || (minor_axis == 0.0) ||
1363 : (major_axis == 0.0))
1364 : {
1365 0 : CPLError(CE_Warning, CPLE_AppDefined,
1366 : "Warning- incomplete"
1367 : " ellipsoid information. Using wgs-84 "
1368 : "parameters.\n");
1369 0 : oLL.SetWellKnownGeogCS("WGS84");
1370 0 : oPrj.SetWellKnownGeogCS("WGS84");
1371 : }
1372 14 : else if (EQUAL(pszEllipsoidName, "WGS84"))
1373 : {
1374 14 : oLL.SetWellKnownGeogCS("WGS84");
1375 14 : oPrj.SetWellKnownGeogCS("WGS84");
1376 : }
1377 : else
1378 : {
1379 0 : const double inv_flattening =
1380 0 : major_axis / (major_axis - minor_axis);
1381 0 : oLL.SetGeogCS("", "", pszEllipsoidName, major_axis,
1382 : inv_flattening);
1383 0 : oPrj.SetGeogCS("", "", pszEllipsoidName, major_axis,
1384 : inv_flattening);
1385 : }
1386 :
1387 14 : oLL.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1388 14 : poDS->m_oGCPSRS = std::move(oLL);
1389 : }
1390 :
1391 : /* --------------------------------------------------------------------
1392 : */
1393 : /* Collect GCPs. */
1394 : /* --------------------------------------------------------------------
1395 : */
1396 17 : CPLXMLNode *psGeoGrid = CPLGetXMLNode(
1397 : psAnnotation.get(),
1398 : "=product.geolocationGrid.geolocationGridPointList");
1399 :
1400 17 : if (psGeoGrid != nullptr)
1401 : {
1402 17 : if (poDS->nGCPCount > 0)
1403 : {
1404 6 : GDALDeinitGCPs(poDS->nGCPCount, poDS->pasGCPList);
1405 6 : CPLFree(poDS->pasGCPList);
1406 : }
1407 :
1408 : /* count GCPs */
1409 17 : poDS->nGCPCount = 0;
1410 :
1411 17 : for (CPLXMLNode *psNode = psGeoGrid->psChild;
1412 52 : psNode != nullptr; psNode = psNode->psNext)
1413 : {
1414 35 : if (EQUAL(psNode->pszValue, "geolocationGridPoint"))
1415 18 : poDS->nGCPCount++;
1416 : }
1417 :
1418 17 : if (poDS->nGCPCount > 0)
1419 : {
1420 34 : poDS->pasGCPList = static_cast<GDAL_GCP *>(
1421 17 : CPLCalloc(sizeof(GDAL_GCP), poDS->nGCPCount));
1422 :
1423 17 : poDS->nGCPCount = 0;
1424 :
1425 17 : for (CPLXMLNode *psNode = psGeoGrid->psChild;
1426 52 : psNode != nullptr; psNode = psNode->psNext)
1427 : {
1428 : GDAL_GCP *psGCP =
1429 35 : poDS->pasGCPList + poDS->nGCPCount;
1430 :
1431 35 : if (!EQUAL(psNode->pszValue,
1432 : "geolocationGridPoint"))
1433 17 : continue;
1434 :
1435 18 : poDS->nGCPCount++;
1436 :
1437 : char szID[32];
1438 18 : snprintf(szID, sizeof(szID), "%d", poDS->nGCPCount);
1439 18 : psGCP->pszId = CPLStrdup(szID);
1440 18 : psGCP->pszInfo = CPLStrdup("");
1441 18 : psGCP->dfGCPPixel =
1442 18 : CPLAtof(CPLGetXMLValue(psNode, "pixel", "0"));
1443 18 : psGCP->dfGCPLine =
1444 18 : CPLAtof(CPLGetXMLValue(psNode, "line", "0"));
1445 18 : psGCP->dfGCPX = CPLAtof(
1446 : CPLGetXMLValue(psNode, "longitude", ""));
1447 18 : psGCP->dfGCPY =
1448 18 : CPLAtof(CPLGetXMLValue(psNode, "latitude", ""));
1449 18 : psGCP->dfGCPZ =
1450 18 : CPLAtof(CPLGetXMLValue(psNode, "height", ""));
1451 : }
1452 : }
1453 : }
1454 :
1455 : // Create bands
1456 17 : if (EQUAL(osProductType.c_str(), "SLC"))
1457 : {
1458 : // only add bands if no subDS or uncalibrated subdataset
1459 : // with complex data. (Calibrated will always be intensity
1460 : // only)
1461 3 : if (!bCalibrated &&
1462 0 : (eRequestType == UNKNOWN || eRequestType == COMPLEX))
1463 : {
1464 1 : poBandFile->MarkAsShared();
1465 : SAFESLCRasterBand *poBand0 = new SAFESLCRasterBand(
1466 1 : poDS.get(), eDataType, osSwath, osPolarization,
1467 1 : std::move(poBandFile), SAFESLCRasterBand::COMPLEX);
1468 1 : poDS->SetBand(poDS->GetRasterCount() + 1, poBand0);
1469 : }
1470 2 : else if (eRequestType == INTENSITY) // Intensity
1471 : {
1472 : SAFESLCRasterBand *poBand1 = new SAFESLCRasterBand(
1473 2 : poDS.get(), eDataType, osSwath, osPolarization,
1474 2 : std::move(poBandFile),
1475 2 : SAFESLCRasterBand::INTENSITY);
1476 2 : poDS->SetBand(poDS->GetRasterCount() + 1, poBand1);
1477 : }
1478 : }
1479 14 : else if (!bCalibrated &&
1480 4 : (eRequestType == UNKNOWN || eRequestType == AMPLITUDE))
1481 : {
1482 : SAFERasterBand *poBand = new SAFERasterBand(
1483 14 : poDS.get(), eDataType, osSwath, osPolarization,
1484 14 : std::move(poBandFile));
1485 14 : poDS->SetBand(poDS->GetRasterCount() + 1, poBand);
1486 : }
1487 0 : else if (bCalibrated &&
1488 0 : (eRequestType == UNKNOWN || eRequestType == COMPLEX))
1489 : {
1490 : auto poBand = std::make_unique<SAFECalibratedRasterBand>(
1491 0 : poDS.get(), eDataType, osSwath, osPolarization,
1492 0 : std::move(poBandFile), osCalibrationFilePath,
1493 0 : eCalibrationType);
1494 0 : if (!poBand->ReadLUT())
1495 : {
1496 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1497 : "Reading calibration LUT(s) failed: %s.",
1498 : osCalibrationFilePath.c_str());
1499 0 : return nullptr;
1500 : }
1501 0 : poDS->SetBand(poDS->GetRasterCount() + 1,
1502 0 : std::move(poBand));
1503 : }
1504 : }
1505 : }
1506 : }
1507 :
1508 : // loop through all Swath/pols to add subdatasets
1509 11 : if (!bIsSubDS)
1510 : {
1511 : const CPLString aosCalibrationValues[4] = {"SIGMA0", "BETA0", "GAMMA",
1512 36 : "UNCALIB"};
1513 : const CPLString aosDataUnitValues[3] = {"AMPLITUDE", "COMPLEX",
1514 30 : "INTENSITY"};
1515 6 : if (!isWave)
1516 : {
1517 10 : for (const auto &iterSwath : oMapSwaths2Pols)
1518 : {
1519 10 : CPLString osSubDS1 = iterSwath.first;
1520 10 : CPLString osSubDS2;
1521 :
1522 15 : for (const auto &pol : iterSwath.second)
1523 : {
1524 10 : if (!osSubDS2.empty())
1525 5 : osSubDS2 += "+";
1526 10 : osSubDS2 += pol;
1527 : // Create single band or multiband complex SubDataset
1528 10 : int i = 0;
1529 10 : if (bIsSLC)
1530 : {
1531 0 : for (i = 0; i < 3; i++)
1532 : {
1533 0 : CPLString osCalibTemp = aosCalibrationValues[i];
1534 0 : poDS->AddSubDataset(
1535 : CPLSPrintf("SENTINEL1_CALIB:%s:%s:%s_%s:%s",
1536 : osCalibTemp.c_str(),
1537 : osMDFilename.c_str(),
1538 : osSubDS1.c_str(), pol.c_str(),
1539 : aosDataUnitValues[2].c_str()),
1540 : CPLSPrintf("Single band with %s swath and %s "
1541 : "polarization and %s calibration",
1542 : osSubDS1.c_str(), pol.c_str(),
1543 : osCalibTemp.c_str()));
1544 : }
1545 :
1546 0 : CPLString osCalibTemp = aosCalibrationValues[i];
1547 0 : poDS->AddSubDataset(
1548 : CPLSPrintf("SENTINEL1_CALIB:%s:%s:%s_%s:%s",
1549 : osCalibTemp.c_str(),
1550 : osMDFilename.c_str(), osSubDS1.c_str(),
1551 : pol.c_str(),
1552 : aosDataUnitValues[1].c_str()),
1553 : CPLSPrintf("Single band with %s swath and %s "
1554 : "polarization and %s calibration",
1555 : osSubDS1.c_str(), pol.c_str(),
1556 : osCalibTemp.c_str()));
1557 0 : poDS->AddSubDataset(
1558 : CPLSPrintf("SENTINEL1_CALIB:%s:%s:%s_%s:%s",
1559 : osCalibTemp.c_str(),
1560 : osMDFilename.c_str(), osSubDS1.c_str(),
1561 : pol.c_str(),
1562 : aosDataUnitValues[2].c_str()),
1563 : CPLSPrintf("Single band with %s swath and %s "
1564 : "polarization and %s calibration",
1565 : osSubDS1.c_str(), pol.c_str(),
1566 : osCalibTemp.c_str()));
1567 : }
1568 : else
1569 : {
1570 10 : i = 3;
1571 10 : CPLString osCalibTemp = aosCalibrationValues[i];
1572 :
1573 10 : poDS->AddSubDataset(
1574 : CPLSPrintf("SENTINEL1_CALIB:%s:%s:%s_%s:%s",
1575 : osCalibTemp.c_str(),
1576 : osMDFilename.c_str(), osSubDS1.c_str(),
1577 : pol.c_str(),
1578 : aosDataUnitValues[0].c_str()),
1579 : CPLSPrintf("Single band with %s swath and %s "
1580 : "polarization and %s calibration",
1581 : osSubDS1.c_str(), pol.c_str(),
1582 : osCalibTemp.c_str()));
1583 : }
1584 : }
1585 :
1586 5 : if (iterSwath.second.size() > 1)
1587 : {
1588 : // Create single band subdataset with all polarizations
1589 5 : int i = 0;
1590 5 : if (bIsSLC)
1591 : {
1592 0 : for (i = 0; i < 3; i++)
1593 : {
1594 0 : CPLString osCalibTemp = aosCalibrationValues[i];
1595 :
1596 0 : poDS->AddSubDataset(
1597 : CPLSPrintf("SENTINEL1_CALIB:%s:%s:%s:%s",
1598 : osCalibTemp.c_str(),
1599 : osMDFilename.c_str(),
1600 : osSubDS1.c_str(),
1601 : aosDataUnitValues[2].c_str()),
1602 : CPLSPrintf(
1603 : "%s swath with all polarizations (%s) as "
1604 : "bands and %s calibration",
1605 : osSubDS1.c_str(), osSubDS2.c_str(),
1606 : osCalibTemp.c_str()));
1607 : }
1608 :
1609 0 : CPLString osCalibTemp = aosCalibrationValues[i];
1610 0 : poDS->AddSubDataset(
1611 : CPLSPrintf("SENTINEL1_CALIB:%s:%s:%s:%s",
1612 : osCalibTemp.c_str(),
1613 : osMDFilename.c_str(), osSubDS1.c_str(),
1614 : aosDataUnitValues[1].c_str()),
1615 : CPLSPrintf(
1616 : "%s swath with all polarizations (%s) as "
1617 : "bands and %s calibration",
1618 : osSubDS1.c_str(), osSubDS2.c_str(),
1619 : osCalibTemp.c_str()));
1620 0 : poDS->AddSubDataset(
1621 : CPLSPrintf("SENTINEL1_CALIB:%s:%s:%s:%s",
1622 : osCalibTemp.c_str(),
1623 : osMDFilename.c_str(), osSubDS1.c_str(),
1624 : aosDataUnitValues[2].c_str()),
1625 : CPLSPrintf(
1626 : "%s swath with all polarizations (%s) as "
1627 : "bands and %s calibration",
1628 : osSubDS1.c_str(), osSubDS2.c_str(),
1629 : osCalibTemp.c_str()));
1630 : }
1631 : else
1632 : {
1633 5 : i = 3;
1634 5 : CPLString osCalibTemp = aosCalibrationValues[i];
1635 5 : poDS->AddSubDataset(
1636 : CPLSPrintf("SENTINEL1_CALIB:%s:%s:%s:%s",
1637 : osCalibTemp.c_str(),
1638 : osMDFilename.c_str(), osSubDS1.c_str(),
1639 : aosDataUnitValues[0].c_str()),
1640 : CPLSPrintf(
1641 : "%s swath with all polarizations (%s) as "
1642 : "bands and %s calibration",
1643 : osSubDS1.c_str(), osSubDS2.c_str(),
1644 : osCalibTemp.c_str()));
1645 : }
1646 : }
1647 : }
1648 : }
1649 : else
1650 : {
1651 3 : for (const CPLString &osImgSwPol : oImageNumberSwPol)
1652 : {
1653 4 : CPLString osImgSwPolTmp = osImgSwPol;
1654 2 : const char *pszImage = strchr(osImgSwPolTmp.c_str(), ' ');
1655 4 : CPLString osImage, osSwath, osPolarization;
1656 2 : if (pszImage != nullptr)
1657 : {
1658 2 : osImage = osImgSwPolTmp;
1659 2 : osImage.resize(pszImage - osImgSwPolTmp.c_str());
1660 2 : osImgSwPolTmp = pszImage + strlen(" ");
1661 2 : const char *pszSwath = strchr(osImgSwPolTmp.c_str(), ' ');
1662 2 : if (pszSwath != nullptr)
1663 : {
1664 2 : osSwath = osImgSwPolTmp;
1665 2 : osSwath.resize(pszSwath - osImgSwPolTmp.c_str());
1666 2 : osPolarization = pszSwath + strlen(" ");
1667 2 : int i = 0;
1668 :
1669 2 : if (bIsSLC)
1670 : {
1671 8 : for (i = 0; i < 3; i++)
1672 : {
1673 6 : CPLString osCalibTemp = aosCalibrationValues[i];
1674 :
1675 6 : poDS->AddSubDataset(
1676 : CPLSPrintf(
1677 : "SENTINEL1_CALIB:%s:%s:%s_%s_%s:%s",
1678 : osCalibTemp.c_str(),
1679 : osMDFilename.c_str(), osSwath.c_str(),
1680 : osPolarization.c_str(), osImage.c_str(),
1681 : aosDataUnitValues[2].c_str()),
1682 : CPLSPrintf(
1683 : "Single band with %s swath and %s "
1684 : "polarization and %s calibration",
1685 : osSwath.c_str(), osPolarization.c_str(),
1686 : osCalibTemp.c_str()));
1687 : }
1688 :
1689 2 : CPLString osCalibTemp = aosCalibrationValues[i];
1690 :
1691 2 : poDS->AddSubDataset(
1692 : CPLSPrintf(
1693 : "SENTINEL1_CALIB:%s:%s:%s_%s_%s:%s",
1694 : osCalibTemp.c_str(), osMDFilename.c_str(),
1695 : osSwath.c_str(), osPolarization.c_str(),
1696 : osImage.c_str(),
1697 : aosDataUnitValues[1].c_str()),
1698 : CPLSPrintf("Single band with %s swath and %s "
1699 : "polarization and %s calibration",
1700 : osSwath.c_str(),
1701 : osPolarization.c_str(),
1702 : osCalibTemp.c_str()));
1703 :
1704 2 : poDS->AddSubDataset(
1705 : CPLSPrintf(
1706 : "SENTINEL1_CALIB:%s:%s:%s_%s_%s:%s",
1707 : osCalibTemp.c_str(), osMDFilename.c_str(),
1708 : osSwath.c_str(), osPolarization.c_str(),
1709 : osImage.c_str(),
1710 : aosDataUnitValues[2].c_str()),
1711 : CPLSPrintf("Single band with %s swath and %s "
1712 : "polarization and %s calibration",
1713 : osSwath.c_str(),
1714 : osPolarization.c_str(),
1715 : osCalibTemp.c_str()));
1716 : }
1717 : else
1718 : {
1719 0 : i = 3;
1720 0 : CPLString osCalibTemp = aosCalibrationValues[i];
1721 :
1722 0 : poDS->AddSubDataset(
1723 : CPLSPrintf(
1724 : "SENTINEL1_CALIB:%s:%s:%s_%s_%s:%s",
1725 : osCalibTemp.c_str(), osMDFilename.c_str(),
1726 : osSwath.c_str(), osPolarization.c_str(),
1727 : osImage.c_str(),
1728 : aosDataUnitValues[0].c_str()),
1729 : CPLSPrintf("Single band with %s swath and %s "
1730 : "polarization and %s calibration",
1731 : osSwath.c_str(),
1732 : osPolarization.c_str(),
1733 : osCalibTemp.c_str()));
1734 : }
1735 : }
1736 : }
1737 : }
1738 : }
1739 : }
1740 11 : if (poDS->GetRasterCount() == 0)
1741 : {
1742 0 : CPLError(CE_Failure, CPLE_OpenFailed, "Measurement bands not found.");
1743 0 : return nullptr;
1744 : }
1745 :
1746 : /* -------------------------------------------------------------------- */
1747 : /* Collect more metadata elements */
1748 : /* -------------------------------------------------------------------- */
1749 :
1750 : /* -------------------------------------------------------------------- */
1751 : /* Platform information */
1752 : /* -------------------------------------------------------------------- */
1753 : const CPLXMLNode *psPlatformAttrs =
1754 11 : SAFEDataset::GetMetaDataObject(psMetaDataObjects, "platform");
1755 :
1756 11 : if (psPlatformAttrs != nullptr)
1757 : {
1758 : const char *pszItem =
1759 11 : CPLGetXMLValue(psPlatformAttrs,
1760 : "metadataWrap.xmlData.safe:platform"
1761 : ".safe:familyName",
1762 : "");
1763 11 : poDS->SetMetadataItem("SATELLITE_IDENTIFIER", pszItem);
1764 :
1765 : pszItem =
1766 11 : CPLGetXMLValue(psPlatformAttrs,
1767 : "metadataWrap.xmlData.safe:platform"
1768 : ".safe:instrument.safe:familyName.abbreviation",
1769 : "");
1770 11 : poDS->SetMetadataItem("SENSOR_IDENTIFIER", pszItem);
1771 :
1772 11 : pszItem = CPLGetXMLValue(psPlatformAttrs,
1773 : "metadataWrap.xmlData.safe:platform"
1774 : ".safe:instrument.safe:extension"
1775 : ".s1sarl1:instrumentMode.s1sarl1:mode",
1776 : "UNK");
1777 11 : poDS->SetMetadataItem("BEAM_MODE", pszItem);
1778 :
1779 11 : pszItem = CPLGetXMLValue(psPlatformAttrs,
1780 : "metadataWrap.xmlData.safe:platform"
1781 : ".safe:instrument.safe:extension"
1782 : ".s1sarl1:instrumentMode.s1sarl1:swath",
1783 : "UNK");
1784 11 : poDS->SetMetadataItem("BEAM_SWATH", pszItem);
1785 : }
1786 :
1787 : /* -------------------------------------------------------------------- */
1788 : /* Acquisition Period information */
1789 : /* -------------------------------------------------------------------- */
1790 : const CPLXMLNode *psAcquisitionAttrs =
1791 11 : SAFEDataset::GetMetaDataObject(psMetaDataObjects, "acquisitionPeriod");
1792 :
1793 11 : if (psAcquisitionAttrs != nullptr)
1794 : {
1795 : const char *pszItem =
1796 11 : CPLGetXMLValue(psAcquisitionAttrs,
1797 : "metadataWrap.xmlData.safe:acquisitionPeriod"
1798 : ".safe:startTime",
1799 : "UNK");
1800 11 : poDS->SetMetadataItem("ACQUISITION_START_TIME", pszItem);
1801 11 : pszItem = CPLGetXMLValue(psAcquisitionAttrs,
1802 : "metadataWrap.xmlData.safe:acquisitionPeriod"
1803 : ".safe:stopTime",
1804 : "UNK");
1805 11 : poDS->SetMetadataItem("ACQUISITION_STOP_TIME", pszItem);
1806 : }
1807 :
1808 : /* -------------------------------------------------------------------- */
1809 : /* Processing information */
1810 : /* -------------------------------------------------------------------- */
1811 : const CPLXMLNode *psProcessingAttrs =
1812 11 : SAFEDataset::GetMetaDataObject(psMetaDataObjects, "processing");
1813 :
1814 11 : if (psProcessingAttrs != nullptr)
1815 : {
1816 11 : const char *pszItem = CPLGetXMLValue(
1817 : psProcessingAttrs,
1818 : "metadataWrap.xmlData.safe:processing.safe:facility.name", "UNK");
1819 11 : poDS->SetMetadataItem("FACILITY_IDENTIFIER", pszItem);
1820 : }
1821 :
1822 : /* -------------------------------------------------------------------- */
1823 : /* Measurement Orbit Reference information */
1824 : /* -------------------------------------------------------------------- */
1825 11 : const CPLXMLNode *psOrbitAttrs = SAFEDataset::GetMetaDataObject(
1826 : psMetaDataObjects, "measurementOrbitReference");
1827 :
1828 11 : if (psOrbitAttrs != nullptr)
1829 : {
1830 : const char *pszItem =
1831 11 : CPLGetXMLValue(psOrbitAttrs,
1832 : "metadataWrap.xmlData.safe:orbitReference"
1833 : ".safe:orbitNumber",
1834 : "UNK");
1835 11 : poDS->SetMetadataItem("ORBIT_NUMBER", pszItem);
1836 11 : pszItem = CPLGetXMLValue(psOrbitAttrs,
1837 : "metadataWrap.xmlData.safe:orbitReference"
1838 : ".safe:extension.s1:orbitProperties.s1:pass",
1839 : "UNK");
1840 11 : poDS->SetMetadataItem("ORBIT_DIRECTION", pszItem);
1841 : }
1842 :
1843 : /* -------------------------------------------------------------------- */
1844 : /* Footprint */
1845 : /* -------------------------------------------------------------------- */
1846 11 : const CPLXMLNode *psFrameSet = SAFEDataset::GetMetaDataObject(
1847 : psMetaDataObjects, "measurementFrameSet");
1848 :
1849 11 : if (psFrameSet)
1850 : {
1851 11 : const auto psFootPrint = CPLGetXMLNode(psFrameSet, "metadataWrap."
1852 : "xmlData."
1853 : "safe:frameSet."
1854 : "safe:frame."
1855 : "safe:footPrint");
1856 11 : if (psFootPrint)
1857 : {
1858 : const char *pszSRSName =
1859 11 : CPLGetXMLValue(psFootPrint, "srsName", nullptr);
1860 : const char *pszCoordinates =
1861 11 : CPLGetXMLValue(psFootPrint, "gml:coordinates", nullptr);
1862 11 : if (pszSRSName &&
1863 11 : EQUAL(pszSRSName,
1864 11 : "http://www.opengis.net/gml/srs/epsg.xml#4326") &&
1865 : pszCoordinates)
1866 : {
1867 : const CPLStringList aosValues(
1868 22 : CSLTokenizeString2(pszCoordinates, " ,", 0));
1869 11 : if (aosValues.size() == 8)
1870 : {
1871 22 : poDS->SetMetadataItem(
1872 : "FOOTPRINT",
1873 : CPLSPrintf("POLYGON((%s %s,%s %s,%s %s,%s %s, %s %s))",
1874 : aosValues[1], aosValues[0], aosValues[3],
1875 : aosValues[2], aosValues[5], aosValues[4],
1876 : aosValues[7], aosValues[6], aosValues[1],
1877 11 : aosValues[0]));
1878 : }
1879 : }
1880 : }
1881 : }
1882 :
1883 : /* -------------------------------------------------------------------- */
1884 : /* Initialize any PAM information. */
1885 : /* -------------------------------------------------------------------- */
1886 22 : const CPLString osDescription = osMDFilename;
1887 :
1888 : /* -------------------------------------------------------------------- */
1889 : /* Initialize any PAM information. */
1890 : /* -------------------------------------------------------------------- */
1891 11 : poDS->SetDescription(osDescription);
1892 :
1893 11 : poDS->SetPhysicalFilename(osMDFilename);
1894 11 : if (!osSubdatasetName.empty())
1895 : {
1896 5 : poDS->SetDescription(poOpenInfo->pszFilename);
1897 5 : poDS->SetSubdatasetName(osSubdatasetName);
1898 : }
1899 :
1900 11 : poDS->TryLoadXML();
1901 :
1902 : /* -------------------------------------------------------------------- */
1903 : /* Check for overviews. */
1904 : /* -------------------------------------------------------------------- */
1905 11 : poDS->oOvManager.Initialize(poDS.get(), ":::VIRTUAL:::");
1906 :
1907 11 : return poDS.release();
1908 : }
1909 :
1910 : /************************************************************************/
1911 : /* AddSubDataset() */
1912 : /************************************************************************/
1913 25 : void SAFEDataset::AddSubDataset(const CPLString &osName,
1914 : const CPLString &osDesc)
1915 : {
1916 25 : ++m_nSubDSNum;
1917 25 : papszSubDatasets = CSLAddNameValue(
1918 : papszSubDatasets, CPLSPrintf("SUBDATASET_%d_NAME", m_nSubDSNum),
1919 : osName.c_str());
1920 25 : papszSubDatasets = CSLAddNameValue(
1921 : papszSubDatasets, CPLSPrintf("SUBDATASET_%d_DESC", m_nSubDSNum),
1922 : osDesc.c_str());
1923 25 : }
1924 :
1925 : /************************************************************************/
1926 : /* GetGCPCount() */
1927 : /************************************************************************/
1928 :
1929 2 : int SAFEDataset::GetGCPCount()
1930 : {
1931 2 : return nGCPCount;
1932 : }
1933 :
1934 : /************************************************************************/
1935 : /* GetGCPSpatialRef() */
1936 : /************************************************************************/
1937 :
1938 0 : const OGRSpatialReference *SAFEDataset::GetGCPSpatialRef() const
1939 :
1940 : {
1941 0 : return m_oGCPSRS.IsEmpty() ? nullptr : &m_oGCPSRS;
1942 : }
1943 :
1944 : /************************************************************************/
1945 : /* GetGCPs() */
1946 : /************************************************************************/
1947 :
1948 2 : const GDAL_GCP *SAFEDataset::GetGCPs()
1949 :
1950 : {
1951 2 : return pasGCPList;
1952 : }
1953 :
1954 : /************************************************************************/
1955 : /* GetMetadataDomainList() */
1956 : /************************************************************************/
1957 :
1958 0 : char **SAFEDataset::GetMetadataDomainList()
1959 : {
1960 0 : return BuildMetadataDomainList(GDALDataset::GetMetadataDomainList(), TRUE,
1961 0 : "SUBDATASETS", nullptr);
1962 : }
1963 :
1964 : /************************************************************************/
1965 : /* GetMetadata() */
1966 : /************************************************************************/
1967 :
1968 3 : char **SAFEDataset::GetMetadata(const char *pszDomain)
1969 : {
1970 3 : if (pszDomain != nullptr && STARTS_WITH_CI(pszDomain, "SUBDATASETS"))
1971 3 : return papszSubDatasets;
1972 :
1973 0 : return GDALDataset::GetMetadata(pszDomain);
1974 : }
1975 :
1976 : /************************************************************************/
1977 : /* GDALRegister_SAFE() */
1978 : /************************************************************************/
1979 :
1980 2024 : void GDALRegister_SAFE()
1981 : {
1982 2024 : if (GDALGetDriverByName("SAFE") != nullptr)
1983 283 : return;
1984 :
1985 1741 : GDALDriver *poDriver = new GDALDriver();
1986 :
1987 1741 : poDriver->SetDescription("SAFE");
1988 1741 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1989 1741 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1990 1741 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Sentinel-1 SAR SAFE Product");
1991 1741 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/safe.html");
1992 :
1993 1741 : poDriver->pfnOpen = SAFEDataset::Open;
1994 1741 : poDriver->pfnIdentify = SAFEDataset::Identify;
1995 :
1996 1741 : GetGDALDriverManager()->RegisterDriver(poDriver);
1997 : }
1998 :
1999 : /************************************************************************/
2000 : /* Azimuth time handling functions */
2001 : /************************************************************************/
2002 :
2003 : SAFECalibratedRasterBand::TimePoint
2004 0 : SAFECalibratedRasterBand::getTimePoint(const char *pszTime)
2005 : {
2006 : int nYear, nMonth, nDay, nHour, nMinute, nSeconds;
2007 : long nMicroSeconds;
2008 0 : sscanf(pszTime, "%d-%d-%dT%d:%d:%d.%ld", &nYear, &nMonth, &nDay, &nHour,
2009 : &nMinute, &nSeconds, &nMicroSeconds);
2010 :
2011 : struct tm oTm;
2012 0 : oTm.tm_sec = nSeconds;
2013 0 : oTm.tm_min = nMinute;
2014 0 : oTm.tm_hour = nHour;
2015 0 : oTm.tm_mday = nDay;
2016 0 : oTm.tm_mon = nMonth - 1;
2017 0 : oTm.tm_year = nYear - 1900;
2018 0 : oTm.tm_isdst = -1;
2019 :
2020 0 : std::time_t oTt = static_cast<std::time_t>(CPLYMDHMSToUnixTime(&oTm));
2021 :
2022 0 : TimePoint oTp = std::chrono::system_clock::from_time_t(oTt);
2023 0 : TimePoint oTp1 = oTp + std::chrono::microseconds(nMicroSeconds);
2024 0 : return oTp1;
2025 : }
2026 :
2027 0 : double SAFECalibratedRasterBand::getTimeDiff(TimePoint oT1, TimePoint oT2)
2028 : {
2029 0 : std::chrono::duration<double> oResult = (oT2 - oT1);
2030 0 : return oResult.count();
2031 : }
2032 :
2033 : SAFECalibratedRasterBand::TimePoint
2034 0 : SAFECalibratedRasterBand::getazTime(TimePoint oStart, TimePoint oStop,
2035 : long nNumOfLines, int nOffset)
2036 : {
2037 0 : double dfTemp = getTimeDiff(oStart, oStop);
2038 0 : dfTemp /= (nNumOfLines - 1);
2039 0 : unsigned long nTimeDiffMicro = static_cast<unsigned long>(dfTemp * 1000000);
2040 : TimePoint oResult =
2041 0 : oStart + (nOffset * std::chrono::microseconds(nTimeDiffMicro));
2042 0 : return oResult;
2043 : }
2044 :
2045 : /************************************************************************/
2046 : /* Utility functions used in interpolation */
2047 : /************************************************************************/
2048 :
2049 0 : int SAFECalibratedRasterBand::getCalibrationVectorIndex(int nLineNo)
2050 : {
2051 0 : for (size_t i = 1; i < m_anLineLUT.size(); i++)
2052 : {
2053 0 : if (nLineNo < m_anLineLUT[i])
2054 0 : return static_cast<int>(i - 1);
2055 : }
2056 0 : return 0;
2057 : }
2058 :
2059 0 : int SAFECalibratedRasterBand::getPixelIndex(int nPixelNo)
2060 : {
2061 0 : for (int i = 1; i < m_nNumPixels; i++)
2062 : {
2063 0 : if (nPixelNo < m_anPixelLUT[i])
2064 0 : return i - 1;
2065 : }
2066 0 : return 0;
2067 : }
|
