Repository navigation
v1.5 — 20m resolution analysis with GPU acceleration (81.2% Wueller agreement) - #1
Merged
Merged
Conversation
Adds a `use_gpu` flag that swaps numpy + scipy.ndimage.map_coordinates
for cupy + cupyx.scipy.ndimage.map_coordinates inside the LOLA ray-march.
Numerically equivalent — max |diff| 2.3e-4 deg over all 231 M finite cells
of the 240 m grid, NaN masks identical, GPU output bit-identical
run-to-run. On the NVIDIA Spark (GB10 Blackwell, cc 12.1, CUDA 13):
250.9 s -> 4.0 s (62x speedup) at 240 m.
Workaround in-line: cupyx.scipy.ndimage.map_coordinates JIT-fails when
cval is `float('nan')` / `np.float32('nan')`; passing `np.nan` (the
module-level singleton) compiles cleanly. Documented at the call site.
Adds `[project.optional-dependencies] gpu` (cupy-cuda13x on aarch64,
cupy-cuda12x on x86_64) and `tests/test_los_to_earth_gpu.py` which
auto-skips off-GPU and asserts GPU/CPU equivalence + determinism on a
synthetic 64x64 grid.
Foundation for the v1.5 per-region 20 m run: each padded tile peaks at
~25 GB horizon array (100 km buffer, ~11k x 11k cells) and is now
GPU-fast.
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
Add `selene preprocess --tiled-per-region` v1.5 mode: skip the global 240 m preprocess loop and instead run the Earth-LOS horizon profile derivation per USGS region tile (polygon bbox + 100 km buffer) on the GPU. Tiles are cached as data/processed/horizon_profile_southpole_<resolution>m_<region_code>.npz with `region_code` lower-cased (sp, mm, mp, hw, n1, n2, g2, cb, ma). Add `selene download lola --resolution 20m` for the v1.5 source product (~1.85 GB IMG/LBL pair, same PDS Geosciences directory as the v1.4 80 m product). The 240 m default codepath is unchanged; tiled mode is opt-in via flag. Validated end-to-end on Slater Plain at the 80 m source: 7.5 M pixels, 31 s wall time on Blackwell GB10, finite output everywhere with the initial -90° sentinel surviving only in the tile buffer ring (0.07 % of pixels, 0 % inside the polygon interior). Adds tests/test_preprocess_tiled.py covering compute_tile_specs, region filtering, end-to-end synthetic roundtrip, and idempotency.
Add `selene rank-per-region --tiled-per-region [--resolution 20]`: for each USGS polygon, derive 20 m slope from a windowed LOLA DEM, derive 20 m Earth-LOS visibility from the v1.5 horizon-profile NPZ produced by `selene preprocess --tiled-per-region`, resample the 240 m illumination and aggregate-score COGs onto the tile grid, apply NASA's four HLS hard filters at 20 m, run greedy NMS within the polygon, write sites + summary to `data/outputs/per_region_tiled/`. The headline change vs v1.4.2 240 m: the HLS 100 m buffer-from-steep constraint becomes meaningful. At 240 m the buffer is 0.42 px (looser than the published 100 m); at 20 m it is 5 px and the constraint applies as written. Smoke-tested on Slater Plain at 80 m: eligibility tightens from 2.66 % (v1.4.2 240 m) to 1.68 % at 80 m, with the same 9 sites under the 2 km NMS cap. Tighter tightening expected at 20 m once the source is downloaded. Adds tests/test_rank_per_region_tiled.py with end-to-end synthetic flow (skipping the LOS filter so the smooth-bowl test surface can pass) and a missing-NPZ regression guard.
Add a use_gpu CuPy dispatch to the libration-cycle visibility sweep — the dominant CPU cost in the v1.5 tiled ranker at 20 m. Output is bitwise-equivalent to the NumPy path on real data (visibility fractions are exact rationals count/24, so no float drift). Measured 28x speedup on a 7.5 M-pixel real-data tile (Slater Plain, 80 m): 15.1 s -> 0.5 s. The tiled ranker's per-region wall time on the same input drops from 21 s to 9.4 s once LOS visibility no longer dominates. Wires use_gpu=True in rank_per_region_tiled.process_tile with a RuntimeError fallback to the NumPy path on hosts without CuPy. Adds tests/test_los_to_earth_gpu.py::test_gpu_matches_cpu_visibility_fraction asserting bitwise equivalence on a small synthetic horizon profile; auto-skips on CPU-only hosts.
End-to-end v1.5 rank-per-region runs were OOM-killed mid-pipeline on the GB10 host because each tile's ~17 GB float32 horizon profile and ~3 GB lat/lon/gamma stayed pinned through the rest of the tile's HLS work, and the cupy memory pool was never drained on the early-return "no HLS-compliant cells" path. Three changes: - Drop horizon_tile and pixel_lat/pixel_lon/gamma right after the Earth-LOS sweep — they are not used past that call, and on GB10 unified memory they otherwise crowd both numpy and cupy out of DRAM before the polygon mask + EDT + score-resample work runs. - Mirror the late-return cleanup on the no-compliant-cells return so an empty-polygon tile (e.g. Malapert Massif) does not leak its cupy + numpy state into the next iteration. - Wrap the per-tile call site in run() in try/except/finally: the finally drains the cupy memory pool unconditionally, and the except catches per-tile failures (e.g. corrupt NPZ cache) so a single bad tile does not throw away the work already done on the other 8 polygons.
… plots Adds the v1.5 result section under "Quantitative comparison against Wueller et al. 2026" plus a roadmap entry. The 20 m per-region tiled run produces 69 sites across 8 of 9 USGS regions; 56 of 69 (81 %) match within 5 km of an in-scope Wueller 2026 site, and the median matched-pair distance tightens from 1.88 km at 240 m to 1.76 km (-6.5 %). Selene-to-Wueller match rate is essentially flat across resolutions (80 % -> 81 %); Wueller-to-selene rate slips 2.7 pp (63 % -> 60 %) because the finer slope/buffer re-ranks sites within MMP, SP, and CB to a different cell within the same polygon, not because the methodology disagrees. The structural outliers from v1.4.2 — Malapert Massif (zero compliant cells) and de Gerlache Rim 2 (only 2 cells, 0/2 within 5 km of dGR2) — persist at fine resolution, confirming they are genuine terrain, not 240 m sampling artefacts. - README: new "Resolution analysis (v1.5 — 20 m Wueller-class)" subsection with the comparison table, per-region detail table, the two new images, the engineering footprint of v1.5, and a roadmap bullet capturing the headline. - docs/img/selene_vs_wueller_20m.png — overlay map regenerated from the 20 m per-region catalog. - docs/img/resolution_sensitivity_v15.png — per-region side-by-side of v1.4.2 (240 m) vs v1.5 (20 m) match percentage and median matched-pair distance. - notebooks/10_resolution_v15.py — script that produces both plots.
…e, Barker reference) Subsumes the v1.4.2 numbers in the opening into v1.5's 81.2% / 1.76 km headline. Promotes selene_vs_wueller_20m.png to the headline image. Adds Barker et al. 2021 reference for the 20 m LOLA south polar product.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Re-runs the per-region pipeline at 20 m (Wueller-class resolution), with GPU-accelerated horizon profile derivation enabling per-region tiled analysis on commodity unified-memory hardware. Methodology tightens at peer-reviewed resolution: 56 of 69 selene sites (81.2%) match within 5 km of an in-scope Wueller 2026 site, median 1.76 km (was 80.0% / 1.88 km at 240 m).
Motivation
v1.4.2 shipped agreement against Wueller's catalog at 240 m resolution. Wueller themselves used 20 m. This PR closes that resolution gap so the comparison is direct, and ships GPU acceleration as the architecture that makes 20 m feasible. The GPU port is also infrastructure for v1.6 (5 m exploratory analysis) without further work.
Checklist
ruff check .passesruff format --check .passespytestpasses locally (350 collected, 346 passed, 4 skipped: LEND only)test_los_to_earth_gpu.py,test_preprocess_tiled.py,test_rank_per_region_tiled.py,test_per_region_validation.py)Scope notes
This PR explicitly does not address:
The GPU memory-cleanup scaffolding in
rank_per_region_tiled.py(try/finally + explicit cupy mempool drain + per-tile exception isolation) is intentionally retained as durable hardening rather than removed once stable; v1.6's 5 m resolution will scale per-tile allocations ~16× and require this defensive structure.Headline result
Six of eight in-scope regions hold or improve their median matched-pair distance at 20 m. Structural outliers (Malapert Massif: 0 compliant cells; de Gerlache Rim 2: 2 cells, 0/2 matched) persist at fine resolution — confirming they are genuine terrain-driven, not 240 m sampling artifacts.
Performance
GPU acceleration via CuPy on Blackwell (DGX Spark GB10):
Citation