diff --git a/config/preview-folders.json b/config/preview-folders.json new file mode 100644 index 00000000..e40d18cf --- /dev/null +++ b/config/preview-folders.json @@ -0,0 +1,37 @@ +{ + "version": 0, + "status": "DRAFT - awaiting the definitive folder list from the IPMA colleague; the entries below are examples taken from survey-extensions-compact_VV01.xlsx and WILL be replaced by that list", + "description": "Survey-template folders where asset previews (issues #145/#154) are to be generated. Folders NOT listed here get no preview. This file does not affect discovery or metadata extraction (#116). Paths are relative to the survey root. 'extensions' optionally narrows which files in the folder qualify. 'preview_kind' hints the derivation type: raster | vector | auto.", + "folders": [ + { + "folder": "s06-mbes/s05-processed-data", + "extensions": [".xyz", ".tif", ".tiff", ".flt", ".asc"], + "preview_kind": "raster" + }, + { + "folder": "s06-mbes/s09-backscatter/s04-processed-data", + "extensions": [".tif", ".xyz"], + "preview_kind": "raster" + }, + { + "folder": "s07-sss/s05-processed-data", + "extensions": [".tif"], + "preview_kind": "raster" + }, + { + "folder": "s08-chirp/s05-processed-data", + "extensions": [".sgy"], + "preview_kind": "auto" + }, + { + "folder": "s13-uhrs/s05-processed-data", + "extensions": [".sgy"], + "preview_kind": "auto" + }, + { + "folder": "s28-tracks/s02-final", + "extensions": [".shp"], + "preview_kind": "vector" + } + ] +} diff --git a/docs/admin-guide.md b/docs/admin-guide.md index 19218c74..3a094f49 100644 --- a/docs/admin-guide.md +++ b/docs/admin-guide.md @@ -98,10 +98,19 @@ whitelisted machine's IP address. ssh seis-lab-data-production -N -D 1080 ``` - And then use a web browser with proxy support. For example, google chrome: + And then use a web browser with proxy support. + + On Linux, for example with Google Chrome: + + ```shell + google-chrome --proxy-server="socks5://localhost:1080" + ``` + + On macOS, invoke the browser binary through `open`, using a separate profile so the proxy is applied even + if the browser is already running. For example with Brave: ```shell - google-chrome --proxy-server="socks5://localhost:1080 + open -na "Brave Browser" --args --proxy-server="socks5://localhost:1080" --user-data-dir="/tmp/brave-proxy" ``` With this connection in place, the production environment can be found by browsing to diff --git a/docs/admin-guide.pt.md b/docs/admin-guide.pt.md index cc7453c4..282705e8 100644 --- a/docs/admin-guide.pt.md +++ b/docs/admin-guide.pt.md @@ -104,10 +104,19 @@ a partir do endereço IP de uma máquina previamente autorizada. ssh seis-lab-data-production -N -D 1080 ``` - E depois usar um navegador web com suporte de proxy. Por exemplo, o Google Chrome: + E depois usar um navegador web com suporte de proxy. + + Em Linux, por exemplo com o Google Chrome: + + ```shell + google-chrome --proxy-server="socks5://localhost:1080" + ``` + + Em macOS, invocar o binário do navegador através do `open`, usando um perfil separado para que o proxy + seja aplicado mesmo que o navegador já esteja aberto. Por exemplo, com o Brave: ```shell - google-chrome --proxy-server="socks5://localhost:1080 + open -na "Brave Browser" --args --proxy-server="socks5://localhost:1080" --user-data-dir="/tmp/brave-proxy" ``` Com esta ligação ativa, o ambiente de produção pode ser acedido navegando para diff --git a/scripts/README.md b/scripts/README.md new file mode 100644 index 00000000..ab2522e0 --- /dev/null +++ b/scripts/README.md @@ -0,0 +1,16 @@ +# scripts + +One-off utility scripts. Run Python scripts with `uv run python scripts/` (or `python3 scripts/`), shell scripts with `bash scripts/.sh`. Python scripts accept `--help`. + +- **`quick-get-sample-data.sh`** — download a sample of the survey archive from the production server (rsync over ssh). Usage: `bash scripts/quick-get-sample-data.sh `. +- **`validate_extractors.py`** — validate the metadata extractors against a whole archive, header-only and read-only, writing one CSV row per file. Pure standard library so it runs on any `python3` without installing anything; the KMALL and SEG-Y parsing mirrors `tasks/extractors/kmall.py` and `segy.py` and must be kept in sync with them. Usage: `python3 scripts/validate_extractors.py --root --extension .kmall --output report.csv`. +- **`run-archive-validation.sh`** — run the validator against the production archive in chunks, over ssh, leaving nothing on the server (the validator is streamed in, the CSV streamed back). One report per chunk, so a dropped connection costs one chunk rather than a whole multi-hour run. Usage: `bash scripts/run-archive-validation.sh `; see the script for the chunk tables and the `SLD_VALIDATION_*` overrides. +- **`analyse_segy_reports.py`** — analyse the `validate_extractors.py` CSV reports for SEG-Y, printing the agreed checks (geographic files by units code, garbage/partial coverage, small-support and mis-scaled boxes, partial-metadata signatures, errors by size, suspect dates). Pure standard library; encodes the analysis rules so a fresh archive scan is re-audited the same way rather than ad-hoc. Usage: `python3 scripts/analyse_segy_reports.py --reports `. +- **`probe_coords.py`** — read-only forensic probe that dumps the sampled src/cdp coordinate distribution of one or more SEG-Y files, to tell apart a real line stretched by a burst of bad-nav traces from a coordinate field that is wholly noise. Pure standard library; stream it to a host that holds the files (`ssh "python3 - /path/a.sgy" < scripts/probe_coords.py`). Usage: `python3 scripts/probe_coords.py [...]`. +- **`recreate_dirs.py`** — recreate an empty directory tree from a `tree -d` listing. Usage: `python3 scripts/recreate_dirs.py ` (e.g. with `sample-file-tree.txt`). +- **`survey_extensions.py`** — aggregate file extensions per `.example-survey` template folder, from a list of file paths. Usage: `python3 scripts/survey_extensions.py --file-list --template-dir sample-data/.example-survey --output survey-extensions.txt`. +- **`emodnet.py`** — download EMODNet bathymetry tiles. +- **`mbtiles.py`** — build an MBTiles file from a `z/x/y` directory of PNG tiles. +- **`screenshotter.py`** — render a MapLibre map to a PNG (Playwright), for use as a project / survey mission / record image. Uses `map-template.j2.html`. + +Data files: `sample-file-tree.txt` (input for `recreate_dirs.py`), `survey-extensions.txt` (output of `survey_extensions.py`). diff --git a/scripts/analyse_segy_reports.py b/scripts/analyse_segy_reports.py new file mode 100644 index 00000000..0033f68b --- /dev/null +++ b/scripts/analyse_segy_reports.py @@ -0,0 +1,261 @@ +#!/usr/bin/env python3 +"""Analyse the CSV reports produced by validate_extractors.py for SEG-Y. + +Reads every CSV in a reports directory (one per chunk, same header) and prints +the checks agreed for the archive scan, so the analysis is reproducible instead +of ad-hoc. Pure standard library, like the validator itself. The rules encoded +here are the ones the extractor's coverage/bbox decisions rest on, so a fresh +archive scan can be re-audited the same way: + + - geographic files (the only ones that draw a map rectangle) by units code, + never by a non-empty min_lon; + - every row with garbage_count > 0, and every bbox built from a sliver of + accepted points; + - metres-box magnitude per family, so a mis-scaled folder stands out; + - partial-metadata rows by signature, errors split by size, suspect dates. + +Usage: python3 scripts/analyse_segy_reports.py --reports [--top N] +""" + +import argparse +import collections +import csv +import datetime as dt +import glob +import os + +# Portugal geographic window: a loose sanity box, deliberately wide enough for +# the whole EEZ (real deep-water lines reach lon -15.7), so "outside" is a +# prompt to look, not a verdict. +GEO_LON = (-16.0, -5.0) +GEO_LAT = (34.0, 44.0) +# Loose projected-metres envelope for an offshore-Portugal survey (TM06-centric; +# UTM rows legitimately fall outside, which the per-family view separates out). +MET_X = (-400_000.0, 400_000.0) +MET_Y = (-500_000.0, 800_000.0) +SMALL_ACCEPTED = 5 +GEO_UNIT_LABELS = {"arc-seconds", "degrees", "dms"} + + +def fnum(s): + if s is None or s == "": + return None + try: + return float(s) + except ValueError: + return None + + +def units_has_geographic(units_counts): + # units_counts like "1:98;2:2" - a code 2 (arc-seconds), 3 (degrees) or + # 4 (dms) token means at least one geographic trace was sampled. + for token in units_counts.split(";"): + if token and token.split(":", 1)[0] in ("2", "3", "4"): + return True + return False + + +def folder_key(path, depth=3): + return "/".join(path.split("/")[:depth]) + + +def load_reports(reports_dir): + rows = [] + for path in sorted(glob.glob(os.path.join(reports_dir, "*.csv"))): + chunk = os.path.basename(path)[:-4] + with open(path) as fh: + for row in csv.DictReader(fh): + row["_chunk"] = chunk + rows.append(row) + return rows + + +def section(title): + print("\n" + "=" * 78 + f"\n{title}\n" + "=" * 78) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--reports", required=True, help="directory of validate_extractors CSVs" + ) + parser.add_argument( + "--top", type=int, default=25, help="max rows to list per finding" + ) + args = parser.parse_args() + + rows = load_reports(args.reports) + if not rows: + parser.error(f"no CSV reports found in {args.reports}") + ok = [r for r in rows if r["status"] == "ok"] + err = [r for r in rows if r["status"] == "error"] + + section("PER-CHUNK COUNTS") + per_chunk = collections.Counter(r["_chunk"] for r in rows) + for chunk in sorted(per_chunk): + chunk_rows = [r for r in rows if r["_chunk"] == chunk] + n_ok = sum(1 for r in chunk_rows if r["status"] == "ok") + n_er = sum(1 for r in chunk_rows if r["status"] == "error") + print(f" {chunk:24s} {len(chunk_rows):6d} rows {n_ok:6d} ok {n_er:3d} err") + print(f" {'TOTAL':24s} {len(rows):6d} rows {len(ok):6d} ok {len(err):3d} err") + + # ---- geographic files (the only ones that draw a map rectangle) ------ + section("GEOGRAPHIC FILES (units code 2/3/4 present, or geographic label)") + geo = [ + r + for r in ok + if units_has_geographic(r.get("units_counts", "")) + or r.get("coordinate_units", "") in GEO_UNIT_LABELS + ] + print(f" {len(geo)} of {len(ok)} ok rows carry any geographic trace") + for r in geo[: args.top]: + lo, la, xo, xa = ( + fnum(r["min_lon"]), + fnum(r["min_lat"]), + fnum(r["max_lon"]), + fnum(r["max_lat"]), + ) + window = "" + if None not in (lo, la, xo, xa): + inside = ( + GEO_LON[0] <= lo <= GEO_LON[1] + and GEO_LON[0] <= xo <= GEO_LON[1] + and GEO_LAT[0] <= la <= GEO_LAT[1] + and GEO_LAT[0] <= xa <= GEO_LAT[1] + ) + window = " window:OK" if inside else " window:OUTSIDE(suspect)" + print( + f" [{r['_chunk']}] units={r['units_counts']!r} cov={r['coverage']}" + f" lonlat=({r['min_lon']},{r['min_lat']})..({r['max_lon']},{r['max_lat']})" + f"{window} {r['path']}" + ) + + # ---- units distribution --------------------------------------------- + section("UNITS: dominant coordinate_units across ok rows") + for label, n in collections.Counter( + r.get("coordinate_units", "") or "(none)" for r in ok + ).most_common(): + print(f" {label:18s} {n:6d}") + + # ---- garbage / partial coverage (junk-navigation net) --------------- + section("GARBAGE / PARTIAL COVERAGE") + garb = [r for r in ok if (fnum(r.get("garbage_count")) or 0) > 0] + partial = [r for r in ok if r.get("coverage") == "partial"] + print(f" garbage_count > 0 : {len(garb)} coverage=partial : {len(partial)}") + hist = collections.Counter(int(fnum(r.get("garbage_count")) or 0) for r in ok) + print(" garbage_count histogram:", {k: hist[k] for k in sorted(hist)}) + for tag, group in (("garbage>0", garb), ("partial", partial)): + clusters = collections.Counter(folder_key(r["path"]) for r in group) + if clusters: + print(f" {tag} by folder (3-seg):") + for folder, n in clusters.most_common(args.top): + print(f" {n:5d} {folder}") + + # ---- small-support bboxes (accepted below the flag floor) ----------- + section(f"SMALL-SUPPORT BBOXES (0 < accepted_count < {SMALL_ACCEPTED}, box kept)") + small = [ + r + for r in ok + if 0 < (fnum(r.get("accepted_count")) or 0) < SMALL_ACCEPTED + and (r.get("native_min_x") or r.get("min_lon")) + ] + print(f" {len(small)} rows") + for r in small[: args.top]: + print( + f" a={r['accepted_count']} p={r['parsed_count']} " + f"({r['native_min_x']},{r['native_min_y']})" + f"..({r['native_max_x']},{r['native_max_y']}) {r['path']}" + ) + + # ---- metres-box magnitude, per family ------------------------------- + section("METRES BBOX MAGNITUDE (per-family centres; envelope outliers)") + metres = [ + r + for r in ok + if r.get("coordinate_units") in ("metres", "unset") and r.get("native_min_x") + ] + outliers = 0 + families = collections.defaultdict(list) + for r in metres: + x0, y0, x1, y1 = ( + fnum(r["native_min_x"]), + fnum(r["native_min_y"]), + fnum(r["native_max_x"]), + fnum(r["native_max_y"]), + ) + if None in (x0, y0, x1, y1): + continue + if not ( + MET_X[0] <= x0 and x1 <= MET_X[1] and MET_Y[0] <= y0 and y1 <= MET_Y[1] + ): + outliers += 1 + families[folder_key(r["path"], 2)].append(((x0 + x1) / 2, (y0 + y1) / 2)) + print(f" {len(metres)} metres/unset boxes; {outliers} outside the loose envelope") + for family in sorted(families): + centres = families[family] + xs = sorted(c[0] for c in centres) + ys = sorted(c[1] for c in centres) + print( + f" {family:28s} n={len(centres):5d} " + f"x[{xs[0]:>12.0f}..{xs[-1]:>12.0f}] y[{ys[0]:>12.0f}..{ys[-1]:>12.0f}]" + ) + + # ---- partial-metadata rows by signature ----------------------------- + section("PARTIAL-METADATA (ok rows with no bbox) by signature") + nobbox = [r for r in ok if not r.get("native_min_x") and not r.get("min_lon")] + sig = collections.Counter() + example = {} + for r in nobbox: + fmt, spt, tc = ( + r.get("sample_format", ""), + r.get("samples_per_trace", ""), + r.get("trace_count", ""), + ) + if fmt.startswith("format code") or fmt.startswith("format "): + key = "unknown-sample-format" + elif spt == "": + key = "samples_per_trace=0" + elif tc == "": + key = "rev2-bailout-or-no-traces" + elif (fnum(tc) or 0) == 0: + key = "trace_count=0" + else: + key = "no-usable-navigation" + sig[key] += 1 + example.setdefault(key, r) + print(f" {len(nobbox)} ok rows without any bbox") + for key, n in sig.most_common(): + print(f" {n:6d} {key} e.g. {example[key]['path']}") + + # ---- errors, split by size ------------------------------------------ + section("ERRORS split by size_bytes") + too_small = [r for r in err if (fnum(r.get("size_bytes")) or 0) < 3840] + unrecognised = [r for r in err if (fnum(r.get("size_bytes")) or 0) >= 3840] + print(f" < 3840 bytes (too small / empty) : {len(too_small)}") + for r in too_small[: args.top]: + print(f" size={r['size_bytes']} {r['path']}") + print(f" >= 3840 (header unrecognised) : {len(unrecognised)}") + for r in unrecognised[: args.top]: + print(f" size={r['size_bytes']} {r['path']}") + + # ---- suspect dates -------------------------------------------------- + section("DATES: date_begin before 2016 or span > 60 days") + suspect = [] + for r in ok: + begin = r.get("date_begin", "") + if not begin: + continue + try: + d0 = dt.date.fromisoformat(begin) + d1 = dt.date.fromisoformat(r["date_end"]) if r.get("date_end") else d0 + except ValueError: + continue + if d0.year < 2016 or (d1 - d0).days > 60: + suspect.append((r, d0, d1)) + print(f" {len(suspect)} suspect") + for r, d0, d1 in suspect[: args.top]: + print(f" {d0}..{d1} ({(d1 - d0).days}d) {r['path']}") + + +if __name__ == "__main__": + main() diff --git a/scripts/probe_coords.py b/scripts/probe_coords.py new file mode 100644 index 00000000..9eb8b396 --- /dev/null +++ b/scripts/probe_coords.py @@ -0,0 +1,157 @@ +#!/usr/bin/env python3 +"""Dump the sampled src/cdp coordinates of one or more SEG-Y files. + +A read-only, standard-library forensic probe. It reads exactly the ~100 sampled +trace headers the extractor uses (segy.py) and prints, separately for the src +and cdp fields, how the scalco-scaled coordinates are distributed. The point is +to tell apart the two junk-navigation modes the archive holds: + + - a tight cluster with a few wild outliers -> a REAL line stretched by a + burst of bad-nav traces (the box should be discarded, the file is real); + - a broad spread across many values -> the coordinate field itself is + noise (true junk navigation). + +Stream it to a host that holds the files, so nothing is written there: + + ssh "python3 - /path/a.sgy /path/b.sgy" < scripts/probe_coords.py + +Usage: python3 scripts/probe_coords.py [ ...] +""" + +import argparse +import os +import statistics +import struct +import sys + +_TEXTUAL = 3200 +_BINHDR = 400 +_TRHDR = 240 +_DATA_START = _TEXTUAL + _BINHDR +_SAMPLE_TARGET = 100 +_INT32_SENTINELS = (-(2**31), 2**31 - 1) +_SAMPLE_FORMATS = {1: 4, 2: 4, 3: 2, 5: 4, 6: 8, 8: 1, 10: 4, 11: 2, 16: 1} + +# binary: sample_interval(16) samples_per_trace(20) sample_format(24) rev(300) +# ext_header_count(304) extra_trace_headers(306) +_BIN_FMT = "16x H 2x H 2x H 274x H 2x h i" +_BIN = {">": struct.Struct(">" + _BIN_FMT), "<": struct.Struct("<" + _BIN_FMT)} +# trace: scalco(70) src_x(72) src_y(76) units(88) ... cdp_x(180) cdp_y(184) +_TR_FMT = "70x h i i 8x H 66x H H H H H 14x i i" +_TR = {">": struct.Struct(">" + _TR_FMT), "<": struct.Struct("<" + _TR_FMT)} + + +def scale(value, scalco): + if scalco < 0: + return value / -scalco + if scalco > 0: + return float(value * scalco) + return float(value) + + +def sample_indices(n): + if n <= _SAMPLE_TARGET: + return list(range(n)) + last = n - 1 + return sorted({s * last // (_SAMPLE_TARGET - 1) for s in range(_SAMPLE_TARGET)}) + + +def summarise(name, xs, ys): + if not xs: + print(f" {name}: no usable points") + return + for axis, vals in (("x", xs), ("y", ys)): + vals = sorted(vals) + med = statistics.median(vals) + near = sum(1 for v in vals if abs(v - med) <= 5000) + print( + f" {name}.{axis}: n={len(vals):3d} min={vals[0]:12.1f} " + f"med={med:12.1f} max={vals[-1]:12.1f} " + f"within5km_of_med={near} outliers={len(vals) - near}" + ) + + +def probe(path): + size = os.path.getsize(path) + print(f"\n{path} ({size} bytes)") + with open(path, "rb") as fh: + fh.seek(_TEXTUAL) + buf = fh.read(_BINHDR) + order = None + for candidate in (">", "<"): + si, spt, fmt, rev, exth, extra = _BIN[candidate].unpack_from(buf) + if 1 <= fmt <= 16: + order = candidate + break + if order is None: + print(" not a recognisable SEG-Y binary header") + return + bps = _SAMPLE_FORMATS.get(fmt) + rev_major = rev >> 8 if rev > 0xFF else rev + print( + f" order={order!r} fmt={fmt} bytes/sample={bps} " + f"samples/trace={spt} rev={rev_major}" + ) + if not bps or spt == 0: + print(" cannot locate traces (unknown format or samples/trace 0)") + return + data_start = _DATA_START + if ( + rev_major >= 1 + and 0 < exth <= 100 + and _DATA_START + exth * _TEXTUAL + _TRHDR <= size + ): + data_start = _DATA_START + exth * _TEXTUAL + tlen = _TRHDR + spt * bps + ntr = (size - data_start) // tlen + print(f" trace_count={ntr}") + srcx, srcy, cdpx, cdpy = [], [], [], [] + scalcos, unit_codes = set(), set() + for index in sample_indices(ntr): + fh.seek(data_start + index * tlen) + b = fh.read(_TRHDR) + if len(b) < _TR[order].size: + break + scalco, sx, sy, units, *_rest, cx, cy = _TR[order].unpack_from(b) + scalcos.add(scalco) + unit_codes.add(units) + if ( + sx not in _INT32_SENTINELS + and sy not in _INT32_SENTINELS + and not (sx == 0 and sy == 0) + ): + srcx.append(scale(sx, scalco)) + srcy.append(scale(sy, scalco)) + if ( + cx not in _INT32_SENTINELS + and cy not in _INT32_SENTINELS + and not (cx == 0 and cy == 0) + ): + cdpx.append(scale(cx, scalco)) + cdpy.append(scale(cy, scalco)) + print( + f" scalco seen={sorted(scalcos)} units codes seen={sorted(unit_codes)}" + ) + summarise("src", srcx, srcy) + summarise("cdp", cdpx, cdpy) + print( + " READING: within5km_of_med high + few outliers => real line + bad traces;" + ) + print( + " values spread across many => the field itself is noise (junk)." + ) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("files", nargs="+", help="SEG-Y files to probe") + args = parser.parse_args() + for path in args.files: + try: + probe(path) + except Exception as err: # a survey - report and carry on + print(f"\n{path}\n ERROR {type(err).__name__}: {err}", file=sys.stderr) + + +if __name__ == "__main__": + main() diff --git a/scripts/quick-get-sample-data.sh b/scripts/quick-get-sample-data.sh old mode 100644 new mode 100755 index 8b73a93e..d4cdb60f --- a/scripts/quick-get-sample-data.sh +++ b/scripts/quick-get-sample-data.sh @@ -8,17 +8,51 @@ fi SIMULATED_ARCHIVE_ROOT=$1 +# This script needs GPL rsync (rsync.samba.org) for resume (--partial) and +# include/exclude filters. +RSYNC=rsync +if ! command -v "${RSYNC}" >/dev/null 2>&1 || "${RSYNC}" --version 2>&1 | grep -qi openrsync; then + RSYNC="" + for candidate in /opt/homebrew/bin/rsync /usr/local/bin/rsync; do + if [[ -x "${candidate}" ]] && ! "${candidate}" --version 2>&1 | grep -qi openrsync; then + RSYNC=${candidate} + break + fi + done +fi +if [[ -z "${RSYNC}" ]]; then + echo "GPL rsync not found. Install it." >&2 + exit 1 +fi + +# Fetch a single file. rsync resumes interrupted transfers (--partial) fetch() { local remote_path=$1 local local_dir=$2 - local filename - filename=$(basename "${remote_path}") - if [[ -f "${local_dir}/${filename}" ]]; then - echo "Skipping ${filename} (already present)" - return - fi mkdir -p "${local_dir}" - scp "seis-lab-data-production:${remote_path}" "${local_dir}/" + # -s (--protect-args): paths sent literally, so spaces in remote paths work + # across rsync versions. --append-verify: resume from the partial offset at + # full speed (no slow delta scan), then verify the whole-file checksum. + "${RSYNC}" -aP -s --append-verify "seis-lab-data-production:${remote_path}" "${local_dir}/" +} + +# Fetch a primary file together with all same-stem siblings in its directory +# (shapefile sidecars, FLT+.hdr, ECW+.eww/.prj, S-57 base .000 + updates, a +# GeoTIFF + its .aux.xml, ...). Pass the path to the primary file; the {stem}.* +# siblings are selected with rsync include/exclude filters (no remote glob), so +# it also works when the remote directory path contains spaces. The stem strips +# only the last extension, so it matches sibling files but not unrelated ones. +fetch_with_sidecars() { + local remote_primary=$1 + local local_dir=$2 + local base + base=$(basename "${remote_primary}") + local stem=${base%.*} + local remote_dir + remote_dir=$(dirname "${remote_primary}") + mkdir -p "${local_dir}" + "${RSYNC}" -aP -s --append-verify --include="${stem}.*" --exclude="*" \ + "seis-lab-data-production:${remote_dir}/" "${local_dir}/" } # raw bathymetry (82MB) @@ -36,11 +70,11 @@ fetch /mnt/seislab_data/surveys/owf-seism-2024/s06-mbes/s05-processed-data/LEI_A "${SIMULATED_ARCHIVE_ROOT}"/surveys/owf-seism-2024/s06-mbes/s05-processed-data/ # HUGE FILE!: seismic uhrs raw data (139GB) -fetch '/mnt/seislab_data/surveys/owf-seism-2024/s13-uhrs/s02-raw-data/BACKED\ UP\ -\ FIGUEIRA/FIG-24-M001/FIG-24-M001_rev1.segy' \ +fetch '/mnt/seislab_data/surveys/owf-seism-2024/s13-uhrs/s02-raw-data/BACKED UP - FIGUEIRA/FIG-24-M001/FIG-24-M001_rev1.segy' \ "${SIMULATED_ARCHIVE_ROOT}/surveys/owf-seism-2024/s13-uhrs/s02-raw-data/BACKED UP - FIGUEIRA/FIG-24-M001" # HUGE FILE!: seismic uhrs raw data (84GB) -fetch '/mnt/seislab_data/surveys/owf-seism-2024/s13-uhrs/s02-raw-data/BACKED\ UP\ -\ FIGUEIRA/FIG-24-M010/FIG-24-M010_rev1.segy' \ +fetch '/mnt/seislab_data/surveys/owf-seism-2024/s13-uhrs/s02-raw-data/BACKED UP - FIGUEIRA/FIG-24-M010/FIG-24-M010_rev1.segy' \ "${SIMULATED_ARCHIVE_ROOT}/surveys/owf-seism-2024/s13-uhrs/s02-raw-data/BACKED UP - FIGUEIRA/FIG-24-M010" # seismic uhrs qc data (~2GB) @@ -48,8 +82,8 @@ fetch /mnt/seislab_data/surveys/owf-seism-2024/s13-uhrs/s03-qc/s03-bstk/FIG-24-M "${SIMULATED_ARCHIVE_ROOT}"/surveys/owf-seism-2024/s13-uhrs/s03-qc/s03-bstk/ # HUGE FILE!: seismic uhrs processed prestack data (~138GB) -# fetch /mnt/seislab_data/surveys/owf-seism-2024/s13-uhrs/s05-processed-data/s01-prestk/s01-mul/FIG-24-M001_CMP_MUL.sgy \ -# "${SIMULATED_ARCHIVE_ROOT}"/surveys/owf-seism-2024/s13-uhrs/s05-processed-data/s01-prestk/s01-mul/ +fetch /mnt/seislab_data/surveys/owf-seism-2024/s13-uhrs/s05-processed-data/s01-prestk/s01-mul/FIG-24-M001_CMP_MUL.sgy \ + "${SIMULATED_ARCHIVE_ROOT}"/surveys/owf-seism-2024/s13-uhrs/s05-processed-data/s01-prestk/s01-mul/ # seismic uhrs processed mul data (~2GB) fetch /mnt/seislab_data/surveys/owf-seism-2024/s13-uhrs/s05-processed-data/s02-poststk/s01-twt/s01-mul/FIG-24-M001_MUL.sgy \ @@ -70,3 +104,90 @@ fetch /mnt/seislab_data/surveys/owf-seism-2024/s13-uhrs/s05-processed-data/s03-v # seismic uhrs interval velocities (~2.2GB) fetch /mnt/seislab_data/surveys/owf-seism-2024/s13-uhrs/s05-processed-data/s03-velocities/s03-vel-int/FIG-24-M001_VEL_INT.sgy \ "${SIMULATED_ARCHIVE_ROOT}"/surveys/owf-seism-2024/s13-uhrs/s05-processed-data/s03-velocities/s03-vel-int/ + +# seismic raster files +fetch /mnt/seislab_data/surveys/owf-seism-2024/s02-preparation/s02-lineplan/s01-final/01-Line_Plan/Reference_bathymetry/EMODnet/F3_2022.tif \ + "${SIMULATED_ARCHIVE_ROOT}"/surveys/owf-seism-2024/s02-preparation/s02-lineplan/s01-final/01-Line_Plan/Reference_bathymetry/ + +fetch /mnt/seislab_data/surveys/owf-seism-2024//s04-gis-master-survey/s01-final/negative_depth_seism2024/negative_depth_FIG.tif \ + "${SIMULATED_ARCHIVE_ROOT}"/surveys/owf-seism-2024//s04-gis-master-survey/s01-final/negative_depth_seism2024/ + +fetch /mnt/seislab_data/surveys/owf-seism-2024/s04-gis-master-survey/s01-final/negative_depth_seism2024/FIG_All_Mainlines_and_Xlines_MBES_Grid_4m.tif \ + "${SIMULATED_ARCHIVE_ROOT}"/surveys/owf-seism-2024/s04-gis-master-survey/s01-final/negative_depth_seism2024/ + +# shapefiles + sidecars +fetch_with_sidecars /mnt/seislab_data/surveys/owf-seism-2024/s02-preparation/s02-lineplan/s01-final/01-Line_Plan/Reference_bathymetry/ingmar.shp \ + "${SIMULATED_ARCHIVE_ROOT}"/surveys/owf-seism-2024/s02-preparation/s02-lineplan/s01-final/01-Line_Plan/Reference_bathymetry/ + +fetch_with_sidecars /mnt/seislab_data/surveys/owf-seism-2024/s09-sbp-other/s05-processed-data/s01-twt/oars/2024-09-29/NAV/ShapeFiles/FIG-LEI-24-T002_20240928184838_000.shp \ + "${SIMULATED_ARCHIVE_ROOT}"/surveys/owf-seism-2024/s09-sbp-other/s05-processed-data/s01-twt/oars/2024-09-29/NAV/ShapeFiles/ + +# =========================================================================== +# GDAL-supported samples for the metadata-stripping / quick-preview assessment. +# One sizeable file per format (the largest sensible instance where one exists), +# broad raster + vector coverage. GeoTIFF DTM, XYZ and small shapefiles are +# already covered above, so they are not repeated here. +# =========================================================================== + +# --- GDAL raster samples --- + +# GeoTIFF, scanned/georeferenced IH chart (~1GB). Carries a .tif.aux.xml sidecar +# (GDAL PAM: stats + georef/metadata) plus embedded TIFF tags, so it is the best +# case for the metadata-stripping assessment. Pulled with its sidecar. +fetch_with_sidecars /mnt/seislab_data/surveys/sat-uhrs-2025/s04-gis-master-survey/s01-final/mapas_georef/IH_SED_SUP_folha05_2005_alterado.tif \ + "${SIMULATED_ARCHIVE_ROOT}"/surveys/sat-uhrs-2025/s04-gis-master-survey/s01-final/mapas_georef/ + +# ESRI ASCII Grid / AAIGrid (~6.4GB). Plain ASCII: large and slow to parse, and +# no embedded CRS (would need a .prj, none present here). Single file. +fetch /mnt/seislab_data/surveys/owf-2025/s06-mbes/s05-processed-data/03_DTM_POINTFILE/04_ASC/FIG25_MBES_AVG_ZH_2m_asc_TM06_v0_EXT.asc \ + "${SIMULATED_ARCHIVE_ROOT}"/surveys/owf-2025/s06-mbes/s05-processed-data/03_DTM_POINTFILE/04_ASC/ + +# ESRI float grid / EHdr (~7.5GB). The .hdr is REQUIRED by GDAL to read the .flt; +# .flt.gi / .flt.xml are auxiliary. Pulled with all same-stem siblings. +fetch_with_sidecars /mnt/seislab_data/surveys/owf-2025/s17-magnetometer/s03-qc/s01-24h-deliverables/s02-flt/20250819/OI_634_A7805_TVG_LEI_25_nT_Residual_LMB04_M081.flt \ + "${SIMULATED_ARCHIVE_ROOT}"/surveys/owf-2025/s17-magnetometer/s03-qc/s01-24h-deliverables/s02-flt/20250819/ + +# ECW (~3MB). Needs the (proprietary) ECW driver/plugin in GDAL to read. Ships a +# .eww world file and a .prj; pulled as sidecars. +fetch_with_sidecars /mnt/seislab_data/surveys/mineplat-sampling-2019/s06-mbes/s10-tracklines/2019-04-04/BACK_20190404.ecw \ + "${SIMULATED_ARCHIVE_ROOT}"/surveys/mineplat-sampling-2019/s06-mbes/s10-tracklines/2019-04-04/ + +# --- GDAL vector samples --- + +# Line shapefile, the vector "monster" (~1GB .shp). Needs .shx/.dbf/.prj sidecars. +fetch_with_sidecars /mnt/seislab_data/surveys/owf-2025/s28-tracks/s02-final/FIG25_SBP_CHIRP_Trackplot_ln_TM06_v1_EXT.shp \ + "${SIMULATED_ARCHIVE_ROOT}"/surveys/owf-2025/s28-tracks/s02-final/ + +# DXF / CAD (~25MB). Single file; no embedded CRS in the DXF itself. +fetch /mnt/seislab_data/surveys/owf-2025/s04-gis-master-survey/Permit_areas.dxf \ + "${SIMULATED_ARCHIVE_ROOT}"/surveys/owf-2025/s04-gis-master-survey/ + +# AutoCAD DWG / CAD, georeferenced survey line plan (~0.05MB). Binary CAD read by the +# OGR DWG driver (a different path from the ASCII DXF above). Carries a CRS (TM06), so +# it exercises bbox/CRS extraction. Single file. +fetch /mnt/seislab_data/surveys/owf-2026/s02-preparation/s02-lineplan/s01-final/2026-04-21_Main_LEI_Survey_Lines/DWG/105634-Line_Plan_TM06_20260421.dwg \ + "${SIMULATED_ARCHIVE_ROOT}"/surveys/owf-2026/s02-preparation/s02-lineplan/s01-final/2026-04-21_Main_LEI_Survey_Lines/DWG/ + +# AutoCAD DWG, largest instance (~15MB; spaces in the remote path). Vessel general- +# arrangement drawing: mechanical CAD, likely no usable georeferencing; included as the +# size / metadata-stripping stress case. Single file. +fetch '/mnt/seislab_data/surveys/owf-seism-2024/s02-preparation/s01-documents/1.Vessel_MarioRuivo_Equipment_info/PLANO DE ARRANJO GERAL FINAL _ GRAFICA.dwg' \ + "${SIMULATED_ARCHIVE_ROOT}/surveys/owf-seism-2024/s02-preparation/s01-documents/1.Vessel_MarioRuivo_Equipment_info/" + +# GeoPackage (~0.3MB). Self-contained SQLite container; single file. +fetch /mnt/seislab_data/surveys/owf-2025/s02-preparation/s02-lineplan/s02-old/20250428_LinePlan/IPMA_2025_Prelim_Survey_Line_Plan_20250428.gpkg \ + "${SIMULATED_ARCHIVE_ROOT}"/surveys/owf-2025/s02-preparation/s02-lineplan/s02-old/20250428_LinePlan/ + +# GeoJSON (~0.2MB). Single file; WGS84 by spec unless a crs member says otherwise. +fetch /mnt/seislab_data/surveys/owf-seism-2024/s02-preparation/s02-lineplan/s02-old/LinePlan_Planeamento_preliminar/Planning_Preliminary_lineplan.geojson \ + "${SIMULATED_ARCHIVE_ROOT}"/surveys/owf-seism-2024/s02-preparation/s02-lineplan/s02-old/LinePlan_Planeamento_preliminar/ + +# KML (~0.7MB; spaces in the remote path). Single file; KML coords are WGS84. +fetch '/mnt/seislab_data/surveys/mineplat-geophysical-2019/s06-mbes/s04-processing-flows/PDS Projects/Testes_Pinger/LinhaCosta_UTM29N.kml' \ + "${SIMULATED_ARCHIVE_ROOT}/surveys/mineplat-geophysical-2019/s06-mbes/s04-processing-flows/PDS Projects/Testes_Pinger/" + +# S-57 ENC chart (~0.9MB; spaces in the remote path). Multi-file dataset: base +# cell .000 plus sequential updates .001..008 (+ .TXT/.met). GDAL OGR S57 driver +# applies the updates onto the base. Pulled as same-stem siblings. +fetch_with_sidecars '/mnt/seislab_data/surveys/mineplat-2017/s06-mbes/s02-raw-data/2017_04_01_MBES/Projects Common Files/S-57/PT324205/PT324205.000' \ + "${SIMULATED_ARCHIVE_ROOT}/surveys/mineplat-2017/s06-mbes/s02-raw-data/2017_04_01_MBES/Projects Common Files/S-57/PT324205/" diff --git a/scripts/run-archive-validation.sh b/scripts/run-archive-validation.sh new file mode 100755 index 00000000..fe82c552 --- /dev/null +++ b/scripts/run-archive-validation.sh @@ -0,0 +1,106 @@ +#!/usr/bin/env bash +# Run scripts/validate_extractors.py against the production archive, in chunks. +# +# The validator is streamed in over stdin and its CSV streamed back, so nothing +# is ever written on the production server. Each chunk produces its own report, +# which is what makes a multi-hour run survivable: a dropped connection costs +# one chunk, not the night. +# +# Usage: bash scripts/run-archive-validation.sh +# Override the defaults with the SLD_VALIDATION_HOST / SLD_VALIDATION_ROOT / +# SLD_VALIDATION_OUT environment variables. + +set -uo pipefail +# NOTE: -e is deliberately absent. A chunk that fails (broken pipe, unreadable +# subtree) must not cancel the chunks after it. + +HOST=${SLD_VALIDATION_HOST:-seis-lab-data-production} +ROOT=${SLD_VALIDATION_ROOT:-/mnt/seislab_data/surveys} +OUT=${SLD_VALIDATION_OUT:-${HOME}/archive-validation-reports} +VALIDATOR=scripts/validate_extractors.py +# The two-hop ssh route to production has no keepalives of its own, so a +# silently dead stream would hang until morning instead of failing. +KEEPALIVE=(-o ServerAliveInterval=30 -o ServerAliveCountMax=6) + +# Chunks are "namesubpathextensionexpected", one per line. The +# expected counts come from the 2026-07-20 archive scan and are a sanity check +# only: os.walk skips an unreadable subtree silently and still exits 0, so a +# short report is otherwise indistinguishable from a complete one. +# +# SEG-Y needs one pass per extension, and the chunk boundaries follow the +# archive's own structure. The "previews" set is everything the IPMA +# preview-folder list covers (plus all of owf-seism-2024 and sat-uhrs-2025); +# "rest" is the raw and QC material, which the catalog still ingests. +read -r -d '' CHUNKS_KMALL <<'EOF' +kmall-2025 owf-2025 .kmall 5222 +kmall-2024 owf-seism-2024 .kmall 442 +EOF + +read -r -d '' CHUNKS_SEGY_PREVIEWS <<'EOF' +2024-sgy owf-seism-2024 .sgy 2814 +2024-segy owf-seism-2024/s13-uhrs .segy 337 +satuhrs-segy sat-uhrs-2025 .segy 88 +satuhrs-sgy sat-uhrs-2025 .sgy 5 +2025-innomar-proc owf-2025/s10-innomar/s05-processed-data .sgy 3348 +2025-uhrs-proc owf-2025/s13-uhrs/s05-processed-data .sgy 3732 +2025-chirp-proc owf-2025/s08-chirp/s05-processed-data .sgy 12206 +EOF + +read -r -d '' CHUNKS_SEGY_REST <<'EOF' +2025-uhrs-qc owf-2025/s13-uhrs/s03-qc .sgy 2352 +2025-uhrs-raw-segy owf-2025/s13-uhrs/s02-raw-data .segy 1797 +2025-uhrs-raw-sgy owf-2025/s13-uhrs/s02-raw-data .sgy 594 +2025-innomar-raw owf-2025/s10-innomar/s02-raw-data .sgy 7169 +2025-innomar-qc owf-2025/s10-innomar/s03-qc .sgy 6791 +2025-chirp-qc owf-2025/s08-chirp/s03-qc .sgy 9693 +2025-chirp-raw owf-2025/s08-chirp/s02-raw-data .sgy 15126 +2025-reporting owf-2025/s30-reporting .sgy 9 +EOF + +case "${1:-}" in + kmall) CHUNKS=${CHUNKS_KMALL} ;; + segy-previews) CHUNKS=${CHUNKS_SEGY_PREVIEWS} ;; + segy-rest) CHUNKS=${CHUNKS_SEGY_REST} ;; + *) + cat >&2 < + + kmall every .kmall in the two priority surveys (~3.5 h) + segy-previews the SEG-Y the IPMA preview-folder list covers, plus all of + owf-seism-2024 and sat-uhrs-2025 (~8.7 h) + segy-rest the remaining SEG-Y: raw and QC material (~8.9 h) + +Run from the repository root. On a laptop, wrap it in caffeinate and stay on +mains power: + + caffeinate -i bash $0 segy-previews +USAGE + exit 1 + ;; +esac + +if [[ ! -f "${VALIDATOR}" ]]; then + echo "Run from the repository root: ${VALIDATOR} not found." >&2 + exit 1 +fi +mkdir -p "${OUT}" + +date "+start: %Y-%m-%d %H:%M:%S (set: ${1}, host: ${HOST})" +while IFS=$'\t' read -r name subpath extension expected; do + [[ -z "${name}" ]] && continue + echo + date "+%H:%M:%S >>> ${name} (${expected} files expected)" + ssh "${KEEPALIVE[@]}" "${HOST}" \ + "nice -n 19 ionice -c3 python3 - --root ${ROOT}/${subpath} --extension ${extension} --output /dev/stdout" \ + < "${VALIDATOR}" > "${OUT}/${name}.csv" + got=$(( $(wc -l < "${OUT}/${name}.csv") - 1 )) + if (( got < expected )); then + echo " WARNING: ${got} rows, expected ${expected} - unreadable subtree or dropped connection" + else + echo " ${got} rows" + fi +done <<< "${CHUNKS}" + +echo +date "+end: %Y-%m-%d %H:%M:%S" +echo "reports in ${OUT}" diff --git a/scripts/survey_extensions.py b/scripts/survey_extensions.py new file mode 100644 index 00000000..66c28c04 --- /dev/null +++ b/scripts/survey_extensions.py @@ -0,0 +1,210 @@ +#!/usr/bin/env python3 +"""Aggregate file extensions per template folder across all surveys. + +Reads a list of file paths (one per line, relative to the surveys root, e.g. +``./owf-seism-2024/s06-mbes/s02-raw-data/file.kmall``) produced on the server with: + + cd /mnt/seislab_data/surveys && find . -type f > survey-files.txt + +For every file it finds the deepest folder of the ``.example-survey`` template that +is an ancestor of (or equal to) the file's directory, and records the file's +extension there. Files that live in extra (non-template) subfolders therefore roll +up to the nearest template folder. The result is the union of extensions per +template folder across all surveys, with no per-survey identification. + +Output (default): a tree of the whole template with the extensions (and file +counts) annotated on each folder, plus a global list of all extensions by count. +Aggregation is kept separate from rendering, so emitting CSV/Excel later is trivial. + +Usage: + python3 scripts/survey_extensions.py --file-list survey-files.txt \\ + --template-dir sample-data/.example-survey --output survey-extensions.txt + + # or read the list from stdin (single command): + ssh seis-lab-data-production 'cd /mnt/seislab_data/surveys && find . -type f' \\ + | python3 scripts/survey_extensions.py --template-dir sample-data/.example-survey +""" + +from __future__ import annotations + +import argparse +import collections +import os +import sys + +# OS noise that is never worth indexing. +JUNK_NAMES = {".ds_store", "thumbs.db"} + + +def load_template_folders(template_dir: str) -> set: + """Return the set of folder paths (relative, posix) inside the template.""" + folders = set() + for dirpath, _, _ in os.walk(template_dir): + rel = os.path.relpath(dirpath, template_dir).replace(os.sep, "/") + if rel != ".": + folders.add(rel) + return folders + + +def extension_of(filename: str) -> str: + """Lowercased extension including the dot, or ``(no-ext)``. + + Uses the last dotted component, so ``x.tif.aux.xml`` -> ``.xml``. A leading dot + with no other dot (hidden file) counts as no extension. + """ + dot = filename.rfind(".") + return filename[dot:].lower() if dot > 0 else "(no-ext)" + + +def aggregate(lines, template_folders): + """Map each file to its nearest template folder and count extensions. + + Returns ``(per_folder, outside, n_files, surveys)`` where ``per_folder`` is + ``folder -> Counter(ext)`` and ``outside`` holds files with no template ancestor. + """ + per_folder = collections.defaultdict(collections.Counter) + outside = collections.Counter() + surveys = set() + n_files = 0 + rollup_cache = {} + + def rollup(reldir): + cached = rollup_cache.get(reldir, 0) + if cached != 0: + return cached + cur = reldir + while cur: + if cur in template_folders: + rollup_cache[reldir] = cur + return cur + cur = cur.rsplit("/", 1)[0] if "/" in cur else "" + rollup_cache[reldir] = None + return None + + for line in lines: + path = line.rstrip("\n") + if path.startswith("./"): + path = path[2:] + sep = path.find("/") + if sep < 0: + continue # a file directly under the surveys root, with no survey - skip + rest = path[sep + 1 :] + slash = rest.rfind("/") + reldir, filename = ( + (rest[:slash], rest[slash + 1 :]) if slash >= 0 else ("", rest) + ) + if filename.lower() in JUNK_NAMES: + continue + n_files += 1 + surveys.add(path[:sep]) + ext = extension_of(filename) + folder = rollup(reldir) if reldir else None + if folder is None: + outside[ext] += 1 + else: + per_folder[folder][ext] += 1 + return per_folder, outside, n_files, surveys + + +def render_tree(template_folders, per_folder): + """Render the whole template as an indented tree annotated with extensions.""" + # Build a nested dict from the template folder paths. + root = {} + for folder in template_folders: + node = root + for part in folder.split("/"): + node = node.setdefault(part, {}) + + lines = [] + + def walk(node, parts, prefix): + names = sorted(node) + for i, name in enumerate(names): + last = i == len(names) - 1 + connector = "└── " if last else "├── " + counter = per_folder.get("/".join(parts + [name])) + exts = ( + " " + ", ".join(f"{e}({c})" for e, c in counter.most_common()) + if counter + else "" + ) + lines.append(f"{prefix}{connector}{name}/{exts}") + walk(node[name], parts + [name], prefix + (" " if last else "│ ")) + + walk(root, [], "") + return lines + + +def main() -> int: + parser = argparse.ArgumentParser( + description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter + ) + parser.add_argument( + "--file-list", help="file with one path per line (default: stdin)" + ) + parser.add_argument( + "--template-dir", required=True, help="the .example-survey template directory" + ) + parser.add_argument("--output", help="output file (default: stdout)") + args = parser.parse_args() + + template_folders = load_template_folders(args.template_dir) + + source = ( + open(args.file_list, encoding="utf-8", errors="replace") + if args.file_list + else sys.stdin + ) + try: + per_folder, outside, n_files, surveys = aggregate(source, template_folders) + finally: + if args.file_list: + source.close() + + # Global extension universe across all folders - handy for spotting noise. + universe = collections.Counter() + for counter in per_folder.values(): + universe.update(counter) + universe.update(outside) + + out = [] + out.append("# Survey file extensions per .example-survey folder") + out.append( + f"# files: {n_files:,} | surveys: {len(surveys)} | " + f"template folders: {len(template_folders)} | folders with files: {len(per_folder)} | " + f"distinct extensions: {len(universe)}" + ) + out.append( + "# Each folder lists ext(file_count), aggregated across all surveys (no per-survey detail)." + ) + out.append( + "# Review and delete the extensions that are noise; what remains is what to index." + ) + out.append("") + out.append("## All extensions, by file count") + out.extend(f" {ext}\t{count}" for ext, count in universe.most_common()) + out.append("") + out.append("## Template tree (.example-survey)") + out.append(".example-survey/") + out.extend(render_tree(template_folders, per_folder)) + if outside: + out.append("") + out.append("## Files outside the template (no template ancestor)") + out.extend(f" {ext}\t{count}" for ext, count in outside.most_common()) + + text = "\n".join(out) + "\n" + if args.output: + with open(args.output, "w", encoding="utf-8") as fh: + fh.write(text) + print( + f"wrote {args.output}: {n_files:,} files, {len(per_folder)} folders with extensions, " + f"{len(universe)} distinct extensions", + file=sys.stderr, + ) + else: + sys.stdout.write(text) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/validate_extractors.py b/scripts/validate_extractors.py new file mode 100644 index 00000000..93b75474 --- /dev/null +++ b/scripts/validate_extractors.py @@ -0,0 +1,580 @@ +#!/usr/bin/env python3 +"""Standalone extractor validation harness for production-side runs. + +Walks archive files with header-only reads and writes a CSV report, one row per +file. Pure standard library on purpose: it must run on any python3 (e.g. the +production server) without installing anything, so it does NOT import the +seis_lab_data package. The KMALL and SEG-Y parsing below mirrors +src/seis_lab_data/tasks/extractors/kmall.py and segy.py - keep them in sync. +The CSV also carries report-only diagnostics (parsed/garbage/accepted counts, +the units and coordinate-source histograms) that the extractors deliberately do +not keep: at scale we need the tallies behind a conclusion, not just the +conclusion. Those are not a divergence from the mirrored parsing. + +Typical production run (read-only, low priority, resumable): + + nice -n 19 ionice -c3 python3 validate_extractors.py \ + --root /mnt/seislab_data/surveys \ + --extension .kmall \ + --output ~/kmall-report.csv \ + --resume + +Interrupt at will; --resume skips files already present in the output. + +SEG-Y takes one pass per extension into the same report, the second with +--resume (which is also what keeps the first pass's rows): + + python3 validate_extractors.py ... --extension .sgy --output ~/segy-report.csv + python3 validate_extractors.py ... --extension .segy --output ~/segy-report.csv --resume +""" + +import argparse +import collections +import csv +import datetime as dt +import math +import os +import struct +import sys +import time +import typing + +# numBytesDgm, dgmType, dgmVersion, systemID, echoSounderID, time_sec, time_nanosec +_DATAGRAM_HEADER = struct.Struct(" size + or not dgm_type.startswith(b"#") + ): + if datagram_count == 0: + raise ValueError("does not look like a KMALL file") + break + datagram_count += 1 + versions.setdefault(dgm_type, set()).add(dgm_version) + timestamp = time_sec + time_nanosec / 1e9 + if echo_sounder_id is None: + echo_sounder_id = sounder_id + min_time = timestamp if min_time is None else min(min_time, timestamp) + max_time = timestamp if max_time is None else max(max_time, timestamp) + if dgm_type in _POSITION_DATAGRAMS: + body = fh.read( + min(num_bytes - _DATAGRAM_HEADER.size, _POSITION_BODY_BYTES) + ) + fix = _parse_position_fix(body) + if fix is not None: + lon, lat = fix + fix_counts[dgm_type] += 1 + box = boxes[dgm_type] + boxes[dgm_type] = ( + (lon, lat, lon, lat) + if box is None + else ( + min(box[0], lon), + min(box[1], lat), + max(box[2], lon), + max(box[3], lat), + ) + ) + position += num_bytes + + if datagram_count == 0: + raise ValueError("does not look like a KMALL file") + + bbox = None + position_count = 0 + position_source = "" + for dgm_type in _POSITION_DATAGRAMS: + if boxes[dgm_type] is not None: + bbox = boxes[dgm_type] + position_count = fix_counts[dgm_type] + position_source = dgm_type.decode() + break + + return { + "datagram_count": datagram_count, + "position_count": position_count, + "position_source": position_source, + "echo_sounder_id": echo_sounder_id, + "dgm_versions": ";".join( + f"{k.decode()}:{','.join(str(v) for v in sorted(vals))}" + for k, vals in sorted(versions.items()) + ), + "min_lon": bbox[0] if bbox else "", + "min_lat": bbox[1] if bbox else "", + "max_lon": bbox[2] if bbox else "", + "max_lat": bbox[3] if bbox else "", + "date_begin": _utc_date(min_time), + "date_end": _utc_date(max_time), + } + + +def _parse_position_fix(body): + if len(body) < 2: + return None + (num_bytes_common,) = struct.unpack_from(" len(body): + return None + latitude, longitude = struct.unpack_from(" 90 or abs(longitude) > 180: + return None + return longitude, latitude + + +def _utc_date(timestamp): + if timestamp is None: + return "" + return dt.datetime.fromtimestamp(timestamp, dt.timezone.utc).date().isoformat() + + +# SEG-Y header sampling, mirroring src/seis_lab_data/tasks/extractors/segy.py. + +_TEXTUAL_HEADER_BYTES = 3200 +_BINARY_HEADER_BYTES = 400 +_TRACE_HEADER_BYTES = 240 +_TRACE_DATA_START = _TEXTUAL_HEADER_BYTES + _BINARY_HEADER_BYTES +_MINIMUM_FILE_BYTES = _TRACE_DATA_START + _TRACE_HEADER_BYTES +_MAX_EXTENDED_HEADERS = 100 +_SAMPLE_TARGET = 100 +_GARBAGE_FRACTION_LIMIT = 0.5 +_MAX_PLAUSIBLE_METRES = 1e7 +# see segy.py: a box spanning more than a real survey line means the coordinate +# fields hold noise, so it is discarded rather than trimmed. 200 km catches both +# junk modes the 66k scan found (origin-noise 3,173-17,303 km, and the A7808 +# ET_SBP burst corruption at 494-496 km with garbage 0) with the 62-450 km range +# empty of real files (real max 61.5 km). +_MAX_PLAUSIBLE_SPAN = {"metres": 200_000.0, "geographic": 5.0} +# see segy.py: a box built from a sliver of the sampled traces is dropped and the +# coverage flagged partial, once enough traces were sampled to trust the ratio. +_MIN_SUPPORT_PARSED = 20 +_MIN_SUPPORT_FLOOR = 5 +_MIN_SUPPORT_FRACTION = 0.10 +_INT32_SENTINELS = (-(2**31), 2**31 - 1) + +_SAMPLE_FORMATS = { + 1: ("ibm-float32", 4), + 2: ("int32", 4), + 3: ("int16", 2), + 5: ("ieee-float32", 4), + 6: ("ieee-float64", 8), + 8: ("int8", 1), + 10: ("uint32", 4), + 11: ("uint16", 2), + 16: ("uint8", 1), +} + +_UNITS_ARC_SECONDS = 2 +_UNITS_DEGREES = 3 +_UNITS_DMS = 4 +_UNITS_LABELS = {0: "unset", 1: "metres", 2: "arc-seconds", 3: "degrees", 4: "dms"} + + +class _BinaryHeader(typing.NamedTuple): + sample_interval: int # offset 16 + samples_per_trace: int # offset 20 + sample_format: int # offset 24 + revision: int # offset 300; major revision, normalised on parsing + extended_header_count: int # offset 304 + extra_trace_headers: int # offset 306; rev2 additional 240-byte trace headers + + +_BINARY_HEADER_FORMAT = "16x H 2x H 2x H 274x H 2x h i" +_BIG_ENDIAN_BINARY_HEADER = struct.Struct(">" + _BINARY_HEADER_FORMAT) +_LITTLE_ENDIAN_BINARY_HEADER = struct.Struct("<" + _BINARY_HEADER_FORMAT) + + +class _TraceHeader(typing.NamedTuple): + scalco: int # offset 70 + src_x: int # offset 72 + src_y: int # offset 76 + units: int # offset 88 + year: int # offset 156 + day: int # offset 158 + hour: int # offset 160 + minute: int # offset 162 + second: int # offset 164 + cdp_x: int # offset 180 + cdp_y: int # offset 184 + + +_TRACE_HEADER_FORMAT = "70x h i i 8x H 66x H H H H H 14x i i" +_BIG_ENDIAN_TRACE_HEADER = struct.Struct(">" + _TRACE_HEADER_FORMAT) +_LITTLE_ENDIAN_TRACE_HEADER = struct.Struct("<" + _TRACE_HEADER_FORMAT) + + +def validate_segy(path): + """Mirror of extract_segy_metadata, returning a plain report dict.""" + size = os.path.getsize(path) + if size < _MINIMUM_FILE_BYTES: + raise ValueError("does not look like a SEG-Y file") + with open(path, "rb") as fh: + fh.seek(_TEXTUAL_HEADER_BYTES) + headers = _parse_binary_header(fh.read(_BINARY_HEADER_BYTES)) + if headers is None: + raise ValueError("does not look like a SEG-Y file") + binary_header, trace_header_struct = headers + format_label, bytes_per_sample = _SAMPLE_FORMATS.get( + binary_header.sample_format, + ("format code %d" % binary_header.sample_format, None), + ) + row = { + "sample_interval": binary_header.sample_interval or "", + "samples_per_trace": binary_header.samples_per_trace or "", + "sample_format": format_label, + } + if ( + bytes_per_sample is None + or binary_header.samples_per_trace == 0 + or (binary_header.revision >= 2 and binary_header.extra_trace_headers) + ): + return row + data_start = _locate_trace_data(binary_header, size) + trace_length = ( + _TRACE_HEADER_BYTES + binary_header.samples_per_trace * bytes_per_sample + ) + trace_count = (size - data_start) // trace_length + row["trace_count"] = trace_count + if trace_count == 0: + return row + + boxes = {"metres": None, "geographic": None} + accepted = {"metres": 0, "geographic": 0} + units_counts = collections.Counter() + source_counts = collections.Counter() + parsed_count = 0 + garbage_count = 0 + min_date = max_date = None + for index in _sample_indices(trace_count): + fh.seek(data_start + index * trace_length) + buffer = fh.read(_TRACE_HEADER_BYTES) + if len(buffer) < trace_header_struct.size: + break + trace = _TraceHeader(*trace_header_struct.unpack_from(buffer)) + parsed_count += 1 + units_counts[trace.units] += 1 + date = _parse_trace_date(trace) + if date is not None: + min_date = date if min_date is None else min(min_date, date) + max_date = date if max_date is None else max(max_date, date) + raw = _select_raw_coordinate(trace) + if raw is None or trace.units == _UNITS_DMS: + continue + # report-only: which of the two coordinate sources answered + source_counts["src" if raw == (trace.src_x, trace.src_y) else "cdp"] += 1 + point = _plausible_point(raw, trace.scalco, trace.units) + if point is None: + garbage_count += 1 + continue + category, x, y = point + accepted[category] += 1 + box = boxes[category] + boxes[category] = ( + (x, y, x, y) + if box is None + else (min(box[0], x), min(box[1], y), max(box[2], x), max(box[3], y)) + ) + + for category, box in boxes.items(): + span = _MAX_PLAUSIBLE_SPAN[category] + if box is not None and (box[2] - box[0] > span or box[3] - box[1] > span): + boxes[category] = None + garbage_count += accepted[category] + accepted[category] = 0 + + # Report-only diagnostics, deliberately NOT in segy.py: the extractor keeps + # only its conclusions, but at scale we need the raw tallies behind them - + # a desynced sampling grid lands around the 0.5 garbage fraction, so + # coverage alone hides it, whereas real files score exactly 0. + row["parsed_count"] = parsed_count + row["garbage_count"] = garbage_count + row["accepted_count"] = accepted["metres"] + accepted["geographic"] + row["units_counts"] = ";".join( + "%d:%d" % (code, n) for code, n in sorted(units_counts.items()) + ) + row["position_source"] = ";".join( + "%s:%d" % (name, n) for name, n in sorted(source_counts.items()) + ) + if parsed_count == 0: + return row + dominant_units = units_counts.most_common(1)[0][0] + row["coordinate_units"] = _UNITS_LABELS.get( + dominant_units, "code %d" % dominant_units + ) + # accepted_count above keeps the pre-drop tally for the report; here a sliver + # of support drops the box and flags the coverage, mirroring segy.py. + low_support = parsed_count >= _MIN_SUPPORT_PARSED and 0 < row[ + "accepted_count" + ] < max(_MIN_SUPPORT_FLOOR, parsed_count * _MIN_SUPPORT_FRACTION) + if low_support: + boxes["metres"] = boxes["geographic"] = None + accepted["metres"] = accepted["geographic"] = 0 + row["coverage"] = ( + "partial" + if low_support or garbage_count > parsed_count * _GARBAGE_FRACTION_LIMIT + else "full" + ) + row["date_begin"] = min_date.isoformat() if min_date else "" + row["date_end"] = max_date.isoformat() if max_date else "" + if accepted["geographic"] > accepted["metres"]: + bbox = boxes["geographic"] + row["min_lon"], row["min_lat"], row["max_lon"], row["max_lat"] = bbox + else: + bbox = boxes["metres"] + if bbox is not None: + ( + row["native_min_x"], + row["native_min_y"], + row["native_max_x"], + row["native_max_y"], + ) = bbox + return row + + +def _parse_binary_header(buffer): + if len(buffer) < _BIG_ENDIAN_BINARY_HEADER.size: + return None + for binary_struct, trace_struct in ( + (_BIG_ENDIAN_BINARY_HEADER, _BIG_ENDIAN_TRACE_HEADER), + (_LITTLE_ENDIAN_BINARY_HEADER, _LITTLE_ENDIAN_TRACE_HEADER), + ): + header = _BinaryHeader(*binary_struct.unpack_from(buffer)) + if 1 <= header.sample_format <= 16: + code = header.revision + return ( + header._replace(revision=code >> 8 if code > 0xFF else code), + trace_struct, + ) + return None + + +def _locate_trace_data(binary_header, size): + count = binary_header.extended_header_count + if ( + binary_header.revision >= 1 + and 0 < count <= _MAX_EXTENDED_HEADERS + and _TRACE_DATA_START + count * _TEXTUAL_HEADER_BYTES + _TRACE_HEADER_BYTES + <= size + ): + return _TRACE_DATA_START + count * _TEXTUAL_HEADER_BYTES + return _TRACE_DATA_START + + +def _sample_indices(trace_count): + if trace_count <= _SAMPLE_TARGET: + return list(range(trace_count)) + last = trace_count - 1 + return sorted( + {step * last // (_SAMPLE_TARGET - 1) for step in range(_SAMPLE_TARGET)} + ) + + +def _select_raw_coordinate(trace): + for raw_x, raw_y in ((trace.src_x, trace.src_y), (trace.cdp_x, trace.cdp_y)): + if raw_x in _INT32_SENTINELS or raw_y in _INT32_SENTINELS: + continue + if raw_x == 0 and raw_y == 0: + continue + return raw_x, raw_y + return None + + +def _plausible_point(raw, scalco, units): + x = _apply_coordinate_scalar(raw[0], scalco) + y = _apply_coordinate_scalar(raw[1], scalco) + if not (math.isfinite(x) and math.isfinite(y)): + return None + if units in (_UNITS_ARC_SECONDS, _UNITS_DEGREES): + if units == _UNITS_ARC_SECONDS: + x /= 3600.0 + y /= 3600.0 + if abs(x) > 180.0 or abs(y) > 90.0: + return None + return "geographic", x, y + if abs(x) >= _MAX_PLAUSIBLE_METRES or abs(y) >= _MAX_PLAUSIBLE_METRES: + return None + return "metres", x, y + + +def _apply_coordinate_scalar(value, scalco): + if scalco < 0: + return value / -scalco + if scalco > 0: + return float(value * scalco) + return float(value) + + +def _parse_trace_date(trace): + if not 1970 <= trace.year <= 2100: + return None + if not 1 <= trace.day <= 366: + return None + if trace.hour > 23 or trace.minute > 59 or trace.second > 59: + return None + date = dt.date(trace.year, 1, 1) + dt.timedelta(days=trace.day - 1) + if date.year != trace.year: + return None + return date + + +_VALIDATORS = { + ".kmall": validate_kmall, + ".segy": validate_segy, + ".sgy": validate_segy, +} + + +def iter_target_files(root, extension, file_list): + if file_list: + with open(file_list) as fh: + for line in fh: + candidate = line.rstrip("\n").split("\t")[-1] + if candidate.lower().endswith(extension): + yield os.path.join(root, candidate.removeprefix("./")) + else: + for dirpath, _dirnames, filenames in os.walk(root): + for name in sorted(filenames): + if name.lower().endswith(extension): + yield os.path.join(dirpath, name) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--root", required=True, help="archive root directory") + parser.add_argument("--extension", default=".kmall", choices=sorted(_VALIDATORS)) + parser.add_argument("--output", required=True, help="CSV report path") + parser.add_argument( + "--file-list", + help="optional listing (plain paths, or sizepath); paths are " + "resolved against --root; when absent the root is walked", + ) + parser.add_argument( + "--resume", + action="store_true", + help="skip files already present in the output CSV", + ) + parser.add_argument("--limit", type=int, help="stop after N files (trial runs)") + args = parser.parse_args() + + validator = _VALIDATORS[args.extension] + done = set() + has_report = os.path.exists(args.output) and os.path.getsize(args.output) > 0 + if has_report: + # SEG-Y needs one pass per extension (.sgy then .segy) into the same + # report, so refuse to silently truncate hours of work: --resume both + # keeps the existing rows and skips the files already in them. + if not args.resume: + sys.exit( + f"{args.output} already holds a report; pass --resume to add to " + "it (or choose another --output)" + ) + with open(args.output) as fh: + done = {row["path"] for row in csv.DictReader(fh)} + + mode = "a" if has_report else "w" + processed = ok = failed = 0 + with open(args.output, mode, newline="") as out: + writer = csv.DictWriter(out, fieldnames=_CSV_COLUMNS) + if mode == "w": + writer.writeheader() + for path in iter_target_files(args.root, args.extension, args.file_list): + relative = os.path.relpath(path, args.root) + if relative in done: + continue + row = {"path": relative} + started = time.perf_counter() + try: + row["size_bytes"] = os.path.getsize(path) + row.update(validator(path)) + row["status"] = "ok" + ok += 1 + except Exception as err: # report and carry on - this is a survey + row["status"] = "error" + row["error"] = f"{type(err).__name__}: {err}" + failed += 1 + row["walk_seconds"] = f"{time.perf_counter() - started:.2f}" + writer.writerow(row) + out.flush() + processed += 1 + if processed % 50 == 0: + print( + f"{processed} files ({ok} ok, {failed} errors)...", + file=sys.stderr, + flush=True, + ) + if args.limit and processed >= args.limit: + break + + print(f"done: {processed} files, {ok} ok, {failed} errors", file=sys.stderr) + + +if __name__ == "__main__": + main() diff --git a/src/seis_lab_data/db/queries/recordassets.py b/src/seis_lab_data/db/queries/recordassets.py index 8a5b2c83..a5bee106 100644 --- a/src/seis_lab_data/db/queries/recordassets.py +++ b/src/seis_lab_data/db/queries/recordassets.py @@ -71,11 +71,21 @@ async def get_record_asset_by_english_name( async def get_record_asset_by_file_path( - session: AsyncSession, file_path: str + session: AsyncSession, + file_path: str, + survey_mission_id: identifiers.SurveyMissionId, ) -> models.RecordAsset | None: + # Scoped per mission: the same relative path may legitimately exist in + # several missions, each deserving its own record. statement = ( select(models.RecordAsset) + .join( + models.SurveyRelatedRecord, + models.RecordAsset.survey_related_record_id + == models.SurveyRelatedRecord.id, + ) .where(models.RecordAsset.relative_path == file_path) + .where(models.SurveyRelatedRecord.survey_mission_id == survey_mission_id) .options(_SELECT_IN_LOAD_OPTIONS) ) return (await session.exec(statement)).first() diff --git a/src/seis_lab_data/operations/discovery.py b/src/seis_lab_data/operations/discovery.py index 8dd24325..a2a01179 100644 --- a/src/seis_lab_data/operations/discovery.py +++ b/src/seis_lab_data/operations/discovery.py @@ -1,10 +1,12 @@ import logging +import math import re import uuid from collections.abc import AsyncIterator from typing import AsyncGenerator -from anyio import Path +import shapely +from anyio import Path, to_thread from sqlmodel.ext.asyncio.session import AsyncSession from .. import ( @@ -32,6 +34,7 @@ user as user_schemas, ) from .. import dispatch +from ..tasks.extractors import dispatch as extractor_dispatch from . import ( surveymissions as mission_ops, @@ -40,6 +43,11 @@ logger = logging.getLogger(__name__) +# Buffer applied to every extracted bbox (~10 m in degrees): noise on multi-km +# footprints, but keeps point/line extents visible on the map and avoids the +# degenerate zero-area polygons that shapely marks invalid. +_BBOX_BUFFER = 1e-4 + async def create_asset_discovery_configuration( *, @@ -316,6 +324,17 @@ async def _discover_mission_records( ) logger.debug(f"{mission_root_path=}") for asset_discovery_conf in asset_discovery_configs: + if ( + asset_discovery_conf.dataset_category_id is None + or asset_discovery_conf.workflow_stage_id is None + ): + logger.warning( + "Skipping asset discovery configuration %s (%r): it has no " + "dataset category or workflow stage", + asset_discovery_conf.id, + asset_discovery_conf.name, + ) + continue logger.debug(f"Searching for asset {asset_discovery_conf=}...") full_path_regexp = "/".join( ( @@ -329,9 +348,28 @@ async def _discover_mission_records( relative_file_path = str(found_path.relative_to(mission_root_path)) if ( await asset_queries.get_record_asset_by_file_path( - session, relative_file_path + session, + relative_file_path, + identifiers.SurveyMissionId(mission.id), ) ) is None: + # best-effort metadata extraction: a failure must never abort + # discovery or record creation + metadata = None + try: + metadata = await to_thread.run_sync( + extractor_dispatch.dispatch_extractor, str(found_path) + ) + except Exception as err: + logger.warning( + "Metadata extraction failed for %s: %s", found_path, err + ) + logger.debug("Extraction failure detail", exc_info=True) + bbox_wkt = ( + _bbox_4326_tuple_to_wkt(metadata.bbox_4326) + if metadata is not None and metadata.bbox_4326 is not None + else None + ) # create a new record and a new asset await record_ops.create_survey_related_record( request_id=request_id, @@ -340,13 +378,22 @@ async def _discover_mission_records( owner_id=identifiers.UserId(user.id), survey_mission_id=identifiers.SurveyMissionId(mission.id), name=common.LocalizableDraftName(en=found_path.name), - description=common.LocalizableDraftDescription(en=""), + description=common.LocalizableDraftDescription( + en=metadata.describe() if metadata is not None else "" + ), dataset_category_id=identifiers.DatasetCategoryId( asset_discovery_conf.dataset_category_id ), workflow_stage_id=identifiers.WorkflowStageId( asset_discovery_conf.workflow_stage_id ), + bbox_4326=bbox_wkt, + temporal_extent_begin=( + metadata.temporal_extent_begin if metadata else None + ), + temporal_extent_end=( + metadata.temporal_extent_end if metadata else None + ), assets=[ record_schemas.RecordAssetCreate( id=identifiers.RecordAssetId(uuid.uuid4()), @@ -366,6 +413,38 @@ async def _discover_mission_records( ) +def _bbox_4326_tuple_to_wkt( + bbox: tuple[float, float, float, float], +) -> str | None: + """Turn an extracted lon/lat bbox into a WKT polygon for the record. + + Garbage coordinates (non-finite or out of lon/lat range, e.g. from a bogus + CRS transform) are discarded with a warning. The surviving bbox is expanded + by _BBOX_BUFFER on every side. Antimeridian-crossing bboxes are not handled + (Portuguese waters). + """ + minx, miny, maxx, maxy = bbox + tolerance = 1e-6 + if not all(math.isfinite(value) for value in bbox) or not ( + -180 - tolerance <= minx <= maxx <= 180 + tolerance + and -90 - tolerance <= miny <= maxy <= 90 + tolerance + ): + logger.warning( + "Discarding extracted bbox with out-of-range or non-finite coordinates: %s", + bbox, + ) + return None + # The envelope of the two corners is always valid (Point, LineString or + # Polygon), unlike shapely.box on a degenerate bbox; buffering with square + # caps and mitre joins keeps the result an axis-aligned rectangle. + envelope = shapely.MultiPoint([(minx, miny), (maxx, maxy)]).envelope + buffered = envelope.buffer(_BBOX_BUFFER, cap_style="square", join_style="mitre") + # Snap to a 1e-5 degree grid (~1 m): finer precision is meaningless under + # the ~10 m buffer, and the frontend's TerraDraw silently rejects + # coordinates with more than 9 decimal places, dropping the map rectangle. + return shapely.set_precision(buffered, 1e-5).wkt + + async def _discover_asset_paths( full_path_regexp: str, ) -> AsyncGenerator[Path, None]: diff --git a/src/seis_lab_data/tasks/extractors/__init__.py b/src/seis_lab_data/tasks/extractors/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/src/seis_lab_data/tasks/extractors/common.py b/src/seis_lab_data/tasks/extractors/common.py new file mode 100644 index 00000000..30db98a7 --- /dev/null +++ b/src/seis_lab_data/tasks/extractors/common.py @@ -0,0 +1,52 @@ +import logging + +from osgeo import gdal, osr + +logger = logging.getLogger(__name__) + +gdal.UseExceptions() + + +def identify_srs( + srs: "osr.SpatialReference | None", +) -> tuple[int | None, str | None, str | None]: + """Best-effort (epsg, crs_name, crs_wkt) for a spatial reference. + + AutoIdentifyEPSG fails on some real archive rasters (e.g. ETRS89/PT-TM06), + so epsg may stay None while crs_name/crs_wkt are still populated. + """ + if srs is None: + return None, None, None + try: + srs.AutoIdentifyEPSG() + except RuntimeError as err: + logger.debug("AutoIdentifyEPSG failed: %s", err) + # Only trust the code when the authority is EPSG and the code is numeric: + # an ESRI code (e.g. ESRI:102629) would be mislabelled as EPSG, and a + # non-numeric code (e.g. IGNF:LAMB93) would crash int(). + code = srs.GetAuthorityCode(None) + if srs.GetAuthorityName(None) == "EPSG" and code and code.isdigit(): + epsg = int(code) + else: + epsg = None + return epsg, srs.GetName(), srs.ExportToWkt() + + +def project_bbox_to_wgs84( + bbox_native: tuple[float, float, float, float] | None, + src_srs: "osr.SpatialReference | None", +) -> tuple[float, float, float, float] | None: + if bbox_native is None or src_srs is None: + return None + wgs84 = osr.SpatialReference() + wgs84.ImportFromEPSG(4326) + wgs84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) + src = src_srs.Clone() + src.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) + transform = osr.CoordinateTransformation(src, wgs84) + minx, miny, maxx, maxy = bbox_native + corners = [(minx, miny), (minx, maxy), (maxx, miny), (maxx, maxy)] + projected = [transform.TransformPoint(x, y) for x, y in corners] + lons = [p[0] for p in projected] + lats = [p[1] for p in projected] + return (min(lons), min(lats), max(lons), max(lats)) diff --git a/src/seis_lab_data/tasks/extractors/dispatch.py b/src/seis_lab_data/tasks/extractors/dispatch.py index b1eb3fb5..dce1fb0d 100644 --- a/src/seis_lab_data/tasks/extractors/dispatch.py +++ b/src/seis_lab_data/tasks/extractors/dispatch.py @@ -1,14 +1,64 @@ +import logging from pathlib import Path -from .schemas import RasterMetadata +from .gdal_raster import extract_raster_metadata +from .gdal_vector import extract_vector_metadata +from .kmall import extract_kmall_metadata +from .schemas import ExtractionResult +from .segy import extract_segy_metadata -_RASTER_EXTENSIONS = frozenset({".tif", ".tiff", ".xyz", ".asc", ".flt", ".nc", ".nc4"}) +logger = logging.getLogger(__name__) +_RASTER_EXTENSIONS = frozenset( + {".tif", ".tiff", ".asc", ".flt", ".grd", ".xyz", ".nc", ".nc4"} +) +# .nc/.nc4 are kept for gridded NetCDF (served by the plain GDAL raster path); the +# ubuntu-small image has no netCDF driver yet, so extraction fails best-effort until the +# base image moves to gdal:ubuntu-full. Production has 0 .nc files today. -def dispatch_extractor(path: Path | str) -> RasterMetadata | None: - p = Path(path) - if p.suffix.lower() in _RASTER_EXTENSIONS: - from .gdal_raster import extract_raster_metadata +# Sidecars (.shx/.dbf/.prj/.cpg/.qmd/.sbn/.sbx ...) are deliberately absent so they never +# trigger extraction; GDAL reads them implicitly when opening the primary file. .kmz is +# excluded (needs LIBKML, absent from the image); .gdb (a directory) and S-57 .000 too. +_VECTOR_EXTENSIONS = frozenset({".shp", ".gpkg", ".geojson", ".kml", ".dxf"}) + +_KMALL_EXTENSIONS = frozenset({".kmall"}) + +# No big-file log for SEG-Y: header sampling reads ~100 trace headers no matter +# the file size, so even a 200 GB file costs a constant number of reads. +_SEGY_EXTENSIONS = frozenset({".sgy", ".segy"}) + +_BIG_FILE_LOG_BYTES = 1024**3 # log-only threshold; no hard size limit + +def dispatch_extractor(path: Path | str) -> ExtractionResult | None: + """Route a file to its metadata extractor by extension. + + Pure sync and potentially slow: GDAL's XYZ driver scans the whole file on open, + so a multi-GB grid can take ~1 minute, and KMALL files get a full datagram-header + walk (seconds per GB); SEG-Y files stay fast at any size (constant number of + header reads). Async callers must run this in a worker thread (e.g. + anyio.to_thread.run_sync). Returns None for unsupported extensions and + directories. + """ + p = Path(path) + if not p.is_file(): + # A real directory named "F3_2022.tif" exists in the archive. + return None + suffix = p.suffix.lower() + if suffix in _RASTER_EXTENSIONS: + size = p.stat().st_size + if size > _BIG_FILE_LOG_BYTES: + logger.info("Extracting metadata from large raster %s (%d bytes)", p, size) return extract_raster_metadata(p) + if suffix in _VECTOR_EXTENSIONS: + return extract_vector_metadata(p) + if suffix in _KMALL_EXTENSIONS: + size = p.stat().st_size + if size > _BIG_FILE_LOG_BYTES: + logger.info( + "Extracting metadata from large KMALL file %s (%d bytes)", p, size + ) + return extract_kmall_metadata(p) + if suffix in _SEGY_EXTENSIONS: + return extract_segy_metadata(p) return None diff --git a/src/seis_lab_data/tasks/extractors/gdal_raster.py b/src/seis_lab_data/tasks/extractors/gdal_raster.py index 8a015f1c..8a333550 100644 --- a/src/seis_lab_data/tasks/extractors/gdal_raster.py +++ b/src/seis_lab_data/tasks/extractors/gdal_raster.py @@ -3,6 +3,7 @@ from osgeo import gdal, osr +from . import common from .schemas import RasterMetadata logger = logging.getLogger(__name__) @@ -25,17 +26,18 @@ def extract_raster_metadata(path: Path | str) -> RasterMetadata: else: pixel_size_x = gt[1] pixel_size_y = gt[5] - x0 = gt[0] - y0 = gt[3] - x1 = x0 + width * gt[1] + height * gt[2] - y1 = y0 + width * gt[4] + height * gt[5] - minx, maxx = sorted((x0, x1)) - miny, maxy = sorted((y0, y1)) - bbox_native = (minx, miny, maxx, maxy) - - epsg, crs_wkt, src_srs = _read_srs(ds) + # Transform all four pixel corners so rotated/sheared geotransforms + # (gt[2] or gt[4] nonzero) get a correct bbox, not just two corners. + corners_px = [(0, 0), (width, 0), (0, height), (width, height)] + xs = [gt[0] + px * gt[1] + py * gt[2] for px, py in corners_px] + ys = [gt[3] + px * gt[4] + py * gt[5] for px, py in corners_px] + bbox_native = (min(xs), min(ys), max(xs), max(ys)) + + wkt = ds.GetProjection() + src_srs = osr.SpatialReference(wkt) if wkt else None + epsg, crs_name, crs_wkt = common.identify_srs(src_srs) nodata = ds.GetRasterBand(1).GetNoDataValue() if band_count else None - bbox_4326 = _project_to_wgs84(bbox_native, src_srs) + bbox_4326 = common.project_bbox_to_wgs84(bbox_native, src_srs) return RasterMetadata( driver=driver, @@ -43,6 +45,7 @@ def extract_raster_metadata(path: Path | str) -> RasterMetadata: height=height, band_count=band_count, epsg=epsg, + crs_name=crs_name, crs_wkt=crs_wkt, pixel_size_x=pixel_size_x, pixel_size_y=pixel_size_y, @@ -52,118 +55,3 @@ def extract_raster_metadata(path: Path | str) -> RasterMetadata: ) finally: ds = None # noqa: F841 - - -def _read_srs(ds) -> tuple[int | None, str | None, "osr.SpatialReference | None"]: - wkt = ds.GetProjection() - if not wkt: - return None, None, None - srs = osr.SpatialReference(wkt) - try: - srs.AutoIdentifyEPSG() - except RuntimeError as err: - logger.debug("AutoIdentifyEPSG failed for raster: %s", err) - code = srs.GetAuthorityCode(None) - epsg = int(code) if code else None - return epsg, wkt, srs - - -def _project_to_wgs84( - bbox_native: tuple[float, float, float, float] | None, - src_srs: "osr.SpatialReference | None", -) -> tuple[float, float, float, float] | None: - if bbox_native is None or src_srs is None: - return None - wgs84 = osr.SpatialReference() - wgs84.ImportFromEPSG(4326) - wgs84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) - src = src_srs.Clone() - src.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) - transform = osr.CoordinateTransformation(src, wgs84) - minx, miny, maxx, maxy = bbox_native - corners = [(minx, miny), (minx, maxy), (maxx, miny), (maxx, maxy)] - projected = [transform.TransformPoint(x, y) for x, y in corners] - lons = [p[0] for p in projected] - lats = [p[1] for p in projected] - return (min(lons), min(lats), max(lons), max(lats)) - - -# todo reverificar se está bem aqui ou não. Como isto passou do discovery_demo/protocols.py ver se não quebrei qualquer coisa. - - -# async def extractor_gdal_raster( -# *, -# survey_mission: models.SurveyMission, -# record_configuration: discovery_schemas.SurveyRecordDiscoveryConfiguration, -# discovered_assets: list[DiscoveredAsset], -# session: AsyncSession, -# ) -> record_schemas.SurveyRelatedRecordCreate: -# if not discovered_assets: -# raise ValueError("extractor1 requires at least one discovered asset") - -# bboxes: list[tuple[float, float, float, float]] = [] -# for _, abs_path in discovered_assets: -# metadata = dispatch_extractor(abs_path) -# if metadata is not None and metadata.bbox_4326 is not None: -# bboxes.append(metadata.bbox_4326) - -# bbox_wkt = _aggregate_bbox_wkt(bboxes) - -# dataset_category = await record_queries.get_dataset_category_by_english_name( -# session, record_configuration.dataset_category -# ) -# domain_type = await record_queries.get_domain_type_by_english_name( -# session, record_configuration.domain_type -# ) -# workflow_stage = await record_queries.get_workflow_stage_by_english_name( -# session, record_configuration.workflow_stage -# ) -# if dataset_category is None: -# raise ValueError( -# f"Unknown dataset_category {record_configuration.dataset_category!r}" -# ) -# if domain_type is None: -# raise ValueError(f"Unknown domain_type {record_configuration.domain_type!r}") -# if workflow_stage is None: -# raise ValueError( -# f"Unknown workflow_stage {record_configuration.workflow_stage!r}" -# ) - -# relative_path = "/".join( -# ( -# record_configuration.domain_type, -# record_configuration.dataset_category, -# record_configuration.workflow_stage, -# ) -# ) - -# return SurveyRelatedRecordCreate( -# id=SurveyRelatedRecordId(uuid.uuid4()), -# owner_id=UserId(survey_mission.owner_id), -# survey_mission_id=SurveyMissionId(survey_mission.id), -# name=LocalizableDraftName(**record_configuration.name), -# description=LocalizableDraftDescription( -# **(record_configuration.description or {}) -# ), -# dataset_category_id=DatasetCategoryId(dataset_category.id), -# domain_type_id=DomainTypeId(domain_type.id), -# workflow_stage_id=WorkflowStageId(workflow_stage.id), -# relative_path=relative_path, -# bbox_4326=bbox_wkt, -# temporal_extent_begin=None, -# temporal_extent_end=None, -# links=list(record_configuration.links), -# assets=[asset for asset, _ in discovered_assets], -# related_records=[], -# ) - - -# def _aggregate_bbox_wkt( -# bboxes: list[tuple[float, float, float, float]], -# ) -> str | None: -# if not bboxes: -# return None -# polys = [shapely.box(minx, miny, maxx, maxy) for minx, miny, maxx, maxy in bboxes] -# union = shapely.unary_union(polys) -# minx, miny, maxx, maxy = union.bounds -# return shapely.box(minx, miny, maxx, maxy).wkt diff --git a/src/seis_lab_data/tasks/extractors/gdal_vector.py b/src/seis_lab_data/tasks/extractors/gdal_vector.py new file mode 100644 index 00000000..8e0f4730 --- /dev/null +++ b/src/seis_lab_data/tasks/extractors/gdal_vector.py @@ -0,0 +1,84 @@ +import logging +from pathlib import Path + +from osgeo import gdal, ogr + +from . import common +from .schemas import VectorMetadata + +logger = logging.getLogger(__name__) + +gdal.UseExceptions() + + +def extract_vector_metadata(path: Path | str) -> VectorMetadata: + ds = gdal.OpenEx(str(path), gdal.OF_VECTOR | gdal.OF_READONLY) + try: + driver = ds.GetDriver().ShortName + layer_count = ds.GetLayerCount() + + feature_count = 0 + geometry_names: set[str] = set() + native_boxes: list[tuple[float, float, float, float]] = [] + native_srs_names: set[str | None] = set() + boxes_4326: list[tuple[float, float, float, float]] = [] + epsg = crs_name = crs_wkt = None + srs_captured = False + + for layer_index in range(layer_count): + lyr = ds.GetLayer(layer_index) + count = max(lyr.GetFeatureCount(), 0) + feature_count += count + geometry_names.add(ogr.GeometryTypeToName(lyr.GetGeomType())) + + if count == 0: + continue + try: + # OGR GetExtent returns (minx, maxx, miny, maxy); reorder it. + minx, maxx, miny, maxy = lyr.GetExtent(force=1) + except RuntimeError: + continue + # Capture the CRS from an extent-bearing layer so the reported CRS + # matches bbox_native's CRS (a 0-feature first layer must not set it). + srs = lyr.GetSpatialRef() + if not srs_captured and srs is not None: + epsg, crs_name, crs_wkt = common.identify_srs(srs) + srs_captured = True + native_boxes.append((minx, miny, maxx, maxy)) + native_srs_names.add(srs.GetName() if srs is not None else None) + projected = common.project_bbox_to_wgs84((minx, miny, maxx, maxy), srs) + if projected is not None: + boxes_4326.append(projected) + + # Native bbox only makes sense when the extent-bearing layers share one CRS. + bbox_native = _union(native_boxes) if len(native_srs_names) <= 1 else None + geometry_type = ( + next(iter(geometry_names)) if len(geometry_names) == 1 else "mixed" + ) + + return VectorMetadata( + driver=driver, + epsg=epsg, + crs_name=crs_name, + crs_wkt=crs_wkt, + bbox_native=bbox_native, + bbox_4326=_union(boxes_4326), + layer_count=layer_count, + feature_count=feature_count, + geometry_type=geometry_type, + ) + finally: + ds = None # noqa: F841 + + +def _union( + boxes: list[tuple[float, float, float, float]], +) -> tuple[float, float, float, float] | None: + if not boxes: + return None + return ( + min(b[0] for b in boxes), + min(b[1] for b in boxes), + max(b[2] for b in boxes), + max(b[3] for b in boxes), + ) diff --git a/src/seis_lab_data/tasks/extractors/kmall.py b/src/seis_lab_data/tasks/extractors/kmall.py new file mode 100644 index 00000000..f0a59bae --- /dev/null +++ b/src/seis_lab_data/tasks/extractors/kmall.py @@ -0,0 +1,156 @@ +import datetime as dt +import logging +import math +import struct +from pathlib import Path + +from .schemas import KmallMetadata + +logger = logging.getLogger(__name__) + +# The standalone prod validator scripts/validate_extractors.py re-implements this +# walk (it cannot import the package). Any change to the datagram/fix parsing +# below must be mirrored there, and vice versa. + +# numBytesDgm, dgmType, dgmVersion, systemID, echoSounderID, time_sec, time_nanosec +_DATAGRAM_HEADER = struct.Struct(" KmallMetadata: + """Extract metadata from a Kongsberg KMALL file by walking datagram headers. + + Reads only the 20-byte datagram headers plus the head of each position + datagram - payloads are never loaded, so multi-GB files take seconds. A + truncated tail (e.g. an acquisition crash) keeps the metadata gathered up + to that point; a file whose first datagram is invalid raises ValueError. + """ + p = Path(path) + size = p.stat().st_size + datagram_count = 0 + # running min/max, not first/last: real files interleave sensor datagrams + # with header times backstepping up to ~3 s, which could shift a date at a + # UTC midnight boundary + min_time = max_time = None + echo_sounder_id = None + boxes: dict[bytes, tuple[float, float, float, float] | None] = { + dgm_type: None for dgm_type in _POSITION_DATAGRAMS + } + fix_counts = {dgm_type: 0 for dgm_type in _POSITION_DATAGRAMS} + + with p.open("rb") as fh: + position = 0 + while position + _DATAGRAM_HEADER.size <= size: + fh.seek(position) + ( + num_bytes, + dgm_type, + _dgm_version, + _system_id, + sounder_id, + time_sec, + time_nanosec, + ) = _DATAGRAM_HEADER.unpack(fh.read(_DATAGRAM_HEADER.size)) + if ( + num_bytes < _MIN_DATAGRAM_BYTES + or position + num_bytes > size + or not dgm_type.startswith(b"#") + ): + if datagram_count == 0: + raise ValueError(f"{p.name} does not look like a KMALL file") + logger.debug( + "Truncated or corrupt datagram at offset %d in %s - keeping " + "the metadata gathered so far", + position, + p, + ) + break + datagram_count += 1 + timestamp = time_sec + time_nanosec / 1e9 + if echo_sounder_id is None: + echo_sounder_id = sounder_id + min_time = timestamp if min_time is None else min(min_time, timestamp) + max_time = timestamp if max_time is None else max(max_time, timestamp) + if dgm_type in _POSITION_DATAGRAMS: + body = fh.read( + min(num_bytes - _DATAGRAM_HEADER.size, _POSITION_BODY_BYTES) + ) + fix = _parse_position_fix(body) + if fix is not None: + lon, lat = fix + fix_counts[dgm_type] += 1 + box = boxes[dgm_type] + boxes[dgm_type] = ( + (lon, lat, lon, lat) + if box is None + else ( + min(box[0], lon), + min(box[1], lat), + max(box[2], lon), + max(box[3], lat), + ) + ) + position += num_bytes + + if datagram_count == 0: + raise ValueError(f"{p.name} does not look like a KMALL file") + + bbox = None + position_count = 0 + for dgm_type in _POSITION_DATAGRAMS: + if boxes[dgm_type] is not None: + bbox = boxes[dgm_type] + position_count = fix_counts[dgm_type] + break + + return KmallMetadata( + driver="KMALL", + # positions are geographic degrees from the GNSS feed; the format + # declares no datum, WGS84-family in practice (sub-metre difference, + # far below the bbox buffer) - documented approximation + epsg=4326 if bbox is not None else None, + crs_name="WGS 84" if bbox is not None else None, + bbox_native=bbox, + bbox_4326=bbox, + temporal_extent_begin=_utc_date(min_time), + temporal_extent_end=_utc_date(max_time), + echo_sounder_id=echo_sounder_id, + datagram_count=datagram_count, + position_count=position_count, + ) + + +def _parse_position_fix(body: bytes) -> tuple[float, float] | None: + """Best-effort (lon, lat) from a #SPO/#CPO body; None when unusable.""" + if len(body) < 2: + return None + (num_bytes_common,) = struct.unpack_from(" len(body): + return None + latitude, longitude = struct.unpack_from(" 90 or abs(longitude) > 180: + # the spec uses out-of-range sentinels (e.g. 200.0) for "unavailable" + return None + return longitude, latitude + + +def _utc_date(timestamp: float | None) -> dt.date | None: + if timestamp is None: + return None + return dt.datetime.fromtimestamp(timestamp, dt.timezone.utc).date() diff --git a/src/seis_lab_data/tasks/extractors/protocols.py b/src/seis_lab_data/tasks/extractors/protocols.py deleted file mode 100644 index 54b90533..00000000 --- a/src/seis_lab_data/tasks/extractors/protocols.py +++ /dev/null @@ -1,13 +0,0 @@ -from typing import Protocol - -from ...schemas.discovery import SurveyRecordDiscoveryConfiguration -from ...schemas.surveyrelatedrecords import SurveyRelatedRecordCreate - - -class NewRecordExtractorProtocol(Protocol): - async def __call__( - self, - initial_record: SurveyRelatedRecordCreate, - record_configuration: SurveyRecordDiscoveryConfiguration, - ) -> SurveyRelatedRecordCreate: - """Extract relevant metadata from the record's assets""" diff --git a/src/seis_lab_data/tasks/extractors/schemas.py b/src/seis_lab_data/tasks/extractors/schemas.py index 38f8949a..4afc15f8 100644 --- a/src/seis_lab_data/tasks/extractors/schemas.py +++ b/src/seis_lab_data/tasks/extractors/schemas.py @@ -1,15 +1,127 @@ +import datetime as dt + import pydantic +from ... import constants + + +def _format_crs(epsg: int | None, crs_name: str | None) -> str: + if crs_name and epsg is not None: + return f"{crs_name} (EPSG:{epsg})" + if epsg is not None: + return f"EPSG:{epsg}" + if crs_name: + return crs_name + return "unknown" + -class RasterMetadata(pydantic.BaseModel): - driver: str +def _native_bbox_clause(bbox: tuple[float, float, float, float] | None) -> str: + if bbox is None: + return "" + minx, miny, maxx, maxy = bbox + return f" Native bbox: {minx:.8g}, {miny:.8g}, {maxx:.8g}, {maxy:.8g}." + + +class ExtractedMetadata(pydantic.BaseModel): + driver: str | None = None + epsg: int | None = None + crs_name: str | None = None + crs_wkt: str | None = None + bbox_native: tuple[float, float, float, float] | None = None + bbox_4326: tuple[float, float, float, float] | None = None + temporal_extent_begin: dt.date | None = None + temporal_extent_end: dt.date | None = None + + def describe(self) -> str: + raise NotImplementedError + + +class RasterMetadata(ExtractedMetadata): width: int height: int band_count: int - epsg: int | None = None - crs_wkt: str | None = None pixel_size_x: float | None = None pixel_size_y: float | None = None nodata: float | None = None - bbox_native: tuple[float, float, float, float] | None = None - bbox_4326: tuple[float, float, float, float] | None = None + + def describe(self) -> str: + summary = ( + f"Auto-extracted: {self.driver} raster, " + f"{self.width} x {self.height} px, {self.band_count} band(s)" + ) + if self.pixel_size_x is not None and self.pixel_size_y is not None: + summary += ( + f", pixel size {abs(self.pixel_size_x):g} x {abs(self.pixel_size_y):g}" + ) + if self.nodata is not None: + summary += f", nodata {self.nodata:g}" + summary += f". CRS: {_format_crs(self.epsg, self.crs_name)}." + summary += _native_bbox_clause(self.bbox_native) + return summary[: constants.DESCRIPTION_MAX_LENGTH] + + +class VectorMetadata(ExtractedMetadata): + layer_count: int + feature_count: int + geometry_type: str + + def describe(self) -> str: + summary = ( + f"Auto-extracted: {self.driver} vector, {self.layer_count} layer(s), " + f"{self.feature_count} feature(s), geometry: {self.geometry_type}. " + f"CRS: {_format_crs(self.epsg, self.crs_name)}." + ) + summary += _native_bbox_clause(self.bbox_native) + return summary[: constants.DESCRIPTION_MAX_LENGTH] + + +class KmallMetadata(ExtractedMetadata): + echo_sounder_id: int | None = None + datagram_count: int + position_count: int + + def describe(self) -> str: + sounder = ( + f"EM {self.echo_sounder_id}" if self.echo_sounder_id else "unknown sounder" + ) + summary = ( + f"Auto-extracted: Kongsberg KMALL, {sounder}, " + f"{self.datagram_count} datagram(s), " + f"{self.position_count} position fix(es). " + f"CRS: {_format_crs(self.epsg, self.crs_name)}." + ) + summary += _native_bbox_clause(self.bbox_native) + return summary[: constants.DESCRIPTION_MAX_LENGTH] + + +class SegyMetadata(ExtractedMetadata): + trace_count: int | None = None + samples_per_trace: int | None = None + # the raw header value, deliberately unitless: time files store microseconds + # here but depth-domain files store a length (50 == 0.05 m) and rev0/1 has + # no flag to tell them apart + sample_interval: int | None = None + sample_format: str | None = None + coordinate_units: str | None = None + coverage: str = "full" + + def describe(self) -> str: + summary = "Auto-extracted: SEG-Y" + if self.trace_count is not None: + summary += f", {self.trace_count} trace(s)" + if self.samples_per_trace is not None: + summary += f", {self.samples_per_trace} sample(s) per trace" + if self.sample_interval is not None: + summary += f", sample interval {self.sample_interval} (raw)" + if self.sample_format is not None: + summary += f", {self.sample_format}" + if self.coordinate_units is not None: + summary += f", coordinates in {self.coordinate_units}" + summary += f". CRS: {_format_crs(self.epsg, self.crs_name)}." + summary += _native_bbox_clause(self.bbox_native) + if self.coverage == "partial": + summary += " Coordinate sampling was partly unreliable." + return summary[: constants.DESCRIPTION_MAX_LENGTH] + + +ExtractionResult = RasterMetadata | VectorMetadata | KmallMetadata | SegyMetadata diff --git a/src/seis_lab_data/tasks/extractors/segy.py b/src/seis_lab_data/tasks/extractors/segy.py new file mode 100644 index 00000000..03fecda9 --- /dev/null +++ b/src/seis_lab_data/tasks/extractors/segy.py @@ -0,0 +1,377 @@ +import collections +import datetime as dt +import math +import struct +import typing +from pathlib import Path + +from .schemas import SegyMetadata + +# The standalone prod validator scripts/validate_extractors.py re-implements this +# sampling walk (it cannot import the package). Any change to the header parsing +# below must be mirrored there, and vice versa. + +_TEXTUAL_HEADER_BYTES = 3200 +_BINARY_HEADER_BYTES = 400 +_TRACE_HEADER_BYTES = 240 +_TRACE_DATA_START = _TEXTUAL_HEADER_BYTES + _BINARY_HEADER_BYTES +# A plausible SEG-Y holds at least the two file headers plus one trace header. +_MINIMUM_FILE_BYTES = _TRACE_DATA_START + _TRACE_HEADER_BYTES +# Extended textual header counts beyond this are treated as garbage. +_MAX_EXTENDED_HEADERS = 100 +# How many trace headers to sample per file; the first and last trace are always +# included, so the cost is constant regardless of file size. +_SAMPLE_TARGET = 100 +# Above this fraction of implausible sampled coordinates the file is flagged as +# partially covered - the tell of a variable-length file whose size happened to +# floor-divide evenly, leaving samples landing mid-trace. +_GARBAGE_FRACTION_LIMIT = 0.5 +# Projected coordinates beyond this magnitude are implausible in any real CRS. +_MAX_PLAUSIBLE_METRES = 1e7 +# A file whose sampled positions span more than one survey line possibly could +# has no usable navigation. Two junk modes are known from the 66k-file archive +# scan: coordinate fields holding noise around the projection origin (3,173 to +# 17,303 km across, the A7808 GEOM family), and a burst of corrupt-nav traces +# inside an otherwise real line (the A7808 ET_SBP files reach 494 to 496 km with +# garbage_count 0, so ONLY the span catches them). Real archive files measure +# 14 m to 61.5 km across, with the 62 to 450 km range completely empty, so the +# 200 km cap discards both junk modes with wide margin and touches no real file. +# Boxes are discarded rather than trimmed: the noise clusters AT the origin, so +# rejecting the outliers instead would leave a small, credible-looking box in +# the middle of Portugal. +_MAX_PLAUSIBLE_SPAN = {"metres": 200_000.0, "geographic": 5.0} +# A bbox built from only a handful of the sampled traces covers a sliver of the +# real line while claiming full coverage. When enough traces were sampled but +# almost none carried usable coordinates, the box is dropped and the coverage +# flagged partial: a family of 95-182 GB GEOM files publishes a ~130 m box from +# 2-3 of 100 sampled fixes. The parsed floor keeps genuinely small files (a few +# traces, a few fixes) from tripping it. +_MIN_SUPPORT_PARSED = 20 +_MIN_SUPPORT_FLOOR = 5 +_MIN_SUPPORT_FRACTION = 0.10 +# INT32 min/max are "unset" fill values, never coordinates (real archive traces +# carry INT32_MIN in the CDP fields). +_INT32_SENTINELS = (-(2**31), 2**31 - 1) + +# Data sample format code -> (label, bytes per sample). Anything else yields +# partial metadata because the trace length cannot be computed. +_SAMPLE_FORMATS = { + 1: ("ibm-float32", 4), + 2: ("int32", 4), + 3: ("int16", 2), + 5: ("ieee-float32", 4), + 6: ("ieee-float64", 8), + 8: ("int8", 1), + 10: ("uint32", 4), + 11: ("uint16", 2), + 16: ("uint8", 1), +} + +# Coordinate units codes: 2 and 3 are geographic (the only ones that yield a +# bbox_4326), 4 is packed DMS which is never converted (storing the packed +# integers would fake a bbox), and everything else is treated as projected +# metres in an UNDECLARED system - the CRS is never guessed. +_UNITS_ARC_SECONDS = 2 +_UNITS_DEGREES = 3 +_UNITS_DMS = 4 +_UNITS_LABELS = {0: "unset", 1: "metres", 2: "arc-seconds", 3: "degrees", 4: "dms"} + + +class _BinaryHeader(typing.NamedTuple): + """The fields used from the 400-byte binary file header (file offset 3200).""" + + sample_interval: int # offset 16; raw and unitless, see SegyMetadata + samples_per_trace: int # offset 20; authoritative (the per-trace copy lies) + sample_format: int # offset 24; data sample format code + revision: int # offset 300; major SEG-Y revision, normalised on parsing + extended_header_count: int # offset 304; int16, -1 means "variable" + extra_trace_headers: int # offset 306; rev2 additional 240-byte trace headers + + +# Pad bytes ("x") skip everything between the _BinaryHeader fields above. +_BINARY_HEADER_FORMAT = "16x H 2x H 2x H 274x H 2x h i" +_BIG_ENDIAN_BINARY_HEADER = struct.Struct(">" + _BINARY_HEADER_FORMAT) +_LITTLE_ENDIAN_BINARY_HEADER = struct.Struct("<" + _BINARY_HEADER_FORMAT) + + +class _TraceHeader(typing.NamedTuple): + """The fields used from each 240-byte trace header.""" + + scalco: int # offset 70; coordinate scalar: negative divides, positive multiplies + src_x: int # offset 72; source coordinates, the preferred position source + src_y: int # offset 76 + units: int # offset 88; coordinate units code, see _UNITS_LABELS + year: int # offset 156; trace time (2-digit years are ambiguous -> skipped) + day: int # offset 158; day of year + hour: int # offset 160 + minute: int # offset 162 + second: int # offset 164 + cdp_x: int # offset 180; CDP coordinates, the fallback position source + cdp_y: int # offset 184 + + +_TRACE_HEADER_FORMAT = "70x h i i 8x H 66x H H H H H 14x i i" +_BIG_ENDIAN_TRACE_HEADER = struct.Struct(">" + _TRACE_HEADER_FORMAT) +_LITTLE_ENDIAN_TRACE_HEADER = struct.Struct("<" + _TRACE_HEADER_FORMAT) + + +def extract_segy_metadata(path: Path | str) -> SegyMetadata: + """Extract metadata from a SEG-Y file by sampling ~100 trace headers. + + Only the two file headers plus a constant number of trace headers are read, + so a multi-GB file costs the same as a tiny one. Sampled coordinates give a + native bbox; bbox_4326 is only set when the traces declare geographic units, + because projected metres in an undeclared CRS must never be guessed at (the + units go into the description instead). A file too small to hold a single + trace, or whose binary header is not recognisable in either byte order, + raises ValueError. + """ + p = Path(path) + size = p.stat().st_size + if size < _MINIMUM_FILE_BYTES: + raise ValueError(f"{p.name} does not look like a SEG-Y file") + with p.open("rb") as fh: + fh.seek(_TEXTUAL_HEADER_BYTES) + headers = _parse_binary_header(fh.read(_BINARY_HEADER_BYTES)) + if headers is None: + raise ValueError(f"{p.name} does not look like a SEG-Y file") + binary_header, trace_header_struct = headers + format_label, bytes_per_sample = _SAMPLE_FORMATS.get( + binary_header.sample_format, + (f"format code {binary_header.sample_format}", None), + ) + header_facts = { + "driver": "SEG-Y", + "sample_interval": binary_header.sample_interval or None, + "samples_per_trace": binary_header.samples_per_trace or None, + "sample_format": format_label, + } + if ( + bytes_per_sample is None + or binary_header.samples_per_trace == 0 + # rev2 lets every trace carry N additional 240-byte headers; no real + # file uses them, so rather than assume a layout we keep the header + # facts - a wrong trace length would desync the whole sampling grid + or (binary_header.revision >= 2 and binary_header.extra_trace_headers) + ): + # without a trustworthy trace length the traces cannot be walked; + # keep the header facts instead of failing + return SegyMetadata(**header_facts) + data_start = _locate_trace_data(binary_header, size) + trace_length = ( + _TRACE_HEADER_BYTES + binary_header.samples_per_trace * bytes_per_sample + ) + # FLOOR division: a non-zero remainder is a trailing chunk (real files + # carry small trailers), never a reason to fail. + trace_count = (size - data_start) // trace_length + if trace_count == 0: + return SegyMetadata(**header_facts, trace_count=0) + + boxes: dict[str, tuple[float, float, float, float] | None] = { + "metres": None, + "geographic": None, + } + accepted = {"metres": 0, "geographic": 0} + units_counts: collections.Counter[int] = collections.Counter() + parsed_count = 0 + garbage_count = 0 + min_date = max_date = None + for index in _sample_indices(trace_count): + fh.seek(data_start + index * trace_length) + buffer = fh.read(_TRACE_HEADER_BYTES) + if len(buffer) < trace_header_struct.size: + break + trace = _TraceHeader(*trace_header_struct.unpack_from(buffer)) + parsed_count += 1 + units_counts[trace.units] += 1 + date = _parse_trace_date(trace) + if date is not None: + min_date = date if min_date is None else min(min_date, date) + max_date = date if max_date is None else max(max_date, date) + raw = _select_raw_coordinate(trace) + if raw is None or trace.units == _UNITS_DMS: + continue + point = _plausible_point(raw, trace.scalco, trace.units) + if point is None: + garbage_count += 1 + continue + category, x, y = point + accepted[category] += 1 + box = boxes[category] + boxes[category] = ( + (x, y, x, y) + if box is None + else (min(box[0], x), min(box[1], y), max(box[2], x), max(box[3], y)) + ) + + if parsed_count == 0: + return SegyMetadata(**header_facts, trace_count=trace_count) + for category, box in boxes.items(): + span = _MAX_PLAUSIBLE_SPAN[category] + if box is not None and (box[2] - box[0] > span or box[3] - box[1] > span): + # the points behind it were garbage that slipped the per-value + # filter, so counting them as such is what flags the coverage + boxes[category] = None + garbage_count += accepted[category] + accepted[category] = 0 + # A box supported by only a sliver of the sampled traces is dropped and the + # coverage flagged, but only once enough traces were sampled to trust the + # ratio - a genuinely small file with a few real fixes must keep its box. + accepted_total = accepted["metres"] + accepted["geographic"] + low_support = parsed_count >= _MIN_SUPPORT_PARSED and 0 < accepted_total < max( + _MIN_SUPPORT_FLOOR, parsed_count * _MIN_SUPPORT_FRACTION + ) + if low_support: + boxes["metres"] = boxes["geographic"] = None + accepted["metres"] = accepted["geographic"] = 0 + dominant_units = units_counts.most_common(1)[0][0] + shared_facts = { + **header_facts, + "trace_count": trace_count, + "coordinate_units": _UNITS_LABELS.get(dominant_units, f"code {dominant_units}"), + "coverage": ( + "partial" + if low_support or garbage_count > parsed_count * _GARBAGE_FRACTION_LIMIT + else "full" + ), + "temporal_extent_begin": min_date, + "temporal_extent_end": max_date, + } + if accepted["geographic"] > accepted["metres"]: + return SegyMetadata( + **shared_facts, + # geographic positions carry no datum declaration; WGS84-family in + # practice (sub-metre difference, far below the bbox buffer) - the + # same documented approximation as KMALL + epsg=4326, + crs_name="WGS 84", + bbox_native=boxes["geographic"], + bbox_4326=boxes["geographic"], + ) + return SegyMetadata(**shared_facts, bbox_native=boxes["metres"]) + + +def _parse_binary_header( + buffer: bytes, +) -> tuple[_BinaryHeader, struct.Struct] | None: + """Read the binary header, detecting endianness via the format code. + + Big-endian (every real file so far) is tried first; a byte order is accepted + when it yields a sample format code in 1..16. Returns the header together + with the matching trace-header struct, or None when neither order works or + the file ended mid-header (it can shrink between the stat and the read). + """ + if len(buffer) < _BIG_ENDIAN_BINARY_HEADER.size: + return None + for binary_struct, trace_struct in ( + (_BIG_ENDIAN_BINARY_HEADER, _BIG_ENDIAN_TRACE_HEADER), + (_LITTLE_ENDIAN_BINARY_HEADER, _LITTLE_ENDIAN_TRACE_HEADER), + ): + header = _BinaryHeader(*binary_struct.unpack_from(buffer)) + if 1 <= header.sample_format <= 16: + # The revision is a uint16 whose major number is the high byte + # (conformant rev1 writers store 0x0100 and rev2 splits it into + # major/minor bytes), except for writers that store a plain 1. + # Reading a single byte instead would miss it on a little-endian file. + code = header.revision + return ( + header._replace(revision=code >> 8 if code > 0xFF else code), + trace_struct, + ) + return None + + +def _locate_trace_data(binary_header: _BinaryHeader, size: int) -> int: + """Offset of the first trace: 3600 plus any trusted extended textual headers. + + The count at offset 304 is only meaningful from revision 1 on (the rev0 + bytes are undefined); -1 means "variable" and a count that is absurd or + leaves no room for a single trace header is ignored too. + """ + count = binary_header.extended_header_count + if ( + binary_header.revision >= 1 + and 0 < count <= _MAX_EXTENDED_HEADERS + and _TRACE_DATA_START + count * _TEXTUAL_HEADER_BYTES + _TRACE_HEADER_BYTES + <= size + ): + return _TRACE_DATA_START + count * _TEXTUAL_HEADER_BYTES + return _TRACE_DATA_START + + +def _sample_indices(trace_count: int) -> list[int]: + """Up to ~_SAMPLE_TARGET evenly spread indices, endpoints always included.""" + if trace_count <= _SAMPLE_TARGET: + return list(range(trace_count)) + last = trace_count - 1 + return sorted( + {step * last // (_SAMPLE_TARGET - 1) for step in range(_SAMPLE_TARGET)} + ) + + +def _select_raw_coordinate(trace: _TraceHeader) -> tuple[int, int] | None: + """Raw (x, y) from the source fields, falling back to CDP; None when absent. + + INT32 min/max fill and all-zero pairs mean "no coordinates", not values. + """ + for raw_x, raw_y in ((trace.src_x, trace.src_y), (trace.cdp_x, trace.cdp_y)): + if raw_x in _INT32_SENTINELS or raw_y in _INT32_SENTINELS: + continue + if raw_x == 0 and raw_y == 0: + continue + return raw_x, raw_y + return None + + +def _plausible_point( + raw: tuple[int, int], scalco: int, units: int +) -> tuple[str, float, float] | None: + """Scaled, plausibility-checked point as (category, x, y); None when garbage. + + Geographic units must land inside the lon/lat range; anything else is + treated as projected metres and rejected only on absurd magnitude. This + filter is what catches samples that landed mid-trace. + """ + x = _apply_coordinate_scalar(raw[0], scalco) + y = _apply_coordinate_scalar(raw[1], scalco) + if not (math.isfinite(x) and math.isfinite(y)): + return None + if units in (_UNITS_ARC_SECONDS, _UNITS_DEGREES): + if units == _UNITS_ARC_SECONDS: + x /= 3600.0 + y /= 3600.0 + if abs(x) > 180.0 or abs(y) > 90.0: + return None + return "geographic", x, y + if abs(x) >= _MAX_PLAUSIBLE_METRES or abs(y) >= _MAX_PLAUSIBLE_METRES: + return None + return "metres", x, y + + +def _apply_coordinate_scalar(value: int, scalco: int) -> float: + if scalco < 0: + return value / -scalco + if scalco > 0: + return float(value * scalco) + return float(value) + + +def _parse_trace_date(trace: _TraceHeader) -> dt.date | None: + """Best-effort date from one trace's time fields; None when implausible. + + A single corrupt time field must never abort the extraction (the same + guarantee as KMALL), so every field is range-checked per trace: a 4-digit + year (2-digit years are ambiguous), a day of year, and a sane time of day + as a corruption tell. + """ + if not 1970 <= trace.year <= 2100: + return None + if not 1 <= trace.day <= 366: + return None + if trace.hour > 23 or trace.minute > 59 or trace.second > 59: + return None + date = dt.date(trace.year, 1, 1) + dt.timedelta(days=trace.day - 1) + if date.year != trace.year: + # day 366 of a non-leap year rolls into January + return None + return date diff --git a/tests/conftest.py b/tests/conftest.py index 22315657..b509984c 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -15,6 +15,7 @@ ) from seis_lab_data.db.commands import ( datasetcategories as category_commands, + discovery as discovery_commands, projects as project_commands, surveymissions as mission_commands, surveyrelatedrecords as record_commands, @@ -119,13 +120,19 @@ async def bootstrap_workflow_stages(db, db_session_maker, bootstrap_data): @pytest_asyncio.fixture async def bootstrap_asset_discovery_configurations( - db, db_session_maker, bootstrap_data + db, + db_session_maker, + bootstrap_data, + bootstrap_dataset_categories, + bootstrap_workflow_stages, ): created = [] async with db_session_maker() as session: for to_create in bootstrap_data[constants.ResourceType.ASSET_DISCOVERY_CONFIG]: created.append( - await stage_commands.create_workflow_stage(session, to_create) + await discovery_commands.create_asset_discovery_configuration( + session, to_create + ) ) yield created diff --git a/tests/test_extractors.py b/tests/test_extractors.py new file mode 100644 index 00000000..21633045 --- /dev/null +++ b/tests/test_extractors.py @@ -0,0 +1,1166 @@ +"""Unit tests for the metadata extractors package. + +Pure functions, no DB: these run by default (no marker). GDAL is required; the +module self-skips where it is unavailable. Everything is exercised against +synthetic files built in-test, so the suite needs no sample-data archive. +""" + +import datetime as dt +import logging +import math +import struct + +import pytest + +pytest.importorskip("osgeo") + +from osgeo import gdal, ogr, osr # noqa: E402 + +from seis_lab_data.tasks.extractors import ( # noqa: E402 + dispatch, + schemas, +) +from seis_lab_data.tasks.extractors.gdal_raster import ( # noqa: E402 + extract_raster_metadata, +) +from seis_lab_data.tasks.extractors.gdal_vector import ( # noqa: E402 + extract_vector_metadata, +) +from seis_lab_data.tasks.extractors.kmall import ( # noqa: E402 + extract_kmall_metadata, +) +from seis_lab_data.tasks.extractors.segy import ( # noqa: E402 + extract_segy_metadata, +) + +gdal.UseExceptions() + +# EPSG:3763 (ETRS89 / Portugal TM06); coordinates near its false origin land in +# central Portugal, so projected bboxes fall in a checkable lon/lat window. +_PT_TM06 = 3763 +_POINTS = [(0.0, 0.0), (1000.0, 1000.0), (2000.0, 500.0)] + + +def _srs(epsg: int) -> osr.SpatialReference: + srs = osr.SpatialReference() + srs.ImportFromEPSG(epsg) + return srs + + +@pytest.fixture(scope="session") +def synthetic_geotiff(tmp_path_factory): + path = tmp_path_factory.mktemp("raster") / "grid.tif" + ds = gdal.GetDriverByName("GTiff").Create(str(path), 10, 10, 1, gdal.GDT_Float32) + ds.SetGeoTransform((0.0, 10.0, 0.0, 2000.0, 0.0, -10.0)) + ds.SetProjection(_srs(_PT_TM06).ExportToWkt()) + band = ds.GetRasterBand(1) + band.SetNoDataValue(-9999.0) + band.Fill(1.0) + ds = None + return path + + +@pytest.fixture(scope="session") +def synthetic_xyz(tmp_path_factory): + path = tmp_path_factory.mktemp("xyz") / "grid.xyz" + lines = [f"{x} {y} {x + y}" for y in range(3) for x in range(3)] + path.write_text("\n".join(lines) + "\n") + return path + + +def _write_shapefile(path, srs): + driver = ogr.GetDriverByName("ESRI Shapefile") + ds = driver.CreateDataSource(str(path)) + layer = ds.CreateLayer("points", srs=srs, geom_type=ogr.wkbPoint) + for x, y in _POINTS: + feature = ogr.Feature(layer.GetLayerDefn()) + geom = ogr.Geometry(ogr.wkbPoint) + geom.AddPoint(x, y) + feature.SetGeometry(geom) + layer.CreateFeature(feature) + feature = None + ds = None + + +@pytest.fixture(scope="session") +def synthetic_shapefile(tmp_path_factory): + path = tmp_path_factory.mktemp("vector") / "points.shp" + _write_shapefile(path, _srs(_PT_TM06)) + return path + + +@pytest.fixture(scope="session") +def synthetic_shapefile_no_prj(tmp_path_factory): + path = tmp_path_factory.mktemp("vector_noprj") / "points.shp" + _write_shapefile(path, None) + return path + + +@pytest.fixture(scope="session") +def broken_vector(tmp_path_factory): + # A valid shapefile with its .shx removed: reproduces the real prr_eolicas + # "missing .shx" condition, which GDAL refuses to open (SHAPE_RESTORE_SHX + # defaults to NO) -> RuntimeError. + path = tmp_path_factory.mktemp("broken") / "points.shp" + _write_shapefile(path, _srs(_PT_TM06)) + path.with_suffix(".shx").unlink() + return path + + +def test_raster_metadata_fields(synthetic_geotiff): + result = extract_raster_metadata(synthetic_geotiff) + + assert isinstance(result, schemas.RasterMetadata) + assert result.driver == "GTiff" + assert result.width == 10 + assert result.height == 10 + assert result.band_count == 1 + assert result.pixel_size_x == 10.0 + assert result.pixel_size_y == -10.0 + assert result.nodata == -9999.0 + assert result.epsg == _PT_TM06 + assert "TM06" in result.crs_name + + assert result.bbox_native is not None + minx, miny, maxx, maxy = result.bbox_native + assert minx < maxx and miny < maxy + + assert result.bbox_4326 is not None + lon_min, lat_min, lon_max, lat_max = result.bbox_4326 + assert -10.0 < lon_min <= lon_max < -6.0 + assert 36.0 < lat_min <= lat_max < 43.0 + + +def test_xyz_raster(synthetic_xyz): + result = extract_raster_metadata(synthetic_xyz) + + assert result.driver == "XYZ" + assert result.width == 3 + assert result.height == 3 + # No .prj alongside an XYZ grid: no CRS, native bbox only. + assert result.epsg is None + assert result.bbox_native is not None + assert result.bbox_4326 is None + + +def test_vector_metadata_fields(synthetic_shapefile): + result = extract_vector_metadata(synthetic_shapefile) + + assert isinstance(result, schemas.VectorMetadata) + assert result.driver == "ESRI Shapefile" + assert result.layer_count == 1 + assert result.feature_count == 3 + assert "Point" in result.geometry_type + assert result.epsg == _PT_TM06 + + assert result.bbox_native is not None + minx, miny, maxx, maxy = result.bbox_native + # Reorder trap: OGR GetExtent returns (minx, maxx, miny, maxy). If it were stored + # verbatim, miny would be 2000 and this assertion would fail. + assert minx < maxx + assert miny < maxy + + assert result.bbox_4326 is not None + lon_min, lat_min, lon_max, lat_max = result.bbox_4326 + assert -10.0 < lon_min <= lon_max < -6.0 + assert 36.0 < lat_min <= lat_max < 43.0 + + +def test_vector_without_prj(synthetic_shapefile_no_prj): + result = extract_vector_metadata(synthetic_shapefile_no_prj) + + assert result.epsg is None + assert result.crs_name is None + assert result.crs_wkt is None + assert result.bbox_4326 is None + assert result.bbox_native is not None + assert result.feature_count == 3 + + +def test_broken_vector_raises(broken_vector): + with pytest.raises(RuntimeError): + extract_vector_metadata(broken_vector) + + +def test_dispatch_routes_raster(synthetic_geotiff): + assert isinstance( + dispatch.dispatch_extractor(synthetic_geotiff), schemas.RasterMetadata + ) + + +def test_dispatch_routes_xyz(synthetic_xyz): + result = dispatch.dispatch_extractor(synthetic_xyz) + assert isinstance(result, schemas.RasterMetadata) + assert result.driver == "XYZ" + + +def test_dispatch_routes_vector(synthetic_shapefile): + assert isinstance( + dispatch.dispatch_extractor(synthetic_shapefile), schemas.VectorMetadata + ) + + +@pytest.mark.parametrize("suffix", [".shx", ".dbf", ".prj", ".cpg", ".qmd", ".sbn"]) +def test_dispatch_ignores_sidecars(tmp_path, suffix): + sidecar = tmp_path / f"points{suffix}" + sidecar.write_bytes(b"\x00") + assert dispatch.dispatch_extractor(sidecar) is None + + +def test_dispatch_unknown_extension(tmp_path): + unknown = tmp_path / "notes.txt" + unknown.write_text("hello") + assert dispatch.dispatch_extractor(unknown) is None + + +def test_dispatch_directory_named_like_raster(tmp_path): + # A real directory named "F3_2022.tif" exists in the archive; it must not be opened. + fake = tmp_path / "F3_2022.tif" + fake.mkdir() + assert dispatch.dispatch_extractor(fake) is None + + +def test_dispatch_big_file_logs(synthetic_geotiff, monkeypatch, caplog): + monkeypatch.setattr(dispatch, "_BIG_FILE_LOG_BYTES", 0) + with caplog.at_level(logging.INFO, logger=dispatch.logger.name): + result = dispatch.dispatch_extractor(synthetic_geotiff) + assert isinstance(result, schemas.RasterMetadata) + assert any("large raster" in record.message for record in caplog.records) + + +def test_raster_describe(synthetic_geotiff): + text = extract_raster_metadata(synthetic_geotiff).describe() + assert text.startswith("Auto-extracted: GTiff raster") + assert "TM06" in text + assert "EPSG:3763" in text + assert len(text) <= 500 + + +def test_vector_describe(synthetic_shapefile): + text = extract_vector_metadata(synthetic_shapefile).describe() + assert text.startswith("Auto-extracted: ESRI Shapefile vector") + assert "Point" in text + assert len(text) <= 500 + + +def test_describe_truncates_to_500(): + absurd = schemas.RasterMetadata( + driver="GTiff", + width=1, + height=1, + band_count=1, + crs_name="X" * 2000, + ) + assert len(absurd.describe()) == 500 + + +def test_describe_crs_unknown(): + text = schemas.VectorMetadata( + driver="ESRI Shapefile", + layer_count=1, + feature_count=0, + geometry_type="Point", + ).describe() + assert "CRS: unknown." in text + + +# 2024-09-28 18:48:31 UTC, matching the archive's EM712 line 0419 +_KMALL_T0 = 1_727_549_311 +_KMALL_HEADER = struct.Struct(" walk stops + + result = extract_kmall_metadata(path) + + assert result.datagram_count == 3 + assert result.position_count == 1 + + +def test_kmall_truncated_tail_beyond_eof(tmp_path): + # the realistic acquisition-crash shape: a well-formed header whose + # num_bytes claims data beyond the end of the file + path = tmp_path / "line.kmall" + _write_kmall(path, [(40.5, -9.3)]) + with path.open("ab") as fh: + fh.write(_KMALL_HEADER.pack(99999, b"#MRZ", 0, 0, 712, _KMALL_T0 + 2, 0)) + + result = extract_kmall_metadata(path) + + assert result.datagram_count == 3 + assert result.position_count == 1 + + +def test_kmall_garbage_raises(tmp_path): + path = tmp_path / "junk.kmall" + path.write_bytes(b"this is definitely not a kmall file") + with pytest.raises(ValueError): + extract_kmall_metadata(path) + + +@pytest.mark.parametrize("content", [b"", b"tiny"]) +def test_kmall_too_small_raises(tmp_path, content): + # empty or sub-header files never enter the walk and must still raise + path = tmp_path / "empty.kmall" + path.write_bytes(content) + with pytest.raises(ValueError): + extract_kmall_metadata(path) + + +def test_dispatch_routes_kmall(tmp_path): + path = tmp_path / "line.kmall" + _write_kmall(path, [(40.5, -9.3)]) + assert isinstance(dispatch.dispatch_extractor(path), schemas.KmallMetadata) + + +def test_dispatch_big_kmall_logs(tmp_path, monkeypatch, caplog): + path = tmp_path / "line.kmall" + _write_kmall(path, [(40.5, -9.3)]) + monkeypatch.setattr(dispatch, "_BIG_FILE_LOG_BYTES", 0) + with caplog.at_level(logging.INFO, logger=dispatch.logger.name): + result = dispatch.dispatch_extractor(path) + assert isinstance(result, schemas.KmallMetadata) + assert any("large KMALL file" in record.message for record in caplog.records) + + +def test_kmall_describe(tmp_path): + path = tmp_path / "line.kmall" + _write_kmall(path, [(40.5, -9.3)]) + text = extract_kmall_metadata(path).describe() + assert text.startswith("Auto-extracted: Kongsberg KMALL, EM 712") + assert "position fix(es)" in text + assert "EPSG:4326" in text + assert len(text) <= 500 + + +# SEG-Y builders: tiny files with the real layout (3200-byte textual header + +# 400-byte binary header + fixed-length traces), big-endian by default. Day 272 +# of 2024 is 2024-09-28, matching the archive's survey dates. +_INT32_MIN = -(2**31) +_INT32_MAX = 2**31 - 1 + + +def _segy_binary_header(ns, sample_format, interval, revision, n_ext, order): + header = bytearray(400) + struct.pack_into(order + "H", header, 16, interval) + struct.pack_into(order + "H", header, 20, ns) + struct.pack_into(order + "H", header, 24, sample_format) + header[300] = revision + struct.pack_into(order + "h", header, 304, n_ext) + return bytes(header) + + +def _segy_trace( + ns=4, + bytes_per_sample=4, + src=(0, 0), + cdp=(0, 0), + scalco=0, + units=1, + year=2024, + day=272, + hour=12, + minute=30, + second=15, + order=">", +): + header = bytearray(240) + struct.pack_into(order + "h", header, 70, scalco) + struct.pack_into(order + "ii", header, 72, *src) + struct.pack_into(order + "H", header, 88, units) + struct.pack_into(order + "HHHHH", header, 156, year, day, hour, minute, second) + struct.pack_into(order + "ii", header, 180, *cdp) + return bytes(header) + b"\x00" * (ns * bytes_per_sample) + + +def _write_segy( + path, + traces, + ns=4, + sample_format=5, + interval=250, + revision=1, + n_ext=0, + order=">", + extended=b"", + trailer=b"", +): + path.write_bytes( + b" " * 3200 + + _segy_binary_header(ns, sample_format, interval, revision, n_ext, order) + + extended + + b"".join(traces) + + trailer + ) + + +def test_segy_metres_path(tmp_path): + path = tmp_path / "line.sgy" + _write_segy( + path, + [ + _segy_trace(src=(-50000, 150000), day=272), + _segy_trace(src=(-48000, 152000), day=273), + _segy_trace(src=(-49000, 151000), day=272), + ], + ) + + result = extract_segy_metadata(path) + + assert isinstance(result, schemas.SegyMetadata) + assert result.driver == "SEG-Y" + assert result.trace_count == 3 + assert result.samples_per_trace == 4 + assert result.sample_interval == 250 + assert result.sample_format == "ieee-float32" + assert result.coordinate_units == "metres" + assert result.coverage == "full" + # undeclared projected metres: native bbox only, the CRS is never guessed + assert result.bbox_native == (-50000.0, 150000.0, -48000.0, 152000.0) + assert result.bbox_4326 is None + assert result.epsg is None + assert result.temporal_extent_begin == dt.date(2024, 9, 28) + assert result.temporal_extent_end == dt.date(2024, 9, 29) + + +def test_segy_geographic_degrees(tmp_path): + # units=3 with scalco=-10000: raw ints are degrees times 10000 + path = tmp_path / "line.sgy" + _write_segy( + path, + [ + _segy_trace(src=(-87500, 422300), scalco=-10000, units=3), + _segy_trace(src=(-87000, 422800), scalco=-10000, units=3), + ], + ) + + result = extract_segy_metadata(path) + + assert result.coordinate_units == "degrees" + assert result.epsg == 4326 + assert result.crs_name == "WGS 84" + assert result.bbox_4326 == (-8.75, 42.23, -8.7, 42.28) + assert result.bbox_native == result.bbox_4326 + + +def test_segy_geographic_arc_seconds(tmp_path): + # the real CW HAT.sgy shape: units=2, stationary, arc-seconds with a + # negative scalar; -8.7385 deg is -31458.6 arc-seconds + path = tmp_path / "line.sgy" + _write_segy( + path, + [_segy_trace(src=(-314586, 1520280), scalco=-10, units=2)], + ) + + result = extract_segy_metadata(path) + + assert result.coordinate_units == "arc-seconds" + assert result.epsg == 4326 + lon_min, lat_min, lon_max, lat_max = result.bbox_4326 + assert lon_min == lon_max == pytest.approx(-8.7385) + assert lat_min == lat_max == pytest.approx(42.23) + + +def test_segy_dms_units_skipped(tmp_path): + # packed DDDMMSS.ss is never converted in v1; storing the packed integers + # verbatim would fake a bbox + path = tmp_path / "line.sgy" + _write_segy(path, [_segy_trace(src=(83030450, 421512), units=4)]) + + result = extract_segy_metadata(path) + + assert result.coordinate_units == "dms" + assert result.bbox_native is None + assert result.bbox_4326 is None + assert result.coverage == "full" + + +def test_segy_scalco_multiplies(tmp_path): + path = tmp_path / "line.sgy" + _write_segy(path, [_segy_trace(src=(-5000, 15000), scalco=10)]) + + result = extract_segy_metadata(path) + + assert result.bbox_native == (-50000.0, 150000.0, -50000.0, 150000.0) + + +def test_segy_cdp_sentinel_fill_ignored(tmp_path): + # real archive traces carry INT32_MIN fill in the CDP fields while src holds + # the good coordinates. The scalar matters: scaled by -10000 a sentinel + # becomes -214748.36, well inside the plausible-metres range, so only the + # sentinel guard itself keeps it out of the bbox + path = tmp_path / "line.sgy" + _write_segy( + path, + [ + _segy_trace(src=(-500000000, 1500000000), scalco=-10000), + _segy_trace(src=(0, 0), cdp=(_INT32_MIN, _INT32_MIN), scalco=-10000), + _segy_trace(src=(_INT32_MAX, _INT32_MAX), scalco=-10000), + ], + ) + + result = extract_segy_metadata(path) + + assert result.bbox_native == (-50000.0, 150000.0, -50000.0, 150000.0) + assert result.coverage == "full" + + +def test_segy_src_takes_precedence_over_cdp(tmp_path): + # src is the primary source; a divergent cdp feed must be IGNORED, never + # merged into the bbox (the KMALL SPO/CPO rule, applied to SEG-Y) + path = tmp_path / "line.sgy" + _write_segy( + path, + [ + _segy_trace(src=(-50000, 150000), cdp=(-90000, 190000)), + _segy_trace(src=(-48000, 152000), cdp=(-91000, 191000)), + ], + ) + + result = extract_segy_metadata(path) + + assert result.bbox_native == (-50000.0, 150000.0, -48000.0, 152000.0) + + +def test_segy_cdp_fallback(tmp_path): + # chirp/Innomar/sbp populate only src, UHRS both; a src-less trace must + # still yield its CDP coordinates + path = tmp_path / "line.sgy" + _write_segy(path, [_segy_trace(src=(0, 0), cdp=(-50000, 150000))]) + + result = extract_segy_metadata(path) + + assert result.bbox_native == (-50000.0, 150000.0, -50000.0, 150000.0) + + +def test_segy_without_coordinates(tmp_path): + path = tmp_path / "line.sgy" + _write_segy(path, [_segy_trace(), _segy_trace()]) + + result = extract_segy_metadata(path) + + assert result.bbox_native is None + assert result.bbox_4326 is None + assert result.coverage == "full" + assert result.temporal_extent_begin == dt.date(2024, 9, 28) + + +def test_segy_mostly_garbage_coordinates_partial(tmp_path): + # implausible coordinate magnitudes are the tell of samples landing + # mid-trace; above half the samples the coverage is flagged partial, but + # the plausible traces still contribute + path = tmp_path / "line.sgy" + garbage = [_segy_trace(src=(2_000_000_000, 7)) for _ in range(8)] + good = [ + _segy_trace(src=(-50000, 150000)), + _segy_trace(src=(-48000, 152000)), + ] + _write_segy(path, garbage + good) + + result = extract_segy_metadata(path) + + assert result.coverage == "partial" + assert result.bbox_native == (-50000.0, 150000.0, -48000.0, 152000.0) + + +def test_segy_two_digit_year_temporal_none(tmp_path): + path = tmp_path / "line.sgy" + _write_segy(path, [_segy_trace(src=(-50000, 150000), year=24)]) + + result = extract_segy_metadata(path) + + assert result.temporal_extent_begin is None + assert result.temporal_extent_end is None + assert result.bbox_native is not None + + +def test_segy_corrupt_time_fields_skipped(tmp_path): + # a single corrupt time field skips only that trace's date and never aborts + # the extraction (the KMALL guarantee) + path = tmp_path / "line.sgy" + _write_segy( + path, + [ + _segy_trace(src=(-50000, 150000), hour=99), + _segy_trace(src=(-48000, 152000), year=2023, day=366), # not a leap year + _segy_trace(src=(-49000, 151000), day=272), + ], + ) + + result = extract_segy_metadata(path) + + assert result.temporal_extent_begin == dt.date(2024, 9, 28) + assert result.temporal_extent_end == dt.date(2024, 9, 28) + assert result.bbox_native == (-50000.0, 150000.0, -48000.0, 152000.0) + + +def test_segy_trailing_bytes_tolerated(tmp_path): + # a real BSTK file carries a 5,056-byte trailer: floor division ignores it + path = tmp_path / "line.sgy" + _write_segy( + path, + [_segy_trace(src=(-50000, 150000)), _segy_trace(src=(-48000, 152000))], + trailer=b"\x87" * 100, + ) + + result = extract_segy_metadata(path) + + assert result.trace_count == 2 + assert result.coverage == "full" + assert result.bbox_native == (-50000.0, 150000.0, -48000.0, 152000.0) + + +def test_segy_zero_samples_partial_metadata(tmp_path): + # ns=0 in the binary header leaves no trustworthy trace length; keep the + # header facts, never fall back to the unreliable per-trace ns field + path = tmp_path / "line.sgy" + _write_segy(path, [b"\x00" * 240], ns=0) + + result = extract_segy_metadata(path) + + assert result.samples_per_trace is None + assert result.trace_count is None + assert result.sample_format == "ieee-float32" + assert result.sample_interval == 250 + assert result.bbox_native is None + + +def test_segy_unknown_format_partial_metadata(tmp_path): + path = tmp_path / "line.sgy" + _write_segy(path, [_segy_trace()], sample_format=9) + + result = extract_segy_metadata(path) + + assert result.sample_format == "format code 9" + assert result.trace_count is None + assert result.bbox_native is None + + +def test_segy_little_endian(tmp_path): + path = tmp_path / "line.sgy" + _write_segy( + path, + [_segy_trace(src=(-50000, 150000), order="<")], + order="<", + ) + + result = extract_segy_metadata(path) + + assert result.sample_format == "ieee-float32" + assert result.samples_per_trace == 4 + assert result.bbox_native == (-50000.0, 150000.0, -50000.0, 150000.0) + + +def test_segy_rev0_extended_count_ignored(tmp_path): + # the rev0 bytes at offset 304 are undefined noise; trusting them would + # push data_start past the traces + path = tmp_path / "line.sgy" + _write_segy( + path, + [_segy_trace(src=(-50000, 150000)), _segy_trace(src=(-48000, 152000))], + revision=0, + n_ext=50, + ) + + result = extract_segy_metadata(path) + + assert result.trace_count == 2 + assert result.bbox_native == (-50000.0, 150000.0, -48000.0, 152000.0) + + +def test_segy_variable_extended_count_untrusted(tmp_path): + # n_ext == -1 means "variable, stanza-terminated" and is not trusted + path = tmp_path / "line.sgy" + _write_segy(path, [_segy_trace(src=(-50000, 150000))], revision=1, n_ext=-1) + + result = extract_segy_metadata(path) + + assert result.trace_count == 1 + assert result.bbox_native == (-50000.0, 150000.0, -50000.0, 150000.0) + + +def test_segy_extended_textual_headers(tmp_path): + path = tmp_path / "line.sgy" + _write_segy( + path, + [_segy_trace(src=(-50000, 150000)), _segy_trace(src=(-48000, 152000))], + revision=1, + n_ext=1, + extended=b" " * 3200, + ) + + result = extract_segy_metadata(path) + + assert result.trace_count == 2 + assert result.bbox_native == (-50000.0, 150000.0, -48000.0, 152000.0) + + +def test_segy_sampling_includes_endpoints(tmp_path): + # 150 single-sample traces: only ~100 are sampled, but the first and last + # trace (the line endpoints) must always contribute to the bbox + path = tmp_path / "line.sgy" + first = _segy_trace(ns=1, bytes_per_sample=1, src=(-52000, 148000)) + middle = [ + _segy_trace(ns=1, bytes_per_sample=1, src=(-49000, 151000)) for _ in range(148) + ] + last = _segy_trace(ns=1, bytes_per_sample=1, src=(-45000, 155000)) + _write_segy(path, [first, *middle, last], ns=1, sample_format=8) + + result = extract_segy_metadata(path) + + assert result.trace_count == 150 + assert result.sample_format == "int8" + assert result.bbox_native == (-52000.0, 148000.0, -45000.0, 155000.0) + + +@pytest.mark.parametrize( + "content", + [b"", b"tiny", b"\x00" * 4000, b"\xff" * 4000], +) +def test_segy_not_a_segy_raises(tmp_path, content): + # too small to hold a trace, or a binary header whose format code is + # invalid in both byte orders + path = tmp_path / "junk.sgy" + path.write_bytes(content) + with pytest.raises(ValueError): + extract_segy_metadata(path) + + +# production holds 925 .SEGY files besides the 68k lowercase ones, so the +# uppercase variants are pinned here too +@pytest.mark.parametrize("suffix", [".sgy", ".segy", ".SGY", ".SEGY"]) +def test_dispatch_routes_segy(tmp_path, suffix): + path = tmp_path / f"line{suffix}" + _write_segy(path, [_segy_trace(src=(-50000, 150000))]) + assert isinstance(dispatch.dispatch_extractor(path), schemas.SegyMetadata) + + +def test_segy_implausible_span_discards_bbox(tmp_path): + # real owf-2025 files exist whose src AND cdp fields hold noise around the + # projection origin with sporadic spikes into the millions; the resulting + # box spanned thousands of km while coverage still said "full". Rejecting + # the outliers instead would leave a small credible box AT the origin, + # which is on land in central Portugal, so the whole box goes. + path = tmp_path / "line.sgy" + _write_segy( + path, + [ + _segy_trace(src=(140, 120), scalco=-100), + _segy_trace(src=(-240, -330), scalco=-100), + _segy_trace(src=(588515700, 600505400), scalco=-100), + _segy_trace(src=(-100, -80), scalco=-100), + ], + ) + + result = extract_segy_metadata(path) + + assert result.bbox_native is None + assert result.bbox_4326 is None + assert result.coverage == "partial" + # the rest of the metadata is unaffected and still worth having + assert result.trace_count == 4 + assert result.sample_format == "ieee-float32" + assert result.temporal_extent_begin == dt.date(2024, 9, 28) + + +def test_segy_implausible_span_without_any_rejected_value(tmp_path): + # the case a garbage-count rule cannot catch: every value stays inside the + # per-value plausibility filter, so only the span betrays the noise + path = tmp_path / "line.sgy" + _write_segy( + path, + [ + _segy_trace(src=(140, 120), scalco=-100), + _segy_trace(src=(-986000000, -493000000), scalco=-100), + _segy_trace(src=(995000000, 497000000), scalco=-100), + ], + ) + + result = extract_segy_metadata(path) + + assert result.bbox_native is None + assert result.coverage == "partial" + + +def test_segy_geographic_implausible_span_discards_bbox(tmp_path): + # the same protection on the only path that draws a map rectangle + path = tmp_path / "line.sgy" + _write_segy( + path, + [ + _segy_trace(src=(-87500, 422300), scalco=-10000, units=3), + _segy_trace(src=(1200000, -600000), scalco=-10000, units=3), + ], + ) + + result = extract_segy_metadata(path) + + assert result.bbox_4326 is None + assert result.bbox_native is None + assert result.epsg is None # no CRS claim without a box to attach it to + assert result.coverage == "partial" + + +def test_segy_survey_wide_span_is_kept(tmp_path): + # a legitimately long line must survive: real files measure up to 61.5 km + # and the cap is 200 km, so a 60 km line is still stored + path = tmp_path / "line.sgy" + _write_segy( + path, + [ + _segy_trace(src=(-13000000, 2200000), scalco=-100), + _segy_trace(src=(-13000000, 8200000), scalco=-100), + ], + ) + + result = extract_segy_metadata(path) + + assert result.bbox_native == (-130000.0, 22000.0, -130000.0, 82000.0) + assert result.coverage == "full" + + +def test_segy_burst_corruption_span_discards_bbox(tmp_path): + # the second junk mode the 66k scan found (A7808 ET_SBP): a real ~2.5 km + # line with a burst of corrupt-nav traces stretched the box to ~494 km with + # garbage_count 0, which the old 500 km cap let through as coverage=full. + # The 200 km cap catches it while leaving every real line (<= 61.5 km) alone. + path = tmp_path / "line.sgy" + _write_segy( + path, + [ + _segy_trace(src=(-12198249, 5017758), scalco=-100), # (-121982, 50177) + _segy_trace(src=(-738249, -44428591), scalco=-100), # (-7382, -444285) + ], + ) + + result = extract_segy_metadata(path) + + assert result.bbox_native is None + assert result.coverage == "partial" + # the header facts are unaffected and still worth keeping + assert result.trace_count == 2 + + +def test_segy_low_support_bbox_dropped(tmp_path): + # the third junk mode the scan found (A7808 GEOM): a 95-182 GB file where + # only 2-3 of 100 sampled traces carry a fix and the rest are coordinate-less. + # A box from 3 points covers a sliver of the real line while claiming full + # coverage, so it is dropped and the coverage flagged. + path = tmp_path / "line.sgy" + traces = [_segy_trace(src=(-50000, 150000))] * 3 + [_segy_trace()] * 97 + _write_segy(path, traces) + + result = extract_segy_metadata(path) + + assert result.trace_count == 100 + assert result.bbox_native is None + assert result.coverage == "partial" + + +def test_segy_low_support_kept_for_small_file(tmp_path): + # the guard: a genuinely small file with only a few traces keeps its box even + # when only some carry a fix - too few samples to trust the support ratio + path = tmp_path / "line.sgy" + traces = [ + _segy_trace(src=(-50000, 150000)), + _segy_trace(src=(-48000, 152000)), + ] + [_segy_trace()] * 6 + _write_segy(path, traces) + + result = extract_segy_metadata(path) + + assert result.trace_count == 8 + assert result.bbox_native == (-50000.0, 150000.0, -48000.0, 152000.0) + assert result.coverage == "full" + + +def test_segy_units_unset_treated_as_metres(tmp_path): + # the real raw-file case: the UHRS rev1.segy files leave units at 0 + path = tmp_path / "line.sgy" + _write_segy(path, [_segy_trace(src=(-50000, 150000), units=0)]) + + result = extract_segy_metadata(path) + + assert result.coordinate_units == "unset" + assert result.bbox_native == (-50000.0, 150000.0, -50000.0, 150000.0) + assert result.bbox_4326 is None + assert result.epsg is None + + +def test_segy_geographic_out_of_range_rejected(tmp_path): + # the plausibility filter guards the only path that draws a map rectangle: + # a trace claiming degrees but carrying projected metres must not become a + # bbox somewhere off the planet + path = tmp_path / "line.sgy" + _write_segy( + path, + [ + _segy_trace(src=(-87500, 422300), scalco=-10000, units=3), + _segy_trace(src=(-500000, 1500000), units=3), # metres in a degrees file + ], + ) + + result = extract_segy_metadata(path) + + assert result.bbox_4326 == (-8.75, 42.23, -8.75, 42.23) + assert result.coverage == "full" # one garbage trace out of two is not "> half" + + +def test_segy_coverage_boundary_at_exactly_half(tmp_path): + # exactly half the samples garbage stays "full"; one more tips it to partial + path = tmp_path / "line.sgy" + good = [_segy_trace(src=(-50000, 150000)), _segy_trace(src=(-48000, 152000))] + garbage = [_segy_trace(src=(2_000_000_000, 7)) for _ in range(2)] + _write_segy(path, good + garbage) + assert extract_segy_metadata(path).coverage == "full" + + _write_segy(path, good + garbage + [_segy_trace(src=(2_000_000_000, 7))]) + assert extract_segy_metadata(path).coverage == "partial" + + +def test_segy_extended_headers_not_fitting_are_ignored(tmp_path): + # a declared extended-header count that would leave no room for a trace is + # garbage; trusting it would push data_start past the end of the file and + # yield a negative trace_count + path = tmp_path / "line.sgy" + _write_segy( + path, + [_segy_trace(src=(-50000, 150000))], + revision=1, + n_ext=99, # 99 * 3200 bytes that are not in the file + ) + + result = extract_segy_metadata(path) + + assert result.trace_count == 1 + assert result.bbox_native == (-50000.0, 150000.0, -50000.0, 150000.0) + + +def test_segy_little_endian_revision_word(tmp_path): + # a little-endian writer storing the rev1-style uint16 0x0100 puts the major + # number in the SECOND byte; reading a single byte would see revision 0 and + # then ignore the extended header, shifting the whole trace grid + path = tmp_path / "line.sgy" + header = bytearray(_segy_binary_header(4, 5, 250, 1, 1, "<")) + header[300:302] = b"\x00\x01" + path.write_bytes( + b" " * 3200 + + bytes(header) + + b" " * 3200 + + _segy_trace(src=(-50000, 150000), order="<") + + _segy_trace(src=(-48000, 152000), order="<") + ) + + result = extract_segy_metadata(path) + + assert result.trace_count == 2 + assert result.bbox_native == (-50000.0, 150000.0, -48000.0, 152000.0) + + +def test_segy_rev2_extra_trace_headers_partial_metadata(tmp_path): + # rev2 lets each trace carry additional 240-byte headers; the trace length + # is then unknown to us, and a desynced grid would store a confidently + # wrong bbox, so only the header facts survive + path = tmp_path / "line.sgy" + header = bytearray(_segy_binary_header(4, 5, 250, 2, 0, ">")) + struct.pack_into(">i", header, 306, 1) + path.write_bytes( + b" " * 3200 + + bytes(header) + + b"".join( + _segy_trace(src=src) + b"\x00" * 240 + for src in ((-50000, 150000), (-48000, 152000)) + ) + ) + + result = extract_segy_metadata(path) + + assert result.sample_format == "ieee-float32" + assert result.trace_count is None + assert result.bbox_native is None + + +def test_segy_truncated_binary_header_raises(tmp_path): + # a file that ends inside the binary header must still raise ValueError + # rather than let a struct error escape + path = tmp_path / "line.sgy" + _write_segy(path, [_segy_trace(src=(-50000, 150000))]) + truncated = tmp_path / "truncated.sgy" + truncated.write_bytes(path.read_bytes()[:3400]) + + with pytest.raises(ValueError): + extract_segy_metadata(truncated) + + +def test_segy_trace_count_zero(tmp_path): + # a file holding the headers but less than one full trace + path = tmp_path / "line.sgy" + path.write_bytes( + b" " * 3200 + _segy_binary_header(1000, 5, 250, 1, 0, ">") + b"\x00" * 300 + ) + + result = extract_segy_metadata(path) + + assert result.trace_count == 0 + assert result.samples_per_trace == 1000 + assert result.bbox_native is None + assert result.temporal_extent_begin is None + + +def test_segy_describe_partial_metadata(tmp_path): + # a file whose traces cannot be walked still gets a useful description + path = tmp_path / "line.sgy" + _write_segy(path, [_segy_trace()], sample_format=9) + text = extract_segy_metadata(path).describe() + assert text.startswith("Auto-extracted: SEG-Y") + assert "format code 9" in text + assert "trace(s)" not in text + assert "CRS: unknown." in text + + +def test_segy_describe_partial_coverage(tmp_path): + path = tmp_path / "line.sgy" + garbage = [_segy_trace(src=(2_000_000_000, 7)) for _ in range(8)] + _write_segy(path, garbage + [_segy_trace(src=(-50000, 150000))]) + text = extract_segy_metadata(path).describe() + assert "partly unreliable" in text + + +def test_segy_describe(tmp_path): + path = tmp_path / "line.sgy" + _write_segy(path, [_segy_trace(src=(-50000, 150000))]) + text = extract_segy_metadata(path).describe() + assert text.startswith("Auto-extracted: SEG-Y") + assert "1 trace(s)" in text + assert "coordinates in metres" in text + assert "CRS: unknown." in text + assert "Native bbox:" in text + assert len(text) <= 500 + + +def test_segy_describe_geographic(tmp_path): + path = tmp_path / "line.sgy" + _write_segy(path, [_segy_trace(src=(-87500, 422300), scalco=-10000, units=3)]) + text = extract_segy_metadata(path).describe() + assert "coordinates in degrees" in text + assert "EPSG:4326" in text diff --git a/tests/test_operations_discovery.py b/tests/test_operations_discovery.py new file mode 100644 index 00000000..dc0483fd --- /dev/null +++ b/tests/test_operations_discovery.py @@ -0,0 +1,446 @@ +"""Integration tests for mission discovery with metadata extraction. + +These run against the test database and a fake archive built under tmp_path. +GDAL is required to build the synthetic fixture files; the module self-skips +where it is unavailable. +""" + +import datetime as dt +import logging +import math +import uuid + +import pytest +import pytest_asyncio +import shapely +import sqlmodel +from geoalchemy2.shape import to_shape + +pytest.importorskip("osgeo") + +from osgeo import gdal, osr # noqa: E402 + +from seis_lab_data import constants # noqa: E402 +from seis_lab_data.db import models # noqa: E402 +from seis_lab_data.db.commands import ( # noqa: E402 + discovery as discovery_commands, + projects as project_commands, + surveymissions as mission_commands, +) +from seis_lab_data.operations import discovery as discovery_ops # noqa: E402 +from seis_lab_data.schemas import ( # noqa: E402 + common as common_schemas, + discovery as discovery_schemas, + events as event_schemas, + identifiers, + projects as project_schemas, + surveymissions as mission_schemas, +) +from seis_lab_data.tasks.extractors import schemas as extractor_schemas # noqa: E402 + +gdal.UseExceptions() + +_MISSION_RELATIVE_PATH = "surveys/test-mission" + + +class _EventCollector: + def __init__(self): + self.events: list[event_schemas.SeisLabDataEvent] = [] + + async def __call__(self, event: event_schemas.SeisLabDataEvent) -> None: + self.events.append(event) + + +def _write_geotiff(path): + path.parent.mkdir(parents=True, exist_ok=True) + srs = osr.SpatialReference() + srs.ImportFromEPSG(3763) + ds = gdal.GetDriverByName("GTiff").Create(str(path), 10, 10, 1, gdal.GDT_Float32) + ds.SetGeoTransform((0.0, 10.0, 0.0, 2000.0, 0.0, -10.0)) + ds.SetProjection(srs.ExportToWkt()) + ds.GetRasterBand(1).SetNoDataValue(-9999.0) + ds.GetRasterBand(1).Fill(1.0) + ds = None + + +async def _create_mission( + session, admin_user, project_id, mission_relative_path, name_suffix="" +): + mission = await mission_commands.create_survey_mission( + session, + mission_schemas.SurveyMissionCreate( + id=identifiers.SurveyMissionId(uuid.uuid4()), + owner_id=identifiers.UserId(admin_user.id), + project_id=project_id, + name=common_schemas.LocalizableDraftName( + en=f"Discovery test mission{name_suffix}" + ), + description=common_schemas.LocalizableDraftDescription(en=""), + relative_path=mission_relative_path, + ), + ) + return mission + + +@pytest_asyncio.fixture +async def discovery_env( + tmp_path, + settings, + db, + db_session_maker, + admin_user, + bootstrap_dataset_categories, + bootstrap_workflow_stages, +): + """A fake archive plus one project, mission and discovery configuration.""" + settings.readonly_archive_root_directory = tmp_path + async with db_session_maker() as session: + project = await project_commands.create_project( + session, + project_schemas.ProjectCreate( + id=identifiers.ProjectId(uuid.uuid4()), + owner_id=identifiers.UserId(admin_user.id), + name=common_schemas.LocalizableDraftName(en="Discovery test project"), + root_path="", + ), + ) + mission = await _create_mission( + session, + admin_user, + identifiers.ProjectId(project.id), + _MISSION_RELATIVE_PATH, + ) + configuration = await discovery_commands.create_asset_discovery_configuration( + session, + discovery_schemas.AssetDiscoveryConfigurationCreate( + id=identifiers.AssetDiscoveryConfId(uuid.uuid4()), + name="test tif", + relative_path_regexp=r"s01/.*\.tif", + workflow_stage_id=identifiers.WorkflowStageId( + bootstrap_workflow_stages[0].id + ), + dataset_category_id=identifiers.DatasetCategoryId( + bootstrap_dataset_categories[0].id + ), + ), + ) + return { + "archive_root": tmp_path, + "settings": settings, + "project": project, + "mission": mission, + "configuration": configuration, + } + + +async def _run_discovery(db_session_maker, mission_id, settings, admin_user): + collector = _EventCollector() + async with db_session_maker() as session: + await discovery_ops.run_mission_discovery( + request_id=identifiers.RequestId(uuid.uuid4()), + mission_id=identifiers.SurveyMissionId(mission_id), + session=session, + event_dispatcher=collector, + settings=settings, + user=admin_user, + ) + return collector + + +async def _get_mission_records(db_session_maker, mission_id): + async with db_session_maker() as session: + statement = sqlmodel.select(models.SurveyRelatedRecord).where( + models.SurveyRelatedRecord.survey_mission_id == mission_id + ) + return (await session.exec(statement)).all() + + +def _ended_events(collector): + return [ + e + for e in collector.events + if isinstance(e, event_schemas.DiscoveryEvent) + and e.modification == constants.DiscoveryStage.ENDED + ] + + +def test_bbox_tuple_to_wkt_buffers_every_side(): + wkt = discovery_ops._bbox_4326_tuple_to_wkt((10.0, 40.0, 10.5, 40.5)) + assert shapely.from_wkt(wkt).bounds == pytest.approx( + (9.9999, 39.9999, 10.5001, 40.5001) + ) + + +def test_bbox_tuple_to_wkt_pads_degenerate_point(): + wkt = discovery_ops._bbox_4326_tuple_to_wkt((-8.5, 39.0, -8.5, 39.0)) + polygon = shapely.from_wkt(wkt) + assert polygon.is_valid + assert polygon.area > 0 + assert polygon.bounds == pytest.approx((-8.5001, 38.9999, -8.4999, 39.0001)) + + +def test_bbox_tuple_to_wkt_discards_out_of_range(caplog): + with caplog.at_level(logging.WARNING, logger=discovery_ops.logger.name): + result = discovery_ops._bbox_4326_tuple_to_wkt((200.0, 40.0, 201.0, 41.0)) + assert result is None + assert any("Discarding extracted bbox" in r.message for r in caplog.records) + + +def test_bbox_tuple_to_wkt_discards_non_finite(): + assert discovery_ops._bbox_4326_tuple_to_wkt((math.nan, 40.0, 10.0, 41.0)) is None + assert discovery_ops._bbox_4326_tuple_to_wkt((10.0, 40.0, math.inf, 41.0)) is None + + +def test_bbox_tuple_to_wkt_discards_inverted_bbox(): + assert discovery_ops._bbox_4326_tuple_to_wkt((10.5, 40.0, 10.0, 41.0)) is None + + +def test_bbox_tuple_to_wkt_rounds_to_five_decimal_places(): + # the frontend's TerraDraw silently drops features with coordinates over 9 + # decimal places, so stored bboxes must come out rounded + wkt = discovery_ops._bbox_4326_tuple_to_wkt( + (-9.370596897483855, 40.51113433598111, -9.159202392776214, 40.805561916939304) + ) + for coordinate in shapely.from_wkt(wkt).exterior.coords: + for value in coordinate: + assert value == round(value, 5) + + +@pytest.mark.integration +@pytest.mark.asyncio +async def test_discovery_extracts_metadata(db_session_maker, admin_user, discovery_env): + _write_geotiff( + discovery_env["archive_root"] / _MISSION_RELATIVE_PATH / "s01/sub/grid.tif" + ) + + collector = await _run_discovery( + db_session_maker, + discovery_env["mission"].id, + discovery_env["settings"], + admin_user, + ) + + ended = _ended_events(collector) + assert len(ended) == 1 + assert ended[0].succeeded is True + + records = await _get_mission_records(db_session_maker, discovery_env["mission"].id) + assert len(records) == 1 + record = records[0] + assert record.name["en"] == "grid.tif" + assert record.description["en"].startswith("Auto-extracted: GTiff raster") + assert record.bbox_4326 is not None + # the synthetic EPSG:3763 grid reprojects to central Portugal + minx, miny, maxx, maxy = to_shape(record.bbox_4326).bounds + assert -8.2 < minx <= maxx < -8.0 + assert 39.6 < miny <= maxy < 39.8 + async with db_session_maker() as session: + statement = sqlmodel.select(models.RecordAsset).where( + models.RecordAsset.survey_related_record_id == record.id + ) + assets = (await session.exec(statement)).all() + assert len(assets) == 1 + assert assets[0].relative_path == "s01/sub/grid.tif" + # the extracted summary goes on the record only; the asset stays empty + assert assets[0].description["en"] == "" + + +@pytest.mark.integration +@pytest.mark.asyncio +async def test_discovery_is_idempotent( + db_session_maker, admin_user, discovery_env, monkeypatch +): + _write_geotiff( + discovery_env["archive_root"] / _MISSION_RELATIVE_PATH / "s01/grid.tif" + ) + extraction_calls = [] + real_dispatch = discovery_ops.extractor_dispatch.dispatch_extractor + + def counting_dispatch(path): + extraction_calls.append(path) + return real_dispatch(path) + + monkeypatch.setattr( + discovery_ops.extractor_dispatch, "dispatch_extractor", counting_dispatch + ) + + await _run_discovery( + db_session_maker, + discovery_env["mission"].id, + discovery_env["settings"], + admin_user, + ) + first_records = await _get_mission_records( + db_session_maker, discovery_env["mission"].id + ) + second_collector = await _run_discovery( + db_session_maker, + discovery_env["mission"].id, + discovery_env["settings"], + admin_user, + ) + + records = await _get_mission_records(db_session_maker, discovery_env["mission"].id) + assert len(records) == 1 + assert records[0].id == first_records[0].id + # dedup runs BEFORE extraction, so the second run must not extract at all + assert len(extraction_calls) == 1 + # and the second run must not even attempt a record creation (a failed + # attempt would surface as a succeeded=False modification event) + assert not any( + isinstance(e, event_schemas.ResourceModificationEvent) and e.succeeded is False + for e in second_collector.events + ) + + +@pytest.mark.integration +@pytest.mark.asyncio +async def test_discovery_dedup_is_scoped_per_mission( + db_session_maker, admin_user, discovery_env +): + # the same mission-relative file path exists in two missions: each mission + # must get its own record (dedup must not be global) + second_mission_relative_path = "surveys/test-mission-2" + _write_geotiff( + discovery_env["archive_root"] / _MISSION_RELATIVE_PATH / "s01/grid.tif" + ) + _write_geotiff( + discovery_env["archive_root"] / second_mission_relative_path / "s01/grid.tif" + ) + async with db_session_maker() as session: + second_mission = await _create_mission( + session, + admin_user, + identifiers.ProjectId(discovery_env["project"].id), + second_mission_relative_path, + name_suffix=" 2", + ) + + for mission_id in (discovery_env["mission"].id, second_mission.id): + await _run_discovery( + db_session_maker, mission_id, discovery_env["settings"], admin_user + ) + + first_records = await _get_mission_records( + db_session_maker, discovery_env["mission"].id + ) + second_records = await _get_mission_records(db_session_maker, second_mission.id) + assert len(first_records) == 1 + assert len(second_records) == 1 + + +@pytest.mark.integration +@pytest.mark.asyncio +async def test_discovery_survives_extraction_failure( + db_session_maker, admin_user, discovery_env, caplog +): + bad_file = discovery_env["archive_root"] / _MISSION_RELATIVE_PATH / "s01/bad.tif" + bad_file.parent.mkdir(parents=True, exist_ok=True) + bad_file.write_bytes(b"this is not a tif") + + with caplog.at_level(logging.WARNING, logger=discovery_ops.logger.name): + collector = await _run_discovery( + db_session_maker, + discovery_env["mission"].id, + discovery_env["settings"], + admin_user, + ) + + ended = _ended_events(collector) + assert len(ended) == 1 + assert ended[0].succeeded is True + assert any( + "Metadata extraction failed" in record.message for record in caplog.records + ) + + records = await _get_mission_records(db_session_maker, discovery_env["mission"].id) + assert len(records) == 1 + record = records[0] + assert record.name["en"] == "bad.tif" + assert record.description["en"] == "" + assert record.bbox_4326 is None + + +@pytest.mark.integration +@pytest.mark.asyncio +async def test_discovery_passes_temporal_extent_through( + db_session_maker, admin_user, discovery_env, monkeypatch +): + # KMALL and SEG-Y do emit temporal dates, but a faked dispatcher keeps this + # test about the passthrough itself rather than about any one extractor + target = discovery_env["archive_root"] / _MISSION_RELATIVE_PATH / "s01/dated.tif" + target.parent.mkdir(parents=True, exist_ok=True) + target.write_bytes(b"content is irrelevant, dispatch is faked") + + def fake_dispatch(path): + return extractor_schemas.RasterMetadata( + driver="GTiff", + width=1, + height=1, + band_count=1, + bbox_4326=(-8.2, 39.6, -8.1, 39.7), + temporal_extent_begin=dt.date(2024, 9, 28), + temporal_extent_end=dt.date(2024, 9, 29), + ) + + monkeypatch.setattr( + discovery_ops.extractor_dispatch, "dispatch_extractor", fake_dispatch + ) + + await _run_discovery( + db_session_maker, + discovery_env["mission"].id, + discovery_env["settings"], + admin_user, + ) + + records = await _get_mission_records(db_session_maker, discovery_env["mission"].id) + assert len(records) == 1 + record = records[0] + assert record.temporal_extent_begin == dt.date(2024, 9, 28) + assert record.temporal_extent_end == dt.date(2024, 9, 29) + # the stored bbox is the metadata bbox expanded by the always-applied buffer + assert to_shape(record.bbox_4326).bounds == pytest.approx( + (-8.2001, 39.5999, -8.0999, 39.7001) + ) + + +@pytest.mark.integration +@pytest.mark.asyncio +async def test_discovery_skips_configuration_with_null_foreign_keys( + db_session_maker, admin_user, discovery_env, caplog +): + # a configuration with a NULL FK exists in real deployments; it must be + # skipped with a warning while the other configurations keep working + async with db_session_maker() as session: + session.add( + models.AssetDiscoveryConfiguration( + id=uuid.uuid4(), + name="broken config", + relative_path_regexp=r"s01/.*\.tif", + workflow_stage_id=None, + dataset_category_id=None, + ) + ) + await session.commit() + _write_geotiff( + discovery_env["archive_root"] / _MISSION_RELATIVE_PATH / "s01/grid.tif" + ) + + with caplog.at_level(logging.WARNING, logger=discovery_ops.logger.name): + collector = await _run_discovery( + db_session_maker, + discovery_env["mission"].id, + discovery_env["settings"], + admin_user, + ) + + ended = _ended_events(collector) + assert len(ended) == 1 + assert ended[0].succeeded is True + assert any( + "Skipping asset discovery configuration" in record.message + for record in caplog.records + ) + records = await _get_mission_records(db_session_maker, discovery_env["mission"].id) + assert len(records) == 1