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