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