Coverage for python/lsst/images/_mask.py: 70%
509 statements
« prev ^ index » next coverage.py v7.16.0, created at 2026-09-19 02:55 -0700
1# This file is part of lsst-images.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# Use of this source code is governed by a 3-clause BSD-style
10# license that can be found in the LICENSE file.
12from __future__ import annotations
14__all__ = (
15 "Mask",
16 "MaskPlane",
17 "MaskPlaneBit",
18 "MaskSchema",
19 "MaskSerializationModel",
20 "get_legacy_deep_coadd_mask_planes",
21 "get_legacy_difference_image_mask_planes",
22 "get_legacy_non_cell_coadd_mask_planes",
23 "get_legacy_visit_image_mask_planes",
24)
26import dataclasses
27import math
28from collections.abc import Callable, Iterable, Iterator, Mapping, Sequence, Set
29from types import EllipsisType
30from typing import TYPE_CHECKING, Any, ClassVar, cast
32import astropy.io.fits
33import astropy.wcs
34import numpy as np
35import numpy.typing as npt
36import pydantic
38from lsst.resources import ResourcePath, ResourcePathExpression
40from . import fits
41from ._generalized_image import GeneralizedImage
42from ._geom import YX, Box, NoOverlapError
43from ._transforms import Frame, SkyProjection, SkyProjectionSerializationModel
44from .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 compare(self, other: Mask) -> dict[str, tuple[int, int]]:
750 """Return a plane-by-plane comparison with another mask.
752 Parameters
753 ----------
754 other
755 The mask to compare against.
757 Returns
758 -------
759 `dict` [`str`, `tuple` [`int`, `int`]]
760 Dictionary mask planes as keys, where the values are:
762 - the number of pixels with that plane set in ``self`` but not
763 ``other``;
764 - the number of pixels with that plane set in ``other`` but not
765 ``self``.
767 Mask planes where the images have the same pixels set are not
768 included, so the result is empty for identical masks. Mask planes
769 present in only one operand are treated as though they were unset
770 for all pixels in the other operand.
772 Notes
773 -----
774 Two masks with different schemas can have an empty `compare` result
775 while still not satisfying an equality check, as long as the planes
776 they have in common have the same pixels set and any planes present in
777 only one operand have no pixels set.
778 """
779 if self.bbox != other.bbox:
780 raise ValueError("masks must have the same bounding box to compute a difference")
781 empty = np.zeros(self._array.shape[:-1], dtype=bool)
782 result: dict[str, tuple[int, int]] = {}
783 for name in sorted(self.schema.names | other.schema.names):
784 a = self.get(name) if name in self.schema.names else empty
785 b = other.get(name) if name in other.schema.names else empty
786 added = int(np.count_nonzero(a & ~b))
787 removed = int(np.count_nonzero(b & ~a))
788 if added or removed:
789 result[name] = (added, removed)
790 return result
792 def set(self, plane: str, boolean_mask: np.ndarray | EllipsisType = ...) -> None:
793 """Set a mask plane.
795 Parameters
796 ----------
797 plane
798 Name of the mask plane to set.
799 boolean_mask
800 A 2-d boolean array with the same shape as `bbox` that is `True`
801 where the bit for ``plane`` should be set and `False` where it
802 should be left unchanged (*not* set to zero). May be ``...`` to
803 set the bit everywhere.
804 """
805 bit = self.schema.bit(plane)
806 if boolean_mask is not ...: 806 ↛ 808line 806 didn't jump to line 808 because the condition on line 806 was always true
807 boolean_mask = boolean_mask.astype(bool)
808 self._array[boolean_mask, bit.index] |= bit.mask
810 def clear(self, plane: str | None = None, boolean_mask: np.ndarray | EllipsisType = ...) -> None:
811 """Clear one or more mask planes.
813 Parameters
814 ----------
815 plane
816 Name of the mask plane to set. If `None` all mask planes are
817 cleared.
818 boolean_mask
819 A 2-d boolean array with the same shape as `bbox` that is `True`
820 where the bit for ``plane`` should be cleared and `False` where it
821 should be left unchanged. May be ``...`` to clear the bit
822 everywhere.
823 """
824 if boolean_mask is not ...: 824 ↛ 825line 824 didn't jump to line 825 because the condition on line 824 was never true
825 boolean_mask = boolean_mask.astype(bool)
826 if plane is None: 826 ↛ 829line 826 didn't jump to line 829 because the condition on line 826 was always true
827 self._array[boolean_mask, :] = 0
828 else:
829 bit = self.schema.bit(plane)
830 self._array[boolean_mask, bit.index] &= ~bit.mask
832 def add_plane(self, name: str, description: str) -> Mask:
833 """Return a new mask with one additional mask plane.
835 This is a convenience wrapper around `add_planes` for the common case
836 of adding a single plane.
838 Parameters
839 ----------
840 name
841 Unique name for the new mask plane.
842 description
843 Human-readable documentation for the new mask plane.
845 Returns
846 -------
847 `Mask`
848 A new mask whose schema includes the new plane; see `add_planes`
849 for the reallocation and view semantics.
851 Raises
852 ------
853 ValueError
854 Raised if a plane named ``name`` already exists.
855 """
856 return self.add_planes([MaskPlane(name, description)])
858 def add_planes(self, planes: Iterable[MaskPlane | None], *, drop: Iterable[str] = ()) -> Mask:
859 """Return a new mask with planes added and/or dropped.
861 Parameters
862 ----------
863 planes
864 New mask planes to append, in order, after the planes retained
865 from this mask. `None` entries reserve unused bits (placeholders),
866 exactly as in `MaskSchema`.
867 drop
868 Names of existing planes to remove from the schema.
870 Returns
871 -------
872 `Mask`
873 A new mask with the updated schema. Retained planes keep their
874 pixel values (copied by name); newly added planes start cleared.
876 Raises
877 ------
878 ValueError
879 Raised if a name in ``drop`` is not an existing plane, or if a
880 plane in ``planes`` collides with a retained plane name.
882 Notes
883 -----
884 Adding or dropping planes always reallocates the backing array and
885 returns a new `Mask`; this mask is left unchanged and any views or
886 subimages of it continue to refer to the original array with the
887 original schema. This is deliberate: there is no way to update the
888 schema of an existing view, and a stale view must never set bits that
889 its now-outdated schema regards as unused. Dropping a plane compacts
890 the schema, so planes after it are reassigned to lower bits and the
891 pixel values are repacked by plane name to match.
892 """
893 drop_set = set(drop)
894 if unknown := drop_set - set(self._schema.names):
895 raise ValueError(f"Cannot drop mask planes that do not exist: {sorted(unknown)}.")
896 retained = [plane for plane in self._schema if plane is None or plane.name not in drop_set]
897 names = {plane.name for plane in retained if plane is not None}
898 new_planes = list(planes)
899 for plane in new_planes:
900 if plane is None:
901 continue
902 if plane.name in names:
903 raise ValueError(f"Mask plane {plane.name!r} already exists.")
904 names.add(plane.name)
905 new_schema = MaskSchema([*retained, *new_planes], dtype=self._schema.dtype)
906 result = Mask(0, schema=new_schema, bbox=self._bbox, sky_projection=self._sky_projection)
907 # The retained planes are exactly the names common to both schemas, and
908 # ``result`` starts cleared and shares this mask's bbox, so ``update``
909 # transfers their pixel values (and nothing else) by name.
910 result.update(self)
911 return self._transfer_metadata(result, copy=True)
913 def serialize[P: pydantic.BaseModel](
914 self,
915 archive: OutputArchive[P],
916 *,
917 update_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
918 save_projection: bool = True,
919 add_offset_wcs: str | None = "A",
920 tile_shape: tuple[int, ...] | None = None,
921 options_name: str | None = None,
922 ) -> MaskSerializationModel[P]:
923 """Serialize the mask to an output archive.
925 Parameters
926 ----------
927 archive
928 Archive to write to.
929 update_header
930 A callback that will be given the FITS header for the HDU
931 containing this mask in order to add keys to it. This callback
932 may be provided but will not be called if the output format is not
933 FITS. As multiple HDUs may be added, this function may be called
934 multiple times.
935 save_projection
936 If `True`, save the `SkyProjection` attached to the image, if there
937 is one. This does not affect whether a FITS WCS corresponding to
938 the projection is written (it always is, if available, and if
939 ``add_offset_wcs`` is not ``" "``).
940 add_offset_wcs
941 A FITS WCS single-character suffix to use when adding a linear
942 WCS that maps the FITS array to the logical pixel coordinates
943 defined by ``bbox.start`` / ``yx0``. Set to `None` to not write
944 this WCS. If this is set to ``" "``, it will prevent the
945 `SkyProjection` from being saved as a FITS WCS.
946 tile_shape
947 The recommended shape of each tile, if the archive will save
948 the array in distinct tiles for faster subarray retrieval.
949 This is a hint; archives are not required to use this value.
950 options_name
951 Use this name to look up archive options.
952 """
953 if _archive_prefers_native_mask_arrays(archive):
954 # HDS presents array dimensions in Fortran order, which is the
955 # reverse of the h5py dataset shape. Store the in-memory trailing
956 # mask-byte axis first in HDF5 so Starlink tools see HDS axes
957 # (x, y, byte), without changing the bit packing within a pixel.
958 array_model = archive.add_array(
959 np.moveaxis(self._array, -1, 0),
960 update_header=update_header,
961 tile_shape=tile_shape,
962 options_name=options_name,
963 )
964 if not isinstance(array_model, ArrayReferenceModel): 964 ↛ 965line 964 didn't jump to line 965 because the condition on line 964 was never true
965 raise RuntimeError("Native mask arrays require reference array storage.")
966 array_model.shape = list(self._array.shape)
967 data: list[ArrayReferenceModel | InlineArrayModel] = [array_model]
968 else:
969 data = []
970 for schema_2d in self.schema.split(np.int32):
971 mask_2d = Mask(0, bbox=self.bbox, schema=schema_2d, sky_projection=self._sky_projection)
972 mask_2d.update(self)
973 data.append(
974 mask_2d._serialize_2d(
975 archive,
976 update_header=update_header,
977 add_offset_wcs=add_offset_wcs,
978 tile_shape=tile_shape,
979 options_name=options_name,
980 )
981 )
982 serialized_projection: SkyProjectionSerializationModel[P] | None = None
983 if save_projection and self.sky_projection is not None: 983 ↛ 984line 983 didn't jump to line 984 because the condition on line 983 was never true
984 serialized_projection = archive.serialize_direct("sky_projection", self.sky_projection.serialize)
985 serialized_dtype = NumberType.from_numpy(self.schema.dtype)
986 assert is_integer(serialized_dtype), "Mask dtypes should always be integers."
987 return MaskSerializationModel.model_construct(
988 data=data,
989 yx0=list(self.bbox.start),
990 planes=list(self.schema),
991 dtype=serialized_dtype,
992 sky_projection=serialized_projection,
993 metadata=self.metadata,
994 )
996 def _serialize_2d[P: pydantic.BaseModel](
997 self,
998 archive: OutputArchive[P],
999 *,
1000 update_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
1001 add_offset_wcs: str | None = "A",
1002 tile_shape: tuple[int, ...] | None = None,
1003 options_name: str | None = None,
1004 ) -> ArrayReferenceModel | InlineArrayModel:
1005 def _update_header(header: astropy.io.fits.Header) -> None:
1006 update_header(header)
1007 self.schema.update_header(header)
1008 if self.sky_projection is not None and add_offset_wcs != " ":
1009 if self.fits_wcs:
1010 header.update(self.fits_wcs.to_header(relax=True))
1011 if add_offset_wcs is not None: 1011 ↛ exitline 1011 didn't return from function '_update_header' because the condition on line 1011 was always true
1012 fits.add_offset_wcs(header, x=self.bbox.x.start, y=self.bbox.y.start, key=add_offset_wcs)
1014 assert self.array.shape[2] == 1, "Mask should be split before calling this method."
1015 return archive.add_array(
1016 self._array[:, :, 0],
1017 update_header=_update_header,
1018 tile_shape=tile_shape,
1019 options_name=options_name,
1020 )
1022 @staticmethod
1023 def _get_archive_tree_type[P: pydantic.BaseModel](
1024 pointer_type: type[P],
1025 ) -> type[MaskSerializationModel[P]]:
1026 """Return the serialization model type for this object for an archive
1027 type that uses the given pointer type.
1028 """
1029 return MaskSerializationModel[pointer_type] # type: ignore
1031 _archive_default_name: ClassVar[str] = "mask"
1032 """The name this object should be serialized with when written as the
1033 top-level object.
1034 """
1036 @staticmethod
1037 def from_legacy(
1038 legacy: Any,
1039 plane_map: Mapping[str, MaskPlane] | None = None,
1040 *,
1041 sky_projection: SkyProjection[Any] | None = None,
1042 ) -> Mask:
1043 """Convert from an `lsst.afw.image.Mask` instance.
1045 Parameters
1046 ----------
1047 legacy
1048 An `lsst.afw.image.Mask` instance. This will not share pixel
1049 data with the new object.
1050 plane_map
1051 A mapping from legacy mask plane name to the new plane name and
1052 description. If not provided, the right legacy mask plane will be
1053 guessed, but this can depend on which mask planes the legacy
1054 mask actually has set.
1055 sky_projection
1056 Projection from pixels to xky.
1057 """
1058 return Mask._from_legacy_array(
1059 legacy.array,
1060 legacy.getMaskPlaneDict(),
1061 yx0=YX(y=legacy.getY0(), x=legacy.getX0()),
1062 plane_map=plane_map,
1063 sky_projection=sky_projection,
1064 )
1066 def to_legacy(self, plane_map: Mapping[str, MaskPlane] | None = None) -> Any:
1067 """Convert to an `lsst.afw.image.Mask` instance.
1069 The pixel data will not be shared between the two objects.
1071 Parameters
1072 ----------
1073 plane_map
1074 A mapping from legacy mask plane name to the new plane name and
1075 description.
1076 """
1077 import lsst.afw.image
1078 import lsst.geom
1080 result = lsst.afw.image.Mask(self.bbox.to_legacy())
1081 if plane_map is None: 1081 ↛ 1083line 1081 didn't jump to line 1083 because the condition on line 1081 was always true
1082 plane_map = {plane.name: plane for plane in self.schema if plane is not None}
1083 for old_name, new_plane in plane_map.items():
1084 old_bit = result.addMaskPlane(old_name)
1085 old_bitmask = 1 << old_bit
1086 if old_bitmask == 2147483648: 1086 ↛ 1089line 1086 didn't jump to line 1089 because the condition on line 1086 was never true
1087 # afw uses int32 masks, but relies on overflow wrapping, which
1088 # numpy doesn't like.
1089 old_bitmask = -2147483648
1090 if new_plane in self.schema: 1090 ↛ 1083line 1090 didn't jump to line 1083 because the condition on line 1090 was always true
1091 result.array[self.get(new_plane.name)] |= old_bitmask
1092 return result
1094 @staticmethod
1095 def _from_legacy_array(
1096 array2d: np.ndarray,
1097 old_planes: Mapping[str, int],
1098 *,
1099 yx0: YX[int],
1100 plane_map: Mapping[str, MaskPlane] | None = None,
1101 sky_projection: SkyProjection | None = None,
1102 ) -> Mask:
1103 if plane_map is None: 1103 ↛ 1104line 1103 didn't jump to line 1104 because the condition on line 1103 was never true
1104 plane_map = _guess_legacy_plane_map(old_planes)
1105 planes: list[MaskPlane] = list(plane_map.values()) if plane_map is not None else []
1106 new_name_to_old_bitmask: dict[str, int] = {}
1107 for old_name, old_bit in old_planes.items():
1108 old_bitmask = 1 << old_bit
1109 if old_bitmask == 2147483648: 1109 ↛ 1112line 1109 didn't jump to line 1112 because the condition on line 1109 was never true
1110 # afw uses int32 masks, but relies on overflow wrapping, which
1111 # numpy doesn't like.
1112 old_bitmask = -2147483648
1113 if new_plane := plane_map.get(old_name):
1114 # Already added to 'planes' at initialization.
1115 new_name_to_old_bitmask[new_plane.name] = old_bitmask
1116 else:
1117 if n_orphaned := np.count_nonzero(array2d & old_bitmask): 1117 ↛ 1118line 1117 didn't jump to line 1118 because the condition on line 1117 was never true
1118 raise RuntimeError(
1119 f"Legacy mask plane {old_name!r} is not remapped, "
1120 f"but {n_orphaned} pixels have this bit set."
1121 )
1122 schema = MaskSchema(planes)
1123 mask = Mask(0, schema=schema, yx0=yx0, shape=array2d.shape, sky_projection=sky_projection)
1124 for new_name, old_bitmask in new_name_to_old_bitmask.items():
1125 mask.set(new_name, array2d & old_bitmask)
1126 return mask
1128 @staticmethod
1129 def read_legacy(
1130 uri: ResourcePathExpression,
1131 *,
1132 plane_map: Mapping[str, MaskPlane] | None = None,
1133 ext: str | int = 1,
1134 fits_wcs_frame: Frame | None = None,
1135 ) -> Mask:
1136 """Read a FITS file written by `lsst.afw.image.Mask.writeFits`.
1138 Parameters
1139 ----------
1140 uri
1141 URI or file name.
1142 plane_map
1143 A mapping from legacy mask plane name to the new plane name and
1144 description. If not provided, the right legacy mask plane will be
1145 guessed, but this can depend on which mask planes the legacy
1146 mask actually has set.
1147 ext
1148 Name or index of the FITS HDU to read.
1149 fits_wcs_frame
1150 If not `None` and the HDU containing the mask has a FITS WCS,
1151 attach a `SkyProjection` to the returned mask by converting that
1152 WCS.
1153 """
1154 opaque_metadata = fits.FitsOpaqueMetadata()
1155 fs, fspath = ResourcePath(uri).to_fsspec()
1156 with fs.open(fspath) as stream, astropy.io.fits.open(stream) as hdu_list:
1157 opaque_metadata.extract_legacy_primary_header(hdu_list[0].header)
1158 result = Mask._read_legacy_hdu(
1159 hdu_list[ext], opaque_metadata, plane_map=plane_map, fits_wcs_frame=fits_wcs_frame
1160 )
1161 result._opaque_metadata = opaque_metadata
1162 return result
1164 @staticmethod
1165 def _read_legacy_hdu(
1166 hdu: astropy.io.fits.ImageHDU | astropy.io.fits.CompImageHDU | astropy.io.fits.BinTableHDU,
1167 opaque_metadata: fits.FitsOpaqueMetadata,
1168 plane_map: Mapping[str, MaskPlane] | None = None,
1169 fits_wcs_frame: Frame | None = None,
1170 strip_legacy_planes: bool = True,
1171 ) -> Mask:
1172 if isinstance(hdu, astropy.io.fits.BinTableHDU): 1172 ↛ 1173line 1172 didn't jump to line 1173 because the condition on line 1172 was never true
1173 hdu = astropy.io.fits.CompImageHDU(bintable=hdu)
1174 yx0 = fits.read_yx0(hdu.header)
1175 hdu.header.remove("LTV1", ignore_missing=True)
1176 hdu.header.remove("LTV2", ignore_missing=True)
1177 sky_projection: SkyProjection | None = None
1178 if fits_wcs_frame is not None: 1178 ↛ 1179line 1178 didn't jump to line 1179 because the condition on line 1178 was never true
1179 try:
1180 fits_wcs = astropy.wcs.WCS(hdu.header)
1181 except KeyError:
1182 pass
1183 else:
1184 sky_projection = SkyProjection.from_fits_wcs(
1185 fits_wcs, pixel_frame=fits_wcs_frame, x0=yx0.x, y0=yx0.y
1186 )
1187 if any(card.keyword.startswith("MSKN") for card in hdu.header.cards):
1188 # New ``lsst.images`` form: plane definitions are self-describing
1189 # via MSKN/MSKM/MSKD cards, so no plane_map is needed. The on-disk
1190 # array packs every plane into one element; ``set`` repacks each
1191 # plane into the (default uint8) in-memory layout by name.
1192 schema = MaskSchema.from_fits_header(hdu.header)
1193 mask = Mask(0, schema=schema, yx0=yx0, shape=hdu.data.shape, sky_projection=sky_projection)
1194 for n, plane in enumerate(schema):
1195 if plane is not None: 1195 ↛ 1194line 1195 didn't jump to line 1194 because the condition on line 1195 was always true
1196 mask.set(plane.name, hdu.data & hdu.header.get(f"MSKM{n:04d}", 1 << n))
1197 schema.strip_header(hdu.header)
1198 else:
1199 # Legacy ``lsst.afw.image`` form: bit indices in MP_* cards are
1200 # mapped to new planes via ``plane_map``.
1201 old_planes = MaskPlane.read_legacy(hdu.header, strip=strip_legacy_planes)
1202 resolved_map = plane_map if plane_map is not None else _guess_legacy_plane_map(old_planes)
1203 mask = Mask._from_legacy_array(
1204 hdu.data, old_planes, yx0=yx0, plane_map=resolved_map, sky_projection=sky_projection
1205 )
1206 if not strip_legacy_planes:
1207 # Keep the MP_ cards for backwards compatibility, but re-index
1208 # them to the (reshuffled) positions of the new schema so a
1209 # legacy reader sees each plane at the bit it is actually
1210 # packed into on disk.
1211 _reindex_legacy_plane_cards(hdu.header, old_planes, resolved_map, mask.schema)
1212 fits.strip_wcs_cards(hdu.header)
1213 hdu.header.strip()
1214 hdu.header.remove("EXTTYPE", ignore_missing=True)
1215 hdu.header.remove("INHERIT", ignore_missing=True)
1216 # afw set BUNIT on masks because of limitations in how FITS
1217 # metadata is handled there.
1218 hdu.header.remove("BUNIT", ignore_missing=True)
1219 opaque_metadata.add_header(hdu.header)
1220 return mask
1223class MaskSerializationModel[P: pydantic.BaseModel](ArchiveTree):
1224 """Pydantic model used to represent the serialized form of a `.Mask`."""
1226 SCHEMA_NAME: ClassVar[str] = "mask"
1227 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
1228 MIN_READ_VERSION: ClassVar[int] = 1
1229 PUBLIC_TYPE: ClassVar[type] = Mask
1231 data: list[ArrayReferenceModel | InlineArrayModel] = pydantic.Field(
1232 description="References to pixel data."
1233 )
1234 yx0: list[int] = pydantic.Field(
1235 description="Coordinate of the first pixels in the array, ordered (y, x)."
1236 )
1237 planes: list[MaskPlane | None] = pydantic.Field(description="Definitions of the bitplanes in the mask.")
1238 dtype: IntegerType = pydantic.Field(description="Data type of the in-memory mask.")
1239 sky_projection: SkyProjectionSerializationModel[P] | None = pydantic.Field(
1240 default=None,
1241 exclude_if=is_none,
1242 description="Projection that maps the logical pixel grid onto the sky.",
1243 )
1245 @property
1246 def bbox(self) -> Box:
1247 """The 2-d bounding box of the mask."""
1248 shape = self.data[0].shape
1249 if len(shape) == 3:
1250 shape = shape[:2]
1251 return Box.from_shape(shape, start=self.yx0)
1253 def deserialize(
1254 self,
1255 archive: InputArchive[Any],
1256 *,
1257 bbox: Box | None = None,
1258 strip_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
1259 **kwargs: Any,
1260 ) -> Mask:
1261 """Deserialize a mask from an input archive.
1263 Parameters
1264 ----------
1265 archive
1266 Archive to read from.
1267 bbox
1268 Bounding box of a subimage to read instead.
1269 strip_header
1270 A callable that strips out any FITS header cards added by the
1271 ``update_header`` argument in the corresponding call to
1272 `Mask.serialize`.
1273 **kwargs
1274 Unsupported keyword arguments are accepted only to provide better
1275 error messages (raising `serialization.InvalidParameterError`).
1276 """
1277 if kwargs: 1277 ↛ 1278line 1277 didn't jump to line 1278 because the condition on line 1277 was never true
1278 raise InvalidParameterError(f"Unrecognized parameters for Mask: {set(kwargs.keys())}.")
1280 def strip_header_and_legacy_planes(header: astropy.io.fits.Header) -> None:
1281 # The authoritative schema comes from the serialized tree, so drop
1282 # any legacy MP_* cards (written only for afw compatibility in the
1283 # legacy-cutout scenario) rather than carrying them as opaque
1284 # metadata, where they could drift out of sync or be re-propagated.
1285 strip_header(header)
1286 _strip_legacy_plane_cards(header)
1288 slices: tuple[slice, ...] | EllipsisType = ...
1289 if bbox is not None:
1290 slices = bbox.slice_within(self.bbox)
1291 else:
1292 bbox = self.bbox
1293 if not is_integer(self.dtype): 1293 ↛ 1294line 1293 didn't jump to line 1294 because the condition on line 1293 was never true
1294 raise ArchiveReadError(f"Mask array has a non-integer dtype: {self.dtype}.")
1295 schema = MaskSchema(self.planes, dtype=self.dtype.to_numpy())
1296 sky_projection = self.sky_projection.deserialize(archive) if self.sky_projection is not None else None
1297 if len(self.data) == 1 and tuple(self.data[0].shape) == tuple(self.bbox.shape) + (schema.mask_size,):
1298 storage_slices = slices if slices is ... else (slice(None),) + slices
1299 array = archive.get_array(
1300 self.data[0], strip_header=strip_header_and_legacy_planes, slices=storage_slices
1301 )
1302 array = np.moveaxis(array, 0, -1)
1303 return Mask(array, schema=schema, bbox=bbox, sky_projection=sky_projection)._finish_deserialize(
1304 self
1305 )
1306 result = Mask(0, schema=schema, bbox=bbox, sky_projection=sky_projection)
1307 schemas_2d = schema.split(np.int32)
1308 if len(schemas_2d) != len(self.data): 1308 ↛ 1309line 1308 didn't jump to line 1309 because the condition on line 1308 was never true
1309 raise ArchiveReadError(
1310 f"Number of mask arrays ({len(self.data)}) does not match expectation ({len(schemas_2d)})."
1311 )
1312 for array_model, schema_2d in zip(self.data, schemas_2d):
1313 mask_2d = self._deserialize_2d(
1314 array_model,
1315 schema_2d,
1316 bbox.start,
1317 archive,
1318 strip_header=strip_header_and_legacy_planes,
1319 slices=slices,
1320 )
1321 result.update(mask_2d)
1322 return result._finish_deserialize(self)
1324 @staticmethod
1325 def _deserialize_2d(
1326 ref: ArrayReferenceModel | InlineArrayModel,
1327 schema_2d: MaskSchema,
1328 yx0: Sequence[int],
1329 archive: InputArchive[Any],
1330 *,
1331 slices: tuple[slice, ...] | EllipsisType = ...,
1332 strip_header: Callable[[astropy.io.fits.Header], None] = no_header_updates,
1333 ) -> Mask:
1334 def _strip_header(header: astropy.io.fits.Header) -> None:
1335 strip_header(header)
1336 schema_2d.strip_header(header)
1337 fits.strip_wcs_cards(header)
1339 array_2d = archive.get_array(ref, strip_header=_strip_header, slices=slices)
1340 return Mask(array_2d[:, :, np.newaxis], schema=schema_2d, yx0=yx0)
1342 def deserialize_component(self, component: str, archive: InputArchive[Any], **kwargs: Any) -> Any:
1343 if kwargs:
1344 raise InvalidParameterError(f"Unsupported parameters for Mask components: {set(kwargs.keys())}.")
1345 return super().deserialize_component(component, archive)
1348def _archive_prefers_native_mask_arrays(archive: OutputArchive[Any]) -> bool:
1349 """Return whether an archive wants masks in their native 3-D layout."""
1350 current: Any = archive
1351 while current is not None:
1352 if getattr(current, "_prefer_native_mask_arrays", False):
1353 return True
1354 current = getattr(current, "_parent", None)
1355 return False
1358def get_legacy_visit_image_mask_planes() -> dict[str, MaskPlane]:
1359 """Return a mapping from legacy mask plane name to `MaskPlane` instance
1360 for LSST visit images, c. DP2.
1361 """
1362 return {
1363 "BAD": MaskPlane("BAD", "Bad pixel in the instrument, including bad amplifiers."),
1364 "SAT": MaskPlane(
1365 "SATURATED", "Pixel was saturated or affected by saturation in a neighboring pixel."
1366 ),
1367 "INTRP": MaskPlane("INTERPOLATED", "Original pixel value was interpolated."),
1368 "CR": MaskPlane("COSMIC_RAY", "A cosmic ray affected this pixel."),
1369 "EDGE": MaskPlane(
1370 "DETECTION_EDGE",
1371 "Pixel was too close to the edge to be considered for detection, "
1372 "due to the finite size of the detection kernel.",
1373 ),
1374 "DETECTED": MaskPlane("DETECTED", "Pixel was part of a detected source."),
1375 "SUSPECT": MaskPlane("SUSPECT", "Pixel was close to the saturation level. "),
1376 "NO_DATA": MaskPlane("NO_DATA", "No data was available for this pixel."),
1377 "VIGNETTED": MaskPlane("VIGNETTED", "Pixel was vignetted by the optics."),
1378 "PARTLY_VIGNETTED": MaskPlane("PARTLY_VIGNETTED", "Pixel was partly vignetted by the optics."),
1379 "CROSSTALK": MaskPlane("CROSSTALK", "Pixel was affected by crosstalk and corrected accordingly."),
1380 "ITL_DIP": MaskPlane(
1381 "ITL_DIP", "Pixel was affected by a dark vertical trail from a bright source, on an ITL CCD."
1382 ),
1383 "NOT_DEBLENDED": MaskPlane(
1384 "NOT_DEBLENDED",
1385 "Pixel belonged to a detection that was not deblended, usually due to size limits.",
1386 ),
1387 "SPIKE": MaskPlane(
1388 "SPIKE", "Pixel is in the neighborhood of a diffraction spike from a bright star."
1389 ),
1390 "UNMASKEDNAN": MaskPlane("UNMASKED_NAN", "Pixel was found to be NaN unexpectedly."),
1391 }
1394def get_legacy_difference_image_mask_planes() -> dict[str, MaskPlane]:
1395 """Return a mapping from legacy mask plane name to `MaskPlane` instance
1396 for LSST difference images, c. DP2.
1397 """
1398 result = get_legacy_visit_image_mask_planes()
1399 result["DETECTED_NEGATIVE"] = MaskPlane(
1400 "DETECTED_NEGATIVE", "Pixel was part of a detected source with negative flux."
1401 )
1402 result["SAT_TEMPLATE"] = MaskPlane("SAT_TEMPLATE", "Template pixel was saturated.")
1403 result["HIGH_VARIANCE"] = MaskPlane(
1404 "HIGH_VARIANCE", "Template pixel had fewer-than-usual input epochs and hence high noise."
1405 )
1406 result["STREAK"] = MaskPlane(
1407 "STREAK", "An extended streak (probably an artificial satellite) affected this pixel."
1408 )
1409 return result
1412def get_legacy_deep_coadd_mask_planes() -> dict[str, MaskPlane]:
1413 """Return a mapping from legacy mask plane name to `MaskPlane` instance
1414 for LSST deep coadds, c. DP2.
1415 """
1416 return {
1417 "NO_DATA": MaskPlane("NO_DATA", "No data was available for this pixel."),
1418 "INTRP": MaskPlane("INTERPOLATED", "Pixel value is the result of interpolating nearby good pixels."),
1419 "CR": MaskPlane(
1420 "COSMIC_RAY",
1421 "A cosmic ray affected this pixel on at least one input image (and was interpolated).",
1422 ),
1423 "SAT": MaskPlane(
1424 "SATURATED",
1425 "More than 10% of the potential input visits had a saturated pixel at this location "
1426 "('potential' because saturated pixel values are not actually propagated to the coadd). "
1427 "SATURATED always implies REJECTED, and is often a reason for NO_DATA.",
1428 ),
1429 "EDGE": MaskPlane(
1430 "DETECTION_EDGE",
1431 "Pixel was too close to the edge of the patch to be considered for detection, "
1432 "due to the finite size of the detection kernel.",
1433 ),
1434 "CLIPPED": MaskPlane(
1435 "CLIPPED",
1436 "Region was identified as a probable artifact when comparing multiple single-visit warps. "
1437 "CLIPPED always implies REJECTED.",
1438 ),
1439 "REJECTED": MaskPlane(
1440 "REJECTED",
1441 "At least one input visit was left out of the coadd for this pixel due to masking. "
1442 "REJECTED always implies INEXACT_PSF.",
1443 ),
1444 "DETECTED": MaskPlane("DETECTED", "Pixel was part of a detected source."),
1445 "INEXACT_PSF": MaskPlane(
1446 "INEXACT_PSF",
1447 "The set of visits contributing to this pixel differs from the set of visits "
1448 "contributing to the PSF model for its cell.",
1449 ),
1450 }
1453def get_legacy_non_cell_coadd_mask_planes() -> dict[str, MaskPlane]:
1454 """Return a mapping from legacy mask plane name to `MaskPlane` instance
1455 for LSST non-cell coadds such as ``template_coadd`` in DP2, and all
1456 DP1 coadds.
1458 These coadds carry the visit-level planes propagated from their input
1459 warps in addition to the coadd-specific planes, and flag chip edges with
1460 ``SENSOR_EDGE`` (cell coadds use ``CELL_EDGE`` instead).
1461 """
1462 result = get_legacy_deep_coadd_mask_planes()
1463 result["BAD"] = MaskPlane("BAD", "Bad pixel in the instrument, including bad amplifiers.")
1464 result["SUSPECT"] = MaskPlane("SUSPECT", "Pixel was close to the saturation level.")
1465 result["CROSSTALK"] = MaskPlane("CROSSTALK", "Pixel was affected by crosstalk and corrected accordingly.")
1466 result["DETECTED_NEGATIVE"] = MaskPlane(
1467 "DETECTED_NEGATIVE", "Pixel was part of a detected source with negative flux."
1468 )
1469 result["NOT_DEBLENDED"] = MaskPlane(
1470 "NOT_DEBLENDED",
1471 "Pixel belonged to a detection that was not deblended, usually due to size limits.",
1472 )
1473 result["UNMASKEDNAN"] = MaskPlane("UNMASKED_NAN", "Pixel was found to be NaN unexpectedly.")
1474 result["SENSOR_EDGE"] = MaskPlane(
1475 "SENSOR_EDGE",
1476 "Pixel is near the edge of a contributing sensor/chip, so the coadd PSF is discontinuous there.",
1477 )
1478 return result
1481def _guess_legacy_plane_map(old_planes: Mapping[str, int]) -> dict[str, MaskPlane]:
1482 """Guess which of the ``get_legacy_*_plane_map`` created the given mask
1483 plane dictionary and call it.
1484 """
1485 if "SAT_TEMPLATE" in old_planes: 1485 ↛ 1486line 1485 didn't jump to line 1486 because the condition on line 1485 was never true
1486 return get_legacy_difference_image_mask_planes()
1487 if "INEXACT_PSF" in old_planes:
1488 # Both cell and non-cell coadds have INEXACT_PSF, but only non-cell
1489 # (assemble_coadd) coadds flag chip edges with SENSOR_EDGE; cell coadds
1490 # use CELL_EDGE.
1491 if "SENSOR_EDGE" in old_planes:
1492 return get_legacy_non_cell_coadd_mask_planes()
1493 return get_legacy_deep_coadd_mask_planes()
1494 return get_legacy_visit_image_mask_planes()
1497def _reindex_legacy_plane_cards(
1498 header: astropy.io.fits.Header,
1499 old_planes: Mapping[str, int],
1500 plane_map: Mapping[str, MaskPlane],
1501 schema: MaskSchema,
1502) -> None:
1503 """Rewrite retained legacy ``MP_`` cards in place to match a reshuffled
1504 schema.
1506 Parameters
1507 ----------
1508 header
1509 Header whose ``MP_`` cards are updated in place.
1510 old_planes
1511 Mapping from legacy mask plane name to its original (on-disk) bit
1512 index, as returned by `MaskPlane.read_legacy`.
1513 plane_map
1514 Mapping from legacy mask plane name to the `MaskPlane` it was remapped
1515 to in ``schema``.
1516 schema
1517 The reconstructed schema that defines the new bit positions.
1519 Notes
1520 -----
1521 Each ``MP_<legacy name>`` card is set to the index that its remapped plane
1522 occupies in ``schema`` (equivalently, the ``MSKN`` index written on
1523 serialization). Cards for legacy planes that are not represented in the
1524 new schema are removed, since they no longer correspond to any stored bit.
1525 Legacy masks have at most 31 planes, so every plane maps to a single bit in
1526 one on-disk element and the index is unambiguous.
1527 """
1528 new_index = {plane.name: n for n, plane in enumerate(schema) if plane is not None}
1529 for old_name in old_planes:
1530 keyword = f"MP_{old_name}"
1531 new_plane = plane_map.get(old_name)
1532 if new_plane is not None and (index := new_index.get(new_plane.name)) is not None:
1533 header[keyword] = index
1534 else:
1535 del header[keyword]
1538def _strip_legacy_plane_cards(header: astropy.io.fits.Header) -> None:
1539 """Remove all legacy ``MP_*`` mask-plane cards from a FITS header.
1541 These are written only so that legacy tooling can read masks reconstructed
1542 from legacy cutouts; the ``lsst.images`` reader uses the serialized schema
1543 instead, so it strips them rather than carrying them as opaque metadata.
1544 """
1545 for keyword in [card.keyword for card in header.cards if card.keyword.startswith("MP_")]:
1546 del header[keyword]