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