specwcs: Discussion and description of the spectral image formats
Package: onedspec
1. types of spectral data
Spectra are stored as one, two, or three dimensional images with one axis being the dispersion axis. A pixel value is the flux over some interval of wavelength and position. The simplest example of a spectrum is a one dimensional image which has pixel values as a function of wavelength.
There are two types of higher dimensional spectral image formats. One type has spatial axes for the other dimensions and the dispersion axis may be along any of the image axes. Typically this type of format is used for long slit (two dimensional) and Fabry-Perot (three dimensional) spectra. This type of spectra is referred to as spatial spectra and the world coordinate system (WCS) format is called ndspec. The details of the world coordinate systems are discussed later.
The second type of higher dimensional spectral image consists of multiple, independent, one dimensional spectra stored in the higher dimensions with the first image axis being the dispersion axis; i.e. each line is a spectrum. This format allows associating many spectra and related parameters in a single data object. This type of spectra is referred to as multispec and the there are two coordinate system formats, equispec and multispec. The equispec format applies to the common case where all spectra have the same linear dispersion relation. The multispec format applies to the general case of spectra with differing dispersion relations or non-linear dispersion functions. These multi-spectrum formats are important since maintaining large numbers of spectra as individual one dimensional images is very unwieldy for the user and inefficient for the software.
Examples of multispec spectral images are spectra extracted from a multi-fiber or multi-aperture spectrograph or orders from an echelle spectrum. The second axis is some arbitrary indexing of the spectra, called apertures, and the third dimension is used for associated quantities. The IRAF apextract package may produce multiple spectra from a CCD image in successive image lines with an optimally weighted spectrum, a simple aperture sum spectrum, a background spectrum, and sigma spectrum as the associated quantities along the third dimension of the image.
Many onedspec package tasks which are designed to operate on individual one dimensional spectra may operate on spatial spectra by summing a number of neighboring spectra across the dispersion axis. This eliminates the need to "extract" one dimensional spectra from the natural format of this type of data in order to use tasks oriented towards the display and analysis of one dimensional spectra. The dispersion axis is either given in the image header by the keyword DISPAXIS or the package dispaxis parameter. The summing factors across the dispersion are specified by the nsum package parameter. See "help onedspec.package" for information on these parmaeters.
One dimensional spectra, whether from multispec spatial spectra, have several associated quantities which may appear in the image header as part of the coordinate system description. The primary identification of a spectrum is an integer aperture number. This number must be unique within a single image. There is also an integer beam number used for various purposes such as discriminating object, sky, and arc spectra in multi-fiber/multi-aperture data or to identifying the order number in echelle data. For spectra summed from spatial spectra the aperture number is the central line, column, or band. In 3D images the aperture index wraps around the lowest non-dispersion axis. Since most one dimensional spectra are derived from an integration over one or more spatial axes, two additional aperture parameters record the aperture limits. These limits refer to the original pixel limits along the spatial axis. This information is primarily for record keeping but in some cases it is used for spatial interpolation during dispersion calibration. These values are set either by the apextract tasks or when summing neighboring vectors in spatial spectra.
An important task to be aware of for manipulating spectra between image formats is scopy. This task allows selecting spectra from multispec images and grouping them in various ways and also "extracts" apertures from long slit and 3D spectra simply and without resort to the more general apextract package.
2. world coordinate systems
IRAF images have three types of coordinate systems. The pixel array coordinates of an image or image section, i.e. the lines and columns, are called the logical coordinates. The logical coordinates of individual pixels change as sections of the image are used or extracted. Pixel coordinates are tied to the data, i.e. are fixed to features in the image, are called physical coordinates. Initially the logical and physical coordinates are the equivalent but differ when image sections or other tasks which modify the sampling of the pixels are applied.
The last type of coordinate system is called the world coordinate system. Like the physical coordinates, the world coordinates are tied to the features in the image and remain unchanged when sections of the image are used or extracted. If a world coordinate system is not defined for an image, the physical coordinate system is considered to be the world coordinate system. In spectral images the world coordinate system includes dispersion coordinates such as wavelengths. In many tasks outside the spectroscopy packages, for example the plot, tv and images packages, one may select the type of coordinate system to be used. To make plots and get coordinates in dispersion units for spectra with these tasks one selects the "world" system. The spectral tasks always use world coordinates.
The coordinate systems are defined in the image headers using a set of reserved keywords which are set, changed, and updated by various tasks. Some of the keywords consist of simple single values following the FITS convention. Others, the WAT keywords, encode long strings of information, one for each coordinate axis and one applying to all axes, into a set of sequential keywords. The values of these keywords must then be pasted together to recover the string. The long strings contain multiple pieces called WCS attributes. In general the WCS keywords should be left to IRAF tasks to modify. However, if one wants modify them directly some tasks which may be used are hedit, hfix, wcsedit, wcsreset, specshift, dopcor, and sapertures. The first two are useful for the simple keywords, the two "wcs" tasks are useful for the linear ndspec and equispec formats, the next two are for the common cases of shifting the coordinate zero point or applying a doppler correction, and the last one is the one to use for the more complex multispec format attributes.
3. physical coordinate system
The physical coordinate system is used by the spectral tasks when there is no dispersion coordinate information (such as before dispersion calibration), to map the physical dispersion axis to the logical dispersion axis, and in the multispec world coordinate system dispersion functions which are defined in terms of physical coordinates.
The transformation between logical and physical coordinates is defined by the header keywords LTVi, LTMi_j (where i and j are axis numbers) through the vector equation
l = |m| * p + v
where l is a logical coordinate vector, p is a physical coordinate vector, v is the origin translation vector specified by the LTV keywords and |m| is the scale/rotation matrix specified by the LTM keywords. For spectra rotation terms (nondiagonal matrix elements) generally do not make sense (in fact many tasks will not work if there is a rotation) so the transformations along each axis are given by the linear equation
where l is a logical coordinate vector, p is a physical coordinate vector, v is the origin translation vector specified by the LTV keywords and |m| is the scale/rotation matrix specified by the LTM keywords. For spectra a rotation term (nondiagonal matrix elements) generally does not make sense (in fact many tasks will not work if there is a rotation) so the transformations along each axis are given by the linear equation
li = LTMi_i * pi + LTVi.
If all the LTM/LTV keywords are missing they are assumed to have zero values except that the diagonal matrix terms, LTMi_i, are assumed to be 1. Note that if some of the keywords are present then a missing LTMi_i will take the value zero which generally causes an arithmetic or matrix inversion error in the IRAF tasks.
The dimensional mapping between logical and physical axes is given by the keywords WCSDIM and WAXMAP01. The WCSDIM keyword gives the dimensionality of the physical and world coordinate system. There must be coordinate information for that many axes in the header (though some may be missing and take their default values). If the WCSDIM keyword is missing it is assumed to be the same as the logical image dimensionality.
The syntax of the WAXMAP keyword are pairs of integer values, one for each physical axis. The first number of each pair indicates which current logical axis corresponds to the original physical axis (in order) or zero if that axis is missing. When the first number is zero the second number gives the offset to the element of the original axis which is missing. As an example consider a three dimensional image in which the second plane is extracted (an IRAF image section of [*,2,*]). The keyword would then appear as WAXMAP01 = '1 0 0 1 2 0'. If this keyword is missing the mapping is 1:1; i.e. the dimensionality and order of the axes are the same.
The dimensional mapping is important because the dispersion axis for the nspec spatial spectra as specified by the DISPAXIS keyword or task parameter, or the axis definitions for the equispec and or multispec formats are always in terms of the original physical axes.
4. linear spectral world coordinate systems
When there is a linear or logarithmic relation between pixels and dispersion coordinates which is the same for all spectra the WCS header format is simple and uses the FITS convention (with the CD matrix keywords proposed by Hanisch and Wells 1992) for the logical pixel to world coordinate transformation. This format applies to one, two, and three dimensional data. The higher dimensional data may have either linear spatial axes or the equispec format where each one dimensional spectrum stored along the lines of the image has the same dispersion.
The FITS image header keywords describing the spectral world coordinates are CTYPEi, CRPIXi, CRVALi, and CDi_j where i and j are axis numbers. As with the physical coordinate transformation the nondiagonal or rotation terms are not expected in the spectral WCS and may cause problems if they are not zero. The CTYPEi keywords will have the value LINEAR to identify the type of coordinate system. The transformation between dispersion coordinate, wi, and logical pixel coordinate, li, along axis i is given by
wi = CRVALi + CDi_i * (li - CRPIXi)
If the keywords are missing then the values are assumed to be zero except for the diagonal elements of the scale/rotation matrix, the CDi_i, which are assumed to be 1. If only some of the keywords are present then any missing CDi_i keywords will take the value 0 which will cause IRAF tasks to fail with arithmetic or matrix inversion errors. If the CTYPEi keyword is missing it is assumed to be "LINEAR".
If the pixel sampling is logarithmic in the dispersion coordinate, as required for radial velocity cross-correlations, the WCS coordinate values are logarithmic and wi (above) is the logarithm of the dispersion coordinate. The spectral tasks (though not other tasks) will recognize this case and automatically apply the anti-log. The two types of pixel sampling are identified by the value of the keyword DC-FLAG. A value of 0 defines a linear sampling of the dispersion and a value of 1 defines a logarithmic sampling of the dispersion. Thus, in all cases the spectral tasks will display and analyze the spectra in the same dispersion units regardless of the pixel sampling.
Other keywords which may be present are DISPAXIS for 2 and 3 dimensional spatial spectra, and the WCS attributes "system", "wtype", "label", and "units". The system attribute will usually have the value "world" for spatial spectra and "equispec" for equispec spectra. The wtype attribute will have the value "linear". Currently the label will be either "Pixel" or "Wavelength" and the units will be "Angstroms" for dispersion corrected spectra. In the future there will be more generality in the units for dispersion calibrated spectra.
Figure 1 shows the WCS keywords for a two dimensional long slit spectrum. The coordinate system is defined to be a generic "world" system and the wtype attributes and CTYPE keywords define the axes to be linear. The other attributes define a label and unit for the second axis, which is the dispersion axis as indicated by the DISPAXIS keyword. The LTM/LTV keywords in this example show that a subsection of the original image has been extracted with a factor of 2 block averaging along the dispersion axis. The dispersion coordinates are given in terms of the logical pixel coordinates by the FITS keywords as defined previously.
Figure 1: Long Slit Spectrum
WAT0_001= 'system=world'
WAT1_001= 'wtype=linear'
WAT2_001= 'wtype=linear label=Wavelength units=Angstroms'
WCSDIM = 2
DISPAXIS= 2
DC-FLAG = 0
CTYPE1 = 'LINEAR '
LTV1 = -10.
LTM1_1 = 1.
CRPIX1 = -9.
CRVAL1 = 19.5743865966797
CD1_1 = 1.01503419876099
CTYPE2 = 'LINEAR '
LTV2 = -49.5
LTM2_2 = 0.5
CRPIX2 = -49.
CRVAL2 = 4204.462890625
CD2_2 = 12.3337936401367
Figure 2 shows the WCS keywords for a three dimensional image where each line is an independent spectrum or associated data but where all spectra have the same linear dispersion. This type of coordinate system has the system name "equispec". The ancillary information about each aperture is found in the APNUM keywords. These give the aperture number, beam number, and extraction limits. In this example the LTM/LTV keywords have their default values; i.e. the logical and physical coordinates are the same.
Figure 2: Equispec Spectrum
WAT0_001= 'system=equispec'
WAT1_001= 'wtype=linear label=Wavelength units=Angstroms'
WAT2_001= 'wtype=linear'
WAT3_001= 'wtype=linear'
WCSDIM = 3
DC-FLAG = 0
APNUM1 = '41 3 7.37 13.48'
APNUM2 = '15 1 28.04 34.15'
APNUM3 = '33 2 43.20 49.32'
CTYPE1 = 'LINEAR '
LTM1_1 = 1.
CRPIX1 = 1.
CRVAL1 = 4204.463
CD1_1 = 6.16689700000001
CTYPE2 = 'LINEAR '
LTM2_2 = 1.
CD2_2 = 1.
CTYPE3 = 'LINEAR '
LTM3_3 = 1.
CD3_3 = 1.
5. multispec spectral world coordinate system
The multispec spectral world coordinate system applies only to one dimensional spectra; i.e. there is no analog for the spatial type spectra. It is used either when there are multiple 1D spectra with differing dispersion functions in a single image or when the dispersion functions are nonlinear.
The multispec coordinate system is always two dimensional though there may be an independent third axis. The two axes are coupled and they both have axis type "multispec". When the image is one dimensional the physical line is given by the dimensional reduction keyword WAXMAP. The second, line axis, has world coordinates of aperture number. The aperture numbers are integer values and need not be in any particular order but do need to be unique. This aspect of the WCS is not of particular user interest but applications use the inverse world to physical transformation to select a spectrum line given a specified aperture.
The dispersion functions are specified by attribute strings with the identifier specN where N is the physical image line. The attribute strings contain a series of numeric fields. The fields are indicated symbolically as follows.
specN = ap beam dtype w1 dw nw z aplow aphigh [functions_i]
where there are zero or more functions having the following fields,
function_i = wt_i w0_i ftype_i [parameters] [coefficients]
The first nine fields in the attribute are common to all the dispersion functions. The first field of the WCS attribute is the aperture number, the second field is the beam number, and the third field is the dispersion type with the same function as DC-FLAG in the nspec and equispec formats. A value of -1 indicates the coordinates are not dispersion coordinates (the spectrum is not dispersion calibrated), a value of 0 indicates linear dispersion sampling, a value of 1 indicates log-linear dispersion sampling, and a value of 2 indicates a nonlinear dispersion.
The next two fields are the dispersion coordinate of the first physical pixel and the average dispersion interval per physical pixel. For linear and log-linear dispersion types the dispersion parameters are exact while for the nonlinear dispersion functions they are approximate. The next field is the number of valid pixels, hence it is possible to have spectra with varying lengths in the same image. In that case the image is as big as the biggest spectrum and the number of pixels selects the actual data in each image line. The next (seventh) field is a doppler factor. This doppler factor is applied to all dispersion coordinates by multiplying by 1/(1+z) (assuming wavelength dispersion units). Thus a value of 0 is no doppler correction. The last two fields are extraction aperture limits as discussed previously.
Following these fields are zero or more function descriptions. For linear or log-linear dispersion coordinate systems there are no function fields. For the nonlinear dispersion systems the function fields specify a weight, a zero point offset, the type of dispersion function, and the parameters and coefficients describing it. The function type codes, ftype_i, are 1 for a chebyshev polynomial, 2 for a legendre polynomial, 3 for a cubic spline, 4 for a linear spline, 5 for a pixel coordinate array, and 6 for a sampled coordinate array. The number of fields before the next function and the number of functions are determined from the parameters of the preceding function until the end of the attribute is reached.
The equation below shows how the final wavelength is computed based on the nfunc individual dispersion functions W_i(p). Note that this is completely general in that different function types may be combined. However, in practice when multiple functions are used they are generally of the same type and represent a calibration before and after the actual object observation with the weights based on the relative time difference between the calibration dispersion functions and the object observation.
w = sum from i=1 to nfunc {wt_i * (w0_i + W_i(p)) / (1 + z)}
The multispec coordinate systems define a transformation between physical pixel, p, and world coordinates, w. Generally there is an intermediate coordinate system used. The following equations define these coordinates. The first one shows the transformation between logical, l, and physical, p, coordinates based on the LTM/LTV keywords. The polynomial functions are defined in terms of a normalized coordinate, n, as shown in the second equation. The normalized coordinates run between -1 and 1 over the range of physical coordinates, pmin and pmax which are parameters of the function, upon which the coefficients were defined. The spline functions map the physical range into an index over the number of evenly divided spline pieces, npieces, which is a parameter of the function. This mapping is shown in the third and fourth equations where s is the continuous spline coordinate and j is the nearest integer less than or equal to s.
p = (l - LTV1) / LTM1_1
n = (p - pmiddle) / (prange / 2)
= (p - (pmax+pmin)/2) / ((pmax-pmin) / 2)
s = (p - pmin) / (pmax - pmin) * npieces
j = int(s)
5.1 linear and log linear dispersion function
The linear and log-linear dispersion functions are described by a wavelength at the first physical pixel and a wavelength increment per physical pixel. A doppler correction may also be applied. The equations below show the two forms. Note that the coordinates returned are always wavelength even though the pixel sampling and the dispersion parameters may be log-linear.
w = (w1 + dw * (p - 1)) / (1 + z)
w = 10 ** {(w1 + dw * (p - 1)) / (1 + z)}
Figure 3 shows an example from a multispec image with independent linear dispersion coordinates. This is a linearized echelle spectrum where each order (identified by the beam number) is stored as a separate image line.
Figure 3: Echelle Spectrum with Linear Dispersion Function
WAT0_001= 'system=multispec'
WAT1_001= 'wtype=multispec label=Wavelength units=Angstroms'
WAT2_001= 'wtype=multispec spec1 = "1 113 0 4955.44287109375 0.05...
WAT2_002= '5 256 0. 23.22 31.27" spec2 = "2 112 0 4999.0810546875...
WAT2_003= '58854293 256 0. 46.09 58.44" spec3 = "3 111 0 5043.505...
WAT2_004= '928358078002 256 0. 69.28 77.89"
WCSDIM = 2
CTYPE1 = 'MULTISPE'
LTM1_1 = 1.
CD1_1 = 1.
CTYPE2 = 'MULTISPE'
LTM2_2 = 1.
CD2_2 = 1.
5.2 chebyshev polynomial dispersion function
The parameters for the chebyshev polynomial dispersion function are the order (number of coefficients) and the normalizing range of physical coordinates, pmin and pmax, over which the function is defined and which are used to compute n. Following the parameters are the order coefficients, ci. The equation below shows how to evaluate the function using an iterative definition where x_1 = 1, x_2 = n, and x_i = 2 * n * x_{i-1} - x_{i-2}.
The parameters for the chebyshev polynomial dispersion function are the order (number of coefficients) and the normalizing range of physical coordinates, pmin and pmax, over which the function is defined and which are used to compute n. Following the parameters are the order coefficients, c_i. The equation below shows how to evaluate the function using an iterative definition where x_1 = 1, x_2 = n, and x_i = 2 * n * x_{i-1} - x_{i-2}.
W = sum from i=1 to order {c_i * x_i}
5.3 legendre polynomial dispersion function
The parameters for the legendre polynomial dispersion function are the order (number of coefficients) and the normalizing range of physical coordinates, pmin and pmax, over which the function is defined and which are used to compute n. Following the parameters are the order coefficients, c_i. The equation below shows how to evaluate the function using an iterative definition where x_1 = 1, x_2 = n, and x_i = ((2i-3)*n*x_{i-1}-(i-2)*x_{i-2})/(i-1).
W = sum from i=1 to order {c_i * x_i}
Figure 4 shows an example from a multispec image with independent nonlinear dispersion coordinates. This is again from an echelle spectrum. Note that the IRAF echelle package determines a two dimensional dispersion function, in this case a bidimensional legendre polynomial, with the independent variables being the order number and the extracted pixel coordinate. To assign and store this function in the image is simply a matter of collapsing the two dimensional dispersion function by fixing the order number and combining all the terms with the same order.
Figure 4: Echelle Spectrum with Legendre Polynomial Function
WAT0_001= 'system=multispec'
WAT1_001= 'wtype=multispec label=Wavelength units=Angstroms'
WAT2_001= 'wtype=multispec spec1 = "1 113 2 4955.442888635351 0.05...
WAT2_002= '83 256 0. 23.22 31.27 1. 0. 2 4 1. 256. 4963.0163112090...
WAT2_003= '976664 -0.3191636898579552 -0.8169352858733255" spec2 =...
WAT2_004= '9.081188912082 0.06387049476832223 256 0. 46.09 58.44 1...
WAT2_005= '56. 5007.401409453303 8.555959076467951 -0.176732458267...
WAT2_006= '09935064388" spec3 = "3 111 2 5043.505764869474 0.07097...
WAT2_007= '256 0. 69.28 77.89 1. 0. 2 4 1. 256. 5052.586239197408 ...
WAT2_008= '271 -0.03173489817897474 -7.190562320405975E-4"
WCSDIM = 2
CTYPE1 = 'MULTISPE'
LTM1_1 = 1.
CD1_1 = 1.
CTYPE2 = 'MULTISPE'
LTM2_2 = 1.
CD2_2 = 1.
5.4 linear spline dispersion function
The parameters for the linear spline dispersion function are the number of spline pieces, npieces, and the range of physical coordinates, pmin and pmax, over which the function is defined and which are used to compute the spline coordinate s. Following the parameters are the npieces+1 coefficients, c_i. The two coefficients used in a linear combination are selected based on the spline coordinate, where a and b are the fractions of the interval in the spline piece between the spline knots, a=(j+1)-s, b=s-j, and x_0=a, and x_1=b.
W = sum from i=0 to 1 {c_(i+j) * x_i}
5.5 cubic spline dispersion function
The parameters for the cubic spline dispersion function are the number of spline pieces, npieces, and the range of physical coordinates, pmin and pmax, over which the function is defined and which are used to compute the spline coordinate s. Following the parameters are the npieces+3 coefficients, c_i. The four coefficients used are selected based on the spline coordinate. The fractions of the interval between the integer spline knots are given by a and b, a=(j+1)-s, b=s-j, and x_0 =a sup 3, x_1 =(1+3*a*(1+a*b)), x_2 =(1+3*b*(1+a*b)), and x_3 =b**3.
The parameters for the cubic spline dispersion function are the number of spline pieces, npieces, and the range of physical coordinates, pmin and pmax, over which the function is defined and which are used to compute the spline coordinate s. Following the parameters are the npieces+3 coefficients, c_i. The four coefficients used are selected based on the spline coordinate. The fractions of the interval between the integer spline knots are given by a and b, a=(j+1)-s, b=s-j, and x_0=a**3, x_1=(1+3*a*(1+a*b)), x_2=(1+3*b*(1+a*b)), and x_3=b**3.
W = sum from i=0 to 3 {c_(i+j) * x_i}
5.6 pixel array dispersion function
The parameters for the pixel array dispersion function consists of just the number of coordinates ncoords. Following this are the wavelengths at integer physical pixel coordinates starting with 1. To evaluate a wavelength at some physical coordinate, not necessarily an integer, a linear interpolation is used between the nearest integer physical coordinates and the desired physical coordinate where a and b are the usual fractional intervals k=int(p), a=(k+1)-p, b=p-k, and x_0=a, and x_1=b.
W = sum from i=0 to 1 {c_(i+j) * x_i}
5.7 sampled array dispersion function
The parameters for the sampled array dispersion function consists of the number of coordinate pairs, ncoords, and a dummy field. Following these are the physical coordinate and wavelength pairs which are in increasing order. The nearest physical coordinates to the desired physical coordinate are located and a linear interpolation is computed between the two sample points.