Coverage for python/lsst/images/_mask.py: 87%
468 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-24 01:50 -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.
12from __future__ import annotations
14__all__ = (
15 "Mask",
16 "MaskPlane",
17 "MaskPlaneBit",
18 "MaskSchema",
19 "MaskSerializationModel",
20 "get_legacy_deep_coadd_mask_planes",
21 "get_legacy_difference_image_mask_planes",
22 "get_legacy_non_cell_coadd_mask_planes",
23 "get_legacy_visit_image_mask_planes",
24)
26import dataclasses
27import math
28from collections.abc import Callable, Iterable, Iterator, Mapping, Sequence, Set
29from types import EllipsisType
30from typing import TYPE_CHECKING, Any, ClassVar, cast
32import astropy.io.fits
33import astropy.wcs
34import numpy as np
35import numpy.typing as npt
36import pydantic
38from lsst.resources import ResourcePath, ResourcePathExpression
40from . import fits
41from ._generalized_image import GeneralizedImage
42from ._geom import YX, Box, NoOverlapError
43from ._transforms import Frame, SkyProjection, SkyProjectionSerializationModel
44from .serialization import (
45 ArchiveReadError,
46 ArchiveTree,
47 ArrayReferenceModel,
48 InlineArrayModel,
49 InputArchive,
50 IntegerType,
51 InvalidParameterError,
52 MetadataValue,
53 NumberType,
54 OutputArchive,
55 is_integer,
56 no_header_updates,
57)
58from .utils import is_none
60if TYPE_CHECKING:
61 try:
62 from lsst.afw.image import Mask as LegacyMask
63 except ImportError:
64 type LegacyMask = Any # type: ignore[no-redef]
67@dataclasses.dataclass(frozen=True)
68class MaskPlane:
69 """Name and description of a single plane in a mask array."""
71 name: str
72 """Unique name for the mask plane (`str`)."""
74 description: str
75 """Human-readable documentation for the mask plane (`str`)."""
77 @classmethod
78 def read_legacy(cls, header: astropy.io.fits.Header, *, strip: bool = True) -> dict[str, int]:
79 """Read mask plane descriptions written by
80 `lsst.afw.image.Mask.writeFits`.
82 Parameters
83 ----------
84 header
85 FITS header.
86 strip
87 If `True` (default), delete the ``MP_`` cards from the header after
88 reading them, as appropriate when the mask is being reinterpreted
89 for new code only. If `False`, leave them in place so they can be
90 propagated for backwards compatibility (re-indexed to the new
91 schema by the caller).
93 Returns
94 -------
95 `dict` [`str`, `int`]
96 A dictionary mapping mask plane name to integer bit index.
97 """
98 result: dict[str, int] = {}
99 for card in list(header.cards):
100 if card.keyword.startswith("MP_"):
101 result[card.keyword.removeprefix("MP_")] = card.value
102 if strip:
103 del header[card.keyword]
104 return result
107@dataclasses.dataclass(frozen=True)
108class MaskPlaneBit:
109 """The nested array index and mask value associated with a single mask
110 plane.
111 """
113 index: int
114 """Index into the last dimension of the mask array where this plane's bit
115 is stored.
116 """
118 mask: np.integer
119 """Bitmask that selects just this plane's bit from a mask array value
120 (`numpy.integer`).
121 """
123 @classmethod
124 def compute(cls, overall_index: int, stride: int, mask_type: type[np.integer]) -> MaskPlaneBit:
125 """Construct a `MaskPlaneBit` from the overall index of a plane in a
126 `MaskSchema` and the stride (number of bits per mask array element).
127 """
128 index, bit = divmod(overall_index, stride)
129 return cls(index, mask_type(1 << bit))
132class MaskSchema:
133 """A schema for a bit-packed mask array.
135 Parameters
136 ----------
137 planes
138 Iterable of `MaskPlane` instances that define the schema. `None`
139 values may be included to reserve bits for future use.
140 dtype
141 The numpy data type of the mask arrays that use this schema.
143 Notes
144 -----
145 A `MaskSchema` is a collection of mask planes, which each correspond to a
146 single bit in a mask array. Mask schemas are immutable and associated with
147 a particular array data type, allowing them to safely precompute the index
148 and bitmask for each plane.
150 `MaskSchema` indexing is by integer (the overall index of a plane in the
151 schema). The `descriptions` attribute may be indexed by plane name to get
152 the description for that plane, and the `bitmask` method can be used to
153 obtain an array that can be used to select one or more planes by name in
154 a mask array that uses this schema.
156 If no mask planes are provided, a `None` placeholder is automatically
157 added.
158 """
160 def __init__(self, planes: Iterable[MaskPlane | None], dtype: npt.DTypeLike = np.uint8) -> None:
161 self._planes: tuple[MaskPlane | None, ...] = tuple(planes) or (None,)
162 self._dtype = cast(np.dtype[np.integer], np.dtype(dtype))
163 stride = self.bits_per_element(self._dtype)
164 self._descriptions = {plane.name: plane.description for plane in self._planes if plane is not None}
165 self._mask_size = math.ceil(len(self._planes) / stride)
166 self._bits: dict[str, MaskPlaneBit] = {
167 plane.name: MaskPlaneBit.compute(n, stride, self._dtype.type)
168 for n, plane in enumerate(self._planes)
169 if plane is not None
170 }
172 @staticmethod
173 def bits_per_element(dtype: npt.DTypeLike) -> int:
174 """Return the number of mask bits per array element for the given
175 data type.
176 """
177 dtype = np.dtype(dtype)
178 match dtype.kind:
179 case "u":
180 return dtype.itemsize * 8
181 case "i":
182 return dtype.itemsize * 8 - 1
183 case _:
184 raise TypeError(f"dtype for masks must be an integer; got {dtype} with kind={dtype.kind}.")
186 def __iter__(self) -> Iterator[MaskPlane | None]:
187 return iter(self._planes)
189 def __len__(self) -> int:
190 return len(self._planes)
192 def __contains__(self, plane: str | MaskPlane) -> bool:
193 return getattr(plane, "name", plane) in self.names
195 def __getitem__(self, i: int) -> MaskPlane | None:
196 return self._planes[i]
198 def __repr__(self) -> str:
199 return f"MaskSchema({list(self._planes)}, dtype={self._dtype!r})"
201 def __str__(self) -> str:
202 return "\n".join(
203 [
204 f"{name} [{bit.index}@{hex(bit.mask)}]: {self._descriptions[name]}"
205 for name, bit in self._bits.items()
206 ]
207 )
209 def __eq__(self, other: object) -> bool:
210 if isinstance(other, MaskSchema): 210 ↛ 212line 210 didn't jump to line 212 because the condition on line 210 was always true
211 return self._planes == other._planes and self._dtype == other._dtype
212 return False
214 @property
215 def dtype(self) -> np.dtype:
216 """The numpy data type of the mask arrays that use this schema."""
217 return self._dtype
219 @property
220 def mask_size(self) -> int:
221 """The number of elements in the last dimension of any mask array that
222 uses this schema.
223 """
224 return self._mask_size
226 @property
227 def names(self) -> Set[str]:
228 """The names of the mask planes, in bit order."""
229 return self._bits.keys()
231 @property
232 def descriptions(self) -> Mapping[str, str]:
233 """A mapping from plane name to description."""
234 return self._descriptions
236 def bit(self, plane: str) -> MaskPlaneBit:
237 """Return the last array index and mask for the given mask plane."""
238 return self._bits[plane]
240 def bitmask(self, *planes: str) -> np.ndarray:
241 """Return a 1-d mask array that represents the union (i.e. bitwise OR)
242 of the planes with the given names.
244 Parameters
245 ----------
246 *planes
247 Mask plane names.
249 Returns
250 -------
251 numpy.ndarray
252 A 1-d array with shape ``(mask_size,)``.
253 """
254 result = np.zeros(self.mask_size, dtype=self._dtype)
255 for plane in planes:
256 bit = self._bits[plane]
257 result[bit.index] |= bit.mask
258 return result
260 def split(self, dtype: npt.DTypeLike) -> list[MaskSchema]:
261 """Split the schema into an equivalent series of schemas that each
262 have a `mask_size` of ``1``, dropping all `None` placeholders.
264 Parameters
265 ----------
266 dtype
267 Data type of the new mask pixels.
269 Returns
270 -------
271 `list` [`MaskSchema`]
272 A list of mask schemas that together include all planes in
273 ``self`` and have `mask_size` equal to ``1``. If there are no
274 mask planes (only `None` placeholders) in ``self``, a single mask
275 schema with a `None` placeholder is returned; otherwise `None`
276 placeholders are returned.
277 """
278 dtype = np.dtype(dtype)
279 planes: list[MaskPlane] = []
280 schemas: list[MaskSchema] = []
281 n_planes_per_schema = self.bits_per_element(dtype)
282 for plane in self._planes:
283 if plane is not None:
284 planes.append(plane)
285 if len(planes) == n_planes_per_schema:
286 schemas.append(MaskSchema(planes, dtype=dtype))
287 planes.clear()
288 if planes: 288 ↛ 290line 288 didn't jump to line 290 because the condition on line 288 was always true
289 schemas.append(MaskSchema(planes, dtype=dtype))
290 if not schemas: 290 ↛ 291line 290 didn't jump to line 291 because the condition on line 290 was never true
291 schemas.append(MaskSchema([None], dtype=dtype))
292 return schemas
294 def update_header(self, header: astropy.io.fits.Header) -> None:
295 """Add a description of this mask schema to a FITS header."""
296 for n, plane in enumerate(self):
297 if plane is not None:
298 bit = self.bit(plane.name)
299 if bit.index != 0: 299 ↛ 300line 299 didn't jump to line 300 because the condition on line 299 was never true
300 raise TypeError("Only mask schemas with mask_size==1 can be described in FITS.")
301 header.set(f"MSKN{n:04d}", plane.name, f"Name for mask plane {n}.")
302 header.set(f"MSKM{n:04d}", bit.mask, f"Bitmask for plane n={n}; always 1<<n.")
303 # We don't add a comment to the description card, because it's
304 # likely to overrun a single card and get the CONTINUE
305 # treatment. That will cause Astropy to warn about the comment
306 # being truncated and that's worse than just leaving it
307 # unexplained; it's pretty obvious from context what it is.
308 header.set(f"MSKD{n:04d}", plane.description)
310 def strip_header(self, header: astropy.io.fits.Header) -> None:
311 """Remove all header cards added by `update_header`."""
312 for n, plane in enumerate(self):
313 if plane is not None: 313 ↛ 312line 313 didn't jump to line 312 because the condition on line 313 was always true
314 header.remove(f"MSKN{n:04d}", ignore_missing=True)
315 header.remove(f"MSKM{n:04d}", ignore_missing=True)
316 header.remove(f"MSKD{n:04d}", ignore_missing=True)
318 @classmethod
319 def from_fits_header(cls, header: astropy.io.fits.Header, dtype: npt.DTypeLike = np.uint8) -> MaskSchema:
320 """Reconstruct a schema from the ``MSKN``/``MSKD`` cards written by
321 `update_header`.
323 Parameters
324 ----------
325 header
326 FITS header containing ``MSKN{n:04d}`` plane-name cards and
327 ``MSKD{n:04d}`` description cards.
328 dtype
329 Data type of the mask arrays that will use this schema. The cards
330 describe a ``mask_size==1`` serialized form and do not record the
331 in-memory dtype, so the caller must supply it; it defaults to the
332 same ``uint8`` used by the `Mask` constructor.
334 Returns
335 -------
336 `MaskSchema`
337 Schema whose planes are ordered by their ``MSKN`` index, with
338 `None` placeholders inserted for any gaps in that numbering.
340 Raises
341 ------
342 ValueError
343 Raised if the header contains no ``MSKN`` cards.
344 """
345 planes_by_index: dict[int, MaskPlane] = {}
346 for card in header.cards:
347 if card.keyword.startswith("MSKN"):
348 n = int(card.keyword.removeprefix("MSKN"))
349 planes_by_index[n] = MaskPlane(card.value, header.get(f"MSKD{n:04d}", ""))
350 if not planes_by_index:
351 raise ValueError("Header has no MSKN cards describing a mask schema.")
352 planes = [planes_by_index.get(n) for n in range(max(planes_by_index) + 1)]
353 return cls(planes, dtype=dtype)
356class Mask(GeneralizedImage):
357 """A 2-d bitmask image backed by a 3-d byte array.
359 Parameters
360 ----------
361 array_or_fill
362 Array or fill value for the mask. If a fill value, ``bbox`` or
363 ``shape`` must be provided.
364 schema
365 Schema that defines the planes and their bit assignments.
366 bbox
367 Bounding box for the mask. This sets the shape of the first two
368 dimensions of the array.
369 yx0
370 Logical coordinates of the first pixel in the array, ordered ``y``,
371 ``x`` (unless an `XY` instance is passed). Ignored if
372 ``bbox`` is provided. Defaults to zeros.
373 shape
374 Leading dimensions of the array, ordered ``y``, ``x`` (unless an `XY`
375 instance is passed). Only needed if ``array_or_fill`` is not an
376 array and ``bbox`` is not provided. Like the bbox, this does not
377 include the last dimension of the array.
378 sky_projection
379 Projection that maps the pixel grid to the sky.
380 metadata
381 Arbitrary flexible metadata to associate with the mask.
383 Notes
384 -----
385 Indexing the `array` attribute of a `Mask` does not take into account its
386 ``yx0`` offset, but accessing a subimage mask by indexing a `Mask` with
387 a `Box` does, and the `bbox` of the subimage is set to match its location
388 within the original mask.
390 A mask's ``bbox`` corresponds to the leading dimensions of its backing
391 `numpy.ndarray`, while the last dimension's size is always equal to the
392 `~MaskSchema.mask_size` of its schema, since a schema can in general
393 require multiple array elements to represent all of its planes.
394 """
396 def __init__(
397 self,
398 array_or_fill: np.ndarray | int = 0,
399 /,
400 *,
401 schema: MaskSchema,
402 bbox: Box | None = None,
403 yx0: Sequence[int] | None = None,
404 shape: Sequence[int] | None = None,
405 sky_projection: SkyProjection | None = None,
406 metadata: dict[str, MetadataValue] | None = None,
407 ) -> None:
408 super().__init__(metadata)
409 if shape is not None:
410 shape = tuple(shape)
411 if isinstance(array_or_fill, np.ndarray):
412 array = np.array(array_or_fill, dtype=schema.dtype, copy=None)
413 if array.ndim != 3:
414 raise ValueError("Mask array must be 3-d.")
415 if bbox is None:
416 bbox = Box.from_shape(array.shape[:-1], start=yx0)
417 elif bbox.shape + (schema.mask_size,) != array.shape:
418 raise ValueError(
419 f"Explicit bbox shape {bbox.shape} and schema of size {schema.mask_size} do not "
420 f"match array with shape {array.shape}."
421 )
422 if shape is not None and shape + (schema.mask_size,) != array.shape:
423 raise ValueError(
424 f"Explicit shape {shape} and schema of size {schema.mask_size} do "
425 f"not match array with shape {array.shape}."
426 )
428 else:
429 if bbox is None:
430 if shape is None:
431 raise TypeError("No bbox, size, or array provided.")
432 bbox = Box.from_shape(shape, start=yx0)
433 array = np.full(bbox.shape + (schema.mask_size,), array_or_fill, dtype=schema.dtype)
434 self._array = array
435 self._bbox: Box = bbox
436 self._schema: MaskSchema = schema
437 self._sky_projection = sky_projection
439 @property
440 def array(self) -> np.ndarray:
441 """The low-level array (`numpy.ndarray`).
443 Assigning to this attribute modifies the existing array in place; the
444 bounding box and underlying data pointer are never changed.
445 """
446 return self._array
448 @array.setter
449 def array(self, value: np.ndarray | int) -> None:
450 self._array[:, :] = value
452 @property
453 def schema(self) -> MaskSchema:
454 """Schema that defines the planes and their bit assignments
455 (`MaskSchema`).
456 """
457 return self._schema
459 @property
460 def bbox(self) -> Box:
461 """2-d bounding box of the mask (`Box`).
463 This sets the shape of the first two dimensions of the array.
464 """
465 return self._bbox
467 @property
468 def sky_projection(self) -> SkyProjection[Any] | None:
469 """The projection that maps this mask's pixel grid to the sky
470 (`SkyProjection` | `None`).
472 Notes
473 -----
474 The pixel coordinates used by this projection account for the bounding
475 box ``start`` (i.e. ``yx0``); they are not just array indices.
476 """
477 return self._sky_projection
479 def __getitem__(self, bbox: Box | EllipsisType) -> Mask:
480 if bbox is ...:
481 return self
482 super().__getitem__(bbox)
483 return self._transfer_metadata(
484 Mask(
485 self.array[bbox.y.slice_within(self._bbox.y), bbox.x.slice_within(self._bbox.x), :],
486 bbox=bbox,
487 schema=self.schema,
488 sky_projection=self._sky_projection,
489 ),
490 bbox=bbox,
491 )
493 def __setitem__(self, bbox: Box | EllipsisType, value: Mask) -> None:
494 subview = self[bbox]
495 subview.clear()
496 subview.update(value)
498 def __str__(self) -> str:
499 return f"Mask({self.bbox!s}, {list(self.schema.names)})"
501 def __repr__(self) -> str:
502 return f"Mask(..., bbox={self.bbox!r}, schema={self.schema!r})"
504 def __eq__(self, other: object) -> bool:
505 if not isinstance(other, Mask):
506 return NotImplemented
507 return (
508 self._bbox == other._bbox
509 and self._schema == other._schema
510 and np.array_equal(self._array, other._array, equal_nan=True)
511 )
513 def copy(self) -> Mask:
514 """Deep-copy the mask and metadata."""
515 return self._transfer_metadata(
516 Mask(
517 self._array.copy(), bbox=self._bbox, schema=self._schema, sky_projection=self._sky_projection
518 ),
519 copy=True,
520 )
522 def view(
523 self,
524 *,
525 schema: MaskSchema | EllipsisType = ...,
526 sky_projection: SkyProjection | None | EllipsisType = ...,
527 yx0: Sequence[int] | EllipsisType = ...,
528 ) -> Mask:
529 """Make a view of the mask, with optional updates.
531 Notes
532 -----
533 This can only be used to make changes to schema descriptions; plane
534 names must remain the same (in the same order).
535 """
536 if schema is ...: 536 ↛ 539line 536 didn't jump to line 539 because the condition on line 536 was always true
537 schema = self._schema
538 else:
539 if list(schema.names) != list(self.schema.names):
540 raise ValueError("Cannot create a mask view with a schema with different names.")
541 if sky_projection is ...: 541 ↛ 542line 541 didn't jump to line 542 because the condition on line 541 was never true
542 sky_projection = self._sky_projection
543 if yx0 is ...: 543 ↛ 545line 543 didn't jump to line 545 because the condition on line 543 was always true
544 yx0 = self._bbox.start
545 return self._transfer_metadata(
546 Mask(self._array, yx0=yx0, schema=schema, sky_projection=sky_projection)
547 )
549 def update(self, other: Mask) -> None:
550 """Update ``self`` to include all common mask values set in ``other``.
552 Notes
553 -----
554 This only operates on the intersection of the two mask bounding boxes
555 and the mask planes that are present in both. Mask bits are only set,
556 not cleared (i.e. this uses ``|=`` updates, not ``=`` assignments).
557 """
558 lhs = self
559 rhs = other
560 if other.bbox != self.bbox: 560 ↛ 561line 560 didn't jump to line 561 because the condition on line 560 was never true
561 try:
562 bbox = self.bbox.intersection(other.bbox)
563 except NoOverlapError:
564 return
565 lhs = self[bbox]
566 rhs = other[bbox]
567 for name in self.schema.names & other.schema.names:
568 lhs.set(name, rhs.get(name))
570 def get(self, plane: str) -> np.ndarray:
571 """Return a 2-d boolean array for the given mask plane.
573 Parameters
574 ----------
575 plane
576 Name of the mask plane.
578 Returns
579 -------
580 numpy.ndarray
581 A 2-d boolean array with the same shape as `bbox` that is `True`
582 where the bit for ``plane`` is set and `False` elsewhere.
583 """
584 bit = self.schema.bit(plane)
585 return (self._array[..., bit.index] & bit.mask).astype(bool)
587 def set(self, plane: str, boolean_mask: np.ndarray | EllipsisType = ...) -> None:
588 """Set a mask plane.
590 Parameters
591 ----------
592 plane
593 Name of the mask plane to set
594 boolean_mask
595 A 2-d boolean array with the same shape as `bbox` that is `True`
596 where the bit for ``plane`` should be set and `False` where it
597 should be left unchanged (*not* set to zero). May be ``...`` to
598 set the bit everywhere.
599 """
600 bit = self.schema.bit(plane)
601 if boolean_mask is not ...: 601 ↛ 603line 601 didn't jump to line 603 because the condition on line 601 was always true
602 boolean_mask = boolean_mask.astype(bool)
603 self._array[boolean_mask, bit.index] |= bit.mask
605 def clear(self, plane: str | None = None, boolean_mask: np.ndarray | EllipsisType = ...) -> None:
606 """Clear one or more mask planes.
608 Parameters
609 ----------
610 plane
611 Name of the mask plane to set. If `None` all mask planes are
612 cleared.
613 boolean_mask
614 A 2-d boolean array with the same shape as `bbox` that is `True`
615 where the bit for ``plane`` should be cleared and `False` where it
616 should be left unchanged. May be ``...`` to clear the bit
617 everywhere.
618 """
619 if boolean_mask is not ...: 619 ↛ 620line 619 didn't jump to line 620 because the condition on line 619 was never true
620 boolean_mask = boolean_mask.astype(bool)
621 if plane is None: 621 ↛ 624line 621 didn't jump to line 624 because the condition on line 621 was always true
622 self._array[boolean_mask, :] = 0
623 else:
624 bit = self.schema.bit(plane)
625 self._array[boolean_mask, bit.index] &= ~bit.mask
627 def add_plane(self, name: str, description: str) -> Mask:
628 """Return a new mask with one additional mask plane.
630 This is a convenience wrapper around `add_planes` for the common case
631 of adding a single plane.
633 Parameters
634 ----------
635 name
636 Unique name for the new mask plane.
637 description
638 Human-readable documentation for the new mask plane.
640 Returns
641 -------
642 `Mask`
643 A new mask whose schema includes the new plane; see `add_planes`
644 for the reallocation and view semantics.
646 Raises
647 ------
648 ValueError
649 Raised if a plane named ``name`` already exists.
650 """
651 return self.add_planes([MaskPlane(name, description)])
653 def add_planes(self, planes: Iterable[MaskPlane | None], *, drop: Iterable[str] = ()) -> Mask:
654 """Return a new mask with planes added and/or dropped.
656 Parameters
657 ----------
658 planes
659 New mask planes to append, in order, after the planes retained
660 from this mask. `None` entries reserve unused bits (placeholders),
661 exactly as in `MaskSchema`.
662 drop
663 Names of existing planes to remove from the schema.
665 Returns
666 -------
667 `Mask`
668 A new mask with the updated schema. Retained planes keep their
669 pixel values (copied by name); newly added planes start cleared.
671 Raises
672 ------
673 ValueError
674 Raised if a name in ``drop`` is not an existing plane, or if a
675 plane in ``planes`` collides with a retained plane name.
677 Notes
678 -----
679 Adding or dropping planes always reallocates the backing array and
680 returns a new `Mask`; this mask is left unchanged and any views or
681 subimages of it continue to refer to the original array with the
682 original schema. This is deliberate: there is no way to update the
683 schema of an existing view, and a stale view must never set bits that
684 its now-outdated schema regards as unused. Dropping a plane compacts
685 the schema, so planes after it are reassigned to lower bits and the
686 pixel values are repacked by plane name to match.
687 """
688 drop_set = set(drop)
689 if unknown := drop_set - set(self._schema.names):
690 raise ValueError(f"Cannot drop mask planes that do not exist: {sorted(unknown)}.")
691 retained = [plane for plane in self._schema if plane is None or plane.name not in drop_set]
692 names = {plane.name for plane in retained if plane is not None}
693 new_planes = list(planes)
694 for plane in new_planes:
695 if plane is None:
696 continue
697 if plane.name in names:
698 raise ValueError(f"Mask plane {plane.name!r} already exists.")
699 names.add(plane.name)
700 new_schema = MaskSchema([*retained, *new_planes], dtype=self._schema.dtype)
701 result = Mask(0, schema=new_schema, bbox=self._bbox, sky_projection=self._sky_projection)
702 # The retained planes are exactly the names common to both schemas, and
703 # ``result`` starts cleared and shares this mask's bbox, so ``update``
704 # transfers their pixel values (and nothing else) by name.
705 result.update(self)
706 return self._transfer_metadata(result, copy=True)
708 def serialize[P: pydantic.BaseModel](
709 self,
710 archive: OutputArchive[P],
711 *,
712 update_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
713 save_projection: bool = True,
714 add_offset_wcs: str | None = "A",
715 tile_shape: tuple[int, ...] | None = None,
716 options_name: str | None = None,
717 ) -> MaskSerializationModel[P]:
718 """Serialize the mask to an output archive.
720 Parameters
721 ----------
722 archive
723 Archive to write to.
724 update_header
725 A callback that will be given the FITS header for the HDU
726 containing this mask in order to add keys to it. This callback
727 may be provided but will not be called if the output format is not
728 FITS. As multiple HDUs may be added, this function may be called
729 multiple times.
730 save_projection
731 If `True`, save the `SkyProjection` attached to the image, if there
732 is one. This does not affect whether a FITS WCS corresponding to
733 the projection is written (it always is, if available, and if
734 ``add_offset_wcs`` is not ``" "``).
735 add_offset_wcs
736 A FITS WCS single-character suffix to use when adding a linear
737 WCS that maps the FITS array to the logical pixel coordinates
738 defined by ``bbox.start`` / ``yx0``. Set to `None` to not write
739 this WCS. If this is set to ``" "``, it will prevent the
740 `SkyProjection` from being saved as a FITS WCS.
741 tile_shape
742 The recommended shape of each tile, if the archive will save
743 the array in distinct tiles for faster subarray retrieval.
744 This is a hint; archives are not required to use this value.
745 options_name
746 Use this name to look up archive options.
747 """
748 if _archive_prefers_native_mask_arrays(archive):
749 # HDS presents array dimensions in Fortran order, which is the
750 # reverse of the h5py dataset shape. Store the in-memory trailing
751 # mask-byte axis first in HDF5 so Starlink tools see HDS axes
752 # (x, y, byte), without changing the bit packing within a pixel.
753 array_model = archive.add_array(
754 np.moveaxis(self._array, -1, 0),
755 update_header=update_header,
756 tile_shape=tile_shape,
757 options_name=options_name,
758 )
759 if not isinstance(array_model, ArrayReferenceModel): 759 ↛ 760line 759 didn't jump to line 760 because the condition on line 759 was never true
760 raise RuntimeError("Native mask arrays require reference array storage.")
761 array_model.shape = list(self._array.shape)
762 data: list[ArrayReferenceModel | InlineArrayModel] = [array_model]
763 else:
764 data = []
765 for schema_2d in self.schema.split(np.int32):
766 mask_2d = Mask(0, bbox=self.bbox, schema=schema_2d, sky_projection=self._sky_projection)
767 mask_2d.update(self)
768 data.append(
769 mask_2d._serialize_2d(
770 archive,
771 update_header=update_header,
772 add_offset_wcs=add_offset_wcs,
773 tile_shape=tile_shape,
774 options_name=options_name,
775 )
776 )
777 serialized_projection: SkyProjectionSerializationModel[P] | None = None
778 if save_projection and self.sky_projection is not None: 778 ↛ 779line 778 didn't jump to line 779 because the condition on line 778 was never true
779 serialized_projection = archive.serialize_direct("sky_projection", self.sky_projection.serialize)
780 serialized_dtype = NumberType.from_numpy(self.schema.dtype)
781 assert is_integer(serialized_dtype), "Mask dtypes should always be integers."
782 return MaskSerializationModel.model_construct(
783 data=data,
784 yx0=list(self.bbox.start),
785 planes=list(self.schema),
786 dtype=serialized_dtype,
787 sky_projection=serialized_projection,
788 metadata=self.metadata,
789 )
791 def _serialize_2d[P: pydantic.BaseModel](
792 self,
793 archive: OutputArchive[P],
794 *,
795 update_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
796 add_offset_wcs: str | None = "A",
797 tile_shape: tuple[int, ...] | None = None,
798 options_name: str | None = None,
799 ) -> ArrayReferenceModel | InlineArrayModel:
800 def _update_header(header: astropy.io.fits.Header) -> None:
801 update_header(header)
802 self.schema.update_header(header)
803 if self.sky_projection is not None and add_offset_wcs != " ":
804 if self.fits_wcs:
805 header.update(self.fits_wcs.to_header(relax=True))
806 if add_offset_wcs is not None: 806 ↛ exitline 806 didn't return from function '_update_header' because the condition on line 806 was always true
807 fits.add_offset_wcs(header, x=self.bbox.x.start, y=self.bbox.y.start, key=add_offset_wcs)
809 assert self.array.shape[2] == 1, "Mask should be split before calling this method."
810 return archive.add_array(
811 self._array[:, :, 0],
812 update_header=_update_header,
813 tile_shape=tile_shape,
814 options_name=options_name,
815 )
817 @staticmethod
818 def _get_archive_tree_type[P: pydantic.BaseModel](
819 pointer_type: type[P],
820 ) -> type[MaskSerializationModel[P]]:
821 """Return the serialization model type for this object for an archive
822 type that uses the given pointer type.
823 """
824 return MaskSerializationModel[pointer_type] # type: ignore
826 _archive_default_name: ClassVar[str] = "mask"
827 """The name this object should be serialized with when written as the
828 top-level object.
829 """
831 @staticmethod
832 def from_legacy(
833 legacy: Any,
834 plane_map: Mapping[str, MaskPlane] | None = None,
835 ) -> Mask:
836 """Convert from an `lsst.afw.image.Mask` instance.
838 Parameters
839 ----------
840 legacy
841 An `lsst.afw.image.Mask` instance. This will not share pixel
842 data with the new object.
843 plane_map
844 A mapping from legacy mask plane name to the new plane name and
845 description. If not provided, the right legacy mask plane will be
846 guessed, but this can depend on which mask planes the legacy
847 mask actually has set.
848 """
849 return Mask._from_legacy_array(
850 legacy.array,
851 legacy.getMaskPlaneDict(),
852 yx0=YX(y=legacy.getY0(), x=legacy.getX0()),
853 plane_map=plane_map,
854 )
856 def to_legacy(self, plane_map: Mapping[str, MaskPlane] | None = None) -> Any:
857 """Convert to an `lsst.afw.image.Mask` instance.
859 The pixel data will not be shared between the two objects.
861 Parameters
862 ----------
863 plane_map
864 A mapping from legacy mask plane name to the new plane name and
865 description.
866 """
867 import lsst.afw.image
868 import lsst.geom
870 result = lsst.afw.image.Mask(self.bbox.to_legacy())
871 if plane_map is None: 871 ↛ 873line 871 didn't jump to line 873 because the condition on line 871 was always true
872 plane_map = {plane.name: plane for plane in self.schema if plane is not None}
873 for old_name, new_plane in plane_map.items():
874 old_bit = result.addMaskPlane(old_name)
875 old_bitmask = 1 << old_bit
876 if old_bitmask == 2147483648: 876 ↛ 879line 876 didn't jump to line 879 because the condition on line 876 was never true
877 # afw uses int32 masks, but relies on overflow wrapping, which
878 # numpy doesn't like.
879 old_bitmask = -2147483648
880 if new_plane in self.schema: 880 ↛ 873line 880 didn't jump to line 873 because the condition on line 880 was always true
881 result.array[self.get(new_plane.name)] |= old_bitmask
882 return result
884 @staticmethod
885 def _from_legacy_array(
886 array2d: np.ndarray,
887 old_planes: Mapping[str, int],
888 *,
889 yx0: YX[int],
890 plane_map: Mapping[str, MaskPlane] | None = None,
891 sky_projection: SkyProjection | None = None,
892 ) -> Mask:
893 if plane_map is None: 893 ↛ 894line 893 didn't jump to line 894 because the condition on line 893 was never true
894 plane_map = _guess_legacy_plane_map(old_planes)
895 planes: list[MaskPlane] = list(plane_map.values()) if plane_map is not None else []
896 new_name_to_old_bitmask: dict[str, int] = {}
897 for old_name, old_bit in old_planes.items():
898 old_bitmask = 1 << old_bit
899 if old_bitmask == 2147483648: 899 ↛ 902line 899 didn't jump to line 902 because the condition on line 899 was never true
900 # afw uses int32 masks, but relies on overflow wrapping, which
901 # numpy doesn't like.
902 old_bitmask = -2147483648
903 if new_plane := plane_map.get(old_name):
904 # Already added to 'planes' at initialization.
905 new_name_to_old_bitmask[new_plane.name] = old_bitmask
906 else:
907 if n_orphaned := np.count_nonzero(array2d & old_bitmask): 907 ↛ 908line 907 didn't jump to line 908 because the condition on line 907 was never true
908 raise RuntimeError(
909 f"Legacy mask plane {old_name!r} is not remapped, "
910 f"but {n_orphaned} pixels have this bit set."
911 )
912 schema = MaskSchema(planes)
913 mask = Mask(0, schema=schema, yx0=yx0, shape=array2d.shape, sky_projection=sky_projection)
914 for new_name, old_bitmask in new_name_to_old_bitmask.items():
915 mask.set(new_name, array2d & old_bitmask)
916 return mask
918 @staticmethod
919 def read_legacy(
920 uri: ResourcePathExpression,
921 *,
922 plane_map: Mapping[str, MaskPlane] | None = None,
923 ext: str | int = 1,
924 fits_wcs_frame: Frame | None = None,
925 ) -> Mask:
926 """Read a FITS file written by `lsst.afw.image.Mask.writeFits`.
928 Parameters
929 ----------
930 uri
931 URI or file name.
932 plane_map
933 A mapping from legacy mask plane name to the new plane name and
934 description. If not provided, the right legacy mask plane will be
935 guessed, but this can depend on which mask planes the legacy
936 mask actually has set.
937 ext
938 Name or index of the FITS HDU to read.
939 fits_wcs_frame
940 If not `None` and the HDU containing the mask has a FITS WCS,
941 attach a `SkyProjection` to the returned mask by converting that
942 WCS.
943 """
944 opaque_metadata = fits.FitsOpaqueMetadata()
945 fs, fspath = ResourcePath(uri).to_fsspec()
946 with fs.open(fspath) as stream, astropy.io.fits.open(stream) as hdu_list:
947 opaque_metadata.extract_legacy_primary_header(hdu_list[0].header)
948 result = Mask._read_legacy_hdu(
949 hdu_list[ext], opaque_metadata, plane_map=plane_map, fits_wcs_frame=fits_wcs_frame
950 )
951 result._opaque_metadata = opaque_metadata
952 return result
954 @staticmethod
955 def _read_legacy_hdu(
956 hdu: astropy.io.fits.ImageHDU | astropy.io.fits.CompImageHDU | astropy.io.fits.BinTableHDU,
957 opaque_metadata: fits.FitsOpaqueMetadata,
958 plane_map: Mapping[str, MaskPlane] | None = None,
959 fits_wcs_frame: Frame | None = None,
960 strip_legacy_planes: bool = True,
961 ) -> Mask:
962 if isinstance(hdu, astropy.io.fits.BinTableHDU): 962 ↛ 963line 962 didn't jump to line 963 because the condition on line 962 was never true
963 hdu = astropy.io.fits.CompImageHDU(bintable=hdu)
964 yx0 = fits.read_yx0(hdu.header)
965 hdu.header.remove("LTV1", ignore_missing=True)
966 hdu.header.remove("LTV2", ignore_missing=True)
967 sky_projection: SkyProjection | None = None
968 if fits_wcs_frame is not None: 968 ↛ 969line 968 didn't jump to line 969 because the condition on line 968 was never true
969 try:
970 fits_wcs = astropy.wcs.WCS(hdu.header)
971 except KeyError:
972 pass
973 else:
974 sky_projection = SkyProjection.from_fits_wcs(
975 fits_wcs, pixel_frame=fits_wcs_frame, x0=yx0.x, y0=yx0.y
976 )
977 if any(card.keyword.startswith("MSKN") for card in hdu.header.cards):
978 # New ``lsst.images`` form: plane definitions are self-describing
979 # via MSKN/MSKM/MSKD cards, so no plane_map is needed. The on-disk
980 # array packs every plane into one element; ``set`` repacks each
981 # plane into the (default uint8) in-memory layout by name.
982 schema = MaskSchema.from_fits_header(hdu.header)
983 mask = Mask(0, schema=schema, yx0=yx0, shape=hdu.data.shape, sky_projection=sky_projection)
984 for n, plane in enumerate(schema):
985 if plane is not None: 985 ↛ 984line 985 didn't jump to line 984 because the condition on line 985 was always true
986 mask.set(plane.name, hdu.data & hdu.header.get(f"MSKM{n:04d}", 1 << n))
987 schema.strip_header(hdu.header)
988 else:
989 # Legacy ``lsst.afw.image`` form: bit indices in MP_* cards are
990 # mapped to new planes via ``plane_map``.
991 old_planes = MaskPlane.read_legacy(hdu.header, strip=strip_legacy_planes)
992 resolved_map = plane_map if plane_map is not None else _guess_legacy_plane_map(old_planes)
993 mask = Mask._from_legacy_array(
994 hdu.data, old_planes, yx0=yx0, plane_map=resolved_map, sky_projection=sky_projection
995 )
996 if not strip_legacy_planes:
997 # Keep the MP_ cards for backwards compatibility, but re-index
998 # them to the (reshuffled) positions of the new schema so a
999 # legacy reader sees each plane at the bit it is actually
1000 # packed into on disk.
1001 _reindex_legacy_plane_cards(hdu.header, old_planes, resolved_map, mask.schema)
1002 fits.strip_wcs_cards(hdu.header)
1003 hdu.header.strip()
1004 hdu.header.remove("EXTTYPE", ignore_missing=True)
1005 hdu.header.remove("INHERIT", ignore_missing=True)
1006 # afw set BUNIT on masks because of limitations in how FITS
1007 # metadata is handled there.
1008 hdu.header.remove("BUNIT", ignore_missing=True)
1009 opaque_metadata.add_header(hdu.header)
1010 return mask
1013class MaskSerializationModel[P: pydantic.BaseModel](ArchiveTree):
1014 """Pydantic model used to represent the serialized form of a `.Mask`."""
1016 SCHEMA_NAME: ClassVar[str] = "mask"
1017 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
1018 MIN_READ_VERSION: ClassVar[int] = 1
1019 PUBLIC_TYPE: ClassVar[type] = Mask
1021 data: list[ArrayReferenceModel | InlineArrayModel] = pydantic.Field(
1022 description="References to pixel data."
1023 )
1024 yx0: list[int] = pydantic.Field(
1025 description="Coordinate of the first pixels in the array, ordered (y, x)."
1026 )
1027 planes: list[MaskPlane | None] = pydantic.Field(description="Definitions of the bitplanes in the mask.")
1028 dtype: IntegerType = pydantic.Field(description="Data type of the in-memory mask.")
1029 sky_projection: SkyProjectionSerializationModel[P] | None = pydantic.Field(
1030 default=None,
1031 exclude_if=is_none,
1032 description="Projection that maps the logical pixel grid onto the sky.",
1033 )
1035 @property
1036 def bbox(self) -> Box:
1037 """The 2-d bounding box of the mask."""
1038 shape = self.data[0].shape
1039 if len(shape) == 3:
1040 shape = shape[:2]
1041 return Box.from_shape(shape, start=self.yx0)
1043 def deserialize(
1044 self,
1045 archive: InputArchive[Any],
1046 *,
1047 bbox: Box | None = None,
1048 strip_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
1049 **kwargs: Any,
1050 ) -> Mask:
1051 """Deserialize a mask from an input archive.
1053 Parameters
1054 ----------
1055 archive
1056 Archive to read from.
1057 bbox
1058 Bounding box of a subimage to read instead.
1059 strip_header
1060 A callable that strips out any FITS header cards added by the
1061 ``update_header`` argument in the corresponding call to
1062 `Mask.serialize`.
1063 **kwargs
1064 Unsupported keyword arguments are accepted only to provide better
1065 error messages (raising `serialization.InvalidParameterError`).
1066 """
1067 if kwargs: 1067 ↛ 1068line 1067 didn't jump to line 1068 because the condition on line 1067 was never true
1068 raise InvalidParameterError(f"Unrecognized parameters for Mask: {set(kwargs.keys())}.")
1070 def strip_header_and_legacy_planes(header: astropy.io.fits.Header) -> None:
1071 # The authoritative schema comes from the serialized tree, so drop
1072 # any legacy MP_* cards (written only for afw compatibility in the
1073 # legacy-cutout scenario) rather than carrying them as opaque
1074 # metadata, where they could drift out of sync or be re-propagated.
1075 strip_header(header)
1076 _strip_legacy_plane_cards(header)
1078 slices: tuple[slice, ...] | EllipsisType = ...
1079 if bbox is not None:
1080 slices = bbox.slice_within(self.bbox)
1081 else:
1082 bbox = self.bbox
1083 if not is_integer(self.dtype): 1083 ↛ 1084line 1083 didn't jump to line 1084 because the condition on line 1083 was never true
1084 raise ArchiveReadError(f"Mask array has a non-integer dtype: {self.dtype}.")
1085 schema = MaskSchema(self.planes, dtype=self.dtype.to_numpy())
1086 sky_projection = self.sky_projection.deserialize(archive) if self.sky_projection is not None else None
1087 if len(self.data) == 1 and tuple(self.data[0].shape) == tuple(self.bbox.shape) + (schema.mask_size,):
1088 storage_slices = slices if slices is ... else (slice(None),) + slices
1089 array = archive.get_array(
1090 self.data[0], strip_header=strip_header_and_legacy_planes, slices=storage_slices
1091 )
1092 array = np.moveaxis(array, 0, -1)
1093 return Mask(array, schema=schema, bbox=bbox, sky_projection=sky_projection)._finish_deserialize(
1094 self
1095 )
1096 result = Mask(0, schema=schema, bbox=bbox, sky_projection=sky_projection)
1097 schemas_2d = schema.split(np.int32)
1098 if len(schemas_2d) != len(self.data): 1098 ↛ 1099line 1098 didn't jump to line 1099 because the condition on line 1098 was never true
1099 raise ArchiveReadError(
1100 f"Number of mask arrays ({len(self.data)}) does not match expectation ({len(schemas_2d)})."
1101 )
1102 for array_model, schema_2d in zip(self.data, schemas_2d):
1103 mask_2d = self._deserialize_2d(
1104 array_model,
1105 schema_2d,
1106 bbox.start,
1107 archive,
1108 strip_header=strip_header_and_legacy_planes,
1109 slices=slices,
1110 )
1111 result.update(mask_2d)
1112 return result._finish_deserialize(self)
1114 @staticmethod
1115 def _deserialize_2d(
1116 ref: ArrayReferenceModel | InlineArrayModel,
1117 schema_2d: MaskSchema,
1118 yx0: Sequence[int],
1119 archive: InputArchive[Any],
1120 *,
1121 slices: tuple[slice, ...] | EllipsisType = ...,
1122 strip_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
1123 ) -> Mask:
1124 def _strip_header(header: astropy.io.fits.Header) -> None:
1125 strip_header(header)
1126 schema_2d.strip_header(header)
1127 fits.strip_wcs_cards(header)
1129 array_2d = archive.get_array(ref, strip_header=_strip_header, slices=slices)
1130 return Mask(array_2d[:, :, np.newaxis], schema=schema_2d, yx0=yx0)
1132 def deserialize_component(self, component: str, archive: InputArchive[Any], **kwargs: Any) -> Any:
1133 if kwargs:
1134 raise InvalidParameterError(f"Unsupported parameters for Mask components: {set(kwargs.keys())}.")
1135 return super().deserialize_component(component, archive)
1138def _archive_prefers_native_mask_arrays(archive: OutputArchive[Any]) -> bool:
1139 """Return whether an archive wants masks in their native 3-D layout."""
1140 current: Any = archive
1141 while current is not None:
1142 if getattr(current, "_prefer_native_mask_arrays", False):
1143 return True
1144 current = getattr(current, "_parent", None)
1145 return False
1148def get_legacy_visit_image_mask_planes() -> dict[str, MaskPlane]:
1149 """Return a mapping from legacy mask plane name to `MaskPlane` instance
1150 for LSST visit images, c. DP2.
1151 """
1152 return {
1153 "BAD": MaskPlane("BAD", "Bad pixel in the instrument, including bad amplifiers."),
1154 "SAT": MaskPlane(
1155 "SATURATED", "Pixel was saturated or affected by saturation in a neighboring pixel."
1156 ),
1157 "INTRP": MaskPlane("INTERPOLATED", "Original pixel value was interpolated."),
1158 "CR": MaskPlane("COSMIC_RAY", "A cosmic ray affected this pixel."),
1159 "EDGE": MaskPlane(
1160 "DETECTION_EDGE",
1161 "Pixel was too close to the edge to be considered for detection, "
1162 "due to the finite size of the detection kernel.",
1163 ),
1164 "DETECTED": MaskPlane("DETECTED", "Pixel was part of a detected source."),
1165 "SUSPECT": MaskPlane("SUSPECT", "Pixel was close to the saturation level. "),
1166 "NO_DATA": MaskPlane("NO_DATA", "No data was available for this pixel."),
1167 "VIGNETTED": MaskPlane("VIGNETTED", "Pixel was vignetted by the optics."),
1168 "PARTLY_VIGNETTED": MaskPlane("PARTLY_VIGNETTED", "Pixel was partly vignetted by the optics."),
1169 "CROSSTALK": MaskPlane("CROSSTALK", "Pixel was affected by crosstalk and corrected accordingly."),
1170 "ITL_DIP": MaskPlane(
1171 "ITL_DIP", "Pixel was affected by a dark vertical trail from a bright source, on an ITL CCD."
1172 ),
1173 "NOT_DEBLENDED": MaskPlane(
1174 "NOT_DEBLENDED",
1175 "Pixel belonged to a detection that was not deblended, usually due to size limits.",
1176 ),
1177 "SPIKE": MaskPlane(
1178 "SPIKE", "Pixel is in the neighborhood of a diffraction spike from a bright star."
1179 ),
1180 "UNMASKEDNAN": MaskPlane("UNMASKED_NAN", "Pixel was found to be NaN unexpectedly."),
1181 }
1184def get_legacy_difference_image_mask_planes() -> dict[str, MaskPlane]:
1185 """Return a mapping from legacy mask plane name to `MaskPlane` instance
1186 for LSST difference images, c. DP2.
1187 """
1188 result = get_legacy_visit_image_mask_planes()
1189 result["DETECTED_NEGATIVE"] = MaskPlane(
1190 "DETECTED_NEGATIVE", "Pixel was part of a detected source with negative flux."
1191 )
1192 result["SAT_TEMPLATE"] = MaskPlane("SAT_TEMPLATE", "Template pixel was saturated.")
1193 result["HIGH_VARIANCE"] = MaskPlane("HIGH_VARIANCE", "TODO[DM-55036]")
1194 result["STREAK"] = MaskPlane(
1195 "STREAK", "An extended streak (probably an artificial satellite) affected this pixel."
1196 )
1197 return result
1200def get_legacy_deep_coadd_mask_planes() -> dict[str, MaskPlane]:
1201 """Return a mapping from legacy mask plane name to `MaskPlane` instance
1202 for LSST deep coadds, c. DP2.
1203 """
1204 return {
1205 "NO_DATA": MaskPlane("NO_DATA", "No data was available for this pixel."),
1206 "INTRP": MaskPlane("INTERPOLATED", "Pixel value is the result of interpolating nearby good pixels."),
1207 "CR": MaskPlane(
1208 "COSMIC_RAY",
1209 "A cosmic ray affected this pixel on at least one input image (and was interpolated).",
1210 ),
1211 "SAT": MaskPlane(
1212 "SATURATED",
1213 "More than 10% of the potential input visits had a saturated pixel at this location "
1214 "('potential' because saturated pixel values are not actually propagated to the coadd). "
1215 "SATURATED always implies REJECTED, and is often a reason for NO_DATA.",
1216 ),
1217 "EDGE": MaskPlane(
1218 "DETECTION_EDGE",
1219 "Pixel was too close to the edge of the patch to be considered for detection, "
1220 "due to the finite size of the detection kernel.",
1221 ),
1222 "CLIPPED": MaskPlane(
1223 "CLIPPED",
1224 "Region was identified as a probable artifact when comparing multiple single-visit warps. "
1225 "CLIPPED always implies REJECTED.",
1226 ),
1227 "REJECTED": MaskPlane(
1228 "REJECTED",
1229 "At least one input visit was left out of the coadd for this pixel due to masking. "
1230 "REJECTED always implies INEXACT_PSF.",
1231 ),
1232 "DETECTED": MaskPlane("DETECTED", "Pixel was part of a detected source."),
1233 "INEXACT_PSF": MaskPlane(
1234 "INEXACT_PSF",
1235 "The set of visits contributing to this pixel differs from the set of visits "
1236 "contributing to the PSF model for its cell.",
1237 ),
1238 }
1241def get_legacy_non_cell_coadd_mask_planes() -> dict[str, MaskPlane]:
1242 """Return a mapping from legacy mask plane name to `MaskPlane` instance
1243 for LSST non-cell coadds such as ``template_coadd`` in DP2, and all
1244 DP1 coadds.
1246 These coadds carry the visit-level planes propagated from their input
1247 warps in addition to the coadd-specific planes, and flag chip edges with
1248 ``SENSOR_EDGE`` (cell coadds use ``CELL_EDGE`` instead).
1249 """
1250 result = get_legacy_deep_coadd_mask_planes()
1251 result["BAD"] = MaskPlane("BAD", "Bad pixel in the instrument, including bad amplifiers.")
1252 result["SUSPECT"] = MaskPlane("SUSPECT", "Pixel was close to the saturation level.")
1253 result["CROSSTALK"] = MaskPlane("CROSSTALK", "Pixel was affected by crosstalk and corrected accordingly.")
1254 result["DETECTED_NEGATIVE"] = MaskPlane(
1255 "DETECTED_NEGATIVE", "Pixel was part of a detected source with negative flux."
1256 )
1257 result["NOT_DEBLENDED"] = MaskPlane(
1258 "NOT_DEBLENDED",
1259 "Pixel belonged to a detection that was not deblended, usually due to size limits.",
1260 )
1261 result["UNMASKEDNAN"] = MaskPlane("UNMASKED_NAN", "Pixel was found to be NaN unexpectedly.")
1262 result["SENSOR_EDGE"] = MaskPlane(
1263 "SENSOR_EDGE",
1264 "Pixel is near the edge of a contributing sensor/chip, so the coadd PSF is discontinuous there.",
1265 )
1266 return result
1269def _guess_legacy_plane_map(old_planes: Mapping[str, int]) -> dict[str, MaskPlane]:
1270 """Guess which of the ``get_legacy_*_plane_map`` created the given mask
1271 plane dictionary and call it.
1272 """
1273 if "SAT_TEMPLATE" in old_planes: 1273 ↛ 1274line 1273 didn't jump to line 1274 because the condition on line 1273 was never true
1274 return get_legacy_difference_image_mask_planes()
1275 if "INEXACT_PSF" in old_planes:
1276 # Both cell and non-cell coadds have INEXACT_PSF, but only non-cell
1277 # (assemble_coadd) coadds flag chip edges with SENSOR_EDGE; cell coadds
1278 # use CELL_EDGE.
1279 if "SENSOR_EDGE" in old_planes:
1280 return get_legacy_non_cell_coadd_mask_planes()
1281 return get_legacy_deep_coadd_mask_planes()
1282 return get_legacy_visit_image_mask_planes()
1285def _reindex_legacy_plane_cards(
1286 header: astropy.io.fits.Header,
1287 old_planes: Mapping[str, int],
1288 plane_map: Mapping[str, MaskPlane],
1289 schema: MaskSchema,
1290) -> None:
1291 """Rewrite retained legacy ``MP_`` cards in place to match a reshuffled
1292 schema.
1294 Parameters
1295 ----------
1296 header
1297 Header whose ``MP_`` cards are updated in place.
1298 old_planes
1299 Mapping from legacy mask plane name to its original (on-disk) bit
1300 index, as returned by `MaskPlane.read_legacy`.
1301 plane_map
1302 Mapping from legacy mask plane name to the `MaskPlane` it was remapped
1303 to in ``schema``.
1304 schema
1305 The reconstructed schema that defines the new bit positions.
1307 Notes
1308 -----
1309 Each ``MP_<legacy name>`` card is set to the index that its remapped plane
1310 occupies in ``schema`` (equivalently, the ``MSKN`` index written on
1311 serialization). Cards for legacy planes that are not represented in the
1312 new schema are removed, since they no longer correspond to any stored bit.
1313 Legacy masks have at most 31 planes, so every plane maps to a single bit in
1314 one on-disk element and the index is unambiguous.
1315 """
1316 new_index = {plane.name: n for n, plane in enumerate(schema) if plane is not None}
1317 for old_name in old_planes:
1318 keyword = f"MP_{old_name}"
1319 new_plane = plane_map.get(old_name)
1320 if new_plane is not None and (index := new_index.get(new_plane.name)) is not None:
1321 header[keyword] = index
1322 else:
1323 del header[keyword]
1326def _strip_legacy_plane_cards(header: astropy.io.fits.Header) -> None:
1327 """Remove all legacy ``MP_*`` mask-plane cards from a FITS header.
1329 These are written only so that legacy tooling can read masks reconstructed
1330 from legacy cutouts; the ``lsst.images`` reader uses the serialized schema
1331 instead, so it strips them rather than carrying them as opaque metadata.
1332 """
1333 for keyword in [card.keyword for card in header.cards if card.keyword.startswith("MP_")]:
1334 del header[keyword]