lunar_toolkit.sfs

Shape-from-shading DEMs from LRO NAC images with ISIS and ASP (Bertone et al. 2026 workflow).

SfsConfig(resolution=5.0, projwin=None, preset='bertone', ba_resolution=1.0, reduce='auto', native_resolution=1.0, csm=True, web_spice=True, registration='stereo', stereo_pairs=30, stereo_convergence=(5.0, 60.0), max_displacement=400.0, refine=True, max_offset_p85=1.8, max_large_residual_fraction=0.85, large_residual_px=1.0, min_matches=1, min_images=3, selection='coverage', max_shadow_fraction=0.95, selection_views=3, selection_azimuth_sep=30.0, max_sfs_images=None, selection_clusters=30, sfs_cameras='csm', cameras_from='ba', published_region=None, max_sfs_slew=5.0, selection_grid=250, region_id=None, ensemble_members=0, ensemble_size=None, ensemble_crop=2500.0, seed=0, processes=None, threads=None, bundle_adjust_options=dict(), refine_options=dict(), sfs_options=dict()) dataclass

Settings for SfsPipeline.

Parameters:
  • resolution (float, default: 5.0 ) –

    SfS DEM grid spacing (m/px).

  • projwin (tuple, default: None ) –

    (ulx, uly, lrx, lry) crop of the prior DEM, in its projected units.

  • preset (('bertone', 'asp'), default: "bertone" ) –

    Option defaults: Bertone et al. (2026) or the ASP SfS recipe (PRESETS).

  • ba_resolution (float, default: 1.0 ) –

    Grid spacing of the images mapprojected for bundle-adjustment matching and stereo registration (default 1 m, NAC native; None = resolution).

  • reduce ((int, 'auto' or None), default: 'auto' ) –

    ISIS reduce factor bringing the cubes to resolution before SfS. "auto": round(resolution / native) per image (native = observation pixel_scale, else native_resolution). Reduced cubes use ISIS cameras (with --use-approx-camera-models); bundle adjustment still uses CSM.

  • native_resolution (float, default: 1.0 ) –

    Assumed image m/px when the observation table does not give it.

  • csm (bool, default: True ) –

    CSM cameras (7–15× faster than ISIS) for mapproject, bundle_adjust and stereo.

  • web_spice (bool, default: True ) –

    spiceinit web=true (USGS SPICE web service) instead of local kernels.

  • registration (('stereo', 'residuals', 'none'), default: "stereo" ) –

    Registering the adjusted cameras to the prior DEM (translation only). "stereo": DEMs from up to stereo_pairs pairs, mosaicked and pc_align-ed (falls back to "residuals" without pairs); "residuals": the bundle-adjustment triangulated points.

  • stereo_pairs (int, default: 30 ) –

    Maximum stereo pairs for "stereo" registration.

  • stereo_convergence ((min, max), default: (5.0, 60.0) ) –

    Convergence-angle limits (deg) for those pairs.

  • max_displacement (float, default: 400.0 ) –

    pc_align --max-displacement (m).

  • refine (bool, default: True ) –

    Re-run bundle adjustment on the kept, registered cameras with --heights-from-dem (reusing the matches).

  • max_offset_p85 (float, default: 1.8 ) –

    Drop images whose 85th-percentile mapprojected match offset (m) is larger. Applied to the DEM-registered, refined cameras (ASP's ba_align_ref report); with refine off, to the unregistered ba/ cameras.

  • max_large_residual_fraction (float, default: 0.85 ) –

    Drop images with more than this fraction of reprojection errors over large_residual_px pixels.

  • min_matches (int, default: 1 ) –

    Drop images with fewer inlier matches.

  • min_images (int, default: 3 ) –

    Stop with PipelineError when a quality filter leaves fewer images.

  • selection (('coverage', 'image_subset', 'all'), default: "coverage" ) –

    SfS image selection: greedy per-pixel illumination coverage (select_by_illumination), ASP image_subset per Sun-azimuth bin, or all kept images.

  • max_shadow_fraction (float, default: 0.95 ) –

    Images at least this shadowed (below the SfS shadow threshold) are dropped.

  • selection_views ((int, float), default: 3 ) –

    Illumination directions wanted per pixel and their minimum separation (deg).

  • selection_azimuth_sep ((int, float), default: 3 ) –

    Illumination directions wanted per pixel and their minimum separation (deg).

  • max_sfs_images (int, default: None ) –

    Cap on the selected images ("coverage" / "kmedoids" selection).

  • selection_clusters (int, default: 30 ) –

    Medoids seeding "kmedoids" selection (Bertone et al. 2026, §2: ~30).

  • sfs_cameras (('csm', 'isis'), default: "csm" ) –

    Camera models for SfS on reduced cubes: CSM cameras re-mapped to the reduced image (lunar_toolkit.csm.isis_reduce_with_csm, as in Bertone et al. 2026, App. A) or ISIS cameras with --use-approx-camera-models.

  • cameras_from (('ba', 'published'), default: "ba" ) –

    "published" uses the bundle-adjusted cameras of Bertone et al. (2026) for published_region (lunar_toolkit.reference.published_cameras) and skips bundle adjustment, registration and refinement; images without one are removed.

  • published_region (str or int, default: None ) –

    Region of the published cameras (default: from region_id A3NN).

  • max_sfs_slew (float, default: 5.0 ) –

    Only images with |slew| < this (deg) are selected for SfS (Bertone et al. 2026, §2: nadir images); slewed images still join bundle adjustment and stereo registration. Unknown slews count as nadir; None disables the cut.

  • selection_grid (int, default: 250 ) –

    Cells along the longer side of the selection grid.

  • region_id (str, default: None ) –

    ROIID for product file names (default: the work directory name).

  • ensemble_members (int, default: 0 ) –

    SfS runs on random subsets of the kept images (sharing at most 40% of their images) over an ensemble_crop m square at the DEM centre, for the ensemble uncertainty maps (Bertone et al. 2026, §3.2); 0 skips this.

  • ensemble_size (int, default: None ) –

    Images per ensemble member (default: as many as selected).

  • ensemble_crop (float, default: 2500.0 ) –

    Side of the ensemble crop (m).

  • seed (int, default: 0 ) –

    Random seed of the ensemble image subsets.

  • processes (int, default: None ) –

    Passed to parallel_bundle_adjust, parallel_stereo, parallel_sfs and mapproject.

  • threads (int, default: None ) –

    Passed to parallel_bundle_adjust, parallel_stereo, parallel_sfs and mapproject.

  • bundle_adjust_options (dict, default: dict() ) –

    Overrides merged onto the preset; a value of None drops a default.

  • refine_options (dict, default: dict() ) –

    Overrides merged onto the preset; a value of None drops a default.

  • sfs_options (dict, default: dict() ) –

    Overrides merged onto the preset; a value of None drops a default.

options(step)

Preset options for "bundle_adjust", "refine" or "sfs" with overrides applied.

PipelineError

Bases: RuntimeError

The SfS pipeline cannot continue (e.g. too few images left).

SfsImage(product_id, cube, camera=None, sun_azimuth=None, sun_elevation=None, pixel_scale=None, mapprojected=None, shadow_fraction=None, slew_angle=None) dataclass

One input image: cube, optional CSM camera, Sun angles and derived files.

RemovedImage(image, stage, reason) dataclass

An image dropped by a quality filter: stage (ba, refine, selection) and why.

SfsPipeline(work_dir, prior_dem, config=None, force=False, *, observations=None, aoi=None, data_dir=None, run_config=None, provenance=True)

Build an SfS DEM over a region from NAC images and a prior DEM.

run follows Bertone et al. (2026, PSJ 7:118) and the ASP large-area SfS recipe:

  1. calibrate the images and write CSM cameras (images/);
  2. crop the prior DEM and resample it to the target grid (dem/);
  3. sort by Sun azimuth (sfs --query);
  4. mapproject at native resolution (prj/) and bundle-adjust with matches between the mapprojected images (ba/);
  5. drop images with no matches or mostly large residuals (selection/removed.csv);
  6. register the cameras to the prior DEM, translation only, from stereo DEMs (stereo/) or the BA points; apply to all cameras (align/, ba_align/);
  7. refine with --heights-from-dem and filter again, now also on mapprojected match offsets (ba_refine/); fewer than min_images left stops the run;
  8. mapproject at the target resolution (prj_sfs/), write max-lit QA mosaics, all images and per 30° Sun-azimuth bin (qa/);
  9. select an illumination-diverse subset (selection/);
  10. reduce the selected cubes to the target resolution and run parallel_sfs (sfs/);
  11. products (products/ via lunar_toolkit.sfs_products, else hillshade
    • diff).

Steps whose outputs exist are skipped, so a rerun resumes.

Parameters:
  • work_dir (Path) –

    Directory for all intermediate and final files.

  • prior_dem (Path) –

    Reference terrain, e.g. a LOLA polar stereographic LDEM GeoTIFF.

  • config (SfsConfig, default: None ) –
  • force (bool, default: False ) –

    Recompute steps even when their outputs exist.

  • observations (DataFrame, default: None ) –

    Observation table (lunar_toolkit.observation schema) with the input images, for native pixel scales and stereo pairs; looked up in the local EDR CumIndex when needed and not given.

  • aoi (optional, default: None ) –

    Region the stereo pairs must overlap in: anything lunar_toolkit.stereo.find_stereo_pairs takes (dict, GeoJSON / WKT).

  • data_dir (Path, default: None ) –

    CumIndex data directory.

  • run_config (RunConfig, default: None ) –

    The full run settings (lunar_toolkit.config) written to config.yaml and provenance.json; default: work_dir, prior_dem, the sources and config.

  • provenance (bool, default: True ) –

    Write config.yaml and append this run to provenance.json (never in a dry run).
Attributes:
  • all_images, images (list of SfsImage) –

    Every input / those passing the quality filters (Sun-azimuth order).

  • removed (list of RemovedImage) –
  • selected, unselected (list of SfsImage) –

    The SfS inputs and the remaining good images (e.g. for an uncertainty ensemble); both carry mapprojected (prj_sfs/) and sun_azimuth.

  • dem (Path) –

    The cropped prior DEM on the target grid.
Example

::

cfg = SfsConfig(resolution=5, projwin=(-10000, 50000, 0, 40000))
p = SfsPipeline("work/malapert", "LDEM_87S_5MPP_ADJ.TIF", cfg)
dem = p.run(["raw/M1234567890LE.IMG", ...])

sfs_camera_prefix property

Adjustments renamed for the reduced ISIS cubes.

camera_prefix property

Adjustments of the final (registered, refined) cameras; None for published cameras.

run(sources)

All steps; returns the final SfS DEM.

Writes config.yaml and appends a run to provenance.json first (see start_provenance), then records each step there.

effective_config(sources=())

run_config with this pipeline's work dir, prior DEM and SfsConfig.

start_provenance(sources)

Write config.yaml and start a provenance.json run; returns the config path.

Inputs are sources and the prior DEM. Nothing is written in a dry run or with provenance=False.

prepare_images(sources)

Calibrate EDRs (.IMG) or take ready cubes; write CSM cameras.

prepare_dem()

Crop the prior DEM to projwin and resample to the target grid.

sort_by_sun()

Fill Sun azimuth / elevation and sort the images by azimuth.

Bundle adjustment matches neighbours in this order, and parallel_sfs expects it.

mapproject(images, subdir, resolution, **options)

Mapproject images onto the prior DEM (cropped to projwin) into <subdir>/.

bundle_adjust()

parallel_bundle_adjust with matches from the mapprojected images.

--mapproj-dem writes the per-image match-offset report used by filter_images.

workspace(read_only=False)

The work dir's DuckDB workspace (lunar.duckdb, lunar_toolkit.workspace).

filter_images(prefix, stage)

Drop images that ASP's bundle-adjustment reports flag (see flag_images).

Loads <prefix>-* reports (residual stats and raw pixels, match offsets, point maps, convergence angles) into the workspace as run <prefix>'s directory name (ba, ba_refine) and removes the images of its image_rejections view, with this config's thresholds. Removed images are logged and appended to selection/removed.csv.

At ba the match-offset rule waits for the refinement when that runs: before registration the network floats tens of metres off the DEM, so convergent pairs mapproject apart by dh * tan(convergence) however well they are tied (ASP: "do not use this option for any initial evaluation of bundle adjustment").

register()

Register the adjusted cameras to the prior DEM (translation only); returns their prefix.

stereo_pairs()

Up to stereo_pairs kept-image pairs ranked by find_stereo_pairs.

refine()

Bundle-adjust the kept, registered cameras again with --heights-from-dem.

Reuses the ba/ matches (--match-files-prefix); --mapproj-dem gives the match-offset report for a second filter_images.

mapproject_for_sfs()

Mapproject the kept images at the target resolution with the final cameras (prj_sfs/).

qa_mosaics(bin_width=30.0)

Max-lit mosaics of all kept images and per Sun-azimuth bin (qa/).

Ghosting in these shows misregistered images (Bertone et al. 2026, App. A). Bins are named by their start, e.g. max_lit_az030.tif for 30–60°.

select_images()

Split the kept images into SfS inputs and the rest (selection/).

Images at least max_shadow_fraction shadowed are removed first. Writes selection.csv (every kept image: azimuth, shadow fraction, mapprojected path, selected) plus selected.txt / unselected.txt (one product ID per line).

reduce_factor(im)

ISIS reduce factor for im (1 = keep native resolution).

sfs_inputs(images)

Cubes, cameras and adjustment prefix for SfS on images.

Cubes needing it are reduce-d to the target resolution (Bertone et al. 2026, App. A: SfS's own resampling leaves artifacts). With sfs_cameras="csm" each CSM camera is re-mapped to its reduced cube (csm.isis_reduce_with_csm); otherwise ISIS cameras are used. The final .adjust files (resolution independent) are copied under sfs_camera_prefix with the new names; published cameras need none. Returns (cubes, cameras, prefix, isis_cameras).

use_published_cameras()

Replace the cameras with the published Bertone et al. (2026) adjusted ones.

Each published camera (a CSM state of the reduced full image) is re-mapped to the image of its cube (reference.published_camera_native) and written to images/<id>.published.json; images without one are removed.

shape_from_shading(images=None, prefix=None, dem=None, **options)

parallel_sfs on the target-grid prior DEM; returns <prefix>-DEM-final.tif.

images default to selected (else all kept images), prefix to sfs/run and dem (input DEM, e.g. a crop) to dem. options override the configured SfS options.

products(dem)

Standard products via lunar_toolkit.sfs_products.make_products when available.

Otherwise (module missing or its dependencies not installed) a hillshade of the SfS DEM and its difference from the prior.

ensemble()

SfS ensemble over a central crop and its spread (ensemble/, products/).

Members are ensemble_subsets of the kept images; returns ensemble_uncertainty's maps.

PolarStereographic(hemisphere, lon_0=0.0, radius=MOON_RADIUS_KM * 1000.0, x_0=0.0, y_0=0.0) dataclass

A polar stereographic projection (sphere, true scale at the pole) from a PROJ string.

from_proj4(proj4) classmethod

Parse e.g. +proj=stere +lat_0=-90 +lon_0=0 +k=1 +x_0=0 +y_0=0 +R=1737400 +units=m.

Raises:
  • ValueError –

    For anything else (other projections, oblique stereographic, scale ≠ 1, non-metre units).

flag_images(residual_stats, residual_norms, offsets, *, max_offset_p85=1.8, max_large_residual_fraction=0.85, large_residual_px=1.0, min_matches=1)

Images to drop after bundle adjustment (Bertone et al. 2026, §2), keyed by file name.

An image is flagged when it has fewer than min_matches residuals or match offsets, an 85th-percentile mapprojected match offset over max_offset_p85 m, or more than max_large_residual_fraction of its reprojection errors over large_residual_px. The value is the reason. The rules are the workspace's image_rejections view (lunar_toolkit.workspace), evaluated in an in-memory database.

select_by_illumination(lit, azimuth, *, views=3, azimuth_sep=30.0, max_images=None, min_gain=1)

Greedy per-pixel illumination-diverse image selection on a raster grid.

The raster analogue of lunar_toolkit.selection.select_for_coverage, after Bertone et al. (2026, §2): a lit cell is helped by an image lighting it if it has fewer than views illumination directions so far and the image's Sun azimuth is at least azimuth_sep from each of them. Each step takes the image helping most cells, until none helps min_gain cells or max_images are chosen.

Parameters:
  • lit ((n_images, n_cells) bool array) –
  • azimuth ((n_images,) array) –

    Sun azimuths (deg); NaN never helps.
Returns:
  • (chosen, count) –

    Image indices in pick order, and the directions reached per cell.

prefilter_candidates(df, max_slew=None, max_resolution=1.5)

Images finer than max_resolution m/px (and, with max_slew, |slew| < it).

Works on CumIndex rows (SLEW_ANGLE, RESOLUTION) or observation rows (slew_angle, pixel_scale). Rows with unknown values are dropped. The nadir cut of Bertone et al. (2026, §2) belongs to SfS selection (SfsConfig.max_sfs_slew), so searches keep slewed images by default: they give the stereo convergence the registration step needs.

cap_candidates(obs, aoi, max_images, *, max_sfs_slew=5.0, stereo_pairs=30, stereo_convergence=(5.0, 60.0), views=3, azimuth_sep=30.0)

At most max_images candidate images, chosen before anything is downloaded.

Keeps both images of the best stereo pairs (up to half the budget; they feed the registration step), then nadir images (|slew| < max_sfs_slew) chosen for illumination coverage over aoi (select_for_coverage: views directions at least azimuth_sep apart per cell), then the finest remaining nadir images. Returns the chosen rows in their original order.

search_images(aoi, *, max_slew=None, max_resolution=1.5, max_images=None, config=None, data_dir=None)

NAC EDRs over aoi from the local CumIndex, pre-filtered and capped; observation schema.

config supplies the SfS slew cut and stereo settings used by cap_candidates.

aoi_projwin(aoi, dem, snap=None)

(ulx, uly, lrx, lry) bounding aoi in the polar stereographic projection of dem.

aoi is anything lunar_toolkit.geometry.load_aoi takes (GeoJSON / WKT file, WKT string). The window is widened to multiples of snap (m) when given.