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