Read footprint polygons from remote ZIP archives

docs/source/user_guide/index.md

@@ -9,6 +9,7 @@ installation
 cli_usage
 modules_and_bundles
 hooks_and_presets
+remote_archive_footprint
 streams
 recipes
 plugins_and_extensions

docs/source/user_guide/remote_archive_footprint.md

@@ -0,0 +1,52 @@
+# Remote archive footprints
+
+`remote_archive_footprint` reads polygon boundaries from a shapefile inside a
+remote ZIP and attaches them to each Fetchez entry before the archive is
+downloaded. It uses HTTP range reads to retrieve the ZIP directory and the
+selected shapefile components, leaving unrelated files in the archive unread.
+The server must support byte ranges. No additional dependencies are required.
+
+The hook accepts direct HTTP ZIP URLs from any module. It does not inspect
+raster values, interpret project names, select newer surveys, or apply TNM
+resolution rules. This first version supports ZIP shapefiles; standalone
+GeoPackages such as WESM are a separate reader/integration task.
+
+Use the hook on archives known to contain footprint polygons. For example:
+
+```bash
+fetchez run --global-hook list-only -R -124.5/-124.0/41.5/44.0 \
+  tnm --datasets 1_9as --hook remote_archive_footprint
+```
+
+If an archive contains multiple shapefiles, select one by its exact member
+path using `layer`, for example `layer=metadata/footprint.shp`. Paths are
+case-sensitive. Without a selection, exactly one shapefile must be present.
+Use `fields=project/year` to retain specific attribute columns, or omit
+`fields` to keep all columns. Python callers can also pass a list of field
+names. Field names are preserved without TNM-specific renaming.
+
+The hook sets:
+
+- `geometry`: the union of the selected layer's polygons, as a Shapely geometry
+  in WGS84 longitude/latitude. It replaces any previous entry geometry and can
+  be consumed by `spatial_cull`.
+- `footprint_features`: individual records containing `fid`, WGS84 `geometry`
+  as WKT, and `properties` with the selected attributes. These records preserve
+  feature boundaries even though the entry-level geometry is combined.
+- `footprint_layer`: the selected ZIP member path. Together with the unchanged
+  entry `url` and each `fid`, this identifies the source of the extracted data.
+
+Polygon holes are preserved as supplied. Unlike a raster-extent hook, this
+reader does not invent a perimeter or fill holes in published vector geometry.
+It transforms existing vertices; it does not densify sparse curved edges.
+It neither clips to an ROI nor turns individual features into duplicate file
+downloads. Downstream code must interpret coverage and provenance semantics.
+
+Missing layers, required shapefile components, selected fields, CRS, empty
+layers, and invalid/non-polygon geometry raise errors. Extents crossing WGS84
+longitude limits are rejected rather than returned as misleading polygons.
+The hook does not change the pipeline's normal error-handling policy.
+
+Only selected shapefile components are written to a temporary directory, using
+fixed local filenames. ZIP member paths are never used as extraction paths.
+Temporary files and HTTP responses are closed after reading, including errors.

src/fetchez/hooks/remote_archive_footprint.py

@@ -0,0 +1,129 @@
+"""Read polygon footprints from remote ZIP shapefiles before downloading."""
+
+import math
+import shutil
+import zipfile
+from pathlib import Path, PurePosixPath
+from tempfile import TemporaryDirectory
+
+import requests
+import shapely
+from pyogrio.raw import read
+from pyproj import Transformer
+from shapely.ops import transform
+
+from fetchez.core import HttpFile
+from fetchez.hooks import FetchHook
+
+
+class RemoteArchiveFootprintHook(FetchHook):
+    """Attach WGS84 polygons and selected source attributes to each entry.
+
+    ``layer`` is the exact path of a .shp member inside the ZIP. It may be
+    omitted when there is only one shapefile. ``fields`` is a slash-separated
+    list of attribute names; omit it to retain all attributes.
+    """
+
+    name = "remote_archive_footprint"
+    meta_stage = "manifest"
+    meta_desc = "Read footprint polygons from a remote ZIP shapefile."
+
+    def __init__(self, layer=None, fields=None, **kwargs):
+        super().__init__(**kwargs)
+        self.layer = layer
+        self.fields = fields.split("/") if isinstance(fields, str) else fields
+
+    def _read_features(self, url, session):
+        with TemporaryDirectory(prefix="fetchez-footprint-") as directory:
+            with HttpFile(url, session=session) as remote:
+                with zipfile.ZipFile(remote) as archive:
+                    names = archive.namelist()
+                    layers = [name for name in names if name.lower().endswith(".shp")]
+                    layer = self.layer
+                    if layer is None:
+                        if len(layers) != 1:
+                            raise ValueError(
+                                f"Select a footprint layer in {url}: {layers}"
+                            )
+                        layer = layers[0]
+                    if layer not in layers:
+                        raise ValueError(f"Footprint layer not found: {layer} in {url}")
+                    stem = str(PurePosixPath(layer).with_suffix(""))
+                    for suffix in (".shp", ".shx", ".dbf", ".prj", ".cpg"):
+                        matches = [
+                            name
+                            for name in names
+                            if name.lower() == (stem + suffix).lower()
+                        ]
+                        if not matches and suffix == ".cpg":
+                            continue
+                        if len(matches) != 1:
+                            raise ValueError(
+                                f"Missing or ambiguous {suffix}: {layer} in {url}"
+                            )
+                        # Never use archive member paths as local extraction paths.
+                        with archive.open(matches[0]) as source:
+                            with (Path(directory) / ("footprint" + suffix)).open(
+                                "wb"
+                            ) as target:
+                                shutil.copyfileobj(source, target)
+            meta, fids, geometries, columns = read(
+                str(Path(directory) / "footprint.shp"),
+                columns=self.fields,
+                return_fids=True,
+            )
+        if not meta.get("crs"):
+            raise ValueError(f"Footprint layer has no CRS: {layer} in {url}")
+        if self.fields and set(self.fields) - set(meta["fields"]):
+            raise ValueError(f"Footprint fields not found: {self.fields} in {url}")
+        if geometries is None or len(geometries) == 0:
+            raise ValueError(f"Footprint layer is empty: {layer} in {url}")
+        project = Transformer.from_crs(meta["crs"], 4326, always_xy=True)
+        features = []
+        footprints = []
+        for offset, wkb in enumerate(geometries):
+            geometry = shapely.from_wkb(wkb)
+            if (
+                geometry is None
+                or geometry.geom_type not in ("Polygon", "MultiPolygon")
+                or geometry.is_empty
+                or not geometry.is_valid
+            ):
+                raise ValueError(f"Invalid footprint polygon: {layer} in {url}")
+            geometry = transform(project.transform, geometry)
+            if not geometry.is_valid or geometry.is_empty or geometry.area <= 0:
+                raise ValueError(f"Invalid projected footprint: {layer} in {url}")
+            coords = shapely.get_coordinates(geometry)
+            if (
+                (abs(coords[:, 0]) > 180).any()
+                or (abs(coords[:, 1]) > 90).any()
+                or geometry.bounds[2] - geometry.bounds[0] > 180
+            ):
+                raise ValueError(f"Footprint crosses WGS84 limits: {layer} in {url}")
+            properties = {}
+            for name, column in zip(meta["fields"], columns, strict=True):
+                value = column[offset]
+                value = value.item() if hasattr(value, "item") else value
+                if hasattr(value, "isoformat"):
+                    value = value.isoformat()
+                elif isinstance(value, float) and not math.isfinite(value):
+                    value = None
+                properties[str(name)] = value
+            footprints.append(geometry)
+            features.append(
+                {
+                    "fid": int(fids[offset]),
+                    "geometry": shapely.to_wkt(geometry, rounding_precision=-1),
+                    "properties": properties,
+                }
+            )
+        return layer, shapely.union_all(footprints), features
+
+    def run(self, entries):
+        with requests.Session() as session:
+            for _, entry in entries:
+                layer, geometry, features = self._read_features(entry["url"], session)
+                entry["geometry"] = geometry
+                entry["footprint_features"] = features
+                entry["footprint_layer"] = layer
+        return entries

tests/test_remote_archive_footprint.py

@@ -0,0 +1,177 @@
+"""Use real archived shapefiles and local HTTP range responses."""
+
+import io
+import json
+import threading
+import zipfile
+from http.server import BaseHTTPRequestHandler, ThreadingHTTPServer
+
+import numpy as np
+import pytest
+import shapely
+from pyogrio.raw import write
+from pyproj import CRS
+
+from fetchez.hooks.remote_archive_footprint import RemoteArchiveFootprintHook
+from fetchez.hooks.spatial_cull import SpatialCullHook
+from fetchez.hooks.audit import Audit
+
+
+@pytest.fixture
+def remote_archive(tmp_path):
+    state = {"payload": b"", "bytes": 0}
+
+    class Handler(BaseHTTPRequestHandler):
+        def do_HEAD(self):
+            self.send_response(200)
+            self.send_header("Content-Length", str(len(state["payload"])))
+            self.end_headers()
+
+        def do_GET(self):
+            start, end = map(
+                int, self.headers["Range"].removeprefix("bytes=").split("-")
+            )
+            self.send_response(206)
+            self.send_header(
+                "Content-Range", f"bytes {start}-{end}/{len(state['payload'])}"
+            )
+            self.send_header("Content-Length", str(end - start + 1))
+            self.end_headers()
+            state["bytes"] += end - start + 1
+            self.wfile.write(state["payload"][start : end + 1])
+
+        def log_message(self, *args):
+            pass
+
+    server = ThreadingHTTPServer(("127.0.0.1", 0), Handler)
+    thread = threading.Thread(target=server.serve_forever, daemon=True)
+    thread.start()
+
+    def make(
+        polygons=None, crs="EPSG:4326", missing=None, second=False, prefix="nested"
+    ):
+        if polygons is None:
+            polygons = [
+                shapely.box(-124, 43, -123, 44),
+                shapely.box(-123, 43, -122, 44),
+            ]
+        write(
+            tmp_path / "source.shp",
+            geometry=np.array([shapely.to_wkb(p) for p in polygons], dtype=object),
+            field_data=[
+                np.array([f"survey-{i}" for i in range(len(polygons))]),
+                np.full(len(polygons), 2020, dtype="int32"),
+                np.full(len(polygons), np.nan),
+            ],
+            fields=["project", "year", "score"],
+            driver="ESRI Shapefile",
+            geometry_type=polygons[0].geom_type,
+            crs=CRS(crs).to_wkt(),
+        )
+        buffer = io.BytesIO()
+        with zipfile.ZipFile(buffer, "w", compression=zipfile.ZIP_DEFLATED) as archive:
+            archive.writestr(
+                "unrelated.bin", b"x" * 2_000_000, compress_type=zipfile.ZIP_STORED
+            )
+            for path in tmp_path.glob("source.*"):
+                if path.suffix == missing:
+                    continue
+                archive.write(path, f"{prefix}/source{path.suffix.upper()}")
+                if second:
+                    archive.write(path, f"other/source{path.suffix.upper()}")
+        state["payload"] = buffer.getvalue()
+        state["bytes"] = 0
+        return f"http://127.0.0.1:{server.server_port}/data.zip", state
+
+    yield make
+    server.shutdown()
+    server.server_close()
+    thread.join()
+
+
+def test_reads_only_footprints_and_preserves_features(remote_archive):
+    url, state = remote_archive()
+    module = object()
+    entry = {"url": url, "year": 2020, "title": "source"}
+    entries = [(module, entry)]
+    hook = RemoteArchiveFootprintHook(fields="project")
+    assert hook.stage == "manifest"
+    assert hook.run(entries) is entries
+    assert entry["url"] == url
+    assert entry["geometry"].bounds == (-124, 43, -122, 44)
+    assert len(entry["footprint_features"]) == 2
+    assert entry["footprint_features"][1]["properties"] == {"project": "survey-1"}
+    assert entry["footprint_features"][0]["fid"] == 0
+    assert entry["footprint_layer"] == "nested/source.SHP"
+    json.dumps(Audit()._sanitize(entry), allow_nan=False)
+    assert 0 < state["bytes"] < len(state["payload"]) / 4
+    newer = (module, {**entry, "year": 2021})
+    assert SpatialCullHook().run([*entries, newer]) == [newer]
+
+
+def test_requires_selection_for_multiple_layers(remote_archive):
+    url, _ = remote_archive(second=True)
+    with pytest.raises(ValueError, match="Select a footprint layer"):
+        RemoteArchiveFootprintHook().run([(None, {"url": url})])
+    entry = {"url": url}
+    RemoteArchiveFootprintHook(layer="other/source.SHP").run([(None, entry)])
+    assert entry["footprint_layer"] == "other/source.SHP"
+
+
+@pytest.mark.parametrize("suffix", [".shx", ".dbf", ".prj"])
+def test_missing_sidecar_does_not_replace_geometry(remote_archive, suffix):
+    url, _ = remote_archive(missing=suffix)
+    entry = {"url": url, "geometry": "original"}
+    with pytest.raises(ValueError, match="Missing or ambiguous"):
+        RemoteArchiveFootprintHook().run([(None, entry)])
+    assert entry["geometry"] == "original"
+
+
+def test_missing_field_is_rejected(remote_archive):
+    url, _ = remote_archive()
+    with pytest.raises(ValueError, match="fields not found"):
+        RemoteArchiveFootprintHook(fields="absent").run([(None, {"url": url})])
+
+
+def test_projected_footprint(remote_archive):
+    url, _ = remote_archive([shapely.box(0, 0, 1000, 1000)], crs="EPSG:3857")
+    entry = {"url": url}
+    RemoteArchiveFootprintHook().run([(None, entry)])
+    assert entry["geometry"].bounds == pytest.approx((0, 0, 0.00898315, 0.00898315))
+
+
+def test_preserves_polygon_holes(remote_archive):
+    polygon = shapely.box(-124, 43, -122, 45).difference(
+        shapely.box(-123.5, 43.5, -123, 44)
+    )
+    url, _ = remote_archive([polygon])
+    entry = {"url": url}
+    RemoteArchiveFootprintHook().run([(None, entry)])
+    assert entry["geometry"].equals(polygon)
+    properties = entry["footprint_features"][0]["properties"]
+    assert properties == {"project": "survey-0", "year": 2020, "score": None}
+    assert shapely.from_wkt(entry["footprint_features"][0]["geometry"]).equals(polygon)
+    json.dumps(Audit()._sanitize(entry), allow_nan=False)
+
+
+def test_archive_paths_are_not_extracted_verbatim(remote_archive):
+    url, _ = remote_archive(prefix="../../outside")
+    entry = {"url": url}
+    RemoteArchiveFootprintHook().run([(None, entry)])
+    assert entry["geometry"].area == 2
+
+
+def test_non_polygon_is_rejected(remote_archive):
+    url, _ = remote_archive([shapely.Point(0, 0)])
+    with pytest.raises(ValueError, match="Invalid footprint polygon"):
+        RemoteArchiveFootprintHook().run([(None, {"url": url})])
+
+
+def test_empty_manifest():
+    assert RemoteArchiveFootprintHook().run([]) == []
+
+
+def test_antimeridian_footprint_is_rejected(remote_archive):
+    url, _ = remote_archive([shapely.box(179, 0, 181, 1)])
+    with pytest.raises(ValueError, match="WGS84 limits"):
+        RemoteArchiveFootprintHook().run([(None, {"url": url})])