Accessing Sentinel-2 L2A Scene Classification Layer (SCL) with GDAL

Question

I have a Sentinel-2 L2A and was able to access it with:

In[1]:
image = gdal.Open(s2l2a_image.zip)
subdatasets = image.GetSubdatasets()

Out[1]:
[('SENTINEL2_L2A:/vsizip/s2l2a_image.zip/S2A_MSIL2A_20190620T100031_N0212_R122_T33UXQ_20190620T131720.SAFE/MTD_MSIL2A.xml:10m:EPSG_32633',
  'Bands B2, B3, B4, B8 with 10m resolution, UTM 33N'),
 ('SENTINEL2_L2A:/vsizip/s2l2a_image.zip/S2A_MSIL2A_20190620T100031_N0212_R122_T33UXQ_20190620T131720.SAFE/MTD_MSIL2A.xml:20m:EPSG_32633',
  'Bands B5, B6, B7, B8A, B11, B12 with 20m resolution, UTM 33N'),
 ('SENTINEL2_L2A:/vsizip/s2l2a_image.zip/S2A_MSIL2A_20190620T100031_N0212_R122_T33UXQ_20190620T131720.SAFE/MTD_MSIL2A.xml:60m:EPSG_32633',
  'Bands B1, B9 with 60m resolution, UTM 33N'),
 ('SENTINEL2_L2A:/vsizip/s2l2a_image.zip/S2A_MSIL2A_20190620T100031_N0212_R122_T33UXQ_20190620T131720.SAFE/MTD_MSIL2A.xml:TCI:EPSG_32633',
  'True color image, UTM 33N')]

This is great, however, official GDAL documentation for Sentinel-2 driver (here) talks about L2A-specific layers, such as Scene Classification Layer (SCL) and others. I tried a lot, inspected subdatasets and such, but found no way to access them.

Answer 1

You can still navigate to the SCL layer opening it as a zipped file, e.g. prepending the filename with /vsizip/ and then providing the full path within the archive.

gdal.Open("/vsizip/S2B_MSIL2A_20200328T113319_N0214_R080_T29TNJ_20200328T142123.zip"
+ "/S2B_MSIL2A_20200328T113319_N0214_R080_T29TNJ_20200328T142123.SAFE/"
+ "GRANULE/L2A_T29TNJ_A015979_20200328T113627/"
+ "IMG_DATA/R20m/T29TNJ_20200328T113319_SCL_20m.jp2")

To find out the actual path(s), you can

import zipfile
fname = "S2B_MSIL2A_20200328T113319_N0214_R080_T29TNJ_20200328T142123.zip"
zz = zipfile.ZipFile(fname)

scls = [f.filename 
       for f in zz.filelist 
       if f.filename.find("SCL") >= 0]
# You may want to restrict the filename match to be e.g. SCL_20m or SCL_60m
# scls contains a list with the two scene masks at 20 and 60m:
#['S2B_MSIL2A_20200328T113319_N0214_R080_T29TNJ_20200328T142123.SAFE/GRANULE/L2A_T29TNJ_A015979_20200328T113627/IMG_DATA/R60m/T29TNJ_20200328T113319_SCL_60m.jp2',
 #'S2B_MSIL2A_20200328T113319_N0214_R080_T29TNJ_20200328T142123.SAFE/GRANULE/L2A_T29TNJ_A015979_20200328T113627/IMG_DATA/R20m/T29TNJ_20200328T113319_SCL_20m.jp2']
for scl in scls:
    g = gdal.Open(f"/vsizip/{fname}/{scl}"))
    # Do something with the data

Answer 2

Oddly, I had an outdated GDAL 2.4.3 version where this was an issue. It has been resolved in GDAL 3.1.0.

Consequently, it is possible to access even L2A-specific scenes through gdal.Open() interface.