Coverage for lunar_toolkit/workspace.py: 93%
460 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 02:58 +0000
1"""DuckDB workspace (``<work_dir>/lunar.duckdb``): image footprints, stereo overlaps and bundle-adjustment statistics."""
3from __future__ import annotations
5import logging
6import math
7import re
8from collections.abc import Mapping, Sequence
9from pathlib import Path
10from typing import TYPE_CHECKING, Any
12import numpy as np
14from .constants import MOON_RADIUS_KM, parse_nac_id
16if TYPE_CHECKING: # pragma: no cover
17 import pandas as pd
19 from .asp import MatchOffsetStats, ResidualStats
20 from .sfs import SfsConfig
22log = logging.getLogger(__name__)
24DB_NAME = "lunar.duckdb"
25# lunar pairs --engine auto uses DuckDB from this many images up
26PAIRS_ENGINE_THRESHOLD = 300
27# Footprint edges are great circles (as in find_stereo_pairs), densified at this step (deg)
28GREAT_CIRCLE_STEP_DEG = 0.1
29# The sphere and projections of the images.footprint geometries (by images.hemi)
30LONLAT_PROJ4 = f"+proj=longlat +R={MOON_RADIUS_KM * 1000:g} +no_defs"
31PS_PROJ4 = {
32 hemi: f"+proj=stere +lat_0={lat} +lon_0=0 +k=1 +x_0=0 +y_0=0 +R={MOON_RADIUS_KM * 1000:g} +units=m +no_defs"
33 for hemi, lat in (("S", -90), ("N", 90))
34}
35BA_TABLES = (
36 "ba_residual_stats",
37 "ba_pixel_residuals",
38 "ba_match_offsets",
39 "ba_points",
40 "ba_convergence",
41 "ba_camera_offsets",
42)
43_R_M = MOON_RADIUS_KM * 1000.0
44_IDENT_RE = re.compile(r"^[A-Za-z_][A-Za-z0-9_]*$")
46_OBS_TYPES = {
47 "product_id": "VARCHAR",
48 "instrument": "VARCHAR",
49 "product_type": "VARCHAR",
50 "start_time": "TIMESTAMP",
51 "center_lat": "DOUBLE",
52 "center_lon": "DOUBLE",
53 "min_lat": "DOUBLE",
54 "max_lat": "DOUBLE",
55 "lon_west": "DOUBLE",
56 "lon_east": "DOUBLE",
57 "footprint_wkt": "VARCHAR",
58 "incidence": "DOUBLE",
59 "emission": "DOUBLE",
60 "phase": "DOUBLE",
61 "sub_solar_azimuth": "DOUBLE",
62 "north_azimuth": "DOUBLE",
63 "sub_spacecraft_lat": "DOUBLE",
64 "sub_spacecraft_lon": "DOUBLE",
65 "spacecraft_altitude": "DOUBLE",
66 "pixel_scale": "DOUBLE",
67 "url": "VARCHAR",
68 "label_url": "VARCHAR",
69 "source": "VARCHAR",
70}
71_IMAGE_EXTRA = {
72 "obs_key": "VARCHAR",
73 "hemi": "VARCHAR",
74 "cap_lat": "DOUBLE",
75 "cap_lon": "DOUBLE",
76 "footprint_km2": "DOUBLE",
77 "xmin": "DOUBLE",
78 "xmax": "DOUBLE",
79 "ymin": "DOUBLE",
80 "ymax": "DOUBLE",
81}
83_BA_DDL = [
84 "CREATE TABLE IF NOT EXISTS ba_runs (run VARCHAR, prefix VARCHAR, loaded_at TIMESTAMP, files VARCHAR[])",
85 "CREATE TABLE IF NOT EXISTS ba_residual_stats (run VARCHAR, stage VARCHAR, image VARCHAR, "
86 "image_path VARCHAR, product_id VARCHAR, mean_px DOUBLE, median_px DOUBLE, residual_count BIGINT)",
87 "CREATE TABLE IF NOT EXISTS ba_pixel_residuals (run VARCHAR, stage VARCHAR, image VARCHAR, "
88 "product_id VARCHAR, norm_px DOUBLE)",
89 "CREATE TABLE IF NOT EXISTS ba_match_offsets (run VARCHAR, image VARCHAR, image_path VARCHAR, "
90 "product_id VARCHAR, p25 DOUBLE, p50 DOUBLE, p75 DOUBLE, p85 DOUBLE, p95 DOUBLE, match_count BIGINT)",
91 "CREATE TABLE IF NOT EXISTS ba_points (run VARCHAR, stage VARCHAR, lon DOUBLE, lat DOUBLE, "
92 "height_m DOUBLE, mean_residual_px DOUBLE, num_observations INTEGER, kind VARCHAR, hemi VARCHAR, "
93 "x DOUBLE, y DOUBLE)",
94 "CREATE TABLE IF NOT EXISTS ba_convergence (run VARCHAR, left_image VARCHAR, right_image VARCHAR, "
95 "left_product_id VARCHAR, right_product_id VARCHAR, p25 DOUBLE, p50 DOUBLE, p75 DOUBLE, matches BIGINT)",
96 "CREATE TABLE IF NOT EXISTS ba_camera_offsets (run VARCHAR, image VARCHAR, product_id VARCHAR, "
97 "horizontal_m DOUBLE, vertical_m DOUBLE)",
98 "CREATE TABLE IF NOT EXISTS rejection_params (max_offset_p85 DOUBLE, "
99 "max_large_residual_fraction DOUBLE, large_residual_px DOUBLE, min_matches BIGINT)",
100]
102_PS_MACROS = [
103 # areal scale of the polar stereographic projection is _ps_k^2 (k = 1 at the pole)
104 f"CREATE OR REPLACE MACRO _ps_k(x, y) AS 1.0 + (x * x + y * y) / {4.0 * _R_M * _R_M!r}",
105 f"CREATE OR REPLACE MACRO _ps_lat(x, y, hemi) AS (CASE WHEN hemi = 'S' THEN 1.0 ELSE -1.0 END) * "
106 f"(degrees(2.0 * atan(sqrt(x * x + y * y) / {2.0 * _R_M!r})) - 90.0)",
107 "CREATE OR REPLACE MACRO _ps_lon(x, y, hemi) AS "
108 "(degrees(atan2(x, CASE WHEN hemi = 'S' THEN y ELSE -y END)) + 360.0) % 360.0",
109 "CREATE OR REPLACE MACRO _az_diff(a, b) AS least(((b - a) % 360.0 + 360.0) % 360.0, "
110 "360.0 - ((b - a) % 360.0 + 360.0) % 360.0)",
111 "CREATE OR REPLACE MACRO _true_km2(g) AS "
112 "ST_Area(g) / pow(_ps_k(ST_X(ST_Centroid(g)), ST_Y(ST_Centroid(g))), 2) / 1e6",
113]
115# Per-image residual stats before / after BA; {params} is a rejection_params-shaped relation.
116_IMAGE_STATS_SQL = """
117WITH names AS (
118 SELECT DISTINCT run, image FROM (
119 SELECT run, image FROM ba_residual_stats
120 UNION ALL SELECT run, image FROM ba_match_offsets
121 UNION ALL SELECT DISTINCT run, image FROM ba_pixel_residuals
122 )
123),
124raw AS (
125 SELECT r.run, r.image,
126 count(*) FILTER (WHERE r.stage = 'final') AS raw_count,
127 avg(CAST(r.norm_px > p.large_residual_px AS DOUBLE)) FILTER (WHERE r.stage = 'final')
128 AS large_fraction,
129 count(*) FILTER (WHERE r.stage = 'initial') AS raw_count_initial,
130 avg(CAST(r.norm_px > p.large_residual_px AS DOUBLE)) FILTER (WHERE r.stage = 'initial')
131 AS large_fraction_initial
132 FROM ba_pixel_residuals r CROSS JOIN ({params}) p GROUP BY r.run, r.image
133)
134SELECT n.run, n.image, upper(split_part(n.image, '.', 1)) AS product_id,
135 i.mean_px AS initial_mean, i.median_px AS initial_median, i.residual_count AS initial_count,
136 f.mean_px AS final_mean, f.median_px AS final_median, f.residual_count AS final_count,
137 f.median_px - i.median_px AS median_change,
138 o.p25 AS offset_p25, o.p50 AS offset_p50, o.p75 AS offset_p75, o.p85 AS offset_p85,
139 o.p95 AS offset_p95, o.match_count AS offset_count,
140 r.raw_count, r.large_fraction, r.raw_count_initial, r.large_fraction_initial,
141 c.horizontal_m AS camera_offset_horizontal_m, c.vertical_m AS camera_offset_vertical_m
142FROM names n
143LEFT JOIN ba_residual_stats i ON i.run = n.run AND i.image = n.image AND i.stage = 'initial'
144LEFT JOIN ba_residual_stats f ON f.run = n.run AND f.image = n.image AND f.stage = 'final'
145LEFT JOIN ba_match_offsets o ON o.run = n.run AND o.image = n.image
146LEFT JOIN raw r ON r.run = n.run AND r.image = n.image
147LEFT JOIN ba_camera_offsets c ON c.run = n.run AND c.image = n.image
148"""
150# Bertone et al. (2026, Sec. 2) rejection rules, first match wins (as sfs.flag_images).
151_REJECTIONS_SQL = """
152SELECT run, image, product_id, rule, value, threshold,
153 CASE rule
154 WHEN 'matches' THEN printf('%d matches', CAST(value AS BIGINT))
155 WHEN 'offset' THEN printf('match offset p85 %.2f m > %g m', value, threshold)
156 ELSE printf('%.0f%% of residuals > %g px', 100.0 * value, large_residual_px)
157 END AS reason
158FROM (
159 SELECT s.*, p.large_residual_px,
160 CASE WHEN s.min_count < p.min_matches THEN 'matches'
161 WHEN s.offset_p85 > p.max_offset_p85 THEN 'offset'
162 WHEN s.raw_count > 0 AND s.large_fraction > p.max_large_residual_fraction THEN 'residuals'
163 END AS rule,
164 CASE rule WHEN 'matches' THEN s.min_count WHEN 'offset' THEN s.offset_p85
165 ELSE s.large_fraction END AS value,
166 CASE rule WHEN 'matches' THEN p.min_matches WHEN 'offset' THEN p.max_offset_p85
167 ELSE p.max_large_residual_fraction END AS threshold
168 FROM (SELECT *, least(final_count, offset_count) AS min_count FROM ({stats})) s CROSS JOIN ({params}) p
169)
170WHERE rule IS NOT NULL
171"""
173_RUN_SUMMARY_SQL = """
174WITH s AS (
175 SELECT run,
176 count(*) AS images,
177 avg(initial_median) AS initial_median_mean,
178 stddev_pop(initial_median) AS initial_median_std,
179 median(initial_median) AS initial_median_median,
180 avg(final_median) AS final_median_mean,
181 stddev_pop(final_median) AS final_median_std,
182 median(final_median) AS final_median_median,
183 median(offset_p85) AS offset_p85_median,
184 max(offset_p85) AS offset_p85_max,
185 count(*) FILTER (WHERE offset_p85 > (SELECT max_offset_p85 FROM rejection_params))
186 AS offset_p85_over_threshold,
187 median(offset_p95) AS offset_p95_median,
188 max(offset_p95) AS offset_p95_max
189 FROM ba_image_stats GROUP BY run
190),
191rej AS (SELECT run, count(*) AS rejected FROM image_rejections GROUP BY run),
192conv AS (SELECT run, count(*) AS convergence_pairs, median(p50) AS convergence_p50_median
193 FROM ba_convergence GROUP BY run),
194pts AS (
195 SELECT run,
196 count(*) FILTER (WHERE stage = 'initial') AS points_initial,
197 count(*) FILTER (WHERE stage = 'final') AS points_final,
198 median(mean_residual_px) FILTER (WHERE stage = 'initial') AS point_residual_median_initial,
199 median(mean_residual_px) FILTER (WHERE stage = 'final') AS point_residual_median_final
200 FROM ba_points GROUP BY run
201)
202SELECT s.*, coalesce(rej.rejected, 0) AS rejected, conv.convergence_pairs, conv.convergence_p50_median,
203 pts.points_initial, pts.points_final, pts.point_residual_median_initial,
204 pts.point_residual_median_final
205FROM s LEFT JOIN rej USING (run) LEFT JOIN conv USING (run) LEFT JOIN pts USING (run)
206ORDER BY s.run
207"""
209_PARAMS_VIEW = "SELECT * FROM rejection_params"
210_PARAMS_MACRO = (
211 "SELECT p85 AS max_offset_p85, frac AS max_large_residual_fraction, px AS large_residual_px, "
212 "min_matches AS min_matches"
213)
216class WorkspaceError(RuntimeError):
217 """The workspace cannot do what was asked (e.g. no spatial extension)."""
220def workspace_path(work_dir: Path | str) -> Path:
221 """``<work_dir>/lunar.duckdb``, or *work_dir* itself when it names a ``.duckdb`` file."""
222 p = Path(work_dir)
223 return p if p.suffix == ".duckdb" else p / DB_NAME
226def spatial_available() -> bool:
227 """Whether DuckDB's spatial extension loads (installing it when needed)."""
228 global _SPATIAL
229 if _SPATIAL is None: 229 ↛ 237line 229 didn't jump to line 237 because the condition on line 229 was always true
230 import duckdb
232 con = duckdb.connect()
233 try:
234 _SPATIAL = _load_spatial(con) is None
235 finally:
236 con.close()
237 return _SPATIAL
240_SPATIAL: bool | None = None
243def _load_spatial(con: Any) -> str | None:
244 """Load (installing if needed) the spatial extension; the error message, or None."""
245 import duckdb
247 try:
248 con.execute("LOAD spatial")
249 return None
250 except duckdb.Error:
251 pass
252 try:
253 con.execute("INSTALL spatial")
254 con.execute("LOAD spatial")
255 return None
256 except duckdb.Error as exc:
257 return str(exc).splitlines()[0]
260def _ident(name: str) -> str:
261 if not _IDENT_RE.match(name):
262 raise ValueError(f"Not a table name: {name!r}")
263 return name
266def _sql_str(value: Any) -> str:
267 return "'" + str(value).replace("'", "''") + "'"
270def _num(v: float) -> float | None:
271 """*v*, None for NaN (SQL NULL)."""
272 return None if v is None or math.isnan(v) else float(v)
275def _product_id(image: str) -> str:
276 return Path(str(image)).name.split(".")[0].upper()
279class Workspace:
280 """A DuckDB database of images, stereo pairs and bundle-adjustment results.
282 Parameters
283 ----------
284 path : Path or str
285 Database file (``":memory:"`` for a scratch one); a directory means
286 ``<path>/lunar.duckdb``. Parent directories are created.
287 memory_limit : str
288 DuckDB ``memory_limit``.
289 threads : int, optional
290 DuckDB ``threads`` (default: DuckDB's, all cores).
291 read_only : bool
292 Open an existing database read-only (no schema changes).
293 spatial : bool
294 Load the spatial extension (``INSTALL spatial`` once needs network; it is
295 cached in ``~/.duckdb``). Without it the bundle-adjustment tables and views
296 work, while [`add_observations`][lunar_toolkit.workspace.Workspace.add_observations] and
297 [`pairs`][lunar_toolkit.workspace.Workspace.pairs] raise
298 [`WorkspaceError`][lunar_toolkit.workspace.WorkspaceError].
300 Notes
301 -----
302 Tables: ``images`` (observation columns plus ``footprint``, a GEOMETRY in metres
303 in the polar stereographic projection of ``hemi``, `PS_PROJ4`),
304 ``pairs`` (the last [`pairs`][lunar_toolkit.workspace.Workspace.pairs] result), ``ba_runs``,
305 ``ba_residual_stats``, ``ba_pixel_residuals``, ``ba_match_offsets``, ``ba_points``, ``ba_convergence``,
306 ``ba_camera_offsets`` (by ``run``) and ``rejection_params``. Views:
307 ``ba_image_stats`` (per image before / after), ``ba_run_summary`` and
308 ``image_rejections`` (Bertone et al. 2026 rules with ``rejection_params``); the
309 table macro ``image_rejections_at(max_offset_p85, max_large_residual_fraction,
310 large_residual_px, min_matches)`` takes other thresholds.
312 Example
313 -------
314 ::
316 with Workspace("work/malapert") as ws:
317 ws.add_observations(obs)
318 pairs = ws.pairs(aoi={"point": (-85.9, 0.0), "radius_km": 5})
319 ws.load_bundle_adjust("work/malapert/ba/run")
320 ws.query("SELECT * FROM image_rejections")
321 """
323 def __init__(
324 self,
325 path: Path | str,
326 *,
327 memory_limit: str = "1GB",
328 threads: int | None = None,
329 read_only: bool = False,
330 spatial: bool = True,
331 ) -> None:
332 import duckdb
334 self.path = path if str(path) == ":memory:" else workspace_path(path)
335 if not read_only and isinstance(self.path, Path):
336 self.path.parent.mkdir(parents=True, exist_ok=True)
337 config: dict[str, Any] = {"memory_limit": memory_limit}
338 if threads:
339 config["threads"] = int(threads)
340 self.read_only = read_only
341 self.con = duckdb.connect(str(self.path), read_only=read_only, config=config)
342 self.con.execute("SET enable_progress_bar = false")
343 self.spatial_error = _load_spatial(self.con) if spatial else "not requested"
344 if not read_only:
345 self._create_schema()
347 @property
348 def has_spatial(self) -> bool:
349 return self.spatial_error is None
351 def __enter__(self) -> Workspace:
352 return self
354 def __exit__(self, *exc: Any) -> None:
355 self.close()
357 def close(self) -> None:
358 self.con.close()
360 def _require_spatial(self) -> None:
361 if self.spatial_error is not None:
362 raise WorkspaceError(
363 f"DuckDB spatial extension unavailable ({self.spatial_error}). Install it once "
364 "with network access: python -c \"import duckdb; duckdb.sql('INSTALL spatial')\""
365 )
367 def _create_schema(self) -> None:
368 from .lroc.cumindex import _MACROS
370 for sql in [*_MACROS, *_BA_DDL]:
371 self.con.execute(sql)
372 if self.has_spatial:
373 for sql in _PS_MACROS:
374 self.con.execute(sql)
375 if not self._scalar("SELECT count(*) FROM rejection_params"):
376 self.set_rejection_params()
377 stats = _IMAGE_STATS_SQL.format(params=_PARAMS_VIEW)
378 rejections = _REJECTIONS_SQL.format(params=_PARAMS_VIEW, stats="SELECT * FROM ba_image_stats")
379 macro = _REJECTIONS_SQL.format(
380 params=_PARAMS_MACRO, stats=_IMAGE_STATS_SQL.format(params=_PARAMS_MACRO)
381 )
382 for sql in (
383 f"CREATE OR REPLACE VIEW ba_image_stats AS {stats}",
384 f"CREATE OR REPLACE VIEW image_rejections AS {rejections}",
385 f"CREATE OR REPLACE MACRO image_rejections_at(p85, frac, px, min_matches) AS TABLE {macro}",
386 f"CREATE OR REPLACE VIEW ba_run_summary AS {_RUN_SUMMARY_SQL}",
387 ):
388 self.con.execute(sql)
390 def query(self, sql: str, params: Sequence[Any] | None = None) -> pd.DataFrame:
391 """Run SQL; returns the result as a DataFrame (empty for statements without one)."""
392 import pandas as pd
394 cur = self.con.execute(sql, params) if params is not None else self.con.execute(sql)
395 return cur.df() if cur.description else pd.DataFrame()
397 def tables(self) -> pd.DataFrame:
398 """``name``, ``type`` (table / view) and ``rows`` of the workspace's tables and views."""
399 import pandas as pd
401 rows = self.con.execute(
402 "SELECT table_name, CASE table_type WHEN 'VIEW' THEN 'view' ELSE 'table' END "
403 "FROM information_schema.tables WHERE table_schema = 'main' ORDER BY table_type, table_name"
404 ).fetchall()
405 out = []
406 for name, kind in rows:
407 n = self._scalar(f'SELECT count(*) FROM "{name}"')
408 out.append({"name": name, "type": kind, "rows": int(n)})
409 return pd.DataFrame(out, columns=["name", "type", "rows"])
411 def _scalar(self, sql: str, params: Sequence[Any] | None = None) -> Any:
412 row = (
413 self.con.execute(sql, params).fetchone()
414 if params is not None
415 else self.con.execute(sql).fetchone()
416 )
417 return None if row is None else row[0]
419 def has_table(self, name: str) -> bool:
420 return bool(
421 self._scalar(
422 "SELECT count(*) FROM information_schema.tables WHERE table_schema = 'main' AND table_name = ?",
423 [name],
424 )
425 )
427 def export(self, table: str, path: Path | str) -> Path:
428 """Write a table or view to Parquet (``.parquet``) or CSV (anything else).
430 GEOMETRY columns are written as WKT in CSV and as WKB (GeoParquet) in Parquet.
431 """
432 name = _ident(table)
433 if not self.has_table(name):
434 raise WorkspaceError(f"No table or view {name!r} in {self.path}")
435 path = Path(path)
436 path.parent.mkdir(parents=True, exist_ok=True)
437 parquet = path.suffix.lower() in (".parquet", ".pq")
438 cols = self.con.execute(f'DESCRIBE "{name}"').fetchall()
439 sel = ", ".join(
440 f'ST_AsText("{c}") AS "{c}"' if t.startswith("GEOMETRY") and not parquet else f'"{c}"'
441 for c, t, *_ in cols
442 )
443 fmt = "(FORMAT parquet)" if parquet else "(FORMAT csv, HEADER true)"
444 self.con.execute(f'COPY (SELECT {sel} FROM "{name}") TO {_sql_str(path)} {fmt}')
445 return path
447 def add_observations(self, df: pd.DataFrame, *, rings: Sequence[Any] | None = None) -> int:
448 """Insert or replace (by ``product_id``) observation rows in ``images``.
450 Parameters
451 ----------
452 df : DataFrame
453 [`lunar_toolkit.observation`][lunar_toolkit.observation] schema rows; missing columns are NULL.
454 rings : sequence, optional
455 Per row, [`lunar_toolkit.geometry.footprint_rings`][lunar_toolkit.geometry.footprint_rings] output
456 to use instead of the rows' own (``footprint_wkt``, else the ``min_lat`` ... ``lon_east`` box).
458 Returns
459 -------
460 int
461 Rows written. ``footprint`` is the footprint polygon in polar stereographic
462 metres ([`lunar_toolkit.geometry.footprint_geometry`][lunar_toolkit.geometry.footprint_geometry]
463 with [`ps_projection`][lunar_toolkit.geometry.ps_projection] of ``hemi``, the hemisphere of
464 the image centre; edges densified along great circles every
465 `GREAT_CIRCLE_STEP_DEG`), ``footprint_km2`` its area corrected for the projection's
466 scale, ``cap_lat`` / ``cap_lon`` the centre of its vertices, ``obs_key`` the
467 observation (NAC LE / RE halves share one). ``idx`` numbers rows in insertion
468 order.
469 """
470 import pandas as pd
472 from .geometry import footprint_geometry, footprint_rings, ps_projection
473 from .observation import OBSERVATION_COLUMNS, normalize_observations
474 from .stereo import _cap
476 self._require_spatial()
477 self._create_images()
478 if df is None or len(df) == 0: 478 ↛ 479line 478 didn't jump to line 479 because the condition on line 478 was never true
479 return 0
480 obs = normalize_observations(df.to_dict("records"))
481 if rings is None:
482 rings = [footprint_rings(r, "footprint_wkt") for r in obs.to_dict("records")]
483 extra: dict[str, list[Any]] = {k: [] for k in (*_IMAGE_EXTRA, "fp_wkb")}
484 wkts = []
485 for row, rr in zip(obs.to_dict("records"), rings):
486 pid = str(row["product_id"])
487 extra["obs_key"].append((parse_nac_id(pid) or (pid.upper(),))[0])
488 geom = None
489 cap = (math.nan, math.nan)
490 if rr is not None and len(rr[0][0]) >= 3: 490 ↛ 496line 490 didn't jump to line 496 because the condition on line 490 was always true
491 c, _ = _cap([poly[0] for poly in rr])
492 cap = (
493 math.degrees(math.asin(max(-1.0, min(1.0, c[2])))),
494 math.degrees(math.atan2(c[1], c[0])) % 360,
495 )
496 lat = row["center_lat"] if pd.notna(row["center_lat"]) else cap[0]
497 hemi = None if pd.isna(lat) else ("S" if lat < 0 else "N")
498 if hemi is not None and rr is not None and len(rr[0][0]) >= 3: 498 ↛ 502line 498 didn't jump to line 502 because the condition on line 498 was always true
499 geom = footprint_geometry(_great_circle_rings(rr), ps_projection(hemi))
500 if geom is not None and (geom.is_empty or geom.area <= 0): 500 ↛ 501line 500 didn't jump to line 501 because the condition on line 500 was never true
501 geom = None
502 extra["hemi"].append(hemi)
503 extra["cap_lat"].append(cap[0])
504 extra["cap_lon"].append(cap[1])
505 extra["fp_wkb"].append(None if geom is None else geom.wkb)
506 if geom is not None: 506 ↛ 512line 506 didn't jump to line 512 because the condition on line 506 was always true
507 cx, cy = geom.centroid.x, geom.centroid.y
508 k = 1.0 + (cx * cx + cy * cy) / (4.0 * _R_M * _R_M)
509 xmin, ymin, xmax, ymax = geom.bounds
510 vals = (geom.area / k**2 / 1e6, xmin, xmax, ymin, ymax)
511 else:
512 vals = (math.nan,) * 5
513 for key, v in zip(("footprint_km2", "xmin", "xmax", "ymin", "ymax"), vals):
514 extra[key].append(v)
515 wkt = row.get("footprint_wkt")
516 if (not isinstance(wkt, str) or not wkt.strip()) and rr is not None and len(rr[0][0]) >= 3:
517 wkt = _rings_wkt(rr)
518 wkts.append(wkt if isinstance(wkt, str) else None)
519 frame = obs[OBSERVATION_COLUMNS].copy()
520 frame["footprint_wkt"] = wkts
521 frame["product_id"] = frame["product_id"].astype(str)
522 for k, v in extra.items():
523 frame[k] = v
524 frame = frame.drop_duplicates("product_id", keep="last").reset_index(drop=True)
525 start = self._scalar("SELECT coalesce(max(idx) + 1, 0) FROM images")
526 frame.insert(0, "idx", np.arange(start, start + len(frame), dtype=np.int64))
527 cols = ["idx", *_OBS_TYPES, *_IMAGE_EXTRA]
528 sel = ", ".join(
529 f"CAST({c} AS {t}) AS {c}" for c, t in {"idx": "BIGINT", **_OBS_TYPES, **_IMAGE_EXTRA}.items()
530 )
531 self.con.register("_new_images", frame)
532 self.con.execute("BEGIN TRANSACTION")
533 try:
534 self.con.execute("DELETE FROM images WHERE product_id IN (SELECT product_id FROM _new_images)")
535 self.con.execute(
536 f"INSERT INTO images ({', '.join(cols)}, footprint, added_at) "
537 f"SELECT {sel}, CASE WHEN fp_wkb IS NULL THEN NULL ELSE ST_GeomFromWKB(fp_wkb) END, "
538 "current_localtimestamp() FROM _new_images"
539 )
540 self.con.execute("COMMIT")
541 except BaseException:
542 self.con.execute("ROLLBACK")
543 raise
544 finally:
545 self.con.unregister("_new_images")
546 return len(frame)
548 def add_index_rows(self, df: pd.DataFrame, dataset: str | None = "edr") -> int:
549 """[`add_observations`][lunar_toolkit.workspace.Workspace.add_observations] of ``CumIndex`` rows
550 (footprints from ``FP_WKT``, corners or ``FP_*``)."""
551 from .geometry import footprint_rings
552 from .observation import observations_from_cumindex
554 if df is None or len(df) == 0: 554 ↛ 555line 554 didn't jump to line 555 because the condition on line 554 was never true
555 self._require_spatial()
556 self._create_images()
557 return 0
558 rings = [footprint_rings(r, "FP_WKT") for r in df.to_dict("records")]
559 return self.add_observations(observations_from_cumindex(df, dataset=dataset), rings=rings)
561 def _create_images(self) -> None:
562 cols = ", ".join(f"{c} {t}" for c, t in {"idx": "BIGINT", **_OBS_TYPES, **_IMAGE_EXTRA}.items())
563 self.con.execute(
564 f"CREATE TABLE IF NOT EXISTS images ({cols}, footprint GEOMETRY, added_at TIMESTAMP)"
565 )
567 def pairs(
568 self,
569 aoi: Any = None,
570 *,
571 convergence: tuple[float, float] = (10.0, 60.0),
572 ideal_convergence: tuple[float, float] = (15.0, 45.0),
573 max_incidence_diff: float = 20.0,
574 max_azimuth_diff: float = 30.0,
575 min_overlap: float = 0.1,
576 min_aoi_coverage: float = 0.0,
577 max_resolution_ratio: float | None = None,
578 same_instrument: bool = True,
579 require_geometry: bool = False,
580 min_time_sep_s: float = 60.0,
581 weights: dict[str, float] | None = None,
582 product_ids: Sequence[str] | None = None,
583 table: str | None = "pairs",
584 ) -> pd.DataFrame:
585 """Rank stereo pairs of ``images`` with a SQL spatial self-join.
587 Parameters are those of
588 [`lunar_toolkit.stereo.find_stereo_pairs`][lunar_toolkit.stereo.find_stereo_pairs], plus *product_ids*
589 (only these images) and *table* (where the result is stored; None keeps it out of the database).
591 Returns
592 -------
593 DataFrame
594 `lunar_toolkit.stereo.PAIR_COLUMNS`; ``index_1`` / ``index_2`` are
595 ``images.idx``.
597 Notes
598 -----
599 Candidate pairs come from a bounding-box spatial join (``ST_Intersects_Extent``,
600 R-tree built by DuckDB's spatial join operator) of the footprints of each
601 hemisphere, clipped to the AOI; then the attribute filters, overlap area / fraction (``ST_Area`` of
602 ``ST_Intersection`` in metres, divided by the projection's areal scale at the
603 centroid), AOI coverage, incidence / azimuth / emission / resolution / time
604 differences and the overlap centroid are SQL. The convergence angle and
605 ``b_h`` are then computed for the candidates with numpy
606 ([`lunar_toolkit.stereo.convergence_angles`][lunar_toolkit.stereo.convergence_angles]: sub-spacecraft
607 point + altitude, pushbroom-corrected), exactly as ``find_stereo_pairs``, and the
608 score is [`lunar_toolkit.stereo.score_pairs`][lunar_toolkit.stereo.score_pairs]. Footprint edges are
609 great circles in both (densified here), so overlap areas and fractions agree with
610 ``find_stereo_pairs`` to ~1e-3 relative.
611 """
612 import pandas as pd
614 from .geometry import _polar_hemi, aoi_geometry, aoi_region, load_aoi, ps_projection
615 from .stereo import PAIR_COLUMNS, convergence_angles, score_pairs
617 self._require_spatial()
618 self._create_images()
619 lo_c, hi_c = convergence
620 aoi_parts: dict[str, tuple[bytes, float]] = {}
621 if aoi is not None:
622 a = load_aoi(aoi)
623 region = aoi_region(a)
624 hemi = _polar_hemi(region) or ("S" if sum(region["lat"]) < 0 else "N")
625 geom = aoi_geometry(a, ps_projection(hemi))
626 k = 1.0 + (geom.centroid.x**2 + geom.centroid.y**2) / (4.0 * _R_M * _R_M)
627 aoi_parts[hemi] = (geom.wkb, geom.area / k**2 / 1e6)
628 hemis = [hemi]
629 else:
630 hemis = [
631 h
632 for (h,) in self.con.execute(
633 "SELECT DISTINCT hemi FROM images WHERE hemi IS NOT NULL"
634 ).fetchall()
635 ]
637 conds = ["a.obs_key <> b.obs_key"]
638 if same_instrument: 638 ↛ 640line 638 didn't jump to line 640 because the condition on line 638 was always true
639 conds.append("a.instrument IS NOT DISTINCT FROM b.instrument")
640 conds += [
641 f"NOT coalesce(abs(b.t - a.t) < {float(min_time_sep_s)!r}, false)",
642 f"NOT coalesce(abs(b.incidence - a.incidence) > {float(max_incidence_diff)!r}, false)",
643 f"NOT coalesce(_az_diff(a.az, b.az) > {float(max_azimuth_diff)!r}, false)",
644 f"NOT coalesce(a.emission + b.emission < {lo_c - 5.0!r}, false)",
645 f"NOT coalesce(abs(b.emission - a.emission) > {hi_c + 5.0!r}, false)",
646 ]
647 if require_geometry:
648 conds.append("a.has_geom AND b.has_geom")
649 if max_resolution_ratio is not None:
650 conds.append(
651 f"NOT coalesce({_ratio_sql('a.ps', 'b.ps')} > {float(max_resolution_ratio)!r}, false)"
652 )
654 frames = []
655 for hemi in hemis:
656 wkb, aoi_km2 = aoi_parts.get(hemi, (None, math.nan))
657 where = [f"hemi = {_sql_str(hemi)}", "footprint IS NOT NULL"]
658 params: list[Any] = []
659 if product_ids is not None:
660 self.con.register("_pair_ids", pd.DataFrame({"product_id": [str(p) for p in product_ids]}))
661 where.append("product_id IN (SELECT product_id FROM _pair_ids)")
662 g = "footprint"
663 if wkb is not None:
664 where.append("ST_Intersects(footprint, ST_GeomFromWKB(?))")
665 g = "ST_Intersection(footprint, ST_GeomFromWKB(?))"
666 params = [wkb, wkb]
667 self.con.execute(
668 "CREATE OR REPLACE TEMP TABLE _pair_fp AS SELECT *, _true_km2(g) AS g_km2 FROM ("
669 "SELECT idx, product_id, obs_key, instrument, start_time, epoch(start_time) AS t, incidence, "
670 "emission, sub_solar_azimuth AS az, pixel_scale AS ps, sub_spacecraft_lat AS sclat, "
671 "sub_spacecraft_lon AS sclon, spacecraft_altitude AS alt, "
672 "(sub_spacecraft_lat IS NOT NULL AND sub_spacecraft_lon IS NOT NULL "
673 "AND spacecraft_altitude IS NOT NULL) AS has_geom, "
674 "CASE WHEN center_lat IS NOT NULL AND center_lon IS NOT NULL THEN center_lat ELSE cap_lat END "
675 "AS vlat, CASE WHEN center_lat IS NOT NULL AND center_lon IS NOT NULL THEN center_lon "
676 f"ELSE cap_lon END AS vlon, {g} AS g FROM images WHERE {' AND '.join(where)}"
677 ") WHERE ST_Area(g) > 0",
678 params or None,
679 )
680 # Bounding-box prefilter: DuckDB plans its R-tree spatial join only for a bare
681 # ST_Intersects[_Extent] (+ inequality) condition, so candidates go to a temp table.
682 self.con.execute(
683 "CREATE OR REPLACE TEMP TABLE _pair_cand AS SELECT a.idx AS i, b.idx AS j "
684 "FROM _pair_fp a JOIN _pair_fp b ON ST_Intersects_Extent(a.g, b.g) AND a.idx < b.idx"
685 )
686 aoi_sql = f"ov_km2 / {aoi_km2!r}" if wkb is not None else "NULL::DOUBLE"
687 cov_filter = f"AND aoi_coverage >= {max(min_aoi_coverage, 1e-9)!r}" if wkb is not None else ""
688 self.con.execute(
689 "CREATE OR REPLACE TEMP TABLE _pair_ab AS SELECT a, b FROM _pair_cand c "
690 "JOIN (SELECT * EXCLUDE (g) FROM _pair_fp) a ON a.idx = c.i "
691 "JOIN (SELECT * EXCLUDE (g) FROM _pair_fp) b ON b.idx = c.j"
692 )
693 # Each step is materialised: otherwise the a / b filters become join conditions
694 # and ST_Intersection runs before them, once per reference to its result.
695 self.con.execute(f"""
696 CREATE OR REPLACE TEMP TABLE _pair_ov AS
697 SELECT a.idx AS index_1, b.idx AS index_2, a.product_id AS product_id_1,
698 b.product_id AS product_id_2, a.instrument,
699 a.start_time AS start_time_1, b.start_time AS start_time_2,
700 abs(b.t - a.t) / 86400.0 AS dt_days,
701 abs(b.emission - a.emission) AS convergence_min,
702 a.emission + b.emission AS convergence_max,
703 a.incidence AS incidence_1, b.incidence AS incidence_2,
704 abs(b.incidence - a.incidence) AS incidence_diff,
705 a.az AS azimuth_1, b.az AS azimuth_2, _az_diff(a.az, b.az) AS azimuth_diff,
706 a.emission AS emission_1, b.emission AS emission_2,
707 a.ps AS pixel_scale_1, b.ps AS pixel_scale_2,
708 {_ratio_sql("a.ps", "b.ps")} AS resolution_ratio,
709 a.g_km2 AS g1_km2, b.g_km2 AS g2_km2,
710 a.has_geom AND b.has_geom AS has_geom,
711 a.sclat AS sclat_1, a.sclon AS sclon_1, a.alt AS alt_1,
712 b.sclat AS sclat_2, b.sclon AS sclon_2, b.alt AS alt_2,
713 a.vlat AS vlat_1, a.vlon AS vlon_1, b.vlat AS vlat_2, b.vlon AS vlon_2
714 FROM _pair_ab
715 WHERE {" AND ".join(conds)}
716 """)
717 self.con.execute(
718 "CREATE OR REPLACE TEMP TABLE _pair_x AS SELECT o.*, ST_Intersection(a.g, b.g) AS ov "
719 "FROM _pair_ov o JOIN _pair_fp a ON a.idx = o.index_1 JOIN _pair_fp b ON b.idx = o.index_2"
720 )
721 self.con.execute(
722 "CREATE OR REPLACE TEMP TABLE _pair_m AS SELECT * EXCLUDE (ov, c), ST_X(c) AS cx, "
723 "ST_Y(c) AS cy FROM (SELECT *, ST_Area(ov) AS ov_m2, ST_Centroid(ov) AS c FROM _pair_x "
724 "WHERE NOT ST_IsEmpty(ov))"
725 )
726 df = self.con.execute(f"""
727 SELECT * FROM (
728 SELECT *, ov_km2 / greatest(least(g1_km2, g2_km2), 1e-12) AS overlap_frac,
729 {aoi_sql} AS aoi_coverage
730 FROM (SELECT * EXCLUDE (ov_m2), ov_m2 / pow(_ps_k(cx, cy), 2) / 1e6 AS ov_km2 FROM _pair_m)
731 WHERE ov_km2 > 0
732 )
733 WHERE overlap_frac >= {float(min_overlap)!r} {cov_filter}
734 ORDER BY index_1, index_2
735 """).df()
736 df["overlap_lat"] = _ps_lat(df["cx"].to_numpy(float), df["cy"].to_numpy(float), hemi)
737 df["overlap_lon"] = _ps_lon(df["cx"].to_numpy(float), df["cy"].to_numpy(float), hemi)
738 frames.append(df)
739 for tmp in ("_pair_m", "_pair_x", "_pair_ov", "_pair_ab", "_pair_cand", "_pair_fp"):
740 self.con.execute(f"DROP TABLE IF EXISTS {tmp}")
741 if product_ids is not None:
742 self.con.unregister("_pair_ids")
744 frames = [f for f in frames if len(f)]
745 if not frames:
746 res = pd.DataFrame(columns=PAIR_COLUMNS)
747 else:
748 res = pd.concat(frames, ignore_index=True).sort_values(["index_1", "index_2"], kind="stable")
749 res = res.reset_index(drop=True)
750 geom = res["has_geom"].fillna(False).to_numpy(bool)
751 conv = np.full(len(res), np.nan)
752 bh = np.full(len(res), np.nan)
753 if geom.any():
754 r = res[geom]
756 def col(c: str) -> np.ndarray:
757 return r[c].to_numpy(float)
759 conv[geom], bh[geom] = convergence_angles(
760 (col("overlap_lat"), col("overlap_lon")),
761 (col("sclat_1"), col("sclon_1"), col("alt_1")),
762 (col("sclat_2"), col("sclon_2"), col("alt_2")),
763 (col("vlat_1"), col("vlon_1")),
764 (col("vlat_2"), col("vlon_2")),
765 )
766 res["convergence"], res["b_h"] = conv, bh
767 res["overlap_km2"] = res["ov_km2"]
768 keep = np.isnan(conv) | ((conv >= lo_c) & (conv <= hi_c))
769 res = res[keep].reset_index(drop=True)
770 for c in ("start_time_1", "start_time_2"):
771 res[c] = pd.to_datetime(res[c])
772 res["score"] = score_pairs(
773 res,
774 convergence=convergence,
775 ideal_convergence=ideal_convergence,
776 max_incidence_diff=max_incidence_diff,
777 max_azimuth_diff=max_azimuth_diff,
778 weights=weights,
779 aoi=aoi is not None,
780 )
781 res = (
782 res[PAIR_COLUMNS].sort_values("score", ascending=False, kind="stable").reset_index(drop=True)
783 )
784 if table is not None:
785 name = _ident(table)
786 self.con.register("_pairs_result", res)
787 try:
788 self.con.execute(f'CREATE OR REPLACE TABLE "{name}" AS SELECT * FROM _pairs_result')
789 finally:
790 self.con.unregister("_pairs_result")
791 return res
793 def set_rejection_params(self, config: SfsConfig | None = None, **overrides: Any) -> None:
794 """Thresholds of the ``image_rejections`` view, from *config* (default
795 [`SfsConfig`][lunar_toolkit.sfs.SfsConfig]).
797 Keyword overrides: ``max_offset_p85``, ``max_large_residual_fraction``,
798 ``large_residual_px``, ``min_matches`` (None disables the first two rules).
799 """
800 from .sfs import SfsConfig
802 cfg = config or SfsConfig()
803 vals = {
804 "max_offset_p85": cfg.max_offset_p85,
805 "max_large_residual_fraction": cfg.max_large_residual_fraction,
806 "large_residual_px": cfg.large_residual_px,
807 "min_matches": cfg.min_matches,
808 }
809 unknown = set(overrides) - set(vals)
810 if unknown:
811 raise TypeError(f"Unknown rejection parameter(s): {', '.join(sorted(unknown))}")
812 vals.update(overrides)
813 self.con.execute("DELETE FROM rejection_params")
814 self.con.execute(
815 "INSERT INTO rejection_params VALUES (?, ?, ?, ?)",
816 [
817 vals["max_offset_p85"],
818 vals["max_large_residual_fraction"],
819 vals["large_residual_px"],
820 int(vals["min_matches"] or 0),
821 ],
822 )
824 def load_bundle_adjust(self, prefix: Path | str, run: str | None = None) -> dict[str, int]:
825 """Load ASP ``bundle_adjust`` outputs ``<prefix>-*`` as run *run* (replacing it).
827 *run* defaults to the prefix's directory name (``ba``, ``ba_refine``). Reads,
828 when present: ``-initial`` / ``-final_residuals_stats.txt`` and
829 ``_raw_pixels.txt`` (``ba_residual_stats``, ``ba_pixel_residuals``),
830 ``-mapproj_match_offset_stats.txt`` (``ba_match_offsets``; -1 becomes NULL),
831 ``-initial`` / ``-final_residuals_pointmap.csv`` (``ba_points``, with polar
832 stereographic ``x`` / ``y``), ``-convergence_angles.txt`` (``ba_convergence``)
833 and ``-camera_offsets.txt`` (``ba_camera_offsets``).
835 Returns
836 -------
837 dict
838 Rows loaded per table (empty when no report exists).
839 """
840 from . import asp
842 prefix = Path(prefix)
843 run = run or prefix.parent.name or prefix.name
844 f = {name: Path(f"{prefix}-{name}") for name in _BA_FILES}
845 found = [p for p in f.values() if p.exists()]
846 if not found:
847 return {}
848 stats = {
849 stage: asp.read_residual_stats(f[f"{stage}_residuals_stats.txt"])
850 for stage in ("initial", "final")
851 if f[f"{stage}_residuals_stats.txt"].exists()
852 }
853 raw = {
854 stage: _read_raw_pixels(f[f"{stage}_residuals_raw_pixels.txt"])
855 for stage in ("initial", "final")
856 if f[f"{stage}_residuals_raw_pixels.txt"].exists()
857 }
858 offsets = (
859 asp.read_mapproj_offset_stats(f["mapproj_match_offset_stats.txt"])
860 if f["mapproj_match_offset_stats.txt"].exists()
861 else None
862 )
863 self.con.execute("BEGIN TRANSACTION")
864 try:
865 self._insert_ba(run, stats=stats, raw=raw, offsets=offsets)
866 for stage in ("initial", "final"):
867 p = f[f"{stage}_residuals_pointmap.csv"]
868 if p.exists():
869 self._load_points(run, stage, p)
870 if f["convergence_angles.txt"].exists():
871 rows = _text_rows(self.con, f["convergence_angles.txt"], " ", 6)
872 self.con.execute(
873 "INSERT INTO ba_convergence SELECT ?, parse_filename(c1), parse_filename(c2), "
874 "upper(split_part(parse_filename(c1), '.', 1)), upper(split_part(parse_filename(c2), '.', 1)), "
875 f"TRY_CAST(c3 AS DOUBLE), TRY_CAST(c4 AS DOUBLE), TRY_CAST(c5 AS DOUBLE), "
876 f"CAST(TRY_CAST(c6 AS DOUBLE) AS BIGINT) FROM ({rows}) WHERE c6 IS NOT NULL",
877 [run],
878 )
879 if f["camera_offsets.txt"].exists():
880 rows = _text_rows(self.con, f["camera_offsets.txt"], " ", 3)
881 self.con.execute(
882 "INSERT INTO ba_camera_offsets SELECT ?, parse_filename(c1), "
883 "upper(split_part(parse_filename(c1), '.', 1)), TRY_CAST(c2 AS DOUBLE), "
884 f"TRY_CAST(c3 AS DOUBLE) FROM ({rows}) WHERE c3 IS NOT NULL",
885 [run],
886 )
887 self.con.execute("DELETE FROM ba_runs WHERE run = ?", [run])
888 self.con.execute(
889 "INSERT INTO ba_runs VALUES (?, ?, current_localtimestamp(), ?)",
890 [run, str(prefix), [p.name for p in found]],
891 )
892 self.con.execute("COMMIT")
893 except BaseException:
894 self.con.execute("ROLLBACK")
895 raise
896 return {t: int(self._scalar(f"SELECT count(*) FROM {t} WHERE run = ?", [run])) for t in BA_TABLES}
898 def _insert_ba(
899 self,
900 run: str,
901 *,
902 stats: Mapping[str, Mapping[str, ResidualStats]] | None = None,
903 raw: Mapping[str, pd.DataFrame] | None = None,
904 offsets: Mapping[str, MatchOffsetStats] | None = None,
905 ) -> None:
906 """Replace run *run* with parsed reports: per stage residual stats / raw norms, offsets."""
907 import pandas as pd
909 for t in BA_TABLES:
910 self.con.execute(f"DELETE FROM {t} WHERE run = ?", [run])
911 rows = [
912 (
913 run,
914 stage,
915 Path(s.image).name,
916 s.image,
917 _product_id(s.image),
918 _num(s.mean),
919 _num(s.median),
920 s.count,
921 )
922 for stage, by_image in (stats or {}).items()
923 for s in by_image.values()
924 ]
925 if rows:
926 self.con.executemany("INSERT INTO ba_residual_stats VALUES (?, ?, ?, ?, ?, ?, ?, ?)", rows)
927 for stage, norms in (raw or {}).items():
928 if len(norms): 928 ↛ 927line 928 didn't jump to line 927 because the condition on line 928 was always true
929 frame = pd.DataFrame(
930 {
931 "run": run,
932 "stage": stage,
933 "image": norms["image"].astype(str),
934 "product_id": [_product_id(i) for i in norms["image"]],
935 "norm_px": norms["norm_px"].astype(float),
936 }
937 )
938 self.con.register("_raw", frame)
939 try:
940 self.con.execute("INSERT INTO ba_pixel_residuals SELECT * FROM _raw")
941 finally:
942 self.con.unregister("_raw")
943 orows = [
944 (
945 run,
946 Path(o.image).name,
947 o.image,
948 _product_id(o.image),
949 *(None if not v >= 0 else v for v in (o.p25, o.p50, o.p75, o.p85, o.p95)),
950 o.count,
951 )
952 for o in (offsets or {}).values()
953 ]
954 if orows:
955 self.con.executemany("INSERT INTO ba_match_offsets VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?)", orows)
957 def _load_points(self, run: str, stage: str, path: Path) -> None:
958 rows = _text_rows(self.con, path, ",", 5)
959 self.con.execute(
960 f"""
961 INSERT INTO ba_points
962 SELECT run, stage, lon, lat, height_m, res, nobs, kind, hemi,
963 _ps_x(lat, lon, hemi), _ps_y(lat, lon, hemi)
964 FROM (
965 SELECT ? AS run, ? AS stage, TRY_CAST(c1 AS DOUBLE) AS lon, TRY_CAST(c2 AS DOUBLE) AS lat,
966 TRY_CAST(c3 AS DOUBLE) AS height_m, TRY_CAST(c4 AS DOUBLE) AS res,
967 TRY_CAST(split_part(trim(c5), ' ', 1) AS INTEGER) AS nobs,
968 CASE WHEN c5 LIKE '%GCP%' THEN 'gcp' WHEN c5 LIKE '%DEM%' THEN 'dem' ELSE 'tie' END AS kind,
969 CASE WHEN TRY_CAST(c2 AS DOUBLE) < 0 THEN 'S' ELSE 'N' END AS hemi
970 FROM ({rows}) WHERE c5 IS NOT NULL
971 ) WHERE lon IS NOT NULL AND lat IS NOT NULL
972 """,
973 [run, stage],
974 )
976 def runs(self) -> list[str]:
977 """Loaded bundle-adjustment runs."""
978 return [r for (r,) in self.con.execute("SELECT run FROM ba_runs ORDER BY loaded_at, run").fetchall()]
980 def rejections(self, run: str) -> dict[str, str]:
981 """``image_rejections`` of *run*: image file name -> reason."""
982 rows = self.con.execute("SELECT image, reason FROM image_rejections WHERE run = ?", [run]).fetchall()
983 return dict(rows)
985 def ba_report(self, run: str) -> pd.DataFrame | None:
986 """Per-image table of *run* in
987 [`lunar_toolkit.sfs_report.read_ba_reports`][lunar_toolkit.sfs_report.read_ba_reports] form, or
988 None."""
989 if run not in self.runs(): 989 ↛ 990line 989 didn't jump to line 990 because the condition on line 989 was never true
990 return None
991 cols = ["image"] + [f"{k}_{s}" for k in ("initial", "final") for s in ("mean", "median", "count")]
992 cols += [f"offset_p{p}" for p in (25, 50, 75, 85, 95)] + ["offset_count"]
993 df = self.con.execute(
994 f"SELECT {', '.join(cols)} FROM ba_image_stats WHERE run = ? AND "
995 "(initial_count IS NOT NULL OR final_count IS NOT NULL OR offset_count IS NOT NULL) "
996 "ORDER BY image",
997 [run],
998 ).df()
999 for c in cols[1:]:
1000 df[c] = df[c].astype(float)
1001 return df
1003 def convergence(self, run: str) -> pd.DataFrame:
1004 """``ba_convergence`` of *run* as ``left``, ``right``, ``p25``, ``p50``, ``p75``, ``matches``."""
1005 return self.con.execute(
1006 "SELECT left_image AS left, right_image AS right, p25, p50, p75, matches "
1007 "FROM ba_convergence WHERE run = ?",
1008 [run],
1009 ).df()
1011 def point_residuals(self, run: str, stage: str = "final") -> np.ndarray:
1012 """``mean_residual_px`` of the *run* / *stage* residual point map."""
1013 rows = self.con.execute(
1014 "SELECT mean_residual_px FROM ba_points WHERE run = ? AND stage = ?", [run, stage]
1015 ).fetchnumpy()
1016 return np.asarray(rows["mean_residual_px"], float)
1018 def summary(self) -> pd.DataFrame:
1019 """The ``ba_run_summary`` view."""
1020 return self.con.execute("SELECT * FROM ba_run_summary").df()
1023_BA_FILES = (
1024 "initial_residuals_stats.txt",
1025 "final_residuals_stats.txt",
1026 "initial_residuals_raw_pixels.txt",
1027 "final_residuals_raw_pixels.txt",
1028 "mapproj_match_offset_stats.txt",
1029 "initial_residuals_pointmap.csv",
1030 "final_residuals_pointmap.csv",
1031 "convergence_angles.txt",
1032 "camera_offsets.txt",
1033)
1036def _great_circle_rings(rings: list, step_deg: float = GREAT_CIRCLE_STEP_DEG) -> list:
1037 """*rings* with vertices inserted along great circles, at most *step_deg* apart."""
1038 from .stereo import unit_vector
1040 out = []
1041 for poly in rings:
1042 dense = []
1043 for ring in poly:
1044 r = list(ring[:-1] if len(ring) > 1 and tuple(ring[0]) == tuple(ring[-1]) else ring)
1045 if len(r) < 3: 1045 ↛ 1046line 1045 didn't jump to line 1046 because the condition on line 1045 was never true
1046 dense.append(r)
1047 continue
1048 a = np.asarray(r, float)
1049 v = unit_vector(a[:, 1], a[:, 0])
1050 w = np.roll(v, -1, axis=0)
1051 theta = np.arctan2(np.linalg.norm(np.cross(v, w), axis=1), np.sum(v * w, axis=1))
1052 pts: list[tuple[float, float]] = []
1053 for k, (vk, wk, th) in enumerate(zip(v, w, theta)):
1054 pts.append((float(a[k, 0]), float(a[k, 1])))
1055 n = int(math.ceil(math.degrees(th) / step_deg))
1056 if n < 2 or th < 1e-12:
1057 continue
1058 f = np.arange(1, n)[:, None] / n
1059 q = (np.sin((1 - f) * th) * vk + np.sin(f * th) * wk) / math.sin(th)
1060 lat = np.degrees(np.arcsin(np.clip(q[:, 2], -1, 1)))
1061 lon = np.degrees(np.arctan2(q[:, 1], q[:, 0]))
1062 # keep the ring's longitude unwrapping (values may lie outside 0..360)
1063 lon = a[k, 0] + (lon - a[k, 0] + 180.0) % 360.0 - 180.0
1064 pts += list(zip(lon.tolist(), lat.tolist()))
1065 dense.append(pts)
1066 out.append(dense)
1067 return out
1070def _ratio_sql(a: str, b: str) -> str:
1071 return f"(CASE WHEN {a} >= {b} THEN {a} / {b} ELSE {b} / {a} END)"
1074def _ps_lat(x: np.ndarray, y: np.ndarray, hemi: str) -> np.ndarray:
1075 sgn = 1.0 if hemi == "S" else -1.0
1076 return sgn * (np.degrees(2.0 * np.arctan(np.hypot(x, y) / (2.0 * _R_M))) - 90.0)
1079def _ps_lon(x: np.ndarray, y: np.ndarray, hemi: str) -> np.ndarray:
1080 return np.degrees(np.arctan2(x, y if hemi == "S" else -y)) % 360.0
1083def _rings_wkt(rings: list) -> str:
1084 from shapely.geometry import MultiPolygon, Polygon
1086 polys = [Polygon(r[0], [h for h in r[1:] if len(h) >= 3]) for r in rings if len(r[0]) >= 3]
1087 return (polys[0] if len(polys) == 1 else MultiPolygon(polys)).wkt
1090def _text_rows(con: Any, path: Path, delim: str, n: int) -> str:
1091 """SQL reading *path* as *n* VARCHAR columns ``c1``.. (``#`` lines dropped, short rows NULL-padded)."""
1092 cols = "{" + ", ".join(f"'c{k}': 'VARCHAR'" for k in range(1, n + 1)) + "}"
1093 return (
1094 f"SELECT * FROM read_csv({_sql_str(path)}, delim={_sql_str(delim)}, header=false, columns={cols}, "
1095 "null_padding=true, quote='', escape='', auto_detect=false, strict_mode=false) "
1096 "WHERE NOT starts_with(trim(c1), '#')"
1097 )
1100def _read_raw_pixels(path: Path) -> pd.DataFrame:
1101 """``image``, ``norm_px`` of a ``*_residuals_raw_pixels.txt`` (``name, count`` then ``ex, ey`` lines)."""
1102 import pandas as pd
1104 df = pd.read_csv(path, header=None, names=["a", "b"], dtype=str, skipinitialspace=True, comment=None)
1105 a = pd.to_numeric(df["a"], errors="coerce")
1106 head = a.isna().to_numpy()
1107 if not head.any(): 1107 ↛ 1108line 1107 didn't jump to line 1108 because the condition on line 1107 was never true
1108 return pd.DataFrame({"image": pd.Series(dtype=str), "norm_px": pd.Series(dtype=float)})
1109 names = df["a"].where(head).ffill()
1110 keep = ~head & names.notna().to_numpy()
1111 ex, ey = a[keep].to_numpy(float), pd.to_numeric(df["b"][keep], errors="coerce").to_numpy(float)
1112 image = names[keep].map(lambda s: Path(str(s).strip()).name)
1113 return pd.DataFrame({"image": image.to_numpy(), "norm_px": np.hypot(ex, ey)})
1116def find_stereo_pairs_duckdb(
1117 obs_df: pd.DataFrame, aoi: Any = None, *, workspace: Workspace | None = None, **kwargs: Any
1118) -> pd.DataFrame:
1119 """[`lunar_toolkit.stereo.find_stereo_pairs`][lunar_toolkit.stereo.find_stereo_pairs] through
1120 [`Workspace.pairs`][lunar_toolkit.workspace.Workspace.pairs].
1122 The observations are written to ``images`` of *workspace* (default: a scratch
1123 in-memory one) and ``index_1`` / ``index_2`` are mapped back to *obs_df* row
1124 labels; keyword arguments are those of ``find_stereo_pairs``.
1125 """
1126 import pandas as pd
1128 from .stereo import PAIR_COLUMNS
1130 if obs_df is None or len(obs_df) < 2: 1130 ↛ 1131line 1130 didn't jump to line 1131 because the condition on line 1130 was never true
1131 return pd.DataFrame(columns=PAIR_COLUMNS)
1132 pids = obs_df["product_id"].astype(str)
1133 if pids.duplicated().any(): 1133 ↛ 1134line 1133 didn't jump to line 1134 because the condition on line 1133 was never true
1134 raise ValueError("find_stereo_pairs_duckdb needs unique product_id values")
1135 label = dict(zip(pids, obs_df.index))
1136 ws = workspace or Workspace(":memory:")
1137 try:
1138 ws.add_observations(obs_df)
1139 res = ws.pairs(
1140 aoi, product_ids=list(pids), table="pairs" if workspace is not None else None, **kwargs
1141 )
1142 finally:
1143 if workspace is None:
1144 ws.close()
1145 res["index_1"] = [label[p] for p in res["product_id_1"]]
1146 res["index_2"] = [label[p] for p in res["product_id_2"]]
1147 return res
1150def flag_images_sql(
1151 residual_stats: Mapping[str, ResidualStats],
1152 residual_norms: Mapping[str, np.ndarray],
1153 offsets: Mapping[str, MatchOffsetStats],
1154 **thresholds: Any,
1155) -> dict[str, str]:
1156 """``image_rejections`` of parsed final-stage reports (see
1157 [`lunar_toolkit.sfs.flag_images`][lunar_toolkit.sfs.flag_images])."""
1158 import pandas as pd
1160 raw = pd.DataFrame(
1161 {
1162 "image": [n for n, v in residual_norms.items() for _ in range(len(v))],
1163 "norm_px": np.concatenate([np.asarray(v, float) for v in residual_norms.values()])
1164 if residual_norms
1165 else np.empty(0),
1166 }
1167 )
1168 with Workspace(":memory:", spatial=False, threads=1) as ws:
1169 ws.set_rejection_params(**thresholds)
1170 stats = {k: _with_image(v, k) for k, v in residual_stats.items()}
1171 offs = {k: _with_image(v, k) for k, v in offsets.items()}
1172 ws._insert_ba("_", stats={"final": stats}, raw={"final": raw}, offsets=offs)
1173 return ws.rejections("_")
1176def _with_image(obj: Any, name: str) -> Any:
1177 """*obj* (a frozen report row) keyed by *name*: its ``image`` replaced when the names differ."""
1178 import dataclasses
1180 return obj if Path(obj.image).name == name else dataclasses.replace(obj, image=name)