Explore EOPF Zarr Sentinel-2 files using GDAL
Explore how to open, visualise and plot Sentinel-2 L1C MSI EOPF Zarr format
Explore EOPF Zarr Sentinel-2 files using GDAL¶
Run this notebook interactively with all dependencies pre-installed
Table of Contents¶
Introduction¶
This notebook demonstrates accessing and visualizing Sentinel-2 satellite imagery stored in EOPF Zarr format using GDAL and the EOPF GDAL driver developed by the EOPF Sentinels Zarr Samples project. Learn how to access remote Zarr files, extract metadata, and visualize multispectral bands.
The source code of the EOPF GDAL driver is available on GitHub at https://
Import Libraries¶
from osgeo import gdal
import matplotlib.pyplot as plt
import json
# Check if EOPF-Zarr driver is available
gdal.AllRegister()
driver = gdal.GetDriverByName("EOPFZARR")
if driver:
print("✅ EOPF-Zarr driver loaded successfully!")
print(f"Driver description: {driver.GetDescription()}")
else:
print("❌ EOPF-Zarr driver not found")✅ EOPF-Zarr driver loaded successfully!
Driver description: EOPFZARR
Define Product Path¶
The following Sentinel-2 Zarr product is available on an object store, accessible via https:
s2_product_url = "https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:sample-data/tutorial_data/cpm_v253/S2B_MSIL1C_20250113T103309_N0511_R108_T32TLQ_20250113T122458.zarr"Exploring the Zarr Product¶
Product Level Metadata - Using EOPF GDAL Zarr Driver¶
The EOPF GDAL Zarr Driver requires a specific path structure to use the select the correct driver. You need to always prepend `EOPFZARR:"vsicurl/` to Zarr path.
Please check the following code with a minimal working example:
eopf_gdal_path = f'EOPFZARR:"/vsicurl/{s2_product_url}"'
ds_eopf = gdal.Open(eopf_gdal_path, gdal.GA_ReadOnly)
if ds_eopf is None:
raise RuntimeError("Failed to open the Zarr dataset with GDAL.")
print("Dataset opened successfully!")
print("Metadata:", json.dumps(ds_eopf.GetMetadata(), indent=2))Dataset opened successfully!
Metadata: {
"EOPF_PRODUCT": "YES",
"EPSG": "32632",
"proj:epsg": "32632",
"spatial_ref": "PROJCS[\"WGS 84 / UTM zone 32N\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",9],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32632\"]]",
"geo_transform": "300000.000000000000,214.453125000000,0.000000000000,5000040.000000000000,0.000000000000,-214.453125000000",
"utm_easting_min": "300000.00000000",
"utm_easting_max": "409800.00000000",
"utm_northing_min": "4890240.00000000",
"utm_northing_max": "5000040.00000000"
}
/opt/conda/lib/python3.11/site-packages/osgeo/gdal.py:311: FutureWarning: Neither gdal.UseExceptions() nor gdal.DontUseExceptions() has been explicitly called. In GDAL 4.0, exceptions will be enabled by default.
warnings.warn(
Product Level Metadata - Using default Zarr GDAL Driver¶
Similarly, the default GDAL Zarr Driver requires a specific path structure. You need to always prepend ZARR:"vsicurl/ to Zarr path.
The default GDAL Zarr driver is not capable to read the product metadata.
default_gdal_path = f'ZARR:"/vsicurl/{s2_product_url}"'
ds = gdal.Open(default_gdal_path, gdal.GA_ReadOnly)
if ds is None:
raise RuntimeError("Failed to open the Zarr dataset with GDAL.")
print("Dataset opened successfully!")
print("Metadata:", ds.GetMetadata())Dataset opened successfully!
Metadata: {}
Group Level Metadata - Using EOPF GDAL Zarr Driver¶
The Zarr groups are identifies as subdatasets in GDAL.
The EOPF GDAL Zarr driver is capable to read and parse more metadata fields compared to the default Zarr driver.
sub_ds_eopf = ds_eopf.GetSubDatasets()
for idx, s in enumerate(sub_ds_eopf[:5]): # Show first 5 subdatasets
print(f"{idx}: {s[0]}")0: EOPFZARR:"/vsicurl/https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:sample-data/tutorial_data/cpm_v253/S2B_MSIL1C_20250113T103309_N0511_R108_T32TLQ_20250113T122458.zarr/conditions/geometry/angle"
1: EOPFZARR:"/vsicurl/https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:sample-data/tutorial_data/cpm_v253/S2B_MSIL1C_20250113T103309_N0511_R108_T32TLQ_20250113T122458.zarr/conditions/geometry/band"
2: EOPFZARR:"/vsicurl/https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:sample-data/tutorial_data/cpm_v253/S2B_MSIL1C_20250113T103309_N0511_R108_T32TLQ_20250113T122458.zarr/conditions/geometry/detector"
3: EOPFZARR:"/vsicurl/https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:sample-data/tutorial_data/cpm_v253/S2B_MSIL1C_20250113T103309_N0511_R108_T32TLQ_20250113T122458.zarr/conditions/geometry/x"
4: EOPFZARR:"/vsicurl/https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:sample-data/tutorial_data/cpm_v253/S2B_MSIL1C_20250113T103309_N0511_R108_T32TLQ_20250113T122458.zarr/conditions/geometry/y"
target_subds_eopf = gdal.Open(sub_ds_eopf[80][0], gdal.GA_ReadOnly)
print("Band metadata:", json.dumps(target_subds_eopf.GetMetadata(), indent=2))Band metadata: {
"EOPF_PRODUCT": "YES",
"EPSG": "32632",
"proj:wkt2": "PROJCS[\"WGS 84 / UTM zone 32N\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",9],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32632\"]]",
"spatial_ref": "PROJCS[\"WGS 84 / UTM zone 32N\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",9],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32632\"]]",
"geo_transform": "300000.000000000000,60.000000000000,0.000000000000,5000040.000000000000,0.000000000000,-60.000000000000",
"utm_easting_min": "300000.00000000",
"utm_easting_max": "409800.00000000",
"utm_northing_min": "4890240.00000000",
"utm_northing_max": "5000040.00000000",
"_eopf_attrs": "{ \"add_offset\": -0.1, \"coordinates\": [ \"x\", \"y\" ], \"dimensions\": [ \"y\", \"x\" ], \"fill_value\": 0, \"scale_factor\": 0.0001, \"units\": \"digital_counts\" }",
"dtype": "<u2",
"fill_value": "0",
"long_name": "TOA reflectance from MSI acquisition at spectral band b01 442.3 nm",
"valid_max": "65535",
"valid_min": "1"
}
ERROR 3: Cannot open file '/vsicurl/https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:sample-data/tutorial_data/cpm_v253/S2B_MSIL1C_20250113T103309_N0511_R108_T32TLQ_20250113T122458.zarr/measurements/reflectance/r60m/b01/.zmetadata'
ERROR 3: Load json file /vsicurl/https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:sample-data/tutorial_data/cpm_v253/S2B_MSIL1C_20250113T103309_N0511_R108_T32TLQ_20250113T122458.zarr/measurements/reflectance/r60m/b01/.zmetadata failed
Warning 1: Unable to save auxiliary information in ZARR:"/vsicurl/https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:sample-data/tutorial_data/cpm_v253/S2B_MSIL1C_20250113T103309_N0511_R108_T32TLQ_20250113T122458.zarr/measurements/reflectance/r60m/b01".aux.xml.
Group Level Metadata - Using default GDAL Zarr Driver¶
sub_ds = ds.GetSubDatasets()
for idx, s in enumerate(sub_ds[:5]): # Show first 5 subdatasets
print(f"{idx}: {s[0]}")0: ZARR:"/vsicurl/https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:sample-data/tutorial_data/cpm_v253/S2B_MSIL1C_20250113T103309_N0511_R108_T32TLQ_20250113T122458.zarr":/conditions/geometry/angle
1: ZARR:"/vsicurl/https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:sample-data/tutorial_data/cpm_v253/S2B_MSIL1C_20250113T103309_N0511_R108_T32TLQ_20250113T122458.zarr":/conditions/geometry/band
2: ZARR:"/vsicurl/https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:sample-data/tutorial_data/cpm_v253/S2B_MSIL1C_20250113T103309_N0511_R108_T32TLQ_20250113T122458.zarr":/conditions/geometry/detector
3: ZARR:"/vsicurl/https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:sample-data/tutorial_data/cpm_v253/S2B_MSIL1C_20250113T103309_N0511_R108_T32TLQ_20250113T122458.zarr":/conditions/geometry/x
4: ZARR:"/vsicurl/https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:sample-data/tutorial_data/cpm_v253/S2B_MSIL1C_20250113T103309_N0511_R108_T32TLQ_20250113T122458.zarr":/conditions/geometry/y
target_subds = gdal.Open(sub_ds[80][0], gdal.GA_ReadOnly)
print("Band metadata:", json.dumps(target_subds.GetMetadata(), indent=2))Band metadata: {
"_eopf_attrs": "{ \"add_offset\": -0.1, \"coordinates\": [ \"x\", \"y\" ], \"dimensions\": [ \"y\", \"x\" ], \"fill_value\": 0, \"scale_factor\": 0.0001, \"units\": \"digital_counts\" }",
"dtype": "<u2",
"fill_value": "0",
"long_name": "TOA reflectance from MSI acquisition at spectral band b01 442.3 nm",
"proj:bbox": "{300000,4890240,409800,5000040}",
"proj:epsg": "32632",
"proj:shape": "{1830,1830}",
"proj:transform": "{60,0,300000,0,-60,5000040,0,0,1}",
"proj:wkt2": "PROJCS[\"WGS 84 / UTM zone 32N\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",9],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"32632\"]]",
"valid_max": "65535",
"valid_min": "1"
}
Data Visualization¶
The way we access the data is the same for both GDAL drivers
band = target_subds_eopf.GetRasterBand(1)
data = band.ReadAsArray()
print(f"Data dimensions: {data.shape}")Data dimensions: (1830, 1830)
plt.figure(figsize=(10, 8))
plt.imshow(data, cmap="viridis")
plt.colorbar(label="Reflectance Value")
plt.title("Sentinel-2 Band B01 (60m Resolution)")
plt.axis("off")
plt.show()
