Computing Leaf Area Index Using Sentinel-2 Data

Introduction

This notebook computes Leaf Area Index (LAI) 🌿 from Sentinel-2 Level-2A satellite data. Input data is read directly from cloud-native Zarr archives published by the Earth Observation Processing Framework (EOPF), so no SAFE archive downloads are required, saving what can take several minutes per Sentinel-2 L2A product. It processes each Sentinel-2 tile in full — consistent with the behaviour of the SNAP BiophysicalOp processor. Each STAC item returned by the query is processed independently and saved as a separate Zarr store under OUTPUT_DIR.

What we will learn

🗂️ The difference between using EOPF Zarr datasets and SAFE archives for processing Sentinel-2 data.
🔎 How to integrate EOPF STAC queries and results into a processing pipeline workflow using a neural network inferrence framework.
🛰️ How to extract sensor metadata, such as sun and view zenith and azimuth angles, for computations.
📚 How to combine extracted metadata with observation data for neural network inference.
⚙️ How to run stacked tiles through neural netowrk inference.
📘 How to save results to the Zarr format.
🗺️ How to visualize results using an interactive map.
⚙️⚙️ How to batch process a full growing season.

Prerequisites

⚠️ This is a more advanced walkthrough, so you should have experience processing data cubes using Python with Xarray and NumPy.

📜 Background

Traditional LAI processing with SNAP’s BiophysicalOp begins with downloading an entire Sentinel-2 L2A SAFE archive — typically 600 MB to 1.2 GB per granule — before a single pixel can be processed. On a 100 Mbps connection that equates to roughly 1–2 minutes of idle waiting per tile; on the slower connections common in research institutions or field deployments it can exceed 10 minutes. This notebook eliminates that bottleneck entirely. By opening EOPF products as lazy xarray DataTree objects backed by cloud-hosted Zarr stores, only the compressed chunks are fetched from object storage.

Scaling to a time series or multi-region study exposes an even larger gap. SNAP requires loading each product in the GUI, configuring the BiophysicalOp operator, and exporting results one at a time — a workflow that does not parallelise or automate easily. With the EOPF STAC API, expanding from a single acquisition to a full growing-season time series is a matter of widening the date window and wrapping the processing step in a loop. A European site typically receives 20–30 cloud-free Sentinel-2 passes between April and October; processing all of them in SNAP demands downloading 12–36 GB and chaining 20–30 manual processing runs — a multi-hour exercise even on fast hardware. The same campaign with this notebook streams only the required chunks per scene and processes each tile in approximately 60-90 seconds (with GPU acceleration for running the LAI neural network), compressing the entire time series into an unattended run measured in minutes rather than hours.

Beyond speed, the workflow is designed for reproducibility and frictionless deployment. SNAP carries a roughly 1 GB installation footprint, requires a Java runtime, and exposes its automation interface through the gpt command-line operator — a non-trivial barrier for containerised cloud workers, CI pipelines, or shared compute environments. This notebook’s entire software stack is declared in a single pixi.toml lockfile and installed with one command on any platform. Results are written to CF-compliant Zarr stores — a format that can be streamed directly by Python or a web mapping service without further conversion. Together, these properties make the workflow naturally suited to reproducible research and cloud-native execution at scale.

🌿 What is LAI?

Leaf Area Index (LAI) is a dimensionless biophysical quantity defined as half the total developed area of photosynthetically active elements of the vegetation per unit horizontal ground surface area (m² m⁻²). It is one of the Essential Climate Variables (ECVs) recognised by the Global Climate Observing System (GCOS) and underpins a wide range of Earth observation applications:

Agriculture — crop monitoring, yield forecasting, and irrigation scheduling
Ecology — forest carbon estimation, phenology studies, and habitat mapping
Climate science — land surface model parameterisation and validation
Environmental monitoring — deforestation detection and drought assessment

Typical LAI values range from 0 (bare soil) to roughly 8–10 m² m⁻² for dense closed-canopy forests or intensively managed crops at peak growth.

🛰️ Origins: The SNAP Sentinel-2 Optical Toolbox

The algorithm implemented here is a Python re-implementation of the BiophysicalOp processor from ESA’s open-source SNAP (Sentinel Application Platform), specifically the Optical Toolbox. The underlying science is described in the S2ToolBox Algorithm Theoretical Basis Document (ATBD):

Weiss, M., Baret, F., & Jay, S. (2020). S2ToolBox Level 2 products: LAI, FAPAR, FCOVER – Version 2.1.

The method employs a shallow two-layer feed-forward neural network trained over a large look-up table generated by the PROSAIL radiative-transfer model (PROSPECT for leaf optics + SAIL for canopy architecture). ESA distributes sensor-specific network weight files for Sentinel-2A, -2B, and -2C at both 10 m and 20 m spatial resolutions; this notebook automatically selects the correct weights from the detected satellite ID.

Differences from the Reference Java Implementation

Aspect	SNAP (Java)	This Notebook
Input data	Local SAFE files	EOPF Zarr via object storage
NN runtime	Java native arrays	Modular MAX engine (compiled graph)
Output format	BEAM-DIMAP	CF-compliant Zarr + Cloud-Optimised GeoTIFF

⚙️ Processing Pipeline

The workflow follows the original SNAP BiophysicalOp processing chain:

STAC Discovery — query the EOPF catalog for scenes intersecting the AOI and date window
Angle Interpolation — bilinearly interpolate the sparse 23 × 23 angle grids to full image resolution, per detector module
Input Assembly — stack 8 reflectance bands with 3 trigonometric angle terms (11 features total)
Input Validation — flag pixels outside the NN training domain (per-band min/max check + convex-hull check)
NN Inference — two-layer tanh network executed by the Modular MAX engine
Output Validation — flag and clip outputs to physically plausible bounds
Export — write LAI and quality flags to Zarr and/or GeoTIFF

The following modules cover each stage of the processing pipeline:

Module	Role
`pystac.client`	Query the EOPF STAC catalog to discover available Sentinel-2 L2A products by spatial extent, time range, and cloud cover
`xarray`	Lazily open and navigate the hierarchical EOPF Zarr stores as `DataTree` objects — only the requested array chunks are fetched from object storage
`numpy`	Efficient array operations for angle computation, input stacking, masking, and result post-processing
`biophysical.angles`	Detector-footprint-aware interpolation and aggregation of per-band view zenith/azimuth angle grids
`biophysical.io`	Save LAI values and quality flags to CF-compliant Zarr or Cloud-Optimised GeoTIFF
`biophysical.models`	Enum types defining supported sensors (`S2A`, `S2B`, `S2C` at 10 m / 20 m) and biophysical variables (`LAI`)
`biophysical.processing`	Full BiophysicalOp pipeline: input validation → MAX-accelerated NN inference → output validation
`biophysical.utils`	Bilinear interpolation helpers and Sentinel-2 product name parser

Note: This notebook takes advantage of the biophysical package included with it due to the complexity of the above algorithm and its accompanying pre- and post-processing requirements.

🛠️ Setup

Environment Setup with uv

This notebook uses Modular’s MAX Inference Engine for neural network inference because of its ability to run the same code on either the CPU or any of several different GPUs from NVIDIA, AMD, and soon Apple, if one is available. The biophysicalop package provides the LAI retrieval implementation.

There are two ways to add it to your project:

Option A — add as a `uv` source dependency (recommended and already completed for this project)

Modify your project’s pyproject.toml to include biophysicalop as a git source and configure access to Modular’s nightly wheel index:

[project]
dependencies = [
    "biophysicalop",
    # ...
]

[tool.uv.sources]
biophysicalop = { git = "https://github.com/Geomatys/biophysicalop.git" }

Then sync your environment:

uv sync

Option B — `pip` install in a notebook cell

For environments where modular version 26.2 is already available, install the package directly from within a notebook cell:

%pip install git+https://github.com/Geomatys/biophysicalop.git

⚙️ Processing Pipeline

import sys
sys.path.insert(0, ".")  # Add current directory to sys.path for local imports

import shutil
from pathlib import Path

import warnings , os, s3fs

import branca.colormap as bcm
import folium
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import zarr
from branca.element import Element
from IPython.display import display
from pyproj import Transformer
from pystac.client import Client

from biophysicalop.angles import (
    add_detector_dimension,
    select_angle_by_detector_mask,
    compute_view_angle_mean,
    downsample_detector_masks,
    )
from biophysicalop.io import save_zarr, epsg_from_tile_id
from biophysicalop.models import BiophysicalVariables, SatelliteSensors
from biophysicalop.processing import BiophysicalOpProcess, process
from biophysicalop.utils import bilinear_interpolate_xy, parse_sentinel2_product_name

from save_zarr_utils import (
    prep_zarr
)

We prepare our credentials for S3 access:

warnings.filterwarnings("ignore")
try:
    bucket = os.environ["BUCKET_NAME"]
    access = os.environ["ACCESS_KEY"]
    secret = os.environ["SECRET_KEY"]
    bucket_endpoint = os.environ["BUCKET_ENDPOINT"]
except KeyError as e:
    raise RuntimeError(
        f"Missing required environment variable: {e}\n"
        "Please set BUCKET_NAME, ACCESS_KEY, and SECRET_KEY"
    )

# S3 filesystem
fs = s3fs.S3FileSystem(
    key=access,
    secret=secret,
    client_kwargs={"endpoint_url": bucket_endpoint},
)

Setup the constants to be used in the processing

The parameters below control what data is fetched and at what resolution. The most important choices are:

Time and quality

START_DATE_TIME / END_DATE_TIME — ISO-8601 UTC bounds for the temporal search.
CLOUD_COVER_MAX — Maximum acceptable cloud cover as a percentage (i.e. 5.0 = 5.0%)
OUTPUT_DIR — Directory where per-tile Zarr stores are written.

Spatial resolution

SPATIAL_RESOLUTION = 20 (recommended) — uses 8 spectral bands (B03–B12, excluding B08) and the 20 m NN weights. Produces a 5 490 × 5 490 pixel grid for a standard 100 km × 100 km Sentinel-2 tile and matches the reference SNAP BiophysicalOp output.
SPATIAL_RESOLUTION = 10 — uses only 3 bands (B03, B04, B08) with a separate, lighter 10 m network. Produces a 10 980 × 10 980 grid at the cost of reduced spectral diversity.

Band and geometry lookups

SENTINEL2_BAND_GROUPS — maps each band name to its native resolution group inside the EOPF Zarr hierarchy (r10m, r20m, or r60m).
SENTINEL2_GEOMETRY_BAND_INDICES — maps spectral band names to the integer index used inside conditions/geometry.
SENTINEL2_LAI_BANDS_20M — the ordered list of 8 reflectance bands that form the first 8 features of the 20 m neural network input vector.

# User-defined parameters
# --------------------------------------------------------------------------------
# [minx, miny, maxx, maxy]
AOI_BBOX = (3.1360, 43.4477, 4.1134, 44.1579)  # an inset of tile 31TEJ, covering Montpellier, France

# Temporal range for Sentinel-2 data filtering (in ISO 8601 format)
# These dates are chosen to cover the growing season in the Northern Hemisphere,
# which is typically from April to October.
START_DATE_TIME = "2025-04-01T00:00:00Z"
END_DATE_TIME = "2025-10-31T23:59:59.999999Z"

CLOUD_COVER_MAX = 5.0  # maximum acceptable cloud cover as a percentage (i.e. 5.0 = 5.0%)

OUTPUT_DIR = Path("lai_tiles")  # per-tile zarr output directory
# --------------------------------------------------------------------------------


TEMPORAL_RANGE = f"{START_DATE_TIME}/{END_DATE_TIME}"
SPATIAL_RESOLUTION = 20  # in meters, must be one of {10, 20} for Sentinel-2
SENTINEL2_BAND_GROUPS = {
    "b01": "r60m",
    "b02": "r10m",
    "b03": "r10m",
    "b04": "r10m",
    "b05": "r20m",
    "b06": "r20m",
    "b07": "r20m",
    "b08": "r10m",
    "b8a": "r20m",
    "b09": "r60m",
    "b10": "r60m",
    "b11": "r20m",
    "b12": "r20m",
}
SENTINEL2_GEOMETRY_BAND_INDICES = {
    "b01": "b00",
    "b02": "b01",
    "b03": "b02",
    "b04": "b03",
    "b05": "b04",
    "b06": "b05",
    "b07": "b06",
    "b08": "b07",
    "b8a": "b08",
    "b09": "b09",
    "b10": "b10",
    "b11": "b11",
    "b12": "b12",
}
SENTINEL2_LAI_BANDS_10M = (
    "b03",
    "b04",
    "b08",
)
SENTINEL2_LAI_BANDS_20M = (
    "b03",
    "b04",
    "b05",
    "b06",
    "b07",
    "b8a",
    "b11",
    "b12",
    )

# Constants for EOPF
SAT_COLLECTION = "sentinel-2-l2a"
EOPF_STAC_API_ENDPOINT = "https://stac.core.eopf.eodc.eu/"

🗂️ The Difference Between Using EOPF Zarr Datasets and SAFE Archives

The Earth Observation Processing Framework (EOPF) repackages ESA Copernicus products as Zarr archives on object storage. Zarr is a cloud-native, chunk-based array format: instead of downloading a multi-gigabyte SAFE archive, only the specific array chunks needed are fetched over HTTP, enabling efficient sub-region and sub-band access. As stated in the background, traditional LAI processing with the SNAP Optical Toolbox BiophysicalOp processor begins with downloading an entire Sentinel-2 L2A SAFE archive, usually requiring 600 MB to 1.2 GB to be downloaded per granule.

🔎 Integrating EOPF STAC Queries and Results into a processing Pipeline

Search EOPF STAC for Sentinel-2 Images

Use the user defined time range from the above constants to query for Sentinel-2 Products.

Product discovery follows the open SpatioTemporal Asset Catalog (STAC) specification. Each STAC item represents one Sentinel-2 granule and includes a "product" asset pointing to the root URL of its Zarr store.

The search below returns Sentinel-2 L2A granules that: - Spatially intersect the AOI_BBOX bounding box - Temporally overlap the TEMPORAL_RANGE window - Have less than CLOUD_COVER_MAX% cloud cover (using the ESA-estimated eo:cloud_cover metadata field)

eopf_catalog = Client.open(url=EOPF_STAC_API_ENDPOINT)
s2_results = eopf_catalog.search(
    collections=[SAT_COLLECTION],
    bbox=AOI_BBOX,
    datetime=TEMPORAL_RANGE,
    query={"eo:cloud_cover": {"lt": CLOUD_COVER_MAX}}  # Filter for low cloud cover
)
items = tuple(s2_results.items())
print(f"Total Sentinel-2 L2A items found: {len(items)}")
for item in items:
    print(f"Item ID: {item.id}")
    print(f"  Geometry: {item.geometry}")
    print(f"  Bounding Box: {item.bbox}")
    print(f"  Datetime: {item.datetime}")
    print(f"  Cloud Cover: {item.properties.get('eo:cloud_cover', 'N/A')}")
    product = item.assets["product"]
    print(f"  Product HREF: {product.href}")  # Use product href to get zarr url

Total Sentinel-2 L2A items found: 15
Item ID: S2A_MSIL2A_20251030T104211_N0511_R008_T31TEJ_20251030T144716
  Geometry: {'type': 'Polygon', 'coordinates': [[[2.999749474120944, 44.25341758107761], [4.374940058749246, 44.24514241182411], [4.352494032483418, 43.256794294201605], [2.999753565511865, 43.26479037737325], [2.999749474120944, 44.25341758107761]]]}
  Bounding Box: [2.999749474120944, 43.256794294201605, 4.374940058749246, 44.25341758107761]
  Datetime: 2025-10-30 10:42:11.024000+00:00
  Cloud Cover: 0.713588
  Product HREF: https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/cpm_v270/S2A_MSIL2A_20251030T104211_N0511_R008_T31TEJ_20251030T144716.zarr
Item ID: S2B_MSIL2A_20251030T103049_N0511_R108_T31TEJ_20251030T130433
  Geometry: {'type': 'Polygon', 'coordinates': [[[4.005694176135717, 43.25884423693857], [4.006666737240778, 43.2618584522954], [4.054214641392845, 43.408714965095214], [4.102001492637269, 43.55556092722918], [4.15009382059985, 43.70238356209406], [4.198619251512779, 43.84908934443938], [4.247012902182092, 43.995914013230525], [4.296322214027034, 44.14255029482407], [4.330147398004033, 44.24541195035435], [4.374940058749246, 44.24514241182411], [4.352494032483418, 43.256794294201605], [4.005694176135717, 43.25884423693857]]]}
  Bounding Box: [4.005694176135717, 43.256794294201605, 4.374940058749246, 44.24541195035435]
  Datetime: 2025-10-30 10:30:49.024000+00:00
  Cloud Cover: 0.082116
  Product HREF: https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/cpm_v270/S2B_MSIL2A_20251030T103049_N0511_R108_T31TEJ_20251030T130433.zarr
Item ID: S2C_MSIL2A_20251028T104151_N0511_R008_T31TEJ_20251028T145122
  Geometry: {'type': 'Polygon', 'coordinates': [[[2.999749474120944, 44.25341758107761], [4.374940058749246, 44.24514241182411], [4.352494032483418, 43.256794294201605], [2.999753565511865, 43.26479037737325], [2.999749474120944, 44.25341758107761]]]}
  Bounding Box: [2.999749474120944, 43.256794294201605, 4.374940058749246, 44.25341758107761]
  Datetime: 2025-10-28 10:41:51.025000+00:00
  Cloud Cover: 1.649756
  Product HREF: https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/cpm_v270/S2C_MSIL2A_20251028T104151_N0511_R008_T31TEJ_20251028T145122.zarr
Item ID: S2C_MSIL2A_20251018T104051_N0511_R008_T31TEJ_20251018T205314
  Geometry: {'type': 'Polygon', 'coordinates': [[[2.999749735219841, 44.190326696329606], [3.080339941079615, 44.1770058700864], [3.080170636302186, 44.17652324772636], [3.080263323696628, 44.176506816370726], [3.07996368494882, 44.175653051885135], [3.35386759094171, 44.126916725834604], [3.354053308066213, 44.12743560297753], [3.354333588729886, 44.12738583299684], [3.354514055378814, 44.12788843271424], [3.354581669773897, 44.12787644622808], [3.354590676309644, 44.12790162287181], [3.358243193930242, 44.127227341720285], [3.372804642699915, 44.12464592981835], [3.372768770285561, 44.12454581534072], [3.644154551226836, 44.074446038723785], [3.644301827954194, 44.074852523185726], [3.646290132854245, 44.07448576793084], [3.661629478014961, 44.071509882675], [3.661172835930347, 44.070251005517264], [3.931814146994537, 44.017803429006285], [3.931932284678499, 44.01812510406819], [3.931939135312465, 44.01812377567094], [3.932226810709366, 44.01890469249662], [3.948922655063657, 44.01548321726738], [3.948769020412453, 44.01506485326411], [4.237104469035296, 43.95606006746774], [4.237037982837612, 43.95588185462314], [4.238085647405527, 43.955659811266806], [4.237886297413936, 43.95512701199173], [4.238034088381697, 43.9550957154162], [4.237925609151044, 43.95480486726842], [4.367723347430212, 43.92737467141782], [4.352494032483418, 43.256794294201605], [2.999753565511865, 43.26479037737325], [2.999749735219841, 44.190326696329606]]]}
  Bounding Box: [2.999749735219841, 43.256794294201605, 4.367723347430212, 44.190326696329606]
  Datetime: 2025-10-18 10:40:51.025000+00:00
  Cloud Cover: 0.409785
  Product HREF: https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/cpm_v270/S2C_MSIL2A_20251018T104051_N0511_R008_T31TEJ_20251018T205314.zarr
Item ID: S2C_MSIL2A_20251018T104051_N0511_R008_T31TEJ_20251018T200113
  Geometry: {'type': 'Polygon', 'coordinates': [[[4.363694039141706, 43.74995532711596], [4.157368162595406, 43.79346341981967], [4.157638319129653, 43.7941901876288], [4.157591188638929, 43.79419986379704], [4.157611894201473, 43.79425557226157], [4.157562147972216, 43.7942657913887], [4.157741992589658, 43.79474893476958], [3.888329331098034, 43.85016310759728], [3.888269674712022, 43.85000068524352], [3.888227552877011, 43.850009342874394], [3.887905751285096, 43.8491322269987], [3.583984502544144, 43.90807962528523], [3.584264495792644, 43.90885361671259], [3.584252941555609, 43.90885575324023], [3.584476459424989, 43.90947277964113], [3.313848366818832, 43.95946280377139], [3.313799872198968, 43.95932740516691], [3.313762211068698, 43.95933436227994], [3.313470629605916, 43.95851939638556], [3.006287053779157, 44.013333526199695], [3.006588778942535, 44.01419031849532], [2.999750459510941, 44.015311915785524], [2.999749474120944, 44.25341758107761], [4.374940058749246, 44.24514241182411], [4.363694039141706, 43.74995532711596]]]}
  Bounding Box: [2.999749474120944, 43.74995532711596, 4.374940058749246, 44.25341758107761]
  Datetime: 2025-10-18 10:40:51.025000+00:00
  Cloud Cover: 0.65563
  Product HREF: https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/cpm_v270/S2C_MSIL2A_20251018T104051_N0511_R008_T31TEJ_20251018T200113.zarr
Item ID: S2A_MSIL2A_20251017T103041_N0511_R108_T31TEJ_20251017T155021
  Geometry: {'type': 'Polygon', 'coordinates': [[[4.009598665597077, 43.258821157402444], [4.015407696185386, 43.2767208010392], [4.06316280730122, 43.42358105442678], [4.110894246531354, 43.57038983482697], [4.158693425647685, 43.71716151753356], [4.206996096092864, 43.863816709106246], [4.25510711981884, 44.01056802547626], [4.304059103877955, 44.15719324523995], [4.33312976289036, 44.24539400406059], [4.374940058749246, 44.24514241182411], [4.352494032483418, 43.256794294201605], [4.009598665597077, 43.258821157402444]]]}
  Bounding Box: [4.009598665597077, 43.256794294201605, 4.374940058749246, 44.24539400406059]
  Datetime: 2025-10-17 10:30:41.024000+00:00
  Cloud Cover: 0.079225
  Product HREF: https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/cpm_v270/S2A_MSIL2A_20251017T103041_N0511_R108_T31TEJ_20251017T155021.zarr
Item ID: S2C_MSIL2A_20251015T103031_N0511_R108_T31TEJ_20251015T161117
  Geometry: {'type': 'Polygon', 'coordinates': [[[4.009840633040255, 43.25881972712676], [4.048973646286205, 43.38025370013829], [4.096295595835079, 43.527178549418565], [4.14381103027854, 43.674062901396816], [4.170348465596649, 43.75558371079492], [4.192183484413271, 43.82252533122775], [4.240110737137591, 43.969365249374015], [4.260182216799283, 44.02997724998347], [4.288269822020025, 44.1144251058632], [4.3313063232805, 44.245404976555285], [4.374940058749246, 44.24514241182411], [4.352494032483418, 43.256794294201605], [4.009840633040255, 43.25881972712676]]]}
  Bounding Box: [4.009840633040255, 43.256794294201605, 4.374940058749246, 44.245404976555285]
  Datetime: 2025-10-15 10:30:31.025000+00:00
  Cloud Cover: 0.027531
  Product HREF: https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/cpm_v270/S2C_MSIL2A_20251015T103031_N0511_R108_T31TEJ_20251015T161117.zarr
Item ID: S2A_MSIL2A_20251010T104001_N0511_R008_T31TEJ_20251010T172518
  Geometry: {'type': 'Polygon', 'coordinates': [[[2.999751647531111, 43.72824351260538], [3.196169918404005, 43.693363313713085], [3.196219759900159, 43.69350318669725], [3.199216067303486, 43.6929703361596], [3.199403482706296, 43.69349686389792], [3.199622798471584, 43.69345793364484], [3.199877147377955, 43.694171491838055], [3.487026245631773, 43.64098954616449], [3.487166955732334, 43.64137999270429], [3.490161997825271, 43.6408255847484], [3.501472739139136, 43.63862606856076], [3.501065580989618, 43.637496909841886], [3.501109162888151, 43.63748843936892], [3.500891007117723, 43.63688283808694], [3.770586220690224, 43.584471752404816], [3.770744522107497, 43.58490579503936], [3.773369769440987, 43.58439513744098], [3.773831933369448, 43.5856623909434], [3.776305013657012, 43.585181469315465], [3.789823927300506, 43.58240661765728], [3.789700888334701, 43.58206987103524], [4.074260310509901, 43.52381547570125], [4.074225034692251, 43.52372015894707], [4.074455899276074, 43.52367105089948], [4.074437013109758, 43.523620117506105], [4.074518786667455, 43.523602737473446], [4.074005415212866, 43.52221593164724], [4.344670655368001, 43.46506504133068], [4.344710603564738, 43.46517165126475], [4.344959885842306, 43.46511905967128], [4.344990861999915, 43.46520182125857], [4.3473700608977, 43.464699065346885], [4.347705994578857, 43.46559640436082], [4.357187636230289, 43.46346403490077], [4.352494032483418, 43.256794294201605], [2.999753565511865, 43.26479037737325], [2.999751647531111, 43.72824351260538]]]}
  Bounding Box: [2.999751647531111, 43.256794294201605, 4.357187636230289, 43.72824351260538]
  Datetime: 2025-10-10 10:40:01.024000+00:00
  Cloud Cover: 0.00881
  Product HREF: https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/cpm_v270/S2A_MSIL2A_20251010T104001_N0511_R008_T31TEJ_20251010T172518.zarr
Item ID: S2C_MSIL2A_20251005T102921_N0511_R108_T31TEJ_20251005T155605
  Geometry: {'type': 'Polygon', 'coordinates': [[[3.996902592802162, 43.25889620421039], [4.002550239623694, 43.2763128635939], [4.050358167958506, 43.42322155121074], [4.098264130761962, 43.57007655526528], [4.146294506040694, 43.71688208047622], [4.194763315665021, 43.86358218775699], [4.243058650581385, 44.010342845215725], [4.292051709437927, 44.15699201427013], [4.321205161369018, 44.24546576000262], [4.374940058749246, 44.24514241182411], [4.352494032483418, 43.256794294201605], [3.996902592802162, 43.25889620421039]]]}
  Bounding Box: [3.996902592802162, 43.256794294201605, 4.374940058749246, 44.24546576000262]
  Datetime: 2025-10-05 10:29:21.025000+00:00
  Cloud Cover: 0.620101
  Product HREF: https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/cpm_v270/S2C_MSIL2A_20251005T102921_N0511_R108_T31TEJ_20251005T155605.zarr
Item ID: S2A_MSIL2A_20250907T103031_N0511_R108_T31TEJ_20250907T142815
  Geometry: {'type': 'Polygon', 'coordinates': [[[4.007933578555628, 43.25883099977413], [4.053914061479624, 43.400130355868015], [4.10200129449421, 43.5470258320004], [4.150286246823151, 43.693872715999696], [4.198977069700992, 43.84056658443097], [4.246877625970667, 43.98740530584433], [4.295831549412578, 44.13392803962985], [4.332295072100586, 44.245399026788085], [4.374940058749246, 44.24514241182411], [4.352494032483418, 43.256794294201605], [4.007933578555628, 43.25883099977413]]]}
  Bounding Box: [4.007933578555628, 43.256794294201605, 4.374940058749246, 44.245399026788085]
  Datetime: 2025-09-07 10:30:31.024000+00:00
  Cloud Cover: 4.397467
  Product HREF: https://data.eodc.eu/collections/EOPF_ZARR/products/cpm_v270/S02MSIL2A/2025/09/07/S2A_MSIL2A_20250907T103031_N0511_R108_T31TEJ_20250907T142815.zarr
Item ID: S2A_MSIL2A_20250530T103041_N0511_R108_T31TEJ_20250530T142711
  Geometry: {'type': 'Polygon', 'coordinates': [[[3.996151030886651, 43.258900646711915], [4.023793990467271, 43.34402002455392], [4.071660137619537, 43.490925621177965], [4.11966554524329, 43.63778121797476], [4.168007854238619, 43.7844957106434], [4.216194567973007, 43.93121158394251], [4.26500308831756, 44.07775724585984], [4.313670484861823, 44.22443496225287], [4.320641676579815, 44.24546915075591], [4.374940058749246, 44.24514241182411], [4.352494032483418, 43.256794294201605], [3.996151030886651, 43.258900646711915]]]}
  Bounding Box: [3.996151030886651, 43.256794294201605, 4.374940058749246, 44.24546915075591]
  Datetime: 2025-05-30 10:30:41.024000+00:00
  Cloud Cover: 0.023985
  Product HREF: https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/cpm_v270/S2A_MSIL2A_20250530T103041_N0511_R108_T31TEJ_20250530T142711.zarr
Item ID: S2A_MSIL2A_20250430T102701_N0511_R108_T31TEJ_20250430T190517
  Geometry: {'type': 'Polygon', 'coordinates': [[[3.999719577844396, 43.258879552940435], [4.01066893050314, 43.292618330840085], [4.058667839661785, 43.43963555941233], [4.106914943696001, 43.586620584359615], [4.15531245648374, 43.73349867183368], [4.203877089687448, 43.88029686606191], [4.25234288801703, 44.02709602620672], [4.301478140176486, 44.17367392313218], [4.325129722814647, 44.24544214406853], [4.374940058749246, 44.24514241182411], [4.352494032483418, 43.256794294201605], [3.999719577844396, 43.258879552940435]]]}
  Bounding Box: [3.999719577844396, 43.256794294201605, 4.374940058749246, 44.24544214406853]
  Datetime: 2025-04-30 10:27:01.024000+00:00
  Cloud Cover: 0.047449
  Product HREF: https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/cpm_v270/S2A_MSIL2A_20250430T102701_N0511_R108_T31TEJ_20250430T190517.zarr
Item ID: S2A_MSIL2A_20250420T103041_N0511_R108_T31TEJ_20250420T202912
  Geometry: {'type': 'Polygon', 'coordinates': [[[3.999054267090983, 43.258883485609154], [4.046727311476494, 43.40550124694311], [4.094877054390124, 43.55241547881358], [4.143268956901982, 43.69933392159373], [4.192035284878529, 43.846133169476616], [4.216811041370139, 43.92123219977511], [4.240519415716781, 43.99294472589098], [4.249841276796993, 44.02070565347856], [4.289983445930977, 44.139607974567234], [4.324815706219906, 44.24544403365424], [4.374940058749246, 44.24514241182411], [4.352494032483418, 43.256794294201605], [3.999054267090983, 43.258883485609154]]]}
  Bounding Box: [3.999054267090983, 43.256794294201605, 4.374940058749246, 44.24544403365424]
  Datetime: 2025-04-20 10:30:41.024000+00:00
  Cloud Cover: 1.786733
  Product HREF: https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/cpm_v270/S2A_MSIL2A_20250420T103041_N0511_R108_T31TEJ_20250420T202912.zarr
Item ID: S2C_MSIL2A_20250418T103041_N0511_R108_T31TEJ_20250418T160655
  Geometry: {'type': 'Polygon', 'coordinates': [[[3.996524135750957, 43.258898441279605], [4.029442036161715, 43.36057738888665], [4.077147593719051, 43.50739796354541], [4.125125045129605, 43.6542448039122], [4.173779370552519, 43.80100107173931], [4.221976447715368, 43.94791608964147], [4.270987202437039, 44.09461300238891], [4.31946199620316, 44.24146793252264], [4.320796341402683, 44.24546822006485], [4.374940058749246, 44.24514241182411], [4.352494032483418, 43.256794294201605], [3.996524135750957, 43.258898441279605]]]}
  Bounding Box: [3.996524135750957, 43.256794294201605, 4.374940058749246, 44.24546822006485]
  Datetime: 2025-04-18 10:30:41.025000+00:00
  Cloud Cover: 0.165078
  Product HREF: https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/cpm_v270/S2C_MSIL2A_20250418T103041_N0511_R108_T31TEJ_20250418T160655.zarr
Item ID: S2B_MSIL2A_20250406T103629_N0511_R008_T31TEJ_20250406T133008
  Geometry: {'type': 'Polygon', 'coordinates': [[[2.999749474120944, 44.25341758107761], [4.374940058749246, 44.24514241182411], [4.352494032483418, 43.256794294201605], [2.999753565511865, 43.26479037737325], [2.999749474120944, 44.25341758107761]]]}
  Bounding Box: [2.999749474120944, 43.256794294201605, 4.374940058749246, 44.25341758107761]
  Datetime: 2025-04-06 10:36:29.024000+00:00
  Cloud Cover: 0.004834
  Product HREF: https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/cpm_v270/S2B_MSIL2A_20250406T103629_N0511_R008_T31TEJ_20250406T133008.zarr

🧭 Walkthrough

The cells below walk through each processing step in detail for the first scene returned by the STAC query. This single-scene walkthrough illustrates the full pipeline — angle interpolation, input assembly, NN inference, and output saving — before the batch processing section repeats the same logic efficiently for every scene.

# get the url to the product of the first item in the search results
product = items[0].assets["product"]
print(f"  Product HREF: {product.href}")

# lazily open the product zarr with xarray
dt = xr.open_datatree(product.href, engine="zarr")
dt

  Product HREF: https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:notebook-data/tutorial_data/cpm_v270/S2A_MSIL2A_20251030T104211_N0511_R008_T31TEJ_20251030T144716.zarr

The EOPF Zarr store is opened lazily as an xarray.DataTree — a hierarchical tree of Dataset nodes that mirrors the internal Zarr group structure. No pixel data is downloaded at this point; arrays are only fetched when explicitly accessed or computed.

The main sub-trees used in this notebook are:

Group path	Contents
`measurements/reflectance/r10m`	Bottom-of-Atmosphere (BOA) reflectance at 10 m — bands B02, B03, B04, B08
`measurements/reflectance/r20m`	BOA reflectance at 20 m — bands B05, B06, B07, B8A, B11, B12
`conditions/geometry`	Sun and per-band, per-detector view angle grids (23 × 23 points at ~5 km spacing)
`conditions/mask/detector_footprint/r10m` and `r20m`	Per-detector binary footprint masks indicating which detector module observed each pixel
`quality/l2a_quicklook/r60m`	True-Colour Image (TCI) at 60 m for quick visual inspection

Reflectance scaling — bands are stored as scaled uint16 integers. The physical BOA reflectance is recovered as: reflectance = (raw_uint16 − 1000) / 10 000 (BOA_ADD_OFFSET = −1000, QUANTIFICATION_VALUE = 10 000, per the Sentinel-2 Product Specification.)

Build a True-Color Composite

The true-colour composite (B04 red, B03 green, B02 blue) at 60 m resolution serves as a quick visual sanity check to confirm:

The correct granule was selected (scene footprint covers the AOI)
Cloud cover is within acceptable limits for the LAI retrieval
The image is correctly georeferenced when overlaid on the CartoDB basemap

The composite is built directly from the measurements/reflectance/r60m group and displayed with a 2–98 % percentile linear stretch for natural-looking contrast. The LAI product is computed at 20 m.

Note: If you see the message Make this Notebook Trusted to load map: File -> Trust Notebook instead of a map after running the next cell, run the code in folium_test.ipynb and then click back on this notebook’s tab. The map should then be displayed.

# ---- Build a true-colour RGB composite from 60 m reflectance bands ----
# R = B04 (red), G = B03 (green), B = B02 (blue)
r60m = dt["measurements/reflectance/r60m"]
band_r = r60m["b04"].values  # red   – BOA reflectance, float32
band_g = r60m["b03"].values  # green
band_b = r60m["b02"].values  # blue

rgb = np.stack([band_r, band_g, band_b], axis=-1)  # (y, x, 3)
del band_r, band_g, band_b

# 2–98 % percentile linear stretch across all three bands for display contrast
p2, p98 = np.nanpercentile(rgb, 2), np.nanpercentile(rgb, 98)
tci_rgb = np.clip((rgb - p2) / (p98 - p2) * 255, 0, 255).astype(np.uint8)
del rgb

# ---- Compute WGS-84 bounding box from the 60 m coordinate arrays ----
ref_da = r60m["b02"]
x_coords = ref_da.coords["x"].values
y_coords = ref_da.coords["y"].values
half_px_x = abs(float(x_coords[1] - x_coords[0])) / 2
half_px_y = abs(float(y_coords[0] - y_coords[1])) / 2
x_min = float(x_coords.min()) - half_px_x
x_max = float(x_coords.max()) + half_px_x
y_min = float(y_coords.min()) - half_px_y
y_max = float(y_coords.max()) + half_px_y

tci_product_id = parse_sentinel2_product_name(Path(product.href))
tci_epsg = epsg_from_tile_id(tci_product_id.tile_id)
tci_transformer = Transformer.from_crs(f"EPSG:{tci_epsg}", "EPSG:4326", always_xy=True)
tci_corners = [
    tci_transformer.transform(x, y)
    for x, y in [(x_min, y_min), (x_min, y_max), (x_max, y_min), (x_max, y_max)]
]
tci_lons, tci_lats = zip(*tci_corners)
tci_bounds = [[min(tci_lats), min(tci_lons)], [max(tci_lats), max(tci_lons)]]
tci_center = [(min(tci_lats) + max(tci_lats)) / 2, (min(tci_lons) + max(tci_lons)) / 2]

# ---- Folium map ----
tci_map = folium.Map(location=tci_center, zoom_start=9, tiles="CartoDB positron")
folium.raster_layers.ImageOverlay(
    image=tci_rgb,
    bounds=tci_bounds,
    opacity=1.0,
    name="True Colour Composite (B4/B3/B2)",
).add_to(tci_map)
folium.LayerControl().add_to(tci_map)
tci_map

Make this Notebook Trusted to load map: File -> Trust Notebook

🛰️ Extract Sensor Metadata for Computations

Map detector indexes based on the data type of the detector array in the source dataset.

match dt["conditions/geometry"]["detector"].dtype:
    case "<U3":
        DETECTORS = ("d01", "d02", "d03", "d04", "d05", "d06", "d07", "d08", "d09", "d10", "d11", "d12")
    case "O":
        DETECTORS = ("d01", "d02", "d03", "d04", "d05", "d06", "d07", "d08", "d09", "d10", "d11", "d12")        
    case "int64":
        DETECTORS = (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12)
    case _:
        raise ValueError(f"Unexpected dtype for detector coordinate: {dt["conditions/geometry"]["detector"].dtype}")
print(f"Using detector labels: {DETECTORS}")

Using detector labels: ('d01', 'd02', 'd03', 'd04', 'd05', 'd06', 'd07', 'd08', 'd09', 'd10', 'd11', 'd12')

Each Sentinel-2 satellite has 12 focal-plane detector modules arranged in a staggered row across the instrument’s focal plane. As the satellite moves along its orbit, different detector modules observe different columns of the swath. Because each detector views the ground from a slightly different angle, the correct per-pixel view geometry depends on which module made the observation.

In the EOPF Zarr, the detector footprint is stored as a label image (uint8 integers 1–12, or string labels "d01"–"d12" depending on the product version). The viewing-angle grids in conditions/geometry are stored separately for each detector. The code block above detects which label format is used in the open store and sets DETECTORS accordingly, so that the angle-selection logic in the next section works correctly regardless of the data version.

Process the Sun and View Zenith and Azimuth Angles

The sun and view zenith and azimuth angles are needed for the Leaf Area Index model. These values are provided in a 23 x 23 grid with a resolution of 5,000 x 5,000 pixels and must be interpolated to the same resolution as the reflectance bands.

The first step is to get the target number of pixels, depending on whether the desired output spatial resolution is 10 meters or 20 meters.

match SPATIAL_RESOLUTION:
    case 10:
        # set the path to the reflectance group for 10m spatial resolution
        reflectance_group_path = "measurements/reflectance/r10m"
        # get the list of bands needed for LAI calculation at 10m resolution
        bands = SENTINEL2_LAI_BANDS_10M
    case 20:
        # set the path to the reflectance group for 20m spatial resolution
        reflectance_group_path = "measurements/reflectance/r20m"
        # get the list of bands needed for LAI calculation at 20m resolution
        bands = SENTINEL2_LAI_BANDS_20M
    case _:
        raise ValueError("SPATIAL_RESOLUTION must be either 10 or 20 meters for Sentinel-2 data.")

target_width = dt[reflectance_group_path].sizes["x"]
target_height = dt[reflectance_group_path].sizes["y"]

target_shape = (target_height, target_width)
print(f"Target shape for interpolation: {target_shape}")

Target shape for interpolation: (5490, 5490)

The SPATIAL_RESOLUTION setting determines both the output pixel spacing and which neural network configuration is used:

Resolution	Spectral bands	NN input features	Standard tile grid	NN weight set
20 m (recommended)	B03, B04, B05, B06, B07, B8A, B11, B12	8 reflectance + 3 angles = 11	5 490 × 5 490	`S2x_20m_LAI`
10 m	B03, B04, B08	3 reflectance + 3 angles = 6	10 980 × 10 980	`S2x_10m_LAI`

The 20 m configuration is the reference setting used by SNAP BiophysicalOp. It benefits from richer spectral coverage (Red-Edge and SWIR bands are particularly sensitive to canopy structure) and is recommended for quantitative analysis. The 10 m option trades spectral richness for higher spatial detail.

Sun Angles

The Sun Zenith Angle (SZA) and Sun Azimuth Angle (SAA) describe the position of the sun at the moment of acquisition. In the EOPF Zarr, they are stored as a sparse 23 × 23 grid (one value roughly every 5 km across the 100 km tile) and must be bilinearly interpolated to the full image resolution before use.

Sun angles are scene-wide — the same grid applies to all spectral bands and all detector modules.

# Source values
da_sun_zenith_src = dt["conditions/geometry"].sel(angle="zenith").sun_angles
da_sun_azimuth_src = dt["conditions/geometry"].sel(angle="azimuth").sun_angles

# Interpolate to target shape
da_sun_zenith = bilinear_interpolate_xy(da_sun_zenith_src, target_shape)
da_sun_azimuth = bilinear_interpolate_xy(da_sun_azimuth_src, target_shape)

# Extract values as float32 numpy arrays (float64 is unnecessary precision)
sun_zenith = da_sun_zenith.values.astype(np.float32)
sun_azimuth = da_sun_azimuth.values.astype(np.float32)

# Free DataArrays no longer needed
del da_sun_zenith_src, da_sun_azimuth_src, da_sun_zenith, da_sun_azimuth

View Angles

The View Zenith Angle (VZA) and View Azimuth Angle (VAA) describe the satellite-to-ground observation direction for each pixel. Unlike the sun angles, viewing geometry in Sentinel-2 is band- and detector-specific: each of the 12 focal-plane detector modules has its own angle grid, and different spectral bands are mapped to slightly different positions on the focal plane.

The processing below handles this complexity in four steps for each of the 8 LAI bands: 1. Extract the per-detector 23 × 23 angle grids from the Zarr store 2. Bilinearly interpolate each detector’s grid to the full image resolution 3. Use the detector footprint masks to select the correct interpolated angle for each pixel 4. Store the merged per-band angle arrays for subsequent averaging

For each of the 8 LAI spectral bands the loop below performs:

Extract source grids — retrieve the per-detector 23 × 23 zenith and azimuth arrays from conditions/geometry.
Bilinear interpolation — resample each detector’s grid from 23 × 23 points up to the full target_shape using a bicubic-extended bilinear scheme.
Build detector masks — derive a binary mask for each of the 12 detectors from the footprint label image; for 10 m bands working at a 20 m target, the masks are spatially downsampled to match.
Select by detector — for every pixel, pick the interpolated angle value from the grid of the detector that covered that pixel (according to the footprint mask). Pixels covered by no detector receive NaN.
Store — save the merged zenith and azimuth arrays, keyed by band name, ready for averaging in the next step.

Note: RuntimeWarning: invalid value encountered in subtract ext[2:, :] = np.where(v, 2 * ext[1:-1, :] - ext[:-2, :], ext[2:, :]) is expected due to np.inf values being present in the data. These values are masked out during processing and does not affect the results.

weighted_view_zeniths = {}
weighted_view_azimuths = {}

for band in bands:
    print(f"Processing view angles for band: {band}...")
    band_group = SENTINEL2_BAND_GROUPS[band]

    # 1. get the viewing angles for the current band
    print("  Extracting source viewing angles...", end="\r")
    try:
        da_view_zenith_src = dt["conditions/geometry"].sel(angle="zenith", band=SENTINEL2_GEOMETRY_BAND_INDICES[band]).viewing_incidence_angles
        da_view_azimuth_src = dt["conditions/geometry"].sel(angle="azimuth", band=SENTINEL2_GEOMETRY_BAND_INDICES[band]).viewing_incidence_angles
    except KeyError:
        print(f"  [!] No viewing geometry for band {band} (index {SENTINEL2_GEOMETRY_BAND_INDICES[band]}), skipping.")
        continue
    print("  [x] Extracted source viewing angles.")

    # 2. interpolate to target shape
    print("  Interpolating viewing angles to target shape...", end="\r")
    da_view_zenith = bilinear_interpolate_xy(da_view_zenith_src, target_shape)
    da_view_azimuth = bilinear_interpolate_xy(da_view_azimuth_src, target_shape)
    print("  [x] Interpolated viewing angles to target shape.")

    # 3. calculate detector masks
    print("  Calculating detector masks...", end="\r")
    da_footprint = dt[f"conditions/mask/detector_footprint/{band_group}"][band]
    da_detector_masks = add_detector_dimension(da_footprint, dt)
    if int(band_group[1:3]) != SPATIAL_RESOLUTION:
        # For 10m bands, we need to select the appropriate band from the detector masks
        da_detector_masks = downsample_detector_masks(da_detector_masks, target_shape)
    print("  [x] Calculated detector masks.")

    # 4. select the viewing angles based on the detector masks
    print("  Merging view angles...", end="\r")
    da_merged_view_zenith = select_angle_by_detector_mask(da_view_zenith, da_detector_masks)
    da_merged_view_azimuth = select_angle_by_detector_mask(da_view_azimuth, da_detector_masks)
    print("  [x] Merged view angles.")

    # 7. store the results in dictionaries keyed by band
    weighted_view_zeniths[band] = da_merged_view_zenith
    weighted_view_azimuths[band] = da_merged_view_azimuth

Processing view angles for band: b03...
  Extracting source viewing angles...
  [x] Extracted source viewing angles.
  Interpolating viewing angles to target shape...
  [x] Interpolated viewing angles to target shape.
  Calculating detector masks...
  [x] Calculated detector masks.
  Merging view angles...
  [x] Merged view angles.
Processing view angles for band: b04...
  Extracting source viewing angles...
  [x] Extracted source viewing angles.
  Interpolating viewing angles to target shape...
  [x] Interpolated viewing angles to target shape.
  Calculating detector masks...
  [x] Calculated detector masks.
  Merging view angles...
  [x] Merged view angles.
Processing view angles for band: b05...
  Extracting source viewing angles...
  [x] Extracted source viewing angles.
  Interpolating viewing angles to target shape...
  [x] Interpolated viewing angles to target shape.
  Calculating detector masks...
  [x] Calculated detector masks.
  Merging view angles...
  [x] Merged view angles.
Processing view angles for band: b06...
  Extracting source viewing angles...
  [x] Extracted source viewing angles.
  Interpolating viewing angles to target shape...
  [x] Interpolated viewing angles to target shape.
  Calculating detector masks...
  [x] Calculated detector masks.
  Merging view angles...
  [x] Merged view angles.
Processing view angles for band: b07...
  Extracting source viewing angles...
  [x] Extracted source viewing angles.
  Interpolating viewing angles to target shape...
  [x] Interpolated viewing angles to target shape.
  Calculating detector masks...
  [x] Calculated detector masks.
  Merging view angles...
  [x] Merged view angles.
Processing view angles for band: b8a...
  Extracting source viewing angles...
  [x] Extracted source viewing angles.
  Interpolating viewing angles to target shape...
  [x] Interpolated viewing angles to target shape.
  Calculating detector masks...
  [x] Calculated detector masks.
  Merging view angles...
  [x] Merged view angles.
Processing view angles for band: b11...
  Extracting source viewing angles...
  [x] Extracted source viewing angles.
  Interpolating viewing angles to target shape...
  [x] Interpolated viewing angles to target shape.
  Calculating detector masks...
  [x] Calculated detector masks.
  Merging view angles...
  [x] Merged view angles.
Processing view angles for band: b12...
  Extracting source viewing angles...
  [x] Extracted source viewing angles.
  Interpolating viewing angles to target shape...
  [x] Interpolated viewing angles to target shape.
  Calculating detector masks...
  [x] Calculated detector masks.
  Merging view angles...
  [x] Merged view angles.

The 8 per-band merged view-angle grids are averaged pixel-by-pixel to produce a single mean view zenith and mean view azimuth map for the scene. This matches the SNAP BiophysicalOp approach, which computes view_zenith_mean across all selected bands before feeding into the neural network.

Pixels with no detector coverage — gaps between adjacent detector swath columns — have NaN view angles. These are tracked via gap_mask and will be set to NaN in the final LAI output after inference.

da_view_zenith_mean = compute_view_angle_mean(weighted_view_zeniths)
da_view_azimuth_mean = compute_view_angle_mean(weighted_view_azimuths)

# extract values as float32 numpy arrays
view_zenith_mean = da_view_zenith_mean.values.astype(np.float32)
view_azimuth_mean = da_view_azimuth_mean.values.astype(np.float32)

# Free intermediate dictionaries and DataArrays no longer needed
del weighted_view_zeniths, weighted_view_azimuths
del da_view_zenith_mean, da_view_azimuth_mean

Relative Azimuth Angle

The neural network uses the Relative Azimuth Angle (RAA), defined as the difference between the sun and view azimuth directions:

RAA = SAA − VAA

RAA captures the illumination–observation geometry independently of the individual azimuth values. An RAA near 0° corresponds to forward-scatter geometry (sun and sensor on opposite sides of nadir); RAA near ±180° corresponds to back-scatter geometry. Both strongly influence canopy reflectance and are important features for the LAI retrieval.

realtive_azimuth = (sun_azimuth - view_azimuth_mean).astype(np.float32)

Trigonometric Angle Features

The neural network receives the cosines of the three angles rather than the angles themselves (in degrees), following the SNAP BiophysicalOp convention:

Feature	Physical meaning
`cos(VZA)`	Inversely related to the optical path length through the canopy; decreases as the sensor looks more obliquely
`cos(SZA)`	Proportional to the solar irradiance on a horizontal surface; decreases toward the terminator
`cos(RAA)`	Encodes the forward/back-scatter geometry; ranges from −1 (forward scatter) to +1 (back scatter)

Taking cosines approximates a more linear relationship between angles and top-of-canopy reflectance, improving NN training convergence. It also provides natural physical bounds (−1 to +1) that complement the [−1, +1] normalisation applied to all features before inference.

Pixels where cos(VZA) is NaN (gap pixels with no detector coverage) are tracked in gap_mask and excluded from the final output.

cos_view_zentith_mean = np.cos(np.radians(view_zenith_mean, dtype=np.float32), dtype=np.float32)
cos_sun_zenith = np.cos(np.radians(sun_zenith, dtype=np.float32), dtype=np.float32)
cos_relative_azimuth = np.cos(np.radians(realtive_azimuth, dtype=np.float32), dtype=np.float32)

# Track pixels with no detector coverage (NaN view angles).
gap_mask = np.isnan(cos_view_zentith_mean)

# Free source angle arrays no longer needed
del view_zenith_mean, sun_zenith, sun_azimuth, view_azimuth_mean, realtive_azimuth

📚 Combine Extracted Metadata with Observation Data for Neural Network Inference

Assemble the Neural Network Input Array

The SNAP BiophysicalOp neural network expects a fixed-length feature vector for each pixel. At 20 m resolution this vector has 11 elements:

Index	Feature	Band / Source	Central Wavelength (nm)
0	B03	Green	560 nm
1	B04	Red	665 nm
2	B05	Red-Edge 1	705 nm
3	B06	Red-Edge 2	740 nm
4	B07	Red-Edge 3	783 nm
5	B8A	NIR narrow	865 nm
6	B11	SWIR 1	1 610 nm
7	B12	SWIR 2	2 190 nm
8	cos(VZA)	Mean view zenith angle	—
9	cos(SZA)	Sun zenith angle	—
10	cos(RAA)	Relative azimuth angle	—

The bands are stacked along the first axis to produce a (11, H, W) array in float32. The Red-Edge and SWIR bands carry most of the LAI sensitivity: Red-Edge bands respond to chlorophyll content and canopy structure, while SWIR bands are sensitive to water content and canopy depth.

data = []
# reflectance bands - load as float32
for band in bands:
    data.append(dt[reflectance_group_path][band].values.astype(np.float32))

# sun and view angles (already float32)
data.append(cos_view_zentith_mean)
data.append(cos_sun_zenith)
data.append(cos_relative_azimuth)

input_array = np.stack(data, axis=0)

# Free intermediate lists and cosine arrays
del data, cos_view_zentith_mean, cos_sun_zenith, cos_relative_azimuth

⚙️ Run Stacked Tiles through LAI Neural Network

Neural Network Architecture

The LAI retrieval uses a two-layer feed-forward network with the following structure:

Input (11)  →  [Linear + tanh]  →  Hidden (5)  →  [Linear]  →  Output (1)

Forward pass — mirroring SNAP’s NeuralNetAlgo:

Input normalisation — each of the 11 features is linearly rescaled from its per-feature training range [x_min, x_max] to [−1, +1] using coefficients from the LAI_Normalization resource file.
Layer 1 (5 hidden neurons) — h = tanh(W₁ · x_norm + b₁), where W₁ is 5 × 11.
Layer 2 (1 output, no activation) — z = W₂ · h + b₂, where W₂ is length 5.
Output denormalisation — z is rescaled back to physical LAI units (m² m⁻²) using coefficients from the LAI_Denormalization resource file.

Sensor-specific weight files (LAI_Weights_Layer1_Neurons, LAI_Weights_Layer1_Bias, LAI_Weights_Layer2_Neurons, LAI_Weights_Layer1_Bias, LAI_Normalization, LAI_Denormalization) are loaded from resources/esa/biophysical/3_0/<sensor>/LAI/. The inference is compiled and executed by the Modular MAX engine using a static computation graph for efficiency.

Note: This graph has to be compiled and takes longer to run the first time the code is executed than on subsequent runs. Also, it uses a GPU for accelerated computation, if a compatible one is available (check Modular’s Documentation for compatible NVIDIA, AMD, and Apple GPUs).

Input Validation (before inference)

check_input_out_of_range runs two checks on the reflectance bands:

Stage 1 — Min/Max bounds: each of the 8 reflectance bands must lie within the per-band training domain bounds defined in the LAI_DefinitionDomain_MinMax resource file. Pixels violating any band bound receive the input_out_of_range flag.
Stage 2 — Convex hull: the combined 8-band vector is projected to a discrete grid and checked against the set of valid combinations in the LAI_DefinitionDomain_Grid resource file. This detects spectral combinations that are individually within bounds but collectively outside the training data distribution.

The network still runs for all pixels, including those flagged as out-of-range, consistent with the SNAP reference implementation.

Output Validation (after inference)

check_output_out_of_range applies physical plausibility bounds to the NN output:

Flag	Meaning
`output_set_to_min`	LAI was below the minimum allowed value and clamped up
`output_set_to_max`	LAI exceeded the maximum allowed value and clamped down
`output_too_low`	Output is below a wider “implausible” threshold
`output_too_high`	Output exceeds a wider “implausible” threshold

s2_product_id = parse_sentinel2_product_name(Path(product.href))
if bands == SENTINEL2_LAI_BANDS_10M:
    match s2_product_id.mission_id:
        case "S2A":
            sat_sensor = SatelliteSensors.S2A_10m
        case "S2B":
            sat_sensor = SatelliteSensors.S2B_10m
        case "S2C":
            sat_sensor = SatelliteSensors.S2C_10m
        case _:
            raise ValueError(f"Unsupported Sentinel-2 mission ID: {s2_product_id.mission_id}")
else:  # bands == SENTINEL2_LAI_BANDS_20M
    match s2_product_id.mission_id:
        case "S2A":
            sat_sensor = SatelliteSensors.S2A
        case "S2B":
            sat_sensor = SatelliteSensors.S2B
        case "S2C":
            sat_sensor = SatelliteSensors.S2C
        case _:
            raise ValueError(f"Unsupported Sentinel-2 mission ID: {s2_product_id.mission_id}")

lai_model = BiophysicalVariables.LAI

biophysical_op_process = BiophysicalOpProcess(
    biophysical_op=lai_model,
    sensor=sat_sensor,
    input_data=input_array,
)

biophysical_op_process = process(biophysical_op_process)

# Overwrite gap-pixel outputs with NaN; clear any spurious threshold flags
# that the NN produced from NaN inputs rather than real out-of-range values.
if gap_mask.any():
    result = biophysical_op_process.result
    if result.data is not None:
        result.data[gap_mask] = np.nan
    for attr in (
        "output_threshold_set_to_min_output_mask",
        "output_threshold_set_to_max_output_mask",
        "output_too_low_mask",
        "output_too_high_mask",
    ):
        mask = getattr(result, attr)
        if mask is not None:
            mask[gap_mask] = False
            if not mask.any():
                setattr(result, attr, None)
                flag_attr = attr.replace("_mask", "")
                setattr(result, flag_attr, False)

After inference, biophysical_op_process.result holds:

Attribute	Shape / type	Description
`result.data`	`float32 (H, W)`	LAI in m² m⁻²; `NaN` for gap pixels (no detector coverage)
`result.input_out_of_range_mask`	`bool (H, W)`	Stage 1 or Stage 2 domain check failed
`result.output_threshold_set_to_min_output_mask`	`bool (H, W)`	LAI clipped to minimum physical bound
`result.output_threshold_set_to_max_output_mask`	`bool (H, W)`	LAI clipped to maximum physical bound
`result.output_too_low_mask`	`bool (H, W)`	LAI below the wider implausible-range threshold
`result.output_too_high_mask`	`bool (H, W)`	LAI above the wider implausible-range threshold

Gap pixels (swath edges with no detector footprint) have their LAI set to NaN and their threshold flags cleared to avoid artefacts at swath boundaries.

📘 Save Results to Zarr Format

We prepare the path to save to our S3 bucket:

output_filepath = f"{lai_model.value.label}_{Path(product.href).stem}.zarr"
print(output_filepath)
base_store = fs.get_mapper(f"{bucket}/{output_filepath}")
protocol = base_store.fs.protocol[0]
full_path = f"{protocol}://{base_store.root}"

LAI_S2A_MSIL2A_20251030T104211_N0511_R008_T31TEJ_20251030T144716.zarr

From the processed data, we bring our Dataset and save it to .zarr

final_store = prep_zarr(biophysical_op_process, dt, full_path, reflectance_group_path)
# final_store.to_zarr(base_store, mode="w") # Run this line to save the results to disc

final_store

🗺️ Visualize Results

Finally, to visualise the generated zarr:

# Derive EPSG
z_vis = final_store
epsg_str = dict(z_vis["crs"].attrs).get("epsg_code", "")
epsg_vis = int(epsg_str.split(":")[-1])

# Bounding box
x_vals = z_vis.x.values
y_vals = z_vis.y.values
dx = int(x_vals[1] - x_vals[0])
dy = int(y_vals[0] - y_vals[1])
x_min = float(x_vals.min()) - dx / 2
x_max = float(x_vals.max()) + dx / 2
y_min = float(y_vals.min()) - dy / 2
y_max = float(y_vals.max()) + dy / 2

transformer_vis = Transformer.from_crs(f"EPSG:{epsg_vis}", "EPSG:4326", always_xy=True)
corners_wgs84 = [
    transformer_vis.transform(x, y)
    for x, y in [(x_min, y_min), (x_min, y_max), (x_max, y_min), (x_max, y_max)]
]
lons, lats = zip(*corners_wgs84)
lat_min, lat_max = min(lats), max(lats)
lon_min, lon_max = min(lons), max(lons)
vis_bounds = [[lat_min, lon_min], [lat_max, lon_max]]
vis_center = [(lat_min + lat_max) / 2, (lon_min + lon_max) / 2]
print(f"WGS-84 bounds: {vis_bounds}")

DOWNSAMPLE = 10

WGS-84 bounds: [[43.25679429374438, 2.9997494741209363], [44.253417580616336, 4.374940058743222]]

lai_vmin, lai_vmax = 0.0, 8.0
cmap_lai = plt.get_cmap("YlGn")
norm_lai  = mcolors.Normalize(vmin=lai_vmin, vmax=lai_vmax)
colormap_vis = bcm.LinearColormap(
    ["#ffffcc", "#c2e699", "#78c679", "#31a354", "#006837"],
    vmin=lai_vmin, vmax=lai_vmax, caption="LAI (m\u00b2 m\u207b\u00b2)",
)

MASK_VARS_VIS = [v for v in ("input_out_of_range", "output_set_to_min", "output_set_to_max", "output_too_low", "output_too_high") if v in z_vis]
MASK_COLORS_VIS = {
    "input_out_of_range": (230,  26,  26),
    "output_set_to_min":  ( 26, 102, 230),
    "output_set_to_max":  (255, 128,   0),
    "output_too_low":     (178,   0, 230),
    "output_too_high":    (255, 217,   0),
}

legend_items_vis = "".join(
    f"<li style='margin:3px 0;display:flex;align-items:center;'>"
    f"<span style='background:#{r:02x}{g:02x}{b:02x};width:14px;height:14px;"
    f"min-width:14px;display:inline-block;margin-right:8px;border:1px solid #555;"
    f"border-radius:2px;'></span>{var.replace('_', ' ').title()}</li>"
    for var, (r, g, b) in ((v, MASK_COLORS_VIS[v]) for v in MASK_VARS_VIS)
)
legend_html_vis = (
    "<div style='position:fixed;bottom:30px;left:30px;z-index:9999;background:white;"
    "padding:10px 14px;border-radius:6px;font-size:13px;line-height:1.5;"
    "box-shadow:2px 2px 6px rgba(0,0,0,0.3);'>"
    "<b style='display:block;margin-bottom:6px;'>Quality Flags</b>"
    f"<ul style='list-style:none;margin:0;padding:0;'>{legend_items_vis}</ul>"
    "</div>"
)

# LAI overlay
lai_d = z_vis["LAI"].values[::DOWNSAMPLE, ::DOWNSAMPLE]
rgba_lai = cmap_lai(norm_lai(np.where(np.isnan(lai_d), 0.0, lai_d)))
rgba_lai[np.isnan(lai_d)] = 0
rgba_lai_u8 = (rgba_lai * 255).astype(np.uint8)

# TCI overlay (from walkthrough)
TCI_DS = 3
tci_d = tci_rgb[::TCI_DS, ::TCI_DS]
tci_rgba = np.dstack([tci_d, np.full(tci_d.shape[:2], 255, dtype=np.uint8)])

# Quality flags overlay
H_d, W_d = lai_d.shape
mask_rgba = np.zeros((H_d, W_d, 4), dtype=np.uint8)
for var in MASK_VARS_VIS:
    flag_mask = z_vis[var].values[::DOWNSAMPLE, ::DOWNSAMPLE].astype(bool)
    r, g, b = MASK_COLORS_VIS[var]
    mask_rgba[flag_mask] = [r, g, b, 217]

fmap_vis = folium.Map(location=vis_center, zoom_start=9, tiles=None)
folium.TileLayer("OpenStreetMap",    name="OpenStreetMap").add_to(fmap_vis)
folium.TileLayer("CartoDB voyager",  name="CartoDB Voyager").add_to(fmap_vis)
folium.TileLayer("CartoDB positron", name="CartoDB Positron").add_to(fmap_vis)
folium.raster_layers.ImageOverlay(
    image=tci_rgba, bounds=tci_bounds, opacity=1.0,
    name="True Colour Composite (B4/B3/B2)", show=True,
).add_to(fmap_vis)
folium.raster_layers.ImageOverlay(
    image=rgba_lai_u8, bounds=vis_bounds, opacity=0.85,
    name="LAI", show=True,
).add_to(fmap_vis)
folium.raster_layers.ImageOverlay(
    image=mask_rgba, bounds=vis_bounds, opacity=1.0,
    name="Quality Flags", show=True,
).add_to(fmap_vis)
colormap_vis.add_to(fmap_vis)
fmap_vis.get_root().html.add_child(Element(legend_html_vis))
folium.LayerControl(collapsed=False).add_to(fmap_vis)
display(fmap_vis)

Make this Notebook Trusted to load map: File -> Trust Notebook

⚙️⚙️ Batch Processing a Full Growing Season — Per-Tile Output

The walkthrough above illustrated the complete pipeline for a single scene. This section applies the same pipeline to every Sentinel-2 tile returned by the STAC query (April – October 2025).

Each STAC item is processed independently: the full tile is run through the BiophysicalOp pipeline and saved as a separate Zarr store under OUTPUT_DIR. No AOI clipping, tile merging, or time-series combination is performed — each output zarr contains one complete Sentinel-2 tile, consistent with the behaviour of the SNAP BiophysicalOp processor.

Existing stores are skipped (if they contain at least one finite LAI pixel), so the loop is safely resumable after interruption.

Batch Processing — Per-Tile Output

The process_scene function runs the full BiophysicalOp LAI pipeline for one STAC item, processing the entire Sentinel-2 tile without any AOI clipping.

Tip

If you would like to replicate the following workflow locally, feel free to activate them, as they produce around 1GB of data! Just remember to adapt the storage paths we set for the S3 Bucket!

# def process_scene(item, output_path: Path) -> Path:
#     """Run the full BiophysicalOp LAI pipeline for one STAC item.

#     Processes the entire Sentinel-2 tile without any AOI clipping,
#     consistent with the behaviour of the SNAP BiophysicalOp processor.

#     Parameters
#     ----------
#     item : pystac.Item
#         A single Sentinel-2 L2A STAC item.
#     output_path : Path
#         Destination Zarr store path for this tile.

#     Returns
#     -------
#     Path | None
#         The path that was written, or None if processing failed.
#     """
#     product = item.assets["product"]
#     try:
#         dt_s = xr.open_datatree(product.href, engine="zarr")
#     except Exception as e:
#         print(f"  [!] Error opening datatree for {product.href}: {e}")
#         return None

#     match SPATIAL_RESOLUTION:
#         case 20:
#             rg_path = "measurements/reflectance/r20m"
#             b_list = SENTINEL2_LAI_BANDS_20M
#         case 10:
#             rg_path = "measurements/reflectance/r10m"
#             b_list = SENTINEL2_LAI_BANDS_10M
#         case _:
#             raise ValueError(f"Unsupported SPATIAL_RESOLUTION: {SPATIAL_RESOLUTION}")

#     th = dt_s[rg_path].sizes["y"]
#     tw = dt_s[rg_path].sizes["x"]
#     tshape = (th, tw)

#     # Sun angles — bilinearly interpolate 23x23 grid to full image resolution
#     sun_zenith = bilinear_interpolate_xy(
#         dt_s["conditions/geometry"].sel(angle="zenith").sun_angles, tshape
#     ).values.astype(np.float32)
#     sun_azimuth = bilinear_interpolate_xy(
#         dt_s["conditions/geometry"].sel(angle="azimuth").sun_angles, tshape
#     ).values.astype(np.float32)

#     # View angles — per-band interpolation + detector-footprint merging
#     weighted_view_zeniths, weighted_view_azimuths = {}, {}
#     for band in b_list:
#         band_group = SENTINEL2_BAND_GROUPS[band]
#         try:
#             da_vz_src = dt_s["conditions/geometry"].sel(
#                 angle="zenith", band=SENTINEL2_GEOMETRY_BAND_INDICES[band]
#             ).viewing_incidence_angles
#             da_va_src = dt_s["conditions/geometry"].sel(
#                 angle="azimuth", band=SENTINEL2_GEOMETRY_BAND_INDICES[band]
#             ).viewing_incidence_angles
#         except KeyError:
#             continue
#         da_vz = bilinear_interpolate_xy(da_vz_src, tshape)
#         da_va = bilinear_interpolate_xy(da_va_src, tshape)
#         da_footprint = dt_s[f"conditions/mask/detector_footprint/{band_group}"][band]
#         da_det_masks = add_detector_dimension(da_footprint, dt_s)
#         if int(band_group[1:3]) != SPATIAL_RESOLUTION:
#             da_det_masks = downsample_detector_masks(da_det_masks, tshape)
#         weighted_view_zeniths[band] = select_angle_by_detector_mask(da_vz, da_det_masks)
#         weighted_view_azimuths[band] = select_angle_by_detector_mask(da_va, da_det_masks)

#     view_zenith_mean = compute_view_angle_mean(weighted_view_zeniths).values.astype(np.float32)
#     view_azimuth_mean = compute_view_angle_mean(weighted_view_azimuths).values.astype(np.float32)
#     del weighted_view_zeniths, weighted_view_azimuths

#     relative_azimuth = (sun_azimuth - view_azimuth_mean).astype(np.float32)
#     cos_vz = np.cos(np.radians(view_zenith_mean, dtype=np.float32), dtype=np.float32)
#     cos_sz = np.cos(np.radians(sun_zenith, dtype=np.float32), dtype=np.float32)
#     cos_ra = np.cos(np.radians(relative_azimuth, dtype=np.float32), dtype=np.float32)
#     gap_mask = np.isnan(cos_vz)
#     del view_zenith_mean, sun_zenith, sun_azimuth, view_azimuth_mean, relative_azimuth

#     # Assemble NN input: [8 reflectance bands | cos_vz | cos_sz | cos_ra]
#     data = [dt_s[rg_path][band].values.astype(np.float32) for band in b_list]
#     data += [cos_vz, cos_sz, cos_ra]
#     input_array = np.stack(data, axis=0)
#     del data, cos_vz, cos_sz, cos_ra

#     # Sensor / NN weights selection
#     pid = parse_sentinel2_product_name(Path(product.href))
#     if b_list == SENTINEL2_LAI_BANDS_10M:
#         match pid.mission_id:
#             case "S2A":
#                 sat_sensor = SatelliteSensors.S2A_10m
#             case "S2B":
#                 sat_sensor = SatelliteSensors.S2B_10m
#             case "S2C":
#                 sat_sensor = SatelliteSensors.S2C_10m
#             case _:
#                 raise ValueError(f"Unknown mission: {pid.mission_id}")
#     else:
#         match pid.mission_id:
#             case "S2A":
#                 sat_sensor = SatelliteSensors.S2A
#             case "S2B":
#                 sat_sensor = SatelliteSensors.S2B
#             case "S2C":
#                 sat_sensor = SatelliteSensors.S2C
#             case _:
#                 raise ValueError(f"Unknown mission: {pid.mission_id}")

#     bop = BiophysicalOpProcess(
#         biophysical_op=BiophysicalVariables.LAI,
#         sensor=sat_sensor,
#         input_data=input_array,
#     )
#     bop = process(bop)
#     del input_array

#     # Overwrite gap-pixel outputs with NaN; clear any spurious threshold flags
#     if gap_mask.any():
#         result = bop.result
#         if result.data is not None:
#             result.data[gap_mask] = np.nan
#         for attr in (
#             "output_threshold_set_to_min_output_mask",
#             "output_threshold_set_to_max_output_mask",
#             "output_too_low_mask",
#             "output_too_high_mask",
#         ):
#             mask = getattr(result, attr)
#             if mask is not None:
#                 mask[gap_mask] = False
#                 if not mask.any():
#                     setattr(result, attr, None)
#                     setattr(result, attr.replace("_mask", ""), False)

#     save_zarr(bop, dt_s, output_path, rg_path)
#     del dt_s, bop, gap_mask
#     return output_path

Per-Tile Batch Loop

The batch loop below iterates over every STAC item returned by the search query. Each item is processed independently by process_scene and saved as a separate Zarr store in OUTPUT_DIR. Existing stores that already contain valid LAI data are skipped, making the loop safely resumable after interruption.

Note: RuntimeWarning: invalid value encountered in subtract is expected due to np.inf values present in the data. These values are masked out during processing and do not affect the results.

# def _scene_has_valid_lai(path: Path) -> bool:
#     """Return True only if the saved zarr has at least one finite LAI pixel."""
#     try:
#         z = zarr.open(str(path), mode="r")
#         if "LAI" not in z:
#             return False
#         return bool(np.any(np.isfinite(z["LAI"][:])))
#     except Exception:
#         return False


# OUTPUT_DIR.mkdir(exist_ok=True)
# tile_zarr_paths = []

# for item in items:
#     product = item.assets["product"]
#     stem = Path(product.href).stem
#     output_path = OUTPUT_DIR / f"{BiophysicalVariables.LAI.value.label}_{stem}.zarr"

#     if output_path.exists():
#         if _scene_has_valid_lai(output_path):
#             print(f"  already saved: {output_path.name}")
#             tile_zarr_paths.append(output_path)
#         else:
#             shutil.rmtree(output_path)
#             print(f"  removed empty: {output_path.name}")
#         continue

#     print(f"  processing: {stem} ...")
#     result = process_scene(item, output_path)
#     if result is not None:
#         tile_zarr_paths.append(output_path)
#         print(f"         -> {output_path.name}")
#     else:
#         print("         [skip] processing failed")

# print(f"\n{len(tile_zarr_paths)} tile(s) saved to {OUTPUT_DIR}/")

💪 Now it is your turn

Task 1: Explore Your Own Area of Interest

Go to http://bboxfinder.com/ and select an area of interest (AOI) (e.g., a research site, and agricultural area, a forest, etc.)
Copy the bounding box coordinates of your area of interest
Change the the coordinates in the AOI_BBOX parameter to your new area of interest.

Task 2: Temporal Analysis

Change the START_DATE_TIME and END_DATE_TIME parameters.
See how the LAI changes for different times of year (e.g., growing season versus the dormant/non-growing season, start of the growing season versus peak growing season).
Repeat for different time spans.

Task 3 (Advanced): Optimize the Calculations

Review the code in the notebook and the biophysical module.
How could the Numpy and SciPy code be refactored into the MAX framework to optimize calculations and leverage an available GPU for faster processing?

Conclusion

This notebook demonstrated an end-to-end workflow for computing Leaf Area Index (LAI) from Sentinel-2 Level-2A data using cloud-native EOPF Zarr archives — without downloading bulky SAFE files.

By combining the EOPF STAC API with lazy Zarr I/O and GPU-accelerated neural network inference, this workflow delivers LAI maps that are numerically equivalent to those produced by SNAP while dramatically reducing both data transfer and processing time. A European growing season of 20–30 Sentinel-2 acquisitions — which would require downloading 12–36 GB and running 20–30 manual SNAP processing steps — can instead be processed in a single unattended run. The declarative pixi.toml environment and CF-compliant Zarr outputs make the results immediately reproducible and cloud-deployment ready, providing a practical template for operationalising any SNAP-based biophysical retrieval at scale.

What’s next?

In the following notebook, we will walk through the implementation of the UN-SPIDER Recommended Practice for SAR flood mapping approach to compare pre-flood and post-flood imagery.
Let’s continue exploring the power of EOPF Zarr.

This notebook was developed as a contribution to the EOPF Zarr Community Notebook Competition.

Created by Geomatys - providing expertise in geodesy, complex data processing and advanced representation modes for over 20 years.

Introduction

What we will learn

Prerequisites

📜 Background

🌿 What is LAI?

🛰️ Origins: The SNAP Sentinel-2 Optical Toolbox

Differences from the Reference Java Implementation

⚙️ Processing Pipeline

🛠️ Setup

Environment Setup with uv

Option A — add as a uv source dependency (recommended and already completed for this project)

Option B — pip install in a notebook cell

⚙️ Processing Pipeline

Setup the constants to be used in the processing

🗂️ The Difference Between Using EOPF Zarr Datasets and SAFE Archives

🔎 Integrating EOPF STAC Queries and Results into a processing Pipeline

Search EOPF STAC for Sentinel-2 Images

🧭 Walkthrough

Build a True-Color Composite

🛰️ Extract Sensor Metadata for Computations

Process the Sun and View Zenith and Azimuth Angles

Sun Angles

View Angles

Relative Azimuth Angle

Trigonometric Angle Features

📚 Combine Extracted Metadata with Observation Data for Neural Network Inference

Assemble the Neural Network Input Array

⚙️ Run Stacked Tiles through LAI Neural Network

Neural Network Architecture

Input Validation (before inference)

Output Validation (after inference)

📘 Save Results to Zarr Format

🗺️ Visualize Results

⚙️⚙️ Batch Processing a Full Growing Season — Per-Tile Output

Batch Processing — Per-Tile Output

Per-Tile Batch Loop

💪 Now it is your turn

Task 1: Explore Your Own Area of Interest

Task 2: Temporal Analysis

Task 3 (Advanced): Optimize the Calculations

Conclusion

What’s next?

Option A — add as a `uv` source dependency (recommended and already completed for this project)

Option B — `pip` install in a notebook cell