Coverage for python/lsst/images/tests/_minify_for_fixtures.py: 0%
175 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-06 08:34 +0000
1# This file is part of lsst-images.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# Use of this source code is governed by a 3-clause BSD-style
10# license that can be found in the LICENSE file.
12"""Minify a real on-disk archive into a small JSON test fixture.
14Reads a FITS or NDF file via the appropriate input archive, takes a
15small subset of the in-memory object, and writes JSON via
16``JsonOutputArchive``. Used to populate ``tests/data/schema_v1/legacy/``
17with derived-from-real test data that exercises the full read path
18including the absence-of-stamp legacy default.
20Per top-level type the subset rule is:
22 VisitImage Crop the image/mask/variance planes to a small (~16x16)
23 corner, keeping the real single-instance structures (PSF
24 such as Piff, detector frames) that synthetic fixtures
25 cannot reproduce -- the whole point of deriving a fixture
26 from real data. Homogeneous repeated collections (detector
27 amplifiers, aperture-correction entries) are trimmed to a
28 representative few, since one entry exercises the schema as
29 well as sixteen. The projection's pixel->sky mapping is
30 replaced by its linear (affine) approximation over the kept
31 box: a real TAN-SIP WCS serializes as a ~100 KB AST
32 polynomial dump, but over a 16x16 box it is linear to far
33 below a pixel, so the affine form is schema-identical and
34 orders of magnitude smaller. A Piff PSF's field
35 interpolation is truncated to a low order (the order-4
36 solution table is ~225 KB; order 0, the field-averaged PSF,
37 is schema-identical and ~13x smaller).
39 CellCoadd Crop to a small block of cells (preferring a block that
40 includes a missing cell so the sparse-grid path is
41 exercised) and then *morph* that block onto a tiny cell
42 grid: each cell's planes are decimated from the native
43 cell size down to a few pixels and re-stitched, and the
44 PSF kernels are cropped to a small odd window. The grid
45 topology (number of cells, the missing-cell set, band,
46 mask schema and provenance shape) is preserved; the pixel
47 values and WCS are *not* physically meaningful. This is
48 the "morph cells in place" fallback: it sidesteps the
49 outer-ring problem (inputs/PSFs that overlap kept cells)
50 by rebuilding a self-consistent miniature coadd rather
51 than trying to carve an accurate subset out of the real
52 one. An accurate per-cell subset would inline several
53 150x150 planes per cell and produce multi-megabyte JSON,
54 which defeats the purpose of a fixture.
56Run interactively (CellCoadd works with just this package installed;
57VisitImage needs a full Rubin environment so the real PSF can be read)::
59 python -c "
60 from lsst.images.tests._minify_for_fixtures import minify
61 minify('cell_example.fits', 'tests/data/schema_v1/legacy/cell_coadd.json')
62 minify('dp1.fits', 'tests/data/schema_v1/legacy/visit_image_dp1.json')
63 minify('dp2.fits', 'tests/data/schema_v1/legacy/visit_image_dp2.json')
64 "
66The helper is invoked manually by developers when they have a real
67on-disk file to derive from; it is not exercised by CI.
68"""
70from __future__ import annotations
72__all__ = ("minify",)
74import json
75import os
76from collections.abc import Callable
77from typing import Any
79import numpy as np
81from .. import DifferenceImage, VisitImage
82from .. import json as images_json
83from .._cell_grid import CellGrid, CellGridBounds, CellIJ, PatchDefinition
84from .._geom import YX, Box
85from .._image import Image
86from .._mask import Mask
87from .._transforms import Projection, TractFrame, Transform
88from .._transforms._ast import PolyMap
89from ..cells import CellCoadd
90from ..cells._provenance import CoaddProvenance
91from ..cells._psf import CellPointSpreadFunction
92from ..psfs import PiffWrapper
93from ..serialization import backend_for_path, read
94from ._creation import make_random_projection
96# Default morph parameters for CellCoadd. ``CELL_SIZE`` should divide the
97# native cell size evenly; ``KERNEL_SIZE`` must be odd. ``MAX_INPUTS`` caps
98# the provenance ``inputs`` table (a real coadd has hundreds of visits); the
99# full provenance schema is already exercised by the ``coadd_provenance``
100# fixture, so here we keep just enough rows to be representative.
101_CELL_SIZE = 6
102_KERNEL_SIZE = 5
103_MAX_INPUTS = 6
105# Default trim parameters for VisitImage. Amplifiers and aperture-correction
106# entries are homogeneous collections, so a couple of each cover the schema
107# just as well as the full set (a real detector has 16 amplifiers and dozens
108# of aperture corrections).
109_MAX_AMPLIFIERS = 2
110_MAX_APERTURE_CORRECTIONS = 2
112# Field-interpolation order to truncate a Piff PSF to (the solution table of a
113# real order-4 PixelGrid PSF dominates the fixture at ~225 KB). Order 0 is the
114# field-averaged PSF; set to `None` to leave the PSF untouched.
115_PSF_INTERP_ORDER = 0
117# Maximum permitted deviation (radians) when approximating a projection's
118# pixel->sky mapping with an affine one. Over a fixture's tiny box the real
119# mapping is linear well below this, so the fit always succeeds.
120_PROJECTION_LINEAR_APPROX_TOL = 1e-8
123def minify(in_path: str, out_path: str, *, schema_name: str | None = None) -> None:
124 """Read a real archive at ``in_path``, take a small subset, and write JSON.
126 Parameters
127 ----------
128 in_path
129 Path to a FITS (``.fits`` / ``.fits.gz``) or NDF (``.sdf`` / ``.ndf``)
130 file to read.
131 out_path
132 Path to the JSON fixture to write. The parent directory is
133 created if it does not exist.
134 schema_name
135 Top-level schema name (e.g. ``"visit_image"`` or ``"cell_coadd"``).
136 If `None`, it is auto-detected from the file.
138 Raises
139 ------
140 ValueError
141 If the file extension is not recognised.
142 NotImplementedError
143 If the top-level type is not one this helper knows how to subset.
144 """
145 backend = backend_for_path(in_path)
146 if schema_name is None:
147 schema_name = backend.input_archive.get_basic_info(in_path).schema_name
149 cls, subsetter = _dispatch(schema_name)
151 obj: Any = read(in_path, cls)
152 subset = subsetter(obj)
154 tree = images_json.write(subset)
155 dumped = tree.model_dump(mode="json")
156 os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True)
157 with open(out_path, "w") as stream:
158 stream.write(json.dumps(dumped, indent=2, sort_keys=False) + "\n")
161def _dispatch(schema_name: str) -> tuple[type, Callable[[Any], Any]]:
162 """Return the ``(class, subsetter)`` pair for a top-level schema name."""
163 registry: dict[str, tuple[type, Callable[[Any], Any]]] = {
164 "visit_image": (VisitImage, _subset_visit_image),
165 "difference_image": (DifferenceImage, _subset_visit_image),
166 "cell_coadd": (CellCoadd, _subset_cell_coadd),
167 }
168 try:
169 return registry[schema_name]
170 except KeyError:
171 raise NotImplementedError(
172 f"No minify rule for schema {schema_name!r}; supported: {sorted(registry)}."
173 ) from None
176# -- VisitImage ------------------------------------------------------------
179def _subset_visit_image(
180 visit_like_image: VisitImage | DifferenceImage,
181 *,
182 size: int = 16,
183 max_amplifiers: int = _MAX_AMPLIFIERS,
184 max_aperture_corrections: int = _MAX_APERTURE_CORRECTIONS,
185 linearize_projection: bool = True,
186 projection_tol: float = _PROJECTION_LINEAR_APPROX_TOL,
187 psf_interp_order: int | None = _PSF_INTERP_ORDER,
188) -> VisitImage | DifferenceImage:
189 """Crop a VisitImage's pixel planes to a small corner and trim its
190 homogeneous collections.
192 The detector frames are a single structure carried through unchanged by
193 ``__getitem__``. The detector's amplifiers and the aperture-correction map
194 are repeated, schema-identical entries, so they are trimmed to a
195 representative few. The projection's pixel->sky mapping is replaced by its
196 affine approximation over the kept box (see ``_linear_approx_projection``)
197 unless ``linearize_projection`` is false. A Piff PSF's field interpolation
198 is truncated to ``psf_interp_order`` (see ``_simplify_piff_psf``) unless
199 that is `None`.
200 """
201 bbox = visit_like_image.bbox
202 y0 = bbox.y.start
203 x0 = bbox.x.start
204 y1 = min(y0 + size, bbox.y.stop)
205 x1 = min(x0 + size, bbox.x.stop)
206 subset = visit_like_image[Box.factory[y0:y1, x0:x1]]
208 # ``subset`` is a fresh throwaway object whose detector amplifier list,
209 # aperture-correction map and PSF are live, mutable components. Trim them
210 # in place through the public accessors rather than reaching for private
211 # attributes.
212 del subset.detector.amplifiers[max_amplifiers:]
213 aperture_corrections = subset.aperture_corrections
214 for key in list(aperture_corrections)[max_aperture_corrections:]:
215 del aperture_corrections[key]
216 if psf_interp_order is not None and isinstance(subset.psf, PiffWrapper):
217 _simplify_piff_psf(subset.psf, order=psf_interp_order)
219 if not linearize_projection or subset.projection is None:
220 return subset
222 # The pixel planes carry the projection immutably (there is no public
223 # setter for it), so install the affine approximation by rebuilding the
224 # VisitImage from its public components with re-viewed planes. Only the
225 # image plane's projection is actually serialized, but keeping all three
226 # consistent avoids surprises.
227 linear = _linear_approx_projection(subset.projection, subset.image.bbox, tol=projection_tol)
228 return type(visit_like_image)(
229 subset.image.view(projection=linear),
230 mask=subset.mask.view(projection=linear),
231 variance=subset.variance.view(projection=linear),
232 projection=linear,
233 psf=subset.psf,
234 obs_info=subset.obs_info,
235 bounds=subset.bounds,
236 summary_stats=subset.summary_stats,
237 detector=subset.detector,
238 photometric_scaling=subset.photometric_scaling,
239 aperture_corrections=subset.aperture_corrections,
240 backgrounds=subset.backgrounds,
241 band=subset.band,
242 metadata=subset.metadata,
243 )
246def _linear_approx_projection(projection: Projection, bbox: Box, *, tol: float) -> Projection:
247 """Return a copy of ``projection`` whose pixel->sky mapping is replaced by
248 its best linear (affine) approximation over ``bbox``.
250 Real WCS mappings (e.g. TAN-SIP) serialize as large AST polynomial dumps.
251 Over the small box of a fixture they are linear to far below a pixel, so
252 an affine approximation is schema-identical but orders of magnitude
253 smaller. The result carries no FITS approximation (the affine is itself
254 trivially FITS-representable).
256 This is written as a self-contained ``projection -> projection`` transform
257 so it can be promoted to a public ``Projection.linear_approx(bbox, tol)``
258 method later with essentially no change. It assumes a 2-D pixel->sky
259 mapping.
261 Parameters
262 ----------
263 projection
264 The projection to approximate.
265 bbox
266 Box (in pixel coordinates) over which the approximation must hold.
267 tol
268 Maximum permitted deviation from linearity, as a Cartesian
269 displacement in the output (sky, radians) coordinates. AST raises
270 ``RuntimeError`` if no fit within ``tol`` exists.
271 """
272 transform = projection.pixel_to_sky_transform
273 mapping = transform._ast_mapping
274 lbnd = [bbox.x.start, bbox.y.start]
275 ubnd = [bbox.x.stop, bbox.y.stop]
276 # linearApprox yields [offsets; Jacobian] as a (1 + n_out, n_in) array on
277 # both AST backends (astshim returns the flat buffer in the same order, so
278 # the reshape recovers the same layout the starlink-pyast bridge returns).
279 fit = np.asarray(mapping.linearApprox(lbnd, ubnd, tol), dtype=float).reshape(3, 2)
280 offset = fit[0] # (lon0, lat0), radians
281 jacobian = fit[1:] # jacobian[i, j] = d(out_i) / d(in_j), in = (x, y)
282 jacobian_inv = np.linalg.inv(jacobian)
283 forward = _affine_polymap_coeffs(jacobian, offset)
284 inverse = _affine_polymap_coeffs(jacobian_inv, -jacobian_inv @ offset)
285 affine = Transform(
286 transform.in_frame,
287 transform.out_frame,
288 PolyMap(forward, inverse),
289 in_bounds=projection.pixel_bounds,
290 )
291 return affine.as_projection()
294def _affine_polymap_coeffs(matrix: np.ndarray, offset: np.ndarray) -> np.ndarray:
295 """Build AST ``PolyMap`` coefficients for ``out = matrix @ in + offset``.
297 Each row is ``[coefficient, output_axis (1-based), power_of_in_1, ...]``;
298 one constant row plus one row per input per output axis. Returned as a
299 float array, which is the form both AST backends require.
300 """
301 n = len(offset)
302 coeffs: list[list[float]] = []
303 for i in range(n):
304 coeffs.append([float(offset[i]), i + 1, *([0] * n)])
305 for j in range(n):
306 powers = [1 if k == j else 0 for k in range(n)]
307 coeffs.append([float(matrix[i][j]), i + 1, *powers])
308 return np.array(coeffs, dtype=float)
311def _simplify_piff_psf(psf: PiffWrapper, *, order: int) -> None:
312 """Truncate a Piff PSF's field interpolation to ``order``, in place.
314 A real Piff PSF interpolates a per-pixel model across the focal plane with
315 a high-order 2-D polynomial; that solution table dominates the serialized
316 size (a 25x25 PixelGrid x order-4 polynomial is ~225 KB). Truncating to
317 ``order`` keeps only the lowest-order field terms -- order 0 is the
318 field-averaged PSF -- which is schema-identical but far smaller, and needs
319 no stars or refit (the fitted ``stars`` are already dropped on serialize).
321 Only ``BasisPolynomial``-interpolated PSFs are handled; anything else (a
322 higher-order model already at/under ``order``, a non-polynomial interp) is
323 left untouched.
325 ``piff`` is imported lazily because it is an optional dependency; this is
326 only ever reached when the PSF being simplified is itself a Piff PSF.
327 """
328 interp = getattr(psf.piff_psf, "interp", None)
329 if interp is None or type(interp).__name__ != "BasisPolynomial" or interp.q is None:
330 return
331 if order >= max(interp._orders):
332 return
334 from piff import BasisPolynomial
336 # ``q`` has one column per active basis term; the terms are the True cells
337 # of ``_mask`` in row-major (i, j) order (see BasisPolynomial.basis). Make
338 # the same ordering for a lower-order interp and copy the shared columns.
339 def _terms(orders: tuple[int, ...], mask: np.ndarray) -> list[tuple[int, ...]]:
340 grids = np.meshgrid(*[np.arange(o + 1) for o in orders], indexing="ij")
341 return list(zip(*(grid[mask].tolist() for grid in grids)))
343 old_terms = _terms(interp._orders, interp._mask)
344 truncated = BasisPolynomial(order, keys=list(interp._keys))
345 new_terms = _terms(truncated._orders, truncated._mask)
346 column_of = {term: index for index, term in enumerate(old_terms)}
347 truncated.q = np.ascontiguousarray(interp.q[:, [column_of[term] for term in new_terms]])
348 psf.piff_psf.interp = truncated
351# -- CellCoadd -------------------------------------------------------------
354def _subset_cell_coadd(
355 cell_coadd: CellCoadd,
356 *,
357 cell_size: int = _CELL_SIZE,
358 kernel_size: int = _KERNEL_SIZE,
359 max_inputs: int = _MAX_INPUTS,
360) -> CellCoadd:
361 """Crop a CellCoadd to a small block of cells and morph it onto a tiny
362 grid (see the module docstring for the rationale).
363 """
364 if kernel_size % 2 == 0:
365 raise ValueError(f"kernel_size must be odd, got {kernel_size}.")
367 # 1. Pick a block of (up to) 2x2 cells, preferring one that contains a
368 # missing cell so the sparse-grid path is exercised. Falls back to the
369 # first available block when the coadd is fully dense.
370 block = cell_coadd[_choose_block_bbox(cell_coadd)]
372 grid = block.grid
373 cs = grid.cell_shape
374 start = block.bounds.grid_start
375 stop = block.bounds.grid_stop
376 n_i = stop.i - start.i
377 n_j = stop.j - start.j
379 # 2. Build a tiny full-patch grid with the same cell *count* as the
380 # original patch but ``cell_size`` pixels per cell, anchored at (0, 0).
381 full_shape = grid.grid_shape
382 new_grid = CellGrid(
383 bbox=Box.factory[0 : full_shape.i * cell_size, 0 : full_shape.j * cell_size],
384 cell_shape=YX(y=cell_size, x=cell_size),
385 )
386 new_block_bbox = _scale_box_to_grid(block.bbox, grid, cell_size)
387 new_bounds = CellGridBounds(grid=new_grid, bbox=new_block_bbox, missing=block.bounds.missing)
389 # 3. Decimate each plane. Because the block's planes tile the kept cells
390 # contiguously, a uniform stride that maps one native cell onto
391 # ``cell_size`` samples is equivalent to per-cell decimation.
392 step_y = max(1, cs.y // cell_size)
393 step_x = max(1, cs.x // cell_size)
394 ny = n_i * cell_size
395 nx = n_j * cell_size
397 def shrink2d(array: np.ndarray) -> np.ndarray:
398 return np.ascontiguousarray(array[::step_y, ::step_x][:ny, :nx])
400 def shrink3d(array: np.ndarray) -> np.ndarray:
401 return np.ascontiguousarray(array[::step_y, ::step_x, :][:ny, :nx, :])
403 # 4. Synthetic-but-valid projection over the tiny tract frame.
404 rng = np.random.default_rng(0)
405 tract_frame = TractFrame(skymap=cell_coadd.skymap, tract=cell_coadd.tract, bbox=new_grid.bbox)
406 projection = make_random_projection(rng, tract_frame, new_block_bbox)
408 unit = cell_coadd.unit
409 image = Image(shrink2d(block.image.array), bbox=new_block_bbox, unit=unit, projection=projection)
410 mask = Mask(shrink3d(block.mask.array), schema=block.mask.schema, bbox=new_block_bbox)
411 variance = Image(shrink2d(block.variance.array), bbox=new_block_bbox, unit=unit**2)
412 mask_fractions = {
413 name: Image(shrink2d(plane.array), bbox=new_block_bbox)
414 for name, plane in block.mask_fractions.items()
415 }
416 noise_realizations = [
417 Image(shrink2d(plane.array), bbox=new_block_bbox) for plane in block.noise_realizations
418 ]
420 # 5. Crop the PSF kernels to a small odd window about their centre,
421 # keeping the (n_i, n_j) per-cell structure and NaN-for-missing cells.
422 psf_array = block.psf._array
423 ky, kx = psf_array.shape[2:]
424 half = kernel_size // 2
425 cy, cx = ky // 2, kx // 2
426 psf_array = np.ascontiguousarray(psf_array[:, :, cy - half : cy + half + 1, cx - half : cx + half + 1])
427 psf = CellPointSpreadFunction(psf_array, bounds=new_bounds)
429 # 6. Patch geometry scaled onto the tiny grid; provenance and backgrounds
430 # are reused as-is (provenance is cell-indexed and already subset).
431 patch = PatchDefinition(
432 id=block.patch.id,
433 index=block.patch.index,
434 inner_bbox=_scale_box_to_grid(block.patch.inner_bbox, grid, cell_size),
435 cells=new_grid,
436 )
438 provenance = block._provenance
439 if provenance is not None:
440 provenance = _trim_provenance(provenance, max_inputs=max_inputs)
442 return CellCoadd(
443 image,
444 mask=mask,
445 variance=variance,
446 mask_fractions=mask_fractions,
447 noise_realizations=noise_realizations,
448 projection=projection,
449 band=block.band,
450 psf=psf,
451 patch=patch,
452 provenance=provenance,
453 backgrounds=block._backgrounds,
454 )
457def _trim_provenance(provenance: CoaddProvenance, *, max_inputs: int) -> CoaddProvenance:
458 """Cap the provenance ``inputs`` table to ``max_inputs`` rows and drop any
459 contributions that reference the removed inputs.
461 The two-table structure, polygon arrays and string dictionary-compression
462 paths are all preserved; only the number of contributing visits shrinks.
463 """
464 inputs = provenance.inputs
465 if len(inputs) <= max_inputs:
466 return provenance
467 kept_inputs = inputs[:max_inputs]
468 keys = {(str(row["instrument"]), int(row["visit"]), int(row["detector"])) for row in kept_inputs}
469 contributions = provenance.contributions
470 mask = np.array(
471 [
472 (str(instrument), int(visit), int(detector)) in keys
473 for instrument, visit, detector in zip(
474 contributions["instrument"], contributions["visit"], contributions["detector"]
475 )
476 ],
477 dtype=bool,
478 )
479 return CoaddProvenance(inputs=kept_inputs, contributions=contributions[mask])
482def _choose_block_bbox(cell_coadd: CellCoadd) -> Box:
483 """Return the pixel bbox of a (up to) 2x2 block of cells to keep.
485 Prefers a block containing a missing cell; otherwise the block anchored at
486 the start of the populated region. Never raises if there is no missing
487 cell.
488 """
489 bounds = cell_coadd.bounds
490 grid = bounds.grid
491 start = bounds.grid_start
492 stop = bounds.grid_stop
493 span_i = min(2, stop.i - start.i)
494 span_j = min(2, stop.j - start.j)
496 target = next(iter(sorted(bounds.missing)), None)
497 if target is not None:
498 # Anchor the block so it includes the missing cell, clamped to the
499 # populated index range.
500 i0 = min(max(target.i, start.i), stop.i - span_i)
501 j0 = min(max(target.j, start.j), stop.j - span_j)
502 else:
503 i0 = start.i
504 j0 = start.j
506 lo = grid.bbox_of(CellIJ(i=i0, j=j0))
507 hi = grid.bbox_of(CellIJ(i=i0 + span_i - 1, j=j0 + span_j - 1))
508 return Box.factory[lo.y.start : hi.y.stop, lo.x.start : hi.x.stop]
511def _scale_box_to_grid(box: Box, grid: CellGrid, cell_size: int) -> Box:
512 """Map a grid-aligned box onto a grid with ``cell_size`` pixels per cell,
513 anchored at the origin.
514 """
515 cs = grid.cell_shape
516 s = grid.bbox.start
517 iy0 = (box.y.start - s.y) // cs.y
518 iy1 = (box.y.stop - s.y) // cs.y
519 ix0 = (box.x.start - s.x) // cs.x
520 ix1 = (box.x.stop - s.x) // cs.x
521 return Box.factory[iy0 * cell_size : iy1 * cell_size, ix0 * cell_size : ix1 * cell_size]