Coadd 3D Spectra
Overview
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.
pypeit_coadd_datacube
The primary script is called pypeit_coadd_datacube
, which takes
an input file to guide the process.
usage
The script usage can be displayed by calling the script with the
-h
option:
$ 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)
coadd3d file
The 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
saved as BB1245p4238.coadd3d
:
# 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, pypeit_coadd_datacube
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
run
Then run the script:
pypeit_coadd_datacube BB1245p4238.coadd3d -o
Combination options
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 spec_subpixel
and spat_subpixel
keywords
as follows:
[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
with the NGP
method. This algorithm is effectively a 3D histogram. This approach is faster
than 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.
Flux calibration
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
Sky Subtraction
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 (all lowercase):
[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.
Grating correction
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
coadd3d
file:
[reduce]
[[cube]]
grating_corr = False
Astrometric correction
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 coadd3d
file:
[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
keyword argument:
[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 _whitelight.fits
.
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 coadd3d
file:
[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.coadd3d 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()