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