Topics

Open Santinel 2 archives with Rasterio


lucavlad50@...
 

Hello,
I saw that rasterio support reading Santinel 2 archives directly, but I can not find any example. How should I access the bands inside the archive after I read it. 
For example I have:

with rio.open("satelite_datasets/S2A_MSIL2A_20201020T092031_N0214_R093_T34TGS_20201020T120739.zip"as f:
    print(...)
What should I put instead of "..." to access for example R10 band 4? 

Thank you in advance!


Sean Gillies
 

Hi,

Please see https://rasterio.readthedocs.io/en/latest/quickstart.html#reading-raster-data for a quick tour of data access features.

On Fri, Feb 12, 2021 at 8:01 AM <lucavlad50@...> wrote:
Hello,
I saw that rasterio support reading Santinel 2 archives directly, but I can not find any example. How should I access the bands inside the archive after I read it. 
For example I have:

with rio.open("satelite_datasets/S2A_MSIL2A_20201020T092031_N0214_R093_T34TGS_20201020T120739.zip"as f:
    print(...)
What should I put instead of "..." to access for example R10 band 4? 

Thank you in advance!
_._,_._,

--
Sean Gillies


lucavlad50@...
 

Hello,
Thanks for the reply. I saw that but there is no example of a zip file being opened. I saw this: 
on: https://gdal.org/drivers/raster/sentinel2.html
and was wondering if it is possible to read directly from the zip file without unzipping it. 
Something like: 
with rio.open("satelite_datasets/S2A_MSIL2A_20201020T092031_N0214_R093_T34TGS_20201020T120739.zip"as f:
    print(f.read(4))
where 4 is band B04 in  R10mm.

Sorry if the question seems stupid, I am new to this. Thank you!


Luke
 

GDALs Sentinel2 driver exposes the data as subdatasets, with (sort of) one for each resolution - https://gdal.org/drivers/raster/sentinel2.html#level-2a

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'}


lucavlad50@...
 

That's it, thank you very much!