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