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