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 52997 : int SENTINEL2Dataset::Identify(GDALOpenInfo *poOpenInfo)
372 : {
373 52997 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1B:"))
374 36 : return TRUE;
375 52961 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1C:"))
376 78 : return TRUE;
377 52883 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1C_TILE:"))
378 26 : return TRUE;
379 52857 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L2A:"))
380 80 : return TRUE;
381 :
382 52777 : 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 52777 : 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 158286 : if ((STARTS_WITH_CI(pszJustFilename, "S2A_MSIL1C_") ||
393 52775 : STARTS_WITH_CI(pszJustFilename, "S2B_MSIL1C_") ||
394 52774 : STARTS_WITH_CI(pszJustFilename, "S2A_MSIL2A_") ||
395 52773 : STARTS_WITH_CI(pszJustFilename, "S2B_MSIL2A_") ||
396 52776 : STARTS_WITH_CI(pszJustFilename, "S2A_OPER_PRD_MSI") ||
397 52776 : STARTS_WITH_CI(pszJustFilename, "S2B_OPER_PRD_MSI") ||
398 52774 : STARTS_WITH_CI(pszJustFilename, "S2A_USER_PRD_MSI") ||
399 105552 : STARTS_WITH_CI(pszJustFilename, "S2B_USER_PRD_MSI")) &&
400 52777 : EQUAL(CPLGetExtensionSafe(pszJustFilename).c_str(), "zip"))
401 : {
402 0 : return TRUE;
403 : }
404 :
405 52734 : if (poOpenInfo->nHeaderBytes < 100)
406 48904 : return FALSE;
407 :
408 3830 : const char *pszHeader =
409 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader);
410 :
411 3830 : 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 3817 : 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 3809 : 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 3791 : 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 3783 : if (strstr(pszHeader, "<n1:Level-2A_User_Product") != nullptr &&
428 12 : strstr(pszHeader, "User_Product_Level-2A") != nullptr)
429 12 : return TRUE;
430 :
431 3771 : if (SENTINEL2isZipped(pszHeader, poOpenInfo->nHeaderBytes))
432 10 : return TRUE;
433 :
434 3792 : 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 37 : EQUAL(CPLGetExtensionSafe(pszJustFilename).c_str(), "zip"))
505 : {
506 0 : const CPLString osBasename(CPLGetBasenameSafe(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, CPLGetBasenameSafe().c_str() 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 37 : EQUAL(CPLGetExtensionSafe(pszJustFilename).c_str(), "zip"))
531 : {
532 0 : const CPLString osBasename(CPLGetBasenameSafe(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(CPLGetExtensionSafe(osSAFE).c_str(), "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 37 : EQUAL(CPLGetExtensionSafe(pszJustFilename).c_str(), "zip"))
550 : {
551 0 : const CPLString osBasename(CPLGetBasenameSafe(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(CPLGetExtensionSafe(osSAFE).c_str(), "SAFE"))
558 0 : osSAFE += ".SAFE";
559 0 : osFilename = osFilename + "/" + osSAFE + "/MTD_MSIL2A.xml";
560 0 : if (strncmp(osFilename, "/vsizip/", strlen("/vsizip/")) != 0)
561 0 : osFilename = "/vsizip/" + osFilename;
562 0 : CPLDebug("SENTINEL2", "Trying %s", osFilename.c_str());
563 0 : GDALOpenInfo oOpenInfo(osFilename, GA_ReadOnly);
564 0 : return Open(&oOpenInfo);
565 : }
566 :
567 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 3807 : static bool SENTINEL2isZipped(const char *pszHeader, int nHeaderBytes)
649 : {
650 3807 : 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 3869 : 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 3854 : 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 CPLFormFilenameSafe()...
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(CPLGetDirnameSafe(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 = CPLGetDirnameSafe(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 52 : CPLString osTopDir(CPLFormFilenameSafe(
1730 52 : CPLFormFilenameSafe(CPLGetDirnameSafe(pszFilename).c_str(), "..",
1731 : nullptr)
1732 : .c_str(),
1733 52 : "..", nullptr));
1734 :
1735 : // Workaround to avoid long filenames on Windows
1736 26 : if (CPLIsFilenameRelative(pszFilename))
1737 : {
1738 : // GRANULE/bla/bla.xml
1739 30 : const std::string osPath = CPLGetPathSafe(pszFilename);
1740 16 : if (osPath.find('/') != std::string::npos ||
1741 1 : osPath.find('\\') != std::string::npos)
1742 : {
1743 14 : osTopDir = CPLGetPathSafe(CPLGetPathSafe(osPath.c_str()).c_str());
1744 14 : if (osTopDir == "")
1745 0 : osTopDir = ".";
1746 : }
1747 : }
1748 :
1749 26 : char **papszContents = VSIReadDir(osTopDir);
1750 26 : CPLString osMainMTD;
1751 221 : for (char **papszIter = papszContents; papszIter && *papszIter; ++papszIter)
1752 : {
1753 211 : if (strlen(*papszIter) >= strlen("S2A_XXXX_MTD") &&
1754 28 : (STARTS_WITH_CI(*papszIter, "S2A_") ||
1755 19 : STARTS_WITH_CI(*papszIter, "S2B_")) &&
1756 16 : EQUALN(*papszIter + strlen("S2A_XXXX"), "_MTD", 4))
1757 : {
1758 16 : osMainMTD = CPLFormFilenameSafe(osTopDir, *papszIter, nullptr);
1759 16 : break;
1760 : }
1761 : }
1762 26 : CSLDestroy(papszContents);
1763 52 : return osMainMTD;
1764 : }
1765 :
1766 : /************************************************************************/
1767 : /* SENTINEL2GetResolutionSetAndMainMDFromGranule() */
1768 : /************************************************************************/
1769 :
1770 26 : static void SENTINEL2GetResolutionSetAndMainMDFromGranule(
1771 : const char *pszFilename, const char *pszRootPathWithoutEqual,
1772 : int nResolutionOfInterest, std::set<int> &oSetResolutions,
1773 : std::map<int, std::set<CPLString>> &oMapResolutionsToBands, char **&papszMD,
1774 : CPLXMLNode **ppsRootMainMTD)
1775 : {
1776 52 : CPLString osMainMTD(SENTINEL2GetMainMTDFilenameFromGranuleMTD(pszFilename));
1777 :
1778 : // Parse product MTD if available
1779 26 : papszMD = nullptr;
1780 42 : if (!osMainMTD.empty() &&
1781 : /* env var for debug only */
1782 16 : CPLTestBool(CPLGetConfigOption("SENTINEL2_USE_MAIN_MTD", "YES")))
1783 : {
1784 14 : CPLXMLNode *psRootMainMTD = CPLParseXMLFile(osMainMTD);
1785 14 : if (psRootMainMTD != nullptr)
1786 : {
1787 14 : CPLStripXMLNamespace(psRootMainMTD, nullptr, TRUE);
1788 :
1789 14 : CPLXMLNode *psProductInfo = CPLGetXMLNode(
1790 : psRootMainMTD, CPLSPrintf("=%s.General_Info.Product_Info",
1791 : pszRootPathWithoutEqual));
1792 14 : if (psProductInfo != nullptr)
1793 : {
1794 14 : SENTINEL2GetResolutionSet(psProductInfo, oSetResolutions,
1795 : oMapResolutionsToBands);
1796 : }
1797 :
1798 14 : papszMD = SENTINEL2GetUserProductMetadata(psRootMainMTD,
1799 : pszRootPathWithoutEqual);
1800 14 : if (ppsRootMainMTD != nullptr)
1801 6 : *ppsRootMainMTD = psRootMainMTD;
1802 : else
1803 8 : CPLDestroyXMLNode(psRootMainMTD);
1804 : }
1805 : }
1806 : else
1807 : {
1808 : // If no main MTD file found, then probe all bands for resolution (of
1809 : // interest if there's one, or all resolutions otherwise)
1810 168 : for (size_t i = 0; i < NB_BANDS; i++)
1811 : {
1812 156 : if (nResolutionOfInterest != 0 &&
1813 117 : asBandDesc[i].nResolution != nResolutionOfInterest)
1814 : {
1815 83 : continue;
1816 : }
1817 : CPLString osBandName =
1818 146 : asBandDesc[i].pszBandName + 1; /* skip B character */
1819 73 : if (atoi(osBandName) < 10)
1820 62 : osBandName = "0" + osBandName;
1821 :
1822 : CPLString osTile(SENTINEL2GetTilename(
1823 146 : CPLGetPathSafe(pszFilename).c_str(),
1824 365 : CPLGetBasenameSafe(pszFilename).c_str(), osBandName));
1825 : VSIStatBufL sStat;
1826 73 : if (VSIStatExL(osTile, &sStat, VSI_STAT_EXISTS_FLAG) == 0)
1827 : {
1828 20 : oMapResolutionsToBands[asBandDesc[i].nResolution].insert(
1829 20 : osBandName);
1830 20 : oSetResolutions.insert(asBandDesc[i].nResolution);
1831 : }
1832 : }
1833 : }
1834 26 : }
1835 :
1836 : /************************************************************************/
1837 : /* OpenL1BGranule() */
1838 : /************************************************************************/
1839 :
1840 18 : GDALDataset *SENTINEL2Dataset::OpenL1BGranule(const char *pszFilename,
1841 : CPLXMLNode **ppsRoot,
1842 : int nResolutionOfInterest,
1843 : std::set<CPLString> *poBandSet)
1844 : {
1845 18 : CPLXMLNode *psRoot = CPLParseXMLFile(pszFilename);
1846 18 : if (psRoot == nullptr)
1847 : {
1848 3 : CPLDebug("SENTINEL2", "Cannot XML parse %s", pszFilename);
1849 3 : return nullptr;
1850 : }
1851 :
1852 15 : char *pszOriginalXML = CPLSerializeXMLTree(psRoot);
1853 30 : CPLString osOriginalXML;
1854 15 : if (pszOriginalXML)
1855 15 : osOriginalXML = pszOriginalXML;
1856 15 : CPLFree(pszOriginalXML);
1857 :
1858 30 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRoot);
1859 15 : CPLStripXMLNamespace(psRoot, nullptr, TRUE);
1860 :
1861 15 : SENTINEL2DatasetContainer *poDS = new SENTINEL2DatasetContainer();
1862 :
1863 15 : if (!osOriginalXML.empty())
1864 : {
1865 : char *apszXMLMD[2];
1866 15 : apszXMLMD[0] = const_cast<char *>(osOriginalXML.c_str());
1867 15 : apszXMLMD[1] = nullptr;
1868 15 : poDS->GDALDataset::SetMetadata(apszXMLMD, "xml:SENTINEL2");
1869 : }
1870 :
1871 30 : std::set<int> oSetResolutions;
1872 15 : std::map<int, std::set<CPLString>> oMapResolutionsToBands;
1873 15 : char **papszMD = nullptr;
1874 15 : SENTINEL2GetResolutionSetAndMainMDFromGranule(
1875 : pszFilename, "Level-1B_User_Product", nResolutionOfInterest,
1876 : oSetResolutions, oMapResolutionsToBands, papszMD, nullptr);
1877 15 : if (poBandSet != nullptr)
1878 12 : *poBandSet = oMapResolutionsToBands[nResolutionOfInterest];
1879 :
1880 15 : char **papszGranuleMD = SENTINEL2GetL1BGranuleMetadata(psRoot);
1881 15 : papszMD = CSLMerge(papszMD, papszGranuleMD);
1882 15 : CSLDestroy(papszGranuleMD);
1883 :
1884 : // Remove CLOUD_COVERAGE_ASSESSMENT that comes from main metadata, if
1885 : // granule CLOUDY_PIXEL_PERCENTAGE is present.
1886 21 : if (CSLFetchNameValue(papszMD, "CLOUDY_PIXEL_PERCENTAGE") != nullptr &&
1887 6 : CSLFetchNameValue(papszMD, "CLOUD_COVERAGE_ASSESSMENT") != nullptr)
1888 : {
1889 6 : papszMD =
1890 6 : CSLSetNameValue(papszMD, "CLOUD_COVERAGE_ASSESSMENT", nullptr);
1891 : }
1892 :
1893 15 : poDS->GDALDataset::SetMetadata(papszMD);
1894 15 : CSLDestroy(papszMD);
1895 :
1896 : // Get the footprint
1897 : const char *pszPosList =
1898 15 : CPLGetXMLValue(psRoot,
1899 : "=Level-1B_Granule_ID.Geometric_Info.Granule_Footprint."
1900 : "Granule_Footprint.Footprint.EXT_POS_LIST",
1901 : nullptr);
1902 15 : if (pszPosList != nullptr)
1903 : {
1904 12 : CPLString osPolygon = SENTINEL2GetPolygonWKTFromPosList(pszPosList);
1905 6 : if (!osPolygon.empty())
1906 6 : poDS->GDALDataset::SetMetadataItem("FOOTPRINT", osPolygon.c_str());
1907 : }
1908 :
1909 : /* Create subdatsets per resolution (10, 20, 60m) */
1910 15 : int iSubDSNum = 1;
1911 37 : for (std::set<int>::const_iterator oIterRes = oSetResolutions.begin();
1912 59 : oIterRes != oSetResolutions.end(); ++oIterRes)
1913 : {
1914 22 : const int nResolution = *oIterRes;
1915 :
1916 22 : poDS->GDALDataset::SetMetadataItem(
1917 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
1918 : CPLSPrintf("SENTINEL2_L1B:%s:%dm", pszFilename, nResolution),
1919 : "SUBDATASETS");
1920 :
1921 : CPLString osBandNames = SENTINEL2GetBandListForResolution(
1922 44 : oMapResolutionsToBands[nResolution]);
1923 :
1924 : CPLString osDesc(CPLSPrintf("Bands %s with %dm resolution",
1925 22 : osBandNames.c_str(), nResolution));
1926 22 : poDS->GDALDataset::SetMetadataItem(
1927 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
1928 : "SUBDATASETS");
1929 :
1930 22 : iSubDSNum++;
1931 : }
1932 :
1933 15 : if (ppsRoot != nullptr)
1934 : {
1935 12 : *ppsRoot = oXMLHolder.Release();
1936 : }
1937 :
1938 15 : return poDS;
1939 : }
1940 :
1941 : /************************************************************************/
1942 : /* SENTINEL2SetBandMetadata() */
1943 : /************************************************************************/
1944 :
1945 209 : static void SENTINEL2SetBandMetadata(GDALRasterBand *poBand,
1946 : const CPLString &osBandName)
1947 : {
1948 418 : CPLString osLookupBandName(osBandName);
1949 209 : if (osLookupBandName[0] == '0')
1950 138 : osLookupBandName = osLookupBandName.substr(1);
1951 209 : if (atoi(osLookupBandName) > 0)
1952 168 : osLookupBandName = "B" + osLookupBandName;
1953 :
1954 418 : CPLString osBandDesc(osLookupBandName);
1955 : const SENTINEL2BandDescription *psBandDesc =
1956 209 : SENTINEL2GetBandDesc(osLookupBandName);
1957 209 : if (psBandDesc != nullptr)
1958 : {
1959 : osBandDesc +=
1960 168 : CPLSPrintf(", central wavelength %d nm", psBandDesc->nWaveLength);
1961 168 : poBand->SetColorInterpretation(psBandDesc->eColorInterp);
1962 168 : poBand->SetMetadataItem("BANDNAME", psBandDesc->pszBandName);
1963 168 : poBand->SetMetadataItem("BANDWIDTH",
1964 168 : CPLSPrintf("%d", psBandDesc->nBandWidth));
1965 168 : poBand->SetMetadataItem("BANDWIDTH_UNIT", "nm");
1966 168 : poBand->SetMetadataItem("WAVELENGTH",
1967 168 : CPLSPrintf("%d", psBandDesc->nWaveLength));
1968 168 : poBand->SetMetadataItem("WAVELENGTH_UNIT", "nm");
1969 :
1970 168 : poBand->SetMetadataItem(
1971 : "CENTRAL_WAVELENGTH_UM",
1972 168 : CPLSPrintf("%.3f", double(psBandDesc->nWaveLength) / 1000),
1973 168 : "IMAGERY");
1974 168 : poBand->SetMetadataItem(
1975 : "FWHM_UM",
1976 168 : CPLSPrintf("%.3f", double(psBandDesc->nBandWidth) / 1000),
1977 168 : "IMAGERY");
1978 : }
1979 : else
1980 : {
1981 : const SENTINEL2_L2A_BandDescription *psL2ABandDesc =
1982 41 : SENTINEL2GetL2ABandDesc(osBandName);
1983 41 : if (psL2ABandDesc != nullptr)
1984 : {
1985 41 : osBandDesc += ", ";
1986 41 : osBandDesc += psL2ABandDesc->pszBandDescription;
1987 : }
1988 :
1989 41 : poBand->SetMetadataItem("BANDNAME", osBandName);
1990 : }
1991 209 : poBand->SetDescription(osBandDesc);
1992 209 : }
1993 :
1994 : /************************************************************************/
1995 : /* OpenL1BSubdataset() */
1996 : /************************************************************************/
1997 :
1998 18 : GDALDataset *SENTINEL2Dataset::OpenL1BSubdataset(GDALOpenInfo *poOpenInfo)
1999 : {
2000 36 : CPLString osFilename;
2001 18 : CPLAssert(STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1B:"));
2002 18 : osFilename = poOpenInfo->pszFilename + strlen("SENTINEL2_L1B:");
2003 18 : const char *pszPrecision = strrchr(osFilename.c_str(), ':');
2004 18 : if (pszPrecision == nullptr || pszPrecision == osFilename.c_str())
2005 : {
2006 2 : CPLError(CE_Failure, CPLE_AppDefined,
2007 : "Invalid syntax for SENTINEL2_L1B:");
2008 2 : return nullptr;
2009 : }
2010 16 : const int nSubDSPrecision = atoi(pszPrecision + 1);
2011 16 : if (nSubDSPrecision != 10 && nSubDSPrecision != 20 && nSubDSPrecision != 60)
2012 : {
2013 2 : CPLError(CE_Failure, CPLE_NotSupported, "Unsupported precision: %d",
2014 : nSubDSPrecision);
2015 2 : return nullptr;
2016 : }
2017 14 : osFilename.resize(pszPrecision - osFilename.c_str());
2018 :
2019 14 : CPLXMLNode *psRoot = nullptr;
2020 28 : std::set<CPLString> oSetBands;
2021 : GDALDataset *poTmpDS =
2022 14 : OpenL1BGranule(osFilename, &psRoot, nSubDSPrecision, &oSetBands);
2023 14 : if (poTmpDS == nullptr)
2024 : {
2025 2 : CPLDebug("SENTINEL2", "Failed to open L1B granule %s",
2026 : osFilename.c_str());
2027 2 : return nullptr;
2028 : }
2029 :
2030 24 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRoot);
2031 :
2032 24 : std::vector<CPLString> aosBands;
2033 31 : for (std::set<CPLString>::const_iterator oIterBandnames = oSetBands.begin();
2034 50 : oIterBandnames != oSetBands.end(); ++oIterBandnames)
2035 : {
2036 19 : aosBands.push_back(*oIterBandnames);
2037 : }
2038 : /* Put 2=Blue, 3=Green, 4=Band bands in RGB order for conveniency */
2039 13 : if (aosBands.size() >= 3 && aosBands[0] == "02" && aosBands[1] == "03" &&
2040 1 : aosBands[2] == "04")
2041 : {
2042 1 : aosBands[0] = "04";
2043 1 : aosBands[2] = "02";
2044 : }
2045 :
2046 12 : int nBits = 0; /* 0 = unknown yet*/
2047 12 : int nValMax = 0; /* 0 = unknown yet*/
2048 12 : int nRows = 0;
2049 12 : int nCols = 0;
2050 12 : CPLXMLNode *psGranuleDimensions = CPLGetXMLNode(
2051 : psRoot, "=Level-1B_Granule_ID.Geometric_Info.Granule_Dimensions");
2052 12 : if (psGranuleDimensions == nullptr)
2053 : {
2054 3 : for (size_t i = 0; i < aosBands.size(); i++)
2055 : {
2056 : CPLString osTile(SENTINEL2GetTilename(
2057 2 : CPLGetPathSafe(osFilename).c_str(),
2058 4 : CPLGetBasenameSafe(osFilename).c_str(), aosBands[i]));
2059 1 : if (SENTINEL2GetTileInfo(osTile, &nCols, &nRows, &nBits))
2060 : {
2061 1 : if (nBits <= 16)
2062 1 : nValMax = (1 << nBits) - 1;
2063 : else
2064 : {
2065 0 : CPLDebug("SENTINEL2", "Unexpected bit depth %d", nBits);
2066 0 : nValMax = 65535;
2067 : }
2068 1 : break;
2069 : }
2070 : }
2071 : }
2072 : else
2073 : {
2074 9 : for (CPLXMLNode *psIter = psGranuleDimensions->psChild;
2075 23 : psIter != nullptr; psIter = psIter->psNext)
2076 : {
2077 22 : if (psIter->eType != CXT_Element)
2078 9 : continue;
2079 26 : if (EQUAL(psIter->pszValue, "Size") &&
2080 13 : atoi(CPLGetXMLValue(psIter, "resolution", "")) ==
2081 : nSubDSPrecision)
2082 : {
2083 8 : const char *pszRows = CPLGetXMLValue(psIter, "NROWS", nullptr);
2084 8 : if (pszRows == nullptr)
2085 : {
2086 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
2087 : "NROWS");
2088 1 : delete poTmpDS;
2089 1 : return nullptr;
2090 : }
2091 7 : const char *pszCols = CPLGetXMLValue(psIter, "NCOLS", nullptr);
2092 7 : if (pszCols == nullptr)
2093 : {
2094 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
2095 : "NCOLS");
2096 1 : delete poTmpDS;
2097 1 : return nullptr;
2098 : }
2099 6 : nRows = atoi(pszRows);
2100 6 : nCols = atoi(pszCols);
2101 6 : break;
2102 : }
2103 : }
2104 : }
2105 10 : if (nRows <= 0 || nCols <= 0)
2106 : {
2107 3 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find granule dimension");
2108 3 : delete poTmpDS;
2109 3 : return nullptr;
2110 : }
2111 :
2112 7 : SENTINEL2Dataset *poDS = new SENTINEL2Dataset(nCols, nRows);
2113 7 : poDS->aosNonJP2Files.push_back(osFilename);
2114 :
2115 : // Transfer metadata
2116 7 : poDS->GDALDataset::SetMetadata(poTmpDS->GetMetadata());
2117 7 : poDS->GDALDataset::SetMetadata(poTmpDS->GetMetadata("xml:SENTINEL2"),
2118 : "xml:SENTINEL2");
2119 :
2120 7 : delete poTmpDS;
2121 :
2122 : /* -------------------------------------------------------------------- */
2123 : /* Initialize bands. */
2124 : /* -------------------------------------------------------------------- */
2125 :
2126 21 : int nSaturatedVal = atoi(CSLFetchNameValueDef(
2127 7 : poDS->GetMetadata(), "SPECIAL_VALUE_SATURATED", "-1"));
2128 7 : int nNodataVal = atoi(CSLFetchNameValueDef(poDS->GetMetadata(),
2129 : "SPECIAL_VALUE_NODATA", "-1"));
2130 :
2131 : const bool bAlpha =
2132 7 : CPLTestBool(SENTINEL2GetOption(poOpenInfo, "ALPHA", "FALSE"));
2133 7 : const int nBands = ((bAlpha) ? 1 : 0) + static_cast<int>(aosBands.size());
2134 7 : const int nAlphaBand = (!bAlpha) ? 0 : nBands;
2135 7 : const GDALDataType eDT = GDT_UInt16;
2136 :
2137 27 : for (int nBand = 1; nBand <= nBands; nBand++)
2138 : {
2139 20 : VRTSourcedRasterBand *poBand = nullptr;
2140 :
2141 20 : if (nBand != nAlphaBand)
2142 : {
2143 19 : poBand = new VRTSourcedRasterBand(
2144 19 : poDS, nBand, eDT, poDS->nRasterXSize, poDS->nRasterYSize);
2145 : }
2146 : else
2147 : {
2148 1 : poBand = new SENTINEL2AlphaBand(
2149 : poDS, nBand, eDT, poDS->nRasterXSize, poDS->nRasterYSize,
2150 1 : nSaturatedVal, nNodataVal);
2151 : }
2152 :
2153 20 : poDS->SetBand(nBand, poBand);
2154 20 : if (nBand == nAlphaBand)
2155 1 : poBand->SetColorInterpretation(GCI_AlphaBand);
2156 :
2157 20 : CPLString osBandName;
2158 20 : if (nBand != nAlphaBand)
2159 : {
2160 19 : osBandName = aosBands[nBand - 1];
2161 19 : SENTINEL2SetBandMetadata(poBand, osBandName);
2162 : }
2163 : else
2164 1 : osBandName = aosBands[0];
2165 :
2166 : CPLString osTile(SENTINEL2GetTilename(
2167 40 : CPLGetPathSafe(osFilename).c_str(),
2168 80 : CPLGetBasenameSafe(osFilename).c_str(), osBandName));
2169 :
2170 20 : bool bTileFound = false;
2171 20 : if (nValMax == 0)
2172 : {
2173 : /* It is supposed to be 12 bits, but some products have 15 bits */
2174 6 : if (SENTINEL2GetTileInfo(osTile, nullptr, nullptr, &nBits))
2175 : {
2176 5 : bTileFound = true;
2177 5 : if (nBits <= 16)
2178 5 : nValMax = (1 << nBits) - 1;
2179 : else
2180 : {
2181 0 : CPLDebug("SENTINEL2", "Unexpected bit depth %d", nBits);
2182 0 : nValMax = 65535;
2183 : }
2184 : }
2185 : }
2186 : else
2187 : {
2188 : VSIStatBufL sStat;
2189 14 : if (VSIStatExL(osTile, &sStat, VSI_STAT_EXISTS_FLAG) == 0)
2190 14 : bTileFound = true;
2191 : }
2192 20 : if (!bTileFound)
2193 : {
2194 1 : CPLError(CE_Warning, CPLE_AppDefined,
2195 : "Tile %s not found on filesystem. Skipping it",
2196 : osTile.c_str());
2197 1 : continue;
2198 : }
2199 :
2200 19 : if (nBand != nAlphaBand)
2201 : {
2202 18 : poBand->AddSimpleSource(osTile, 1, 0, 0, poDS->nRasterXSize,
2203 18 : poDS->nRasterYSize, 0, 0,
2204 18 : poDS->nRasterXSize, poDS->nRasterYSize);
2205 : }
2206 : else
2207 : {
2208 1 : poBand->AddComplexSource(osTile, 1, 0, 0, poDS->nRasterXSize,
2209 1 : poDS->nRasterYSize, 0, 0,
2210 1 : poDS->nRasterXSize, poDS->nRasterYSize,
2211 : nValMax /* offset */, 0 /* scale */);
2212 : }
2213 :
2214 19 : if ((nBits % 8) != 0)
2215 : {
2216 19 : poBand->SetMetadataItem("NBITS", CPLSPrintf("%d", nBits),
2217 19 : "IMAGE_STRUCTURE");
2218 : }
2219 : }
2220 :
2221 : /* -------------------------------------------------------------------- */
2222 : /* Add georeferencing. */
2223 : /* -------------------------------------------------------------------- */
2224 : // const char* pszDirection =
2225 : // poDS->GetMetadataItem("DATATAKE_1_SENSING_ORBIT_DIRECTION");
2226 7 : const char *pszFootprint = poDS->GetMetadataItem("FOOTPRINT");
2227 7 : if (pszFootprint != nullptr)
2228 : {
2229 : // For descending orbits, we have observed that the order of points in
2230 : // the polygon is UL, LL, LR, UR. That might not be true for ascending
2231 : // orbits but let's assume it...
2232 4 : OGRGeometry *poGeom = nullptr;
2233 4 : if (OGRGeometryFactory::createFromWkt(pszFootprint, nullptr, &poGeom) ==
2234 4 : OGRERR_NONE &&
2235 8 : poGeom != nullptr &&
2236 4 : wkbFlatten(poGeom->getGeometryType()) == wkbPolygon)
2237 : {
2238 : OGRLinearRing *poRing =
2239 4 : reinterpret_cast<OGRPolygon *>(poGeom)->getExteriorRing();
2240 4 : if (poRing != nullptr && poRing->getNumPoints() == 5)
2241 : {
2242 : GDAL_GCP asGCPList[5];
2243 4 : memset(asGCPList, 0, sizeof(asGCPList));
2244 20 : for (int i = 0; i < 4; i++)
2245 : {
2246 16 : asGCPList[i].dfGCPX = poRing->getX(i);
2247 16 : asGCPList[i].dfGCPY = poRing->getY(i);
2248 16 : asGCPList[i].dfGCPZ = poRing->getZ(i);
2249 : }
2250 4 : asGCPList[0].dfGCPPixel = 0;
2251 4 : asGCPList[0].dfGCPLine = 0;
2252 4 : asGCPList[1].dfGCPPixel = 0;
2253 4 : asGCPList[1].dfGCPLine = poDS->nRasterYSize;
2254 4 : asGCPList[2].dfGCPPixel = poDS->nRasterXSize;
2255 4 : asGCPList[2].dfGCPLine = poDS->nRasterYSize;
2256 4 : asGCPList[3].dfGCPPixel = poDS->nRasterXSize;
2257 4 : asGCPList[3].dfGCPLine = 0;
2258 :
2259 4 : int nGCPCount = 4;
2260 : CPLXMLNode *psGeometryHeader =
2261 4 : CPLGetXMLNode(psRoot, "=Level-1B_Granule_ID.Geometric_Info."
2262 : "Granule_Position.Geometric_Header");
2263 4 : if (psGeometryHeader != nullptr)
2264 : {
2265 4 : const char *pszGC = CPLGetXMLValue(
2266 : psGeometryHeader, "GROUND_CENTER", nullptr);
2267 : const char *pszQLCenter =
2268 4 : CPLGetXMLValue(psGeometryHeader, "QL_CENTER", nullptr);
2269 4 : if (pszGC != nullptr && pszQLCenter != nullptr &&
2270 4 : EQUAL(pszQLCenter, "0 0"))
2271 : {
2272 4 : char **papszTokens = CSLTokenizeString(pszGC);
2273 4 : if (CSLCount(papszTokens) >= 2)
2274 : {
2275 4 : nGCPCount = 5;
2276 4 : asGCPList[4].dfGCPX = CPLAtof(papszTokens[1]);
2277 4 : asGCPList[4].dfGCPY = CPLAtof(papszTokens[0]);
2278 4 : if (CSLCount(papszTokens) >= 3)
2279 4 : asGCPList[4].dfGCPZ = CPLAtof(papszTokens[2]);
2280 4 : asGCPList[4].dfGCPLine = poDS->nRasterYSize / 2.0;
2281 4 : asGCPList[4].dfGCPPixel = poDS->nRasterXSize / 2.0;
2282 : }
2283 4 : CSLDestroy(papszTokens);
2284 : }
2285 : }
2286 :
2287 4 : poDS->SetGCPs(nGCPCount, asGCPList, SRS_WKT_WGS84_LAT_LONG);
2288 4 : GDALDeinitGCPs(nGCPCount, asGCPList);
2289 : }
2290 : }
2291 4 : delete poGeom;
2292 : }
2293 :
2294 : /* -------------------------------------------------------------------- */
2295 : /* Initialize overview information. */
2296 : /* -------------------------------------------------------------------- */
2297 7 : poDS->SetDescription(poOpenInfo->pszFilename);
2298 7 : CPLString osOverviewFile;
2299 : osOverviewFile =
2300 7 : CPLSPrintf("%s_%dm.tif.ovr", osFilename.c_str(), nSubDSPrecision);
2301 7 : poDS->SetMetadataItem("OVERVIEW_FILE", osOverviewFile, "OVERVIEWS");
2302 7 : poDS->oOvManager.Initialize(poDS, ":::VIRTUAL:::");
2303 :
2304 7 : return poDS;
2305 : }
2306 :
2307 : /************************************************************************/
2308 : /* SENTINEL2GetGranuleList_L1CSafeCompact() */
2309 : /************************************************************************/
2310 :
2311 12 : static bool SENTINEL2GetGranuleList_L1CSafeCompact(
2312 : CPLXMLNode *psMainMTD, const char *pszFilename,
2313 : std::vector<L1CSafeCompatGranuleDescription> &osList)
2314 : {
2315 12 : CPLXMLNode *psProductInfo = CPLGetXMLNode(
2316 : psMainMTD, "=Level-1C_User_Product.General_Info.Product_Info");
2317 12 : if (psProductInfo == nullptr)
2318 : {
2319 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
2320 : "=Level-1C_User_Product.General_Info.Product_Info");
2321 0 : return false;
2322 : }
2323 :
2324 : CPLXMLNode *psProductOrganisation =
2325 12 : CPLGetXMLNode(psProductInfo, "Product_Organisation");
2326 12 : if (psProductOrganisation == nullptr)
2327 : {
2328 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
2329 : "Product_Organisation");
2330 0 : return false;
2331 : }
2332 :
2333 12 : CPLString osDirname(CPLGetDirnameSafe(pszFilename));
2334 : #ifdef HAVE_READLINK
2335 : char szPointerFilename[2048];
2336 : int nBytes = static_cast<int>(
2337 12 : readlink(pszFilename, szPointerFilename, sizeof(szPointerFilename)));
2338 12 : if (nBytes != -1)
2339 : {
2340 : const int nOffset =
2341 0 : std::min(nBytes, static_cast<int>(sizeof(szPointerFilename) - 1));
2342 0 : szPointerFilename[nOffset] = '\0';
2343 0 : osDirname = CPLGetDirnameSafe(szPointerFilename);
2344 : }
2345 : #endif
2346 :
2347 12 : const char chSeparator = SENTINEL2GetPathSeparator(osDirname);
2348 24 : for (CPLXMLNode *psIter = psProductOrganisation->psChild; psIter != nullptr;
2349 12 : psIter = psIter->psNext)
2350 : {
2351 12 : if (psIter->eType != CXT_Element ||
2352 12 : !EQUAL(psIter->pszValue, "Granule_List"))
2353 : {
2354 0 : continue;
2355 : }
2356 24 : for (CPLXMLNode *psIter2 = psIter->psChild; psIter2 != nullptr;
2357 12 : psIter2 = psIter2->psNext)
2358 : {
2359 12 : if (psIter2->eType != CXT_Element ||
2360 12 : !EQUAL(psIter2->pszValue, "Granule"))
2361 : {
2362 0 : continue;
2363 : }
2364 :
2365 : const char *pszImageFile =
2366 12 : CPLGetXMLValue(psIter2, "IMAGE_FILE", nullptr);
2367 12 : if (pszImageFile == nullptr || strlen(pszImageFile) < 3)
2368 : {
2369 0 : CPLDebug("SENTINEL2", "Missing IMAGE_FILE element");
2370 0 : continue;
2371 : }
2372 24 : L1CSafeCompatGranuleDescription oDesc;
2373 12 : oDesc.osBandPrefixPath = osDirname + chSeparator + pszImageFile;
2374 12 : oDesc.osBandPrefixPath.resize(oDesc.osBandPrefixPath.size() -
2375 : 3); // strip B12
2376 : // GRANULE/L1C_T30TXT_A007999_20170102T111441/IMG_DATA/T30TXT_20170102T111442_B12
2377 : // --> GRANULE/L1C_T30TXT_A007999_20170102T111441/MTD_TL.xml
2378 : oDesc.osMTDTLPath =
2379 24 : osDirname + chSeparator +
2380 48 : CPLGetDirnameSafe(CPLGetDirnameSafe(pszImageFile).c_str()) +
2381 12 : chSeparator + "MTD_TL.xml";
2382 12 : osList.push_back(oDesc);
2383 : }
2384 : }
2385 :
2386 12 : return true;
2387 : }
2388 :
2389 : /************************************************************************/
2390 : /* SENTINEL2GetGranuleList_L2ASafeCompact() */
2391 : /************************************************************************/
2392 :
2393 15 : static bool SENTINEL2GetGranuleList_L2ASafeCompact(
2394 : CPLXMLNode *psMainMTD, const char *pszFilename,
2395 : std::vector<L1CSafeCompatGranuleDescription> &osList)
2396 : {
2397 15 : const char *pszNodePath =
2398 : "=Level-2A_User_Product.General_Info.Product_Info";
2399 15 : CPLXMLNode *psProductInfo = CPLGetXMLNode(psMainMTD, pszNodePath);
2400 15 : if (psProductInfo == nullptr)
2401 : {
2402 6 : pszNodePath = "=Level-2A_User_Product.General_Info.L2A_Product_Info";
2403 6 : psProductInfo = CPLGetXMLNode(psMainMTD, pszNodePath);
2404 : }
2405 15 : if (psProductInfo == nullptr)
2406 : {
2407 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s", pszNodePath);
2408 0 : return false;
2409 : }
2410 :
2411 : CPLXMLNode *psProductOrganisation =
2412 15 : CPLGetXMLNode(psProductInfo, "Product_Organisation");
2413 15 : if (psProductOrganisation == nullptr)
2414 : {
2415 : psProductOrganisation =
2416 6 : CPLGetXMLNode(psProductInfo, "L2A_Product_Organisation");
2417 6 : if (psProductOrganisation == nullptr)
2418 : {
2419 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
2420 : "Product_Organisation");
2421 0 : return false;
2422 : }
2423 : }
2424 :
2425 15 : CPLString osDirname(CPLGetDirnameSafe(pszFilename));
2426 : #ifdef HAVE_READLINK
2427 : char szPointerFilename[2048];
2428 : int nBytes = static_cast<int>(
2429 15 : readlink(pszFilename, szPointerFilename, sizeof(szPointerFilename)));
2430 15 : if (nBytes != -1)
2431 : {
2432 : const int nOffset =
2433 0 : std::min(nBytes, static_cast<int>(sizeof(szPointerFilename) - 1));
2434 0 : szPointerFilename[nOffset] = '\0';
2435 0 : osDirname = CPLGetDirnameSafe(szPointerFilename);
2436 : }
2437 : #endif
2438 :
2439 15 : const char chSeparator = SENTINEL2GetPathSeparator(osDirname);
2440 42 : for (CPLXMLNode *psIter = psProductOrganisation->psChild; psIter != nullptr;
2441 27 : psIter = psIter->psNext)
2442 : {
2443 27 : if (psIter->eType != CXT_Element ||
2444 27 : !EQUAL(psIter->pszValue, "Granule_List"))
2445 : {
2446 0 : continue;
2447 : }
2448 54 : for (CPLXMLNode *psIter2 = psIter->psChild; psIter2 != nullptr;
2449 27 : psIter2 = psIter2->psNext)
2450 : {
2451 27 : if (psIter2->eType != CXT_Element ||
2452 27 : !EQUAL(psIter2->pszValue, "Granule"))
2453 : {
2454 0 : continue;
2455 : }
2456 :
2457 : const char *pszImageFile =
2458 27 : CPLGetXMLValue(psIter2, "IMAGE_FILE", nullptr);
2459 27 : if (pszImageFile == nullptr)
2460 : {
2461 : pszImageFile =
2462 18 : CPLGetXMLValue(psIter2, "IMAGE_FILE_2A", nullptr);
2463 18 : if (pszImageFile == nullptr || strlen(pszImageFile) < 3)
2464 : {
2465 0 : CPLDebug("SENTINEL2", "Missing IMAGE_FILE element");
2466 0 : continue;
2467 : }
2468 : }
2469 27 : L1CSafeCompatGranuleDescription oDesc;
2470 27 : oDesc.osBandPrefixPath = osDirname + chSeparator + pszImageFile;
2471 27 : if (oDesc.osBandPrefixPath.size() < 36)
2472 : {
2473 0 : CPLDebug("SENTINEL2", "Band prefix path too short");
2474 0 : continue;
2475 : }
2476 27 : oDesc.osBandPrefixPath.resize(oDesc.osBandPrefixPath.size() - 36);
2477 : // GRANULE/L1C_T30TXT_A007999_20170102T111441/IMG_DATA/T30TXT_20170102T111442_B12_60m
2478 : // --> GRANULE/L1C_T30TXT_A007999_20170102T111441/MTD_TL.xml
2479 : oDesc.osMTDTLPath =
2480 54 : osDirname + chSeparator +
2481 81 : CPLGetDirnameSafe(CPLGetDirnameSafe(pszImageFile).c_str());
2482 27 : if (oDesc.osMTDTLPath.size() < 9)
2483 : {
2484 0 : CPLDebug("SENTINEL2", "MTDTL path too short");
2485 0 : continue;
2486 : }
2487 27 : oDesc.osMTDTLPath.resize(oDesc.osMTDTLPath.size() - 9);
2488 27 : oDesc.osMTDTLPath = oDesc.osMTDTLPath + chSeparator + "MTD_TL.xml";
2489 27 : osList.push_back(oDesc);
2490 : }
2491 : }
2492 :
2493 15 : return true;
2494 : }
2495 :
2496 : /************************************************************************/
2497 : /* OpenL1C_L2A() */
2498 : /************************************************************************/
2499 :
2500 17 : GDALDataset *SENTINEL2Dataset::OpenL1C_L2A(const char *pszFilename,
2501 : SENTINEL2Level eLevel)
2502 : {
2503 17 : CPLXMLNode *psRoot = CPLParseXMLFile(pszFilename);
2504 17 : if (psRoot == nullptr)
2505 : {
2506 2 : CPLDebug("SENTINEL2", "Cannot XML parse %s", pszFilename);
2507 2 : return nullptr;
2508 : }
2509 :
2510 15 : char *pszOriginalXML = CPLSerializeXMLTree(psRoot);
2511 30 : CPLString osOriginalXML;
2512 15 : if (pszOriginalXML)
2513 15 : osOriginalXML = pszOriginalXML;
2514 15 : CPLFree(pszOriginalXML);
2515 :
2516 30 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRoot);
2517 15 : CPLStripXMLNamespace(psRoot, nullptr, TRUE);
2518 :
2519 15 : const char *pszNodePath =
2520 : (eLevel == SENTINEL2_L1C)
2521 : ? "=Level-1C_User_Product.General_Info.Product_Info"
2522 : : "=Level-2A_User_Product.General_Info.Product_Info";
2523 15 : CPLXMLNode *psProductInfo = CPLGetXMLNode(psRoot, pszNodePath);
2524 15 : if (psProductInfo == nullptr && eLevel == SENTINEL2_L2A)
2525 : {
2526 4 : pszNodePath = "=Level-2A_User_Product.General_Info.L2A_Product_Info";
2527 4 : psProductInfo = CPLGetXMLNode(psRoot, pszNodePath);
2528 : }
2529 15 : if (psProductInfo == nullptr)
2530 : {
2531 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s", pszNodePath);
2532 1 : return nullptr;
2533 : }
2534 :
2535 : const bool bIsSafeCompact =
2536 14 : EQUAL(CPLGetXMLValue(psProductInfo, "Query_Options.PRODUCT_FORMAT", ""),
2537 : "SAFE_COMPACT");
2538 :
2539 28 : std::set<int> oSetResolutions;
2540 28 : std::map<int, std::set<CPLString>> oMapResolutionsToBands;
2541 14 : if (bIsSafeCompact)
2542 : {
2543 70 : for (unsigned int i = 0; i < NB_BANDS; ++i)
2544 : {
2545 : // L2 does not contain B10
2546 65 : if (i == 10 && eLevel == SENTINEL2_L2A)
2547 3 : continue;
2548 62 : const SENTINEL2BandDescription *psBandDesc = &asBandDesc[i];
2549 62 : oSetResolutions.insert(psBandDesc->nResolution);
2550 : CPLString osName =
2551 124 : psBandDesc->pszBandName + 1; /* skip B character */
2552 62 : if (atoi(osName) < 10)
2553 50 : osName = "0" + osName;
2554 62 : oMapResolutionsToBands[psBandDesc->nResolution].insert(osName);
2555 : }
2556 5 : if (eLevel == SENTINEL2_L2A)
2557 : {
2558 39 : for (const auto &sL2ABandDesc : asL2ABandDesc)
2559 : {
2560 36 : oSetResolutions.insert(sL2ABandDesc.nResolution);
2561 36 : oMapResolutionsToBands[sL2ABandDesc.nResolution].insert(
2562 36 : sL2ABandDesc.pszBandName);
2563 : }
2564 : }
2565 : }
2566 15 : else if (eLevel == SENTINEL2_L1C &&
2567 6 : !SENTINEL2GetResolutionSet(psProductInfo, oSetResolutions,
2568 : oMapResolutionsToBands))
2569 : {
2570 3 : CPLDebug("SENTINEL2", "Failed to get resolution set");
2571 3 : return nullptr;
2572 : }
2573 :
2574 22 : std::vector<CPLString> aosGranuleList;
2575 11 : if (bIsSafeCompact)
2576 : {
2577 : std::vector<L1CSafeCompatGranuleDescription>
2578 5 : aoL1CSafeCompactGranuleList;
2579 7 : if (eLevel == SENTINEL2_L1C &&
2580 2 : !SENTINEL2GetGranuleList_L1CSafeCompact(
2581 : psRoot, pszFilename, aoL1CSafeCompactGranuleList))
2582 : {
2583 0 : CPLDebug("SENTINEL2", "Failed to get granule list");
2584 0 : return nullptr;
2585 : }
2586 8 : else if (eLevel == SENTINEL2_L2A &&
2587 3 : !SENTINEL2GetGranuleList_L2ASafeCompact(
2588 : psRoot, pszFilename, aoL1CSafeCompactGranuleList))
2589 : {
2590 0 : CPLDebug("SENTINEL2", "Failed to get granule list");
2591 0 : return nullptr;
2592 : }
2593 12 : for (size_t i = 0; i < aoL1CSafeCompactGranuleList.size(); ++i)
2594 : {
2595 7 : aosGranuleList.push_back(
2596 7 : aoL1CSafeCompactGranuleList[i].osMTDTLPath);
2597 : }
2598 : }
2599 6 : else if (!SENTINEL2GetGranuleList(
2600 : psRoot, eLevel, pszFilename, aosGranuleList,
2601 : (eLevel == SENTINEL2_L1C) ? nullptr : &oSetResolutions,
2602 : (eLevel == SENTINEL2_L1C) ? nullptr : &oMapResolutionsToBands))
2603 : {
2604 0 : CPLDebug("SENTINEL2", "Failed to get granule list");
2605 0 : return nullptr;
2606 : }
2607 11 : if (oSetResolutions.empty())
2608 : {
2609 0 : CPLDebug("SENTINEL2", "Resolution set is empty");
2610 0 : return nullptr;
2611 : }
2612 :
2613 11 : std::set<int> oSetEPSGCodes;
2614 27 : for (size_t i = 0; i < aosGranuleList.size(); i++)
2615 : {
2616 16 : int nEPSGCode = 0;
2617 16 : if (SENTINEL2GetGranuleInfo(eLevel, aosGranuleList[i],
2618 32 : *(oSetResolutions.begin()), &nEPSGCode))
2619 : {
2620 13 : oSetEPSGCodes.insert(nEPSGCode);
2621 : }
2622 : }
2623 :
2624 11 : SENTINEL2DatasetContainer *poDS = new SENTINEL2DatasetContainer();
2625 11 : char **papszMD = SENTINEL2GetUserProductMetadata(
2626 : psRoot, (eLevel == SENTINEL2_L1C) ? "Level-1C_User_Product"
2627 : : "Level-2A_User_Product");
2628 11 : poDS->GDALDataset::SetMetadata(papszMD);
2629 11 : CSLDestroy(papszMD);
2630 :
2631 11 : if (!osOriginalXML.empty())
2632 : {
2633 : char *apszXMLMD[2];
2634 11 : apszXMLMD[0] = const_cast<char *>(osOriginalXML.c_str());
2635 11 : apszXMLMD[1] = nullptr;
2636 11 : poDS->GDALDataset::SetMetadata(apszXMLMD, "xml:SENTINEL2");
2637 : }
2638 :
2639 11 : const char *pszPrefix =
2640 : (eLevel == SENTINEL2_L1C) ? "SENTINEL2_L1C" : "SENTINEL2_L2A";
2641 :
2642 : /* Create subdatsets per resolution (10, 20, 60m) and EPSG codes */
2643 11 : int iSubDSNum = 1;
2644 38 : for (std::set<int>::const_iterator oIterRes = oSetResolutions.begin();
2645 65 : oIterRes != oSetResolutions.end(); ++oIterRes)
2646 : {
2647 27 : const int nResolution = *oIterRes;
2648 :
2649 50 : for (std::set<int>::const_iterator oIterEPSG = oSetEPSGCodes.begin();
2650 73 : oIterEPSG != oSetEPSGCodes.end(); ++oIterEPSG)
2651 : {
2652 23 : const int nEPSGCode = *oIterEPSG;
2653 23 : poDS->GDALDataset::SetMetadataItem(
2654 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
2655 : CPLSPrintf("%s:%s:%dm:EPSG_%d", pszPrefix, pszFilename,
2656 : nResolution, nEPSGCode),
2657 : "SUBDATASETS");
2658 :
2659 : CPLString osBandNames = SENTINEL2GetBandListForResolution(
2660 46 : oMapResolutionsToBands[nResolution]);
2661 :
2662 : CPLString osDesc(CPLSPrintf("Bands %s with %dm resolution",
2663 23 : osBandNames.c_str(), nResolution));
2664 23 : if (nEPSGCode >= 32601 && nEPSGCode <= 32660)
2665 23 : osDesc += CPLSPrintf(", UTM %dN", nEPSGCode - 32600);
2666 0 : else if (nEPSGCode >= 32701 && nEPSGCode <= 32760)
2667 0 : osDesc += CPLSPrintf(", UTM %dS", nEPSGCode - 32700);
2668 : else
2669 0 : osDesc += CPLSPrintf(", EPSG:%d", nEPSGCode);
2670 23 : poDS->GDALDataset::SetMetadataItem(
2671 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
2672 : "SUBDATASETS");
2673 :
2674 23 : iSubDSNum++;
2675 : }
2676 : }
2677 :
2678 : /* Expose TCI or PREVIEW subdatasets */
2679 11 : if (bIsSafeCompact)
2680 : {
2681 10 : for (std::set<int>::const_iterator oIterEPSG = oSetEPSGCodes.begin();
2682 15 : oIterEPSG != oSetEPSGCodes.end(); ++oIterEPSG)
2683 : {
2684 5 : const int nEPSGCode = *oIterEPSG;
2685 5 : poDS->GDALDataset::SetMetadataItem(
2686 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
2687 : CPLSPrintf("%s:%s:TCI:EPSG_%d", pszPrefix, pszFilename,
2688 : nEPSGCode),
2689 : "SUBDATASETS");
2690 :
2691 5 : CPLString osDesc("True color image");
2692 5 : if (nEPSGCode >= 32601 && nEPSGCode <= 32660)
2693 5 : osDesc += CPLSPrintf(", UTM %dN", nEPSGCode - 32600);
2694 0 : else if (nEPSGCode >= 32701 && nEPSGCode <= 32760)
2695 0 : osDesc += CPLSPrintf(", UTM %dS", nEPSGCode - 32700);
2696 : else
2697 0 : osDesc += CPLSPrintf(", EPSG:%d", nEPSGCode);
2698 5 : poDS->GDALDataset::SetMetadataItem(
2699 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
2700 : "SUBDATASETS");
2701 :
2702 5 : iSubDSNum++;
2703 : }
2704 : }
2705 : else
2706 : {
2707 10 : for (std::set<int>::const_iterator oIterEPSG = oSetEPSGCodes.begin();
2708 14 : oIterEPSG != oSetEPSGCodes.end(); ++oIterEPSG)
2709 : {
2710 4 : const int nEPSGCode = *oIterEPSG;
2711 4 : poDS->GDALDataset::SetMetadataItem(
2712 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
2713 : CPLSPrintf("%s:%s:PREVIEW:EPSG_%d", pszPrefix, pszFilename,
2714 : nEPSGCode),
2715 : "SUBDATASETS");
2716 :
2717 4 : CPLString osDesc("RGB preview");
2718 4 : if (nEPSGCode >= 32601 && nEPSGCode <= 32660)
2719 4 : osDesc += CPLSPrintf(", UTM %dN", nEPSGCode - 32600);
2720 0 : else if (nEPSGCode >= 32701 && nEPSGCode <= 32760)
2721 0 : osDesc += CPLSPrintf(", UTM %dS", nEPSGCode - 32700);
2722 : else
2723 0 : osDesc += CPLSPrintf(", EPSG:%d", nEPSGCode);
2724 4 : poDS->GDALDataset::SetMetadataItem(
2725 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
2726 : "SUBDATASETS");
2727 :
2728 4 : iSubDSNum++;
2729 : }
2730 : }
2731 :
2732 11 : pszNodePath =
2733 : (eLevel == SENTINEL2_L1C)
2734 : ? "=Level-1C_User_Product.Geometric_Info.Product_Footprint."
2735 : "Product_Footprint.Global_Footprint.EXT_POS_LIST"
2736 : : "=Level-2A_User_Product.Geometric_Info.Product_Footprint."
2737 : "Product_Footprint.Global_Footprint.EXT_POS_LIST";
2738 11 : const char *pszPosList = CPLGetXMLValue(psRoot, pszNodePath, nullptr);
2739 11 : if (pszPosList != nullptr)
2740 : {
2741 18 : CPLString osPolygon = SENTINEL2GetPolygonWKTFromPosList(pszPosList);
2742 9 : if (!osPolygon.empty())
2743 9 : poDS->GDALDataset::SetMetadataItem("FOOTPRINT", osPolygon.c_str());
2744 : }
2745 :
2746 11 : return poDS;
2747 : }
2748 :
2749 : /************************************************************************/
2750 : /* SENTINEL2GetL1BCTileMetadata() */
2751 : /************************************************************************/
2752 :
2753 11 : static char **SENTINEL2GetL1BCTileMetadata(CPLXMLNode *psMainMTD)
2754 : {
2755 22 : CPLStringList aosList;
2756 :
2757 11 : CPLXMLNode *psRoot = CPLGetXMLNode(psMainMTD, "=Level-1C_Tile_ID");
2758 11 : if (psRoot == nullptr)
2759 : {
2760 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find =Level-1C_Tile_ID");
2761 0 : return nullptr;
2762 : }
2763 11 : CPLXMLNode *psGeneralInfo = CPLGetXMLNode(psRoot, "General_Info");
2764 11 : for (CPLXMLNode *psIter =
2765 11 : (psGeneralInfo ? psGeneralInfo->psChild : nullptr);
2766 58 : psIter != nullptr; psIter = psIter->psNext)
2767 : {
2768 47 : if (psIter->eType != CXT_Element)
2769 0 : continue;
2770 47 : const char *pszValue = CPLGetXMLValue(psIter, nullptr, nullptr);
2771 47 : if (pszValue != nullptr)
2772 : {
2773 38 : aosList.AddNameValue(psIter->pszValue, pszValue);
2774 : }
2775 : }
2776 :
2777 11 : CPLXMLNode *psQII = CPLGetXMLNode(psRoot, "Quality_Indicators_Info");
2778 11 : if (psQII != nullptr)
2779 : {
2780 9 : CPLXMLNode *psICCQI = CPLGetXMLNode(psQII, "Image_Content_QI");
2781 9 : for (CPLXMLNode *psIter = (psICCQI ? psICCQI->psChild : nullptr);
2782 27 : psIter != nullptr; psIter = psIter->psNext)
2783 : {
2784 18 : if (psIter->eType != CXT_Element)
2785 0 : continue;
2786 18 : if (psIter->psChild != nullptr &&
2787 18 : psIter->psChild->eType == CXT_Text)
2788 : {
2789 18 : aosList.AddNameValue(psIter->pszValue,
2790 18 : psIter->psChild->pszValue);
2791 : }
2792 : }
2793 : }
2794 :
2795 11 : return aosList.StealList();
2796 : }
2797 :
2798 : /************************************************************************/
2799 : /* OpenL1CTile() */
2800 : /************************************************************************/
2801 :
2802 14 : GDALDataset *SENTINEL2Dataset::OpenL1CTile(const char *pszFilename,
2803 : CPLXMLNode **ppsRootMainMTD,
2804 : int nResolutionOfInterest,
2805 : std::set<CPLString> *poBandSet)
2806 : {
2807 14 : CPLXMLNode *psRoot = CPLParseXMLFile(pszFilename);
2808 14 : if (psRoot == nullptr)
2809 : {
2810 3 : CPLDebug("SENTINEL2", "Cannot XML parse %s", pszFilename);
2811 3 : return nullptr;
2812 : }
2813 :
2814 11 : char *pszOriginalXML = CPLSerializeXMLTree(psRoot);
2815 22 : CPLString osOriginalXML;
2816 11 : if (pszOriginalXML)
2817 11 : osOriginalXML = pszOriginalXML;
2818 11 : CPLFree(pszOriginalXML);
2819 :
2820 22 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRoot);
2821 11 : CPLStripXMLNamespace(psRoot, nullptr, TRUE);
2822 :
2823 22 : std::set<int> oSetResolutions;
2824 22 : std::map<int, std::set<CPLString>> oMapResolutionsToBands;
2825 11 : char **papszMD = nullptr;
2826 11 : SENTINEL2GetResolutionSetAndMainMDFromGranule(
2827 : pszFilename, "Level-1C_User_Product", nResolutionOfInterest,
2828 : oSetResolutions, oMapResolutionsToBands, papszMD, ppsRootMainMTD);
2829 11 : if (poBandSet != nullptr)
2830 8 : *poBandSet = oMapResolutionsToBands[nResolutionOfInterest];
2831 :
2832 11 : SENTINEL2DatasetContainer *poDS = new SENTINEL2DatasetContainer();
2833 :
2834 11 : char **papszGranuleMD = SENTINEL2GetL1BCTileMetadata(psRoot);
2835 11 : papszMD = CSLMerge(papszMD, papszGranuleMD);
2836 11 : CSLDestroy(papszGranuleMD);
2837 :
2838 : // Remove CLOUD_COVERAGE_ASSESSMENT that comes from main metadata, if
2839 : // granule CLOUDY_PIXEL_PERCENTAGE is present.
2840 20 : if (CSLFetchNameValue(papszMD, "CLOUDY_PIXEL_PERCENTAGE") != nullptr &&
2841 9 : CSLFetchNameValue(papszMD, "CLOUD_COVERAGE_ASSESSMENT") != nullptr)
2842 : {
2843 7 : papszMD =
2844 7 : CSLSetNameValue(papszMD, "CLOUD_COVERAGE_ASSESSMENT", nullptr);
2845 : }
2846 :
2847 11 : poDS->GDALDataset::SetMetadata(papszMD);
2848 11 : CSLDestroy(papszMD);
2849 :
2850 11 : if (!osOriginalXML.empty())
2851 : {
2852 : char *apszXMLMD[2];
2853 11 : apszXMLMD[0] = const_cast<char *>(osOriginalXML.c_str());
2854 11 : apszXMLMD[1] = nullptr;
2855 11 : poDS->GDALDataset::SetMetadata(apszXMLMD, "xml:SENTINEL2");
2856 : }
2857 :
2858 : /* Create subdatsets per resolution (10, 20, 60m) */
2859 11 : int iSubDSNum = 1;
2860 36 : for (std::set<int>::const_iterator oIterRes = oSetResolutions.begin();
2861 61 : oIterRes != oSetResolutions.end(); ++oIterRes)
2862 : {
2863 25 : const int nResolution = *oIterRes;
2864 :
2865 25 : poDS->GDALDataset::SetMetadataItem(
2866 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
2867 : CPLSPrintf("%s:%s:%dm", "SENTINEL2_L1C_TILE", pszFilename,
2868 : nResolution),
2869 : "SUBDATASETS");
2870 :
2871 : CPLString osBandNames = SENTINEL2GetBandListForResolution(
2872 50 : oMapResolutionsToBands[nResolution]);
2873 :
2874 : CPLString osDesc(CPLSPrintf("Bands %s with %dm resolution",
2875 25 : osBandNames.c_str(), nResolution));
2876 25 : poDS->GDALDataset::SetMetadataItem(
2877 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
2878 : "SUBDATASETS");
2879 :
2880 25 : iSubDSNum++;
2881 : }
2882 :
2883 : /* Expose PREVIEW subdataset */
2884 11 : poDS->GDALDataset::SetMetadataItem(
2885 : CPLSPrintf("SUBDATASET_%d_NAME", iSubDSNum),
2886 : CPLSPrintf("%s:%s:PREVIEW", "SENTINEL2_L1C_TILE", pszFilename),
2887 : "SUBDATASETS");
2888 :
2889 11 : CPLString osDesc("RGB preview");
2890 11 : poDS->GDALDataset::SetMetadataItem(
2891 : CPLSPrintf("SUBDATASET_%d_DESC", iSubDSNum), osDesc.c_str(),
2892 : "SUBDATASETS");
2893 :
2894 11 : return poDS;
2895 : }
2896 :
2897 : /************************************************************************/
2898 : /* SENTINEL2GetOption() */
2899 : /************************************************************************/
2900 :
2901 55 : static const char *SENTINEL2GetOption(GDALOpenInfo *poOpenInfo,
2902 : const char *pszName,
2903 : const char *pszDefaultVal)
2904 : {
2905 : #ifdef GDAL_DMD_OPENOPTIONLIST
2906 : const char *pszVal =
2907 55 : CSLFetchNameValue(poOpenInfo->papszOpenOptions, pszName);
2908 55 : if (pszVal != nullptr)
2909 3 : return pszVal;
2910 : #else
2911 : (void)poOpenInfo;
2912 : #endif
2913 52 : return CPLGetConfigOption(CPLSPrintf("SENTINEL2_%s", pszName),
2914 52 : pszDefaultVal);
2915 : }
2916 :
2917 : /************************************************************************/
2918 : /* LaunderUnit() */
2919 : /************************************************************************/
2920 :
2921 143 : static CPLString LaunderUnit(const char *pszUnit)
2922 : {
2923 143 : CPLString osUnit;
2924 1144 : for (int i = 0; pszUnit[i] != '\0';)
2925 : {
2926 1001 : if (strncmp(pszUnit + i,
2927 : "\xC2"
2928 : "\xB2",
2929 : 2) == 0) /* square / 2 */
2930 : {
2931 143 : i += 2;
2932 143 : osUnit += "2";
2933 : }
2934 858 : else if (strncmp(pszUnit + i,
2935 : "\xC2"
2936 : "\xB5",
2937 : 2) == 0) /* micro */
2938 : {
2939 143 : i += 2;
2940 143 : osUnit += "u";
2941 : }
2942 : else
2943 : {
2944 715 : osUnit += pszUnit[i];
2945 715 : i++;
2946 : }
2947 : }
2948 143 : return osUnit;
2949 : }
2950 :
2951 : /************************************************************************/
2952 : /* SENTINEL2GetTileInfo() */
2953 : /************************************************************************/
2954 :
2955 35 : static bool SENTINEL2GetTileInfo(const char *pszFilename, int *pnWidth,
2956 : int *pnHeight, int *pnBits)
2957 : {
2958 : static const unsigned char jp2_box_jp[] = {0x6a, 0x50, 0x20,
2959 : 0x20}; /* 'jP ' */
2960 35 : VSILFILE *fp = VSIFOpenL(pszFilename, "rb");
2961 35 : if (fp == nullptr)
2962 2 : return false;
2963 : GByte abyHeader[8];
2964 33 : if (VSIFReadL(abyHeader, 8, 1, fp) != 1)
2965 : {
2966 0 : VSIFCloseL(fp);
2967 0 : return false;
2968 : }
2969 33 : if (memcmp(abyHeader + 4, jp2_box_jp, 4) == 0)
2970 : {
2971 3 : bool bRet = false;
2972 : /* Just parse the ihdr box instead of doing a full dataset opening */
2973 3 : GDALJP2Box oBox(fp);
2974 3 : if (oBox.ReadFirst())
2975 : {
2976 9 : while (strlen(oBox.GetType()) > 0)
2977 : {
2978 9 : if (EQUAL(oBox.GetType(), "jp2h"))
2979 : {
2980 6 : GDALJP2Box oChildBox(fp);
2981 3 : if (!oChildBox.ReadFirstChild(&oBox))
2982 0 : break;
2983 :
2984 3 : while (strlen(oChildBox.GetType()) > 0)
2985 : {
2986 3 : if (EQUAL(oChildBox.GetType(), "ihdr"))
2987 : {
2988 3 : GByte *pabyData = oChildBox.ReadBoxData();
2989 3 : GIntBig nLength = oChildBox.GetDataLength();
2990 3 : if (pabyData != nullptr && nLength >= 4 + 4 + 2 + 1)
2991 : {
2992 3 : bRet = true;
2993 3 : if (pnHeight)
2994 : {
2995 1 : memcpy(pnHeight, pabyData, 4);
2996 1 : CPL_MSBPTR32(pnHeight);
2997 : }
2998 3 : if (pnWidth != nullptr)
2999 : {
3000 : // cppcheck-suppress nullPointer
3001 1 : memcpy(pnWidth, pabyData + 4, 4);
3002 1 : CPL_MSBPTR32(pnWidth);
3003 : }
3004 3 : if (pnBits)
3005 : {
3006 3 : GByte byPBC = pabyData[4 + 4 + 2];
3007 3 : if (byPBC != 255)
3008 : {
3009 3 : *pnBits = 1 + (byPBC & 0x7f);
3010 : }
3011 : else
3012 0 : *pnBits = 0;
3013 : }
3014 : }
3015 3 : CPLFree(pabyData);
3016 3 : break;
3017 : }
3018 0 : if (!oChildBox.ReadNextChild(&oBox))
3019 0 : break;
3020 : }
3021 3 : break;
3022 : }
3023 :
3024 6 : if (!oBox.ReadNext())
3025 0 : break;
3026 : }
3027 : }
3028 3 : VSIFCloseL(fp);
3029 3 : return bRet;
3030 : }
3031 : else /* for unit tests, we use TIFF */
3032 : {
3033 30 : VSIFCloseL(fp);
3034 30 : GDALDataset *poDS = (GDALDataset *)GDALOpen(pszFilename, GA_ReadOnly);
3035 30 : bool bRet = false;
3036 30 : if (poDS != nullptr)
3037 : {
3038 30 : if (poDS->GetRasterCount() != 0)
3039 : {
3040 30 : bRet = true;
3041 30 : if (pnWidth)
3042 0 : *pnWidth = poDS->GetRasterXSize();
3043 30 : if (pnHeight)
3044 0 : *pnHeight = poDS->GetRasterYSize();
3045 30 : if (pnBits)
3046 : {
3047 : const char *pszNBits =
3048 30 : poDS->GetRasterBand(1)->GetMetadataItem(
3049 30 : "NBITS", "IMAGE_STRUCTURE");
3050 30 : if (pszNBits == nullptr)
3051 : {
3052 : GDALDataType eDT =
3053 4 : poDS->GetRasterBand(1)->GetRasterDataType();
3054 4 : pszNBits = CPLSPrintf("%d", GDALGetDataTypeSize(eDT));
3055 : }
3056 30 : *pnBits = atoi(pszNBits);
3057 : }
3058 : }
3059 30 : GDALClose(poDS);
3060 : }
3061 30 : return bRet;
3062 : }
3063 : }
3064 :
3065 : /************************************************************************/
3066 : /* OpenL1C_L2ASubdataset() */
3067 : /************************************************************************/
3068 :
3069 79 : GDALDataset *SENTINEL2Dataset::OpenL1C_L2ASubdataset(GDALOpenInfo *poOpenInfo,
3070 : SENTINEL2Level eLevel)
3071 : {
3072 158 : CPLString osFilename;
3073 79 : const char *pszPrefix =
3074 : (eLevel == SENTINEL2_L1C) ? "SENTINEL2_L1C" : "SENTINEL2_L2A";
3075 79 : CPLAssert(STARTS_WITH_CI(poOpenInfo->pszFilename, pszPrefix));
3076 79 : osFilename = poOpenInfo->pszFilename + strlen(pszPrefix) + 1;
3077 79 : const char *pszEPSGCode = strrchr(osFilename.c_str(), ':');
3078 79 : if (pszEPSGCode == nullptr || pszEPSGCode == osFilename.c_str())
3079 : {
3080 10 : CPLError(CE_Failure, CPLE_AppDefined,
3081 : "Invalid syntax for %s:", pszPrefix);
3082 10 : return nullptr;
3083 : }
3084 69 : osFilename[(size_t)(pszEPSGCode - osFilename.c_str())] = '\0';
3085 69 : const char *pszPrecision = strrchr(osFilename.c_str(), ':');
3086 69 : if (pszPrecision == nullptr)
3087 : {
3088 10 : CPLError(CE_Failure, CPLE_AppDefined,
3089 : "Invalid syntax for %s:", pszPrefix);
3090 10 : return nullptr;
3091 : }
3092 :
3093 59 : if (!STARTS_WITH_CI(pszEPSGCode + 1, "EPSG_"))
3094 : {
3095 5 : CPLError(CE_Failure, CPLE_AppDefined,
3096 : "Invalid syntax for %s:", pszPrefix);
3097 5 : return nullptr;
3098 : }
3099 :
3100 54 : const int nSubDSEPSGCode = atoi(pszEPSGCode + 1 + strlen("EPSG_"));
3101 54 : const bool bIsPreview = STARTS_WITH_CI(pszPrecision + 1, "PREVIEW");
3102 54 : const bool bIsTCI = STARTS_WITH_CI(pszPrecision + 1, "TCI");
3103 105 : const int nSubDSPrecision = (bIsPreview) ? 320
3104 51 : : (bIsTCI) ? 10
3105 49 : : atoi(pszPrecision + 1);
3106 54 : if (!bIsTCI && !bIsPreview && nSubDSPrecision != 10 &&
3107 20 : nSubDSPrecision != 20 && nSubDSPrecision != 60)
3108 : {
3109 5 : CPLError(CE_Failure, CPLE_NotSupported, "Unsupported precision: %d",
3110 : nSubDSPrecision);
3111 5 : return nullptr;
3112 : }
3113 :
3114 49 : osFilename.resize(pszPrecision - osFilename.c_str());
3115 98 : std::vector<CPLString> aosNonJP2Files;
3116 49 : aosNonJP2Files.push_back(osFilename);
3117 :
3118 49 : CPLXMLNode *psRoot = CPLParseXMLFile(osFilename);
3119 49 : if (psRoot == nullptr)
3120 : {
3121 7 : CPLDebug("SENTINEL2", "Cannot XML parse %s", osFilename.c_str());
3122 7 : return nullptr;
3123 : }
3124 :
3125 42 : char *pszOriginalXML = CPLSerializeXMLTree(psRoot);
3126 84 : CPLString osOriginalXML;
3127 42 : if (pszOriginalXML)
3128 42 : osOriginalXML = pszOriginalXML;
3129 42 : CPLFree(pszOriginalXML);
3130 :
3131 84 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRoot);
3132 42 : CPLStripXMLNamespace(psRoot, nullptr, TRUE);
3133 :
3134 : CPLXMLNode *psProductInfo =
3135 : eLevel == SENTINEL2_L1C
3136 42 : ? CPLGetXMLNode(psRoot,
3137 : "=Level-1C_User_Product.General_Info.Product_Info")
3138 18 : : CPLGetXMLNode(psRoot,
3139 42 : "=Level-2A_User_Product.General_Info.Product_Info");
3140 42 : if (psProductInfo == nullptr && eLevel == SENTINEL2_L2A)
3141 : {
3142 11 : psProductInfo = CPLGetXMLNode(
3143 : psRoot, "=Level-2A_User_Product.General_Info.L2A_Product_Info");
3144 : }
3145 42 : if (psProductInfo == nullptr)
3146 : {
3147 1 : CPLDebug("SENTINEL2", "Product Info not found");
3148 1 : return nullptr;
3149 : }
3150 :
3151 : const bool bIsSafeCompact =
3152 41 : EQUAL(CPLGetXMLValue(psProductInfo, "Query_Options.PRODUCT_FORMAT", ""),
3153 : "SAFE_COMPACT");
3154 :
3155 : const char *pszProductURI =
3156 41 : CPLGetXMLValue(psProductInfo, "PRODUCT_URI", nullptr);
3157 41 : SENTINEL2ProductType pType = MSI2A;
3158 41 : if (pszProductURI == nullptr)
3159 : {
3160 : pszProductURI =
3161 32 : CPLGetXMLValue(psProductInfo, "PRODUCT_URI_2A", nullptr);
3162 32 : pType = MSI2Ap;
3163 : }
3164 41 : if (pszProductURI == nullptr)
3165 27 : pszProductURI = "";
3166 :
3167 82 : std::vector<CPLString> aosGranuleList;
3168 82 : std::map<int, std::set<CPLString>> oMapResolutionsToBands;
3169 82 : std::vector<L1CSafeCompatGranuleDescription> aoL1CSafeCompactGranuleList;
3170 41 : if (bIsSafeCompact)
3171 : {
3172 308 : for (unsigned int i = 0; i < NB_BANDS; ++i)
3173 : {
3174 : // L2 does not contain B10
3175 286 : if (i == 10 && eLevel == SENTINEL2_L2A)
3176 12 : continue;
3177 274 : const SENTINEL2BandDescription *psBandDesc = &asBandDesc[i];
3178 : CPLString osName =
3179 548 : psBandDesc->pszBandName + 1; /* skip B character */
3180 274 : if (atoi(osName) < 10)
3181 220 : osName = "0" + osName;
3182 274 : oMapResolutionsToBands[psBandDesc->nResolution].insert(osName);
3183 : }
3184 22 : if (eLevel == SENTINEL2_L2A)
3185 : {
3186 156 : for (const auto &sL2ABandDesc : asL2ABandDesc)
3187 : {
3188 144 : oMapResolutionsToBands[sL2ABandDesc.nResolution].insert(
3189 144 : sL2ABandDesc.pszBandName);
3190 : }
3191 : }
3192 32 : if (eLevel == SENTINEL2_L1C &&
3193 10 : !SENTINEL2GetGranuleList_L1CSafeCompact(
3194 : psRoot, osFilename, aoL1CSafeCompactGranuleList))
3195 : {
3196 0 : CPLDebug("SENTINEL2", "Failed to get granule list");
3197 0 : return nullptr;
3198 : }
3199 34 : if (eLevel == SENTINEL2_L2A &&
3200 12 : !SENTINEL2GetGranuleList_L2ASafeCompact(
3201 : psRoot, osFilename, aoL1CSafeCompactGranuleList))
3202 : {
3203 0 : CPLDebug("SENTINEL2", "Failed to get granule list");
3204 0 : return nullptr;
3205 : }
3206 54 : for (size_t i = 0; i < aoL1CSafeCompactGranuleList.size(); ++i)
3207 : {
3208 32 : aosGranuleList.push_back(
3209 32 : aoL1CSafeCompactGranuleList[i].osMTDTLPath);
3210 : }
3211 : }
3212 19 : else if (!SENTINEL2GetGranuleList(
3213 : psRoot, eLevel, osFilename, aosGranuleList, nullptr,
3214 : (eLevel == SENTINEL2_L1C) ? nullptr : &oMapResolutionsToBands))
3215 : {
3216 1 : CPLDebug("SENTINEL2", "Failed to get granule list");
3217 1 : return nullptr;
3218 : }
3219 :
3220 80 : std::vector<CPLString> aosBands;
3221 80 : std::set<CPLString> oSetBands;
3222 40 : if (bIsPreview || bIsTCI)
3223 : {
3224 5 : aosBands.push_back("04");
3225 5 : aosBands.push_back("03");
3226 5 : aosBands.push_back("02");
3227 : }
3228 35 : else if (eLevel == SENTINEL2_L1C && !bIsSafeCompact)
3229 : {
3230 : CPLXMLNode *psBandList =
3231 10 : CPLGetXMLNode(psRoot, "=Level-1C_User_Product.General_Info.Product_"
3232 : "Info.Query_Options.Band_List");
3233 10 : if (psBandList == nullptr)
3234 : {
3235 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s",
3236 : "Query_Options.Band_List");
3237 0 : return nullptr;
3238 : }
3239 :
3240 116 : for (CPLXMLNode *psIter = psBandList->psChild; psIter != nullptr;
3241 106 : psIter = psIter->psNext)
3242 : {
3243 106 : if (psIter->eType != CXT_Element ||
3244 106 : !EQUAL(psIter->pszValue, "BAND_NAME"))
3245 72 : continue;
3246 106 : const char *pszBandName = CPLGetXMLValue(psIter, nullptr, "");
3247 : const SENTINEL2BandDescription *psBandDesc =
3248 106 : SENTINEL2GetBandDesc(pszBandName);
3249 106 : if (psBandDesc == nullptr)
3250 : {
3251 0 : CPLDebug("SENTINEL2", "Unknown band name %s", pszBandName);
3252 0 : continue;
3253 : }
3254 106 : if (psBandDesc->nResolution != nSubDSPrecision)
3255 72 : continue;
3256 : CPLString osName =
3257 68 : psBandDesc->pszBandName + 1; /* skip B character */
3258 34 : if (atoi(osName) < 10)
3259 30 : osName = "0" + osName;
3260 34 : oSetBands.insert(osName);
3261 : }
3262 10 : if (oSetBands.empty())
3263 : {
3264 0 : CPLDebug("SENTINEL2", "Band set is empty");
3265 0 : return nullptr;
3266 10 : }
3267 : }
3268 : else
3269 : {
3270 25 : oSetBands = oMapResolutionsToBands[nSubDSPrecision];
3271 : }
3272 :
3273 40 : if (aosBands.empty())
3274 : {
3275 192 : for (std::set<CPLString>::const_iterator oIterBandnames =
3276 35 : oSetBands.begin();
3277 419 : oIterBandnames != oSetBands.end(); ++oIterBandnames)
3278 : {
3279 192 : aosBands.push_back(*oIterBandnames);
3280 : }
3281 : /* Put 2=Blue, 3=Green, 4=Band bands in RGB order for conveniency */
3282 82 : if (aosBands.size() >= 3 && aosBands[0] == "02" &&
3283 82 : aosBands[1] == "03" && aosBands[2] == "04")
3284 : {
3285 17 : aosBands[0] = "04";
3286 17 : aosBands[2] = "02";
3287 : }
3288 : }
3289 :
3290 : /* -------------------------------------------------------------------- */
3291 : /* Create dataset. */
3292 : /* -------------------------------------------------------------------- */
3293 :
3294 40 : char **papszMD = SENTINEL2GetUserProductMetadata(
3295 : psRoot, (eLevel == SENTINEL2_L1C) ? "Level-1C_User_Product"
3296 : : "Level-2A_User_Product");
3297 :
3298 : const int nSaturatedVal =
3299 40 : atoi(CSLFetchNameValueDef(papszMD, "SPECIAL_VALUE_SATURATED", "-1"));
3300 : const int nNodataVal =
3301 40 : atoi(CSLFetchNameValueDef(papszMD, "SPECIAL_VALUE_NODATA", "-1"));
3302 :
3303 : const bool bAlpha =
3304 40 : CPLTestBool(SENTINEL2GetOption(poOpenInfo, "ALPHA", "FALSE"));
3305 :
3306 40 : SENTINEL2Dataset *poDS = CreateL1CL2ADataset(
3307 : eLevel, pType, bIsSafeCompact, aosGranuleList,
3308 : aoL1CSafeCompactGranuleList, aosNonJP2Files, nSubDSPrecision,
3309 : bIsPreview, bIsTCI, nSubDSEPSGCode, bAlpha, aosBands, nSaturatedVal,
3310 80 : nNodataVal, CPLString(pszProductURI));
3311 40 : if (poDS == nullptr)
3312 : {
3313 12 : CSLDestroy(papszMD);
3314 12 : return nullptr;
3315 : }
3316 :
3317 28 : if (!osOriginalXML.empty())
3318 : {
3319 28 : char *apszXMLMD[2] = {nullptr};
3320 28 : apszXMLMD[0] = const_cast<char *>(osOriginalXML.c_str());
3321 28 : apszXMLMD[1] = nullptr;
3322 28 : poDS->GDALDataset::SetMetadata(apszXMLMD, "xml:SENTINEL2");
3323 : }
3324 :
3325 28 : poDS->GDALDataset::SetMetadata(papszMD);
3326 28 : CSLDestroy(papszMD);
3327 :
3328 : /* -------------------------------------------------------------------- */
3329 : /* Add extra band metadata. */
3330 : /* -------------------------------------------------------------------- */
3331 28 : poDS->AddL1CL2ABandMetadata(eLevel, psRoot, aosBands);
3332 :
3333 : /* -------------------------------------------------------------------- */
3334 : /* Initialize overview information. */
3335 : /* -------------------------------------------------------------------- */
3336 28 : poDS->SetDescription(poOpenInfo->pszFilename);
3337 28 : CPLString osOverviewFile;
3338 28 : if (bIsPreview)
3339 : osOverviewFile = CPLSPrintf("%s_PREVIEW_EPSG_%d.tif.ovr",
3340 3 : osFilename.c_str(), nSubDSEPSGCode);
3341 25 : else if (bIsTCI)
3342 : osOverviewFile = CPLSPrintf("%s_TCI_EPSG_%d.tif.ovr",
3343 2 : osFilename.c_str(), nSubDSEPSGCode);
3344 : else
3345 : osOverviewFile =
3346 : CPLSPrintf("%s_%dm_EPSG_%d.tif.ovr", osFilename.c_str(),
3347 23 : nSubDSPrecision, nSubDSEPSGCode);
3348 28 : poDS->SetMetadataItem("OVERVIEW_FILE", osOverviewFile, "OVERVIEWS");
3349 28 : poDS->oOvManager.Initialize(poDS, ":::VIRTUAL:::");
3350 :
3351 28 : return poDS;
3352 : }
3353 :
3354 : /************************************************************************/
3355 : /* AddL1CL2ABandMetadata() */
3356 : /************************************************************************/
3357 :
3358 34 : void SENTINEL2Dataset::AddL1CL2ABandMetadata(
3359 : SENTINEL2Level eLevel, CPLXMLNode *psRoot,
3360 : const std::vector<CPLString> &aosBands)
3361 : {
3362 : CPLXMLNode *psIC =
3363 34 : CPLGetXMLNode(psRoot, (eLevel == SENTINEL2_L1C)
3364 : ? "=Level-1C_User_Product.General_Info."
3365 : "Product_Image_Characteristics"
3366 : : "=Level-2A_User_Product.General_Info."
3367 : "Product_Image_Characteristics");
3368 34 : if (psIC == nullptr)
3369 : {
3370 8 : psIC = CPLGetXMLNode(psRoot, "=Level-2A_User_Product.General_Info.L2A_"
3371 : "Product_Image_Characteristics");
3372 : }
3373 34 : if (psIC != nullptr)
3374 : {
3375 : CPLXMLNode *psSIL =
3376 32 : CPLGetXMLNode(psIC, "Reflectance_Conversion.Solar_Irradiance_List");
3377 32 : if (psSIL != nullptr)
3378 : {
3379 448 : for (CPLXMLNode *psIter = psSIL->psChild; psIter != nullptr;
3380 416 : psIter = psIter->psNext)
3381 : {
3382 416 : if (psIter->eType != CXT_Element ||
3383 416 : !EQUAL(psIter->pszValue, "SOLAR_IRRADIANCE"))
3384 : {
3385 0 : continue;
3386 : }
3387 : const char *pszBandId =
3388 416 : CPLGetXMLValue(psIter, "bandId", nullptr);
3389 416 : const char *pszUnit = CPLGetXMLValue(psIter, "unit", nullptr);
3390 416 : const char *pszValue = CPLGetXMLValue(psIter, nullptr, nullptr);
3391 416 : if (pszBandId != nullptr && pszUnit != nullptr &&
3392 : pszValue != nullptr)
3393 : {
3394 416 : int nIdx = atoi(pszBandId);
3395 416 : if (nIdx >= 0 && nIdx < (int)NB_BANDS)
3396 : {
3397 2089 : for (int i = 0; i < nBands; i++)
3398 : {
3399 1816 : GDALRasterBand *poBand = GetRasterBand(i + 1);
3400 : const char *pszBandName =
3401 1816 : poBand->GetMetadataItem("BANDNAME");
3402 1816 : if (pszBandName != nullptr &&
3403 1806 : EQUAL(asBandDesc[nIdx].pszBandName,
3404 : pszBandName))
3405 : {
3406 143 : poBand->GDALRasterBand::SetMetadataItem(
3407 : "SOLAR_IRRADIANCE", pszValue);
3408 143 : poBand->GDALRasterBand::SetMetadataItem(
3409 : "SOLAR_IRRADIANCE_UNIT",
3410 286 : LaunderUnit(pszUnit));
3411 143 : break;
3412 : }
3413 : }
3414 : }
3415 : }
3416 : }
3417 : }
3418 :
3419 32 : CPLXMLNode *psOL = CPLGetXMLNode(
3420 : psIC, (eLevel == SENTINEL2_L1C) ? "Radiometric_Offset_List"
3421 : : "BOA_ADD_OFFSET_VALUES_LIST");
3422 32 : if (psOL != nullptr)
3423 : {
3424 56 : for (CPLXMLNode *psIter = psOL->psChild; psIter != nullptr;
3425 52 : psIter = psIter->psNext)
3426 : {
3427 52 : if (psIter->eType != CXT_Element ||
3428 52 : !EQUAL(psIter->pszValue, (eLevel == SENTINEL2_L1C)
3429 : ? "RADIO_ADD_OFFSET"
3430 : : "BOA_ADD_OFFSET"))
3431 : {
3432 0 : continue;
3433 : }
3434 : const char *pszBandId =
3435 52 : CPLGetXMLValue(psIter, "band_id", nullptr);
3436 52 : const char *pszValue = CPLGetXMLValue(psIter, nullptr, nullptr);
3437 52 : if (pszBandId != nullptr && pszValue != nullptr)
3438 : {
3439 52 : int nIdx = atoi(pszBandId);
3440 52 : if (nIdx >= 0 && nIdx < (int)NB_BANDS)
3441 : {
3442 303 : for (int i = 0; i < nBands; i++)
3443 : {
3444 271 : GDALRasterBand *poBand = GetRasterBand(i + 1);
3445 : const char *pszBandName =
3446 271 : poBand->GetMetadataItem("BANDNAME");
3447 271 : if (pszBandName != nullptr &&
3448 271 : EQUAL(asBandDesc[nIdx].pszBandName,
3449 : pszBandName))
3450 : {
3451 20 : poBand->GDALRasterBand::SetMetadataItem(
3452 : (eLevel == SENTINEL2_L1C)
3453 : ? "RADIO_ADD_OFFSET"
3454 : : "BOA_ADD_OFFSET",
3455 : pszValue);
3456 20 : break;
3457 : }
3458 : }
3459 : }
3460 : }
3461 : }
3462 : }
3463 : }
3464 :
3465 : /* -------------------------------------------------------------------- */
3466 : /* Add Scene Classification category values (L2A) */
3467 : /* -------------------------------------------------------------------- */
3468 34 : CPLXMLNode *psSCL = CPLGetXMLNode(
3469 : psRoot, "=Level-2A_User_Product.General_Info."
3470 : "Product_Image_Characteristics.Scene_Classification_List");
3471 34 : if (psSCL == nullptr)
3472 : {
3473 29 : psSCL = CPLGetXMLNode(
3474 : psRoot,
3475 : "=Level-2A_User_Product.General_Info."
3476 : "L2A_Product_Image_Characteristics.L2A_Scene_Classification_List");
3477 : }
3478 34 : int nSCLBand = 0;
3479 199 : for (int nBand = 1; nBand <= static_cast<int>(aosBands.size()); nBand++)
3480 : {
3481 172 : if (EQUAL(aosBands[nBand - 1], "SCL"))
3482 : {
3483 7 : nSCLBand = nBand;
3484 7 : break;
3485 : }
3486 : }
3487 34 : if (psSCL != nullptr && nSCLBand > 0)
3488 : {
3489 14 : std::vector<CPLString> osCategories;
3490 91 : for (CPLXMLNode *psIter = psSCL->psChild; psIter != nullptr;
3491 84 : psIter = psIter->psNext)
3492 : {
3493 84 : if (psIter->eType != CXT_Element ||
3494 84 : (!EQUAL(psIter->pszValue, "L2A_Scene_Classification_ID") &&
3495 36 : !EQUAL(psIter->pszValue, "Scene_Classification_ID")))
3496 : {
3497 0 : continue;
3498 : }
3499 : const char *pszText =
3500 84 : CPLGetXMLValue(psIter, "SCENE_CLASSIFICATION_TEXT", nullptr);
3501 84 : if (pszText == nullptr)
3502 48 : pszText = CPLGetXMLValue(
3503 : psIter, "L2A_SCENE_CLASSIFICATION_TEXT", nullptr);
3504 : const char *pszIdx =
3505 84 : CPLGetXMLValue(psIter, "SCENE_CLASSIFICATION_INDEX", nullptr);
3506 84 : if (pszIdx == nullptr)
3507 48 : pszIdx = CPLGetXMLValue(
3508 : psIter, "L2A_SCENE_CLASSIFICATION_INDEX", nullptr);
3509 84 : if (pszText && pszIdx && atoi(pszIdx) >= 0 && atoi(pszIdx) < 100)
3510 : {
3511 84 : int nIdx = atoi(pszIdx);
3512 84 : if (nIdx >= (int)osCategories.size())
3513 84 : osCategories.resize(nIdx + 1);
3514 84 : if (STARTS_WITH_CI(pszText, "SC_"))
3515 84 : osCategories[nIdx] = pszText + 3;
3516 : else
3517 0 : osCategories[nIdx] = pszText;
3518 : }
3519 : }
3520 : char **papszCategories =
3521 7 : (char **)CPLCalloc(osCategories.size() + 1, sizeof(char *));
3522 91 : for (size_t i = 0; i < osCategories.size(); i++)
3523 84 : papszCategories[i] = CPLStrdup(osCategories[i]);
3524 7 : GetRasterBand(nSCLBand)->SetCategoryNames(papszCategories);
3525 7 : CSLDestroy(papszCategories);
3526 : }
3527 34 : }
3528 :
3529 : /************************************************************************/
3530 : /* CreateL1CL2ADataset() */
3531 : /************************************************************************/
3532 :
3533 48 : SENTINEL2Dataset *SENTINEL2Dataset::CreateL1CL2ADataset(
3534 : SENTINEL2Level eLevel, SENTINEL2ProductType pType, bool bIsSafeCompact,
3535 : const std::vector<CPLString> &aosGranuleList,
3536 : const std::vector<L1CSafeCompatGranuleDescription>
3537 : &aoL1CSafeCompactGranuleList,
3538 : std::vector<CPLString> &aosNonJP2Files, int nSubDSPrecision,
3539 : bool bIsPreview, bool bIsTCI,
3540 : int nSubDSEPSGCode, /* or -1 if not known at this point */
3541 : bool bAlpha, const std::vector<CPLString> &aosBands, int nSaturatedVal,
3542 : int nNodataVal, const CPLString &osProductURI)
3543 : {
3544 :
3545 : /* Iterate over granule metadata to know the layer extent */
3546 : /* and the location of each granule */
3547 48 : double dfMinX = 1.0e20;
3548 48 : double dfMinY = 1.0e20;
3549 48 : double dfMaxX = -1.0e20;
3550 48 : double dfMaxY = -1.0e20;
3551 96 : std::vector<SENTINEL2GranuleInfo> aosGranuleInfoList;
3552 48 : const int nDesiredResolution = (bIsPreview || bIsTCI) ? 0 : nSubDSPrecision;
3553 :
3554 48 : if (bIsSafeCompact)
3555 : {
3556 22 : CPLAssert(aosGranuleList.size() == aoL1CSafeCompactGranuleList.size());
3557 : }
3558 :
3559 116 : for (size_t i = 0; i < aosGranuleList.size(); i++)
3560 : {
3561 68 : int nEPSGCode = 0;
3562 68 : double dfULX = 0.0;
3563 68 : double dfULY = 0.0;
3564 68 : int nResolution = 0;
3565 68 : int nWidth = 0;
3566 68 : int nHeight = 0;
3567 68 : if (SENTINEL2GetGranuleInfo(eLevel, aosGranuleList[i],
3568 : nDesiredResolution, &nEPSGCode, &dfULX,
3569 65 : &dfULY, &nResolution, &nWidth, &nHeight) &&
3570 117 : (nSubDSEPSGCode == nEPSGCode || nSubDSEPSGCode < 0) &&
3571 49 : nResolution != 0)
3572 : {
3573 48 : nSubDSEPSGCode = nEPSGCode;
3574 48 : aosNonJP2Files.push_back(aosGranuleList[i]);
3575 :
3576 48 : if (dfULX < dfMinX)
3577 35 : dfMinX = dfULX;
3578 48 : if (dfULY > dfMaxY)
3579 35 : dfMaxY = dfULY;
3580 48 : const double dfLRX = dfULX + nResolution * nWidth;
3581 48 : const double dfLRY = dfULY - nResolution * nHeight;
3582 48 : if (dfLRX > dfMaxX)
3583 42 : dfMaxX = dfLRX;
3584 48 : if (dfLRY < dfMinY)
3585 42 : dfMinY = dfLRY;
3586 :
3587 96 : SENTINEL2GranuleInfo oGranuleInfo;
3588 48 : oGranuleInfo.osPath = CPLGetPathSafe(aosGranuleList[i]);
3589 48 : if (bIsSafeCompact)
3590 : {
3591 : oGranuleInfo.osBandPrefixPath =
3592 22 : aoL1CSafeCompactGranuleList[i].osBandPrefixPath;
3593 : }
3594 48 : oGranuleInfo.dfMinX = dfULX;
3595 48 : oGranuleInfo.dfMinY = dfLRY;
3596 48 : oGranuleInfo.dfMaxX = dfLRX;
3597 48 : oGranuleInfo.dfMaxY = dfULY;
3598 48 : oGranuleInfo.nWidth = nWidth / (nSubDSPrecision / nResolution);
3599 48 : oGranuleInfo.nHeight = nHeight / (nSubDSPrecision / nResolution);
3600 48 : aosGranuleInfoList.push_back(oGranuleInfo);
3601 : }
3602 : }
3603 48 : if (dfMinX > dfMaxX)
3604 : {
3605 13 : CPLError(CE_Failure, CPLE_NotSupported,
3606 : "No granule found for EPSG code %d", nSubDSEPSGCode);
3607 13 : return nullptr;
3608 : }
3609 :
3610 35 : const int nRasterXSize =
3611 35 : static_cast<int>((dfMaxX - dfMinX) / nSubDSPrecision + 0.5);
3612 35 : const int nRasterYSize =
3613 35 : static_cast<int>((dfMaxY - dfMinY) / nSubDSPrecision + 0.5);
3614 35 : SENTINEL2Dataset *poDS = new SENTINEL2Dataset(nRasterXSize, nRasterYSize);
3615 :
3616 : #if defined(__GNUC__)
3617 : #pragma GCC diagnostic push
3618 : #pragma GCC diagnostic ignored "-Wnull-dereference"
3619 : #endif
3620 35 : poDS->aosNonJP2Files = aosNonJP2Files;
3621 : #if defined(__GNUC__)
3622 : #pragma GCC diagnostic pop
3623 : #endif
3624 :
3625 35 : OGRSpatialReference oSRS;
3626 35 : char *pszProjection = nullptr;
3627 70 : if (oSRS.importFromEPSG(nSubDSEPSGCode) == OGRERR_NONE &&
3628 35 : oSRS.exportToWkt(&pszProjection) == OGRERR_NONE)
3629 : {
3630 35 : poDS->SetProjection(pszProjection);
3631 35 : CPLFree(pszProjection);
3632 : }
3633 : else
3634 : {
3635 0 : CPLDebug("SENTINEL2", "Invalid EPSG code %d", nSubDSEPSGCode);
3636 : }
3637 :
3638 35 : double adfGeoTransform[6] = {0};
3639 35 : adfGeoTransform[0] = dfMinX;
3640 35 : adfGeoTransform[1] = nSubDSPrecision;
3641 35 : adfGeoTransform[2] = 0;
3642 35 : adfGeoTransform[3] = dfMaxY;
3643 35 : adfGeoTransform[4] = 0;
3644 35 : adfGeoTransform[5] = -nSubDSPrecision;
3645 35 : poDS->SetGeoTransform(adfGeoTransform);
3646 35 : poDS->GDALDataset::SetMetadataItem("COMPRESSION", "JPEG2000",
3647 : "IMAGE_STRUCTURE");
3648 35 : if (bIsPreview || bIsTCI)
3649 7 : poDS->GDALDataset::SetMetadataItem("INTERLEAVE", "PIXEL",
3650 : "IMAGE_STRUCTURE");
3651 :
3652 35 : int nBits = (bIsPreview || bIsTCI) ? 8 : 0 /* 0 = unknown yet*/;
3653 35 : int nValMax = (bIsPreview || bIsTCI) ? 255 : 0 /* 0 = unknown yet*/;
3654 : const int nBands =
3655 30 : (bIsPreview || bIsTCI)
3656 37 : ? 3
3657 28 : : ((bAlpha) ? 1 : 0) + static_cast<int>(aosBands.size());
3658 35 : const int nAlphaBand = (bIsPreview || bIsTCI || !bAlpha) ? 0 : nBands;
3659 35 : const GDALDataType eDT = (bIsPreview || bIsTCI) ? GDT_Byte : GDT_UInt16;
3660 :
3661 227 : for (int nBand = 1; nBand <= nBands; nBand++)
3662 : {
3663 192 : VRTSourcedRasterBand *poBand = nullptr;
3664 :
3665 192 : if (nBand != nAlphaBand)
3666 : {
3667 190 : poBand = new VRTSourcedRasterBand(
3668 190 : poDS, nBand, eDT, poDS->nRasterXSize, poDS->nRasterYSize);
3669 : }
3670 : else
3671 : {
3672 2 : poBand = new SENTINEL2AlphaBand(
3673 : poDS, nBand, eDT, poDS->nRasterXSize, poDS->nRasterYSize,
3674 2 : nSaturatedVal, nNodataVal);
3675 : }
3676 :
3677 192 : poDS->SetBand(nBand, poBand);
3678 192 : if (nBand == nAlphaBand)
3679 2 : poBand->SetColorInterpretation(GCI_AlphaBand);
3680 :
3681 384 : CPLString osBandName;
3682 192 : if (nBand != nAlphaBand)
3683 : {
3684 190 : osBandName = aosBands[nBand - 1];
3685 190 : SENTINEL2SetBandMetadata(poBand, osBandName);
3686 : }
3687 : else
3688 2 : osBandName = aosBands[0];
3689 :
3690 459 : for (size_t iSrc = 0; iSrc < aosGranuleInfoList.size(); iSrc++)
3691 : {
3692 267 : const SENTINEL2GranuleInfo &oGranuleInfo = aosGranuleInfoList[iSrc];
3693 267 : CPLString osTile;
3694 :
3695 267 : if (bIsSafeCompact && eLevel != SENTINEL2_L2A)
3696 : {
3697 33 : if (bIsTCI)
3698 : {
3699 6 : osTile = oGranuleInfo.osBandPrefixPath + "TCI.jp2";
3700 : }
3701 : else
3702 : {
3703 27 : osTile = oGranuleInfo.osBandPrefixPath + "B";
3704 27 : if (osBandName.size() == 1)
3705 0 : osTile += "0" + osBandName;
3706 27 : else if (osBandName.size() == 3)
3707 2 : osTile += osBandName.substr(1);
3708 : else
3709 25 : osTile += osBandName;
3710 27 : osTile += ".jp2";
3711 : }
3712 : }
3713 : else
3714 : {
3715 468 : osTile = SENTINEL2GetTilename(
3716 234 : oGranuleInfo.osPath, CPLGetFilename(oGranuleInfo.osPath),
3717 : osBandName, osProductURI, bIsPreview,
3718 234 : (eLevel == SENTINEL2_L1C) ? 0 : nSubDSPrecision);
3719 113 : if (bIsSafeCompact && eLevel == SENTINEL2_L2A &&
3720 419 : pType == MSI2Ap && osTile.size() >= 34 &&
3721 306 : osTile.substr(osTile.size() - 18, 3) != "MSK")
3722 : {
3723 60 : osTile.insert(osTile.size() - 34, "L2A_");
3724 : }
3725 234 : if (bIsTCI && osTile.size() >= 14)
3726 : {
3727 0 : osTile.replace(osTile.size() - 11, 3, "TCI");
3728 : }
3729 : }
3730 :
3731 267 : bool bTileFound = false;
3732 267 : if (nValMax == 0)
3733 : {
3734 : /* It is supposed to be 12 bits, but some products have 15 bits
3735 : */
3736 28 : if (SENTINEL2GetTileInfo(osTile, nullptr, nullptr, &nBits))
3737 : {
3738 27 : bTileFound = true;
3739 27 : if (nBits <= 16)
3740 27 : nValMax = (1 << nBits) - 1;
3741 : else
3742 : {
3743 0 : CPLDebug("SENTINEL2", "Unexpected bit depth %d", nBits);
3744 0 : nValMax = 65535;
3745 : }
3746 : }
3747 : }
3748 : else
3749 : {
3750 : VSIStatBufL sStat;
3751 239 : if (VSIStatExL(osTile, &sStat, VSI_STAT_EXISTS_FLAG) == 0)
3752 239 : bTileFound = true;
3753 : }
3754 267 : if (!bTileFound)
3755 : {
3756 1 : CPLError(CE_Warning, CPLE_AppDefined,
3757 : "Tile %s not found on filesystem. Skipping it",
3758 : osTile.c_str());
3759 1 : continue;
3760 : }
3761 :
3762 266 : const int nDstXOff = static_cast<int>(
3763 266 : (oGranuleInfo.dfMinX - dfMinX) / nSubDSPrecision + 0.5);
3764 266 : const int nDstYOff = static_cast<int>(
3765 266 : (dfMaxY - oGranuleInfo.dfMaxY) / nSubDSPrecision + 0.5);
3766 :
3767 266 : if (nBand != nAlphaBand)
3768 : {
3769 263 : poBand->AddSimpleSource(
3770 242 : osTile, (bIsPreview || bIsTCI) ? nBand : 1, 0, 0,
3771 263 : oGranuleInfo.nWidth, oGranuleInfo.nHeight, nDstXOff,
3772 263 : nDstYOff, oGranuleInfo.nWidth, oGranuleInfo.nHeight);
3773 : }
3774 : else
3775 : {
3776 3 : poBand->AddComplexSource(
3777 3 : osTile, 1, 0, 0, oGranuleInfo.nWidth, oGranuleInfo.nHeight,
3778 3 : nDstXOff, nDstYOff, oGranuleInfo.nWidth,
3779 3 : oGranuleInfo.nHeight, nValMax /* offset */, 0 /* scale */);
3780 : }
3781 : }
3782 :
3783 192 : if ((nBits % 8) != 0)
3784 : {
3785 143 : poBand->SetMetadataItem("NBITS", CPLSPrintf("%d", nBits),
3786 143 : "IMAGE_STRUCTURE");
3787 : }
3788 : }
3789 :
3790 35 : return poDS;
3791 : }
3792 :
3793 : /************************************************************************/
3794 : /* OpenL1CTileSubdataset() */
3795 : /************************************************************************/
3796 :
3797 13 : GDALDataset *SENTINEL2Dataset::OpenL1CTileSubdataset(GDALOpenInfo *poOpenInfo)
3798 : {
3799 26 : CPLString osFilename;
3800 13 : CPLAssert(STARTS_WITH_CI(poOpenInfo->pszFilename, "SENTINEL2_L1C_TILE:"));
3801 13 : osFilename = poOpenInfo->pszFilename + strlen("SENTINEL2_L1C_TILE:");
3802 13 : const char *pszPrecision = strrchr(osFilename.c_str(), ':');
3803 13 : if (pszPrecision == nullptr || pszPrecision == osFilename.c_str())
3804 : {
3805 2 : CPLError(CE_Failure, CPLE_AppDefined,
3806 : "Invalid syntax for SENTINEL2_L1C_TILE:");
3807 2 : return nullptr;
3808 : }
3809 11 : const bool bIsPreview = STARTS_WITH_CI(pszPrecision + 1, "PREVIEW");
3810 11 : const int nSubDSPrecision = (bIsPreview) ? 320 : atoi(pszPrecision + 1);
3811 11 : if (!bIsPreview && nSubDSPrecision != 10 && nSubDSPrecision != 20 &&
3812 : nSubDSPrecision != 60)
3813 : {
3814 1 : CPLError(CE_Failure, CPLE_NotSupported, "Unsupported precision: %d",
3815 : nSubDSPrecision);
3816 1 : return nullptr;
3817 : }
3818 10 : osFilename.resize(pszPrecision - osFilename.c_str());
3819 :
3820 20 : std::set<CPLString> oSetBands;
3821 10 : CPLXMLNode *psRootMainMTD = nullptr;
3822 : GDALDataset *poTmpDS =
3823 10 : OpenL1CTile(osFilename, &psRootMainMTD, nSubDSPrecision, &oSetBands);
3824 20 : SENTINEL2_CPLXMLNodeHolder oXMLHolder(psRootMainMTD);
3825 10 : if (poTmpDS == nullptr)
3826 2 : return nullptr;
3827 :
3828 16 : std::vector<CPLString> aosBands;
3829 8 : if (bIsPreview)
3830 : {
3831 2 : aosBands.push_back("04");
3832 2 : aosBands.push_back("03");
3833 2 : aosBands.push_back("02");
3834 : }
3835 : else
3836 : {
3837 21 : for (std::set<CPLString>::const_iterator oIterBandnames =
3838 6 : oSetBands.begin();
3839 48 : oIterBandnames != oSetBands.end(); ++oIterBandnames)
3840 : {
3841 21 : aosBands.push_back(*oIterBandnames);
3842 : }
3843 : /* Put 2=Blue, 3=Green, 4=Band bands in RGB order for conveniency */
3844 14 : if (aosBands.size() >= 3 && aosBands[0] == "02" &&
3845 14 : aosBands[1] == "03" && aosBands[2] == "04")
3846 : {
3847 3 : aosBands[0] = "04";
3848 3 : aosBands[2] = "02";
3849 : }
3850 : }
3851 :
3852 : /* -------------------------------------------------------------------- */
3853 : /* Create dataset. */
3854 : /* -------------------------------------------------------------------- */
3855 :
3856 16 : std::vector<CPLString> aosGranuleList;
3857 8 : aosGranuleList.push_back(osFilename);
3858 :
3859 16 : const int nSaturatedVal = atoi(CSLFetchNameValueDef(
3860 8 : poTmpDS->GetMetadata(), "SPECIAL_VALUE_SATURATED", "-1"));
3861 16 : const int nNodataVal = atoi(CSLFetchNameValueDef(
3862 8 : poTmpDS->GetMetadata(), "SPECIAL_VALUE_NODATA", "-1"));
3863 :
3864 : const bool bAlpha =
3865 8 : CPLTestBool(SENTINEL2GetOption(poOpenInfo, "ALPHA", "FALSE"));
3866 :
3867 16 : std::vector<CPLString> aosNonJP2Files;
3868 8 : SENTINEL2Dataset *poDS = CreateL1CL2ADataset(
3869 : SENTINEL2_L1C, MSI2A,
3870 : false, // bIsSafeCompact
3871 16 : aosGranuleList, std::vector<L1CSafeCompatGranuleDescription>(),
3872 : aosNonJP2Files, nSubDSPrecision, bIsPreview,
3873 : false, // bIsTCI
3874 : -1 /*nSubDSEPSGCode*/, bAlpha, aosBands, nSaturatedVal, nNodataVal,
3875 16 : CPLString());
3876 8 : if (poDS == nullptr)
3877 : {
3878 1 : delete poTmpDS;
3879 1 : return nullptr;
3880 : }
3881 :
3882 : // Transfer metadata
3883 7 : poDS->GDALDataset::SetMetadata(poTmpDS->GetMetadata());
3884 7 : poDS->GDALDataset::SetMetadata(poTmpDS->GetMetadata("xml:SENTINEL2"),
3885 : "xml:SENTINEL2");
3886 :
3887 7 : delete poTmpDS;
3888 :
3889 : /* -------------------------------------------------------------------- */
3890 : /* Add extra band metadata. */
3891 : /* -------------------------------------------------------------------- */
3892 7 : if (psRootMainMTD != nullptr)
3893 6 : poDS->AddL1CL2ABandMetadata(SENTINEL2_L1C, psRootMainMTD, aosBands);
3894 :
3895 : /* -------------------------------------------------------------------- */
3896 : /* Initialize overview information. */
3897 : /* -------------------------------------------------------------------- */
3898 7 : poDS->SetDescription(poOpenInfo->pszFilename);
3899 7 : CPLString osOverviewFile;
3900 7 : if (bIsPreview)
3901 2 : osOverviewFile = CPLSPrintf("%s_PREVIEW.tif.ovr", osFilename.c_str());
3902 : else
3903 : osOverviewFile =
3904 5 : CPLSPrintf("%s_%dm.tif.ovr", osFilename.c_str(), nSubDSPrecision);
3905 7 : poDS->SetMetadataItem("OVERVIEW_FILE", osOverviewFile, "OVERVIEWS");
3906 7 : poDS->oOvManager.Initialize(poDS, ":::VIRTUAL:::");
3907 :
3908 7 : return poDS;
3909 : }
3910 :
3911 : /************************************************************************/
3912 : /* GDALRegister_SENTINEL2() */
3913 : /************************************************************************/
3914 :
3915 1682 : void GDALRegister_SENTINEL2()
3916 : {
3917 1682 : if (GDALGetDriverByName("SENTINEL2") != nullptr)
3918 301 : return;
3919 :
3920 1381 : GDALDriver *poDriver = new GDALDriver();
3921 :
3922 1381 : poDriver->SetDescription("SENTINEL2");
3923 : #ifdef GDAL_DCAP_RASTER
3924 1381 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
3925 : #endif
3926 1381 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
3927 1381 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Sentinel 2");
3928 1381 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
3929 1381 : "drivers/raster/sentinel2.html");
3930 1381 : poDriver->SetMetadataItem(GDAL_DMD_SUBDATASETS, "YES");
3931 :
3932 : #ifdef GDAL_DMD_OPENOPTIONLIST
3933 1381 : poDriver->SetMetadataItem(
3934 : GDAL_DMD_OPENOPTIONLIST,
3935 : "<OpenOptionList>"
3936 : " <Option name='ALPHA' type='boolean' description='Whether to expose "
3937 : "an alpha band' default='NO'/>"
3938 1381 : "</OpenOptionList>");
3939 : #endif
3940 :
3941 1381 : poDriver->pfnOpen = SENTINEL2Dataset::Open;
3942 1381 : poDriver->pfnIdentify = SENTINEL2Dataset::Identify;
3943 :
3944 1381 : GetGDALDriverManager()->RegisterDriver(poDriver);
3945 : }
|