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