lunar_toolkit.sfs_products

Standard SfS DEM products and uncertainty estimates (Bertone et al. 2026, PSJ 7:118).

Grid(crs, transform, width, height) dataclass

A georeferenced raster grid (CRS, affine transform, size).

res property

Pixel size (map units); pixels are assumed square.

bounds property

(xmin, ymin, xmax, ymax).

ProductImage(path, sun_azimuth, resolution=None) dataclass

A mapprojected SfS input image, its Sun azimuth (deg from north) and native GSD (m/px).

IlluminationMaps(count, solar_bins, best_resolution, grid) dataclass

Per-pixel illumination statistics of the SfS inputs on the DEM grid.

Attributes:
  • count (ndarray(uint16)) –

    Number of lit images (value above the shadow threshold).

  • solar_bins (ndarray(uint8)) –

    Number of Sun-azimuth bins with at least one lit image.

  • best_resolution (ndarray(float32) or None) –

    Finest native GSD among the lit images (NaN where none); None when some image has no resolution.

  • grid (Grid) –

EnsembleStats(elevation_std, slope_std, count, grid) dataclass

Per-pixel spread of an ensemble of SfS DEMs on a common grid.

Attributes:
  • elevation_std, slope_std (ndarray) –

    Sample standard deviation (ddof=1) of elevation (m) and pixel-scale slope (deg); NaN where fewer than two members have data.

  • count (ndarray) –

    Members with a valid elevation at each pixel.

  • grid (Grid) –

PointComparison(lon, lat, diff) dataclass

DEM minus reference points, as written by geodiff (lon, lat in degrees; diff in m).

product_name(region_id, code, resolution, *suffix)

Table 1 file name, e.g. product_name("SR01", "TSLP", 5, 10) → SR01_GLDTSLP_5_10.tif.

read_raster(path, grid=None)

Band 1 of path as float64 with NaN for nodata, optionally warped onto grid.

Warping averages when the source is finer than grid and is bilinear otherwise.

write_raster(path, array, grid, nodata=np.nan)

Write array as a one-band LZW-compressed GeoTIFF on grid.

baseline_step(baseline, res)

Neighbour offset k (pixels) for a slope / roughness baseline of 2 k res; at least 1.

gradient(z, res, step=1)

Horn (1981) gradient (dz/dx, dz/dy) with neighbours step pixels away (gdaldem at 1).

Rows run north to south, so dz/dy is positive uphill to the north. Edge pixels are NaN.

slope_deg(z, res, step=1)

Slope (degrees) over a baseline of 2 step res; step=1 matches gdaldem slope.

vrm(z, res, step=1)

Vector ruggedness measure (Sappington et al. 2007) over a (2 step + 1)² window.

Unit surface normals come from the pixel-scale Horn gradient; VRM is 1 - |Σ n| / N over the window: 0 for a plane, up to 1 for rugged terrain.

as_product_images(images)

Accept ProductImage or (path, sun_azimuth[, resolution]) tuples.

sun_bin(azimuth, bin_width=30.0)

Index of the Sun-azimuth bin ([0, w), [w, 2w), …) holding azimuth.

illumination_maps(images, grid, *, shadow_threshold=0.005, bin_width=30.0, block_rows=1024)

Lit-image count, solar-bin count and best input resolution per pixel of grid.

Each image is warped onto grid in blocks of block_rows rows, so memory stays at a few arrays of the grid size.

copy_elevation(sfs_dem, out, *, force=False)

GLD01: the SfS DEM, LZW-compressed under its product name.

copy_height_error(height_error, out, *, force=False)

GLD03: ASP's sfs --estimate-height-errors map under its product name.

difference(sfs_dem, prior_dem, out, *, force=False)

GLD04: sfs_dem - prior_dem on the SfS grid (geodiff).

hillshade(dem, out, *, azimuth=315, elevation=45, force=False)

GLD07: ASP hillshade (Table 1 default: Sun azimuth 315°, elevation 45°).

orthomosaic(images, out, *, grid=None, resolution=None, work_dir=None, force=False)

GLD08: max-lit dem_mosaic --max of the mapprojected images, warped onto grid.

dem_mosaic cannot reproduce an arbitrary target grid exactly, so the mosaic is made on the inputs' grid and gdalwarp -r average puts it on the DEM grid (or at resolution when grid is unknown).

terrain_products(sfs_dem, out_dir, *, region_id, resolution, baselines=(None,), force=False)

GLD05 slope (deg) and GLD06 VRM roughness at each baseline.

None is the pixel scale (3×3 kernel), named with resolution as the baseline. A baseline b (m) uses neighbours k = round(b / 2 res) pixels away (baseline_step), so baselines up to 2 res equal the pixel scale. Keys are slope / roughness for the pixel scale and slope_<b> / roughness_<b> otherwise.

illumination_products(images, grid, out_dir, *, region_id, resolution, shadow_threshold=0.005, bin_width=30.0, force=False)

GLD02 lit-image count, GLD10 solar-bin count and GLD09 best input resolution (if known).

make_products(sfs_dem, prior_dem, out_dir, images, *, region_id, resolution, height_error=None, shadow_threshold=0.005, bin_width=30.0, baselines=(None,), hillshade_sun=(315, 45), force=False)

Table 1 products GLD01–GLD10 of an SfS DEM, named <region_id>_GLD<code>_<resolution>….

Parameters:
  • sfs_dem (Path) –

    The SfS DEM (e.g. sfs/run-DEM-final.tif).

  • prior_dem (Path) –

    Reference DEM the SfS DEM is differenced against (e.g. the LOLA LDEM).

  • out_dir (Path) –

    Destination; intermediates go to out_dir/work.

  • images (sequence) –

    Mapprojected SfS input images, as ProductImage or (path, sun_azimuth[, native_resolution]).

  • region_id (str) –

    ROIID of the file names.

  • resolution (float) –

    RES of the file names (m/px).

  • height_error (Path, default: None ) –

    <prefix>-height-error.tif from sfs --estimate-height-errors; GLD03 is skipped without it.

  • shadow_threshold (float, default: 0.005 ) –

    Image values at or below this are shadow (as sfs --shadow-threshold).

  • bin_width (float, default: 30.0 ) –

    Sun-azimuth bin width (deg) for the solar-bin count.

  • baselines (sequence, default: (None,) ) –

    Slope / roughness baselines (m); None is the pixel scale.

  • hillshade_sun ((azimuth, elevation), default: (315, 45) ) –

    Hillshade illumination (deg).

  • force (bool, default: False ) –

    Recompute products that already exist.
Returns:
  • dict –

    Product key (PRODUCT_CODES, plus slope_<b> / roughness_<b>) → path. Under dry_run only the CLI-tool products are planned; the raster-math ones need the real DEM.

crop_window(x, y, size=2500.0)

(ulx, uly, lrx, lry) of a size × size square centred on (x, y).

ensemble_subsets(pool, n_subsets, size, *, max_overlap=0.4, seed=0, max_tries=10000)

Random subsets of pool sharing at most max_overlap * size items pairwise.

Bertone et al. (2026, §3.2) limit the overlap of any two SfS ensemble subsets to 40% of their images. Subsets keep the order of pool (e.g. sorted by Sun azimuth). Fewer than n_subsets are returned, with a warning, if none more fit within max_tries draws.

ensemble_stats(members, grid=None)

Elevation and slope spread of ensemble DEMs, warped onto grid (default: the first's).

ensemble_uncertainty(members, out_dir, *, region_id, resolution, reference=None, force=False)

Ensemble elevation / slope uncertainty maps (Bertone et al. 2026, §3.2).

Parameters:
  • members (sequence of Path) –

    SfS DEMs of the same crop (e.g. 2.5 km, crop_window) from image subsets (ensemble_subsets).

  • out_dir (Path) –
  • region_id (str) –
  • resolution (float) –

    RES of the file names (m/px).

  • reference (Path, default: None ) –

    Raster whose grid the members are warped onto (e.g. the cropped prior DEM); default the first member's grid.

  • force (bool, default: False ) –
Returns:
  • dict –

    elevation_std → <ROI>_GLDISGM_<RES>_ENS.tif (m) and slope_std → <ROI>_GLDSSGM_<RES>_ENS.tif (deg).

compare_to_points(dem, points_csv, out_prefix, *, csv_format='1:lon 2:lat 3:height_above_datum', csv_srs=None, force=False)

Validate dem against altimetry points (e.g. LOLA RDR) with geodiff --csv-format.

Points off the DEM are dropped by geodiff. Returns an empty comparison under dry_run.