LCOV - code coverage report
Current view: top level - frmts/sentinel2 - sentinel2dataset.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 1724 1912 90.2 %
Date: 2025-01-18 12:42:00 Functions: 43 43 100.0 %

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

Generated by: LCOV version 1.14