Lookup-table WCS for irregular spectral sampling doesn't look up #5462
Description
Bearing in mind this may be a) not implemented or b) user error, I seem to be having difficulties expressing the -TAB coordinate type. I have a (4096, 4096, 19) (in FITS axes order) datacube with 19 planes representing simulated observations through different filters.
I want to assign a bogus WCS so that I can visualize angular dimensions on each plane in ds9 or similar. It occurred to me that the third dimension actually is wavelength, and I think I can represent it in FITS WCS using a lookup table. However, when I try to do mywcs.all_pix2world I get unexpected behavior that might even be uninitialized memory (i.e. calling the same function twice gives two different results for the wavelength axis).
I call mywcs.all_pix2world([[2048, 2048, i],], 0) where i is in [0, 19). The first two dimensions come out as I expect (that is, 0, 0) but the last is usually 0 and sometimes a very small float value (like 1e-322):
>>> mywcs.all_pix2world([[2048, 2048, 2],], 0)
array([[ 5.55555555e-005, 5.55555556e-005, 1.58101007e-322]])Here's my WCS header:
WCSAXES = 3 / Number of coordinate axes
CRPIX1 = 2048.0 / Pixel coordinate of reference point
CRPIX2 = 2048.0 / Pixel coordinate of reference point
CRPIX3 = 0.0 / Pixel coordinate of reference point
CDELT1 = 5.5555555555556E-05 / [deg] Coordinate increment at reference point
CDELT2 = 5.5555555555556E-05 / [deg] Coordinate increment at reference point
CDELT3 = 1.0 / [m] Coordinate increment at reference point
CUNIT1 = 'deg' / Units of coordinate increment and value
CUNIT2 = 'deg' / Units of coordinate increment and value
CUNIT3 = 'm' / Units of coordinate increment and value
CTYPE1 = 'RA---TAN' / Right ascension, gnomonic projection
CTYPE2 = 'DEC--TAN' / Declination, gnomonic projection
CTYPE3 = 'WAVE-TAB' / Vacuum wavelength
CRVAL1 = 0.0 / [deg] Coordinate value at reference point
CRVAL2 = 0.0 / [deg] Coordinate value at reference point
CRVAL3 = 0.0 / [m] Coordinate value at reference point
PS3_0 = 'FILTERS' / Coordinate transformation parameter
PS3_1 = 'lambda_eff' / Coordinate transformation parameter
LONPOLE = 180.0 / [deg] Native longitude of celestial pole
LATPOLE = 0.0 / [deg] Native latitude of celestial pole
RADESYS = 'ICRS' / Equatorial coordinate system
Here's the code to build it up:
mywcs = wcs.WCS(fobj=hdulist, naxis=3)
pixel_scale = 0.2
mywcs.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'WAVE-TAB']
mywcs.wcs.cunit = ['deg', 'deg', 'm']
mywcs.wcs.cdelt = [
pixel_scale / 60 / 60,
pixel_scale / 60 / 60,
1
]
mywcs.wcs.crpix = [2048, 2048, 0]
mywcs.wcs.crval = [0, 0, 0]
mywcs.wcs.set_ps([(3, 0, 'FILTERS'),
(3, 1, 'lambda_eff')])
mywcs.to_header()
Here's the header of the extension FILTERS:
XTENSION= 'BINTABLE' / binary table extension
BITPIX = 8 / 8-bit bytes
NAXIS = 2 / 2-dimensional binary table
NAXIS1 = 86 / width of table in bytes
NAXIS2 = 19 / number of rows in table
PCOUNT = 0 / size of special data area
GCOUNT = 1 / one data group (required keyword)
TFIELDS = 8 / number of fields in each row
EXTNAME = 'FILTERS ' / name of this binary table extension
TTYPE1 = 'filter ' / label for field
TFORM1 = '30A ' / format of field
TUNIT1 = 'Filter name'
TTYPE2 = 'lambda_eff' / label for field
TFORM2 = 'D ' / format of field
TUNIT2 = 'm '
TTYPE3 = 'ewidth_lambda' / label for field
TFORM3 = 'D ' / format of field
TUNIT3 = 'm '
TTYPE4 = 'ewidth_nu' / label for field
TFORM4 = 'D ' / format of field
TUNIT4 = 'Hz '
TTYPE5 = 'L_lambda_eff0' / label for field
TFORM5 = 'D ' / format of field
TUNIT5 = 'W/m '
TTYPE6 = 'AB_mag0 ' / label for field
TFORM6 = 'D ' / format of field
TTYPE7 = 'L_lambda_eff_nonscatter0' / label for field
TFORM7 = 'D ' / format of field
TUNIT7 = 'W/m '
TTYPE8 = 'AB_mag_nonscatter0' / label for field
TFORM8 = 'D ' / format of field
I'm following (or trying to follow) the explanation of the -TAB standard in Griesen et al. 2006 (PDF). (more links)