Coverage for lunar_toolkit/sfs_report.py: 95%
777 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 02:58 +0000
1"""Validation report of an SfS work dir: metrics, figures and HTML (Bertone et al. 2026, Sec. 3)."""
3from __future__ import annotations
5import base64
6import html
7import json
8import logging
9import math
10import warnings
11from collections.abc import Callable, Iterable, Sequence
12from dataclasses import dataclass, field, fields
13from pathlib import Path
14from typing import Any
16import numpy as np
17import pandas as pd
19from . import asp
20from .sfs_products import Grid, illumination_maps, product_name, read_raster, slope_deg
21from .tools import ToolError
23log = logging.getLogger(__name__)
25SERIES = ("#2a78d6", "#eb6834", "#1baf7a")
26COVERAGE_COLORS = ("#000000", "#e34948", "#eda100", "#b8b8b8")
27COVERAGE_LABELS = ("0", "1", "2", ">=3")
28DIFF_EDGES = np.linspace(-5.0, 5.0, 101)
29NORM_EDGES = np.linspace(-5.0, 5.0, 101)
30SLOPE_EDGES = (0, 5, 10, 15, 20, 25, 30, 35, 40, 90)
31BA_STAGES = (("ba", "ba/run"), ("refine", "ba_refine/run"))
33SECTIONS = {
34 "ba": "Bundle adjustment residuals and match offsets (Fig. 5)",
35 "removed": "Removed images (Sec. 2)",
36 "coverage": "Illumination coverage (Figs. 3, 6)",
37 "prior": "SfS minus prior DEM (Table 2, Fig. 7 left)",
38 "lola": "SfS minus LOLA points vs slope (Table 2, Fig. 7)",
39 "ensemble": "Ensemble uncertainty (Fig. 8)",
40 "uncertainty": "Residuals normalised by the LDEM uncertainty maps (Fig. 9)",
41 "published": "Comparison with the published SDEM (GLDELEV)",
42}
45@dataclass
46class ReportInputs:
47 """Files and settings of a report; anything None is discovered or its section skipped.
49 Parameters
50 ----------
51 work_dir : Path
52 An [`SfsPipeline`][lunar_toolkit.sfs.SfsPipeline] work directory.
53 sfs_dem, prior_dem : Path, optional
54 Default ``products/<ROI>_GLDELEV_<RES>.tif`` (else ``sfs/run-DEM-final.tif``)
55 and ``dem/prior-<RES>m.tif``.
56 region_id, resolution : optional
57 Product name parts (default: ``config.yaml``, else the work dir name and 5).
58 shadow_threshold, max_offset_p85 : float, optional
59 Default from ``config.yaml`` (or the Bertone preset).
60 ldec, toterr, slperr, ldem_slope : Path, optional
61 PGDA LOLA count, elevation / slope uncertainty and slope rasters.
62 lola_points : Path, optional
63 ASP CSV of LOLA points
64 ([`lunar_toolkit.reference.write_asp_csv`][lunar_toolkit.reference.write_asp_csv]).
65 lola_xyzi : Path, optional
66 PGDA ``.xyzi`` point cloud, converted to *lola_points* when those are not given.
67 published_dem : Path, optional
68 Published Bertone et al. SDEM (``<ROI>_GLDELEV_005.tif``).
69 slope_edges : sequence of float
70 Slope bins (deg) for the differences vs slope.
71 slope_split : float
72 Slope (deg) separating the low / high-slope ensemble regimes.
73 preview_size : int
74 Longer side (px) of the grid maps and coverage comparisons are drawn on.
75 """
77 work_dir: Path
78 sfs_dem: Path | None = None
79 prior_dem: Path | None = None
80 region_id: str | None = None
81 resolution: float | None = None
82 shadow_threshold: float | None = None
83 max_offset_p85: float | None = None
84 ldec: Path | None = None
85 toterr: Path | None = None
86 slperr: Path | None = None
87 ldem_slope: Path | None = None
88 lola_points: Path | None = None
89 lola_xyzi: Path | None = None
90 published_dem: Path | None = None
91 slope_edges: Sequence[float] = SLOPE_EDGES
92 slope_split: float = 20.0
93 preview_size: int = 800
95 @classmethod
96 def from_work_dir(cls, work_dir: Path, **overrides: Any) -> ReportInputs:
97 """Inputs with defaults read from ``<work_dir>/config.yaml`` and the standard paths."""
98 work_dir = Path(work_dir)
99 kw: dict[str, Any] = {k: v for k, v in overrides.items() if v is not None}
100 cfg = _work_dir_config(work_dir)
101 if cfg is not None:
102 kw.setdefault("resolution", cfg.resolution)
103 kw.setdefault("region_id", cfg.region_id)
104 kw.setdefault("shadow_threshold", cfg.options("sfs").get("shadow_threshold"))
105 kw.setdefault("max_offset_p85", cfg.max_offset_p85)
106 inp = cls(work_dir, **{k: v for k, v in kw.items() if v is not None})
107 inp.resolution = inp.resolution or 5.0
108 inp.region_id = inp.region_id or work_dir.resolve().name
109 inp.shadow_threshold = 0.002 if inp.shadow_threshold is None else inp.shadow_threshold
110 inp.max_offset_p85 = 1.8 if inp.max_offset_p85 is None else inp.max_offset_p85
111 if inp.sfs_dem is None: 111 ↛ 113line 111 didn't jump to line 113 because the condition on line 111 was always true
112 inp.sfs_dem = _first_existing(inp.product("ELEV"), work_dir / "sfs" / "run-DEM-final.tif")
113 if inp.prior_dem is None: 113 ↛ 115line 113 didn't jump to line 115 because the condition on line 113 was always true
114 inp.prior_dem = _first_existing(work_dir / "dem" / f"prior-{inp.resolution:g}m.tif")
115 return inp
117 def product(self, code: str, *suffix: Any) -> Path:
118 """``<work_dir>/products/<ROI>_GLD<code>_<RES>[_suffix].tif``."""
119 return (
120 self.work_dir
121 / "products"
122 / product_name(self.region_id or "", code, self.resolution or 5, *suffix)
123 )
126def _work_dir_config(work_dir: Path) -> Any:
127 path = work_dir / "config.yaml"
128 if not path.exists():
129 return None
130 from .config import ConfigError, load_config
132 try:
133 return load_config(path).sfs
134 except ConfigError as exc:
135 log.warning("Ignoring %s: %s", path, exc)
136 return None
139def _first_existing(*paths: Path) -> Path | None:
140 return next((p for p in paths if p.exists()), None)
143@dataclass
144class Histogram:
145 """Counts of *values* in *edges*; values beyond the edges go to the outer bins."""
147 label: str
148 edges: np.ndarray
149 counts: np.ndarray
151 @classmethod
152 def of(cls, label: str, values: np.ndarray, edges: np.ndarray | Sequence[float]) -> Histogram:
153 e = np.asarray(edges, float)
154 v = np.asarray(values, float)
155 v = np.clip(v[np.isfinite(v)], e[0], e[-1])
156 return cls(label, e, np.histogram(v, e)[0])
158 @property
159 def total(self) -> int:
160 return int(self.counts.sum())
163@dataclass
164class Metrics:
165 """Everything a report shows: JSON summary, tables and histograms per section."""
167 inputs: ReportInputs
168 summary: dict[str, Any] = field(default_factory=dict)
169 tables: dict[str, pd.DataFrame] = field(default_factory=dict)
170 histograms: dict[str, list[Histogram]] = field(default_factory=dict)
171 hist2d: dict[str, tuple[np.ndarray, np.ndarray, np.ndarray]] = field(default_factory=dict)
172 skipped: dict[str, str] = field(default_factory=dict)
175def stats(values: Any) -> dict[str, float]:
176 """Count, mean, median, std, NMAD, RMS, range, 5/95th percentiles and the same of ``|values|``."""
177 v = np.asarray(values, float).ravel()
178 v = v[np.isfinite(v)]
179 if v.size == 0:
180 return {"count": 0}
181 med = float(np.median(v))
182 a = np.abs(v)
183 return {
184 "count": int(v.size),
185 "mean": float(v.mean()),
186 "median": med,
187 "std": float(v.std()),
188 "nmad": float(1.4826 * np.median(np.abs(v - med))),
189 "rms": float(np.sqrt(np.mean(v * v))),
190 "min": float(v.min()),
191 "max": float(v.max()),
192 "p05": float(np.percentile(v, 5)),
193 "p95": float(np.percentile(v, 95)),
194 "abs_mean": float(a.mean()),
195 "abs_median": float(np.median(a)),
196 "abs_std": float(a.std()),
197 }
200def split_stats(values: np.ndarray, ldec: np.ndarray | None = None) -> dict[str, dict[str, float]]:
201 """[`stats`][lunar_toolkit.sfs_report.stats] of all values and, with a LOLA count raster, of ``LDEC > 0``
202 / ``LDEC = 0`` pixels."""
203 out = {"all": stats(values)}
204 if ldec is not None: 204 ↛ 208line 204 didn't jump to line 208 because the condition on line 204 was always true
205 ok = np.isfinite(ldec)
206 out["ldec>0"] = stats(values[ok & (ldec > 0)])
207 out["ldec=0"] = stats(values[ok & (ldec <= 0)])
208 return out
211def binned(values: np.ndarray, by: np.ndarray, edges: Sequence[float]) -> pd.DataFrame:
212 """Per bin of *by*: count, median, quartiles, mean ``|values|`` and RMS of *values*."""
213 v, b = np.asarray(values, float).ravel(), np.asarray(by, float).ravel()
214 ok = np.isfinite(v) & np.isfinite(b)
215 v, b = v[ok], b[ok]
216 rows = []
217 for lo, hi in zip(edges[:-1], edges[1:]):
218 x = v[(b >= lo) & (b < hi)]
219 row: dict[str, Any] = {"bin_min": lo, "bin_max": hi, "count": int(x.size)}
220 if x.size:
221 p25, p50, p75 = np.percentile(x, [25, 50, 75])
222 row.update(
223 median=p50,
224 p25=p25,
225 p75=p75,
226 abs_mean=float(np.abs(x).mean()),
227 rms=float(np.sqrt(np.mean(x * x))),
228 )
229 rows.append(row)
230 return pd.DataFrame(
231 rows, columns=["bin_min", "bin_max", "count", "median", "p25", "p75", "abs_mean", "rms"]
232 )
235def normal_summary(z: np.ndarray) -> dict[str, float]:
236 """How close normalised residuals are to a unit Gaussian: std, NMAD and fractions within 1-3 sigma."""
237 z = np.asarray(z, float).ravel()
238 z = z[np.isfinite(z)]
239 if z.size == 0: 239 ↛ 240line 239 didn't jump to line 240 because the condition on line 239 was never true
240 return {"count": 0}
241 a = np.abs(z)
242 return {
243 "count": int(z.size),
244 "mean": float(z.mean()),
245 "std": float(z.std()),
246 "nmad": float(1.4826 * np.median(np.abs(z - np.median(z)))),
247 "within_1": float(np.mean(a < 1)),
248 "within_2": float(np.mean(a < 2)),
249 "within_3": float(np.mean(a < 3)),
250 }
253def coverage_fractions(solar_bins: np.ndarray, valid: np.ndarray | None = None) -> dict[str, float]:
254 """Fractions of pixels lit from 0, 1, 2 and >=3 non-redundant Sun directions (and >=1, >=2)."""
255 b = np.asarray(solar_bins)
256 b = b[valid] if valid is not None else b.ravel()
257 n = b.size
258 if n == 0: 258 ↛ 259line 258 didn't jump to line 259 because the condition on line 258 was never true
259 return {"pixels": 0}
260 f = [float(np.mean(b == k)) for k in (0, 1, 2)] + [float(np.mean(b >= 3))]
261 return {
262 "pixels": int(n),
263 "frac_0": f[0],
264 "frac_1": f[1],
265 "frac_2": f[2],
266 "frac_3plus": f[3],
267 "frac_ge1": f[1] + f[2] + f[3],
268 "frac_ge2": f[2] + f[3],
269 }
272def read_ba_reports(prefix: Path) -> pd.DataFrame | None:
273 """Per-image bundle-adjustment report of ``<prefix>-*``: residuals before / after and match offsets.
275 Columns ``image``, ``initial_mean`` / ``initial_median`` / ``initial_count``,
276 ``final_*`` (pixels) and ``offset_p25`` ... ``offset_p95``, ``offset_count`` (m).
277 None when none of the reports exist.
278 """
279 files = {
280 "initial": Path(f"{prefix}-initial_residuals_stats.txt"),
281 "final": Path(f"{prefix}-final_residuals_stats.txt"),
282 "offsets": Path(f"{prefix}-mapproj_match_offset_stats.txt"),
283 }
284 if not any(p.exists() for p in files.values()):
285 return None
286 rows: dict[str, dict[str, Any]] = {}
287 for key in ("initial", "final"):
288 if files[key].exists(): 288 ↛ 287line 288 didn't jump to line 287 because the condition on line 288 was always true
289 for name, s in asp.read_residual_stats(files[key]).items():
290 rows.setdefault(name, {}).update(
291 {f"{key}_mean": s.mean, f"{key}_median": s.median, f"{key}_count": s.count}
292 )
293 if files["offsets"].exists(): 293 ↛ 298line 293 didn't jump to line 298 because the condition on line 293 was always true
294 for name, o in asp.read_mapproj_offset_stats(files["offsets"]).items():
295 rows.setdefault(name, {}).update(
296 {f"offset_p{p}": getattr(o, f"p{p}") for p in (25, 50, 75, 85, 95)}, offset_count=o.count
297 )
298 cols = ["image"] + [f"{k}_{s}" for k in ("initial", "final") for s in ("mean", "median", "count")]
299 cols += [f"offset_p{p}" for p in (25, 50, 75, 85, 95)] + ["offset_count"]
300 df = pd.DataFrame([{"image": k, **v} for k, v in sorted(rows.items())], columns=cols)
301 offsets = [c for c in cols if c.startswith("offset_p")]
302 df[offsets] = df[offsets].where(df[offsets] >= 0)
303 return df
306def ba_summary(table: pd.DataFrame, resolution: float, max_offset_p85: float) -> dict[str, Any]:
307 """Fig. 5 numbers: per-image median residual mean +- std before / after, offset percentiles."""
308 out: dict[str, Any] = {"images": int(len(table))}
309 for key in ("initial", "final"):
310 med = table[f"{key}_median"].dropna()
311 if len(med): 311 ↛ 309line 311 didn't jump to line 309 because the condition on line 311 was always true
312 out[f"{key}_median_px"] = {
313 "mean": float(med.mean()),
314 "std": float(med.std(ddof=0)),
315 "median": float(med.median()),
316 }
317 p85, p95 = table["offset_p85"].dropna(), table["offset_p95"].dropna()
318 if len(p85): 318 ↛ 322line 318 didn't jump to line 322 because the condition on line 318 was always true
319 out["offset_p85_m"] = {"median": float(p85.median()), "max": float(p85.max())}
320 out["offset_p85_over_threshold"] = int((p85 > max_offset_p85).sum())
321 out["offset_threshold_m"] = max_offset_p85
322 if len(p95): 322 ↛ 325line 322 didn't jump to line 325 because the condition on line 322 was always true
323 out["offset_p95_m"] = {"median": float(p95.median()), "max": float(p95.max())}
324 out["frac_p95_below_resolution"] = float((p95 < resolution).mean())
325 return out
328def read_convergence(path: Path) -> pd.DataFrame:
329 """``<prefix>-convergence_angles.txt``: ``left``, ``right``, ``p25``, ``p50``, ``p75`` (deg), ``matches``."""
330 rows = []
331 for line in Path(path).read_text().splitlines():
332 parts = line.split()
333 if len(parts) >= 6 and not parts[0].startswith("#"):
334 rows.append(
335 [Path(parts[0]).name, Path(parts[1]).name, *map(float, parts[2:5]), int(float(parts[5]))]
336 )
337 return pd.DataFrame(rows, columns=["left", "right", "p25", "p50", "p75", "matches"])
340def read_pointmap_residuals(path: Path) -> np.ndarray:
341 """``mean_residual`` column (pixels) of a ``*_residuals_pointmap.csv`` (GCP rows included)."""
342 df = pd.read_csv(path, comment="#", header=None, usecols=[3], skipinitialspace=True)
343 return pd.to_numeric(df[3], errors="coerce").to_numpy(float)
346def removed_images(work_dir: Path) -> pd.DataFrame:
347 """``selection/removed.csv`` with a ``cause`` column (matches, offset, residuals, shadow, other)."""
348 path = Path(work_dir) / "selection" / "removed.csv"
349 cols = ["product_id", "stage", "reason", "cause"]
350 if not path.exists():
351 return pd.DataFrame(columns=cols)
352 df = pd.read_csv(path, dtype=str).fillna("")
353 df["cause"] = df["reason"].map(_cause)
354 return df[cols]
357def _cause(reason: str) -> str:
358 for key, cause in ( 358 ↛ 366line 358 didn't jump to line 366 because the loop on line 358 didn't complete
359 ("matches", "matches"),
360 ("offset", "offset"),
361 ("residuals", "residuals"),
362 ("shadow", "shadow"),
363 ):
364 if key in reason:
365 return cause
366 return "other"
369def selection_table(work_dir: Path) -> pd.DataFrame | None:
370 """``selection/selection.csv`` (every kept image; ``selected`` as bool), or None."""
371 path = Path(work_dir) / "selection" / "selection.csv"
372 if not path.exists():
373 return None
374 df = pd.read_csv(path)
375 df["selected"] = df["selected"].astype(str).str.lower().isin(("true", "1"))
376 return df
379def preview_grid(grid: Grid, max_size: int) -> Grid:
380 """*grid* coarsened by an integer factor so its longer side is at most *max_size*."""
381 f = max(1, math.ceil(max(grid.width, grid.height) / max_size))
382 if f == 1:
383 return grid
384 from rasterio.transform import Affine
386 t = grid.transform
387 coarse = Affine(t.a * f, t.b * f, t.c, t.d * f, t.e * f, t.f)
388 return Grid(grid.crs, coarse, math.ceil(grid.width / f), math.ceil(grid.height / f))
391def solar_bins_from_images(
392 paths: Sequence[Path],
393 azimuths: Sequence[float],
394 grid: Grid,
395 shadow_threshold: float,
396 bin_width: float = 30.0,
397) -> np.ndarray:
398 """Non-redundant Sun directions per pixel of *grid* from mapprojected images (``illumination_maps``)."""
399 from .sfs_products import ProductImage
401 ims = [ProductImage(Path(p), float(a)) for p, a in zip(paths, azimuths)]
402 return illumination_maps(ims, grid, shadow_threshold=shadow_threshold, bin_width=bin_width).solar_bins
405def sample(array: np.ndarray, grid: Grid, x: np.ndarray, y: np.ndarray) -> np.ndarray:
406 """Nearest-pixel values of *array* on *grid* at map coordinates; NaN outside."""
407 inv = ~grid.transform
408 x, y = np.asarray(x, float), np.asarray(y, float)
409 col, row = inv.a * x + inv.b * y + inv.c, inv.d * x + inv.e * y + inv.f
410 c, r = np.floor(col).astype(int), np.floor(row).astype(int)
411 ok = (c >= 0) & (c < grid.width) & (r >= 0) & (r < grid.height)
412 out = np.full(c.shape, np.nan)
413 out[ok] = array[r[ok], c[ok]]
414 return out
417def lola_differences(
418 dem: Path, points_csv: Path, out_prefix: Path, z: np.ndarray | None = None, grid: Grid | None = None
419) -> pd.DataFrame:
420 """DEM minus LOLA points (``x``, ``y`` polar stereographic m; ``diff`` m).
422 Uses ``geodiff`` ([`compare_to_points`][lunar_toolkit.sfs_products.compare_to_points]);
423 without ASP, the nearest DEM pixel of each point (*z* on *grid*, read when not given).
424 """
425 from . import reference as ref
426 from .geometry import ps_xy
427 from .sfs_products import compare_to_points
429 try:
430 cmp = compare_to_points(
431 dem, points_csv, out_prefix, csv_format=ref.ASP_CSV_FORMAT, csv_srs=ref.PS_SRS
432 )
433 except ToolError as exc:
434 log.info("geodiff unavailable (%s); sampling the DEM at the points", exc)
435 else:
436 xy = [ps_xy(la, lo, "S") for lo, la in zip(cmp.lon, cmp.lat)]
437 x = np.array([p[0] for p in xy], float)
438 y = np.array([p[1] for p in xy], float)
439 return pd.DataFrame({"x": x, "y": y, "diff": cmp.diff})
440 pts = np.loadtxt(points_csv, delimiter=",", comments="#", ndmin=2)
441 if z is None or grid is None:
442 z, grid = read_raster(dem)
443 return pd.DataFrame(
444 {"x": pts[:, 0], "y": pts[:, 1], "diff": sample(z, grid, pts[:, 0], pts[:, 1]) - pts[:, 2]}
445 )
448class _Terrain:
449 """The SfS DEM, prior and derived arrays on the SfS grid, read once (float32)."""
451 def __init__(self, inputs: ReportInputs):
452 self.inputs = inputs
453 self._cache: dict[str, np.ndarray | None] = {}
454 self.grid = Grid.from_raster(inputs.sfs_dem) if inputs.sfs_dem else None
456 def get(self, name: str, make: Callable[[], np.ndarray | None]) -> np.ndarray | None:
457 if name not in self._cache:
458 self._cache[name] = make()
459 return self._cache[name]
461 def raster(self, path: Path | None) -> np.ndarray | None:
462 if path is None or self.grid is None or not Path(path).exists():
463 return None
464 return read_raster(path, self.grid)[0].astype(np.float32)
466 @property
467 def z(self) -> np.ndarray | None:
468 return self.get("z", lambda: self.raster(self.inputs.sfs_dem))
470 @property
471 def prior(self) -> np.ndarray | None:
472 return self.get("prior", lambda: self.raster(self.inputs.prior_dem))
474 @property
475 def diff(self) -> np.ndarray | None:
476 return self.get("diff", lambda: None if self.z is None or self.prior is None else self.z - self.prior)
478 @property
479 def slope(self) -> np.ndarray | None:
480 def make() -> np.ndarray | None:
481 res = self.inputs.resolution or 5
482 product = self.raster(_first_existing(self.inputs.product("TSLP", res)))
483 if product is not None: 483 ↛ 484line 483 didn't jump to line 484 because the condition on line 483 was never true
484 return product
485 return None if self.z is None or self.grid is None else _slope(self.z, self.grid.res)
487 return self.get("slope", make)
490def _slope(z: np.ndarray, res: float) -> np.ndarray:
491 return slope_deg(z.astype(np.float64), res).astype(np.float32)
494def _quiet(func: Callable[..., Any], *args: Any) -> Any:
495 with warnings.catch_warnings():
496 warnings.simplefilter("ignore", RuntimeWarning)
497 return func(*args)
500def _ba_section(m: Metrics, t: _Terrain) -> None:
501 inp = m.inputs
502 ws = open_workspace(inp.work_dir)
503 try:
504 for stage, rel in BA_STAGES:
505 prefix = inp.work_dir / rel
506 run = Path(rel).parent.name
507 from_ws = ws is not None and run in ws.runs()
508 table = ws.ba_report(run) if ws is not None and from_ws else read_ba_reports(prefix)
509 if table is None:
510 continue
511 m.tables[f"ba_{stage}"] = table
512 summary = ba_summary(table, inp.resolution or 5, inp.max_offset_p85 or 1.8)
513 summary["source"] = "workspace" if from_ws else "files"
514 conv = Path(f"{prefix}-convergence_angles.txt")
515 angles = None
516 if ws is not None and from_ws:
517 angles = ws.convergence(run)
518 elif conv.exists(): 518 ↛ 520line 518 didn't jump to line 520 because the condition on line 518 was always true
519 angles = read_convergence(conv)
520 if angles is not None and len(angles): 520 ↛ 522line 520 didn't jump to line 522 because the condition on line 520 was always true
521 summary["convergence_deg"] = {"pairs": int(len(angles)), **_med(angles["p50"])}
522 m.summary.setdefault("ba", {})[stage] = summary
523 for key in ("initial", "final"):
524 pm = Path(f"{prefix}-{key}_residuals_pointmap.csv")
525 values = None
526 if ws is not None and from_ws:
527 values = ws.point_residuals(run, key)
528 values = values if len(values) else None
529 elif pm.exists(): 529 ↛ 531line 529 didn't jump to line 531 because the condition on line 529 was always true
530 values = read_pointmap_residuals(pm)
531 if values is not None: 531 ↛ 523line 531 didn't jump to line 523 because the condition on line 531 was always true
532 m.histograms.setdefault(f"pointmap_{stage}", []).append(
533 Histogram.of(key, values, np.logspace(-2, 2, 81))
534 )
535 finally:
536 if ws is not None:
537 ws.close()
538 if "ba" not in m.summary:
539 m.skipped["ba"] = "no bundle-adjustment reports under ba/ or ba_refine/"
542def open_workspace(work_dir: Path) -> Any:
543 """The work dir's ``lunar.duckdb`` [`Workspace`][lunar_toolkit.workspace.Workspace] (read-only), or None.
545 None when there is none or it cannot be opened (e.g. locked by a running pipeline).
546 """
547 import duckdb
549 from .workspace import Workspace, workspace_path
551 path = workspace_path(work_dir)
552 if not path.exists():
553 return None
554 try:
555 return Workspace(path, read_only=True, spatial=False)
556 except duckdb.Error as exc:
557 log.warning("Not reading %s: %s", path, exc)
558 return None
561def _med(s: pd.Series) -> dict[str, float]:
562 s = s.dropna()
563 return {"median": float(s.median()), "max": float(s.max())} if len(s) else {}
566def _removed_section(m: Metrics, t: _Terrain) -> None:
567 wd = m.inputs.work_dir
568 removed = removed_images(wd)
569 sel = selection_table(wd)
570 ba = m.tables.get("ba_ba")
571 counts: dict[str, Any] = {}
572 if ba is not None:
573 counts["bundle_adjusted"] = int(len(ba))
574 if sel is not None:
575 counts["kept"] = int(len(sel))
576 counts["selected"] = int(sel["selected"].sum())
577 if removed.empty and not counts:
578 m.skipped["removed"] = "no selection/removed.csv or selection.csv"
579 return
580 m.tables["removed"] = removed
581 counts["removed"] = int(len(removed))
582 counts["by_stage"] = {str(k): int(v) for k, v in removed["stage"].value_counts().items()}
583 counts["by_cause"] = {str(k): int(v) for k, v in removed["cause"].value_counts().items()}
584 m.summary["removed"] = counts
587def _coverage_section(m: Metrics, t: _Terrain) -> None:
588 inp = m.inputs
589 out: dict[str, Any] = {}
590 product = inp.product("SBCT")
591 if product.exists():
592 bins = read_raster(product)[0]
593 out["selected_full_res"] = coverage_fractions(np.nan_to_num(bins, nan=0.0))
594 sel = selection_table(inp.work_dir)
595 if sel is not None and t.grid is not None:
596 sel = sel[sel["mapprojected"].map(lambda p: isinstance(p, str) and Path(p).exists())]
597 sel = sel[np.isfinite(pd.to_numeric(sel["sun_azimuth"], errors="coerce"))]
598 coarse = preview_grid(t.grid, inp.preview_size)
599 groups = {"kept": sel, "selected": sel[sel["selected"]]}
600 for name, df in groups.items():
601 if len(df): 601 ↛ 600line 601 didn't jump to line 600 because the condition on line 601 was always true
602 bins = solar_bins_from_images(
603 df["mapprojected"].tolist(),
604 df["sun_azimuth"].astype(float).tolist(),
605 coarse,
606 inp.shadow_threshold or 0.0,
607 )
608 out[name] = {"images": int(len(df)), **coverage_fractions(bins)}
609 if not out:
610 m.skipped["coverage"] = "no solar-bin product (GLDSBCT) or mapprojected SfS images"
611 return
612 m.summary["coverage"] = out
615def _prior_section(m: Metrics, t: _Terrain) -> None:
616 inp = m.inputs
617 if t.diff is None:
618 m.skipped["prior"] = "SfS DEM or prior DEM missing"
619 return
620 ldec = t.get("ldec", lambda: t.raster(inp.ldec))
621 m.summary["prior"] = split_stats(t.diff, ldec)
622 hists = [Histogram.of("SfS - prior (all)", t.diff, DIFF_EDGES)]
623 if ldec is not None: 623 ↛ 628line 623 didn't jump to line 628 because the condition on line 623 was always true
624 hists += [
625 Histogram.of("LDEC > 0", t.diff[np.isfinite(ldec) & (ldec > 0)], DIFF_EDGES),
626 Histogram.of("LDEC = 0", t.diff[np.isfinite(ldec) & (ldec <= 0)], DIFF_EDGES),
627 ]
628 m.histograms["prior_diff"] = hists
629 if t.slope is not None: 629 ↛ exitline 629 didn't return from function '_prior_section' because the condition on line 629 was always true
630 m.tables["prior_vs_slope"] = binned(t.diff, t.slope, inp.slope_edges)
633def _lola_section(m: Metrics, t: _Terrain, work: Path) -> None:
634 inp = m.inputs
635 points = inp.lola_points
636 if points is None and inp.lola_xyzi is not None and t.grid is not None and Path(inp.lola_xyzi).exists(): 636 ↛ 637line 636 didn't jump to line 637 because the condition on line 636 was never true
637 from . import reference as ref
639 x0, y0, x1, y1 = t.grid.bounds
640 points = ref.write_asp_csv(Path(inp.lola_xyzi), work / "lola_points.csv", ref.PsBox(x0, y0, x1, y1))
641 if points is None or not Path(points).exists() or inp.sfs_dem is None or t.z is None or t.grid is None:
642 m.skipped["lola"] = "no LOLA points (--lola-points / --lola-xyzi) or SfS DEM"
643 return
644 pts = lola_differences(inp.sfs_dem, Path(points), work / "lola", t.z, t.grid)
645 m.summary["lola"] = stats(pts["diff"])
646 m.histograms["lola_diff"] = [Histogram.of("SfS - LOLA", pts["diff"].to_numpy(), DIFF_EDGES)]
647 if t.slope is not None: 647 ↛ exitline 647 didn't return from function '_lola_section' because the condition on line 647 was always true
648 slope = sample(t.slope, t.grid, pts["x"].to_numpy(), pts["y"].to_numpy())
649 m.tables["lola_vs_slope"] = binned(pts["diff"].to_numpy(), slope, inp.slope_edges)
650 ok = np.isfinite(slope) & np.isfinite(pts["diff"].to_numpy())
651 counts, xe, ye = np.histogram2d(
652 slope[ok],
653 np.clip(pts["diff"].to_numpy()[ok], DIFF_EDGES[0], DIFF_EDGES[-1]),
654 [np.linspace(0, 45, 46), DIFF_EDGES],
655 )
656 m.hist2d["lola_vs_slope"] = (xe, ye, counts)
659def _ensemble_section(m: Metrics, t: _Terrain) -> None:
660 inp = m.inputs
661 paths = {k: inp.product(c, "ENS") for k, c in (("elevation_std", "ISGM"), ("slope_std", "SSGM"))}
662 if not paths["elevation_std"].exists():
663 m.skipped["ensemble"] = "no ensemble products (*_GLDISGM_*_ENS.tif)"
664 return
665 esd, grid = read_raster(paths["elevation_std"])
666 ssd = read_raster(paths["slope_std"], grid)[0] if paths["slope_std"].exists() else None
667 out: dict[str, Any] = {"elevation_std_m": stats(esd)}
668 if ssd is not None: 668 ↛ 670line 668 didn't jump to line 670 because the condition on line 668 was always true
669 out["slope_std_deg"] = stats(ssd)
670 m.histograms["ensemble"] = [Histogram.of("elevation std (m)", esd, np.linspace(0, 5, 101))]
671 slope = None
672 if inp.sfs_dem is not None: 672 ↛ 675line 672 didn't jump to line 675 because the condition on line 672 was always true
673 z = read_raster(inp.sfs_dem, grid)[0]
674 slope = slope_deg(z, grid.res)
675 if slope is not None: 675 ↛ 682line 675 didn't jump to line 682 because the condition on line 675 was always true
676 lo, hi = slope < inp.slope_split, slope >= inp.slope_split
677 out["by_slope"] = {
678 f"slope<{inp.slope_split:g}": stats(esd[lo]),
679 f"slope>={inp.slope_split:g}": stats(esd[hi]),
680 }
681 m.tables["ensemble_vs_slope"] = binned(esd, slope, inp.slope_edges)
682 height_error = inp.product("ISGM")
683 if height_error.exists(): 683 ↛ 688line 683 didn't jump to line 688 because the condition on line 683 was always true
684 he = read_raster(height_error, grid)[0]
685 ratio = _quiet(np.divide, esd, np.where(he > 0, he, np.nan))
686 out["ensemble_over_height_error"] = stats(ratio)
687 m.histograms["ensemble"].append(Histogram.of("ASP height error (m)", he, np.linspace(0, 5, 101)))
688 m.summary["ensemble"] = out
691def _uncertainty_section(m: Metrics, t: _Terrain) -> None:
692 inp = m.inputs
693 out: dict[str, Any] = {}
694 hists = []
695 toterr = t.raster(inp.toterr)
696 if toterr is not None and t.diff is not None:
697 z = _quiet(np.divide, t.diff, np.where(toterr > 0, toterr, np.nan))
698 out["elevation"] = normal_summary(z)
699 hists.append(Histogram.of("elevation", z, NORM_EDGES))
700 slperr = t.raster(inp.slperr)
701 if slperr is not None and t.slope is not None:
702 ldem_slope = t.raster(inp.ldem_slope)
703 if ldem_slope is None and t.prior is not None and t.grid is not None: 703 ↛ 705line 703 didn't jump to line 705 because the condition on line 703 was always true
704 ldem_slope = _slope(t.prior, t.grid.res)
705 if ldem_slope is not None: 705 ↛ 709line 705 didn't jump to line 709 because the condition on line 705 was always true
706 z = _quiet(np.divide, t.slope - ldem_slope, np.where(slperr > 0, slperr, np.nan))
707 out["slope"] = normal_summary(z)
708 hists.append(Histogram.of("slope", z, NORM_EDGES))
709 if not out:
710 m.skipped["uncertainty"] = "no PGDA toterr / slperr maps (--toterr / --slperr / --region)"
711 return
712 m.summary["uncertainty"] = out
713 m.histograms["normalised"] = hists
716def _published_section(m: Metrics, t: _Terrain) -> None:
717 inp = m.inputs
718 pub = t.raster(inp.published_dem)
719 if pub is None or t.z is None or t.grid is None:
720 m.skipped["published"] = "no published SDEM (--published-dem)"
721 return
722 diff = t.z - pub
723 t._cache["published_diff"] = diff
724 out = {"elevation": stats(diff)}
725 if t.slope is not None: 725 ↛ 727line 725 didn't jump to line 727 because the condition on line 725 was always true
726 out["slope"] = stats(t.slope - _slope(pub, t.grid.res))
727 m.summary["published"] = out
728 m.histograms["published_diff"] = [Histogram.of("SfS - published", diff, DIFF_EDGES)]
731def metrics(
732 inputs: ReportInputs | Path, work: Path | None = None, sections: Iterable[str] | None = None
733) -> Metrics:
734 """All report metrics; sections whose inputs are missing are listed in ``skipped``.
736 Parameters
737 ----------
738 inputs : ReportInputs or Path
739 Inputs, or a work dir
740 ([`ReportInputs.from_work_dir`][lunar_toolkit.sfs_report.ReportInputs.from_work_dir]).
741 work : Path, optional
742 Scratch directory for converted LOLA points and ``geodiff`` output
743 (default: ``<work_dir>/report/work``).
744 sections : iterable of str, optional
745 Subset of `SECTIONS` (default: all).
746 """
747 inp = inputs if isinstance(inputs, ReportInputs) else ReportInputs.from_work_dir(inputs)
748 work = Path(work) if work else inp.work_dir / "report" / "work"
749 m = Metrics(inp)
750 t = _Terrain(inp)
751 runners: dict[str, Callable[[], None]] = {
752 "ba": lambda: _ba_section(m, t),
753 "removed": lambda: _removed_section(m, t),
754 "coverage": lambda: _coverage_section(m, t),
755 "prior": lambda: _prior_section(m, t),
756 "lola": lambda: _lola_section(m, t, work),
757 "ensemble": lambda: _ensemble_section(m, t),
758 "uncertainty": lambda: _uncertainty_section(m, t),
759 "published": lambda: _published_section(m, t),
760 }
761 wanted = list(sections) if sections is not None else list(SECTIONS)
762 unknown = set(wanted) - set(runners)
763 if unknown:
764 raise ValueError(f"Unknown report sections {sorted(unknown)}; choose from {', '.join(SECTIONS)}")
765 for name in wanted:
766 log.info("Report: %s", SECTIONS[name])
767 runners[name]()
768 m.summary["inputs"] = {f.name: _jsonable(getattr(inp, f.name)) for f in fields(inp)}
769 m.summary["skipped"] = dict(m.skipped)
770 return m
773def _jsonable(v: Any) -> Any:
774 if isinstance(v, Path):
775 return str(v)
776 if isinstance(v, (tuple, list)):
777 return [_jsonable(x) for x in v]
778 if isinstance(v, (np.floating, np.integer)): 778 ↛ 779line 778 didn't jump to line 779 because the condition on line 778 was never true
779 return v.item()
780 if isinstance(v, dict):
781 return {str(k): _jsonable(x) for k, x in v.items()}
782 if isinstance(v, float) and not math.isfinite(v): 782 ↛ 783line 782 didn't jump to line 783 because the condition on line 782 was never true
783 return None
784 return v
787def write_tables(m: Metrics, out_dir: Path) -> dict[str, Path]:
788 """``tables/<name>.csv`` for every table, plus ``metrics.json`` (summary and histograms)."""
789 out_dir = Path(out_dir)
790 (out_dir / "tables").mkdir(parents=True, exist_ok=True)
791 paths = {}
792 for name, df in m.tables.items():
793 paths[name] = out_dir / "tables" / f"{name}.csv"
794 df.to_csv(paths[name], index=False)
795 for name, df in summary_tables(m).items():
796 paths[name] = out_dir / "tables" / f"{name}.csv"
797 df.to_csv(paths[name], index=False)
798 doc = {
799 **_jsonable(m.summary),
800 "histograms": {
801 k: [{"label": h.label, "edges": h.edges.tolist(), "counts": h.counts.tolist()} for h in v]
802 for k, v in m.histograms.items()
803 },
804 }
805 paths["metrics"] = out_dir / "metrics.json"
806 paths["metrics"].write_text(json.dumps(doc, indent=1))
807 return paths
810def summary_tables(m: Metrics) -> dict[str, pd.DataFrame]:
811 """Table 2-style rows (one per comparison) and the coverage fractions as tables."""
812 out = {}
813 rows = []
814 for label, s in (
815 *((f"SfS - prior ({k})", v) for k, v in m.summary.get("prior", {}).items()),
816 *((("SfS - LOLA points", m.summary["lola"]),) if "lola" in m.summary else ()),
817 *((("SfS - published", m.summary["published"]["elevation"]),) if "published" in m.summary else ()),
818 ):
819 rows.append({"comparison": label, **s})
820 if rows:
821 out["elevation_differences"] = pd.DataFrame(rows)
822 cov = m.summary.get("coverage", {})
823 if cov:
824 out["coverage"] = pd.DataFrame([{"images_set": k, **v} for k, v in cov.items()])
825 return out
828def _figure(nrows: int = 1, ncols: int = 1, size: tuple[float, float] = (6.4, 4.2)) -> tuple[Any, Any]:
829 import matplotlib
831 matplotlib.use("Agg")
832 import matplotlib.pyplot as plt
834 fig, ax = plt.subplots(nrows, ncols, figsize=size, constrained_layout=True)
835 return fig, ax
838def _save(fig: Any, path: Path) -> Path:
839 import matplotlib.pyplot as plt
841 path.parent.mkdir(parents=True, exist_ok=True)
842 fig.savefig(path, dpi=110)
843 plt.close(fig)
844 return path
847def _style(ax: Any) -> None:
848 ax.grid(True, color="#e0e0e0", linewidth=0.6)
849 ax.set_axisbelow(True)
850 for side in ("top", "right"):
851 ax.spines[side].set_visible(False)
854def _hist_figure(
855 hists: Sequence[Histogram],
856 xlabel: str,
857 path: Path,
858 *,
859 density: bool = True,
860 logx: bool = False,
861 gaussians: Sequence[tuple[float, str]] = (),
862) -> Path:
863 fig, ax = _figure()
864 for h, color in zip(hists, SERIES):
865 y = h.counts / max(h.total, 1) / np.diff(h.edges) if density else h.counts
866 ax.stairs(y, h.edges, color=color, linewidth=1.6, label=f"{h.label} (n={h.total:,})")
867 for sigma, color in gaussians:
868 x = np.linspace(hists[0].edges[0], hists[0].edges[-1], 400)
869 ax.plot(
870 x,
871 np.exp(-0.5 * (x / sigma) ** 2) / (sigma * math.sqrt(2 * math.pi)),
872 "--",
873 color=color,
874 linewidth=1.2,
875 label=f"Gaussian sigma={sigma:.2f}",
876 )
877 if logx:
878 ax.set_xscale("log")
879 ax.set_xlabel(xlabel)
880 ax.set_ylabel("density" if density else "count")
881 _style(ax)
882 if len(hists) + len(gaussians) > 1:
883 ax.legend(frameon=False, fontsize=8)
884 return _save(fig, path)
887def _map_figure(
888 array: np.ndarray,
889 grid: Grid,
890 path: Path,
891 label: str,
892 *,
893 cmap: str = "cividis",
894 symmetric: bool = False,
895 vmin: float | None = None,
896 vmax: float | None = None,
897) -> Path:
898 fig, ax = _figure(size=(6.0, 5.0))
899 x0, y0, x1, y1 = grid.bounds
900 finite = array[np.isfinite(array)]
901 if symmetric and finite.size:
902 lim = float(np.percentile(np.abs(finite), 98)) or 1.0
903 vmin, vmax = -lim, lim
904 elif finite.size and vmin is None:
905 vmin, vmax = (float(v) for v in np.percentile(finite, [2, 98]))
906 im = ax.imshow(
907 array,
908 extent=(x0 / 1e3, x1 / 1e3, y0 / 1e3, y1 / 1e3),
909 cmap=cmap,
910 vmin=vmin,
911 vmax=vmax,
912 interpolation="nearest",
913 )
914 fig.colorbar(im, ax=ax, label=label, shrink=0.85)
915 _map_axes(ax)
916 return _save(fig, path)
919def _map_axes(ax: Any) -> None:
920 from matplotlib.ticker import MaxNLocator
922 ax.xaxis.set_major_locator(MaxNLocator(5))
923 ax.yaxis.set_major_locator(MaxNLocator(5))
924 ax.set_xlabel("x (km, polar stereographic)")
925 ax.set_ylabel("y (km)")
928def hillshade_array(z: np.ndarray, res: float, azimuth: float = 315.0, elevation: float = 45.0) -> np.ndarray:
929 """Lambertian hillshade (0-1) of *z* with the Sun at *azimuth* / *elevation* (deg)."""
930 from .sfs_products import gradient
932 p, q = gradient(z, res)
933 az, el = math.radians(azimuth), math.radians(elevation)
934 sx, sy, sz = math.sin(az) * math.cos(el), math.cos(az) * math.cos(el), math.sin(el)
935 return np.clip((-p * sx - q * sy + sz) / np.sqrt(1 + p * p + q * q), 0, 1)
938def figures(m: Metrics, out_dir: Path) -> dict[str, Path]:
939 """PNG figures of every section with data, in ``<out_dir>/figures``; returns name -> path."""
940 import matplotlib.colors as mcolors
942 inp = m.inputs
943 fig_dir = Path(out_dir) / "figures"
944 out: dict[str, Path] = {}
945 if inp.sfs_dem is not None and Path(inp.sfs_dem).exists(): 945 ↛ 981line 945 didn't jump to line 981 because the condition on line 945 was always true
946 grid = preview_grid(Grid.from_raster(inp.sfs_dem), inp.preview_size)
947 z = read_raster(inp.sfs_dem, grid)[0]
948 out["hillshade"] = _map_figure(
949 hillshade_array(z, grid.res),
950 grid,
951 fig_dir / "hillshade.png",
952 "hillshade (Sun az 315, el 45 deg)",
953 cmap="gray",
954 vmin=0,
955 vmax=1,
956 )
957 out["slope_map"] = _map_figure(
958 slope_deg(z, grid.res),
959 grid,
960 fig_dir / "slope_map.png",
961 "slope (deg)",
962 cmap="cividis",
963 vmin=0,
964 vmax=40,
965 )
966 if inp.prior_dem is not None and Path(inp.prior_dem).exists(): 966 ↛ 971line 966 didn't jump to line 971 because the condition on line 966 was always true
967 d = z - read_raster(inp.prior_dem, grid)[0]
968 out["prior_diff_map"] = _map_figure(
969 d, grid, fig_dir / "prior_diff_map.png", "SfS - prior (m)", cmap="RdBu_r", symmetric=True
970 )
971 if inp.published_dem is not None and Path(inp.published_dem).exists():
972 d = z - read_raster(inp.published_dem, grid)[0]
973 out["published_diff_map"] = _map_figure(
974 d,
975 grid,
976 fig_dir / "published_diff_map.png",
977 "SfS - published SDEM (m)",
978 cmap="RdBu_r",
979 symmetric=True,
980 )
981 sbct = inp.product("SBCT")
982 if sbct.exists(): 982 ↛ 999line 982 didn't jump to line 999 because the condition on line 982 was always true
983 g = preview_grid(Grid.from_raster(sbct), inp.preview_size)
984 bins = np.nan_to_num(_nearest(sbct, g), nan=0.0)
985 cmap = mcolors.ListedColormap(COVERAGE_COLORS)
986 fig, ax = _figure(size=(6.0, 5.0))
987 x0, y0, x1, y1 = g.bounds
988 im = ax.imshow(
989 np.minimum(bins, 3),
990 extent=(x0 / 1e3, x1 / 1e3, y0 / 1e3, y1 / 1e3),
991 cmap=cmap,
992 norm=mcolors.BoundaryNorm([-0.5, 0.5, 1.5, 2.5, 3.5], 4),
993 interpolation="nearest",
994 )
995 cb = fig.colorbar(im, ax=ax, ticks=[0, 1, 2, 3], shrink=0.85, label="non-redundant Sun directions")
996 cb.ax.set_yticklabels(COVERAGE_LABELS)
997 _map_axes(ax)
998 out["coverage_map"] = _save(fig, fig_dir / "coverage_map.png")
999 cov = {k: v for k, v in m.summary.get("coverage", {}).items() if v.get("pixels")}
1000 if cov:
1001 out["coverage_bars"] = _coverage_bars(cov, fig_dir / "coverage_bars.png")
1002 for stage in ("ba", "refine"):
1003 table = m.tables.get(f"ba_{stage}")
1004 if table is not None:
1005 out.update(_ba_figures(table, stage, inp, fig_dir))
1006 if f"pointmap_{stage}" in m.histograms:
1007 out[f"pointmap_{stage}"] = _hist_figure(
1008 m.histograms[f"pointmap_{stage}"],
1009 "triangulated point mean residual (px)",
1010 fig_dir / f"pointmap_residuals_{stage}.png",
1011 logx=True,
1012 )
1013 diffs = m.histograms.get("prior_diff", []) + m.histograms.get("lola_diff", [])
1014 if diffs:
1015 out["diff_histograms"] = _hist_figure(
1016 diffs[:3], "elevation difference (m)", fig_dir / "diff_histograms.png"
1017 )
1018 if "lola_diff" in m.histograms and len(diffs) > 3:
1019 out["lola_histogram"] = _hist_figure(
1020 m.histograms["lola_diff"], "SfS - LOLA (m)", fig_dir / "lola_histogram.png"
1021 )
1022 for name in ("lola_vs_slope", "prior_vs_slope", "ensemble_vs_slope"):
1023 if name in m.tables:
1024 out[name] = _slope_figure(m, name, fig_dir / f"{name}.png")
1025 if "ensemble" in m.histograms:
1026 out["ensemble_histogram"] = _hist_figure(
1027 m.histograms["ensemble"], "elevation uncertainty (m)", fig_dir / "ensemble_histogram.png"
1028 )
1029 esd = inp.product("ISGM", "ENS")
1030 g = preview_grid(Grid.from_raster(esd), inp.preview_size)
1031 out["ensemble_map"] = _map_figure(
1032 read_raster(esd, g)[0],
1033 g,
1034 fig_dir / "ensemble_map.png",
1035 "ensemble elevation std (m)",
1036 cmap="cividis",
1037 )
1038 for h in m.histograms.get("normalised", []):
1039 sigma = m.summary["uncertainty"][h.label]["std"]
1040 out[f"normalised_{h.label}"] = _hist_figure(
1041 [h],
1042 f"(SfS - LDEM) / LDEM uncertainty, {h.label}",
1043 fig_dir / f"normalised_{h.label}.png",
1044 gaussians=[(1.0, "#e34948"), (sigma, "#008300")],
1045 )
1046 return out
1049def _nearest(path: Path, grid: Grid) -> np.ndarray:
1050 import rasterio
1051 from rasterio.enums import Resampling
1053 with rasterio.open(path) as src:
1054 data = src.read(1, out_shape=(grid.height, grid.width), resampling=Resampling.nearest, masked=True)
1055 return data.astype(np.float64).filled(np.nan)
1058def _coverage_bars(cov: dict[str, dict[str, float]], path: Path) -> Path:
1059 fig, ax = _figure(size=(6.4, 2.2 + 0.5 * len(cov)))
1060 names = list(cov)
1061 left = np.zeros(len(names))
1062 keys = ("frac_3plus", "frac_2", "frac_1", "frac_0")
1063 for key, color, label in zip(keys, COVERAGE_COLORS[::-1], ("3+", "2", "1", "0")):
1064 vals = np.array([100 * cov[n].get(key, 0.0) for n in names])
1065 ax.barh(
1066 names, vals, left=left, color=color, edgecolor="white", linewidth=1, label=f"{label} directions"
1067 )
1068 left += vals
1069 ax.set_xlim(0, 100)
1070 ax.set_xlabel("pixels (%)")
1071 ax.legend(frameon=False, fontsize=8, ncol=4, loc="upper center", bbox_to_anchor=(0.5, -0.35))
1072 _style(ax)
1073 return _save(fig, path)
1076def _ba_figures(table: pd.DataFrame, stage: str, inp: ReportInputs, fig_dir: Path) -> dict[str, Path]:
1077 out = {}
1078 med = [table[c].dropna().to_numpy() for c in ("initial_median", "final_median")]
1079 if all(len(x) for x in med): 1079 ↛ 1098line 1079 didn't jump to line 1098 because the condition on line 1079 was always true
1080 edges = np.logspace(-2, 2, 61)
1081 hists = [Histogram.of(k, x, edges) for k, x in zip(("before BA", "after BA"), med)]
1082 out[f"ba_residuals_{stage}"] = _hist_figure(
1083 hists,
1084 "per-image median reprojection error (px)",
1085 fig_dir / f"ba_residuals_{stage}.png",
1086 density=False,
1087 logx=True,
1088 )
1089 both = table.dropna(subset=["initial_median", "final_median"])
1090 fig, ax = _figure(size=(5.0, 4.6))
1091 ax.scatter(both["initial_median"], both["final_median"], s=10, color=SERIES[0])
1092 ax.set_xscale("log")
1093 ax.set_yscale("log")
1094 ax.set_xlabel("median residual before BA (px)")
1095 ax.set_ylabel("median residual after BA (px)")
1096 _style(ax)
1097 out[f"ba_before_after_{stage}"] = _save(fig, fig_dir / f"ba_before_after_{stage}.png")
1098 p85, p95 = table["offset_p85"].dropna(), table["offset_p95"].dropna()
1099 if len(p85): 1099 ↛ 1126line 1099 didn't jump to line 1126 because the condition on line 1099 was always true
1100 hi = max(float(np.nanmax(p95)) if len(p95) else 0.0, (inp.resolution or 5) * 1.2, 2.0)
1101 edges = np.linspace(0, hi, 61)
1102 fig, ax = _figure()
1103 for s, color, label in ((p85, SERIES[0], "p85"), (p95, SERIES[1], "p95")):
1104 ax.stairs(
1105 np.histogram(np.clip(s, 0, hi), edges)[0], edges, color=color, linewidth=1.6, label=label
1106 )
1107 ax.axvline(
1108 inp.max_offset_p85 or 1.8,
1109 color="#52514e",
1110 linestyle="--",
1111 linewidth=1,
1112 label=f"p85 threshold {inp.max_offset_p85:g} m",
1113 )
1114 ax.axvline(
1115 inp.resolution or 5,
1116 color="#52514e",
1117 linestyle=":",
1118 linewidth=1,
1119 label=f"SfS resolution {inp.resolution:g} m",
1120 )
1121 ax.set_xlabel("per-image mapprojected match offset (m)")
1122 ax.set_ylabel("images")
1123 ax.legend(frameon=False, fontsize=8)
1124 _style(ax)
1125 out[f"match_offsets_{stage}"] = _save(fig, fig_dir / f"match_offsets_{stage}.png")
1126 return out
1129def _slope_figure(m: Metrics, name: str, path: Path) -> Path:
1130 t = m.tables[name].dropna(subset=["median"])
1131 fig, ax = _figure()
1132 if name in m.hist2d:
1133 xe, ye, counts = m.hist2d[name]
1134 ax.pcolormesh(
1135 xe, ye, np.where(counts.T > 0, counts.T, np.nan), cmap="Greys", norm="log", shading="flat"
1136 )
1137 mid = (t["bin_min"] + np.minimum(t["bin_max"], 45)) / 2
1138 ax.fill_between(mid, t["p25"], t["p75"], color=SERIES[0], alpha=0.25, linewidth=0, label="interquartile")
1139 ax.plot(mid, t["median"], "-o", color=SERIES[0], markersize=4, linewidth=1.6, label="median")
1140 ylabel = {
1141 "lola_vs_slope": "SfS - LOLA (m)",
1142 "prior_vs_slope": "SfS - prior (m)",
1143 "ensemble_vs_slope": "ensemble elevation std (m)",
1144 }[name]
1145 ax.set_xlabel("slope (deg)")
1146 ax.set_ylabel(ylabel)
1147 ax.legend(frameon=False, fontsize=8)
1148 _style(ax)
1149 return _save(fig, path)
1152_CSS = """
1153body { background: #ffffff; color: #111111; font-family: Georgia, 'Times New Roman', serif;
1154 max-width: 60em; margin: 2em auto; padding: 0 16px; line-height: 1.45; }
1155h1, h2, h3 { font-weight: normal; }
1156h2 { border-bottom: 1px solid #999999; padding-bottom: 0.2em; margin-top: 2em; }
1157table { border-collapse: collapse; margin: 0.8em 0; font-size: 0.85em; display: block; overflow-x: auto; }
1158th, td { border-top: 1px solid #cccccc; border-bottom: 1px solid #cccccc; padding: 0.25em 0.6em; text-align: right; }
1159th { background: #f2f2f2; }
1160td:first-child, th:first-child { text-align: left; }
1161figure { margin: 1em 0; }
1162figure img { max-width: 100%; height: auto; border: 1px solid #dddddd; }
1163figcaption, .note { font-size: 0.85em; color: #444444; }
1164code { font-size: 0.9em; }
1165"""
1168def _fmt(v: Any) -> str:
1169 if isinstance(v, float):
1170 if not math.isfinite(v):
1171 return ""
1172 return f"{v:.4g}" if abs(v) < 1e5 else f"{v:,.0f}"
1173 if isinstance(v, (int, np.integer)) and not isinstance(v, bool):
1174 return f"{v:,}"
1175 return html.escape(str(v))
1178def _table_html(df: pd.DataFrame, max_rows: int = 40) -> str:
1179 head = "".join(f"<th>{html.escape(str(c))}</th>" for c in df.columns)
1180 body = "".join(
1181 "<tr>" + "".join(f"<td>{_fmt(v)}</td>" for v in row) + "</tr>"
1182 for row in df.head(max_rows).itertuples(index=False)
1183 )
1184 more = (
1185 f'<p class="note">First {max_rows} of {len(df)} rows; see the CSV.</p>' if len(df) > max_rows else ""
1186 )
1187 return f"<table><tr>{head}</tr>{body}</table>{more}"
1190def _dict_table(d: dict[str, Any]) -> str:
1191 rows = []
1192 for k, v in d.items():
1193 if isinstance(v, dict) and all(not isinstance(x, dict) for x in v.values()):
1194 v = ", ".join(f"{a}={_fmt(b)}" for a, b in v.items())
1195 rows.append(f"<tr><td>{html.escape(k)}</td><td>{v}</td></tr>")
1196 elif not isinstance(v, dict): 1196 ↛ 1192line 1196 didn't jump to line 1192 because the condition on line 1196 was always true
1197 rows.append(f"<tr><td>{html.escape(k)}</td><td>{_fmt(v)}</td></tr>")
1198 return f"<table>{''.join(rows)}</table>" if rows else ""
1201def _stats_rows(d: dict[str, dict[str, Any]]) -> str:
1202 rows = [{"subset": k, **v} for k, v in d.items() if isinstance(v, dict)]
1203 return _table_html(pd.DataFrame(rows)) if rows else ""
1206def _img(path: Path, caption: str) -> str:
1207 data = base64.b64encode(Path(path).read_bytes()).decode("ascii")
1208 return f'<figure><img src="data:image/png;base64,{data}" alt="{html.escape(caption)}"><figcaption>{html.escape(caption)}</figcaption></figure>'
1211_SECTION_FIGURES = {
1212 "ba": (
1213 "ba_residuals_ba",
1214 "ba_before_after_ba",
1215 "match_offsets_ba",
1216 "pointmap_ba",
1217 "ba_residuals_refine",
1218 "ba_before_after_refine",
1219 "match_offsets_refine",
1220 "pointmap_refine",
1221 ),
1222 "coverage": ("coverage_map", "coverage_bars"),
1223 "prior": ("hillshade", "slope_map", "prior_diff_map", "diff_histograms", "prior_vs_slope"),
1224 "lola": ("lola_histogram", "lola_vs_slope"),
1225 "ensemble": ("ensemble_map", "ensemble_histogram", "ensemble_vs_slope"),
1226 "uncertainty": ("normalised_elevation", "normalised_slope"),
1227 "published": ("published_diff_map",),
1228}
1229_SECTION_TABLES = {
1230 "ba": ("ba_ba", "ba_refine"),
1231 "removed": ("removed",),
1232 "prior": ("prior_vs_slope",),
1233 "lola": ("lola_vs_slope",),
1234 "ensemble": ("ensemble_vs_slope",),
1235}
1236_CAPTIONS = {
1237 "ba_residuals": "Per-image median reprojection error before and after bundle adjustment (cf. Fig. 5 top right).",
1238 "ba_before_after": "Median reprojection error per image, before vs after bundle adjustment.",
1239 "match_offsets": "Per-image 85th / 95th percentile mapprojected match offsets (cf. Fig. 5 bottom right).",
1240 "pointmap": "Mean reprojection error of the triangulated points, before and after bundle adjustment.",
1241 "coverage_map": "Non-redundant Sun directions (30 deg bins) per pixel of the SfS inputs (cf. Fig. 6).",
1242 "coverage_bars": "Pixels lit from 0, 1, 2 and >=3 directions, kept vs selected images (cf. Fig. 3).",
1243 "hillshade": "Hillshade of the SfS DEM (preview grid).",
1244 "slope_map": "Slope of the SfS DEM on the preview grid (baseline = 2 preview pixels).",
1245 "prior_diff_map": "SfS DEM minus prior DEM.",
1246 "diff_histograms": "Elevation differences (cf. Fig. 7 left); values beyond the axis go to the outer bins.",
1247 "lola_histogram": "SfS DEM minus LOLA points.",
1248 "lola_vs_slope": "SfS minus LOLA points vs slope: median and interquartile range (cf. Fig. 7 right).",
1249 "prior_vs_slope": "SfS minus prior DEM vs slope: median and interquartile range.",
1250 "ensemble_map": "Ensemble elevation standard deviation (cf. Fig. 8 top left).",
1251 "ensemble_histogram": "Ensemble elevation std and ASP height error (cf. Fig. 8 bottom right).",
1252 "ensemble_vs_slope": "Ensemble elevation std vs slope (cf. Fig. 8 bottom left).",
1253 "normalised": "SfS - LDEM residuals over the LDEM clone uncertainty, with a unit Gaussian (red) and "
1254 "a Gaussian of the measured sigma (green) (cf. Fig. 9).",
1255 "published_diff_map": "SfS DEM minus the published Bertone et al. SDEM.",
1256}
1259def _caption(name: str) -> str:
1260 for key in sorted(_CAPTIONS, key=len, reverse=True): 1260 ↛ 1263line 1260 didn't jump to line 1263 because the loop on line 1260 didn't complete
1261 if name.startswith(key):
1262 return _CAPTIONS[key]
1263 return name
1266def write_html(m: Metrics, figs: dict[str, Path], out: Path, tables: dict[str, Path] | None = None) -> Path:
1267 """Self-contained ``report.html``: summary tables, embedded figures, skipped sections."""
1268 inp = m.inputs
1269 out = Path(out)
1270 parts = [
1271 "<!doctype html><html lang='en'><head><meta charset='utf-8'>",
1272 "<meta name='viewport' content='width=device-width, initial-scale=1'>",
1273 f"<title>SfS report {html.escape(inp.region_id or '')}</title><style>{_CSS}</style></head><body>",
1274 f"<h1>SfS DEM report: {html.escape(inp.region_id or '')}</h1>",
1275 f"<p class='note'>Work dir <code>{html.escape(str(inp.work_dir))}</code>; SfS DEM "
1276 f"<code>{html.escape(str(inp.sfs_dem))}</code>; prior <code>{html.escape(str(inp.prior_dem))}</code>; "
1277 f"{inp.resolution:g} m/px. Metrics follow Bertone et al. (2026, PSJ 7:118); "
1278 "tables are in <code>tables/</code> and <code>metrics.json</code>.</p>",
1279 ]
1280 summary_tabs = summary_tables(m)
1281 if "elevation_differences" in summary_tabs:
1282 parts += [
1283 "<h2>Elevation differences (Table 2)</h2>",
1284 _table_html(summary_tabs["elevation_differences"]),
1285 ]
1286 for key, title in SECTIONS.items():
1287 parts.append(f"<h2>{html.escape(title)}</h2>")
1288 if key in m.skipped:
1289 parts.append(f"<p class='note'>Skipped: {html.escape(m.skipped[key])}.</p>")
1290 continue
1291 s = m.summary.get(key, {})
1292 if key == "coverage" and "coverage" in summary_tabs:
1293 parts.append(_table_html(summary_tabs["coverage"]))
1294 elif key in ("prior", "uncertainty", "published") or (key == "ensemble" and "by_slope" in s):
1295 parts.append(_stats_rows({k: v for k, v in s.items() if k != "by_slope"}))
1296 if "by_slope" in s:
1297 parts.append(_stats_rows(s["by_slope"]))
1298 elif key == "ba":
1299 for stage, st in s.items():
1300 parts += [f"<h3>{html.escape(stage)}</h3>", _dict_table(st)]
1301 else:
1302 parts.append(_dict_table(s) if key != "lola" else _stats_rows({"SfS - LOLA": s}))
1303 for name in _SECTION_FIGURES.get(key, ()):
1304 if name in figs:
1305 parts.append(_img(figs[name], _caption(name)))
1306 for name in _SECTION_TABLES.get(key, ()):
1307 if name in m.tables and len(m.tables[name]):
1308 parts += [f"<h3>{html.escape(name)}</h3>", _table_html(m.tables[name])]
1309 if tables: 1309 ↛ 1316line 1309 didn't jump to line 1316 because the condition on line 1309 was always true
1310 links = "".join(
1311 f"<li><a href='{html.escape(p.relative_to(out.parent).as_posix())}'>{html.escape(p.name)}</a></li>"
1312 for p in tables.values()
1313 if p.is_relative_to(out.parent)
1314 )
1315 parts += ["<h2>Files</h2>", f"<ul>{links}</ul>"]
1316 parts.append("</body></html>")
1317 out.parent.mkdir(parents=True, exist_ok=True)
1318 out.write_text("\n".join(parts), encoding="ascii", errors="xmlcharrefreplace")
1319 return out
1322def make_report(
1323 inputs: ReportInputs | Path, out_dir: Path | None = None, sections: Iterable[str] | None = None
1324) -> Path:
1325 """Metrics, tables, figures and ``report.html`` in *out_dir* (default ``<work_dir>/report``)."""
1326 inp = inputs if isinstance(inputs, ReportInputs) else ReportInputs.from_work_dir(inputs)
1327 out_dir = Path(out_dir) if out_dir else inp.work_dir / "report"
1328 m = metrics(inp, out_dir / "work", sections)
1329 tables = write_tables(m, out_dir)
1330 figs = figures(m, out_dir)
1331 return write_html(m, figs, out_dir / "report.html", tables)