diff --git a/CHANGES.rst b/CHANGES.rst index f29a830a9..98ecbe356 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,6 +1,10 @@ 0.13.0 (unreleased) =================== +outlier_detection +----------------- +- Implemented ``outlier-detection step``. [#981] + associations ------------ diff --git a/docs/roman/references_general/dq_flags.inc b/docs/roman/references_general/dq_flags.inc index 6ab5b24e1..284a09f40 100644 --- a/docs/roman/references_general/dq_flags.inc +++ b/docs/roman/references_general/dq_flags.inc @@ -10,7 +10,7 @@ Bit Value Name Description 4 16 GW_AFFECTED_DATA Data affected by the GW read window 5 32 PERSISTENCE High persistence (was RESERVED_2) 6 64 AD_FLOOR Below A/D floor (0 DN, was RESERVED_3) -7 128 RESERVED_4 +7 128 OUTLIER Detected as outlier in coadded image 8 256 UNRELIABLE_ERROR Uncertainty exceeds quoted error 9 512 NON_SCIENCE Pixel not on science portion of detector 10 1024 DEAD Dead pixel diff --git a/romancal/datamodels/container.py b/romancal/datamodels/container.py index 4d6652ffe..65affce76 100644 --- a/romancal/datamodels/container.py +++ b/romancal/datamodels/container.py @@ -8,6 +8,7 @@ from collections.abc import Iterable, Sequence from pathlib import Path +import numpy as np from roman_datamodels import datamodels as rdm from romancal.lib.basic_utils import is_association @@ -18,6 +19,7 @@ "ModelContainer", ] +_ONE_MB = 1 << 20 RECOGNIZED_MEMBER_FIELDS = [ "tweakreg_catalog", ] @@ -529,6 +531,82 @@ def get_crds_parameters(self): return crds_header + def set_buffer(self, buffer_size, overlap=None): + """Set buffer size for scrolling section-by-section access. + + Parameters + ---------- + buffer_size : float, None + Define size of buffer in MB for each section. + If `None`, a default buffer size of 1MB will be used. + + overlap : int, optional + Define the number of rows of overlaps between sections. + If `None`, no overlap will be used. + """ + self.overlap = 0 if overlap is None else overlap + self.grow = 0 + + with rdm.open(self._models[0]) as model: + imrows, imcols = model.data.shape + data_item_size = model.data.itemsize + data_item_type = model.data.dtype + del model + + min_buffer_size = imcols * data_item_size + + self.buffer_size = ( + min_buffer_size if buffer_size is None else (buffer_size * _ONE_MB) + ) + + section_nrows = min(imrows, int(self.buffer_size // min_buffer_size)) + + if section_nrows == 0: + self.buffer_size = min_buffer_size + logger.warning( + "WARNING: Buffer size is too small to hold a single row." + f"Increasing buffer size to {self.buffer_size / _ONE_MB}MB" + ) + section_nrows = 1 + + nbr = section_nrows - self.overlap + nsec = (imrows - self.overlap) // nbr + if (imrows - self.overlap) % nbr > 0: + nsec += 1 + + self.n_sections = nsec + self.nbr = nbr + self.section_nrows = section_nrows + self.imrows = imrows + self.imcols = imcols + self.imtype = data_item_type + + def get_sections(self): + """Iterator to return the sections from all members of the container.""" + + for k in range(self.n_sections): + e1 = k * self.nbr + e2 = e1 + self.section_nrows + + if k == self.n_sections - 1: # last section + e2 = min(e2, self.imrows) + e1 = min(e1, e2 - self.overlap - 1) + + data_list = np.empty( + (len(self._models), e2 - e1, self.imcols), dtype=self.imtype + ) + wht_list = np.empty( + (len(self._models), e2 - e1, self.imcols), dtype=self.imtype + ) + for i, model in enumerate(self._models): + model = rdm.open(model, memmap=self._memmap) + + data_list[i, :, :] = model.data[e1:e2].copy() + wht_list[i, :, :] = model.weight[e1:e2].copy() + del model + + yield (data_list, wht_list, (e1, e2)) + def make_file_with_index(file_path, idx): """Append an index to a filename diff --git a/romancal/lib/dqflags.py b/romancal/lib/dqflags.py index 78f84da93..50b3f0837 100644 --- a/romancal/lib/dqflags.py +++ b/romancal/lib/dqflags.py @@ -30,7 +30,7 @@ "GW_AFFECTED_DATA": 2**4, # Data affected by the GW read window "PERSISTENCE": 2**5, # High persistence (was RESERVED_2) "AD_FLOOR": 2**6, # Below A/D floor (0 DN, was RESERVED_3) - "RESERVED_4": 2**7, # + "OUTLIER": 2**7, # Flagged by outlier detection (was RESERVED_4) "UNRELIABLE_ERROR": 2**8, # Uncertainty exceeds quoted error "NON_SCIENCE": 2**9, # Pixel not on science portion of detector "DEAD": 2**10, # Dead pixel diff --git a/romancal/outlier_detection/outlier_detection.py b/romancal/outlier_detection/outlier_detection.py new file mode 100644 index 000000000..484342fbd --- /dev/null +++ b/romancal/outlier_detection/outlier_detection.py @@ -0,0 +1,467 @@ +"""Primary code for performing outlier detection on JWST observations.""" + +import logging +import warnings +from functools import partial + +import numpy as np +from astropy.stats import sigma_clip +from astropy.units import Quantity +from drizzle.cdrizzle import tblot +from roman_datamodels import datamodels as rdm +from scipy import ndimage + +from romancal.datamodels import ModelContainer +from romancal.lib import dqflags +from romancal.resample import resample +from romancal.resample.resample_utils import build_driz_weight, calc_gwcs_pixmap + +from ..stpipe import RomanStep + +log = logging.getLogger(__name__) +log.setLevel(logging.DEBUG) + +DO_NOT_USE = dqflags.pixel["DO_NOT_USE"] +OUTLIER = dqflags.pixel["OUTLIER"] + + +__all__ = ["OutlierDetection", "flag_cr", "abs_deriv"] + + +class OutlierDetection: + """Main class for performing outlier detection. + + This is the controlling routine for the outlier detection process. + It loads and sets the various input data and parameters needed by + the various functions and then controls the operation of this process + through all the steps used for the detection. + + Notes + ----- + This routine performs the following operations:: + + 1. Extracts parameter settings from input model and merges + them with any user-provided values + 2. Resamples all input images into grouped observation mosaics. + 3. Creates a median image from all grouped observation mosaics. + 4. Blot median image to match each original input image. + 5. Perform statistical comparison between blotted image and original + image to identify outliers. + 6. Updates input data model DQ arrays with mask of detected outliers. + + """ + + default_suffix = "i2d" + + def __init__(self, input_models, **pars): + """ + Initialize the class with input ModelContainers. + + Parameters + ---------- + input_models : list of DataModels, str + list of data models as ModelContainer or ASN file, + one data model for each input image + + pars : dict, optional + Optional user-specified parameters to modify how outlier_detection + will operate. Valid parameters include: + - resample_suffix + + """ + self.input_models = input_models + + self.outlierpars = {} + self.outlierpars.update(pars) + + self.resample_suffix = f"_outlier_{self.default_suffix}.asdf" + log.debug(f"Defined output product suffix as: {self.resample_suffix}") + + # Define how file names are created + self.make_output_path = pars.get( + "make_output_path", partial(RomanStep._make_output_path, None) + ) + + def do_detection(self): + """Flag outlier pixels in DQ of input images.""" # self._convert_inputs() + pars = self.outlierpars + + if pars["resample_data"]: + # Start by creating resampled/mosaic images for + # each group of exposures + resamp = resample.ResampleData( + self.input_models, single=True, blendheaders=False, **pars + ) + drizzled_models = resamp.do_drizzle() + + else: + # for non-dithered data, the resampled image is just the original image + drizzled_models = self.input_models + for i in range(len(self.input_models)): + drizzled_models[i].weight = build_driz_weight( + self.input_models[i], + weight_type="ivm", + good_bits=pars["good_bits"], + ) + + # Initialize intermediate products used in the outlier detection + median_model = rdm.open(drizzled_models[0]).copy() + + # Perform median combination on set of drizzled mosaics + median_model.data = Quantity( + self.create_median(drizzled_models), unit=median_model.data.unit + ) + median_model_output_path = self.make_output_path( + basepath=median_model.meta.filename.replace(self.resample_suffix, ".asdf"), + suffix="median", + ) + median_model.save(median_model_output_path) + log.info(f"Saved model in {median_model_output_path}") + + if pars["resample_data"]: + # Blot the median image back to recreate each input image specified + # in the original input list/ASN/ModelContainer + blot_models = self.blot_median(median_model) + + else: + # Median image will serve as blot image + blot_models = ModelContainer(return_open=False) + for i in range(len(self.input_models)): + blot_models.append(median_model) + + # Perform outlier detection using statistical comparisons between + # each original input image and its blotted version of the median image + self.detect_outliers(blot_models) + + # clean-up (just to be explicit about being finished with + # these results) + del median_model, blot_models + + def create_median(self, resampled_models): + """Create a median image from the singly resampled images. + + NOTES + ----- + This version is simplified from astrodrizzle's version in the + following ways: + - type of combination: fixed to 'median' + - 'minmed' not implemented as an option + """ + maskpt = self.outlierpars.get("maskpt", 0.7) + + log.info("Computing median") + + # Compute weight means without keeping DataModel for eacn input open + # Start by insuring that the ModelContainer does NOT open and keep each datamodel + ropen_orig = resampled_models._return_open + resampled_models._return_open = False # turn off auto-opening of models + # keep track of resulting computation for each input resampled datamodel + weight_thresholds = [] + # For each model, compute the bad-pixel threshold from the weight arrays + for resampled in resampled_models: + m = rdm.open(resampled) + weight = m.weight + # necessary in order to assure that mask gets applied correctly + if hasattr(weight, "_mask"): + del weight._mask + mask_zero_weight = np.equal(weight, 0.0) + mask_nans = np.isnan(weight) + # Combine the masks + weight_masked = np.ma.array( + weight, mask=np.logical_or(mask_zero_weight, mask_nans) + ) + # Sigma-clip the unmasked data + weight_masked = sigma_clip(weight_masked, sigma=3, maxiters=5) + mean_weight = np.mean(weight_masked) + # Mask pixels where weight falls below maskpt percent + weight_threshold = mean_weight * maskpt + weight_thresholds.append(weight_threshold) + # close and delete the model, just to explicitly try to keep the memory as clean as possible + m.close() + del m + # Reset ModelContainer attribute to original value + resampled_models._return_open = ropen_orig + + # Now, set up buffered access to all input models + resampled_models.set_buffer(1.0) # Set buffer at 1Mb + resampled_sections = resampled_models.get_sections() + median_image = np.empty( + (resampled_models.imrows, resampled_models.imcols), + resampled_models.imtype, + ) + median_image[:] = np.nan # initialize with NaNs + + for resampled_sci, resampled_weight, (row1, row2) in resampled_sections: + # Create a mask for each input image, masking out areas where there is + # no data or the data has very low weight + badmasks = [] + for weight, weight_threshold in zip(resampled_weight, weight_thresholds): + badmask = np.less(weight, weight_threshold) + log.debug( + "Percentage of pixels with low weight: {}".format( + np.sum(badmask) / len(weight.flat) * 100 + ) + ) + badmasks.append(badmask) + + # Fill resampled_sci array with nan's where mask values are True + for f1, f2 in zip(resampled_sci, badmasks): + for elem1, elem2 in zip(f1, f2): + elem1[elem2] = np.nan + + del badmasks + + # For a stack of images with "bad" data replaced with Nan + # use np.nanmedian to compute the median. + with warnings.catch_warnings(): + warnings.filterwarnings( + action="ignore", + message="All-NaN slice encountered", + category=RuntimeWarning, + ) + median_image[row1:row2] = np.nanmedian(resampled_sci, axis=0) + del resampled_sci, resampled_weight + + return median_image + + def blot_median(self, median_model): + """Blot resampled median image back to the detector images.""" + interp = self.outlierpars.get("interp", "linear") + sinscl = self.outlierpars.get("sinscl", 1.0) + + # Initialize container for output blot images + blot_models = [] + + log.info("Blotting median") + for model in self.input_models: + blotted_median = model.copy() + + # clean out extra data not related to blot result + blotted_median.err *= 0.0 # None + blotted_median.dq *= 0 # None + + # apply blot to re-create model.data from median image + blotted_median.data = Quantity( + gwcs_blot(median_model, model, interp=interp, sinscl=sinscl), + unit=blotted_median.data.unit, + ) + + model_path = self.make_output_path( + basepath=model.meta.filename, suffix="blot" + ) + blotted_median.save(model_path) + log.info(f"Saved model in {model_path}") + + # Append model name to the ModelContainer so it is not passed in memory + blot_models.append(model_path) + + return ModelContainer(blot_models, return_open=False) + + def detect_outliers(self, blot_models): + """Flag DQ array for cosmic rays in input images. + + The science frame in each ImageModel in input_models is compared to + the corresponding blotted median image in blot_models. The result is + an updated DQ array in each ImageModel in input_models. + + Parameters + ---------- + input_models: JWST ModelContainer object + data model container holding science ImageModels, modified in place + + blot_models : JWST ModelContainer object + data model container holding ImageModels of the median output frame + blotted back to the wcs and frame of the ImageModels in + input_models + + Returns + ------- + None + The dq array in each input model is modified in place + + """ + log.info("Flagging outliers") + for i, (image, blot) in enumerate(zip(self.input_models, blot_models)): + blot = rdm.open(blot) + flag_cr(image, blot, **self.outlierpars) + self.input_models[i] = image + + +def flag_cr( + sci_image, + blot_image, + snr="5.0 4.0", + scale="1.2 0.7", + backg=0, + resample_data=True, + **kwargs, +): + """Masks outliers in science image by updating DQ in-place + + Mask blemishes in dithered data by comparing a science image + with a model image and the derivative of the model image. + + Parameters + ---------- + sci_image : ~romancal.DataModel.ImageModel + the science data + + blot_image : ~romancal.DataModel.ImageModel + the blotted median image of the dithered science frames + + snr : str + Signal-to-noise ratio + + scale : str + scaling factor applied to the derivative + + backg : float + Background value (scalar) to subtract + + resample_data : bool + Boolean to indicate whether blot_image is created from resampled, + dithered data or not + """ + snr1, snr2 = (float(val) for val in snr.split()) + scale1, scale2 = (float(val) for val in scale.split()) + + # Get background level of science data if it has not been subtracted, so it + # can be added into the level of the blotted data, which has been + # background-subtracted + if ( + hasattr(sci_image.meta, "background") + and sci_image.meta.background.subtracted is False + and sci_image.meta.background.level is not None + ): + subtracted_background = sci_image.meta.background.level + log.debug(f"Adding background level {subtracted_background} to blotted image") + else: + # No subtracted background. Allow user-set value, which defaults to 0 + subtracted_background = backg + + sci_data = sci_image.data + blot_data = blot_image.data + blot_deriv = abs_deriv(blot_data.value) + err_data = np.nan_to_num(sci_image.err) + + # create the outlier mask + if resample_data: # dithered outlier detection + blot_data += subtracted_background + diff_noise = np.abs(sci_data - blot_data) + + # Create a boolean mask based on a scaled version of + # the derivative image (dealing with interpolating issues?) + # and the standard n*sigma above the noise + threshold1 = scale1 * blot_deriv + snr1 * err_data.value + mask1 = np.greater(diff_noise.value, threshold1) + + # Smooth the boolean mask with a 3x3 boxcar kernel + kernel = np.ones((3, 3), dtype=int) + mask1_smoothed = ndimage.convolve(mask1, kernel, mode="nearest") + + # Create a 2nd boolean mask based on the 2nd set of + # scale and threshold values + threshold2 = scale2 * blot_deriv + snr2 * err_data.value + mask2 = np.greater(diff_noise.value, threshold2) + + # Final boolean mask + cr_mask = mask1_smoothed & mask2 + + else: # stack outlier detection + diff_noise = np.abs(sci_data - blot_data) + + # straightforward detection of outliers for non-dithered data since + # err_data includes all noise sources (photon, read, and flat for baseline) + cr_mask = np.greater(diff_noise.value, snr1 * err_data.value) + + # Count existing DO_NOT_USE pixels + count_existing = np.count_nonzero(sci_image.dq & DO_NOT_USE) + + # Update the DQ array in the input image. + sci_image.dq = np.bitwise_or(sci_image.dq, cr_mask * (DO_NOT_USE | OUTLIER)) + + # Report number (and percent) of new DO_NOT_USE pixels found + count_outlier = np.count_nonzero(sci_image.dq & DO_NOT_USE) + count_added = count_outlier - count_existing + percent_cr = count_added / (sci_image.shape[0] * sci_image.shape[1]) * 100 + log.info(f"New pixels flagged as outliers: {count_added} ({percent_cr:.2f}%)") + + +def abs_deriv(array): + """Take the absolute derivate of a numpy array.""" + tmp = np.zeros(array.shape, dtype=np.float64) + out = np.zeros(array.shape, dtype=np.float64) + + tmp[1:, :] = array[:-1, :] + tmp, out = _absolute_subtract(array, tmp, out) + tmp[:-1, :] = array[1:, :] + tmp, out = _absolute_subtract(array, tmp, out) + + tmp[:, 1:] = array[:, :-1] + tmp, out = _absolute_subtract(array, tmp, out) + tmp[:, :-1] = array[:, 1:] + tmp, out = _absolute_subtract(array, tmp, out) + + return out + + +def _absolute_subtract(array, tmp, out): + tmp = np.abs(array - tmp) + out = np.maximum(tmp, out) + tmp = tmp * 0.0 + return tmp, out + + +def gwcs_blot(median_model, blot_img, interp="poly5", sinscl=1.0): + """ + Resample the output/resampled image to recreate an input image based on + the input image's world coordinate system + + Parameters + ---------- + median_model : `~roman_datamodels.datamodels.MosaicModel` + + blot_img : datamodel + Datamodel containing header and WCS to define the 'blotted' image + + interp : str, optional + The type of interpolation used in the resampling. The + possible values are "nearest" (nearest neighbor interpolation), + "linear" (bilinear interpolation), "poly3" (cubic polynomial + interpolation), "poly5" (quintic polynomial interpolation), + "sinc" (sinc interpolation), "lan3" (3rd order Lanczos + interpolation), and "lan5" (5th order Lanczos interpolation). + + sincscl : float, optional + The scaling factor for sinc interpolation. + """ + blot_wcs = blot_img.meta.wcs + + # Compute the mapping between the input and output pixel coordinates + pixmap = calc_gwcs_pixmap(blot_wcs, median_model.meta.wcs, blot_img.data.shape) + log.debug(f"Pixmap shape: {pixmap[:, :, 0].shape}") + log.debug(f"Sci shape: {blot_img.data.shape}") + + pix_ratio = 1 + log.info(f"Blotting {blot_img.data.shape} <-- {median_model.data.shape}") + + outsci = np.zeros(blot_img.shape, dtype=np.float32) + + # Currently tblot cannot handle nans in the pixmap, so we need to give some + # other value. -1 is not optimal and may have side effects. But this is + # what we've been doing up until now, so more investigation is needed + # before a change is made. Preferably, fix tblot in drizzle. + pixmap[np.isnan(pixmap)] = -1 + tblot( + median_model.data, + pixmap, + outsci, + scale=pix_ratio, + kscale=1.0, + interp=interp, + exptime=1.0, + misval=0.0, + sinscl=sinscl, + ) + + return outsci diff --git a/romancal/outlier_detection/outlier_detection_step.py b/romancal/outlier_detection/outlier_detection_step.py index f746d1a03..5f9072796 100755 --- a/romancal/outlier_detection/outlier_detection_step.py +++ b/romancal/outlier_detection/outlier_detection_step.py @@ -1,8 +1,11 @@ -""" -Detect and flag outlier in a science image -""" +"""Public common step definition for OutlierDetection processing.""" -import roman_datamodels as rdm +import contextlib +import os +from functools import partial + +from romancal.datamodels import ModelContainer +from romancal.outlier_detection import outlier_detection from ..stpipe import RomanStep @@ -10,28 +13,161 @@ class OutlierDetectionStep(RomanStep): - """Detect and flag outliers in a science image.""" + """Flag outlier bad pixels and cosmic rays in DQ array of each input image. + + Input images can be listed in an input association file or already opened + with a ModelContainer. DQ arrays are modified in place. + + Parameters + ----------- + input_data : asn file or ~romancal.datamodels.ModelContainer + Single filename association table, or a datamodels.ModelContainer. + + """ + + class_alias = "outlier_detection" + + # The members of spec needs to be a super-set of all parameters needed + # by the various versions of the outlier_detection algorithms, and each + # version will pick and choose what they need while ignoring the rest. + spec = """ + weight_type = option('ivm','exptime',default='ivm') + pixfrac = float(default=1.0) + kernel = string(default='square') # drizzle kernel + fillval = string(default='INDEF') + nlow = integer(default=0) + nhigh = integer(default=0) + maskpt = float(default=0.7) + grow = integer(default=1) + snr = string(default='5.0 4.0') + scale = string(default='1.2 0.7') + backg = float(default=0.0) + kernel_size = string(default='7 7') + threshold_percent = float(default=99.8) + ifu_second_check = boolean(default=False) + save_intermediate_results = boolean(default=False) + resample_data = boolean(default=True) + good_bits = string(default="~DO_NOT_USE") # DQ flags to allow + scale_detection = boolean(default=False) + search_output_file = boolean(default=False) + allowed_memory = float(default=None) # Fraction of memory to use for the combined image + in_memory = boolean(default=False) + """ + + def process(self, input_models): + """Perform outlier detection processing on input data.""" + + self.input_models = input_models + self.input_container = isinstance(self.input_models, ModelContainer) + # Setup output path naming if associations are involved. + asn_id = None + with contextlib.suppress(AttributeError, KeyError): + asn_id = self.input_models.meta.asn_table.asn_id + if asn_id is None: + asn_id = self.search_attr("asn_id") + if asn_id is not None: + # TODO: revisit this once we have a better idea of the associations + _make_output_path = self.search_attr("_make_output_path", parent_first=True) + + self._make_output_path = partial(_make_output_path, asn_id=asn_id) + + # Setup outlier detection parameters + pars = { + "weight_type": self.weight_type, + "pixfrac": self.pixfrac, + "kernel": self.kernel, + "fillval": self.fillval, + "nlow": self.nlow, + "nhigh": self.nhigh, + "maskpt": self.maskpt, + "grow": self.grow, + "snr": self.snr, + "scale": self.scale, + "backg": self.backg, + "kernel_size": self.kernel_size, + "threshold_percent": self.threshold_percent, + "ifu_second_check": self.ifu_second_check, + "allowed_memory": self.allowed_memory, + "in_memory": self.in_memory, + "save_intermediate_results": self.save_intermediate_results, + "resample_data": self.resample_data, + "good_bits": self.good_bits, + "make_output_path": self.make_output_path, + } + + # Add logic here to select which version of OutlierDetection + # needs to be used depending on the input data + if self.input_container: + single_model = self.input_models[0] + else: + single_model = self.input_models + exptype = single_model.meta.exposure.type + self.check_input() + + if exptype == "WFI_IMAGE": + # imaging with resampling + detection_step = outlier_detection.OutlierDetection + pars["resample_suffix"] = "i2d" + else: + self.log.error( + "Outlier detection failed for unknown/unsupported ", + f"exposure type: {exptype}", + ) + self.valid_input = False - def process(self, input): - input_model = rdm.open(input, lazy_load=False) + if not self.valid_input: + if self.input_container: + for model in self.input_models: + model.meta.cal_step.outlier_detection = "SKIPPED" + else: + self.input_models.meta.cal_step.outlier_detection = "SKIPPED" + self.skip = True + return self.input_models - # No reference files - # reference_file_model = {} - # Do the outlier detection - # output_model = outlier_detection.do_correction( - # input_model, - # reference_file_model, - # ) + self.log.debug(f"Using {detection_step.__name__} class for outlier_detection") - output_model = input_model + # Set up outlier detection, then do detection + step = detection_step(self.input_models, **pars) + step.do_detection() - # Close the input and reference files - input_model.close() + state = "COMPLETE" + if self.input_container: + if not self.save_intermediate_results: + self.log.debug( + "The following files will be deleted since save_intermediate_results=False:" + ) + for model in self.input_models: + model.meta.cal_step.outlier_detection = state + if not self.save_intermediate_results: + # Remove unwanted files + crf_path = self.make_output_path(basepath=model.meta.filename) + # These lines to be used when/if outlier_i2d files follow output_dir + # crf_file = os.path.basename(crf_path) + # outlr_path = crf_path.replace(crf_file, outlr_file) + outlr_file = model.meta.filename.replace("cal", "outlier_i2d") + blot_path = crf_path.replace("crf", "blot") + median_path = blot_path.replace("blot", "median") - if self.save_results: - try: - self.suffix = "outlier_detection" - except AttributeError: - self["suffix"] = "outlier_detection" + for fle in [outlr_file, blot_path, median_path]: + if os.path.isfile(fle): + os.remove(fle) + self.log.debug(f" {fle}") + else: + self.input_models.meta.cal_step.outlier_detection = state + return self.input_models - return output_model + def check_input(self): + """Use this method to determine whether input is valid or not.""" + if self.input_container: + ninputs = len(self.input_models) + if not isinstance(self.input_models, ModelContainer): + self.log.warning("Input is not a ModelContainer") + self.log.warning("Outlier detection step will be skipped") + self.valid_input = False + elif ninputs < 2: + self.log.warning(f"Input only contains {ninputs} exposure") + self.log.warning("Outlier detection step will be skipped") + self.valid_input = False + else: + self.valid_input = True + self.log.info(f"Performing outlier detection on {ninputs} inputs") diff --git a/romancal/resample/resample.py b/romancal/resample/resample.py index 25206370e..aa71aecf4 100644 --- a/romancal/resample/resample.py +++ b/romancal/resample/resample.py @@ -186,6 +186,7 @@ def resample_many_to_many(self): """ for exposure in self.input_models.models_grouped: output_model = self.blank_output + output_model.meta["resample"] = mk_resample() # Determine output file type from input exposure filenames # Use this for defining the output filename indx = exposure[0].meta.filename.rfind(".") @@ -204,6 +205,7 @@ def resample_many_to_many(self): ) log.info(f"{len(exposure)} exposures to drizzle together") + output_list = [] for img in exposure: img = datamodels.open(img) # TODO: should weight_type=None here? @@ -212,6 +214,8 @@ def resample_many_to_many(self): ) # apply sky subtraction + if not hasattr(img.meta, "background"): + self._create_background_attribute(img) blevel = img.meta.background.level if not img.meta.background.subtracted and blevel is not None: data = img.data - blevel @@ -237,13 +241,17 @@ def resample_many_to_many(self): if not self.in_memory: # Write out model to disk, then return filename output_name = output_model.meta.filename + # cast context array to uint32 + output_model.context = output_model.context.astype("uint32") output_model.save(output_name) log.info(f"Exposure {output_name} saved to file") - self.output_models.append(output_name) + output_list.append(output_name) else: - self.output_models.append(output_model.copy()) + output_list.append(output_model.copy()) + + self.output_models = ModelContainer(output_list, return_open=self.in_memory) output_model.data *= 0.0 - output_model.wht *= 0.0 + output_model.weight *= 0.0 return self.output_models @@ -274,11 +282,8 @@ def resample_many_to_one(self): inwht = resample_utils.build_driz_weight( img, weight_type=self.weight_type, good_bits=self.good_bits ) - # apply sky subtraction - # NOTE: mocking a sky-subtracted image (remove this later on) - img.meta["background"] = {} - img.meta.background["level"] = 0 - img.meta.background["subtracted"] = True + if not hasattr(img.meta, "background"): + self._create_background_attribute(img) blevel = img.meta.background.level if not img.meta.background.subtracted and blevel is not None: data = img.data - blevel @@ -333,6 +338,11 @@ def resample_many_to_one(self): return self.output_models + def _create_background_attribute(self, img): + img.meta["background"] = {} + img.meta.background["level"] = 0 + img.meta.background["subtracted"] = True + def resample_variance_array(self, name, output_model): """Resample variance arrays from ``self.input_models`` to the ``output_model``. @@ -413,7 +423,8 @@ def resample_variance_array(self, name, output_model): setattr(output_model, name, output_variance) def resample_exposure_time(self, output_model): - """Resample the exposure time from ``self.input_models`` to the ``output_model``. + """Resample the exposure time from ``self.input_models`` to the + ``output_model``. Create an exposure time image that is the drizzled sum of the input images. @@ -466,10 +477,7 @@ def resample_exposure_time(self, output_model): def update_exposure_times(self, output_model, exptime_tot): """Update exposure time metadata (in-place).""" m = exptime_tot > 0 - if np.any(m): - total_exposure_time = np.median(exptime_tot[m]) - else: - total_exposure_time = 0 + total_exposure_time = np.median(exptime_tot[m]) if np.any(m) else 0 max_exposure_time = np.max(exptime_tot) log.info( f"Mean, max exposure times: {total_exposure_time:.1f}, " diff --git a/romancal/stpipe/integration.py b/romancal/stpipe/integration.py index e74507734..72ee451e2 100644 --- a/romancal/stpipe/integration.py +++ b/romancal/stpipe/integration.py @@ -32,7 +32,7 @@ def get_steps(): ("romancal.step.RefPixStep", None, False), ("romancal.step.SaturationStep", None, False), ("romancal.step.AssignWcsStep", None, False), - ("romancal.step.OutlierDetectionStep", None, False), + ("romancal.step.OutlierDetectionStep", "outlier_detection", False), ("romancal.step.SourceDetectionStep", None, False), ("romancal.step.SkyMatchStep", "skymatch", False), ("romancal.step.TweakRegStep", "tweakreg", False),