Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: Sentinel2 products
5 : * Author: Even Rouault, <even.rouault at spatialys.com>
6 : * Funded by: Centre National d'Etudes Spatiales (CNES)
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2015, Even Rouault, <even.rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_minixml.h"
15 : #include "cpl_string.h"
16 : #include "gdal_pam.h"
17 : #include "ogr_spatialref.h"
18 : #include "ogr_geometry.h"
19 : #include "gdaljp2metadata.h"
20 : #include "../vrt/vrtdataset.h"
21 :
22 : #include <algorithm>
23 : #include <map>
24 : #include <set>
25 : #include <vector>
26 :
27 : #ifdef HAVE_UNISTD_H
28 : #include <unistd.h>
29 : #endif
30 :
31 : #ifndef STARTS_WITH_CI
32 : #define STARTS_WITH_CI(a, b) EQUALN(a, b, strlen(b))
33 : #endif
34 :
35 : #define DIGIT_ZERO '0'
36 :
37 : typedef enum
38 : {
39 : SENTINEL2_L1B,
40 : SENTINEL2_L1C,
41 : SENTINEL2_L2A
42 : } SENTINEL2Level;
43 :
44 : typedef enum
45 : {
46 : MSI2Ap,
47 : MSI2A
48 : } SENTINEL2ProductType;
49 :
50 : typedef struct
51 : {
52 : const char *pszBandName;
53 : int nResolution; /* meters */
54 : int nWaveLength; /* nanometers */
55 : int nBandWidth; /* nanometers */
56 : GDALColorInterp eColorInterp;
57 : } SENTINEL2BandDescription;
58 :
59 : constexpr int RES_10M = 10;
60 : constexpr int RES_20M = 20;
61 : constexpr int RES_60M = 60;
62 :
63 : /* clang-format off */
64 : static const SENTINEL2BandDescription asBandDesc[] = {
65 : {"B1", RES_60M, 443, 20, GCI_CoastalBand},
66 : {"B2", RES_10M, 490, 65, GCI_BlueBand},
67 : {"B3", RES_10M, 560, 35, GCI_GreenBand},
68 : {"B4", RES_10M, 665, 30, GCI_RedBand},
69 : {"B5", RES_20M, 705, 15, GCI_RedEdgeBand}, // rededge071
70 : {"B6", RES_20M, 740, 15, GCI_RedEdgeBand}, // rededge075
71 : {"B7", RES_20M, 783, 20, GCI_RedEdgeBand}, // rededge078
72 : {"B8", RES_10M, 842, 115, GCI_NIRBand}, // nir
73 : {"B8A", RES_20M, 865, 20, GCI_NIRBand}, // nir08
74 : {"B9", RES_60M, 945, 20, GCI_NIRBand}, // nir09
75 : {"B10", RES_60M, 1375, 30, GCI_OtherIRBand}, // cirrus
76 : {"B11", RES_20M, 1610, 90, GCI_SWIRBand}, // swir16
77 : {"B12", RES_20M, 2190, 180, GCI_SWIRBand}, // swir11
78 : };
79 :
80 : /* clang-format on */
81 :
82 : typedef enum
83 : {
84 : TL_IMG_DATA, /* Tile is located in IMG_DATA/ */
85 : TL_IMG_DATA_Rxxm, /* Tile is located in IMG_DATA/Rxxm/ */
86 : TL_QI_DATA /* Tile is located in QI_DATA/ */
87 : } SENTINEL2_L2A_Tilelocation;
88 :
89 : struct SENTINEL2_L2A_BandDescription
90 : {
91 : const char *pszBandName = nullptr;
92 : const char *pszBandDescription = nullptr;
93 : int nResolution = 0; /* meters */
94 : SENTINEL2_L2A_Tilelocation eLocation = TL_IMG_DATA;
95 : };
96 :
97 : class L1CSafeCompatGranuleDescription
98 : {
99 : public:
100 : // GRANULE/L1C_T30TXT_A007999_20170102T111441/MTD_TL.xml
101 : CPLString osMTDTLPath{};
102 : // GRANULE/L1C_T30TXT_A007999_20170102T111441/IMG_DATA/T30TXT_20170102T111442_
103 : CPLString osBandPrefixPath{};
104 : };
105 :
106 : static const char *L2A_BandDescription_AOT =
107 : "Aerosol Optical Thickness map (at 550nm)";
108 : static const char *L2A_BandDescription_WVP = "Scene-average Water Vapour map";
109 : static const char *L2A_BandDescription_SCL = "Scene Classification";
110 : static const char *L2A_BandDescription_CLD =
111 : "Raster mask values range from 0 for high confidence clear sky to 100 for "
112 : "high confidence cloudy";
113 : static const char *L2A_BandDescription_SNW =
114 : "Raster mask values range from 0 for high confidence NO snow/ice to 100 "
115 : "for high confidence snow/ice";
116 :
117 : static const SENTINEL2_L2A_BandDescription asL2ABandDesc[] = {
118 : {"AOT", L2A_BandDescription_AOT, RES_10M, TL_IMG_DATA_Rxxm},
119 : {"AOT", L2A_BandDescription_AOT, RES_20M, TL_IMG_DATA_Rxxm},
120 : {"AOT", L2A_BandDescription_AOT, RES_60M, TL_IMG_DATA_Rxxm},
121 : {"WVP", L2A_BandDescription_WVP, RES_10M, TL_IMG_DATA_Rxxm},
122 : {"WVP", L2A_BandDescription_WVP, RES_20M, TL_IMG_DATA_Rxxm},
123 : {"WVP", L2A_BandDescription_WVP, RES_60M, TL_IMG_DATA_Rxxm},
124 : {"SCL", L2A_BandDescription_SCL, RES_20M, TL_IMG_DATA_Rxxm},
125 : {"SCL", L2A_BandDescription_SCL, RES_60M, TL_IMG_DATA_Rxxm},
126 : {"CLD", L2A_BandDescription_CLD, RES_20M, TL_QI_DATA},
127 : {"CLD", L2A_BandDescription_CLD, RES_60M, TL_QI_DATA},
128 : {"SNW", L2A_BandDescription_SNW, RES_20M, TL_QI_DATA},
129 : {"SNW", L2A_BandDescription_SNW, RES_60M, TL_QI_DATA},
130 : };
131 :
132 : static bool SENTINEL2isZipped(const char *pszHeader, int nHeaderBytes);
133 : static const char *SENTINEL2GetOption(GDALOpenInfo *poOpenInfo,
134 : const char *pszName,
135 : const char *pszDefaultVal = nullptr);
136 : static bool SENTINEL2GetTileInfo(const char *pszFilename, int *pnWidth,
137 : int *pnHeight, int *pnBits);
138 :
139 : /************************************************************************/
140 : /* IsS2Prefixed() */
141 : /************************************************************************/
142 :
143 : // IsS2Prefixed(pszStr, "foo") checks that pszStr starts with
144 : // "S2x_foo" where x=A/B/C/D
145 235213 : static bool IsS2Prefixed(const char *pszStr, const char *pszPrefixAfterS2X)
146 : {
147 282 : return pszStr[0] == 'S' && pszStr[1] == '2' && pszStr[2] >= 'A' &&
148 235669 : pszStr[2] <= 'Z' &&
149 174 : (*pszPrefixAfterS2X == 0 ||
150 235368 : STARTS_WITH_CI(pszStr + 3, pszPrefixAfterS2X));
151 : }
152 :
153 : /************************************************************************/
154 : /* SENTINEL2GranuleInfo */
155 : /************************************************************************/
156 :
157 : class SENTINEL2GranuleInfo
158 : {
159 : public:
160 : CPLString osPath{};
161 : CPLString osBandPrefixPath{}; // for Sentinel 2C SafeCompact
162 : double dfMinX = 0, dfMinY = 0, dfMaxX = 0, dfMaxY = 0;
163 : int nWidth = 0, nHeight = 0;
164 : };
165 :
166 : /************************************************************************/
167 : /* ==================================================================== */
168 : /* SENTINEL2Dataset */
169 : /* ==================================================================== */
170 : /************************************************************************/
171 :
172 86 : class SENTINEL2DatasetContainer final : public GDALPamDataset
173 : {
174 : public:
175 43 : SENTINEL2DatasetContainer() = default;
176 : ~SENTINEL2DatasetContainer() override;
177 : };
178 :
179 : SENTINEL2DatasetContainer::~SENTINEL2DatasetContainer() = default;
180 :
181 : class SENTINEL2Dataset final : public VRTDataset
182 : {
183 : std::vector<CPLString> aosNonJP2Files{};
184 :
185 : void AddL1CL2ABandMetadata(SENTINEL2Level eLevel, CPLXMLNode *psRoot,
186 : const std::vector<CPLString> &aosBands);
187 :
188 : static SENTINEL2Dataset *CreateL1CL2ADataset(
189 : SENTINEL2Level eLevel, SENTINEL2ProductType pType, bool bIsSafeCompact,
190 : const std::vector<CPLString> &aosGranuleList,
191 : const std::vector<L1CSafeCompatGranuleDescription>
192 : &aoL1CSafeCompactGranuleList,
193 : std::vector<CPLString> &aosNonJP2Files, int nSubDSPrecision,
194 : bool bIsPreview, bool bIsTCI, int nSubDSEPSGCode, bool bAlpha,
195 : const std::vector<CPLString> &aosBands, int nSaturatedVal,
196 : int nNodataVal, const CPLString &osProductURI);
197 :
198 : public:
199 : SENTINEL2Dataset(int nXSize, int nYSize);
200 : ~SENTINEL2Dataset() override;
201 :
202 : char **GetFileList() override;
203 :
204 : static GDALDataset *Open(GDALOpenInfo *);
205 : static GDALDataset *OpenL1BUserProduct(GDALOpenInfo *);
206 : static GDALDataset *
207 : OpenL1BGranule(const char *pszFilename, CPLXMLNode **ppsRoot = nullptr,
208 : int nResolutionOfInterest = 0,
209 : std::set<CPLString> *poBandSet = nullptr);
210 : static GDALDataset *OpenL1BSubdataset(GDALOpenInfo *);
211 : static GDALDataset *OpenL1BSubdatasetWithGeoloc(GDALOpenInfo *);
212 : static GDALDataset *OpenL1C_L2A(const char *pszFilename,
213 : SENTINEL2Level eLevel);
214 : static GDALDataset *OpenL1CTile(const char *pszFilename,
215 : CPLXMLNode **ppsRootMainMTD = nullptr,
216 : int nResolutionOfInterest = 0,
217 : std::set<CPLString> *poBandSet = nullptr);
218 : static GDALDataset *OpenL1CTileSubdataset(GDALOpenInfo *);
219 : static GDALDataset *OpenL1C_L2ASubdataset(GDALOpenInfo *,
220 : SENTINEL2Level eLevel);
221 :
222 : static int Identify(GDALOpenInfo *);
223 : };
224 :
225 : /************************************************************************/
226 : /* ==================================================================== */
227 : /* SENTINEL2AlphaBand */
228 : /* ==================================================================== */
229 : /************************************************************************/
230 :
231 : class SENTINEL2AlphaBand final : public VRTSourcedRasterBand
232 : {
233 : int m_nSaturatedVal;
234 : int m_nNodataVal;
235 :
236 : public:
237 : SENTINEL2AlphaBand(GDALDataset *poDS, int nBand, GDALDataType eType,
238 : int nXSize, int nYSize, int nSaturatedVal,
239 : int nNodataVal);
240 :
241 : CPLErr IRasterIO(GDALRWFlag, int, int, int, int, void *, int, int,
242 : GDALDataType, GSpacing nPixelSpace, GSpacing nLineSpace,
243 : GDALRasterIOExtraArg *psExtraArg) override;
244 : };
245 :
246 : /************************************************************************/
247 : /* SENTINEL2AlphaBand() */
248 : /************************************************************************/
249 :
250 3 : SENTINEL2AlphaBand::SENTINEL2AlphaBand(GDALDataset *poDSIn, int nBandIn,
251 : GDALDataType eType, int nXSize,
252 : int nYSize, int nSaturatedVal,
253 3 : int nNodataVal)
254 : : VRTSourcedRasterBand(poDSIn, nBandIn, eType, nXSize, nYSize),
255 3 : m_nSaturatedVal(nSaturatedVal), m_nNodataVal(nNodataVal)
256 : {
257 3 : }
258 :
259 : /************************************************************************/
260 : /* IRasterIO() */
261 : /************************************************************************/
262 :
263 29 : CPLErr SENTINEL2AlphaBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
264 : int nXSize, int nYSize, void *pData,
265 : int nBufXSize, int nBufYSize,
266 : GDALDataType eBufType,
267 : GSpacing nPixelSpace, GSpacing nLineSpace,
268 : GDALRasterIOExtraArg *psExtraArg)
269 : {
270 : // Query the first band. Quite arbitrary, but hopefully all bands have
271 : // the same nodata/saturated pixels.
272 29 : CPLErr eErr = poDS->GetRasterBand(1)->RasterIO(
273 : eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
274 : eBufType, nPixelSpace, nLineSpace, psExtraArg);
275 29 : if (eErr == CE_None)
276 : {
277 29 : const char *pszNBITS = GetMetadataItem("NBITS", "IMAGE_STRUCTURE");
278 29 : const int nBits = (pszNBITS) ? atoi(pszNBITS) : 16;
279 29 : const GUInt16 nMaxVal = (nBits > 8 && nBits <= 16)
280 29 : ? static_cast<GUInt16>((1 << nBits) - 1)
281 : : 65535;
282 :
283 : // Replace pixels matching m_nSaturatedVal and m_nNodataVal by 0
284 : // and others by the maxVal.
285 7023 : for (int iY = 0; iY < nBufYSize; iY++)
286 : {
287 24465000 : for (int iX = 0; iX < nBufXSize; iX++)
288 : {
289 : // Optimized path for likely most common case
290 24458000 : if (eBufType == GDT_UInt16)
291 : {
292 12229000 : GUInt16 *panPtr = reinterpret_cast<GUInt16 *>(
293 12229000 : static_cast<GByte *>(pData) + iY * nLineSpace +
294 12229000 : iX * nPixelSpace);
295 12229000 : if (*panPtr == 0 || *panPtr == m_nSaturatedVal ||
296 0 : *panPtr == m_nNodataVal)
297 : {
298 12229000 : *panPtr = 0;
299 : }
300 : else
301 0 : *panPtr = nMaxVal;
302 : }
303 : // Generic path for other datatypes
304 : else
305 : {
306 : double dfVal;
307 12229000 : GDALCopyWords(static_cast<GByte *>(pData) +
308 12229000 : iY * nLineSpace + iX * nPixelSpace,
309 : eBufType, 0, &dfVal, GDT_Float64, 0, 1);
310 12229000 : if (dfVal == 0.0 || dfVal == m_nSaturatedVal ||
311 0 : dfVal == m_nNodataVal)
312 : {
313 12229000 : dfVal = 0;
314 : }
315 : else
316 0 : dfVal = nMaxVal;
317 12229000 : GDALCopyWords(&dfVal, GDT_Float64, 0,
318 : static_cast<GByte *>(pData) +
319 12229000 : iY * nLineSpace + iX * nPixelSpace,
320 : eBufType, 0, 1);
321 : }
322 : }
323 : }
324 : }
325 :
326 29 : return eErr;
327 : }
328 :
329 : /************************************************************************/
330 : /* SENTINEL2Dataset() */
331 : /************************************************************************/
332 :
333 45 : SENTINEL2Dataset::SENTINEL2Dataset(int nXSize, int nYSize)
334 45 : : VRTDataset(nXSize, nYSize)
335 : {
336 45 : poDriver = nullptr;
337 45 : SetWritable(FALSE);
338 45 : }
339 :
340 : /************************************************************************/
341 : /* ~SENTINEL2Dataset() */
342 : /************************************************************************/
343 :
344 90 : SENTINEL2Dataset::~SENTINEL2Dataset()
345 : {
346 90 : }
347 :
348 : /************************************************************************/
349 : /* GetFileList() */
350 : /************************************************************************/
351 :
352 2 : char **SENTINEL2Dataset::GetFileList()
353 : {
354 4 : CPLStringList aosList;
355 7 : for (size_t i = 0; i < aosNonJP2Files.size(); i++)
356 5 : aosList.AddString(aosNonJP2Files[i]);
357 2 : char **papszFileList = VRTDataset::GetFileList();
358 5 : for (char **papszIter = papszFileList; papszIter && *papszIter; ++papszIter)
359 3 : aosList.AddString(*papszIter);
360 2 : CSLDestroy(papszFileList);
361 4 : return aosList.StealList();
362 : }
363 :
364 : /************************************************************************/
365 : /* Identify() */
366 : /************************************************************************/
367 :
368 58994 : int SENTINEL2Dataset::Identify(GDALOpenInfo *poOpenInfo)
369 : {
370 58994 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1B:"))
371 36 : return TRUE;
372 58958 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1B_WITH_GEOLOC:"))
373 18 : return TRUE;
374 58940 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1C:"))
375 78 : return TRUE;
376 58862 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1C_TILE:"))
377 26 : return TRUE;
378 58836 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L2A:"))
379 80 : return TRUE;
380 :
381 58756 : const char *pszJustFilename = CPLGetFilename(poOpenInfo->pszFilename);
382 :
383 : // We don't handle direct tile access for L1C SafeCompact products
384 : // We could, but this isn't just done yet.
385 58757 : if (EQUAL(pszJustFilename, "MTD_TL.xml"))
386 0 : return FALSE;
387 :
388 : /* Accept directly .zip as provided by https://scihub.esa.int/
389 : * First we check just by file name as it is faster than looking
390 : * inside to detect content. */
391 58754 : if ((IsS2Prefixed(pszJustFilename, "_MSIL1C_") ||
392 117510 : IsS2Prefixed(pszJustFilename, "_MSIL2A_") ||
393 117507 : IsS2Prefixed(pszJustFilename, "_OPER_PRD_MSI") ||
394 176268 : IsS2Prefixed(pszJustFilename, "_USER_PRD_MSI")) &&
395 58754 : EQUAL(CPLGetExtensionSafe(pszJustFilename).c_str(), "zip"))
396 : {
397 0 : return TRUE;
398 : }
399 :
400 58755 : if (poOpenInfo->nHeaderBytes < 100)
401 54778 : return FALSE;
402 :
403 3977 : const char *pszHeader =
404 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader);
405 :
406 3977 : if (strstr(pszHeader, "<n1:Level-1B_User_Product") != nullptr &&
407 15 : strstr(pszHeader, "User_Product_Level-1B.xsd") != nullptr)
408 15 : return TRUE;
409 :
410 3962 : if (strstr(pszHeader, "<n1:Level-1B_Granule_ID") != nullptr &&
411 8 : strstr(pszHeader, "S2_PDI_Level-1B_Granule_Metadata.xsd") != nullptr)
412 8 : return TRUE;
413 :
414 3954 : if (strstr(pszHeader, "<n1:Level-1C_User_Product") != nullptr &&
415 18 : strstr(pszHeader, "User_Product_Level-1C.xsd") != nullptr)
416 18 : return TRUE;
417 :
418 3936 : if (strstr(pszHeader, "<n1:Level-1C_Tile_ID") != nullptr &&
419 8 : strstr(pszHeader, "S2_PDI_Level-1C_Tile_Metadata.xsd") != nullptr)
420 8 : return TRUE;
421 :
422 3928 : if (strstr(pszHeader, "<n1:Level-2A_User_Product") != nullptr &&
423 12 : strstr(pszHeader, "User_Product_Level-2A") != nullptr)
424 12 : return TRUE;
425 :
426 3916 : if (SENTINEL2isZipped(pszHeader, poOpenInfo->nHeaderBytes))
427 10 : return TRUE;
428 :
429 3903 : return FALSE;
430 : }
431 :
432 : /************************************************************************/
433 : /* SENTINEL2_CPLXMLNodeHolder */
434 : /************************************************************************/
435 :
436 : class SENTINEL2_CPLXMLNodeHolder
437 : {
438 : CPLXMLNode *m_psNode = nullptr;
439 :
440 : CPL_DISALLOW_COPY_ASSIGN(SENTINEL2_CPLXMLNodeHolder)
441 :
442 : public:
443 190 : explicit SENTINEL2_CPLXMLNodeHolder(CPLXMLNode *psNode) : m_psNode(psNode)
444 : {
445 190 : }
446 :
447 190 : ~SENTINEL2_CPLXMLNodeHolder()
448 190 : {
449 190 : if (m_psNode)
450 174 : CPLDestroyXMLNode(m_psNode);
451 190 : }
452 :
453 12 : CPLXMLNode *Release()
454 : {
455 12 : CPLXMLNode *psRet = m_psNode;
456 12 : m_psNode = nullptr;
457 12 : return psRet;
458 : }
459 : };
460 :
461 : /************************************************************************/
462 : /* Open() */
463 : /************************************************************************/
464 :
465 157 : GDALDataset *SENTINEL2Dataset::Open(GDALOpenInfo *poOpenInfo)
466 : {
467 157 : if (!Identify(poOpenInfo))
468 : {
469 0 : return nullptr;
470 : }
471 :
472 157 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1B:"))
473 : {
474 18 : CPLDebug("SENTINEL2", "Using OpenL1BSubdataset");
475 18 : return OpenL1BSubdataset(poOpenInfo);
476 : }
477 :
478 139 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1B_WITH_GEOLOC:"))
479 : {
480 9 : CPLDebug("SENTINEL2", "Using OpenL1BSubdatasetWithGeoloc");
481 9 : return OpenL1BSubdatasetWithGeoloc(poOpenInfo);
482 : }
483 :
484 130 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1C:"))
485 : {
486 39 : CPLDebug("SENTINEL2", "Using OpenL1C_L2ASubdataset");
487 39 : return OpenL1C_L2ASubdataset(poOpenInfo, SENTINEL2_L1C);
488 : }
489 :
490 91 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1C_TILE:"))
491 : {
492 13 : CPLDebug("SENTINEL2", "Using OpenL1CTileSubdataset");
493 13 : return OpenL1CTileSubdataset(poOpenInfo);
494 : }
495 :
496 78 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L2A:"))
497 : {
498 40 : CPLDebug("SENTINEL2", "Using OpenL1C_L2ASubdataset");
499 40 : return OpenL1C_L2ASubdataset(poOpenInfo, SENTINEL2_L2A);
500 : }
501 :
502 38 : const char *pszJustFilename = CPLGetFilename(poOpenInfo->pszFilename);
503 38 : if ((IsS2Prefixed(pszJustFilename, "_OPER_PRD_MSI") ||
504 76 : IsS2Prefixed(pszJustFilename, "_USER_PRD_MSI")) &&
505 38 : EQUAL(CPLGetExtensionSafe(pszJustFilename).c_str(), "zip"))
506 : {
507 0 : const CPLString osBasename(CPLGetBasenameSafe(pszJustFilename));
508 0 : CPLString osFilename(poOpenInfo->pszFilename);
509 0 : CPLString osMTD(osBasename);
510 : // Normally given above constraints, osMTD.size() should be >= 16
511 : // but if pszJustFilename is too long, CPLGetBasenameSafe().c_str() will return
512 : // an empty string.
513 0 : if (osMTD.size() < 16)
514 0 : return nullptr;
515 0 : osMTD[9] = 'M';
516 0 : osMTD[10] = 'T';
517 0 : osMTD[11] = 'D';
518 0 : osMTD[13] = 'S';
519 0 : osMTD[14] = 'A';
520 0 : osMTD[15] = 'F';
521 0 : CPLString osSAFE(CPLString(osBasename) + ".SAFE");
522 0 : CPL_IGNORE_RET_VAL(osBasename);
523 0 : osFilename = osFilename + "/" + osSAFE + "/" + osMTD + ".xml";
524 0 : if (strncmp(osFilename, "/vsizip/", strlen("/vsizip/")) != 0)
525 0 : osFilename = "/vsizip/" + osFilename;
526 0 : CPLDebug("SENTINEL2", "Trying %s", osFilename.c_str());
527 0 : GDALOpenInfo oOpenInfo(osFilename, GA_ReadOnly);
528 0 : return Open(&oOpenInfo);
529 : }
530 38 : else if (IsS2Prefixed(pszJustFilename, "_MSIL1C_") &&
531 38 : EQUAL(CPLGetExtensionSafe(pszJustFilename).c_str(), "zip"))
532 : {
533 0 : CPLString osFilename(poOpenInfo->pszFilename);
534 0 : CPLString osSAFE(CPLGetBasenameSafe(pszJustFilename));
535 : // S2B_MSIL1C_20171004T233419_N0206_R001_T54DWM_20171005T001811.SAFE.zip
536 : // has .SAFE.zip extension, but other products have just a .zip
537 : // extension. So for the subdir in the zip only add .SAFE when needed
538 0 : if (!EQUAL(CPLGetExtensionSafe(osSAFE).c_str(), "SAFE"))
539 0 : osSAFE += ".SAFE";
540 0 : osFilename = osFilename + "/" + osSAFE + "/MTD_MSIL1C.xml";
541 0 : if (strncmp(osFilename, "/vsizip/", strlen("/vsizip/")) != 0)
542 0 : osFilename = "/vsizip/" + osFilename;
543 0 : CPLDebug("SENTINEL2", "Trying %s", osFilename.c_str());
544 0 : GDALOpenInfo oOpenInfo(osFilename, GA_ReadOnly);
545 0 : return Open(&oOpenInfo);
546 : }
547 38 : else if (IsS2Prefixed(pszJustFilename, "_MSIL2A_") &&
548 38 : EQUAL(CPLGetExtensionSafe(pszJustFilename).c_str(), "zip"))
549 : {
550 0 : CPLString osFilename(poOpenInfo->pszFilename);
551 0 : CPLString osSAFE(CPLGetBasenameSafe(pszJustFilename));
552 : // S2B_MSIL1C_20171004T233419_N0206_R001_T54DWM_20171005T001811.SAFE.zip
553 : // has .SAFE.zip extension, but other products have just a .zip
554 : // extension. So for the subdir in the zip only add .SAFE when needed
555 0 : if (!EQUAL(CPLGetExtensionSafe(osSAFE).c_str(), "SAFE"))
556 0 : osSAFE += ".SAFE";
557 0 : osFilename = osFilename + "/" + osSAFE + "/MTD_MSIL2A.xml";
558 0 : if (strncmp(osFilename, "/vsizip/", strlen("/vsizip/")) != 0)
559 0 : osFilename = "/vsizip/" + osFilename;
560 0 : CPLDebug("SENTINEL2", "Trying %s", osFilename.c_str());
561 0 : GDALOpenInfo oOpenInfo(osFilename, GA_ReadOnly);
562 0 : return Open(&oOpenInfo);
563 : }
564 :
565 38 : const char *pszHeader =
566 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader);
567 :
568 38 : if (strstr(pszHeader, "<n1:Level-1B_User_Product") != nullptr &&
569 8 : strstr(pszHeader, "User_Product_Level-1B.xsd") != nullptr)
570 : {
571 8 : CPLDebug("SENTINEL2", "Trying OpenL1BUserProduct");
572 8 : return OpenL1BUserProduct(poOpenInfo);
573 : }
574 :
575 30 : if (strstr(pszHeader, "<n1:Level-1B_Granule_ID") != nullptr &&
576 4 : strstr(pszHeader, "S2_PDI_Level-1B_Granule_Metadata.xsd") != nullptr)
577 : {
578 4 : CPLDebug("SENTINEL2", "Trying OpenL1BGranule");
579 4 : return OpenL1BGranule(poOpenInfo->pszFilename);
580 : }
581 :
582 26 : if (strstr(pszHeader, "<n1:Level-1C_User_Product") != nullptr &&
583 10 : strstr(pszHeader, "User_Product_Level-1C.xsd") != nullptr)
584 : {
585 10 : CPLDebug("SENTINEL2", "Trying OpenL1C_L2A");
586 10 : return OpenL1C_L2A(poOpenInfo->pszFilename, SENTINEL2_L1C);
587 : }
588 :
589 16 : if (strstr(pszHeader, "<n1:Level-1C_Tile_ID") != nullptr &&
590 4 : strstr(pszHeader, "S2_PDI_Level-1C_Tile_Metadata.xsd") != nullptr)
591 : {
592 4 : CPLDebug("SENTINEL2", "Trying OpenL1CTile");
593 4 : return OpenL1CTile(poOpenInfo->pszFilename);
594 : }
595 :
596 12 : if (strstr(pszHeader, "<n1:Level-2A_User_Product") != nullptr &&
597 7 : strstr(pszHeader, "User_Product_Level-2A") != nullptr)
598 : {
599 7 : CPLDebug("SENTINEL2", "Trying OpenL1C_L2A");
600 7 : return OpenL1C_L2A(poOpenInfo->pszFilename, SENTINEL2_L2A);
601 : }
602 :
603 5 : if (SENTINEL2isZipped(pszHeader, poOpenInfo->nHeaderBytes))
604 : {
605 5 : CPLString osFilename(poOpenInfo->pszFilename);
606 5 : if (strncmp(osFilename, "/vsizip/", strlen("/vsizip/")) != 0)
607 5 : osFilename = "/vsizip/" + osFilename;
608 :
609 5 : auto psDir = VSIOpenDir(osFilename.c_str(), 1, nullptr);
610 5 : if (psDir == nullptr)
611 : {
612 0 : CPLError(CE_Failure, CPLE_AppDefined,
613 : "SENTINEL2: Cannot open ZIP file %s", osFilename.c_str());
614 0 : return nullptr;
615 : }
616 10 : while (const VSIDIREntry *psEntry = VSIGetNextDirEntry(psDir))
617 : {
618 10 : const char *pszInsideFilename = CPLGetFilename(psEntry->pszName);
619 15 : if (VSI_ISREG(psEntry->nMode) &&
620 5 : (STARTS_WITH_CI(pszInsideFilename, "MTD_MSIL2A") ||
621 7 : STARTS_WITH_CI(pszInsideFilename, "MTD_MSIL1C") ||
622 5 : IsS2Prefixed(pszInsideFilename, "_OPER_MTD_SAFL1B") ||
623 3 : IsS2Prefixed(pszInsideFilename, "_OPER_MTD_SAFL1C") ||
624 1 : IsS2Prefixed(pszInsideFilename, "_USER_MTD_SAFL2A") ||
625 0 : IsS2Prefixed(pszInsideFilename, "_USER_MTD_SAFL2A")))
626 : {
627 5 : osFilename = osFilename + "/" + psEntry->pszName;
628 5 : CPLDebug("SENTINEL2", "Trying %s", osFilename.c_str());
629 10 : GDALOpenInfo oOpenInfo(osFilename, GA_ReadOnly);
630 5 : VSICloseDir(psDir);
631 5 : return Open(&oOpenInfo);
632 : }
633 5 : }
634 0 : VSICloseDir(psDir);
635 : }
636 :
637 0 : return nullptr;
638 : }
639 :
640 : /************************************************************************/
641 : /* SENTINEL2isZipped() */
642 : /************************************************************************/
643 :
644 3918 : static bool SENTINEL2isZipped(const char *pszHeader, int nHeaderBytes)
645 : {
646 3918 : if (nHeaderBytes < 50)
647 0 : return FALSE;
648 :
649 : /* According to Sentinel-2 Products Specification Document,
650 : * all files are located inside a folder with a specific name pattern
651 : * Ref: S2-PDGS-TAS-DI-PSD Issue: 14.6.
652 : */
653 : return (
654 : // inside a ZIP file
655 3980 : memcmp(pszHeader, "\x50\x4b", 2) == 0 &&
656 : (
657 : // a "4.2.1 Compact Naming Convention" confirming file
658 62 : (memcmp(pszHeader + 34, "MSIL2A", 6) == 0 ||
659 59 : memcmp(pszHeader + 34, "MSIL1C", 6) == 0) ||
660 : // a "4.2 S2 User Product Naming Convention" confirming file
661 56 : (memcmp(pszHeader + 34, "OPER_PRD_MSIL2A", 15) == 0 ||
662 56 : memcmp(pszHeader + 34, "OPER_PRD_MSIL1B", 15) == 0 ||
663 53 : memcmp(pszHeader + 34, "OPER_PRD_MSIL1C", 15) == 0) ||
664 : // some old / validation naming convention
665 50 : (memcmp(pszHeader + 34, "USER_PRD_MSIL2A", 15) == 0 ||
666 47 : memcmp(pszHeader + 34, "USER_PRD_MSIL1B", 15) == 0 ||
667 3965 : memcmp(pszHeader + 34, "USER_PRD_MSIL1C", 15) == 0)));
668 : }
669 :
670 : /************************************************************************/
671 : /* SENTINEL2GetBandDesc() */
672 : /************************************************************************/
673 :
674 : static const SENTINEL2BandDescription *
675 598 : SENTINEL2GetBandDesc(const char *pszBandName)
676 : {
677 4282 : for (const auto &sBandDesc : asBandDesc)
678 : {
679 4238 : if (EQUAL(sBandDesc.pszBandName, pszBandName))
680 554 : return &sBandDesc;
681 : }
682 44 : return nullptr;
683 : }
684 :
685 : /************************************************************************/
686 : /* SENTINEL2GetL2ABandDesc() */
687 : /************************************************************************/
688 :
689 : static const SENTINEL2_L2A_BandDescription *
690 191 : SENTINEL2GetL2ABandDesc(const char *pszBandName)
691 : {
692 1721 : for (const auto &sBandDesc : asL2ABandDesc)
693 : {
694 1636 : if (EQUAL(sBandDesc.pszBandName, pszBandName))
695 106 : return &sBandDesc;
696 : }
697 85 : return nullptr;
698 : }
699 :
700 : /************************************************************************/
701 : /* SENTINEL2GetGranuleInfo() */
702 : /************************************************************************/
703 :
704 84 : static bool SENTINEL2GetGranuleInfo(
705 : SENTINEL2Level eLevel, const CPLString &osGranuleMTDPath,
706 : int nDesiredResolution, int *pnEPSGCode = nullptr, double *pdfULX = nullptr,
707 : double *pdfULY = nullptr, int *pnResolution = nullptr,
708 : int *pnWidth = nullptr, int *pnHeight = nullptr)
709 : {
710 : static bool bTryOptimization = true;
711 84 : CPLXMLNode *psRoot = nullptr;
712 :
713 84 : if (bTryOptimization)
714 : {
715 : /* Small optimization: in practice the interesting info are in the */
716 : /* first bytes of the Granule MTD, which can be very long sometimes */
717 : /* so only read them, and hack the buffer a bit to form a valid XML */
718 : char szBuffer[3072];
719 84 : VSILFILE *fp = VSIFOpenL(osGranuleMTDPath, "rb");
720 84 : size_t nRead = 0;
721 162 : if (fp == nullptr ||
722 78 : (nRead = VSIFReadL(szBuffer, 1, sizeof(szBuffer) - 1, fp)) == 0)
723 : {
724 6 : if (fp)
725 0 : VSIFCloseL(fp);
726 6 : CPLError(CE_Failure, CPLE_AppDefined,
727 : "SENTINEL2GetGranuleInfo: Cannot read %s",
728 : osGranuleMTDPath.c_str());
729 6 : return false;
730 : }
731 78 : szBuffer[nRead] = 0;
732 78 : VSIFCloseL(fp);
733 78 : char *pszTileGeocoding = strstr(szBuffer, "</Tile_Geocoding>");
734 78 : if (eLevel == SENTINEL2_L1C && pszTileGeocoding != nullptr &&
735 44 : strstr(szBuffer, "<n1:Level-1C_Tile_ID") != nullptr &&
736 44 : strstr(szBuffer, "<n1:Geometric_Info") != nullptr &&
737 44 : static_cast<size_t>(pszTileGeocoding - szBuffer) <
738 : sizeof(szBuffer) -
739 : strlen("</Tile_Geocoding></n1:Geometric_Info></"
740 : "n1:Level-1C_Tile_ID>") -
741 : 1)
742 : {
743 44 : strcpy(
744 : pszTileGeocoding,
745 : "</Tile_Geocoding></n1:Geometric_Info></n1:Level-1C_Tile_ID>");
746 44 : psRoot = CPLParseXMLString(szBuffer);
747 : }
748 34 : else if (eLevel == SENTINEL2_L2A && pszTileGeocoding != nullptr &&
749 34 : strstr(szBuffer, "<n1:Level-2A_Tile_ID") != nullptr &&
750 34 : strstr(szBuffer, "<n1:Geometric_Info") != nullptr &&
751 34 : static_cast<size_t>(pszTileGeocoding - szBuffer) <
752 : sizeof(szBuffer) -
753 : strlen("</Tile_Geocoding></n1:Geometric_Info></"
754 : "n1:Level-2A_Tile_ID>") -
755 : 1)
756 : {
757 34 : strcpy(
758 : pszTileGeocoding,
759 : "</Tile_Geocoding></n1:Geometric_Info></n1:Level-2A_Tile_ID>");
760 34 : psRoot = CPLParseXMLString(szBuffer);
761 : }
762 : else
763 0 : bTryOptimization = false;
764 : }
765 :
766 : // If the above doesn't work, then read the whole file...
767 78 : if (psRoot == nullptr)
768 0 : psRoot = CPLParseXMLFile(osGranuleMTDPath);
769 78 : if (psRoot == nullptr)
770 : {
771 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot XML parse %s",
772 : osGranuleMTDPath.c_str());
773 0 : return false;
774 : }
775 156 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRoot);
776 78 : CPLStripXMLNamespace(psRoot, nullptr, TRUE);
777 :
778 78 : const char *pszNodePath =
779 : (eLevel == SENTINEL2_L1C)
780 : ? "=Level-1C_Tile_ID.Geometric_Info.Tile_Geocoding"
781 : : "=Level-2A_Tile_ID.Geometric_Info.Tile_Geocoding";
782 78 : CPLXMLNode *psTileGeocoding = CPLGetXMLNode(psRoot, pszNodePath);
783 78 : if (psTileGeocoding == nullptr)
784 : {
785 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s in %s",
786 : pszNodePath, osGranuleMTDPath.c_str());
787 0 : return false;
788 : }
789 :
790 : const char *pszCSCode =
791 78 : CPLGetXMLValue(psTileGeocoding, "HORIZONTAL_CS_CODE", nullptr);
792 78 : if (pszCSCode == nullptr)
793 : {
794 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s in %s",
795 : "HORIZONTAL_CS_CODE", osGranuleMTDPath.c_str());
796 0 : return false;
797 : }
798 78 : if (!STARTS_WITH_CI(pszCSCode, "EPSG:"))
799 : {
800 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid CS code (%s) for %s",
801 : pszCSCode, osGranuleMTDPath.c_str());
802 0 : return false;
803 : }
804 78 : int nEPSGCode = atoi(pszCSCode + strlen("EPSG:"));
805 78 : if (pnEPSGCode != nullptr)
806 78 : *pnEPSGCode = nEPSGCode;
807 :
808 768 : for (CPLXMLNode *psIter = psTileGeocoding->psChild; psIter != nullptr;
809 690 : psIter = psIter->psNext)
810 : {
811 690 : if (psIter->eType != CXT_Element)
812 78 : continue;
813 831 : if (EQUAL(psIter->pszValue, "Size") &&
814 219 : (nDesiredResolution == 0 ||
815 219 : atoi(CPLGetXMLValue(psIter, "resolution", "")) ==
816 : nDesiredResolution))
817 : {
818 77 : nDesiredResolution = atoi(CPLGetXMLValue(psIter, "resolution", ""));
819 77 : const char *pszRows = CPLGetXMLValue(psIter, "NROWS", nullptr);
820 77 : if (pszRows == nullptr)
821 : {
822 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s in %s",
823 : "NROWS", osGranuleMTDPath.c_str());
824 0 : return false;
825 : }
826 77 : const char *pszCols = CPLGetXMLValue(psIter, "NCOLS", nullptr);
827 77 : if (pszCols == nullptr)
828 : {
829 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s in %s",
830 : "NCOLS", osGranuleMTDPath.c_str());
831 0 : return false;
832 : }
833 77 : if (pnResolution)
834 64 : *pnResolution = nDesiredResolution;
835 77 : if (pnWidth)
836 64 : *pnWidth = atoi(pszCols);
837 77 : if (pnHeight)
838 64 : *pnHeight = atoi(pszRows);
839 : }
840 763 : else if (EQUAL(psIter->pszValue, "Geoposition") &&
841 228 : (nDesiredResolution == 0 ||
842 228 : atoi(CPLGetXMLValue(psIter, "resolution", "")) ==
843 : nDesiredResolution))
844 : {
845 77 : nDesiredResolution = atoi(CPLGetXMLValue(psIter, "resolution", ""));
846 77 : const char *pszULX = CPLGetXMLValue(psIter, "ULX", nullptr);
847 77 : if (pszULX == nullptr)
848 : {
849 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s in %s",
850 : "ULX", osGranuleMTDPath.c_str());
851 0 : return false;
852 : }
853 77 : const char *pszULY = CPLGetXMLValue(psIter, "ULY", nullptr);
854 77 : if (pszULY == nullptr)
855 : {
856 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s in %s",
857 : "ULY", osGranuleMTDPath.c_str());
858 0 : return false;
859 : }
860 77 : if (pnResolution)
861 64 : *pnResolution = nDesiredResolution;
862 77 : if (pdfULX)
863 64 : *pdfULX = CPLAtof(pszULX);
864 77 : if (pdfULY)
865 64 : *pdfULY = CPLAtof(pszULY);
866 : }
867 : }
868 :
869 78 : return true;
870 : }
871 :
872 : /************************************************************************/
873 : /* SENTINEL2GetPathSeparator() */
874 : /************************************************************************/
875 :
876 : // For the sake of simplifying our unit tests, we limit the use of \\ to when
877 : // it is strictly necessary. Otherwise we could use CPLFormFilenameSafe()...
878 416 : static char SENTINEL2GetPathSeparator(const char *pszBasename)
879 : {
880 416 : if (STARTS_WITH_CI(pszBasename, "\\\\?\\"))
881 0 : return '\\';
882 : else
883 416 : return '/';
884 : }
885 :
886 : /************************************************************************/
887 : /* SENTINEL2GetGranuleList() */
888 : /************************************************************************/
889 :
890 27 : static bool SENTINEL2GetGranuleList(
891 : CPLXMLNode *psMainMTD, SENTINEL2Level eLevel, const char *pszFilename,
892 : std::vector<CPLString> &osList, std::set<int> *poSetResolutions = nullptr,
893 : std::map<int, std::set<CPLString>> *poMapResolutionsToBands = nullptr)
894 : {
895 27 : const char *pszNodePath =
896 : (eLevel == SENTINEL2_L1B) ? "Level-1B_User_Product"
897 : : (eLevel == SENTINEL2_L1C) ? "Level-1C_User_Product"
898 : : "Level-2A_User_Product";
899 :
900 : CPLXMLNode *psRoot =
901 27 : CPLGetXMLNode(psMainMTD, CPLSPrintf("=%s", pszNodePath));
902 27 : if (psRoot == nullptr)
903 : {
904 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find =%s", pszNodePath);
905 0 : return false;
906 : }
907 27 : pszNodePath = "General_Info.Product_Info";
908 27 : CPLXMLNode *psProductInfo = CPLGetXMLNode(psRoot, pszNodePath);
909 27 : if (psProductInfo == nullptr && eLevel == SENTINEL2_L2A)
910 : {
911 9 : pszNodePath = "General_Info.L2A_Product_Info";
912 9 : psProductInfo = CPLGetXMLNode(psRoot, pszNodePath);
913 : }
914 27 : if (psProductInfo == nullptr)
915 : {
916 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s", pszNodePath);
917 0 : return false;
918 : }
919 :
920 27 : pszNodePath = "Product_Organisation";
921 : CPLXMLNode *psProductOrganisation =
922 27 : CPLGetXMLNode(psProductInfo, pszNodePath);
923 27 : if (psProductOrganisation == nullptr && eLevel == SENTINEL2_L2A)
924 : {
925 9 : pszNodePath = "L2A_Product_Organisation";
926 9 : psProductOrganisation = CPLGetXMLNode(psProductInfo, pszNodePath);
927 : }
928 27 : if (psProductOrganisation == nullptr)
929 : {
930 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s", pszNodePath);
931 1 : return false;
932 : }
933 :
934 52 : CPLString osDirname(CPLGetDirnameSafe(pszFilename));
935 : #if !defined(_WIN32)
936 : char szPointerFilename[2048];
937 : int nBytes = static_cast<int>(
938 26 : readlink(pszFilename, szPointerFilename, sizeof(szPointerFilename)));
939 26 : if (nBytes != -1)
940 : {
941 : const int nOffset =
942 0 : std::min(nBytes, static_cast<int>(sizeof(szPointerFilename) - 1));
943 0 : szPointerFilename[nOffset] = '\0';
944 0 : osDirname = CPLGetDirnameSafe(szPointerFilename);
945 : }
946 : #endif
947 :
948 : const bool bIsMSI2Ap =
949 26 : EQUAL(CPLGetXMLValue(psProductInfo, "PRODUCT_TYPE", ""), "S2MSI2Ap");
950 26 : const bool bIsCompact = EQUAL(
951 : CPLGetXMLValue(psProductInfo, "PRODUCT_FORMAT", ""), "SAFE_COMPACT");
952 52 : CPLString oGranuleId("L2A_");
953 52 : std::set<CPLString> aoSetGranuleId;
954 73 : for (CPLXMLNode *psIter = psProductOrganisation->psChild; psIter != nullptr;
955 47 : psIter = psIter->psNext)
956 : {
957 47 : if (psIter->eType != CXT_Element ||
958 47 : !EQUAL(psIter->pszValue, "Granule_List"))
959 : {
960 0 : continue;
961 : }
962 102 : for (CPLXMLNode *psIter2 = psIter->psChild; psIter2 != nullptr;
963 55 : psIter2 = psIter2->psNext)
964 : {
965 55 : if (psIter2->eType != CXT_Element ||
966 47 : (!EQUAL(psIter2->pszValue, "Granule") &&
967 47 : !EQUAL(psIter2->pszValue, "Granules")))
968 : {
969 16 : continue;
970 : }
971 : const char *pszGranuleId =
972 47 : CPLGetXMLValue(psIter2, "granuleIdentifier", nullptr);
973 47 : if (pszGranuleId == nullptr)
974 : {
975 4 : CPLDebug("SENTINEL2", "Missing granuleIdentifier attribute");
976 4 : continue;
977 : }
978 :
979 43 : if (eLevel == SENTINEL2_L2A)
980 : {
981 161 : for (CPLXMLNode *psIter3 = psIter2->psChild; psIter3 != nullptr;
982 150 : psIter3 = psIter3->psNext)
983 : {
984 150 : if (psIter3->eType != CXT_Element ||
985 121 : (!EQUAL(psIter3->pszValue, "IMAGE_ID_2A") &&
986 0 : !EQUAL(psIter3->pszValue, "IMAGE_FILE") &&
987 0 : !EQUAL(psIter3->pszValue, "IMAGE_FILE_2A")))
988 : {
989 29 : continue;
990 : }
991 : const char *pszTileName =
992 121 : CPLGetXMLValue(psIter3, nullptr, "");
993 121 : size_t nLen = strlen(pszTileName);
994 : // If granule name ends with resolution: _60m
995 121 : if (nLen > 4 && pszTileName[nLen - 4] == '_' &&
996 121 : pszTileName[nLen - 1] == 'm')
997 : {
998 121 : int nResolution = atoi(pszTileName + nLen - 3);
999 121 : if (poSetResolutions != nullptr)
1000 35 : (*poSetResolutions).insert(nResolution);
1001 121 : if (poMapResolutionsToBands != nullptr)
1002 : {
1003 121 : nLen -= 4;
1004 121 : if (nLen > 4 && pszTileName[nLen - 4] == '_' &&
1005 86 : pszTileName[nLen - 3] == 'B')
1006 : {
1007 86 : (*poMapResolutionsToBands)[nResolution].insert(
1008 86 : CPLString(pszTileName).substr(nLen - 2, 2));
1009 : }
1010 35 : else if (nLen > strlen("S2A_USER_MSI_") &&
1011 35 : pszTileName[8] == '_' &&
1012 35 : pszTileName[12] == '_' &&
1013 35 : !EQUALN(pszTileName + 9, "MSI", 3))
1014 : {
1015 35 : (*poMapResolutionsToBands)[nResolution].insert(
1016 35 : CPLString(pszTileName).substr(9, 3));
1017 : }
1018 : }
1019 : }
1020 : }
1021 : }
1022 :
1023 : // For L2A we can have several time the same granuleIdentifier
1024 : // for the different resolutions
1025 43 : if (aoSetGranuleId.find(pszGranuleId) != aoSetGranuleId.end())
1026 0 : continue;
1027 43 : aoSetGranuleId.insert(pszGranuleId);
1028 :
1029 : /* S2A_OPER_MSI_L1C_TL_SGS__20151024T023555_A001758_T53JLJ_N01.04
1030 : * --> */
1031 : /* S2A_OPER_MTD_L1C_TL_SGS__20151024T023555_A001758_T53JLJ */
1032 : // S2B_OPER_MSI_L2A_TL_MPS__20180823T122014_A007641_T34VFJ_N02.08
1033 43 : CPLString osGranuleMTD = pszGranuleId;
1034 43 : if (bIsCompact == 0 &&
1035 43 : osGranuleMTD.size() > strlen("S2A_OPER_MSI_") &&
1036 39 : osGranuleMTD[8] == '_' && osGranuleMTD[12] == '_' &&
1037 39 : osGranuleMTD[osGranuleMTD.size() - 7] == '_' &&
1038 125 : osGranuleMTD[osGranuleMTD.size() - 6] == 'N' &&
1039 39 : osGranuleMTD[7] == 'R')
1040 : {
1041 39 : osGranuleMTD[9] = 'M';
1042 39 : osGranuleMTD[10] = 'T';
1043 39 : osGranuleMTD[11] = 'D';
1044 39 : osGranuleMTD.resize(osGranuleMTD.size() - 7);
1045 : }
1046 4 : else if (bIsMSI2Ap)
1047 : {
1048 0 : osGranuleMTD = "MTD_TL";
1049 0 : oGranuleId = "L2A_";
1050 : // S2A_MSIL2A_20170823T094031_N0205_R036_T34VFJ_20170823T094252.SAFE
1051 : // S2A_USER_MSI_L2A_TL_SGS__20170823T133142_A011330_T34VFJ_N02.05
1052 : // --> L2A_T34VFJ_A011330_20170823T094252
1053 : const char *pszProductURI =
1054 0 : CPLGetXMLValue(psProductInfo, "PRODUCT_URI_2A", nullptr);
1055 0 : if (pszProductURI != nullptr)
1056 : {
1057 0 : CPLString psProductURI(pszProductURI);
1058 0 : if (psProductURI.size() < 60)
1059 : {
1060 0 : CPLDebug("SENTINEL2", "Invalid PRODUCT_URI_2A");
1061 0 : continue;
1062 : }
1063 0 : oGranuleId += psProductURI.substr(38, 7);
1064 0 : oGranuleId += CPLString(pszGranuleId).substr(41, 8).c_str();
1065 0 : oGranuleId += psProductURI.substr(45, 15);
1066 0 : pszGranuleId = oGranuleId.c_str();
1067 : }
1068 : }
1069 : else
1070 : {
1071 4 : CPLDebug("SENTINEL2", "Invalid granule ID: %s", pszGranuleId);
1072 4 : continue;
1073 : }
1074 39 : osGranuleMTD += ".xml";
1075 :
1076 39 : const char chSeparator = SENTINEL2GetPathSeparator(osDirname);
1077 39 : CPLString osGranuleMTDPath = osDirname;
1078 39 : osGranuleMTDPath += chSeparator;
1079 39 : osGranuleMTDPath += "GRANULE";
1080 39 : osGranuleMTDPath += chSeparator;
1081 39 : osGranuleMTDPath += pszGranuleId;
1082 39 : osGranuleMTDPath += chSeparator;
1083 39 : osGranuleMTDPath += osGranuleMTD;
1084 39 : if (CPLHasPathTraversal(osGranuleMTDPath.c_str()))
1085 : {
1086 0 : CPLError(CE_Failure, CPLE_NotSupported,
1087 : "Path traversal detected in %s",
1088 : osGranuleMTDPath.c_str());
1089 0 : return false;
1090 : }
1091 39 : osList.push_back(std::move(osGranuleMTDPath));
1092 : }
1093 : }
1094 :
1095 26 : return true;
1096 : }
1097 :
1098 : /************************************************************************/
1099 : /* SENTINEL2GetUserProductMetadata() */
1100 : /************************************************************************/
1101 :
1102 74 : static char **SENTINEL2GetUserProductMetadata(CPLXMLNode *psMainMTD,
1103 : const char *pszRootNode)
1104 : {
1105 148 : CPLStringList aosList;
1106 :
1107 : CPLXMLNode *psRoot =
1108 74 : CPLGetXMLNode(psMainMTD, CPLSPrintf("=%s", pszRootNode));
1109 74 : if (psRoot == nullptr)
1110 : {
1111 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find =%s", pszRootNode);
1112 0 : return nullptr;
1113 : }
1114 74 : const char *psPIPath = "General_Info.Product_Info";
1115 74 : CPLXMLNode *psProductInfo = CPLGetXMLNode(psRoot, psPIPath);
1116 74 : if (psProductInfo == nullptr && EQUAL(pszRootNode, "Level-2A_User_Product"))
1117 : {
1118 15 : psPIPath = "General_Info.L2A_Product_Info";
1119 15 : psProductInfo = CPLGetXMLNode(psRoot, psPIPath);
1120 : }
1121 74 : if (psProductInfo == nullptr)
1122 : {
1123 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find =%s", psPIPath);
1124 0 : return nullptr;
1125 : }
1126 74 : int nDataTakeCounter = 1;
1127 869 : for (CPLXMLNode *psIter = psProductInfo->psChild; psIter != nullptr;
1128 795 : psIter = psIter->psNext)
1129 : {
1130 795 : if (psIter->eType != CXT_Element)
1131 0 : continue;
1132 795 : if (psIter->psChild != nullptr && psIter->psChild->eType == CXT_Text)
1133 : {
1134 539 : aosList.AddNameValue(psIter->pszValue, psIter->psChild->pszValue);
1135 : }
1136 256 : else if (EQUAL(psIter->pszValue, "Datatake"))
1137 : {
1138 128 : CPLString osPrefix(CPLSPrintf("DATATAKE_%d_", nDataTakeCounter));
1139 64 : nDataTakeCounter++;
1140 : const char *pszId =
1141 64 : CPLGetXMLValue(psIter, "datatakeIdentifier", nullptr);
1142 64 : if (pszId)
1143 64 : aosList.AddNameValue((osPrefix + "ID").c_str(), pszId);
1144 448 : for (CPLXMLNode *psIter2 = psIter->psChild; psIter2 != nullptr;
1145 384 : psIter2 = psIter2->psNext)
1146 : {
1147 384 : if (psIter2->eType != CXT_Element)
1148 64 : continue;
1149 320 : if (psIter2->psChild != nullptr &&
1150 320 : psIter2->psChild->eType == CXT_Text)
1151 : {
1152 320 : aosList.AddNameValue((osPrefix + psIter2->pszValue).c_str(),
1153 320 : psIter2->psChild->pszValue);
1154 : }
1155 : }
1156 : }
1157 : }
1158 :
1159 74 : const char *psICPath = "General_Info.Product_Image_Characteristics";
1160 74 : CPLXMLNode *psIC = CPLGetXMLNode(psRoot, psICPath);
1161 74 : if (psIC == nullptr)
1162 : {
1163 23 : psICPath = "General_Info.L2A_Product_Image_Characteristics";
1164 23 : psIC = CPLGetXMLNode(psRoot, psICPath);
1165 : }
1166 74 : if (psIC != nullptr)
1167 : {
1168 1091 : for (CPLXMLNode *psIter = psIC->psChild; psIter != nullptr;
1169 1027 : psIter = psIter->psNext)
1170 : {
1171 1027 : if (psIter->eType != CXT_Element ||
1172 1020 : !EQUAL(psIter->pszValue, "Special_Values"))
1173 : {
1174 899 : continue;
1175 : }
1176 : const char *pszText =
1177 128 : CPLGetXMLValue(psIter, "SPECIAL_VALUE_TEXT", nullptr);
1178 : const char *pszIndex =
1179 128 : CPLGetXMLValue(psIter, "SPECIAL_VALUE_INDEX", nullptr);
1180 128 : if (pszText && pszIndex)
1181 : {
1182 : aosList.AddNameValue(
1183 128 : (CPLString("SPECIAL_VALUE_") + pszText).c_str(), pszIndex);
1184 : }
1185 : }
1186 :
1187 : const char *pszQuantValue =
1188 64 : CPLGetXMLValue(psIC, "QUANTIFICATION_VALUE", nullptr);
1189 64 : if (pszQuantValue != nullptr)
1190 30 : aosList.AddNameValue("QUANTIFICATION_VALUE", pszQuantValue);
1191 :
1192 : const char *pszRCU =
1193 64 : CPLGetXMLValue(psIC, "Reflectance_Conversion.U", nullptr);
1194 64 : if (pszRCU != nullptr)
1195 52 : aosList.AddNameValue("REFLECTANCE_CONVERSION_U", pszRCU);
1196 :
1197 : // L2A specific
1198 : CPLXMLNode *psQVL =
1199 64 : CPLGetXMLNode(psIC, "L1C_L2A_Quantification_Values_List");
1200 64 : if (psQVL == nullptr)
1201 : {
1202 51 : psQVL = CPLGetXMLNode(psIC, "Quantification_Values_List");
1203 : }
1204 64 : for (CPLXMLNode *psIter = psQVL ? psQVL->psChild : nullptr;
1205 143 : psIter != nullptr; psIter = psIter->psNext)
1206 : {
1207 79 : if (psIter->eType != CXT_Element)
1208 : {
1209 0 : continue;
1210 : }
1211 79 : aosList.AddNameValue(psIter->pszValue,
1212 79 : CPLGetXMLValue(psIter, nullptr, nullptr));
1213 79 : const char *pszUnit = CPLGetXMLValue(psIter, "unit", nullptr);
1214 79 : if (pszUnit)
1215 : aosList.AddNameValue(CPLSPrintf("%s_UNIT", psIter->pszValue),
1216 79 : pszUnit);
1217 : }
1218 :
1219 : const char *pszRefBand =
1220 64 : CPLGetXMLValue(psIC, "REFERENCE_BAND", nullptr);
1221 64 : if (pszRefBand != nullptr)
1222 : {
1223 41 : const unsigned int nIdx = atoi(pszRefBand);
1224 41 : if (nIdx < CPL_ARRAYSIZE(asBandDesc))
1225 : aosList.AddNameValue("REFERENCE_BAND",
1226 41 : asBandDesc[nIdx].pszBandName);
1227 : }
1228 : }
1229 :
1230 74 : CPLXMLNode *psQII = CPLGetXMLNode(psRoot, "Quality_Indicators_Info");
1231 74 : if (psQII != nullptr)
1232 : {
1233 : const char *pszCC =
1234 64 : CPLGetXMLValue(psQII, "Cloud_Coverage_Assessment", nullptr);
1235 64 : if (pszCC)
1236 64 : aosList.AddNameValue("CLOUD_COVERAGE_ASSESSMENT", pszCC);
1237 :
1238 64 : const char *pszDegradedAnc = CPLGetXMLValue(
1239 : psQII, "Technical_Quality_Assessment.DEGRADED_ANC_DATA_PERCENTAGE",
1240 : nullptr);
1241 64 : if (pszDegradedAnc)
1242 : aosList.AddNameValue("DEGRADED_ANC_DATA_PERCENTAGE",
1243 64 : pszDegradedAnc);
1244 :
1245 64 : const char *pszDegradedMSI = CPLGetXMLValue(
1246 : psQII, "Technical_Quality_Assessment.DEGRADED_MSI_DATA_PERCENTAGE",
1247 : nullptr);
1248 64 : if (pszDegradedMSI)
1249 : aosList.AddNameValue("DEGRADED_MSI_DATA_PERCENTAGE",
1250 64 : pszDegradedMSI);
1251 :
1252 : CPLXMLNode *psQualInspect =
1253 64 : CPLGetXMLNode(psQII, "Quality_Control_Checks.Quality_Inspections");
1254 64 : for (CPLXMLNode *psIter =
1255 64 : (psQualInspect ? psQualInspect->psChild : nullptr);
1256 386 : psIter != nullptr; psIter = psIter->psNext)
1257 : {
1258 : // MSIL2A approach
1259 322 : if (psIter->psChild != nullptr &&
1260 322 : psIter->psChild->psChild != nullptr &&
1261 57 : psIter->psChild->psNext != nullptr &&
1262 57 : psIter->psChild->psChild->eType == CXT_Text &&
1263 57 : psIter->psChild->psNext->eType == CXT_Text)
1264 : {
1265 57 : aosList.AddNameValue(psIter->psChild->psChild->pszValue,
1266 57 : psIter->psChild->psNext->pszValue);
1267 57 : continue;
1268 : }
1269 :
1270 265 : if (psIter->eType != CXT_Element)
1271 0 : continue;
1272 265 : if (psIter->psChild != nullptr &&
1273 265 : psIter->psChild->eType == CXT_Text)
1274 : {
1275 265 : aosList.AddNameValue(psIter->pszValue,
1276 265 : psIter->psChild->pszValue);
1277 : }
1278 : }
1279 :
1280 64 : CPLXMLNode *psICCQI = CPLGetXMLNode(psQII, "Image_Content_QI");
1281 64 : if (psICCQI == nullptr)
1282 : {
1283 : CPLXMLNode *psL2A_QII =
1284 55 : CPLGetXMLNode(psRoot, "L2A_Quality_Indicators_Info");
1285 55 : if (psL2A_QII != nullptr)
1286 : {
1287 13 : psICCQI = CPLGetXMLNode(psL2A_QII, "Image_Content_QI");
1288 : }
1289 : }
1290 64 : if (psICCQI != nullptr)
1291 : {
1292 370 : for (CPLXMLNode *psIter = psICCQI->psChild; psIter != nullptr;
1293 348 : psIter = psIter->psNext)
1294 : {
1295 348 : if (psIter->eType != CXT_Element)
1296 0 : continue;
1297 348 : if (psIter->psChild != nullptr &&
1298 348 : psIter->psChild->eType == CXT_Text)
1299 : {
1300 348 : aosList.AddNameValue(psIter->pszValue,
1301 348 : psIter->psChild->pszValue);
1302 : }
1303 : }
1304 : }
1305 : }
1306 :
1307 74 : return aosList.StealList();
1308 : }
1309 :
1310 : /************************************************************************/
1311 : /* SENTINEL2GetResolutionSet() */
1312 : /************************************************************************/
1313 :
1314 28 : static bool SENTINEL2GetResolutionSet(
1315 : CPLXMLNode *psProductInfo, std::set<int> &oSetResolutions,
1316 : std::map<int, std::set<CPLString>> &oMapResolutionsToBands)
1317 : {
1318 :
1319 : CPLXMLNode *psBandList =
1320 28 : CPLGetXMLNode(psProductInfo, "Query_Options.Band_List");
1321 28 : if (psBandList == nullptr)
1322 : {
1323 4 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
1324 : "Query_Options.Band_List");
1325 4 : return false;
1326 : }
1327 :
1328 300 : for (CPLXMLNode *psIter = psBandList->psChild; psIter != nullptr;
1329 276 : psIter = psIter->psNext)
1330 : {
1331 276 : if (psIter->eType != CXT_Element ||
1332 276 : !EQUAL(psIter->pszValue, "BAND_NAME"))
1333 : {
1334 2 : continue;
1335 : }
1336 276 : const char *pszBandName = CPLGetXMLValue(psIter, nullptr, "");
1337 : const SENTINEL2BandDescription *psBandDesc =
1338 276 : SENTINEL2GetBandDesc(pszBandName);
1339 276 : if (psBandDesc == nullptr)
1340 : {
1341 2 : CPLDebug("SENTINEL2", "Unknown band name %s", pszBandName);
1342 2 : continue;
1343 : }
1344 274 : oSetResolutions.insert(psBandDesc->nResolution);
1345 548 : CPLString osName = psBandDesc->pszBandName + 1; /* skip B character */
1346 274 : if (atoi(osName) < 10)
1347 211 : osName = "0" + osName;
1348 274 : oMapResolutionsToBands[psBandDesc->nResolution].insert(
1349 274 : std::move(osName));
1350 : }
1351 24 : if (oSetResolutions.empty())
1352 : {
1353 2 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find any band");
1354 2 : return false;
1355 : }
1356 22 : return true;
1357 : }
1358 :
1359 : /************************************************************************/
1360 : /* SENTINEL2GetPolygonWKTFromPosList() */
1361 : /************************************************************************/
1362 :
1363 18 : static CPLString SENTINEL2GetPolygonWKTFromPosList(const char *pszPosList)
1364 : {
1365 18 : CPLString osPolygon;
1366 18 : char **papszTokens = CSLTokenizeString(pszPosList);
1367 18 : int nTokens = CSLCount(papszTokens);
1368 18 : int nDim = 2;
1369 18 : if ((nTokens % 3) == 0 && nTokens >= 3 * 4 &&
1370 6 : EQUAL(papszTokens[0], papszTokens[nTokens - 3]) &&
1371 6 : EQUAL(papszTokens[1], papszTokens[nTokens - 2]) &&
1372 6 : EQUAL(papszTokens[2], papszTokens[nTokens - 1]))
1373 : {
1374 6 : nDim = 3;
1375 : }
1376 18 : if ((nTokens % nDim) == 0)
1377 : {
1378 18 : osPolygon = "POLYGON((";
1379 108 : for (char **papszIter = papszTokens; *papszIter; papszIter += nDim)
1380 : {
1381 90 : if (papszIter != papszTokens)
1382 72 : osPolygon += ", ";
1383 90 : osPolygon += papszIter[1];
1384 90 : osPolygon += " ";
1385 90 : osPolygon += papszIter[0];
1386 90 : if (nDim == 3)
1387 : {
1388 30 : osPolygon += " ";
1389 30 : osPolygon += papszIter[2];
1390 : }
1391 : }
1392 18 : osPolygon += "))";
1393 : }
1394 18 : CSLDestroy(papszTokens);
1395 18 : return osPolygon;
1396 : }
1397 :
1398 : /************************************************************************/
1399 : /* SENTINEL2GetBandListForResolution() */
1400 : /************************************************************************/
1401 :
1402 : static CPLString
1403 76 : SENTINEL2GetBandListForResolution(const std::set<CPLString> &oBandnames)
1404 : {
1405 76 : CPLString osBandNames;
1406 374 : for (std::set<CPLString>::const_iterator oIterBandnames =
1407 76 : oBandnames.begin();
1408 824 : oIterBandnames != oBandnames.end(); ++oIterBandnames)
1409 : {
1410 374 : if (!osBandNames.empty())
1411 298 : osBandNames += ", ";
1412 374 : const char *pszName = *oIterBandnames;
1413 374 : if (*pszName == DIGIT_ZERO)
1414 254 : pszName++;
1415 374 : if (atoi(pszName) > 0)
1416 328 : osBandNames += "B" + CPLString(pszName);
1417 : else
1418 46 : osBandNames += pszName;
1419 : }
1420 76 : return osBandNames;
1421 : }
1422 :
1423 : /************************************************************************/
1424 : /* OpenL1BUserProduct() */
1425 : /************************************************************************/
1426 :
1427 8 : GDALDataset *SENTINEL2Dataset::OpenL1BUserProduct(GDALOpenInfo *poOpenInfo)
1428 : {
1429 8 : CPLXMLNode *psRoot = CPLParseXMLFile(poOpenInfo->pszFilename);
1430 8 : if (psRoot == nullptr)
1431 : {
1432 1 : CPLDebug("SENTINEL2", "Cannot XML parse %s", poOpenInfo->pszFilename);
1433 1 : return nullptr;
1434 : }
1435 :
1436 7 : char *pszOriginalXML = CPLSerializeXMLTree(psRoot);
1437 14 : CPLString osOriginalXML;
1438 7 : if (pszOriginalXML)
1439 7 : osOriginalXML = pszOriginalXML;
1440 7 : CPLFree(pszOriginalXML);
1441 :
1442 14 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRoot);
1443 7 : CPLStripXMLNamespace(psRoot, nullptr, TRUE);
1444 :
1445 7 : CPLXMLNode *psProductInfo = CPLGetXMLNode(
1446 : psRoot, "=Level-1B_User_Product.General_Info.Product_Info");
1447 7 : if (psProductInfo == nullptr)
1448 : {
1449 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
1450 : "=Level-1B_User_Product.General_Info.Product_Info");
1451 1 : return nullptr;
1452 : }
1453 :
1454 12 : auto poDS = std::make_unique<SENTINEL2DatasetContainer>();
1455 : char **papszMD =
1456 6 : SENTINEL2GetUserProductMetadata(psRoot, "Level-1B_User_Product");
1457 6 : poDS->GDALDataset::SetMetadata(papszMD);
1458 6 : CSLDestroy(papszMD);
1459 :
1460 6 : if (!osOriginalXML.empty())
1461 : {
1462 6 : char *apszXMLMD[2] = {nullptr};
1463 6 : apszXMLMD[0] = const_cast<char *>(osOriginalXML.c_str());
1464 6 : apszXMLMD[1] = nullptr;
1465 6 : poDS->GDALDataset::SetMetadata(apszXMLMD, "xml:SENTINEL2");
1466 : }
1467 :
1468 6 : const char *pszPosList = CPLGetXMLValue(
1469 : psRoot,
1470 : "=Level-1B_User_Product.Geometric_Info.Product_Footprint."
1471 : "Product_Footprint.Global_Footprint.EXT_POS_LIST",
1472 : nullptr);
1473 6 : if (pszPosList != nullptr)
1474 : {
1475 6 : CPLString osPolygon = SENTINEL2GetPolygonWKTFromPosList(pszPosList);
1476 3 : if (!osPolygon.empty())
1477 3 : poDS->GDALDataset::SetMetadataItem("FOOTPRINT", osPolygon.c_str());
1478 : }
1479 :
1480 6 : const char *pszMode = CSLFetchNameValueDef(poOpenInfo->papszOpenOptions,
1481 : "L1B_MODE", "DEFAULT");
1482 :
1483 : // Detect if this L1B product has associated geolocation arrays, as provided
1484 : // by the Sen2VM MPC project.
1485 : const std::string osDatastripRoot = CPLFormFilenameSafe(
1486 12 : CPLGetPathSafe(poOpenInfo->pszFilename).c_str(), "DATASTRIP", nullptr);
1487 : VSIStatBufL sStat;
1488 6 : if (VSIStatL(osDatastripRoot.c_str(), &sStat) == 0 &&
1489 6 : VSI_ISDIR(sStat.st_mode) && !EQUAL(pszMode, "GRANULE"))
1490 : {
1491 1 : int iSubDSNum = 1;
1492 1 : const CPLStringList aosList(VSIReadDir(osDatastripRoot.c_str()));
1493 4 : for (const char *pszDatastripId : aosList)
1494 : {
1495 3 : if (IsS2Prefixed(pszDatastripId, "_OPER_MSI_L1B_"))
1496 : {
1497 : const std::string osGEO_DATA = CPLFormFilenameSafe(
1498 1 : CPLFormFilenameSafe(osDatastripRoot.c_str(), pszDatastripId,
1499 : nullptr)
1500 : .c_str(),
1501 2 : "GEO_DATA", nullptr);
1502 2 : const CPLStringList aosListVRT(VSIReadDir(osGEO_DATA.c_str()));
1503 2 : std::vector<std::string> aosVRTs;
1504 6 : for (const char *pszVRT : aosListVRT)
1505 : {
1506 5 : if (EQUAL(CPLGetExtensionSafe(pszVRT).c_str(), "VRT"))
1507 : {
1508 3 : aosVRTs.push_back(pszVRT);
1509 : }
1510 : }
1511 1 : std::sort(aosVRTs.begin(), aosVRTs.end());
1512 4 : for (const std::string &osVRT : aosVRTs)
1513 : {
1514 6 : poDS->GDALDataset::SetMetadataItem(
1515 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
1516 : CPLSPrintf("SENTINEL2_L1B_WITH_GEOLOC:\"%s\":%s",
1517 : poOpenInfo->pszFilename,
1518 6 : CPLGetBasenameSafe(osVRT.c_str()).c_str()),
1519 : "SUBDATASETS");
1520 3 : ++iSubDSNum;
1521 : }
1522 : }
1523 : }
1524 1 : if (iSubDSNum > 1)
1525 : {
1526 1 : return poDS.release();
1527 : }
1528 0 : if (EQUAL(pszMode, "DATASTRIP"))
1529 : {
1530 0 : CPLError(CE_Failure, CPLE_AppDefined,
1531 : "Could not find VRT geolocation files");
1532 0 : return nullptr;
1533 : }
1534 : }
1535 :
1536 10 : std::set<int> oSetResolutions;
1537 10 : std::map<int, std::set<CPLString>> oMapResolutionsToBands;
1538 5 : if (!SENTINEL2GetResolutionSet(psProductInfo, oSetResolutions,
1539 : oMapResolutionsToBands))
1540 : {
1541 3 : CPLDebug("SENTINEL2", "Failed to get resolution set");
1542 3 : return nullptr;
1543 : }
1544 :
1545 4 : std::vector<CPLString> aosGranuleList;
1546 2 : if (!SENTINEL2GetGranuleList(psRoot, SENTINEL2_L1B, poOpenInfo->pszFilename,
1547 : aosGranuleList))
1548 : {
1549 0 : CPLDebug("SENTINEL2", "Failed to get granule list");
1550 0 : return nullptr;
1551 : }
1552 :
1553 : /* Create subdatsets per granules and resolution (10, 20, 60m) */
1554 2 : int iSubDSNum = 1;
1555 4 : for (size_t i = 0; i < aosGranuleList.size(); i++)
1556 : {
1557 8 : for (std::set<int>::const_iterator oIterRes = oSetResolutions.begin();
1558 14 : oIterRes != oSetResolutions.end(); ++oIterRes)
1559 : {
1560 6 : const int nResolution = *oIterRes;
1561 :
1562 12 : poDS->GDALDataset::SetMetadataItem(
1563 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
1564 6 : CPLSPrintf("SENTINEL2_L1B:%s:%dm", aosGranuleList[i].c_str(),
1565 : nResolution),
1566 : "SUBDATASETS");
1567 :
1568 : CPLString osBandNames = SENTINEL2GetBandListForResolution(
1569 12 : oMapResolutionsToBands[nResolution]);
1570 :
1571 : CPLString osDesc(
1572 : CPLSPrintf("Bands %s of granule %s with %dm resolution",
1573 : osBandNames.c_str(),
1574 6 : CPLGetFilename(aosGranuleList[i]), nResolution));
1575 6 : poDS->GDALDataset::SetMetadataItem(
1576 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
1577 : "SUBDATASETS");
1578 :
1579 6 : iSubDSNum++;
1580 : }
1581 : }
1582 :
1583 2 : return poDS.release();
1584 : }
1585 :
1586 : /************************************************************************/
1587 : /* SENTINEL2GetL1BGranuleMetadata() */
1588 : /************************************************************************/
1589 :
1590 15 : static char **SENTINEL2GetL1BGranuleMetadata(CPLXMLNode *psMainMTD)
1591 : {
1592 30 : CPLStringList aosList;
1593 :
1594 15 : CPLXMLNode *psRoot = CPLGetXMLNode(psMainMTD, "=Level-1B_Granule_ID");
1595 15 : if (psRoot == nullptr)
1596 : {
1597 1 : CPLError(CE_Failure, CPLE_AppDefined,
1598 : "Cannot find =Level-1B_Granule_ID");
1599 1 : return nullptr;
1600 : }
1601 14 : CPLXMLNode *psGeneralInfo = CPLGetXMLNode(psRoot, "General_Info");
1602 14 : for (CPLXMLNode *psIter =
1603 14 : (psGeneralInfo ? psGeneralInfo->psChild : nullptr);
1604 54 : psIter != nullptr; psIter = psIter->psNext)
1605 : {
1606 40 : if (psIter->eType != CXT_Element)
1607 0 : continue;
1608 40 : const char *pszValue = CPLGetXMLValue(psIter, nullptr, nullptr);
1609 40 : if (pszValue != nullptr)
1610 : {
1611 34 : aosList.AddNameValue(psIter->pszValue, pszValue);
1612 : }
1613 : }
1614 :
1615 14 : CPLXMLNode *psGeometryHeader = CPLGetXMLNode(
1616 : psRoot, "Geometric_Info.Granule_Position.Geometric_Header");
1617 14 : if (psGeometryHeader != nullptr)
1618 : {
1619 6 : const char *pszVal = CPLGetXMLValue(
1620 : psGeometryHeader, "Incidence_Angles.ZENITH_ANGLE", nullptr);
1621 6 : if (pszVal)
1622 6 : aosList.AddNameValue("INCIDENCE_ZENITH_ANGLE", pszVal);
1623 :
1624 6 : pszVal = CPLGetXMLValue(psGeometryHeader,
1625 : "Incidence_Angles.AZIMUTH_ANGLE", nullptr);
1626 6 : if (pszVal)
1627 6 : aosList.AddNameValue("INCIDENCE_AZIMUTH_ANGLE", pszVal);
1628 :
1629 6 : pszVal = CPLGetXMLValue(psGeometryHeader, "Solar_Angles.ZENITH_ANGLE",
1630 : nullptr);
1631 6 : if (pszVal)
1632 6 : aosList.AddNameValue("SOLAR_ZENITH_ANGLE", pszVal);
1633 :
1634 6 : pszVal = CPLGetXMLValue(psGeometryHeader, "Solar_Angles.AZIMUTH_ANGLE",
1635 : nullptr);
1636 6 : if (pszVal)
1637 6 : aosList.AddNameValue("SOLAR_AZIMUTH_ANGLE", pszVal);
1638 : }
1639 :
1640 14 : CPLXMLNode *psQII = CPLGetXMLNode(psRoot, "Quality_Indicators_Info");
1641 14 : if (psQII != nullptr)
1642 : {
1643 6 : CPLXMLNode *psICCQI = CPLGetXMLNode(psQII, "Image_Content_QI");
1644 6 : for (CPLXMLNode *psIter = (psICCQI ? psICCQI->psChild : nullptr);
1645 18 : psIter != nullptr; psIter = psIter->psNext)
1646 : {
1647 12 : if (psIter->eType != CXT_Element)
1648 0 : continue;
1649 12 : if (psIter->psChild != nullptr &&
1650 12 : psIter->psChild->eType == CXT_Text)
1651 : {
1652 12 : aosList.AddNameValue(psIter->pszValue,
1653 12 : psIter->psChild->pszValue);
1654 : }
1655 : }
1656 : }
1657 :
1658 14 : return aosList.StealList();
1659 : }
1660 :
1661 : /************************************************************************/
1662 : /* SENTINEL2GetTilename() */
1663 : /************************************************************************/
1664 :
1665 : static CPLString
1666 343 : SENTINEL2GetTilename(const std::string &osGranulePath,
1667 : const std::string &osGranuleName,
1668 : const std::string &osBandName,
1669 : const std::string &osProductURI = CPLString(),
1670 : bool bIsPreview = false, int nPrecisionL2A = 0)
1671 : {
1672 343 : bool granuleNameMatchTilename = true;
1673 686 : CPLString osJPEG2000Name(osGranuleName);
1674 343 : if (osJPEG2000Name.size() > 7 &&
1675 496 : osJPEG2000Name[osJPEG2000Name.size() - 7] == '_' &&
1676 153 : osJPEG2000Name[osJPEG2000Name.size() - 6] == 'N')
1677 : {
1678 136 : osJPEG2000Name.resize(osJPEG2000Name.size() - 7);
1679 : }
1680 :
1681 : const SENTINEL2_L2A_BandDescription *psL2ABandDesc =
1682 343 : (nPrecisionL2A) ? SENTINEL2GetL2ABandDesc(osBandName.c_str()) : nullptr;
1683 :
1684 343 : CPLString osTile(osGranulePath);
1685 343 : const char chSeparator = SENTINEL2GetPathSeparator(osTile);
1686 343 : if (!osTile.empty())
1687 343 : osTile += chSeparator;
1688 343 : bool procBaseLineIs1 = false;
1689 553 : if (osJPEG2000Name.size() > 12 && osJPEG2000Name[8] == '_' &&
1690 210 : osJPEG2000Name[12] == '_')
1691 210 : procBaseLineIs1 = true;
1692 343 : if (bIsPreview ||
1693 65 : (psL2ABandDesc != nullptr && (psL2ABandDesc->eLocation == TL_QI_DATA)))
1694 : {
1695 43 : osTile += "QI_DATA";
1696 43 : osTile += chSeparator;
1697 43 : if (procBaseLineIs1)
1698 : {
1699 25 : if (atoi(osBandName.c_str()) > 0)
1700 : {
1701 21 : osJPEG2000Name[9] = 'P';
1702 21 : osJPEG2000Name[10] = 'V';
1703 21 : osJPEG2000Name[11] = 'I';
1704 : }
1705 4 : else if (nPrecisionL2A && osBandName.size() == 3)
1706 : {
1707 4 : osJPEG2000Name[9] = osBandName[0];
1708 4 : osJPEG2000Name[10] = osBandName[1];
1709 4 : osJPEG2000Name[11] = osBandName[2];
1710 : }
1711 25 : osTile += osJPEG2000Name;
1712 : }
1713 : else
1714 : {
1715 18 : osTile += "MSK_";
1716 18 : osTile += osBandName;
1717 18 : osTile += "PRB";
1718 : }
1719 43 : if (nPrecisionL2A && !bIsPreview)
1720 22 : osTile += CPLSPrintf("_%02dm", nPrecisionL2A);
1721 : }
1722 : else
1723 : {
1724 300 : osTile += "IMG_DATA";
1725 300 : osTile += chSeparator;
1726 343 : if (((psL2ABandDesc != nullptr &&
1727 300 : psL2ABandDesc->eLocation == TL_IMG_DATA_Rxxm) ||
1728 682 : (psL2ABandDesc == nullptr && nPrecisionL2A != 0)) &&
1729 125 : (!procBaseLineIs1 || osBandName != "SCL"))
1730 : {
1731 123 : osTile += CPLSPrintf("R%02dm", nPrecisionL2A);
1732 123 : osTile += chSeparator;
1733 : }
1734 300 : if (procBaseLineIs1)
1735 : {
1736 185 : if (atoi(osBandName.c_str()) > 0)
1737 : {
1738 179 : osJPEG2000Name[9] = 'M';
1739 179 : osJPEG2000Name[10] = 'S';
1740 179 : osJPEG2000Name[11] = 'I';
1741 : }
1742 6 : else if (nPrecisionL2A && osBandName.size() == 3)
1743 : {
1744 6 : osJPEG2000Name[9] = osBandName[0];
1745 6 : osJPEG2000Name[10] = osBandName[1];
1746 6 : osJPEG2000Name[11] = osBandName[2];
1747 : }
1748 : }
1749 210 : else if (osProductURI.size() > 44 &&
1750 210 : osProductURI.substr(3, 8) == "_MSIL2A_")
1751 : {
1752 95 : osTile += osProductURI.substr(38, 6);
1753 95 : osTile += osProductURI.substr(10, 16);
1754 95 : granuleNameMatchTilename = false;
1755 : }
1756 : else
1757 : {
1758 20 : CPLDebug("SENTINEL2", "Invalid granule path: %s",
1759 : osGranulePath.c_str());
1760 : }
1761 300 : if (granuleNameMatchTilename)
1762 205 : osTile += osJPEG2000Name;
1763 300 : if (atoi(osBandName.c_str()) > 0)
1764 : {
1765 257 : osTile += "_B";
1766 257 : if (osBandName.size() == 3 && osBandName[0] == '0')
1767 12 : osTile += osBandName.substr(1);
1768 : else
1769 245 : osTile += osBandName;
1770 : }
1771 43 : else if (!procBaseLineIs1)
1772 : {
1773 37 : osTile += "_";
1774 37 : osTile += osBandName;
1775 : }
1776 300 : if (nPrecisionL2A)
1777 125 : osTile += CPLSPrintf("_%02dm", nPrecisionL2A);
1778 : }
1779 343 : osTile += ".jp2";
1780 686 : return osTile;
1781 : }
1782 :
1783 : /************************************************************************/
1784 : /* SENTINEL2GetMainMTDFilenameFromGranuleMTD() */
1785 : /************************************************************************/
1786 :
1787 : static CPLString
1788 29 : SENTINEL2GetMainMTDFilenameFromGranuleMTD(const char *pszFilename)
1789 : {
1790 : // Look for product MTD file
1791 58 : CPLString osTopDir(CPLFormFilenameSafe(
1792 58 : CPLFormFilenameSafe(CPLGetDirnameSafe(pszFilename).c_str(), "..",
1793 : nullptr)
1794 : .c_str(),
1795 58 : "..", nullptr));
1796 :
1797 : // Workaround to avoid long filenames on Windows
1798 29 : if (CPLIsFilenameRelative(pszFilename))
1799 : {
1800 : // GRANULE/bla/bla.xml
1801 36 : const std::string osPath = CPLGetPathSafe(pszFilename);
1802 19 : if (osPath.find('/') != std::string::npos ||
1803 1 : osPath.find('\\') != std::string::npos)
1804 : {
1805 17 : osTopDir = CPLGetPathSafe(CPLGetPathSafe(osPath.c_str()).c_str());
1806 17 : if (osTopDir == "")
1807 0 : osTopDir = ".";
1808 : }
1809 : }
1810 :
1811 29 : char **papszContents = VSIReadDir(osTopDir);
1812 29 : CPLString osMainMTD;
1813 228 : for (char **papszIter = papszContents; papszIter && *papszIter; ++papszIter)
1814 : {
1815 34 : if (strlen(*papszIter) >= strlen("S2A_XXXX_MTD") &&
1816 252 : IsS2Prefixed(*papszIter, "") &&
1817 19 : EQUALN(*papszIter + strlen("S2A_XXXX"), "_MTD", 4))
1818 : {
1819 19 : osMainMTD = CPLFormFilenameSafe(osTopDir, *papszIter, nullptr);
1820 19 : break;
1821 : }
1822 : }
1823 29 : CSLDestroy(papszContents);
1824 58 : return osMainMTD;
1825 : }
1826 :
1827 : /************************************************************************/
1828 : /* SENTINEL2GetResolutionSetAndMainMDFromGranule() */
1829 : /************************************************************************/
1830 :
1831 29 : static void SENTINEL2GetResolutionSetAndMainMDFromGranule(
1832 : const char *pszFilename, const char *pszRootPathWithoutEqual,
1833 : int nResolutionOfInterest, std::set<int> &oSetResolutions,
1834 : std::map<int, std::set<CPLString>> &oMapResolutionsToBands, char **&papszMD,
1835 : CPLXMLNode **ppsRootMainMTD)
1836 : {
1837 58 : CPLString osMainMTD(SENTINEL2GetMainMTDFilenameFromGranuleMTD(pszFilename));
1838 :
1839 : // Parse product MTD if available
1840 29 : papszMD = nullptr;
1841 48 : if (!osMainMTD.empty() &&
1842 : /* env var for debug only */
1843 19 : CPLTestBool(CPLGetConfigOption("SENTINEL2_USE_MAIN_MTD", "YES")))
1844 : {
1845 17 : CPLXMLNode *psRootMainMTD = CPLParseXMLFile(osMainMTD);
1846 17 : if (psRootMainMTD != nullptr)
1847 : {
1848 17 : CPLStripXMLNamespace(psRootMainMTD, nullptr, TRUE);
1849 :
1850 17 : CPLXMLNode *psProductInfo = CPLGetXMLNode(
1851 : psRootMainMTD, CPLSPrintf("=%s.General_Info.Product_Info",
1852 : pszRootPathWithoutEqual));
1853 17 : if (psProductInfo != nullptr)
1854 : {
1855 17 : SENTINEL2GetResolutionSet(psProductInfo, oSetResolutions,
1856 : oMapResolutionsToBands);
1857 : }
1858 :
1859 17 : papszMD = SENTINEL2GetUserProductMetadata(psRootMainMTD,
1860 : pszRootPathWithoutEqual);
1861 17 : if (ppsRootMainMTD != nullptr)
1862 6 : *ppsRootMainMTD = psRootMainMTD;
1863 : else
1864 11 : CPLDestroyXMLNode(psRootMainMTD);
1865 : }
1866 : }
1867 : else
1868 : {
1869 : // If no main MTD file found, then probe all bands for resolution (of
1870 : // interest if there's one, or all resolutions otherwise)
1871 168 : for (const auto &sBandDesc : asBandDesc)
1872 : {
1873 156 : if (nResolutionOfInterest != 0 &&
1874 117 : sBandDesc.nResolution != nResolutionOfInterest)
1875 : {
1876 83 : continue;
1877 : }
1878 : CPLString osBandName =
1879 146 : sBandDesc.pszBandName + 1; /* skip B character */
1880 73 : if (atoi(osBandName) < 10)
1881 62 : osBandName = "0" + osBandName;
1882 :
1883 : CPLString osTile(SENTINEL2GetTilename(
1884 146 : CPLGetPathSafe(pszFilename).c_str(),
1885 365 : CPLGetBasenameSafe(pszFilename).c_str(), osBandName));
1886 : VSIStatBufL sStat;
1887 73 : if (VSIStatExL(osTile, &sStat, VSI_STAT_EXISTS_FLAG) == 0)
1888 : {
1889 20 : oMapResolutionsToBands[sBandDesc.nResolution].insert(
1890 20 : std::move(osBandName));
1891 20 : oSetResolutions.insert(sBandDesc.nResolution);
1892 : }
1893 : }
1894 : }
1895 29 : }
1896 :
1897 : /************************************************************************/
1898 : /* OpenL1BGranule() */
1899 : /************************************************************************/
1900 :
1901 18 : GDALDataset *SENTINEL2Dataset::OpenL1BGranule(const char *pszFilename,
1902 : CPLXMLNode **ppsRoot,
1903 : int nResolutionOfInterest,
1904 : std::set<CPLString> *poBandSet)
1905 : {
1906 18 : CPLXMLNode *psRoot = CPLParseXMLFile(pszFilename);
1907 18 : if (psRoot == nullptr)
1908 : {
1909 3 : CPLDebug("SENTINEL2", "Cannot XML parse %s", pszFilename);
1910 3 : return nullptr;
1911 : }
1912 :
1913 15 : char *pszOriginalXML = CPLSerializeXMLTree(psRoot);
1914 30 : CPLString osOriginalXML;
1915 15 : if (pszOriginalXML)
1916 15 : osOriginalXML = pszOriginalXML;
1917 15 : CPLFree(pszOriginalXML);
1918 :
1919 30 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRoot);
1920 15 : CPLStripXMLNamespace(psRoot, nullptr, TRUE);
1921 :
1922 15 : SENTINEL2DatasetContainer *poDS = new SENTINEL2DatasetContainer();
1923 :
1924 15 : if (!osOriginalXML.empty())
1925 : {
1926 : char *apszXMLMD[2];
1927 15 : apszXMLMD[0] = const_cast<char *>(osOriginalXML.c_str());
1928 15 : apszXMLMD[1] = nullptr;
1929 15 : poDS->GDALDataset::SetMetadata(apszXMLMD, "xml:SENTINEL2");
1930 : }
1931 :
1932 30 : std::set<int> oSetResolutions;
1933 15 : std::map<int, std::set<CPLString>> oMapResolutionsToBands;
1934 15 : char **papszMD = nullptr;
1935 15 : SENTINEL2GetResolutionSetAndMainMDFromGranule(
1936 : pszFilename, "Level-1B_User_Product", nResolutionOfInterest,
1937 : oSetResolutions, oMapResolutionsToBands, papszMD, nullptr);
1938 15 : if (poBandSet != nullptr)
1939 12 : *poBandSet = oMapResolutionsToBands[nResolutionOfInterest];
1940 :
1941 15 : char **papszGranuleMD = SENTINEL2GetL1BGranuleMetadata(psRoot);
1942 15 : papszMD = CSLMerge(papszMD, papszGranuleMD);
1943 15 : CSLDestroy(papszGranuleMD);
1944 :
1945 : // Remove CLOUD_COVERAGE_ASSESSMENT that comes from main metadata, if
1946 : // granule CLOUDY_PIXEL_PERCENTAGE is present.
1947 21 : if (CSLFetchNameValue(papszMD, "CLOUDY_PIXEL_PERCENTAGE") != nullptr &&
1948 6 : CSLFetchNameValue(papszMD, "CLOUD_COVERAGE_ASSESSMENT") != nullptr)
1949 : {
1950 6 : papszMD =
1951 6 : CSLSetNameValue(papszMD, "CLOUD_COVERAGE_ASSESSMENT", nullptr);
1952 : }
1953 :
1954 15 : poDS->GDALDataset::SetMetadata(papszMD);
1955 15 : CSLDestroy(papszMD);
1956 :
1957 : // Get the footprint
1958 : const char *pszPosList =
1959 15 : CPLGetXMLValue(psRoot,
1960 : "=Level-1B_Granule_ID.Geometric_Info.Granule_Footprint."
1961 : "Granule_Footprint.Footprint.EXT_POS_LIST",
1962 : nullptr);
1963 15 : if (pszPosList != nullptr)
1964 : {
1965 12 : CPLString osPolygon = SENTINEL2GetPolygonWKTFromPosList(pszPosList);
1966 6 : if (!osPolygon.empty())
1967 6 : poDS->GDALDataset::SetMetadataItem("FOOTPRINT", osPolygon.c_str());
1968 : }
1969 :
1970 : /* Create subdatsets per resolution (10, 20, 60m) */
1971 15 : int iSubDSNum = 1;
1972 37 : for (std::set<int>::const_iterator oIterRes = oSetResolutions.begin();
1973 59 : oIterRes != oSetResolutions.end(); ++oIterRes)
1974 : {
1975 22 : const int nResolution = *oIterRes;
1976 :
1977 22 : poDS->GDALDataset::SetMetadataItem(
1978 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
1979 : CPLSPrintf("SENTINEL2_L1B:%s:%dm", pszFilename, nResolution),
1980 : "SUBDATASETS");
1981 :
1982 : CPLString osBandNames = SENTINEL2GetBandListForResolution(
1983 44 : oMapResolutionsToBands[nResolution]);
1984 :
1985 : CPLString osDesc(CPLSPrintf("Bands %s with %dm resolution",
1986 22 : osBandNames.c_str(), nResolution));
1987 22 : poDS->GDALDataset::SetMetadataItem(
1988 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
1989 : "SUBDATASETS");
1990 :
1991 22 : iSubDSNum++;
1992 : }
1993 :
1994 15 : if (ppsRoot != nullptr)
1995 : {
1996 12 : *ppsRoot = oXMLHolder.Release();
1997 : }
1998 :
1999 15 : return poDS;
2000 : }
2001 :
2002 : /************************************************************************/
2003 : /* SENTINEL2SetBandMetadata() */
2004 : /************************************************************************/
2005 :
2006 212 : static void SENTINEL2SetBandMetadata(GDALRasterBand *poBand,
2007 : const std::string &osBandName)
2008 : {
2009 424 : CPLString osLookupBandName(osBandName);
2010 212 : if (osLookupBandName[0] == '0')
2011 141 : osLookupBandName = osLookupBandName.substr(1);
2012 212 : if (atoi(osLookupBandName) > 0)
2013 171 : osLookupBandName = "B" + osLookupBandName;
2014 :
2015 424 : CPLString osBandDesc(osLookupBandName);
2016 : const SENTINEL2BandDescription *psBandDesc =
2017 212 : SENTINEL2GetBandDesc(osLookupBandName);
2018 212 : if (psBandDesc != nullptr)
2019 : {
2020 : osBandDesc +=
2021 171 : CPLSPrintf(", central wavelength %d nm", psBandDesc->nWaveLength);
2022 171 : poBand->SetColorInterpretation(psBandDesc->eColorInterp);
2023 171 : poBand->SetMetadataItem("BANDNAME", psBandDesc->pszBandName);
2024 171 : poBand->SetMetadataItem("BANDWIDTH",
2025 171 : CPLSPrintf("%d", psBandDesc->nBandWidth));
2026 171 : poBand->SetMetadataItem("BANDWIDTH_UNIT", "nm");
2027 171 : poBand->SetMetadataItem("WAVELENGTH",
2028 171 : CPLSPrintf("%d", psBandDesc->nWaveLength));
2029 171 : poBand->SetMetadataItem("WAVELENGTH_UNIT", "nm");
2030 :
2031 171 : poBand->SetMetadataItem(
2032 : "CENTRAL_WAVELENGTH_UM",
2033 171 : CPLSPrintf("%.3f", double(psBandDesc->nWaveLength) / 1000),
2034 171 : "IMAGERY");
2035 171 : poBand->SetMetadataItem(
2036 : "FWHM_UM",
2037 171 : CPLSPrintf("%.3f", double(psBandDesc->nBandWidth) / 1000),
2038 171 : "IMAGERY");
2039 : }
2040 : else
2041 : {
2042 : const SENTINEL2_L2A_BandDescription *psL2ABandDesc =
2043 41 : SENTINEL2GetL2ABandDesc(osBandName.c_str());
2044 41 : if (psL2ABandDesc != nullptr)
2045 : {
2046 41 : osBandDesc += ", ";
2047 41 : osBandDesc += psL2ABandDesc->pszBandDescription;
2048 : }
2049 :
2050 41 : poBand->SetMetadataItem("BANDNAME", osBandName.c_str());
2051 : }
2052 212 : poBand->SetDescription(osBandDesc);
2053 212 : }
2054 :
2055 : /************************************************************************/
2056 : /* OpenL1BSubdataset() */
2057 : /************************************************************************/
2058 :
2059 18 : GDALDataset *SENTINEL2Dataset::OpenL1BSubdataset(GDALOpenInfo *poOpenInfo)
2060 : {
2061 36 : CPLString osFilename;
2062 18 : CPLAssert(STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1B:"));
2063 18 : osFilename = poOpenInfo->pszFilename + strlen("SENTINEL2_L1B:");
2064 18 : const char *pszPrecision = strrchr(osFilename.c_str(), ':');
2065 18 : if (pszPrecision == nullptr || pszPrecision == osFilename.c_str())
2066 : {
2067 2 : CPLError(CE_Failure, CPLE_AppDefined,
2068 : "Invalid syntax for SENTINEL2_L1B:");
2069 2 : return nullptr;
2070 : }
2071 16 : const int nSubDSPrecision = atoi(pszPrecision + 1);
2072 16 : if (nSubDSPrecision != RES_10M && nSubDSPrecision != RES_20M &&
2073 : nSubDSPrecision != RES_60M)
2074 : {
2075 2 : CPLError(CE_Failure, CPLE_NotSupported, "Unsupported precision: %d",
2076 : nSubDSPrecision);
2077 2 : return nullptr;
2078 : }
2079 14 : osFilename.resize(pszPrecision - osFilename.c_str());
2080 :
2081 14 : CPLXMLNode *psRoot = nullptr;
2082 28 : std::set<CPLString> oSetBands;
2083 : GDALDataset *poTmpDS =
2084 14 : OpenL1BGranule(osFilename, &psRoot, nSubDSPrecision, &oSetBands);
2085 14 : if (poTmpDS == nullptr)
2086 : {
2087 2 : CPLDebug("SENTINEL2", "Failed to open L1B granule %s",
2088 : osFilename.c_str());
2089 2 : return nullptr;
2090 : }
2091 :
2092 24 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRoot);
2093 :
2094 24 : std::vector<CPLString> aosBands;
2095 31 : for (std::set<CPLString>::const_iterator oIterBandnames = oSetBands.begin();
2096 50 : oIterBandnames != oSetBands.end(); ++oIterBandnames)
2097 : {
2098 19 : aosBands.push_back(*oIterBandnames);
2099 : }
2100 : /* Put 2=Blue, 3=Green, 4=Band bands in RGB order for conveniency */
2101 13 : if (aosBands.size() >= 3 && aosBands[0] == "02" && aosBands[1] == "03" &&
2102 1 : aosBands[2] == "04")
2103 : {
2104 1 : aosBands[0] = "04";
2105 1 : aosBands[2] = "02";
2106 : }
2107 :
2108 12 : int nBits = 0; /* 0 = unknown yet*/
2109 12 : int nValMax = 0; /* 0 = unknown yet*/
2110 12 : int nRows = 0;
2111 12 : int nCols = 0;
2112 12 : CPLXMLNode *psGranuleDimensions = CPLGetXMLNode(
2113 : psRoot, "=Level-1B_Granule_ID.Geometric_Info.Granule_Dimensions");
2114 12 : if (psGranuleDimensions == nullptr)
2115 : {
2116 3 : for (size_t i = 0; i < aosBands.size(); i++)
2117 : {
2118 : CPLString osTile(SENTINEL2GetTilename(
2119 2 : CPLGetPathSafe(osFilename).c_str(),
2120 4 : CPLGetBasenameSafe(osFilename).c_str(), aosBands[i]));
2121 1 : if (SENTINEL2GetTileInfo(osTile, &nCols, &nRows, &nBits))
2122 : {
2123 1 : if (nBits <= 16)
2124 1 : nValMax = (1 << nBits) - 1;
2125 : else
2126 : {
2127 0 : CPLDebug("SENTINEL2", "Unexpected bit depth %d", nBits);
2128 0 : nValMax = 65535;
2129 : }
2130 1 : break;
2131 : }
2132 : }
2133 : }
2134 : else
2135 : {
2136 9 : for (CPLXMLNode *psIter = psGranuleDimensions->psChild;
2137 23 : psIter != nullptr; psIter = psIter->psNext)
2138 : {
2139 22 : if (psIter->eType != CXT_Element)
2140 9 : continue;
2141 26 : if (EQUAL(psIter->pszValue, "Size") &&
2142 13 : atoi(CPLGetXMLValue(psIter, "resolution", "")) ==
2143 : nSubDSPrecision)
2144 : {
2145 8 : const char *pszRows = CPLGetXMLValue(psIter, "NROWS", nullptr);
2146 8 : if (pszRows == nullptr)
2147 : {
2148 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
2149 : "NROWS");
2150 1 : delete poTmpDS;
2151 1 : return nullptr;
2152 : }
2153 7 : const char *pszCols = CPLGetXMLValue(psIter, "NCOLS", nullptr);
2154 7 : if (pszCols == nullptr)
2155 : {
2156 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
2157 : "NCOLS");
2158 1 : delete poTmpDS;
2159 1 : return nullptr;
2160 : }
2161 6 : nRows = atoi(pszRows);
2162 6 : nCols = atoi(pszCols);
2163 6 : break;
2164 : }
2165 : }
2166 : }
2167 10 : if (nRows <= 0 || nCols <= 0)
2168 : {
2169 3 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find granule dimension");
2170 3 : delete poTmpDS;
2171 3 : return nullptr;
2172 : }
2173 :
2174 7 : SENTINEL2Dataset *poDS = new SENTINEL2Dataset(nCols, nRows);
2175 7 : poDS->aosNonJP2Files.push_back(osFilename);
2176 :
2177 : // Transfer metadata
2178 7 : poDS->GDALDataset::SetMetadata(poTmpDS->GetMetadata());
2179 7 : poDS->GDALDataset::SetMetadata(poTmpDS->GetMetadata("xml:SENTINEL2"),
2180 : "xml:SENTINEL2");
2181 :
2182 7 : delete poTmpDS;
2183 :
2184 : /* -------------------------------------------------------------------- */
2185 : /* Initialize bands. */
2186 : /* -------------------------------------------------------------------- */
2187 :
2188 21 : int nSaturatedVal = atoi(CSLFetchNameValueDef(
2189 7 : poDS->GetMetadata(), "SPECIAL_VALUE_SATURATED", "-1"));
2190 7 : int nNodataVal = atoi(CSLFetchNameValueDef(poDS->GetMetadata(),
2191 : "SPECIAL_VALUE_NODATA", "-1"));
2192 :
2193 : const bool bAlpha =
2194 7 : CPLTestBool(SENTINEL2GetOption(poOpenInfo, "ALPHA", "FALSE"));
2195 7 : const int nBands = ((bAlpha) ? 1 : 0) + static_cast<int>(aosBands.size());
2196 7 : const int nAlphaBand = (!bAlpha) ? 0 : nBands;
2197 7 : const GDALDataType eDT = GDT_UInt16;
2198 :
2199 27 : for (int nBand = 1; nBand <= nBands; nBand++)
2200 : {
2201 20 : VRTSourcedRasterBand *poBand = nullptr;
2202 :
2203 20 : if (nBand != nAlphaBand)
2204 : {
2205 19 : poBand = new VRTSourcedRasterBand(
2206 19 : poDS, nBand, eDT, poDS->nRasterXSize, poDS->nRasterYSize);
2207 : }
2208 : else
2209 : {
2210 1 : poBand = new SENTINEL2AlphaBand(
2211 : poDS, nBand, eDT, poDS->nRasterXSize, poDS->nRasterYSize,
2212 1 : nSaturatedVal, nNodataVal);
2213 : }
2214 :
2215 20 : poDS->SetBand(nBand, poBand);
2216 20 : if (nBand == nAlphaBand)
2217 1 : poBand->SetColorInterpretation(GCI_AlphaBand);
2218 :
2219 20 : CPLString osBandName;
2220 20 : if (nBand != nAlphaBand)
2221 : {
2222 19 : osBandName = aosBands[nBand - 1];
2223 19 : SENTINEL2SetBandMetadata(poBand, osBandName);
2224 : }
2225 : else
2226 1 : osBandName = aosBands[0];
2227 :
2228 : CPLString osTile(SENTINEL2GetTilename(
2229 40 : CPLGetPathSafe(osFilename).c_str(),
2230 80 : CPLGetBasenameSafe(osFilename).c_str(), osBandName));
2231 :
2232 20 : bool bTileFound = false;
2233 20 : if (nValMax == 0)
2234 : {
2235 : /* It is supposed to be 12 bits, but some products have 15 bits */
2236 6 : if (SENTINEL2GetTileInfo(osTile, nullptr, nullptr, &nBits))
2237 : {
2238 5 : bTileFound = true;
2239 5 : if (nBits <= 16)
2240 5 : nValMax = (1 << nBits) - 1;
2241 : else
2242 : {
2243 0 : CPLDebug("SENTINEL2", "Unexpected bit depth %d", nBits);
2244 0 : nValMax = 65535;
2245 : }
2246 : }
2247 : }
2248 : else
2249 : {
2250 : VSIStatBufL sStat;
2251 14 : if (VSIStatExL(osTile, &sStat, VSI_STAT_EXISTS_FLAG) == 0)
2252 14 : bTileFound = true;
2253 : }
2254 20 : if (!bTileFound)
2255 : {
2256 1 : CPLError(CE_Warning, CPLE_AppDefined,
2257 : "Tile %s not found on filesystem. Skipping it",
2258 : osTile.c_str());
2259 1 : continue;
2260 : }
2261 :
2262 19 : if (nBand != nAlphaBand)
2263 : {
2264 18 : poBand->AddSimpleSource(osTile, 1, 0, 0, poDS->nRasterXSize,
2265 18 : poDS->nRasterYSize, 0, 0,
2266 18 : poDS->nRasterXSize, poDS->nRasterYSize);
2267 : }
2268 : else
2269 : {
2270 1 : poBand->AddComplexSource(osTile, 1, 0, 0, poDS->nRasterXSize,
2271 1 : poDS->nRasterYSize, 0, 0,
2272 1 : poDS->nRasterXSize, poDS->nRasterYSize,
2273 : nValMax /* offset */, 0 /* scale */);
2274 : }
2275 :
2276 19 : if ((nBits % 8) != 0)
2277 : {
2278 19 : poBand->SetMetadataItem("NBITS", CPLSPrintf("%d", nBits),
2279 19 : "IMAGE_STRUCTURE");
2280 : }
2281 : }
2282 :
2283 : /* -------------------------------------------------------------------- */
2284 : /* Add georeferencing. */
2285 : /* -------------------------------------------------------------------- */
2286 : // const char* pszDirection =
2287 : // poDS->GetMetadataItem("DATATAKE_1_SENSING_ORBIT_DIRECTION");
2288 7 : const char *pszFootprint = poDS->GetMetadataItem("FOOTPRINT");
2289 7 : if (pszFootprint != nullptr)
2290 : {
2291 : // For descending orbits, we have observed that the order of points in
2292 : // the polygon is UL, LL, LR, UR. That might not be true for ascending
2293 : // orbits but let's assume it...
2294 4 : OGRGeometry *poGeom = nullptr;
2295 4 : if (OGRGeometryFactory::createFromWkt(pszFootprint, nullptr, &poGeom) ==
2296 4 : OGRERR_NONE &&
2297 8 : poGeom != nullptr &&
2298 4 : wkbFlatten(poGeom->getGeometryType()) == wkbPolygon)
2299 : {
2300 : OGRLinearRing *poRing =
2301 4 : reinterpret_cast<OGRPolygon *>(poGeom)->getExteriorRing();
2302 4 : if (poRing != nullptr && poRing->getNumPoints() == 5)
2303 : {
2304 : GDAL_GCP asGCPList[5];
2305 4 : memset(asGCPList, 0, sizeof(asGCPList));
2306 20 : for (int i = 0; i < 4; i++)
2307 : {
2308 16 : asGCPList[i].dfGCPX = poRing->getX(i);
2309 16 : asGCPList[i].dfGCPY = poRing->getY(i);
2310 16 : asGCPList[i].dfGCPZ = poRing->getZ(i);
2311 : }
2312 4 : asGCPList[0].dfGCPPixel = 0;
2313 4 : asGCPList[0].dfGCPLine = 0;
2314 4 : asGCPList[1].dfGCPPixel = 0;
2315 4 : asGCPList[1].dfGCPLine = poDS->nRasterYSize;
2316 4 : asGCPList[2].dfGCPPixel = poDS->nRasterXSize;
2317 4 : asGCPList[2].dfGCPLine = poDS->nRasterYSize;
2318 4 : asGCPList[3].dfGCPPixel = poDS->nRasterXSize;
2319 4 : asGCPList[3].dfGCPLine = 0;
2320 :
2321 4 : int nGCPCount = 4;
2322 : CPLXMLNode *psGeometryHeader =
2323 4 : CPLGetXMLNode(psRoot, "=Level-1B_Granule_ID.Geometric_Info."
2324 : "Granule_Position.Geometric_Header");
2325 4 : if (psGeometryHeader != nullptr)
2326 : {
2327 4 : const char *pszGC = CPLGetXMLValue(
2328 : psGeometryHeader, "GROUND_CENTER", nullptr);
2329 : const char *pszQLCenter =
2330 4 : CPLGetXMLValue(psGeometryHeader, "QL_CENTER", nullptr);
2331 4 : if (pszGC != nullptr && pszQLCenter != nullptr &&
2332 4 : EQUAL(pszQLCenter, "0 0"))
2333 : {
2334 4 : char **papszTokens = CSLTokenizeString(pszGC);
2335 4 : if (CSLCount(papszTokens) >= 2)
2336 : {
2337 4 : nGCPCount = 5;
2338 4 : asGCPList[4].dfGCPX = CPLAtof(papszTokens[1]);
2339 4 : asGCPList[4].dfGCPY = CPLAtof(papszTokens[0]);
2340 4 : if (CSLCount(papszTokens) >= 3)
2341 4 : asGCPList[4].dfGCPZ = CPLAtof(papszTokens[2]);
2342 4 : asGCPList[4].dfGCPLine = poDS->nRasterYSize / 2.0;
2343 4 : asGCPList[4].dfGCPPixel = poDS->nRasterXSize / 2.0;
2344 : }
2345 4 : CSLDestroy(papszTokens);
2346 : }
2347 : }
2348 :
2349 4 : poDS->SetGCPs(nGCPCount, asGCPList, SRS_WKT_WGS84_LAT_LONG);
2350 4 : GDALDeinitGCPs(nGCPCount, asGCPList);
2351 : }
2352 : }
2353 4 : delete poGeom;
2354 : }
2355 :
2356 : /* -------------------------------------------------------------------- */
2357 : /* Initialize overview information. */
2358 : /* -------------------------------------------------------------------- */
2359 7 : poDS->SetDescription(poOpenInfo->pszFilename);
2360 7 : CPLString osOverviewFile;
2361 : osOverviewFile =
2362 7 : CPLSPrintf("%s_%dm.tif.ovr", osFilename.c_str(), nSubDSPrecision);
2363 7 : poDS->SetMetadataItem("OVERVIEW_FILE", osOverviewFile, "OVERVIEWS");
2364 7 : poDS->oOvManager.Initialize(poDS, ":::VIRTUAL:::");
2365 :
2366 7 : return poDS;
2367 : }
2368 :
2369 : /************************************************************************/
2370 : /* OpenL1BSubdatasetWithGeoloc() */
2371 : /************************************************************************/
2372 :
2373 : GDALDataset *
2374 9 : SENTINEL2Dataset::OpenL1BSubdatasetWithGeoloc(GDALOpenInfo *poOpenInfo)
2375 : {
2376 18 : CPLString osFilename;
2377 9 : CPLAssert(
2378 : STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1B_WITH_GEOLOC:"));
2379 : const CPLStringList aosTokens(
2380 9 : CSLTokenizeString2(poOpenInfo->pszFilename, ":",
2381 18 : CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES));
2382 9 : if (aosTokens.size() != 3)
2383 : {
2384 1 : CPLError(CE_Failure, CPLE_AppDefined,
2385 : "OpenL1BSubdatasetWithGeoloc(): invalid number of tokens in "
2386 : "subdataset name");
2387 1 : return nullptr;
2388 : }
2389 8 : const char *pszMainXMLFilename = aosTokens[1];
2390 8 : const char *pszGeolocVRT = aosTokens[2];
2391 :
2392 8 : const size_t nLen = strlen(pszGeolocVRT);
2393 8 : if (nLen <= strlen("_Dxx_Byy") || pszGeolocVRT[nLen - 7] != 'D' ||
2394 7 : pszGeolocVRT[nLen - 3] != 'B')
2395 : {
2396 1 : CPLError(CE_Failure, CPLE_AppDefined,
2397 : "Invalid subdataset component name");
2398 1 : return nullptr;
2399 : }
2400 :
2401 : // Find in which datastrip we are
2402 : const std::string osDatastrip(
2403 7 : CPLString(pszGeolocVRT, nLen - strlen("_Dxx_Byy"))
2404 21 : .replaceAll("_GEO_", "_MSI_"));
2405 14 : const std::string osDetectorId(pszGeolocVRT + nLen - 6, 2);
2406 14 : const std::string osBandId(pszGeolocVRT + nLen - 2, 2);
2407 :
2408 7 : const char chSeparator = SENTINEL2GetPathSeparator(pszMainXMLFilename);
2409 7 : const char achSeparator[] = {chSeparator, 0};
2410 :
2411 : const std::string osDatastripRoot(
2412 7 : std::string(CPLGetDirnameSafe(pszMainXMLFilename))
2413 7 : .append(achSeparator)
2414 14 : .append("DATASTRIP"));
2415 14 : const CPLStringList aosList(VSIReadDir(osDatastripRoot.c_str()));
2416 14 : std::string osDatastripSubdir;
2417 10 : for (const char *pszDatastripId : aosList)
2418 : {
2419 8 : if (STARTS_WITH_CI(pszDatastripId, osDatastrip.c_str()))
2420 : {
2421 5 : osDatastripSubdir = pszDatastripId;
2422 5 : break;
2423 : }
2424 : }
2425 7 : if (osDatastripSubdir.empty())
2426 : {
2427 2 : CPLError(CE_Failure, CPLE_AppDefined,
2428 : "Cannot find a file in %s starting with %s",
2429 : osDatastripRoot.c_str(), osDatastrip.c_str());
2430 2 : return nullptr;
2431 : }
2432 :
2433 5 : const std::string osDatastripSubdirFull = std::string(osDatastripRoot)
2434 5 : .append(achSeparator)
2435 10 : .append(osDatastripSubdir);
2436 5 : CPL_IGNORE_RET_VAL(osDatastripRoot);
2437 : const std::string osXMLDatastrip =
2438 5 : std::string(osDatastripSubdirFull)
2439 5 : .append(achSeparator)
2440 10 : .append(CPLString(osDatastrip).replaceAll("_MSI_", "_MTD_"))
2441 10 : .append(".xml");
2442 5 : if (CPLHasPathTraversal(osXMLDatastrip.c_str()))
2443 : {
2444 0 : CPLError(CE_Failure, CPLE_NotSupported, "Path traversal detected in %s",
2445 : osXMLDatastrip.c_str());
2446 0 : return nullptr;
2447 : }
2448 10 : CPLXMLTreeCloser poXML(CPLParseXMLFile(osXMLDatastrip.c_str()));
2449 5 : if (!poXML)
2450 : {
2451 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot XML parse %s",
2452 : osXMLDatastrip.c_str());
2453 0 : return nullptr;
2454 : }
2455 5 : CPLStripXMLNamespace(poXML.get(), nullptr, TRUE);
2456 :
2457 : const CPLXMLNode *psDetectorList =
2458 5 : CPLGetXMLNode(poXML.get(), "=Level-1B_DataStrip_ID.Image_Data_Info."
2459 : "Granules_Information.Detector_List");
2460 5 : if (!psDetectorList)
2461 : {
2462 0 : CPLError(CE_Failure, CPLE_AppDefined,
2463 : "Cannot find <Detector_List> in %s", osXMLDatastrip.c_str());
2464 0 : return nullptr;
2465 : }
2466 :
2467 : struct GranuleDesc
2468 : {
2469 : const char *pszGranuleId = nullptr;
2470 : int nPosition = 0;
2471 : };
2472 :
2473 10 : std::vector<GranuleDesc> asGranules;
2474 :
2475 : // Get the list of granules for the current detector and band
2476 6 : for (const CPLXMLNode *psDetector = psDetectorList->psChild; psDetector;
2477 1 : psDetector = psDetector->psNext)
2478 : {
2479 15 : if (psDetector->eType == CXT_Element &&
2480 10 : EQUAL(psDetector->pszValue, "Detector") &&
2481 5 : CPLGetXMLValue(psDetector, "detectorId", "") == osDetectorId)
2482 : {
2483 : const CPLXMLNode *psGranuleList =
2484 4 : CPLGetXMLNode(psDetector, "Granule_List");
2485 4 : if (psGranuleList)
2486 : {
2487 4 : for (const CPLXMLNode *psGranule = psGranuleList->psChild;
2488 24 : psGranule; psGranule = psGranule->psNext)
2489 : {
2490 20 : if (psGranule->eType == CXT_Element &&
2491 20 : EQUAL(psGranule->pszValue, "Granule"))
2492 : {
2493 : const char *pszGranuleId =
2494 20 : CPLGetXMLValue(psGranule, "granuleId", "");
2495 : const char *pszPosition =
2496 20 : CPLGetXMLValue(psGranule, "POSITION", "");
2497 20 : GranuleDesc sDesc;
2498 20 : sDesc.pszGranuleId = pszGranuleId;
2499 20 : sDesc.nPosition = atoi(pszPosition);
2500 20 : asGranules.push_back(sDesc);
2501 : }
2502 : }
2503 : }
2504 4 : break;
2505 : }
2506 : }
2507 5 : if (asGranules.empty())
2508 : {
2509 1 : CPLError(CE_Failure, CPLE_AppDefined,
2510 : "Cannot find granules for detector %s in %s",
2511 : osDetectorId.c_str(), osXMLDatastrip.c_str());
2512 1 : return nullptr;
2513 : }
2514 4 : std::sort(asGranules.begin(), asGranules.end(),
2515 32 : [](const GranuleDesc &sDescA, const GranuleDesc &sDescB)
2516 32 : { return sDescA.nPosition < sDescB.nPosition; });
2517 4 : const int nGranuleCount = static_cast<int>(asGranules.size());
2518 4 : constexpr int ROWS_PER_10M_GRANULE = 2304;
2519 24 : for (int i = 0; i < nGranuleCount; ++i)
2520 : {
2521 20 : if (asGranules[i].nPosition != 1 + i * ROWS_PER_10M_GRANULE)
2522 : {
2523 0 : CPLError(
2524 : CE_Failure, CPLE_NotSupported,
2525 : "Granule %s is declared to be at position %d, was expecting %d",
2526 0 : asGranules[i].pszGranuleId, asGranules[i].nPosition,
2527 0 : 1 + i * ROWS_PER_10M_GRANULE);
2528 0 : return nullptr;
2529 : }
2530 : }
2531 :
2532 8 : const std::string osBandName = std::string("B").append(
2533 4 : osBandId == "8A"
2534 12 : ? osBandId
2535 8 : : std::string(CPLSPrintf("%d", atoi(osBandId.c_str()))));
2536 : const SENTINEL2BandDescription *psBandDesc =
2537 4 : SENTINEL2GetBandDesc(osBandName.c_str());
2538 4 : if (psBandDesc == nullptr)
2539 : {
2540 1 : CPLError(CE_Failure, CPLE_NotSupported, "Unknown band id %s",
2541 : osBandId.c_str());
2542 1 : return nullptr;
2543 : }
2544 :
2545 3 : const int nResolution = psBandDesc->nResolution;
2546 5 : const int nRowsPerGranule = nResolution == RES_10M ? ROWS_PER_10M_GRANULE
2547 2 : : nResolution == RES_20M ? 1152
2548 : : 384;
2549 5 : const int nColsPerGranule = nResolution == RES_10M ? 2552
2550 2 : : nResolution == RES_20M ? 1276
2551 : : 425;
2552 : auto poDS = std::make_unique<SENTINEL2Dataset>(
2553 6 : nColsPerGranule, nRowsPerGranule * nGranuleCount);
2554 :
2555 3 : constexpr GDALDataType eDT = GDT_UInt16;
2556 : auto poBand = new VRTSourcedRasterBand(
2557 3 : poDS.get(), 1, eDT, poDS->nRasterXSize, poDS->nRasterYSize);
2558 3 : poDS->SetBand(1, poBand);
2559 :
2560 3 : SENTINEL2SetBandMetadata(poBand, osBandId);
2561 :
2562 : // Create the virtual raster by adding each granule
2563 : const std::string osGranuleRoot =
2564 3 : std::string(CPLGetDirnameSafe(pszMainXMLFilename))
2565 3 : .append(achSeparator)
2566 6 : .append("GRANULE");
2567 3 : int nValMax = 0;
2568 3 : int nBits = 0;
2569 18 : for (int i = 0; i < nGranuleCount; ++i)
2570 : {
2571 15 : const auto &sGranule = asGranules[i];
2572 : const std::string osTile(
2573 45 : SENTINEL2GetTilename(std::string(osGranuleRoot)
2574 15 : .append(achSeparator)
2575 15 : .append(sGranule.pszGranuleId),
2576 30 : sGranule.pszGranuleId, osBandId.c_str()));
2577 :
2578 15 : bool bTileFound = false;
2579 15 : if (nValMax == 0)
2580 : {
2581 : /* It is supposed to be 12 bits, but some products have 15 bits */
2582 3 : int nGranuleWidth = 0;
2583 3 : int nGranuleHeight = 0;
2584 3 : if (SENTINEL2GetTileInfo(osTile.c_str(), &nGranuleWidth,
2585 : &nGranuleHeight, &nBits))
2586 : {
2587 3 : bTileFound = true;
2588 3 : if (nBits <= 16)
2589 3 : nValMax = (1 << nBits) - 1;
2590 : else
2591 : {
2592 0 : CPLDebug("SENTINEL2", "Unexpected bit depth %d", nBits);
2593 0 : nValMax = 65535;
2594 : }
2595 3 : if (nGranuleWidth != nColsPerGranule ||
2596 3 : nGranuleHeight != nRowsPerGranule)
2597 : {
2598 0 : CPLError(CE_Failure, CPLE_AppDefined,
2599 : "Tile %s has not expected dimensions. "
2600 : "Got %dx%d, expected %dx%d",
2601 : osTile.c_str(), nGranuleWidth, nGranuleHeight,
2602 : nColsPerGranule, nRowsPerGranule);
2603 0 : return nullptr;
2604 : }
2605 : }
2606 : }
2607 : else
2608 : {
2609 : VSIStatBufL sStat;
2610 12 : if (VSIStatExL(osTile.c_str(), &sStat, VSI_STAT_EXISTS_FLAG) == 0)
2611 12 : bTileFound = true;
2612 : }
2613 15 : if (!bTileFound)
2614 : {
2615 0 : CPLError(CE_Failure, CPLE_AppDefined, "Tile %s not found.",
2616 : osTile.c_str());
2617 0 : return nullptr;
2618 : }
2619 :
2620 15 : poBand->AddSimpleSource(osTile.c_str(), 1, 0, 0, nColsPerGranule,
2621 15 : nRowsPerGranule, 0, i * nRowsPerGranule,
2622 : nColsPerGranule, nRowsPerGranule);
2623 : }
2624 :
2625 3 : if ((nBits % 8) != 0)
2626 : {
2627 0 : poBand->SetMetadataItem("NBITS", CPLSPrintf("%d", nBits),
2628 0 : "IMAGE_STRUCTURE");
2629 : }
2630 :
2631 : // Get metadata from top MTD XML filename
2632 6 : std::set<int> oSetResolutions;
2633 6 : std::map<int, std::set<CPLString>> oMapResolutionsToBands;
2634 3 : char **papszMD = nullptr;
2635 : std::string osGranuleMTD = CPLFormFilenameSafe(
2636 3 : CPLFormFilenameSafe(osGranuleRoot.c_str(), asGranules[0].pszGranuleId,
2637 : nullptr)
2638 : .c_str(),
2639 9 : asGranules[0].pszGranuleId, nullptr);
2640 3 : if (osGranuleMTD.size() > strlen("_NXX.YY") &&
2641 6 : osGranuleMTD[osGranuleMTD.size() - 7] == '_' &&
2642 3 : osGranuleMTD[osGranuleMTD.size() - 6] == 'N')
2643 : {
2644 3 : osGranuleMTD.resize(osGranuleMTD.size() - strlen("_NXX.YY"));
2645 : }
2646 3 : SENTINEL2GetResolutionSetAndMainMDFromGranule(
2647 : osGranuleMTD.c_str(), "Level-1B_User_Product", nResolution,
2648 : oSetResolutions, oMapResolutionsToBands, papszMD, nullptr);
2649 144 : for (const auto [pszKey, pszValue] :
2650 147 : cpl::IterateNameValue(const_cast<CSLConstList>(papszMD)))
2651 : {
2652 72 : if (!EQUAL(pszKey, "FOOTPRINT"))
2653 72 : poDS->SetMetadataItem(pszKey, pszValue);
2654 : }
2655 3 : CSLDestroy(papszMD);
2656 :
2657 : // Attach geolocation array
2658 3 : const std::string osGeolocVRT = std::string(osDatastripSubdirFull)
2659 3 : .append(achSeparator)
2660 3 : .append("GEO_DATA")
2661 3 : .append(achSeparator)
2662 3 : .append(pszGeolocVRT)
2663 6 : .append(".vrt");
2664 3 : CPL_IGNORE_RET_VAL(osDatastripSubdirFull);
2665 : auto poGeolocDS = std::unique_ptr<GDALDataset>(
2666 6 : GDALDataset::Open(osGeolocVRT.c_str(), GDAL_OF_RASTER));
2667 3 : if (poGeolocDS)
2668 : {
2669 6 : CPLStringList osMD(CSLDuplicate(poGeolocDS->GetMetadata()));
2670 3 : osMD.SetNameValue("X_DATASET", osGeolocVRT.c_str());
2671 3 : osMD.SetNameValue("Y_DATASET", osGeolocVRT.c_str());
2672 3 : const char *pszSRS = osMD.FetchNameValue("SRS");
2673 6 : OGRSpatialReference oSRS;
2674 3 : if (oSRS.SetFromUserInput(
2675 : pszSRS,
2676 3 : OGRSpatialReference::SET_FROM_USER_INPUT_LIMITATIONS_get()) ==
2677 : OGRERR_NONE)
2678 : {
2679 3 : osMD.SetNameValue("SRS", oSRS.exportToWkt().c_str());
2680 : }
2681 :
2682 3 : poDS->SetMetadata(osMD.List(), "GEOLOCATION");
2683 : }
2684 :
2685 3 : return poDS.release();
2686 : }
2687 :
2688 : /************************************************************************/
2689 : /* SENTINEL2GetGranuleList_L1CSafeCompact() */
2690 : /************************************************************************/
2691 :
2692 12 : static bool SENTINEL2GetGranuleList_L1CSafeCompact(
2693 : CPLXMLNode *psMainMTD, const char *pszFilename,
2694 : std::vector<L1CSafeCompatGranuleDescription> &osList)
2695 : {
2696 12 : CPLXMLNode *psProductInfo = CPLGetXMLNode(
2697 : psMainMTD, "=Level-1C_User_Product.General_Info.Product_Info");
2698 12 : if (psProductInfo == nullptr)
2699 : {
2700 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
2701 : "=Level-1C_User_Product.General_Info.Product_Info");
2702 0 : return false;
2703 : }
2704 :
2705 : CPLXMLNode *psProductOrganisation =
2706 12 : CPLGetXMLNode(psProductInfo, "Product_Organisation");
2707 12 : if (psProductOrganisation == nullptr)
2708 : {
2709 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
2710 : "Product_Organisation");
2711 0 : return false;
2712 : }
2713 :
2714 24 : CPLString osDirname(CPLGetDirnameSafe(pszFilename));
2715 : #if !defined(_WIN32)
2716 : char szPointerFilename[2048];
2717 : int nBytes = static_cast<int>(
2718 12 : readlink(pszFilename, szPointerFilename, sizeof(szPointerFilename)));
2719 12 : if (nBytes != -1)
2720 : {
2721 : const int nOffset =
2722 0 : std::min(nBytes, static_cast<int>(sizeof(szPointerFilename) - 1));
2723 0 : szPointerFilename[nOffset] = '\0';
2724 0 : osDirname = CPLGetDirnameSafe(szPointerFilename);
2725 : }
2726 : #endif
2727 :
2728 12 : const char chSeparator = SENTINEL2GetPathSeparator(osDirname);
2729 24 : for (CPLXMLNode *psIter = psProductOrganisation->psChild; psIter != nullptr;
2730 12 : psIter = psIter->psNext)
2731 : {
2732 12 : if (psIter->eType != CXT_Element ||
2733 12 : !EQUAL(psIter->pszValue, "Granule_List"))
2734 : {
2735 0 : continue;
2736 : }
2737 24 : for (CPLXMLNode *psIter2 = psIter->psChild; psIter2 != nullptr;
2738 12 : psIter2 = psIter2->psNext)
2739 : {
2740 12 : if (psIter2->eType != CXT_Element ||
2741 12 : !EQUAL(psIter2->pszValue, "Granule"))
2742 : {
2743 0 : continue;
2744 : }
2745 :
2746 : const char *pszImageFile =
2747 12 : CPLGetXMLValue(psIter2, "IMAGE_FILE", nullptr);
2748 12 : if (pszImageFile == nullptr || strlen(pszImageFile) < 3)
2749 : {
2750 0 : CPLDebug("SENTINEL2", "Missing IMAGE_FILE element");
2751 0 : continue;
2752 : }
2753 12 : L1CSafeCompatGranuleDescription oDesc;
2754 12 : oDesc.osBandPrefixPath = osDirname + chSeparator + pszImageFile;
2755 12 : if (CPLHasPathTraversal(oDesc.osBandPrefixPath.c_str()))
2756 : {
2757 0 : CPLError(CE_Failure, CPLE_NotSupported,
2758 : "Path traversal detected in %s",
2759 : oDesc.osBandPrefixPath.c_str());
2760 0 : return false;
2761 : }
2762 12 : oDesc.osBandPrefixPath.resize(oDesc.osBandPrefixPath.size() -
2763 : 3); // strip B12
2764 : // GRANULE/L1C_T30TXT_A007999_20170102T111441/IMG_DATA/T30TXT_20170102T111442_B12
2765 : // --> GRANULE/L1C_T30TXT_A007999_20170102T111441/MTD_TL.xml
2766 : oDesc.osMTDTLPath =
2767 24 : osDirname + chSeparator +
2768 48 : CPLGetDirnameSafe(CPLGetDirnameSafe(pszImageFile).c_str()) +
2769 12 : chSeparator + "MTD_TL.xml";
2770 12 : osList.push_back(std::move(oDesc));
2771 : }
2772 : }
2773 :
2774 12 : return true;
2775 : }
2776 :
2777 : /************************************************************************/
2778 : /* SENTINEL2GetGranuleList_L2ASafeCompact() */
2779 : /************************************************************************/
2780 :
2781 15 : static bool SENTINEL2GetGranuleList_L2ASafeCompact(
2782 : CPLXMLNode *psMainMTD, const char *pszFilename,
2783 : std::vector<L1CSafeCompatGranuleDescription> &osList)
2784 : {
2785 15 : const char *pszNodePath =
2786 : "=Level-2A_User_Product.General_Info.Product_Info";
2787 15 : CPLXMLNode *psProductInfo = CPLGetXMLNode(psMainMTD, pszNodePath);
2788 15 : if (psProductInfo == nullptr)
2789 : {
2790 6 : pszNodePath = "=Level-2A_User_Product.General_Info.L2A_Product_Info";
2791 6 : psProductInfo = CPLGetXMLNode(psMainMTD, pszNodePath);
2792 : }
2793 15 : if (psProductInfo == nullptr)
2794 : {
2795 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s", pszNodePath);
2796 0 : return false;
2797 : }
2798 :
2799 : CPLXMLNode *psProductOrganisation =
2800 15 : CPLGetXMLNode(psProductInfo, "Product_Organisation");
2801 15 : if (psProductOrganisation == nullptr)
2802 : {
2803 : psProductOrganisation =
2804 6 : CPLGetXMLNode(psProductInfo, "L2A_Product_Organisation");
2805 6 : if (psProductOrganisation == nullptr)
2806 : {
2807 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
2808 : "Product_Organisation");
2809 0 : return false;
2810 : }
2811 : }
2812 :
2813 30 : CPLString osDirname(CPLGetDirnameSafe(pszFilename));
2814 : #if !defined(_WIN32)
2815 : char szPointerFilename[2048];
2816 : int nBytes = static_cast<int>(
2817 15 : readlink(pszFilename, szPointerFilename, sizeof(szPointerFilename)));
2818 15 : if (nBytes != -1)
2819 : {
2820 : const int nOffset =
2821 0 : std::min(nBytes, static_cast<int>(sizeof(szPointerFilename) - 1));
2822 0 : szPointerFilename[nOffset] = '\0';
2823 0 : osDirname = CPLGetDirnameSafe(szPointerFilename);
2824 : }
2825 : #endif
2826 :
2827 15 : const char chSeparator = SENTINEL2GetPathSeparator(osDirname);
2828 42 : for (CPLXMLNode *psIter = psProductOrganisation->psChild; psIter != nullptr;
2829 27 : psIter = psIter->psNext)
2830 : {
2831 27 : if (psIter->eType != CXT_Element ||
2832 27 : !EQUAL(psIter->pszValue, "Granule_List"))
2833 : {
2834 0 : continue;
2835 : }
2836 54 : for (CPLXMLNode *psIter2 = psIter->psChild; psIter2 != nullptr;
2837 27 : psIter2 = psIter2->psNext)
2838 : {
2839 27 : if (psIter2->eType != CXT_Element ||
2840 27 : !EQUAL(psIter2->pszValue, "Granule"))
2841 : {
2842 0 : continue;
2843 : }
2844 :
2845 : const char *pszImageFile =
2846 27 : CPLGetXMLValue(psIter2, "IMAGE_FILE", nullptr);
2847 27 : if (pszImageFile == nullptr)
2848 : {
2849 : pszImageFile =
2850 18 : CPLGetXMLValue(psIter2, "IMAGE_FILE_2A", nullptr);
2851 18 : if (pszImageFile == nullptr || strlen(pszImageFile) < 3)
2852 : {
2853 0 : CPLDebug("SENTINEL2", "Missing IMAGE_FILE element");
2854 0 : continue;
2855 : }
2856 : }
2857 27 : L1CSafeCompatGranuleDescription oDesc;
2858 27 : oDesc.osBandPrefixPath = osDirname + chSeparator + pszImageFile;
2859 27 : if (oDesc.osBandPrefixPath.size() < 36)
2860 : {
2861 0 : CPLDebug("SENTINEL2", "Band prefix path too short");
2862 0 : continue;
2863 : }
2864 27 : if (CPLHasPathTraversal(oDesc.osBandPrefixPath.c_str()))
2865 : {
2866 0 : CPLError(CE_Failure, CPLE_NotSupported,
2867 : "Path traversal detected in %s",
2868 : oDesc.osBandPrefixPath.c_str());
2869 0 : return false;
2870 : }
2871 27 : oDesc.osBandPrefixPath.resize(oDesc.osBandPrefixPath.size() - 36);
2872 : // GRANULE/L1C_T30TXT_A007999_20170102T111441/IMG_DATA/T30TXT_20170102T111442_B12_60m
2873 : // --> GRANULE/L1C_T30TXT_A007999_20170102T111441/MTD_TL.xml
2874 : oDesc.osMTDTLPath =
2875 54 : osDirname + chSeparator +
2876 81 : CPLGetDirnameSafe(CPLGetDirnameSafe(pszImageFile).c_str());
2877 27 : if (oDesc.osMTDTLPath.size() < 9)
2878 : {
2879 0 : CPLDebug("SENTINEL2", "MTDTL path too short");
2880 0 : continue;
2881 : }
2882 27 : oDesc.osMTDTLPath.resize(oDesc.osMTDTLPath.size() - 9);
2883 27 : oDesc.osMTDTLPath = oDesc.osMTDTLPath + chSeparator + "MTD_TL.xml";
2884 27 : osList.push_back(std::move(oDesc));
2885 : }
2886 : }
2887 :
2888 15 : return true;
2889 : }
2890 :
2891 : /************************************************************************/
2892 : /* OpenL1C_L2A() */
2893 : /************************************************************************/
2894 :
2895 17 : GDALDataset *SENTINEL2Dataset::OpenL1C_L2A(const char *pszFilename,
2896 : SENTINEL2Level eLevel)
2897 : {
2898 17 : CPLXMLNode *psRoot = CPLParseXMLFile(pszFilename);
2899 17 : if (psRoot == nullptr)
2900 : {
2901 2 : CPLDebug("SENTINEL2", "Cannot XML parse %s", pszFilename);
2902 2 : return nullptr;
2903 : }
2904 :
2905 15 : char *pszOriginalXML = CPLSerializeXMLTree(psRoot);
2906 30 : CPLString osOriginalXML;
2907 15 : if (pszOriginalXML)
2908 15 : osOriginalXML = pszOriginalXML;
2909 15 : CPLFree(pszOriginalXML);
2910 :
2911 30 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRoot);
2912 15 : CPLStripXMLNamespace(psRoot, nullptr, TRUE);
2913 :
2914 15 : const char *pszNodePath =
2915 : (eLevel == SENTINEL2_L1C)
2916 : ? "=Level-1C_User_Product.General_Info.Product_Info"
2917 : : "=Level-2A_User_Product.General_Info.Product_Info";
2918 15 : CPLXMLNode *psProductInfo = CPLGetXMLNode(psRoot, pszNodePath);
2919 15 : if (psProductInfo == nullptr && eLevel == SENTINEL2_L2A)
2920 : {
2921 4 : pszNodePath = "=Level-2A_User_Product.General_Info.L2A_Product_Info";
2922 4 : psProductInfo = CPLGetXMLNode(psRoot, pszNodePath);
2923 : }
2924 15 : if (psProductInfo == nullptr)
2925 : {
2926 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s", pszNodePath);
2927 1 : return nullptr;
2928 : }
2929 :
2930 : const bool bIsSafeCompact =
2931 14 : EQUAL(CPLGetXMLValue(psProductInfo, "Query_Options.PRODUCT_FORMAT", ""),
2932 : "SAFE_COMPACT");
2933 :
2934 28 : std::set<int> oSetResolutions;
2935 28 : std::map<int, std::set<CPLString>> oMapResolutionsToBands;
2936 14 : if (bIsSafeCompact)
2937 : {
2938 70 : for (const auto &sBandDesc : asBandDesc)
2939 : {
2940 : // L2A does not contain B10
2941 65 : if (eLevel == SENTINEL2_L2A &&
2942 39 : strcmp(sBandDesc.pszBandName, "B10") == 0)
2943 3 : continue;
2944 62 : oSetResolutions.insert(sBandDesc.nResolution);
2945 124 : CPLString osName = sBandDesc.pszBandName + 1; /* skip B character */
2946 62 : if (atoi(osName) < 10)
2947 50 : osName = "0" + osName;
2948 62 : oMapResolutionsToBands[sBandDesc.nResolution].insert(
2949 62 : std::move(osName));
2950 : }
2951 5 : if (eLevel == SENTINEL2_L2A)
2952 : {
2953 39 : for (const auto &sL2ABandDesc : asL2ABandDesc)
2954 : {
2955 36 : oSetResolutions.insert(sL2ABandDesc.nResolution);
2956 36 : oMapResolutionsToBands[sL2ABandDesc.nResolution].insert(
2957 36 : sL2ABandDesc.pszBandName);
2958 : }
2959 : }
2960 : }
2961 15 : else if (eLevel == SENTINEL2_L1C &&
2962 6 : !SENTINEL2GetResolutionSet(psProductInfo, oSetResolutions,
2963 : oMapResolutionsToBands))
2964 : {
2965 3 : CPLDebug("SENTINEL2", "Failed to get resolution set");
2966 3 : return nullptr;
2967 : }
2968 :
2969 22 : std::vector<CPLString> aosGranuleList;
2970 11 : if (bIsSafeCompact)
2971 : {
2972 : std::vector<L1CSafeCompatGranuleDescription>
2973 5 : aoL1CSafeCompactGranuleList;
2974 7 : if (eLevel == SENTINEL2_L1C &&
2975 2 : !SENTINEL2GetGranuleList_L1CSafeCompact(
2976 : psRoot, pszFilename, aoL1CSafeCompactGranuleList))
2977 : {
2978 0 : CPLDebug("SENTINEL2", "Failed to get granule list");
2979 0 : return nullptr;
2980 : }
2981 8 : else if (eLevel == SENTINEL2_L2A &&
2982 3 : !SENTINEL2GetGranuleList_L2ASafeCompact(
2983 : psRoot, pszFilename, aoL1CSafeCompactGranuleList))
2984 : {
2985 0 : CPLDebug("SENTINEL2", "Failed to get granule list");
2986 0 : return nullptr;
2987 : }
2988 12 : for (size_t i = 0; i < aoL1CSafeCompactGranuleList.size(); ++i)
2989 : {
2990 7 : aosGranuleList.push_back(
2991 7 : aoL1CSafeCompactGranuleList[i].osMTDTLPath);
2992 : }
2993 : }
2994 6 : else if (!SENTINEL2GetGranuleList(
2995 : psRoot, eLevel, pszFilename, aosGranuleList,
2996 : (eLevel == SENTINEL2_L1C) ? nullptr : &oSetResolutions,
2997 : (eLevel == SENTINEL2_L1C) ? nullptr : &oMapResolutionsToBands))
2998 : {
2999 0 : CPLDebug("SENTINEL2", "Failed to get granule list");
3000 0 : return nullptr;
3001 : }
3002 11 : if (oSetResolutions.empty())
3003 : {
3004 0 : CPLDebug("SENTINEL2", "Resolution set is empty");
3005 0 : return nullptr;
3006 : }
3007 :
3008 11 : std::set<int> oSetEPSGCodes;
3009 27 : for (size_t i = 0; i < aosGranuleList.size(); i++)
3010 : {
3011 16 : int nEPSGCode = 0;
3012 16 : if (SENTINEL2GetGranuleInfo(eLevel, aosGranuleList[i],
3013 32 : *(oSetResolutions.begin()), &nEPSGCode))
3014 : {
3015 13 : oSetEPSGCodes.insert(nEPSGCode);
3016 : }
3017 : }
3018 :
3019 11 : SENTINEL2DatasetContainer *poDS = new SENTINEL2DatasetContainer();
3020 11 : char **papszMD = SENTINEL2GetUserProductMetadata(
3021 : psRoot, (eLevel == SENTINEL2_L1C) ? "Level-1C_User_Product"
3022 : : "Level-2A_User_Product");
3023 11 : poDS->GDALDataset::SetMetadata(papszMD);
3024 11 : CSLDestroy(papszMD);
3025 :
3026 11 : if (!osOriginalXML.empty())
3027 : {
3028 : char *apszXMLMD[2];
3029 11 : apszXMLMD[0] = const_cast<char *>(osOriginalXML.c_str());
3030 11 : apszXMLMD[1] = nullptr;
3031 11 : poDS->GDALDataset::SetMetadata(apszXMLMD, "xml:SENTINEL2");
3032 : }
3033 :
3034 11 : const char *pszPrefix =
3035 : (eLevel == SENTINEL2_L1C) ? "SENTINEL2_L1C" : "SENTINEL2_L2A";
3036 :
3037 : /* Create subdatsets per resolution (10, 20, 60m) and EPSG codes */
3038 11 : int iSubDSNum = 1;
3039 38 : for (std::set<int>::const_iterator oIterRes = oSetResolutions.begin();
3040 65 : oIterRes != oSetResolutions.end(); ++oIterRes)
3041 : {
3042 27 : const int nResolution = *oIterRes;
3043 :
3044 50 : for (std::set<int>::const_iterator oIterEPSG = oSetEPSGCodes.begin();
3045 73 : oIterEPSG != oSetEPSGCodes.end(); ++oIterEPSG)
3046 : {
3047 23 : const int nEPSGCode = *oIterEPSG;
3048 23 : poDS->GDALDataset::SetMetadataItem(
3049 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
3050 : CPLSPrintf("%s:%s:%dm:EPSG_%d", pszPrefix, pszFilename,
3051 : nResolution, nEPSGCode),
3052 : "SUBDATASETS");
3053 :
3054 : CPLString osBandNames = SENTINEL2GetBandListForResolution(
3055 46 : oMapResolutionsToBands[nResolution]);
3056 :
3057 : CPLString osDesc(CPLSPrintf("Bands %s with %dm resolution",
3058 23 : osBandNames.c_str(), nResolution));
3059 23 : if (nEPSGCode >= 32601 && nEPSGCode <= 32660)
3060 23 : osDesc += CPLSPrintf(", UTM %dN", nEPSGCode - 32600);
3061 0 : else if (nEPSGCode >= 32701 && nEPSGCode <= 32760)
3062 0 : osDesc += CPLSPrintf(", UTM %dS", nEPSGCode - 32700);
3063 : else
3064 0 : osDesc += CPLSPrintf(", EPSG:%d", nEPSGCode);
3065 23 : poDS->GDALDataset::SetMetadataItem(
3066 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
3067 : "SUBDATASETS");
3068 :
3069 23 : iSubDSNum++;
3070 : }
3071 : }
3072 :
3073 : /* Expose TCI or PREVIEW subdatasets */
3074 11 : if (bIsSafeCompact)
3075 : {
3076 10 : for (std::set<int>::const_iterator oIterEPSG = oSetEPSGCodes.begin();
3077 15 : oIterEPSG != oSetEPSGCodes.end(); ++oIterEPSG)
3078 : {
3079 5 : const int nEPSGCode = *oIterEPSG;
3080 5 : poDS->GDALDataset::SetMetadataItem(
3081 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
3082 : CPLSPrintf("%s:%s:TCI:EPSG_%d", pszPrefix, pszFilename,
3083 : nEPSGCode),
3084 : "SUBDATASETS");
3085 :
3086 5 : CPLString osDesc("True color image");
3087 5 : if (nEPSGCode >= 32601 && nEPSGCode <= 32660)
3088 5 : osDesc += CPLSPrintf(", UTM %dN", nEPSGCode - 32600);
3089 0 : else if (nEPSGCode >= 32701 && nEPSGCode <= 32760)
3090 0 : osDesc += CPLSPrintf(", UTM %dS", nEPSGCode - 32700);
3091 : else
3092 0 : osDesc += CPLSPrintf(", EPSG:%d", nEPSGCode);
3093 5 : poDS->GDALDataset::SetMetadataItem(
3094 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
3095 : "SUBDATASETS");
3096 :
3097 5 : iSubDSNum++;
3098 : }
3099 : }
3100 : else
3101 : {
3102 10 : for (std::set<int>::const_iterator oIterEPSG = oSetEPSGCodes.begin();
3103 14 : oIterEPSG != oSetEPSGCodes.end(); ++oIterEPSG)
3104 : {
3105 4 : const int nEPSGCode = *oIterEPSG;
3106 4 : poDS->GDALDataset::SetMetadataItem(
3107 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
3108 : CPLSPrintf("%s:%s:PREVIEW:EPSG_%d", pszPrefix, pszFilename,
3109 : nEPSGCode),
3110 : "SUBDATASETS");
3111 :
3112 4 : CPLString osDesc("RGB preview");
3113 4 : if (nEPSGCode >= 32601 && nEPSGCode <= 32660)
3114 4 : osDesc += CPLSPrintf(", UTM %dN", nEPSGCode - 32600);
3115 0 : else if (nEPSGCode >= 32701 && nEPSGCode <= 32760)
3116 0 : osDesc += CPLSPrintf(", UTM %dS", nEPSGCode - 32700);
3117 : else
3118 0 : osDesc += CPLSPrintf(", EPSG:%d", nEPSGCode);
3119 4 : poDS->GDALDataset::SetMetadataItem(
3120 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
3121 : "SUBDATASETS");
3122 :
3123 4 : iSubDSNum++;
3124 : }
3125 : }
3126 :
3127 11 : pszNodePath =
3128 : (eLevel == SENTINEL2_L1C)
3129 : ? "=Level-1C_User_Product.Geometric_Info.Product_Footprint."
3130 : "Product_Footprint.Global_Footprint.EXT_POS_LIST"
3131 : : "=Level-2A_User_Product.Geometric_Info.Product_Footprint."
3132 : "Product_Footprint.Global_Footprint.EXT_POS_LIST";
3133 11 : const char *pszPosList = CPLGetXMLValue(psRoot, pszNodePath, nullptr);
3134 11 : if (pszPosList != nullptr)
3135 : {
3136 18 : CPLString osPolygon = SENTINEL2GetPolygonWKTFromPosList(pszPosList);
3137 9 : if (!osPolygon.empty())
3138 9 : poDS->GDALDataset::SetMetadataItem("FOOTPRINT", osPolygon.c_str());
3139 : }
3140 :
3141 11 : return poDS;
3142 : }
3143 :
3144 : /************************************************************************/
3145 : /* SENTINEL2GetL1BCTileMetadata() */
3146 : /************************************************************************/
3147 :
3148 11 : static char **SENTINEL2GetL1BCTileMetadata(CPLXMLNode *psMainMTD)
3149 : {
3150 22 : CPLStringList aosList;
3151 :
3152 11 : CPLXMLNode *psRoot = CPLGetXMLNode(psMainMTD, "=Level-1C_Tile_ID");
3153 11 : if (psRoot == nullptr)
3154 : {
3155 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find =Level-1C_Tile_ID");
3156 0 : return nullptr;
3157 : }
3158 11 : CPLXMLNode *psGeneralInfo = CPLGetXMLNode(psRoot, "General_Info");
3159 11 : for (CPLXMLNode *psIter =
3160 11 : (psGeneralInfo ? psGeneralInfo->psChild : nullptr);
3161 58 : psIter != nullptr; psIter = psIter->psNext)
3162 : {
3163 47 : if (psIter->eType != CXT_Element)
3164 0 : continue;
3165 47 : const char *pszValue = CPLGetXMLValue(psIter, nullptr, nullptr);
3166 47 : if (pszValue != nullptr)
3167 : {
3168 38 : aosList.AddNameValue(psIter->pszValue, pszValue);
3169 : }
3170 : }
3171 :
3172 11 : CPLXMLNode *psQII = CPLGetXMLNode(psRoot, "Quality_Indicators_Info");
3173 11 : if (psQII != nullptr)
3174 : {
3175 9 : CPLXMLNode *psICCQI = CPLGetXMLNode(psQII, "Image_Content_QI");
3176 9 : for (CPLXMLNode *psIter = (psICCQI ? psICCQI->psChild : nullptr);
3177 27 : psIter != nullptr; psIter = psIter->psNext)
3178 : {
3179 18 : if (psIter->eType != CXT_Element)
3180 0 : continue;
3181 18 : if (psIter->psChild != nullptr &&
3182 18 : psIter->psChild->eType == CXT_Text)
3183 : {
3184 18 : aosList.AddNameValue(psIter->pszValue,
3185 18 : psIter->psChild->pszValue);
3186 : }
3187 : }
3188 : }
3189 :
3190 11 : return aosList.StealList();
3191 : }
3192 :
3193 : /************************************************************************/
3194 : /* OpenL1CTile() */
3195 : /************************************************************************/
3196 :
3197 14 : GDALDataset *SENTINEL2Dataset::OpenL1CTile(const char *pszFilename,
3198 : CPLXMLNode **ppsRootMainMTD,
3199 : int nResolutionOfInterest,
3200 : std::set<CPLString> *poBandSet)
3201 : {
3202 14 : CPLXMLNode *psRoot = CPLParseXMLFile(pszFilename);
3203 14 : if (psRoot == nullptr)
3204 : {
3205 3 : CPLDebug("SENTINEL2", "Cannot XML parse %s", pszFilename);
3206 3 : return nullptr;
3207 : }
3208 :
3209 11 : char *pszOriginalXML = CPLSerializeXMLTree(psRoot);
3210 22 : CPLString osOriginalXML;
3211 11 : if (pszOriginalXML)
3212 11 : osOriginalXML = pszOriginalXML;
3213 11 : CPLFree(pszOriginalXML);
3214 :
3215 22 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRoot);
3216 11 : CPLStripXMLNamespace(psRoot, nullptr, TRUE);
3217 :
3218 22 : std::set<int> oSetResolutions;
3219 22 : std::map<int, std::set<CPLString>> oMapResolutionsToBands;
3220 11 : char **papszMD = nullptr;
3221 11 : SENTINEL2GetResolutionSetAndMainMDFromGranule(
3222 : pszFilename, "Level-1C_User_Product", nResolutionOfInterest,
3223 : oSetResolutions, oMapResolutionsToBands, papszMD, ppsRootMainMTD);
3224 11 : if (poBandSet != nullptr)
3225 8 : *poBandSet = oMapResolutionsToBands[nResolutionOfInterest];
3226 :
3227 11 : SENTINEL2DatasetContainer *poDS = new SENTINEL2DatasetContainer();
3228 :
3229 11 : char **papszGranuleMD = SENTINEL2GetL1BCTileMetadata(psRoot);
3230 11 : papszMD = CSLMerge(papszMD, papszGranuleMD);
3231 11 : CSLDestroy(papszGranuleMD);
3232 :
3233 : // Remove CLOUD_COVERAGE_ASSESSMENT that comes from main metadata, if
3234 : // granule CLOUDY_PIXEL_PERCENTAGE is present.
3235 20 : if (CSLFetchNameValue(papszMD, "CLOUDY_PIXEL_PERCENTAGE") != nullptr &&
3236 9 : CSLFetchNameValue(papszMD, "CLOUD_COVERAGE_ASSESSMENT") != nullptr)
3237 : {
3238 7 : papszMD =
3239 7 : CSLSetNameValue(papszMD, "CLOUD_COVERAGE_ASSESSMENT", nullptr);
3240 : }
3241 :
3242 11 : poDS->GDALDataset::SetMetadata(papszMD);
3243 11 : CSLDestroy(papszMD);
3244 :
3245 11 : if (!osOriginalXML.empty())
3246 : {
3247 : char *apszXMLMD[2];
3248 11 : apszXMLMD[0] = const_cast<char *>(osOriginalXML.c_str());
3249 11 : apszXMLMD[1] = nullptr;
3250 11 : poDS->GDALDataset::SetMetadata(apszXMLMD, "xml:SENTINEL2");
3251 : }
3252 :
3253 : /* Create subdatsets per resolution (10, 20, 60m) */
3254 11 : int iSubDSNum = 1;
3255 36 : for (std::set<int>::const_iterator oIterRes = oSetResolutions.begin();
3256 61 : oIterRes != oSetResolutions.end(); ++oIterRes)
3257 : {
3258 25 : const int nResolution = *oIterRes;
3259 :
3260 25 : poDS->GDALDataset::SetMetadataItem(
3261 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
3262 : CPLSPrintf("%s:%s:%dm", "SENTINEL2_L1C_TILE", pszFilename,
3263 : nResolution),
3264 : "SUBDATASETS");
3265 :
3266 : CPLString osBandNames = SENTINEL2GetBandListForResolution(
3267 50 : oMapResolutionsToBands[nResolution]);
3268 :
3269 : CPLString osDesc(CPLSPrintf("Bands %s with %dm resolution",
3270 25 : osBandNames.c_str(), nResolution));
3271 25 : poDS->GDALDataset::SetMetadataItem(
3272 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
3273 : "SUBDATASETS");
3274 :
3275 25 : iSubDSNum++;
3276 : }
3277 :
3278 : /* Expose PREVIEW subdataset */
3279 11 : poDS->GDALDataset::SetMetadataItem(
3280 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
3281 : CPLSPrintf("%s:%s:PREVIEW", "SENTINEL2_L1C_TILE", pszFilename),
3282 : "SUBDATASETS");
3283 :
3284 11 : CPLString osDesc("RGB preview");
3285 11 : poDS->GDALDataset::SetMetadataItem(
3286 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
3287 : "SUBDATASETS");
3288 :
3289 11 : return poDS;
3290 : }
3291 :
3292 : /************************************************************************/
3293 : /* SENTINEL2GetOption() */
3294 : /************************************************************************/
3295 :
3296 55 : static const char *SENTINEL2GetOption(GDALOpenInfo *poOpenInfo,
3297 : const char *pszName,
3298 : const char *pszDefaultVal)
3299 : {
3300 : #ifdef GDAL_DMD_OPENOPTIONLIST
3301 : const char *pszVal =
3302 55 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, pszName);
3303 55 : if (pszVal != nullptr)
3304 3 : return pszVal;
3305 : #else
3306 : (void)poOpenInfo;
3307 : #endif
3308 52 : return CPLGetConfigOption(CPLSPrintf("SENTINEL2_%s", pszName),
3309 52 : pszDefaultVal);
3310 : }
3311 :
3312 : /************************************************************************/
3313 : /* LaunderUnit() */
3314 : /************************************************************************/
3315 :
3316 143 : static CPLString LaunderUnit(const char *pszUnit)
3317 : {
3318 143 : CPLString osUnit;
3319 1144 : for (int i = 0; pszUnit[i] != '\0';)
3320 : {
3321 1001 : if (strncmp(pszUnit + i,
3322 : "\xC2"
3323 : "\xB2",
3324 : 2) == 0) /* square / 2 */
3325 : {
3326 143 : i += 2;
3327 143 : osUnit += "2";
3328 : }
3329 858 : else if (strncmp(pszUnit + i,
3330 : "\xC2"
3331 : "\xB5",
3332 : 2) == 0) /* micro */
3333 : {
3334 143 : i += 2;
3335 143 : osUnit += "u";
3336 : }
3337 : else
3338 : {
3339 715 : osUnit += pszUnit[i];
3340 715 : i++;
3341 : }
3342 : }
3343 143 : return osUnit;
3344 : }
3345 :
3346 : /************************************************************************/
3347 : /* SENTINEL2GetTileInfo() */
3348 : /************************************************************************/
3349 :
3350 38 : static bool SENTINEL2GetTileInfo(const char *pszFilename, int *pnWidth,
3351 : int *pnHeight, int *pnBits)
3352 : {
3353 : static const unsigned char jp2_box_jp[] = {0x6a, 0x50, 0x20,
3354 : 0x20}; /* 'jP ' */
3355 38 : VSILFILE *fp = VSIFOpenL(pszFilename, "rb");
3356 38 : if (fp == nullptr)
3357 2 : return false;
3358 : GByte abyHeader[8];
3359 36 : if (VSIFReadL(abyHeader, 8, 1, fp) != 1)
3360 : {
3361 0 : VSIFCloseL(fp);
3362 0 : return false;
3363 : }
3364 36 : if (memcmp(abyHeader + 4, jp2_box_jp, 4) == 0)
3365 : {
3366 3 : bool bRet = false;
3367 : /* Just parse the ihdr box instead of doing a full dataset opening */
3368 3 : GDALJP2Box oBox(fp);
3369 3 : if (oBox.ReadFirst())
3370 : {
3371 9 : while (strlen(oBox.GetType()) > 0)
3372 : {
3373 9 : if (EQUAL(oBox.GetType(), "jp2h"))
3374 : {
3375 6 : GDALJP2Box oChildBox(fp);
3376 3 : if (!oChildBox.ReadFirstChild(&oBox))
3377 0 : break;
3378 :
3379 3 : while (strlen(oChildBox.GetType()) > 0)
3380 : {
3381 3 : if (EQUAL(oChildBox.GetType(), "ihdr"))
3382 : {
3383 3 : GByte *pabyData = oChildBox.ReadBoxData();
3384 3 : GIntBig nLength = oChildBox.GetDataLength();
3385 3 : if (pabyData != nullptr && nLength >= 4 + 4 + 2 + 1)
3386 : {
3387 3 : bRet = true;
3388 3 : if (pnHeight)
3389 : {
3390 1 : memcpy(pnHeight, pabyData, 4);
3391 1 : CPL_MSBPTR32(pnHeight);
3392 : }
3393 3 : if (pnWidth != nullptr)
3394 : {
3395 : // cppcheck-suppress nullPointer
3396 1 : memcpy(pnWidth, pabyData + 4, 4);
3397 1 : CPL_MSBPTR32(pnWidth);
3398 : }
3399 3 : if (pnBits)
3400 : {
3401 3 : GByte byPBC = pabyData[4 + 4 + 2];
3402 3 : if (byPBC != 255)
3403 : {
3404 3 : *pnBits = 1 + (byPBC & 0x7f);
3405 : }
3406 : else
3407 0 : *pnBits = 0;
3408 : }
3409 : }
3410 3 : CPLFree(pabyData);
3411 3 : break;
3412 : }
3413 0 : if (!oChildBox.ReadNextChild(&oBox))
3414 0 : break;
3415 : }
3416 3 : break;
3417 : }
3418 :
3419 6 : if (!oBox.ReadNext())
3420 0 : break;
3421 : }
3422 : }
3423 3 : VSIFCloseL(fp);
3424 3 : return bRet;
3425 : }
3426 : else /* for unit tests, we use TIFF */
3427 : {
3428 33 : VSIFCloseL(fp);
3429 : auto poDS = std::unique_ptr<GDALDataset>(GDALDataset::Open(
3430 33 : pszFilename, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR));
3431 33 : bool bRet = false;
3432 33 : if (poDS != nullptr)
3433 : {
3434 33 : if (poDS->GetRasterCount() != 0)
3435 : {
3436 33 : bRet = true;
3437 33 : if (pnWidth)
3438 3 : *pnWidth = poDS->GetRasterXSize();
3439 33 : if (pnHeight)
3440 3 : *pnHeight = poDS->GetRasterYSize();
3441 33 : if (pnBits)
3442 : {
3443 : const char *pszNBits =
3444 33 : poDS->GetRasterBand(1)->GetMetadataItem(
3445 33 : "NBITS", "IMAGE_STRUCTURE");
3446 33 : if (pszNBits == nullptr)
3447 : {
3448 : GDALDataType eDT =
3449 7 : poDS->GetRasterBand(1)->GetRasterDataType();
3450 : pszNBits =
3451 7 : CPLSPrintf("%d", GDALGetDataTypeSizeBits(eDT));
3452 : }
3453 33 : *pnBits = atoi(pszNBits);
3454 : }
3455 : }
3456 : }
3457 33 : return bRet;
3458 : }
3459 : }
3460 :
3461 : /************************************************************************/
3462 : /* OpenL1C_L2ASubdataset() */
3463 : /************************************************************************/
3464 :
3465 79 : GDALDataset *SENTINEL2Dataset::OpenL1C_L2ASubdataset(GDALOpenInfo *poOpenInfo,
3466 : SENTINEL2Level eLevel)
3467 : {
3468 158 : CPLString osFilename;
3469 79 : const char *pszPrefix =
3470 : (eLevel == SENTINEL2_L1C) ? "SENTINEL2_L1C" : "SENTINEL2_L2A";
3471 79 : CPLAssert(STARTS_WITH_CI(poOpenInfo->pszFilename, pszPrefix));
3472 79 : osFilename = poOpenInfo->pszFilename + strlen(pszPrefix) + 1;
3473 79 : const char *pszColumn = strrchr(osFilename.c_str(), ':');
3474 79 : if (pszColumn == nullptr || pszColumn == osFilename.c_str())
3475 : {
3476 10 : CPLError(CE_Failure, CPLE_AppDefined,
3477 : "Invalid syntax for %s:", pszPrefix);
3478 10 : return nullptr;
3479 : }
3480 69 : const auto nPos = static_cast<size_t>(pszColumn - osFilename.c_str());
3481 138 : const std::string osEPSGCode = osFilename.substr(nPos + 1);
3482 69 : osFilename.resize(nPos);
3483 69 : const char *pszPrecision = strrchr(osFilename.c_str(), ':');
3484 69 : if (pszPrecision == nullptr)
3485 : {
3486 10 : CPLError(CE_Failure, CPLE_AppDefined,
3487 : "Invalid syntax for %s:", pszPrefix);
3488 10 : return nullptr;
3489 : }
3490 :
3491 59 : if (!STARTS_WITH_CI(osEPSGCode.c_str(), "EPSG_"))
3492 : {
3493 5 : CPLError(CE_Failure, CPLE_AppDefined,
3494 : "Invalid syntax for %s:", pszPrefix);
3495 5 : return nullptr;
3496 : }
3497 :
3498 54 : const int nSubDSEPSGCode = atoi(osEPSGCode.c_str() + strlen("EPSG_"));
3499 54 : const bool bIsPreview = STARTS_WITH_CI(pszPrecision + 1, "PREVIEW");
3500 54 : const bool bIsTCI = STARTS_WITH_CI(pszPrecision + 1, "TCI");
3501 105 : const int nSubDSPrecision = (bIsPreview) ? 320
3502 51 : : (bIsTCI) ? RES_10M
3503 49 : : atoi(pszPrecision + 1);
3504 54 : if (!bIsTCI && !bIsPreview && nSubDSPrecision != RES_10M &&
3505 20 : nSubDSPrecision != RES_20M && nSubDSPrecision != RES_60M)
3506 : {
3507 5 : CPLError(CE_Failure, CPLE_NotSupported, "Unsupported precision: %d",
3508 : nSubDSPrecision);
3509 5 : return nullptr;
3510 : }
3511 :
3512 49 : osFilename.resize(pszPrecision - osFilename.c_str());
3513 98 : std::vector<CPLString> aosNonJP2Files;
3514 49 : aosNonJP2Files.push_back(osFilename);
3515 :
3516 49 : CPLXMLNode *psRoot = CPLParseXMLFile(osFilename);
3517 49 : if (psRoot == nullptr)
3518 : {
3519 7 : CPLDebug("SENTINEL2", "Cannot XML parse %s", osFilename.c_str());
3520 7 : return nullptr;
3521 : }
3522 :
3523 42 : char *pszOriginalXML = CPLSerializeXMLTree(psRoot);
3524 84 : CPLString osOriginalXML;
3525 42 : if (pszOriginalXML)
3526 42 : osOriginalXML = pszOriginalXML;
3527 42 : CPLFree(pszOriginalXML);
3528 :
3529 84 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRoot);
3530 42 : CPLStripXMLNamespace(psRoot, nullptr, TRUE);
3531 :
3532 : CPLXMLNode *psProductInfo =
3533 : eLevel == SENTINEL2_L1C
3534 42 : ? CPLGetXMLNode(psRoot,
3535 : "=Level-1C_User_Product.General_Info.Product_Info")
3536 18 : : CPLGetXMLNode(psRoot,
3537 42 : "=Level-2A_User_Product.General_Info.Product_Info");
3538 42 : if (psProductInfo == nullptr && eLevel == SENTINEL2_L2A)
3539 : {
3540 11 : psProductInfo = CPLGetXMLNode(
3541 : psRoot, "=Level-2A_User_Product.General_Info.L2A_Product_Info");
3542 : }
3543 42 : if (psProductInfo == nullptr)
3544 : {
3545 1 : CPLDebug("SENTINEL2", "Product Info not found");
3546 1 : return nullptr;
3547 : }
3548 :
3549 : const bool bIsSafeCompact =
3550 41 : EQUAL(CPLGetXMLValue(psProductInfo, "Query_Options.PRODUCT_FORMAT", ""),
3551 : "SAFE_COMPACT");
3552 :
3553 : const char *pszProductURI =
3554 41 : CPLGetXMLValue(psProductInfo, "PRODUCT_URI", nullptr);
3555 41 : SENTINEL2ProductType pType = MSI2A;
3556 41 : if (pszProductURI == nullptr)
3557 : {
3558 : pszProductURI =
3559 32 : CPLGetXMLValue(psProductInfo, "PRODUCT_URI_2A", nullptr);
3560 32 : pType = MSI2Ap;
3561 : }
3562 41 : if (pszProductURI == nullptr)
3563 27 : pszProductURI = "";
3564 :
3565 82 : std::vector<CPLString> aosGranuleList;
3566 82 : std::map<int, std::set<CPLString>> oMapResolutionsToBands;
3567 82 : std::vector<L1CSafeCompatGranuleDescription> aoL1CSafeCompactGranuleList;
3568 41 : if (bIsSafeCompact)
3569 : {
3570 308 : for (const auto &sBandDesc : asBandDesc)
3571 : {
3572 : // L2 does not contain B10
3573 286 : if (eLevel == SENTINEL2_L2A &&
3574 156 : strcmp(sBandDesc.pszBandName, "B10") == 0)
3575 12 : continue;
3576 548 : CPLString osName = sBandDesc.pszBandName + 1; /* skip B character */
3577 274 : if (atoi(osName) < 10)
3578 220 : osName = "0" + osName;
3579 274 : oMapResolutionsToBands[sBandDesc.nResolution].insert(
3580 274 : std::move(osName));
3581 : }
3582 22 : if (eLevel == SENTINEL2_L2A)
3583 : {
3584 156 : for (const auto &sL2ABandDesc : asL2ABandDesc)
3585 : {
3586 144 : oMapResolutionsToBands[sL2ABandDesc.nResolution].insert(
3587 144 : sL2ABandDesc.pszBandName);
3588 : }
3589 : }
3590 32 : if (eLevel == SENTINEL2_L1C &&
3591 10 : !SENTINEL2GetGranuleList_L1CSafeCompact(
3592 : psRoot, osFilename, aoL1CSafeCompactGranuleList))
3593 : {
3594 0 : CPLDebug("SENTINEL2", "Failed to get granule list");
3595 0 : return nullptr;
3596 : }
3597 34 : if (eLevel == SENTINEL2_L2A &&
3598 12 : !SENTINEL2GetGranuleList_L2ASafeCompact(
3599 : psRoot, osFilename, aoL1CSafeCompactGranuleList))
3600 : {
3601 0 : CPLDebug("SENTINEL2", "Failed to get granule list");
3602 0 : return nullptr;
3603 : }
3604 54 : for (size_t i = 0; i < aoL1CSafeCompactGranuleList.size(); ++i)
3605 : {
3606 32 : aosGranuleList.push_back(
3607 32 : aoL1CSafeCompactGranuleList[i].osMTDTLPath);
3608 : }
3609 : }
3610 19 : else if (!SENTINEL2GetGranuleList(
3611 : psRoot, eLevel, osFilename, aosGranuleList, nullptr,
3612 : (eLevel == SENTINEL2_L1C) ? nullptr : &oMapResolutionsToBands))
3613 : {
3614 1 : CPLDebug("SENTINEL2", "Failed to get granule list");
3615 1 : return nullptr;
3616 : }
3617 :
3618 80 : std::vector<CPLString> aosBands;
3619 80 : std::set<CPLString> oSetBands;
3620 40 : if (bIsPreview || bIsTCI)
3621 : {
3622 5 : aosBands.push_back("04");
3623 5 : aosBands.push_back("03");
3624 5 : aosBands.push_back("02");
3625 : }
3626 35 : else if (eLevel == SENTINEL2_L1C && !bIsSafeCompact)
3627 : {
3628 : CPLXMLNode *psBandList =
3629 10 : CPLGetXMLNode(psRoot, "=Level-1C_User_Product.General_Info.Product_"
3630 : "Info.Query_Options.Band_List");
3631 10 : if (psBandList == nullptr)
3632 : {
3633 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
3634 : "Query_Options.Band_List");
3635 0 : return nullptr;
3636 : }
3637 :
3638 116 : for (CPLXMLNode *psIter = psBandList->psChild; psIter != nullptr;
3639 106 : psIter = psIter->psNext)
3640 : {
3641 106 : if (psIter->eType != CXT_Element ||
3642 106 : !EQUAL(psIter->pszValue, "BAND_NAME"))
3643 72 : continue;
3644 106 : const char *pszBandName = CPLGetXMLValue(psIter, nullptr, "");
3645 : const SENTINEL2BandDescription *psBandDesc =
3646 106 : SENTINEL2GetBandDesc(pszBandName);
3647 106 : if (psBandDesc == nullptr)
3648 : {
3649 0 : CPLDebug("SENTINEL2", "Unknown band name %s", pszBandName);
3650 0 : continue;
3651 : }
3652 106 : if (psBandDesc->nResolution != nSubDSPrecision)
3653 72 : continue;
3654 : CPLString osName =
3655 68 : psBandDesc->pszBandName + 1; /* skip B character */
3656 34 : if (atoi(osName) < 10)
3657 30 : osName = "0" + osName;
3658 34 : oSetBands.insert(std::move(osName));
3659 : }
3660 10 : if (oSetBands.empty())
3661 : {
3662 0 : CPLDebug("SENTINEL2", "Band set is empty");
3663 0 : return nullptr;
3664 10 : }
3665 : }
3666 : else
3667 : {
3668 25 : oSetBands = oMapResolutionsToBands[nSubDSPrecision];
3669 : }
3670 :
3671 40 : if (aosBands.empty())
3672 : {
3673 192 : for (std::set<CPLString>::const_iterator oIterBandnames =
3674 35 : oSetBands.begin();
3675 419 : oIterBandnames != oSetBands.end(); ++oIterBandnames)
3676 : {
3677 192 : aosBands.push_back(*oIterBandnames);
3678 : }
3679 : /* Put 2=Blue, 3=Green, 4=Band bands in RGB order for conveniency */
3680 82 : if (aosBands.size() >= 3 && aosBands[0] == "02" &&
3681 82 : aosBands[1] == "03" && aosBands[2] == "04")
3682 : {
3683 17 : aosBands[0] = "04";
3684 17 : aosBands[2] = "02";
3685 : }
3686 : }
3687 :
3688 : /* -------------------------------------------------------------------- */
3689 : /* Create dataset. */
3690 : /* -------------------------------------------------------------------- */
3691 :
3692 40 : char **papszMD = SENTINEL2GetUserProductMetadata(
3693 : psRoot, (eLevel == SENTINEL2_L1C) ? "Level-1C_User_Product"
3694 : : "Level-2A_User_Product");
3695 :
3696 : const int nSaturatedVal =
3697 40 : atoi(CSLFetchNameValueDef(papszMD, "SPECIAL_VALUE_SATURATED", "-1"));
3698 : const int nNodataVal =
3699 40 : atoi(CSLFetchNameValueDef(papszMD, "SPECIAL_VALUE_NODATA", "-1"));
3700 :
3701 : const bool bAlpha =
3702 40 : CPLTestBool(SENTINEL2GetOption(poOpenInfo, "ALPHA", "FALSE"));
3703 :
3704 40 : SENTINEL2Dataset *poDS = CreateL1CL2ADataset(
3705 : eLevel, pType, bIsSafeCompact, aosGranuleList,
3706 : aoL1CSafeCompactGranuleList, aosNonJP2Files, nSubDSPrecision,
3707 : bIsPreview, bIsTCI, nSubDSEPSGCode, bAlpha, aosBands, nSaturatedVal,
3708 80 : nNodataVal, CPLString(pszProductURI));
3709 40 : if (poDS == nullptr)
3710 : {
3711 12 : CSLDestroy(papszMD);
3712 12 : return nullptr;
3713 : }
3714 :
3715 28 : if (!osOriginalXML.empty())
3716 : {
3717 28 : char *apszXMLMD[2] = {nullptr};
3718 28 : apszXMLMD[0] = const_cast<char *>(osOriginalXML.c_str());
3719 28 : apszXMLMD[1] = nullptr;
3720 28 : poDS->GDALDataset::SetMetadata(apszXMLMD, "xml:SENTINEL2");
3721 : }
3722 :
3723 28 : poDS->GDALDataset::SetMetadata(papszMD);
3724 28 : CSLDestroy(papszMD);
3725 :
3726 : /* -------------------------------------------------------------------- */
3727 : /* Add extra band metadata. */
3728 : /* -------------------------------------------------------------------- */
3729 28 : poDS->AddL1CL2ABandMetadata(eLevel, psRoot, aosBands);
3730 :
3731 : /* -------------------------------------------------------------------- */
3732 : /* Initialize overview information. */
3733 : /* -------------------------------------------------------------------- */
3734 28 : poDS->SetDescription(poOpenInfo->pszFilename);
3735 28 : CPLString osOverviewFile;
3736 28 : if (bIsPreview)
3737 : osOverviewFile = CPLSPrintf("%s_PREVIEW_EPSG_%d.tif.ovr",
3738 3 : osFilename.c_str(), nSubDSEPSGCode);
3739 25 : else if (bIsTCI)
3740 : osOverviewFile = CPLSPrintf("%s_TCI_EPSG_%d.tif.ovr",
3741 2 : osFilename.c_str(), nSubDSEPSGCode);
3742 : else
3743 : osOverviewFile =
3744 : CPLSPrintf("%s_%dm_EPSG_%d.tif.ovr", osFilename.c_str(),
3745 23 : nSubDSPrecision, nSubDSEPSGCode);
3746 28 : poDS->SetMetadataItem("OVERVIEW_FILE", osOverviewFile, "OVERVIEWS");
3747 28 : poDS->oOvManager.Initialize(poDS, ":::VIRTUAL:::");
3748 :
3749 28 : return poDS;
3750 : }
3751 :
3752 : /************************************************************************/
3753 : /* AddL1CL2ABandMetadata() */
3754 : /************************************************************************/
3755 :
3756 34 : void SENTINEL2Dataset::AddL1CL2ABandMetadata(
3757 : SENTINEL2Level eLevel, CPLXMLNode *psRoot,
3758 : const std::vector<CPLString> &aosBands)
3759 : {
3760 : CPLXMLNode *psIC =
3761 34 : CPLGetXMLNode(psRoot, (eLevel == SENTINEL2_L1C)
3762 : ? "=Level-1C_User_Product.General_Info."
3763 : "Product_Image_Characteristics"
3764 : : "=Level-2A_User_Product.General_Info."
3765 : "Product_Image_Characteristics");
3766 34 : if (psIC == nullptr)
3767 : {
3768 8 : psIC = CPLGetXMLNode(psRoot, "=Level-2A_User_Product.General_Info.L2A_"
3769 : "Product_Image_Characteristics");
3770 : }
3771 34 : if (psIC != nullptr)
3772 : {
3773 : CPLXMLNode *psSIL =
3774 32 : CPLGetXMLNode(psIC, "Reflectance_Conversion.Solar_Irradiance_List");
3775 32 : if (psSIL != nullptr)
3776 : {
3777 448 : for (CPLXMLNode *psIter = psSIL->psChild; psIter != nullptr;
3778 416 : psIter = psIter->psNext)
3779 : {
3780 416 : if (psIter->eType != CXT_Element ||
3781 416 : !EQUAL(psIter->pszValue, "SOLAR_IRRADIANCE"))
3782 : {
3783 0 : continue;
3784 : }
3785 : const char *pszBandId =
3786 416 : CPLGetXMLValue(psIter, "bandId", nullptr);
3787 416 : const char *pszUnit = CPLGetXMLValue(psIter, "unit", nullptr);
3788 416 : const char *pszValue = CPLGetXMLValue(psIter, nullptr, nullptr);
3789 416 : if (pszBandId != nullptr && pszUnit != nullptr &&
3790 : pszValue != nullptr)
3791 : {
3792 416 : const unsigned int nIdx = atoi(pszBandId);
3793 416 : if (nIdx < CPL_ARRAYSIZE(asBandDesc))
3794 : {
3795 2089 : for (int i = 0; i < nBands; i++)
3796 : {
3797 1816 : GDALRasterBand *poBand = GetRasterBand(i + 1);
3798 : const char *pszBandName =
3799 1816 : poBand->GetMetadataItem("BANDNAME");
3800 1816 : if (pszBandName != nullptr &&
3801 1806 : EQUAL(asBandDesc[nIdx].pszBandName,
3802 : pszBandName))
3803 : {
3804 143 : poBand->GDALRasterBand::SetMetadataItem(
3805 : "SOLAR_IRRADIANCE", pszValue);
3806 143 : poBand->GDALRasterBand::SetMetadataItem(
3807 : "SOLAR_IRRADIANCE_UNIT",
3808 286 : LaunderUnit(pszUnit));
3809 143 : break;
3810 : }
3811 : }
3812 : }
3813 : }
3814 : }
3815 : }
3816 :
3817 32 : CPLXMLNode *psOL = CPLGetXMLNode(
3818 : psIC, (eLevel == SENTINEL2_L1C) ? "Radiometric_Offset_List"
3819 : : "BOA_ADD_OFFSET_VALUES_LIST");
3820 32 : if (psOL != nullptr)
3821 : {
3822 56 : for (CPLXMLNode *psIter = psOL->psChild; psIter != nullptr;
3823 52 : psIter = psIter->psNext)
3824 : {
3825 52 : if (psIter->eType != CXT_Element ||
3826 52 : !EQUAL(psIter->pszValue, (eLevel == SENTINEL2_L1C)
3827 : ? "RADIO_ADD_OFFSET"
3828 : : "BOA_ADD_OFFSET"))
3829 : {
3830 0 : continue;
3831 : }
3832 : const char *pszBandId =
3833 52 : CPLGetXMLValue(psIter, "band_id", nullptr);
3834 52 : const char *pszValue = CPLGetXMLValue(psIter, nullptr, nullptr);
3835 52 : if (pszBandId != nullptr && pszValue != nullptr)
3836 : {
3837 52 : const unsigned int nIdx = atoi(pszBandId);
3838 52 : if (nIdx < CPL_ARRAYSIZE(asBandDesc))
3839 : {
3840 303 : for (int i = 0; i < nBands; i++)
3841 : {
3842 271 : GDALRasterBand *poBand = GetRasterBand(i + 1);
3843 : const char *pszBandName =
3844 271 : poBand->GetMetadataItem("BANDNAME");
3845 271 : if (pszBandName != nullptr &&
3846 271 : EQUAL(asBandDesc[nIdx].pszBandName,
3847 : pszBandName))
3848 : {
3849 20 : poBand->GDALRasterBand::SetMetadataItem(
3850 : (eLevel == SENTINEL2_L1C)
3851 : ? "RADIO_ADD_OFFSET"
3852 : : "BOA_ADD_OFFSET",
3853 : pszValue);
3854 20 : break;
3855 : }
3856 : }
3857 : }
3858 : }
3859 : }
3860 : }
3861 : }
3862 :
3863 : /* -------------------------------------------------------------------- */
3864 : /* Add Scene Classification category values (L2A) */
3865 : /* -------------------------------------------------------------------- */
3866 34 : CPLXMLNode *psSCL = CPLGetXMLNode(
3867 : psRoot, "=Level-2A_User_Product.General_Info."
3868 : "Product_Image_Characteristics.Scene_Classification_List");
3869 34 : if (psSCL == nullptr)
3870 : {
3871 29 : psSCL = CPLGetXMLNode(
3872 : psRoot,
3873 : "=Level-2A_User_Product.General_Info."
3874 : "L2A_Product_Image_Characteristics.L2A_Scene_Classification_List");
3875 : }
3876 34 : int nSCLBand = 0;
3877 199 : for (int nBand = 1; nBand <= static_cast<int>(aosBands.size()); nBand++)
3878 : {
3879 172 : if (EQUAL(aosBands[nBand - 1], "SCL"))
3880 : {
3881 7 : nSCLBand = nBand;
3882 7 : break;
3883 : }
3884 : }
3885 34 : if (psSCL != nullptr && nSCLBand > 0)
3886 : {
3887 14 : std::vector<std::string> aosCategories(100);
3888 7 : size_t nMaxIdx = 0;
3889 7 : bool bFirst = true;
3890 91 : for (CPLXMLNode *psIter = psSCL->psChild; psIter != nullptr;
3891 84 : psIter = psIter->psNext)
3892 : {
3893 84 : if (psIter->eType != CXT_Element ||
3894 84 : (!EQUAL(psIter->pszValue, "L2A_Scene_Classification_ID") &&
3895 36 : !EQUAL(psIter->pszValue, "Scene_Classification_ID")))
3896 : {
3897 0 : continue;
3898 : }
3899 : const char *pszText =
3900 84 : CPLGetXMLValue(psIter, "SCENE_CLASSIFICATION_TEXT", nullptr);
3901 84 : if (pszText == nullptr)
3902 48 : pszText = CPLGetXMLValue(
3903 : psIter, "L2A_SCENE_CLASSIFICATION_TEXT", nullptr);
3904 : const char *pszIdx =
3905 84 : CPLGetXMLValue(psIter, "SCENE_CLASSIFICATION_INDEX", nullptr);
3906 84 : if (pszIdx == nullptr)
3907 48 : pszIdx = CPLGetXMLValue(
3908 : psIter, "L2A_SCENE_CLASSIFICATION_INDEX", nullptr);
3909 84 : if (pszText && pszIdx)
3910 : {
3911 84 : const size_t nIdx = atoi(pszIdx);
3912 84 : if (nIdx < aosCategories.size())
3913 : {
3914 84 : if (bFirst)
3915 7 : nMaxIdx = nIdx;
3916 : else
3917 77 : nMaxIdx = std::max(nMaxIdx, nIdx);
3918 84 : bFirst = false;
3919 84 : if (STARTS_WITH_CI(pszText, "SC_"))
3920 84 : aosCategories[nIdx] = pszText + 3;
3921 : else
3922 0 : aosCategories[nIdx] = pszText;
3923 : }
3924 : }
3925 : }
3926 7 : if (!bFirst)
3927 : {
3928 7 : aosCategories.resize(nMaxIdx + 1);
3929 14 : GetRasterBand(nSCLBand)->SetCategoryNames(
3930 14 : CPLStringList(aosCategories).List());
3931 : }
3932 : }
3933 34 : }
3934 :
3935 : /************************************************************************/
3936 : /* CreateL1CL2ADataset() */
3937 : /************************************************************************/
3938 :
3939 48 : SENTINEL2Dataset *SENTINEL2Dataset::CreateL1CL2ADataset(
3940 : SENTINEL2Level eLevel, SENTINEL2ProductType pType, bool bIsSafeCompact,
3941 : const std::vector<CPLString> &aosGranuleList,
3942 : const std::vector<L1CSafeCompatGranuleDescription>
3943 : &aoL1CSafeCompactGranuleList,
3944 : std::vector<CPLString> &aosNonJP2Files, int nSubDSPrecision,
3945 : bool bIsPreview, bool bIsTCI,
3946 : int nSubDSEPSGCode, /* or -1 if not known at this point */
3947 : bool bAlpha, const std::vector<CPLString> &aosBands, int nSaturatedVal,
3948 : int nNodataVal, const CPLString &osProductURI)
3949 : {
3950 :
3951 : /* Iterate over granule metadata to know the layer extent */
3952 : /* and the location of each granule */
3953 48 : double dfMinX = 1.0e20;
3954 48 : double dfMinY = 1.0e20;
3955 48 : double dfMaxX = -1.0e20;
3956 48 : double dfMaxY = -1.0e20;
3957 96 : std::vector<SENTINEL2GranuleInfo> aosGranuleInfoList;
3958 48 : const int nDesiredResolution = (bIsPreview || bIsTCI) ? 0 : nSubDSPrecision;
3959 :
3960 48 : if (bIsSafeCompact)
3961 : {
3962 22 : CPLAssert(aosGranuleList.size() == aoL1CSafeCompactGranuleList.size());
3963 : }
3964 :
3965 116 : for (size_t i = 0; i < aosGranuleList.size(); i++)
3966 : {
3967 68 : int nEPSGCode = 0;
3968 68 : double dfULX = 0.0;
3969 68 : double dfULY = 0.0;
3970 68 : int nResolution = 0;
3971 68 : int nWidth = 0;
3972 68 : int nHeight = 0;
3973 68 : if (SENTINEL2GetGranuleInfo(eLevel, aosGranuleList[i],
3974 : nDesiredResolution, &nEPSGCode, &dfULX,
3975 65 : &dfULY, &nResolution, &nWidth, &nHeight) &&
3976 117 : (nSubDSEPSGCode == nEPSGCode || nSubDSEPSGCode < 0) &&
3977 49 : nResolution != 0)
3978 : {
3979 48 : nSubDSEPSGCode = nEPSGCode;
3980 48 : aosNonJP2Files.push_back(aosGranuleList[i]);
3981 :
3982 48 : if (dfULX < dfMinX)
3983 35 : dfMinX = dfULX;
3984 48 : if (dfULY > dfMaxY)
3985 35 : dfMaxY = dfULY;
3986 48 : const double dfLRX = dfULX + nResolution * nWidth;
3987 48 : const double dfLRY = dfULY - nResolution * nHeight;
3988 48 : if (dfLRX > dfMaxX)
3989 42 : dfMaxX = dfLRX;
3990 48 : if (dfLRY < dfMinY)
3991 42 : dfMinY = dfLRY;
3992 :
3993 96 : SENTINEL2GranuleInfo oGranuleInfo;
3994 48 : oGranuleInfo.osPath = CPLGetPathSafe(aosGranuleList[i]);
3995 48 : if (bIsSafeCompact)
3996 : {
3997 : oGranuleInfo.osBandPrefixPath =
3998 22 : aoL1CSafeCompactGranuleList[i].osBandPrefixPath;
3999 : }
4000 48 : oGranuleInfo.dfMinX = dfULX;
4001 48 : oGranuleInfo.dfMinY = dfLRY;
4002 48 : oGranuleInfo.dfMaxX = dfLRX;
4003 48 : oGranuleInfo.dfMaxY = dfULY;
4004 48 : oGranuleInfo.nWidth = nWidth / (nSubDSPrecision / nResolution);
4005 48 : oGranuleInfo.nHeight = nHeight / (nSubDSPrecision / nResolution);
4006 48 : aosGranuleInfoList.push_back(std::move(oGranuleInfo));
4007 : }
4008 : }
4009 48 : if (dfMinX > dfMaxX)
4010 : {
4011 13 : CPLError(CE_Failure, CPLE_NotSupported,
4012 : "No granule found for EPSG code %d", nSubDSEPSGCode);
4013 13 : return nullptr;
4014 : }
4015 :
4016 35 : const int nRasterXSize =
4017 35 : static_cast<int>((dfMaxX - dfMinX) / nSubDSPrecision + 0.5);
4018 35 : const int nRasterYSize =
4019 35 : static_cast<int>((dfMaxY - dfMinY) / nSubDSPrecision + 0.5);
4020 35 : SENTINEL2Dataset *poDS = new SENTINEL2Dataset(nRasterXSize, nRasterYSize);
4021 :
4022 : #if defined(__GNUC__)
4023 : #pragma GCC diagnostic push
4024 : #pragma GCC diagnostic ignored "-Wnull-dereference"
4025 : #endif
4026 35 : poDS->aosNonJP2Files = aosNonJP2Files;
4027 : #if defined(__GNUC__)
4028 : #pragma GCC diagnostic pop
4029 : #endif
4030 :
4031 35 : OGRSpatialReference oSRS;
4032 35 : char *pszProjection = nullptr;
4033 70 : if (oSRS.importFromEPSG(nSubDSEPSGCode) == OGRERR_NONE &&
4034 35 : oSRS.exportToWkt(&pszProjection) == OGRERR_NONE)
4035 : {
4036 35 : poDS->SetProjection(pszProjection);
4037 35 : CPLFree(pszProjection);
4038 : }
4039 : else
4040 : {
4041 0 : CPLDebug("SENTINEL2", "Invalid EPSG code %d", nSubDSEPSGCode);
4042 : }
4043 :
4044 35 : poDS->SetGeoTransform(GDALGeoTransform(dfMinX, nSubDSPrecision, 0, dfMaxY,
4045 35 : 0, -nSubDSPrecision));
4046 35 : poDS->GDALDataset::SetMetadataItem("COMPRESSION", "JPEG2000",
4047 : "IMAGE_STRUCTURE");
4048 35 : if (bIsPreview || bIsTCI)
4049 7 : poDS->GDALDataset::SetMetadataItem("INTERLEAVE", "PIXEL",
4050 : "IMAGE_STRUCTURE");
4051 :
4052 35 : int nBits = (bIsPreview || bIsTCI) ? 8 : 0 /* 0 = unknown yet*/;
4053 35 : int nValMax = (bIsPreview || bIsTCI) ? 255 : 0 /* 0 = unknown yet*/;
4054 : const int nBands =
4055 30 : (bIsPreview || bIsTCI)
4056 37 : ? 3
4057 28 : : ((bAlpha) ? 1 : 0) + static_cast<int>(aosBands.size());
4058 35 : const int nAlphaBand = (bIsPreview || bIsTCI || !bAlpha) ? 0 : nBands;
4059 35 : const GDALDataType eDT = (bIsPreview || bIsTCI) ? GDT_Byte : GDT_UInt16;
4060 :
4061 227 : for (int nBand = 1; nBand <= nBands; nBand++)
4062 : {
4063 192 : VRTSourcedRasterBand *poBand = nullptr;
4064 :
4065 192 : if (nBand != nAlphaBand)
4066 : {
4067 190 : poBand = new VRTSourcedRasterBand(
4068 190 : poDS, nBand, eDT, poDS->nRasterXSize, poDS->nRasterYSize);
4069 : }
4070 : else
4071 : {
4072 2 : poBand = new SENTINEL2AlphaBand(
4073 : poDS, nBand, eDT, poDS->nRasterXSize, poDS->nRasterYSize,
4074 2 : nSaturatedVal, nNodataVal);
4075 : }
4076 :
4077 192 : poDS->SetBand(nBand, poBand);
4078 192 : if (nBand == nAlphaBand)
4079 2 : poBand->SetColorInterpretation(GCI_AlphaBand);
4080 :
4081 384 : CPLString osBandName;
4082 192 : if (nBand != nAlphaBand)
4083 : {
4084 190 : osBandName = aosBands[nBand - 1];
4085 190 : SENTINEL2SetBandMetadata(poBand, osBandName);
4086 : }
4087 : else
4088 2 : osBandName = aosBands[0];
4089 :
4090 459 : for (size_t iSrc = 0; iSrc < aosGranuleInfoList.size(); iSrc++)
4091 : {
4092 267 : const SENTINEL2GranuleInfo &oGranuleInfo = aosGranuleInfoList[iSrc];
4093 267 : CPLString osTile;
4094 :
4095 267 : if (bIsSafeCompact && eLevel != SENTINEL2_L2A)
4096 : {
4097 33 : if (bIsTCI)
4098 : {
4099 6 : osTile = oGranuleInfo.osBandPrefixPath + "TCI.jp2";
4100 : }
4101 : else
4102 : {
4103 27 : osTile = oGranuleInfo.osBandPrefixPath + "B";
4104 27 : if (osBandName.size() == 1)
4105 0 : osTile += "0" + osBandName;
4106 27 : else if (osBandName.size() == 3)
4107 2 : osTile += osBandName.substr(1);
4108 : else
4109 25 : osTile += osBandName;
4110 27 : osTile += ".jp2";
4111 : }
4112 : }
4113 : else
4114 : {
4115 468 : osTile = SENTINEL2GetTilename(
4116 : oGranuleInfo.osPath, CPLGetFilename(oGranuleInfo.osPath),
4117 : osBandName, osProductURI, bIsPreview,
4118 234 : (eLevel == SENTINEL2_L1C) ? 0 : nSubDSPrecision);
4119 113 : if (bIsSafeCompact && eLevel == SENTINEL2_L2A &&
4120 419 : pType == MSI2Ap && osTile.size() >= 34 &&
4121 306 : osTile.substr(osTile.size() - 18, 3) != "MSK")
4122 : {
4123 60 : osTile.insert(osTile.size() - 34, "L2A_");
4124 : }
4125 234 : if (bIsTCI && osTile.size() >= 14)
4126 : {
4127 0 : osTile.replace(osTile.size() - 11, 3, "TCI");
4128 : }
4129 : }
4130 :
4131 267 : bool bTileFound = false;
4132 267 : if (nValMax == 0)
4133 : {
4134 : /* It is supposed to be 12 bits, but some products have 15 bits
4135 : */
4136 28 : if (SENTINEL2GetTileInfo(osTile, nullptr, nullptr, &nBits))
4137 : {
4138 27 : bTileFound = true;
4139 27 : if (nBits <= 16)
4140 27 : nValMax = (1 << nBits) - 1;
4141 : else
4142 : {
4143 0 : CPLDebug("SENTINEL2", "Unexpected bit depth %d", nBits);
4144 0 : nValMax = 65535;
4145 : }
4146 : }
4147 : }
4148 : else
4149 : {
4150 : VSIStatBufL sStat;
4151 239 : if (VSIStatExL(osTile, &sStat, VSI_STAT_EXISTS_FLAG) == 0)
4152 239 : bTileFound = true;
4153 : }
4154 267 : if (!bTileFound)
4155 : {
4156 1 : CPLError(CE_Warning, CPLE_AppDefined,
4157 : "Tile %s not found on filesystem. Skipping it",
4158 : osTile.c_str());
4159 1 : continue;
4160 : }
4161 :
4162 266 : const int nDstXOff = static_cast<int>(
4163 266 : (oGranuleInfo.dfMinX - dfMinX) / nSubDSPrecision + 0.5);
4164 266 : const int nDstYOff = static_cast<int>(
4165 266 : (dfMaxY - oGranuleInfo.dfMaxY) / nSubDSPrecision + 0.5);
4166 :
4167 266 : if (nBand != nAlphaBand)
4168 : {
4169 263 : poBand->AddSimpleSource(
4170 242 : osTile, (bIsPreview || bIsTCI) ? nBand : 1, 0, 0,
4171 263 : oGranuleInfo.nWidth, oGranuleInfo.nHeight, nDstXOff,
4172 263 : nDstYOff, oGranuleInfo.nWidth, oGranuleInfo.nHeight);
4173 : }
4174 : else
4175 : {
4176 3 : poBand->AddComplexSource(
4177 3 : osTile, 1, 0, 0, oGranuleInfo.nWidth, oGranuleInfo.nHeight,
4178 3 : nDstXOff, nDstYOff, oGranuleInfo.nWidth,
4179 3 : oGranuleInfo.nHeight, nValMax /* offset */, 0 /* scale */);
4180 : }
4181 : }
4182 :
4183 192 : if ((nBits % 8) != 0)
4184 : {
4185 143 : poBand->SetMetadataItem("NBITS", CPLSPrintf("%d", nBits),
4186 143 : "IMAGE_STRUCTURE");
4187 : }
4188 : }
4189 :
4190 35 : return poDS;
4191 : }
4192 :
4193 : /************************************************************************/
4194 : /* OpenL1CTileSubdataset() */
4195 : /************************************************************************/
4196 :
4197 13 : GDALDataset *SENTINEL2Dataset::OpenL1CTileSubdataset(GDALOpenInfo *poOpenInfo)
4198 : {
4199 26 : CPLString osFilename;
4200 13 : CPLAssert(STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1C_TILE:"));
4201 13 : osFilename = poOpenInfo->pszFilename + strlen("SENTINEL2_L1C_TILE:");
4202 13 : const char *pszPrecision = strrchr(osFilename.c_str(), ':');
4203 13 : if (pszPrecision == nullptr || pszPrecision == osFilename.c_str())
4204 : {
4205 2 : CPLError(CE_Failure, CPLE_AppDefined,
4206 : "Invalid syntax for SENTINEL2_L1C_TILE:");
4207 2 : return nullptr;
4208 : }
4209 11 : const bool bIsPreview = STARTS_WITH_CI(pszPrecision + 1, "PREVIEW");
4210 11 : const int nSubDSPrecision = (bIsPreview) ? 320 : atoi(pszPrecision + 1);
4211 11 : if (!bIsPreview && nSubDSPrecision != RES_10M &&
4212 2 : nSubDSPrecision != RES_20M && nSubDSPrecision != RES_60M)
4213 : {
4214 1 : CPLError(CE_Failure, CPLE_NotSupported, "Unsupported precision: %d",
4215 : nSubDSPrecision);
4216 1 : return nullptr;
4217 : }
4218 10 : osFilename.resize(pszPrecision - osFilename.c_str());
4219 :
4220 20 : std::set<CPLString> oSetBands;
4221 10 : CPLXMLNode *psRootMainMTD = nullptr;
4222 : GDALDataset *poTmpDS =
4223 10 : OpenL1CTile(osFilename, &psRootMainMTD, nSubDSPrecision, &oSetBands);
4224 20 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRootMainMTD);
4225 10 : if (poTmpDS == nullptr)
4226 2 : return nullptr;
4227 :
4228 16 : std::vector<CPLString> aosBands;
4229 8 : if (bIsPreview)
4230 : {
4231 2 : aosBands.push_back("04");
4232 2 : aosBands.push_back("03");
4233 2 : aosBands.push_back("02");
4234 : }
4235 : else
4236 : {
4237 21 : for (std::set<CPLString>::const_iterator oIterBandnames =
4238 6 : oSetBands.begin();
4239 48 : oIterBandnames != oSetBands.end(); ++oIterBandnames)
4240 : {
4241 21 : aosBands.push_back(*oIterBandnames);
4242 : }
4243 : /* Put 2=Blue, 3=Green, 4=Band bands in RGB order for conveniency */
4244 14 : if (aosBands.size() >= 3 && aosBands[0] == "02" &&
4245 14 : aosBands[1] == "03" && aosBands[2] == "04")
4246 : {
4247 3 : aosBands[0] = "04";
4248 3 : aosBands[2] = "02";
4249 : }
4250 : }
4251 :
4252 : /* -------------------------------------------------------------------- */
4253 : /* Create dataset. */
4254 : /* -------------------------------------------------------------------- */
4255 :
4256 16 : std::vector<CPLString> aosGranuleList;
4257 8 : aosGranuleList.push_back(osFilename);
4258 :
4259 16 : const int nSaturatedVal = atoi(CSLFetchNameValueDef(
4260 8 : poTmpDS->GetMetadata(), "SPECIAL_VALUE_SATURATED", "-1"));
4261 16 : const int nNodataVal = atoi(CSLFetchNameValueDef(
4262 8 : poTmpDS->GetMetadata(), "SPECIAL_VALUE_NODATA", "-1"));
4263 :
4264 : const bool bAlpha =
4265 8 : CPLTestBool(SENTINEL2GetOption(poOpenInfo, "ALPHA", "FALSE"));
4266 :
4267 16 : std::vector<CPLString> aosNonJP2Files;
4268 8 : SENTINEL2Dataset *poDS = CreateL1CL2ADataset(
4269 : SENTINEL2_L1C, MSI2A,
4270 : false, // bIsSafeCompact
4271 16 : aosGranuleList, std::vector<L1CSafeCompatGranuleDescription>(),
4272 : aosNonJP2Files, nSubDSPrecision, bIsPreview,
4273 : false, // bIsTCI
4274 : -1 /*nSubDSEPSGCode*/, bAlpha, aosBands, nSaturatedVal, nNodataVal,
4275 16 : CPLString());
4276 8 : if (poDS == nullptr)
4277 : {
4278 1 : delete poTmpDS;
4279 1 : return nullptr;
4280 : }
4281 :
4282 : // Transfer metadata
4283 7 : poDS->GDALDataset::SetMetadata(poTmpDS->GetMetadata());
4284 7 : poDS->GDALDataset::SetMetadata(poTmpDS->GetMetadata("xml:SENTINEL2"),
4285 : "xml:SENTINEL2");
4286 :
4287 7 : delete poTmpDS;
4288 :
4289 : /* -------------------------------------------------------------------- */
4290 : /* Add extra band metadata. */
4291 : /* -------------------------------------------------------------------- */
4292 7 : if (psRootMainMTD != nullptr)
4293 6 : poDS->AddL1CL2ABandMetadata(SENTINEL2_L1C, psRootMainMTD, aosBands);
4294 :
4295 : /* -------------------------------------------------------------------- */
4296 : /* Initialize overview information. */
4297 : /* -------------------------------------------------------------------- */
4298 7 : poDS->SetDescription(poOpenInfo->pszFilename);
4299 7 : CPLString osOverviewFile;
4300 7 : if (bIsPreview)
4301 2 : osOverviewFile = CPLSPrintf("%s_PREVIEW.tif.ovr", osFilename.c_str());
4302 : else
4303 : osOverviewFile =
4304 5 : CPLSPrintf("%s_%dm.tif.ovr", osFilename.c_str(), nSubDSPrecision);
4305 7 : poDS->SetMetadataItem("OVERVIEW_FILE", osOverviewFile, "OVERVIEWS");
4306 7 : poDS->oOvManager.Initialize(poDS, ":::VIRTUAL:::");
4307 :
4308 7 : return poDS;
4309 : }
4310 :
4311 : /************************************************************************/
4312 : /* GDALRegister_SENTINEL2() */
4313 : /************************************************************************/
4314 :
4315 2024 : void GDALRegister_SENTINEL2()
4316 : {
4317 2024 : if (GDALGetDriverByName("SENTINEL2") != nullptr)
4318 283 : return;
4319 :
4320 1741 : GDALDriver *poDriver = new GDALDriver();
4321 :
4322 1741 : poDriver->SetDescription("SENTINEL2");
4323 1741 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
4324 1741 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
4325 1741 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Sentinel 2");
4326 1741 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
4327 1741 : "drivers/raster/sentinel2.html");
4328 1741 : poDriver->SetMetadataItem(GDAL_DMD_SUBDATASETS, "YES");
4329 :
4330 1741 : poDriver->SetMetadataItem(
4331 : GDAL_DMD_OPENOPTIONLIST,
4332 : "<OpenOptionList>"
4333 : " <Option name='ALPHA' type='boolean' description='Whether to expose "
4334 : "an alpha band' default='NO'/>"
4335 : " <Option name='L1B_MODE' type='string-select' "
4336 : "default='DEFAULT'>\n"
4337 : " <Value>DEFAULT</Value>\n"
4338 : " <Value>GRANULE</Value>\n"
4339 : " <Value>DATASTRIP</Value>\n"
4340 : " </Option>\n"
4341 1741 : "</OpenOptionList>");
4342 :
4343 1741 : poDriver->pfnOpen = SENTINEL2Dataset::Open;
4344 1741 : poDriver->pfnIdentify = SENTINEL2Dataset::Identify;
4345 :
4346 1741 : GetGDALDriverManager()->RegisterDriver(poDriver);
4347 : }
|
