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