Coverage for python/lsst/images/tests/_minify_for_fixtures.py: 0%
181 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-12 07:51 +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 SkyProjection, TractFrame, Transform
88from .._transforms._ast import PolyMap
89from ..cells import CellCoadd, CellField, CellPointSpreadFunction, CoaddProvenance
90from ..psfs import PiffWrapper
91from ..serialization import backend_for_path, read
92from ._creation import make_random_sky_projection
94# Default morph parameters for CellCoadd. ``CELL_SIZE`` should divide the
95# native cell size evenly; ``KERNEL_SIZE`` must be odd. ``MAX_INPUTS`` caps
96# the provenance ``inputs`` table (a real coadd has hundreds of visits); the
97# full provenance schema is already exercised by the ``coadd_provenance``
98# fixture, so here we keep just enough rows to be representative.
99_CELL_SIZE = 6
100_KERNEL_SIZE = 5
101_MAX_INPUTS = 6
103# Default trim parameters for VisitImage. Amplifiers and aperture-correction
104# entries are homogeneous collections, so a couple of each cover the schema
105# just as well as the full set (a real detector has 16 amplifiers and dozens
106# of aperture corrections).
107_MAX_AMPLIFIERS = 2
108_MAX_APERTURE_CORRECTIONS = 2
110# Field-interpolation order to truncate a Piff PSF to (the solution table of a
111# real order-4 PixelGrid PSF dominates the fixture at ~225 KB). Order 0 is the
112# field-averaged PSF; set to `None` to leave the PSF untouched.
113_PSF_INTERP_ORDER = 0
115# Maximum permitted deviation (radians) when approximating a projection's
116# pixel->sky mapping with an affine one. Over a fixture's tiny box the real
117# mapping is linear well below this, so the fit always succeeds.
118_PROJECTION_LINEAR_APPROX_TOL = 1e-8
121def minify(in_path: str, out_path: str, *, schema_name: str | None = None) -> None:
122 """Read a real archive at ``in_path``, take a small subset, and write JSON.
124 Parameters
125 ----------
126 in_path
127 Path to a FITS (``.fits`` / ``.fits.gz``) or NDF (``.sdf`` / ``.ndf``)
128 file to read.
129 out_path
130 Path to the JSON fixture to write. The parent directory is
131 created if it does not exist.
132 schema_name
133 Top-level schema name (e.g. ``"visit_image"`` or ``"cell_coadd"``).
134 If `None`, it is auto-detected from the file.
136 Raises
137 ------
138 ValueError
139 If the file extension is not recognised.
140 NotImplementedError
141 If the top-level type is not one this helper knows how to subset.
142 """
143 backend = backend_for_path(in_path)
144 if schema_name is None:
145 schema_name = backend.input_archive.get_basic_info(in_path).schema_name
147 cls, subsetter = _dispatch(schema_name)
149 obj: Any = read(in_path, cls)
150 subset = subsetter(obj)
152 tree = images_json.write(subset)
153 dumped = tree.model_dump(mode="json")
154 os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True)
155 with open(out_path, "w") as stream:
156 stream.write(json.dumps(dumped, indent=2, sort_keys=False) + "\n")
159def _dispatch(schema_name: str) -> tuple[type, Callable[[Any], Any]]:
160 """Return the ``(class, subsetter)`` pair for a top-level schema name."""
161 registry: dict[str, tuple[type, Callable[[Any], Any]]] = {
162 "visit_image": (VisitImage, _subset_visit_image),
163 "difference_image": (DifferenceImage, _subset_visit_image),
164 "cell_coadd": (CellCoadd, _subset_cell_coadd),
165 }
166 try:
167 return registry[schema_name]
168 except KeyError:
169 raise NotImplementedError(
170 f"No minify rule for schema {schema_name!r}; supported: {sorted(registry)}."
171 ) from None
174# -- VisitImage ------------------------------------------------------------
177def _subset_visit_image(
178 visit_like_image: VisitImage | DifferenceImage,
179 *,
180 size: int = 16,
181 max_amplifiers: int = _MAX_AMPLIFIERS,
182 max_aperture_corrections: int = _MAX_APERTURE_CORRECTIONS,
183 linearize_projection: bool = True,
184 projection_tol: float = _PROJECTION_LINEAR_APPROX_TOL,
185 psf_interp_order: int | None = _PSF_INTERP_ORDER,
186) -> VisitImage | DifferenceImage:
187 """Crop a VisitImage's pixel planes to a small corner and trim its
188 homogeneous collections.
190 The detector frames are a single structure carried through unchanged by
191 ``__getitem__``. The detector's amplifiers and the aperture-correction map
192 are repeated, schema-identical entries, so they are trimmed to a
193 representative few. The projection's pixel->sky mapping is replaced by its
194 affine approximation over the kept box (see ``_linear_approx_projection``)
195 unless ``linearize_projection`` is false. A Piff PSF's field interpolation
196 is truncated to ``psf_interp_order`` (see ``_simplify_piff_psf``) unless
197 that is `None`.
198 """
199 bbox = visit_like_image.bbox
200 y0 = bbox.y.start
201 x0 = bbox.x.start
202 y1 = min(y0 + size, bbox.y.stop)
203 x1 = min(x0 + size, bbox.x.stop)
204 subset = visit_like_image[Box.factory[y0:y1, x0:x1]]
206 # ``subset`` is a fresh throwaway object whose detector amplifier list,
207 # aperture-correction map and PSF are live, mutable components. Trim them
208 # in place through the public accessors rather than reaching for private
209 # attributes.
210 del subset.detector.amplifiers[max_amplifiers:]
211 aperture_corrections = subset.aperture_corrections
212 for key in list(aperture_corrections)[max_aperture_corrections:]:
213 del aperture_corrections[key]
214 if psf_interp_order is not None and isinstance(subset.psf, PiffWrapper):
215 _simplify_piff_psf(subset.psf, order=psf_interp_order)
217 if not linearize_projection or subset.sky_projection is None:
218 return subset
220 # The pixel planes carry the projection immutably (there is no public
221 # setter for it), so install the affine approximation by rebuilding the
222 # VisitImage from its public components with re-viewed planes. Only the
223 # image plane's projection is actually serialized, but keeping all three
224 # consistent avoids surprises.
225 linear = _linear_approx_projection(subset.sky_projection, subset.image.bbox, tol=projection_tol)
226 return type(visit_like_image)(
227 subset.image.view(sky_projection=linear),
228 mask=subset.mask.view(sky_projection=linear),
229 variance=subset.variance.view(sky_projection=linear),
230 sky_projection=linear,
231 psf=subset.psf,
232 obs_info=subset.obs_info,
233 bounds=subset.bounds,
234 summary_stats=subset.summary_stats,
235 detector=subset.detector,
236 photometric_scaling=subset.photometric_scaling,
237 aperture_corrections=subset.aperture_corrections,
238 backgrounds=subset.backgrounds,
239 band=subset.band,
240 metadata=subset.metadata,
241 )
244def _linear_approx_projection(sky_projection: SkyProjection, bbox: Box, *, tol: float) -> SkyProjection:
245 """Return a copy of ``sky_projection`` whose pixel->sky mapping is replaced
246 by its best linear (affine) approximation over ``bbox``.
248 Real WCS mappings (e.g. TAN-SIP) serialize as large AST polynomial dumps.
249 Over the small box of a fixture they are linear to far below a pixel, so
250 an affine approximation is schema-identical but orders of magnitude
251 smaller. The result carries no FITS approximation (the affine is itself
252 trivially FITS-representable).
254 This is written as a self-contained ``sky_projection -> sky_projection``
255 transform so it can be promoted to a public
256 ``SkyProjection.linear_approx(bbox, tol)`` method later with essentially
257 no change. It assumes a 2-D pixel->sky
258 mapping.
260 Parameters
261 ----------
262 sky_projection
263 The projection to approximate.
264 bbox
265 Box (in pixel coordinates) over which the approximation must hold.
266 tol
267 Maximum permitted deviation from linearity, as a Cartesian
268 displacement in the output (sky, radians) coordinates. AST raises
269 ``RuntimeError`` if no fit within ``tol`` exists.
270 """
271 transform = sky_projection.pixel_to_sky_transform
272 mapping = transform._ast_mapping
273 lbnd = [bbox.x.start, bbox.y.start]
274 ubnd = [bbox.x.stop, bbox.y.stop]
275 # linearApprox yields [offsets; Jacobian] as a (1 + n_out, n_in) array on
276 # both AST backends (astshim returns the flat buffer in the same order, so
277 # the reshape recovers the same layout the starlink-pyast bridge returns).
278 fit = np.asarray(mapping.linearApprox(lbnd, ubnd, tol), dtype=float).reshape(3, 2)
279 offset = fit[0] # (lon0, lat0), radians
280 jacobian = fit[1:] # jacobian[i, j] = d(out_i) / d(in_j), in = (x, y)
281 jacobian_inv = np.linalg.inv(jacobian)
282 forward = _affine_polymap_coeffs(jacobian, offset)
283 inverse = _affine_polymap_coeffs(jacobian_inv, -jacobian_inv @ offset)
284 affine = Transform(
285 transform.in_frame,
286 transform.out_frame,
287 PolyMap(forward, inverse),
288 in_bounds=sky_projection.pixel_bounds,
289 )
290 return SkyProjection(affine)
293def _affine_polymap_coeffs(matrix: np.ndarray, offset: np.ndarray) -> np.ndarray:
294 """Build AST ``PolyMap`` coefficients for ``out = matrix @ in + offset``.
296 Each row is ``[coefficient, output_axis (1-based), power_of_in_1, ...]``;
297 one constant row plus one row per input per output axis. Returned as a
298 float array, which is the form both AST backends require.
299 """
300 n = len(offset)
301 coeffs: list[list[float]] = []
302 for i in range(n):
303 coeffs.append([float(offset[i]), i + 1, *([0] * n)])
304 for j in range(n):
305 powers = [1 if k == j else 0 for k in range(n)]
306 coeffs.append([float(matrix[i][j]), i + 1, *powers])
307 return np.array(coeffs, dtype=float)
310def _simplify_piff_psf(psf: PiffWrapper, *, order: int) -> None:
311 """Truncate a Piff PSF's field interpolation to ``order``, in place.
313 A real Piff PSF interpolates a per-pixel model across the focal plane with
314 a high-order 2-D polynomial; that solution table dominates the serialized
315 size (a 25x25 PixelGrid x order-4 polynomial is ~225 KB). Truncating to
316 ``order`` keeps only the lowest-order field terms -- order 0 is the
317 field-averaged PSF -- which is schema-identical but far smaller, and needs
318 no stars or refit (the fitted ``stars`` are already dropped on serialize).
320 Only ``BasisPolynomial``-interpolated PSFs are handled; anything else (a
321 higher-order model already at/under ``order``, a non-polynomial interp) is
322 left untouched.
324 ``piff`` is imported lazily because it is an optional dependency; this is
325 only ever reached when the PSF being simplified is itself a Piff PSF.
326 """
327 interp = getattr(psf.piff_psf, "interp", None)
328 if interp is None or type(interp).__name__ != "BasisPolynomial" or interp.q is None:
329 return
330 if order >= max(interp._orders):
331 return
333 from piff import BasisPolynomial
335 # ``q`` has one column per active basis term; the terms are the True cells
336 # of ``_mask`` in row-major (i, j) order (see BasisPolynomial.basis). Make
337 # the same ordering for a lower-order interp and copy the shared columns.
338 def _terms(orders: tuple[int, ...], mask: np.ndarray) -> list[tuple[int, ...]]:
339 grids = np.meshgrid(*[np.arange(o + 1) for o in orders], indexing="ij")
340 return list(zip(*(grid[mask].tolist() for grid in grids)))
342 old_terms = _terms(interp._orders, interp._mask)
343 truncated = BasisPolynomial(order, keys=list(interp._keys))
344 new_terms = _terms(truncated._orders, truncated._mask)
345 column_of = {term: index for index, term in enumerate(old_terms)}
346 truncated.q = np.ascontiguousarray(interp.q[:, [column_of[term] for term in new_terms]])
347 psf.piff_psf.interp = truncated
350# -- CellCoadd -------------------------------------------------------------
353def _subset_cell_coadd(
354 cell_coadd: CellCoadd,
355 *,
356 cell_size: int = _CELL_SIZE,
357 kernel_size: int = _KERNEL_SIZE,
358 max_inputs: int = _MAX_INPUTS,
359) -> CellCoadd:
360 """Crop a CellCoadd to a small block of cells and morph it onto a tiny
361 grid (see the module docstring for the rationale).
362 """
363 if kernel_size % 2 == 0:
364 raise ValueError(f"kernel_size must be odd, got {kernel_size}.")
366 # 1. Pick a block of (up to) 2x2 cells, preferring one that contains a
367 # missing cell so the sparse-grid path is exercised. Falls back to the
368 # first available block when the coadd is fully dense.
369 block = cell_coadd[_choose_block_bbox(cell_coadd)]
371 grid = block.grid
372 cs = grid.cell_shape
373 start = block.bounds.subgrid_start
374 stop = block.bounds.subgrid_stop
375 n_i = stop.i - start.i
376 n_j = stop.j - start.j
378 # 2. Build a tiny full-patch grid with the same cell *count* as the
379 # original patch but ``cell_size`` pixels per cell, anchored at (0, 0).
380 full_shape = grid.grid_size
381 new_grid = CellGrid(
382 bbox=Box.factory[0 : full_shape.i * cell_size, 0 : full_shape.j * cell_size],
383 cell_shape=YX(y=cell_size, x=cell_size),
384 )
385 new_block_bbox = _scale_box_to_grid(block.bbox, grid, cell_size)
386 new_bounds = CellGridBounds(grid=new_grid, bbox=new_block_bbox, missing=block.bounds.missing)
388 # 3. Decimate each plane. Because the block's planes tile the kept cells
389 # contiguously, a uniform stride that maps one native cell onto
390 # ``cell_size`` samples is equivalent to per-cell decimation.
391 step_y = max(1, cs.y // cell_size)
392 step_x = max(1, cs.x // cell_size)
393 ny = n_i * cell_size
394 nx = n_j * cell_size
396 def shrink2d(array: np.ndarray) -> np.ndarray:
397 return np.ascontiguousarray(array[::step_y, ::step_x][:ny, :nx])
399 def shrink3d(array: np.ndarray) -> np.ndarray:
400 return np.ascontiguousarray(array[::step_y, ::step_x, :][:ny, :nx, :])
402 # 4. Synthetic-but-valid sky_projection over the tiny tract frame.
403 rng = np.random.default_rng(0)
404 tract_frame = TractFrame(skymap=cell_coadd.skymap, tract=cell_coadd.tract, bbox=new_grid.bbox)
405 sky_projection = make_random_sky_projection(rng, tract_frame, new_block_bbox)
407 unit = cell_coadd.unit
408 image = Image(shrink2d(block.image.array), bbox=new_block_bbox, unit=unit, sky_projection=sky_projection)
409 mask = Mask(shrink3d(block.mask.array), schema=block.mask.schema, bbox=new_block_bbox)
410 variance = Image(shrink2d(block.variance.array), bbox=new_block_bbox, unit=unit**2)
411 mask_fractions = {
412 name: Image(shrink2d(plane.array), bbox=new_block_bbox)
413 for name, plane in block.mask_fractions.items()
414 }
415 noise_realizations = [
416 Image(shrink2d(plane.array), bbox=new_block_bbox) for plane in block.noise_realizations
417 ]
419 # 5. Crop the PSF kernels to a small odd window about their centre,
420 # keeping the (n_i, n_j) per-cell structure and NaN-for-missing cells.
421 psf_array = block.psf._array
422 ky, kx = psf_array.shape[2:]
423 half = kernel_size // 2
424 cy, cx = ky // 2, kx // 2
425 psf_array = np.ascontiguousarray(psf_array[:, :, cy - half : cy + half + 1, cx - half : cx + half + 1])
426 psf = CellPointSpreadFunction(psf_array, bounds=new_bounds)
428 # 6. Patch geometry scaled onto the tiny grid; provenance and backgrounds
429 # are reused as-is (provenance is cell-indexed and already subset).
430 patch = PatchDefinition(
431 id=block.patch.id,
432 index=block.patch.index,
433 inner_bbox=_scale_box_to_grid(block.patch.inner_bbox, grid, cell_size),
434 cells=new_grid,
435 )
437 provenance = block._provenance
438 if provenance is not None:
439 provenance = _trim_provenance(provenance, max_inputs=max_inputs)
441 # Aperture corrections are not subset when CellCoadd is subset with a
442 # bounding box, because they're always tiny. But that makes setting
443 # up a consistent new grid for them tricky.
444 aperture_corrections = {}
445 new_apcorr_bounds = None
446 for i, (name, field) in enumerate(block.aperture_corrections.items()):
447 if new_apcorr_bounds is None:
448 new_apcorr_bounds = CellGridBounds(
449 grid=new_grid,
450 bbox=_scale_box_to_grid(field.bounds.bbox, grid, cell_size),
451 missing=cell_coadd.bounds.missing,
452 )
453 aperture_corrections[name] = CellField(new_apcorr_bounds, field._array)
454 if i >= 2:
455 break
457 return CellCoadd(
458 image,
459 mask=mask,
460 variance=variance,
461 mask_fractions=mask_fractions,
462 noise_realizations=noise_realizations,
463 sky_projection=sky_projection,
464 band=block.band,
465 psf=psf,
466 patch=patch,
467 provenance=provenance,
468 backgrounds=block._backgrounds,
469 aperture_corrections=aperture_corrections,
470 )
473def _trim_provenance(provenance: CoaddProvenance, *, max_inputs: int) -> CoaddProvenance:
474 """Cap the provenance ``inputs`` table to ``max_inputs`` rows and drop any
475 contributions that reference the removed inputs.
477 The two-table structure, polygon arrays and string dictionary-compression
478 paths are all preserved; only the number of contributing visits shrinks.
479 """
480 inputs = provenance.inputs
481 if len(inputs) <= max_inputs:
482 return provenance
483 kept_inputs = inputs[:max_inputs]
484 keys = {(str(row["instrument"]), int(row["visit"]), int(row["detector"])) for row in kept_inputs}
485 contributions = provenance.contributions
486 mask = np.array(
487 [
488 (str(instrument), int(visit), int(detector)) in keys
489 for instrument, visit, detector in zip(
490 contributions["instrument"], contributions["visit"], contributions["detector"]
491 )
492 ],
493 dtype=bool,
494 )
495 return CoaddProvenance(inputs=kept_inputs, contributions=contributions[mask])
498def _choose_block_bbox(cell_coadd: CellCoadd) -> Box:
499 """Return the pixel bbox of a (up to) 2x2 block of cells to keep.
501 Prefers a block containing a missing cell; otherwise the block anchored at
502 the start of the populated region. Never raises if there is no missing
503 cell.
504 """
505 bounds = cell_coadd.bounds
506 grid = bounds.grid
507 start = bounds.subgrid_start
508 stop = bounds.subgrid_stop
509 span_i = min(2, stop.i - start.i)
510 span_j = min(2, stop.j - start.j)
512 target = next(iter(sorted(bounds.missing)), None)
513 if target is not None:
514 # Anchor the block so it includes the missing cell, clamped to the
515 # populated index range.
516 i0 = min(max(target.i, start.i), stop.i - span_i)
517 j0 = min(max(target.j, start.j), stop.j - span_j)
518 else:
519 i0 = start.i
520 j0 = start.j
522 lo = grid.bbox_of(CellIJ(i=i0, j=j0))
523 hi = grid.bbox_of(CellIJ(i=i0 + span_i - 1, j=j0 + span_j - 1))
524 return Box.factory[lo.y.start : hi.y.stop, lo.x.start : hi.x.stop]
527def _scale_box_to_grid(box: Box, grid: CellGrid, cell_size: int) -> Box:
528 """Map a grid-aligned box onto a grid with ``cell_size`` pixels per cell,
529 anchored at the origin.
530 """
531 cs = grid.cell_shape
532 s = grid.bbox.start
533 iy0 = (box.y.start - s.y) // cs.y
534 iy1 = (box.y.stop - s.y) // cs.y
535 ix0 = (box.x.start - s.x) // cs.x
536 ix1 = (box.x.stop - s.x) // cs.x
537 return Box.factory[iy0 * cell_size : iy1 * cell_size, ix0 * cell_size : ix1 * cell_size]