diff --git a/changes/1357.outlier_detection.0.rst b/changes/1357.outlier_detection.0.rst new file mode 100644 index 000000000..3c3667401 --- /dev/null +++ b/changes/1357.outlier_detection.0.rst @@ -0,0 +1 @@ +Remove unused arguments to outlier detection. diff --git a/changes/1357.outlier_detection.1.rst b/changes/1357.outlier_detection.1.rst new file mode 100644 index 000000000..798d48c35 --- /dev/null +++ b/changes/1357.outlier_detection.1.rst @@ -0,0 +1 @@ +Update input handling to raise an exception on an invalid input instead of issuing a warning and skipping the step. diff --git a/changes/1357.outlier_detection.2.rst b/changes/1357.outlier_detection.2.rst new file mode 100644 index 000000000..dc0be573d --- /dev/null +++ b/changes/1357.outlier_detection.2.rst @@ -0,0 +1 @@ +Use stcal common code in outlier detection. diff --git a/docs/roman/outlier_detection/arguments.rst b/docs/roman/outlier_detection/arguments.rst index 1ecbe63ba..0f8e6446f 100644 --- a/docs/roman/outlier_detection/arguments.rst +++ b/docs/roman/outlier_detection/arguments.rst @@ -1,16 +1,19 @@ .. _outlier_detection_step_args: +For more details about step arguments (including datatypes, possible values +and defaults) see :py:obj:`romancal.outlier_detection.OutlierDetectionStep.spec`. + Step Arguments ============== The ``outlier_detection`` step has the following optional arguments that control the behavior of the processing: -``--weight_type`` (string, default='exptime') +``--weight_type`` The type of data weighting to use during resampling the images for creating the median image used for detecting outliers; options are `'ivm'`, `'exptime'`, and `None` (see :ref:`weight_type_options_details_section` for details). -``--pixfrac`` (float, default=1.0) +``--pixfrac`` Fraction by which input pixels are “shrunk” before being drizzled onto the output image grid, given as a real number between 0 and 1. This specifies the size of the footprint, or “dropsize”, of a pixel in units of the input pixel size. If `pixfrac` @@ -20,7 +23,7 @@ behavior of the processing: output drizzled image is fully populated with pixels from the input image. Valid values range from 0.0 to 1.0. -``--kernel`` (string, default='square') +``--kernel`` This parameter specifies the form of the kernel function used to distribute flux onto the separate output images, for the initial separate drizzling operation only. The value options for this parameter include: @@ -43,7 +46,7 @@ behavior of the processing: should never be used for ``pixfrac != 1.0``, and is not recommended for ``scale!=1.0``. -``--fillval`` (string, default='INDEF') +``--fillval`` The value for this parameter is to be assigned to the output pixels that have zero weight or which do not receive flux from any input pixels during drizzling. This parameter corresponds to the ``fillval`` parameter of the @@ -55,77 +58,53 @@ behavior of the processing: Any floating-point value, given as a string, is valid. A value of 'INDEF' will use the last zero weight flux. -``--nlow`` (integer, default=0) - The number of low values in each pixel stack to ignore when computing the median - value. - -``--nhigh`` (integer, default=0) - The number of high values in each pixel stack to ignore when computing the median - value. - -``--maskpt`` (float, default=0.7) +``--maskpt`` Percentage of weight image values below which they are flagged as bad and rejected from the median image. Valid values range from 0.0 to 1.0. -``--grow`` (integer, default=1) - The distance, in pixels, beyond the limit set by the rejection algorithm being - used, for additional pixels to be rejected in an image. - -``--snr`` (string, default='4.0 3.0') +``--snr`` The signal-to-noise values to use for bad pixel identification. Since cosmic rays often extend across several pixels the user must specify two cut-off values for determining whether a pixel should be masked: the first for detecting the primary cosmic ray, and the second (typically lower threshold) for masking lower-level bad pixels adjacent to those found in the first pass. Valid values are a pair of - floating-point values in a single string. + floating-point values in a single string (for example "5.0 4.0"). -``--scale`` (string, default='0.5 0.4') +``--scale`` The scaling factor applied to derivative used to identify bad pixels. Since cosmic rays often extend across several pixels the user must specify two cut-off values for determining whether a pixel should be masked: the first for detecting the primary cosmic ray, and the second (typically lower threshold) for masking lower-level bad pixels adjacent to those found in the first pass. Valid values are a pair of - floating-point values in a single string. + floating-point values in a single string (for example "1.2 0.7"). -``--backg`` (float, default=0.0) +``--backg`` User-specified background value (scalar) to subtract during final identification step of outliers in `driz_cr` computation. -``--kernel_size`` (string, default='7 7') - Size of kernel to be used during resampling of the data - (i.e. when `resample_data=True`). - -``--save_intermediate_results`` (boolean, default=False) - Specifies whether or not to write out intermediate products such as median image or +``--save_intermediate_results`` + Boolean specifying whether or not to write out intermediate products such as median image or resampled individual input exposures to disk. Typically, only used to track down problems with final results when too many or too few pixels are flagged as outliers. -``--resample_data`` (boolean, default=True) - Specifies whether or not to resample the input images when performing outlier +``--resample_data`` + Boolean specifying whether or not to resample the input images when performing outlier detection. -``--good_bits`` (string, default=0) +``--good_bits`` The DQ bit values from the input image DQ arrays that should be considered 'good' when creating masks of bad pixels during outlier detection when resampling the data. See `Roman's Data Quality Flags `_ for details. -``--allowed_memory`` (float, default=None) - Specifies the fractional amount of free memory to allow when creating the resampled - image. If ``None``, the environment variable ``DMODEL_ALLOWED_MEMORY`` is used. If - not defined, no check is made. If the resampled image would be larger than specified, - an ``OutputTooLargeError`` exception will be generated. For example, if set to - ``0.5``, only resampled images that use less than half the available memory can be - created. - -``--in_memory`` (boolean, default=False) - Specifies whether or not to keep all intermediate products and datamodels in +``--in_memory`` + Boolean specifying whether or not to keep all intermediate products and datamodels in memory at the same time during the processing of this step. If set to `False`, - all input and output data will be written to disk at the start of the step - (as much as `roman_datamodels` will allow, anyway), then read in to memory only when - accessed. This results in a much lower memory profile at the expense of file I/O, - which can allow large mosaics to process in more limited amounts of memory. + any `ModelLibrary` opened by this step will use ``on_disk=True`` and use temporary + files to store model modifications. Additionally any resampled images will + be kept in memory (as long as needed). This can result in much lower memory + usage (at the expense of file I/O) to process large associations. .. _weight_type_options_details_section: diff --git a/docs/roman/outlier_detection/outlier_detection.rst b/docs/roman/outlier_detection/outlier_detection.rst index c11aa482c..7fee83327 100644 --- a/docs/roman/outlier_detection/outlier_detection.rst +++ b/docs/roman/outlier_detection/outlier_detection.rst @@ -55,20 +55,13 @@ Specifically, this routine performs the following operations: * The median image is created by combining all grouped mosaic images or non-resampled input data pixel-by-pixel. - * The ``nlow`` and ``nhigh`` parameters specify how many low and high values - to ignore when computing the median for any given pixel. * The ``maskpt`` parameter sets the percentage of the weight image values to use, and any pixel with a weight below this value gets flagged as "bad" and ignored when resampled. - * The ``grow`` parameter sets the width, in pixels, beyond the limit set by - the rejection algorithm being used, for additional pixels to be rejected in - an image. - * The median image is written out to disk as `__median` by default. #. By default, the median image is blotted back (inverse of resampling) to match each original input image. - * Blotted images are written out to disk as `__blot` by default. * **If resampling is turned off**, the median image is compared directly to each input image. @@ -136,26 +129,16 @@ memory usage at the expense of file I/O. The control over this memory model hap with the use of the ``in_memory`` parameter. The full impact of this parameter during processing includes: -#. The ``save_open`` parameter gets set to `False` +#. The ``on_disk`` parameter gets set to `True` when opening the input :py:class:`~romancal.datamodels.library.ModelLibrary` - object. This forces all input models in the input - :py:class:`~romancal.datamodels.library.ModelLibrary` to get written out to disk. - It then uses the filename of the input model during subsequent processing. + object. This causes modified models to be written to temporary files. -#. The ``in_memory`` parameter gets passed to the :py:class:`~romancal.resample.ResampleStep` - to set whether or not to keep the resampled images in memory or not. By default, - the outlier detection processing sets this parameter to `False` so that each resampled - image gets written out to disk. - -#. Computing the median image works section-by-section by only keeping 1Mb of each input - in memory at a time. As a result, only the final output product array for the final - median image along with a stack of 1Mb image sections are kept in memory. - -#. The final resampling step also avoids keeping all inputs in memory by only reading - each input into memory 1 at a time as it gets resampled onto the final output product. +#. Computing the median image uses temporary files. Each resampled group + is split into sections (1 per "row") and each section is appended to a different + temporary file. After resampling all groups, each temporary file is read and a + median is computed for all sections in that file (yielding a median for that + section across all resampled groups). Finally, these median sections are + combined into a final median image. These changes result in a minimum amount of memory usage during processing at the obvious expense of reading and writing the products from disk. - - -.. automodapi:: romancal.outlier_detection.outlier_detection diff --git a/docs/roman/outlier_detection/outlier_detection_step.rst b/docs/roman/outlier_detection/outlier_detection_step.rst index 0be906290..353bca28d 100644 --- a/docs/roman/outlier_detection/outlier_detection_step.rst +++ b/docs/roman/outlier_detection/outlier_detection_step.rst @@ -4,9 +4,8 @@ OutlierDetectionStep -------------------- This module provides the sole interface to all methods of performing outlier detection -on Roman observations. The outlier detection algorithm used for WFI data is implemented -in :py:class:`~romancal.outlier_detection.outlier_detection.OutlierDetection` -and described in :ref:`outlier-detection-imaging`. +on Roman observations. The outlier detection algorithm used for WFI data is +described in :ref:`outlier-detection-imaging`. .. note:: Whether the data are being provided in an `association file`_ or as a list of ASDF filenames, diff --git a/romancal/outlier_detection/_fileio.py b/romancal/outlier_detection/_fileio.py new file mode 100644 index 000000000..d1dde8489 --- /dev/null +++ b/romancal/outlier_detection/_fileio.py @@ -0,0 +1,45 @@ +import logging + +from astropy.units import Quantity + +log = logging.getLogger(__name__) +log.setLevel(logging.DEBUG) + + +def save_median(example_model, median_data, median_wcs, make_output_path): + _save_intermediate_output( + _make_median_model(example_model, median_data, median_wcs), + "median", + make_output_path, + ) + + +def save_drizzled(drizzled_model, make_output_path): + _save_intermediate_output(drizzled_model, "outlier_i2d", make_output_path) + + +def _make_median_model(example_model, data, wcs): + model = example_model.copy() + model.data = Quantity(data, unit=model.data.unit) + model.meta.filename = "drizzled_median.asdf" + model.meta.wcs = wcs + return model + + +def _save_intermediate_output(model, suffix, make_output_path): + """ + Ensure all intermediate outputs from OutlierDetectionStep have consistent file naming conventions + + Notes + ----- + self.make_output_path() is updated globally for the step in the main pipeline + to include the asn_id in the output path, so no need to handle it here. + """ + + # outlier_?2d is not a known suffix, and make_output_path cannot handle an + # underscore in an unknown suffix, so do a manual string replacement + input_path = model.meta.filename.replace("_outlier_", "_") + + output_path = make_output_path(input_path, suffix=suffix) + model.save(output_path) + log.info(f"Saved {suffix} model in {output_path}") diff --git a/romancal/outlier_detection/outlier_detection.py b/romancal/outlier_detection/outlier_detection.py deleted file mode 100644 index 47daeddd5..000000000 --- a/romancal/outlier_detection/outlier_detection.py +++ /dev/null @@ -1,404 +0,0 @@ -"""Primary code for performing outlier detection on Roman observations.""" - -import copy -import logging -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.dqflags import pixel -from scipy import ndimage - -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) - -__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 ModelLibrary. - - Parameters - ---------- - input_models : ~romancal.datamodels.ModelLibrary - A `~romancal.datamodels.ModelLibrary` object containing the data - to be processed. - - 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 = dict(pars) - self.resample_suffix = f"_outlier_{self.default_suffix if pars.get('resample_suffix') is None else pars.get('resample_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 - with drizzled_models: - for i, model in enumerate(drizzled_models): - model["weight"] = build_driz_weight( - model, - weight_type="ivm", - good_bits=pars["good_bits"], - ) - drizzled_models.shelve(model, i) - - # Perform median combination on set of drizzled mosaics - median_data = self.create_median(drizzled_models) - - # Initialize intermediate products used in the outlier detection - with drizzled_models: - example_model = drizzled_models.borrow(0) - median_wcs = copy.deepcopy(example_model.meta.wcs) - if pars["save_intermediate_results"]: - median_model = example_model.copy() - median_model.data = Quantity(median_data, unit=median_model.data.unit) - median_model.meta.filename = "drizzled_median.asdf" - median_model_output_path = self.make_output_path( - basepath=median_model.meta.filename, - suffix="median", - ) - median_model.save(median_model_output_path) - log.info(f"Saved model in {median_model_output_path}") - drizzled_models.shelve(example_model, 0, modify=False) - - # Perform outlier detection using statistical comparisons between - # each original input image and its blotted version of the median image - self.detect_outliers(median_data, median_wcs, pars["resample_data"]) - - # clean-up (just to be explicit about being finished with - # these results) - del median_data, median_wcs - - 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") - - data = [] - - # Compute weight means without keeping DataModel for eacn input open - # keep track of resulting computation for each input resampled datamodel - weight_thresholds = [] - # For each model, compute the bad-pixel threshold from the weight arrays - with resampled_models: - for i, model in enumerate(resampled_models): - weight = model.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) - this_data = model.data.copy() - this_data[model.weight < weight_threshold] = np.nan - data.append(this_data) - - resampled_models.shelve(model, i, modify=False) - - median_image = np.nanmedian(data, axis=0) - return median_image - - def detect_outliers(self, median_data, median_wcs, resampled): - """Flag DQ array for cosmic rays in input images. - - The science frame in each ImageModel in self.input_models is compared to - the a blotted median image (generated with median_data and median_wcs). - The result is an updated DQ array in each ImageModel in input_models. - - Parameters - ---------- - median_data : numpy.ndarray - Median array that will be used as the "reference" for detecting - outliers. - - median_wcs : gwcs.WCS - WCS for the median data - - resampled : bool - True if the median data was generated from resampling the input - images. - - Returns - ------- - None - The dq array in each input model is modified in place - - """ - interp = self.outlierpars.get("interp", "linear") - sinscl = self.outlierpars.get("sinscl", 1.0) - log.info("Flagging outliers") - with self.input_models: - for i, image in enumerate(self.input_models): - # make blot_data Quantity (same unit as image.data) - if resampled: - # blot back onto image - blot_data = Quantity( - gwcs_blot( - median_data, median_wcs, image, interp=interp, sinscl=sinscl - ), - unit=image.data.unit, - ) - else: - # use median - blot_data = Quantity(median_data, unit=image.data.unit, copy=True) - flag_cr(image, blot_data, **self.outlierpars) - self.input_models.shelve(image, i) - - -def flag_cr( - sci_image, - blot_data, - 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_data : Quantity - 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_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 & pixel.DO_NOT_USE) - - # Update the DQ array values in the input image but preserve datatype. - sci_image.dq = np.bitwise_or( - sci_image.dq, cr_mask * (pixel.DO_NOT_USE | pixel.OUTLIER) - ).astype(np.uint32) - - # Report number (and percent) of new DO_NOT_USE pixels found - count_outlier = np.count_nonzero(sci_image.dq & pixel.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 derivative 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_data, median_wcs, blot_img, interp="poly5", sinscl=1.0): - """ - Resample the median_data to recreate an input image based on - the blot_img's WCS. - - Parameters - ---------- - median_data : numpy.ndarray - Median data used as the source data for blotting. - - median_wcs : gwcs.WCS - WCS for median_data. - - 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). - - sinscl : 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_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_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_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 332791e4b..1029aad87 100644 --- a/romancal/outlier_detection/outlier_detection_step.py +++ b/romancal/outlier_detection/outlier_detection_step.py @@ -1,10 +1,9 @@ """Public common step definition for OutlierDetection processing.""" from functools import partial -from pathlib import Path from romancal.datamodels import ModelLibrary -from romancal.outlier_detection import outlier_detection +from romancal.outlier_detection.utils import detect_outliers from ..stpipe import RomanStep @@ -34,37 +33,23 @@ class OutlierDetectionStep(RomanStep): pixfrac = float(default=1.0) # Fraction by which input pixels are shrunk before being drizzled onto the output image grid kernel = string(default='square') # Shape of the kernel used for flux distribution onto output images fillval = string(default='INDEF') # Value assigned to output pixels that have zero weight or no flux during drizzling - nlow = integer(default=0) # The number of low values in each pixel stack to ignore when computing the median value - nhigh = integer(default=0) # The number of high values in each pixel stack to ignore when computing the median value maskpt = float(default=0.7) # Percentage of weight image values below which they are flagged as bad pixels - grow = integer(default=1) # The distance beyond the rejection limit for additional pixels to be rejected in an image snr = string(default='5.0 4.0') # The signal-to-noise values to use for bad pixel identification scale = string(default='1.2 0.7') # The scaling factor applied to derivative used to identify bad pixels backg = float(default=0.0) # User-specified background value to subtract during final identification step - kernel_size = string(default='7 7') # Size of kernel to be used during resampling of the data save_intermediate_results = boolean(default=False) # Specifies whether or not to write out intermediate products to disk resample_data = boolean(default=True) # Specifies whether or not to resample the input images when performing outlier detection good_bits = string(default="~DO_NOT_USE+NON_SCIENCE") # DQ bit value to be considered 'good' - allowed_memory = float(default=None) # Fraction of memory to use for the combined image in_memory = boolean(default=False) # Specifies whether or not to keep all intermediate products and datamodels in memory """ # noqa: E501 def process(self, input_models): """Perform outlier detection processing on input data.""" - self.skip = False - if isinstance(input_models, ModelLibrary): library = input_models else: - try: - library = ModelLibrary(input_models) - except Exception: - self.log.warning( - "Skipping outlier_detection - input cannot be parsed into a ModelLibrary." - ) - self.skip = True - return input_models + library = ModelLibrary(input_models) # check number of input models if len(library) < 2: @@ -74,90 +59,50 @@ def process(self, input_models): self.log.warning( "Skipping outlier_detection - at least two imaging observations are needed." ) - self.skip = True - # check that all inputs are WFI_IMAGE - if not self.skip: - with library: - for i, model in enumerate(library): - if model.meta.exposure.type != "WFI_IMAGE": - self.skip = True - library.shelve(model, i, modify=False) - if self.skip: - self.log.warning( - "Skipping outlier_detection - all WFI_IMAGE exposures are required." - ) - - # if skipping for any reason above... - if self.skip: - # set meta.cal_step.outlier_detection to SKIPPED - with library: - for i, model in enumerate(library): - model.meta.cal_step["outlier_detection"] = "SKIPPED" - library.shelve(model, i) + def set_skip(model, index): + model.meta.cal_step["outlier_detection"] = "SKIPPED" + + list(library.map_function(set_skip)) + return library + # check that all inputs are WFI_IMAGE + def get_exptype(model, index): + return model.meta.exposure.type + + exptypes = list(library.map_function(get_exptype, modify=False)) + if any(exptype != "WFI_IMAGE" for exptype in exptypes): + raise ValueError( + f"outlier_detection only supports WFI_IMAGE exposure types: {set(exptypes)}" + ) + # Setup output path naming if associations are involved. asn_id = library.asn.get("asn_id", None) if asn_id is not None: _make_output_path = self.search_attr("_make_output_path", parent_first=True) self._make_output_path = partial(_make_output_path, asn_id=asn_id) - detection_step = outlier_detection.OutlierDetection - 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, - "save_intermediate_results": self.save_intermediate_results, - "resample_data": self.resample_data, - "good_bits": self.good_bits, - "allowed_memory": self.allowed_memory, - "in_memory": self.in_memory, - "make_output_path": self.make_output_path, - "resample_suffix": "i2d", - } - - self.log.debug(f"Using {detection_step.__name__} class for outlier_detection") + snr1, snr2 = (float(v) for v in self.snr.split()) + scale1, scale2 = (float(v) for v in self.scale.split()) # Set up outlier detection, then do detection - step = detection_step(library, **pars) - step.do_detection() - - state = "COMPLETE" - - if not self.save_intermediate_results: - self.log.debug( - "The following files will be deleted since \ - save_intermediate_results=False:" - ) - with library: - for i, model in enumerate(library): - model.meta.cal_step["outlier_detection"] = state - if not self.save_intermediate_results: - # remove intermediate files found in - # make_output_path() and the local dir - intermediate_files_paths = [ - Path(self.make_output_path()).parent, - Path().cwd(), - ] - intermediate_files_suffixes = ( - "*blot.asdf", - "*median.asdf", - f'*outlier_{pars.get("resample_suffix")}.asdf', - ) - for current_path in intermediate_files_paths: - for suffix in intermediate_files_suffixes: - for filename in current_path.glob(suffix): - filename.unlink() - self.log.debug(f" {filename}") - library.shelve(model, i) + detect_outliers( + library, + self.weight_type, + self.pixfrac, + self.kernel, + self.fillval, + self.maskpt, + snr1, + snr2, + scale1, + scale2, + self.backg, + self.save_intermediate_results, + self.resample_data, + self.good_bits, + self.in_memory, + self.make_output_path, + ) return library diff --git a/romancal/outlier_detection/tests/test_outlier_detection.py b/romancal/outlier_detection/tests/test_outlier_detection.py index d27ed9ed9..c859a1d52 100644 --- a/romancal/outlier_detection/tests/test_outlier_detection.py +++ b/romancal/outlier_detection/tests/test_outlier_detection.py @@ -6,22 +6,20 @@ from astropy.units import Quantity from romancal.datamodels import ModelLibrary -from romancal.outlier_detection import OutlierDetectionStep, outlier_detection +from romancal.outlier_detection import OutlierDetectionStep @pytest.mark.parametrize( "input_models", [ - list(), "", ], ) -def test_outlier_raises_error_on_invalid_input_models(input_models, caplog): - """Test that OutlierDetection logs out a WARNING if input is invalid.""" +def test_outlier_raises_error_on_invalid_input_models(input_models): + """Test that OutlierDetection raises an Exception if input is invalid.""" - OutlierDetectionStep.call(input_models) - - assert "WARNING" in [x.levelname for x in caplog.records] + with pytest.raises(IsADirectoryError): + OutlierDetectionStep.call(input_models) def test_outlier_skips_step_on_invalid_number_of_elements_in_input(base_image): @@ -37,21 +35,17 @@ def test_outlier_skips_step_on_invalid_number_of_elements_in_input(base_image): res.shelve(m, i, modify=False) -def test_outlier_skips_step_on_exposure_type_different_from_wfi_image(base_image): +def test_outlier_raises_exception_on_exposure_type_different_from_wfi_image(base_image): """ - Test if the outlier detection step is skipped when the exposure type is different from WFI image. + Test if the outlier detection step raises an Exception for non-image inputs. """ img_1 = base_image() img_1.meta.exposure.type = "WFI_PRISM" img_2 = base_image() img_2.meta.exposure.type = "WFI_PRISM" - res = OutlierDetectionStep.call(ModelLibrary([img_1, img_2])) - - with res: - for i, m in enumerate(res): - assert m.meta.cal_step.outlier_detection == "SKIPPED" - res.shelve(m, i, modify=False) + with pytest.raises(ValueError, match="only supports WFI_IMAGE"): + OutlierDetectionStep.call(ModelLibrary([img_1, img_2])) def test_outlier_valid_input_asn(tmp_path, base_image, create_mock_asn_file): @@ -104,91 +98,26 @@ def test_outlier_valid_input_modelcontainer(tmp_path, base_image): res.shelve(m, i, modify=False) -@pytest.mark.parametrize( - "pars", - [ - { - "weight_type": "exptime", - "pixfrac": 1.0, - "kernel": "square", - "fillval": "INDEF", - "nlow": 0, - "nhigh": 0, - "maskpt": 0.7, - "grow": 1, - "snr": "4.0 3.0", - "scale": "0.5 0.4", - "backg": 0.0, - "kernel_size": "7 7", - "save_intermediate_results": False, - "resample_data": True, - "good_bits": 0, - "allowed_memory": None, - "in_memory": True, - "make_output_path": None, - "resample_suffix": "i2d", - }, - { - "weight_type": "exptime", - "save_intermediate_results": True, - "make_output_path": None, - "resample_suffix": "some_other_suffix", - }, - ], -) -def test_outlier_init_default_parameters(pars, base_image): - """ - Test parameter setting on initialization for OutlierDetection. - """ - img_1 = base_image() - img_1.meta.filename = "img_1.asdf" - input_models = ModelLibrary([img_1]) - - step = outlier_detection.OutlierDetection(input_models, **pars) - - assert step.input_models == input_models - assert step.outlierpars == pars - assert step.make_output_path == pars["make_output_path"] - assert step.resample_suffix == f"_outlier_{pars['resample_suffix']}.asdf" - - def test_outlier_do_detection_write_files_to_custom_location(tmp_path, base_image): """ Test that OutlierDetection can create files on disk in a custom location. """ img_1 = base_image() - img_1.meta.filename = "img_1.asdf" + img_1.meta.filename = "img1_cal.asdf" + img_1.meta.background.level = 0 * u.DN / u.s img_2 = base_image() - img_2.meta.filename = "img_2.asdf" + img_2.meta.filename = "img2_cal.asdf" + img_2.meta.background.level = 0 * u.DN / u.s input_models = ModelLibrary([img_1, img_2]) - outlier_step = OutlierDetectionStep() + outlier_step = OutlierDetectionStep( + in_memory=False, + save_intermediate_results=True, + ) # set output dir for all files created by the step outlier_step.output_dir = tmp_path.as_posix() - # make sure files are written out to disk - outlier_step.in_memory = False - pars = { - "weight_type": "exptime", - "pixfrac": 1.0, - "kernel": "square", - "fillval": "INDEF", - "nlow": 0, - "nhigh": 0, - "maskpt": 0.7, - "grow": 1, - "snr": "4.0 3.0", - "scale": "0.5 0.4", - "backg": 0.0, - "kernel_size": "7 7", - "save_intermediate_results": True, - "resample_data": False, - "good_bits": 0, - "allowed_memory": None, - "in_memory": outlier_step.in_memory, - "make_output_path": outlier_step.make_output_path, - "resample_suffix": "i2d", - } + outlier_step(input_models) # meta.filename for the median image created by OutlierDetection.do_detection() median_path = tmp_path / "drizzled_median.asdf" @@ -197,9 +126,6 @@ def test_outlier_do_detection_write_files_to_custom_location(tmp_path, base_imag median_path, ] - step = outlier_detection.OutlierDetection(input_models, **pars) - step.do_detection() - assert all(x.exists() for x in outlier_files_path) @@ -284,9 +210,7 @@ def test_identical_images(tmp_path, base_image, caplog): result = outlier_step(input_models) # assert that log shows no new outliers detected - assert "New pixels flagged as outliers: 0 (0.00%)" in { - x.message for x in caplog.records - } + assert "0 pixels marked as outliers" in {x.message for x in caplog.records} # assert that DQ array has nothing flagged as outliers with result: for i, model in enumerate(result): diff --git a/romancal/outlier_detection/utils.py b/romancal/outlier_detection/utils.py new file mode 100644 index 000000000..56d3bbcce --- /dev/null +++ b/romancal/outlier_detection/utils.py @@ -0,0 +1,315 @@ +import copy +import logging +from functools import partial + +import numpy as np +from roman_datamodels.dqflags import pixel +from stcal.outlier_detection.median import MedianComputer +from stcal.outlier_detection.utils import ( + compute_weight_threshold, + flag_crs, + flag_resampled_crs, + gwcs_blot, +) + +from romancal.resample.resample import ResampleData +from romancal.resample.resample_utils import build_driz_weight + +from . import _fileio + +__all__ = ["detect_outliers"] + + +log = logging.getLogger(__name__) +log.setLevel(logging.DEBUG) + + +def _median_with_resampling( + input_models, + resamp, + maskpt, + save_intermediate_results, + make_output_path, + buffer_size=None, +): + """ + Compute median of resampled data from models in a library. + + Parameters + ---------- + input_models : ModelLibrary + The input datamodels. + + resamp : resample.resample.ResampleData object + The controlling object for the resampling process. + + maskpt : float + The weight threshold for masking out low weight pixels. + + save_intermediate_results : bool + if True, save the drizzled models and median model to fits. + + make_output_path : function + The functools.partial instance to pass to save_median. + + buffer_size : int + The size of chunk in bytes that will be read into memory when computing the median. + This parameter has no effect if the input library has its on_disk attribute + set to False. If None or 0 the buffer size will be set to the size of one + resampled image. + """ + if not resamp.single: + raise ValueError( + "median_with_resampling should only be used for resample_many_to_many" + ) + + in_memory = not input_models._on_disk + indices_by_group = list(input_models.group_indices.values()) + ngroups = len(indices_by_group) + example_model = None + median_wcs = resamp.output_wcs + + with input_models: + for i, indices in enumerate(indices_by_group): + + drizzled_model = resamp.resample_group(input_models, indices) + + if save_intermediate_results: + # write the drizzled model to file + _fileio.save_drizzled(drizzled_model, make_output_path) + + if i == 0: + input_shape = (ngroups,) + drizzled_model.data.shape + dtype = drizzled_model.data.dtype + computer = MedianComputer(input_shape, in_memory, buffer_size, dtype) + example_model = drizzled_model + + weight_threshold = compute_weight_threshold(drizzled_model.weight, maskpt) + drizzled_model.data[drizzled_model.weight < weight_threshold] = np.nan + computer.append(drizzled_model.data, i) + del drizzled_model + + # Perform median combination on set of drizzled mosaics + median_data = computer.evaluate() + + if save_intermediate_results: + # drizzled model already contains asn_id + _fileio.save_median( + example_model, + median_data, + median_wcs, + partial(make_output_path, asn_id=None), + ) + + return median_data, median_wcs + + +def _median_without_resampling( + input_models, + maskpt, + weight_type, + good_bits, + save_intermediate_results, + make_output_path, + buffer_size=None, +): + """ + Compute median of data from models in a library. + + Parameters + ---------- + input_models : ModelLibrary + The input datamodels. + + maskpt : float + The weight threshold for masking out low weight pixels. + + weight_type : str + The type of weighting to use when combining images. Options are: + 'ivm' (inverse variance) or 'exptime' (exposure time). + + good_bits : int + The bit values that are considered good when determining the + data quality of the input. + + save_intermediate_results : bool + if True, save the models and median model to fits. + + make_output_path : function + The functools.partial instance to pass to save_median. + + buffer_size : int + The size of chunk in bytes that will be read into memory when computing the median. + This parameter has no effect if the input library has its on_disk attribute + set to False. If None or 0 the buffer size will be set to the size of one + input image. + + """ + in_memory = not input_models._on_disk + ngroups = len(input_models) + example_model = None + + with input_models: + for i in range(len(input_models)): + + model = input_models.borrow(i) + wht = build_driz_weight( + model, + # FIXME this was hard-coded to "ivm" + weight_type="ivm", + good_bits=good_bits, + ) + + if save_intermediate_results: + # write the model to file + _fileio.save_drizzled(model, make_output_path) + + if i == 0: + input_shape = (ngroups,) + model.data.shape + dtype = model.data.dtype + computer = MedianComputer(input_shape, in_memory, buffer_size, dtype) + example_model = model + median_wcs = copy.deepcopy(model.meta.wcs) + + weight_threshold = compute_weight_threshold(wht, maskpt) + + data_copy = model.data.copy() + data_copy[wht < weight_threshold] = np.nan + computer.append(data_copy, i) + + input_models.shelve(model, i, modify=True) + del data_copy, model + + # Perform median combination on set of drizzled mosaics + median_data = computer.evaluate() + + if save_intermediate_results: + _fileio.save_median(example_model, median_data, median_wcs, make_output_path) + + return median_data, median_wcs + + +def _flag_resampled_model_crs( + image, + median_data, + median_wcs, + snr1, + snr2, + scale1, + scale2, + backg, + save_intermediate_results, + make_output_path, +): + blot = gwcs_blot(median_data, median_wcs, image.data.shape, image.meta.wcs, 1.0) + + # 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(image.meta, "background") + and image.meta.background.subtracted is False + and image.meta.background.level is not None + ): + backg = image.meta.background.level.value + log.debug( + f"Adding background level {image.meta.background.level} to blotted image" + ) + + cr_mask = flag_resampled_crs( + image.data.value, image.err.value, blot, snr1, snr2, scale1, scale2, backg + ) + + # update the dq flags in-place + image.dq |= cr_mask * np.uint32(pixel.DO_NOT_USE | pixel.OUTLIER) + log.info(f"{np.count_nonzero(cr_mask)} pixels marked as outliers") + + +def _flag_model_crs(image, median_data, snr): + cr_mask = flag_crs(image.data.value, image.err.value, median_data, snr) + + # Update dq array in-place + image.dq |= cr_mask * np.uint32(pixel.DO_NOT_USE | pixel.OUTLIER) + + log.info(f"{np.count_nonzero(cr_mask)} pixels marked as outliers") + + +def detect_outliers( + library, + weight_type, + pixfrac, + kernel, + fillval, + maskpt, + snr1, + snr2, + scale1, + scale2, + backg, + save_intermediate_results, + resample_data, + good_bits, + in_memory, + make_output_path, +): + # setup ResampleData + # call + if resample_data: + resamp = ResampleData( + library, + single=True, + blendheaders=False, + # FIXME prior code provided weight_type when only wht_type is understood + # both default to 'ivm' but tests that set this to something else did + # not change the resampling weight type. For now, disabling it to match + # the previous behavior. + # wht_type=weight_type + pixfrac=pixfrac, + kernel=kernel, + fillval=fillval, + good_bits=good_bits, + in_memory=in_memory, + ) + median_data, median_wcs = _median_with_resampling( + library, + resamp, + maskpt, + save_intermediate_results=save_intermediate_results, + make_output_path=make_output_path, + ) + else: + median_data, median_wcs = _median_without_resampling( + library, + maskpt, + weight_type, + good_bits, + save_intermediate_results=save_intermediate_results, + make_output_path=make_output_path, + ) + + # Perform outlier detection using statistical comparisons between + # each original input image and its blotted version of the median image + with library: + for image in library: + if resample_data: + _flag_resampled_model_crs( + image, + median_data, + median_wcs, + snr1, + snr2, + scale1, + scale2, + backg, + save_intermediate_results, + make_output_path, + ) + else: + _flag_model_crs(image, median_data, snr1) + + # mark step as complete + image.meta.cal_step["outlier_detection"] = "COMPLETE" + + library.shelve(image, modify=True) + + return library diff --git a/romancal/resample/resample.py b/romancal/resample/resample.py index ca6363f1d..55bc277be 100644 --- a/romancal/resample/resample.py +++ b/romancal/resample/resample.py @@ -1,4 +1,6 @@ +import json import logging +import os from typing import List import numpy as np @@ -8,6 +10,8 @@ from roman_datamodels import datamodels, maker_utils, stnode from stcal.alignment.util import compute_scale +from romancal.associations.asn_from_list import asn_from_list + from ..assign_wcs import utils from ..datamodels import ModelLibrary from . import gwcs_drizzle, resample_utils @@ -191,101 +195,118 @@ def do_drizzle(self): else: return self.resample_many_to_one() - def resample_many_to_many(self): - """Resample many inputs to many outputs where outputs have a common frame. + def resample_group(self, input_models, indices): + """ + Resample models at the provided indices. - Coadd only different detectors of the same exposure (e.g. map SCA 1 and - 10 onto the same output image), as they image different areas of the - sky. + This is used by outlier detection and will not blend metadata + and will not resample variance, error or exposure time. - Used for outlier detection + Parameters + ---------- + input_models : ModelLibrary + + indices : list + List of model indices to include in this resampling """ - output_list = [] - for group_id, indices in self.input_models.group_indices.items(): - output_model = self.blank_output - output_model.meta["resample"] = maker_utils.mk_resample() + output_model = self.blank_output.copy() + output_model.meta["resample"] = maker_utils.mk_resample() - # copy over asn information - if (asn_pool := self.input_models.asn.get("asn_pool", None)) is not None: - output_model.meta.asn.pool_name = asn_pool - if ( - asn_table_name := self.input_models.asn.get("table_name", None) - ) is not None: - output_model.meta.asn.table_name = asn_table_name + copy_asn_info_from_library(input_models, output_model) - with self.input_models: - example_image = self.input_models.borrow(indices[0]) - - # Determine output file type from input exposure filenames - # Use this for defining the output filename - indx = example_image.meta.filename.rfind(".") - output_type = example_image.meta.filename[indx:] - output_root = "_".join( - example_image.meta.filename.replace(output_type, "").split("_")[:-1] - ) - output_model.meta.filename = f"{output_root}_outlier_i2d{output_type}" + with input_models: + example_image = input_models.borrow(indices[0]) - self.input_models.shelve(example_image, indices[0], modify=False) + # Determine output file type from input exposure filenames + # Use this for defining the output filename + indx = example_image.meta.filename.rfind(".") + output_type = example_image.meta.filename[indx:] + output_root = "_".join( + example_image.meta.filename.replace(output_type, "").split("_")[:-1] + ) + output_model.meta.filename = f"{output_root}_outlier_i2d{output_type}" + input_models.shelve(example_image, indices[0], modify=False) + del example_image + + # Initialize the output with the wcs + driz = gwcs_drizzle.GWCSDrizzle( + output_model, + pixfrac=self.pixfrac, + kernel=self.kernel, + fillval=self.fillval, + ) - # Initialize the output with the wcs - driz = gwcs_drizzle.GWCSDrizzle( - output_model, - pixfrac=self.pixfrac, - kernel=self.kernel, - fillval=self.fillval, + log.info(f"{len(indices)} exposures to drizzle together") + for index in indices: + img = input_models.borrow(index) + # TODO: should weight_type=None here? + inwht = resample_utils.build_driz_weight( + img, weight_type=self.weight_type, good_bits=self.good_bits ) - log.info(f"{len(indices)} exposures to drizzle together") - for index in indices: - img = self.input_models.borrow(index) - # TODO: should weight_type=None here? - inwht = resample_utils.build_driz_weight( - img, weight_type=self.weight_type, good_bits=self.good_bits - ) + # apply sky subtraction + if ( + hasattr(img.meta, "background") + and img.meta.background.subtracted is False + and img.meta.background.level is not None + ): + data = img.data - img.meta.background.level + else: + data = img.data - # apply sky subtraction - if ( - hasattr(img.meta, "background") - and img.meta.background.subtracted is False - and img.meta.background.level is not None - ): - data = img.data - img.meta.background.level - else: - data = img.data - - xmin, xmax, ymin, ymax = resample_utils.resample_range( - data.shape, img.meta.wcs.bounding_box - ) + xmin, xmax, ymin, ymax = resample_utils.resample_range( + data.shape, img.meta.wcs.bounding_box + ) - driz.add_image( - data, - img.meta.wcs, - inwht=inwht, - xmin=xmin, - xmax=xmax, - ymin=ymin, - ymax=ymax, - ) - del data - self.input_models.shelve(img, index) - - # cast context array to uint32 - output_model.context = output_model.context.astype("uint32") - - # copy over asn information - if not self.in_memory: - # Write out model to disk, then return filename - output_name = output_model.meta.filename - output_model.save(output_name) - log.info(f"Exposure {output_name} saved to file") - output_list.append(output_name) - else: - output_list.append(output_model.copy()) + driz.add_image( + data, + img.meta.wcs, + inwht=inwht, + xmin=xmin, + xmax=xmax, + ymin=ymin, + ymax=ymax, + ) + del data + self.input_models.shelve(img, index, modify=False) + del img - output_model.data *= 0.0 - output_model.weight *= 0.0 + # cast context array to uint32 + output_model.context = output_model.context.astype("uint32") + return output_model - return ModelLibrary(output_list) + def resample_many_to_many(self): + """Resample many inputs to many outputs where outputs have a common frame. + + Coadd only different detectors of the same exposure (e.g. map SCA 1 and + 10 onto the same output image), as they image different areas of the + sky. + """ + output_models = [] + for group_id, indices in self.input_models.group_indices.items(): + output_model = self.resample_group(self.input_models, indices) + + if not self.in_memory: + # Write out model to disk, then return filename + output_name = output_model.meta.filename + if self.output_dir is not None: + output_name = os.path.join(self.output_dir, output_name) + output_model.save(output_name) + log.info(f"Saved model in {output_name}") + output_models.append(output_name) + else: + output_models.append(output_model) + + if not self.in_memory: + # build ModelLibrary as an association from the output files + # this saves memory if there are multiple groups + asn = asn_from_list(output_models, product_name="outlier_i2d") + asn_dict = json.loads( + asn.dump()[1] + ) # serializes the asn and converts to dict + return ModelLibrary(asn_dict, on_disk=True) + # otherwise just build it as a list of in-memory models + return ModelLibrary(output_models, on_disk=False) def resample_many_to_one(self): """Resample and coadd many inputs to a single output. @@ -946,3 +967,11 @@ def populate_mosaic_individual( input_metas = [datamodel.meta for datamodel in input_models] for input_meta in input_metas: output_model.append_individual_image_meta(input_meta) + + +def copy_asn_info_from_library(input_models, output_model): + # copy over asn information + if (asn_pool := input_models.asn.get("asn_pool", None)) is not None: + output_model.meta.asn.pool_name = asn_pool + if (asn_table_name := input_models.asn.get("table_name", None)) is not None: + output_model.meta.asn.table_name = asn_table_name