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