Coverage for lunar_toolkit/stereo_dem.py: 93%
50 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 02:58 +0000
1"""Stereo DEM from one LRO NAC pair with ISIS and ASP."""
3from __future__ import annotations
5from pathlib import Path
6from typing import Any
8from . import asp, isis
9from .tools import reuse_existing
11# ASP's recommended correlator settings for NAC stereo.
12STEREO_DEFAULTS: dict[str, Any] = {
13 "stereo_algorithm": "asp_mgm",
14 "subpixel_mode": 9,
15}
17# bundle_adjust of a lone pair before stereo.
18PAIR_BUNDLE_ADJUST_DEFAULTS: dict[str, Any] = {
19 "ip_per_image": 20000,
20 "min_matches": 1,
21 "camera_weight": 0,
22 "tri_weight": 0.1,
23 "tri_robust_threshold": 0.1,
24 "remove_outliers_params": "75.0 3.0 3.0 5.0",
25 "datum": "D_MOON",
26}
29def stereo_dem(
30 left: Path,
31 right: Path,
32 work_dir: Path,
33 *,
34 resolution: float = 2.0,
35 prior_dem: Path | None = None,
36 image_resolution: float = 1.0,
37 align: bool = True,
38 max_displacement: float = 200.0,
39 csm: bool = True,
40 web_spice: bool = True,
41 cameras: tuple[Path, Path] | None = None,
42 mapprojected: tuple[Path, Path] | None = None,
43 bundle_adjust_prefix: Path | None = None,
44 processes: int | None = None,
45 threads: int | None = None,
46 stereo_options: dict[str, Any] | None = None,
47 force: bool = False,
48) -> Path:
49 """Build a DEM from a NAC stereo pair; returns the final DEM path.
51 Parameters
52 ----------
53 left, right : Path
54 NAC EDR ``.IMG`` files or calibrated ``.cub`` files.
55 resolution : float
56 DEM grid spacing (m/px); about 3× the image resolution is typical.
57 prior_dem : Path, optional
58 Reference DEM: inputs are mapprojected onto it at *image_resolution*
59 before correlation, and the result is aligned to it when *align*.
60 cameras : (Path, Path), optional
61 Existing CSM cameras for *left* / *right* (with *csm*; default: written by
62 ``isd_generate`` into ``<work_dir>/images``).
63 mapprojected : (Path, Path), optional
64 Existing images mapprojected onto *prior_dem* (any camera / adjustments:
65 ASP reads those from the geoheader), instead of mapprojecting here.
66 bundle_adjust_prefix : Path, optional
67 Adjustments to triangulate with, instead of bundle-adjusting the pair here.
68 """
69 work_dir = Path(work_dir)
70 img_dir = work_dir / "images"
72 cubes: list[Path] = []
73 cams: list[Path] = []
74 for k, src in enumerate((Path(left), Path(right))):
75 cube, cam = isis.prepare_lronac(src, img_dir, csm=csm and not cameras, web=web_spice, force=force)
76 if csm and cameras:
77 cam = Path(cameras[k])
78 cubes.append(cube)
79 if cam:
80 cams.append(cam)
82 ba = Path(bundle_adjust_prefix) if bundle_adjust_prefix else work_dir / "ba" / "run"
83 if not bundle_adjust_prefix and not reuse_existing(Path(f"{ba}-final_residuals_pointmap.csv"), force):
84 asp.bundle_adjust(cubes, ba, cams, threads=threads, **PAIR_BUNDLE_ADJUST_DEFAULTS)
86 opts = {
87 **STEREO_DEFAULTS,
88 "processes": processes,
89 "threads": threads,
90 "bundle_adjust_prefix": ba,
91 **(stereo_options or {}),
92 }
93 out = work_dir / "stereo" / "run"
94 pc = Path(f"{out}-PC.tif")
95 pair_cams: list[Path | None] = [*cams] if cams else [None, None]
96 if not reuse_existing(pc, force): 96 ↛ 125line 96 didn't jump to line 125 because the condition on line 96 was always true
97 if prior_dem:
98 maps = [Path(m) for m in mapprojected] if mapprojected else []
99 if not maps:
100 for cube, camera in zip(cubes, pair_cams):
101 m = work_dir / "prj" / f"{cube.name.split('.')[0]}.map.tif"
102 if not reuse_existing(m, force): 102 ↛ 112line 102 didn't jump to line 112 because the condition on line 102 was always true
103 asp.mapproject(
104 prior_dem,
105 cube,
106 m,
107 camera,
108 tr=image_resolution,
109 bundle_adjust_prefix=ba,
110 threads=threads,
111 )
112 maps.append(m)
113 asp.stereo(
114 maps[0],
115 maps[1],
116 out,
117 *(cubes if not cams else cams),
118 dem=prior_dem,
119 alignment_method="none",
120 **opts,
121 )
122 else:
123 asp.stereo(cubes[0], cubes[1], out, *pair_cams, alignment_method="local_epipolar", **opts)
125 dem = Path(f"{out}-DEM.tif")
126 if not reuse_existing(dem, force): 126 ↛ 128line 126 didn't jump to line 128 because the condition on line 126 was always true
127 asp.point2dem(pc, tr=resolution, datum="D_MOON", errorimage=True)
128 if not (prior_dem and align):
129 return dem
131 align_prefix = work_dir / "align" / "run"
132 trans = Path(f"{align_prefix}-trans_source.tif")
133 if not reuse_existing(trans, force): 133 ↛ 142line 133 didn't jump to line 142 because the condition on line 133 was always true
134 asp.pc_align(
135 prior_dem,
136 pc,
137 align_prefix,
138 max_displacement=max_displacement,
139 save_transformed_source_points=True,
140 highest_accuracy=True,
141 )
142 aligned = Path(f"{align_prefix}-trans_source-DEM.tif")
143 if not reuse_existing(aligned, force): 143 ↛ 145line 143 didn't jump to line 145 because the condition on line 143 was always true
144 asp.point2dem(trans, tr=resolution, datum="D_MOON")
145 return aligned