Cpm geozarr eodc changes - #247

is there documentation somewhere for this attribute? is it part of a published specification?

no this was implemented test wise to reflect an internal discussion at eodc at an eopf meeting. During the L1C multiscale calculation not all bands are present at 60m, as they are only featured at their native resolution. As we need all bands for the multiscales they have to be provided at a uniform resolution of 60m before being resampled to 120m, ...

Thus, we need to do the resampling during processing (for which the inject_bands functionality is already there! nice!). But we want to differentiate between native resolution bands and our resampled ones, which is why we decided to add additional metadata in the attributes. These are jsut there to inform users.. Naming and content can be adapted

@shollendEODC It is good to be explicit, but I would suggest to keep it more concise and not overwhelm users with too much metadata. Users should decide themselves how they want to handle it. Often you have some compact metadata attributes and then a detailed explanation in a product user manual (here it does not exist as far as I know, but maybe they need to create a new one anyway for explaining the GeoZarr structure for each product). I would suggest to remove it, since you have anyway the '_resampling_mode'.

I removed '_processing_mode' but left the other parameters

@shollendEODC It is good to be explicit, but I would suggest to keep it more concise and not overwhelm users with too much metadata. Users should decide themselves how they want to handle it. Often you have some compact metadata attributes and then a detailed explanation in a product user manual (here it does not exist as far as I know, but maybe they need to create a new one anyway for explaining the GeoZarr structure for each product). I would suggest to remove it, since you have anyway the '_resampling_mode'.

Would be nice to put this into an Enum, since then you are less prone to typos and it is easier to compare sth to it.

put it in an enum (for now in s2_multiscale.py) but would move to own utils module once the other modes are implemented to have a unified enum somewhere.
its also likely a similar designation exists within cpm, didnt check if it is accessible but for now its alrigt

Outdated description? The maximum is taken, not the first (which will be always the first)?

What happens if a 2D metadata array has a size larger than the spatial_chunk? Because the rule here seems to be that if the chunking is not set and it is not 1D, there is one chunk for the whole dataset. Why doing it only for the 1D case?

logic to calculate zarr encoding chukns moved to its own function get_chunking_for_encoding to differentiate between xarray chunksizes (a tuple of all chunks: ((1024, 806), (1024, 806)) and zarr chunks (tuple of maximal chunk (1024, 1024).

shards would be a better name to be aligned with chunks.

Please remove your personal remarks/comments to have less noise in the code.

Try to put your personal things either to an already existing ignored folder or at least into one single folder. Otherwise .gitignore gets spammed if every contributor has special folders. And .vscode should be ignored at a higher level, maybe you could uncomment line 192?

should be done

But why not taking the max for each dimension?

(spatial_chunk, spatial_chunk)

I would raise a warning here, because this should never happen if both metadata and measurement arrays are rechunked before?

implemented an Error, as it would lead to an error later on during zarr writing anyway

Just a bit confused here why there is this clause - shouldn't we get a 1 element list from the call before: max_chunksizes = [...], which we just can convert to a tuple? Could you double-check please?

yeah, it needs to be a tuple for encoding[chunks] = (chuknsize, ) but a tuple cast works nicer and removes the if clause, many thanks

tin .vscode, there is a useful launch.json that we ahre to make E2E integration tests

attrs["fill_value"] resolves here - the fixtures this PR adds carry fill_value in the array attributes of measurements/reflectance/* (e.g. attrs: {add_offset, dtype, fill_value: 0, long_name, scale_factor, ...}), as do real CPM 3.0.0 products. So the masking does what it intends for the variables this branch covers.

Two things worth considering anyway.

1. The lookup is more brittle than it looks. attrs["fill_value"] is a CPM-specific attribute - not CF _FillValue, and not the zarr-level fill_value. Where a zarr v3 fill value actually surfaces once xarray opens the array:

For the uint16 reflectance case, neither the CF encoding key nor the zarr-level fill value is reachable through xarray - the CPM attribute is the only thing that works. Fine today, but if CPM stops writing that attribute, or sanitize_array_attrs strips it, this branch silently reverts to unmasked means, with no error and no failing test. A fallback chain (attrs["fill_value"] -> encoding["_FillValue"] -> the zarr array's own fill_value) removes that failure mode cheaply. Prior art on the same confusion: cogeotiff/rio-tiler#914 (comment)

2. This branch has no test. _coarsen_variable has dedicated tests for the classification and quality_mask branches, but not for the reflectance masking path - which is the actual fix for #246. Worth adding one: an array with fill-valued border pixels, asserting the coarsened mean ignores them and the fill value is restored afterwards. That would also catch point 1 if the attribute ever disappears.

will add a fallback to encoding["_FillValue"]

These three attributes should be removed. The multiscales convention already declares all of them, and duplicating them here creates a second source of truth that can drift.

Mapping them one for one onto ScaleLevel in data_api/geozarr/multiscales/zcm.py:

The first two are already written for every level by add_multiscales_metadata_to_parent() - see the scale_level_data["derived_from"] and scale_level_data["transform"] assignments around line 1100. So _derived_from and _factor are pure duplication today.

resampling_method is the one genuinely missing piece: the model supports it but we never populate it. If the intent is to record that overviews are mean-resampled, the right fix is to set ScaleLevel.resampling_method there, not to stamp a private attribute on each array.

these are purposely filled in with a '_' before them to highlight that they are not real pyramid overviews but ingested bands. Open for discussion if this nomenclature makes sense or not, but I think its important to highlight it in some way

They are kept for now, mutliscale metadata is not written for these injected bands, so some sort of metadata should be added to differentiate between the native resolution and injected bands.

Multiscale metdata are written via the multiscale zarr convention at the root of the pyramid group. This is an unecessary duplication. See https://github.com/shollendEODC/data-model/blob/b629cd284d5ed911b59cd88144b1e3600eefe7be/src/eopf_geozarr/s2_optimization/s2_multiscale.py#L957

Alright, I agree. Will be removed.

as per our discussion, these attributes are removed. Some other form of attribute/annotation will be added but this will be shifted to a new PR

For context, this branch is copy-pasted from create_measurements_encoding, where the same code exists without this note - so it isn't a regression introduced here.

But --experimental-scale-offset-codec is a real opt-in flag on both the CLI and the CPM writer, so a shipped feature currently carries an in-code note saying it doesn't work. Could you clarify which it is? If it works now, drop the line. If it doesn't, better to raise an issue and have the flag fail loudly than leave the note in the source.

yeah i will test it again and remove the code if it fails

Not blocking this PR, but flagging it here so it is on record.

This copies the input product's root attributes through verbatim, which in practice means stac_discovery lands in the output unchanged. That object describes the input layout. EODC's 31 August product (built from this branch) shows the effect: the output root carries stac_discovery with 83 assets, every href pointing at input paths like /conditions/geometry/sun_angles, and not one of them references the pyramid levels this conversion just created (r120m, r360m, r720m).

So the store ships a STAC object that no longer matches the store it lives in, and a consumer has no way to tell it is stale. Arguably worse than not copying it at all.

Two smaller points on the mechanics:

• The function is a generic root.attrs.update(...) with a STAC-specific name, and it re-opens the store in r+ immediately after write_store_root_bbox() closed it, duplicating the get_storage_options boilerplate. Between this, write_store_root_geo_metadata and write_store_root_bbox (itself already a thin wrapper), there are now three root-writing entry points.
• Because the update() is unconditional and runs after the minispec writer, any key collision silently wins. I checked: current CPM inputs carry only other_metadata, processing_history and stac_discovery at the root, so nothing collides today. But we are actively asking CPM to start declaring zarr_conventions at the root (https://gitlab.eopf.copernicus.eu/cpm/eopf-cpm/-/issues/1189), which would turn this into a silent overwrite of the computed spatial:bbox / proj:code footprint. Worth making the merge explicit about which side wins.

The broader question of what unspecified CPM metadata we should carry through at all is out of scope here and tracked in #257.

will include:

• the omitted pyramid levels in the stac metadata
• key checks for update()

will not consolidate root writing scripts now -> can be done in general refactoring for inclusion of all modes

Pyramid levels will NOT be added to the STAC metadata to avoid its cluttering and prevent users from accessing them. Functionality remains but is not used.

sure this is understandable. Nevertheless, the root pyramid asset "reflectance" should be referenced. Multiscale consumers expect to discover the multiscale attributes in the reference group.

To clarify: Do you mean the exposure in the root/zarr.json under

"attributes"/"stac_discovery"/"assets"/ ->
        "b01_120m": {
          "href": "/measurements/reflectance/r120m/b01",
          "title": "b01_120m"
        },

Not that, no - those are per-band, per-level entries (b01_120m, b02_120m, ...), which is exactly the clutter you want to avoid. I agree with avoiding them.

What I mean is a single asset pointing at the pyramid group, shaped like a real STAC asset. The reference is what our own catalogue already publishes for this product - see the reflectance asset in this item:

"reflectance": {
  "href": ".../measurements/reflectance",
  "type": "application/vnd.zarr; version=3; profile=multiscales",
  "title": "Surface Reflectance",
  "roles": ["data", "reflectance"],
  "gsd": 10,
  "proj:code": "EPSG:32626",
  "proj:shape": [10980, 10980],
  "bands": [ ... ],
  "cube:dimensions": { ... },
  "cube:variables": { ... }
}

The profile=multiscales media-type parameter is the important part: it tells a consumer this asset is a pyramid root, so it opens that one node, reads multiscales, and resolves r10m ... r720m from the convention. No per-level assets needed - that is the point of the convention.

The gap today: the input carries 25 reflectance assets, all pointing at leaf arrays (/measurements/reflectance/r20m/b03), and none at /measurements/reflectance. A consumer has no entry point to the pyramid at all.

Suggested minimum below - href, type, title, roles, gsd, plus proj:code / proj:shape, all available at this call site. I have deliberately left out bands, cube:dimensions and cube:variables: they are buildable from BAND_INFO, but that needs eo:common_name, descriptions and a nm-to-micrometre conversion which do not exist in the repo yet, so better as a follow-up than scope creep here.

This also lets add_multiscale_pyramids_to_stac_metadata() be deleted rather than left commented out.

One note on hrefs: every asset in this object uses root-absolute paths (/measurements/...), so I kept that for consistency. The STAC Zarr Convention v0.1 requires hrefs relative to the group carrying the object (measurements/...), but converting one asset alone would make the object inconsistent - that belongs with the stac:item migration in https://gitlab.eopf.copernicus.eu/cpm/eopf-cpm/-/issues/1189.

Not that, those are per-band, per-level entries (b01_120m, b02_120m, …), which is exactly the clutter you want to avoid. I agree with avoiding them.

What I mean is a single asset pointing at the pyramid group, shaped like a real STAC asset. The reference is what our own catalogue already publishes for this product. See the reflectance asset in this item:

"reflectance": {
  "href": ".../measurements/reflectance",
  "type": "application/vnd.zarr; version=3; profile=multiscales",
  "title": "Surface Reflectance",
  "roles": ["data", "reflectance"],
  "gsd": 10,
  "proj:code": "EPSG:32626",
  "proj:shape": [10980, 10980],
  "bands": [ ... ],
  "cube:dimensions": { ... },
  "cube:variables": { ... }
}

The profile=multiscales media-type parameter is the important part: it tells a consumer this asset is a pyramid root, so it opens that one node, reads multiscales, and discovers r10m … r720m from the convention. No per-level assets needed, that is the point of the convention.

The gap today: the input carries 25 reflectance assets, all pointing at leaf arrays (/measurements/reflectance/r20m/b03), and none at /measurements/reflectance. A consumer has no entry point to the pyramid at all.

Suggested minimum below (href, type, title, roles, plus proj:code / proj:shape / gsd), which are all available at this call site. I have deliberately left out bands, cube:dimensions and cube:variables: they are buildable from BAND_INFO, but that needs eo:common_name, descriptions and a nm→µm conversion which does not exist in the repo yet, so it is better as a follow-up than as scope creep here.

This also lets add_multiscale_pyramids_to_stac_metadata() be deleted rather than left commented out.

okay, I will take a look and let you know

As discussed, this PR will be merged without this feature, but will be implemented in a new PR

Fine as a temporary measure - the docstring is honest about why it exists, and hand-rolling this beats guessing at a V3 model under time pressure.

Worth recording though: this is a schema expressed as imperative code. Sentinel2Root is V2-only, and pyz/v3.py already provides GroupSpec / ArraySpec over pydantic_zarr.v3 but nothing in src/ imports it yet - so the V3 scaffolding exists, the models just haven't been ported.

Tracked in #258 to fold this back into the pydantic model rather than leave detection split between a model and a function.