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