import logging
import numpy as np
from stcal.outlier_detection.utils import compute_weight_threshold
from jwst import datamodels as dm
from jwst.outlier_detection._fileio import save_median
from jwst.outlier_detection.utils import flag_model_crs, nanmedian3D
from jwst.resample.resample_utils import build_mask
log = logging.getLogger(__name__)
__all__ = ["detect_outliers"]
[docs]
def detect_outliers(
input_model,
save_intermediate_results,
good_bits,
maskpt,
rolling_window_width,
snr,
make_output_path,
):
"""
Flag outliers in tso data.
Parameters
----------
input_model : `~stdatamodels.jwst.datamodels.CubeModel`
The input cube model.
save_intermediate_results : bool
If `True`, save the rolling median model as a
`~stdatamodels.jwst.datamodels.CubeModel`.
good_bits : int
DQ flag bit values indicating good pixels.
maskpt : float
The percentage of the mean weight to use as a threshold for masking.
rolling_window_width : int
The width of the rolling median window.
snr : float
The signal-to-noise ratio threshold for flagging outliers.
make_output_path : callable
A function that generates a path for saving intermediate results.
Returns
-------
`~stdatamodels.jwst.datamodels.CubeModel`
The input model with outliers flagged.
"""
if isinstance(input_model, dm.ModelContainer):
raise TypeError("OutlierDetectionTSO does not support ModelContainer input.")
weighted_cube = weight_no_resample(input_model, good_bits)
weight_threshold = compute_weight_threshold(weighted_cube.wht, maskpt)
if (rolling_window_width > 1) and (rolling_window_width < weighted_cube.shape[0]):
medians = compute_rolling_median(weighted_cube, weight_threshold, w=rolling_window_width)
else:
medians = nanmedian3D(weighted_cube.data, overwrite_input=False)
# this is a 2-D array, need to repeat it into the time axis
# for consistent shape with rolling median case
medians = np.broadcast_to(medians, weighted_cube.shape)
# Save median model if pars['save_intermediate_results'] is True
# this will be a CubeModel with rolling median values.
if save_intermediate_results:
median_model = dm.CubeModel(data=medians) # type: ignore[name-defined]
with dm.open(weighted_cube) as dm0:
median_model.update(dm0)
save_median(median_model, make_output_path)
del median_model
# no need for blotting, resample is turned off for TSO
# go straight to outlier detection
log.info("Flagging outliers")
flag_model_crs(
input_model,
medians,
snr,
)
return input_model
def weight_no_resample(input_model, good_bits):
"""
Give weights to model without resampling.
Parameters
----------
input_model : `~stdatamodels.jwst.datamodels.CubeModel`
The input cube model.
good_bits : int
DQ flag bit values indicating good pixels.
Returns
-------
`~stdatamodels.jwst.datamodels.CubeModel`
A copy of the input cube model with weights assigned in the ``wht`` extension.
"""
# Notes
# -----
# Prior to PR #8473, the `build_driz_weight` function was used to
# create the weights for the input models for TSO data. However, that
# function was simply returning a copy of the DQ array because the
# var_noise was not being passed in by calwebb_tso3. As of PR #8473,
# a cube model that includes the var_noise is passed into TSO
# outlier detection, so `build_driz_weight` would weight the cube model
# by the variance. Therefore `build_driz_weight` was removed in order to
# preserve the original behavior. If it is determined later that exposure
# time or inverse variance weighting should be used here, build_driz_weight
# should be re-implemented.
weighted_cube = input_model.copy()
dqmask = build_mask(input_model.dq, good_bits)
weighted_cube.wht = dqmask.astype(np.float32)
return weighted_cube
def compute_rolling_median(
model: dm.CubeModel, # type: ignore[name-defined]
weight_threshold: np.ndarray,
w: int = 25,
) -> np.ndarray:
"""
Set bad and low-weight data to NaN, then compute the rolling median over the time axis.
Parameters
----------
model : `~stdatamodels.jwst.datamodels.CubeModel`
The input cube model.
weight_threshold : ndarray
The weight thresholds for each integration.
w : int
The window size for the rolling median.
Returns
-------
ndarray
The rolling median of the input data. Same dimensions as input.
"""
sci = model.data
weight = model.wht
badmask = np.less(weight, weight_threshold)
log.debug(f"Percentage of pixels with low weight: {np.sum(badmask) / weight.size * 100}")
# Fill resampled_sci array with nan's where mask values are True
sci[badmask] = np.nan
del badmask
if w > sci.shape[0]:
raise ValueError("Window size must be less than the number of integrations.")
meds = moving_median_over_zeroth_axis(sci, w)
del sci
return meds
def moving_median_over_zeroth_axis(x: np.ndarray, w: int) -> np.ndarray:
"""
Calculate the median of a moving window over the zeroth axis of an N-d array.
Slide a window of size ``w`` over the array along axis 0, and for each position,
calculate the median of the values inside that window (across axis 0 only).
The result at each step is stored in the center position of the window,
producing an output array with the same shape as the input.
Because the window cannot fully overlap the data at the beginning and end,
those edge positions are filled with the nearest computed median value to
avoid missing data.
Parameters
----------
x : ndarray
The input array.
w : int
The window size.
Returns
-------
ndarray
The rolling median of the input array. Same dimensions as input.
"""
if w <= 1:
raise ValueError("Rolling median window size must be greater than 1.")
out = np.full(x.shape, np.nan)
hw, odd_window = divmod(w, 2)
for start_index in range(x.shape[0] - w + 1):
end_index = start_index + w
np.median(x[start_index:end_index], axis=0, out=out[start_index + hw])
# Fill in the edges with the nearest valid value
out[:hw] = out[hw]
if odd_window:
out[-hw:] = out[-hw - 1]
else:
out[-hw + 1 :] = out[-hw]
return out