fix(s2): write geo metadata for conditions/mask groups so SCL is readable by lhoupert · Pull Request #262 · EOPF-Explorer/data-model

lhoupert · 2026-09-10T11:58:15Z

conditions/mask/l2a_classification had no spatial/proj attrs, no spatial_ref and no grid_mapping, so a GeoZarr reader could not open it — the cause of the SCL_20m failure in EOPF-Explorer/titiler-eopf#163. Geo metadata is now written for any group gridded on x/y rather than for a hardcoded list of paths. Verified by opening a freshly converted store with titiler's GeoZarrReader: the group lists /r20m:scl and /r60m:scl, where before it raised ValueError: not enough values to unpack.

For reviewers: this also changes spatial:bbox and spatial:transform for every geo-stamped group — they were computed from pixel centres and are now pixel edges, a half-pixel shift that propagates to the store-root footprint. Two tests that encoded the old values were updated, and the three structure fixtures were regenerated (the only new nodes are spatial_ref under conditions/mask/**). The pixel-centre premise behind that outset is verified against a shipped store (S2C_MSIL2A_20260427T101021_N0512_R022_T33UWT_20260427T151616_converted.zarr): every geo-stamped group at every level implies the same left edge 499980.00 — measurements/reflectance/r10m x[0]=499985, r60m 500010, r720m 500340, conditions/mask/l2a_classification/r20m 499990. The half-pixel shift is therefore not a convention change: it makes spatial:bbox/spatial:transform agree with coordinates that were already centres. Happy to split that half out if you would rather ship the SCL fix alone.

A second commit here widens the CF standard-name-table except to (OSError, HTTPException) and warns when validation degrades to off — the same block, verbatim, as the s1-tiling branch (#216), which rewrote this except from the other direction. Identical text on both sides means whichever of the two merges second merges cleanly; simulated both orders locally and they produce the same tree.

Author attestation

I am a human, these are my changes, and I have reviewed and understood every change and can explain why each is correct.

AI-assisted: implementation and tests drafted with Claude Code, then reviewed against a converted store and titiler's reader; every behavioural claim above was checked by running it.

🤖 Generated with Claude Code

…able conditions/mask/l2a_classification had no spatial/proj attrs, no spatial_ref and no grid_mapping, so a GeoZarr reader could not open it — the cause of the SCL_20m failure in EOPF-Explorer/titiler-eopf#163. Geo metadata is now written for any group gridded on x/y rather than for a hardcoded list of paths, and spatial:bbox / spatial:transform now describe pixel edges instead of pixel centres, consistent with the spatial:registration: "pixel" these groups declare. Co-Authored-By: Claude Opus 5 <[email protected]> Claude-Session: https://claude.ai/code/session_01JEbpfA8CPcBQJzMpVJBViF

github-actions · 2026-09-10T11:59:01Z

Security Audit Report

View workflow run

Bandit — Static Security Analysis (Security tab)

12 issue(s) found: 12 low

✅ No issues at or above HIGH severity.

12 low issue(s) below threshold not shown in table.

pip-audit — Dependency Vulnerabilities (Security tab)

Package	Version	ID	Fix Versions	Description
tornado	6.5.7	PYSEC-2026-3928	6.5.8	## Summary Tornado parses `application/x-www-form-urlencoded` bodies with `urllib.parse.parse_qs` and does not pass `ma
tornado	6.5.7	GHSA-wwv5-g3v4-889x	6.5.8	## Summary The CVE-2026-35536 fix added a validation loop that rejects `[\x00-\x20\x3b\x7f]`, but only for the hardcoded
tornado	6.5.7	GHSA-8423-8fgw-73vq	6.5.8	## Description ### Summary `parse_multipart_form_data` (httputil.py:34) calls `data.split(b"--"+boundary+b"\r\n")` **b

3 vulnerability/vulnerabilities found (3 fixable) across 1 package(s).

Result: ❌ Blocking issues found — see details above.

lhoupert

@emmanuelmathot I would like your feedback on a few implemenation details

lhoupert · 2026-09-10T14:25:27Z

+    Returns:
+        True if geographic metadata should be written for this group
+    """
+    return {"x", "y"} <= set(dataset.dims)


@emmanuelmathot instead of listing which group paths get CRS metadata, the converter now checks whether a group actually has x/y dimensions. One side effect: /conditions/geometry (the sun and viewing angle grids same UTM grid as the bands) now gets CRS metadata too. It never had any before. Is it ok?

lhoupert · 2026-09-10T14:25:27Z

+            half_y = _half_pixel(y_coords)
+            x_min, x_max = float(x_coords.min()) - half_x, float(x_coords.max()) + half_x
+            y_min, y_max = float(y_coords.min()) - half_y, float(y_coords.max()) + half_y
            spatial_data["spatial:bbox"] = [x_min, y_min, x_max, y_max]


The bounding box here was calculated from the centres of the corner pixels, but the group says its coordinates describe pixel areas, should the box reach the outer pixel edges?

lhoupert · 2026-09-10T14:25:27Z

+    # Coordinates label pixel centres, but an affine transform maps pixel *edges*,
+    # so the origin sits half a pixel outside the first centre. Without this the
+    # transform can never agree with rioxarray's and is half a cell off.
+    x_min = float(x_coords.min()) - pixel_size_x / 2


Same issue as the bounding box. This function built the transform starting from the centre of the first pixel, while rioxarray builds it from the pixel corner.

lhoupert · 2026-09-10T14:25:27Z

-    if "/measurements/" in path or "/quality/" in path:
+    # Add the geo metadata before writing for the groups on a projected grid
+    if needs_geo_metadata(dataset):
        write_geo_metadata(dataset, crs=crs)


Left out on purpose: SCL exists at 20 m and 60 m but the pair is not declared as a pyramid, so clients will not find the 60 m overview but display works without it.

Should I create a separate issue to declare one?

Yes, please open the issue (text below). Also note that conditions/mask/** goes through create_original_encoding. In the fixture, scl/r20m is one 5490×5490 chunk with no sharding, while reflectance is sharded with spatial_chunk. A tile read of SCL must decompress the full array. Raster groups on x/y should use the same chunking and sharding as reflectance. This can go into the same follow-up.

lhoupert · 2026-09-10T14:25:27Z

    CF_STANDARD_NAMES = get_cf_standard_names(url=CF_STANDARD_NAME_URL)
    DO_CF_NAME_VALIDATION = True
-except URLError:
+except (URLError, HTTPException, OSError):


Unrelated to SCL, but it blocked Claude. This code downloads the CF standard-names table as soon as the package is imported, and only handled the connection being refused. If the download is cut off midway a different error escapes, and import eopf_geozarr fails completely, which took down the whole test suite ... :-(

tests/test_cli_e2e.py compares against tests/_test_data/geozarr_examples/, a second golden set distinct from optimized_geozarr_examples/. It was missed when the first set was regenerated, so CI failed on all three products. Via the CLI the converter passes a real CRS down, so /conditions/geometry is now geo-stamped too — visible here and not in the other fixture set. Co-Authored-By: Claude Opus 5 <[email protected]> Claude-Session: https://claude.ai/code/session_01JEbpfA8CPcBQJzMpVJBViF

`URLError`, `TimeoutError`, `ConnectionResetError` and `ssl.SSLError` are all `OSError` subclasses, so `(OSError, HTTPException)` subsumes the previous tuple and covers the rest of the failure modes of an import-time download. Also warn instead of degrading silently: with validation off, `check_standard_name` accepts anything, so a caller that believes it is validating is not. The same block lands verbatim on the `s1-tiling` branch (#216), which rewrote this `except` from the other direction: identical text on both sides means whichever PR merges second merges cleanly instead of conflicting. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01CHpM7QqgsF4jzhPxUxdU7V

`URLError`, `TimeoutError`, `ConnectionResetError` and `ssl.SSLError` are all `OSError` subclasses, so `(OSError, HTTPException)` subsumes the previous tuple and covers the rest of the failure modes of an import-time download. The same block lands verbatim on #262's branch, which rewrote this `except` from the other direction: identical text on both sides means whichever PR merges second merges cleanly instead of conflicting. Co-Authored-By: Claude Opus 5 (1M context) <[email protected]> Claude-Session: https://claude.ai/code/session_01CHpM7QqgsF4jzhPxUxdU7V

@lhoupert

…e different pixel offset definitions, bbox was off by half a pixel, adressed by @lhoupert in EOPF-Explorer#262, impleemted here as well as a merge would be too complciated

emmanuelmathot

Thanks @lhoupert. The SCL fix works, and needs_geo_metadata() is a good change. The half-pixel correction is also correct.

Before merge, please do these steps:

Check the registration of /conditions/geometry (see https://github.com/EOPF-Explorer/data-model/pull/262/changes#r4091279766) before we stamp it with a pixel-edge bbox.
Move the CF standard-name commit (58d4766) to a separate PR. It is not related to SCL, and it is also in #216.
Record the parent bbox limit in the PR description (see https://github.com/EOPF-Explorer/data-model/pull/262/changes#r4091258938). #266 tracks the fix.
Clean up the small items in the inline comments

About your question on a separate issue for the SCL pyramid: yes, I opened #274.

emmanuelmathot · 2026-09-24T08:00:49Z

We need #266 fixed first. See PR #275

emmanuelmathot · 2026-09-24T08:03:24Z

L1156–1171 (needs_geo_metadata)

The CRS on the angle grids is correct, so it is OK to stamp them. But the half-pixel outset is only correct when the coordinates are cell centres. The S2 angle grids are tie-point grids with a 5 km step. I think that their coordinates are grid nodes, starting at the tile corner. If that is true, the outset adds 2.5 km on each side, and spatial:registration: "pixel" is incorrect for these groups.

Can you check this on the shipped store that you used in the description? Compare conditions/geometry x[0] with the tile left edge (499980). If x[0] == 499980, the grid is node-registered. In that case, do one of these:

write spatial:registration: "node" without an outset, or

leave geometry out of this PR and track it separately.

emmanuelmathot · 2026-09-24T08:05:17Z

+def _half_pixel(coords: np.ndarray) -> float:
+    """Half the grid spacing of a coordinate array, or 0.0 when it has no spacing."""
+    if len(coords) < 2:
+        return 0.0
+    return float(np.abs(coords[1] - coords[0])) / 2


The pixel size is now calculated in three places: get_grid_spacing (L59), _transform_from_coordinates (L85–86) and _half_pixel. If we calculate the bbox from the selected transform (see https://github.com/EOPF-Explorer/data-model/pull/262/changes#r4091303319), we do not need this helper, and the bbox always agrees with spatial:transform.

Suggested change

def _half_pixel(coords: np.ndarray) -> float:

"""Half the grid spacing of a coordinate array, or 0.0 when it has no spacing."""

if len(coords) < 2:

return 0.0

return float(np.abs(coords[1] - coords[0])) / 2

def _bbox_from_transform(

transform: tuple[float, float, float, float, float, float], height: int, width: int

) -> list[float]:

"""Outer pixel edges ``[xmin, ymin, xmax, ymax]`` of a north-up grid."""

a, _, c, _, e, f = transform

return [c, f + height * e, c + width * a, f]

emmanuelmathot · 2026-09-24T08:06:23Z

@@ -1199,12 +1235,13 @@ def write_geo_metadata(
            if spatial_transform is not None and not all(t == 0 for t in spatial_transform):
                spatial_data["spatial:transform"] = list(spatial_transform)

-            # Add spatial shape if data variables exist
-            if dataset.data_vars:
-                first_var = next(iter(dataset.data_vars.values()))
-                if first_var.ndim >= 2:
-                    height, width = first_var.shape[-2:]
-                    spatial_data["spatial:shape"] = [height, width]
+            # Shape comes from the raster dims themselves: the first data variable
+            # need not be 2-D, nor have (y, x) as its trailing dims.
+            if "y" in dataset.sizes and "x" in dataset.sizes:
+                spatial_data["spatial:shape"] = [
+                    int(dataset.sizes["y"]),
+                    int(dataset.sizes["x"]),
+                ]


This block calculates the bbox from the coordinates and calculates the transform separately. The two values can go out of sync. Please calculate the bbox from the transform and the shape. If there is no transform (fewer than 2 pixels), keep the old behaviour.

Suggested change

if "x" in dataset.coords and "y" in dataset.coords:

height, width = int(dataset.sizes["y"]), int(dataset.sizes["x"])

spatial_data["spatial:shape"] = [height, width]

spatial_transform = _preferred_spatial_transform(dataset)

if spatial_transform is not None and not all(t == 0 for t in spatial_transform):

spatial_data["spatial:transform"] = list(spatial_transform)

spatial_data["spatial:bbox"] = _bbox_from_transform(

spatial_transform, height, width

)

else:

x_coords = dataset.coords["x"].values

y_coords = dataset.coords["y"].values

spatial_data["spatial:bbox"] = [

float(x_coords.min()),

float(y_coords.min()),

float(x_coords.max()),

float(y_coords.max()),

]

emmanuelmathot · 2026-09-24T08:07:18Z

+        if "spatial_ref" in dataset.coords and "spatial_ref" not in encoding:
+            encoding["spatial_ref"] = {"compressors": None}


This line writes into the encoding dict that the caller gave. It is a hidden side effect. Please make a copy.

Suggested change

if "spatial_ref" in dataset.coords and "spatial_ref" not in encoding:

encoding["spatial_ref"] = {"compressors": None}

if "spatial_ref" in dataset.coords and "spatial_ref" not in encoding:

encoding = {**encoding, "spatial_ref": {"compressors": None}}

emmanuelmathot · 2026-09-24T08:09:30Z

+        y_dim, x_dim = dims
+        return xr.Dataset(
+            {"var": ([y_dim, x_dim], np.zeros((4, 4)))},
+            coords={y_dim: (([y_dim]), np.arange(4.0)), x_dim: (([x_dim]), np.arange(4.0))},


Suggested change

coords={y_dim: (([y_dim]), np.arange(4.0)), x_dim: (([x_dim]), np.arange(4.0))},

coords={y_dim: np.arange(4.0), x_dim: np.arange(4.0)},

emmanuelmathot · 2026-09-24T08:09:59Z

+        assert written_ds.rio.crs is not None
+        assert written_ds.rio.crs.to_epsg() == 32632
+        assert "spatial_ref" in written_ds.coords
+        assert "proj:code" in written_ds.attrs or "proj:wkt2" in written_ds.attrs


This assertion passes for any CRS. Please assert the expected value, for example written_ds.attrs["proj:code"] == "EPSG:32632" (or the WKT2 equivalent, if that is what build_convention_attrs writes).

emmanuelmathot · 2026-09-24T08:10:17Z

        assert "zarr_conventions" not in sample_dataset_no_crs.attrs
+
+
+class TestNeedsGeoMetadata:


test_classification_mask_is_readable_after_write and test_lat_lon_group_is_left_alone test stream_write_dataset from end to end. They do not test needs_geo_metadata. Please move them to a separate class, for example TestStreamWriteGeoMetadata, so that a failure points to the correct function.

emmanuelmathot · 2026-09-24T08:10:51Z

+            # `spatial:registration` below declares "pixel", so the bbox covers the
+            # pixel edges. Coordinates are centres, hence the half-pixel outset —
+            # without it the footprint is a half-pixel narrower than the raster and
+            # disagrees with the arrays' own `proj:bbox`.


Please make these comments shorter. One line is sufficient, for example: # Coordinates are pixel centres; the bbox covers pixel edges. The reason for the change belongs in the PR description, not in the code.

emmanuelmathot · 2026-09-24T08:11:10Z

    pixel_size_y = float(np.abs(y_coords[1] - y_coords[0]))
-    x_min = float(x_coords.min())
-    y_max = float(y_coords.max())
+    # Coordinates label pixel centres, but an affine transform maps pixel *edges*,


Please make these comments shorter. One line is sufficient, for example: # Coordinates are pixel centres; the bbox covers pixel edges. The reason for the change belongs in the PR description, not in the code.

lhoupert · 2026-09-29T09:53:14Z

Superseded : #277 made SCL readable, #275 has the bbox change, and #285 removes the CF check.

lhoupert marked this pull request as draft September 10, 2026 12:28

lhoupert mentioned this pull request Sep 10, 2026

fix(stac): point S2 AOT/WVP assets at the store root (titiler-eopf#163) EOPF-Explorer/data-pipeline#412

Merged

1 task done

lhoupert commented Sep 10, 2026

View reviewed changes

lhoupert force-pushed the fix--conditions-mask-geo-metadata branch from 529e040 to dd77bbb Compare September 10, 2026 14:40

lhoupert marked this pull request as ready for review September 10, 2026 14:44

lhoupert requested a review from emmanuelmathot September 10, 2026 14:44

shollendEODC mentioned this pull request Sep 15, 2026

Fix correct bounds for spatial:bbox attributes due to mismatch/offset between coords and bounds #266

Closed

lhoupert mentioned this pull request Sep 22, 2026

No way to tell which assets are renderable: AOT/SCL/WVP 500 on one route, empty 200 on the other EOPF-Explorer/titiler-eopf#163

Closed

emmanuelmathot requested changes Sep 24, 2026

View reviewed changes

This was referenced Sep 24, 2026

fix(s2): compute spatial:bbox and spatial:transform from pixel edges #275

Merged

feat!: reintegrate cpm-geozarr into main #277

Merged

Extension of functionalities to all other modes #248

Open

lhoupert closed this Sep 29, 2026

-def _half_pixel(coords: np.ndarray) -> float:
-    """Half the grid spacing of a coordinate array, or 0.0 when it has no spacing."""
-    if len(coords) < 2:
-        return 0.0
-    return float(np.abs(coords[1] - coords[0])) / 2
+def _bbox_from_transform(
+    transform: tuple[float, float, float, float, float, float], height: int, width: int
+) -> list[float]:
+    """Outer pixel edges ``[xmin, ymin, xmax, ymax]`` of a north-up grid."""
+    a, _, c, _, e, f = transform
+    return [c, f + height * e, c + width * a, f]

+        if "x" in dataset.coords and "y" in dataset.coords:
+            height, width = int(dataset.sizes["y"]), int(dataset.sizes["x"])
+            spatial_data["spatial:shape"] = [height, width]
+            spatial_transform = _preferred_spatial_transform(dataset)
+            if spatial_transform is not None and not all(t == 0 for t in spatial_transform):
+                spatial_data["spatial:transform"] = list(spatial_transform)
+                spatial_data["spatial:bbox"] = _bbox_from_transform(
+                    spatial_transform, height, width
+                )
+            else:
+                x_coords = dataset.coords["x"].values
+                y_coords = dataset.coords["y"].values
+                spatial_data["spatial:bbox"] = [
+                    float(x_coords.min()),
+                    float(y_coords.min()),
+                    float(x_coords.max()),
+                    float(y_coords.max()),
+                ]

		if "spatial_ref" in dataset.coords and "spatial_ref" not in encoding:
		encoding["spatial_ref"] = {"compressors": None}

	coords={y_dim: (([y_dim]), np.arange(4.0)), x_dim: (([x_dim]), np.arange(4.0))},
	coords={y_dim: np.arange(4.0), x_dim: np.arange(4.0)},

		assert "zarr_conventions" not in sample_dataset_no_crs.attrs


		class TestNeedsGeoMetadata:

Conversation

lhoupert commented Sep 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Author attestation

Uh oh!

github-actions Bot commented Sep 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Security Audit Report

Bandit — Static Security Analysis (Security tab)

pip-audit — Dependency Vulnerabilities (Security tab)

Uh oh!

lhoupert left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

lhoupert Sep 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

lhoupert Sep 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

lhoupert Sep 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

lhoupert Sep 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

lhoupert Sep 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

emmanuelmathot left a comment

Choose a reason for hiding this comment

Uh oh!

emmanuelmathot Sep 24, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

lhoupert commented Sep 29, 2026

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

lhoupert commented Sep 10, 2026 •

edited

Loading

github-actions Bot commented Sep 10, 2026 •

edited

Loading

lhoupert left a comment •

edited

Loading

lhoupert Sep 10, 2026 •

edited

Loading

lhoupert Sep 10, 2026 •

edited

Loading

lhoupert Sep 10, 2026 •

edited

Loading

lhoupert Sep 10, 2026 •

edited

Loading

lhoupert Sep 10, 2026 •

edited

Loading

emmanuelmathot Sep 24, 2026 •

edited

Loading