Coverage for lunar_toolkit/geometry.py: 92%
351 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 02:58 +0000
1"""Search regions and AOI polygons (lon/lat, 0–360 east; polar stereographic past ±60°)."""
3from __future__ import annotations
5import json
6import math
7import re
8from collections.abc import Callable
9from dataclasses import dataclass
10from pathlib import Path
11from typing import TYPE_CHECKING, TypeAlias
13import shapely.geometry
14import shapely.wkt
15from shapely.affinity import translate
16from shapely.geometry import MultiPoint, Polygon
17from shapely.ops import unary_union
18from shapely.prepared import prep
19from shapely.validation import make_valid
21from .constants import FOOTPRINT_CORNERS, MOON_RADIUS_KM, POLAR_LAT_THRESHOLD
23if TYPE_CHECKING: # pragma: no cover
24 import pandas as pd
27def ps_xy(lat: float, lon: float, hemi: str) -> tuple[float, float]:
28 """Polar stereographic (x, y) in metres; mirrors the SQL ``_ps_*`` macros."""
29 sgn = 1.0 if hemi == "S" else -1.0
30 rho = 2.0 * MOON_RADIUS_KM * 1000.0 * math.tan(math.pi / 4 + sgn * math.radians(lat) / 2)
31 return rho * math.sin(math.radians(lon)), sgn * rho * math.cos(math.radians(lon))
34def _ps_hemis(s: float, n: float) -> list[str]:
35 """Hemispheres (``'S'`` / ``'N'``) a latitude range reaches poleward of ±60°."""
36 return [h for h, past in (("S", s < -POLAR_LAT_THRESHOLD), ("N", n > POLAR_LAT_THRESHOLD)) if past]
39def query_region(
40 lat: tuple[float, float] | None = None,
41 lon: tuple[float, float] | None = None,
42 point: tuple[float, float] | None = None,
43 radius_km: float = 0.0,
44) -> dict:
45 """Normalise a bbox or point+radius into lat range, lon arc and PS boxes."""
46 if point is not None:
47 plat, plon = float(point[0]), float(point[1])
48 delta = math.degrees(radius_km / MOON_RADIUS_KM)
49 s, n = plat - delta, plat + delta
50 if n >= 90 or s <= -90 or radius_km <= 0 and abs(plat) >= 90:
51 west, span = 0.0, 360.0
52 else:
53 ratio = math.sin(math.radians(delta)) / math.cos(math.radians(plat))
54 dlon = 90.0 if ratio >= 1 else math.degrees(math.asin(ratio))
55 west, span = (plon - dlon) % 360.0, min(2 * dlon, 360.0)
56 s, n = max(s, -90.0), min(n, 90.0)
57 ps: dict[str, tuple[float, float, float, float]] = {}
58 for hemi in _ps_hemis(s, n):
59 x, y = ps_xy(plat, plon, hemi)
60 # PS scale k = 2 / (1 + sin|φ|); use the largest k in the circle
61 edge = min(abs(s), abs(n)) if s * n > 0 else 0.0
62 k = 2.0 / (1.0 + math.sin(math.radians(edge)))
63 r = radius_km * 1000.0 * k * 1.01 + 1.0
64 ps[hemi] = (x - r, x + r, y - r, y + r)
65 return {"lat": (s, n), "west": west, "span": span, "ps": ps}
67 s, n = (float(lat[0]), float(lat[1])) if lat is not None else (-90.0, 90.0)
68 if s > n: 68 ↛ 69line 68 didn't jump to line 69 because the condition on line 68 was never true
69 s, n = n, s
70 if lon is None or float(lon[1]) - float(lon[0]) >= 360.0:
71 west, span = 0.0, 360.0
72 else:
73 west = float(lon[0]) % 360.0
74 span = (float(lon[1]) - float(lon[0])) % 360.0
75 ps = {hemi: _ps_bbox_of_box(s, n, west, span, hemi) for hemi in _ps_hemis(s, n)}
76 return {"lat": (s, n), "west": west, "span": span, "ps": ps}
79def _ps_bbox_of_box(
80 s: float, n: float, west: float, span: float, hemi: str, samples: int = 181
81) -> tuple[float, float, float, float]:
82 """Projected bbox of a lat/lon box by sampling its boundary."""
83 # Only the part of the box on this hemisphere's side of ±60° matters for
84 # rows that carry PS coordinates, but clipping at the equator is enough.
85 if hemi == "S":
86 s2, n2 = s, min(n, 0.0)
87 else:
88 s2, n2 = max(s, 0.0), n
89 pts = []
90 for i in range(samples):
91 f = i / (samples - 1)
92 lo = west + span * f
93 la = s2 + (n2 - s2) * f
94 pts += [(s2, lo), (n2, lo), (la, west), (la, west + span)]
95 xy = [ps_xy(a, b, hemi) for a, b in pts]
96 xs, ys = [p[0] for p in xy], [p[1] for p in xy]
97 pole_lat = -90.0 if hemi == "S" else 90.0
98 if span >= 360.0 and s2 <= pole_lat <= n2:
99 xs.append(0.0)
100 ys.append(0.0)
101 pad = 0.005 * max(max(xs) - min(xs), max(ys) - min(ys), 1.0)
102 return min(xs) - pad, max(xs) + pad, min(ys) - pad, max(ys) + pad
105def region_bbox(region: dict) -> tuple[float, float, float, float]:
106 """``(west, south, east, north)`` of a [`query_region`][lunar_toolkit.geometry.query_region] (east may be
107 < west)."""
108 s, n = region["lat"]
109 if region["span"] >= 360.0:
110 return 0.0, s, 360.0, n
111 west = region["west"]
112 east = (west + region["span"]) % 360.0
113 return west, s, east, n
116Ring = list[tuple[float, float]] # open ring of (lon, lat) vertices, degrees
117AoiLike: TypeAlias = "Aoi | dict | str | Path" # anything load_aoi accepts
118Projection = Callable[[Ring], list[tuple[float, float]]] # lon/lat ring -> planar (x, y) ring
120_WKT_COLUMNS = ("FP_WKT", "footprint_wkt")
123@dataclass
124class Aoi:
125 """Area of interest: one or more lon/lat polygons (degrees, positive east).
127 ``polygons`` holds, per polygon, its exterior ring followed by any holes.
128 Rings are open (no repeated closing vertex), densified and longitude
129 unwrapped so consecutive vertices never jump across the 0/360 seam
130 (values may lie outside 0..360). A ring whose longitudes wind once
131 around the globe encloses a pole.
132 """
134 polygons: list[list[Ring]]
135 source: str = ""
137 @property
138 def wkt(self) -> str:
139 """``MULTIPOLYGON`` WKT of the (unwrapped) rings."""
141 def ring_wkt(r: Ring) -> str:
142 return "(" + ", ".join(f"{lo:.6f} {la:.6f}" for lo, la in r + r[:1]) + ")"
144 return (
145 "MULTIPOLYGON("
146 + ", ".join("(" + ", ".join(ring_wkt(r) for r in poly) + ")" for poly in self.polygons)
147 + ")"
148 )
151_WKT_RE = re.compile(
152 r"^\s*((MULTI)?(POINT|LINESTRING|POLYGON)|GEOMETRYCOLLECTION)\s*(Z|M|ZM)?\s*(\(|EMPTY)", re.IGNORECASE
153)
156def _is_file(text: str) -> bool:
157 """True when *text* names an existing file (WKT / GeoJSON strings never do)."""
158 if not text or "\n" in text or text.startswith("{") or _WKT_RE.match(text):
159 return False
160 try:
161 return Path(text).is_file()
162 except (OSError, ValueError):
163 return False
166def load_aoi(source: AoiLike, densify_deg: float | None = 1.0) -> Aoi:
167 """Load an area of interest from GeoJSON, WKT or a box / circle region.
169 Parameters
170 ----------
171 source : str, Path, dict or Aoi
172 A ``.geojson`` / ``.json`` file (Feature, FeatureCollection or bare
173 Geometry), a ``.wkt`` (or any text) file, a raw WKT string, a GeoJSON
174 mapping, a region dict ``{"point": (lat, lon), "radius_km": r}``
175 (spherical cap) or ``{"lat": (south, north), "lon": (west, east)}``
176 (``west > east`` crosses 0/360; omit ``lon`` for all longitudes), or
177 an existing [`Aoi`][lunar_toolkit.geometry.Aoi] (returned unchanged). Only
178 Polygon / MultiPolygon parts are used (also inside
179 GeometryCollections); coordinates are lon/lat degrees, positive east,
180 in either −180..180 or 0..360. Polygons may cross the 0/360 (or
181 ±180) seam or wind around a pole.
182 densify_deg : float, optional
183 Insert vertices so no edge spans more than this many degrees of
184 latitude or longitude. Edges are taken as straight in lon/lat (as in
185 GeoJSON); densifying keeps them that way once projected to polar
186 stereographic. ``None`` keeps the input vertices.
188 Returns
189 -------
190 Aoi
191 """
192 if isinstance(source, Aoi):
193 return source
194 label = ""
195 if isinstance(source, dict) and "type" not in source:
196 raw = [[_region_ring(source)]]
197 label = repr(source)
198 elif isinstance(source, dict):
199 raw = _shapely_rings(_geojson_geometries(source))
200 else:
201 text = str(source).strip()
202 if _is_file(text):
203 label = text
204 text = Path(text).read_text(encoding="utf-8").strip()
205 if text.startswith("{"):
206 geoms = _geojson_geometries(json.loads(text))
207 else:
208 try:
209 geoms = [shapely.wkt.loads(text)]
210 except Exception as exc:
211 raise ValueError(f"Not a GeoJSON/WKT file or WKT string: {text[:80]!r}") from exc
212 raw = _shapely_rings(geoms)
213 label = label or text[:80]
215 polygons: list[list[Ring]] = []
216 for rings in raw:
217 rings = [r for r in (_open_ring(r) for r in rings) if len(r) >= 3]
218 if not rings: 218 ↛ 219line 218 didn't jump to line 219 because the condition on line 218 was never true
219 continue
220 if densify_deg:
221 rings = [_densify_ring(r, densify_deg) for r in rings]
222 rings = [_unwrap_ring(r) for r in rings]
223 out = _shift_to_0_360(rings[:1])
224 ext_mid = _lon_mid(out[0])
225 for h in rings[1:]: 225 ↛ 226line 225 didn't jump to line 226 because the loop on line 225 never started
226 kh = round((_lon_mid(h) - ext_mid) / 360.0)
227 out.append([(lo - 360.0 * kh, la) for lo, la in h])
228 polygons.append(out)
229 if not polygons:
230 raise ValueError("AOI contains no Polygon or MultiPolygon geometry")
231 for poly in polygons:
232 for r in poly:
233 if any(abs(la) > 90.0 for _, la in r):
234 raise ValueError("AOI latitude outside ±90° — coordinates must be (lon, lat)")
235 return Aoi(polygons=polygons, source=label)
238def _region_ring(region: dict, n: int = 72) -> Ring:
239 """Boundary ring of a ``{"point", "radius_km"}`` cap or ``{"lat", "lon"}`` box."""
240 if "point" in region:
241 lat, lon = map(math.radians, region["point"])
242 d = float(region.get("radius_km", 0.0)) / MOON_RADIUS_KM
243 if d <= 0:
244 raise ValueError("point AOI needs radius_km > 0")
245 ring = []
246 for k in range(n):
247 b = 2 * math.pi * k / n
248 la = math.asin(math.sin(lat) * math.cos(d) + math.cos(lat) * math.sin(d) * math.cos(b))
249 dlo = math.atan2(
250 math.sin(b) * math.sin(d) * math.cos(lat), math.cos(d) - math.sin(lat) * math.sin(la)
251 )
252 ring.append((math.degrees(lon + dlo), math.degrees(la)))
253 return ring
254 if "lat" not in region and "lon" not in region:
255 raise ValueError(f"Region dict needs 'point' + 'radius_km' or 'lat' / 'lon', got {sorted(region)}")
256 r = query_region(lat=region.get("lat"), lon=region.get("lon"))
257 (s, n), west, span = r["lat"], r["west"], r["span"]
258 return _box_ring(s, n, west, span)
261def _shapely_rings(geoms: list) -> list[list[Ring]]:
262 """Per polygon part of *geoms*: exterior then hole rings, as (lon, lat) lists."""
263 return [
264 [[(float(p[0]), float(p[1])) for p in r.coords] for r in (poly.exterior, *poly.interiors)]
265 for g in geoms
266 for poly in _polygon_parts(g)
267 ]
270def _geojson_geometries(obj: dict) -> list:
271 kind = obj.get("type")
272 if kind == "FeatureCollection":
273 return [g for f in obj.get("features", []) for g in _geojson_geometries(f)]
274 if kind == "Feature":
275 return _geojson_geometries(obj["geometry"]) if obj.get("geometry") else []
276 if kind == "GeometryCollection": 276 ↛ 277line 276 didn't jump to line 277 because the condition on line 276 was never true
277 return [g for sub in obj.get("geometries", []) for g in _geojson_geometries(sub)]
278 return [shapely.geometry.shape(obj)]
281def _polygon_parts(geom) -> list:
282 kind = geom.geom_type
283 if kind == "Polygon":
284 return [] if geom.is_empty else [geom]
285 if kind in ("MultiPolygon", "GeometryCollection"):
286 return [p for g in geom.geoms for p in _polygon_parts(g)]
287 return []
290def _wrap180(d: float) -> float:
291 return (d + 180.0) % 360.0 - 180.0
294def _open_ring(ring: Ring) -> Ring:
295 if len(ring) > 1 and ring[0] == ring[-1]:
296 return ring[:-1]
297 return ring
300def _densify_ring(ring: Ring, step: float) -> Ring:
301 """Subdivide edges (seam-aware) so none exceeds ``step`` degrees."""
302 out: Ring = []
303 for i, (lo0, la0) in enumerate(ring):
304 lo1, la1 = ring[(i + 1) % len(ring)]
305 dlo, dla = _wrap180(lo1 - lo0), la1 - la0
306 n = max(1, math.ceil(max(abs(dlo), abs(dla)) / step))
307 out += [(lo0 + dlo * j / n, la0 + dla * j / n) for j in range(n)]
308 return out
311def _unwrap_ring(ring: Ring) -> Ring:
312 """Make longitudes continuous along the ring (no ±360 jumps)."""
313 out = [ring[0]]
314 for lo, la in ring[1:]:
315 out.append((out[-1][0] + _wrap180(lo - out[-1][0]), la))
316 return out
319def _winding(ring: Ring) -> float:
320 """Net longitude swept by the closed ring (≈ ±360 if it encloses a pole)."""
321 return sum(_wrap180(ring[(i + 1) % len(ring)][0] - ring[i][0]) for i in range(len(ring)))
324def _encloses_pole(ring: Ring) -> bool:
325 return abs(_winding(ring)) > 180.0
328def _shift_to_0_360(rings: list[Ring]) -> list[Ring]:
329 """Shift all rings by the multiple of 360° that puts the first one's west edge in [0, 360)."""
330 k = math.floor(min(lo for lo, _ in rings[0]) / 360.0)
331 return [[(lo - 360.0 * k, la) for lo, la in r] for r in rings]
334def _lon_mid(ring: Ring) -> float:
335 los = [lo for lo, _ in ring]
336 return (min(los) + max(los)) / 2.0
339def _covering_arc(arcs: list[tuple[float, float]]) -> tuple[float, float]:
340 """Smallest ``(west, span)`` arc covering all ``(west, span)`` arcs."""
341 pieces = []
342 for w, sp in arcs:
343 if sp >= 360.0:
344 return 0.0, 360.0
345 w %= 360.0
346 if w + sp > 360.0:
347 pieces += [(w, 360.0), (0.0, w + sp - 360.0)]
348 else:
349 pieces.append((w, w + sp))
350 pieces.sort()
351 merged = [list(pieces[0])]
352 for a, b in pieces[1:]:
353 if a <= merged[-1][1]:
354 merged[-1][1] = max(merged[-1][1], b)
355 else:
356 merged.append([a, b])
357 # Largest uncovered gap (including the one wrapping past 360)
358 gaps = [(merged[i + 1][0] - merged[i][1], merged[i + 1][0]) for i in range(len(merged) - 1)]
359 gaps.append((merged[0][0] + 360.0 - merged[-1][1], merged[0][0]))
360 size, east_of_gap = max(gaps)
361 if size <= 0.0: 361 ↛ 362line 361 didn't jump to line 362 because the condition on line 361 was never true
362 return 0.0, 360.0
363 return east_of_gap % 360.0, 360.0 - size
366def aoi_region(aoi: AoiLike) -> dict:
367 """Bounding region of an AOI, in [`query_region`][lunar_toolkit.geometry.query_region] form.
369 Returns ``{"lat": (s, n), "west", "span", "ps"}`` so the bbox prefilter
370 used for boxes and circles applies unchanged; refine the prefiltered
371 rows with [`filter_aoi`][lunar_toolkit.geometry.filter_aoi].
372 """
373 aoi = load_aoi(aoi)
374 lats = [la for poly in aoi.polygons for la in (p[1] for p in poly[0])]
375 s, n = min(lats), max(lats)
376 arcs = []
377 for poly in aoi.polygons:
378 ext = poly[0]
379 if _encloses_pole(ext):
380 if sum(la for _, la in ext) < 0: 380 ↛ 383line 380 didn't jump to line 383 because the condition on line 380 was always true
381 s = -90.0
382 else:
383 n = 90.0
384 arcs.append((0.0, 360.0))
385 else:
386 los = [lo for lo, _ in ext]
387 arcs.append((min(los), max(los) - min(los)))
388 west, span = _covering_arc(arcs)
389 ps = {}
390 for hemi in _ps_hemis(s, n):
391 same_side = all((la <= 0.0) if hemi == "S" else (la >= 0.0) for la in lats)
392 if not same_side:
393 ps[hemi] = _ps_bbox_of_box(s, n, west, span, hemi)
394 continue
395 xy = [ps_xy(la, lo, hemi) for poly in aoi.polygons for lo, la in poly[0]]
396 xs, ys = [p[0] for p in xy], [p[1] for p in xy]
397 pad = 0.005 * max(max(xs) - min(xs), max(ys) - min(ys)) + 1.0
398 ps[hemi] = (min(xs) - pad, max(xs) + pad, min(ys) - pad, max(ys) + pad)
399 return {"lat": (s, n), "west": west, "span": span, "ps": ps}
402def _polar_hemi(region: dict) -> str | None:
403 """``'S'`` / ``'N'`` when the whole region is poleward of ±60°."""
404 s, n = region["lat"]
405 if n < -POLAR_LAT_THRESHOLD:
406 return "S"
407 if s > POLAR_LAT_THRESHOLD: 407 ↛ 408line 407 didn't jump to line 408 because the condition on line 407 was never true
408 return "N"
409 return None
412def ps_projection(hemi: str) -> Projection:
413 """Ring projection to polar stereographic ([`ps_xy`][lunar_toolkit.geometry.ps_xy]) for hemisphere
414 *hemi*."""
415 # Clamp away from the opposite pole, where the projection diverges
416 lim = 89.0
417 return lambda ring: [ps_xy(min(la, lim) if hemi == "S" else max(la, -lim), lo, hemi) for lo, la in ring]
420def _ll_ring(ring: Ring) -> Ring:
421 """Unwrapped lon/lat ring; pole-enclosing rings are closed via the pole."""
422 r = _unwrap_ring(ring)
423 if _encloses_pole(ring): 423 ↛ 424line 423 didn't jump to line 424 because the condition on line 423 was never true
424 pole = -90.0 if sum(la for _, la in r) < 0 else 90.0
425 lo_end = r[0][0] + _winding(ring)
426 r = r + [(lo_end, r[0][1]), (lo_end, pole), (r[0][0], pole)]
427 return r
430def _make_polygon(rings: list[list[tuple[float, float]]]):
431 """Valid shapely geometry from rings (degenerate rings become hulls)."""
432 if len(rings[0]) < 3: 432 ↛ 433line 432 didn't jump to line 433 because the condition on line 432 was never true
433 return MultiPoint(rings[0]).convex_hull
434 poly = Polygon(rings[0], [h for h in rings[1:] if len(h) >= 3])
435 if not poly.is_valid: 435 ↛ 436line 435 didn't jump to line 436 because the condition on line 435 was never true
436 poly = make_valid(poly)
437 if poly.is_empty or poly.area == 0.0: 437 ↛ 438line 437 didn't jump to line 438 because the condition on line 437 was never true
438 return MultiPoint(rings[0]).convex_hull
439 return poly
442def aoi_geometry(aoi: AoiLike, proj: Projection | None = None):
443 """AOI as one shapely geometry: rings projected by *proj*, else seam-replicated lon/lat."""
444 aoi = load_aoi(aoi)
445 if proj is not None:
446 return unary_union([_make_polygon([proj(r) for r in poly]) for poly in aoi.polygons])
447 g = unary_union([_make_polygon(_shift_to_0_360([_ll_ring(r) for r in poly])) for poly in aoi.polygons])
448 return unary_union([g, translate(g, xoff=360.0), translate(g, xoff=-360.0)])
451def footprint_geometry(rings_list: list[list[Ring]], proj: Projection | None):
452 """Footprint polygons (per polygon: lon/lat exterior, *holes) as shapely geometry.
454 Rings are projected by *proj* (e.g. [`ps_projection`][lunar_toolkit.geometry.ps_projection]), else kept in
455 lon/lat (pole-enclosing rings closed via the pole) starting in 0..360.
456 """
457 parts = []
458 for rings in rings_list:
459 rings = [_open_ring(r) for r in rings]
460 if not rings or len(rings[0]) < 1: 460 ↛ 461line 460 didn't jump to line 461 because the condition on line 460 was never true
461 continue
462 if proj is not None:
463 parts.append(_make_polygon([proj(r) for r in rings]))
464 continue
465 parts.append(_make_polygon(_shift_to_0_360([_ll_ring(r) if len(r) >= 3 else r for r in rings])))
466 return unary_union(parts) if parts else None
469def _box_ring(s: float, n: float, west: float, span: float, step: float = 1.0) -> Ring:
470 """Densified lon/lat box ring (already unwrapped; span may be 360)."""
471 m = max(1, math.ceil(span / step))
472 k = max(1, math.ceil((n - s) / step))
473 bottom = [(west + span * i / m, s) for i in range(m + 1)]
474 right = [(west + span, s + (n - s) * j / k) for j in range(1, k + 1)]
475 top = [(west + span * i / m, n) for i in range(m - 1, -1, -1)]
476 left = [(west, n - (n - s) * j / k) for j in range(1, k)]
477 return bottom + right + top + left
480def _num(v) -> float | None:
481 try:
482 f = float(v)
483 except (TypeError, ValueError):
484 return None
485 return None if math.isnan(f) else f
488def footprint_rings(row: dict, wkt_column: str | None = None) -> list[list[Ring]] | None:
489 """Footprint of a result row as lon/lat rings (per polygon: exterior, *holes).
491 Source: the *wkt_column* WKT, else the four corner lat/lon columns, else
492 the ``FP_*`` / ``min_lat``…``lon_east`` bounding box; ``None`` if none.
493 """
494 wkt = row.get(wkt_column) if wkt_column else None
495 if isinstance(wkt, str) and wkt.strip():
496 try:
497 rings = _shapely_rings([shapely.wkt.loads(wkt)])
498 except Exception:
499 rings = [] # unparsable WKT: fall back to corners / bbox
500 if rings: 500 ↛ 502line 500 didn't jump to line 502 because the condition on line 500 was always true
501 return rings
502 corners = [(_num(row.get(f"{c}_LONGITUDE")), _num(row.get(f"{c}_LATITUDE"))) for c in FOOTPRINT_CORNERS]
503 ring_ll = [(lo, la) for lo, la in corners if lo is not None and la is not None]
504 if len(ring_ll) == len(corners):
505 return [[ring_ll]]
507 s, n = _num(row.get("FP_MIN_LAT")), _num(row.get("FP_MAX_LAT"))
508 west, span = _num(row.get("FP_LON_WEST")), _num(row.get("FP_LON_SPAN"))
509 if s is None or n is None or west is None or span is None:
510 s, n = _num(row.get("min_lat")), _num(row.get("max_lat"))
511 west, east = _num(row.get("lon_west")), _num(row.get("lon_east"))
512 if s is None or n is None or west is None or east is None: 512 ↛ 513line 512 didn't jump to line 513 because the condition on line 512 was never true
513 return None
514 span = east - west if east >= west else east - west + 360.0
515 span = min(span, 360.0)
516 return [[_box_ring(s, n, west, span)]]
519def row_footprint(row: dict, wkt_column: str | None = None, proj: Projection | None = None):
520 """Shapely footprint ([`footprint_geometry`][lunar_toolkit.geometry.footprint_geometry]) of a result row,
521 or ``None``."""
522 rings = footprint_rings(row, wkt_column)
523 return None if rings is None else footprint_geometry(rings, proj)
526def aoi_mask(df: pd.DataFrame, aoi: AoiLike, wkt_column: str | None = None) -> list[bool]:
527 """Per-row flag: does the row's footprint intersect the AOI exactly?
529 Footprint source, per row: the WKT column (``FP_WKT`` for CumIndex,
530 ``footprint_wkt`` for observations, or ``wkt_column``), else the four
531 corner lat/lon columns, else the ``FP_*`` / ``min_lat``…``lon_east``
532 bounding box. Rows with no geometry at all are kept.
534 *aoi* is anything [`load_aoi`][lunar_toolkit.geometry.load_aoi] accepts. Regions entirely poleward
535 of ±60° (``POLAR_LAT_THRESHOLD``) are intersected in polar
536 stereographic ([`ps_xy`][lunar_toolkit.geometry.ps_xy], footprint edges straight in the
537 projection); others in lon/lat with the AOI replicated at ±360° so
538 footprints across the 0/360 seam match.
539 """
540 aoi = load_aoi(aoi)
541 if df is None or len(df) == 0: 541 ↛ 542line 541 didn't jump to line 542 because the condition on line 541 was never true
542 return []
543 hemi = _polar_hemi(aoi_region(aoi))
544 proj = ps_projection(hemi) if hemi else None
545 target = prep(aoi_geometry(aoi, proj))
546 if wkt_column is None:
547 wkt_column = next((c for c in _WKT_COLUMNS if c in df.columns), None)
548 mask = []
549 for row in df.to_dict("records"):
550 geom = row_footprint(row, wkt_column, proj)
551 mask.append(True if geom is None or geom.is_empty else bool(target.intersects(geom)))
552 return mask
555def filter_aoi(df: pd.DataFrame, aoi: AoiLike, wkt_column: str | None = None) -> pd.DataFrame:
556 """Keep only rows whose footprint intersects the AOI (see
557 [`aoi_mask`][lunar_toolkit.geometry.aoi_mask])."""
558 if df is None or len(df) == 0: 558 ↛ 559line 558 didn't jump to line 559 because the condition on line 558 was never true
559 return df
560 return df[aoi_mask(df, aoi, wkt_column=wkt_column)].reset_index(drop=True)