Skip to article frontmatterSkip to article content

Explore EOPF Zarr Sentinel-2 files using GDAL

Explore how to open, visualise and plot Sentinel-2 L1C MSI EOPF Zarr format

Authors
Affiliations
Eurac Research
Eurac Research
Eurac Research
ESA EOPF Zarr Logo

Explore EOPF Zarr Sentinel-2 files using GDAL

🚀 Launch in JupyterHub

Run this notebook interactively with all dependencies pre-installed

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://github.com/EOPF-Sample-Service/GDAL-ZARR-EOPF .

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()
<Figure size 1000x800 with 2 Axes>