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