Pre-processing the data

The codes in STDPipe expect as an input the science-ready images, cleaned as much as possible from instrumental signatures and imaging artefacts. In practice, it means that the image should be

  • bias and dark subtracted

  • flat-fielded.

Also, the artefacts such as saturated stars, bleeding charges, cosmic ray hits etc have to be masked.

All these tasks are outside of STDPipe per se, as they are highly instrument and site specific. On the other hand, they may usually be performed easily using standard Python/NumPy/AstroPy routines and libraries like Astro-SCRAPPY etc.

E.g. to pre-process the raw image using pre-computed master dark and flat, mask some common problems, and cleanup the cosmic rays you may do something like that:

image = fits.getdata(filename).astype(np.double)
header = fits.getheader(filename)

dark = fits.getdata(darkname)
flat = fits.getdata(flatname)

image -= dark
image *= np.median(flat)/flat

saturation = header.get('SATURATE') or 50000 # Guess saturation level from FITS header

mask = np.isnan(image) # mask NaNs in the input image
mask |= image > saturation # mask saturated pixels
mask |= flat < 0.5*np.median(flat) # mask underilluminated/vignetted regions

from astropy.stats import mad_std
mask |= dark > np.median(dark) + 10.0*mad_std(dark) # mask hotter pixels

gain = header.get('GAIN') or 1.0 # Guess gain from FITS header

import astroscrappy
# mask cosmic rays using LACosmic algorithm
cmask,cimage = astroscrappy.detect_cosmics(image, mask, gain=gain, verbose=True)
mask |= cmask

We have a simple routine that implements these steps, stdpipe.pipeline.make_mask(), which is documented below.

stdpipe.pipeline.make_mask(image, header=None, saturation=None, external_mask=None, mask_cosmics=False, gain=None, verbose=True)[source]

Make a basic mask for an image.

The mask is a boolean bitmap with True values marking regions to be excluded from further processing. The following regions are masked:

  • pixels with undefined or non-finite (inf or nan) values

  • regions outside the usable area defined by the DATASEC or TRIMSEC header keyword

  • pixels above the saturation limit, if provided

  • pixels set in the external mask, if provided

  • cosmic rays, if requested

If saturation=True, the saturation level is estimated from the image as median + 0.95 * (max - median).

Parameters:
imagendarray

2D image array to mask.

headerastropy.io.fits.Header, optional

FITS header; used to read DATASEC / TRIMSEC keywords.

saturationfloat or bool, optional

Saturation level in ADU. If True, estimated automatically from the image.

external_maskndarray, optional

Boolean mask to OR with the created mask.

mask_cosmicsbool, optional

If True, detect and mask cosmic rays using the LACosmic algorithm.

gainfloat, optional

Detector gain in e-/ADU, used for cosmic-ray masking.

verbosebool or callable, optional

Whether to show verbose messages. May be boolean or a print-like callable.

Returns:
ndarray of bool

Boolean mask with True marking excluded pixels.

Removing fringes

Images taken in red and near-infrared bands with thinned CCDs often show fringes - a pattern of curved, wave-like ripples caused by interference of night-sky emission lines inside the detector. The pattern is additive and varies smoothly across the field on scales of tens of pixels, sitting between the size of the stars (a few pixels) and of the large-scale sky background (hundreds of pixels). As such, it is too fine to be captured by the usual grid-based background estimators, and too coarse and irregular to be removed by frequency-domain (Fourier) filtering, as real fringes are typically curved and non-stationary rather than strictly periodic.

STDPipe provides stdpipe.fringe_removal.remove_fringes() that models the fringe pattern directly in image space. It masks the stars and other significant sources, estimates the smooth intermediate-scale residual by mask-aware Gaussian smoothing (so the fringe structure is transparently inpainted under the masked positions), and refines this estimate over a few iterations. The resulting fringe map is then subtracted from the image, while the large-scale sky background is left intact.

from stdpipe import fringe_removal

# Remove the fringes, automatically selecting the smoothing scale
cleaned = fringe_removal.remove_fringes(image, mask=mask)

# Also return the model that was subtracted, e.g. for inspection
cleaned, fringe_model = fringe_removal.remove_fringes(
    image, mask=mask, get_fringe_model=True)

The smoothing kernel size is the main parameter. By default (scale='auto') it is selected automatically from the measured fringe-to-noise ratio: strong, sharp fringes get a small kernel that follows them closely, while faint, broad fringes get a larger one that averages the pixel noise down instead of leaking it into the model. You may set scale to a number to fix it manually, and use noise_frac to trade model cleanliness against fringe-tracking fidelity when using the automatic mode.

# Fix the smoothing scale manually (pixels)
cleaned = fringe_removal.remove_fringes(image, mask=mask, scale=10)

# Make the automatic selection less aggressive (cleaner model, smaller kernel)
cleaned = fringe_removal.remove_fringes(image, mask=mask, noise_frac=0.4)

Bright stars carry azimuthally symmetric structures - extended halos, ghost reflection rings, and bowls left by the coarse background estimation - that are locally indistinguishable from fringes. By default (halo_sn=20) these are modelled with radial profiles around bright stars and removed together with the fringes; pass halo_sn=None to disable this and keep the faint stellar wings in the image.

Attention

The fringe model is an additive correction estimated from the science image itself. It assumes the fringes are faint compared to the stars and vary on scales well above the stellar size. For fields dominated by extended low-surface-brightness sources (large galaxies, nebulosity), part of their flux may be absorbed into the fringe model - mask such regions explicitly via the mask argument in that case.

The standalone script examples/fringe_removal_example.py runs the routine on a FITS image and visualises the original image, the reconstructed fringe pattern and the corrected result side by side, exposing all parameters as command-line options - useful for inspecting the correction and tuning it for a given instrument.

stdpipe.fringe_removal.remove_fringes(image, mask=None, scale='auto', noise_frac=0.25, bg_size=256, threshold=2.0, dilate=2, iterations=3, halo_sn=20.0, get_fringe_model=False, verbose=False)[source]

Remove fringes by estimating a smooth intermediate-scale background map with sources masked and inpainted.

The algorithm exploits the scale separation between stars (~FWHM, a few pixels), fringes (tens of pixels) and the sky background (hundreds of pixels):

  1. Subtract a coarse sky background (SEP mesh of bg_size pixels)

  2. Around bright stars (peaks above halo_sn sigmas with footprints wider than the smoothing kernel), model the non-monotone azimuthally symmetric circumstellar structure - ghost reflection rings, background estimation bowls - as azimuthal-median radial profiles with their smooth monotone envelope (the stellar wings and halo) removed, and subtract it from the working residual. This captures circular artifacts with full circular fidelity, which isotropic smoothing cannot do

  3. Mask pixels deviating by more than threshold sigma (sources and outliers), grow the mask by dilate pixels

  4. Smooth the masked residual with a Gaussian of sigma scale pixels using mask-aware normalized convolution - masked pixels are filled by the weighted average of the surrounding fringe structure (inpainting), each region by the finest scale of a multi-scale ladder that still has unmasked support there. The smoothing also absorbs the smooth sub-threshold envelope of stellar halos and wings

  5. Iterate: re-mask on the fringe-subtracted residual (so that fringe crests are not mistaken for sources), smooth it again and add to the fringe map. The additive refinement also recovers the part of the fringe structure attenuated by the smoothing, so narrower fringes are progressively captured with each iteration

The radial circumstellar fields are included in the returned fringe map and thus subtracted from the image along with the fringes.

The resulting fringe map is defined everywhere, including under the stars, and is subtracted from the original image. The coarse sky background is left in the image.

Parameters:
imagendarray

Input image

maskndarray, optional

Boolean mask (True = masked pixels to exclude from fringe estimation)

scalefloat or ‘auto’

Sigma of the Gaussian smoothing kernel in pixels. Should be of order the stellar FWHM or somewhat larger, and well below the narrowest fringe period. Smaller values track narrower fringes at the cost of slightly more noise and source flux absorbed into the map. If 'auto' (default), the scale is selected from the measured fringe-to-noise ratio (6 px for strong, sharp fringes up to a few tens of px for faint, broad ones), so that noise does not leak into the fringe model. Pass a number to set it explicitly.

noise_fracfloat

Aggressiveness of the automatic scale selection (used only when scale='auto'). The smallest scale whose predicted noise floor falls below noise_frac times the fringe amplitude is chosen. Smaller values give cleaner (less noisy) models but larger kernels, which can over-smooth very faint, broad fringes; larger values keep smaller kernels at the cost of more noise in the model. Default: 0.25

bg_sizeint or None

Mesh size in pixels for the coarse sky background subtracted before fringe estimation. Should be much larger than the fringe scale. If None, the global median is used instead. Default: 256

thresholdfloat

Source masking threshold in sigmas of the local background RMS. Applied symmetrically (both positive and negative outliers are masked) to avoid biasing the fringe map. Default: 2.0

dilateint

Number of binary dilation iterations applied to the source mask to cover the wings of the sources. Default: 2

iterationsint

Number of refinement iterations. Each iteration re-derives the source mask from the fringe-subtracted residual (the first one has to mask on the raw residual where fringe crests may exceed the threshold) and accumulates the smoothed residual into the fringe map. More iterations capture narrower fringes but absorb more noise. Default: 3

halo_snfloat or None

Peak significance (in sigmas) above which a source is considered bright enough to have halos and gets the radial halo treatment. Should be well above the peak significance of the fringe crests. Only sources whose thresholded footprint radius also exceeds scale are treated (ordinary stars have compact footprints and are not affected). None disables the halo modeling. Default: 20

get_fringe_modelbool

If True, return (corrected_image, fringe_model). Default: False

verbosebool or callable

Whether to show verbose messages during the run. May be either boolean, or a print-like function.

Returns:
corrected_imagendarray

Image with the fringe map subtracted

fringe_modelndarray, optional

The subtracted fringe map (if get_fringe_model=True)

Notes

Circumstellar structures - halos, ghost rings, background bowls - are locally indistinguishable from fringes (smooth structures at similar scales). The smooth ones are absorbed by the masked multi-scale smoothing along with the fringes; sharply circular ones (ghost donuts, bowls) are additionally captured with full circular fidelity by azimuthal-median radial profiles around bright stars (the halo_sn mechanism), which average the oscillating fringes passing through to nearly zero. Everything smooth or circular around bright stars is thus subtracted, including the sub-threshold part of their PSF wings; the star cores and everything above the masking threshold stay in the image. Use halo_sn=None to disable the radial modeling (isotropic smoothing only).

Faint extended sources (low surface brightness galaxies, nebulae) whose peaks stay below halo_sn remain partially degenerate with the fringe pattern: the part below the masking threshold gets absorbed into the fringe map and over-subtracted. For fields dominated by such sources, lower halo_sn, increase scale, or mask the sources explicitly via the mask argument.

Examples

>>> corrected = remove_fringes(image, mask=mask, scale=6)
>>> corrected, fringes = remove_fringes(image, get_fringe_model=True)

Stacking the images

You may want to stack/coadd or mosaic some images before processing them. While there are dedicated large-scape packages like Montage that handle it properly, it still may be done manually with relatively little efforts using e.g. Python reproject package.

You may check simple example notebook that shows how to do it.

Attention

The stacking modify the statistical properties of resulting image! The reasons are both averaging (or especially median averaging!) of the images that modify effective gain value (typically increasing it by the factor equal to number of averaged images), and pixel interpolation when re-projecting the images onto the same pixel grid.

STDPipe provides two methods for image reprojection:

SWarp reprojection

Alternatively, stdpipe.reproject.reproject_swarp() wraps the SWarp external binary. It is implemented to resemble the calling conventions of the reproject package - i.e. allows directly stacking image files without loading them to memory first. Requires SWarp to be installed.

stdpipe.reproject.reproject_swarp(input=None, wcs=None, shape=None, width=None, height=None, header=None, extra=None, is_flags=False, use_nans=True, get_weights=False, _workdir=None, _tmpdir=None, _exe=None, verbose=False)[source]

Wrapper for running SWarp for re-projecting and mosaicking of images onto target WCS grid.

It accepts as input either list of filenames, or list of tuples where first element is an image, and second one - either FITS header or WCS.

If the input images are integer flags, set is_flags=True so that it will be handled by passing RESAMPLING_TYPE=FLAGS and COMBINE_TYPE=AND.

If use_nans=True, the regions with zero weights will be filled with NaNs (or 0xFFFF).

Any additional configuration parameter may be passed to SWarp through extra argument which should be the dictionary with parameter names as keys.