|
| 1 | +"""Robust ERT file loading for the desktop workbench (Qt-free). |
| 2 | +
|
| 3 | +A single source of truth for turning an ERT data file into a pygimli |
| 4 | +``DataContainerERT`` with correct geometry, topography, and apparent |
| 5 | +resistivity. Both the single-inversion ERT module and the time-lapse pipeline |
| 6 | +use this so they behave identically. |
| 7 | +
|
| 8 | +Why this exists: pygimli's native ``ert.load`` cannot parse several common |
| 9 | +field formats. The E4D survey export (used here with a ``.ohm`` extension), for |
| 10 | +example, carries a leading index column and no ``# a b m n`` token header, so |
| 11 | +``ert.load`` misreads the index column as coordinates, drops the topography, |
| 12 | +and discards every measurement (``size() == 0``). The device-specific parsers in |
| 13 | +:mod:`PyHydroGeophysX.data_processing.ert_data_agent` (resipy or the embedded |
| 14 | +fallback) handle those layouts; this module wires them to pygimli and falls back |
| 15 | +across loaders so a file never silently loads empty. |
| 16 | +""" |
| 17 | + |
| 18 | +from __future__ import annotations |
| 19 | + |
| 20 | +import datetime as _dt |
| 21 | +import re |
| 22 | +import tempfile |
| 23 | +from pathlib import Path |
| 24 | +from typing import Any, Callable, List, Optional, Sequence, Tuple |
| 25 | + |
| 26 | +import numpy as np |
| 27 | + |
| 28 | +LogFn = Callable[[str], None] |
| 29 | + |
| 30 | + |
| 31 | +def _noop(_msg: str) -> None: |
| 32 | + return None |
| 33 | + |
| 34 | + |
| 35 | +# --------------------------------------------------------------------------- |
| 36 | +# StandardERT -> pygimli DataContainerERT |
| 37 | +# --------------------------------------------------------------------------- |
| 38 | +def electrode_elevation(electrodes) -> np.ndarray: |
| 39 | + """Per-electrode elevation aligned with ``electrodes``. |
| 40 | +
|
| 41 | + ``load_ert_resipy`` carries 2D elevation in ``y`` (its |
| 42 | + ``_normalize_elevation_axis`` moves a flat-y/varying-z profile so elevation |
| 43 | + lives in ``y``); fall back to ``z`` if ``y`` is flat. |
| 44 | + """ |
| 45 | + if not electrodes: |
| 46 | + return np.zeros(0) |
| 47 | + ys = np.array([float(e.y) for e in electrodes]) |
| 48 | + zs = np.array([float(e.z) for e in electrodes]) |
| 49 | + return ys if ys.std() >= zs.std() else zs |
| 50 | + |
| 51 | + |
| 52 | +def standard_to_pg(std): |
| 53 | + """Build a pygimli ``DataContainerERT`` from a ``StandardERT`` for inversion. |
| 54 | +
|
| 55 | + Most instrument loaders report transfer resistance (V/I), not apparent |
| 56 | + resistivity, so apparent resistivity is recovered with geometric factors: |
| 57 | + ``rhoa = R * k``. When the source already provides apparent resistivity it is |
| 58 | + used as-is. pygimli's forward operator uses the same ``data["k"]``, so the |
| 59 | + inversion stays self-consistent. Returns ``None`` if pygimli is unavailable |
| 60 | + or no measurement maps onto the electrode set. |
| 61 | + """ |
| 62 | + try: |
| 63 | + import pygimli as pg |
| 64 | + from pygimli.physics import ert as pg_ert |
| 65 | + except Exception: # noqa: BLE001 |
| 66 | + return None |
| 67 | + data = pg.DataContainerERT() |
| 68 | + id_to_idx = {} |
| 69 | + elecs = std.electrodes or [] |
| 70 | + elev = electrode_elevation(elecs) |
| 71 | + for i, e in enumerate(elecs): |
| 72 | + data.createSensor(pg.Pos(float(e.x), 0.0, float(elev[i]))) |
| 73 | + id_to_idx[int(e.id)] = i |
| 74 | + keys = ("A", "B", "M", "N") |
| 75 | + valid = [ |
| 76 | + o for o in (std.observations or []) |
| 77 | + if o.app_res is not None and all(int(getattr(o.quad, k)) in id_to_idx for k in keys) |
| 78 | + ] |
| 79 | + if not valid: |
| 80 | + return None |
| 81 | + data.resize(len(valid)) |
| 82 | + for name, qk in zip(("a", "b", "m", "n"), keys): |
| 83 | + data.set(name, [id_to_idx[int(getattr(o.quad, qk))] for o in valid]) |
| 84 | + vals = np.array([float(o.app_res) for o in valid], dtype=float) |
| 85 | + data.set("r", vals) |
| 86 | + data.set("err", [float(o.rel_err) if o.rel_err else 0.05 for o in valid]) |
| 87 | + try: |
| 88 | + data["k"] = pg_ert.createGeometricFactors(data, numerical=False) |
| 89 | + k = np.asarray(data["k"], dtype=float) |
| 90 | + except Exception: # noqa: BLE001 |
| 91 | + k = np.ones(len(valid)) |
| 92 | + source = str((std.metadata or {}).get("app_res_source", "")).lower() |
| 93 | + rhoa = vals * k if source == "resistance" else vals |
| 94 | + data.set("rhoa", rhoa) |
| 95 | + data.markValid(data("rhoa") > 0) |
| 96 | + return data |
| 97 | + |
| 98 | + |
| 99 | +# --------------------------------------------------------------------------- |
| 100 | +# Robust single-file loader |
| 101 | +# --------------------------------------------------------------------------- |
| 102 | +_AUTO_NAMES = ("", "auto", "none", "auto-detect", "auto-detect (pygimli)") |
| 103 | +#: Count-prefixed / device formats tried when native pygimli load comes back empty. |
| 104 | +_RECOVERY_INSTRUMENTS = ("E4D", "BERT", "Syscal") |
| 105 | + |
| 106 | + |
| 107 | +def _usable(data) -> bool: |
| 108 | + try: |
| 109 | + return data is not None and int(data.size()) > 0 |
| 110 | + except Exception: # noqa: BLE001 |
| 111 | + return False |
| 112 | + |
| 113 | + |
| 114 | +def _ensure_rhoa(data, log: LogFn = _noop) -> None: |
| 115 | + """Make sure a natively-loaded container has positive apparent resistivity.""" |
| 116 | + try: |
| 117 | + from pygimli.physics import ert as pg_ert |
| 118 | + if data is None: |
| 119 | + return |
| 120 | + have_rhoa = data.haveData("rhoa") and np.any(np.asarray(data["rhoa"], dtype=float) > 0) |
| 121 | + if have_rhoa: |
| 122 | + return |
| 123 | + if not data.haveData("k") or not np.any(np.asarray(data["k"], dtype=float) != 0): |
| 124 | + data["k"] = pg_ert.createGeometricFactors(data, numerical=False) |
| 125 | + k = np.asarray(data["k"], dtype=float) |
| 126 | + if data.haveData("r"): |
| 127 | + data["rhoa"] = np.asarray(data["r"], dtype=float) * k |
| 128 | + elif data.haveData("u") and data.haveData("i"): |
| 129 | + data["rhoa"] = np.asarray(data["u"], dtype=float) / np.asarray(data["i"], dtype=float) * k |
| 130 | + except Exception as exc: # noqa: BLE001 |
| 131 | + log(f"Apparent-resistivity computation skipped: {exc}") |
| 132 | + |
| 133 | + |
| 134 | +def _via_instrument(path: str, instrument: str, electrode_file: Optional[str], |
| 135 | + spacing: Optional[float], log: LogFn): |
| 136 | + try: |
| 137 | + from PyHydroGeophysX.data_processing.ert_data_agent import load_ert_resipy |
| 138 | + except Exception as exc: # noqa: BLE001 |
| 139 | + log(f"ert_data_agent unavailable ({exc}).") |
| 140 | + return None |
| 141 | + proj = tempfile.mkdtemp(prefix="phgx_resipy_") |
| 142 | + try: |
| 143 | + std = load_ert_resipy(project_dir=proj, data_file=str(path), instrument=instrument, |
| 144 | + spacing=spacing, electrode_file=electrode_file) |
| 145 | + except Exception as exc: # noqa: BLE001 |
| 146 | + log(f"Instrument '{instrument}' loader error: {exc}") |
| 147 | + return None |
| 148 | + return standard_to_pg(std) |
| 149 | + |
| 150 | + |
| 151 | +def _via_native(path: str, log: LogFn): |
| 152 | + try: |
| 153 | + from pygimli.physics import ert as pg_ert |
| 154 | + data = pg_ert.load(str(path), verbose=False) |
| 155 | + except Exception as exc: # noqa: BLE001 |
| 156 | + log(f"pygimli native load failed: {exc}") |
| 157 | + return None |
| 158 | + _ensure_rhoa(data, log) |
| 159 | + return data |
| 160 | + |
| 161 | + |
| 162 | +def load_ert_container(path: str, instrument: Optional[str] = None, |
| 163 | + electrode_file: Optional[str] = None, |
| 164 | + spacing: Optional[float] = None, log: LogFn = _noop): |
| 165 | + """Load one ERT file into a pygimli ``DataContainerERT``, robustly. |
| 166 | +
|
| 167 | + An explicit ``instrument`` uses the device parsers (handles index-prefixed |
| 168 | + E4D, BERT topography, Syscal, etc.); ``None``/``"auto"`` uses pygimli's own |
| 169 | + reader. If the chosen path yields no measurements the loader falls back: an |
| 170 | + explicit instrument retries native pygimli; a still-empty result triggers a |
| 171 | + short sweep of count-prefixed device formats so an E4D-style file never |
| 172 | + loads empty. Raises ``ValueError`` if nothing parses. |
| 173 | + """ |
| 174 | + path = str(path) |
| 175 | + inst = str(instrument).strip() if instrument else "" |
| 176 | + is_auto = inst.lower() in _AUTO_NAMES |
| 177 | + |
| 178 | + # 1. explicit instrument |
| 179 | + if not is_auto: |
| 180 | + data = _via_instrument(path, inst, electrode_file, spacing, log) |
| 181 | + if _usable(data): |
| 182 | + return data |
| 183 | + log(f"Instrument '{inst}' parsed no usable measurements; trying pygimli auto-detect.") |
| 184 | + |
| 185 | + # 2. native pygimli |
| 186 | + data = _via_native(path, log) |
| 187 | + if _usable(data): |
| 188 | + return data |
| 189 | + |
| 190 | + # 3. recovery sweep across count-prefixed device formats |
| 191 | + for cand in _RECOVERY_INSTRUMENTS: |
| 192 | + if cand.lower() == inst.lower(): |
| 193 | + continue |
| 194 | + data = _via_instrument(path, cand, electrode_file, spacing, log) |
| 195 | + if _usable(data): |
| 196 | + log(f"Auto-recovered '{Path(path).name}' using instrument='{cand}'.") |
| 197 | + return data |
| 198 | + |
| 199 | + raise ValueError( |
| 200 | + f"No ERT measurements could be parsed from '{Path(path).name}'. " |
| 201 | + f"Pick the matching instrument/format in the loader.") |
| 202 | + |
| 203 | + |
| 204 | +# --------------------------------------------------------------------------- |
| 205 | +# Measurement times from filenames |
| 206 | +# --------------------------------------------------------------------------- |
| 207 | +_DATE_PATTERNS: List[Tuple[re.Pattern, bool]] = [ |
| 208 | + (re.compile(r"(\d{4})[-_]?(\d{2})[-_]?(\d{2})[-_ ]?(\d{2})(\d{2})"), True), # ...YYYY-MM-DD_HHMM |
| 209 | + (re.compile(r"(\d{4})[-_]?(\d{2})[-_]?(\d{2})"), False), # ...YYYY-MM-DD |
| 210 | +] |
| 211 | + |
| 212 | + |
| 213 | +def _parse_date(stem: str) -> Optional[_dt.datetime]: |
| 214 | + for pattern, has_time in _DATE_PATTERNS: |
| 215 | + m = pattern.search(stem) |
| 216 | + if not m: |
| 217 | + continue |
| 218 | + try: |
| 219 | + if has_time: |
| 220 | + return _dt.datetime(int(m.group(1)), int(m.group(2)), int(m.group(3)), |
| 221 | + int(m.group(4)), int(m.group(5))) |
| 222 | + return _dt.datetime(int(m.group(1)), int(m.group(2)), int(m.group(3))) |
| 223 | + except ValueError: |
| 224 | + continue |
| 225 | + return None |
| 226 | + |
| 227 | + |
| 228 | +def measurement_times_for(files: Sequence[str]) -> Tuple[List[float], List[str]]: |
| 229 | + """Derive numeric measurement times + display labels from filenames. |
| 230 | +
|
| 231 | + When every filename embeds a distinct date, times are elapsed days from the |
| 232 | + first acquisition and labels are ``YYYY-MM-DD``. Otherwise falls back to a |
| 233 | + sequential ``1..n`` with index labels. |
| 234 | + """ |
| 235 | + dates = [_parse_date(Path(f).stem) for f in files] |
| 236 | + if files and all(d is not None for d in dates): |
| 237 | + t0 = min(dates) |
| 238 | + times = [round((d - t0).total_seconds() / 86400.0, 4) for d in dates] |
| 239 | + labels = [d.strftime("%Y-%m-%d") for d in dates] |
| 240 | + if len(set(times)) == len(times): # distinct -> usable as a time axis |
| 241 | + return times, labels |
| 242 | + n = len(files) |
| 243 | + return [float(i + 1) for i in range(n)], [str(i + 1) for i in range(n)] |
| 244 | + |
| 245 | + |
| 246 | +# --------------------------------------------------------------------------- |
| 247 | +# Normalize a sequence into clean pygimli files for the core inversion |
| 248 | +# --------------------------------------------------------------------------- |
| 249 | +def normalize_for_timelapse(files: Sequence[str], instrument: Optional[str], |
| 250 | + out_dir: str, log: LogFn = _noop): |
| 251 | + """Load each file robustly and write a clean pygimli ``.dat`` into one folder. |
| 252 | +
|
| 253 | + The core :class:`TimeLapseERTInversion` reloads files with ``ert.load``; the |
| 254 | + normalized files are written in pygimli's native unified format (proper token |
| 255 | + headers, geometric factors, ``rhoa = R*k``, topography in ``z``) so that |
| 256 | + reload is correct. Returns ``(clean_dir, basenames, containers)``; all clean |
| 257 | + files share one folder so the windowed inversion (which takes a directory + |
| 258 | + filenames) works directly. |
| 259 | + """ |
| 260 | + base = Path(out_dir) / "qt_ert_timelapse" / "normalized" |
| 261 | + base.mkdir(parents=True, exist_ok=True) |
| 262 | + for stale in base.glob("step_*.dat"): |
| 263 | + try: |
| 264 | + stale.unlink() |
| 265 | + except OSError: |
| 266 | + pass |
| 267 | + basenames: List[str] = [] |
| 268 | + containers: List[Any] = [] |
| 269 | + for i, f in enumerate(files): |
| 270 | + data = load_ert_container(f, instrument=instrument, log=log) |
| 271 | + name = f"step_{i:03d}.dat" |
| 272 | + data.save(str(base / name)) |
| 273 | + basenames.append(name) |
| 274 | + containers.append(data) |
| 275 | + log(f"Prepared {i + 1}/{len(files)}: {Path(f).name} -> " |
| 276 | + f"{int(data.size())} data, {int(data.sensorCount())} electrodes") |
| 277 | + return str(base), basenames, containers |
0 commit comments