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
Truevalues 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
DATASECorTRIMSECheader keywordpixels 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 asmedian + 0.95 * (max - median).- Parameters:
- imagendarray
2D image array to mask.
- headerastropy.io.fits.Header, optional
FITS header; used to read
DATASEC/TRIMSECkeywords.- 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
Truemarking excluded pixels.
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:
Lanczos reprojection (recommended)¶
The default method stdpipe.reproject.reproject_lanczos() is a pure Python implementation of Lanczos interpolation with two features borrowed from SWarp that are critical for photometric accuracy:
Automatic oversampling - when output pixels are larger than input pixels (downscaling), the interpolation is evaluated at multiple sub-pixel positions and averaged. This prevents aliasing of undersampled stars.
Jacobian flux conservation - the output is multiplied by the pixel area ratio (Jacobian determinant) so that total flux is conserved. Set
conserve_flux=Falseto preserve surface brightness instead.
It also supports reprojection of integer flag/mask images (is_flags=True) using nearest-neighbor resampling and bitwise AND combining, matching SWarp’s behavior.
No external binaries are required.
from stdpipe.reproject import reproject_lanczos
# Reproject and coadd a list of FITS files
coadd = reproject_lanczos(filenames, wcs=target_wcs, shape=(1024, 1024))
# Reproject from (image, WCS) tuples
coadd = reproject_lanczos([(image1, wcs1), (image2, wcs2)],
wcs=target_wcs, shape=(1024, 1024))
# Preserve surface brightness instead of total flux
coadd = reproject_lanczos([(image, wcs_in)], wcs=wcs_out,
shape=(512, 512), conserve_flux=False)
# Reproject integer flag/mask image
mask = reproject_lanczos([(flags, wcs_in)], wcs=wcs_out,
shape=(512, 512), is_flags=True)
- stdpipe.reproject.reproject_lanczos(input=None, wcs=None, shape=None, width=None, height=None, header=None, order=3, conserve_flux=True, oversamp=None, is_flags=False, use_nans=True, parallel=False, return_footprint=False, verbose=False)[source]
Reproject images using Lanczos interpolation with automatic oversampling.
Implements SWarp-style oversampling (sub-pixel averaging when output pixels are larger than input pixels) and Jacobian area scaling for flux conservation.
Accepts the same input format as
reproject_swarp(): a list of(image, header/WCS)tuples or a list of FITS filenames. For multiple inputs the reprojected frames are averaged (simple coadd).- Parameters:
- inputlist or tuple
List of
(image, header_or_wcs)tuples or FITS filenames. A single(image, header_or_wcs)tuple is also accepted (wrapped into a list automatically for reproject compatibility).- wcs
~astropy.wcs.WCS, optional Output WCS. Ignored if header already contains WCS.
- shapetuple, optional
Output
(height, width).- width, heightint, optional
Output dimensions (alternative to shape).
- header
~astropy.io.fits.Header, optional Output FITS header (overrides wcs/shape/width/height).
- orderint
Lanczos kernel order (default 3).
- conserve_fluxbool
If True (default), multiply by the Jacobian area ratio so that total flux is conserved. If False, surface brightness is conserved instead.
- oversampint or None
Sub-pixel oversampling factor per axis.
None(default) selects automatically:max(1, round(output_scale / input_scale)).- is_flagsbool
If True, treat input as integer flag/mask images: use nearest-neighbor resampling (no interpolation) and bitwise AND for combining multiple inputs. Overrides order, conserve_flux and oversamp.
- use_nansbool
If True (default), fill regions with no input coverage with NaN (floating-point images) or
0xFFFF(integer flag images).- parallelbool or int
If True, use threads for parallel interpolation (number chosen automatically). If int > 1, use that many threads. Gives ~3-4x speedup on multi-core machines.
- return_footprintbool
If True, return
(coadd, footprint)where footprint is a float array with values between 0.0 (no coverage) and 1.0 (full coverage). When oversampling is active, fractional values indicate partial sub-pixel coverage. Default is False for backward compatibility.- verbosebool or callable
Logging control.
- Returns:
- coadd2D
~numpy.ndarrayor None Reprojected (and optionally coadded) image.
- footprint2D
~numpy.ndarray Coverage map (only returned when
return_footprint=True).
- coadd2D
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=Trueso that it will be handled by passingRESAMPLING_TYPE=FLAGSandCOMBINE_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
extraargument which should be the dictionary with parameter names as keys.