Coverage for lunar_toolkit/sfs.py: 91%
719 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 02:58 +0000
1"""Shape-from-shading DEMs from LRO NAC images with ISIS and ASP (Bertone et al. 2026 workflow)."""
3from __future__ import annotations
5import csv
6import dataclasses
7import importlib
8import logging
9import math
10import re
11import shutil
12from collections.abc import Sequence
13from dataclasses import dataclass, field
14from pathlib import Path
15from typing import TYPE_CHECKING, Any
17import numpy as np
19from . import asp, isis
20from .constants import (
21 CAMERA_SOURCES,
22 MOON_RADIUS_KM,
23 REGISTRATION_MODES,
24 SELECTION_METHODS,
25 SFS_CAMERA_MODELS,
26)
27from .geometry import load_aoi, ps_xy
28from .selection import select_kmedoids
29from .stereo import azimuth_difference
30from .stereo_dem import stereo_dem
31from .tools import ToolError, is_dry_run, reuse_existing
33if TYPE_CHECKING: # pragma: no cover
34 import pandas as pd
36 from .config import RunConfig
37 from .provenance import Provenance
38 from .workspace import Workspace
40log = logging.getLogger(__name__)
42# ASP "Bundle adjustment" recipe for many LRO NAC images under varied illumination
43# (docs/sfs_usage.rst; --match-first-to-last: near the poles the Sun loops the sky).
44ASP_BUNDLE_ADJUST: dict[str, Any] = {
45 "ip_per_image": 30000,
46 "overlap_limit": 100,
47 "num_iterations": 100,
48 "num_passes": 2,
49 "min_matches": 1,
50 "max_pairwise_matches": 2000,
51 "camera_weight": 0,
52 "robust_threshold": 2,
53 "tri_weight": 0.05,
54 "tri_robust_threshold": 0.05,
55 "remove_outliers_params": "75.0 3.0 100 100",
56 "save_intermediate_cameras": True,
57 "match_first_to_last": True,
58 "forced_triangulation_distance": 100000,
59 "min_triangulation_angle": 1e-10,
60 "datum": "D_MOON",
61}
63# ASP "Registration refinement" recipe (bundle_adjust --heights-from-dem, reusing the matches).
64ASP_REFINE: dict[str, Any] = {
65 "max_pairwise_matches": 2000,
66 "match_first_to_last": True,
67 "min_matches": 1,
68 "num_iterations": 100,
69 "num_passes": 2,
70 "camera_weight": 0,
71 "save_intermediate_cameras": True,
72 "heights_from_dem_uncertainty": 20.0,
73 "forced_triangulation_distance": 100000,
74 "min_triangulation_angle": 1e-10,
75 "remove_outliers_params": "75.0 3.0 100 100",
76 "parameter_tolerance": 1e-20,
77 "datum": "D_MOON",
78}
80# ASP "Running SfS in parallel" recipe.
81ASP_SFS: dict[str, Any] = {
82 "tile_size": 200,
83 "padding": 50,
84 "shadow_threshold": 0.005,
85 "crop_input_images": True,
86 "blending_dist": 10,
87 "min_blend_size": 50,
88 "allow_borderline_data": True,
89 "low_light_threshold": 0.03,
90 "smoothness_weight": 0.08,
91 "initial_dem_constraint_weight": 0.0025,
92 "reflectance_type": 1,
93 "max_iterations": 5,
94 "save_sparingly": True,
95}
97# Bertone et al. (2026): match each image against the next 200 in Sun-azimuth order whose
98# footprints overlap (``auto_overlap_percent`` → ``--auto-overlap-params "<dem> 15"``), and
99# soft-anchor cameras at 100 m horizontal / vertical (App. A); intrinsics stay fixed
100# (ASP's default without --solve-intrinsics).
101_UNCERTAINTY = {"camera_position_uncertainty": "100,100"}
102BERTONE_BUNDLE_ADJUST: dict[str, Any] = {
103 **ASP_BUNDLE_ADJUST,
104 # The next 200 images in Sun-azimuth order (Bertone et al. 2026, §2). Footprint-based
105 # pairs (`auto_overlap_percent`) need cameras that fall mostly on the prior DEM, which
106 # a crop does not satisfy; set it as an override for full-region runs.
107 "overlap_limit": 200,
108 **_UNCERTAINTY,
109}
110BERTONE_REFINE: dict[str, Any] = {**ASP_REFINE, **_UNCERTAINTY}
111# §2: mu = 0.05–0.1, lambda = 1e-3, five iterations; App. A: shadow threshold 0.002, albedo fixed.
112BERTONE_SFS: dict[str, Any] = {
113 **ASP_SFS,
114 "shadow_threshold": 0.002,
115 "smoothness_weight": 0.08,
116 "initial_dem_constraint_weight": 1e-3,
117 "max_iterations": 5,
118}
120PRESETS: dict[str, dict[str, dict[str, Any]]] = {
121 "bertone": {"bundle_adjust": BERTONE_BUNDLE_ADJUST, "refine": BERTONE_REFINE, "sfs": BERTONE_SFS},
122 "asp": {"bundle_adjust": ASP_BUNDLE_ADJUST, "refine": ASP_REFINE, "sfs": ASP_SFS},
123}
126_OVERLAP_OPTIONS = ("overlap_list", "auto_overlap_params", "auto_overlap_buffer", "overlap_limit")
129@dataclass
130class SfsConfig:
131 """Settings for [`SfsPipeline`][lunar_toolkit.sfs.SfsPipeline].
133 Parameters
134 ----------
135 resolution : float
136 SfS DEM grid spacing (m/px).
137 projwin : tuple, optional
138 ``(ulx, uly, lrx, lry)`` crop of the prior DEM, in its projected units.
139 preset : {"bertone", "asp"}
140 Option defaults: Bertone et al. (2026) or the ASP SfS recipe (`PRESETS`).
141 ba_resolution : float, optional
142 Grid spacing of the images mapprojected for bundle-adjustment matching and
143 stereo registration (default 1 m, NAC native; None = *resolution*).
144 reduce : int, "auto" or None
145 ISIS ``reduce`` factor bringing the cubes to *resolution* before SfS.
146 ``"auto"``: ``round(resolution / native)`` per image (native = observation
147 ``pixel_scale``, else *native_resolution*). Reduced cubes use ISIS cameras
148 (with ``--use-approx-camera-models``); bundle adjustment still uses CSM.
149 native_resolution : float
150 Assumed image m/px when the observation table does not give it.
151 csm : bool
152 CSM cameras (7–15× faster than ISIS) for mapproject, bundle_adjust and stereo.
153 web_spice : bool
154 ``spiceinit web=true`` (USGS SPICE web service) instead of local kernels.
155 registration : {"stereo", "residuals", "none"}
156 Registering the adjusted cameras to the prior DEM (translation only).
157 ``"stereo"``: DEMs from up to *stereo_pairs* pairs, mosaicked and
158 ``pc_align``-ed (falls back to ``"residuals"`` without pairs);
159 ``"residuals"``: the bundle-adjustment triangulated points.
160 stereo_pairs : int
161 Maximum stereo pairs for ``"stereo"`` registration.
162 stereo_convergence : (min, max)
163 Convergence-angle limits (deg) for those pairs.
164 max_displacement : float
165 ``pc_align --max-displacement`` (m).
166 refine : bool
167 Re-run bundle adjustment on the kept, registered cameras with
168 ``--heights-from-dem`` (reusing the matches).
169 max_offset_p85 : float, optional
170 Drop images whose 85th-percentile mapprojected match offset (m) is larger.
171 Applied to the DEM-registered, refined cameras (ASP's ``ba_align_ref``
172 report); with *refine* off, to the unregistered ``ba/`` cameras.
173 max_large_residual_fraction : float, optional
174 Drop images with more than this fraction of reprojection errors over
175 *large_residual_px* pixels.
176 min_matches : int
177 Drop images with fewer inlier matches.
178 min_images : int
179 Stop with [`PipelineError`][lunar_toolkit.sfs.PipelineError] when a quality filter leaves fewer
180 images.
181 selection : {"coverage", "image_subset", "all"}
182 SfS image selection: greedy per-pixel illumination coverage
183 ([`select_by_illumination`][lunar_toolkit.sfs.select_by_illumination]), ASP ``image_subset`` per
184 Sun-azimuth bin, or all kept images.
185 max_shadow_fraction : float, optional
186 Images at least this shadowed (below the SfS shadow threshold) are dropped.
187 selection_views, selection_azimuth_sep : int, float
188 Illumination directions wanted per pixel and their minimum separation (deg).
189 max_sfs_images : int, optional
190 Cap on the selected images (``"coverage"`` / ``"kmedoids"`` selection).
191 selection_clusters : int
192 Medoids seeding ``"kmedoids"`` selection (Bertone et al. 2026, §2: ~30).
193 sfs_cameras : {"csm", "isis"}
194 Camera models for SfS on reduced cubes: CSM cameras re-mapped to the reduced
195 image ([`lunar_toolkit.csm.isis_reduce_with_csm`][lunar_toolkit.csm.isis_reduce_with_csm], as in
196 Bertone et al. 2026, App. A) or ISIS cameras with ``--use-approx-camera-models``.
197 cameras_from : {"ba", "published"}
198 ``"published"`` uses the bundle-adjusted cameras of Bertone et al. (2026) for
199 *published_region*
200 ([`lunar_toolkit.reference.published_cameras`][lunar_toolkit.reference.published_cameras]) and skips
201 bundle adjustment, registration and refinement; images without one are removed.
202 published_region : str or int, optional
203 Region of the published cameras (default: from ``region_id`` ``A3NN``).
204 max_sfs_slew : float, optional
205 Only images with ``|slew| <`` this (deg) are selected for SfS (Bertone et al.
206 2026, §2: nadir images); slewed images still join bundle adjustment and stereo
207 registration. Unknown slews count as nadir; ``None`` disables the cut.
208 selection_grid : int
209 Cells along the longer side of the selection grid.
210 region_id : str, optional
211 ``ROIID`` for product file names (default: the work directory name).
212 ensemble_members : int
213 SfS runs on random subsets of the kept images (sharing at most 40% of their
214 images) over an *ensemble_crop* m square at the DEM centre, for the ensemble
215 uncertainty maps (Bertone et al. 2026, §3.2); 0 skips this.
216 ensemble_size : int, optional
217 Images per ensemble member (default: as many as selected).
218 ensemble_crop : float
219 Side of the ensemble crop (m).
220 seed : int
221 Random seed of the ensemble image subsets.
222 processes, threads : int, optional
223 Passed to ``parallel_bundle_adjust``, ``parallel_stereo``, ``parallel_sfs`` and ``mapproject``.
224 bundle_adjust_options, refine_options, sfs_options : dict
225 Overrides merged onto the preset; a value of ``None`` drops a default.
226 """
228 resolution: float = 5.0
229 projwin: tuple[float, float, float, float] | None = None
230 preset: str = "bertone"
231 ba_resolution: float | None = 1.0
232 reduce: int | str | None = "auto"
233 native_resolution: float = 1.0
234 csm: bool = True
235 web_spice: bool = True
236 registration: str = "stereo"
237 stereo_pairs: int = 30
238 stereo_convergence: tuple[float, float] = (5.0, 60.0)
239 max_displacement: float = 400.0
240 refine: bool = True
241 max_offset_p85: float | None = 1.8
242 max_large_residual_fraction: float | None = 0.85
243 large_residual_px: float = 1.0
244 min_matches: int = 1
245 min_images: int = 3
246 selection: str = "coverage"
247 max_shadow_fraction: float | None = 0.95
248 selection_views: int = 3
249 selection_azimuth_sep: float = 30.0
250 max_sfs_images: int | None = None
251 selection_clusters: int = 30
252 sfs_cameras: str = "csm"
253 cameras_from: str = "ba"
254 published_region: str | int | None = None
255 max_sfs_slew: float | None = 5.0
256 selection_grid: int = 250
257 region_id: str | None = None
258 ensemble_members: int = 0
259 ensemble_size: int | None = None
260 ensemble_crop: float = 2500.0
261 seed: int = 0
262 processes: int | None = None
263 threads: int | None = None
264 bundle_adjust_options: dict[str, Any] = field(default_factory=dict)
265 refine_options: dict[str, Any] = field(default_factory=dict)
266 sfs_options: dict[str, Any] = field(default_factory=dict)
268 def __post_init__(self) -> None:
269 for name, value, allowed in (
270 ("preset", self.preset, tuple(PRESETS)),
271 ("registration", self.registration, REGISTRATION_MODES),
272 ("selection", self.selection, SELECTION_METHODS),
273 ("sfs_cameras", self.sfs_cameras, SFS_CAMERA_MODELS),
274 ("cameras_from", self.cameras_from, CAMERA_SOURCES),
275 ):
276 if value not in allowed:
277 raise ValueError(f"{name} must be one of {', '.join(allowed)}, not {value!r}")
278 if isinstance(self.reduce, str) and self.reduce != "auto":
279 raise ValueError(f"reduce must be an integer, 'auto' or None, not {self.reduce!r}")
281 def options(self, step: str) -> dict[str, Any]:
282 """Preset options for ``"bundle_adjust"``, ``"refine"`` or ``"sfs"`` with overrides applied."""
283 extra = {
284 "bundle_adjust": self.bundle_adjust_options,
285 "refine": self.refine_options,
286 "sfs": self.sfs_options,
287 }[step]
288 return {**PRESETS[self.preset][step], **extra}
291class PipelineError(RuntimeError):
292 """The SfS pipeline cannot continue (e.g. too few images left)."""
295@dataclass
296class SfsImage:
297 """One input image: cube, optional CSM camera, Sun angles and derived files."""
299 product_id: str
300 cube: Path
301 camera: Path | None = None
302 sun_azimuth: float | None = None
303 sun_elevation: float | None = None
304 pixel_scale: float | None = None
305 mapprojected: Path | None = None
306 shadow_fraction: float | None = None
307 slew_angle: float | None = None
310@dataclass(frozen=True)
311class RemovedImage:
312 """An image dropped by a quality filter: *stage* (``ba``, ``refine``, ``selection``) and why."""
314 image: SfsImage
315 stage: str
316 reason: str
319class SfsPipeline:
320 """Build an SfS DEM over a region from NAC images and a prior DEM.
322 [`run`][lunar_toolkit.sfs.SfsPipeline.run] follows Bertone et al. (2026, PSJ 7:118) and the ASP large-area
323 SfS recipe:
325 1. calibrate the images and write CSM cameras (``images/``);
326 2. crop the prior DEM and resample it to the target grid (``dem/``);
327 3. sort by Sun azimuth (``sfs --query``);
328 4. mapproject at native resolution (``prj/``) and bundle-adjust with matches between
329 the mapprojected images (``ba/``);
330 5. drop images with no matches or mostly large residuals (``selection/removed.csv``);
331 6. register the cameras to the prior DEM, translation only, from stereo DEMs
332 (``stereo/``) or the BA points; apply to all cameras (``align/``, ``ba_align/``);
333 7. refine with ``--heights-from-dem`` and filter again, now also on mapprojected
334 match offsets (``ba_refine/``); fewer than ``min_images`` left stops the run;
335 8. mapproject at the target resolution (``prj_sfs/``), write max-lit QA mosaics,
336 all images and per 30° Sun-azimuth bin (``qa/``);
337 9. select an illumination-diverse subset (``selection/``);
338 10. reduce the selected cubes to the target resolution and run ``parallel_sfs`` (``sfs/``);
339 11. products (``products/`` via [`lunar_toolkit.sfs_products`][lunar_toolkit.sfs_products], else hillshade
340 + diff).
342 Steps whose outputs exist are skipped, so a rerun resumes.
344 Parameters
345 ----------
346 work_dir : Path
347 Directory for all intermediate and final files.
348 prior_dem : Path
349 Reference terrain, e.g. a LOLA polar stereographic LDEM GeoTIFF.
350 config : SfsConfig, optional
351 force : bool
352 Recompute steps even when their outputs exist.
353 observations : DataFrame, optional
354 Observation table ([`lunar_toolkit.observation`][lunar_toolkit.observation] schema) with the input
355 images, for native pixel scales and stereo pairs; looked up in the local
356 EDR ``CumIndex`` when needed and not given.
357 aoi : optional
358 Region the stereo pairs must overlap in: anything
359 [`lunar_toolkit.stereo.find_stereo_pairs`][lunar_toolkit.stereo.find_stereo_pairs] takes (dict,
360 GeoJSON / WKT).
361 data_dir : Path, optional
362 ``CumIndex`` data directory.
363 run_config : RunConfig, optional
364 The full run settings ([`lunar_toolkit.config`][lunar_toolkit.config]) written to ``config.yaml``
365 and ``provenance.json``; default: *work_dir*, *prior_dem*, the sources and *config*.
366 provenance : bool
367 Write ``config.yaml`` and append this run to ``provenance.json`` (never in a dry run).
369 Attributes
370 ----------
371 all_images, images : list of SfsImage
372 Every input / those passing the quality filters (Sun-azimuth order).
373 removed : list of RemovedImage
374 selected, unselected : list of SfsImage
375 The SfS inputs and the remaining good images (e.g. for an uncertainty ensemble);
376 both carry ``mapprojected`` (``prj_sfs/``) and ``sun_azimuth``.
377 dem : Path
378 The cropped prior DEM on the target grid.
380 Example
381 -------
382 ::
384 cfg = SfsConfig(resolution=5, projwin=(-10000, 50000, 0, 40000))
385 p = SfsPipeline("work/malapert", "LDEM_87S_5MPP_ADJ.TIF", cfg)
386 dem = p.run(["raw/M1234567890LE.IMG", ...])
387 """
389 def __init__(
390 self,
391 work_dir: Path,
392 prior_dem: Path,
393 config: SfsConfig | None = None,
394 force: bool = False,
395 *,
396 observations: pd.DataFrame | None = None,
397 aoi: Any = None,
398 data_dir: Path | None = None,
399 run_config: RunConfig | None = None,
400 provenance: bool = True,
401 ) -> None:
402 self.work_dir = Path(work_dir)
403 self.prior_dem = Path(prior_dem)
404 self.config = config or SfsConfig()
405 self.force = force
406 self.observations = observations
407 self.aoi = aoi
408 self.data_dir = data_dir
409 self.all_images: list[SfsImage] = []
410 self.images: list[SfsImage] = []
411 self.removed: list[RemovedImage] = []
412 self.selected: list[SfsImage] = []
413 self.unselected: list[SfsImage] = []
414 self.dem = self.work_dir / "dem" / f"prior-{self.config.resolution:g}m.tif"
415 self._registered = False
416 self.run_config = run_config
417 self.provenance = provenance
418 self._prov: Provenance | None = None
420 def run(self, sources: Sequence[Path]) -> Path:
421 """All steps; returns the final SfS DEM.
423 Writes ``config.yaml`` and appends a run to ``provenance.json`` first (see
424 [`start_provenance`][lunar_toolkit.sfs.SfsPipeline.start_provenance]), then records each step there.
425 """
426 sources = [Path(s) for s in sources]
427 self.start_provenance(sources)
428 step = self._step
429 try:
430 step("prepare_images", self.prepare_images, sources)
431 step("prepare_dem", self.prepare_dem)
432 step("sort_by_sun", self.sort_by_sun)
433 step(
434 "mapproject",
435 self.mapproject,
436 self.images,
437 "prj",
438 self.config.ba_resolution or self.config.resolution,
439 )
440 if self.config.cameras_from == "published":
441 step("published_cameras", self.use_published_cameras)
442 else:
443 step("bundle_adjust", self.bundle_adjust)
444 step("filter_ba", self.filter_images, self.ba_prefix, "ba")
445 step("register", self.register)
446 if self.config.refine:
447 step("refine", self.refine)
448 step("filter_refine", self.filter_images, self.refine_prefix, "refine")
449 step("mapproject_for_sfs", self.mapproject_for_sfs)
450 step("qa_mosaics", self.qa_mosaics)
451 step("select_images", self.select_images)
452 dem = step("shape_from_shading", self.shape_from_shading)
453 step("products", self.products, dem)
454 if self.config.ensemble_members:
455 step("ensemble", self.ensemble)
456 except BaseException:
457 if self._prov:
458 self._prov.finish("error")
459 raise
460 if self._prov: 460 ↛ 461line 460 didn't jump to line 461 because the condition on line 460 was never true
461 self._prov.finish("ok")
462 return dem
464 def effective_config(self, sources: Sequence[Path] = ()) -> RunConfig:
465 """`run_config` with this pipeline's work dir, prior DEM and
466 [`SfsConfig`][lunar_toolkit.sfs.SfsConfig]."""
467 from .config import RunConfig
469 if self.run_config is None: 469 ↛ 470line 469 didn't jump to line 470 because the condition on line 469 was never true
470 return RunConfig(
471 work_dir=self.work_dir,
472 dem=str(self.prior_dem),
473 sources=[str(s) for s in sources],
474 sfs=self.config,
475 )
476 return dataclasses.replace(
477 self.run_config, work_dir=self.work_dir, dem=str(self.prior_dem), sfs=self.config
478 )
480 def start_provenance(self, sources: Sequence[Path]) -> Path | None:
481 """Write ``config.yaml`` and start a ``provenance.json`` run; returns the config path.
483 Inputs are *sources* and the prior DEM. Nothing is written in a dry run
484 or with ``provenance=False``.
485 """
486 if is_dry_run() or not self.provenance:
487 return None
488 from .config import dump_config
489 from .provenance import Provenance
491 cfg = self.effective_config(sources)
492 path = self.work_dir / "config.yaml"
493 dump_config(cfg, path)
494 self._prov = Provenance(self.work_dir, cfg, [*sources, self.prior_dem], checksums=cfg.checksums)
495 return path
497 def _step(self, name: str, func: Any, *args: Any) -> Any:
498 if self._prov is None:
499 return func(*args)
500 with self._prov.step(name) as entry:
501 entry["result"] = result = func(*args)
502 return result
504 @property
505 def ba_prefix(self) -> Path:
506 return self.work_dir / "ba" / "run"
508 @property
509 def align_prefix(self) -> Path:
510 return self.work_dir / "align" / "run"
512 @property
513 def ba_align_prefix(self) -> Path:
514 return self.work_dir / "ba_align" / "run"
516 @property
517 def refine_prefix(self) -> Path:
518 return self.work_dir / "ba_refine" / "run"
520 @property
521 def sfs_camera_prefix(self) -> Path:
522 """Adjustments renamed for the reduced ISIS cubes."""
523 return self.work_dir / "sfs_cameras" / "run"
525 @property
526 def sfs_prefix(self) -> Path:
527 return self.work_dir / "sfs" / "run"
529 @property
530 def sfs_dem(self) -> Path:
531 return Path(f"{self.sfs_prefix}-DEM-final.tif")
533 @property
534 def camera_prefix(self) -> Path | None:
535 """Adjustments of the final (registered, refined) cameras; None for published cameras."""
536 if self.config.cameras_from == "published":
537 return None
538 if self.config.refine:
539 return self.refine_prefix
540 return self.ba_align_prefix if self.config.registration != "none" else self.ba_prefix
542 @property
543 def shadow_threshold(self) -> float:
544 return float(self.config.options("sfs").get("shadow_threshold") or 0.0)
546 def _skip(self, path: Path) -> bool:
547 return reuse_existing(path, self.force)
549 def _parallel_opts(self) -> dict[str, Any]:
550 return {"processes": self.config.processes, "threads": self.config.threads}
552 def _cameras(self, images: Sequence[SfsImage]) -> list[Path]:
553 return [im.camera for im in images if im.camera] if self.config.csm else []
555 def _map_path(self, subdir: str, im: SfsImage) -> Path:
556 return self.work_dir / subdir / f"{im.product_id}.map.tif"
558 def prepare_images(self, sources: Sequence[Path]) -> list[SfsImage]:
559 """Calibrate EDRs (``.IMG``) or take ready cubes; write CSM cameras."""
560 cfg = self.config
561 scales = _by_id(self.observations, "pixel_scale")
562 slews = _by_id(self.observations, "slew_angle")
563 self.all_images = []
564 for src in map(Path, sources):
565 cube, camera = isis.prepare_lronac(
566 src, self.work_dir / "images", csm=cfg.csm, web=cfg.web_spice, force=self.force
567 )
568 pid = src.name.split(".")[0]
569 self.all_images.append(
570 SfsImage(
571 pid, cube, camera, pixel_scale=scales.get(pid.upper()), slew_angle=slews.get(pid.upper())
572 )
573 )
574 self.images = list(self.all_images)
575 self._record_observations()
576 return self.images
578 def _record_observations(self) -> None:
579 """Write the given observation rows to the workspace ``images`` table (not in a dry run)."""
580 if is_dry_run() or self.observations is None or not len(self.observations):
581 return
582 from .workspace import WorkspaceError
584 try:
585 with self.workspace() as ws:
586 ws.add_observations(self.observations)
587 except WorkspaceError as exc:
588 log.warning("Workspace images not written: %s", exc)
590 def prepare_dem(self) -> Path:
591 """Crop the prior DEM to ``projwin`` and resample to the target grid."""
592 cfg = self.config
593 crop = self.work_dir / "dem" / "prior.tif"
594 if not self._skip(crop):
595 asp.gdal_translate(self.prior_dem, crop, projwin=cfg.projwin)
596 if not self._skip(self.dem):
597 asp.gdal_translate(crop, self.dem, tr=cfg.resolution, r="cubicspline")
598 return self.dem
600 def sort_by_sun(self) -> list[SfsImage]:
601 """Fill Sun azimuth / elevation and sort the images by azimuth.
603 Bundle adjustment matches neighbours in this order, and ``parallel_sfs``
604 expects it.
605 """
606 angles = asp.sun_angles(
607 self.dem,
608 [im.cube for im in self.images],
609 self._cameras(self.images),
610 out_prefix=self.work_dir / "query" / "run",
611 )
612 by_image = {Path(a.image).name: a for a in angles}
613 for im in self.images:
614 a = by_image.get(im.cube.name)
615 if a:
616 im.sun_azimuth, im.sun_elevation = a.azimuth, a.elevation
617 if angles:
618 self.images.sort(key=lambda im: (im.sun_azimuth is None, im.sun_azimuth or 0.0))
619 elif not is_dry_run(): 619 ↛ 620line 619 didn't jump to line 620 because the condition on line 619 was never true
620 log.warning("sfs --query reported no Sun angles; keeping input order")
621 return self.images
623 def mapproject(
624 self, images: Sequence[SfsImage], subdir: str, resolution: float, **options: Any
625 ) -> list[Path]:
626 """Mapproject *images* onto the prior DEM (cropped to ``projwin``) into ``<subdir>/``."""
627 cfg = self.config
628 win = cfg.projwin
629 outs = []
630 for im in images:
631 out = self._map_path(subdir, im)
632 if not self._skip(out):
633 asp.mapproject(
634 self.dem,
635 im.cube,
636 out,
637 im.camera if cfg.csm else None,
638 tr=resolution,
639 t_projwin=(win[0], win[3], win[2], win[1]) if win else None,
640 threads=cfg.threads,
641 **options,
642 )
643 outs.append(out)
644 return outs
646 def bundle_adjust(self) -> Path:
647 """``parallel_bundle_adjust`` with matches from the mapprojected images.
649 ``--mapproj-dem`` writes the per-image match-offset report used by
650 [`filter_images`][lunar_toolkit.sfs.SfsPipeline.filter_images].
651 """
652 done = Path(f"{self.ba_prefix}-final_residuals_stats.txt")
653 if self._skip(done):
654 return self.ba_prefix
655 opts = self._ba_options(self.config.options("bundle_adjust"))
656 return asp.bundle_adjust(
657 [im.cube for im in self.images],
658 self.ba_prefix,
659 self._cameras(self.images),
660 mapprojected=[self._map_path("prj", im) for im in self.images],
661 dem=self.dem,
662 parallel=True,
663 mapproj_dem=self.dem,
664 **{**self._parallel_opts(), **opts},
665 )
667 def _ba_options(self, opts: dict[str, Any]) -> dict[str, Any]:
668 """Preset bundle-adjust options made valid for ASP.
670 ``auto_overlap_percent`` becomes ``--auto-overlap-params "<dem> <pct>"``, and of
671 ASP's mutually exclusive pair-selection options only the most specific is kept
672 (overlap list, footprint overlap, overlap buffer, then consecutive-image limit).
673 """
674 opts = dict(opts)
675 pct = opts.pop("auto_overlap_percent", None)
676 if pct is not None and "auto_overlap_params" not in opts:
677 opts["auto_overlap_params"] = f"{self.dem} {pct:g}"
678 given = [k for k in _OVERLAP_OPTIONS if opts.get(k) not in (None, False)]
679 for key in given[1:]:
680 log.info(
681 "bundle_adjust: --%s overrides --%s; dropping the latter",
682 given[0].replace("_", "-"),
683 key.replace("_", "-"),
684 )
685 del opts[key]
686 return opts
688 def workspace(self, read_only: bool = False) -> Workspace:
689 """The work dir's DuckDB workspace (``lunar.duckdb``,
690 [`lunar_toolkit.workspace`][lunar_toolkit.workspace])."""
691 from .workspace import Workspace
693 return Workspace(self.work_dir, threads=self.config.threads, read_only=read_only)
695 def filter_images(self, prefix: Path, stage: str) -> list[RemovedImage]:
696 """Drop images that ASP's bundle-adjustment reports flag (see
697 [`flag_images`][lunar_toolkit.sfs.flag_images]).
699 Loads ``<prefix>-*`` reports (residual stats and raw pixels, match offsets, point
700 maps, convergence angles) into the workspace as run ``<prefix>``'s directory name
701 (``ba``, ``ba_refine``) and removes the images of its ``image_rejections`` view,
702 with this config's thresholds. Removed images are logged and appended to
703 ``selection/removed.csv``.
705 At ``ba`` the match-offset rule waits for the refinement when that runs: before
706 registration the network floats tens of metres off the DEM, so convergent pairs
707 mapproject apart by ``dh * tan(convergence)`` however well they are tied (ASP:
708 "do not use this option for any initial evaluation of bundle adjustment").
709 """
710 cfg = self.config
711 files = {
712 k: Path(f"{prefix}-{name}")
713 for k, name in (
714 ("stats", "final_residuals_stats.txt"),
715 ("raw", "final_residuals_raw_pixels.txt"),
716 ("offsets", "mapproj_match_offset_stats.txt"),
717 )
718 }
719 if not any(p.exists() for p in files.values()):
720 if is_dry_run():
721 return []
722 raise PipelineError(f"No bundle-adjustment reports under {prefix}")
723 overrides: dict[str, Any] = {}
724 if stage == "ba" and cfg.max_offset_p85 is not None:
725 if cfg.refine:
726 overrides["max_offset_p85"] = None
727 else:
728 log.warning("Match offsets of unregistered cameras overstate misalignment; consider refine")
729 with self.workspace() as ws:
730 ws.set_rejection_params(cfg, **overrides)
731 ws.load_bundle_adjust(prefix, run=Path(prefix).parent.name)
732 flags = ws.rejections(Path(prefix).parent.name)
733 removed = [
734 RemovedImage(im, stage, flags[im.cube.name]) for im in self.images if im.cube.name in flags
735 ]
736 return self._remove(removed)
738 def _remove(self, removed: list[RemovedImage]) -> list[RemovedImage]:
739 if not removed:
740 return []
741 gone = {id(r.image) for r in removed}
742 self.images = [im for im in self.images if id(im) not in gone]
743 self.removed += removed
744 for r in removed:
745 log.info("Removed %s (%s): %s", r.image.product_id, r.stage, r.reason)
746 log.info("%d images removed at %s, %d remain", len(removed), removed[0].stage, len(self.images))
747 _write_csv(
748 self.work_dir / "selection" / "removed.csv",
749 ["product_id", "stage", "reason"],
750 [[r.image.product_id, r.stage, r.reason] for r in self.removed],
751 )
752 if len(self.images) < self.config.min_images:
753 raise PipelineError(
754 f"{len(self.images)} images left after the {removed[0].stage} filter, fewer than "
755 f"min_images={self.config.min_images}; see {self.work_dir / 'selection' / 'removed.csv'}"
756 )
757 return removed
759 def register(self) -> Path:
760 """Register the adjusted cameras to the prior DEM (translation only); returns their prefix."""
761 mode = self.config.registration
762 if mode == "none":
763 return self.ba_prefix
764 transform = Path(f"{self.align_prefix}-transform.txt")
765 if not self._skip(transform):
766 if mode != "stereo" or self._stereo_mosaic() is None:
767 self._align_residuals()
768 else:
769 self._align_stereo()
770 marker = Path(f"{self.ba_align_prefix}-applied-transform.txt")
771 if not self._skip(marker):
772 asp.apply_transform(
773 [im.cube for im in self.images],
774 transform,
775 self.ba_align_prefix,
776 self._cameras(self.images),
777 input_adjustments_prefix=self.ba_prefix,
778 )
779 if not is_dry_run(): 779 ↛ 780line 779 didn't jump to line 780 because the condition on line 779 was never true
780 marker.write_text(transform.read_text())
781 return self.ba_align_prefix
783 def _align_residuals(self) -> Path:
784 return asp.pc_align(
785 self.dem,
786 Path(f"{self.ba_prefix}-final_residuals_pointmap.csv"),
787 self.align_prefix,
788 max_displacement=self.config.max_displacement,
789 compute_translation_only=True,
790 csv_format="1:lon 2:lat 3:height_above_datum",
791 save_transformed_source_points=True,
792 )
794 def _align_stereo(self) -> Path:
795 return asp.pc_align(
796 self.dem,
797 self.work_dir / "stereo" / "mosaic.tif",
798 self.align_prefix,
799 max_displacement=self.config.max_displacement,
800 compute_translation_only=True,
801 )
803 def stereo_pairs(self) -> pd.DataFrame:
804 """Up to ``stereo_pairs`` kept-image pairs ranked by
805 [`find_stereo_pairs`][lunar_toolkit.stereo.find_stereo_pairs]."""
806 from .stereo import find_stereo_pairs
808 obs = self._observations()
809 kept = {im.product_id.upper() for im in self.images}
810 obs = obs[obs["product_id"].astype(str).str.upper().isin(kept)]
811 pairs = find_stereo_pairs(obs, aoi=self.aoi, convergence=self.config.stereo_convergence)
812 return pairs.head(self.config.stereo_pairs)
814 def _stereo_mosaic(self) -> Path | None:
815 """Stereo DEMs (triangulated with the ``ba/`` cameras) mosaicked; None without any."""
816 cfg = self.config
817 mosaic = self.work_dir / "stereo" / "mosaic.tif"
818 if self._skip(mosaic): 818 ↛ 819line 818 didn't jump to line 819 because the condition on line 818 was never true
819 return mosaic
820 try:
821 pairs = self.stereo_pairs()
822 except (OSError, KeyError, ValueError) as exc:
823 log.warning("No observations for stereo registration (%s); using BA points", exc)
824 return None
825 by_id = {im.product_id.upper(): im for im in self.images}
826 dems = []
827 for _, pair in pairs.iterrows():
828 a, b = by_id[str(pair["product_id_1"]).upper()], by_id[str(pair["product_id_2"]).upper()]
829 try:
830 dem = stereo_dem(
831 a.cube,
832 b.cube,
833 self.work_dir / "stereo" / f"{a.product_id}__{b.product_id}",
834 resolution=cfg.resolution,
835 prior_dem=self.dem,
836 align=False,
837 csm=cfg.csm,
838 cameras=(a.camera, b.camera) if cfg.csm and a.camera and b.camera else None,
839 mapprojected=(self._map_path("prj", a), self._map_path("prj", b)),
840 bundle_adjust_prefix=self.ba_prefix,
841 processes=cfg.processes,
842 threads=cfg.threads,
843 force=self.force,
844 )
845 except ToolError as exc:
846 log.warning("Stereo %s / %s failed: %s", a.product_id, b.product_id, exc)
847 continue
848 if is_dry_run() or dem.exists(): 848 ↛ 827line 848 didn't jump to line 827 because the condition on line 848 was always true
849 dems.append(dem)
850 if not dems:
851 log.warning("No stereo DEMs for registration; using BA points")
852 return None
853 return asp.dem_mosaic(dems, mosaic)
855 def _observations(self) -> pd.DataFrame:
856 if self.observations is None:
857 from .lroc.cumindex import CumIndex
859 ci = CumIndex("edr", data_dir=self.data_dir)
860 self.observations = ci.to_observations(ci.lookup([im.product_id for im in self.all_images]))
861 return self.observations
863 def refine(self) -> Path:
864 """Bundle-adjust the kept, registered cameras again with ``--heights-from-dem``.
866 Reuses the ``ba/`` matches (``--match-files-prefix``); ``--mapproj-dem`` gives
867 the match-offset report for a second [`filter_images`][lunar_toolkit.sfs.SfsPipeline.filter_images].
868 """
869 done = Path(f"{self.refine_prefix}-final_residuals_stats.txt")
870 if self._skip(done):
871 return self.refine_prefix
872 start = self.ba_align_prefix if self.config.registration != "none" else self.ba_prefix
873 opts = self._ba_options(self.config.options("refine"))
874 return asp.bundle_adjust(
875 [im.cube for im in self.images],
876 self.refine_prefix,
877 self._cameras(self.images),
878 match_files_prefix=self.ba_prefix,
879 input_adjustments_prefix=start,
880 heights_from_dem=self.dem,
881 mapproj_dem=self.dem,
882 threads=self.config.threads,
883 **opts,
884 )
886 def mapproject_for_sfs(self) -> list[Path]:
887 """Mapproject the kept images at the target resolution with the final cameras (``prj_sfs/``)."""
888 outs = self.mapproject(
889 self.images, "prj_sfs", self.config.resolution, bundle_adjust_prefix=self.camera_prefix
890 )
891 for im, out in zip(self.images, outs):
892 im.mapprojected = out
893 return outs
895 def qa_mosaics(self, bin_width: float = 30.0) -> dict[str, Path]:
896 """Max-lit mosaics of all kept images and per Sun-azimuth bin (``qa/``).
898 Ghosting in these shows misregistered images (Bertone et al. 2026, App. A).
899 Bins are named by their start, e.g. ``max_lit_az030.tif`` for 30–60°.
900 """
901 qa = self.work_dir / "qa"
902 maps = [im.mapprojected or self._map_path("prj_sfs", im) for im in self.images]
903 groups: dict[str, list[Path]] = {"max_lit": maps}
904 for im, m in zip(self.images, maps):
905 if im.sun_azimuth is not None:
906 start = int((im.sun_azimuth % 360.0) // bin_width * bin_width)
907 groups.setdefault(f"max_lit_az{start:03d}", []).append(m)
908 out = {}
909 for name, paths in groups.items():
910 out[name] = qa / f"{name}.tif"
911 if paths and not self._skip(out[name]): 911 ↛ 909line 911 didn't jump to line 909 because the condition on line 911 was always true
912 asp.dem_mosaic(paths, out[name], max=True)
913 return out
915 def select_images(self) -> list[SfsImage]:
916 """Split the kept images into SfS inputs and the rest (``selection/``).
918 Images at least ``max_shadow_fraction`` shadowed are removed first. Writes
919 ``selection.csv`` (every kept image: azimuth, shadow fraction, mapprojected path,
920 selected) plus ``selected.txt`` / ``unselected.txt`` (one product ID per line).
921 """
922 cfg = self.config
923 out_dir = self.work_dir / "selection"
924 lit = self._lit_masks()
925 if lit is not None and cfg.max_shadow_fraction is not None:
926 shadowed = [
927 RemovedImage(im, "selection", f"{im.shadow_fraction:.1%} shadowed")
928 for im in self.images
929 if im.shadow_fraction is not None and im.shadow_fraction >= cfg.max_shadow_fraction
930 ]
931 gone = {id(r.image) for r in shadowed}
932 lit = lit[[k for k, im in enumerate(self.images) if id(im) not in gone]]
933 self._remove(shadowed)
934 eligible = self._sfs_eligible()
935 if cfg.selection == "all" or lit is None:
936 chosen = eligible
937 elif cfg.selection == "coverage":
938 az = np.array([np.nan if im.sun_azimuth is None else im.sun_azimuth for im in self.images])
939 picked, views = select_by_illumination(
940 lit[eligible],
941 az[eligible],
942 views=cfg.selection_views,
943 azimuth_sep=cfg.selection_azimuth_sep,
944 max_images=cfg.max_sfs_images,
945 )
946 chosen = [eligible[k] for k in picked]
947 covered = lit[eligible].any(axis=0)
948 if covered.any(): 948 ↛ 978line 948 didn't jump to line 978 because the condition on line 948 was always true
949 log.info(
950 "Selected %d / %d nadir images; lit cells with >=%d directions: %.1f%%",
951 len(chosen),
952 len(eligible),
953 cfg.selection_views,
954 100 * np.mean(views[covered] >= cfg.selection_views),
955 )
956 elif cfg.selection == "kmedoids":
957 az = np.array([np.nan if im.sun_azimuth is None else im.sun_azimuth for im in self.images])
958 res = select_kmedoids(
959 lit[eligible],
960 az[eligible],
961 cfg.selection_clusters,
962 views=cfg.selection_views,
963 azimuth_sep=cfg.selection_azimuth_sep,
964 max_images=cfg.max_sfs_images,
965 seed=cfg.seed,
966 )
967 chosen = sorted(eligible[k] for k in res.chosen)
968 log.info(
969 "k-medoids selected %d / %d nadir images (%d medoids); uncovered %.2f%%, poorly lit %.2f%%",
970 len(chosen),
971 len(eligible),
972 len(res.medoids),
973 100 * res.uncovered,
974 100 * res.poorly_lit,
975 )
976 else:
977 chosen = sorted(set(self._image_subset()) & set(eligible))
978 chosen_set = set(chosen)
979 self.selected = [im for k, im in enumerate(self.images) if k in chosen_set]
980 self.unselected = [im for k, im in enumerate(self.images) if k not in chosen_set]
981 _write_csv(
982 out_dir / "selection.csv",
983 ["product_id", "sun_azimuth", "sun_elevation", "shadow_fraction", "mapprojected", "selected"],
984 [
985 [
986 im.product_id,
987 im.sun_azimuth,
988 im.sun_elevation,
989 im.shadow_fraction,
990 im.mapprojected,
991 k in chosen_set,
992 ]
993 for k, im in enumerate(self.images)
994 ],
995 )
996 for name, ims in (("selected", self.selected), ("unselected", self.unselected)):
997 _write_lines(out_dir / f"{name}.txt", [im.product_id for im in ims])
998 return self.selected
1000 def _sfs_eligible(self) -> list[int]:
1001 """Indices of the kept images allowed into SfS (``|slew| < max_sfs_slew``)."""
1002 limit = self.config.max_sfs_slew
1003 if limit is None:
1004 return list(range(len(self.images)))
1005 if any(im.slew_angle is None for im in self.images):
1006 try:
1007 obs = self._observations()
1008 except (OSError, KeyError, ValueError) as exc:
1009 log.warning("No slew angles (%s); all images count as nadir", exc)
1010 obs = None
1011 if obs is not None and "slew_angle" in obs:
1012 slews = dict(zip(obs["product_id"].astype(str).str.upper(), obs["slew_angle"]))
1013 for im in self.images:
1014 v = slews.get(im.product_id.upper())
1015 if im.slew_angle is None and v is not None and np.isfinite(v): 1015 ↛ 1016line 1015 didn't jump to line 1016 because the condition on line 1015 was never true
1016 im.slew_angle = float(v)
1017 eligible = [
1018 k for k, im in enumerate(self.images) if im.slew_angle is None or abs(im.slew_angle) < limit
1019 ]
1020 if len(eligible) < len(self.images):
1021 log.info(
1022 "%d slewed images (|slew| >= %g deg) kept out of SfS", len(self.images) - len(eligible), limit
1023 )
1024 return eligible
1026 def _lit_masks(self) -> np.ndarray | None:
1027 """``(images, cells)`` lit mask on a common coarse grid; fills ``shadow_fraction``."""
1028 bounds = asp.raster_bounds(self.dem)
1029 if bounds is None or not self.images:
1030 return None
1031 ulx, uly, lrx, lry = bounds
1032 n = self.config.selection_grid
1033 w, h = abs(lrx - ulx), abs(uly - lry)
1034 size = (n, max(1, round(n * h / w))) if w >= h else (max(1, round(n * w / h)), n)
1035 thr = self.shadow_threshold
1036 masks = []
1037 for im in self.images:
1038 r = asp.read_raster(
1039 im.mapprojected or self._map_path("prj_sfs", im), projwin=bounds, outsize=size
1040 )
1041 v = r.data.ravel() if r is not None else np.full(size[0] * size[1], np.nan)
1042 valid = np.isfinite(v) & (v > -1e30)
1043 lit = valid & (v > thr)
1044 im.shadow_fraction = 1.0 - float(lit.sum()) / float(valid.sum()) if valid.any() else 1.0
1045 masks.append(lit)
1046 return np.array(masks)
1048 def _image_subset(self, bin_width: float = 30.0) -> list[int]:
1049 """ASP ``image_subset`` in each Sun-azimuth bin; returns indices of the union."""
1050 bounds = asp.raster_bounds(self.dem)
1051 win = (bounds[0], bounds[3], bounds[2], bounds[1]) if bounds else None
1052 bins: dict[int, list[int]] = {}
1053 for k, im in enumerate(self.images):
1054 bins.setdefault(int(((im.sun_azimuth or 0.0) % 360.0) // bin_width), []).append(k)
1055 by_path = {str(im.mapprojected): k for k, im in enumerate(self.images)}
1056 chosen: list[int] = []
1057 for b, idx in sorted(bins.items()):
1058 subset = asp.image_subset(
1059 [self.images[k].mapprojected or self._map_path("prj_sfs", self.images[k]) for k in idx],
1060 self.work_dir / "selection" / f"image_subset_az{int(b * bin_width):03d}.txt",
1061 threshold=self.shadow_threshold,
1062 t_projwin=win,
1063 threads=self.config.threads,
1064 )
1065 chosen += [by_path[p] for p, _ in subset if p in by_path]
1066 return sorted(set(chosen))
1068 def reduce_factor(self, im: SfsImage) -> int:
1069 """ISIS ``reduce`` factor for *im* (1 = keep native resolution)."""
1070 cfg = self.config
1071 if cfg.reduce is None:
1072 return 1
1073 if cfg.reduce == "auto":
1074 native = im.pixel_scale or cfg.native_resolution
1075 return max(1, int(round(cfg.resolution / native)))
1076 return max(1, int(cfg.reduce))
1078 def sfs_inputs(self, images: Sequence[SfsImage]) -> tuple[list[Path], list[Path], Path | None, bool]:
1079 """Cubes, cameras and adjustment prefix for SfS on *images*.
1081 Cubes needing it are ``reduce``-d to the target resolution (Bertone et al. 2026,
1082 App. A: SfS's own resampling leaves artifacts). With ``sfs_cameras="csm"`` each
1083 CSM camera is re-mapped to its reduced cube
1084 ([`csm.isis_reduce_with_csm`][lunar_toolkit.csm.isis_reduce_with_csm]); otherwise ISIS cameras are
1085 used. The final ``.adjust`` files (resolution independent) are copied under
1086 [`sfs_camera_prefix`][lunar_toolkit.sfs.SfsPipeline.sfs_camera_prefix] with the new names; published
1087 cameras need none. Returns ``(cubes, cameras, prefix, isis_cameras)``.
1088 """
1089 from . import csm
1091 cfg = self.config
1092 factors = [self.reduce_factor(im) for im in images]
1093 prefix = self.camera_prefix
1094 if cfg.csm and all(f == 1 for f in factors):
1095 return [im.cube for im in images], self._cameras(images), prefix, False
1096 use_csm = cfg.csm and cfg.sfs_cameras == "csm"
1097 cubes: list[Path] = []
1098 cameras: list[Path] = []
1099 for im, f in zip(images, factors):
1100 cube, camera = im.cube, im.camera if use_csm else None
1101 if f > 1 and use_csm and im.camera is not None:
1102 cube, camera = csm.isis_reduce_with_csm(
1103 im.cube, im.camera, f, self.work_dir / "images", force=self.force
1104 )
1105 elif f > 1: 1105 ↛ 1109line 1105 didn't jump to line 1109 because the condition on line 1105 was always true
1106 cube = self.work_dir / "images" / f"{im.product_id}.red{f}.cub"
1107 if not self._skip(cube):
1108 isis.reduce(im.cube, cube, f)
1109 if prefix is not None:
1110 src = asp.bundle_adjust_file_name(prefix, im.cube, im.camera if cfg.csm else None)
1111 dst = asp.bundle_adjust_file_name(self.sfs_camera_prefix, cube, camera)
1112 if not is_dry_run() and not self._skip(dst):
1113 dst.parent.mkdir(parents=True, exist_ok=True)
1114 shutil.copyfile(src, dst)
1115 cubes.append(cube)
1116 if camera is not None:
1117 cameras.append(camera)
1118 out_prefix = self.sfs_camera_prefix if prefix is not None else None
1119 return cubes, cameras, out_prefix, not use_csm
1121 def use_published_cameras(self) -> list[SfsImage]:
1122 """Replace the cameras with the published Bertone et al. (2026) adjusted ones.
1124 Each published camera (a CSM state of the reduced full image) is re-mapped to
1125 the image of its cube
1126 ([`reference.published_camera_native`][lunar_toolkit.reference.published_camera_native]) and written
1127 to ``images/<id>.published.json``; images without one are removed.
1128 """
1129 from . import csm
1130 from . import reference as ref
1132 region = self.config.published_region
1133 if region is None: 1133 ↛ 1138line 1133 didn't jump to line 1138 because the condition on line 1133 was always true
1134 rid = self.config.region_id or ""
1135 if not (rid.startswith("A3") and rid[2:].isdigit()): 1135 ↛ 1136line 1135 didn't jump to line 1136 because the condition on line 1135 was never true
1136 raise ValueError("cameras_from: published needs published_region (or region_id A3NN)")
1137 region = int(rid[2:])
1138 out_dir = self.work_dir / "images"
1139 if is_dry_run(): 1139 ↛ 1143line 1139 didn't jump to line 1143 because the condition on line 1139 was always true
1140 for im in self.images:
1141 im.camera = out_dir / f"{im.product_id}.published.json"
1142 return self.images
1143 published = ref.published_cameras(
1144 region, self.work_dir / "published", images=[im.product_id for im in self.images]
1145 )
1146 missing = [
1147 RemovedImage(im, "published", "no published camera")
1148 for im in self.images
1149 if im.product_id not in published
1150 ]
1151 self._remove(missing)
1152 for im in self.images:
1153 camera = out_dir / f"{im.product_id}.published.json"
1154 if not self._skip(camera):
1155 geom = csm.CubeGeometry.from_cube(im.cube)
1156 native = (geom.alpha_samples, geom.alpha_lines)
1157 csm.write_camera(ref.published_camera_native(published[im.product_id], native), camera)
1158 im.camera = camera
1159 return self.images
1161 def shape_from_shading(
1162 self,
1163 images: Sequence[SfsImage] | None = None,
1164 prefix: Path | None = None,
1165 dem: Path | None = None,
1166 **options: Any,
1167 ) -> Path:
1168 """``parallel_sfs`` on the target-grid prior DEM; returns ``<prefix>-DEM-final.tif``.
1170 *images* default to `selected` (else all kept images), *prefix* to
1171 ``sfs/run`` and *dem* (input DEM, e.g. a crop) to `dem`. *options*
1172 override the configured SfS options.
1173 """
1174 images = list(images) if images is not None else (self.selected or self.images)
1175 prefix = Path(prefix) if prefix else self.sfs_prefix
1176 out = Path(f"{prefix}-DEM-final.tif")
1177 if self._skip(out): 1177 ↛ 1178line 1177 didn't jump to line 1178 because the condition on line 1177 was never true
1178 return out
1179 cubes, cameras, ba_prefix, isis_cameras = self.sfs_inputs(images)
1180 opts: dict[str, Any] = {**self.config.options("sfs"), **self._parallel_opts()}
1181 if isis_cameras: 1181 ↛ 1182line 1181 didn't jump to line 1182 because the condition on line 1181 was never true
1182 opts["use_approx_camera_models"] = True
1183 opts.update(options)
1184 return asp.sfs(
1185 dem or self.dem, cubes, prefix, cameras, parallel=True, bundle_adjust_prefix=ba_prefix, **opts
1186 )
1188 def products(self, dem: Path) -> dict[str, Path]:
1189 """Standard products via
1190 [`lunar_toolkit.sfs_products.make_products`][lunar_toolkit.sfs_products.make_products] when available.
1192 Otherwise (module missing or its dependencies not installed) a hillshade of the
1193 SfS DEM and its difference from the prior.
1194 """
1195 try:
1196 products = importlib.import_module(f"{__package__}.sfs_products")
1197 except ImportError as exc:
1198 log.info("sfs_products unavailable (%s); writing hillshade and diff only", exc)
1199 else:
1200 height_error = Path(f"{self.sfs_prefix}-height-error.tif")
1201 return products.make_products(
1202 dem,
1203 self.dem,
1204 self.work_dir / "products",
1205 [
1206 (im.mapprojected, im.sun_azimuth, im.pixel_scale)
1207 for im in self.selected
1208 if im.mapprojected and im.sun_azimuth is not None
1209 ],
1210 region_id=self.config.region_id or self.work_dir.name,
1211 resolution=self.config.resolution,
1212 height_error=height_error if height_error.exists() else None,
1213 shadow_threshold=self.shadow_threshold,
1214 force=self.force,
1215 )
1216 out = {"dem": dem}
1217 shade = Path(f"{self.sfs_prefix}-DEM-final_hillshade.tif")
1218 out["hillshade"] = shade if self._skip(shade) else asp.hillshade(dem, shade)
1219 diff_prefix = Path(f"{self.sfs_prefix}-prior")
1220 diff = Path(f"{diff_prefix}-diff.tif")
1221 out["diff"] = diff if self._skip(diff) else asp.geodiff(dem, self.dem, diff_prefix)
1222 return out
1224 def ensemble(self) -> dict[str, Path]:
1225 """SfS ensemble over a central crop and its spread (``ensemble/``, ``products/``).
1227 Members are [`ensemble_subsets`][lunar_toolkit.sfs_products.ensemble_subsets] of the kept
1228 images; returns [`ensemble_uncertainty`][lunar_toolkit.sfs_products.ensemble_uncertainty]'s maps.
1229 """
1230 from . import sfs_products
1232 cfg = self.config
1233 bounds = asp.raster_bounds(self.dem) or cfg.projwin
1234 if bounds is None: 1234 ↛ 1235line 1234 didn't jump to line 1235 because the condition on line 1234 was never true
1235 log.warning("Ensemble skipped: DEM extent unknown (dry run without projwin)")
1236 return {}
1237 ens = self.work_dir / "ensemble"
1238 crop = ens / "prior.tif"
1239 if not self._skip(crop): 1239 ↛ 1242line 1239 didn't jump to line 1242 because the condition on line 1239 was always true
1240 centre = ((bounds[0] + bounds[2]) / 2, (bounds[1] + bounds[3]) / 2)
1241 asp.gdal_translate(self.dem, crop, projwin=sfs_products.crop_window(*centre, cfg.ensemble_crop))
1242 size = min(cfg.ensemble_size or len(self.selected) or len(self.images), len(self.images))
1243 subsets = sfs_products.ensemble_subsets(self.images, cfg.ensemble_members, size, seed=cfg.seed)
1244 members = [
1245 self.shape_from_shading(subset, ens / f"m{k:02d}" / "run", dem=crop)
1246 for k, subset in enumerate(subsets)
1247 ]
1248 _write_lines(ens / "members.txt", [" ".join(im.product_id for im in subset) for subset in subsets])
1249 if is_dry_run(): 1249 ↛ 1251line 1249 didn't jump to line 1251 because the condition on line 1249 was always true
1250 return {}
1251 return sfs_products.ensemble_uncertainty(
1252 members,
1253 self.work_dir / "products",
1254 region_id=cfg.region_id or self.work_dir.name,
1255 resolution=cfg.resolution,
1256 reference=crop,
1257 force=self.force,
1258 )
1261def flag_images(
1262 residual_stats: dict[str, asp.ResidualStats],
1263 residual_norms: dict[str, np.ndarray],
1264 offsets: dict[str, asp.MatchOffsetStats],
1265 *,
1266 max_offset_p85: float | None = 1.8,
1267 max_large_residual_fraction: float | None = 0.85,
1268 large_residual_px: float = 1.0,
1269 min_matches: int = 1,
1270) -> dict[str, str]:
1271 """Images to drop after bundle adjustment (Bertone et al. 2026, §2), keyed by file name.
1273 An image is flagged when it has fewer than *min_matches* residuals or match offsets,
1274 an 85th-percentile mapprojected match offset over *max_offset_p85* m, or more than
1275 *max_large_residual_fraction* of its reprojection errors over *large_residual_px*.
1276 The value is the reason. The rules are the workspace's ``image_rejections`` view
1277 ([`lunar_toolkit.workspace`][lunar_toolkit.workspace]), evaluated in an in-memory database.
1278 """
1279 from .workspace import flag_images_sql
1281 return flag_images_sql(
1282 residual_stats,
1283 residual_norms,
1284 offsets,
1285 max_offset_p85=max_offset_p85,
1286 max_large_residual_fraction=max_large_residual_fraction,
1287 large_residual_px=large_residual_px,
1288 min_matches=min_matches,
1289 )
1292def select_by_illumination(
1293 lit: np.ndarray,
1294 azimuth: np.ndarray,
1295 *,
1296 views: int = 3,
1297 azimuth_sep: float = 30.0,
1298 max_images: int | None = None,
1299 min_gain: int = 1,
1300) -> tuple[list[int], np.ndarray]:
1301 """Greedy per-pixel illumination-diverse image selection on a raster grid.
1303 The raster analogue of
1304 [`lunar_toolkit.selection.select_for_coverage`][lunar_toolkit.selection.select_for_coverage], after
1305 Bertone et al. (2026, §2): a lit cell is *helped* by an image lighting it if it has fewer than *views*
1306 illumination directions so far and the image's Sun azimuth is at least *azimuth_sep* from each of them.
1307 Each step takes the image helping most cells, until none helps *min_gain* cells or *max_images* are
1308 chosen.
1310 Parameters
1311 ----------
1312 lit : (n_images, n_cells) bool array
1313 azimuth : (n_images,) array
1314 Sun azimuths (deg); NaN never helps.
1316 Returns
1317 -------
1318 (chosen, count)
1319 Image indices in pick order, and the directions reached per cell.
1320 """
1321 lit = np.asarray(lit, bool)
1322 az = np.asarray(azimuth, float)
1323 n, cells = lit.shape
1324 dirs = np.full((cells, views), np.nan)
1325 count = np.zeros(cells, int)
1326 available = np.isfinite(az)
1327 chosen: list[int] = []
1328 while available.any() and (max_images is None or len(chosen) < max_images):
1329 open_cells = count < views
1330 best, best_cells = -1, np.empty(0, int)
1331 for k in np.flatnonzero(available):
1332 idx = np.flatnonzero(lit[k] & open_cells)
1333 if len(idx) <= len(best_cells):
1334 continue
1335 d = azimuth_difference(dirs[idx], az[k])
1336 idx = idx[np.all(np.isnan(d) | (d >= azimuth_sep), axis=1)]
1337 if len(idx) > len(best_cells):
1338 best, best_cells = int(k), idx
1339 if best < 0 or len(best_cells) < min_gain:
1340 break
1341 dirs[best_cells, count[best_cells]] = az[best]
1342 count[best_cells] += 1
1343 available[best] = False
1344 chosen.append(best)
1345 return chosen, count
1348def prefilter_candidates(
1349 df: pd.DataFrame, max_slew: float | None = None, max_resolution: float | None = 1.5
1350) -> pd.DataFrame:
1351 """Images finer than *max_resolution* m/px (and, with *max_slew*, ``|slew| <`` it).
1353 Works on ``CumIndex`` rows (``SLEW_ANGLE``, ``RESOLUTION``) or observation rows
1354 (``slew_angle``, ``pixel_scale``). Rows with unknown values are dropped. The
1355 nadir cut of Bertone et al. (2026, §2) belongs to SfS selection
1356 (`SfsConfig.max_sfs_slew`), so searches keep slewed images by default: they
1357 give the stereo convergence the registration step needs.
1358 """
1359 import pandas as pd
1361 keep = pd.Series(True, index=df.index)
1362 slew_col = next((c for c in ("SLEW_ANGLE", "slew_angle") if c in df), None)
1363 if max_slew is not None:
1364 if slew_col: 1364 ↛ 1367line 1364 didn't jump to line 1367 because the condition on line 1364 was always true
1365 keep &= pd.to_numeric(df[slew_col], errors="coerce").abs() < max_slew
1366 else:
1367 log.warning("No slew column; slew filter skipped")
1368 res_col = next((c for c in ("RESOLUTION", "pixel_scale", "PIXEL_SCALE") if c in df), None)
1369 if max_resolution is not None and res_col: 1369 ↛ 1371line 1369 didn't jump to line 1371 because the condition on line 1369 was always true
1370 keep &= pd.to_numeric(df[res_col], errors="coerce") < max_resolution
1371 return df[keep]
1374def cap_candidates(
1375 obs: pd.DataFrame,
1376 aoi: Any,
1377 max_images: int,
1378 *,
1379 max_sfs_slew: float | None = 5.0,
1380 stereo_pairs: int = 30,
1381 stereo_convergence: tuple[float, float] = (5.0, 60.0),
1382 views: int = 3,
1383 azimuth_sep: float = 30.0,
1384) -> pd.DataFrame:
1385 """At most *max_images* candidate images, chosen before anything is downloaded.
1387 Keeps both images of the best stereo pairs (up to half the budget; they feed the
1388 registration step), then nadir images (``|slew| < max_sfs_slew``) chosen for
1389 illumination coverage over *aoi* ([`select_for_coverage`][lunar_toolkit.selection.select_for_coverage]:
1390 *views* directions at least *azimuth_sep* apart per cell), then the finest remaining
1391 nadir images. Returns the chosen rows in their original order.
1392 """
1393 from .selection import select_for_coverage
1395 if len(obs) <= max_images:
1396 return obs
1397 ids = obs["product_id"].astype(str).str.upper()
1398 chosen: list[str] = []
1400 def take(new: Sequence[str], limit: int) -> None:
1401 for pid in new:
1402 if len(chosen) >= limit:
1403 return
1404 if pid not in chosen:
1405 chosen.append(pid)
1407 if stereo_pairs: 1407 ↛ 1415line 1407 didn't jump to line 1415 because the condition on line 1407 was always true
1408 pairs = _candidate_pairs(obs, aoi, stereo_convergence)
1409 stereo_ids = [
1410 str(v).upper()
1411 for pair in pairs.head(stereo_pairs).itertuples()
1412 for v in (pair.product_id_1, pair.product_id_2)
1413 ]
1414 take(stereo_ids, max_images // 2)
1415 slew = _numeric(obs.get("slew_angle"))
1416 nadir = obs[(slew.abs() < max_sfs_slew) | slew.isna()] if max_sfs_slew is not None else obs
1417 nadir = nadir[~nadir["product_id"].astype(str).str.upper().isin(chosen)]
1418 if len(nadir): 1418 ↛ 1425line 1418 didn't jump to line 1425 because the condition on line 1418 was always true
1419 cover = select_for_coverage(
1420 nadir, aoi, azimuth_sep_deg=azimuth_sep, views_per_cell=views, max_images=max_images - len(chosen)
1421 )
1422 take(cover["product_id"].astype(str).str.upper().tolist(), max_images)
1423 rest = nadir.assign(_res=_numeric(nadir.get("pixel_scale"))).sort_values("_res")
1424 take(rest["product_id"].astype(str).str.upper().tolist(), max_images)
1425 log.info("Capped %d candidates to %d (stereo pairs + illumination coverage)", len(obs), len(chosen))
1426 return obs[ids.isin(set(chosen)).to_numpy()]
1429def _candidate_pairs(obs: pd.DataFrame, aoi: Any, convergence: tuple[float, float]) -> pd.DataFrame:
1430 from .stereo import find_stereo_pairs
1431 from .workspace import PAIRS_ENGINE_THRESHOLD, find_stereo_pairs_duckdb, spatial_available
1433 if len(obs) >= PAIRS_ENGINE_THRESHOLD and spatial_available():
1434 return find_stereo_pairs_duckdb(obs, aoi, convergence=convergence)
1435 return find_stereo_pairs(obs, aoi=aoi, convergence=convergence)
1438def _numeric(values: Any) -> pd.Series:
1439 import pandas as pd
1441 if values is None: 1441 ↛ 1442line 1441 didn't jump to line 1442 because the condition on line 1441 was never true
1442 return pd.Series(dtype=float)
1443 return pd.to_numeric(values, errors="coerce")
1446def search_images(
1447 aoi: Any,
1448 *,
1449 max_slew: float | None = None,
1450 max_resolution: float | None = 1.5,
1451 max_images: int | None = None,
1452 config: SfsConfig | None = None,
1453 data_dir: Path | None = None,
1454) -> pd.DataFrame:
1455 """NAC EDRs over *aoi* from the local ``CumIndex``, pre-filtered and capped; observation schema.
1457 *config* supplies the SfS slew cut and stereo settings used by
1458 [`cap_candidates`][lunar_toolkit.sfs.cap_candidates].
1459 """
1460 from .lroc.cumindex import CumIndex
1462 cfg = config or SfsConfig()
1463 ci = CumIndex("edr", data_dir=data_dir)
1464 df = ci.search(aoi=aoi, nac_only=True, wkt=True)
1465 kept = prefilter_candidates(df, max_slew, max_resolution)
1466 log.info("%d NAC images over the AOI, %d after the resolution cut", len(df), len(kept))
1467 obs = ci.to_observations(kept)
1468 if max_images is not None:
1469 obs = cap_candidates(
1470 obs,
1471 aoi,
1472 max_images,
1473 max_sfs_slew=cfg.max_sfs_slew,
1474 stereo_pairs=cfg.stereo_pairs if cfg.registration == "stereo" else 0,
1475 stereo_convergence=cfg.stereo_convergence,
1476 views=cfg.selection_views,
1477 azimuth_sep=cfg.selection_azimuth_sep,
1478 )
1479 return obs
1482@dataclass(frozen=True)
1483class PolarStereographic:
1484 """A polar stereographic projection (sphere, true scale at the pole) from a PROJ string."""
1486 hemisphere: str
1487 lon_0: float = 0.0
1488 radius: float = MOON_RADIUS_KM * 1000.0
1489 x_0: float = 0.0
1490 y_0: float = 0.0
1492 @classmethod
1493 def from_proj4(cls, proj4: str) -> PolarStereographic:
1494 """Parse e.g. ``+proj=stere +lat_0=-90 +lon_0=0 +k=1 +x_0=0 +y_0=0 +R=1737400 +units=m``.
1496 Raises
1497 ------
1498 ValueError
1499 For anything else (other projections, oblique stereographic, scale ≠ 1, non-metre units).
1500 """
1501 p: dict[str, str] = dict(re.findall(r"\+(\w+)(?:=(\S+))?", proj4))
1502 lat_0 = float(p.get("lat_0", "nan"))
1503 k = float(p.get("k", p.get("k_0", "1")))
1504 lat_ts = float(p["lat_ts"]) if "lat_ts" in p else lat_0
1505 sphere_a = p.get("a") if p.get("a") == p.get("b", p.get("a")) else None
1506 radius = float(p.get("R") or sphere_a or "nan")
1507 if (
1508 p.get("proj") not in ("stere", "ups")
1509 or not math.isclose(abs(lat_0), 90.0)
1510 or not math.isclose(k, 1.0)
1511 or not math.isclose(lat_ts, lat_0)
1512 or p.get("units", "m") != "m"
1513 or not math.isfinite(radius)
1514 ):
1515 raise ValueError(
1516 f"Expected a spherical polar stereographic projection in metres, got {proj4!r}; "
1517 "pass --projwin in the DEM's own units instead"
1518 )
1519 return cls(
1520 "S" if lat_0 < 0 else "N",
1521 float(p.get("lon_0", 0.0)),
1522 radius,
1523 float(p.get("x_0", 0.0)),
1524 float(p.get("y_0", 0.0)),
1525 )
1527 def xy(self, lat: float, lon: float) -> tuple[float, float]:
1528 x, y = ps_xy(lat, lon - self.lon_0, self.hemisphere)
1529 scale = self.radius / (MOON_RADIUS_KM * 1000.0)
1530 return x * scale + self.x_0, y * scale + self.y_0
1533def aoi_projwin(aoi: Any, dem: Path, snap: float | None = None) -> tuple[float, float, float, float]:
1534 """``(ulx, uly, lrx, lry)`` bounding *aoi* in the polar stereographic projection of *dem*.
1536 *aoi* is anything [`lunar_toolkit.geometry.load_aoi`][lunar_toolkit.geometry.load_aoi] takes (GeoJSON /
1537 WKT file, WKT string). The window is widened to multiples of *snap* (m) when given.
1538 """
1539 proj = PolarStereographic.from_proj4(asp.srs_proj4(dem))
1540 pts = [proj.xy(lat, lon) for poly in load_aoi(aoi).polygons for lon, lat in poly[0]]
1541 xs, ys = [p[0] for p in pts], [p[1] for p in pts]
1542 x0, x1, y0, y1 = min(xs), max(xs), min(ys), max(ys)
1543 if snap: 1543 ↛ 1546line 1543 didn't jump to line 1546 because the condition on line 1543 was always true
1544 x0, y0 = math.floor(x0 / snap) * snap, math.floor(y0 / snap) * snap
1545 x1, y1 = math.ceil(x1 / snap) * snap, math.ceil(y1 / snap) * snap
1546 return x0, y1, x1, y0
1549def _by_id(observations: pd.DataFrame | None, column: str) -> dict[str, float]:
1550 """Finite values of a numeric observation *column* keyed by upper-case product ID."""
1551 if observations is None or column not in observations:
1552 return {}
1553 import pandas as pd
1555 values = pd.to_numeric(observations[column], errors="coerce")
1556 ids = observations["product_id"].astype(str).str.upper()
1557 return {i: float(v) for i, v in zip(ids, values) if np.isfinite(v)}
1560def _write_csv(path: Path, header: list[str], rows: list[list[Any]]) -> None:
1561 if is_dry_run():
1562 return
1563 path.parent.mkdir(parents=True, exist_ok=True)
1564 with open(path, "w", newline="") as f:
1565 w = csv.writer(f)
1566 w.writerow(header)
1567 w.writerows([["" if v is None else v for v in r] for r in rows])
1570def _write_lines(path: Path, lines: list[str]) -> None:
1571 if not is_dry_run():
1572 path.parent.mkdir(parents=True, exist_ok=True)
1573 path.write_text("".join(f"{ln}\n" for ln in lines))