astropy: Lookup-table WCS for irregular spectral sampling doesn't look up

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)