LCOV - code coverage report
Current view: top level - frmts/sentinel2 - sentinel2dataset.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 1489 1601 93.0 %
Date: 2017-12-13 Functions: 52 53 98.1 %

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

Generated by: LCOV version 1.10