Coverage for lunar_toolkit/csm.py: 94%
162 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 02:58 +0000
1"""Read, write and resample USGS CSM linescan cameras (ISD and model state) to match ISIS ``reduce`` / ``crop``."""
3from __future__ import annotations
5import json
6import math
7import re
8from dataclasses import dataclass
9from pathlib import Path
10from typing import Any
12from . import isis
13from .tools import is_dry_run
15LINESCAN = "USGS_ASTRO_LINE_SCANNER_SENSOR_MODEL"
18@dataclass
19class CsmCamera:
20 """A CSM camera file: an ALE / ``isd_generate`` ISD or a USGS CSM model state.
22 Attributes
23 ----------
24 data : dict
25 The parsed JSON (ISD keys such as ``image_lines`` or state keys such
26 as ``m_nLines``).
27 kind : {"isd", "state"}
28 header : str or None
29 For a state, the model-name line that precedes the JSON in
30 ``getModelState()`` output and ASP ``.adjusted_state.json`` files
31 (``None`` when the file is bare JSON).
32 """
34 data: dict[str, Any]
35 kind: str
36 header: str | None = None
38 @property
39 def model(self) -> str:
40 key = "name_model" if self.kind == "isd" else "m_modelName"
41 return str(self.data.get(key, self.header or ""))
43 @property
44 def size(self) -> tuple[int, int]:
45 """``(samples, lines)`` of the image the camera describes."""
46 if self.kind == "isd":
47 return int(self.data["image_samples"]), int(self.data["image_lines"])
48 return int(self.data["m_nSamples"]), int(self.data["m_nLines"])
51def parse_camera(text: str) -> CsmCamera:
52 """Parse a CSM camera from text (ISD JSON, or a model state with or without its name line)."""
53 stripped = text.lstrip()
54 header = None
55 if not stripped.startswith("{"):
56 first, _, rest = stripped.partition("\n")
57 header, stripped = first.strip(), rest
58 data = json.loads(stripped)
59 if not isinstance(data, dict): 59 ↛ 60line 59 didn't jump to line 60 because the condition on line 59 was never true
60 raise ValueError("CSM camera JSON must be an object")
61 if "m_nLines" in data:
62 return CsmCamera(data, "state", header)
63 if "image_lines" in data: 63 ↛ 64line 63 didn't jump to line 64 because the condition on line 63 was never true
64 return CsmCamera(data, "isd", None)
65 raise ValueError("Not a CSM ISD (image_lines) or model state (m_nLines)")
68def read_camera(path: Path) -> CsmCamera:
69 """Read a CSM camera file (``.json`` ISD, model state or ASP ``.adjusted_state.json``)."""
70 return parse_camera(Path(path).read_text())
73def format_camera(cam: CsmCamera) -> str:
74 """Text of *cam* in its original layout (state name line kept)."""
75 body = json.dumps(cam.data, indent=2)
76 return f"{cam.header}\n{body}\n" if cam.header else body + "\n"
79def write_camera(cam: CsmCamera, path: Path) -> Path:
80 path = Path(path)
81 path.parent.mkdir(parents=True, exist_ok=True)
82 path.write_text(format_camera(cam))
83 return path
86def transform_camera(
87 cam: CsmCamera,
88 sample_scale: float,
89 sample_offset: float,
90 line_scale: float,
91 line_offset: float,
92 size: tuple[int, int],
93) -> CsmCamera:
94 """Camera for an image whose pixels map affinely onto *cam*'s image.
96 CSM image coordinates put the top-left pixel edge at 0 (first pixel
97 centre at 0.5). A point at ``(x, y)`` of the new image is
98 ``(sample_scale * x + sample_offset, line_scale * y + line_offset)`` of
99 the old one. Across track this scales the detector sample summing and
100 shifts the starting detector sample; along track it scales the
101 integration times (and the line summing, which converts the residual
102 detector-line offset into image lines in ``groundToImage``) and moves the
103 integration start lines so that line ``y`` is exposed when old line
104 ``line_scale * y + line_offset`` was.
106 Parameters
107 ----------
108 cam : CsmCamera
109 A USGS linescan ISD or state.
110 sample_scale, sample_offset, line_scale, line_offset : float
111 The mapping above (ISIS ``reduce`` by *f*: scales *f*, offsets 0;
112 ISIS ``crop`` from sample *s*, line *l*: scales 1, offsets
113 ``s - 1``, ``l - 1``).
114 size : (samples, lines)
115 Size of the new image.
117 Returns
118 -------
119 CsmCamera
120 A new camera of the same kind.
121 """
122 if cam.model != LINESCAN:
123 raise ValueError(f"Only {LINESCAN} cameras are supported, not {cam.model or 'unknown'}")
124 if sample_scale <= 0 or line_scale <= 0:
125 raise ValueError("Scales must be positive")
126 d = json.loads(json.dumps(cam.data))
127 ns, nl = int(size[0]), int(size[1])
129 def new_lines(lines: list[float]) -> list[float]:
130 return [(v - 0.5 - line_offset) / line_scale + 0.5 for v in lines]
132 if cam.kind == "isd":
133 ss = float(d.get("detector_sample_summing", 1.0))
134 d["starting_detector_sample"] = float(d.get("starting_detector_sample", 0.0)) + sample_offset * ss
135 _isd_summing(d, "sample", ss * sample_scale)
136 _isd_summing(d, "line", float(d.get("detector_line_summing", 1.0)) * line_scale)
137 rates = d["line_scan_rate"]
138 d["line_scan_rate"] = [
139 [new_lines([float(r[0])])[0], float(r[1]), float(r[2]) * line_scale] for r in rates
140 ]
141 d["image_samples"], d["image_lines"] = ns, nl
142 return CsmCamera(d, "isd")
144 ss = float(d["m_detectorSampleSumming"])
145 d["m_startingDetectorSample"] = float(d["m_startingDetectorSample"]) + sample_offset * ss
146 d["m_detectorSampleSumming"] = ss * sample_scale
147 d["m_detectorLineSumming"] = float(d["m_detectorLineSumming"]) * line_scale
148 d["m_intTimeLines"] = new_lines([float(v) for v in d["m_intTimeLines"]])
149 d["m_intTimes"] = [float(v) * line_scale for v in d["m_intTimes"]]
150 d["m_nSamples"], d["m_nLines"] = ns, nl
151 if "m_gsd" in d: 151 ↛ 155line 151 didn't jump to line 155 because the condition on line 151 was always true
152 d["m_gsd"] = float(d["m_gsd"]) * math.sqrt(sample_scale * line_scale)
153 if "m_halfSwath" in d: 153 ↛ 155line 153 didn't jump to line 155 because the condition on line 153 was always true
154 d["m_halfSwath"] = d["m_gsd"] * ns / 2.0
155 starts, times, lines = d["m_intTimeStartTimes"], d["m_intTimes"], d["m_intTimeLines"]
156 if "m_halfTime" in d and starts: 156 ↛ 158line 156 didn't jump to line 158 because the condition on line 156 was always true
157 d["m_halfTime"] = (starts[-1] - starts[0] + times[-1] * (nl - lines[-1])) / 2.0
158 return CsmCamera(d, "state", cam.header)
161def _isd_summing(d: dict[str, Any], axis: str, summing: float) -> None:
162 """Set an ISD detector summing; usgscsm reads it as an integer, so a fractional
163 remainder ``k`` goes into that axis' focal2pixel row, detector centre and
164 starting detector pixel (all divided by ``k``, which leaves focal-plane
165 coordinates unchanged)."""
166 n = max(1, round(summing))
167 d[f"detector_{axis}_summing"] = n
168 k = summing / n
169 if abs(k - 1.0) < 1e-12:
170 return
171 key = "lines" if axis == "line" else "samples"
172 d[f"focal2pixel_{key}"] = [float(v) / k for v in d[f"focal2pixel_{key}"]]
173 centre = dict(d.get("detector_center", {}))
174 centre[axis] = float(centre.get(axis, 0.0)) / k
175 d["detector_center"] = centre
176 d[f"starting_detector_{axis}"] = float(d.get(f"starting_detector_{axis}", 0.0)) / k
179def reduce_camera(
180 camera: Path | CsmCamera,
181 factor: float,
182 out_json: Path | None = None,
183 size: tuple[int, int] | None = None,
184) -> CsmCamera:
185 """CSM camera of an image averaged by *factor* (ISIS ``reduce``, ``mode=scale``).
187 ``reduce`` output pixel *k* (0-based) averages input pixels
188 ``k*f .. k*f + f - 1``, so in CSM coordinates (pixel edges at integers)
189 the mapping is a pure scale by *f*: summing and integration times grow
190 by *f*, sizes shrink by *f*.
192 Parameters
193 ----------
194 camera : Path or CsmCamera
195 ISD or state of the full-resolution image.
196 factor : float
197 ``sscale`` = ``lscale``.
198 out_json : Path, optional
199 Where to write the result (same format as the input).
200 size : (samples, lines), optional
201 Size of the reduced image; default ``ceil(size / factor)``.
202 """
203 cam = camera if isinstance(camera, CsmCamera) else read_camera(camera)
204 if size is None:
205 ns, nl = cam.size
206 size = (math.ceil(ns / factor), math.ceil(nl / factor))
207 out = transform_camera(cam, factor, 0.0, factor, 0.0, size)
208 if out_json is not None: 208 ↛ 210line 208 didn't jump to line 210 because the condition on line 208 was always true
209 write_camera(out, out_json)
210 return out
213def read_isis_label(cube: Path, max_bytes: int = 1 << 22) -> str:
214 """The attached PVL label of an ISIS cube (text up to the closing ``End``)."""
215 with open(cube, "rb") as f:
216 head = f.read(max_bytes).decode("latin-1")
217 m = re.search(r"^End\s*$", head, flags=re.M)
218 return head[: m.end()] if m else head
221def _group(label: str, name: str) -> dict[str, str] | None:
222 m = re.search(rf"^\s*Group\s*=\s*{name}\s*$(.*?)^\s*End_Group", label, flags=re.M | re.S)
223 if not m:
224 return None
225 return dict(re.findall(r"^\s*(\w+)\s*=\s*(\S+)", m.group(1), flags=re.M))
228@dataclass(frozen=True)
229class CubeGeometry:
230 """Size of an ISIS cube and its pixel mapping onto the original (alpha) image.
232 ``alpha = scale * beta + offset`` per axis in CSM coordinates (pixel
233 edges at integers), from the cube's ``AlphaCube`` group (identity when
234 absent).
235 """
237 samples: int
238 lines: int
239 alpha_samples: int
240 alpha_lines: int
241 sample_scale: float = 1.0
242 sample_offset: float = 0.0
243 line_scale: float = 1.0
244 line_offset: float = 0.0
246 @classmethod
247 def from_label(cls, label: str) -> CubeGeometry:
248 dims = _group(label, "Dimensions")
249 if not dims:
250 raise ValueError("No Dimensions group in the ISIS label")
251 ns, nl = int(dims["Samples"]), int(dims["Lines"])
252 a = _group(label, "AlphaCube")
253 if not a:
254 return cls(ns, nl, ns, nl)
255 s0, s1 = float(a["AlphaStartingSample"]), float(a["AlphaEndingSample"])
256 l0, l1 = float(a["AlphaStartingLine"]), float(a["AlphaEndingLine"])
257 bs, bl = int(a["BetaSamples"]), int(a["BetaLines"])
258 return cls(
259 ns,
260 nl,
261 int(a["AlphaSamples"]),
262 int(a["AlphaLines"]),
263 (s1 - s0) / bs,
264 s0 - 0.5,
265 (l1 - l0) / bl,
266 l0 - 0.5,
267 )
269 @classmethod
270 def from_cube(cls, cube: Path) -> CubeGeometry:
271 return cls.from_label(read_isis_label(cube))
274def camera_for_cube(cam: CsmCamera, cube: CubeGeometry, source: CubeGeometry | None = None) -> CsmCamera:
275 """Re-map *cam* so it describes *cube* pixel for pixel.
277 *cam* is taken to describe the original (alpha) image when its size
278 matches ``cube.alpha_*`` (cameras from ``isd_generate`` on the
279 uncropped cube, as in the ASP LRONAC example), else the *source* cube
280 that *cube* was derived from (``reduce`` / ``crop`` of a cube whose own
281 camera was written for it).
282 """
283 size = cam.size
284 ss, so, ls, lo = cube.sample_scale, cube.sample_offset, cube.line_scale, cube.line_offset
285 if size != (cube.alpha_samples, cube.alpha_lines):
286 if source is None or size != (source.samples, source.lines):
287 raise ValueError(
288 f"Camera size {size} matches neither the alpha image "
289 f"{(cube.alpha_samples, cube.alpha_lines)} nor the source cube"
290 )
291 # cube -> alpha, then alpha -> source (inverse of the source's AlphaCube)
292 ss, so = ss / source.sample_scale, (so - source.sample_offset) / source.sample_scale
293 ls, lo = ls / source.line_scale, (lo - source.line_offset) / source.line_scale
294 return transform_camera(cam, ss, so, ls, lo, (cube.samples, cube.lines))
297def isis_reduce_with_csm(
298 cube: Path,
299 camera: Path,
300 factor: int,
301 out_dir: Path,
302 algorithm: str = "average",
303 force: bool = False,
304) -> tuple[Path, Path]:
305 """ISIS ``reduce`` a cube by *factor* and write the matching CSM camera.
307 Parameters
308 ----------
309 cube : Path
310 Input ISIS cube (may be a crop: its ``AlphaCube`` is honoured).
311 camera : Path
312 CSM ISD or state of the original image or of *cube*.
313 factor : int
314 ``sscale`` = ``lscale``.
315 out_dir : Path
316 Receives ``<cube stem>.red<factor>.cub`` and ``<camera stem>.red<factor>.json``
317 (same CSM format as *camera*).
318 algorithm : str
319 ``reduce`` algorithm (``average`` as in Bertone et al. 2026, App. A).
320 force : bool
321 Redo existing outputs.
323 Returns
324 -------
325 (reduced_cube, reduced_camera)
326 """
327 out_dir = Path(out_dir)
328 stem = Path(cube).stem
329 red = out_dir / f"{stem}.red{factor}.cub"
330 out_cam = out_dir / f"{Path(camera).stem}.red{factor}.json" # named after its source camera
331 if force or not red.exists():
332 if not is_dry_run(): 332 ↛ 333line 332 didn't jump to line 333 because the condition on line 332 was never true
333 out_dir.mkdir(parents=True, exist_ok=True)
334 isis.reduce(cube, red, factor, algorithm=algorithm)
335 if is_dry_run():
336 return red, out_cam
337 if force or not out_cam.exists(): 337 ↛ 338line 337 didn't jump to line 338 because the condition on line 337 was never true
338 geom = CubeGeometry.from_cube(red)
339 write_camera(camera_for_cube(read_camera(camera), geom, CubeGeometry.from_cube(cube)), out_cam)
340 return red, out_cam