Coverage for lunar_toolkit/sfs_products.py: 94%
295 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 02:58 +0000
1"""Standard SfS DEM products and uncertainty estimates (Bertone et al. 2026, PSJ 7:118)."""
3from __future__ import annotations
5import logging
6import math
7import os
8import random
9import warnings
10from collections.abc import Hashable, Iterable, Sequence
11from dataclasses import dataclass
12from pathlib import Path
13from typing import Any, TypeVar
15import numpy as np
16import rasterio
18from . import asp
19from .tools import is_dry_run, reuse_existing, run
21log = logging.getLogger(__name__)
23T = TypeVar("T", bound=Hashable)
25# Table 1 of Bertone et al. (2026): product key -> GLD code.
26PRODUCT_CODES: dict[str, str] = {
27 "elevation": "ELEV",
28 "mask": "MASK",
29 "uncertainty": "ISGM",
30 "diff": "DIFF",
31 "slope": "TSLP",
32 "roughness": "TRGH",
33 "hillshade": "HILL",
34 "orthomosaic": "OMOS",
35 "best_resolution": "BRES",
36 "solar_bins": "SBCT",
37}
40def product_name(region_id: str, code: str, resolution: float, *suffix: Any) -> str:
41 """Table 1 file name, e.g. ``product_name("SR01", "TSLP", 5, 10)`` → ``SR01_GLDTSLP_5_10.tif``."""
42 parts = [region_id, f"GLD{code}", _label(resolution), *(_label(s) for s in suffix)]
43 return "_".join(parts) + ".tif"
46def _label(value: Any) -> str:
47 return f"{value:g}" if isinstance(value, (int, float)) else str(value)
50@dataclass(frozen=True)
51class Grid:
52 """A georeferenced raster grid (CRS, affine transform, size)."""
54 crs: Any
55 transform: Any
56 width: int
57 height: int
59 @classmethod
60 def from_raster(cls, path: Path) -> Grid:
61 with rasterio.open(path) as src:
62 return cls(src.crs, src.transform, src.width, src.height)
64 @property
65 def res(self) -> float:
66 """Pixel size (map units); pixels are assumed square."""
67 return abs(float(self.transform.a))
69 @property
70 def bounds(self) -> tuple[float, float, float, float]:
71 """``(xmin, ymin, xmax, ymax)``."""
72 t = self.transform
73 x0, y0 = t.c, t.f
74 x1, y1 = x0 + t.a * self.width, y0 + t.e * self.height
75 return min(x0, x1), min(y0, y1), max(x0, x1), max(y0, y1)
78def read_raster(path: Path, grid: Grid | None = None) -> tuple[np.ndarray, Grid]:
79 """Band 1 of *path* as float64 with NaN for nodata, optionally warped onto *grid*.
81 Warping averages when the source is finer than *grid* and is bilinear otherwise.
82 """
83 with rasterio.open(path) as src:
84 if grid is None:
85 data = src.read(1, masked=True)
86 return data.astype(np.float64).filled(np.nan), Grid(src.crs, src.transform, src.width, src.height)
87 with _warped(src, grid) as vrt:
88 data = vrt.read(1, masked=True)
89 return data.astype(np.float64).filled(np.nan), grid
92def write_raster(path: Path, array: np.ndarray, grid: Grid, nodata: float | None = np.nan) -> Path:
93 """Write *array* as a one-band LZW-compressed GeoTIFF on *grid*."""
94 path = Path(path)
95 path.parent.mkdir(parents=True, exist_ok=True)
96 profile = {
97 "driver": "GTiff",
98 "width": grid.width,
99 "height": grid.height,
100 "count": 1,
101 "dtype": array.dtype.name,
102 "crs": grid.crs,
103 "transform": grid.transform,
104 "compress": "lzw",
105 "nodata": nodata,
106 }
107 with rasterio.open(path, "w", **profile) as dst:
108 dst.write(array, 1)
109 return path
112def _warped(src: Any, grid: Grid) -> Any:
113 finer = abs(src.transform.a) < grid.res
114 return rasterio.vrt.WarpedVRT(
115 src,
116 crs=grid.crs,
117 transform=grid.transform,
118 width=grid.width,
119 height=grid.height,
120 resampling=rasterio.enums.Resampling.average if finer else rasterio.enums.Resampling.bilinear,
121 )
124def baseline_step(baseline: float | None, res: float) -> int:
125 """Neighbour offset *k* (pixels) for a slope / roughness baseline of ``2 k res``; at least 1."""
126 if baseline is None:
127 return 1
128 return max(1, round(baseline / (2 * res)))
131def _shift(a: np.ndarray, dy: int, dx: int) -> np.ndarray:
132 """``out[i, j] = a[i + dy, j + dx]``, NaN outside."""
133 out = np.full_like(a, np.nan)
134 h, w = a.shape
135 out[max(0, -dy) : h - max(0, dy), max(0, -dx) : w - max(0, dx)] = a[
136 max(0, dy) : h - max(0, -dy), max(0, dx) : w - max(0, -dx)
137 ]
138 return out
141def gradient(z: np.ndarray, res: float, step: int = 1) -> tuple[np.ndarray, np.ndarray]:
142 """Horn (1981) gradient ``(dz/dx, dz/dy)`` with neighbours *step* pixels away (``gdaldem`` at 1).
144 Rows run north to south, so ``dz/dy`` is positive uphill to the north. Edge pixels are NaN.
145 """
146 k = step
148 def s(dy: int, dx: int) -> np.ndarray:
149 return _shift(z, dy, dx)
151 dzdx = ((s(-k, k) + 2 * s(0, k) + s(k, k)) - (s(-k, -k) + 2 * s(0, -k) + s(k, -k))) / (8 * k * res)
152 dzdy = ((s(-k, -k) + 2 * s(-k, 0) + s(-k, k)) - (s(k, -k) + 2 * s(k, 0) + s(k, k))) / (8 * k * res)
153 return dzdx, dzdy
156def slope_deg(z: np.ndarray, res: float, step: int = 1) -> np.ndarray:
157 """Slope (degrees) over a baseline of ``2 step res``; ``step=1`` matches ``gdaldem slope``."""
158 dzdx, dzdy = gradient(z, res, step)
159 return np.degrees(np.arctan(np.hypot(dzdx, dzdy)))
162def _box_sum(a: np.ndarray, k: int) -> np.ndarray:
163 """Sum over the ``(2k+1)²`` window around each pixel; NaN if any value in it is NaN."""
164 n = 2 * k + 1
166 def window(x: np.ndarray) -> np.ndarray:
167 c = np.pad(np.pad(x, k).cumsum(0).cumsum(1), ((1, 0), (1, 0)))
168 return c[n:, n:] - c[:-n, n:] - c[n:, :-n] + c[:-n, :-n]
170 valid = np.isfinite(a)
171 total = window(np.where(valid, a, 0.0))
172 full = window(valid.astype(np.float64)) > n * n - 0.5
173 return np.where(full, total, np.nan)
176def vrm(z: np.ndarray, res: float, step: int = 1) -> np.ndarray:
177 """Vector ruggedness measure (Sappington et al. 2007) over a ``(2 step + 1)²`` window.
179 Unit surface normals come from the pixel-scale Horn gradient; VRM is
180 ``1 - |Σ n| / N`` over the window: 0 for a plane, up to 1 for rugged terrain.
181 """
182 p, q = gradient(z, res, 1)
183 norm = np.sqrt(1 + p * p + q * q)
184 sums = [_box_sum(c, step) for c in (-p / norm, -q / norm, 1 / norm)]
185 n = (2 * step + 1) ** 2
186 return 1 - np.sqrt(sum(s * s for s in sums)) / n
189@dataclass(frozen=True)
190class ProductImage:
191 """A mapprojected SfS input image, its Sun azimuth (deg from north) and native GSD (m/px)."""
193 path: Path
194 sun_azimuth: float
195 resolution: float | None = None
198def as_product_images(images: Iterable[ProductImage | Sequence[Any]]) -> list[ProductImage]:
199 """Accept [`ProductImage`][lunar_toolkit.sfs_products.ProductImage] or
200 ``(path, sun_azimuth[, resolution])`` tuples."""
201 out = []
202 for im in images:
203 if isinstance(im, ProductImage):
204 out.append(im)
205 else:
206 path, az, *rest = im
207 out.append(ProductImage(Path(path), float(az), float(rest[0]) if rest and rest[0] else None))
208 return out
211def sun_bin(azimuth: float, bin_width: float = 30.0) -> int:
212 """Index of the Sun-azimuth bin (``[0, w)``, ``[w, 2w)``, …) holding *azimuth*."""
213 return int((azimuth % 360.0) // bin_width)
216@dataclass
217class IlluminationMaps:
218 """Per-pixel illumination statistics of the SfS inputs on the DEM grid.
220 Attributes
221 ----------
222 count : ndarray (uint16)
223 Number of lit images (value above the shadow threshold).
224 solar_bins : ndarray (uint8)
225 Number of Sun-azimuth bins with at least one lit image.
226 best_resolution : ndarray (float32) or None
227 Finest native GSD among the lit images (NaN where none); None when
228 some image has no resolution.
229 grid : Grid
230 """
232 count: np.ndarray
233 solar_bins: np.ndarray
234 best_resolution: np.ndarray | None
235 grid: Grid
238def illumination_maps(
239 images: Sequence[ProductImage],
240 grid: Grid,
241 *,
242 shadow_threshold: float = 0.005,
243 bin_width: float = 30.0,
244 block_rows: int = 1024,
245) -> IlluminationMaps:
246 """Lit-image count, solar-bin count and best input resolution per pixel of *grid*.
248 Each image is warped onto *grid* in blocks of *block_rows* rows, so memory
249 stays at a few arrays of the grid size.
250 """
251 nbins = math.ceil(360.0 / bin_width)
252 if nbins > 64: 252 ↛ 253line 252 didn't jump to line 253 because the condition on line 252 was never true
253 raise ValueError(f"bin_width {bin_width} gives {nbins} bins; at most 64 are supported")
254 shape = (grid.height, grid.width)
255 count = np.zeros(shape, np.uint16)
256 bits = np.zeros(shape, np.uint64)
257 with_res = all(im.resolution for im in images)
258 best = np.full(shape, np.inf, np.float32) if with_res else None
259 for im in images:
260 bit = np.uint64(1) << np.uint64(sun_bin(im.sun_azimuth, bin_width))
261 with rasterio.open(im.path) as src, _warped(src, grid) as vrt:
262 for row in range(0, grid.height, block_rows):
263 rows = slice(row, min(row + block_rows, grid.height))
264 win = rasterio.windows.Window(0, row, grid.width, rows.stop - row)
265 data = vrt.read(1, window=win, masked=True)
266 lit = np.ma.filled(data > shadow_threshold, False)
267 count[rows] += lit
268 bits[rows] |= np.where(lit, bit, np.uint64(0))
269 if best is not None:
270 best[rows] = np.where(lit, np.minimum(best[rows], np.float32(im.resolution)), best[rows])
271 bins = np.zeros(shape, np.uint8)
272 for b in range(nbins):
273 bins += ((bits >> np.uint64(b)) & np.uint64(1)).astype(np.uint8)
274 if best is not None:
275 best[np.isinf(best)] = np.nan
276 return IlluminationMaps(count, bins, best, grid)
279def copy_elevation(sfs_dem: Path, out: Path, *, force: bool = False) -> Path:
280 """GLD01: the SfS DEM, LZW-compressed under its product name."""
281 if not reuse_existing(out, force):
282 asp.gdal_translate(sfs_dem, out)
283 return Path(out)
286def copy_height_error(height_error: Path, out: Path, *, force: bool = False) -> Path:
287 """GLD03: ASP's ``sfs --estimate-height-errors`` map under its product name."""
288 if not reuse_existing(out, force): 288 ↛ 290line 288 didn't jump to line 290 because the condition on line 288 was always true
289 asp.gdal_translate(height_error, out)
290 return Path(out)
293def difference(sfs_dem: Path, prior_dem: Path, out: Path, *, force: bool = False) -> Path:
294 """GLD04: ``sfs_dem - prior_dem`` on the SfS grid (``geodiff``)."""
295 out = Path(out)
296 if reuse_existing(out, force):
297 return out
298 prefix = out.with_suffix("")
299 diff = asp.geodiff(sfs_dem, prior_dem, prefix)
300 if not is_dry_run(): 300 ↛ 301line 300 didn't jump to line 301 because the condition on line 300 was never true
301 os.replace(diff, out)
302 return out
305def hillshade(
306 dem: Path, out: Path, *, azimuth: float = 315, elevation: float = 45, force: bool = False
307) -> Path:
308 """GLD07: ASP ``hillshade`` (Table 1 default: Sun azimuth 315°, elevation 45°)."""
309 if not reuse_existing(out, force):
310 asp.hillshade(dem, out, azimuth=azimuth, elevation=elevation)
311 return Path(out)
314def orthomosaic(
315 images: Sequence[Path],
316 out: Path,
317 *,
318 grid: Grid | None = None,
319 resolution: float | None = None,
320 work_dir: Path | None = None,
321 force: bool = False,
322) -> Path:
323 """GLD08: max-lit ``dem_mosaic --max`` of the mapprojected images, warped onto *grid*.
325 ``dem_mosaic`` cannot reproduce an arbitrary target grid exactly, so the
326 mosaic is made on the inputs' grid and ``gdalwarp -r average`` puts it on
327 the DEM grid (or at *resolution* when *grid* is unknown).
328 """
329 out = Path(out)
330 if reuse_existing(out, force):
331 return out
332 work_dir = Path(work_dir) if work_dir else out.parent / "work"
333 max_lit = work_dir / "max_lit.tif"
334 if not reuse_existing(max_lit, force): 334 ↛ 336line 334 didn't jump to line 336 because the condition on line 334 was always true
335 asp.dem_mosaic(images, max_lit, max=True)
336 if grid is not None:
337 target = ["-te", *grid.bounds, "-ts", grid.width, grid.height]
338 elif resolution is not None: 338 ↛ 341line 338 didn't jump to line 341 because the condition on line 338 was always true
339 target = ["-tr", resolution, resolution]
340 else:
341 raise ValueError("orthomosaic needs a grid or a resolution")
342 run("gdalwarp", ["-overwrite", "-r", "average", *target, "-co", "COMPRESS=LZW", max_lit, out])
343 return out
346def terrain_products(
347 sfs_dem: Path,
348 out_dir: Path,
349 *,
350 region_id: str,
351 resolution: float,
352 baselines: Sequence[float | None] = (None,),
353 force: bool = False,
354) -> dict[str, Path]:
355 """GLD05 slope (deg) and GLD06 VRM roughness at each baseline.
357 ``None`` is the pixel scale (3×3 kernel), named with *resolution* as the
358 baseline. A baseline *b* (m) uses neighbours ``k = round(b / 2 res)``
359 pixels away ([`baseline_step`][lunar_toolkit.sfs_products.baseline_step]), so baselines up to ``2 res``
360 equal the pixel scale. Keys are ``slope`` / ``roughness`` for the pixel scale and
361 ``slope_<b>`` / ``roughness_<b>`` otherwise.
362 """
363 out: dict[str, Path] = {}
364 grid = Grid.from_raster(sfs_dem)
365 z: np.ndarray | None = None
366 for baseline in baselines:
367 k = baseline_step(baseline, grid.res)
368 key = "" if baseline is None else f"_{baseline:g}"
369 label = resolution if baseline is None else baseline
370 for name, func in (("slope", slope_deg), ("roughness", vrm)):
371 path = Path(out_dir) / product_name(region_id, PRODUCT_CODES[name], resolution, label)
372 out[name + key] = path
373 if reuse_existing(path, force):
374 continue
375 if z is None:
376 z = read_raster(sfs_dem)[0]
377 write_raster(path, func(z, grid.res, k).astype(np.float32), grid)
378 return out
381def illumination_products(
382 images: Sequence[ProductImage],
383 grid: Grid,
384 out_dir: Path,
385 *,
386 region_id: str,
387 resolution: float,
388 shadow_threshold: float = 0.005,
389 bin_width: float = 30.0,
390 force: bool = False,
391) -> dict[str, Path]:
392 """GLD02 lit-image count, GLD10 solar-bin count and GLD09 best input resolution (if known)."""
393 names = {
394 k: Path(out_dir) / product_name(region_id, PRODUCT_CODES[k], resolution)
395 for k in ("mask", "solar_bins", "best_resolution")
396 }
397 with_res = all(im.resolution for im in images)
398 if not with_res:
399 log.info("Skipping best-resolution product: not every image has a native resolution")
400 del names["best_resolution"]
401 if all(reuse_existing(p, force) for p in names.values()): 401 ↛ 402line 401 didn't jump to line 402 because the condition on line 401 was never true
402 return names
403 maps = illumination_maps(images, grid, shadow_threshold=shadow_threshold, bin_width=bin_width)
404 write_raster(names["mask"], maps.count, grid, nodata=None)
405 write_raster(names["solar_bins"], maps.solar_bins, grid, nodata=None)
406 if maps.best_resolution is not None:
407 write_raster(names["best_resolution"], maps.best_resolution, grid)
408 return names
411def make_products(
412 sfs_dem: Path,
413 prior_dem: Path,
414 out_dir: Path,
415 images: Sequence[ProductImage | Sequence[Any]],
416 *,
417 region_id: str,
418 resolution: float,
419 height_error: Path | None = None,
420 shadow_threshold: float = 0.005,
421 bin_width: float = 30.0,
422 baselines: Sequence[float | None] = (None,),
423 hillshade_sun: tuple[float, float] = (315, 45),
424 force: bool = False,
425) -> dict[str, Path]:
426 """Table 1 products GLD01–GLD10 of an SfS DEM, named ``<region_id>_GLD<code>_<resolution>…``.
428 Parameters
429 ----------
430 sfs_dem : Path
431 The SfS DEM (e.g. ``sfs/run-DEM-final.tif``).
432 prior_dem : Path
433 Reference DEM the SfS DEM is differenced against (e.g. the LOLA LDEM).
434 out_dir : Path
435 Destination; intermediates go to ``out_dir/work``.
436 images : sequence
437 Mapprojected SfS input images, as [`ProductImage`][lunar_toolkit.sfs_products.ProductImage] or
438 ``(path, sun_azimuth[, native_resolution])``.
439 region_id : str
440 ``ROIID`` of the file names.
441 resolution : float
442 ``RES`` of the file names (m/px).
443 height_error : Path, optional
444 ``<prefix>-height-error.tif`` from ``sfs --estimate-height-errors``;
445 GLD03 is skipped without it.
446 shadow_threshold : float
447 Image values at or below this are shadow (as ``sfs --shadow-threshold``).
448 bin_width : float
449 Sun-azimuth bin width (deg) for the solar-bin count.
450 baselines : sequence
451 Slope / roughness baselines (m); ``None`` is the pixel scale.
452 hillshade_sun : (azimuth, elevation)
453 Hillshade illumination (deg).
454 force : bool
455 Recompute products that already exist.
457 Returns
458 -------
459 dict
460 Product key (`PRODUCT_CODES`, plus ``slope_<b>`` / ``roughness_<b>``)
461 → path. Under [`dry_run`][lunar_toolkit.tools.dry_run] only the CLI-tool
462 products are planned; the raster-math ones need the real DEM.
463 """
464 out_dir = Path(out_dir)
465 if not is_dry_run(): 465 ↛ 466line 465 didn't jump to line 466 because the condition on line 465 was never true
466 out_dir.mkdir(parents=True, exist_ok=True)
467 ims = as_product_images(images)
469 def name(key: str, *suffix: Any) -> Path:
470 return out_dir / product_name(region_id, PRODUCT_CODES[key], resolution, *suffix)
472 out: dict[str, Path] = {"elevation": copy_elevation(sfs_dem, name("elevation"), force=force)}
473 if height_error is not None and (is_dry_run() or Path(height_error).exists()):
474 out["uncertainty"] = copy_height_error(height_error, name("uncertainty"), force=force)
475 else:
476 log.info("Skipping elevation uncertainty: no SfS height-error map")
477 out["diff"] = difference(sfs_dem, prior_dem, name("diff"), force=force)
478 az, el = hillshade_sun
479 out["hillshade"] = hillshade(sfs_dem, name("hillshade", az, el), azimuth=az, elevation=el, force=force)
480 grid = Grid.from_raster(sfs_dem) if Path(sfs_dem).exists() else None
481 if ims:
482 out["orthomosaic"] = orthomosaic(
483 [im.path for im in ims],
484 name("orthomosaic"),
485 grid=grid,
486 resolution=resolution,
487 work_dir=out_dir / "work",
488 force=force,
489 )
490 if is_dry_run() or grid is None: 490 ↛ 493line 490 didn't jump to line 493 because the condition on line 490 was always true
491 log.info("Skipping raster-math products (dry run or missing SfS DEM)")
492 return out
493 out.update(
494 terrain_products(
495 sfs_dem, out_dir, region_id=region_id, resolution=resolution, baselines=baselines, force=force
496 )
497 )
498 if ims:
499 out.update(
500 illumination_products(
501 ims,
502 grid,
503 out_dir,
504 region_id=region_id,
505 resolution=resolution,
506 shadow_threshold=shadow_threshold,
507 bin_width=bin_width,
508 force=force,
509 )
510 )
511 return out
514def crop_window(x: float, y: float, size: float = 2500.0) -> tuple[float, float, float, float]:
515 """``(ulx, uly, lrx, lry)`` of a *size* × *size* square centred on ``(x, y)``."""
516 h = size / 2
517 return x - h, y + h, x + h, y - h
520def ensemble_subsets(
521 pool: Sequence[T],
522 n_subsets: int,
523 size: int,
524 *,
525 max_overlap: float = 0.4,
526 seed: int = 0,
527 max_tries: int = 10000,
528) -> list[list[T]]:
529 """Random subsets of *pool* sharing at most ``max_overlap * size`` items pairwise.
531 Bertone et al. (2026, §3.2) limit the overlap of any two SfS ensemble
532 subsets to 40% of their images. Subsets keep the order of *pool* (e.g.
533 sorted by Sun azimuth). Fewer than *n_subsets* are returned, with a
534 warning, if none more fit within *max_tries* draws.
535 """
536 if size > len(pool):
537 raise ValueError(f"subset size {size} exceeds the pool of {len(pool)}")
538 rng = random.Random(seed)
539 limit = math.floor(max_overlap * size + 1e-9)
540 chosen: list[set[int]] = []
541 for _ in range(max_tries):
542 if len(chosen) == n_subsets:
543 break
544 cand = set(rng.sample(range(len(pool)), size))
545 if all(len(cand & c) <= limit for c in chosen):
546 chosen.append(cand)
547 if len(chosen) < n_subsets:
548 log.warning("Only %d of %d subsets with <= %d shared images", len(chosen), n_subsets, limit)
549 return [[pool[i] for i in sorted(c)] for c in chosen]
552@dataclass
553class EnsembleStats:
554 """Per-pixel spread of an ensemble of SfS DEMs on a common grid.
556 Attributes
557 ----------
558 elevation_std, slope_std : ndarray
559 Sample standard deviation (``ddof=1``) of elevation (m) and pixel-scale
560 slope (deg); NaN where fewer than two members have data.
561 count : ndarray
562 Members with a valid elevation at each pixel.
563 grid : Grid
564 """
566 elevation_std: np.ndarray
567 slope_std: np.ndarray
568 count: np.ndarray
569 grid: Grid
572def ensemble_stats(members: Sequence[Path], grid: Grid | None = None) -> EnsembleStats:
573 """Elevation and slope spread of ensemble DEMs, warped onto *grid* (default: the first's)."""
574 if len(members) < 2:
575 raise ValueError("an ensemble needs at least two DEMs")
576 grid = grid or Grid.from_raster(members[0])
577 z = np.stack([read_raster(m, grid)[0] for m in members])
578 s = np.stack([slope_deg(zi, grid.res) for zi in z])
580 def std(a: np.ndarray) -> np.ndarray:
581 with warnings.catch_warnings():
582 warnings.simplefilter("ignore", RuntimeWarning)
583 return np.nanstd(a, axis=0, ddof=1)
585 return EnsembleStats(std(z), std(s), np.isfinite(z).sum(0).astype(np.uint16), grid)
588def ensemble_uncertainty(
589 members: Sequence[Path],
590 out_dir: Path,
591 *,
592 region_id: str,
593 resolution: float,
594 reference: Path | None = None,
595 force: bool = False,
596) -> dict[str, Path]:
597 """Ensemble elevation / slope uncertainty maps (Bertone et al. 2026, §3.2).
599 Parameters
600 ----------
601 members : sequence of Path
602 SfS DEMs of the same crop (e.g. 2.5 km, [`crop_window`][lunar_toolkit.sfs_products.crop_window]) from
603 image subsets ([`ensemble_subsets`][lunar_toolkit.sfs_products.ensemble_subsets]).
604 out_dir : Path
605 region_id : str
606 resolution : float
607 ``RES`` of the file names (m/px).
608 reference : Path, optional
609 Raster whose grid the members are warped onto (e.g. the cropped prior
610 DEM); default the first member's grid.
611 force : bool
613 Returns
614 -------
615 dict
616 ``elevation_std`` → ``<ROI>_GLDISGM_<RES>_ENS.tif`` (m) and
617 ``slope_std`` → ``<ROI>_GLDSSGM_<RES>_ENS.tif`` (deg).
618 """
619 out = {
620 "elevation_std": Path(out_dir) / product_name(region_id, "ISGM", resolution, "ENS"),
621 "slope_std": Path(out_dir) / product_name(region_id, "SSGM", resolution, "ENS"),
622 }
623 if all(reuse_existing(p, force) for p in out.values()): 623 ↛ 624line 623 didn't jump to line 624 because the condition on line 623 was never true
624 return out
625 stats = ensemble_stats(members, Grid.from_raster(reference) if reference else None)
626 write_raster(out["elevation_std"], stats.elevation_std.astype(np.float32), stats.grid)
627 write_raster(out["slope_std"], stats.slope_std.astype(np.float32), stats.grid)
628 return out
631@dataclass
632class PointComparison:
633 """DEM minus reference points, as written by ``geodiff`` (lon, lat in degrees; diff in m)."""
635 lon: np.ndarray
636 lat: np.ndarray
637 diff: np.ndarray
639 def stats(self) -> dict[str, float]:
640 d = self.diff[np.isfinite(self.diff)]
641 if d.size == 0:
642 return {"count": 0}
643 med = float(np.median(d))
644 return {
645 "count": int(d.size),
646 "mean": float(d.mean()),
647 "median": med,
648 "std": float(d.std()),
649 "rms": float(np.sqrt(np.mean(d * d))),
650 "nmad": float(1.4826 * np.median(np.abs(d - med))),
651 "min": float(d.min()),
652 "max": float(d.max()),
653 }
656def compare_to_points(
657 dem: Path,
658 points_csv: Path,
659 out_prefix: Path,
660 *,
661 csv_format: str = "1:lon 2:lat 3:height_above_datum",
662 csv_srs: str | None = None,
663 force: bool = False,
664) -> PointComparison:
665 """Validate *dem* against altimetry points (e.g. LOLA RDR) with ``geodiff --csv-format``.
667 Points off the DEM are dropped by ``geodiff``. Returns an empty
668 comparison under [`dry_run`][lunar_toolkit.tools.dry_run].
669 """
670 diff_csv = Path(f"{out_prefix}-diff.csv")
671 if not reuse_existing(diff_csv, force):
672 asp.geodiff(dem, points_csv, out_prefix, csv_format=csv_format, csv_srs=csv_srs)
673 if is_dry_run():
674 empty = np.empty(0)
675 return PointComparison(empty, empty, empty)
676 data = np.loadtxt(diff_csv, delimiter=",", comments="#", ndmin=2)
677 return PointComparison(data[:, 0], data[:, 1], data[:, 2])