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