Coverage for lunar_toolkit/stereo.py: 96%
222 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 02:58 +0000
1"""Stereo pair ranking by convergence, illumination similarity and overlap."""
3from __future__ import annotations
5import math
6from typing import TYPE_CHECKING, Any
8import numpy as np
10from .constants import MOON_RADIUS_KM, parse_nac_id
11from .geometry import AoiLike, Ring, aoi_geometry, footprint_geometry, footprint_rings, load_aoi
13if TYPE_CHECKING: # pragma: no cover
14 import pandas as pd
16DEFAULT_WEIGHTS: dict[str, float] = {
17 "convergence": 0.35,
18 "overlap": 0.25,
19 "incidence": 0.20,
20 "azimuth": 0.10,
21 "resolution": 0.10,
22}
24PAIR_COLUMNS: list[str] = [
25 "index_1",
26 "index_2",
27 "product_id_1",
28 "product_id_2",
29 "instrument",
30 "start_time_1",
31 "start_time_2",
32 "dt_days",
33 "convergence",
34 "convergence_min",
35 "convergence_max",
36 "b_h",
37 "incidence_1",
38 "incidence_2",
39 "incidence_diff",
40 "azimuth_1",
41 "azimuth_2",
42 "azimuth_diff",
43 "emission_1",
44 "emission_2",
45 "pixel_scale_1",
46 "pixel_scale_2",
47 "resolution_ratio",
48 "overlap_km2",
49 "overlap_frac",
50 "aoi_coverage",
51 "overlap_lat",
52 "overlap_lon",
53 "score",
54]
57def unit_vector(lat, lon) -> np.ndarray:
58 """Unit vector(s) ``(..., 3)`` of planetocentric latitude / longitude (deg)."""
59 la, lo = np.radians(np.asarray(lat, float)), np.radians(np.asarray(lon, float))
60 return np.stack([np.cos(la) * np.cos(lo), np.cos(la) * np.sin(lo), np.sin(la)], axis=-1)
63def _latlon(v: np.ndarray) -> tuple[float, float]:
64 v = v / np.linalg.norm(v)
65 return math.degrees(math.asin(max(-1.0, min(1.0, v[2])))), math.degrees(math.atan2(v[1], v[0])) % 360.0
68def angle_between(a: np.ndarray, b: np.ndarray) -> np.ndarray:
69 """Angle (deg) between vectors along the last axis (numerically stable)."""
70 cross = np.linalg.norm(np.cross(a, b), axis=-1)
71 return np.degrees(np.arctan2(cross, np.sum(a * b, axis=-1)))
74def azimuth_difference(a, b) -> np.ndarray:
75 """Smallest circular difference (deg, 0–180) between azimuths."""
76 d = np.abs((np.asarray(a, float) - np.asarray(b, float)) % 360.0)
77 return np.minimum(d, 360.0 - d)
80def _rotate_onto(v: np.ndarray, a: np.ndarray, b: np.ndarray) -> np.ndarray:
81 """Apply to *v* the minimal rotation taking unit vector *a* onto *b*."""
82 k = np.cross(a, b)
83 s, c = np.linalg.norm(k), float(np.dot(a, b))
84 if s < 1e-15:
85 return v
86 k = k / s
87 return v * c + np.cross(k, v) * s + k * np.dot(k, v) * (1.0 - c)
90def convergence_angle(
91 ground: tuple[float, float],
92 sc1: tuple[float, float, float],
93 sc2: tuple[float, float, float],
94 centre1: tuple[float, float] | None = None,
95 centre2: tuple[float, float] | None = None,
96 radius_km: float = MOON_RADIUS_KM,
97) -> tuple[float, float]:
98 """Convergence angle and base/height ratio for two views of one point.
100 Parameters
101 ----------
102 ground : (lat, lon)
103 Ground point (deg) on the sphere of radius *radius_km*.
104 sc1, sc2 : (sub_spacecraft_lat, sub_spacecraft_lon, altitude_km)
105 Spacecraft positions as sub-spacecraft point + altitude above the
106 sphere, at the time the image centre was acquired.
107 centre1, centre2 : (lat, lon), optional
108 Image centres. Pushbroom images (NAC, ShadowCam) are tens of km
109 long, so the spacecraft was elsewhere when it imaged a ground point
110 away from the centre. When given, each spacecraft position is moved
111 by the rotation that takes the image centre onto *ground*, i.e. the
112 look direction is held fixed in the local frame (constant pointing
113 along a circular orbit).
115 Returns
116 -------
117 (convergence_deg, b_h)
118 Angle between the look vectors spacecraft → ground, and the baseline
119 ``|S1 - S2|`` divided by the mean slant range to the ground point.
120 """
121 gu = unit_vector(*ground)
122 g = radius_km * gu
123 s1 = (radius_km + sc1[2]) * unit_vector(sc1[0], sc1[1])
124 s2 = (radius_km + sc2[2]) * unit_vector(sc2[0], sc2[1])
125 if centre1 is not None:
126 s1 = _rotate_onto(s1, unit_vector(*centre1), gu)
127 if centre2 is not None:
128 s2 = _rotate_onto(s2, unit_vector(*centre2), gu)
129 conv = float(angle_between(g - s1, g - s2))
130 h = 0.5 * (np.linalg.norm(s1 - g) + np.linalg.norm(s2 - g))
131 return conv, float(np.linalg.norm(s1 - s2) / h)
134def _rotate_onto_many(v: np.ndarray, a: np.ndarray, b: np.ndarray) -> np.ndarray:
135 """Row-wise `_rotate_onto` of ``(n, 3)`` arrays."""
136 k = np.cross(a, b)
137 s = np.linalg.norm(k, axis=-1)
138 c = np.sum(a * b, axis=-1)
139 still = s < 1e-15
140 k = k / np.where(still, 1.0, s)[:, None]
141 out = (
142 v * c[:, None]
143 + np.cross(k, v) * s[:, None]
144 + k * np.sum(k * v, axis=-1)[:, None] * (1.0 - c)[:, None]
145 )
146 return np.where(still[:, None], v, out)
149def convergence_angles(
150 ground: tuple[Any, Any],
151 sc1: tuple[Any, Any, Any],
152 sc2: tuple[Any, Any, Any],
153 centre1: tuple[Any, Any] | None = None,
154 centre2: tuple[Any, Any] | None = None,
155 radius_km: float = MOON_RADIUS_KM,
156) -> tuple[np.ndarray, np.ndarray]:
157 """[`convergence_angle`][lunar_toolkit.stereo.convergence_angle] for arrays: each tuple element is an
158 array over pairs.
160 Returns
161 -------
162 (convergence_deg, b_h) : ndarrays
163 NaN where an input is NaN.
164 """
165 gu = unit_vector(*ground).reshape(-1, 3)
166 g = radius_km * gu
167 s1 = (radius_km + np.asarray(sc1[2], float)).reshape(-1, 1) * unit_vector(sc1[0], sc1[1]).reshape(-1, 3)
168 s2 = (radius_km + np.asarray(sc2[2], float)).reshape(-1, 1) * unit_vector(sc2[0], sc2[1]).reshape(-1, 3)
169 if centre1 is not None: 169 ↛ 171line 169 didn't jump to line 171 because the condition on line 169 was always true
170 s1 = _rotate_onto_many(s1, unit_vector(*centre1).reshape(-1, 3), gu)
171 if centre2 is not None: 171 ↛ 173line 171 didn't jump to line 173 because the condition on line 171 was always true
172 s2 = _rotate_onto_many(s2, unit_vector(*centre2).reshape(-1, 3), gu)
173 conv = angle_between(g - s1, g - s2)
174 h = 0.5 * (np.linalg.norm(s1 - g, axis=-1) + np.linalg.norm(s2 - g, axis=-1))
175 return conv, np.linalg.norm(s1 - s2, axis=-1) / h
178class TangentPlane:
179 """Gnomonic projection (km) about a centre point; great circles → lines."""
181 def __init__(self, lat: float, lon: float, radius_km: float = MOON_RADIUS_KM) -> None:
182 self.c = unit_vector(lat, lon)
183 ref = np.array([0.0, 0.0, 1.0]) if abs(self.c[2]) < 0.999 else np.array([1.0, 0.0, 0.0])
184 self.e = np.cross(ref, self.c)
185 self.e /= np.linalg.norm(self.e)
186 self.n = np.cross(self.c, self.e)
187 self.r = radius_km
189 @classmethod
190 def from_vector(cls, v: np.ndarray) -> TangentPlane:
191 return cls(*_latlon(np.asarray(v, float)))
193 def ring(self, ring: Ring) -> list[tuple[float, float]]:
194 """Project a lon/lat ring (deg) to plane (x, y) km (a `Projection`)."""
195 a = np.asarray(ring, float).reshape(-1, 2)
196 v = unit_vector(a[:, 1], a[:, 0])
197 d = np.maximum(v @ self.c, 0.05) # far side (callers prefilter): clamp
198 return list(zip((self.r * (v @ self.e) / d).tolist(), (self.r * (v @ self.n) / d).tolist()))
200 def unproject(self, x: float, y: float) -> tuple[float, float]:
201 return _latlon(self.c + (x / self.r) * self.e + (y / self.r) * self.n)
204def observation_footprints(df: pd.DataFrame) -> tuple[list, np.ndarray, np.ndarray]:
205 """Footprint rings, centre unit vectors and angular radii (deg) of observation rows.
207 Rings are [`lunar_toolkit.geometry.footprint_rings`][lunar_toolkit.geometry.footprint_rings] of each row
208 (``footprint_wkt``, else the ``min_lat``…``lon_east`` box), ``None``
209 (centre 0, radius 0) where a row has no footprint polygon.
210 """
211 rings = [footprint_rings(r, "footprint_wkt") for r in df.to_dict("records")]
212 centres, radii = np.zeros((len(rings), 3)), np.zeros(len(rings))
213 for k, rr in enumerate(rings):
214 if rr is None or len(rr[0][0]) < 3: 214 ↛ 215line 214 didn't jump to line 215 because the condition on line 214 was never true
215 rings[k] = None
216 continue
217 centres[k], radii[k] = _cap([poly[0] for poly in rr])
218 return rings, centres, radii
221def _cap(rings: list) -> tuple[np.ndarray, float]:
222 """Centre unit vector and angular radius (deg) of lon/lat rings' vertices."""
223 ll = np.array([p for r in rings for p in r], float)
224 v = unit_vector(ll[:, 1], ll[:, 0])
225 c = v.sum(axis=0)
226 c /= np.linalg.norm(c) or 1.0
227 return c, float(np.max(angle_between(v, c[None, :])))
230def aoi_plane(aoi: AoiLike) -> tuple[TangentPlane, Any, np.ndarray, float]:
231 """Tangent plane about an AOI, the AOI in it (km), and its centre / angular radius.
233 *aoi* is anything [`lunar_toolkit.geometry.load_aoi`][lunar_toolkit.geometry.load_aoi] accepts,
234 e.g. ``{"point": (lat, lon), "radius_km": r}``,
235 ``{"lat": (s, n), "lon": (w, e)}``, a GeoJSON / WKT path or string.
237 Returns
238 -------
239 (TangentPlane, shapely geometry, ndarray, float)
240 """
241 a = load_aoi(aoi)
242 centre, radius = _cap([poly[0] for poly in a.polygons])
243 plane = TangentPlane.from_vector(centre)
244 return plane, aoi_geometry(a, plane.ring), centre, radius
247def _column(df: pd.DataFrame, col: str) -> np.ndarray:
248 import pandas as pd
250 if col not in df.columns: 250 ↛ 251line 250 didn't jump to line 251 because the condition on line 250 was never true
251 return np.full(len(df), np.nan)
252 return pd.to_numeric(df[col], errors="coerce").to_numpy(float)
255def _ramp(x: np.ndarray, lo: float, ideal_lo: float, ideal_hi: float, hi: float) -> np.ndarray:
256 """1 inside [ideal_lo, ideal_hi], linear to 0 at lo / hi."""
257 up = np.clip((x - lo) / max(ideal_lo - lo, 1e-9), 0.0, 1.0)
258 down = np.clip((hi - x) / max(hi - ideal_hi, 1e-9), 0.0, 1.0)
259 return np.minimum(up, down)
262def find_stereo_pairs(
263 obs_df: pd.DataFrame,
264 aoi: AoiLike | None = None,
265 convergence: tuple[float, float] = (10.0, 60.0),
266 ideal_convergence: tuple[float, float] = (15.0, 45.0),
267 max_incidence_diff: float = 20.0,
268 max_azimuth_diff: float = 30.0,
269 min_overlap: float = 0.1,
270 min_aoi_coverage: float = 0.0,
271 max_resolution_ratio: float | None = None,
272 same_instrument: bool = True,
273 require_geometry: bool = False,
274 min_time_sep_s: float = 60.0,
275 weights: dict[str, float] | None = None,
276) -> pd.DataFrame:
277 """Rank candidate stereo pairs in an observation table.
279 Parameters
280 ----------
281 obs_df : DataFrame
282 Observations in the common schema (extra columns are ignored).
283 aoi : dict, str, Path or Aoi, optional
284 Area of interest: ``{"point": (lat, lon), "radius_km": r}``,
285 ``{"lat": (s, n), "lon": (w, e)}`` or anything
286 [`lunar_toolkit.geometry.load_aoi`][lunar_toolkit.geometry.load_aoi] accepts (GeoJSON / WKT file
287 or string). Overlaps are then clipped to the AOI, pairs not
288 overlapping inside it are dropped and the convergence angle is
289 evaluated at the centroid of the clipped overlap.
290 convergence : (min, max)
291 Hard limits on the convergence angle (deg). Pairs without spacecraft
292 geometry are kept only if the emission bounds ``[|e1-e2|, e1+e2]``
293 (from centre emission angles, so approximate; 5° slack is allowed)
294 intersect this range and ``require_geometry`` is False.
295 ideal_convergence : (min, max)
296 Range scored 1; the score falls linearly to 0 at *convergence* limits.
297 max_incidence_diff, max_azimuth_diff : float
298 Hard limits (deg) on ``|i1 - i2|`` and the circular sub-solar azimuth
299 difference (azimuth is ignored where unknown).
300 min_overlap : float
301 Minimum overlap area as a fraction of the smaller footprint (with
302 an *aoi*: of the smaller footprint-inside-AOI).
303 min_aoi_coverage : float
304 Minimum fraction of the AOI covered by the overlap (with *aoi*).
305 max_resolution_ratio : float, optional
306 Maximum ``pixel_scale`` ratio (coarser / finer).
307 same_instrument : bool
308 Only pair images from the same ``instrument``.
309 require_geometry : bool
310 Drop pairs whose convergence angle cannot be computed.
311 min_time_sep_s : float
312 Minimum start-time separation; also, NAC LE/RE halves of one
313 observation are never paired (their convergence is ~0).
314 weights : dict, optional
315 Score weights, keys as `DEFAULT_WEIGHTS`.
317 Returns
318 -------
319 DataFrame
320 One row per pair, best first, with ``index_1`` / ``index_2`` (row
321 labels in *obs_df*), ``product_id_1`` / ``_2``, ``instrument``,
322 ``start_time_1`` / ``_2``, ``dt_days``, ``convergence`` (deg, NaN if
323 unknown), ``convergence_min`` / ``convergence_max`` (emission bounds),
324 ``b_h`` (baseline / mean slant range, ≈ 2 sin(convergence / 2)),
325 ``incidence_1`` / ``_2``, ``incidence_diff``, ``azimuth_1`` / ``_2``, ``azimuth_diff``,
326 ``emission_1`` / ``_2``, ``pixel_scale_1`` / ``_2``,
327 ``resolution_ratio``, ``overlap_km2``, ``overlap_frac``,
328 ``aoi_coverage`` (NaN without *aoi*), ``overlap_lat`` /
329 ``overlap_lon`` (where convergence is evaluated) and ``score``.
331 Notes
332 -----
333 ``score`` is a weighted mean of terms in [0, 1]:
335 * convergence — 1 inside *ideal_convergence*, linear to 0 at the
336 *convergence* limits (0.5 when the angle is unknown);
337 * overlap — ``aoi_coverage`` with an AOI, else ``overlap_frac``;
338 * incidence — ``1 - incidence_diff / max_incidence_diff``;
339 * azimuth — ``1 - azimuth_diff / max_azimuth_diff`` (0.5 if unknown);
340 * resolution — finer / coarser ``pixel_scale`` (0.5 if unknown).
341 """
342 import pandas as pd
344 cols = PAIR_COLUMNS
345 if obs_df is None or len(obs_df) < 2: 345 ↛ 346line 345 didn't jump to line 346 because the condition on line 345 was never true
346 return pd.DataFrame(columns=cols)
348 df = obs_df
349 labels = list(df.index)
350 n = len(df)
351 pids = (
352 df["product_id"].astype(str).to_numpy() if "product_id" in df else np.array([str(i) for i in labels])
353 )
354 keys = np.array([(parse_nac_id(p) or (p.upper(),))[0] for p in pids]) # NAC L/R halves share one
355 inst = df["instrument"].astype(str).to_numpy() if "instrument" in df else np.full(n, "")
356 t = pd.to_datetime(df["start_time"], errors="coerce") if "start_time" in df else pd.Series([pd.NaT] * n)
357 tsec = (t - pd.Timestamp("2000-01-01")).dt.total_seconds().to_numpy(float)
358 inc, emi = _column(df, "incidence"), _column(df, "emission")
359 az, ps = _column(df, "sub_solar_azimuth"), _column(df, "pixel_scale")
360 sclat, sclon = _column(df, "sub_spacecraft_lat"), _column(df, "sub_spacecraft_lon")
361 alt = _column(df, "spacecraft_altitude")
362 has_geom = ~(np.isnan(sclat) | np.isnan(sclon) | np.isnan(alt))
364 rings, centres, radii = observation_footprints(df)
365 valid = np.array([r is not None for r in rings], bool)
366 clat, clon = _column(df, "center_lat"), _column(df, "center_lon")
367 img_centres = centres.copy()
368 known = ~(np.isnan(clat) | np.isnan(clon))
369 img_centres[known] = unit_vector(clat[known], clon[known])
371 plane = None
372 clipped: list = []
373 aoi_area = np.nan
374 if aoi is not None:
375 plane, target, ac, aoi_radius = aoi_plane(aoi)
376 aoi_area = target.area
377 near = valid & (angle_between(centres, ac[None, :]) < radii + aoi_radius + 1e-6)
378 clipped = [
379 footprint_geometry(r, plane.ring).intersection(target) if near[k] else None
380 for k, r in enumerate(rings)
381 ]
382 valid &= np.array([g is not None and g.area > 0 for g in clipped], bool)
384 out = []
385 lo_c, hi_c = convergence
386 for i in range(n - 1):
387 if not valid[i]:
388 continue
389 j = np.arange(i + 1, n)
390 m = valid[j] & (keys[j] != keys[i])
391 if same_instrument: 391 ↛ 393line 391 didn't jump to line 393 because the condition on line 391 was always true
392 m &= inst[j] == inst[i]
393 dt = np.abs(tsec[j] - tsec[i])
394 m &= ~(dt < min_time_sep_s)
395 dinc = np.abs(inc[j] - inc[i])
396 m &= ~(dinc > max_incidence_diff)
397 daz = azimuth_difference(az[j], az[i])
398 m &= ~(daz > max_azimuth_diff)
399 cmin, cmax = np.abs(emi[j] - emi[i]), emi[j] + emi[i]
400 # Centre emissions only bound the convergence approximately: allow slack.
401 m &= ~(cmax < lo_c - 5.0) & ~(cmin > hi_c + 5.0)
402 if require_geometry:
403 m &= has_geom[j] & has_geom[i]
404 if max_resolution_ratio is not None:
405 rr = np.maximum(ps[j], ps[i]) / np.minimum(ps[j], ps[i])
406 m &= ~(rr > max_resolution_ratio)
407 sep = angle_between(centres[j], centres[i][None, :])
408 m &= sep <= radii[j] + radii[i] + 1e-6
409 for jj in j[m]:
410 rec = _pair_record(
411 i, jj, rings, centres, img_centres, plane, clipped, aoi_area, has_geom, sclat, sclon, alt
412 )
413 if rec is None:
414 continue
415 frac = rec["overlap_frac"]
416 if frac < min_overlap or (aoi is not None and rec["aoi_coverage"] < max(min_aoi_coverage, 1e-9)):
417 continue
418 c = rec["convergence"]
419 if not math.isnan(c) and not (lo_c <= c <= hi_c):
420 continue
421 rec.update(i=i, j=int(jj))
422 out.append(rec)
424 if not out:
425 return pd.DataFrame(columns=cols)
427 res = pd.DataFrame(out)
428 i_, j_ = res.pop("i").to_numpy(), res.pop("j").to_numpy()
429 res.insert(0, "index_1", [labels[k] for k in i_])
430 res.insert(1, "index_2", [labels[k] for k in j_])
431 res["product_id_1"], res["product_id_2"] = pids[i_], pids[j_]
432 res["instrument"] = inst[i_]
433 res["start_time_1"], res["start_time_2"] = t.to_numpy()[i_], t.to_numpy()[j_]
434 res["dt_days"] = np.abs(tsec[j_] - tsec[i_]) / 86400.0
435 res["convergence_min"] = np.abs(emi[j_] - emi[i_])
436 res["convergence_max"] = emi[j_] + emi[i_]
437 res["incidence_1"], res["incidence_2"] = inc[i_], inc[j_]
438 res["incidence_diff"] = np.abs(inc[j_] - inc[i_])
439 res["azimuth_1"], res["azimuth_2"] = az[i_], az[j_]
440 res["azimuth_diff"] = azimuth_difference(az[j_], az[i_])
441 res["emission_1"], res["emission_2"] = emi[i_], emi[j_]
442 res["pixel_scale_1"], res["pixel_scale_2"] = ps[i_], ps[j_]
443 res["resolution_ratio"] = np.maximum(ps[i_], ps[j_]) / np.minimum(ps[i_], ps[j_])
445 res["score"] = score_pairs(
446 res,
447 convergence=convergence,
448 ideal_convergence=ideal_convergence,
449 max_incidence_diff=max_incidence_diff,
450 max_azimuth_diff=max_azimuth_diff,
451 weights=weights,
452 aoi=aoi is not None,
453 )
454 return res[cols].sort_values("score", ascending=False, kind="stable").reset_index(drop=True)
457def score_pairs(
458 pairs: pd.DataFrame,
459 *,
460 convergence: tuple[float, float] = (10.0, 60.0),
461 ideal_convergence: tuple[float, float] = (15.0, 45.0),
462 max_incidence_diff: float = 20.0,
463 max_azimuth_diff: float = 30.0,
464 weights: dict[str, float] | None = None,
465 aoi: bool = False,
466) -> np.ndarray:
467 """The [`find_stereo_pairs`][lunar_toolkit.stereo.find_stereo_pairs] ``score`` of pair rows
468 (``convergence``, ``overlap_frac`` or ``aoi_coverage`` with *aoi*, ``incidence_diff``, ``azimuth_diff``,
469 ``resolution_ratio``)."""
470 w = dict(DEFAULT_WEIGHTS, **(weights or {}))
471 conv = pairs["convergence"].to_numpy(float)
472 s_conv = np.where(np.isnan(conv), 0.5, _ramp(conv, convergence[0], *ideal_convergence, convergence[1]))
473 s_ovl = pairs["aoi_coverage" if aoi else "overlap_frac"].to_numpy(float)
474 s_inc = np.nan_to_num(
475 1.0 - pairs["incidence_diff"].to_numpy(float) / max(max_incidence_diff, 1e-9), nan=0.5
476 )
477 s_az = np.nan_to_num(1.0 - pairs["azimuth_diff"].to_numpy(float) / max(max_azimuth_diff, 1e-9), nan=0.5)
478 s_res = np.nan_to_num(1.0 / pairs["resolution_ratio"].to_numpy(float), nan=0.5)
479 tot = sum(w.values()) or 1.0
480 return (
481 w["convergence"] * s_conv
482 + w["overlap"] * np.clip(s_ovl, 0, 1)
483 + w["incidence"] * np.clip(s_inc, 0, 1)
484 + w["azimuth"] * np.clip(s_az, 0, 1)
485 + w["resolution"] * s_res
486 ) / tot
489def _pair_record(
490 i, j, rings, centres, img_centres, plane, clipped, aoi_area, has_geom, sclat, sclon, alt
491) -> dict | None:
492 if plane is None:
493 pl = TangentPlane.from_vector(centres[i] + centres[j])
494 g1, g2 = footprint_geometry(rings[i], pl.ring), footprint_geometry(rings[j], pl.ring)
495 else:
496 pl, g1, g2 = plane, clipped[i], clipped[j] # already clipped to the AOI
497 inter = g1.intersection(g2)
498 area = inter.area
499 if area <= 0:
500 return None
501 aoi_cov = area / aoi_area if plane is not None and aoi_area > 0 else np.nan
502 c = inter.centroid
503 glat, glon = pl.unproject(c.x, c.y)
504 conv = bh = np.nan
505 if has_geom[i] and has_geom[j]:
506 conv, bh = convergence_angle(
507 (glat, glon),
508 (sclat[i], sclon[i], alt[i]),
509 (sclat[j], sclon[j], alt[j]),
510 _latlon(img_centres[i]),
511 _latlon(img_centres[j]),
512 )
513 return {
514 "convergence": conv,
515 "b_h": bh,
516 "overlap_km2": area,
517 "overlap_frac": area / max(min(g1.area, g2.area), 1e-12),
518 "aoi_coverage": aoi_cov,
519 "overlap_lat": glat,
520 "overlap_lon": glon,
521 }