fix(fci): bounded FDHSI entry selection

docs/fci-eumetsat.md

@@ -85,13 +85,36 @@ No rectangular query bbox is presented as the provider footprint.
 ## Chunk downloads and native grids
 
 EUMDAC supplies `AccessToken`, `DataStore.get_collection().search()`,
-`Product.entries` and `Product.open(entry=...)`. L1c fetch inspects HDF5
-metadata and coordinate vectors with EUMDAC byte-range requests across body
-entries, then downloads the complete entries containing requested-channel
-pixel centres inside the bbox, plus the trailer. Chunk heights can vary;
-selection uses their published scan coordinates. Pixel arrays in other body
-entries are not decoded for selection. A server that does not honor byte ranges
-raises an error. This transport behavior still needs a credentialed live run.
+`Product.entries` and `Product.open(entry=...)`. On a cold L1c layout, fetch
+reads the provider's HDF5 projection and requested-channel x/y coordinate
+vectors across BODY entries. It downloads complete entries containing
+requested-channel pixel centres inside the bbox, plus the trailer. Chunk
+heights can vary; selection uses measured coordinates rather than a table
+of inferred chunk rows. Pixel arrays are not decoded for selection.
+
+Reuse the same `FCIConnector` for successive repeat cycles. It keeps up to eight
+measured layouts in memory, partitioned by collection, published product and
+processor versions, platform/product/encoding/processing modes, exact BODY
+chunk-number set (including count), and requested channel set. A matching later
+cycle evaluates its bbox against those grids with **zero metadata range calls**;
+whole-entry downloads and catalogue requests still occur. Different bboxes
+reuse the coordinates. Changing channels measures a separate layout. Missing
+version/identity metadata or a changed version/chunk set triggers measured
+selection for that product. No cache survives a connector/process restart.
+
+`FCIResult.entry_selection`, Dataset attrs and canonical `entry_selection`
+record `mode` (`measured` or `cached`), fallback `reason`, BODY count, worker
+limit and coordinate basis. Reuse assumes provider grids remain stable within
+that versioned identity. The synthetic fixtures prove the cache mechanism;
+real-header stability and the 25-cycle/30-minute cloud budget still require
+RefCal's cloud-side proof. No hardcoded PUG chunk layout is used.
+
+Cold measurements fetch only demanded 4 KiB metadata pages (plus 64-byte
+superblocks), using a bounded pool: `FCIConnector(metadata_workers=4)` by
+default. Local HDF5 inspection pauses on a missing page; network reads run
+outside HDF5's global lock so separate BODY entries can progress concurrently.
+A server that ignores byte ranges raises an error. This transport behavior
+still needs a credentialed cloud run.
 L2 products are single NetCDF entries and are downloaded before local clipping.
 Short cache filenames avoid Windows NetCDF path limits; result metadata retains
 the original provider entry names.

src/spectraccess/connectors/fci_eumetsat/connector.py

@@ -2,10 +2,14 @@
 from __future__ import annotations
 
 import io
+import gzip
+import json
 import itertools
 import os
 import re
 import shutil
+from collections import OrderedDict
+from concurrent.futures import ThreadPoolExecutor
 from copy import deepcopy
 from contextlib import contextmanager
 from dataclasses import dataclass, field
@@ -117,6 +121,7 @@ class FCIResult:
     bbox: tuple[float, float, float, float]
     entries: tuple[str, ...]
     local_files: Mapping[str, str] = field(default_factory=dict)
+    entry_selection: Mapping[str, Any] = field(default_factory=dict)
 
 
 def _path(result, entry):
@@ -137,16 +142,24 @@ def _provider_failure(exc, stage, collection=None, *, values=()):
                                     provider_error=FCIProviderError, authorization_error=FCIAuthorizationError)
 
 
+class _RangeNeeded(Exception):
+    """Suspend HDF5 inspection until a missing metadata page is fetched."""
+
+    def __init__(self, stream, block):
+        self.stream, self.block = stream, block
+
+
 class _EntryRange(io.RawIOBase):
     """Seekable HDF5 metadata access using EUMDAC's maintained byte-range API.
 
     The HDF5 superblock publishes the EOF address. Small blocks are cached;
     pixel arrays are never accessed while choosing entries. A server that
     ignores Range is rejected rather than downloading a complete repeat cycle.
     """
-    block_size = 65536
+    block_size = 4096
 
-    def __init__(self, product, entry):
+    def __init__(self, product, entry, *, deferred=False):
+        self.deferred = deferred
         self.product, self.entry, self.position = product, entry, 0
         self.blocks = {}
         header = self._range(0, 64)
@@ -193,12 +206,18 @@ def read(self, size=-1):
             block = self.position // self.block_size
             start = block * self.block_size
             if block not in self.blocks:
-                self.blocks[block] = self._range(start, min(start + self.block_size, self.length))
+                if self.deferred:
+                    raise _RangeNeeded(self, block)
+                self.load_block(block)
             stop = min(end, start + len(self.blocks[block]))
             parts.append(self.blocks[block][self.position - start:stop - start])
             self.position = stop
         return b"".join(parts)
 
+    def load_block(self, block):
+        start = block * self.block_size
+        self.blocks[block] = self._range(start, min(start + self.block_size, self.length))
+
     def readinto(self, buffer):
         data = self.read(len(buffer))
         buffer[:len(data)] = data
@@ -240,20 +259,100 @@ def _window(x, y, projection, bbox):
     return slice(rows[0], rows[-1] + 1), slice(min(cols), max(cols) + 1)
 
 
-def _entry_intersects(product, entry, bbox, channels):
-    with _EntryRange(product, entry) as stream, h5netcdf.File(stream, "r") as nc:
-        data = nc.groups["data"]
-        projection = dict(data.variables["mtg_geos_projection"].attrs)
+def _entry_geometry(stream, channels):
+    # Retain measured vectors, including variable chunk heights and channel grids.
+    # Never infer chunk rows from their ordinal number.
+    stream.seek(0)
+    with h5py.File(stream, "r") as nc:
+        data = nc["data"]
+        projection = {k: np.asarray(v).item() if np.asarray(v).size == 1 else v
+                      for k, v in data["mtg_geos_projection"].attrs.items()}
+        grids = {}
         for channel in channels:
-            if channel not in data.groups:
-                continue
-            group = data.groups[channel].groups["measured"]
-            rows, cols = _window(_angles(group.variables["x"]), _angles(group.variables["y"]), projection, bbox)
-            if rows.stop and cols.stop:
-                return True
+            if channel in data:
+                group = data[channel]["measured"]
+                grids[channel] = (_angles(group["x"]), _angles(group["y"]))
+        return projection, grids
+
+
+def _geometry_intersects(geometry, bbox):
+    projection, grids = geometry
+    for x, y in grids.values():
+        rows, cols = _window(x, y, projection, bbox)
+        if rows.stop and cols.stop:
+            return True
     return False
 
 
+def _entry_intersects(product, entry, bbox, channels):
+    with _EntryRange(product, entry) as stream:
+        return _geometry_intersects(_entry_geometry(stream, channels), bbox)
+
+
+def _measure_entries(product, entries, channels, token, workers):
+    """Read metadata pages concurrently, without holding HDF5's global lock.
+
+    HDF5 inspection is local and serial. A cache miss unwinds the file handle;
+    the pool fetches one missing page per entry, then inspection resumes from
+    cached pages. This avoids serial network waits inside HDF5's Python lock.
+    Auth context must be propagated explicitly to the transport threads.
+    """
+    def create(entry):
+        with _eumdac_context(token=token):
+            return _EntryRange(product, entry, deferred=True)
+
+    def load(needed):
+        with _eumdac_context(token=token):
+            needed.stream.load_block(needed.block)
+
+    geometry = {}
+    with ThreadPoolExecutor(max_workers=workers) as pool:
+        streams = list(pool.map(create, entries))
+        try:
+            pending = list(zip(entries, streams))
+            while pending:
+                missing, remaining = [], []
+                for entry, stream in pending:
+                    try:
+                        geometry[entry] = _entry_geometry(stream, channels)
+                    except _RangeNeeded as needed:
+                        missing.append(needed)
+                        remaining.append((entry, stream))
+                list(pool.map(load, missing))
+                pending = remaining
+        finally:
+            for stream in streams:
+                stream.close()
+    return geometry
+
+
+def _layout_identity(target, entries, channels):
+    info = target.raw.get("download_properties", {}).get("productInformation", {})
+    if not info:
+        info = target.raw.get("properties", {}).get("productInformation", {})
+    version = info.get("productVersion", info.get("version"))
+    processing = info.get("processingInformation", {}).get("processorVersion", info.get("processingVersion"))
+    identities = []
+    for entry in entries:
+        match = _SUFFIX.search(entry)
+        if match is None or "-FDHSI-FD-" not in entry:
+            return None
+        # Platform, product kind, distribution format and processing modes also
+        # partition layouts. Timestamps and repeat-cycle number do not.
+        identities.append((entry.split("-CHK-")[0], entry.split("-BODY-")[1].split("_")[0],
+                           tuple(entry[:match.start()].split("_")[-2:]),
+                           match.groups()[2:5], int(match[7])))
+    if not identities or len(set(identities)) != len(identities):
+        return None
+    if not version or not processing:
+        checksum = target.raw.get("properties", {}).get("md5")
+        # Incomplete processing metadata cannot establish cross-cycle identity.
+        # A published archive checksum can still bind this exact product.
+        return ((target.collection, "exact-product-md5", target.product_id, str(checksum),
+                 tuple(sorted(identities)), tuple(sorted(channels))) if checksum else None)
+    return (target.collection, str(version), str(processing), tuple(sorted(identities)), tuple(sorted(channels)))
+
+
 def _file_info(path):
     match = _SUFFIX.search(path.name)
     if match is None or "-FDHSI-FD-" not in path.name:
@@ -316,8 +415,61 @@ class FCIConnector(Connector):
     """EUMDAC discovery and chunk fetch, Satpy values, native-window Dataset."""
     credential_provider = "eumetsat"
 
-    def __init__(self, *, credentials: CredentialSource | None = None):
+    def __init__(self, *, credentials: CredentialSource | None = None, metadata_workers=4):
+        if not isinstance(metadata_workers, int) or metadata_workers < 1:
+            raise ValueError("metadata_workers must be a positive integer")
         self.credentials = credentials
+        self.metadata_workers = metadata_workers
+        self._entry_layouts = OrderedDict()
+
+    def export_layout_cache(self):
+        """Portable measured coordinates; contains no credentials or pixels."""
+        rows = [[key, [[projection, {channel: [x.tolist(), y.tolist()]
+                  for channel, (x, y) in grids.items()}] for projection, grids in layout]]
+                for key, layout in self._entry_layouts.items()]
+        return gzip.compress(json.dumps({"schema": "fci-entry-layout/1", "layouts": rows}).encode())
+
+    def import_layout_cache(self, data):
+        """Restore layouts; normal version/chunk identity checks still apply."""
+        def frozen(value):
+            return tuple(frozen(v) for v in value) if isinstance(value, list) else value
+        body = json.loads(gzip.decompress(data))
+        if body.get("schema") != "fci-entry-layout/1":
+            raise ValueError("unsupported FCI entry layout cache")
+        layouts = OrderedDict()
+        for key, layout in body["layouts"][-8:]:
+            layouts[frozen(key)] = tuple((projection, {channel: (np.asarray(x), np.asarray(y))
+                for channel, (x, y) in grids.items()}) for projection, grids in layout)
+        self._entry_layouts.update(layouts)
+        while len(self._entry_layouts) > 8:
+            self._entry_layouts.popitem(last=False)
+
+    def _select_bodies(self, product, target, entries, channels, bbox, token):
+        key = _layout_identity(target, entries, channels)
+        cached = self._entry_layouts.get(key) if key is not None else None
+        reason = None
+        if cached is not None:
+            self._entry_layouts.move_to_end(key)
+            mode = "cached"
+        else:
+            reason = ("missing_version_or_chunk_identity" if key is None else
+                      "version_or_chunk_set_mismatch" if self._entry_layouts else "cold_cache")
+            geometry = _measure_entries(product, entries, channels, token, self.metadata_workers)
+            mode = "measured"
+            if key is not None:
+                # Store by chunk identity rather than provider entry order.
+                ordered = sorted(entries, key=lambda e: int(_SUFFIX.search(e)[7]))
+                self._entry_layouts[key] = tuple(geometry[e] for e in ordered)
+                while len(self._entry_layouts) > 8:
+                    self._entry_layouts.popitem(last=False)
+        if cached is not None:
+            ordered = sorted(entries, key=lambda e: int(_SUFFIX.search(e)[7]))
+            geometry = dict(zip(ordered, cached))
+        selected = {e for e in entries if _geometry_intersects(geometry[e], bbox)}
+        return selected, {"mode": mode, "reason": reason, "body_count": len(entries),
+                          "metadata_workers": self.metadata_workers,
+                          "basis": "measured_provider_coordinates"}
+
 
     def _store(self, credentials=None):
         credential = resolve("eumetsat", credentials if credentials is not None else self.credentials)
@@ -366,19 +518,24 @@ def fetch(self, target, *, dest=None, bbox=None, channels=CHANNELS, credentials=
         import hashlib
         path = path / hashlib.sha256((target.collection + target.product_id).encode()).hexdigest()[:20]
         path.mkdir(parents=True, exist_ok=True)
-        downloaded, local_files = [], {}
+        downloaded, local_files, selection = [], {}, {}
         with _eumdac_context(token=store.token):
             try:
                 product = store.get_product(target.collection, target.product_id)
                 available = tuple(product.entries)
                 l1c = target.collection == COLLECTIONS["l1c"]
                 selected = []
+                bodies = tuple(e for e in available if e.endswith(".nc") and "-BODY-" in e)
+                selected_bodies = set()
+                if l1c:
+                    selected_bodies, selection = self._select_bodies(
+                        product, target, bodies, channels, bbox, store.token)
                 for entry in available:
                     if not entry.endswith(".nc"):
                         continue
                     if not l1c or "-TRAIL-" in entry:
                         selected.append(entry)
-                    elif "-BODY-" in entry and _entry_intersects(product, entry, bbox, channels):
+                    elif entry in selected_bodies:
                         selected.append(entry)
                 if l1c and any("-BODY-" in e for e in selected) and not any("-TRAIL-" in e for e in selected):
                     raise ValueError("FCI repeat cycle has no trailer entry")
@@ -399,7 +556,7 @@ def fetch(self, target, *, dest=None, bbox=None, channels=CHANNELS, credentials=
                     local_files[filename] = local_name
             except Exception as exc:
                 raise _provider_failure(exc, "entry fetch", target.collection, values=_auth_values(token=store.token)) from exc
-        return FCIResult(path, target, datetime.now(timezone.utc), bbox, tuple(downloaded), local_files)
+        return FCIResult(path, target, datetime.now(timezone.utc), bbox, tuple(downloaded), local_files, selection)
 
     def parse(self, result, **kwargs):
         return self.parse_canonical(result, **kwargs)
@@ -411,6 +568,7 @@ def parse_canonical(self, result, *, bbox=None, channels=CHANNELS):
         frame = target_to_canonical(result.target, retrieved_at=result.retrieved_at,
                                     provider_attrs=ds.attrs.get("provider_attrs", {}),
                                     qa=ds.attrs.get("provider_product_qa"))
+        frame["entry_selection"] = [dict(result.entry_selection)]
         if any("radiometric_noise_lut_noise" in name for name in ds):
             frame["unc_definition"] = NOISE_DEFINITION
         return frame
@@ -427,7 +585,7 @@ def read(self, result, *, bbox=None, channels=CHANNELS):
                              integration_start=result.target.start.isoformat(),
                              integration_end=result.target.end.isoformat(),
                              fetched_bbox=result.bbox, requested_bbox=bbox, source_entries=result.entries,
-                             provider_coverage=result.target.coverage)
+                             provider_coverage=result.target.coverage, entry_selection=dict(result.entry_selection))
         if result.target.footprint:
             dataset.attrs["footprint_geometry"] = result.target.footprint
         return dataset
@@ -544,7 +702,7 @@ def _read_l1(result, bbox, channels):
     output = xr.Dataset()
     for channel, chunks in pieces.items():
         if chunks:
-            merged = xr.concat(chunks, dim=f"y_{channel}", join="exact", combine_attrs="drop_conflicts")
+            merged = xr.concat(chunks, dim=f"y_{channel}", join="exact", combine_attrs="drop_conflicts").sortby(f"y_{channel}")
             # Concatenation contains only existing segment rows; gaps are never padded.
             if len(np.unique(merged[f"y_{channel}"])) != merged.sizes[f"y_{channel}"]:
                 raise ValueError("overlapping FCI body rows in repeat cycle")