Sentinel-2 L1C MSI Zarr Product Exploration

Table of Contents¶

Introduction
Import modules
Open the product
Overview of the product content
Plot a RGB image
Explore product geolocation
A. Interactive map
B. Non-interactive map
Compute radiances
Visualization

Introduction¶

This notebook demonstrates how to open, explore, and visualize Sentinel-2 Level 1C MSI products stored in EOPF Zarr format, including plotting RGB images and analyzing geolocation data.

Import modules¶

from eopf.product.eo_variable import EOVariable
import cartopy.crs as ccrs
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr

Open the product¶

remote_product_path = "https://objectstore.eodc.eu:2222/e05ab01a9d56408d82ac32d69a5aae2a:sample-data/tutorial_data/cpm_v253/S2B_MSIL1C_20250113T103309_N0511_R108_T32TLQ_20250113T122458.zarr"
dt = xr.open_datatree(remote_product_path, engine="zarr", chunks={})

Overview of the product content¶

dt["measurements/reflectance/r60m"]["b01"].plot()
# Alternative notation
# dt.measurements.reflectance.r60m.b01.plot()

The original data was stored in uint16, with scale and offset values to be applied.
In this case, the Zarr contains already reflectance values, which means that the scale and offset are already applied to the data.
We can get some info about the encoding with the following code:

dt["measurements/reflectance/r60m"]["b01"].encoding

{'chunks': (305, 305),
 'preferred_chunks': {'y': 305, 'x': 305},
 'compressor': Blosc(cname='zstd', clevel=3, shuffle=BITSHUFFLE, blocksize=0),
 'filters': None,
 '_FillValue': 0,
 'scale_factor': 0.0001,
 'add_offset': -0.1,
 'dtype': dtype('uint16')}

If we are interested in some band statistics, they can easily be computed using the xarray methods.
Remember that calling .compute() triggers the data download and processing of the result.

min, max, mean = (
    dt["measurements/reflectance/r60m"]["b01"].data.max(),
    dt["measurements/reflectance/r60m"]["b01"].data.min(),
    dt["measurements/reflectance/r60m"]["b01"].data.mean(),
)
min.compute(), max.compute(), mean.compute()

(1.6631, 0.0746, 0.27504428200304576)

Plot a RGB image¶

We can compose a new object for visualization, taking only the red (b04), green (b03) and blue (b02) bands.

Decimation is also applied, which reduces the data we consider for visualization, speeding up the process. Decimation is performed using the .isel xarray method.

.isel(x=slice(None, None, 10), y=slice(None, None, 10))

Finally, we rescale the data for a more meaningful visualization.

rgb_band_paths = (
    "measurements/reflectance/r10m/b04",
    "measurements/reflectance/r10m/b03",
    "measurements/reflectance/r10m/b02",
)
# Concatenate the selected bands in a single xarray DataArray, applying decimation with a factor of 10
concat = xr.concat(
    [
        dt[str(p)].isel(x=slice(None, None, 10), y=slice(None, None, 10))
        for p in rgb_band_paths
    ],
    dim="band",
)

# Rescale the data for a more meaningful visualization
concat = (concat / 0.4).clip(0, 1)

ax = concat.plot.imshow()
ax.axes.set_aspect("equal")
plt.gca().invert_yaxis()

Explore product geolocation¶

The Zarr product metadata contains information about the product geolocation we can explore in various ways.

Interactive map¶

The following code snippet shows an interactive map with the tile’s footprint

import geopandas as gpd

# Create a GeoDataFrame from the feature
gdf = gpd.GeoDataFrame.from_features(
    [
        {
            "type": "Feature",
            "geometry": dt.attrs["stac_discovery"][
                "geometry"
            ],  # Use the actual geometry, not the bbox
            "properties": dt.attrs["stac_discovery"][
                "properties"
            ],  # Include all properties
        }
    ],
    crs="EPSG:4326",  # Set the CRS explicitly (adjust EPSG code if needed)
)

gdf.explore()

Non-interactive map¶

The following code snippet shows the location of the tile on a global map.

# Reproject to a projected CRS
gdf_projected = gdf.to_crs(epsg=3857)
centroids = gdf_projected.centroid

fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Robinson())

# Make the map global rather than have it zoom in to the extents of any plotted data
ax.set_global()

# ax.stock_img()
ax.coastlines()

# Transform back to geographic CRS for plotting
gdf_centroids_geographic = centroids.to_crs(epsg=4326)

ax.plot(
    gdf_centroids_geographic.x,
    gdf_centroids_geographic.y,
    "ro",
    transform=ccrs.PlateCarree(),
)

plt.show()

Compute radiances¶

From: https://s2.pages.eopf.copernicus.eu/pdfs-adfs/MSI/L1/PDFS_S2_MSI_L1C.myst.html#product-overview

$\text{radiance} = \text{reflectance} \times \cos\left(\text{radians}(\text{SunZenithAngle})\right) \times \text{solarIrradiance} \times \frac{U}{\pi}$ $ ]

U: float = dt.attrs["other_metadata"][
    "reflectance_correction_factor_from_the_Sun-Earth_distance_variation_computed_using_the_acquisition_date"
]
U

1.03411047670495

Be carefull, Sun Zenith Angle is expressed on the angles grid (5km), it needs to be reprojected on the 10m grid for computing radiances cosinus is applied now because we can not interpolate angles using a linear interpolation (discontinuity at 0°).
On the other hand, cosines can be interpolated:

cos_sza_5km: EOVariable = np.cos(np.deg2rad(dt["conditions/geometry/sun_angles"])).sel(
    angle="zenith"
)
cos_sza_5km

# We will convert reflectances from band b02 to radiances
BAND: int = 2  # b02

# Band - 1 because Python list index starts at 0
solarIrradiance: float = np.float64(
    dt.attrs["stac_discovery"]["properties"]["bands"][BAND - 1]["solar_illumination"]
)

reflectance_b02_10m: EOVariable = dt["measurements/reflectance/r10m"]["b02"]
reflectance_b02_10m

Interpolate sza on the angles grid to the 10m grid, using the xarray .interp_like() method, which interpolates our input data to match the given target.

cos_sza_10m: xr.DataArray = cos_sza_5km.interp_like(reflectance_b02_10m)
cos_sza_10m.compute()

For simplicity, radiance computation assume that reflectances equal numerical counts

radiance = reflectance_b02_10m * cos_sza_10m * solarIrradiance * U / np.pi
radiance

Visualization¶

Visualize Dask computational graph to compute radiance values. For more info, see the Dask docs: https://docs.dask.org/en/stable/graphviz.html#visualize-task-graphs

radiance.data.visualize()

Visualize the resulting radiance values

radiance.plot()

Jupyter Notebooks

SAFE vs Zarr Metadata

Jupyter Notebooks

SENTINEL-1 L1 GRD Product Format Prototype