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