Coverage for lunar_toolkit/reference.py: 90%
538 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 02:58 +0000
1"""Reference data for the Artemis III candidate regions: PGDA LOLA 5 m site DEMs and Bertone et al. SfS DEMs."""
3from __future__ import annotations
5import fnmatch
6import io
7import json
8import logging
9import math
10import re
11import shutil
12import time
13import urllib.error
14import urllib.parse
15import urllib.request
16import zipfile
17from collections.abc import Iterable, Sequence
18from dataclasses import dataclass
19from pathlib import Path
20from typing import IO, Any
22import numpy as np
23import pandas as pd
25from .constants import APP_NAME, APP_VERSION, DEFAULT_CACHE_DIR, MOON_RADIUS_KM
26from .csm import CsmCamera, read_camera, transform_camera
27from .downloader import download_file, url_filename
28from .geometry import Aoi, load_aoi, ps_xy
29from .tools import reuse_existing
31log = logging.getLogger(__name__)
33PS_SRS = "+proj=stere +lat_0=-90 +lat_ts=-90 +lon_0=0 +x_0=0 +y_0=0 +R=1737400 +units=m +no_defs"
34"""South polar stereographic on the 1737.4 km sphere (PGDA and Bertone et al. rasters)."""
36ASP_CSV_FORMAT = "1:easting 2:northing 3:height_above_datum"
37"""``--csv-format`` for [`write_asp_csv`][lunar_toolkit.reference.write_asp_csv] output (use with
38``--csv-srs`` [`PS_SRS`][lunar_toolkit.reference.PS_SRS])."""
40ZENODO_RECORD = 17954508
41ZENODO_API_URL = f"https://zenodo.org/api/records/{ZENODO_RECORD}"
42ZENODO_FILE_URL = f"https://zenodo.org/records/{ZENODO_RECORD}/files/{{name}}?download=1"
44PGDA_PRODUCT_PAGE = "https://pgda.gsfc.nasa.gov/products/78"
45PGDA_ROOT = "https://pgda.gsfc.nasa.gov/data/LOLA_5mpp/"
46PGDA_KINDS = ("surf", "ldec", "toterr", "slp", "slperr", "xyzi", "clones", "rois")
48REFERENCE_CACHE_DIR = DEFAULT_CACHE_DIR / "reference"
50_USER_AGENT = f"{APP_NAME}/{APP_VERSION} (python-urllib)"
51_RETRY_CODES = (403, 429, 500, 502, 503, 504) # PGDA answers 403 to bursts of requests
54@dataclass(frozen=True)
55class PsBox:
56 """Axis-aligned box in south polar stereographic metres ([`PS_SRS`][lunar_toolkit.reference.PS_SRS]).
58 Parameters
59 ----------
60 x0, y0, x1, y1 : float
61 West, south, east and north edges (``x0 < x1``, ``y0 < y1``).
62 """
64 x0: float
65 y0: float
66 x1: float
67 y1: float
69 @property
70 def bounds(self) -> tuple[float, float, float, float]:
71 """``(xmin, ymin, xmax, ymax)``."""
72 return self.x0, self.y0, self.x1, self.y1
74 @property
75 def projwin(self) -> tuple[float, float, float, float]:
76 """GDAL / ASP ``projwin`` order: ``(ulx, uly, lrx, lry)``."""
77 return self.x0, self.y1, self.x1, self.y0
79 @property
80 def width(self) -> float:
81 return self.x1 - self.x0
83 @property
84 def height(self) -> float:
85 return self.y1 - self.y0
87 @property
88 def center_xy(self) -> tuple[float, float]:
89 return (self.x0 + self.x1) / 2.0, (self.y0 + self.y1) / 2.0
91 @property
92 def center_latlon(self) -> tuple[float, float]:
93 """Centre as ``(lat, lon)`` degrees, lon in −180..180."""
94 return ps_latlon(*self.center_xy)
96 def contains(self, other: PsBox) -> bool:
97 return self.x0 <= other.x0 and self.y0 <= other.y0 and other.x1 <= self.x1 and other.y1 <= self.y1
99 def intersects(self, other: PsBox) -> bool:
100 return self.x0 < other.x1 and other.x0 < self.x1 and self.y0 < other.y1 and other.y0 < self.y1
102 def ring_lonlat(self, step_m: float = 100.0) -> list[tuple[float, float]]:
103 """Closed ``(lon, lat)`` ring of the box edges, sampled every *step_m* metres."""
104 if self.x0 < 0.0 < self.x1 and self.y0 < 0.0 < self.y1:
105 raise ValueError("Box contains the pole; it has no simple lon/lat outline")
106 corners = [(self.x0, self.y0), (self.x1, self.y0), (self.x1, self.y1), (self.x0, self.y1)]
107 pts: list[tuple[float, float]] = []
108 for (xa, ya), (xb, yb) in zip(corners, corners[1:] + corners[:1]):
109 n = max(1, math.ceil(math.hypot(xb - xa, yb - ya) / step_m))
110 pts += [(xa + (xb - xa) * i / n, ya + (yb - ya) * i / n) for i in range(n)]
111 ring = [ps_latlon(x, y)[::-1] for x, y in pts]
112 return ring + ring[:1]
114 def aoi(self, step_m: float = 100.0) -> Aoi:
115 """The box as a lon/lat [`Aoi`][lunar_toolkit.geometry.Aoi] (edges straight in PS).
117 Parameters
118 ----------
119 step_m : float
120 Vertex spacing along the edges; the AOI matches the box to well
121 under a metre once [`geometry`][lunar_toolkit.geometry] projects it back.
122 """
123 geom = {"type": "Polygon", "coordinates": [self.ring_lonlat(step_m)]}
124 aoi = load_aoi(geom, densify_deg=None)
125 aoi.source = f"PS box {self.bounds}"
126 return aoi
128 def wkt(self, step_m: float = 100.0) -> str:
129 """``POLYGON`` WKT of [`ring_lonlat`][lunar_toolkit.reference.PsBox.ring_lonlat] (lon/lat degrees)."""
130 return "POLYGON((" + ", ".join(f"{lo:.7f} {la:.7f}" for lo, la in self.ring_lonlat(step_m)) + "))"
132 def wkt_ps(self) -> str:
133 """``POLYGON`` WKT of the corners in PS metres."""
134 x0, y0, x1, y1 = self.bounds
135 return f"POLYGON(({x0:g} {y0:g}, {x1:g} {y0:g}, {x1:g} {y1:g}, {x0:g} {y1:g}, {x0:g} {y0:g}))"
138def ps_latlon(x: float, y: float) -> tuple[float, float]:
139 """Inverse of south-hemisphere [`ps_xy`][lunar_toolkit.geometry.ps_xy]: ``(lat, lon)`` degrees."""
140 rho = math.hypot(x, y)
141 lat = math.degrees(2.0 * math.atan(rho / (2.0 * MOON_RADIUS_KM * 1000.0))) - 90.0
142 return lat, math.degrees(math.atan2(x, y))
145@dataclass(frozen=True)
146class PgdaSite:
147 """A PGDA LOLA 5 m/px site DEM (Barker et al. 2021): name, description and PS extent."""
149 name: str
150 description: str
151 extent: PsBox
153 @property
154 def url(self) -> str:
155 return f"{PGDA_ROOT}{self.name}/"
158# Extents read from each site's *_surf.tif (all 5 m/px, PS_SRS, float32, NaN nodata).
159PGDA_SITES: dict[str, PgdaSite] = {
160 s[0]: PgdaSite(s[0], s[1], PsBox(*s[2]))
161 for s in [
162 ("Site01", "Connecting ridge", (-19000, -20000, -3000, -4000)),
163 ("Site04", "Shackleton rim", (-9000, -15000, 7000, 1000)),
164 ("Site06", "Nobile rim 1", (74000, 100000, 94000, 120000)),
165 ("Site07", "Peak near Shackleton", (22000, -28000, 38000, -12000)),
166 ("Site11", "de Gerlache rim", (-45000, 7000, -29000, 23000)),
167 ("Site20", "Leibnitz beta plateau", (65000, 110000, 81000, 126000)),
168 ("Site20v2", "Leibnitz beta plateau, extended boundaries", (60000, 105000, 86000, 131000)),
169 ("Site23", "Malapert massif", (-11000, 111000, 10000, 132000)),
170 ("Site42", "de Gerlache-Kocher massif", (-123400, -66100, -103400, -46100)),
171 ("Haworth", "Haworth", (-52900, 75600, -23100, 105400)),
172 ("Shoemaker", "Shoemaker", (66000, 29000, 86000, 49000)),
173 ("DM1", "Amundsen rim", (154000, 51500, 174000, 71500)),
174 ("DM2", "Nobile rim 2", (147000, 85000, 167000, 105000)),
175 ("SL2", "de Gerlache rim 2", (-58500, 13000, -38500, 33000)),
176 ("SL3", "Connecting ridge extension", (-38700, -16100, -18700, 3900)),
177 ("NPA", "Cabeus exterior wall 1", (-98500, 19500, -78500, 39500)),
178 ("NPB", "Amundsen 1", (175000, -8500, 195000, 11500)),
179 ("NPC", "Idel'son L crater 1", (122500, -109500, 142500, -89500)),
180 ("NPD", "Malapert crater 1", (13000, 146000, 33000, 166000)),
181 ("LM1", "Shackleton Rim B", (2000, -2000, 22000, 18000)),
182 ("LM2", "Shoemaker Rim A", (2000, 47000, 22000, 67000)),
183 ("LM3", "Shoemaker Rim B", (4000, 37000, 24000, 57000)),
184 ("LM4", "Shoemaker Rim C", (9000, 15000, 29000, 35000)),
185 ("LM5", "Shoemaker Rim D", (30000, 1000, 50000, 21000)),
186 ("LM6", "Shoemaker Rim E", (49000, 10000, 69000, 30000)),
187 ("LM7", "Faustini Rim A", (54000, -10000, 74000, 10000)),
188 ("LM8", "Shoemaker Rim F", (60000, 30000, 80000, 50000)),
189 ]
190}
193@dataclass(frozen=True)
194class Region:
195 """A 2022 Artemis III candidate landing region as processed by Bertone et al. (2026).
197 ``bbox`` (``min_lon, min_lat, max_lon, max_lat``) and ``center``
198 (``lon, lat``, the bbox midpoint) come from the Zenodo README; ``roi`` is
199 the polar stereographic square the SfS DEM covers, which is the extent of
200 the matching PGDA LOLA 5 m site DEM (``pgda_site``).
201 """
203 id: int
204 name: str
205 bbox: tuple[float, float, float, float]
206 center: tuple[float, float]
207 pgda_site: str | None
208 zenodo_file: str
209 zenodo_size: int
210 zenodo_md5: str
212 @property
213 def slug(self) -> str:
214 return re.sub(r"[^a-z0-9]+", "-", self.name.lower()).strip("-")
216 @property
217 def roi(self) -> PsBox:
218 if self.pgda_site is None: 218 ↛ 219line 218 didn't jump to line 219 because the condition on line 218 was never true
219 return bbox_ps_box(self.bbox)
220 return PGDA_SITES[self.pgda_site].extent
222 @property
223 def projwin(self) -> tuple[float, float, float, float]:
224 """``(ulx, uly, lrx, lry)`` of `roi` in PS metres."""
225 return self.roi.projwin
227 @property
228 def bbox_projwin(self) -> tuple[float, float, float, float]:
229 """``projwin`` of the PS envelope of the lon/lat `bbox` (a wedge; larger than `roi`)."""
230 return bbox_ps_box(self.bbox).projwin
232 @property
233 def zenodo_url(self) -> str:
234 return ZENODO_FILE_URL.format(name=self.zenodo_file)
236 def aoi(self, step_m: float = 100.0) -> Aoi:
237 """`roi` as an [`Aoi`][lunar_toolkit.geometry.Aoi]."""
238 aoi = self.roi.aoi(step_m)
239 aoi.source = f"A3CLR22 #{self.id} {self.name}"
240 return aoi
242 def bbox_aoi(self, densify_deg: float = 0.05) -> Aoi:
243 """The lon/lat `bbox` as an AOI, densified so it stays a wedge in PS."""
244 w, s, e, n = self.bbox
245 aoi = load_aoi(f"POLYGON(({w} {s}, {e} {s}, {e} {n}, {w} {n}, {w} {s}))", densify_deg=densify_deg)
246 aoi.source = f"A3CLR22 #{self.id} {self.name} bbox"
247 return aoi
250# Zenodo README table + file sizes / MD5s from the record API (id, label, center lon/lat,
251# bbox min lon/lat max lon/lat, PGDA site, zip, size, md5).
252_REGION_ROWS = [
253 (0, "Amundsen Rim", (69.303844, -84.225586), (65.095231, -84.648844, 73.512457, -83.802328), "DM1",
254 "A3CLR22_0_Amundsen_Rim.zip", 1775016989, "d1e12877277cf9715381b8795c0873cb"),
255 (1, "Nobile Rim 2", (58.743523, -83.953082), (54.462322, -84.404616, 63.024724, -83.501549), "DM2",
256 "A3CLR22_1_Nobile_Rim_2.zip", 1839070131, "964310cb7f2a5033e2e12cd34e7aa7b1"),
257 (2, "Haworth", (-23.671814, -86.752966), (-34.981846, -87.393534, -12.361782, -86.112398), "Haworth",
258 "A3CLR22_2_Haworth.zip", 4353812518, "e83a923c41c00be51396d06ab7cf6e9e"),
259 (3, "Faustini Rim A", (90.0, -87.863451), (79.508523, -88.189067, 100.491477, -87.537834), "LM7",
260 "A3CLR22_3_Faustini_Rim_A.zip", 1936729230, "161db71236036f7c9b98b501d353482a"),
261 (4, "de Gerlache Rim 2", (-63.434949, -88.222637), (-77.471192, -88.659986, -49.398705, -87.785288), "SL2",
262 "A3CLR22_4_deGerlache_Rim_2.zip", 2187384204, "52921a134793a6b35b0eea6f98ad79d8"),
263 (5, "Connecting Ridge Extension", (-104.473348, -88.993919), (-130.727219, -89.370051, -78.219476, -88.617786),
264 "SL3", "A3CLR22_5_Connecting_Ridge_Extension.zip", 2225897760, "9141ec5d298a3aebab4d3671f7f38a79"),
265 (6, "Connecting Ridge", (-136.678946, -89.462697), (-171.469234, -89.835111, -101.888658, -89.090284),
266 "Site01", "A3CLR22_6_Connecting_Ridge.zip", 2206485514, "0ad982cf42319a300f2a2cff69d39322"),
267 (7, "Nobile Rim 1", (37.444641, -85.43775), (31.660751, -85.899218, 43.22853, -84.976283), "Site06",
268 "A3CLR22_7_Nobile_Rim_1.zip", 1962014939, "a8bcb60f432253f5fa97ff07899bd02a"),
269 (8, "Peak Near Shackleton", (124.684171, -88.808537), (107.525568, -89.173591, 141.842773, -88.443484),
270 "Site07", "A3CLR22_8_Peak_Near_Shackleton.zip", 1360387336, "9ac25cc18b1290e679dfad9b1eaabad7"),
271 (9, "de Gerlache Rim", (-66.370065, -88.674855), (-81.158185, -89.016199, -51.581945, -88.33351), "Site11",
272 "A3CLR22_9_deGerlache_Rim.zip", 1340667387, "a34e603e9d2c2c9a2796fc237d24b1b4"),
273 (10, "Leibnitz Beta Plateau", (31.827284, -85.425819), (27.287988, -85.788335, 36.36658, -85.063304),
274 "Site20", "A3CLR22_10_Leibnitz_Beta_Plateau.zip", 1050441742, "b0a606e6a4ce73faf2ea68c91ce0b93e"),
275 (11, "Malapert Massif", (-0.255799, -85.979917), (-5.659482, -86.325889, 5.147885, -85.633944), "Site23",
276 "A3CLR22_11_Malapert_Massif.zip", 2213137537, "772a578b8f2ebdb078c5a04a3c07182a"),
277 (12, "de Gerlache-Kocher Massif", (-116.537066, -85.826913), (-122.589357, -86.267855, -110.484776, -85.385972),
278 "Site42", "A3CLR22_12_deGerlache_Kocher_Massif.zip", 2108633565, "78b8a042d0c6be2923abdd288e436c23"),
279] # fmt: skip
281REGIONS: tuple[Region, ...] = tuple(
282 Region(id=i, name=name, center=c, bbox=b, pgda_site=site, zenodo_file=f, zenodo_size=size, zenodo_md5=md5)
283 for i, name, c, b, site, f, size, md5 in _REGION_ROWS
284)
287def _norm(s: str) -> str:
288 return re.sub(r"[^a-z0-9]", "", s.lower())
291def get_region(key: Region | int | str) -> Region:
292 """Look up a region by id, name, slug or PGDA site (case-insensitive).
294 ``11``, ``"11"``, ``"Malapert Massif"``, ``"malapert-massif"``,
295 ``"malapert"`` and ``"Site23"`` all give region 11. A partial name must
296 match exactly one region.
297 """
298 if isinstance(key, Region):
299 return key
300 if isinstance(key, int) or str(key).strip().isdigit():
301 i = int(key)
302 for r in REGIONS:
303 if r.id == i:
304 return r
305 raise KeyError(f"No region with id {i} (0–{len(REGIONS) - 1})")
306 k = _norm(str(key))
307 if not k: 307 ↛ 308line 307 didn't jump to line 308 because the condition on line 307 was never true
308 raise KeyError("Empty region key")
309 for r in REGIONS:
310 if k in (
311 _norm(r.name),
312 _norm(r.slug),
313 _norm(r.pgda_site or ""),
314 _norm(r.zenodo_file.removesuffix(".zip")),
315 ):
316 return r
317 hits = [r for r in REGIONS if k and k in _norm(r.name)]
318 if len(hits) == 1:
319 return hits[0]
320 if hits:
321 raise KeyError(f"{key!r} is ambiguous: " + ", ".join(f"{r.id} {r.name}" for r in hits))
322 raise KeyError(f"No region matches {key!r}; known: " + ", ".join(r.name for r in REGIONS))
325def regions_table() -> pd.DataFrame:
326 """All regions as a table (one row per region)."""
327 rows = []
328 for r in REGIONS:
329 rows.append(
330 {
331 "id": r.id,
332 "name": r.name,
333 "slug": r.slug,
334 "center_lon": r.center[0],
335 "center_lat": r.center[1],
336 "min_lon": r.bbox[0],
337 "min_lat": r.bbox[1],
338 "max_lon": r.bbox[2],
339 "max_lat": r.bbox[3],
340 "pgda_site": r.pgda_site,
341 "roi_projwin": r.projwin,
342 "zenodo_file": r.zenodo_file,
343 "zenodo_size": r.zenodo_size,
344 "zenodo_md5": r.zenodo_md5,
345 }
346 )
347 return pd.DataFrame(rows)
350def bbox_ps_box(bbox: Sequence[float], samples: int = 200) -> PsBox:
351 """PS envelope of a lon/lat ``(min_lon, min_lat, max_lon, max_lat)`` box (south)."""
352 w, s, e, n = bbox
353 pts = []
354 for i in range(samples + 1):
355 f = i / samples
356 lo, la = w + (e - w) * f, s + (n - s) * f
357 pts += [ps_xy(s, lo, "S"), ps_xy(n, lo, "S"), ps_xy(la, w, "S"), ps_xy(la, e, "S")]
358 xs, ys = [p[0] for p in pts], [p[1] for p in pts]
359 return PsBox(min(xs), min(ys), max(xs), max(ys))
362def crop_window(
363 region: Region | int | str,
364 size_km: float = 2.5,
365 center: tuple[float, float] | None = None,
366 snap_m: float | None = 5.0,
367) -> PsBox:
368 """Square PS window of *size_km* centred on a region.
370 Parameters
371 ----------
372 region : Region, int or str
373 Anything [`get_region`][lunar_toolkit.reference.get_region] accepts.
374 size_km : float
375 Side length (2.5 km is the Bertone et al. uncertainty-ensemble crop).
376 center : (lat, lon), optional
377 Centre in degrees; default the centre of the region's `Region.roi`.
378 snap_m : float, optional
379 Snap the edges to multiples of this many metres so the window lines
380 up with the 5 m PGDA / SfS pixel grid; ``None`` keeps them exact.
382 Returns
383 -------
384 PsBox
385 Use ``.projwin`` for ASP / GDAL and ``.aoi()`` for searches.
386 """
387 reg = get_region(region)
388 cx, cy = ps_xy(center[0], center[1], "S") if center is not None else reg.roi.center_xy
389 side = size_km * 1000.0
390 x0, y0 = cx - side / 2.0, cy - side / 2.0
391 if snap_m:
392 x0, y0 = round(x0 / snap_m) * snap_m, round(y0 / snap_m) * snap_m
393 box = PsBox(x0, y0, x0 + side, y0 + side)
394 if not reg.roi.contains(box):
395 log.warning("Crop window %s extends outside the %s ROI %s", box.bounds, reg.name, reg.roi.bounds)
396 return box
399def _site_name(site_or_region: Region | int | str) -> str:
400 if isinstance(site_or_region, str) and site_or_region in PGDA_SITES:
401 return site_or_region
402 reg = get_region(site_or_region)
403 if reg.pgda_site is None: 403 ↛ 404line 403 didn't jump to line 404 because the condition on line 403 was never true
404 raise KeyError(f"Region {reg.id} {reg.name} has no PGDA LOLA 5 m site")
405 return reg.pgda_site
408def _classify_pgda(path: str) -> str | None:
409 name = path.rsplit("/", 1)[-1]
410 if "/Clones/" in path:
411 return "clones"
412 if name.endswith("_ROIs.tgz"):
413 return "rois"
414 if name.endswith(".xyzi"):
415 return "xyzi"
416 m = re.search(r"_(surf|ldec|toterr|slp|slperr)\.tif$", name)
417 return m.group(1) if m else None
420def parse_pgda_index(html: str) -> dict[str, dict[str, list[str]]]:
421 """``site → kind → [url, ...]`` from the PGDA product page HTML.
423 Links may be absolute, root-relative or relative to `PGDA_ROOT`
424 (``DM2/Clones/DM2_final_adj_5mpp_0001_err.tif``); clones are sorted by
425 number.
426 """
427 index: dict[str, dict[str, list[str]]] = {}
428 host = "{0.scheme}://{0.netloc}".format(urllib.parse.urlparse(PGDA_ROOT))
429 for href in re.findall(r"""href\s*=\s*["']?([^"'\s>]+)""", html, flags=re.I):
430 if href.startswith(("http://", "https://")):
431 url = href.replace("http://", "https://", 1)
432 elif href.startswith("/"):
433 url = host + href
434 else:
435 url = PGDA_ROOT + href
436 if not url.startswith(PGDA_ROOT):
437 continue
438 rel = url[len(PGDA_ROOT) :]
439 kind = _classify_pgda(rel)
440 if kind is None or "/" not in rel:
441 continue
442 urls = index.setdefault(rel.split("/", 1)[0], {}).setdefault(kind, [])
443 if url not in urls:
444 urls.append(url)
445 for kinds in index.values():
446 kinds.get("clones", []).sort()
447 return index
450def _get(url: str, headers: dict[str, str] | None = None, timeout: float = 60.0, retries: int = 4) -> bytes:
451 """GET *url* with retries on throttling / server errors."""
452 req = urllib.request.Request(url, headers={"User-Agent": _USER_AGENT, **(headers or {})})
453 for attempt in range(retries + 1): 453 ↛ 464line 453 didn't jump to line 464 because the loop on line 453 didn't complete
454 try:
455 with urllib.request.urlopen(req, timeout=timeout) as resp:
456 return resp.read()
457 except urllib.error.HTTPError as exc:
458 if exc.code not in _RETRY_CODES or attempt == retries:
459 raise
460 except urllib.error.URLError:
461 if attempt == retries:
462 raise
463 time.sleep(2.0 * 2**attempt)
464 raise AssertionError("unreachable")
467def pgda_index(
468 refresh: bool = False, cache_dir: Path = REFERENCE_CACHE_DIR, html: str | None = None
469) -> dict[str, dict[str, list[str]]]:
470 """PGDA LOLA 5 m file index, parsed from the product page and cached.
472 Parameters
473 ----------
474 refresh : bool
475 Re-fetch `PGDA_PRODUCT_PAGE` even when a cached index exists.
476 cache_dir : Path
477 Where ``pgda_lola_5mpp_index.json`` is kept.
478 html : str, optional
479 Parse this page text instead of fetching it (nothing is cached).
480 """
481 if html is not None:
482 return parse_pgda_index(html)
483 cache = Path(cache_dir) / "pgda_lola_5mpp_index.json"
484 if cache.exists() and not refresh: 484 ↛ 486line 484 didn't jump to line 486 because the condition on line 484 was always true
485 return json.loads(cache.read_text())
486 index = parse_pgda_index(_get(PGDA_PRODUCT_PAGE).decode("utf-8", "replace"))
487 if not index:
488 raise RuntimeError(f"No LOLA 5 m files found on {PGDA_PRODUCT_PAGE}")
489 cache.parent.mkdir(parents=True, exist_ok=True)
490 cache.write_text(json.dumps(index, indent=1))
491 return index
494def _clone_number(url: str) -> int:
495 m = re.search(r"_(\d{4})_err\.tif$", url)
496 return int(m.group(1)) if m else 0
499def pgda_urls(
500 site_or_region: Region | int | str,
501 kinds: Iterable[str] = ("surf",),
502 clones: int | Iterable[int] | None = None,
503 index: dict[str, dict[str, list[str]]] | None = None,
504) -> dict[str, list[str]]:
505 """``kind → [url, ...]`` for one PGDA site.
507 Parameters
508 ----------
509 site_or_region : Region, int or str
510 A PGDA site name (``"Site23"``) or anything [`get_region`][lunar_toolkit.reference.get_region]
511 accepts.
512 kinds : iterable of str
513 Members of `PGDA_KINDS`.
514 clones : int or iterable of int, optional
515 With ``"clones"``: the first *n* clones, or these clone numbers
516 (1-based); default all 100.
517 index : dict, optional
518 A [`pgda_index`][lunar_toolkit.reference.pgda_index] result (default: the cached index).
519 """
520 site = _site_name(site_or_region)
521 kinds = list(kinds)
522 bad = [k for k in kinds if k not in PGDA_KINDS]
523 if bad:
524 raise ValueError(f"Unknown PGDA kind(s) {bad}; choose from {', '.join(PGDA_KINDS)}")
525 files = (index if index is not None else pgda_index()).get(site)
526 if not files: 526 ↛ 527line 526 didn't jump to line 527 because the condition on line 526 was never true
527 raise KeyError(f"PGDA site {site} not in the file index")
528 out: dict[str, list[str]] = {}
529 for kind in kinds:
530 urls = list(files.get(kind, []))
531 if kind == "clones" and clones is not None:
532 if isinstance(clones, int):
533 urls = urls[:clones]
534 else:
535 wanted = set(clones)
536 urls = [u for u in urls if _clone_number(u) in wanted]
537 if not urls: 537 ↛ 538line 537 didn't jump to line 538 because the condition on line 537 was never true
538 log.warning("PGDA site %s has no %s files", site, kind)
539 out[kind] = urls
540 return out
543def fetch(
544 url: str,
545 dest: Path,
546 expected_md5: str | None = None,
547 force: bool = False,
548 show_progress: bool = True,
549 retries: int = 4,
550) -> Path:
551 """Resumable download of *url* to *dest*, skipped when *dest* exists.
553 Throttling replies (403 / 429) and server errors are retried with
554 back-off; the ``.part`` file is resumed each time.
555 """
556 dest = Path(dest)
557 if reuse_existing(dest, force):
558 return dest
559 part = dest.with_suffix(dest.suffix + ".part")
560 if force and part.exists(): 560 ↛ 561line 560 didn't jump to line 561 because the condition on line 560 was never true
561 part.unlink()
562 for attempt in range(retries + 1): 562 ↛ 573line 562 didn't jump to line 573 because the loop on line 562 didn't complete
563 try:
564 return download_file(url, dest, show_progress=show_progress, expected_md5=expected_md5)
565 except urllib.error.HTTPError as exc:
566 if exc.code not in _RETRY_CODES or attempt == retries: 566 ↛ 567line 566 didn't jump to line 567 because the condition on line 566 was never true
567 raise
568 except (urllib.error.URLError, TimeoutError, ConnectionError):
569 if attempt == retries:
570 raise
571 log.info("Retrying %s", url)
572 time.sleep(2.0 * 2**attempt)
573 raise AssertionError("unreachable")
576def download_pgda(
577 site_or_region: Region | int | str,
578 kinds: Iterable[str] = ("surf",),
579 out_dir: Path = Path("."),
580 clones: int | Iterable[int] | None = None,
581 force: bool = False,
582 show_progress: bool = True,
583 index: dict[str, dict[str, list[str]]] | None = None,
584) -> dict[str, list[Path]]:
585 """Download PGDA LOLA 5 m files of one site by kind (sequentially; resumable).
587 Files land in ``out_dir/<site>/`` (clones in ``out_dir/<site>/Clones/``);
588 existing ones are kept unless *force*. Arguments as in [`pgda_urls`][lunar_toolkit.reference.pgda_urls].
590 Returns
591 -------
592 dict
593 ``kind → [local path, ...]``.
594 """
595 site = _site_name(site_or_region)
596 out: dict[str, list[Path]] = {}
597 for kind, urls in pgda_urls(site, kinds, clones, index).items():
598 sub = Path(out_dir) / site / ("Clones" if kind == "clones" else "")
599 out[kind] = [fetch(u, sub / url_filename(u), force=force, show_progress=show_progress) for u in urls]
600 return out
603def read_xyzi(path: Path, box: PsBox | None = None) -> pd.DataFrame:
604 """Read a PGDA ``.xyzi`` LOLA point cloud.
606 The file is little-endian float64 rows of ``X, Y, Z, RDRid`` with X/Y the
607 polar stereographic coordinates and Z the height above the 1737.4 km
608 sphere, all in km (the same Z as the ``*_surf.tif`` metres), and RDRid the
609 LOLA RDR ``YYDOYHHMM``.
611 Parameters
612 ----------
613 path : Path
614 The ``.xyzi`` file.
615 box : PsBox, optional
616 Keep only points inside this window.
618 Returns
619 -------
620 pandas.DataFrame
621 Columns ``x``, ``y``, ``height`` (metres) and ``rdr_id`` (int64).
622 """
623 raw = np.fromfile(path, dtype="<f8")
624 if raw.size % 4:
625 raise ValueError(f"{path}: size is not a multiple of 4 float64 values")
626 pts = raw.reshape(-1, 4)
627 if box is not None:
628 x, y = pts[:, 0] * 1000.0, pts[:, 1] * 1000.0
629 pts = pts[(x >= box.x0) & (x <= box.x1) & (y >= box.y0) & (y <= box.y1)]
630 return pd.DataFrame(
631 {
632 "x": pts[:, 0] * 1000.0,
633 "y": pts[:, 1] * 1000.0,
634 "height": pts[:, 2] * 1000.0,
635 "rdr_id": pts[:, 3].astype(np.int64),
636 }
637 )
640def write_asp_csv(points: pd.DataFrame | Path, out: Path, box: PsBox | None = None) -> Path:
641 """Write LOLA points as an ASP CSV of ``easting, northing, height_above_datum``.
643 Use with ``--csv-format`` [`ASP_CSV_FORMAT`][lunar_toolkit.reference.ASP_CSV_FORMAT] and ``--csv-srs``
644 [`PS_SRS`][lunar_toolkit.reference.PS_SRS] (the 1737.4 km sphere is then the datum), e.g.
645 ``pc_align --csv-format '1:easting 2:northing 3:height_above_datum'
646 --csv-srs '+proj=stere +lat_0=-90 +lon_0=0 +R=1737400 +units=m' dem.tif lola.csv``.
648 Parameters
649 ----------
650 points : DataFrame or Path
651 A [`read_xyzi`][lunar_toolkit.reference.read_xyzi] table or an ``.xyzi`` path.
652 out : Path
653 CSV to write; a ``#`` header line (ignored by ASP) names the columns.
654 box : PsBox, optional
655 Keep only points inside this window.
656 """
657 df = read_xyzi(Path(points), box) if not isinstance(points, pd.DataFrame) else points
658 if box is not None and isinstance(points, pd.DataFrame): 658 ↛ 659line 658 didn't jump to line 659 because the condition on line 658 was never true
659 df = df[(df.x >= box.x0) & (df.x <= box.x1) & (df.y >= box.y0) & (df.y <= box.y1)]
660 out = Path(out)
661 out.parent.mkdir(parents=True, exist_ok=True)
662 with open(out, "w") as f:
663 f.write(f"# easting_m,northing_m,height_above_datum_m,rdr_id csv-srs: {PS_SRS}\n")
664 np.savetxt(
665 f,
666 df[["x", "y", "height", "rdr_id"]].to_numpy(),
667 fmt=["%.3f", "%.3f", "%.3f", "%d"],
668 delimiter=",",
669 )
670 return out
673def parse_zenodo_files(record: dict[str, Any]) -> list[dict[str, Any]]:
674 """``[{"key", "size", "md5", "url"}, ...]`` from a Zenodo record API response."""
675 files = []
676 for f in record.get("files", []):
677 algo, _, digest = str(f.get("checksum", "")).partition(":")
678 files.append(
679 {
680 "key": f["key"],
681 "size": int(f.get("size", 0)),
682 "md5": digest if algo == "md5" else None,
683 "url": ZENODO_FILE_URL.format(name=f["key"]),
684 }
685 )
686 return sorted(files, key=lambda f: f["key"])
689def zenodo_files(refresh: bool = False, cache_dir: Path = REFERENCE_CACHE_DIR) -> list[dict[str, Any]]:
690 """Files of the Bertone et al. Zenodo record (fetched once, then cached)."""
691 cache = Path(cache_dir) / f"zenodo_{ZENODO_RECORD}.json"
692 if cache.exists() and not refresh: 692 ↛ 695line 692 didn't jump to line 695 because the condition on line 692 was always true
693 record = json.loads(cache.read_text())
694 else:
695 record = json.loads(_get(ZENODO_API_URL))
696 cache.parent.mkdir(parents=True, exist_ok=True)
697 cache.write_text(json.dumps({"id": record.get("id"), "files": record.get("files", [])}, indent=1))
698 return parse_zenodo_files(record)
701def download_sdem(
702 region: Region | int | str,
703 out_dir: Path = Path("."),
704 force: bool = False,
705 show_progress: bool = True,
706) -> Path:
707 """Download a region's Bertone et al. zip (1–4 GB; resumable, MD5-checked)."""
708 reg = get_region(region)
709 return fetch(
710 reg.zenodo_url,
711 Path(out_dir) / reg.zenodo_file,
712 expected_md5=reg.zenodo_md5,
713 force=force,
714 show_progress=show_progress,
715 )
718class HttpRangeFile(io.RawIOBase):
719 """Read-only, seekable view of a remote file through HTTP Range requests.
721 Lets `zipfile.ZipFile` read a remote zip's central directory or a
722 single member without downloading the archive.
723 """
725 def __init__(self, url: str, timeout: float = 60.0) -> None:
726 self.url, self.timeout, self.pos = url, timeout, 0
727 self.requests = 0
728 req = urllib.request.Request(url, headers={"User-Agent": _USER_AGENT, "Range": "bytes=0-0"})
729 with urllib.request.urlopen(req, timeout=timeout) as resp:
730 rng = resp.headers.get("Content-Range", "")
731 if resp.status != 206 or "/" not in rng: 731 ↛ 732line 731 didn't jump to line 732 because the condition on line 731 was never true
732 raise OSError(f"{url} does not support HTTP Range requests")
733 self.size = int(rng.rsplit("/", 1)[1])
735 def readable(self) -> bool:
736 return True
738 def seekable(self) -> bool:
739 return True
741 def tell(self) -> int:
742 return self.pos
744 def seek(self, offset: int, whence: int = io.SEEK_SET) -> int:
745 base = {io.SEEK_SET: 0, io.SEEK_CUR: self.pos, io.SEEK_END: self.size}[whence]
746 self.pos = max(0, base + offset)
747 return self.pos
749 def readinto(self, b: Any) -> int:
750 n = min(len(b), self.size - self.pos)
751 if n <= 0: 751 ↛ 752line 751 didn't jump to line 752 because the condition on line 751 was never true
752 return 0
753 data = _get(self.url, {"Range": f"bytes={self.pos}-{self.pos + n - 1}"}, timeout=self.timeout)[:n]
754 self.requests += 1
755 b[: len(data)] = data
756 self.pos += len(data)
757 return len(data)
760def _open_zip(source: Region | int | str | Path, buffer_size: int = 4 << 20) -> zipfile.ZipFile:
761 """Local zip path, or a region's remote Zenodo zip (read via HTTP Range, *buffer_size* per request)."""
762 if isinstance(source, Path) or (isinstance(source, str) and Path(source).is_file()):
763 return zipfile.ZipFile(source)
764 fh: IO[bytes] = io.BufferedReader(HttpRangeFile(get_region(source).zenodo_url), buffer_size=buffer_size)
765 return zipfile.ZipFile(fh)
768def _is_junk(name: str) -> bool:
769 return name.startswith("__MACOSX/") or name.endswith((".DS_Store", "/"))
772def sdem_members(source: Region | int | str | Path) -> pd.DataFrame:
773 """Members of a Bertone et al. zip: a local path, or a region (remote, central directory only).
775 Returns
776 -------
777 pandas.DataFrame
778 ``name`` (path inside the zip), ``size`` and ``compressed_size``
779 (bytes); macOS metadata and directories are left out.
780 """
781 with _open_zip(source) as zf:
782 rows = [
783 {"name": i.filename, "size": i.file_size, "compressed_size": i.compress_size}
784 for i in zf.infolist()
785 if not _is_junk(i.filename)
786 ]
787 return pd.DataFrame(rows, columns=["name", "size", "compressed_size"])
790def select_members(names: Iterable[str], patterns: Iterable[str]) -> list[str]:
791 """Zip member names matching any pattern.
793 A pattern with ``*``, ``?`` or ``[`` is a glob against the full name and
794 the base name (``"SLOPE_MAPS/*_100.tif"``); otherwise it is a
795 case-insensitive substring of the base name (``"GLDELEV"``, ``"README"``).
796 """
797 pats = list(patterns)
798 out = []
799 for name in names:
800 if _is_junk(name):
801 continue
802 base = name.rsplit("/", 1)[-1]
803 for p in pats:
804 if any(c in p for c in "*?["):
805 hit = fnmatch.fnmatch(name, p) or fnmatch.fnmatch(base, p) or fnmatch.fnmatch(name, f"*/{p}")
806 else:
807 hit = p.lower() in base.lower()
808 if hit:
809 out.append(name)
810 break
811 return out
814def extract_sdem(
815 source: Region | int | str | Path,
816 members: Iterable[str] = ("GLDELEV",),
817 out_dir: Path = Path("."),
818 force: bool = False,
819) -> list[Path]:
820 """Extract selected members of a Bertone et al. zip.
822 Parameters
823 ----------
824 source : Path, Region, int or str
825 A downloaded zip, or a region whose Zenodo zip is then read remotely
826 member by member (HTTP Range; only the selected members are fetched).
827 members : iterable of str
828 [`select_members`][lunar_toolkit.reference.select_members] patterns, e.g.
829 ``("GLDELEV", "GLDSIGM", "README")``.
830 out_dir : Path
831 Destination; paths below the zip's top-level folder are kept
832 (``SLOPE_MAPS/A311_GLDTSLP_005_100.tif``).
833 force : bool
834 Re-extract members that already exist.
836 Returns
837 -------
838 list of Path
839 """
840 out_dir = Path(out_dir)
841 paths = []
842 with _open_zip(source) as zf:
843 names = select_members(zf.namelist(), members)
844 if not names:
845 raise KeyError(f"No zip member matches {list(members)}")
846 for name in names:
847 rel = name.split("/", 1)[1] if "/" in name else name
848 dest = out_dir / rel
849 if not reuse_existing(dest, force):
850 dest.parent.mkdir(parents=True, exist_ok=True)
851 tmp = dest.with_suffix(dest.suffix + ".part")
852 with zf.open(name) as src, open(tmp, "wb") as dst:
853 shutil.copyfileobj(src, dst, 1 << 20)
854 tmp.replace(dest)
855 log.info("Extracted %s", dest)
856 paths.append(dest)
857 return paths
860NAC_FOCAL_MM = {"LEFT": 699.62, "RIGHT": 701.57}
861"""LROC NAC focal lengths (mm) of the ISIS / ALE camera models, by ``name_sensor`` suffix."""
863_PID_RE = re.compile(r"_(M\d+[LR]E)\.(?:json|tif)$")
866def _member_pid(name: str) -> str | None:
867 m = _PID_RE.search(name)
868 return m.group(1) if m else None
871def published_members(source: Region | int | str | Path, kind: str = "GLDACAM") -> pd.DataFrame:
872 """Per-image members of a Bertone et al. zip: cameras (``GLDACAM``) or mapprojected images (``GLDIMGP``).
874 Returns
875 -------
876 pandas.DataFrame
877 ``product_id``, ``name`` and ``size``; regions 0, 7, 8 and 10 have
878 no ``LROC-NAC-adjust/`` cameras (the table is then empty).
879 """
880 members = sdem_members(source)
881 hit = members[members["name"].str.contains(f"_{kind}_", regex=False)].copy()
882 hit.insert(0, "product_id", hit["name"].map(_member_pid))
883 return hit.dropna(subset=["product_id"])[["product_id", "name", "size"]].reset_index(drop=True)
886def parse_image_status(text: str) -> pd.DataFrame:
887 """Parse ``extra/image_status_*.csv``: ``product_id`` and ``status``.
889 Statuses are ``used_for_sfs`` (the SfS inputs, one ``GLDACAM`` camera
890 each), ``not_aligned`` (dropped by bundle adjustment) and
891 ``available_not_selected``. The ``# key: value`` header lines go to
892 ``attrs["summary"]`` (counts as ints).
893 """
894 summary: dict[str, int] = {}
895 rows = []
896 for line in text.splitlines():
897 line = line.strip()
898 if not line:
899 continue
900 if line.startswith("#"):
901 key, sep, val = line[1:].partition(":")
902 if sep and val.strip().isdigit():
903 summary[key.strip()] = int(val)
904 continue
905 pid, _, status = line.partition(",")
906 if pid != "img_name":
907 rows.append((pid.strip(), status.strip()))
908 df = pd.DataFrame(rows, columns=["product_id", "status"])
909 df.attrs["summary"] = summary
910 return df
913def image_status(
914 source: Region | int | str | Path, out_dir: Path = REFERENCE_CACHE_DIR / "published", force: bool = False
915) -> pd.DataFrame:
916 """A region's image status table ([`parse_image_status`][lunar_toolkit.reference.parse_image_status]),
917 extracted once to *out_dir*."""
918 (path,) = extract_sdem(source, ["extra/image_status_*.csv"], Path(out_dir), force=force)
919 return parse_image_status(path.read_text())
922def parse_tiff_extent(head: bytes) -> PsBox:
923 """Extent of a north-up GeoTIFF from its first bytes (header, IFD and georeferencing tags).
925 Reads ``ImageWidth`` / ``ImageLength`` with ``ModelPixelScale`` and
926 ``ModelTiepoint`` (or ``ModelTransformation``) from a classic TIFF whose
927 IFD and tag values lie inside *head* (GDAL writes them first, so 16 KB
928 is plenty for the Bertone et al. ``GLDIMGP`` maps).
929 """
930 import struct
932 bo = {b"II": "<", b"MM": ">"}.get(head[:2])
933 if bo is None or struct.unpack(bo + "H", head[2:4])[0] != 42:
934 raise ValueError("Not a classic TIFF")
935 sizes = {1: 1, 2: 1, 3: 2, 4: 4, 5: 8, 6: 1, 7: 1, 8: 2, 9: 4, 10: 8, 11: 4, 12: 8}
936 codes = {3: "H", 4: "I", 8: "h", 9: "i", 11: "f", 12: "d"}
937 (ifd,) = struct.unpack(bo + "I", head[4:8])
938 (n,) = struct.unpack(bo + "H", head[ifd : ifd + 2])
939 tags: dict[int, tuple[Any, ...]] = {}
940 for i in range(n):
941 tag, typ, cnt = struct.unpack(bo + "HHI", head[ifd + 2 + 12 * i : ifd + 10 + 12 * i])
942 if typ not in codes:
943 continue
944 size = sizes[typ] * cnt
945 raw = head[ifd + 10 + 12 * i : ifd + 14 + 12 * i]
946 if size > 4:
947 (off,) = struct.unpack(bo + "I", raw)
948 if off + size > len(head):
949 raise ValueError(f"TIFF tag {tag} lies beyond the first {len(head)} bytes")
950 raw = head[off : off + size]
951 tags[tag] = struct.unpack(bo + codes[typ] * cnt, raw[:size])
952 w, h = tags[256][0], tags[257][0]
953 if 34264 in tags: 953 ↛ 954line 953 didn't jump to line 954 because the condition on line 953 was never true
954 t = tags[34264]
955 sx, sy, x0, y0 = t[0], -t[5], t[3], t[7]
956 elif 33550 in tags and 33922 in tags: 956 ↛ 961line 956 didn't jump to line 961 because the condition on line 956 was always true
957 sx, sy = tags[33550][:2]
958 i, j, _, x, y, _ = tags[33922][:6]
959 x0, y0 = x - i * sx, y + j * sy
960 else:
961 raise ValueError("TIFF has no GeoTIFF georeferencing tags")
962 return PsBox(x0, y0 - h * sy, x0 + w * sx, y0)
965def _member_head(zf: zipfile.ZipFile, info: zipfile.ZipInfo, n: int) -> bytes:
966 """First *n* uncompressed bytes of a stored or deflated member, in one read of the archive."""
967 import struct
968 import zlib
970 fp = zf.fp
971 assert fp is not None
972 fp.seek(info.header_offset)
973 raw = fp.read(30 + len(info.filename.encode()) + 1024 + 2 * n)
974 if raw[:4] != b"PK\x03\x04": 974 ↛ 975line 974 didn't jump to line 975 because the condition on line 974 was never true
975 raise zipfile.BadZipFile(f"Bad local header for {info.filename}")
976 name_len, extra_len = struct.unpack("<HH", raw[26:30])
977 data = raw[30 + name_len + extra_len :][: info.compress_size]
978 if info.compress_type == zipfile.ZIP_STORED: 978 ↛ 979line 978 didn't jump to line 979 because the condition on line 978 was never true
979 return data[:n]
980 if info.compress_type != zipfile.ZIP_DEFLATED: 980 ↛ 981line 980 didn't jump to line 981 because the condition on line 980 was never true
981 raise NotImplementedError(f"Compression {info.compress_type} of {info.filename}")
982 return zlib.decompressobj(-15).decompress(data, n)
985def published_extents(
986 source: Region | int | str | Path,
987 cache: Path | None = None,
988 head_bytes: int = 16384,
989) -> pd.DataFrame:
990 """PS extents of a region's mapprojected input images (``LROC-NAC-maps/*_GLDIMGP_*.tif``).
992 Only the start of each member is read (remotely: about 2 small HTTP
993 Range requests per image), so this is a cheap way to find the images
994 over a crop window. The maps are cropped to the valid pixels' bounding
995 box within the ROI, so a diagonal strip's extent over-covers.
997 Parameters
998 ----------
999 source : Path, Region, int or str
1000 A downloaded zip or a region (remote).
1001 cache : Path, optional
1002 JSON file of the result, written once and reused.
1003 head_bytes : int
1004 Bytes decompressed from each member.
1006 Returns
1007 -------
1008 pandas.DataFrame
1009 ``product_id``, ``x0``, ``y0``, ``x1``, ``y1`` (PS metres).
1010 """
1011 cols = ["product_id", "x0", "y0", "x1", "y1"]
1012 if cache is not None and Path(cache).is_file():
1013 return pd.DataFrame(json.loads(Path(cache).read_text()), columns=cols)
1014 rows = []
1015 with _open_zip(source, buffer_size=4 << 10) as zf:
1016 for info in zf.infolist():
1017 pid = _member_pid(info.filename)
1018 if pid is None or "_GLDIMGP_" not in info.filename or _is_junk(info.filename):
1019 continue
1020 rows.append([pid, *parse_tiff_extent(_member_head(zf, info, head_bytes)).bounds])
1021 if cache is not None: 1021 ↛ 1024line 1021 didn't jump to line 1024 because the condition on line 1021 was always true
1022 Path(cache).parent.mkdir(parents=True, exist_ok=True)
1023 Path(cache).write_text(json.dumps(rows))
1024 return pd.DataFrame(rows, columns=cols)
1027def images_covering(extents: pd.DataFrame, box: PsBox, min_fraction: float = 0.0) -> list[str]:
1028 """Product IDs whose extent overlaps *box* by more than *min_fraction* of the box area."""
1029 ix = np.clip(np.minimum(extents.x1, box.x1) - np.maximum(extents.x0, box.x0), 0, None)
1030 iy = np.clip(np.minimum(extents.y1, box.y1) - np.maximum(extents.y0, box.y0), 0, None)
1031 frac = ix * iy / (box.width * box.height)
1032 return list(extents.product_id[(frac > min_fraction) & (ix * iy > 0)])
1035def published_cameras(
1036 region: Region | int | str,
1037 out_dir: Path,
1038 images: Iterable[str] | None = None,
1039 window: PsBox | None = None,
1040 source: Path | None = None,
1041 force: bool = False,
1042) -> dict[str, Path]:
1043 """Extract the Bertone et al. bundle-adjusted NAC cameras (``LROC-NAC-adjust/*_GLDACAM_*.json``).
1045 Each file is an ASP / USGS CSM linescan model state (model-name line,
1046 then JSON) with the bundle adjustment and LOLA registration baked in.
1047 It describes the *full* image ISIS-``reduce``-d by an integer factor to
1048 about 5 m/px (see [`published_reduce_factor`][lunar_toolkit.reference.published_reduce_factor] and
1049 [`published_camera_native`][lunar_toolkit.reference.published_camera_native]).
1051 Parameters
1052 ----------
1053 region : Region, int or str
1054 Anything [`get_region`][lunar_toolkit.reference.get_region] accepts.
1055 out_dir : Path
1056 Destination; files land in ``out_dir/LROC-NAC-adjust/``.
1057 images : iterable of str, optional
1058 Only these product IDs.
1059 window : PsBox, optional
1060 Only images whose mapprojected extent overlaps this box
1061 ([`published_extents`][lunar_toolkit.reference.published_extents], cached in ``out_dir``).
1062 source : Path, optional
1063 A downloaded zip instead of the remote Zenodo file.
1064 force : bool
1065 Re-extract existing files.
1067 Returns
1068 -------
1069 dict
1070 ``product_id -> camera path`` (sorted by product ID).
1071 """
1072 reg = get_region(region)
1073 src: Region | Path = Path(source) if source is not None else reg
1074 cams = published_members(src)
1075 if cams.empty: 1075 ↛ 1076line 1075 didn't jump to line 1076 because the condition on line 1075 was never true
1076 raise KeyError(f"Region {reg.id} {reg.name} has no published cameras (LROC-NAC-adjust/)")
1077 wanted = set(cams.product_id)
1078 if images is not None:
1079 wanted &= {str(p).upper() for p in images}
1080 if window is not None:
1081 ext = published_extents(src, cache=Path(out_dir) / f"gldimgp_extents_{reg.id}.json")
1082 wanted &= set(images_covering(ext, window))
1083 names = sorted(n.rsplit("/", 1)[-1] for n in cams.name[cams.product_id.isin(wanted)])
1084 if not names:
1085 return {}
1086 paths = extract_sdem(src, names, Path(out_dir), force=force)
1087 return dict(sorted((p.name.rsplit("_", 1)[-1].removesuffix(".json"), p) for p in paths))
1090def published_reduce_factor(camera: CsmCamera | Path) -> int:
1091 """``reduce`` factor of a published camera: NAC focal length / its (scaled) ``m_focalLength``."""
1092 cam = camera if isinstance(camera, CsmCamera) else read_camera(camera)
1093 side = "LEFT" if str(cam.data.get("m_sensorName", "")).upper().endswith("LEFT") else "RIGHT"
1094 return max(1, round(NAC_FOCAL_MM[side] / float(cam.data["m_focalLength"])))
1097def published_camera_native(camera: CsmCamera | Path, native_size: tuple[int, int]) -> CsmCamera:
1098 """Published camera re-mapped to the full-resolution image (``(samples, lines)`` of the EDR / cube).
1100 The published states keep the detector summing and divide the focal
1101 length and detector-sample origin by the reduce factor *f*; their
1102 ``m_nLines`` / ``m_nSamples`` are the native size divided by ``f**2``
1103 (not ``f``), which only affects usgscsm's initial ground-to-image guess.
1104 The result is a state for the native image with consistent sizes; use
1105 [`lunar_toolkit.csm.camera_for_cube`][lunar_toolkit.csm.camera_for_cube] or
1106 [`lunar_toolkit.csm.isis_reduce_with_csm`][lunar_toolkit.csm.isis_reduce_with_csm] to get any crop or
1107 reduction of it.
1108 """
1109 cam = camera if isinstance(camera, CsmCamera) else read_camera(camera)
1110 f = published_reduce_factor(cam)
1111 out = transform_camera(cam, 1.0 / f, 0.0, 1.0 / f, 0.0, native_size)
1112 if "m_gsd" in cam.data: # already the native ground sample distance 1112 ↛ 1115line 1112 didn't jump to line 1115 because the condition on line 1112 was always true
1113 out.data["m_gsd"] = cam.data["m_gsd"]
1114 out.data["m_halfSwath"] = cam.data["m_gsd"] * native_size[0] / 2.0
1115 return out