Coadd 3D Spectra
This document describes how to combine a set of fully reduced 2D spectra from multiple exposures into a single 3D datacube for IFU spectrographs.
This must be done outside of the data reduction pipeline (run_pypeit); i.e., PypeIt will not coadd your spectra as part of the data reduction process.
The primary script is called
pypeit_coadd_datacube, which takes
an input file to guide the process.
The script usage can be displayed by calling the script with the
$ pypeit_coadd_datacube -h usage: pypeit_coadd_datacube [-h] [--det DET] [-o] [-v VERBOSITY] file Read in an array of spec2D files and convert them into a datacube positional arguments: file filename.coadd3d file options: -h, --help show this help message and exit --det DET Detector (default: 1) -o, --overwrite Overwrite any existing files/directories (default: False) -v VERBOSITY, --verbosity VERBOSITY Verbosity level between 0 [none] and 2 [all]. Default: 1. Level 2 writes a log with filename coadd_datacube_YYYYMMDD-HHMM.log (default: 1)
pypeit_coadd_datacube script requires an
input file to guide the process.
The format of this type of Input File Format
includes a Parameter Block (required)
and a Data Block (required).
In the following example for
keck_kcwi, the coadd3d file will be
# User-defined execution parameters [rdx] spectrograph = keck_kcwi detnum = 1 [reduce] [[cube]] combine = True output_filename = BB1245p4238_datacube.fits save_whitelight = True # Read in the data spec2d read filename | scale_corr Science/spec2d_scienceframe_01.fits | Science/spec2d_scalecorr.fits Science/spec2d_scienceframe_02.fits | Science/spec2d_scalecorr.fits spec2d end
The opening block sets parameters for the reduction steps. Note, by default,
will convert all spec2d files into a spec3d file (i.e. individual datacubes for each exposure).
If you want to combine all exposures into a single datacube, you need to set
combine = True,
as in the above example, and provide an
output_filename. This is very useful if you want to
combine several standard star exposures into a single datacube for flux calibration, for example.
The spec2d block provides a list of Spec2D Output files. You can also specify an optional scale correction
as part of the spec2d block. This relative scale correction ensures that the relative spectral sensitivity of the
datacube is constant across the field of view. The spec2d file used for the
scale_corr column should either be a
twilight or dome flat reduced as a
science frame (see KECK KCWI for a description of what you need to do).
In order to use this functionality, you should not reduce your science data with a spectral illumination correction.
In other words, in your PypeIt Reduction File file, set the following when you execute run_pypeit:
[scienceframe] [[process]] use_specillum = False
Then run the script:
pypeit_coadd_datacube BB1245p4238.coadd3d -o
PypeIt currently supports two different methods to convert an spec2d frame into a datacube;
these options are called
subpixel (default) and
NGP (which is short for, nearest grid point),
and can be set using the following keyword arguments:
[reduce] [[cube]] method = ngp
The default option is called
subpixel, which divides each pixel in the spec2d frame
into many subpixels, and assigns each subpixel to a voxel of the datacube. Flux is conserved,
but voxels are correlated, and the error spectrum does not account for covariance between
adjacent voxels. The subpixellation scale can be separately set in the spatial and spectral
direction on the 2D detector. If you would like to change the subpixellation factors from
the default values (5), you can set the
[reduce] [[cube]] method = subpixel spec_subpixel = 8 spat_subpixel = 10
The total number of subpixels generated for each detector pixel on the spec2d frame is
spec_subpixel x spat_subpixel. The default values (5) divide each spec2d pixel into 25 subpixels
during datacube creation. As an alternative, you can convert the spec2d frames into a datacube
NGP method. This algorithm is effectively a 3D histogram. This approach is faster
subpixel, flux is conserved, and voxels are not correlated. However, this option suffers
the same downsides as any histogram; the choice of bin sizes can change how the datacube appears.
This algorithm takes each pixel on the spec2d frame and puts the flux of this pixel into one voxel
in the datacube. Depending on the binning used, some voxels may be empty (zero flux) while a
neighbouring voxel might contain the flux from two spec2d pixels.
If you would like to flux calibrate your datacube, you need to
produce your standard star datacube first, and when generating
the datacube of the science frame you must pass in the name of
the standard star cube in your
coadd3d file as follows:
[reduce] [[cube]] standard_cube = standard_star_cube.fits
The default behaviour of PypeIt is to subtract the model sky that is derived from the science frame during the reduction. If you would like to turn off sky subtraction, set the following keyword argument:
[reduce] [[cube]] skysub_frame = None
If you would like to use a dedicated sky frame for sky subtraction
that is separate from the science frame, then you need to provide
the relative path+file of the spec2d file that you would like to
use. If you need a different sky frame for different science frames,
then you can specify the
skysub_frame in the
spec2d block of the
.coadd3d file, similar to the way
scale_corr is set in the example
above. If you have dedicated sky frames, then it is generally
recommended to reduce these frames as if they are regular science
frames, but add the following keyword arguments at the top of your
PypeIt Reduction File:
[reduce] [[skysub]] joint_fit = True user_regions = : [flexure] spec_method = slitcen
This ensures that all pixels in the slit are used to generate a complete model of the sky.
The grating correction is needed if any of the data are recorded
with even a very slightly different setup (e.g. data taken on two
different nights with the same intended wavelength coverage,
but the grating angle of the two nights were slightly different).
This is also needed if your standard star observations were taken
with a slightly different setup. This correction requires that you
have taken calibrations (i.e. flatfields) with the two different
setups. By default, the grating correction will be applied, but it
can be disabled by setting the following keyword argument in your
[reduce] [[cube]] grating_corr = False
If you would like to perform an astrometric correction, you
need to install scikit-image (version > 0.17;
see Install via pip or simply install scikit-image with pip directly). The default
option is to perform the astrometric correction, if a Alignment
frame has been computed. To disable the astrometric
correction, set the following keyword argument in your
[reduce] [[cube]] astrometric = False
White light image
A white light image can be generated for the combined frame, or
for each individual frame if
combine=False, by setting the following
[reduce] [[cube]] save_whitelight = True
White light images are not produced by default. The output filename for
the white light images are given the suffix
Spatial alignment with different setups
If you have multiple setups that you want to align so that all
pixels are spatially coincident, you must first produce the
datacube that you wish to use as a reference. Then, define the
WCS parameters using the keyword arguments in your
[reduce] [[cube]] reference_image = reference_cube_whitelight.fits ra_min = 191.398441 ra_max = 191.401419 dec_min = 42.634352 dec_max = 42.639988 spatial_delta = 0.339462
where these values are printed as terminal output after
reference_cube.fits is generated.
Note that PypeIt is not currently setup to stitch together cubes covering different wavelength range, but it can coadd multiple spec2D files into a single datacube if the wavelength setup overlaps, and the spatial positions are very similar.
Combining multiple datacubes
PypeIt is able to combine standard star frames for flux calibration, and should not have any difficulty with this. If your science observations are designed so that there is very little overlap between exposures, you should not assume that the automatic combination algorithm will perform well. Instead, you may prefer to output individual data cubes and manually combine the cubes with some other purpose-built software. If you know the relative offsets very well, then you can specify these, and PypeIt can combine all frames into a single combined datacube. This is the recommended approach, provided that you know the relative offsets of each frame. In the following example, the first cube is assumed to be the reference cube (0.0 offset in both RA and Dec), and the second science frame is offset relative to the first by:
Delta RA x cos(Dec) = 1.0" W Delta Dec = 2.0" N
The offset convention used in PypeIt is that positive offsets translate the RA and Dec of a frame to higher RA (i.e. more East) and higher Dec (i.e. more North). In the above example, frame 2 is 1” to the West of frame 1, meaning that we need to move frame 2 by 1” to the East (i.e. a correction of +1”). Similarly, we need to more frame 2 by 2” South (i.e. a correction of -2”). Therefore, in the above example, the coadd3d file would look like the following:
# User-defined execution parameters [rdx] spectrograph = keck_kcwi detnum = 1 [reduce] [[cube]] combine = True output_filename = BB1245p4238_datacube.fits align = True # Read in the data spec2d read filename | ra_offset | dec_offset Science/spec2d_scienceframe_01.fits | 0.0 | 0.0 Science/spec2d_scienceframe_02.fits | 1.0 | -2.0 spec2d end
Current Coadd3D Data Model
The output is a single fits file that contains a datacube, and
a cube with the same shape that stores the variance. The units
are stored in the
FLUXUNIT header keyword.
Here is a short python script that will allow you to read in and plot a wavelength slice of the cube:
from matplotlib import pyplot as plt from astropy.visualization import ZScaleInterval, ImageNormalize from pypeit.core.datacube import DataCube filename = "datacube.fits" cube = DataCube.from_file(filename) flux_cube = cube.flux # Flux datacube error_cube = cube.sig # Errors associated with each voxel of the flux datacube ivar_cube = cube.ivar # Inverse variance cube wcs = cube.wcs wave_slice = 1000 norm = ImageNormalize(flux_cube[wave_slice,:,:], interval=ZScaleInterval()) fig = plt.figure() fig.add_subplot(111, projection=wcs, slices=('x', 'y', wave_slice)) plt.imshow(flux_cube[wave_slice,:,:], origin='lower', cmap=plt.cm.viridis, norm=norm) plt.xlabel('RA') plt.ylabel('Dec') plt.show()