imcombine

imutil.imcombine.imcombine(
    inputs,
    mask=None,
    extension=None,
    extension_uncertainty=None,
    extension_mask=None,
    uncertainty_type='stddev',
    trimsec=None,
    blank=np.nan,
    offsets=None,
    thresholds=None,
    zero=None,
    zero_to_0th=True,
    zero_section=None,
    scale=None,
    scale_to_0th=True,
    scale_section=None,
    zero_kw=None,
    scale_kw=None,
    weight=None,
    reject=None,
    sigma=None,
    cenfunc='median',
    maxiters=50,
    ddof=1,
    nkeep=1,
    maxrej=None,
    n_minmax=None,
    rdnoise=0.0,
    gain=1.0,
    snoise=0.0,
    pclip=-0.5,
    combine='average',
    dtype='float32',
    dtype_std='float32',
    dtype_low=None,
    dtype_upp=None,
    memlimit=2500000000.0,
    verbose=False,
    full=False,
    imcmb_key='$I',
    exposure_key='EXPTIME',
    output=None,
    output_mask=None,
    output_nrej=None,
    output_std=None,
    output_low=None,
    output_upp=None,
    output_flags=None,
    return_dict=False,
    output_verify='exception',
    overwrite=False,
    checksum=False,
    *,
    diagnostics=None,
    grow=None,
    output_sample_flags=None,
)

A FITS-file helper for imcombiners.ndcombine.

Warning

A few options are not implemented yet: blank, statsec with input/output/overlap, scale_sample, zero_sample, "mode" for scale/zero/weight.

Parameters

inputs : glob pattern, list-like of path-like, list-like of ~astropy.nddata.CCDData-like

The ~glob pattern for files (e.g., "2020*[012].fits") or list of files (each element must be path-like or ~astropy.nddata.CCDData). Although it is not a good idea, a mixed list of ~astropy.nddata.CCDData and paths to the files is also acceptable. For the purpose of imred.imcombine, the best use is to use the ~glob pattern or list of paths.

mask : ~numpy.ndarray, optional. = None

The mask of bad pixels. If given, it must satisfy mask.shape[0] identical to the number of images. Each plane uses the corresponding input’s coordinates after trimsec, before offsets. It is combined with CCDData masks or the selected FITS mask extension for normalization statistics as well as rejection and combination.

Note

If the user ever want to use masking, it’s more convenient to use 'MASK' extension to the FITS files or replace bad pixel to very large or small numbers and use thresholds.

extension : int, str, (str, int) = None

The extension of FITS, uncertainty, and mask to be used. It can be given as integer (0-indexing) of the extension, EXTNAME (single str), or a tuple of str and int: (EXTNAME, EXTVER). If None (default), the first extension with data will be used. If extension_uncertainty or extension_mask is None (default), uncertainty and mask extensions are ignored (turned off). A non-None extension_uncertainty raises NotImplementedError before input I/O: uncertainty propagation is not supported. Attached CCDData uncertainties are ignored without copying. output_std is clipping spread, not propagated measurement uncertainty.

extension_uncertainty : int, str, (str, int) = None

extension_mask : int, str, (str, int) = None

offsets : str or (n, m)-d array = None

If array, it must have shape such that n is the number of images and m is the dimension of the images (if m=3, offsets in x, y, z, … order, not pythonic order), and it is directly regarded as the raw offsets. If str, the raw offsets are obtained by the followings:

- ``CRPIX`` values in the header if ``"wcs"|"world"``
- ``LTV`` values in the header if ``"physical"|"phys"|"phy"``

Warning

The physical coordinate system is defined by the IRAF-like LTM/LTV keywords define the offsets. Currently, only the cases when LTMi_j is 0 or 1 can be managed. Otherwise, we need scaling and it is not supported now.

For both wcs or physical cases, the raw offsets for each frame is nothing but an m-D tuple consists of offset_raw[i] = CRPIX{m-i} or LTV{m-i}[_{m-i}]. The reason to subtract i is because python has z, y, x order of indexing while ~astropy.wcs.WCS information is in x, y, z order. If it is a j-th image, offsets[j, :] = offset_raw, and offsets has shape of (n, m).

This raw offsets are then modified such that the minimum offsets in each axis becomes zero (in pythonic way, np.max(offsets, axis=0) - offsets). The offsets are used to determine the final output image’s shape.

Note

Though the IRAF IMCOMBINE manual says it calculates offsets from the 0-th image center if offsets="wcs", it seems it acutally uses CRPIX from the header… I couldn’t find how IRAF does offset calculation for astropy.wcs.WCS, it’s not reproducible using rounding. Even using astropy.wcs.WCS info correctly, it’s not reproducible. Also, if we only use CRPIX, the offset calculations are completely wrong if CRPIX is not centered at the identical world coordinate (e.g., RA/DEC). IRAF indeed wrongly combines images if this happens.

thresholds : 2-float list-like, optional. = None

The thresholds (lower, upper) applied to all images before any rejection/combination. Default is no thresholding, (-np.inf, +np.inf). One possible usage is to replace bad pixel to very large or small numbers and use this thresholding.

zero : str or 1-d array = None

The zero value to subtract from each image after thresholding, but before scaling/offset shifting/rejection/combination. If an array, it is directly subtracted from each image, (so it must have size identical to the number of images). If str, the zero-level is:

- **Average** if ``'avg'|'average'|'mean'``
- **Median** if ``'med'|'medi'|'median'``
- **Sigma-clipped average** if
  ``'avg_sc'|'average_sc'|'mean_sc'``
- **Sigma-clipped median** if
  ``'med_sc'|'medi_sc'|'median_sc'``

For options for sigma-clipped statistics, see zero_kw.

Note

By using zero of "med_sc", the user can crudely subtract sky value from each frame before combining.

scale : str or 1-d array = None

The way to scale each image after thresholding/zeroing, but before offset shifting/rejection/combination. If an array, it is directly understood as the raw scales, and it must have size identical to the number of images. If str, the raw scale is:

- **Exposure time** (``exposure_key`` in header of each FITS)
  if ``'exp'|'expos'|'exposure'|'exptime'``
- **Average** if ``'avg'|'average'|'mean'``
- **Median** if ``'med'|'medi'|'median'``
- **Sigma-clipped average** if
  ``'avg_sc'|'average_sc'|'mean_sc'``
- **Sigma-clipped median** if
  ``'med_sc'|'medi_sc'|'median_sc'``

The true scale is obtained by scales / scales[0] if scale_to_0th is True, following IRAF’s convention. Otherwise, the absolute value from the raw scale will be used. For options for sigma-clipped statistics, see scale_kw.

Note

Using scale="avg_sc", scale_to_0th=False is useful for flat combining.

zero_to_0th : bool, optional. = True

Whether to re-base the zero values such that all images have identical zero values as that of the 0-th image (in python, zero - zero[0]). This is the behavior of IRAF, so zero_to_0th is True by default.

scale_to_0th : bool, optional. = True

Whether to re-scale the scales such that scale[0] is unity (in python, scale/scale[0]). This is the behavior of IRAF, so scale_to_0th is True by default.

zero_section : str, optional. = None

The sections used for zeroing and scaling. These must be in FITS section format, and are the sections AFTER trimming based on trimsec. Statistical zero/scale/weight values exclude input-masked and nonfinite pixels. A section with no finite unmasked pixels raises ValueError naming the parameter and input index. Numeric calibration values are unaffected. Thresholds and algorithm rejections do not select normalization samples. Callable statistics receive the selected pixels as an (M, 1) array with excluded samples replaced by NaN; use a NaN-aware statistic and return one finite value (nonzero for scale or weight).

scale_section : str, optional. = None

zero_kw : dict = None

Used only if scale or zero are sigma-clipped mean, median, etc (ending with _sc such as median_sc, avg_sc). The keyword arguments for astropy.stats.sigma_clipped_stats. By default, std_ddof=1 (note that ~astropy.stats.sigma_clipped_stats has default std_ddof=0.)

Note

If axis is specified, it will be ignored.
Sigma-clipping in astropy has cumulative rejection algorithm by default. This may give unwanted results especially when sigma in this zero_kw or scale_kw is small (around 1), but this is not likely a problem since both zero and scale will be determined from all the pixels, which is usually more than an order of a million.

scale_kw : dict = None

Note

If axis is specified, it will be ignored.
Sigma-clipping in astropy has cumulative rejection algorithm by default. This may give unwanted results especially when sigma in this zero_kw or scale_kw is small (around 1), but this is not likely a problem since both zero and scale will be determined from all the pixels, which is usually more than an order of a million.

sigma : 2-float list-like, optional. = None

The sigma-factors to be used for sigma-clip rejeciton in (sigma_lower, sigma_upper). Defaults to (3, 3), which means 3-sigma clipping from the “sigma” values determined by the method specified by reject. If a single float, it will be used for both the lower and upper values.

maxiters : int, optional. = 50

The maximum number of iterations to do the rejection (for sigma-clipping). It is silently converted to int if it is not.

ddof : int, optional. = 1

The delta-degrees of freedom (see numpy.std). It is silently converted to int if it is not.

nkeep : float or int, optional. = 1

The minimum number of pixels that should be left after rejection. If nkeep < 1, it is regarded as fraction of the total number of pixels along the axis to combine. This corresponds to positive nkeep parameter of IRAF IMCOMBINE. If number of remaining non-nan value is fewer than nkeep, the masks at that position will be reverted to the previous iteration, and rejection code will be added by number 4.

maxrej : float or int, optional. = None

The maximum number of pixels that can be rejected during the rejection. If maxrej < 1, it is regarded as fraction of the total number of pixels along the axis to combine. This corresponds to negative nkeep parameter of IRAF IMCOMBINE. In IRAF, only one of nkeep and maxrej can be set. If number of rejected pixels at a position exceeds maxrej, the masks at that position will be reverted to the previous iteration, and rejection code will be added by number 8.

cenfunc : str, optional. = 'median'

The centering function to be used in rejection algorithm.

- median if  ``'med'|'medi'|'median'``
- average if ``'avg'|'average'|'mean'``
- lower median if ``'lmed'|'lmd'|'lmedian'``

For lower median, see note in combine.

n_minmax : 2-float or 2-int list-like, optional. = None

The number of low and high pixels to be rejected by the “minmax” algorithm. These numbers are converted to fractions of the total number of input images so that if no rejections have taken place the specified number of pixels are rejected while if pixels have been rejected by masking, thresholding, or non-overlap, then the fraction of the remaining pixels, truncated to an integer, is used.

rdnoise : float, optional. = 0.0

The readnoise of the detector in the unit of electrons, electron gain of the detector in the unit of elctrons/DN (or electrons/ADU), and sensitivity noise as a fraction. Used only if reject="ccdclip" and/or combine="nmodel".

The variance of a single pixel in an image when these are used,

\[ V_\mathrm{DN} = ( \mathtt{rdnoise}/\mathtt{gain} )^2 + \mathrm{DN}/\mathtt{gain} + ( \mathtt{snoise} * \mathrm{DN} )^2 \]

\[ V_\mathrm{electron} = (\mathtt{rdnoise})^2 + (\mathtt{gain} * \mathrm{DN})^2 + (\mathtt{snoise} * \mathtt{gain} * \mathrm{DN})^2 \]

gain : float, optional. = 0.0

The variance of a single pixel in an image when these are used,

\[ V_\mathrm{DN} = ( \mathtt{rdnoise}/\mathtt{gain} )^2 + \mathrm{DN}/\mathtt{gain} + ( \mathtt{snoise} * \mathrm{DN} )^2 \]

\[ V_\mathrm{electron} = (\mathtt{rdnoise})^2 + (\mathtt{gain} * \mathrm{DN})^2 + (\mathtt{snoise} * \mathtt{gain} * \mathrm{DN})^2 \]

snoise : float, optional. = 0.0

The variance of a single pixel in an image when these are used,

\[ V_\mathrm{DN} = ( \mathtt{rdnoise}/\mathtt{gain} )^2 + \mathrm{DN}/\mathtt{gain} + ( \mathtt{snoise} * \mathrm{DN} )^2 \]

\[ V_\mathrm{electron} = (\mathtt{rdnoise})^2 + (\mathtt{gain} * \mathrm{DN})^2 + (\mathtt{snoise} * \mathtt{gain} * \mathrm{DN})^2 \]

pclip : float, optional. = -0.5

The parameter for reject="pclip". If abs(pclip) >= 1, then it specifies a number of pixels above or below the median to use for computing the clipping sigma. If abs(pclip) < 1, then it specifies the fraction of the pixels above or below the median to use. A positive value selects a point above the median and a negative value selects a point below the median. The default of -0.5 selects approximately the quartile point. Better to use negative value to avoid cosmic-ray contamination.

combine : str = 'average'

The function to be used for the final combining after thresholding, zeroing, scaling, rejection, and offset shifting.

- median if  ``'med'|'medi'|'median'``
- average if ``'avg'|'average'|'mean'``
- lower median if ``'lmed'|'lmd'|'lmedian'``
- sum if ``'sum'``
- min/max if ``'min'`` or ``'max'``
- variance if ``'variance'|'var'``
- weighted average if ``'weighted_average'|'wvg'``

Note

The lower median means the median which takes the lower value when even number of data is left. This is suggested to be robust against cosmic-ray hit according to the IRAF IMCOMBINE manual. The combine/rejection kernels are provided by imcombiners.

weight : float, array-like, str, or None = None

Finite, nonzero weights for mean/average combination; signs are retained. Use a shared scalar, one value per image, or a statistic such as "mean" computed within scale_section on each trimmed input. Statistics exclude input-masked and nonfinite pixels, as for zero and scale. A callable receives an (M, 1) array with excluded samples set to NaN and must return a finite, nonzero value. Vectors follow sorted file-path order; CCDData lists retain input order. Each pixel uses only finite, unmasked, unrejected samples; None gives equal weights. No valid samples gives NaN; valid weights summing to zero raise ZeroDivisionError. Weights do not affect rejection.

diagnostics : str or None = None

None preserves the existing full and output-path behavior. "simple" returns the same diagnostics as full=True; "full" also returns sample_flags, costing one retained byte per input sample plus temporary workspace. Flags do not change pixel values.

grow : float or None = None

Non-negative Euclidean radius in pixels around algorithm rejections, within each input image. Input masks, thresholds, and nonfinite pixels are not growth seeds. None disables growth; radii below 1 add nothing. Chunked reads include neighbors across boundaries. Growth uses extra temporary diagnostics and a second combination pass. With no reject, it has no effect. std remains clipping spread, not propagated error.

imcmb_key : str = '$I'

The thing to add as IMCMBnnn in the output FITS file header. If "$I", following the default of IRAF, the file’s name will be added. Otherwise, it should be a header keyword. If the key does not exist in nnn-th file, a null string will be added. If a null string (imcmb_key=""), it does not set the IMCMBnnn keywords nor deletes any existing keyword.

Warning

If more than 999 files are combined, only the first 999 files will be recorded in the header.

exposure_key : str, optional. = 'EXPTIME'

The header keyword which contains the information about the exposure time of each FITS file. This is used only if scaling is done for exposure time (see scale).

uncertainty_type : str = 'stddev'

Reserved for future uncertainty propagation. Only the default "stddev" is accepted; other values raise NotImplementedError.

memlimit : float or None = 2500000000.0

Approximate byte budget for input data, working buffers, and outputs. If the estimate exceeds it, read and combine spatial chunks. The final image and requested diagnostics remain in memory. Statistical normalization may still read a full input image. None or non-positive values disable chunking.

output : path - like = None

Write the final combined image here. dtype controls the output data type; offsets can change the output shape.

output_mask : path - like = None

FITS paths for the combined mask, rejection count, spread, rejection bounds, output flags, and per-sample flags, respectively. Each None default skips writing that diagnostic. See Returns for their meanings.

output_nrej : path - like = None

output_std : path - like = None

output_low : path - like = None

output_upp : path - like = None

output_flags : path - like = None

output_sample_flags : path - like = None

Write the stack-shaped uint8 sample flags as a FITS image and enable detailed diagnostic returns. Uses overwrite, checksum, and output_verify like the other output files.

return_dict : bool, optional. = False

When diagnostics are enabled, True returns a dict and False a tuple. Without diagnostics, the result is always a single CCDData.

Returns

comb : astropy.nddata.CCDData (dtype dtype): The combined data. Returned alone without diagnostics; otherwise it is the first tuple entry or the "comb" dictionary value.
std : ~numpy.ndarray or None: The spread diagnostic returned by imcombiners. It is available for sigma/CCD clipping and is None for rejection algorithms without a spread diagnostic.
: mask_rej, mask_thresh : ~numpy.ndarray or None: The masks from the rejection process and thresholding process (thresholds). Thresholding is done prior to rejection or zero/scale normalization. None means that diagnostic category was not requested or not applicable.
: low, upp : ~numpy.ndarray or None: The lower and upper retained-value bounds reported by imcombiners.
nit : ~numpy.ndarray or None: The number of rejection iterations.
output_flags : ~numpy.ndarray or None: Integer output flags reported by imcombiners.
sample_flags : ~numpy.ndarray of uint8: Returned only for diagnostics="full" or output_sample_flags: appended as the ninth tuple entry or the "sample_flags" dict key. Shape is (number_of_inputs, *combined_shape). Bits follow imcombiners.SampleFlags: 1=input mask, 2=nonfinite, 4=threshold, 8=algorithm rejection, 16=added by growth, 32=previous exclusion, 64=restored by nkeep, and 128=restored by maxrej. Multiple bits may be set; this is not a boolean exclusion mask. output_flags uses a separate per-output-pixel flag namespace; its bit 16 indicates that growth added at least one rejection there.