ZAP (the Zurich Atmosphere Purge) is a high precision sky subtraction tool which can be used as complete sky subtraction solution, or as an enhancement to previously sky-subtracted MUSE integral field spectroscopic data. The method uses PCA to isolate the residual sky subtraction features and remove them from the observed datacube. Though the operation of ZAP is not dependent on perfect flatfielding of the data in a MUSE exposure, better results are obtained when these corrections are made ahead of time. Future development will include expansion to more instruments.
ZAP requires the following packages:
Many linear algebra operations are performed in ZAP, so it can be beneficial to use an alternative BLAS package. In the Anaconda distribution, the default BLAS comes with Numpy linked to OpenBlas, which can amount to a 20% speedup of ZAP.
In its most hands-off form, ZAP can take an input fits datacube, operate on it, and output a final fits datacube:
import zap
zap.process('INPUT.fits', outcubefits='OUTPUT.fits')
Care should be taken, however, since this case assumes a sparse field, and better results can be obtained by applying masks.
The main function is zap.process:
There are a number of options that can be passed to the code which we describe here:
Performs the entire ZAP sky subtraction algorithm.
Work on an input FITS file and optionally writes the product to an output FITS file.
| Parameters: | musecubefits : str
outcubefits : str
clean : bool
zlevel : str
cftype : str
cfwidthSVD : int or float
cfwidthSP : int or float
optimizeType : str
pevals : list
nevals : list
extSVD : str
skycubefits : str
svdoutputfits : str
mask : str
interactive : bool
|
|---|
The code can handle datacubes trimmed in wavelength space. Since the code uses the correlation of segments of the emission line spectrum, it is best to trim the cube at specific wavelengths. The cube can include any connected subset of these segments. (for example 6400 - 8200 Angstroms)
[0, 5400]
[5400, 5850]
[5850, 6440]
[6440, 6750]
[6750, 7200]
[7200, 7700]
[7700, 8265]
[8265, 8602]
[8602, 8731]
[8731, 9275]
[9275, 10000]
This case specifically refers to the case where the sky can be measured in the sky frame itself, using:
zap.process('INPUT.fits', outcubefits='OUTPUT.fits')
In both cases, the code will create a resulting processed datacube named DATACUBE_ZAP.fits and an SVD file named ZAP_SVD.fits in the current directory. While this can work well in the case of very faint sources, masks can improve the results.
For the sparse field case, a mask file can be included, which is a 2d fits image matching the spatial dimensions of the input datacube. Masks are defined to be >= 1 on astronomical sources and 0 at the position of the sky. Set this parameter with the mask keyword
zap.process('INPUT.fits', outcubefits='OUTPUT.fits', mask='mask.fits')
This approach also can address the saturated field case and is robust in the case of strong emission lines, in this case the input is an offset sky observation. To achieve this, we calculate the SVD on an external sky frame using the function zap.SVDoutput
Performs the SVD decomposition of a datacube.
This allows to use the SVD for a different datacube.
| Parameters: | musecubefits : str
svdoutputfits : str
clean : bool
zlevel : str
cftype : str
cfwidth : int or float
mask : str
|
|---|
An example of running the code in this way is as follows:
zap.SVDoutput('Offset_Field_CUBE.fits', svdfn='ZAP_SVD.fits', mask='mask.fits')
zap.process('Source_cube.fits', outcubefits='OUTPUT.fits', extSVD='ZAP_SVD.fits', cfwidthSP=50)
The integration time of this frame does not need to be the same as the object exposure, but rather just a 2-3 minute exposure. Often residuals can be further reduced by changing cfwidthSP to a smaller value. However, this parameter should not be reduced to smaller than 15 pixels.
Aside from the main process, two functions are included that can be run outside of the entire zap process to facilitate some investigations.
Detects NaN values in cube and removes them by replacing them with an interpolation of the nearest neighbors in the data cube. The positions in the cube are retained in nancube for later remasking.
| Parameters: | musecubefits : str
outfn : str
rejectratio : float
boxsz : int
|
|---|
Performs a weighted median filtering of a 1d spectrum
Operates using a cumulative sum curve
| Parameters: | spec : numpy.ndarray
wt : numpy.ndarray
cfwidth : int or float
|
|---|
ZAP can also be used from the command line:
python -m zap INPUT_CUBE.fits
More information use of the command line interface can be found with the command
python -m zap -h
ZAP can also be used interactively from within IPython
import zap
zobj = zap.process('INPUT.fits', interactive=True)
The run method operates on the datacube, and retains all of the data and methods necessary to process a final data cube in a python class named zclass. You can elect to investigate the data product via the zclass, and even reprocess the cube with a different number of eigenspectra per region. A workflow may go as follows:
import zap
from matplotlib import pyplot as plt
# allow ZAP to run the optimize routine
zobj = zap.process('INPUT.fits', optimization='normal', interactive=True)
# plot the variance curves and the selection of the number of eigenspectra used
zobj.plotvarcurve(5)
# plot a spectrum extracted from the original cube
plt.figure()
plt.plot(zobj.cube[:,50:100,50:100].sum(axis=(1,2)), 'b', alpha=0.3)
# plot a spectrum of the cleaned ZAP dataproduct
plt.plot(zobj.cleancube[:,50:100,50:100].sum(axis=(1,2)), 'g')
# choose just the first 3 spectra for all segmments
zobj.reprocess(nevals=3)
# plot a spectrum extracted from the original cube
plt.plot(zobj.cube[:,50:100,50:100].sum(axis=(1,2)), 'b', alpha=0.3)
# plot a spectrum of the cleaned ZAP dataproduct
plt.plot(zobj.cleancube[:,50:100,50:100].sum(axis=(1,2))), 'g')
# choose some number of modes by hand
zobj.reprocess(nevals=[2,5,2,4,6,7,9,8,5,3,5])
# plot a spectrum
plt.plot(zobj.cleancube[:,50:100,50:100].sum(axis=(1,2))), 'k')
# Use the optimization algorithm to identify the best number of modes per segment
zobj.optimize()
# compare to the previous versions
plt.plot(zobj.cleancube[:,50:100,50:100].sum(axis=(1,2))), 'r')
# identify a pixel in the dispersion axis that shows a residual feature in the original
plt.figure()
plt.matshow(zobj.cube[2903,:,:])
# compare this to the zap dataproduct
plt.figure()
plt.matshow(zobj.cleancube[2903,:,:])
# write the processed cube as a single extension fits
zobj.writecube('DATACUBE_ZAP.fits')
# or merge the zap datacube into the original input datacube, replacing the data extension
zobj.writefits(outcubefits='DATACUBE_FINAL_ZAP.fits')
Main class to run each of the steps of ZAP.
Attributes
| cleancube | (numpy.ndarray) The final datacube after removing all of the residual features. |
| contarray | (numpy.ndarray) A 2D array containing the subtracted continuum per spaxel. |
| cube | (numpy.ndarray) The original cube with the zlevel subtraction performed per spaxel. |
| especeval | (list of (eigenspectra, eval)) A list containing the full set of eigenspectra and eigenvalues generated by the SVD calculation that is used toy reconstruct the entire datacube. |
| laxis | (numpy.ndarray) A 1d array containing the wavelength solution generated from the header parameters. |
| wcs | (astropy.wcs.WCS) WCS object with the wavelength solution. |
| lranges | (list) A list of the wavelength bin limits used in segmenting the sepctrum for SVD. |
| nancube | (numpy.ndarray) A 3d boolean datacube containing True in voxels where a NaN value was replaced with an interpolation. |
| nevals | (numpy.ndarray) A 1d array containing the number of eigenvalues used per segment to reconstruct the residuals. |
| normstack | (numpy.ndarray) A normalized version of the datacube decunstructed into a 2d array. |
| varlist | (numpy.ndarray) An array for each segment with the variance curve, calculated for the optimize method. |
| pranges | (numpy.ndarray) The pixel indices of the bounding regions for each spectral segment. |
| recon | (numpy.ndarray) A 2d array containing the reconstructed emission line residuals. |
| run_clean | (bool) Boolean that indicates that the NaN cleaning method was used. |
| run_zlevel | (bool) Boolean indicating that the zero level correction was used. |
| stack | (numpy.ndarray) The datacube deconstructed into a 2d array for use in the the SVD. |
| subespeceval | (list of (eigenspectra, eval)) The subset of eigenvalues and eigenspectra used to reconstruct the sky residuals. |
| variancearray | (numpy.ndarray) A list of length nsegments containing variances calculated per spaxel used for normalization |
| y,x | (numpy.ndarray) The position in the cube of the spaxels that are in the 2d deconstructed stack |
| zlsky | (numpy.ndarray) A 1d array containing the result of the zero level subtraction |
Methods
| chooseevals | |
| make_contcube | |
| mergefits | |
| optimize | |
| plotvarcurve | |
| reconstruct | |
| remold | |
| reprocess | |
| writeSVD | |
| writecube | |
| writeskycube |
Choose the number of eigenspectra/evals to use for reconstruction.
User supplies the number of eigen spectra to be used (neval) or the percentage of the eigenspectra that were calculated (peval) from each spectral segment to be used.
The user can either provide a single value to be used for all segments, or provide an array that defines neval or peval per segment.
Remold the continuum array so it can be investigated.
Takes the continuum stack and returns it into a familiar cube form.
Merge the ZAP cube into the full muse datacube and write.
Function to optimize the number of components used to characterize the residuals.
This function calculates the variance per segment with an increasing number of eigenspectra/eigenvalues. It then deterimines the point at which the second derivative of this variance curve reaches zero. When this occurs, the linear reduction in variance is attributable to the removal of astronomical features rather than emission line residuals.
Reconstruct the residuals from a given set of eigenspectra and eigenvalues
Subtracts the reconstructed residuals and places the cleaned spectra into the duplicated datacube.
A method that redoes the eigenvalue selection, reconstruction, and remolding of the data.
Write the SVD to an individual fits file.
Write the processed datacube to an individual fits file.
Write the processed datacube to an individual fits file.