Coverage for python/lsst/images/_mask.py: 67%
470 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-26 09:32 +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))
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 *,
903 sky_projection: SkyProjection[Any] | None = None,
904 ) -> Mask:
905 """Convert from an `lsst.afw.image.Mask` instance.
907 Parameters
908 ----------
909 legacy
910 An `lsst.afw.image.Mask` instance. This will not share pixel
911 data with the new object.
912 plane_map
913 A mapping from legacy mask plane name to the new plane name and
914 description. If not provided, the right legacy mask plane will be
915 guessed, but this can depend on which mask planes the legacy
916 mask actually has set.
917 sky_projection
918 Projection from pixels to xky.
919 """
920 return Mask._from_legacy_array(
921 legacy.array,
922 legacy.getMaskPlaneDict(),
923 yx0=YX(y=legacy.getY0(), x=legacy.getX0()),
924 plane_map=plane_map,
925 sky_projection=sky_projection,
926 )
928 def to_legacy(self, plane_map: Mapping[str, MaskPlane] | None = None) -> Any:
929 """Convert to an `lsst.afw.image.Mask` instance.
931 The pixel data will not be shared between the two objects.
933 Parameters
934 ----------
935 plane_map
936 A mapping from legacy mask plane name to the new plane name and
937 description.
938 """
939 import lsst.afw.image
940 import lsst.geom
942 result = lsst.afw.image.Mask(self.bbox.to_legacy())
943 if plane_map is None: 943 ↛ 945line 943 didn't jump to line 945 because the condition on line 943 was always true
944 plane_map = {plane.name: plane for plane in self.schema if plane is not None}
945 for old_name, new_plane in plane_map.items():
946 old_bit = result.addMaskPlane(old_name)
947 old_bitmask = 1 << old_bit
948 if old_bitmask == 2147483648: 948 ↛ 951line 948 didn't jump to line 951 because the condition on line 948 was never true
949 # afw uses int32 masks, but relies on overflow wrapping, which
950 # numpy doesn't like.
951 old_bitmask = -2147483648
952 if new_plane in self.schema: 952 ↛ 945line 952 didn't jump to line 945 because the condition on line 952 was always true
953 result.array[self.get(new_plane.name)] |= old_bitmask
954 return result
956 @staticmethod
957 def _from_legacy_array(
958 array2d: np.ndarray,
959 old_planes: Mapping[str, int],
960 *,
961 yx0: YX[int],
962 plane_map: Mapping[str, MaskPlane] | None = None,
963 sky_projection: SkyProjection | None = None,
964 ) -> Mask:
965 if plane_map is None: 965 ↛ 966line 965 didn't jump to line 966 because the condition on line 965 was never true
966 plane_map = _guess_legacy_plane_map(old_planes)
967 planes: list[MaskPlane] = list(plane_map.values()) if plane_map is not None else []
968 new_name_to_old_bitmask: dict[str, int] = {}
969 for old_name, old_bit in old_planes.items():
970 old_bitmask = 1 << old_bit
971 if old_bitmask == 2147483648: 971 ↛ 974line 971 didn't jump to line 974 because the condition on line 971 was never true
972 # afw uses int32 masks, but relies on overflow wrapping, which
973 # numpy doesn't like.
974 old_bitmask = -2147483648
975 if new_plane := plane_map.get(old_name):
976 # Already added to 'planes' at initialization.
977 new_name_to_old_bitmask[new_plane.name] = old_bitmask
978 else:
979 if n_orphaned := np.count_nonzero(array2d & old_bitmask): 979 ↛ 980line 979 didn't jump to line 980 because the condition on line 979 was never true
980 raise RuntimeError(
981 f"Legacy mask plane {old_name!r} is not remapped, "
982 f"but {n_orphaned} pixels have this bit set."
983 )
984 schema = MaskSchema(planes)
985 mask = Mask(0, schema=schema, yx0=yx0, shape=array2d.shape, sky_projection=sky_projection)
986 for new_name, old_bitmask in new_name_to_old_bitmask.items():
987 mask.set(new_name, array2d & old_bitmask)
988 return mask
990 @staticmethod
991 def read_legacy(
992 uri: ResourcePathExpression,
993 *,
994 plane_map: Mapping[str, MaskPlane] | None = None,
995 ext: str | int = 1,
996 fits_wcs_frame: Frame | None = None,
997 ) -> Mask:
998 """Read a FITS file written by `lsst.afw.image.Mask.writeFits`.
1000 Parameters
1001 ----------
1002 uri
1003 URI or file name.
1004 plane_map
1005 A mapping from legacy mask plane name to the new plane name and
1006 description. If not provided, the right legacy mask plane will be
1007 guessed, but this can depend on which mask planes the legacy
1008 mask actually has set.
1009 ext
1010 Name or index of the FITS HDU to read.
1011 fits_wcs_frame
1012 If not `None` and the HDU containing the mask has a FITS WCS,
1013 attach a `SkyProjection` to the returned mask by converting that
1014 WCS.
1015 """
1016 opaque_metadata = fits.FitsOpaqueMetadata()
1017 fs, fspath = ResourcePath(uri).to_fsspec()
1018 with fs.open(fspath) as stream, astropy.io.fits.open(stream) as hdu_list:
1019 opaque_metadata.extract_legacy_primary_header(hdu_list[0].header)
1020 result = Mask._read_legacy_hdu(
1021 hdu_list[ext], opaque_metadata, plane_map=plane_map, fits_wcs_frame=fits_wcs_frame
1022 )
1023 result._opaque_metadata = opaque_metadata
1024 return result
1026 @staticmethod
1027 def _read_legacy_hdu(
1028 hdu: astropy.io.fits.ImageHDU | astropy.io.fits.CompImageHDU | astropy.io.fits.BinTableHDU,
1029 opaque_metadata: fits.FitsOpaqueMetadata,
1030 plane_map: Mapping[str, MaskPlane] | None = None,
1031 fits_wcs_frame: Frame | None = None,
1032 strip_legacy_planes: bool = True,
1033 ) -> Mask:
1034 if isinstance(hdu, astropy.io.fits.BinTableHDU): 1034 ↛ 1035line 1034 didn't jump to line 1035 because the condition on line 1034 was never true
1035 hdu = astropy.io.fits.CompImageHDU(bintable=hdu)
1036 yx0 = fits.read_yx0(hdu.header)
1037 hdu.header.remove("LTV1", ignore_missing=True)
1038 hdu.header.remove("LTV2", ignore_missing=True)
1039 sky_projection: SkyProjection | None = None
1040 if fits_wcs_frame is not None: 1040 ↛ 1041line 1040 didn't jump to line 1041 because the condition on line 1040 was never true
1041 try:
1042 fits_wcs = astropy.wcs.WCS(hdu.header)
1043 except KeyError:
1044 pass
1045 else:
1046 sky_projection = SkyProjection.from_fits_wcs(
1047 fits_wcs, pixel_frame=fits_wcs_frame, x0=yx0.x, y0=yx0.y
1048 )
1049 if any(card.keyword.startswith("MSKN") for card in hdu.header.cards):
1050 # New ``lsst.images`` form: plane definitions are self-describing
1051 # via MSKN/MSKM/MSKD cards, so no plane_map is needed. The on-disk
1052 # array packs every plane into one element; ``set`` repacks each
1053 # plane into the (default uint8) in-memory layout by name.
1054 schema = MaskSchema.from_fits_header(hdu.header)
1055 mask = Mask(0, schema=schema, yx0=yx0, shape=hdu.data.shape, sky_projection=sky_projection)
1056 for n, plane in enumerate(schema):
1057 if plane is not None: 1057 ↛ 1056line 1057 didn't jump to line 1056 because the condition on line 1057 was always true
1058 mask.set(plane.name, hdu.data & hdu.header.get(f"MSKM{n:04d}", 1 << n))
1059 schema.strip_header(hdu.header)
1060 else:
1061 # Legacy ``lsst.afw.image`` form: bit indices in MP_* cards are
1062 # mapped to new planes via ``plane_map``.
1063 old_planes = MaskPlane.read_legacy(hdu.header, strip=strip_legacy_planes)
1064 resolved_map = plane_map if plane_map is not None else _guess_legacy_plane_map(old_planes)
1065 mask = Mask._from_legacy_array(
1066 hdu.data, old_planes, yx0=yx0, plane_map=resolved_map, sky_projection=sky_projection
1067 )
1068 if not strip_legacy_planes:
1069 # Keep the MP_ cards for backwards compatibility, but re-index
1070 # them to the (reshuffled) positions of the new schema so a
1071 # legacy reader sees each plane at the bit it is actually
1072 # packed into on disk.
1073 _reindex_legacy_plane_cards(hdu.header, old_planes, resolved_map, mask.schema)
1074 fits.strip_wcs_cards(hdu.header)
1075 hdu.header.strip()
1076 hdu.header.remove("EXTTYPE", ignore_missing=True)
1077 hdu.header.remove("INHERIT", ignore_missing=True)
1078 # afw set BUNIT on masks because of limitations in how FITS
1079 # metadata is handled there.
1080 hdu.header.remove("BUNIT", ignore_missing=True)
1081 opaque_metadata.add_header(hdu.header)
1082 return mask
1085class MaskSerializationModel[P: pydantic.BaseModel](ArchiveTree):
1086 """Pydantic model used to represent the serialized form of a `.Mask`."""
1088 SCHEMA_NAME: ClassVar[str] = "mask"
1089 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
1090 MIN_READ_VERSION: ClassVar[int] = 1
1091 PUBLIC_TYPE: ClassVar[type] = Mask
1093 data: list[ArrayReferenceModel | InlineArrayModel] = pydantic.Field(
1094 description="References to pixel data."
1095 )
1096 yx0: list[int] = pydantic.Field(
1097 description="Coordinate of the first pixels in the array, ordered (y, x)."
1098 )
1099 planes: list[MaskPlane | None] = pydantic.Field(description="Definitions of the bitplanes in the mask.")
1100 dtype: IntegerType = pydantic.Field(description="Data type of the in-memory mask.")
1101 sky_projection: SkyProjectionSerializationModel[P] | None = pydantic.Field(
1102 default=None,
1103 exclude_if=is_none,
1104 description="Projection that maps the logical pixel grid onto the sky.",
1105 )
1107 @property
1108 def bbox(self) -> Box:
1109 """The 2-d bounding box of the mask."""
1110 shape = self.data[0].shape
1111 if len(shape) == 3:
1112 shape = shape[:2]
1113 return Box.from_shape(shape, start=self.yx0)
1115 def deserialize(
1116 self,
1117 archive: InputArchive[Any],
1118 *,
1119 bbox: Box | None = None,
1120 strip_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
1121 **kwargs: Any,
1122 ) -> Mask:
1123 """Deserialize a mask from an input archive.
1125 Parameters
1126 ----------
1127 archive
1128 Archive to read from.
1129 bbox
1130 Bounding box of a subimage to read instead.
1131 strip_header
1132 A callable that strips out any FITS header cards added by the
1133 ``update_header`` argument in the corresponding call to
1134 `Mask.serialize`.
1135 **kwargs
1136 Unsupported keyword arguments are accepted only to provide better
1137 error messages (raising `serialization.InvalidParameterError`).
1138 """
1139 if kwargs: 1139 ↛ 1140line 1139 didn't jump to line 1140 because the condition on line 1139 was never true
1140 raise InvalidParameterError(f"Unrecognized parameters for Mask: {set(kwargs.keys())}.")
1142 def strip_header_and_legacy_planes(header: astropy.io.fits.Header) -> None:
1143 # The authoritative schema comes from the serialized tree, so drop
1144 # any legacy MP_* cards (written only for afw compatibility in the
1145 # legacy-cutout scenario) rather than carrying them as opaque
1146 # metadata, where they could drift out of sync or be re-propagated.
1147 strip_header(header)
1148 _strip_legacy_plane_cards(header)
1150 slices: tuple[slice, ...] | EllipsisType = ...
1151 if bbox is not None:
1152 slices = bbox.slice_within(self.bbox)
1153 else:
1154 bbox = self.bbox
1155 if not is_integer(self.dtype): 1155 ↛ 1156line 1155 didn't jump to line 1156 because the condition on line 1155 was never true
1156 raise ArchiveReadError(f"Mask array has a non-integer dtype: {self.dtype}.")
1157 schema = MaskSchema(self.planes, dtype=self.dtype.to_numpy())
1158 sky_projection = self.sky_projection.deserialize(archive) if self.sky_projection is not None else None
1159 if len(self.data) == 1 and tuple(self.data[0].shape) == tuple(self.bbox.shape) + (schema.mask_size,):
1160 storage_slices = slices if slices is ... else (slice(None),) + slices
1161 array = archive.get_array(
1162 self.data[0], strip_header=strip_header_and_legacy_planes, slices=storage_slices
1163 )
1164 array = np.moveaxis(array, 0, -1)
1165 return Mask(array, schema=schema, bbox=bbox, sky_projection=sky_projection)._finish_deserialize(
1166 self
1167 )
1168 result = Mask(0, schema=schema, bbox=bbox, sky_projection=sky_projection)
1169 schemas_2d = schema.split(np.int32)
1170 if len(schemas_2d) != len(self.data): 1170 ↛ 1171line 1170 didn't jump to line 1171 because the condition on line 1170 was never true
1171 raise ArchiveReadError(
1172 f"Number of mask arrays ({len(self.data)}) does not match expectation ({len(schemas_2d)})."
1173 )
1174 for array_model, schema_2d in zip(self.data, schemas_2d):
1175 mask_2d = self._deserialize_2d(
1176 array_model,
1177 schema_2d,
1178 bbox.start,
1179 archive,
1180 strip_header=strip_header_and_legacy_planes,
1181 slices=slices,
1182 )
1183 result.update(mask_2d)
1184 return result._finish_deserialize(self)
1186 @staticmethod
1187 def _deserialize_2d(
1188 ref: ArrayReferenceModel | InlineArrayModel,
1189 schema_2d: MaskSchema,
1190 yx0: Sequence[int],
1191 archive: InputArchive[Any],
1192 *,
1193 slices: tuple[slice, ...] | EllipsisType = ...,
1194 strip_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
1195 ) -> Mask:
1196 def _strip_header(header: astropy.io.fits.Header) -> None:
1197 strip_header(header)
1198 schema_2d.strip_header(header)
1199 fits.strip_wcs_cards(header)
1201 array_2d = archive.get_array(ref, strip_header=_strip_header, slices=slices)
1202 return Mask(array_2d[:, :, np.newaxis], schema=schema_2d, yx0=yx0)
1204 def deserialize_component(self, component: str, archive: InputArchive[Any], **kwargs: Any) -> Any:
1205 if kwargs:
1206 raise InvalidParameterError(f"Unsupported parameters for Mask components: {set(kwargs.keys())}.")
1207 return super().deserialize_component(component, archive)
1210def _archive_prefers_native_mask_arrays(archive: OutputArchive[Any]) -> bool:
1211 """Return whether an archive wants masks in their native 3-D layout."""
1212 current: Any = archive
1213 while current is not None:
1214 if getattr(current, "_prefer_native_mask_arrays", False):
1215 return True
1216 current = getattr(current, "_parent", None)
1217 return False
1220def get_legacy_visit_image_mask_planes() -> dict[str, MaskPlane]:
1221 """Return a mapping from legacy mask plane name to `MaskPlane` instance
1222 for LSST visit images, c. DP2.
1223 """
1224 return {
1225 "BAD": MaskPlane("BAD", "Bad pixel in the instrument, including bad amplifiers."),
1226 "SAT": MaskPlane(
1227 "SATURATED", "Pixel was saturated or affected by saturation in a neighboring pixel."
1228 ),
1229 "INTRP": MaskPlane("INTERPOLATED", "Original pixel value was interpolated."),
1230 "CR": MaskPlane("COSMIC_RAY", "A cosmic ray affected this pixel."),
1231 "EDGE": MaskPlane(
1232 "DETECTION_EDGE",
1233 "Pixel was too close to the edge to be considered for detection, "
1234 "due to the finite size of the detection kernel.",
1235 ),
1236 "DETECTED": MaskPlane("DETECTED", "Pixel was part of a detected source."),
1237 "SUSPECT": MaskPlane("SUSPECT", "Pixel was close to the saturation level. "),
1238 "NO_DATA": MaskPlane("NO_DATA", "No data was available for this pixel."),
1239 "VIGNETTED": MaskPlane("VIGNETTED", "Pixel was vignetted by the optics."),
1240 "PARTLY_VIGNETTED": MaskPlane("PARTLY_VIGNETTED", "Pixel was partly vignetted by the optics."),
1241 "CROSSTALK": MaskPlane("CROSSTALK", "Pixel was affected by crosstalk and corrected accordingly."),
1242 "ITL_DIP": MaskPlane(
1243 "ITL_DIP", "Pixel was affected by a dark vertical trail from a bright source, on an ITL CCD."
1244 ),
1245 "NOT_DEBLENDED": MaskPlane(
1246 "NOT_DEBLENDED",
1247 "Pixel belonged to a detection that was not deblended, usually due to size limits.",
1248 ),
1249 "SPIKE": MaskPlane(
1250 "SPIKE", "Pixel is in the neighborhood of a diffraction spike from a bright star."
1251 ),
1252 "UNMASKEDNAN": MaskPlane("UNMASKED_NAN", "Pixel was found to be NaN unexpectedly."),
1253 }
1256def get_legacy_difference_image_mask_planes() -> dict[str, MaskPlane]:
1257 """Return a mapping from legacy mask plane name to `MaskPlane` instance
1258 for LSST difference images, c. DP2.
1259 """
1260 result = get_legacy_visit_image_mask_planes()
1261 result["DETECTED_NEGATIVE"] = MaskPlane(
1262 "DETECTED_NEGATIVE", "Pixel was part of a detected source with negative flux."
1263 )
1264 result["SAT_TEMPLATE"] = MaskPlane("SAT_TEMPLATE", "Template pixel was saturated.")
1265 result["HIGH_VARIANCE"] = MaskPlane(
1266 "HIGH_VARIANCE", "Template pixel had fewer-than-usual input epochs and hence high noise."
1267 )
1268 result["STREAK"] = MaskPlane(
1269 "STREAK", "An extended streak (probably an artificial satellite) affected this pixel."
1270 )
1271 return result
1274def get_legacy_deep_coadd_mask_planes() -> dict[str, MaskPlane]:
1275 """Return a mapping from legacy mask plane name to `MaskPlane` instance
1276 for LSST deep coadds, c. DP2.
1277 """
1278 return {
1279 "NO_DATA": MaskPlane("NO_DATA", "No data was available for this pixel."),
1280 "INTRP": MaskPlane("INTERPOLATED", "Pixel value is the result of interpolating nearby good pixels."),
1281 "CR": MaskPlane(
1282 "COSMIC_RAY",
1283 "A cosmic ray affected this pixel on at least one input image (and was interpolated).",
1284 ),
1285 "SAT": MaskPlane(
1286 "SATURATED",
1287 "More than 10% of the potential input visits had a saturated pixel at this location "
1288 "('potential' because saturated pixel values are not actually propagated to the coadd). "
1289 "SATURATED always implies REJECTED, and is often a reason for NO_DATA.",
1290 ),
1291 "EDGE": MaskPlane(
1292 "DETECTION_EDGE",
1293 "Pixel was too close to the edge of the patch to be considered for detection, "
1294 "due to the finite size of the detection kernel.",
1295 ),
1296 "CLIPPED": MaskPlane(
1297 "CLIPPED",
1298 "Region was identified as a probable artifact when comparing multiple single-visit warps. "
1299 "CLIPPED always implies REJECTED.",
1300 ),
1301 "REJECTED": MaskPlane(
1302 "REJECTED",
1303 "At least one input visit was left out of the coadd for this pixel due to masking. "
1304 "REJECTED always implies INEXACT_PSF.",
1305 ),
1306 "DETECTED": MaskPlane("DETECTED", "Pixel was part of a detected source."),
1307 "INEXACT_PSF": MaskPlane(
1308 "INEXACT_PSF",
1309 "The set of visits contributing to this pixel differs from the set of visits "
1310 "contributing to the PSF model for its cell.",
1311 ),
1312 }
1315def get_legacy_non_cell_coadd_mask_planes() -> dict[str, MaskPlane]:
1316 """Return a mapping from legacy mask plane name to `MaskPlane` instance
1317 for LSST non-cell coadds such as ``template_coadd`` in DP2, and all
1318 DP1 coadds.
1320 These coadds carry the visit-level planes propagated from their input
1321 warps in addition to the coadd-specific planes, and flag chip edges with
1322 ``SENSOR_EDGE`` (cell coadds use ``CELL_EDGE`` instead).
1323 """
1324 result = get_legacy_deep_coadd_mask_planes()
1325 result["BAD"] = MaskPlane("BAD", "Bad pixel in the instrument, including bad amplifiers.")
1326 result["SUSPECT"] = MaskPlane("SUSPECT", "Pixel was close to the saturation level.")
1327 result["CROSSTALK"] = MaskPlane("CROSSTALK", "Pixel was affected by crosstalk and corrected accordingly.")
1328 result["DETECTED_NEGATIVE"] = MaskPlane(
1329 "DETECTED_NEGATIVE", "Pixel was part of a detected source with negative flux."
1330 )
1331 result["NOT_DEBLENDED"] = MaskPlane(
1332 "NOT_DEBLENDED",
1333 "Pixel belonged to a detection that was not deblended, usually due to size limits.",
1334 )
1335 result["UNMASKEDNAN"] = MaskPlane("UNMASKED_NAN", "Pixel was found to be NaN unexpectedly.")
1336 result["SENSOR_EDGE"] = MaskPlane(
1337 "SENSOR_EDGE",
1338 "Pixel is near the edge of a contributing sensor/chip, so the coadd PSF is discontinuous there.",
1339 )
1340 return result
1343def _guess_legacy_plane_map(old_planes: Mapping[str, int]) -> dict[str, MaskPlane]:
1344 """Guess which of the ``get_legacy_*_plane_map`` created the given mask
1345 plane dictionary and call it.
1346 """
1347 if "SAT_TEMPLATE" in old_planes: 1347 ↛ 1348line 1347 didn't jump to line 1348 because the condition on line 1347 was never true
1348 return get_legacy_difference_image_mask_planes()
1349 if "INEXACT_PSF" in old_planes:
1350 # Both cell and non-cell coadds have INEXACT_PSF, but only non-cell
1351 # (assemble_coadd) coadds flag chip edges with SENSOR_EDGE; cell coadds
1352 # use CELL_EDGE.
1353 if "SENSOR_EDGE" in old_planes:
1354 return get_legacy_non_cell_coadd_mask_planes()
1355 return get_legacy_deep_coadd_mask_planes()
1356 return get_legacy_visit_image_mask_planes()
1359def _reindex_legacy_plane_cards(
1360 header: astropy.io.fits.Header,
1361 old_planes: Mapping[str, int],
1362 plane_map: Mapping[str, MaskPlane],
1363 schema: MaskSchema,
1364) -> None:
1365 """Rewrite retained legacy ``MP_`` cards in place to match a reshuffled
1366 schema.
1368 Parameters
1369 ----------
1370 header
1371 Header whose ``MP_`` cards are updated in place.
1372 old_planes
1373 Mapping from legacy mask plane name to its original (on-disk) bit
1374 index, as returned by `MaskPlane.read_legacy`.
1375 plane_map
1376 Mapping from legacy mask plane name to the `MaskPlane` it was remapped
1377 to in ``schema``.
1378 schema
1379 The reconstructed schema that defines the new bit positions.
1381 Notes
1382 -----
1383 Each ``MP_<legacy name>`` card is set to the index that its remapped plane
1384 occupies in ``schema`` (equivalently, the ``MSKN`` index written on
1385 serialization). Cards for legacy planes that are not represented in the
1386 new schema are removed, since they no longer correspond to any stored bit.
1387 Legacy masks have at most 31 planes, so every plane maps to a single bit in
1388 one on-disk element and the index is unambiguous.
1389 """
1390 new_index = {plane.name: n for n, plane in enumerate(schema) if plane is not None}
1391 for old_name in old_planes:
1392 keyword = f"MP_{old_name}"
1393 new_plane = plane_map.get(old_name)
1394 if new_plane is not None and (index := new_index.get(new_plane.name)) is not None:
1395 header[keyword] = index
1396 else:
1397 del header[keyword]
1400def _strip_legacy_plane_cards(header: astropy.io.fits.Header) -> None:
1401 """Remove all legacy ``MP_*`` mask-plane cards from a FITS header.
1403 These are written only so that legacy tooling can read masks reconstructed
1404 from legacy cutouts; the ``lsst.images`` reader uses the serialized schema
1405 instead, so it strips them rather than carrying them as opaque metadata.
1406 """
1407 for keyword in [card.keyword for card in header.cards if card.keyword.startswith("MP_")]:
1408 del header[keyword]