"""Class for standardizing Data Products of Vera C. Rubin Science Pipelines
via the Rubin Data Butler.
"""
import copy
import importlib
import logging
import uuid
from astropy.coordinates import SkyCoord
import astropy.nddata as bitmask
import astropy.time
from astropy.wcs import WCS
from astropy.wcs.utils import fit_wcs_from_points
import astropy.units as u
from networkx import center
import numpy as np
from .standardizer import Standardizer, StandardizerConfig
from kbmod.wcs_utils import max_sky_separation
from kbmod.core.psf import PSF
from kbmod.core.image_stack_py import LayeredImagePy
# Set up logger for this module
logger = logging.getLogger(__name__)
__all__ = [
"ButlerStandardizer",
"ButlerStandardizerConfig",
]
def deferred_import(module, name=None):
"""Imports module/class/function/name as ``name`` into global modules.
If module/class/function/name already exists does nothing.
Intended for importing a large module or a library only when needed, as to
avoid the cost of the import if the functionality depending on the imported
library is not being used.
Used to defer the import of Vera C. Rubin Middleware component only during
`ButlerStandardizer` initialization so that we may be able to import KBMOD
before the heat death of the universe.
Parameters
----------
module : `str`
Fully qualified name of the module/submodule/class/function/name.
name : `str` or `None`, optional
Name the target is imported as.
Raises
------
ImportError :
Target is not reachable. Ensure package is
installed and visible int the environment.
"""
if name is None:
name = module
if globals().get(name, False):
return
try:
globals()[name] = importlib.import_module(module)
except ImportError as e:
raise ImportError(f"Module {module} or name {name} not found.") from e
[docs]class ButlerStandardizerConfig(StandardizerConfig):
do_mask = True
"""Perform masking if ``True``, otherwise return an empty mask."""
do_bitmask = True
"""Mask ``mask_flags`` from the mask plane in the FITS file."""
do_threshold = False
"""Mask all pixels above the given count threshold."""
grow_mask = True
"""Grow mask footprint by ``grow_kernel_shape``"""
brightness_threshold = 10
"""Pixels with value greater than this threshold will be masked."""
grow_kernel_shape = (10, 10)
"""Size of the symmetric square kernel by which mask footprints will be
increased by."""
mask_flags = ["BAD", "CLIPPED", "CR", "CROSSTALK", "EDGE", "NO_DATA", "SAT", "SENSOR_EDGE", "SUSPECT"]
"""List of flags that will be masked."""
psf_std = 1
"""Standard deviation of the Point Spread Function. Used as a fallback
when ``psf_std_from_summary`` is disabled, or when the measured value is
unavailable for an image."""
psf_std_from_summary = True
"""Build each image's PSF kernel from the measured ``psfSigma`` in the
exposure summary statistics. Set False to use the fixed ``psf_std``."""
standardize_metadata = True
"""Fetch and include values from Rubin's Exposure.metadata
PropertyList such as 'OBSID', 'AIRMASS', 'DIMM2SEE' ec. Typically
corresponds to ingested header values, ergo the capitalized column
names.
"""
standardize_effective_summary_stats = False
"""Include the "effective" fit metrics from SummaryStats"""
standardize_uri = False
"""Include an URL-like path to the file"""
wcs_fallback_points = 1000
"""Number of random points to sample across the detector when
an astropy WCS cannot be constructed from the Rubin SkyWCS metadata."""
wcs_fallback_sips_degree = 4
"""Degree of the SIP distortion to fit when creating a fallback WCS when
an astropy WCS cannot be constructed from the Rubin SkyWCS metadata.
If ``None``, no SIP distortion is fitted."""
zero_point = 31
"""Target photometric zero point for science and variance. The input
calibration comes from the Exposure's attached PhotoCalib, not its summary
statistics. This is a unit conversion, not a source-color correction."""
greedy_export = False
"""If True, the standardizer will keep the Exposure object in memory
after the LayeredImage is created. This is useful for large datasets."""
[docs]class ButlerStandardizer(Standardizer):
"""Standardizer for Vera C. Rubin Data Products, namely the underlying
``Exposure`` objects of ``calexp``, ``difim`` and various ``*coadd``
dataset types.
This standardizer will volunteer to process ``DatasetRef`` or ``DatasetId``
objects. The Rubin Data Butler is expected to be provided at
instantiation time, in addition to the target we want to standardize.
Parameters
----------
dataId : `lsst.daf.butler.core.DatasetId`, `lsst.daf.butler.core.DatasetRef` or `int`
Target to standardize.
butler : `lsst.daf.butler.Butler` or `list[lsst.daf.butler.Butler]`
Vera C. Rubin Data Butler or a list of butlers. The butlers are queried
to resolve the given target in the given order.
config : `StandardizerConfig`, `dict` or `None`, optional
Configuration key-values used when standardizing the file.
Attributes
----------
butler : `lsst.daf.butler.Butler`
Vera C. Rubin Data Butler.
ref : `lsst.daf.butler.core.DatasetRef`
Dataset reference to the given target
exp : `lsst.afw.image.exposure.Exposure`
The `Exposure` object targeted by the ``ref``
processable : `list[lsst.afw.image.exposure.Exposure]`
Items marked as processable by the standardizer. See `Standardizer`.
"""
name = "ButlerStandardizer"
priority = 2
configClass = ButlerStandardizerConfig
[docs] @classmethod
def resolveTarget(self, tgt):
# DatasetId is a type alias for UUID's so it'll be like a large int or
# hex or string of int/hex etc. We try to cast to UUID to check if str
# is compliant
# https://github.com/lsst/daf_butler/blob/main/python/lsst/daf/butler/_dataset_ref.py#L265
if isinstance(tgt, uuid.UUID):
return True
if isinstance(tgt, str):
try:
uuid.UUID(tgt)
except ValueError:
return False
else:
return True
# kinda hacky, but I don't want to import the entire Stack before I
# absolutely know I need/have it.
tgttype = str(type(tgt)).lower()
if "datasetref" in tgttype or "datasetid" in tgttype:
return True
return False
@classmethod
def __query_butler(self, tgt, butler):
"""Given a target and a butler, which might not contain the target
queries the butler to resolve it. Butler failures are silenced.
Has to be called after deferred_import.
Parameters
----------
tgt : `lsst.daf.butler.core.DatasetId`, `lsst.daf.butler.core.DatasetRef` or `int`
Target to standardize.
butler : `lsst.daf.butler.Butler` or `list[lsst.daf.butler.Butler]`
Vera C. Rubin Data Butler or a list of butlers. The butlers are queried
to resolve the given target in the given order.
Raises
------
TypeError : When given target is not a DatasetRef, DatasetId, or unique integer ID"
"""
# including records expands the dataId to include
# key pieces of information such as filter and band
# loading datastore_records could be a shortcut to
# relative path inside the repository
if isinstance(tgt, dafButler.DatasetRef):
ref = tgt
elif isinstance(tgt, dafButler.DatasetId):
ref = butler.get_dataset(tgt, dimension_records=True)
elif isinstance(tgt, (uuid.UUID, str)):
did = dafButler.DatasetId(tgt)
ref = butler.get_dataset(did, dimension_records=True)
else:
raise TypeError("Expected DatasetRef, DatasetId or an unique integer ID, " f"got {tgt} instead.")
return ref, butler
def __init__(self, dataId, butler, config=None, **kwargs):
deferred_import("lsst.daf.butler", "dafButler")
# Sometimes we find ourselves having to process data that is
# in the process of migration between multiple repositories.
# We want to prioritize one of these repos as the preffered
# source of data, but it does not yet contain all data.
# So we check all the given butlers to resolve a target in
# order and then skip once we get a hit. To cover the more
# the plain single-butler use case just promote it to a list.
if isinstance(butler, dafButler.Butler):
butlers = [
butler,
]
else:
butlers = butler
for b in butlers:
self.ref, self.butler = self.__query_butler(dataId, b)
if self.ref is not None:
break
if self.ref is None:
raise ValueError(f"Unable to resolve target {dataId} for any butler.")
# Now that target was upgraded to a ref and the correct butler
# is know we can get the info we need from them.
# Somewhere around w_2024_ builds the datastore.root
# was removed as an attribute of the datastore, not sure
# it was ever replaced with anything back-compatible. We simply
# check for the which _datastore attribute is available for this
# butler and then check wherther it has a root or roots attribute.
if hasattr(self.butler, "datastore"):
datastore_root = self.butler.datastore.root
elif hasattr(butler, "_datastore"):
if hasattr(self.butler._datastore, "root"):
datastore_root = self.butler._datastore.root
elif hasattr(butler._datastore, "roots"):
datastore_root = self.butler._datastore.roots
else:
raise AttributeError("Butler does not have a valid datastore root attribute.")
else:
raise AttributeError("Butler does not have a valid datastore attribute.")
super().__init__(str(datastore_root), config=config)
self.exp = None
self.processable = [self.exp]
# This is for lazy loading
self._bbox = None
self._wcs = None
# This is set for internal use because often we can't share
# between methods and they all need these values. Centralizing
# metadatada fetching via butler also helps performance
self._metadata = None
self._naxis1 = None
self._naxis2 = None
def __deepcopy__(self, memo):
# Share self.butler across the copy. DirectButler pickles to a config
# and unpickles by calling create_from_config, which opens a fresh
# psycopg2 connection to the registry. Deep-copying many standardizers
# (e.g. in kbmod.injection.inject_sources_into_ic) would otherwise fan
# out into one new DB connection per image and overwhelm the registry
# pooler under parallel workers.
cls = self.__class__
new = cls.__new__(cls)
memo[id(self)] = new
for k, v in self.__dict__.items():
if k == "butler":
new.butler = v
else:
setattr(new, k, copy.deepcopy(v, memo))
return new
def _fitWCSFallback(self, wcs, naxis1, naxis2, n_rand_pts, sip_degree, sample_outside_chip=True):
"""Create a simple TAN WCS centered on the detector through sampling random points.
Parameters
----------
wcs : `SkyWCS`
Rubin SkyWCS object.
naxis1: `int`
Detector width in pixels (FITS NAXIS1, the WCS x-axis).
naxis2: `int`
Detector height in pixels (FITS NAXIS2, the WCS y-axis).
n_rand_pts : `int`
Number of random points to sample across the detector.
sip_degree : `int` or `None`, optional
Degree of the SIP distortion to fit. If ``None``, no SIP distortion
is fitted.
sample_outside_chip: `bool`, optional
Whether to sample points outside of the chip bounds to interpolate rather
than extrapolate the WCS fit. Default is `True`.
Returns
-------
wcs : `astropy.wcs.WCS`
Fitted WCS object.
"""
if n_rand_pts <= 0:
raise ValueError("Number of random points must be positive.")
if sip_degree is not None and sip_degree <= 0:
raise ValueError("SIP degree must be non-negative or None.")
if not sample_outside_chip:
# Sample random X, Y points across this detector
rand_xy = np.random.rand(n_rand_pts, 2) * [naxis1, naxis2]
rand_x, rand_y = rand_xy[:, 0], rand_xy[:, 1]
else:
# Expand our X, Y grid slightly beyond the bounds of the detector.
grid_offset_percent = 0.1
rand_xy = np.random.rand(n_rand_pts, 2) * np.array(
[[naxis1 * (1 + 2 * grid_offset_percent), naxis2 * (1 + 2 * grid_offset_percent)]]
)
# Subtract our offset from each dimension to ensure we do not begin sampling at "0"
rand_x = rand_xy[:, 0] - (naxis1 * grid_offset_percent)
rand_y = rand_xy[:, 1] - (naxis2 * grid_offset_percent)
# Turn our random X, Y points into an RA, Dec SkyCoord
rand_ra, rand_dec = wcs.pixelToSkyArray(rand_x, rand_y, degrees=True)
world_coords = SkyCoord(ra=rand_ra * u.deg, dec=rand_dec * u.deg, frame="icrs")
# For our center point, we just use the center of the detector
detector_center = wcs.pixelToSky(int(naxis1 // 2), int(naxis2 // 2))
detector_center_coord = SkyCoord(
ra=detector_center.getRa().asDegrees() * u.deg,
dec=detector_center.getDec().asDegrees() * u.deg,
frame="icrs",
)
# Fit a TAN WCS to these points, with optional SIP distortion
fitted_wcs = fit_wcs_from_points(
(rand_x, rand_y), world_coords, detector_center_coord, sip_degree=sip_degree
)
# fit_wcs_from_points sets pixel_shape from the sampled points' bounding
# box, which extends beyond the chip; pin it to the true (width, height).
# Astropy's array_shape is the reverse: (height, width).
fitted_wcs.pixel_shape = (naxis1, naxis2)
return fitted_wcs
def _computeSkyBBox(self, wcs, height, width):
"""Given an Rubin SkyWCS object and the dimensions of an image
calculates the values of world coordinates image center and
image corners.
The corners are given by the following indices:
topleft topright
(0, height) ------------- (width, height)
| |
| x |
| (width/2, height/2) |
| center |
| |
(0, 0) ----------------- (width, 0)
botleft botright
Parameters
----------
wcs : `object`
World coordinate system object, must support standard WCS API.
height : `int`
Image height (number of rows).
width : `int`
Image width (number of columns).
return_type : `str`, optional
A 'dict' or an 'array', the type the result is returned as.
Returns
-------
standardizedBBox : `dict` or `array`
An array of shape ``(5, 2)`` starting with center coordinate,
then bottom left and progressing clockwise around the detector.
When a dictionary, ``ra`` and ``dec`` keys mark center coordinates.
The ``'ra_bl', 'dec_bl'`` mark bottom left. Progressing clocwise
again, ``tl``, ``tr`` and ``br`` mark top left, top right and bottom
right mark the edge position, and ``ra_`` and ``dec_`` prefix the
coordinate.
Notes
-----
Dimensions are passed in (height, width) order. WCS pixel
coordinates use (x, y), so the center is (width/2, height/2), rounded down.
Bottom left corner is taken to be the (0,0)-th pixel and image lies
in the first quadrant of a unit circle to match Astropy's convention.
"""
center = wcs.pixelToSky(int(width // 2), int(height // 2))
botleft = wcs.pixelToSky(0, 0)
topleft = wcs.pixelToSky(0, height)
topright = wcs.pixelToSky(width, height)
botright = wcs.pixelToSky(width, 0)
pts = np.array(
[
[center.getRa().asDegrees(), center.getDec().asDegrees()],
[botleft.getRa().asDegrees(), botleft.getDec().asDegrees()],
[topleft.getRa().asDegrees(), topleft.getDec().asDegrees()],
[topright.getRa().asDegrees(), topright.getDec().asDegrees()],
[botright.getRa().asDegrees(), botright.getDec().asDegrees()],
]
)
return pts
@staticmethod
def _visit_midpoint_utc(visit):
"""Return the native VisitInfo midpoint as an Astropy UTC Time.
VisitInfo.date is already the exposure midpoint. Keep the Time object
until serialization so nominal-start arithmetic retains its precision.
"""
return visit.date.toAstropy().utc
@staticmethod
def _mjd_to_obs_day(mjd_mid):
"""Convert UTC MJD to observing day in YYYYMMDD format.
Uses the Rubin/DECam noon-TAI boundary to group an observing night;
this is a fixed convention, not a local sunset calculation.
Parameters
----------
mjd_mid : `float`
UTC Modified Julian Date at the middle of the exposure.
Returns
-------
obs_day : `int`
Observing day in YYYYMMDD format.
"""
observing_date = astropy.time.Time(mjd_mid, format="mjd", scale="utc").tai
offset = astropy.time.TimeDelta(12 * 3600, format="sec", scale="tai")
observing_date -= offset
return int(observing_date.strftime("%Y%m%d"))
def _fetch_meta(self):
"""Fetch metadata and any dataset components that do not
load the image or large amount of data.
This resolves the majority of the metadata, temporal and
spatial information required by KBMOD, which is stored
in the attributes and then returned in the standardize
mehtods. This prevents having to evaluate queries to the
Butler Registry and Datastore repeatedly, and enables the
evalution of information that spans multiple datasets without
forcing the Butler.get of a larger dataset (f.e. the Exposure
object).
"""
self._metadata = {}
# First we standardize the required metadata. Most of this can
# be extracted from the dataset reference and a visitInfo.
# This includes i) data required to roundtrip the standardizer
# ii) timestamp and iii) pointing information
self._metadata["dataId"] = str(self.ref.id)
self._metadata["collection"] = self.ref.run
self._metadata["datasetType"] = self.ref.datasetType.name
self._metadata["visit"] = self.ref.dataId["visit"]
self._metadata["detector"] = self.ref.dataId["detector"]
self._metadata["band"] = self.ref.dataId["band"]
self._metadata["filter"] = self.ref.dataId["physical_filter"]
visit_ref = self.ref.makeComponentRef("visitInfo")
visit = self.butler.get(visit_ref)
expt = visit.exposureTime
# VisitInfo.date already represents the exposure midpoint, in TAI.
# Derive the nominal start from the duration; do not advance the midpoint.
midpoint = self._visit_midpoint_utc(visit)
mjd_start = midpoint - (expt / 2) * u.s
self._metadata["exposureTime"] = expt
# Note the timescales for MJD. The Butler uses TAI, but we convert
# time stamps to UTC for consistency.
# Name mjd into mjd_mid - make it obvious it's middle of exposure.
self._metadata["mjd_start"] = mjd_start.utc.mjd
self._metadata["mjd_mid"] = midpoint.mjd
self._metadata["obs_day"] = ButlerStandardizer._mjd_to_obs_day(midpoint.mjd)
self._metadata["object"] = visit.object
self._metadata["pointing_ra"] = visit.boresightRaDec.getRa().asDegrees()
self._metadata["pointing_dec"] = visit.boresightRaDec.getDec().asDegrees()
self._metadata["airmass"] = visit.boresightAirmass
obs = visit.getObservatory()
self._metadata["obs_lon"] = obs.getLongitude().asDegrees()
self._metadata["obs_lat"] = obs.getLatitude().asDegrees()
self._metadata["obs_elev"] = obs.getElevation()
# Pointing information is hard to standardize because the
# dimensions of the detector are not easily availible. We get
# those from BBox (in-pixel bounding box). NAXIS values are
# required if we reproject, so we must extract them if we can
bbox_ref = self.ref.makeComponentRef("bbox")
bbox = self.butler.get(bbox_ref)
self._naxis1 = bbox.getWidth()
self._naxis2 = bbox.getHeight()
# If the standardizer is re-used, many generators will be
# depleted, returning None as values. Cast to dict to make
# a copy.
wcs_ref = self.ref.makeComponentRef("wcs")
wcs = self.butler.get(wcs_ref)
# Always resample the WCS via SIP fit to the Rubin SkyWCS; the lossy
# getFitsMetadata() TAN-SIP export is not used. _fitWCSFallback sets
# pixel_shape=(naxis1,naxis2), so serialized NAXIS is correct.
n_rand_pts = self.config["wcs_fallback_points"]
sip_degree = self.config["wcs_fallback_sips_degree"]
self._wcs = self._fitWCSFallback(wcs, self._naxis1, self._naxis2, n_rand_pts, sip_degree)
center_pt = bbox.getCenter()
self._metadata["pixel_scale"] = wcs.getPixelScale(center_pt).asArcseconds()
# calculate the WCS "error" (max on-sky separation between edge
# coordinates from Rubin's more powerful SkyWCS and Astropy's
# Fits-WCS), in degrees
# Both bbox helpers take (height, width), unlike WCS.pixel_shape.
skyBBox = self._computeSkyBBox(wcs, self._naxis2, self._naxis1)
apyBBox = self._computeBBoxArray(self._wcs, self._naxis2, self._naxis1)
self._metadata["wcs_err"] = max_sky_separation(skyBBox, apyBBox)
# TODO: see issue #666
# this will unroll the entire bbox into columns
# make sky bbox the default, since that one is guaranteed to be correct
self._bbox = self._bboxArrayToDict(skyBBox)
self._metadata.update(self._bbox)
# Preserve the summary statistics as upstream provenance. In particular,
# zeroPoint can describe instrumental counts BEFORE Rubin calibrated the
# delivered pixels to nJy; it must not be used to scale those pixels.
summary_ref = self.ref.makeComponentRef("summaryStats")
summary = self.butler.get(summary_ref)
self._metadata["psfSigma"] = summary.psfSigma
self._metadata["psfArea"] = summary.psfArea
self._metadata["nPsfStar"] = summary.nPsfStar
self._metadata["zeroPoint"] = summary.zeroPoint
self._metadata["skyBg"] = summary.skyBg
self._metadata["skyNoise"] = summary.skyNoise
self._metadata["meanVar"] = summary.meanVar
self._metadata["astromOffsetMean"] = summary.astromOffsetMean
self._metadata["astromOffsetStd"] = summary.astromOffsetStd
# The rest of the data here is optional, generally metadata
# is nice to standardize, but keys may change between
# different instruments, summary stats are useful for
# photometric analysis of the results, while the effective
# values are too often NaN. The URI location itself is
# ultimately not very useful, but helpful for data inspection.
if self.config["standardize_metadata"]:
meta_ref = self.ref.makeComponentRef("metadata")
meta = self.butler.get(meta_ref)
# dataId sometimes doesn't have a filter or a band,
# depending on the way the initial ref is resolved? Why
# is middleware so complicated! Best-effort attempt,
# 90% cases?
self._metadata["OBSID"] = meta["OBSID"]
# Note that the following metadata keys may not be present in all
# Rubin butlers, and we only extract them if available.
if "DTNSANAM" in meta:
self._metadata["DTNSANAM"] = meta["DTNSANAM"]
if "AIRMASS" in meta:
self._metadata["AIRMASS"] = meta["AIRMASS"]
d2s = 0.0
if "DIMM2SEE" in meta and meta["DIMM2SEE"] != "NaN":
self._metadata["DIMM2SEE"] = d2s
if "GAINA" in meta:
self._metadata["GAINA"] = meta["GAINA"]
if "GAINB" in meta:
self._metadata["GAINB"] = meta["GAINB"]
# Will be nan for VR filter so it's optional
if self.config["standardize_effective_summary_stats"]:
self._metadata["effTime"] = summary.effTime
self._metadata["effTimePsfSigmaScale"] = summary.effTimePsfSigmaScale
self._metadata["effTimeSkyBgScale"] = summary.effTimeSkyBgScale
self._metadata["effTimeZeroPointScale"] = summary.effTimeZeroPointScale
if self.config["standardize_uri"]:
self._metadata["location"] = self.butler.getURI(
self.ref,
collections=[
self.ref.run,
],
).geturl()
else:
# Save the full string representation of ref.
self._metadata["location"] = str(self.ref)
@property
def wcs(self):
if self._wcs is None:
self._fetch_meta()
return [
self._wcs,
]
@property
def bbox(self):
if self._bbox is None:
self._fetch_meta()
return [
self._bbox,
]
def _get_photometric_scale(self):
"""Return the scalar flux conversion from delivered pixels to our zero point.
As before, this standardizer uses an exposure-wide calibration; calling
PhotoCalib without a position uses its mean calibration. Applying a full
spatially varying calibration is outside this scalar conversion.
"""
self.exp = self.butler.get(self.ref) if self.exp is None else self.exp
photo_calib = self.exp.photoCalib
if photo_calib is None:
raise ValueError("Cannot standardize pixels without an attached PhotoCalib.")
# PhotoCalib describes the pixels we actually have, including already
# calibrated visit/difference images. For nJy pixels it is the identity,
# and the magnitude of one pixel flux unit is 31.4 AB mag. The older
# summary zero point would apply the instrumental calibration AGAIN.
input_zero_point = float(photo_calib.instFluxToMagnitude(1.0))
# Read the mapping: attribute access would return the class default and
# silently ignore a configured target (for example, 31.4 to retain nJy).
target_zero_point = float(self.config["zero_point"])
if not np.isfinite(input_zero_point) or not np.isfinite(target_zero_point):
raise ValueError("Input PhotoCalib and target zero points must be finite.")
# m = ZP - 2.5 log10(flux). Thus nJy -> ZP 31 multiplies science by
# 10**((31 - 31.4)/2.5), and variance by the SQUARE of that factor.
# Do not cache this factor: injection can replace self.exp, and callers
# can change the target zero point between exports.
with np.errstate(over="ignore", under="ignore", invalid="ignore"):
scale = np.power(10.0, (target_zero_point - input_zero_point) / 2.5)
variance_scale = scale * scale
# A finite, positive flux scale can still overflow/underflow when squared.
if not np.isfinite(variance_scale) or variance_scale <= 0:
raise ValueError("Photometric scale and its square must be finite and positive.")
# A Python scalar preserves float32 image dtype under NumPy 2 promotion.
return float(scale)
[docs] def standardizeScienceImage(self):
scale = self._get_photometric_scale()
return [
self.exp.image.array * scale,
]
[docs] def standardizeVarianceImage(self):
scale = self._get_photometric_scale()
return [
self.exp.variance.array * scale**2,
]
[docs] def standardizeMaskImage(self):
self.exp = self.butler.get(self.ref) if self.exp is None else self.exp
if self._naxis1 is None or self._naxis2 is None:
self._fetch_meta()
# Return empty masks if no masking is done. We only return a single
# mask with NumPy shape (height, width), the reverse of WCS.pixel_shape.
if not self.config["do_mask"]:
return [np.zeros((self._naxis2, self._naxis1))]
# Otherwise load the mask extension and process it
mask = self.exp.mask.array.astype(int)
if self.config["do_bitmask"]:
# flip_bits makes ignore_flags into mask_these_flags
bit_flag_map = self.exp.mask.getMaskPlaneDict()
bit_flag_map = {key: int(2**val) for key, val in bit_flag_map.items()}
# Filter out any configured mask flags not present in this
# exposure's mask plane to avoid KeyError at lookup time.
available_flags = [f for f in self.config["mask_flags"] if f in bit_flag_map]
missing_flags = set(self.config["mask_flags"]) - set(available_flags)
if missing_flags:
logger.debug(f"Mask flags {missing_flags} not found in exposure mask plane, skipping.")
mask = bitmask.bitfield_to_boolean_mask(
bitfield=mask,
ignore_flags=available_flags,
flag_name_map=bit_flag_map,
flip_bits=True,
)
if self.config["do_threshold"]:
bmask = self.exp.image.array > self.config["brightness_threshold"]
mask = mask | bmask
if self.config["grow_mask"]:
# Only import the scipy module if we actually need it.
from scipy.signal import convolve2d
grow_kernel = np.ones(self.config["grow_kernel_shape"])
mask = convolve2d(mask, grow_kernel, mode="same").astype(bool)
return [
mask,
]
[docs] def standardizePSF(self):
# A user-supplied fallback value for when the measured PSF width is unavailable.
std = self.config["psf_std"]
if not self.config["psf_std_from_summary"]:
logger.debug(f"Using fallback psf_std: {std}.")
else:
# Use the measured PSF width from the summary statistics if available, otherwise fall back to the user-supplied value.
if self._metadata is None:
# Metadata is lazily fetched so explicitly check and fetch in case this is a fresh standardizer.
self._fetch_meta()
# Keep the raw value around so the warning below can report what the
# butler actually handed back, not the None the failed cast leaves.
raw = self._metadata.get("psfSigma", None)
try:
measured = float(raw)
except (TypeError, ValueError):
measured = None
if measured is None or not np.isfinite(measured) or measured <= 0:
logger.warning(f"Unusable psfSigma {raw!r}. Using fallback psf_std: {std}.")
else:
std = measured
return [PSF.make_gaussian_kernel(std)]
# These exist because standardizers promise to return lists
# for compatiblity for single-data and multi-data sources
[docs] def standardizeWCS(self):
if self._wcs is None:
self._fetch_meta()
return [
self._wcs,
]
[docs] def standardizeBBox(self):
if self._bbox is None:
self._fetch_meta()
return [
self._bbox,
]
[docs] def toLayeredImage(self):
masks = self.standardizeMaskImage()
mask = masks[0].astype(np.float32)
imgs = [
LayeredImagePy(
self.standardizeScienceImage()[0],
self.standardizeVarianceImage()[0],
mask=mask,
psf=self.standardizePSF()[0],
time=self._metadata["mjd_mid"],
),
]
if not self.config["greedy_export"]:
self.exp = None
return imgs