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