Opening Sentinel 2 data with Python [closed]

Question

I'm just trying to learn how to process Sentinel 2 data with Python. All of the tutorials seem to use the data in GeoTIFF format. So I opened the data with SNAP and tried to export as a GeoTIFF but it says that I would need to resample to do that due to the differing pixel sizes. I don't think there should be too much of a problem doing that as long as I'm resampling down the the smallest pixel size (10m) but I'm not sure I'm going about this the right way.

How do I take raw Sentinel (2) data and open it with Python/rasterio?

I am using Level 2A.

Answer 1

GDALs Sentinel2 driver exposes the data as subdatasets, with (sort of) one for each resolution.

To access the subdataset that contains the 10m bands, you could use something like the following:

import rasterio

with rasterio.open('S2A_MSIL2A_20210214T022811_N0214_R046_T51TXL_20210214T044441.zip') as s2a:
    subdatasets = s2a.subdatasets

with rasterio.open(subdatasets[0]) as b10m:
    print(b10m.profile)
    
{'driver': 'SENTINEL2', 'dtype': 'uint16', 'nodata': None, 'width': 10980, 'height': 10980, 'count': 4, 'crs': CRS.from_epsg(32651), 'transform': Affine(10.0, 0.0, 600000.0,
       0.0, -10.0, 5100000.0), 'blockxsize': 128, 'blockysize': 128, 'tiled': True, 'compress': 'jpeg2000'}

Answer 2

If you have obtained your Sentinel-2 data from the ESA in SAFE format (or from other providers such as the DIAS that also distribute it in the same format), you should have several JPEG2000 files (.jp2 extension), one file per band.

You should encounter no issue when reading these files with rasterio, because it uses GDAL for I/O operations and GDAL handles JPEG2000 files (source). There should be no need to convert your files to GeoTIFF.

See this answer to Merging Sentinel 2 RGB bands with rasterio for an example of how to use rasterio to open several bands and merge them into a single file.

Answer 3

Actually, you could use libraries designed to help you to handle satellite products such as EOReader

from eoreader.reader import Reader
from eoreader.bands import GREEN, NDVI, CLOUDS

prod = Reader().open("S2B_MSIL2A_20200114T065229_N0213_R020_T40REQ_20200114T094749.SAFE")

# Load those bands as a dict of xarray.DataArray
band_dict = prod.load([GREEN, NDVI, CLOUDS])
green = band_dict[GREEN]
ndvi = band_dict[NDVI]
clouds = band_dict[CLOUDS]

More details in this notebook.

Disclaimer: I am the maintainer of EOReader.