Coverage for lunar_toolkit/lroc/cumindex.py: 75%
563 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 02:58 +0000
1"""LROC PDS CUMINDEX as a local DuckDB / Parquet database."""
3from __future__ import annotations
5import datetime as _dt
6import json
7import logging
8import re
9import urllib.error
10import urllib.parse
11import urllib.request
12import xml.etree.ElementTree as ET
13from collections.abc import Callable, Iterable
14from pathlib import Path
15from typing import TYPE_CHECKING, Any
17from ..constants import (
18 DATASETS,
19 DEFAULT_CACHE_DIR,
20 FOOTPRINT_CORNERS,
21 INDEX_DIR,
22 MOON_RADIUS_KM,
23 NAC_PRODUCT_TYPES,
24 NAC_SUFFIXES,
25 PDS_STORE_URL,
26 POLAR_LAT_THRESHOLD,
27 VOLINDEX_LBL_FILENAME,
28 VOLINDEX_TAB_FILENAME,
29 find_parquet,
30 is_nac_id,
31 parse_nac_id,
32 pds_url_from_file_spec,
33)
34from ..downloader import _USER_AGENT, PdsDownloader
35from ..geometry import aoi_region, filter_aoi, load_aoi, query_region
37if TYPE_CHECKING: # pragma: no cover
38 import pandas as pd
40 from .nac import NacImageCollection
42log = logging.getLogger(__name__)
44# Bump when the Parquet layout changes in a way readers should know about.
45SCHEMA_VERSION = 2
47# Footprint columns derived at build time (all prefixed FP_).
48FOOTPRINT_COLUMNS = [
49 "FP_MIN_LAT",
50 "FP_MAX_LAT",
51 "FP_LON_WEST",
52 "FP_LON_EAST",
53 "FP_LON_SPAN",
54 "FP_CONTAINS_POLE",
55 "FP_PS_HEMI",
56 "FP_PS_XMIN",
57 "FP_PS_XMAX",
58 "FP_PS_YMIN",
59 "FP_PS_YMAX",
60]
63class CumIndex:
64 """Query the LROC PDS CUMINDEX via DuckDB on a local Parquet database.
66 Parameters
67 ----------
68 dataset : str
69 Which PDS dataset to use: ``"edr"`` (default), ``"cdr"``, or ``"rdr"``.
70 data_dir : Path, optional
71 Directory containing ``cumindex_<dataset>.parquet``. Defaults to
72 ``LUNAR_DATA_DIR``, else the user data dir.
73 cache_dir : Path
74 Where raw CUMINDEX/INDEX .TAB / .LBL downloads are cached.
75 Defaults to ``~/.cache/lunar_toolkit/``.
76 release : str, optional
77 PDS volume to read labels / index tables from (e.g. ``"LROLRC_0066C"``).
78 Defaults to the release recorded in the build sidecar, else the
79 fallback constant. ``build_db()`` auto-discovers the newest volume.
80 """
82 def __init__(
83 self,
84 dataset: str = "edr",
85 data_dir: Path | None = None,
86 cache_dir: Path = DEFAULT_CACHE_DIR,
87 release: str | None = None,
88 ) -> None:
89 if dataset not in DATASETS: 89 ↛ 90line 89 didn't jump to line 90 because the condition on line 89 was never true
90 raise ValueError(f"Unknown dataset {dataset!r}. Valid options: {list(DATASETS)}")
91 self._dataset = dataset
92 self._ds = DATASETS[dataset]
93 self._explicit_data_dir = Path(data_dir) if data_dir is not None else None
94 self._cache_dir = Path(cache_dir)
95 self._cache = PdsDownloader(self._cache_dir)
96 self._release = _norm_volume(release) if release else None
97 self._schema: pd.DataFrame | None = None
98 self._meta: dict | None = None
100 # Dataset URLs
101 @property
102 def dataset(self) -> str:
103 return self._dataset
105 @property
106 def release(self) -> str:
107 """Volume used for label / index URLs (no trailing slash)."""
108 if self._release: 108 ↛ 110line 108 didn't jump to line 110 because the condition on line 108 was always true
109 return self._release
110 side = self.build_info()
111 if side and side.get("release"):
112 return side["release"]
113 return _norm_volume(self._ds["release"])
115 def _index_base_url(self, volume: str | None = None) -> str:
116 vol = _norm_volume(volume) if volume else self.release
117 return self._ds["base_url"] + vol + "/" + self._ds["index_dir"]
119 @property
120 def lbl_url(self) -> str:
121 return self._index_base_url() + self._ds["lbl_filename"]
123 @property
124 def tab_url(self) -> str:
125 return self._index_base_url() + self._ds["tab_filename"]
127 def volume_index_urls(self, volume: str) -> tuple[str, str]:
128 """``(INDEX.LBL, INDEX.TAB)`` URLs for a single PDS volume."""
129 base = self._index_base_url(volume)
130 return base + VOLINDEX_LBL_FILENAME, base + VOLINDEX_TAB_FILENAME
132 # Parquet path / sidecar
133 @property
134 def parquet_path(self) -> Path:
135 """Path to the Parquet database for this dataset."""
136 return find_parquet(self._ds["parquet"], self._explicit_data_dir)
138 @property
139 def meta_path(self) -> Path:
140 """Build-metadata sidecar next to the Parquet (``*.meta.json``)."""
141 p = self.parquet_path
142 return p.with_name(p.stem + ".meta.json")
144 def build_info(self) -> dict | None:
145 """Return the build sidecar as a dict, or None if absent."""
146 mp = self.meta_path
147 if not mp.exists(): 147 ↛ 148line 147 didn't jump to line 148 because the condition on line 147 was never true
148 return None
149 try:
150 return json.loads(mp.read_text())
151 except (OSError, json.JSONDecodeError):
152 return None
154 def _ensure_parquet(self) -> Path:
155 p = self.parquet_path
156 if not p.exists():
157 raise FileNotFoundError(
158 f"Database not found: {p}\n"
159 f"Build it with: lunar build-db {self._dataset}\n"
160 f"Or set LUNAR_DATA_DIR to a directory that contains {self._ds['parquet']}"
161 )
162 return p
164 # Schema / metadata (from .LBL)
165 @property
166 def schema(self) -> pd.DataFrame:
167 """Column definitions from CUMINDEX.LBL (lazy, downloads to cache)."""
168 if self._schema is None:
169 lbl_text = self._read_lbl(self.lbl_url)
170 self._schema = _parse_pds_lbl_columns(lbl_text)
171 self._meta = _parse_lbl_meta(lbl_text)
172 return self._schema
174 @property
175 def meta(self) -> dict:
176 """Top-level scalar fields from CUMINDEX.LBL (lazy)."""
177 if self._meta is None:
178 self._meta = _parse_lbl_meta(self._read_lbl(self.lbl_url))
179 return self._meta
181 def _read_lbl(self, url: str) -> str:
182 return self._cache.ensure_cached(url).read_text(encoding="latin-1")
184 # Volume discovery
185 def discover_volumes(self) -> list[str]:
186 """List this dataset's PDS volumes (sorted oldest → newest).
188 Uses the S3 ``ListObjectsV2`` listing exposed by the PDS store; falls
189 back to probing ``INDEX/INDEX.LBL`` of successive volume numbers.
190 """
191 try:
192 vols = list_store_volumes(self._ds["base_url"], self._ds["volume_re"])
193 if vols:
194 return vols
195 except (OSError, ET.ParseError, ValueError) as exc:
196 log.warning("Volume listing failed (%s); probing instead", exc)
197 start = self.build_info() or {}
198 return probe_volumes(
199 self._ds["base_url"],
200 start.get("release") or _norm_volume(self._ds["release"]),
201 )
203 def latest_release(self) -> str:
204 """Newest volume on the server (fallback: the constant release)."""
205 try:
206 vols = self.discover_volumes()
207 except Exception as exc: # network down etc.
208 log.warning("Could not discover volumes (%s)", exc)
209 vols = []
210 return vols[-1] if vols else _norm_volume(self._ds["release"])
212 # Build database
213 def build_db(
214 self,
215 force: bool = False,
216 release: str | None = None,
217 volume_only: bool = False,
218 memory_limit: str = "4GB",
219 threads: int | None = None,
220 wkt: bool = False,
221 on_stage: Callable[[str], None] | None = None,
222 ) -> Path:
223 """Download the PDS index table and build the Parquet database.
225 Parameters
226 ----------
227 force : bool
228 Rebuild even if the Parquet already exists.
229 release : str, optional
230 Volume to build from (e.g. ``"LROLRC_0066C"``). Default: newest
231 volume discovered on the server, else the fallback constant.
232 volume_only : bool
233 Build from that volume's own ``INDEX/INDEX.TAB`` (one release,
234 ~20 k rows) instead of the cumulative ``CUMINDEX.TAB`` (~3.8 GB).
235 memory_limit : str
236 DuckDB ``memory_limit``; sorts spill to a temp dir under the cache.
237 threads : int, optional
238 DuckDB worker threads (default: all cores).
239 wkt : bool
240 Also store a ``FP_WKT`` footprint polygon column.
241 on_stage : callable, optional
242 Called with ``"download"`` and ``"build"`` as those stages start.
244 Returns
245 -------
246 Path to the written Parquet file.
247 """
248 p = self.parquet_path
249 if p.exists() and not force: 249 ↛ 250line 249 didn't jump to line 250 because the condition on line 249 was never true
250 log.info("Database already exists: %s (force=True rebuilds, update_db() appends new volumes)", p)
251 return p
253 p.parent.mkdir(parents=True, exist_ok=True)
254 clean_stale_build_files(p.parent)
256 if on_stage: 256 ↛ 257line 256 didn't jump to line 257 because the condition on line 256 was never true
257 on_stage("download")
258 vol = _norm_volume(release) if release else self.latest_release()
259 self._release = vol
260 self._schema = self._meta = None
261 if volume_only: 261 ↛ 264line 261 didn't jump to line 264 because the condition on line 261 was always true
262 lbl_url, tab_url = self.volume_index_urls(vol)
263 else:
264 lbl_url, tab_url = self.lbl_url, self.tab_url
266 lbl_text = self._read_lbl(lbl_url)
267 meta = _parse_lbl_meta(lbl_text)
268 columns = parse_lbl_columns(lbl_text)
269 tab_path = self._cache.ensure_cached(tab_url, expected_md5=meta.get("md5_checksum"))
271 if on_stage: 271 ↛ 272line 271 didn't jump to line 272 because the condition on line 271 was never true
272 on_stage("build")
273 n = build_parquet(
274 [(tab_path, columns)],
275 p,
276 temp_dir=self._cache_dir / "duckdb_tmp",
277 memory_limit=memory_limit,
278 threads=threads,
279 wkt=wkt,
280 expected_rows=meta.get("file_records") or meta.get("rows"),
281 )
282 now = _utcnow()
283 self._write_sidecar(
284 {
285 "schema_version": SCHEMA_VERSION,
286 "dataset": self._dataset,
287 "data_set_id": self._ds["data_set_id"],
288 "source": "volume_index" if volume_only else "cumindex",
289 "release": vol,
290 "lbl_url": lbl_url,
291 "tab_url": tab_url,
292 "md5": meta.get("md5_checksum"),
293 "rows": n,
294 "wkt": wkt,
295 "volumes": {
296 vol: {
297 "kind": "volume_index" if volume_only else "cumindex",
298 "rows": n,
299 "md5": meta.get("md5_checksum"),
300 "added_at": now,
301 }
302 },
303 "built_at": now,
304 "updated_at": now,
305 }
306 )
307 return p
309 def update_db(
310 self,
311 volumes: Iterable[str] | None = None,
312 memory_limit: str = "4GB",
313 threads: int | None = None,
314 on_stage: Callable[[str], None] | None = None,
315 ) -> list[str]:
316 """Append per-volume ``INDEX.TAB`` files newer than the built release.
318 Parameters
319 ----------
320 volumes : iterable of str, optional
321 Explicit volumes to add. Default: every discovered volume newer
322 than the sidecar's ``release`` (or ``max(VOLUME_ID)`` in the
323 Parquet when no sidecar exists).
324 on_stage : callable, optional
325 As in [`build_db`][lunar_toolkit.lroc.cumindex.CumIndex.build_db].
327 Returns
328 -------
329 The list of volumes that were appended (empty if already up to date).
330 Rows are de-duplicated on ``PRODUCT_ID`` (newest volume wins).
331 """
332 p = self._ensure_parquet()
333 if on_stage: 333 ↛ 334line 333 didn't jump to line 334 because the condition on line 333 was never true
334 on_stage("download")
335 side = self.build_info() or {}
336 current = side.get("release") or self._max_volume_in_parquet()
337 if volumes is None: 337 ↛ 344line 337 didn't jump to line 344 because the condition on line 337 was always true
338 if current is None: 338 ↛ 339line 338 didn't jump to line 339 because the condition on line 338 was never true
339 raise RuntimeError("Cannot determine the built release; rebuild with --force")
340 found = self.discover_volumes()
341 ck = _volume_key(current)
342 new_vols = [v for v in found if _volume_key(v) > ck]
343 else:
344 new_vols = sorted({_norm_volume(v) for v in volumes}, key=_volume_key)
346 if not new_vols:
347 log.info("%s is up to date (release %s)", p.name, current)
348 return []
350 log.info("Appending %d volume(s) to %s: %s", len(new_vols), p.name, ", ".join(new_vols))
351 clean_stale_build_files(p.parent)
352 sources: list[tuple[Path, list[dict]]] = []
353 added: dict[str, dict] = {}
354 for vol in new_vols:
355 lbl_url, tab_url = self.volume_index_urls(vol)
356 lbl_text = self._read_lbl(lbl_url)
357 meta = _parse_lbl_meta(lbl_text)
358 tab = self._cache.ensure_cached(tab_url, expected_md5=meta.get("md5_checksum"))
359 sources.append((tab, parse_lbl_columns(lbl_text)))
360 added[vol] = {
361 "kind": "volume_index",
362 "rows": meta.get("file_records") or meta.get("rows"),
363 "md5": meta.get("md5_checksum"),
364 "added_at": _utcnow(),
365 }
367 if on_stage: 367 ↛ 368line 367 didn't jump to line 368 because the condition on line 367 was never true
368 on_stage("build")
369 tmp = p.with_name(p.stem + ".update.parquet")
370 n = build_parquet(
371 sources,
372 tmp,
373 base_parquet=p,
374 temp_dir=self._cache_dir / "duckdb_tmp",
375 memory_limit=memory_limit,
376 threads=threads,
377 wkt=bool(side.get("wkt")),
378 )
379 tmp.replace(p)
381 vols = dict(side.get("volumes") or {})
382 vols.update(added)
383 side.update(
384 {
385 "schema_version": SCHEMA_VERSION,
386 "dataset": self._dataset,
387 "data_set_id": self._ds["data_set_id"],
388 "release": new_vols[-1],
389 "rows": n,
390 "volumes": vols,
391 "updated_at": _utcnow(),
392 }
393 )
394 side.setdefault("built_at", side["updated_at"])
395 self._write_sidecar(side)
396 self._release = new_vols[-1]
397 return new_vols
399 def _write_sidecar(self, info: dict) -> None:
400 self.meta_path.write_text(json.dumps(info, indent=2, sort_keys=True))
402 def _max_volume_in_parquet(self) -> str | None:
403 import duckdb
405 row = duckdb.sql(
406 f"SELECT max(VOLUME_ID) FROM read_parquet({_sql_str(self._ensure_parquet())})"
407 ).fetchone()
408 return row[0] if row and row[0] else None
410 # DuckDB query interface
411 def _connect(self):
412 import duckdb
414 con = duckdb.connect()
415 _register_macros(con)
416 con.execute(f"CREATE VIEW cumindex AS SELECT * FROM read_parquet({_sql_str(self._ensure_parquet())})")
417 return con
419 def query(self, sql: str) -> pd.DataFrame:
420 """Run arbitrary DuckDB SQL against the Parquet database.
422 Use the alias ``cumindex`` as the table name:
424 ci.query("SELECT PRODUCT_ID, EMISSION_ANGLE FROM cumindex LIMIT 10")
425 """
426 con = self._connect()
427 try:
428 return con.execute(sql).df()
429 finally:
430 con.close()
432 def lookup(
433 self,
434 product_ids: list[str],
435 columns: list[str] | None = None,
436 ) -> pd.DataFrame:
437 """Return metadata rows for the given product IDs.
439 Queries the Parquet via DuckDB — no full dataset loaded into RAM.
441 Parameters
442 ----------
443 product_ids : list[str]
444 e.g. ``["M1524360231LE", "M1512659523RE"]``
445 columns : list[str], optional
446 Subset of columns to return. ``PRODUCT_ID`` and
447 ``FILE_SPECIFICATION_NAME`` are always included.
448 """
449 ids_clean = [p.strip().upper() for p in product_ids]
450 if not ids_clean:
451 return self._enrich(self.query("SELECT * FROM cumindex LIMIT 0"))
452 ids_sql = ", ".join(_sql_str(i) for i in ids_clean)
453 col_sql = _col_select(columns)
454 sql = f"SELECT {col_sql} FROM cumindex WHERE PRODUCT_ID IN ({ids_sql})"
455 return self._enrich(self.query(sql))
457 def search(
458 self,
459 lat: tuple[float, float] | None = None,
460 lon: tuple[float, float] | None = None,
461 point: tuple[float, float] | None = None,
462 radius_km: float = 0.0,
463 emission_angle: tuple[float, float] | None = None,
464 incidence_angle: tuple[float, float] | None = None,
465 phase_angle: tuple[float, float] | None = None,
466 nac_only: bool = False,
467 columns: list[str] | None = None,
468 wkt: bool = False,
469 limit: int | None = None,
470 aoi=None,
471 ) -> pd.DataFrame:
472 """Return images whose *footprint* intersects a region.
474 Region: either a lat/lon box (``lat=(min, max)``, ``lon=(west, east)``;
475 ``west > east`` crosses the 0/360 seam, longitudes may be given in
476 −180..180 or 0..360, omitted ``lon`` = all longitudes), a circle
477 (``point=(lat, lon)``, ``radius_km``) or an area of interest
478 (``aoi=``: GeoJSON / WKT path, WKT string or
479 [`Aoi`][lunar_toolkit.geometry.Aoi]).
481 Matching is bounding-box based and conservative (may include a few
482 false positives, never drops a true hit for convex-ish footprints):
483 the image's lat range and seam-aware longitude arc must overlap the
484 query's, and for images poleward of ±60° their polar stereographic
485 x/y box must overlap the query's projected box. Footprints that
486 enclose a pole are treated as spanning all longitudes. With ``aoi``
487 the bbox matches are then refined exactly: each corner footprint
488 polygon must intersect the AOI polygon (in polar stereographic when
489 the AOI lies poleward of ±60°), and ``limit`` applies after that.
491 Example
492 -------
493 ::
495 df = ci.search(lat=(-90, -85), nac_only=True, incidence_angle=(80, 95))
496 df = ci.search(lat=(-10, 10), lon=(355, 5)) # across 0°
497 df = ci.search(point=(-89.67, 129.8), radius_km=15) # Shackleton
498 df = ci.search(aoi="aoi.geojson") # QuickMap polygon
499 """
500 if aoi is not None:
501 if lat is not None or lon is not None or point is not None:
502 raise ValueError("Pass either aoi or lat/lon/point, not both")
503 aoi = load_aoi(aoi)
504 region = aoi_region(aoi)
505 else:
506 region = query_region(lat=lat, lon=lon, point=point, radius_km=radius_km)
507 conditions = [region_sql(region)]
508 conditions += _angle_conditions(emission_angle, incidence_angle, phase_angle)
509 if nac_only:
510 conditions.append(_NAC_ONLY_SQL)
512 cols = _col_select(columns, extra=FOOTPRINT_COLUMNS)
513 if wkt or aoi is not None:
514 cols += f", {footprint_wkt_sql(self._parquet_columns())} AS FP_WKT_QUERY"
515 sql = f"SELECT {cols} FROM cumindex WHERE " + " AND ".join(conditions)
516 if limit and aoi is None: 516 ↛ 517line 516 didn't jump to line 517 because the condition on line 516 was never true
517 sql += f" LIMIT {int(limit)}"
518 df = self._enrich(self.query(sql))
519 if aoi is not None:
520 df = filter_aoi(df, aoi, wkt_column="FP_WKT_QUERY")
521 if limit:
522 df = df.head(int(limit))
523 if "FP_WKT_QUERY" in df.columns:
524 query_wkt = df.pop("FP_WKT_QUERY")
525 if wkt and "FP_WKT" not in df.columns:
526 df["FP_WKT"] = query_wkt
527 return df
529 def _parquet_columns(self) -> list[str]:
530 import duckdb
532 rel = duckdb.sql(f"SELECT * FROM read_parquet({_sql_str(self._ensure_parquet())}) LIMIT 0")
533 return list(rel.columns)
535 def get_images(self, product_ids: list[str]) -> NacImageCollection:
536 """Look up product IDs and return a NacImageCollection.
538 Single DuckDB query; each NacImage is pre-populated with its metadata
539 row so subsequent property accesses require zero additional DB calls.
541 Parameters
542 ----------
543 product_ids : list[str]
544 e.g. ``["M1524360231LE", "M1512659523RE"]``
546 Returns
547 -------
548 NacImageCollection
549 Contains only the IDs that were found in the index, in input order.
550 """
551 from .nac import NacImage, NacImageCollection # deferred — avoids circular import
553 df = self.lookup(product_ids)
554 row_map = {str(row["PRODUCT_ID"]): row for _, row in df.iterrows()}
555 ids = [pid.strip().upper() for pid in product_ids]
556 return NacImageCollection(
557 [NacImage(pid, cumindex=self, _row=row_map[pid]) for pid in ids if pid in row_map]
558 )
560 def load(self, nrows: int | None = None) -> pd.DataFrame:
561 """Load the full index into memory as a DataFrame.
563 For most workflows prefer ``lookup()`` or ``search()``.
564 Pass ``nrows`` for a quick preview that needs no Parquet: only the
565 label and the first ``nrows`` records of the TAB are fetched (HTTP
566 Range request), then parsed with the same DuckDB extraction as
567 ``build_db()``.
568 """
569 if nrows is not None: 569 ↛ 591line 569 didn't jump to line 591 because the condition on line 569 was always true
570 import duckdb
572 lbl_text = self._read_lbl(self.lbl_url)
573 meta = _parse_lbl_meta(lbl_text)
574 columns = parse_lbl_columns(lbl_text)
575 row_bytes = meta.get("record_bytes") or meta.get("row_bytes") or 903
576 full = self._cache.cache_path(self.tab_url)
577 if full.exists(): 577 ↛ 580line 577 didn't jump to line 580 because the condition on line 577 was always true
578 tab = full
579 else:
580 tab = full.with_name(full.stem + f".head{int(nrows)}" + full.suffix)
581 if not tab.exists():
582 _fetch_head(self.tab_url, tab, int(nrows) * int(row_bytes))
583 con = duckdb.connect()
584 try:
585 _register_macros(con)
586 sql = _select_sql([(tab, columns)], wkt=False)
587 return con.execute(f"{sql} LIMIT {int(nrows)}").df()
588 finally:
589 con.close()
591 import pyarrow.parquet as pq
593 return pq.read_table(self._ensure_parquet()).to_pandas()
595 # Enrichment
596 def _enrich(self, df: pd.DataFrame) -> pd.DataFrame:
597 """Add PDS_URL and derived convenience columns to a query result."""
598 import pandas as pd
600 if df.empty:
601 return df
603 df = df.copy()
605 if "FILE_SPECIFICATION_NAME" in df.columns:
606 df["PDS_URL"] = df["FILE_SPECIFICATION_NAME"].map(
607 lambda p: pds_url_from_file_spec(p) if isinstance(p, str) and p else None
608 )
609 df["LABEL_URL"] = df["PDS_URL"].map(
610 lambda u: re.sub(r"\.[A-Za-z0-9]+$", ".xml", u) if u else None
611 )
612 if "CENTER_LATITUDE" in df.columns: 612 ↛ 614line 612 didn't jump to line 614 because the condition on line 612 was always true
613 df["LAT"] = pd.to_numeric(df["CENTER_LATITUDE"], errors="coerce")
614 if "CENTER_LONGITUDE" in df.columns: 614 ↛ 616line 614 didn't jump to line 616 because the condition on line 614 was always true
615 df["LON"] = pd.to_numeric(df["CENTER_LONGITUDE"], errors="coerce")
616 if "START_TIME" in df.columns:
617 df["DATE"] = pd.to_datetime(df["START_TIME"], errors="coerce").dt.date
618 if "PRODUCT_ID" in df.columns: 618 ↛ 626line 618 didn't jump to line 626 because the condition on line 618 was always true
619 pid = df["PRODUCT_ID"].astype(str)
620 nac = [parse_nac_id(p) if is_nac_id(p) else None for p in pid]
621 df["IS_NAC"] = [t is not None for t in nac]
622 df["IS_WAC"] = pid.str.contains(r"[CMUV][EC]$") & ~df["IS_NAC"]
623 df["NAC_SIDE"] = [t[1] if t else None for t in nac]
624 df["PRODUCT_TYPE"] = pid.str[-1:].map(NAC_PRODUCT_TYPES).where(df["IS_NAC"] | df["IS_WAC"])
626 return df.reset_index(drop=True)
628 def to_observations(self, df: pd.DataFrame) -> pd.DataFrame:
629 """Convert query results to the common multi-instrument schema."""
630 from ..observation import observations_from_cumindex
632 return observations_from_cumindex(df, dataset=self._dataset)
635# Parquet builder (DuckDB-native)
636def build_parquet(
637 sources: list[tuple[Path, list[dict]]],
638 parquet_path: Path,
639 base_parquet: Path | None = None,
640 temp_dir: Path | None = None,
641 memory_limit: str = "4GB",
642 threads: int | None = None,
643 wkt: bool = False,
644 expected_rows: int | None = None,
645) -> int:
646 """Build a sorted, ZSTD-compressed Parquet from fixed-width PDS TAB file(s).
648 Parameters
649 ----------
650 sources : list of (tab_path, lbl_columns)
651 Each TAB is parsed with its own label layout (record lengths differ
652 between volumes) and the results are combined ``UNION ALL BY NAME``.
653 base_parquet : Path, optional
654 Existing database to merge with (incremental update). Rows are then
655 de-duplicated on ``PRODUCT_ID``, keeping the newest ``VOLUME_ID``.
657 Returns
658 -------
659 int
660 Number of rows written.
662 Notes
663 -----
664 Everything happens in one ``COPY (...) TO`` statement; DuckDB spills
665 the sort to *temp_dir*.
666 """
667 import duckdb
669 parquet_path = Path(parquet_path)
670 parquet_path.parent.mkdir(parents=True, exist_ok=True)
671 tmp_out = parquet_path.with_name(parquet_path.name + ".tmp")
673 con = duckdb.connect()
674 try:
675 if temp_dir is not None: 675 ↛ 678line 675 didn't jump to line 678 because the condition on line 675 was always true
676 Path(temp_dir).mkdir(parents=True, exist_ok=True)
677 con.execute(f"SET temp_directory={_sql_str(temp_dir)}")
678 con.execute(f"SET memory_limit={_sql_str(memory_limit)}")
679 con.execute("SET preserve_insertion_order=false")
680 if threads: 680 ↛ 681line 680 didn't jump to line 681 because the condition on line 680 was never true
681 con.execute(f"SET threads={int(threads)}")
682 _register_macros(con)
684 select = _select_sql(sources, wkt=wkt, base_parquet=base_parquet)
685 names = ", ".join(Path(t).name for t, _ in sources)
686 base = f" + {Path(base_parquet).name}" if base_parquet else ""
687 log.info("Building %s from %s%s", parquet_path.name, names, base)
688 if expected_rows:
689 log.info("Expected rows: %s", f"{expected_rows:,}")
691 con.execute(
692 f"COPY ({select} ORDER BY CENTER_LATITUDE NULLS LAST, PRODUCT_ID) "
693 f"TO {_sql_str(tmp_out)} (FORMAT PARQUET, COMPRESSION ZSTD)"
694 )
695 row = con.execute(f"SELECT count(*) FROM read_parquet({_sql_str(tmp_out)})").fetchone()
696 n = row[0] if row else 0
697 except BaseException:
698 if tmp_out.exists():
699 tmp_out.unlink()
700 raise
701 finally:
702 con.close()
704 tmp_out.replace(parquet_path)
705 size_mb = parquet_path.stat().st_size / 1e6
706 log.info("Done — %s rows → %s (%.1f MB)", f"{n:,}", parquet_path.name, size_mb)
707 if expected_rows and n != expected_rows and base_parquet is None: 707 ↛ 708line 707 didn't jump to line 708 because the condition on line 707 was never true
708 log.warning("Expected %s rows, got %s", f"{expected_rows:,}", f"{n:,}")
709 return int(n)
712def _tab_source_sql(tab_path: Path) -> str:
713 """DuckDB table function reading a fixed-width TAB as one VARCHAR column."""
714 return (
715 f"read_csv({_sql_str(tab_path)}, columns={{'line': 'VARCHAR'}}, "
716 f"delim=chr(31), quote='', escape='', header=false, auto_detect=false, "
717 f"encoding='latin-1')"
718 )
721def _select_sql(
722 sources: list[tuple[Path, list[dict]]],
723 wkt: bool = False,
724 base_parquet: Path | None = None,
725) -> str:
726 """Full SELECT (extraction → union → dedup → footprint columns)."""
727 parts = []
728 all_names: list[str] = []
729 for tab, cols in sources:
730 exprs = ",\n ".join(f"{column_sql(c)} AS {_ident(c['name'])}" for c in cols)
731 parts.append(f"SELECT {exprs}\nFROM {_tab_source_sql(tab)}\nWHERE length(trim(line)) > 0")
732 all_names += [c["name"] for c in cols if c["name"] not in all_names]
733 raw = "\nUNION ALL BY NAME\n".join(f"({p})" for p in parts)
734 if base_parquet is not None:
735 # Incremental update: new rows win. An anti-join against the (small)
736 # set of new product IDs is far cheaper than a window over all rows.
737 raw = (
738 f"WITH new_rows AS (SELECT * FROM ({raw}) QUALIFY row_number() OVER "
739 f"(PARTITION BY PRODUCT_ID ORDER BY VOLUME_ID DESC NULLS LAST) = 1)\n"
740 f"SELECT * FROM (SELECT COLUMNS(c -> NOT starts_with(c, 'FP_')) "
741 f"FROM read_parquet({_sql_str(base_parquet)}) "
742 f"WHERE PRODUCT_ID NOT IN (SELECT PRODUCT_ID FROM new_rows)) "
743 f"UNION ALL BY NAME SELECT * FROM new_rows"
744 )
745 fp = footprint_sql(all_names, wkt=wkt)
746 return f"SELECT *{', ' + fp if fp else ''} FROM ({raw}) AS t"
749def column_sql(col: dict) -> str:
750 """SQL expression extracting one LBL column from the ``line`` VARCHAR.
752 Bytes → ``substr(line, START_BYTE, BYTES)``, trimmed of blanks and
753 quotes, typed by ``DATA_TYPE`` (falling back to name hints) with
754 ``TRY_CAST``; ``NOT_APPLICABLE_CONSTANT`` / ``MISSING_CONSTANT`` → NULL.
755 """
756 start, nbytes = int(col["start_byte"]), int(col["bytes"])
757 raw = f"nullif(trim(substr(line, {start}, {nbytes}), ' \"'), '')"
758 kind = _column_kind(col)
759 sentinels = [
760 str(col[k]).strip().strip('"')
761 for k in ("not_applicable_constant", "missing_constant", "unknown_constant")
762 if col.get(k) not in (None, "")
763 ]
765 if kind in ("BIGINT", "DOUBLE"):
766 expr = f"TRY_CAST({raw} AS {kind})"
767 for s in sentinels:
768 try:
769 v = float(s)
770 except ValueError:
771 continue
772 lit = str(int(v)) if kind == "BIGINT" and v.is_integer() else repr(v)
773 expr = f"nullif({expr}, {lit})"
774 return expr
775 if kind == "TIMESTAMP":
776 return f"TRY_CAST({raw} AS TIMESTAMP)"
778 expr = raw
779 for s in sentinels:
780 expr = f"nullif({expr}, {_sql_str(s)})"
781 if col["name"].upper() == "PRODUCT_ID":
782 expr = f"upper({expr})"
783 return expr
786_NUMERIC_HINTS = [
787 # Schema columns
788 "COLUMN_NUMBER",
789 "START_BYTE",
790 "BYTES",
791 # Image dimensions
792 "LINES",
793 "LINE_SAMPLES",
794 "SAMPLE_BITS",
795 # Angles / azimuths
796 "_ANGLE",
797 "_AZIMUTH",
798 "SOLAR_LONGITUDE",
799 # Spatial coordinates (all corners + center + subsolar + subspacecraft)
800 "_LATITUDE",
801 "_LONGITUDE",
802 # Distance / altitude
803 "ALTITUDE",
804 "DISTANCE",
805 # Resolution
806 "PIXEL_SCALE",
807 "SCALED_PIXEL",
808 "RESOLUTION",
809 # Timing
810 "EXPOSURE_DURATION",
811 # NAC instrument / temperatures
812 "TEMPERATURE",
813 "DAC_RESET",
814 "CHANNEL_A_OFFSET",
815 "CHANNEL_B_OFFSET",
816 "MODE_CODE",
817 "SELECT_CODE",
818 # Orbit
819 "ORBIT_NUMBER",
820]
822_TIME_COLUMNS = {"START_TIME", "STOP_TIME", "NAC_PREROLL_START_TIME", "PRODUCT_CREATION_TIME"}
825def _column_kind(col: dict) -> str:
826 """DuckDB type for an LBL column: LBL DATA_TYPE first, then name hints."""
827 name = col["name"].upper()
828 dtype = str(col.get("data_type") or "").upper()
829 if dtype in (
830 "ASCII_INTEGER",
831 "INTEGER",
832 "MSB_INTEGER",
833 "LSB_INTEGER",
834 "UNSIGNED_INTEGER",
835 "ASCII_UNSIGNED_INTEGER",
836 ):
837 return "BIGINT"
838 if dtype in ("ASCII_REAL", "REAL", "FLOAT", "IEEE_REAL", "PC_REAL"):
839 return "DOUBLE"
840 if dtype in ("TIME", "DATE", "ASCII_TIME", "ASCII_DATE_TIME") or name in _TIME_COLUMNS:
841 return "TIMESTAMP"
842 if dtype in ("", "UNKNOWN") and any(h in name for h in _NUMERIC_HINTS):
843 return "DOUBLE"
844 return "VARCHAR"
847# Footprints
848_MACROS = [
849 # longitude → [0, 360)
850 "CREATE OR REPLACE MACRO _lon360(x) AS ((x % 360.0) + 360.0) % 360.0",
851 # signed shortest longitude step a → b in [-180, 180)
852 "CREATE OR REPLACE MACRO _dlon(a, b) AS ((_lon360(b) - _lon360(a) + 540.0) % 360.0) - 180.0",
853 # polar stereographic radius (m) for hemisphere 'S' / 'N', true scale at the pole
854 f"CREATE OR REPLACE MACRO _ps_rho(lat, hemi) AS 2.0 * {MOON_RADIUS_KM * 1000.0!r} * "
855 "tan(pi() / 4.0 + (CASE WHEN hemi = 'S' THEN 1.0 ELSE -1.0 END) * radians(lat) / 2.0)",
856 "CREATE OR REPLACE MACRO _ps_x(lat, lon, hemi) AS _ps_rho(lat, hemi) * sin(radians(lon))",
857 "CREATE OR REPLACE MACRO _ps_y(lat, lon, hemi) AS "
858 "(CASE WHEN hemi = 'S' THEN 1.0 ELSE -1.0 END) * _ps_rho(lat, hemi) * cos(radians(lon))",
859 # does arc [w1, w1+s1] overlap arc [w2, w2+s2]? (degrees, s in [0, 360])
860 "CREATE OR REPLACE MACRO _arc_overlap(w1, s1, w2, s2) AS "
861 "(s1 >= 360 OR s2 >= 360 OR _lon360(w2 - w1) <= s1 OR _lon360(w1 - w2) <= s2)",
862]
865def _register_macros(con) -> None:
866 for m in _MACROS:
867 con.execute(m)
870def _corner_pairs(names: Iterable[str]) -> list[tuple[str, str]] | None:
871 have = set(names)
872 pairs = [(f"{c}_LATITUDE", f"{c}_LONGITUDE") for c in FOOTPRINT_CORNERS]
873 if all(a in have and b in have for a, b in pairs): 873 ↛ 875line 873 didn't jump to line 875 because the condition on line 873 was always true
874 return pairs
875 return None
878def footprint_sql(names: Iterable[str], wkt: bool = False) -> str:
879 """Derived ``FP_*`` column expressions over an extracted-row relation.
881 * ``FP_MIN_LAT`` / ``FP_MAX_LAT`` — corner (and centre) latitude range,
882 extended to ±90 when the footprint encloses a pole.
883 * ``FP_LON_WEST`` / ``FP_LON_EAST`` / ``FP_LON_SPAN`` — the shortest
884 longitude arc covering all corners, in [0, 360); ``WEST > EAST`` means
885 the arc crosses the 0/360 seam; ``SPAN = 360`` means all longitudes.
886 * ``FP_CONTAINS_POLE`` — winding number of the corner ring ≠ 0.
887 * ``FP_PS_*`` — polar stereographic bbox (m, R = 1737.4 km) for images
888 whose centre is poleward of ±60° (``FP_PS_HEMI`` = ``'S'`` / ``'N'``).
889 """
890 names = list(names)
891 if "CENTER_LATITUDE" not in names or "CENTER_LONGITUDE" not in names: 891 ↛ 892line 891 didn't jump to line 892 because the condition on line 891 was never true
892 return ""
893 clat, clon = "CENTER_LATITUDE", "CENTER_LONGITUDE"
894 pairs = _corner_pairs(names)
895 if pairs: 895 ↛ 899line 895 didn't jump to line 899 because the condition on line 895 was always true
896 lats = [f"coalesce({a}, {clat})" for a, _ in pairs]
897 lons = [f"coalesce({b}, {clon})" for _, b in pairs]
898 else:
899 lats, lons = [clat], [clon]
901 l360 = "[" + ", ".join(f"_lon360({x})" for x in lons) + "]"
902 l180 = "[" + ", ".join(f"(_lon360({x}) + 180.0) % 360.0 - 180.0" for x in lons) + "]"
903 span0 = f"(list_max({l360}) - list_min({l360}))"
904 span180 = f"(list_max({l180}) - list_min({l180}))"
905 if pairs: 905 ↛ 910line 905 didn't jump to line 910 because the condition on line 905 was always true
906 ring = lons + [lons[0]]
907 winding = " + ".join(f"_dlon({ring[i]}, {ring[i + 1]})" for i in range(len(lons)))
908 pole = f"(abs({winding}) > 180.0)"
909 else:
910 pole = "false"
911 west = (
912 f"CASE WHEN {pole} THEN 0.0 WHEN {span0} <= {span180} THEN list_min({l360}) "
913 f"ELSE _lon360(list_min({l180})) END"
914 )
915 span = f"CASE WHEN {pole} THEN 360.0 WHEN {span0} <= {span180} THEN {span0} ELSE {span180} END"
916 lat_list = "[" + ", ".join(lats + [clat]) + "]"
918 hemi = (
919 f"CASE WHEN {clat} < -{POLAR_LAT_THRESHOLD} THEN 'S' WHEN {clat} > {POLAR_LAT_THRESHOLD} THEN 'N' END"
920 )
921 xs = "[" + ", ".join(f"_ps_x({a}, {b}, FP_PS_HEMI)" for a, b in zip(lats, lons)) + "]"
922 ys = "[" + ", ".join(f"_ps_y({a}, {b}, FP_PS_HEMI)" for a, b in zip(lats, lons)) + "]"
924 # Two-level: the inner block defines pole/hemi aliases, reused below via
925 # DuckDB's lateral column aliases.
926 exprs = [
927 f"{pole} AS FP_CONTAINS_POLE",
928 f"CASE WHEN FP_CONTAINS_POLE AND {clat} < 0 THEN -90.0 ELSE list_min({lat_list}) END AS FP_MIN_LAT",
929 f"CASE WHEN FP_CONTAINS_POLE AND {clat} >= 0 THEN 90.0 ELSE list_max({lat_list}) END AS FP_MAX_LAT",
930 f"{west} AS FP_LON_WEST",
931 f"{span} AS FP_LON_SPAN",
932 "CASE WHEN FP_LON_SPAN >= 360 THEN 360.0 ELSE _lon360(FP_LON_WEST + FP_LON_SPAN) END AS FP_LON_EAST",
933 f"{hemi} AS FP_PS_HEMI",
934 f"CASE WHEN FP_PS_HEMI IS NOT NULL THEN list_min({xs}) END AS FP_PS_XMIN",
935 f"CASE WHEN FP_PS_HEMI IS NOT NULL THEN list_max({xs}) END AS FP_PS_XMAX",
936 f"CASE WHEN FP_PS_HEMI IS NOT NULL THEN list_min({ys}) END AS FP_PS_YMIN",
937 f"CASE WHEN FP_PS_HEMI IS NOT NULL THEN list_max({ys}) END AS FP_PS_YMAX",
938 ]
939 if wkt:
940 exprs.append(f"{footprint_wkt_sql(names, west_expr='FP_LON_WEST')} AS FP_WKT")
941 return ",\n ".join(exprs)
944def footprint_wkt_sql(names: Iterable[str], west_expr: str = "FP_LON_WEST") -> str:
945 """SQL for a ``POLYGON((lon lat, ...))`` WKT of the corner ring.
947 Longitudes are unwrapped relative to the footprint's western edge so the
948 ring never jumps at the seam (values may exceed 360). For pole-enclosing
949 footprints the lat/lon ring is only indicative — use the ``FP_PS_*`` box.
950 """
951 pairs = _corner_pairs(names)
952 if not pairs: 952 ↛ 953line 952 didn't jump to line 953 because the condition on line 952 was never true
953 return "NULL"
954 pts = []
955 for a, b in pairs + [pairs[0]]:
956 lon = f"round({west_expr} + _lon360(round({b} - {west_expr}, 6)), 4)"
957 pts.append(f"CAST({lon} AS VARCHAR) || ' ' || CAST(round({a}, 4) AS VARCHAR)")
958 return "'POLYGON((' || " + " || ', ' || ".join(pts) + " || '))'"
961# Search regions
962def region_sql(region: dict) -> str:
963 """WHERE clause matching ``FP_*`` columns against a
964 [`query_region`][lunar_toolkit.geometry.query_region]."""
965 s, n = region["lat"]
966 conds = [f"FP_MAX_LAT >= {s!r}", f"FP_MIN_LAT <= {n!r}"]
967 if region["span"] < 360.0:
968 conds.append(f"_arc_overlap(FP_LON_WEST, FP_LON_SPAN, {region['west']!r}, {region['span']!r})")
969 for hemi, (x0, x1, y0, y1) in region["ps"].items():
970 conds.append(
971 f"(FP_PS_HEMI IS DISTINCT FROM '{hemi}' OR (FP_PS_XMAX >= {x0!r} AND FP_PS_XMIN <= {x1!r} "
972 f"AND FP_PS_YMAX >= {y0!r} AND FP_PS_YMIN <= {y1!r}))"
973 )
974 return "(" + " AND ".join(conds) + ")"
977# Volume discovery / housekeeping
978_VOLUME_RE = re.compile(r"^LROLRC_(\d{4})([A-Z]?)$")
981def _norm_volume(v: str) -> str:
982 return str(v).strip().strip("/").upper()
985def _volume_key(v: str) -> tuple[int, str]:
986 m = _VOLUME_RE.match(_norm_volume(v))
987 if not m: 987 ↛ 988line 987 didn't jump to line 988 because the condition on line 987 was never true
988 return (-1, _norm_volume(v))
989 return (int(m.group(1)), m.group(2))
992def parse_s3_listing(xml_text: str) -> tuple[list[str], str | None]:
993 """Return ``(common_prefixes, next_continuation_token)`` from ListObjectsV2 XML."""
994 root = ET.fromstring(xml_text)
995 ns = ""
996 if root.tag.startswith("{"): 996 ↛ 998line 996 didn't jump to line 998 because the condition on line 996 was always true
997 ns = root.tag.split("}")[0] + "}"
998 prefixes = [el.text for el in root.iter(f"{ns}Prefix") if el.text and el is not root.find(f"{ns}Prefix")]
999 token = None
1000 if (root.findtext(f"{ns}IsTruncated") or "").lower() == "true": 1000 ↛ 1001line 1000 didn't jump to line 1001 because the condition on line 1000 was never true
1001 token = root.findtext(f"{ns}NextContinuationToken")
1002 return prefixes, token
1005def list_store_volumes(base_url: str, volume_re: str, timeout: float = 30.0) -> list[str]:
1006 """List volume directories under a dataset on the PDS S3 store."""
1007 if not base_url.startswith(PDS_STORE_URL):
1008 raise ValueError(f"Not a PDS store URL: {base_url}")
1009 prefix = base_url[len(PDS_STORE_URL) :]
1010 pat = re.compile(rf"^{re.escape(prefix)}({volume_re})/$")
1011 vols: list[str] = []
1012 token: str | None = None
1013 while True:
1014 q = {"list-type": "2", "prefix": prefix, "delimiter": "/"}
1015 if token:
1016 q["continuation-token"] = token
1017 url = PDS_STORE_URL + "?" + urllib.parse.urlencode(q)
1018 req = urllib.request.Request(url, headers={"User-Agent": _USER_AGENT})
1019 with urllib.request.urlopen(req, timeout=timeout) as resp:
1020 prefixes, token = parse_s3_listing(resp.read().decode("utf-8", "replace"))
1021 vols += [m.group(1) for p in prefixes if (m := pat.match(p))]
1022 if not token:
1023 break
1024 return sorted(set(vols), key=_volume_key)
1027def probe_volumes(
1028 base_url: str,
1029 start: str,
1030 max_numbers: int = 12,
1031 timeout: float = 15.0,
1032) -> list[str]:
1033 """Probe ``<vol>/INDEX/INDEX.LBL`` for volumes at/after *start*.
1035 Volume names look like ``LROLRC_0066C`` (number + optional A/B/C… part
1036 letter). Stops after two consecutive numbers without any hit.
1037 """
1038 num, _ = _volume_key(start)
1039 m = re.match(r"^(LROLRC_)(\d{4})", _norm_volume(start))
1040 if not m or num < 0:
1041 return [_norm_volume(start)]
1042 prefix = m.group(1)
1043 found = [_norm_volume(start)]
1044 misses = 0
1045 for n in range(num, num + max_numbers):
1046 hit = False
1047 for suffix in ("", "A", "B", "C", "D", "E"):
1048 vol = f"{prefix}{n:04d}{suffix}"
1049 if _volume_key(vol) <= _volume_key(start):
1050 continue
1051 if _url_exists(base_url + vol + "/" + INDEX_DIR + VOLINDEX_LBL_FILENAME, timeout):
1052 found.append(vol)
1053 hit = True
1054 if n > num:
1055 misses = 0 if hit else misses + 1
1056 if misses >= 2:
1057 break
1058 return sorted(set(found), key=_volume_key)
1061def _url_exists(url: str, timeout: float = 15.0) -> bool:
1062 req = urllib.request.Request(url, method="HEAD", headers={"User-Agent": _USER_AGENT})
1063 try:
1064 with urllib.request.urlopen(req, timeout=timeout) as resp:
1065 return 200 <= getattr(resp, "status", 200) < 300
1066 except urllib.error.HTTPError:
1067 return False
1070def _fetch_head(url: str, dest: Path, nbytes: int, timeout: float = 60.0) -> Path:
1071 """Fetch the first *nbytes* of *url* (HTTP Range) into *dest*."""
1072 dest.parent.mkdir(parents=True, exist_ok=True)
1073 req = urllib.request.Request(url, headers={"User-Agent": _USER_AGENT, "Range": f"bytes=0-{nbytes - 1}"})
1074 with urllib.request.urlopen(req, timeout=timeout) as resp:
1075 data = resp.read(nbytes)
1076 tmp = dest.with_suffix(dest.suffix + ".part")
1077 tmp.write_bytes(data)
1078 tmp.replace(dest)
1079 return dest
1082def clean_stale_build_files(directory: Path) -> list[Path]:
1083 """Delete leftovers of interrupted builds (``*.build.duckdb[.wal]``, ``*.tmp``)."""
1084 removed: list[Path] = []
1085 directory = Path(directory)
1086 if not directory.is_dir(): 1086 ↛ 1087line 1086 didn't jump to line 1087 because the condition on line 1086 was never true
1087 return removed
1088 for pattern in ("*.build.duckdb", "*.build.duckdb.wal", "*.parquet.tmp", "*.update.parquet"):
1089 for f in directory.glob(pattern):
1090 try:
1091 f.unlink()
1092 removed.append(f)
1093 log.info("Removed stale build file %s", f)
1094 except OSError:
1095 pass
1096 return removed
1099# LBL parsing
1100def parse_lbl_columns(lbl_content: str) -> list[dict]:
1101 """Parse ``OBJECT = COLUMN`` blocks of a PDS3 label (pure Python).
1103 Keys are lower-cased LBL keywords; ``column_number``, ``start_byte`` and
1104 ``bytes`` are ints. Sorted by ``column_number``.
1105 """
1106 blocks = re.findall(
1107 r"OBJECT\s*=\s*COLUMN(.*?)END_OBJECT\s*=\s*COLUMN",
1108 lbl_content,
1109 re.DOTALL,
1110 )
1112 def _parse_value(val: str) -> str:
1113 val = val.strip().strip('"')
1114 return re.sub(r"\s*\n\s*", " ", val).strip()
1116 pattern = re.compile(r"([A-Z_]+)\s*=\s*(.*?)(?=\n\s*[A-Z_]+\s*=|\Z)", re.DOTALL)
1117 records = []
1118 for block in blocks:
1119 rec: dict[str, Any] = {
1120 m.group(1).strip().lower(): _parse_value(m.group(2)) for m in pattern.finditer(block)
1121 }
1122 for k in ("column_number", "start_byte", "bytes"):
1123 if k in rec: 1123 ↛ 1122line 1123 didn't jump to line 1122 because the condition on line 1123 was always true
1124 try:
1125 rec[k] = int(rec[k])
1126 except ValueError:
1127 rec[k] = None
1128 if rec.get("name") and rec.get("start_byte") and rec.get("bytes"): 1128 ↛ 1118line 1128 didn't jump to line 1118 because the condition on line 1128 was always true
1129 records.append(rec)
1130 return sorted(records, key=lambda r: r.get("column_number") or 0)
1133def _parse_pds_lbl_columns(lbl_content: str) -> pd.DataFrame:
1134 """``parse_lbl_columns`` as a DataFrame (the ``CumIndex.schema`` format)."""
1135 import pandas as pd
1137 df = pd.DataFrame(parse_lbl_columns(lbl_content))
1138 for col in ("column_number", "start_byte", "bytes"):
1139 if col in df.columns: 1139 ↛ 1138line 1139 didn't jump to line 1138 because the condition on line 1139 was always true
1140 df[col] = pd.to_numeric(df[col], errors="coerce").astype("Int64")
1141 preferred = ["name", "column_number", "data_type", "start_byte", "bytes", "format", "description"]
1142 ordered = [c for c in preferred if c in df.columns]
1143 extras = [c for c in df.columns if c not in preferred]
1144 return df[ordered + extras].reset_index(drop=True)
1147def _parse_lbl_meta(lbl_content: str) -> dict:
1148 """Top-level label fields (ints for sizes, strings for ids / checksum)."""
1149 meta: dict = {}
1150 for key in ("RECORD_BYTES", "FILE_RECORDS", "ROWS", "COLUMNS", "ROW_BYTES"):
1151 m = re.search(rf"^\s*{key}\s*=\s*(\d+)", lbl_content, re.MULTILINE)
1152 if m: 1152 ↛ 1150line 1152 didn't jump to line 1150 because the condition on line 1152 was always true
1153 meta[key.lower()] = int(m.group(1))
1154 for key in ("MD5_CHECKSUM", "VOLUME_ID", "DATA_SET_ID", "PRODUCT_CREATION_TIME"):
1155 m = re.search(rf'^\s*{key}\s*=\s*"?([^"\r\n]+?)"?\s*$', lbl_content, re.MULTILINE)
1156 if m:
1157 meta[key.lower()] = m.group(1).strip()
1158 return meta
1161# Private helpers
1162_NAC_ONLY_SQL = "(" + " OR ".join(f"PRODUCT_ID LIKE '%{s}'" for s in NAC_SUFFIXES) + ")"
1165def _angle_conditions(
1166 emission_angle: tuple[float, float] | None,
1167 incidence_angle: tuple[float, float] | None,
1168 phase_angle: tuple[float, float] | None,
1169) -> list[str]:
1170 out = []
1171 for col, rng in (
1172 ("EMISSION_ANGLE", emission_angle),
1173 ("INCIDENCE_ANGLE", incidence_angle),
1174 ("PHASE_ANGLE", phase_angle),
1175 ):
1176 if rng is not None: 1176 ↛ 1177line 1176 didn't jump to line 1177 because the condition on line 1176 was never true
1177 out.append(f"{col} BETWEEN {float(rng[0])} AND {float(rng[1])}")
1178 return out
1181def _col_select(columns: list[str] | None, extra: list[str] | None = None) -> str:
1182 """Build a SELECT column list, always including the key identity columns."""
1183 if not columns:
1184 return "*"
1185 always = {"PRODUCT_ID", "FILE_SPECIFICATION_NAME", "CENTER_LATITUDE", "CENTER_LONGITUDE", "START_TIME"}
1186 return ", ".join(_ident(c) for c in sorted(always | set(columns) | set(extra or [])))
1189def _ident(name: str) -> str:
1190 return '"' + str(name).replace('"', '""') + '"'
1193def _sql_str(value) -> str:
1194 return "'" + str(value).replace("'", "''") + "'"
1197def _utcnow() -> str:
1198 return _dt.datetime.now(_dt.timezone.utc).replace(microsecond=0).isoformat()