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