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 236756 : static bool IsS2Prefixed(const char *pszStr, const char *pszPrefixAfterS2X)
147 : {
148 282 : return pszStr[0] == 'S' && pszStr[1] == '2' && pszStr[2] >= 'A' &&
149 237212 : pszStr[2] <= 'Z' &&
150 174 : (*pszPrefixAfterS2X == 0 ||
151 236911 : 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 59378 : int SENTINEL2Dataset::Identify(GDALOpenInfo *poOpenInfo)
370 : {
371 59378 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1B:"))
372 36 : return TRUE;
373 59342 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1B_WITH_GEOLOC:"))
374 18 : return TRUE;
375 59324 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1C:"))
376 78 : return TRUE;
377 59246 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1C_TILE:"))
378 26 : return TRUE;
379 59220 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L2A:"))
380 80 : return TRUE;
381 :
382 59140 : 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 59142 : 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 59141 : if ((IsS2Prefixed(pszJustFilename, "_MSIL1C_") ||
393 118283 : IsS2Prefixed(pszJustFilename, "_MSIL2A_") ||
394 118283 : IsS2Prefixed(pszJustFilename, "_OPER_PRD_MSI") ||
395 177425 : IsS2Prefixed(pszJustFilename, "_USER_PRD_MSI")) &&
396 59140 : EQUAL(CPLGetExtensionSafe(pszJustFilename).c_str(), "zip"))
397 : {
398 0 : return TRUE;
399 : }
400 :
401 59142 : if (poOpenInfo->nHeaderBytes < 100)
402 55136 : return FALSE;
403 :
404 4006 : const char *pszHeader =
405 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader);
406 :
407 4006 : 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 3991 : 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 3983 : 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 3965 : 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 3957 : if (strstr(pszHeader, "<n1:Level-2A_User_Product") != nullptr &&
424 12 : strstr(pszHeader, "User_Product_Level-2A") != nullptr)
425 12 : return TRUE;
426 :
427 3945 : if (SENTINEL2isZipped(pszHeader, poOpenInfo->nHeaderBytes))
428 10 : return TRUE;
429 :
430 3934 : 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 190 : explicit SENTINEL2_CPLXMLNodeHolder(CPLXMLNode *psNode) : m_psNode(psNode)
445 : {
446 190 : }
447 :
448 190 : ~SENTINEL2_CPLXMLNodeHolder()
449 190 : {
450 190 : if (m_psNode)
451 174 : CPLDestroyXMLNode(m_psNode);
452 190 : }
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 3949 : static bool SENTINEL2isZipped(const char *pszHeader, int nHeaderBytes)
646 : {
647 3949 : 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 4011 : 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 3996 : 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 XML parse %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 416 : static char SENTINEL2GetPathSeparator(const char *pszBasename)
880 : {
881 416 : if (STARTS_WITH_CI(pszBasename, "\\\\?\\"))
882 0 : return '\\';
883 : else
884 416 : return '/';
885 : }
886 :
887 : /************************************************************************/
888 : /* SENTINEL2GetGranuleList() */
889 : /************************************************************************/
890 :
891 27 : 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 27 : 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 27 : CPLGetXMLNode(psMainMTD, CPLSPrintf("=%s", pszNodePath));
903 27 : if (psRoot == nullptr)
904 : {
905 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find =%s", pszNodePath);
906 0 : return false;
907 : }
908 27 : pszNodePath = "General_Info.Product_Info";
909 27 : CPLXMLNode *psProductInfo = CPLGetXMLNode(psRoot, pszNodePath);
910 27 : if (psProductInfo == nullptr && eLevel == SENTINEL2_L2A)
911 : {
912 9 : pszNodePath = "General_Info.L2A_Product_Info";
913 9 : psProductInfo = CPLGetXMLNode(psRoot, pszNodePath);
914 : }
915 27 : if (psProductInfo == nullptr)
916 : {
917 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s", pszNodePath);
918 0 : return false;
919 : }
920 :
921 27 : pszNodePath = "Product_Organisation";
922 : CPLXMLNode *psProductOrganisation =
923 27 : CPLGetXMLNode(psProductInfo, pszNodePath);
924 27 : if (psProductOrganisation == nullptr && eLevel == SENTINEL2_L2A)
925 : {
926 9 : pszNodePath = "L2A_Product_Organisation";
927 9 : psProductOrganisation = CPLGetXMLNode(psProductInfo, pszNodePath);
928 : }
929 27 : if (psProductOrganisation == nullptr)
930 : {
931 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s", pszNodePath);
932 1 : return false;
933 : }
934 :
935 52 : CPLString osDirname(CPLGetDirnameSafe(pszFilename));
936 : #if !defined(_WIN32)
937 : char szPointerFilename[2048];
938 : int nBytes = static_cast<int>(
939 26 : readlink(pszFilename, szPointerFilename, sizeof(szPointerFilename)));
940 26 : 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 26 : EQUAL(CPLGetXMLValue(psProductInfo, "PRODUCT_TYPE", ""), "S2MSI2Ap");
951 26 : const bool bIsCompact = EQUAL(
952 : CPLGetXMLValue(psProductInfo, "PRODUCT_FORMAT", ""), "SAFE_COMPACT");
953 52 : CPLString oGranuleId("L2A_");
954 52 : std::set<CPLString> aoSetGranuleId;
955 73 : for (CPLXMLNode *psIter = psProductOrganisation->psChild; psIter != nullptr;
956 47 : psIter = psIter->psNext)
957 : {
958 47 : if (psIter->eType != CXT_Element ||
959 47 : !EQUAL(psIter->pszValue, "Granule_List"))
960 : {
961 0 : continue;
962 : }
963 102 : for (CPLXMLNode *psIter2 = psIter->psChild; psIter2 != nullptr;
964 55 : psIter2 = psIter2->psNext)
965 : {
966 55 : if (psIter2->eType != CXT_Element ||
967 47 : (!EQUAL(psIter2->pszValue, "Granule") &&
968 47 : !EQUAL(psIter2->pszValue, "Granules")))
969 : {
970 16 : continue;
971 : }
972 : const char *pszGranuleId =
973 47 : CPLGetXMLValue(psIter2, "granuleIdentifier", nullptr);
974 47 : if (pszGranuleId == nullptr)
975 : {
976 4 : CPLDebug("SENTINEL2", "Missing granuleIdentifier attribute");
977 4 : continue;
978 : }
979 :
980 43 : 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 43 : if (aoSetGranuleId.find(pszGranuleId) != aoSetGranuleId.end())
1027 0 : continue;
1028 43 : 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 43 : CPLString osGranuleMTD = pszGranuleId;
1035 43 : if (bIsCompact == 0 &&
1036 43 : osGranuleMTD.size() > strlen("S2A_OPER_MSI_") &&
1037 39 : osGranuleMTD[8] == '_' && osGranuleMTD[12] == '_' &&
1038 39 : osGranuleMTD[osGranuleMTD.size() - 7] == '_' &&
1039 125 : osGranuleMTD[osGranuleMTD.size() - 6] == 'N' &&
1040 39 : osGranuleMTD[7] == 'R')
1041 : {
1042 39 : osGranuleMTD[9] = 'M';
1043 39 : osGranuleMTD[10] = 'T';
1044 39 : osGranuleMTD[11] = 'D';
1045 39 : 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 39 : osGranuleMTD += ".xml";
1076 :
1077 39 : const char chSeparator = SENTINEL2GetPathSeparator(osDirname);
1078 39 : CPLString osGranuleMTDPath = osDirname;
1079 39 : osGranuleMTDPath += chSeparator;
1080 39 : osGranuleMTDPath += "GRANULE";
1081 39 : osGranuleMTDPath += chSeparator;
1082 39 : osGranuleMTDPath += pszGranuleId;
1083 39 : osGranuleMTDPath += chSeparator;
1084 39 : osGranuleMTDPath += osGranuleMTD;
1085 39 : 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 39 : osList.push_back(std::move(osGranuleMTDPath));
1093 : }
1094 : }
1095 :
1096 26 : 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 XML parse %s", poOpenInfo->pszFilename);
1434 1 : return nullptr;
1435 : }
1436 :
1437 7 : char *pszOriginalXML = CPLSerializeXMLTree(psRoot);
1438 14 : CPLString osOriginalXML;
1439 7 : if (pszOriginalXML)
1440 7 : osOriginalXML = pszOriginalXML;
1441 7 : CPLFree(pszOriginalXML);
1442 :
1443 14 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRoot);
1444 7 : CPLStripXMLNamespace(psRoot, nullptr, TRUE);
1445 :
1446 7 : CPLXMLNode *psProductInfo = CPLGetXMLNode(
1447 : psRoot, "=Level-1B_User_Product.General_Info.Product_Info");
1448 7 : if (psProductInfo == nullptr)
1449 : {
1450 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
1451 : "=Level-1B_User_Product.General_Info.Product_Info");
1452 1 : return nullptr;
1453 : }
1454 :
1455 12 : auto poDS = std::make_unique<SENTINEL2DatasetContainer>();
1456 : char **papszMD =
1457 6 : SENTINEL2GetUserProductMetadata(psRoot, "Level-1B_User_Product");
1458 6 : poDS->GDALDataset::SetMetadata(papszMD);
1459 6 : CSLDestroy(papszMD);
1460 :
1461 6 : if (!osOriginalXML.empty())
1462 : {
1463 6 : char *apszXMLMD[2] = {nullptr};
1464 6 : apszXMLMD[0] = const_cast<char *>(osOriginalXML.c_str());
1465 6 : apszXMLMD[1] = nullptr;
1466 6 : poDS->GDALDataset::SetMetadata(apszXMLMD, "xml:SENTINEL2");
1467 : }
1468 :
1469 6 : const char *pszPosList = CPLGetXMLValue(
1470 : psRoot,
1471 : "=Level-1B_User_Product.Geometric_Info.Product_Footprint."
1472 : "Product_Footprint.Global_Footprint.EXT_POS_LIST",
1473 : nullptr);
1474 6 : if (pszPosList != nullptr)
1475 : {
1476 6 : CPLString osPolygon = SENTINEL2GetPolygonWKTFromPosList(pszPosList);
1477 3 : if (!osPolygon.empty())
1478 3 : poDS->GDALDataset::SetMetadataItem("FOOTPRINT", osPolygon.c_str());
1479 : }
1480 :
1481 6 : const char *pszMode = CSLFetchNameValueDef(poOpenInfo->papszOpenOptions,
1482 : "L1B_MODE", "DEFAULT");
1483 :
1484 : // Detect if this L1B product has associated geolocation arrays, as provided
1485 : // by the Sen2VM MPC project.
1486 : const std::string osDatastripRoot = CPLFormFilenameSafe(
1487 12 : CPLGetPathSafe(poOpenInfo->pszFilename).c_str(), "DATASTRIP", nullptr);
1488 : VSIStatBufL sStat;
1489 6 : if (VSIStatL(osDatastripRoot.c_str(), &sStat) == 0 &&
1490 6 : VSI_ISDIR(sStat.st_mode) && !EQUAL(pszMode, "GRANULE"))
1491 : {
1492 1 : int iSubDSNum = 1;
1493 1 : const CPLStringList aosList(VSIReadDir(osDatastripRoot.c_str()));
1494 4 : for (const char *pszDatastripId : aosList)
1495 : {
1496 3 : if (IsS2Prefixed(pszDatastripId, "_OPER_MSI_L1B_"))
1497 : {
1498 : const std::string osGEO_DATA = CPLFormFilenameSafe(
1499 1 : CPLFormFilenameSafe(osDatastripRoot.c_str(), pszDatastripId,
1500 : nullptr)
1501 : .c_str(),
1502 2 : "GEO_DATA", nullptr);
1503 2 : const CPLStringList aosListVRT(VSIReadDir(osGEO_DATA.c_str()));
1504 2 : std::vector<std::string> aosVRTs;
1505 6 : for (const char *pszVRT : aosListVRT)
1506 : {
1507 5 : if (EQUAL(CPLGetExtensionSafe(pszVRT).c_str(), "VRT"))
1508 : {
1509 3 : aosVRTs.push_back(pszVRT);
1510 : }
1511 : }
1512 1 : std::sort(aosVRTs.begin(), aosVRTs.end());
1513 4 : for (const std::string &osVRT : aosVRTs)
1514 : {
1515 6 : poDS->GDALDataset::SetMetadataItem(
1516 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
1517 : CPLSPrintf("SENTINEL2_L1B_WITH_GEOLOC:\"%s\":%s",
1518 : poOpenInfo->pszFilename,
1519 6 : CPLGetBasenameSafe(osVRT.c_str()).c_str()),
1520 : "SUBDATASETS");
1521 3 : ++iSubDSNum;
1522 : }
1523 : }
1524 : }
1525 1 : if (iSubDSNum > 1)
1526 : {
1527 1 : return poDS.release();
1528 : }
1529 0 : if (EQUAL(pszMode, "DATASTRIP"))
1530 : {
1531 0 : CPLError(CE_Failure, CPLE_AppDefined,
1532 : "Could not find VRT geolocation files");
1533 0 : return nullptr;
1534 : }
1535 : }
1536 :
1537 10 : std::set<int> oSetResolutions;
1538 10 : std::map<int, std::set<CPLString>> oMapResolutionsToBands;
1539 5 : if (!SENTINEL2GetResolutionSet(psProductInfo, oSetResolutions,
1540 : oMapResolutionsToBands))
1541 : {
1542 3 : CPLDebug("SENTINEL2", "Failed to get resolution set");
1543 3 : return nullptr;
1544 : }
1545 :
1546 4 : std::vector<CPLString> aosGranuleList;
1547 2 : if (!SENTINEL2GetGranuleList(psRoot, SENTINEL2_L1B, poOpenInfo->pszFilename,
1548 : aosGranuleList))
1549 : {
1550 0 : CPLDebug("SENTINEL2", "Failed to get granule list");
1551 0 : return nullptr;
1552 : }
1553 :
1554 : /* Create subdatsets per granules and resolution (10, 20, 60m) */
1555 2 : int iSubDSNum = 1;
1556 4 : for (size_t i = 0; i < aosGranuleList.size(); i++)
1557 : {
1558 8 : for (std::set<int>::const_iterator oIterRes = oSetResolutions.begin();
1559 14 : oIterRes != oSetResolutions.end(); ++oIterRes)
1560 : {
1561 6 : const int nResolution = *oIterRes;
1562 :
1563 12 : poDS->GDALDataset::SetMetadataItem(
1564 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
1565 6 : CPLSPrintf("SENTINEL2_L1B:%s:%dm", aosGranuleList[i].c_str(),
1566 : nResolution),
1567 : "SUBDATASETS");
1568 :
1569 : CPLString osBandNames = SENTINEL2GetBandListForResolution(
1570 12 : oMapResolutionsToBands[nResolution]);
1571 :
1572 : CPLString osDesc(
1573 : CPLSPrintf("Bands %s of granule %s with %dm resolution",
1574 : osBandNames.c_str(),
1575 6 : CPLGetFilename(aosGranuleList[i]), nResolution));
1576 6 : poDS->GDALDataset::SetMetadataItem(
1577 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
1578 : "SUBDATASETS");
1579 :
1580 6 : iSubDSNum++;
1581 : }
1582 : }
1583 :
1584 2 : return poDS.release();
1585 : }
1586 :
1587 : /************************************************************************/
1588 : /* SENTINEL2GetL1BGranuleMetadata() */
1589 : /************************************************************************/
1590 :
1591 15 : static char **SENTINEL2GetL1BGranuleMetadata(CPLXMLNode *psMainMTD)
1592 : {
1593 30 : CPLStringList aosList;
1594 :
1595 15 : CPLXMLNode *psRoot = CPLGetXMLNode(psMainMTD, "=Level-1B_Granule_ID");
1596 15 : if (psRoot == nullptr)
1597 : {
1598 1 : CPLError(CE_Failure, CPLE_AppDefined,
1599 : "Cannot find =Level-1B_Granule_ID");
1600 1 : return nullptr;
1601 : }
1602 14 : CPLXMLNode *psGeneralInfo = CPLGetXMLNode(psRoot, "General_Info");
1603 14 : for (CPLXMLNode *psIter =
1604 14 : (psGeneralInfo ? psGeneralInfo->psChild : nullptr);
1605 54 : psIter != nullptr; psIter = psIter->psNext)
1606 : {
1607 40 : if (psIter->eType != CXT_Element)
1608 0 : continue;
1609 40 : const char *pszValue = CPLGetXMLValue(psIter, nullptr, nullptr);
1610 40 : if (pszValue != nullptr)
1611 : {
1612 34 : aosList.AddNameValue(psIter->pszValue, pszValue);
1613 : }
1614 : }
1615 :
1616 14 : CPLXMLNode *psGeometryHeader = CPLGetXMLNode(
1617 : psRoot, "Geometric_Info.Granule_Position.Geometric_Header");
1618 14 : if (psGeometryHeader != nullptr)
1619 : {
1620 6 : const char *pszVal = CPLGetXMLValue(
1621 : psGeometryHeader, "Incidence_Angles.ZENITH_ANGLE", nullptr);
1622 6 : if (pszVal)
1623 6 : aosList.AddNameValue("INCIDENCE_ZENITH_ANGLE", pszVal);
1624 :
1625 6 : pszVal = CPLGetXMLValue(psGeometryHeader,
1626 : "Incidence_Angles.AZIMUTH_ANGLE", nullptr);
1627 6 : if (pszVal)
1628 6 : aosList.AddNameValue("INCIDENCE_AZIMUTH_ANGLE", pszVal);
1629 :
1630 6 : pszVal = CPLGetXMLValue(psGeometryHeader, "Solar_Angles.ZENITH_ANGLE",
1631 : nullptr);
1632 6 : if (pszVal)
1633 6 : aosList.AddNameValue("SOLAR_ZENITH_ANGLE", pszVal);
1634 :
1635 6 : pszVal = CPLGetXMLValue(psGeometryHeader, "Solar_Angles.AZIMUTH_ANGLE",
1636 : nullptr);
1637 6 : if (pszVal)
1638 6 : aosList.AddNameValue("SOLAR_AZIMUTH_ANGLE", pszVal);
1639 : }
1640 :
1641 14 : CPLXMLNode *psQII = CPLGetXMLNode(psRoot, "Quality_Indicators_Info");
1642 14 : if (psQII != nullptr)
1643 : {
1644 6 : CPLXMLNode *psICCQI = CPLGetXMLNode(psQII, "Image_Content_QI");
1645 6 : for (CPLXMLNode *psIter = (psICCQI ? psICCQI->psChild : nullptr);
1646 18 : psIter != nullptr; psIter = psIter->psNext)
1647 : {
1648 12 : if (psIter->eType != CXT_Element)
1649 0 : continue;
1650 12 : if (psIter->psChild != nullptr &&
1651 12 : psIter->psChild->eType == CXT_Text)
1652 : {
1653 12 : aosList.AddNameValue(psIter->pszValue,
1654 12 : psIter->psChild->pszValue);
1655 : }
1656 : }
1657 : }
1658 :
1659 14 : return aosList.StealList();
1660 : }
1661 :
1662 : /************************************************************************/
1663 : /* SENTINEL2GetTilename() */
1664 : /************************************************************************/
1665 :
1666 : static CPLString
1667 343 : SENTINEL2GetTilename(const std::string &osGranulePath,
1668 : const std::string &osGranuleName,
1669 : const std::string &osBandName,
1670 : const std::string &osProductURI = CPLString(),
1671 : bool bIsPreview = false, int nPrecisionL2A = 0)
1672 : {
1673 343 : bool granuleNameMatchTilename = true;
1674 686 : CPLString osJPEG2000Name(osGranuleName);
1675 343 : if (osJPEG2000Name.size() > 7 &&
1676 496 : osJPEG2000Name[osJPEG2000Name.size() - 7] == '_' &&
1677 153 : osJPEG2000Name[osJPEG2000Name.size() - 6] == 'N')
1678 : {
1679 136 : osJPEG2000Name.resize(osJPEG2000Name.size() - 7);
1680 : }
1681 :
1682 : const SENTINEL2_L2A_BandDescription *psL2ABandDesc =
1683 343 : (nPrecisionL2A) ? SENTINEL2GetL2ABandDesc(osBandName.c_str()) : nullptr;
1684 :
1685 343 : CPLString osTile(osGranulePath);
1686 343 : const char chSeparator = SENTINEL2GetPathSeparator(osTile);
1687 343 : if (!osTile.empty())
1688 343 : osTile += chSeparator;
1689 343 : bool procBaseLineIs1 = false;
1690 553 : if (osJPEG2000Name.size() > 12 && osJPEG2000Name[8] == '_' &&
1691 210 : osJPEG2000Name[12] == '_')
1692 210 : procBaseLineIs1 = true;
1693 343 : if (bIsPreview ||
1694 65 : (psL2ABandDesc != nullptr && (psL2ABandDesc->eLocation == TL_QI_DATA)))
1695 : {
1696 43 : osTile += "QI_DATA";
1697 43 : osTile += chSeparator;
1698 43 : if (procBaseLineIs1)
1699 : {
1700 25 : if (atoi(osBandName.c_str()) > 0)
1701 : {
1702 21 : osJPEG2000Name[9] = 'P';
1703 21 : osJPEG2000Name[10] = 'V';
1704 21 : osJPEG2000Name[11] = 'I';
1705 : }
1706 4 : else if (nPrecisionL2A && osBandName.size() == 3)
1707 : {
1708 4 : osJPEG2000Name[9] = osBandName[0];
1709 4 : osJPEG2000Name[10] = osBandName[1];
1710 4 : osJPEG2000Name[11] = osBandName[2];
1711 : }
1712 25 : osTile += osJPEG2000Name;
1713 : }
1714 : else
1715 : {
1716 18 : osTile += "MSK_";
1717 18 : osTile += osBandName;
1718 18 : osTile += "PRB";
1719 : }
1720 43 : if (nPrecisionL2A && !bIsPreview)
1721 22 : osTile += CPLSPrintf("_%02dm", nPrecisionL2A);
1722 : }
1723 : else
1724 : {
1725 300 : osTile += "IMG_DATA";
1726 300 : osTile += chSeparator;
1727 343 : if (((psL2ABandDesc != nullptr &&
1728 300 : psL2ABandDesc->eLocation == TL_IMG_DATA_Rxxm) ||
1729 682 : (psL2ABandDesc == nullptr && nPrecisionL2A != 0)) &&
1730 125 : (!procBaseLineIs1 || osBandName != "SCL"))
1731 : {
1732 123 : osTile += CPLSPrintf("R%02dm", nPrecisionL2A);
1733 123 : osTile += chSeparator;
1734 : }
1735 300 : if (procBaseLineIs1)
1736 : {
1737 185 : if (atoi(osBandName.c_str()) > 0)
1738 : {
1739 179 : osJPEG2000Name[9] = 'M';
1740 179 : osJPEG2000Name[10] = 'S';
1741 179 : osJPEG2000Name[11] = 'I';
1742 : }
1743 6 : else if (nPrecisionL2A && osBandName.size() == 3)
1744 : {
1745 6 : osJPEG2000Name[9] = osBandName[0];
1746 6 : osJPEG2000Name[10] = osBandName[1];
1747 6 : osJPEG2000Name[11] = osBandName[2];
1748 : }
1749 : }
1750 210 : else if (osProductURI.size() > 44 &&
1751 210 : osProductURI.substr(3, 8) == "_MSIL2A_")
1752 : {
1753 95 : osTile += osProductURI.substr(38, 6);
1754 95 : osTile += osProductURI.substr(10, 16);
1755 95 : granuleNameMatchTilename = false;
1756 : }
1757 : else
1758 : {
1759 20 : CPLDebug("SENTINEL2", "Invalid granule path: %s",
1760 : osGranulePath.c_str());
1761 : }
1762 300 : if (granuleNameMatchTilename)
1763 205 : osTile += osJPEG2000Name;
1764 300 : if (atoi(osBandName.c_str()) > 0)
1765 : {
1766 257 : osTile += "_B";
1767 257 : if (osBandName.size() == 3 && osBandName[0] == '0')
1768 12 : osTile += osBandName.substr(1);
1769 : else
1770 245 : osTile += osBandName;
1771 : }
1772 43 : else if (!procBaseLineIs1)
1773 : {
1774 37 : osTile += "_";
1775 37 : osTile += osBandName;
1776 : }
1777 300 : if (nPrecisionL2A)
1778 125 : osTile += CPLSPrintf("_%02dm", nPrecisionL2A);
1779 : }
1780 343 : osTile += ".jp2";
1781 686 : return osTile;
1782 : }
1783 :
1784 : /************************************************************************/
1785 : /* SENTINEL2GetMainMTDFilenameFromGranuleMTD() */
1786 : /************************************************************************/
1787 :
1788 : static CPLString
1789 29 : SENTINEL2GetMainMTDFilenameFromGranuleMTD(const char *pszFilename)
1790 : {
1791 : // Look for product MTD file
1792 58 : CPLString osTopDir(CPLFormFilenameSafe(
1793 58 : CPLFormFilenameSafe(CPLGetDirnameSafe(pszFilename).c_str(), "..",
1794 : nullptr)
1795 : .c_str(),
1796 58 : "..", nullptr));
1797 :
1798 : // Workaround to avoid long filenames on Windows
1799 29 : if (CPLIsFilenameRelative(pszFilename))
1800 : {
1801 : // GRANULE/bla/bla.xml
1802 36 : const std::string osPath = CPLGetPathSafe(pszFilename);
1803 19 : if (osPath.find('/') != std::string::npos ||
1804 1 : osPath.find('\\') != std::string::npos)
1805 : {
1806 17 : osTopDir = CPLGetPathSafe(CPLGetPathSafe(osPath.c_str()).c_str());
1807 17 : if (osTopDir == "")
1808 0 : osTopDir = ".";
1809 : }
1810 : }
1811 :
1812 29 : char **papszContents = VSIReadDir(osTopDir);
1813 29 : CPLString osMainMTD;
1814 226 : for (char **papszIter = papszContents; papszIter && *papszIter; ++papszIter)
1815 : {
1816 32 : if (strlen(*papszIter) >= strlen("S2A_XXXX_MTD") &&
1817 248 : IsS2Prefixed(*papszIter, "") &&
1818 19 : EQUALN(*papszIter + strlen("S2A_XXXX"), "_MTD", 4))
1819 : {
1820 19 : osMainMTD = CPLFormFilenameSafe(osTopDir, *papszIter, nullptr);
1821 19 : break;
1822 : }
1823 : }
1824 29 : CSLDestroy(papszContents);
1825 58 : return osMainMTD;
1826 : }
1827 :
1828 : /************************************************************************/
1829 : /* SENTINEL2GetResolutionSetAndMainMDFromGranule() */
1830 : /************************************************************************/
1831 :
1832 29 : static void SENTINEL2GetResolutionSetAndMainMDFromGranule(
1833 : const char *pszFilename, const char *pszRootPathWithoutEqual,
1834 : int nResolutionOfInterest, std::set<int> &oSetResolutions,
1835 : std::map<int, std::set<CPLString>> &oMapResolutionsToBands, char **&papszMD,
1836 : CPLXMLNode **ppsRootMainMTD)
1837 : {
1838 58 : CPLString osMainMTD(SENTINEL2GetMainMTDFilenameFromGranuleMTD(pszFilename));
1839 :
1840 : // Parse product MTD if available
1841 29 : papszMD = nullptr;
1842 48 : if (!osMainMTD.empty() &&
1843 : /* env var for debug only */
1844 19 : CPLTestBool(CPLGetConfigOption("SENTINEL2_USE_MAIN_MTD", "YES")))
1845 : {
1846 17 : CPLXMLNode *psRootMainMTD = CPLParseXMLFile(osMainMTD);
1847 17 : if (psRootMainMTD != nullptr)
1848 : {
1849 17 : CPLStripXMLNamespace(psRootMainMTD, nullptr, TRUE);
1850 :
1851 17 : CPLXMLNode *psProductInfo = CPLGetXMLNode(
1852 : psRootMainMTD, CPLSPrintf("=%s.General_Info.Product_Info",
1853 : pszRootPathWithoutEqual));
1854 17 : if (psProductInfo != nullptr)
1855 : {
1856 17 : SENTINEL2GetResolutionSet(psProductInfo, oSetResolutions,
1857 : oMapResolutionsToBands);
1858 : }
1859 :
1860 17 : papszMD = SENTINEL2GetUserProductMetadata(psRootMainMTD,
1861 : pszRootPathWithoutEqual);
1862 17 : if (ppsRootMainMTD != nullptr)
1863 6 : *ppsRootMainMTD = psRootMainMTD;
1864 : else
1865 11 : CPLDestroyXMLNode(psRootMainMTD);
1866 : }
1867 : }
1868 : else
1869 : {
1870 : // If no main MTD file found, then probe all bands for resolution (of
1871 : // interest if there's one, or all resolutions otherwise)
1872 168 : for (const auto &sBandDesc : asBandDesc)
1873 : {
1874 156 : if (nResolutionOfInterest != 0 &&
1875 117 : sBandDesc.nResolution != nResolutionOfInterest)
1876 : {
1877 83 : continue;
1878 : }
1879 : CPLString osBandName =
1880 146 : sBandDesc.pszBandName + 1; /* skip B character */
1881 73 : if (atoi(osBandName) < 10)
1882 62 : osBandName = "0" + osBandName;
1883 :
1884 : CPLString osTile(SENTINEL2GetTilename(
1885 146 : CPLGetPathSafe(pszFilename).c_str(),
1886 365 : CPLGetBasenameSafe(pszFilename).c_str(), osBandName));
1887 : VSIStatBufL sStat;
1888 73 : if (VSIStatExL(osTile, &sStat, VSI_STAT_EXISTS_FLAG) == 0)
1889 : {
1890 20 : oMapResolutionsToBands[sBandDesc.nResolution].insert(
1891 20 : std::move(osBandName));
1892 20 : oSetResolutions.insert(sBandDesc.nResolution);
1893 : }
1894 : }
1895 : }
1896 29 : }
1897 :
1898 : /************************************************************************/
1899 : /* OpenL1BGranule() */
1900 : /************************************************************************/
1901 :
1902 18 : GDALDataset *SENTINEL2Dataset::OpenL1BGranule(const char *pszFilename,
1903 : CPLXMLNode **ppsRoot,
1904 : int nResolutionOfInterest,
1905 : std::set<CPLString> *poBandSet)
1906 : {
1907 18 : CPLXMLNode *psRoot = CPLParseXMLFile(pszFilename);
1908 18 : if (psRoot == nullptr)
1909 : {
1910 3 : CPLDebug("SENTINEL2", "Cannot XML parse %s", pszFilename);
1911 3 : return nullptr;
1912 : }
1913 :
1914 15 : char *pszOriginalXML = CPLSerializeXMLTree(psRoot);
1915 30 : CPLString osOriginalXML;
1916 15 : if (pszOriginalXML)
1917 15 : osOriginalXML = pszOriginalXML;
1918 15 : CPLFree(pszOriginalXML);
1919 :
1920 30 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRoot);
1921 15 : CPLStripXMLNamespace(psRoot, nullptr, TRUE);
1922 :
1923 15 : SENTINEL2DatasetContainer *poDS = new SENTINEL2DatasetContainer();
1924 :
1925 15 : if (!osOriginalXML.empty())
1926 : {
1927 : char *apszXMLMD[2];
1928 15 : apszXMLMD[0] = const_cast<char *>(osOriginalXML.c_str());
1929 15 : apszXMLMD[1] = nullptr;
1930 15 : poDS->GDALDataset::SetMetadata(apszXMLMD, "xml:SENTINEL2");
1931 : }
1932 :
1933 30 : std::set<int> oSetResolutions;
1934 15 : std::map<int, std::set<CPLString>> oMapResolutionsToBands;
1935 15 : char **papszMD = nullptr;
1936 15 : SENTINEL2GetResolutionSetAndMainMDFromGranule(
1937 : pszFilename, "Level-1B_User_Product", nResolutionOfInterest,
1938 : oSetResolutions, oMapResolutionsToBands, papszMD, nullptr);
1939 15 : if (poBandSet != nullptr)
1940 12 : *poBandSet = oMapResolutionsToBands[nResolutionOfInterest];
1941 :
1942 15 : char **papszGranuleMD = SENTINEL2GetL1BGranuleMetadata(psRoot);
1943 15 : papszMD = CSLMerge(papszMD, papszGranuleMD);
1944 15 : CSLDestroy(papszGranuleMD);
1945 :
1946 : // Remove CLOUD_COVERAGE_ASSESSMENT that comes from main metadata, if
1947 : // granule CLOUDY_PIXEL_PERCENTAGE is present.
1948 21 : if (CSLFetchNameValue(papszMD, "CLOUDY_PIXEL_PERCENTAGE") != nullptr &&
1949 6 : CSLFetchNameValue(papszMD, "CLOUD_COVERAGE_ASSESSMENT") != nullptr)
1950 : {
1951 6 : papszMD =
1952 6 : CSLSetNameValue(papszMD, "CLOUD_COVERAGE_ASSESSMENT", nullptr);
1953 : }
1954 :
1955 15 : poDS->GDALDataset::SetMetadata(papszMD);
1956 15 : CSLDestroy(papszMD);
1957 :
1958 : // Get the footprint
1959 : const char *pszPosList =
1960 15 : CPLGetXMLValue(psRoot,
1961 : "=Level-1B_Granule_ID.Geometric_Info.Granule_Footprint."
1962 : "Granule_Footprint.Footprint.EXT_POS_LIST",
1963 : nullptr);
1964 15 : if (pszPosList != nullptr)
1965 : {
1966 12 : CPLString osPolygon = SENTINEL2GetPolygonWKTFromPosList(pszPosList);
1967 6 : if (!osPolygon.empty())
1968 6 : poDS->GDALDataset::SetMetadataItem("FOOTPRINT", osPolygon.c_str());
1969 : }
1970 :
1971 : /* Create subdatsets per resolution (10, 20, 60m) */
1972 15 : int iSubDSNum = 1;
1973 37 : for (std::set<int>::const_iterator oIterRes = oSetResolutions.begin();
1974 59 : oIterRes != oSetResolutions.end(); ++oIterRes)
1975 : {
1976 22 : const int nResolution = *oIterRes;
1977 :
1978 22 : poDS->GDALDataset::SetMetadataItem(
1979 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
1980 : CPLSPrintf("SENTINEL2_L1B:%s:%dm", pszFilename, nResolution),
1981 : "SUBDATASETS");
1982 :
1983 : CPLString osBandNames = SENTINEL2GetBandListForResolution(
1984 44 : oMapResolutionsToBands[nResolution]);
1985 :
1986 : CPLString osDesc(CPLSPrintf("Bands %s with %dm resolution",
1987 22 : osBandNames.c_str(), nResolution));
1988 22 : poDS->GDALDataset::SetMetadataItem(
1989 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
1990 : "SUBDATASETS");
1991 :
1992 22 : iSubDSNum++;
1993 : }
1994 :
1995 15 : if (ppsRoot != nullptr)
1996 : {
1997 12 : *ppsRoot = oXMLHolder.Release();
1998 : }
1999 :
2000 15 : return poDS;
2001 : }
2002 :
2003 : /************************************************************************/
2004 : /* SENTINEL2SetBandMetadata() */
2005 : /************************************************************************/
2006 :
2007 212 : static void SENTINEL2SetBandMetadata(GDALRasterBand *poBand,
2008 : const std::string &osBandName)
2009 : {
2010 424 : CPLString osLookupBandName(osBandName);
2011 212 : if (osLookupBandName[0] == '0')
2012 141 : osLookupBandName = osLookupBandName.substr(1);
2013 212 : if (atoi(osLookupBandName) > 0)
2014 171 : osLookupBandName = "B" + osLookupBandName;
2015 :
2016 424 : CPLString osBandDesc(osLookupBandName);
2017 : const SENTINEL2BandDescription *psBandDesc =
2018 212 : SENTINEL2GetBandDesc(osLookupBandName);
2019 212 : if (psBandDesc != nullptr)
2020 : {
2021 : osBandDesc +=
2022 171 : CPLSPrintf(", central wavelength %d nm", psBandDesc->nWaveLength);
2023 171 : poBand->SetColorInterpretation(psBandDesc->eColorInterp);
2024 171 : poBand->SetMetadataItem("BANDNAME", psBandDesc->pszBandName);
2025 171 : poBand->SetMetadataItem("BANDWIDTH",
2026 171 : CPLSPrintf("%d", psBandDesc->nBandWidth));
2027 171 : poBand->SetMetadataItem("BANDWIDTH_UNIT", "nm");
2028 171 : poBand->SetMetadataItem("WAVELENGTH",
2029 171 : CPLSPrintf("%d", psBandDesc->nWaveLength));
2030 171 : poBand->SetMetadataItem("WAVELENGTH_UNIT", "nm");
2031 :
2032 171 : poBand->SetMetadataItem(
2033 : "CENTRAL_WAVELENGTH_UM",
2034 171 : CPLSPrintf("%.3f", double(psBandDesc->nWaveLength) / 1000),
2035 171 : "IMAGERY");
2036 171 : poBand->SetMetadataItem(
2037 : "FWHM_UM",
2038 171 : CPLSPrintf("%.3f", double(psBandDesc->nBandWidth) / 1000),
2039 171 : "IMAGERY");
2040 : }
2041 : else
2042 : {
2043 : const SENTINEL2_L2A_BandDescription *psL2ABandDesc =
2044 41 : SENTINEL2GetL2ABandDesc(osBandName.c_str());
2045 41 : if (psL2ABandDesc != nullptr)
2046 : {
2047 41 : osBandDesc += ", ";
2048 41 : osBandDesc += psL2ABandDesc->pszBandDescription;
2049 : }
2050 :
2051 41 : poBand->SetMetadataItem("BANDNAME", osBandName.c_str());
2052 : }
2053 212 : poBand->SetDescription(osBandDesc);
2054 212 : }
2055 :
2056 : /************************************************************************/
2057 : /* OpenL1BSubdataset() */
2058 : /************************************************************************/
2059 :
2060 18 : GDALDataset *SENTINEL2Dataset::OpenL1BSubdataset(GDALOpenInfo *poOpenInfo)
2061 : {
2062 36 : CPLString osFilename;
2063 18 : CPLAssert(STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1B:"));
2064 18 : osFilename = poOpenInfo->pszFilename + strlen("SENTINEL2_L1B:");
2065 18 : const char *pszPrecision = strrchr(osFilename.c_str(), ':');
2066 18 : if (pszPrecision == nullptr || pszPrecision == osFilename.c_str())
2067 : {
2068 2 : CPLError(CE_Failure, CPLE_AppDefined,
2069 : "Invalid syntax for SENTINEL2_L1B:");
2070 2 : return nullptr;
2071 : }
2072 16 : const int nSubDSPrecision = atoi(pszPrecision + 1);
2073 16 : if (nSubDSPrecision != RES_10M && nSubDSPrecision != RES_20M &&
2074 : nSubDSPrecision != RES_60M)
2075 : {
2076 2 : CPLError(CE_Failure, CPLE_NotSupported, "Unsupported precision: %d",
2077 : nSubDSPrecision);
2078 2 : return nullptr;
2079 : }
2080 14 : osFilename.resize(pszPrecision - osFilename.c_str());
2081 :
2082 14 : CPLXMLNode *psRoot = nullptr;
2083 28 : std::set<CPLString> oSetBands;
2084 : GDALDataset *poTmpDS =
2085 14 : OpenL1BGranule(osFilename, &psRoot, nSubDSPrecision, &oSetBands);
2086 14 : if (poTmpDS == nullptr)
2087 : {
2088 2 : CPLDebug("SENTINEL2", "Failed to open L1B granule %s",
2089 : osFilename.c_str());
2090 2 : return nullptr;
2091 : }
2092 :
2093 24 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRoot);
2094 :
2095 24 : std::vector<CPLString> aosBands;
2096 31 : for (std::set<CPLString>::const_iterator oIterBandnames = oSetBands.begin();
2097 50 : oIterBandnames != oSetBands.end(); ++oIterBandnames)
2098 : {
2099 19 : aosBands.push_back(*oIterBandnames);
2100 : }
2101 : /* Put 2=Blue, 3=Green, 4=Band bands in RGB order for conveniency */
2102 13 : if (aosBands.size() >= 3 && aosBands[0] == "02" && aosBands[1] == "03" &&
2103 1 : aosBands[2] == "04")
2104 : {
2105 1 : aosBands[0] = "04";
2106 1 : aosBands[2] = "02";
2107 : }
2108 :
2109 12 : int nBits = 0; /* 0 = unknown yet*/
2110 12 : int nValMax = 0; /* 0 = unknown yet*/
2111 12 : int nRows = 0;
2112 12 : int nCols = 0;
2113 12 : CPLXMLNode *psGranuleDimensions = CPLGetXMLNode(
2114 : psRoot, "=Level-1B_Granule_ID.Geometric_Info.Granule_Dimensions");
2115 12 : if (psGranuleDimensions == nullptr)
2116 : {
2117 3 : for (size_t i = 0; i < aosBands.size(); i++)
2118 : {
2119 : CPLString osTile(SENTINEL2GetTilename(
2120 2 : CPLGetPathSafe(osFilename).c_str(),
2121 4 : CPLGetBasenameSafe(osFilename).c_str(), aosBands[i]));
2122 1 : if (SENTINEL2GetTileInfo(osTile, &nCols, &nRows, &nBits))
2123 : {
2124 1 : if (nBits <= 16)
2125 1 : nValMax = (1 << nBits) - 1;
2126 : else
2127 : {
2128 0 : CPLDebug("SENTINEL2", "Unexpected bit depth %d", nBits);
2129 0 : nValMax = 65535;
2130 : }
2131 1 : break;
2132 : }
2133 : }
2134 : }
2135 : else
2136 : {
2137 9 : for (CPLXMLNode *psIter = psGranuleDimensions->psChild;
2138 23 : psIter != nullptr; psIter = psIter->psNext)
2139 : {
2140 22 : if (psIter->eType != CXT_Element)
2141 9 : continue;
2142 26 : if (EQUAL(psIter->pszValue, "Size") &&
2143 13 : atoi(CPLGetXMLValue(psIter, "resolution", "")) ==
2144 : nSubDSPrecision)
2145 : {
2146 8 : const char *pszRows = CPLGetXMLValue(psIter, "NROWS", nullptr);
2147 8 : if (pszRows == nullptr)
2148 : {
2149 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
2150 : "NROWS");
2151 1 : delete poTmpDS;
2152 1 : return nullptr;
2153 : }
2154 7 : const char *pszCols = CPLGetXMLValue(psIter, "NCOLS", nullptr);
2155 7 : if (pszCols == nullptr)
2156 : {
2157 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
2158 : "NCOLS");
2159 1 : delete poTmpDS;
2160 1 : return nullptr;
2161 : }
2162 6 : nRows = atoi(pszRows);
2163 6 : nCols = atoi(pszCols);
2164 6 : break;
2165 : }
2166 : }
2167 : }
2168 10 : if (nRows <= 0 || nCols <= 0)
2169 : {
2170 3 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find granule dimension");
2171 3 : delete poTmpDS;
2172 3 : return nullptr;
2173 : }
2174 :
2175 7 : SENTINEL2Dataset *poDS = new SENTINEL2Dataset(nCols, nRows);
2176 7 : poDS->aosNonJP2Files.push_back(osFilename);
2177 :
2178 : // Transfer metadata
2179 7 : poDS->GDALDataset::SetMetadata(poTmpDS->GetMetadata());
2180 7 : poDS->GDALDataset::SetMetadata(poTmpDS->GetMetadata("xml:SENTINEL2"),
2181 : "xml:SENTINEL2");
2182 :
2183 7 : delete poTmpDS;
2184 :
2185 : /* -------------------------------------------------------------------- */
2186 : /* Initialize bands. */
2187 : /* -------------------------------------------------------------------- */
2188 :
2189 21 : int nSaturatedVal = atoi(CSLFetchNameValueDef(
2190 7 : poDS->GetMetadata(), "SPECIAL_VALUE_SATURATED", "-1"));
2191 7 : int nNodataVal = atoi(CSLFetchNameValueDef(poDS->GetMetadata(),
2192 : "SPECIAL_VALUE_NODATA", "-1"));
2193 :
2194 : const bool bAlpha =
2195 7 : CPLTestBool(SENTINEL2GetOption(poOpenInfo, "ALPHA", "FALSE"));
2196 7 : const int nBands = ((bAlpha) ? 1 : 0) + static_cast<int>(aosBands.size());
2197 7 : const int nAlphaBand = (!bAlpha) ? 0 : nBands;
2198 7 : const GDALDataType eDT = GDT_UInt16;
2199 :
2200 27 : for (int nBand = 1; nBand <= nBands; nBand++)
2201 : {
2202 20 : VRTSourcedRasterBand *poBand = nullptr;
2203 :
2204 20 : if (nBand != nAlphaBand)
2205 : {
2206 19 : poBand = new VRTSourcedRasterBand(
2207 19 : poDS, nBand, eDT, poDS->nRasterXSize, poDS->nRasterYSize);
2208 : }
2209 : else
2210 : {
2211 1 : poBand = new SENTINEL2AlphaBand(
2212 : poDS, nBand, eDT, poDS->nRasterXSize, poDS->nRasterYSize,
2213 1 : nSaturatedVal, nNodataVal);
2214 : }
2215 :
2216 20 : poDS->SetBand(nBand, poBand);
2217 20 : if (nBand == nAlphaBand)
2218 1 : poBand->SetColorInterpretation(GCI_AlphaBand);
2219 :
2220 20 : CPLString osBandName;
2221 20 : if (nBand != nAlphaBand)
2222 : {
2223 19 : osBandName = aosBands[nBand - 1];
2224 19 : SENTINEL2SetBandMetadata(poBand, osBandName);
2225 : }
2226 : else
2227 1 : osBandName = aosBands[0];
2228 :
2229 : CPLString osTile(SENTINEL2GetTilename(
2230 40 : CPLGetPathSafe(osFilename).c_str(),
2231 80 : CPLGetBasenameSafe(osFilename).c_str(), osBandName));
2232 :
2233 20 : bool bTileFound = false;
2234 20 : if (nValMax == 0)
2235 : {
2236 : /* It is supposed to be 12 bits, but some products have 15 bits */
2237 6 : if (SENTINEL2GetTileInfo(osTile, nullptr, nullptr, &nBits))
2238 : {
2239 5 : bTileFound = true;
2240 5 : if (nBits <= 16)
2241 5 : nValMax = (1 << nBits) - 1;
2242 : else
2243 : {
2244 0 : CPLDebug("SENTINEL2", "Unexpected bit depth %d", nBits);
2245 0 : nValMax = 65535;
2246 : }
2247 : }
2248 : }
2249 : else
2250 : {
2251 : VSIStatBufL sStat;
2252 14 : if (VSIStatExL(osTile, &sStat, VSI_STAT_EXISTS_FLAG) == 0)
2253 14 : bTileFound = true;
2254 : }
2255 20 : if (!bTileFound)
2256 : {
2257 1 : CPLError(CE_Warning, CPLE_AppDefined,
2258 : "Tile %s not found on filesystem. Skipping it",
2259 : osTile.c_str());
2260 1 : continue;
2261 : }
2262 :
2263 19 : if (nBand != nAlphaBand)
2264 : {
2265 18 : poBand->AddSimpleSource(osTile, 1, 0, 0, poDS->nRasterXSize,
2266 18 : poDS->nRasterYSize, 0, 0,
2267 18 : poDS->nRasterXSize, poDS->nRasterYSize);
2268 : }
2269 : else
2270 : {
2271 1 : poBand->AddComplexSource(osTile, 1, 0, 0, poDS->nRasterXSize,
2272 1 : poDS->nRasterYSize, 0, 0,
2273 1 : poDS->nRasterXSize, poDS->nRasterYSize,
2274 : nValMax /* offset */, 0 /* scale */);
2275 : }
2276 :
2277 19 : if ((nBits % 8) != 0)
2278 : {
2279 19 : poBand->SetMetadataItem("NBITS", CPLSPrintf("%d", nBits),
2280 19 : "IMAGE_STRUCTURE");
2281 : }
2282 : }
2283 :
2284 : /* -------------------------------------------------------------------- */
2285 : /* Add georeferencing. */
2286 : /* -------------------------------------------------------------------- */
2287 : // const char* pszDirection =
2288 : // poDS->GetMetadataItem("DATATAKE_1_SENSING_ORBIT_DIRECTION");
2289 7 : const char *pszFootprint = poDS->GetMetadataItem("FOOTPRINT");
2290 7 : if (pszFootprint != nullptr)
2291 : {
2292 : // For descending orbits, we have observed that the order of points in
2293 : // the polygon is UL, LL, LR, UR. That might not be true for ascending
2294 : // orbits but let's assume it...
2295 4 : OGRGeometry *poGeom = nullptr;
2296 4 : if (OGRGeometryFactory::createFromWkt(pszFootprint, nullptr, &poGeom) ==
2297 4 : OGRERR_NONE &&
2298 8 : poGeom != nullptr &&
2299 4 : wkbFlatten(poGeom->getGeometryType()) == wkbPolygon)
2300 : {
2301 : OGRLinearRing *poRing =
2302 4 : reinterpret_cast<OGRPolygon *>(poGeom)->getExteriorRing();
2303 4 : if (poRing != nullptr && poRing->getNumPoints() == 5)
2304 : {
2305 : GDAL_GCP asGCPList[5];
2306 4 : memset(asGCPList, 0, sizeof(asGCPList));
2307 20 : for (int i = 0; i < 4; i++)
2308 : {
2309 16 : asGCPList[i].dfGCPX = poRing->getX(i);
2310 16 : asGCPList[i].dfGCPY = poRing->getY(i);
2311 16 : asGCPList[i].dfGCPZ = poRing->getZ(i);
2312 : }
2313 4 : asGCPList[0].dfGCPPixel = 0;
2314 4 : asGCPList[0].dfGCPLine = 0;
2315 4 : asGCPList[1].dfGCPPixel = 0;
2316 4 : asGCPList[1].dfGCPLine = poDS->nRasterYSize;
2317 4 : asGCPList[2].dfGCPPixel = poDS->nRasterXSize;
2318 4 : asGCPList[2].dfGCPLine = poDS->nRasterYSize;
2319 4 : asGCPList[3].dfGCPPixel = poDS->nRasterXSize;
2320 4 : asGCPList[3].dfGCPLine = 0;
2321 :
2322 4 : int nGCPCount = 4;
2323 : CPLXMLNode *psGeometryHeader =
2324 4 : CPLGetXMLNode(psRoot, "=Level-1B_Granule_ID.Geometric_Info."
2325 : "Granule_Position.Geometric_Header");
2326 4 : if (psGeometryHeader != nullptr)
2327 : {
2328 4 : const char *pszGC = CPLGetXMLValue(
2329 : psGeometryHeader, "GROUND_CENTER", nullptr);
2330 : const char *pszQLCenter =
2331 4 : CPLGetXMLValue(psGeometryHeader, "QL_CENTER", nullptr);
2332 4 : if (pszGC != nullptr && pszQLCenter != nullptr &&
2333 4 : EQUAL(pszQLCenter, "0 0"))
2334 : {
2335 4 : char **papszTokens = CSLTokenizeString(pszGC);
2336 4 : if (CSLCount(papszTokens) >= 2)
2337 : {
2338 4 : nGCPCount = 5;
2339 4 : asGCPList[4].dfGCPX = CPLAtof(papszTokens[1]);
2340 4 : asGCPList[4].dfGCPY = CPLAtof(papszTokens[0]);
2341 4 : if (CSLCount(papszTokens) >= 3)
2342 4 : asGCPList[4].dfGCPZ = CPLAtof(papszTokens[2]);
2343 4 : asGCPList[4].dfGCPLine = poDS->nRasterYSize / 2.0;
2344 4 : asGCPList[4].dfGCPPixel = poDS->nRasterXSize / 2.0;
2345 : }
2346 4 : CSLDestroy(papszTokens);
2347 : }
2348 : }
2349 :
2350 4 : poDS->SetGCPs(nGCPCount, asGCPList, SRS_WKT_WGS84_LAT_LONG);
2351 4 : GDALDeinitGCPs(nGCPCount, asGCPList);
2352 : }
2353 : }
2354 4 : delete poGeom;
2355 : }
2356 :
2357 : /* -------------------------------------------------------------------- */
2358 : /* Initialize overview information. */
2359 : /* -------------------------------------------------------------------- */
2360 7 : poDS->SetDescription(poOpenInfo->pszFilename);
2361 7 : CPLString osOverviewFile;
2362 : osOverviewFile =
2363 7 : CPLSPrintf("%s_%dm.tif.ovr", osFilename.c_str(), nSubDSPrecision);
2364 7 : poDS->SetMetadataItem("OVERVIEW_FILE", osOverviewFile, "OVERVIEWS");
2365 7 : poDS->oOvManager.Initialize(poDS, ":::VIRTUAL:::");
2366 :
2367 7 : return poDS;
2368 : }
2369 :
2370 : /************************************************************************/
2371 : /* OpenL1BSubdatasetWithGeoloc() */
2372 : /************************************************************************/
2373 :
2374 : GDALDataset *
2375 9 : SENTINEL2Dataset::OpenL1BSubdatasetWithGeoloc(GDALOpenInfo *poOpenInfo)
2376 : {
2377 18 : CPLString osFilename;
2378 9 : CPLAssert(
2379 : STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1B_WITH_GEOLOC:"));
2380 : const CPLStringList aosTokens(
2381 9 : CSLTokenizeString2(poOpenInfo->pszFilename, ":",
2382 18 : CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES));
2383 9 : if (aosTokens.size() != 3)
2384 : {
2385 1 : CPLError(CE_Failure, CPLE_AppDefined,
2386 : "OpenL1BSubdatasetWithGeoloc(): invalid number of tokens in "
2387 : "subdataset name");
2388 1 : return nullptr;
2389 : }
2390 8 : const char *pszMainXMLFilename = aosTokens[1];
2391 8 : const char *pszGeolocVRT = aosTokens[2];
2392 :
2393 8 : const size_t nLen = strlen(pszGeolocVRT);
2394 8 : if (nLen <= strlen("_Dxx_Byy") || pszGeolocVRT[nLen - 7] != 'D' ||
2395 7 : pszGeolocVRT[nLen - 3] != 'B')
2396 : {
2397 1 : CPLError(CE_Failure, CPLE_AppDefined,
2398 : "Invalid subdataset component name");
2399 1 : return nullptr;
2400 : }
2401 :
2402 : // Find in which datastrip we are
2403 : const std::string osDatastrip(
2404 7 : CPLString(pszGeolocVRT, nLen - strlen("_Dxx_Byy"))
2405 21 : .replaceAll("_GEO_", "_MSI_"));
2406 14 : const std::string osDetectorId(pszGeolocVRT + nLen - 6, 2);
2407 14 : const std::string osBandId(pszGeolocVRT + nLen - 2, 2);
2408 :
2409 7 : const char chSeparator = SENTINEL2GetPathSeparator(pszMainXMLFilename);
2410 7 : const char achSeparator[] = {chSeparator, 0};
2411 :
2412 : const std::string osDatastripRoot(
2413 7 : std::string(CPLGetDirnameSafe(pszMainXMLFilename))
2414 7 : .append(achSeparator)
2415 14 : .append("DATASTRIP"));
2416 14 : const CPLStringList aosList(VSIReadDir(osDatastripRoot.c_str()));
2417 14 : std::string osDatastripSubdir;
2418 10 : for (const char *pszDatastripId : aosList)
2419 : {
2420 8 : if (STARTS_WITH_CI(pszDatastripId, osDatastrip.c_str()))
2421 : {
2422 5 : osDatastripSubdir = pszDatastripId;
2423 5 : break;
2424 : }
2425 : }
2426 7 : if (osDatastripSubdir.empty())
2427 : {
2428 2 : CPLError(CE_Failure, CPLE_AppDefined,
2429 : "Cannot find a file in %s starting with %s",
2430 : osDatastripRoot.c_str(), osDatastrip.c_str());
2431 2 : return nullptr;
2432 : }
2433 :
2434 5 : const std::string osDatastripSubdirFull = std::string(osDatastripRoot)
2435 5 : .append(achSeparator)
2436 10 : .append(osDatastripSubdir);
2437 5 : CPL_IGNORE_RET_VAL(osDatastripRoot);
2438 : const std::string osXMLDatastrip =
2439 5 : std::string(osDatastripSubdirFull)
2440 5 : .append(achSeparator)
2441 10 : .append(CPLString(osDatastrip).replaceAll("_MSI_", "_MTD_"))
2442 10 : .append(".xml");
2443 5 : if (CPLHasPathTraversal(osXMLDatastrip.c_str()))
2444 : {
2445 0 : CPLError(CE_Failure, CPLE_NotSupported, "Path traversal detected in %s",
2446 : osXMLDatastrip.c_str());
2447 0 : return nullptr;
2448 : }
2449 10 : CPLXMLTreeCloser poXML(CPLParseXMLFile(osXMLDatastrip.c_str()));
2450 5 : if (!poXML)
2451 : {
2452 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot XML parse %s",
2453 : osXMLDatastrip.c_str());
2454 0 : return nullptr;
2455 : }
2456 5 : CPLStripXMLNamespace(poXML.get(), nullptr, TRUE);
2457 :
2458 : const CPLXMLNode *psDetectorList =
2459 5 : CPLGetXMLNode(poXML.get(), "=Level-1B_DataStrip_ID.Image_Data_Info."
2460 : "Granules_Information.Detector_List");
2461 5 : if (!psDetectorList)
2462 : {
2463 0 : CPLError(CE_Failure, CPLE_AppDefined,
2464 : "Cannot find <Detector_List> in %s", osXMLDatastrip.c_str());
2465 0 : return nullptr;
2466 : }
2467 :
2468 : struct GranuleDesc
2469 : {
2470 : const char *pszGranuleId = nullptr;
2471 : int nPosition = 0;
2472 : };
2473 :
2474 10 : std::vector<GranuleDesc> asGranules;
2475 :
2476 : // Get the list of granules for the current detector and band
2477 6 : for (const CPLXMLNode *psDetector = psDetectorList->psChild; psDetector;
2478 1 : psDetector = psDetector->psNext)
2479 : {
2480 15 : if (psDetector->eType == CXT_Element &&
2481 10 : EQUAL(psDetector->pszValue, "Detector") &&
2482 5 : CPLGetXMLValue(psDetector, "detectorId", "") == osDetectorId)
2483 : {
2484 : const CPLXMLNode *psGranuleList =
2485 4 : CPLGetXMLNode(psDetector, "Granule_List");
2486 4 : if (psGranuleList)
2487 : {
2488 4 : for (const CPLXMLNode *psGranule = psGranuleList->psChild;
2489 24 : psGranule; psGranule = psGranule->psNext)
2490 : {
2491 20 : if (psGranule->eType == CXT_Element &&
2492 20 : EQUAL(psGranule->pszValue, "Granule"))
2493 : {
2494 : const char *pszGranuleId =
2495 20 : CPLGetXMLValue(psGranule, "granuleId", "");
2496 : const char *pszPosition =
2497 20 : CPLGetXMLValue(psGranule, "POSITION", "");
2498 20 : GranuleDesc sDesc;
2499 20 : sDesc.pszGranuleId = pszGranuleId;
2500 20 : sDesc.nPosition = atoi(pszPosition);
2501 20 : asGranules.push_back(sDesc);
2502 : }
2503 : }
2504 : }
2505 4 : break;
2506 : }
2507 : }
2508 5 : if (asGranules.empty())
2509 : {
2510 1 : CPLError(CE_Failure, CPLE_AppDefined,
2511 : "Cannot find granules for detector %s in %s",
2512 : osDetectorId.c_str(), osXMLDatastrip.c_str());
2513 1 : return nullptr;
2514 : }
2515 4 : std::sort(asGranules.begin(), asGranules.end(),
2516 32 : [](const GranuleDesc &sDescA, const GranuleDesc &sDescB)
2517 32 : { return sDescA.nPosition < sDescB.nPosition; });
2518 4 : const int nGranuleCount = static_cast<int>(asGranules.size());
2519 4 : constexpr int ROWS_PER_10M_GRANULE = 2304;
2520 24 : for (int i = 0; i < nGranuleCount; ++i)
2521 : {
2522 20 : if (asGranules[i].nPosition != 1 + i * ROWS_PER_10M_GRANULE)
2523 : {
2524 0 : CPLError(
2525 : CE_Failure, CPLE_NotSupported,
2526 : "Granule %s is declared to be at position %d, was expecting %d",
2527 0 : asGranules[i].pszGranuleId, asGranules[i].nPosition,
2528 0 : 1 + i * ROWS_PER_10M_GRANULE);
2529 0 : return nullptr;
2530 : }
2531 : }
2532 :
2533 8 : const std::string osBandName = std::string("B").append(
2534 4 : osBandId == "8A"
2535 12 : ? osBandId
2536 8 : : std::string(CPLSPrintf("%d", atoi(osBandId.c_str()))));
2537 : const SENTINEL2BandDescription *psBandDesc =
2538 4 : SENTINEL2GetBandDesc(osBandName.c_str());
2539 4 : if (psBandDesc == nullptr)
2540 : {
2541 1 : CPLError(CE_Failure, CPLE_NotSupported, "Unknown band id %s",
2542 : osBandId.c_str());
2543 1 : return nullptr;
2544 : }
2545 :
2546 3 : const int nResolution = psBandDesc->nResolution;
2547 5 : const int nRowsPerGranule = nResolution == RES_10M ? ROWS_PER_10M_GRANULE
2548 2 : : nResolution == RES_20M ? 1152
2549 : : 384;
2550 5 : const int nColsPerGranule = nResolution == RES_10M ? 2552
2551 2 : : nResolution == RES_20M ? 1276
2552 : : 425;
2553 : auto poDS = std::make_unique<SENTINEL2Dataset>(
2554 6 : nColsPerGranule, nRowsPerGranule * nGranuleCount);
2555 :
2556 3 : constexpr GDALDataType eDT = GDT_UInt16;
2557 : auto poBand = new VRTSourcedRasterBand(
2558 3 : poDS.get(), 1, eDT, poDS->nRasterXSize, poDS->nRasterYSize);
2559 3 : poDS->SetBand(1, poBand);
2560 :
2561 3 : SENTINEL2SetBandMetadata(poBand, osBandId);
2562 :
2563 : // Create the virtual raster by adding each granule
2564 : const std::string osGranuleRoot =
2565 3 : std::string(CPLGetDirnameSafe(pszMainXMLFilename))
2566 3 : .append(achSeparator)
2567 6 : .append("GRANULE");
2568 3 : int nValMax = 0;
2569 3 : int nBits = 0;
2570 18 : for (int i = 0; i < nGranuleCount; ++i)
2571 : {
2572 15 : const auto &sGranule = asGranules[i];
2573 : const std::string osTile(
2574 45 : SENTINEL2GetTilename(std::string(osGranuleRoot)
2575 15 : .append(achSeparator)
2576 15 : .append(sGranule.pszGranuleId),
2577 30 : sGranule.pszGranuleId, osBandId.c_str()));
2578 :
2579 15 : bool bTileFound = false;
2580 15 : if (nValMax == 0)
2581 : {
2582 : /* It is supposed to be 12 bits, but some products have 15 bits */
2583 3 : int nGranuleWidth = 0;
2584 3 : int nGranuleHeight = 0;
2585 3 : if (SENTINEL2GetTileInfo(osTile.c_str(), &nGranuleWidth,
2586 : &nGranuleHeight, &nBits))
2587 : {
2588 3 : bTileFound = true;
2589 3 : if (nBits <= 16)
2590 3 : nValMax = (1 << nBits) - 1;
2591 : else
2592 : {
2593 0 : CPLDebug("SENTINEL2", "Unexpected bit depth %d", nBits);
2594 0 : nValMax = 65535;
2595 : }
2596 3 : if (nGranuleWidth != nColsPerGranule ||
2597 3 : nGranuleHeight != nRowsPerGranule)
2598 : {
2599 0 : CPLError(CE_Failure, CPLE_AppDefined,
2600 : "Tile %s has not expected dimensions. "
2601 : "Got %dx%d, expected %dx%d",
2602 : osTile.c_str(), nGranuleWidth, nGranuleHeight,
2603 : nColsPerGranule, nRowsPerGranule);
2604 0 : return nullptr;
2605 : }
2606 : }
2607 : }
2608 : else
2609 : {
2610 : VSIStatBufL sStat;
2611 12 : if (VSIStatExL(osTile.c_str(), &sStat, VSI_STAT_EXISTS_FLAG) == 0)
2612 12 : bTileFound = true;
2613 : }
2614 15 : if (!bTileFound)
2615 : {
2616 0 : CPLError(CE_Failure, CPLE_AppDefined, "Tile %s not found.",
2617 : osTile.c_str());
2618 0 : return nullptr;
2619 : }
2620 :
2621 15 : poBand->AddSimpleSource(osTile.c_str(), 1, 0, 0, nColsPerGranule,
2622 15 : nRowsPerGranule, 0, i * nRowsPerGranule,
2623 : nColsPerGranule, nRowsPerGranule);
2624 : }
2625 :
2626 3 : if ((nBits % 8) != 0)
2627 : {
2628 0 : poBand->SetMetadataItem("NBITS", CPLSPrintf("%d", nBits),
2629 0 : "IMAGE_STRUCTURE");
2630 : }
2631 :
2632 : // Get metadata from top MTD XML filename
2633 6 : std::set<int> oSetResolutions;
2634 6 : std::map<int, std::set<CPLString>> oMapResolutionsToBands;
2635 3 : char **papszMD = nullptr;
2636 : std::string osGranuleMTD = CPLFormFilenameSafe(
2637 3 : CPLFormFilenameSafe(osGranuleRoot.c_str(), asGranules[0].pszGranuleId,
2638 : nullptr)
2639 : .c_str(),
2640 9 : asGranules[0].pszGranuleId, nullptr);
2641 3 : if (osGranuleMTD.size() > strlen("_NXX.YY") &&
2642 6 : osGranuleMTD[osGranuleMTD.size() - 7] == '_' &&
2643 3 : osGranuleMTD[osGranuleMTD.size() - 6] == 'N')
2644 : {
2645 3 : osGranuleMTD.resize(osGranuleMTD.size() - strlen("_NXX.YY"));
2646 : }
2647 3 : SENTINEL2GetResolutionSetAndMainMDFromGranule(
2648 : osGranuleMTD.c_str(), "Level-1B_User_Product", nResolution,
2649 : oSetResolutions, oMapResolutionsToBands, papszMD, nullptr);
2650 144 : for (const auto [pszKey, pszValue] :
2651 147 : cpl::IterateNameValue(const_cast<CSLConstList>(papszMD)))
2652 : {
2653 72 : if (!EQUAL(pszKey, "FOOTPRINT"))
2654 72 : poDS->SetMetadataItem(pszKey, pszValue);
2655 : }
2656 3 : CSLDestroy(papszMD);
2657 :
2658 : // Attach geolocation array
2659 3 : const std::string osGeolocVRT = std::string(osDatastripSubdirFull)
2660 3 : .append(achSeparator)
2661 3 : .append("GEO_DATA")
2662 3 : .append(achSeparator)
2663 3 : .append(pszGeolocVRT)
2664 6 : .append(".vrt");
2665 3 : CPL_IGNORE_RET_VAL(osDatastripSubdirFull);
2666 : auto poGeolocDS = std::unique_ptr<GDALDataset>(
2667 6 : GDALDataset::Open(osGeolocVRT.c_str(), GDAL_OF_RASTER));
2668 3 : if (poGeolocDS)
2669 : {
2670 6 : CPLStringList osMD(CSLDuplicate(poGeolocDS->GetMetadata()));
2671 3 : osMD.SetNameValue("X_DATASET", osGeolocVRT.c_str());
2672 3 : osMD.SetNameValue("Y_DATASET", osGeolocVRT.c_str());
2673 3 : const char *pszSRS = osMD.FetchNameValue("SRS");
2674 6 : OGRSpatialReference oSRS;
2675 3 : if (oSRS.SetFromUserInput(
2676 : pszSRS,
2677 3 : OGRSpatialReference::SET_FROM_USER_INPUT_LIMITATIONS_get()) ==
2678 : OGRERR_NONE)
2679 : {
2680 3 : osMD.SetNameValue("SRS", oSRS.exportToWkt().c_str());
2681 : }
2682 :
2683 3 : poDS->SetMetadata(osMD.List(), "GEOLOCATION");
2684 : }
2685 :
2686 3 : return poDS.release();
2687 : }
2688 :
2689 : /************************************************************************/
2690 : /* SENTINEL2GetGranuleList_L1CSafeCompact() */
2691 : /************************************************************************/
2692 :
2693 12 : static bool SENTINEL2GetGranuleList_L1CSafeCompact(
2694 : CPLXMLNode *psMainMTD, const char *pszFilename,
2695 : std::vector<L1CSafeCompatGranuleDescription> &osList)
2696 : {
2697 12 : CPLXMLNode *psProductInfo = CPLGetXMLNode(
2698 : psMainMTD, "=Level-1C_User_Product.General_Info.Product_Info");
2699 12 : if (psProductInfo == nullptr)
2700 : {
2701 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
2702 : "=Level-1C_User_Product.General_Info.Product_Info");
2703 0 : return false;
2704 : }
2705 :
2706 : CPLXMLNode *psProductOrganisation =
2707 12 : CPLGetXMLNode(psProductInfo, "Product_Organisation");
2708 12 : if (psProductOrganisation == nullptr)
2709 : {
2710 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
2711 : "Product_Organisation");
2712 0 : return false;
2713 : }
2714 :
2715 24 : CPLString osDirname(CPLGetDirnameSafe(pszFilename));
2716 : #if !defined(_WIN32)
2717 : char szPointerFilename[2048];
2718 : int nBytes = static_cast<int>(
2719 12 : readlink(pszFilename, szPointerFilename, sizeof(szPointerFilename)));
2720 12 : if (nBytes != -1)
2721 : {
2722 : const int nOffset =
2723 0 : std::min(nBytes, static_cast<int>(sizeof(szPointerFilename) - 1));
2724 0 : szPointerFilename[nOffset] = '\0';
2725 0 : osDirname = CPLGetDirnameSafe(szPointerFilename);
2726 : }
2727 : #endif
2728 :
2729 12 : const char chSeparator = SENTINEL2GetPathSeparator(osDirname);
2730 24 : for (CPLXMLNode *psIter = psProductOrganisation->psChild; psIter != nullptr;
2731 12 : psIter = psIter->psNext)
2732 : {
2733 12 : if (psIter->eType != CXT_Element ||
2734 12 : !EQUAL(psIter->pszValue, "Granule_List"))
2735 : {
2736 0 : continue;
2737 : }
2738 24 : for (CPLXMLNode *psIter2 = psIter->psChild; psIter2 != nullptr;
2739 12 : psIter2 = psIter2->psNext)
2740 : {
2741 12 : if (psIter2->eType != CXT_Element ||
2742 12 : !EQUAL(psIter2->pszValue, "Granule"))
2743 : {
2744 0 : continue;
2745 : }
2746 :
2747 : const char *pszImageFile =
2748 12 : CPLGetXMLValue(psIter2, "IMAGE_FILE", nullptr);
2749 12 : if (pszImageFile == nullptr || strlen(pszImageFile) < 3)
2750 : {
2751 0 : CPLDebug("SENTINEL2", "Missing IMAGE_FILE element");
2752 0 : continue;
2753 : }
2754 12 : L1CSafeCompatGranuleDescription oDesc;
2755 12 : oDesc.osBandPrefixPath = osDirname + chSeparator + pszImageFile;
2756 12 : if (CPLHasPathTraversal(oDesc.osBandPrefixPath.c_str()))
2757 : {
2758 0 : CPLError(CE_Failure, CPLE_NotSupported,
2759 : "Path traversal detected in %s",
2760 : oDesc.osBandPrefixPath.c_str());
2761 0 : return false;
2762 : }
2763 12 : oDesc.osBandPrefixPath.resize(oDesc.osBandPrefixPath.size() -
2764 : 3); // strip B12
2765 : // GRANULE/L1C_T30TXT_A007999_20170102T111441/IMG_DATA/T30TXT_20170102T111442_B12
2766 : // --> GRANULE/L1C_T30TXT_A007999_20170102T111441/MTD_TL.xml
2767 : oDesc.osMTDTLPath =
2768 24 : osDirname + chSeparator +
2769 48 : CPLGetDirnameSafe(CPLGetDirnameSafe(pszImageFile).c_str()) +
2770 12 : chSeparator + "MTD_TL.xml";
2771 12 : osList.push_back(std::move(oDesc));
2772 : }
2773 : }
2774 :
2775 12 : return true;
2776 : }
2777 :
2778 : /************************************************************************/
2779 : /* SENTINEL2GetGranuleList_L2ASafeCompact() */
2780 : /************************************************************************/
2781 :
2782 15 : static bool SENTINEL2GetGranuleList_L2ASafeCompact(
2783 : CPLXMLNode *psMainMTD, const char *pszFilename,
2784 : std::vector<L1CSafeCompatGranuleDescription> &osList)
2785 : {
2786 15 : const char *pszNodePath =
2787 : "=Level-2A_User_Product.General_Info.Product_Info";
2788 15 : CPLXMLNode *psProductInfo = CPLGetXMLNode(psMainMTD, pszNodePath);
2789 15 : if (psProductInfo == nullptr)
2790 : {
2791 6 : pszNodePath = "=Level-2A_User_Product.General_Info.L2A_Product_Info";
2792 6 : psProductInfo = CPLGetXMLNode(psMainMTD, pszNodePath);
2793 : }
2794 15 : if (psProductInfo == nullptr)
2795 : {
2796 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s", pszNodePath);
2797 0 : return false;
2798 : }
2799 :
2800 : CPLXMLNode *psProductOrganisation =
2801 15 : CPLGetXMLNode(psProductInfo, "Product_Organisation");
2802 15 : if (psProductOrganisation == nullptr)
2803 : {
2804 : psProductOrganisation =
2805 6 : CPLGetXMLNode(psProductInfo, "L2A_Product_Organisation");
2806 6 : if (psProductOrganisation == nullptr)
2807 : {
2808 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
2809 : "Product_Organisation");
2810 0 : return false;
2811 : }
2812 : }
2813 :
2814 30 : CPLString osDirname(CPLGetDirnameSafe(pszFilename));
2815 : #if !defined(_WIN32)
2816 : char szPointerFilename[2048];
2817 : int nBytes = static_cast<int>(
2818 15 : readlink(pszFilename, szPointerFilename, sizeof(szPointerFilename)));
2819 15 : if (nBytes != -1)
2820 : {
2821 : const int nOffset =
2822 0 : std::min(nBytes, static_cast<int>(sizeof(szPointerFilename) - 1));
2823 0 : szPointerFilename[nOffset] = '\0';
2824 0 : osDirname = CPLGetDirnameSafe(szPointerFilename);
2825 : }
2826 : #endif
2827 :
2828 15 : const char chSeparator = SENTINEL2GetPathSeparator(osDirname);
2829 42 : for (CPLXMLNode *psIter = psProductOrganisation->psChild; psIter != nullptr;
2830 27 : psIter = psIter->psNext)
2831 : {
2832 27 : if (psIter->eType != CXT_Element ||
2833 27 : !EQUAL(psIter->pszValue, "Granule_List"))
2834 : {
2835 0 : continue;
2836 : }
2837 54 : for (CPLXMLNode *psIter2 = psIter->psChild; psIter2 != nullptr;
2838 27 : psIter2 = psIter2->psNext)
2839 : {
2840 27 : if (psIter2->eType != CXT_Element ||
2841 27 : !EQUAL(psIter2->pszValue, "Granule"))
2842 : {
2843 0 : continue;
2844 : }
2845 :
2846 : const char *pszImageFile =
2847 27 : CPLGetXMLValue(psIter2, "IMAGE_FILE", nullptr);
2848 27 : if (pszImageFile == nullptr)
2849 : {
2850 : pszImageFile =
2851 18 : CPLGetXMLValue(psIter2, "IMAGE_FILE_2A", nullptr);
2852 18 : if (pszImageFile == nullptr || strlen(pszImageFile) < 3)
2853 : {
2854 0 : CPLDebug("SENTINEL2", "Missing IMAGE_FILE element");
2855 0 : continue;
2856 : }
2857 : }
2858 27 : L1CSafeCompatGranuleDescription oDesc;
2859 27 : oDesc.osBandPrefixPath = osDirname + chSeparator + pszImageFile;
2860 27 : if (oDesc.osBandPrefixPath.size() < 36)
2861 : {
2862 0 : CPLDebug("SENTINEL2", "Band prefix path too short");
2863 0 : continue;
2864 : }
2865 27 : if (CPLHasPathTraversal(oDesc.osBandPrefixPath.c_str()))
2866 : {
2867 0 : CPLError(CE_Failure, CPLE_NotSupported,
2868 : "Path traversal detected in %s",
2869 : oDesc.osBandPrefixPath.c_str());
2870 0 : return false;
2871 : }
2872 27 : oDesc.osBandPrefixPath.resize(oDesc.osBandPrefixPath.size() - 36);
2873 : // GRANULE/L1C_T30TXT_A007999_20170102T111441/IMG_DATA/T30TXT_20170102T111442_B12_60m
2874 : // --> GRANULE/L1C_T30TXT_A007999_20170102T111441/MTD_TL.xml
2875 : oDesc.osMTDTLPath =
2876 54 : osDirname + chSeparator +
2877 81 : CPLGetDirnameSafe(CPLGetDirnameSafe(pszImageFile).c_str());
2878 27 : if (oDesc.osMTDTLPath.size() < 9)
2879 : {
2880 0 : CPLDebug("SENTINEL2", "MTDTL path too short");
2881 0 : continue;
2882 : }
2883 27 : oDesc.osMTDTLPath.resize(oDesc.osMTDTLPath.size() - 9);
2884 27 : oDesc.osMTDTLPath = oDesc.osMTDTLPath + chSeparator + "MTD_TL.xml";
2885 27 : osList.push_back(std::move(oDesc));
2886 : }
2887 : }
2888 :
2889 15 : return true;
2890 : }
2891 :
2892 : /************************************************************************/
2893 : /* OpenL1C_L2A() */
2894 : /************************************************************************/
2895 :
2896 17 : GDALDataset *SENTINEL2Dataset::OpenL1C_L2A(const char *pszFilename,
2897 : SENTINEL2Level eLevel)
2898 : {
2899 17 : CPLXMLNode *psRoot = CPLParseXMLFile(pszFilename);
2900 17 : if (psRoot == nullptr)
2901 : {
2902 2 : CPLDebug("SENTINEL2", "Cannot XML parse %s", pszFilename);
2903 2 : return nullptr;
2904 : }
2905 :
2906 15 : char *pszOriginalXML = CPLSerializeXMLTree(psRoot);
2907 30 : CPLString osOriginalXML;
2908 15 : if (pszOriginalXML)
2909 15 : osOriginalXML = pszOriginalXML;
2910 15 : CPLFree(pszOriginalXML);
2911 :
2912 30 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRoot);
2913 15 : CPLStripXMLNamespace(psRoot, nullptr, TRUE);
2914 :
2915 15 : const char *pszNodePath =
2916 : (eLevel == SENTINEL2_L1C)
2917 : ? "=Level-1C_User_Product.General_Info.Product_Info"
2918 : : "=Level-2A_User_Product.General_Info.Product_Info";
2919 15 : CPLXMLNode *psProductInfo = CPLGetXMLNode(psRoot, pszNodePath);
2920 15 : if (psProductInfo == nullptr && eLevel == SENTINEL2_L2A)
2921 : {
2922 4 : pszNodePath = "=Level-2A_User_Product.General_Info.L2A_Product_Info";
2923 4 : psProductInfo = CPLGetXMLNode(psRoot, pszNodePath);
2924 : }
2925 15 : if (psProductInfo == nullptr)
2926 : {
2927 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s", pszNodePath);
2928 1 : return nullptr;
2929 : }
2930 :
2931 : const bool bIsSafeCompact =
2932 14 : EQUAL(CPLGetXMLValue(psProductInfo, "Query_Options.PRODUCT_FORMAT", ""),
2933 : "SAFE_COMPACT");
2934 :
2935 28 : std::set<int> oSetResolutions;
2936 28 : std::map<int, std::set<CPLString>> oMapResolutionsToBands;
2937 14 : if (bIsSafeCompact)
2938 : {
2939 70 : for (const auto &sBandDesc : asBandDesc)
2940 : {
2941 : // L2A does not contain B10
2942 65 : if (eLevel == SENTINEL2_L2A &&
2943 39 : strcmp(sBandDesc.pszBandName, "B10") == 0)
2944 3 : continue;
2945 62 : oSetResolutions.insert(sBandDesc.nResolution);
2946 124 : CPLString osName = sBandDesc.pszBandName + 1; /* skip B character */
2947 62 : if (atoi(osName) < 10)
2948 50 : osName = "0" + osName;
2949 62 : oMapResolutionsToBands[sBandDesc.nResolution].insert(
2950 62 : std::move(osName));
2951 : }
2952 5 : if (eLevel == SENTINEL2_L2A)
2953 : {
2954 39 : for (const auto &sL2ABandDesc : asL2ABandDesc)
2955 : {
2956 36 : oSetResolutions.insert(sL2ABandDesc.nResolution);
2957 36 : oMapResolutionsToBands[sL2ABandDesc.nResolution].insert(
2958 36 : sL2ABandDesc.pszBandName);
2959 : }
2960 : }
2961 : }
2962 15 : else if (eLevel == SENTINEL2_L1C &&
2963 6 : !SENTINEL2GetResolutionSet(psProductInfo, oSetResolutions,
2964 : oMapResolutionsToBands))
2965 : {
2966 3 : CPLDebug("SENTINEL2", "Failed to get resolution set");
2967 3 : return nullptr;
2968 : }
2969 :
2970 22 : std::vector<CPLString> aosGranuleList;
2971 11 : if (bIsSafeCompact)
2972 : {
2973 : std::vector<L1CSafeCompatGranuleDescription>
2974 5 : aoL1CSafeCompactGranuleList;
2975 7 : if (eLevel == SENTINEL2_L1C &&
2976 2 : !SENTINEL2GetGranuleList_L1CSafeCompact(
2977 : psRoot, pszFilename, aoL1CSafeCompactGranuleList))
2978 : {
2979 0 : CPLDebug("SENTINEL2", "Failed to get granule list");
2980 0 : return nullptr;
2981 : }
2982 8 : else if (eLevel == SENTINEL2_L2A &&
2983 3 : !SENTINEL2GetGranuleList_L2ASafeCompact(
2984 : psRoot, pszFilename, aoL1CSafeCompactGranuleList))
2985 : {
2986 0 : CPLDebug("SENTINEL2", "Failed to get granule list");
2987 0 : return nullptr;
2988 : }
2989 12 : for (size_t i = 0; i < aoL1CSafeCompactGranuleList.size(); ++i)
2990 : {
2991 7 : aosGranuleList.push_back(
2992 7 : aoL1CSafeCompactGranuleList[i].osMTDTLPath);
2993 : }
2994 : }
2995 6 : else if (!SENTINEL2GetGranuleList(
2996 : psRoot, eLevel, pszFilename, aosGranuleList,
2997 : (eLevel == SENTINEL2_L1C) ? nullptr : &oSetResolutions,
2998 : (eLevel == SENTINEL2_L1C) ? nullptr : &oMapResolutionsToBands))
2999 : {
3000 0 : CPLDebug("SENTINEL2", "Failed to get granule list");
3001 0 : return nullptr;
3002 : }
3003 11 : if (oSetResolutions.empty())
3004 : {
3005 0 : CPLDebug("SENTINEL2", "Resolution set is empty");
3006 0 : return nullptr;
3007 : }
3008 :
3009 11 : std::set<int> oSetEPSGCodes;
3010 27 : for (size_t i = 0; i < aosGranuleList.size(); i++)
3011 : {
3012 16 : int nEPSGCode = 0;
3013 16 : if (SENTINEL2GetGranuleInfo(eLevel, aosGranuleList[i],
3014 32 : *(oSetResolutions.begin()), &nEPSGCode))
3015 : {
3016 13 : oSetEPSGCodes.insert(nEPSGCode);
3017 : }
3018 : }
3019 :
3020 11 : SENTINEL2DatasetContainer *poDS = new SENTINEL2DatasetContainer();
3021 11 : char **papszMD = SENTINEL2GetUserProductMetadata(
3022 : psRoot, (eLevel == SENTINEL2_L1C) ? "Level-1C_User_Product"
3023 : : "Level-2A_User_Product");
3024 11 : poDS->GDALDataset::SetMetadata(papszMD);
3025 11 : CSLDestroy(papszMD);
3026 :
3027 11 : if (!osOriginalXML.empty())
3028 : {
3029 : char *apszXMLMD[2];
3030 11 : apszXMLMD[0] = const_cast<char *>(osOriginalXML.c_str());
3031 11 : apszXMLMD[1] = nullptr;
3032 11 : poDS->GDALDataset::SetMetadata(apszXMLMD, "xml:SENTINEL2");
3033 : }
3034 :
3035 11 : const char *pszPrefix =
3036 : (eLevel == SENTINEL2_L1C) ? "SENTINEL2_L1C" : "SENTINEL2_L2A";
3037 :
3038 : /* Create subdatsets per resolution (10, 20, 60m) and EPSG codes */
3039 11 : int iSubDSNum = 1;
3040 38 : for (std::set<int>::const_iterator oIterRes = oSetResolutions.begin();
3041 65 : oIterRes != oSetResolutions.end(); ++oIterRes)
3042 : {
3043 27 : const int nResolution = *oIterRes;
3044 :
3045 50 : for (std::set<int>::const_iterator oIterEPSG = oSetEPSGCodes.begin();
3046 73 : oIterEPSG != oSetEPSGCodes.end(); ++oIterEPSG)
3047 : {
3048 23 : const int nEPSGCode = *oIterEPSG;
3049 23 : poDS->GDALDataset::SetMetadataItem(
3050 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
3051 : CPLSPrintf("%s:%s:%dm:EPSG_%d", pszPrefix, pszFilename,
3052 : nResolution, nEPSGCode),
3053 : "SUBDATASETS");
3054 :
3055 : CPLString osBandNames = SENTINEL2GetBandListForResolution(
3056 46 : oMapResolutionsToBands[nResolution]);
3057 :
3058 : CPLString osDesc(CPLSPrintf("Bands %s with %dm resolution",
3059 23 : osBandNames.c_str(), nResolution));
3060 23 : if (nEPSGCode >= 32601 && nEPSGCode <= 32660)
3061 23 : osDesc += CPLSPrintf(", UTM %dN", nEPSGCode - 32600);
3062 0 : else if (nEPSGCode >= 32701 && nEPSGCode <= 32760)
3063 0 : osDesc += CPLSPrintf(", UTM %dS", nEPSGCode - 32700);
3064 : else
3065 0 : osDesc += CPLSPrintf(", EPSG:%d", nEPSGCode);
3066 23 : poDS->GDALDataset::SetMetadataItem(
3067 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
3068 : "SUBDATASETS");
3069 :
3070 23 : iSubDSNum++;
3071 : }
3072 : }
3073 :
3074 : /* Expose TCI or PREVIEW subdatasets */
3075 11 : if (bIsSafeCompact)
3076 : {
3077 10 : for (std::set<int>::const_iterator oIterEPSG = oSetEPSGCodes.begin();
3078 15 : oIterEPSG != oSetEPSGCodes.end(); ++oIterEPSG)
3079 : {
3080 5 : const int nEPSGCode = *oIterEPSG;
3081 5 : poDS->GDALDataset::SetMetadataItem(
3082 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
3083 : CPLSPrintf("%s:%s:TCI:EPSG_%d", pszPrefix, pszFilename,
3084 : nEPSGCode),
3085 : "SUBDATASETS");
3086 :
3087 5 : CPLString osDesc("True color image");
3088 5 : if (nEPSGCode >= 32601 && nEPSGCode <= 32660)
3089 5 : osDesc += CPLSPrintf(", UTM %dN", nEPSGCode - 32600);
3090 0 : else if (nEPSGCode >= 32701 && nEPSGCode <= 32760)
3091 0 : osDesc += CPLSPrintf(", UTM %dS", nEPSGCode - 32700);
3092 : else
3093 0 : osDesc += CPLSPrintf(", EPSG:%d", nEPSGCode);
3094 5 : poDS->GDALDataset::SetMetadataItem(
3095 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
3096 : "SUBDATASETS");
3097 :
3098 5 : iSubDSNum++;
3099 : }
3100 : }
3101 : else
3102 : {
3103 10 : for (std::set<int>::const_iterator oIterEPSG = oSetEPSGCodes.begin();
3104 14 : oIterEPSG != oSetEPSGCodes.end(); ++oIterEPSG)
3105 : {
3106 4 : const int nEPSGCode = *oIterEPSG;
3107 4 : poDS->GDALDataset::SetMetadataItem(
3108 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
3109 : CPLSPrintf("%s:%s:PREVIEW:EPSG_%d", pszPrefix, pszFilename,
3110 : nEPSGCode),
3111 : "SUBDATASETS");
3112 :
3113 4 : CPLString osDesc("RGB preview");
3114 4 : if (nEPSGCode >= 32601 && nEPSGCode <= 32660)
3115 4 : osDesc += CPLSPrintf(", UTM %dN", nEPSGCode - 32600);
3116 0 : else if (nEPSGCode >= 32701 && nEPSGCode <= 32760)
3117 0 : osDesc += CPLSPrintf(", UTM %dS", nEPSGCode - 32700);
3118 : else
3119 0 : osDesc += CPLSPrintf(", EPSG:%d", nEPSGCode);
3120 4 : poDS->GDALDataset::SetMetadataItem(
3121 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
3122 : "SUBDATASETS");
3123 :
3124 4 : iSubDSNum++;
3125 : }
3126 : }
3127 :
3128 11 : pszNodePath =
3129 : (eLevel == SENTINEL2_L1C)
3130 : ? "=Level-1C_User_Product.Geometric_Info.Product_Footprint."
3131 : "Product_Footprint.Global_Footprint.EXT_POS_LIST"
3132 : : "=Level-2A_User_Product.Geometric_Info.Product_Footprint."
3133 : "Product_Footprint.Global_Footprint.EXT_POS_LIST";
3134 11 : const char *pszPosList = CPLGetXMLValue(psRoot, pszNodePath, nullptr);
3135 11 : if (pszPosList != nullptr)
3136 : {
3137 18 : CPLString osPolygon = SENTINEL2GetPolygonWKTFromPosList(pszPosList);
3138 9 : if (!osPolygon.empty())
3139 9 : poDS->GDALDataset::SetMetadataItem("FOOTPRINT", osPolygon.c_str());
3140 : }
3141 :
3142 11 : return poDS;
3143 : }
3144 :
3145 : /************************************************************************/
3146 : /* SENTINEL2GetL1BCTileMetadata() */
3147 : /************************************************************************/
3148 :
3149 11 : static char **SENTINEL2GetL1BCTileMetadata(CPLXMLNode *psMainMTD)
3150 : {
3151 22 : CPLStringList aosList;
3152 :
3153 11 : CPLXMLNode *psRoot = CPLGetXMLNode(psMainMTD, "=Level-1C_Tile_ID");
3154 11 : if (psRoot == nullptr)
3155 : {
3156 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find =Level-1C_Tile_ID");
3157 0 : return nullptr;
3158 : }
3159 11 : CPLXMLNode *psGeneralInfo = CPLGetXMLNode(psRoot, "General_Info");
3160 11 : for (CPLXMLNode *psIter =
3161 11 : (psGeneralInfo ? psGeneralInfo->psChild : nullptr);
3162 58 : psIter != nullptr; psIter = psIter->psNext)
3163 : {
3164 47 : if (psIter->eType != CXT_Element)
3165 0 : continue;
3166 47 : const char *pszValue = CPLGetXMLValue(psIter, nullptr, nullptr);
3167 47 : if (pszValue != nullptr)
3168 : {
3169 38 : aosList.AddNameValue(psIter->pszValue, pszValue);
3170 : }
3171 : }
3172 :
3173 11 : CPLXMLNode *psQII = CPLGetXMLNode(psRoot, "Quality_Indicators_Info");
3174 11 : if (psQII != nullptr)
3175 : {
3176 9 : CPLXMLNode *psICCQI = CPLGetXMLNode(psQII, "Image_Content_QI");
3177 9 : for (CPLXMLNode *psIter = (psICCQI ? psICCQI->psChild : nullptr);
3178 27 : psIter != nullptr; psIter = psIter->psNext)
3179 : {
3180 18 : if (psIter->eType != CXT_Element)
3181 0 : continue;
3182 18 : if (psIter->psChild != nullptr &&
3183 18 : psIter->psChild->eType == CXT_Text)
3184 : {
3185 18 : aosList.AddNameValue(psIter->pszValue,
3186 18 : psIter->psChild->pszValue);
3187 : }
3188 : }
3189 : }
3190 :
3191 11 : return aosList.StealList();
3192 : }
3193 :
3194 : /************************************************************************/
3195 : /* OpenL1CTile() */
3196 : /************************************************************************/
3197 :
3198 14 : GDALDataset *SENTINEL2Dataset::OpenL1CTile(const char *pszFilename,
3199 : CPLXMLNode **ppsRootMainMTD,
3200 : int nResolutionOfInterest,
3201 : std::set<CPLString> *poBandSet)
3202 : {
3203 14 : CPLXMLNode *psRoot = CPLParseXMLFile(pszFilename);
3204 14 : if (psRoot == nullptr)
3205 : {
3206 3 : CPLDebug("SENTINEL2", "Cannot XML parse %s", pszFilename);
3207 3 : return nullptr;
3208 : }
3209 :
3210 11 : char *pszOriginalXML = CPLSerializeXMLTree(psRoot);
3211 22 : CPLString osOriginalXML;
3212 11 : if (pszOriginalXML)
3213 11 : osOriginalXML = pszOriginalXML;
3214 11 : CPLFree(pszOriginalXML);
3215 :
3216 22 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRoot);
3217 11 : CPLStripXMLNamespace(psRoot, nullptr, TRUE);
3218 :
3219 22 : std::set<int> oSetResolutions;
3220 22 : std::map<int, std::set<CPLString>> oMapResolutionsToBands;
3221 11 : char **papszMD = nullptr;
3222 11 : SENTINEL2GetResolutionSetAndMainMDFromGranule(
3223 : pszFilename, "Level-1C_User_Product", nResolutionOfInterest,
3224 : oSetResolutions, oMapResolutionsToBands, papszMD, ppsRootMainMTD);
3225 11 : if (poBandSet != nullptr)
3226 8 : *poBandSet = oMapResolutionsToBands[nResolutionOfInterest];
3227 :
3228 11 : SENTINEL2DatasetContainer *poDS = new SENTINEL2DatasetContainer();
3229 :
3230 11 : char **papszGranuleMD = SENTINEL2GetL1BCTileMetadata(psRoot);
3231 11 : papszMD = CSLMerge(papszMD, papszGranuleMD);
3232 11 : CSLDestroy(papszGranuleMD);
3233 :
3234 : // Remove CLOUD_COVERAGE_ASSESSMENT that comes from main metadata, if
3235 : // granule CLOUDY_PIXEL_PERCENTAGE is present.
3236 20 : if (CSLFetchNameValue(papszMD, "CLOUDY_PIXEL_PERCENTAGE") != nullptr &&
3237 9 : CSLFetchNameValue(papszMD, "CLOUD_COVERAGE_ASSESSMENT") != nullptr)
3238 : {
3239 7 : papszMD =
3240 7 : CSLSetNameValue(papszMD, "CLOUD_COVERAGE_ASSESSMENT", nullptr);
3241 : }
3242 :
3243 11 : poDS->GDALDataset::SetMetadata(papszMD);
3244 11 : CSLDestroy(papszMD);
3245 :
3246 11 : if (!osOriginalXML.empty())
3247 : {
3248 : char *apszXMLMD[2];
3249 11 : apszXMLMD[0] = const_cast<char *>(osOriginalXML.c_str());
3250 11 : apszXMLMD[1] = nullptr;
3251 11 : poDS->GDALDataset::SetMetadata(apszXMLMD, "xml:SENTINEL2");
3252 : }
3253 :
3254 : /* Create subdatsets per resolution (10, 20, 60m) */
3255 11 : int iSubDSNum = 1;
3256 36 : for (std::set<int>::const_iterator oIterRes = oSetResolutions.begin();
3257 61 : oIterRes != oSetResolutions.end(); ++oIterRes)
3258 : {
3259 25 : const int nResolution = *oIterRes;
3260 :
3261 25 : poDS->GDALDataset::SetMetadataItem(
3262 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
3263 : CPLSPrintf("%s:%s:%dm", "SENTINEL2_L1C_TILE", pszFilename,
3264 : nResolution),
3265 : "SUBDATASETS");
3266 :
3267 : CPLString osBandNames = SENTINEL2GetBandListForResolution(
3268 50 : oMapResolutionsToBands[nResolution]);
3269 :
3270 : CPLString osDesc(CPLSPrintf("Bands %s with %dm resolution",
3271 25 : osBandNames.c_str(), nResolution));
3272 25 : poDS->GDALDataset::SetMetadataItem(
3273 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
3274 : "SUBDATASETS");
3275 :
3276 25 : iSubDSNum++;
3277 : }
3278 :
3279 : /* Expose PREVIEW subdataset */
3280 11 : poDS->GDALDataset::SetMetadataItem(
3281 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
3282 : CPLSPrintf("%s:%s:PREVIEW", "SENTINEL2_L1C_TILE", pszFilename),
3283 : "SUBDATASETS");
3284 :
3285 11 : CPLString osDesc("RGB preview");
3286 11 : poDS->GDALDataset::SetMetadataItem(
3287 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
3288 : "SUBDATASETS");
3289 :
3290 11 : return poDS;
3291 : }
3292 :
3293 : /************************************************************************/
3294 : /* SENTINEL2GetOption() */
3295 : /************************************************************************/
3296 :
3297 55 : static const char *SENTINEL2GetOption(GDALOpenInfo *poOpenInfo,
3298 : const char *pszName,
3299 : const char *pszDefaultVal)
3300 : {
3301 : #ifdef GDAL_DMD_OPENOPTIONLIST
3302 : const char *pszVal =
3303 55 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, pszName);
3304 55 : if (pszVal != nullptr)
3305 3 : return pszVal;
3306 : #else
3307 : (void)poOpenInfo;
3308 : #endif
3309 52 : return CPLGetConfigOption(CPLSPrintf("SENTINEL2_%s", pszName),
3310 52 : pszDefaultVal);
3311 : }
3312 :
3313 : /************************************************************************/
3314 : /* LaunderUnit() */
3315 : /************************************************************************/
3316 :
3317 143 : static CPLString LaunderUnit(const char *pszUnit)
3318 : {
3319 143 : CPLString osUnit;
3320 1144 : for (int i = 0; pszUnit[i] != '\0';)
3321 : {
3322 1001 : if (strncmp(pszUnit + i,
3323 : "\xC2"
3324 : "\xB2",
3325 : 2) == 0) /* square / 2 */
3326 : {
3327 143 : i += 2;
3328 143 : osUnit += "2";
3329 : }
3330 858 : else if (strncmp(pszUnit + i,
3331 : "\xC2"
3332 : "\xB5",
3333 : 2) == 0) /* micro */
3334 : {
3335 143 : i += 2;
3336 143 : osUnit += "u";
3337 : }
3338 : else
3339 : {
3340 715 : osUnit += pszUnit[i];
3341 715 : i++;
3342 : }
3343 : }
3344 143 : return osUnit;
3345 : }
3346 :
3347 : /************************************************************************/
3348 : /* SENTINEL2GetTileInfo() */
3349 : /************************************************************************/
3350 :
3351 38 : static bool SENTINEL2GetTileInfo(const char *pszFilename, int *pnWidth,
3352 : int *pnHeight, int *pnBits)
3353 : {
3354 : static const unsigned char jp2_box_jp[] = {0x6a, 0x50, 0x20,
3355 : 0x20}; /* 'jP ' */
3356 38 : VSILFILE *fp = VSIFOpenL(pszFilename, "rb");
3357 38 : if (fp == nullptr)
3358 2 : return false;
3359 : GByte abyHeader[8];
3360 36 : if (VSIFReadL(abyHeader, 8, 1, fp) != 1)
3361 : {
3362 0 : VSIFCloseL(fp);
3363 0 : return false;
3364 : }
3365 36 : if (memcmp(abyHeader + 4, jp2_box_jp, 4) == 0)
3366 : {
3367 3 : bool bRet = false;
3368 : /* Just parse the ihdr box instead of doing a full dataset opening */
3369 3 : GDALJP2Box oBox(fp);
3370 3 : if (oBox.ReadFirst())
3371 : {
3372 9 : while (strlen(oBox.GetType()) > 0)
3373 : {
3374 9 : if (EQUAL(oBox.GetType(), "jp2h"))
3375 : {
3376 6 : GDALJP2Box oChildBox(fp);
3377 3 : if (!oChildBox.ReadFirstChild(&oBox))
3378 0 : break;
3379 :
3380 3 : while (strlen(oChildBox.GetType()) > 0)
3381 : {
3382 3 : if (EQUAL(oChildBox.GetType(), "ihdr"))
3383 : {
3384 3 : GByte *pabyData = oChildBox.ReadBoxData();
3385 3 : GIntBig nLength = oChildBox.GetDataLength();
3386 3 : if (pabyData != nullptr && nLength >= 4 + 4 + 2 + 1)
3387 : {
3388 3 : bRet = true;
3389 3 : if (pnHeight)
3390 : {
3391 1 : memcpy(pnHeight, pabyData, 4);
3392 1 : CPL_MSBPTR32(pnHeight);
3393 : }
3394 3 : if (pnWidth != nullptr)
3395 : {
3396 : // cppcheck-suppress nullPointer
3397 1 : memcpy(pnWidth, pabyData + 4, 4);
3398 1 : CPL_MSBPTR32(pnWidth);
3399 : }
3400 3 : if (pnBits)
3401 : {
3402 3 : GByte byPBC = pabyData[4 + 4 + 2];
3403 3 : if (byPBC != 255)
3404 : {
3405 3 : *pnBits = 1 + (byPBC & 0x7f);
3406 : }
3407 : else
3408 0 : *pnBits = 0;
3409 : }
3410 : }
3411 3 : CPLFree(pabyData);
3412 3 : break;
3413 : }
3414 0 : if (!oChildBox.ReadNextChild(&oBox))
3415 0 : break;
3416 : }
3417 3 : break;
3418 : }
3419 :
3420 6 : if (!oBox.ReadNext())
3421 0 : break;
3422 : }
3423 : }
3424 3 : VSIFCloseL(fp);
3425 3 : return bRet;
3426 : }
3427 : else /* for unit tests, we use TIFF */
3428 : {
3429 33 : VSIFCloseL(fp);
3430 : auto poDS = std::unique_ptr<GDALDataset>(GDALDataset::Open(
3431 33 : pszFilename, GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR));
3432 33 : bool bRet = false;
3433 33 : if (poDS != nullptr)
3434 : {
3435 33 : if (poDS->GetRasterCount() != 0)
3436 : {
3437 33 : bRet = true;
3438 33 : if (pnWidth)
3439 3 : *pnWidth = poDS->GetRasterXSize();
3440 33 : if (pnHeight)
3441 3 : *pnHeight = poDS->GetRasterYSize();
3442 33 : if (pnBits)
3443 : {
3444 : const char *pszNBits =
3445 33 : poDS->GetRasterBand(1)->GetMetadataItem(
3446 33 : "NBITS", "IMAGE_STRUCTURE");
3447 33 : if (pszNBits == nullptr)
3448 : {
3449 : GDALDataType eDT =
3450 7 : poDS->GetRasterBand(1)->GetRasterDataType();
3451 : pszNBits =
3452 7 : CPLSPrintf("%d", GDALGetDataTypeSizeBits(eDT));
3453 : }
3454 33 : *pnBits = atoi(pszNBits);
3455 : }
3456 : }
3457 : }
3458 33 : return bRet;
3459 : }
3460 : }
3461 :
3462 : /************************************************************************/
3463 : /* OpenL1C_L2ASubdataset() */
3464 : /************************************************************************/
3465 :
3466 79 : GDALDataset *SENTINEL2Dataset::OpenL1C_L2ASubdataset(GDALOpenInfo *poOpenInfo,
3467 : SENTINEL2Level eLevel)
3468 : {
3469 158 : CPLString osFilename;
3470 79 : const char *pszPrefix =
3471 : (eLevel == SENTINEL2_L1C) ? "SENTINEL2_L1C" : "SENTINEL2_L2A";
3472 79 : CPLAssert(STARTS_WITH_CI(poOpenInfo->pszFilename, pszPrefix));
3473 79 : osFilename = poOpenInfo->pszFilename + strlen(pszPrefix) + 1;
3474 79 : const char *pszColumn = strrchr(osFilename.c_str(), ':');
3475 79 : if (pszColumn == nullptr || pszColumn == osFilename.c_str())
3476 : {
3477 10 : CPLError(CE_Failure, CPLE_AppDefined,
3478 : "Invalid syntax for %s:", pszPrefix);
3479 10 : return nullptr;
3480 : }
3481 69 : const auto nPos = static_cast<size_t>(pszColumn - osFilename.c_str());
3482 138 : const std::string osEPSGCode = osFilename.substr(nPos + 1);
3483 69 : osFilename.resize(nPos);
3484 69 : const char *pszPrecision = strrchr(osFilename.c_str(), ':');
3485 69 : if (pszPrecision == nullptr)
3486 : {
3487 10 : CPLError(CE_Failure, CPLE_AppDefined,
3488 : "Invalid syntax for %s:", pszPrefix);
3489 10 : return nullptr;
3490 : }
3491 :
3492 59 : if (!STARTS_WITH_CI(osEPSGCode.c_str(), "EPSG_"))
3493 : {
3494 5 : CPLError(CE_Failure, CPLE_AppDefined,
3495 : "Invalid syntax for %s:", pszPrefix);
3496 5 : return nullptr;
3497 : }
3498 :
3499 54 : const int nSubDSEPSGCode = atoi(osEPSGCode.c_str() + strlen("EPSG_"));
3500 54 : const bool bIsPreview = STARTS_WITH_CI(pszPrecision + 1, "PREVIEW");
3501 54 : const bool bIsTCI = STARTS_WITH_CI(pszPrecision + 1, "TCI");
3502 105 : const int nSubDSPrecision = (bIsPreview) ? 320
3503 51 : : (bIsTCI) ? RES_10M
3504 49 : : atoi(pszPrecision + 1);
3505 54 : if (!bIsTCI && !bIsPreview && nSubDSPrecision != RES_10M &&
3506 20 : nSubDSPrecision != RES_20M && nSubDSPrecision != RES_60M)
3507 : {
3508 5 : CPLError(CE_Failure, CPLE_NotSupported, "Unsupported precision: %d",
3509 : nSubDSPrecision);
3510 5 : return nullptr;
3511 : }
3512 :
3513 49 : osFilename.resize(pszPrecision - osFilename.c_str());
3514 98 : std::vector<CPLString> aosNonJP2Files;
3515 49 : aosNonJP2Files.push_back(osFilename);
3516 :
3517 49 : CPLXMLNode *psRoot = CPLParseXMLFile(osFilename);
3518 49 : if (psRoot == nullptr)
3519 : {
3520 7 : CPLDebug("SENTINEL2", "Cannot XML parse %s", osFilename.c_str());
3521 7 : return nullptr;
3522 : }
3523 :
3524 42 : char *pszOriginalXML = CPLSerializeXMLTree(psRoot);
3525 84 : CPLString osOriginalXML;
3526 42 : if (pszOriginalXML)
3527 42 : osOriginalXML = pszOriginalXML;
3528 42 : CPLFree(pszOriginalXML);
3529 :
3530 84 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRoot);
3531 42 : CPLStripXMLNamespace(psRoot, nullptr, TRUE);
3532 :
3533 : CPLXMLNode *psProductInfo =
3534 : eLevel == SENTINEL2_L1C
3535 42 : ? CPLGetXMLNode(psRoot,
3536 : "=Level-1C_User_Product.General_Info.Product_Info")
3537 18 : : CPLGetXMLNode(psRoot,
3538 42 : "=Level-2A_User_Product.General_Info.Product_Info");
3539 42 : if (psProductInfo == nullptr && eLevel == SENTINEL2_L2A)
3540 : {
3541 11 : psProductInfo = CPLGetXMLNode(
3542 : psRoot, "=Level-2A_User_Product.General_Info.L2A_Product_Info");
3543 : }
3544 42 : if (psProductInfo == nullptr)
3545 : {
3546 1 : CPLDebug("SENTINEL2", "Product Info not found");
3547 1 : return nullptr;
3548 : }
3549 :
3550 : const bool bIsSafeCompact =
3551 41 : EQUAL(CPLGetXMLValue(psProductInfo, "Query_Options.PRODUCT_FORMAT", ""),
3552 : "SAFE_COMPACT");
3553 :
3554 : const char *pszProductURI =
3555 41 : CPLGetXMLValue(psProductInfo, "PRODUCT_URI", nullptr);
3556 41 : SENTINEL2ProductType pType = MSI2A;
3557 41 : if (pszProductURI == nullptr)
3558 : {
3559 : pszProductURI =
3560 32 : CPLGetXMLValue(psProductInfo, "PRODUCT_URI_2A", nullptr);
3561 32 : pType = MSI2Ap;
3562 : }
3563 41 : if (pszProductURI == nullptr)
3564 27 : pszProductURI = "";
3565 :
3566 82 : std::vector<CPLString> aosGranuleList;
3567 82 : std::map<int, std::set<CPLString>> oMapResolutionsToBands;
3568 82 : std::vector<L1CSafeCompatGranuleDescription> aoL1CSafeCompactGranuleList;
3569 41 : if (bIsSafeCompact)
3570 : {
3571 308 : for (const auto &sBandDesc : asBandDesc)
3572 : {
3573 : // L2 does not contain B10
3574 286 : if (eLevel == SENTINEL2_L2A &&
3575 156 : strcmp(sBandDesc.pszBandName, "B10") == 0)
3576 12 : continue;
3577 548 : CPLString osName = sBandDesc.pszBandName + 1; /* skip B character */
3578 274 : if (atoi(osName) < 10)
3579 220 : osName = "0" + osName;
3580 274 : oMapResolutionsToBands[sBandDesc.nResolution].insert(
3581 274 : std::move(osName));
3582 : }
3583 22 : if (eLevel == SENTINEL2_L2A)
3584 : {
3585 156 : for (const auto &sL2ABandDesc : asL2ABandDesc)
3586 : {
3587 144 : oMapResolutionsToBands[sL2ABandDesc.nResolution].insert(
3588 144 : sL2ABandDesc.pszBandName);
3589 : }
3590 : }
3591 32 : if (eLevel == SENTINEL2_L1C &&
3592 10 : !SENTINEL2GetGranuleList_L1CSafeCompact(
3593 : psRoot, osFilename, aoL1CSafeCompactGranuleList))
3594 : {
3595 0 : CPLDebug("SENTINEL2", "Failed to get granule list");
3596 0 : return nullptr;
3597 : }
3598 34 : if (eLevel == SENTINEL2_L2A &&
3599 12 : !SENTINEL2GetGranuleList_L2ASafeCompact(
3600 : psRoot, osFilename, aoL1CSafeCompactGranuleList))
3601 : {
3602 0 : CPLDebug("SENTINEL2", "Failed to get granule list");
3603 0 : return nullptr;
3604 : }
3605 54 : for (size_t i = 0; i < aoL1CSafeCompactGranuleList.size(); ++i)
3606 : {
3607 32 : aosGranuleList.push_back(
3608 32 : aoL1CSafeCompactGranuleList[i].osMTDTLPath);
3609 : }
3610 : }
3611 19 : else if (!SENTINEL2GetGranuleList(
3612 : psRoot, eLevel, osFilename, aosGranuleList, nullptr,
3613 : (eLevel == SENTINEL2_L1C) ? nullptr : &oMapResolutionsToBands))
3614 : {
3615 1 : CPLDebug("SENTINEL2", "Failed to get granule list");
3616 1 : return nullptr;
3617 : }
3618 :
3619 80 : std::vector<CPLString> aosBands;
3620 80 : std::set<CPLString> oSetBands;
3621 40 : if (bIsPreview || bIsTCI)
3622 : {
3623 5 : aosBands.push_back("04");
3624 5 : aosBands.push_back("03");
3625 5 : aosBands.push_back("02");
3626 : }
3627 35 : else if (eLevel == SENTINEL2_L1C && !bIsSafeCompact)
3628 : {
3629 : CPLXMLNode *psBandList =
3630 10 : CPLGetXMLNode(psRoot, "=Level-1C_User_Product.General_Info.Product_"
3631 : "Info.Query_Options.Band_List");
3632 10 : if (psBandList == nullptr)
3633 : {
3634 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
3635 : "Query_Options.Band_List");
3636 0 : return nullptr;
3637 : }
3638 :
3639 116 : for (CPLXMLNode *psIter = psBandList->psChild; psIter != nullptr;
3640 106 : psIter = psIter->psNext)
3641 : {
3642 106 : if (psIter->eType != CXT_Element ||
3643 106 : !EQUAL(psIter->pszValue, "BAND_NAME"))
3644 72 : continue;
3645 106 : const char *pszBandName = CPLGetXMLValue(psIter, nullptr, "");
3646 : const SENTINEL2BandDescription *psBandDesc =
3647 106 : SENTINEL2GetBandDesc(pszBandName);
3648 106 : if (psBandDesc == nullptr)
3649 : {
3650 0 : CPLDebug("SENTINEL2", "Unknown band name %s", pszBandName);
3651 0 : continue;
3652 : }
3653 106 : if (psBandDesc->nResolution != nSubDSPrecision)
3654 72 : continue;
3655 : CPLString osName =
3656 68 : psBandDesc->pszBandName + 1; /* skip B character */
3657 34 : if (atoi(osName) < 10)
3658 30 : osName = "0" + osName;
3659 34 : oSetBands.insert(std::move(osName));
3660 : }
3661 10 : if (oSetBands.empty())
3662 : {
3663 0 : CPLDebug("SENTINEL2", "Band set is empty");
3664 0 : return nullptr;
3665 10 : }
3666 : }
3667 : else
3668 : {
3669 25 : oSetBands = oMapResolutionsToBands[nSubDSPrecision];
3670 : }
3671 :
3672 40 : if (aosBands.empty())
3673 : {
3674 192 : for (std::set<CPLString>::const_iterator oIterBandnames =
3675 35 : oSetBands.begin();
3676 419 : oIterBandnames != oSetBands.end(); ++oIterBandnames)
3677 : {
3678 192 : aosBands.push_back(*oIterBandnames);
3679 : }
3680 : /* Put 2=Blue, 3=Green, 4=Band bands in RGB order for conveniency */
3681 82 : if (aosBands.size() >= 3 && aosBands[0] == "02" &&
3682 82 : aosBands[1] == "03" && aosBands[2] == "04")
3683 : {
3684 17 : aosBands[0] = "04";
3685 17 : aosBands[2] = "02";
3686 : }
3687 : }
3688 :
3689 : /* -------------------------------------------------------------------- */
3690 : /* Create dataset. */
3691 : /* -------------------------------------------------------------------- */
3692 :
3693 40 : char **papszMD = SENTINEL2GetUserProductMetadata(
3694 : psRoot, (eLevel == SENTINEL2_L1C) ? "Level-1C_User_Product"
3695 : : "Level-2A_User_Product");
3696 :
3697 : const int nSaturatedVal =
3698 40 : atoi(CSLFetchNameValueDef(papszMD, "SPECIAL_VALUE_SATURATED", "-1"));
3699 : const int nNodataVal =
3700 40 : atoi(CSLFetchNameValueDef(papszMD, "SPECIAL_VALUE_NODATA", "-1"));
3701 :
3702 : const bool bAlpha =
3703 40 : CPLTestBool(SENTINEL2GetOption(poOpenInfo, "ALPHA", "FALSE"));
3704 :
3705 40 : SENTINEL2Dataset *poDS = CreateL1CL2ADataset(
3706 : eLevel, pType, bIsSafeCompact, aosGranuleList,
3707 : aoL1CSafeCompactGranuleList, aosNonJP2Files, nSubDSPrecision,
3708 : bIsPreview, bIsTCI, nSubDSEPSGCode, bAlpha, aosBands, nSaturatedVal,
3709 80 : nNodataVal, CPLString(pszProductURI));
3710 40 : if (poDS == nullptr)
3711 : {
3712 12 : CSLDestroy(papszMD);
3713 12 : return nullptr;
3714 : }
3715 :
3716 28 : if (!osOriginalXML.empty())
3717 : {
3718 28 : char *apszXMLMD[2] = {nullptr};
3719 28 : apszXMLMD[0] = const_cast<char *>(osOriginalXML.c_str());
3720 28 : apszXMLMD[1] = nullptr;
3721 28 : poDS->GDALDataset::SetMetadata(apszXMLMD, "xml:SENTINEL2");
3722 : }
3723 :
3724 28 : poDS->GDALDataset::SetMetadata(papszMD);
3725 28 : CSLDestroy(papszMD);
3726 :
3727 : /* -------------------------------------------------------------------- */
3728 : /* Add extra band metadata. */
3729 : /* -------------------------------------------------------------------- */
3730 28 : poDS->AddL1CL2ABandMetadata(eLevel, psRoot, aosBands);
3731 :
3732 : /* -------------------------------------------------------------------- */
3733 : /* Initialize overview information. */
3734 : /* -------------------------------------------------------------------- */
3735 28 : poDS->SetDescription(poOpenInfo->pszFilename);
3736 28 : CPLString osOverviewFile;
3737 28 : if (bIsPreview)
3738 : osOverviewFile = CPLSPrintf("%s_PREVIEW_EPSG_%d.tif.ovr",
3739 3 : osFilename.c_str(), nSubDSEPSGCode);
3740 25 : else if (bIsTCI)
3741 : osOverviewFile = CPLSPrintf("%s_TCI_EPSG_%d.tif.ovr",
3742 2 : osFilename.c_str(), nSubDSEPSGCode);
3743 : else
3744 : osOverviewFile =
3745 : CPLSPrintf("%s_%dm_EPSG_%d.tif.ovr", osFilename.c_str(),
3746 23 : nSubDSPrecision, nSubDSEPSGCode);
3747 28 : poDS->SetMetadataItem("OVERVIEW_FILE", osOverviewFile, "OVERVIEWS");
3748 28 : poDS->oOvManager.Initialize(poDS, ":::VIRTUAL:::");
3749 :
3750 28 : return poDS;
3751 : }
3752 :
3753 : /************************************************************************/
3754 : /* AddL1CL2ABandMetadata() */
3755 : /************************************************************************/
3756 :
3757 34 : void SENTINEL2Dataset::AddL1CL2ABandMetadata(
3758 : SENTINEL2Level eLevel, CPLXMLNode *psRoot,
3759 : const std::vector<CPLString> &aosBands)
3760 : {
3761 : CPLXMLNode *psIC =
3762 34 : CPLGetXMLNode(psRoot, (eLevel == SENTINEL2_L1C)
3763 : ? "=Level-1C_User_Product.General_Info."
3764 : "Product_Image_Characteristics"
3765 : : "=Level-2A_User_Product.General_Info."
3766 : "Product_Image_Characteristics");
3767 34 : if (psIC == nullptr)
3768 : {
3769 8 : psIC = CPLGetXMLNode(psRoot, "=Level-2A_User_Product.General_Info.L2A_"
3770 : "Product_Image_Characteristics");
3771 : }
3772 34 : if (psIC != nullptr)
3773 : {
3774 : CPLXMLNode *psSIL =
3775 32 : CPLGetXMLNode(psIC, "Reflectance_Conversion.Solar_Irradiance_List");
3776 32 : if (psSIL != nullptr)
3777 : {
3778 448 : for (CPLXMLNode *psIter = psSIL->psChild; psIter != nullptr;
3779 416 : psIter = psIter->psNext)
3780 : {
3781 416 : if (psIter->eType != CXT_Element ||
3782 416 : !EQUAL(psIter->pszValue, "SOLAR_IRRADIANCE"))
3783 : {
3784 0 : continue;
3785 : }
3786 : const char *pszBandId =
3787 416 : CPLGetXMLValue(psIter, "bandId", nullptr);
3788 416 : const char *pszUnit = CPLGetXMLValue(psIter, "unit", nullptr);
3789 416 : const char *pszValue = CPLGetXMLValue(psIter, nullptr, nullptr);
3790 416 : if (pszBandId != nullptr && pszUnit != nullptr &&
3791 : pszValue != nullptr)
3792 : {
3793 416 : const unsigned int nIdx = atoi(pszBandId);
3794 416 : if (nIdx < CPL_ARRAYSIZE(asBandDesc))
3795 : {
3796 2089 : for (int i = 0; i < nBands; i++)
3797 : {
3798 1816 : GDALRasterBand *poBand = GetRasterBand(i + 1);
3799 : const char *pszBandName =
3800 1816 : poBand->GetMetadataItem("BANDNAME");
3801 1816 : if (pszBandName != nullptr &&
3802 1806 : EQUAL(asBandDesc[nIdx].pszBandName,
3803 : pszBandName))
3804 : {
3805 143 : poBand->GDALRasterBand::SetMetadataItem(
3806 : "SOLAR_IRRADIANCE", pszValue);
3807 143 : poBand->GDALRasterBand::SetMetadataItem(
3808 : "SOLAR_IRRADIANCE_UNIT",
3809 286 : LaunderUnit(pszUnit));
3810 143 : break;
3811 : }
3812 : }
3813 : }
3814 : }
3815 : }
3816 : }
3817 :
3818 32 : CPLXMLNode *psOL = CPLGetXMLNode(
3819 : psIC, (eLevel == SENTINEL2_L1C) ? "Radiometric_Offset_List"
3820 : : "BOA_ADD_OFFSET_VALUES_LIST");
3821 32 : if (psOL != nullptr)
3822 : {
3823 56 : for (CPLXMLNode *psIter = psOL->psChild; psIter != nullptr;
3824 52 : psIter = psIter->psNext)
3825 : {
3826 52 : if (psIter->eType != CXT_Element ||
3827 52 : !EQUAL(psIter->pszValue, (eLevel == SENTINEL2_L1C)
3828 : ? "RADIO_ADD_OFFSET"
3829 : : "BOA_ADD_OFFSET"))
3830 : {
3831 0 : continue;
3832 : }
3833 : const char *pszBandId =
3834 52 : CPLGetXMLValue(psIter, "band_id", nullptr);
3835 52 : const char *pszValue = CPLGetXMLValue(psIter, nullptr, nullptr);
3836 52 : if (pszBandId != nullptr && pszValue != nullptr)
3837 : {
3838 52 : const unsigned int nIdx = atoi(pszBandId);
3839 52 : if (nIdx < CPL_ARRAYSIZE(asBandDesc))
3840 : {
3841 303 : for (int i = 0; i < nBands; i++)
3842 : {
3843 271 : GDALRasterBand *poBand = GetRasterBand(i + 1);
3844 : const char *pszBandName =
3845 271 : poBand->GetMetadataItem("BANDNAME");
3846 271 : if (pszBandName != nullptr &&
3847 271 : EQUAL(asBandDesc[nIdx].pszBandName,
3848 : pszBandName))
3849 : {
3850 20 : poBand->GDALRasterBand::SetMetadataItem(
3851 : (eLevel == SENTINEL2_L1C)
3852 : ? "RADIO_ADD_OFFSET"
3853 : : "BOA_ADD_OFFSET",
3854 : pszValue);
3855 20 : break;
3856 : }
3857 : }
3858 : }
3859 : }
3860 : }
3861 : }
3862 : }
3863 :
3864 : /* -------------------------------------------------------------------- */
3865 : /* Add Scene Classification category values (L2A) */
3866 : /* -------------------------------------------------------------------- */
3867 34 : CPLXMLNode *psSCL = CPLGetXMLNode(
3868 : psRoot, "=Level-2A_User_Product.General_Info."
3869 : "Product_Image_Characteristics.Scene_Classification_List");
3870 34 : if (psSCL == nullptr)
3871 : {
3872 29 : psSCL = CPLGetXMLNode(
3873 : psRoot,
3874 : "=Level-2A_User_Product.General_Info."
3875 : "L2A_Product_Image_Characteristics.L2A_Scene_Classification_List");
3876 : }
3877 34 : int nSCLBand = 0;
3878 199 : for (int nBand = 1; nBand <= static_cast<int>(aosBands.size()); nBand++)
3879 : {
3880 172 : if (EQUAL(aosBands[nBand - 1], "SCL"))
3881 : {
3882 7 : nSCLBand = nBand;
3883 7 : break;
3884 : }
3885 : }
3886 34 : if (psSCL != nullptr && nSCLBand > 0)
3887 : {
3888 14 : std::vector<std::string> aosCategories(100);
3889 7 : size_t nMaxIdx = 0;
3890 7 : bool bFirst = true;
3891 91 : for (CPLXMLNode *psIter = psSCL->psChild; psIter != nullptr;
3892 84 : psIter = psIter->psNext)
3893 : {
3894 84 : if (psIter->eType != CXT_Element ||
3895 84 : (!EQUAL(psIter->pszValue, "L2A_Scene_Classification_ID") &&
3896 36 : !EQUAL(psIter->pszValue, "Scene_Classification_ID")))
3897 : {
3898 0 : continue;
3899 : }
3900 : const char *pszText =
3901 84 : CPLGetXMLValue(psIter, "SCENE_CLASSIFICATION_TEXT", nullptr);
3902 84 : if (pszText == nullptr)
3903 48 : pszText = CPLGetXMLValue(
3904 : psIter, "L2A_SCENE_CLASSIFICATION_TEXT", nullptr);
3905 : const char *pszIdx =
3906 84 : CPLGetXMLValue(psIter, "SCENE_CLASSIFICATION_INDEX", nullptr);
3907 84 : if (pszIdx == nullptr)
3908 48 : pszIdx = CPLGetXMLValue(
3909 : psIter, "L2A_SCENE_CLASSIFICATION_INDEX", nullptr);
3910 84 : if (pszText && pszIdx)
3911 : {
3912 84 : const size_t nIdx = atoi(pszIdx);
3913 84 : if (nIdx < aosCategories.size())
3914 : {
3915 84 : if (bFirst)
3916 7 : nMaxIdx = nIdx;
3917 : else
3918 77 : nMaxIdx = std::max(nMaxIdx, nIdx);
3919 84 : bFirst = false;
3920 84 : if (STARTS_WITH_CI(pszText, "SC_"))
3921 84 : aosCategories[nIdx] = pszText + 3;
3922 : else
3923 0 : aosCategories[nIdx] = pszText;
3924 : }
3925 : }
3926 : }
3927 7 : if (!bFirst)
3928 : {
3929 7 : aosCategories.resize(nMaxIdx + 1);
3930 14 : GetRasterBand(nSCLBand)->SetCategoryNames(
3931 14 : CPLStringList(aosCategories).List());
3932 : }
3933 : }
3934 34 : }
3935 :
3936 : /************************************************************************/
3937 : /* CreateL1CL2ADataset() */
3938 : /************************************************************************/
3939 :
3940 48 : SENTINEL2Dataset *SENTINEL2Dataset::CreateL1CL2ADataset(
3941 : SENTINEL2Level eLevel, SENTINEL2ProductType pType, bool bIsSafeCompact,
3942 : const std::vector<CPLString> &aosGranuleList,
3943 : const std::vector<L1CSafeCompatGranuleDescription>
3944 : &aoL1CSafeCompactGranuleList,
3945 : std::vector<CPLString> &aosNonJP2Files, int nSubDSPrecision,
3946 : bool bIsPreview, bool bIsTCI,
3947 : int nSubDSEPSGCode, /* or -1 if not known at this point */
3948 : bool bAlpha, const std::vector<CPLString> &aosBands, int nSaturatedVal,
3949 : int nNodataVal, const CPLString &osProductURI)
3950 : {
3951 :
3952 : /* Iterate over granule metadata to know the layer extent */
3953 : /* and the location of each granule */
3954 48 : double dfMinX = 1.0e20;
3955 48 : double dfMinY = 1.0e20;
3956 48 : double dfMaxX = -1.0e20;
3957 48 : double dfMaxY = -1.0e20;
3958 96 : std::vector<SENTINEL2GranuleInfo> aosGranuleInfoList;
3959 48 : const int nDesiredResolution = (bIsPreview || bIsTCI) ? 0 : nSubDSPrecision;
3960 :
3961 48 : if (bIsSafeCompact)
3962 : {
3963 22 : CPLAssert(aosGranuleList.size() == aoL1CSafeCompactGranuleList.size());
3964 : }
3965 :
3966 116 : for (size_t i = 0; i < aosGranuleList.size(); i++)
3967 : {
3968 68 : int nEPSGCode = 0;
3969 68 : double dfULX = 0.0;
3970 68 : double dfULY = 0.0;
3971 68 : int nResolution = 0;
3972 68 : int nWidth = 0;
3973 68 : int nHeight = 0;
3974 68 : if (SENTINEL2GetGranuleInfo(eLevel, aosGranuleList[i],
3975 : nDesiredResolution, &nEPSGCode, &dfULX,
3976 65 : &dfULY, &nResolution, &nWidth, &nHeight) &&
3977 117 : (nSubDSEPSGCode == nEPSGCode || nSubDSEPSGCode < 0) &&
3978 49 : nResolution != 0)
3979 : {
3980 48 : nSubDSEPSGCode = nEPSGCode;
3981 48 : aosNonJP2Files.push_back(aosGranuleList[i]);
3982 :
3983 48 : if (dfULX < dfMinX)
3984 35 : dfMinX = dfULX;
3985 48 : if (dfULY > dfMaxY)
3986 35 : dfMaxY = dfULY;
3987 48 : const double dfLRX = dfULX + nResolution * nWidth;
3988 48 : const double dfLRY = dfULY - nResolution * nHeight;
3989 48 : if (dfLRX > dfMaxX)
3990 42 : dfMaxX = dfLRX;
3991 48 : if (dfLRY < dfMinY)
3992 42 : dfMinY = dfLRY;
3993 :
3994 96 : SENTINEL2GranuleInfo oGranuleInfo;
3995 48 : oGranuleInfo.osPath = CPLGetPathSafe(aosGranuleList[i]);
3996 48 : if (bIsSafeCompact)
3997 : {
3998 : oGranuleInfo.osBandPrefixPath =
3999 22 : aoL1CSafeCompactGranuleList[i].osBandPrefixPath;
4000 : }
4001 48 : oGranuleInfo.dfMinX = dfULX;
4002 48 : oGranuleInfo.dfMinY = dfLRY;
4003 48 : oGranuleInfo.dfMaxX = dfLRX;
4004 48 : oGranuleInfo.dfMaxY = dfULY;
4005 48 : oGranuleInfo.nWidth = nWidth / (nSubDSPrecision / nResolution);
4006 48 : oGranuleInfo.nHeight = nHeight / (nSubDSPrecision / nResolution);
4007 48 : aosGranuleInfoList.push_back(std::move(oGranuleInfo));
4008 : }
4009 : }
4010 48 : if (dfMinX > dfMaxX)
4011 : {
4012 13 : CPLError(CE_Failure, CPLE_NotSupported,
4013 : "No granule found for EPSG code %d", nSubDSEPSGCode);
4014 13 : return nullptr;
4015 : }
4016 :
4017 35 : const int nRasterXSize =
4018 35 : static_cast<int>((dfMaxX - dfMinX) / nSubDSPrecision + 0.5);
4019 35 : const int nRasterYSize =
4020 35 : static_cast<int>((dfMaxY - dfMinY) / nSubDSPrecision + 0.5);
4021 35 : SENTINEL2Dataset *poDS = new SENTINEL2Dataset(nRasterXSize, nRasterYSize);
4022 :
4023 : #if defined(__GNUC__)
4024 : #pragma GCC diagnostic push
4025 : #pragma GCC diagnostic ignored "-Wnull-dereference"
4026 : #endif
4027 35 : poDS->aosNonJP2Files = aosNonJP2Files;
4028 : #if defined(__GNUC__)
4029 : #pragma GCC diagnostic pop
4030 : #endif
4031 :
4032 35 : OGRSpatialReference oSRS;
4033 35 : char *pszProjection = nullptr;
4034 70 : if (oSRS.importFromEPSG(nSubDSEPSGCode) == OGRERR_NONE &&
4035 35 : oSRS.exportToWkt(&pszProjection) == OGRERR_NONE)
4036 : {
4037 35 : poDS->SetProjection(pszProjection);
4038 35 : CPLFree(pszProjection);
4039 : }
4040 : else
4041 : {
4042 0 : CPLDebug("SENTINEL2", "Invalid EPSG code %d", nSubDSEPSGCode);
4043 : }
4044 :
4045 35 : poDS->SetGeoTransform(GDALGeoTransform(dfMinX, nSubDSPrecision, 0, dfMaxY,
4046 35 : 0, -nSubDSPrecision));
4047 35 : poDS->GDALDataset::SetMetadataItem("COMPRESSION", "JPEG2000",
4048 : "IMAGE_STRUCTURE");
4049 35 : if (bIsPreview || bIsTCI)
4050 7 : poDS->GDALDataset::SetMetadataItem("INTERLEAVE", "PIXEL",
4051 : "IMAGE_STRUCTURE");
4052 :
4053 35 : int nBits = (bIsPreview || bIsTCI) ? 8 : 0 /* 0 = unknown yet*/;
4054 35 : int nValMax = (bIsPreview || bIsTCI) ? 255 : 0 /* 0 = unknown yet*/;
4055 : const int nBands =
4056 30 : (bIsPreview || bIsTCI)
4057 37 : ? 3
4058 28 : : ((bAlpha) ? 1 : 0) + static_cast<int>(aosBands.size());
4059 35 : const int nAlphaBand = (bIsPreview || bIsTCI || !bAlpha) ? 0 : nBands;
4060 35 : const GDALDataType eDT = (bIsPreview || bIsTCI) ? GDT_Byte : GDT_UInt16;
4061 :
4062 227 : for (int nBand = 1; nBand <= nBands; nBand++)
4063 : {
4064 192 : VRTSourcedRasterBand *poBand = nullptr;
4065 :
4066 192 : if (nBand != nAlphaBand)
4067 : {
4068 190 : poBand = new VRTSourcedRasterBand(
4069 190 : poDS, nBand, eDT, poDS->nRasterXSize, poDS->nRasterYSize);
4070 : }
4071 : else
4072 : {
4073 2 : poBand = new SENTINEL2AlphaBand(
4074 : poDS, nBand, eDT, poDS->nRasterXSize, poDS->nRasterYSize,
4075 2 : nSaturatedVal, nNodataVal);
4076 : }
4077 :
4078 192 : poDS->SetBand(nBand, poBand);
4079 192 : if (nBand == nAlphaBand)
4080 2 : poBand->SetColorInterpretation(GCI_AlphaBand);
4081 :
4082 384 : CPLString osBandName;
4083 192 : if (nBand != nAlphaBand)
4084 : {
4085 190 : osBandName = aosBands[nBand - 1];
4086 190 : SENTINEL2SetBandMetadata(poBand, osBandName);
4087 : }
4088 : else
4089 2 : osBandName = aosBands[0];
4090 :
4091 459 : for (size_t iSrc = 0; iSrc < aosGranuleInfoList.size(); iSrc++)
4092 : {
4093 267 : const SENTINEL2GranuleInfo &oGranuleInfo = aosGranuleInfoList[iSrc];
4094 267 : CPLString osTile;
4095 :
4096 267 : if (bIsSafeCompact && eLevel != SENTINEL2_L2A)
4097 : {
4098 33 : if (bIsTCI)
4099 : {
4100 6 : osTile = oGranuleInfo.osBandPrefixPath + "TCI.jp2";
4101 : }
4102 : else
4103 : {
4104 27 : osTile = oGranuleInfo.osBandPrefixPath + "B";
4105 27 : if (osBandName.size() == 1)
4106 0 : osTile += "0" + osBandName;
4107 27 : else if (osBandName.size() == 3)
4108 2 : osTile += osBandName.substr(1);
4109 : else
4110 25 : osTile += osBandName;
4111 27 : osTile += ".jp2";
4112 : }
4113 : }
4114 : else
4115 : {
4116 468 : osTile = SENTINEL2GetTilename(
4117 : oGranuleInfo.osPath, CPLGetFilename(oGranuleInfo.osPath),
4118 : osBandName, osProductURI, bIsPreview,
4119 234 : (eLevel == SENTINEL2_L1C) ? 0 : nSubDSPrecision);
4120 113 : if (bIsSafeCompact && eLevel == SENTINEL2_L2A &&
4121 419 : pType == MSI2Ap && osTile.size() >= 34 &&
4122 306 : osTile.substr(osTile.size() - 18, 3) != "MSK")
4123 : {
4124 60 : osTile.insert(osTile.size() - 34, "L2A_");
4125 : }
4126 234 : if (bIsTCI && osTile.size() >= 14)
4127 : {
4128 0 : osTile.replace(osTile.size() - 11, 3, "TCI");
4129 : }
4130 : }
4131 :
4132 267 : bool bTileFound = false;
4133 267 : if (nValMax == 0)
4134 : {
4135 : /* It is supposed to be 12 bits, but some products have 15 bits
4136 : */
4137 28 : if (SENTINEL2GetTileInfo(osTile, nullptr, nullptr, &nBits))
4138 : {
4139 27 : bTileFound = true;
4140 27 : if (nBits <= 16)
4141 27 : nValMax = (1 << nBits) - 1;
4142 : else
4143 : {
4144 0 : CPLDebug("SENTINEL2", "Unexpected bit depth %d", nBits);
4145 0 : nValMax = 65535;
4146 : }
4147 : }
4148 : }
4149 : else
4150 : {
4151 : VSIStatBufL sStat;
4152 239 : if (VSIStatExL(osTile, &sStat, VSI_STAT_EXISTS_FLAG) == 0)
4153 239 : bTileFound = true;
4154 : }
4155 267 : if (!bTileFound)
4156 : {
4157 1 : CPLError(CE_Warning, CPLE_AppDefined,
4158 : "Tile %s not found on filesystem. Skipping it",
4159 : osTile.c_str());
4160 1 : continue;
4161 : }
4162 :
4163 266 : const int nDstXOff = static_cast<int>(
4164 266 : (oGranuleInfo.dfMinX - dfMinX) / nSubDSPrecision + 0.5);
4165 266 : const int nDstYOff = static_cast<int>(
4166 266 : (dfMaxY - oGranuleInfo.dfMaxY) / nSubDSPrecision + 0.5);
4167 :
4168 266 : if (nBand != nAlphaBand)
4169 : {
4170 263 : poBand->AddSimpleSource(
4171 242 : osTile, (bIsPreview || bIsTCI) ? nBand : 1, 0, 0,
4172 263 : oGranuleInfo.nWidth, oGranuleInfo.nHeight, nDstXOff,
4173 263 : nDstYOff, oGranuleInfo.nWidth, oGranuleInfo.nHeight);
4174 : }
4175 : else
4176 : {
4177 3 : poBand->AddComplexSource(
4178 3 : osTile, 1, 0, 0, oGranuleInfo.nWidth, oGranuleInfo.nHeight,
4179 3 : nDstXOff, nDstYOff, oGranuleInfo.nWidth,
4180 3 : oGranuleInfo.nHeight, nValMax /* offset */, 0 /* scale */);
4181 : }
4182 : }
4183 :
4184 192 : if ((nBits % 8) != 0)
4185 : {
4186 143 : poBand->SetMetadataItem("NBITS", CPLSPrintf("%d", nBits),
4187 143 : "IMAGE_STRUCTURE");
4188 : }
4189 : }
4190 :
4191 35 : return poDS;
4192 : }
4193 :
4194 : /************************************************************************/
4195 : /* OpenL1CTileSubdataset() */
4196 : /************************************************************************/
4197 :
4198 13 : GDALDataset *SENTINEL2Dataset::OpenL1CTileSubdataset(GDALOpenInfo *poOpenInfo)
4199 : {
4200 26 : CPLString osFilename;
4201 13 : CPLAssert(STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1C_TILE:"));
4202 13 : osFilename = poOpenInfo->pszFilename + strlen("SENTINEL2_L1C_TILE:");
4203 13 : const char *pszPrecision = strrchr(osFilename.c_str(), ':');
4204 13 : if (pszPrecision == nullptr || pszPrecision == osFilename.c_str())
4205 : {
4206 2 : CPLError(CE_Failure, CPLE_AppDefined,
4207 : "Invalid syntax for SENTINEL2_L1C_TILE:");
4208 2 : return nullptr;
4209 : }
4210 11 : const bool bIsPreview = STARTS_WITH_CI(pszPrecision + 1, "PREVIEW");
4211 11 : const int nSubDSPrecision = (bIsPreview) ? 320 : atoi(pszPrecision + 1);
4212 11 : if (!bIsPreview && nSubDSPrecision != RES_10M &&
4213 2 : nSubDSPrecision != RES_20M && nSubDSPrecision != RES_60M)
4214 : {
4215 1 : CPLError(CE_Failure, CPLE_NotSupported, "Unsupported precision: %d",
4216 : nSubDSPrecision);
4217 1 : return nullptr;
4218 : }
4219 10 : osFilename.resize(pszPrecision - osFilename.c_str());
4220 :
4221 20 : std::set<CPLString> oSetBands;
4222 10 : CPLXMLNode *psRootMainMTD = nullptr;
4223 : GDALDataset *poTmpDS =
4224 10 : OpenL1CTile(osFilename, &psRootMainMTD, nSubDSPrecision, &oSetBands);
4225 20 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRootMainMTD);
4226 10 : if (poTmpDS == nullptr)
4227 2 : return nullptr;
4228 :
4229 16 : std::vector<CPLString> aosBands;
4230 8 : if (bIsPreview)
4231 : {
4232 2 : aosBands.push_back("04");
4233 2 : aosBands.push_back("03");
4234 2 : aosBands.push_back("02");
4235 : }
4236 : else
4237 : {
4238 21 : for (std::set<CPLString>::const_iterator oIterBandnames =
4239 6 : oSetBands.begin();
4240 48 : oIterBandnames != oSetBands.end(); ++oIterBandnames)
4241 : {
4242 21 : aosBands.push_back(*oIterBandnames);
4243 : }
4244 : /* Put 2=Blue, 3=Green, 4=Band bands in RGB order for conveniency */
4245 14 : if (aosBands.size() >= 3 && aosBands[0] == "02" &&
4246 14 : aosBands[1] == "03" && aosBands[2] == "04")
4247 : {
4248 3 : aosBands[0] = "04";
4249 3 : aosBands[2] = "02";
4250 : }
4251 : }
4252 :
4253 : /* -------------------------------------------------------------------- */
4254 : /* Create dataset. */
4255 : /* -------------------------------------------------------------------- */
4256 :
4257 16 : std::vector<CPLString> aosGranuleList;
4258 8 : aosGranuleList.push_back(osFilename);
4259 :
4260 16 : const int nSaturatedVal = atoi(CSLFetchNameValueDef(
4261 8 : poTmpDS->GetMetadata(), "SPECIAL_VALUE_SATURATED", "-1"));
4262 16 : const int nNodataVal = atoi(CSLFetchNameValueDef(
4263 8 : poTmpDS->GetMetadata(), "SPECIAL_VALUE_NODATA", "-1"));
4264 :
4265 : const bool bAlpha =
4266 8 : CPLTestBool(SENTINEL2GetOption(poOpenInfo, "ALPHA", "FALSE"));
4267 :
4268 16 : std::vector<CPLString> aosNonJP2Files;
4269 8 : SENTINEL2Dataset *poDS = CreateL1CL2ADataset(
4270 : SENTINEL2_L1C, MSI2A,
4271 : false, // bIsSafeCompact
4272 16 : aosGranuleList, std::vector<L1CSafeCompatGranuleDescription>(),
4273 : aosNonJP2Files, nSubDSPrecision, bIsPreview,
4274 : false, // bIsTCI
4275 : -1 /*nSubDSEPSGCode*/, bAlpha, aosBands, nSaturatedVal, nNodataVal,
4276 16 : CPLString());
4277 8 : if (poDS == nullptr)
4278 : {
4279 1 : delete poTmpDS;
4280 1 : return nullptr;
4281 : }
4282 :
4283 : // Transfer metadata
4284 7 : poDS->GDALDataset::SetMetadata(poTmpDS->GetMetadata());
4285 7 : poDS->GDALDataset::SetMetadata(poTmpDS->GetMetadata("xml:SENTINEL2"),
4286 : "xml:SENTINEL2");
4287 :
4288 7 : delete poTmpDS;
4289 :
4290 : /* -------------------------------------------------------------------- */
4291 : /* Add extra band metadata. */
4292 : /* -------------------------------------------------------------------- */
4293 7 : if (psRootMainMTD != nullptr)
4294 6 : poDS->AddL1CL2ABandMetadata(SENTINEL2_L1C, psRootMainMTD, aosBands);
4295 :
4296 : /* -------------------------------------------------------------------- */
4297 : /* Initialize overview information. */
4298 : /* -------------------------------------------------------------------- */
4299 7 : poDS->SetDescription(poOpenInfo->pszFilename);
4300 7 : CPLString osOverviewFile;
4301 7 : if (bIsPreview)
4302 2 : osOverviewFile = CPLSPrintf("%s_PREVIEW.tif.ovr", osFilename.c_str());
4303 : else
4304 : osOverviewFile =
4305 5 : CPLSPrintf("%s_%dm.tif.ovr", osFilename.c_str(), nSubDSPrecision);
4306 7 : poDS->SetMetadataItem("OVERVIEW_FILE", osOverviewFile, "OVERVIEWS");
4307 7 : poDS->oOvManager.Initialize(poDS, ":::VIRTUAL:::");
4308 :
4309 7 : return poDS;
4310 : }
4311 :
4312 : /************************************************************************/
4313 : /* GDALRegister_SENTINEL2() */
4314 : /************************************************************************/
4315 :
4316 2033 : void GDALRegister_SENTINEL2()
4317 : {
4318 2033 : if (GDALGetDriverByName("SENTINEL2") != nullptr)
4319 283 : return;
4320 :
4321 1750 : GDALDriver *poDriver = new GDALDriver();
4322 :
4323 1750 : poDriver->SetDescription("SENTINEL2");
4324 1750 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
4325 1750 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
4326 1750 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Sentinel 2");
4327 1750 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
4328 1750 : "drivers/raster/sentinel2.html");
4329 1750 : poDriver->SetMetadataItem(GDAL_DMD_SUBDATASETS, "YES");
4330 :
4331 1750 : poDriver->SetMetadataItem(
4332 : GDAL_DMD_OPENOPTIONLIST,
4333 : "<OpenOptionList>"
4334 : " <Option name='ALPHA' type='boolean' description='Whether to expose "
4335 : "an alpha band' default='NO'/>"
4336 : " <Option name='L1B_MODE' type='string-select' "
4337 : "default='DEFAULT'>\n"
4338 : " <Value>DEFAULT</Value>\n"
4339 : " <Value>GRANULE</Value>\n"
4340 : " <Value>DATASTRIP</Value>\n"
4341 : " </Option>\n"
4342 1750 : "</OpenOptionList>");
4343 :
4344 1750 : poDriver->pfnOpen = SENTINEL2Dataset::Open;
4345 1750 : poDriver->pfnIdentify = SENTINEL2Dataset::Identify;
4346 :
4347 1750 : GetGDALDriverManager()->RegisterDriver(poDriver);
4348 : }
|
