Intuitively manage column-wise WCS coordinates

[matteobachetti]

In the X-ray data I'm analyzing I have single photons and their X and Y coordinates. The WCS information that allows the conversion to sky coordinates is contained in the header in the form of TXXXXNN keywords, where NN indicates the column and XXXX the wcs parameter (CDLT for cdelta, CRVL for cval etc.). If I read the FITS data structure, I see that this information is assigned to the FITS column

In [15]: data.columns
Out[15]: 
ColDefs(
(...)
    name = 'X'; format = '1I'; unit = 'pixel'; null = -1; coord_type = 'RA---TAN'; coord_unit = 'deg'; coord_ref_point = 500.5; coord_ref_value = 148.9584; coord_inc = -0.0006828076
    name = 'Y'; format = '1I'; unit = 'pixel'; null = -1; coord_type = 'DEC--TAN'; coord_unit = 'deg'; coord_ref_point = 500.5; coord_ref_value = 69.67944; coord_inc = 0.0006828076
)

but

1. apparently it is lost if I actually try to access the fits column (e.g. data['X'] is a numpy array, it has no attribute coord_type or any other)
2. there is no apparent easy way to transform these columns into coordinates (they are just treated as numpy arrays by SkyCoord et al.)

It would be great to have an easy way to get these columns and transform them to SkyCoord objects, but I couldn't find one.
For reference, what I'm doing right now as a workaround is the following, which probably duplicates some functionality already implemented in Astropy but not easily accessible.

def get_wcs_from_col(hdu, col):
    import re
    num_re = re.compile('^[a-zA-Z]+([0-9]+)$')
    header = hdu.header
    ttypes = [k for k in header.keys() if k.startswith('TTYP')]
    coltypes = [header[k].strip() for k in ttypes]
    typ = [ttyp for ttyp, colt in zip(ttypes, coltypes) if colt == col.strip()][0]
    num = num_re.match(typ).group(1)
    res = type('wcsinfo', (), {})()
    res.form = header['TFORM' + num]
    res.lmin = header['TLMIN' + num]
    res.lmax = header['TLMAX' + num]
    res.crval = header['TCRVL' + num]
    res.crpix = header['TCRPX' + num]
    res.cdelt = header['TCDLT' + num]
    res.ctype = header['TCTYP' + num]
    res.cunit = header['TCUNI' + num]
    return res


def get_wcs_from_bintable(hdu, xcol, ycol):
    from astropy import wcs
    header = hdu.header
    xwcs = get_wcs_from_col(hdu, xcol)
    ywcs = get_wcs_from_col(hdu, ycol)

    w = wcs.WCS(naxis=2)

    w.wcs.crpix = [xwcs.crpix, ywcs.crpix]
    w.wcs.cdelt = np.array([xwcs.cdelt, ywcs.cdelt])
    w.wcs.crval = [xwcs.crval, ywcs.crval]
    w.wcs.ctype = [xwcs.ctype, ywcs.ctype]

    return w

Thanks in advance for your help

[saimn]

To answer question 1., you can access these attributes with the io.fits.Column object, e.g. data.columns['X'] (not data['X'] which is a column of the .data recarray). The list of keywords mapped to attributes is here:

In [9]: from astropy.io.fits.column import KEYWORD_TO_ATTRIBUTE                                

In [10]: KEYWORD_TO_ATTRIBUTE                                                                  
Out[10]: 
OrderedDict([('TTYPE', 'name'),
             ('TFORM', 'format'),
             ('TUNIT', 'unit'),
             ('TNULL', 'null'),
             ('TSCAL', 'bscale'),
             ('TZERO', 'bzero'),
             ('TDISP', 'disp'),
             ('TBCOL', 'start'),
             ('TDIM', 'dim'),
             ('TCTYP', 'coord_type'),
             ('TCUNI', 'coord_unit'),
             ('TCRPX', 'coord_ref_point'),
             ('TCRVL', 'coord_ref_value'),
             ('TCDLT', 'coord_inc'),
             ('TRPOS', 'time_ref_pos')])

[eteq]

Aha! I figured out a way to do this in a way that uses the wcs directly (but admittedly lets to an awkward result):

>>> fitsobj[1]
<astropy.io.fits.hdu.table.BinTableHDU at ...>
>>> w = astropy.wcs.WCS(fitsobj[1].header, keysel=['pixel'])
>>> w
WCS Keywords

Number of WCS axes: 6
CTYPE : 'DET1X'  'DET1Y'  'DET2X'  'DET2Y'  'RA---TAN'  'DEC--TAN'  
CRVAL : 0.0  0.0  0.0  0.0  266.4168  -29.00781  
CRPIX : 0.0  0.0  0.0  0.0  500.5  500.5  
PC1_1 PC1_2 PC1_3 PC1_4 PC1_5 PC1_6  : 1.0  0.0  0.0  0.0  0.0  0.0  
PC2_1 PC2_2 PC2_3 PC2_4 PC2_5 PC2_6  : 0.0  1.0  0.0  0.0  0.0  0.0  
PC3_1 PC3_2 PC3_3 PC3_4 PC3_5 PC3_6  : 0.0  0.0  1.0  0.0  0.0  0.0  
PC4_1 PC4_2 PC4_3 PC4_4 PC4_5 PC4_6  : 0.0  0.0  0.0  1.0  0.0  0.0  
PC5_1 PC5_2 PC5_3 PC5_4 PC5_5 PC5_6  : 0.0  0.0  0.0  0.0  1.0  0.0  
PC6_1 PC6_2 PC6_3 PC6_4 PC6_5 PC6_6  : 0.0  0.0  0.0  0.0  0.0  1.0  
CDELT : 0.12096  0.12096  0.12096  0.12096  -0.0006828076  0.0006828076  
NAXIS : 286  49415

so this creates a 6D WCS where the last two dimensions are the actual desired dimensions. That's pretty user hostile in that you then can't e.g. get a SkyCoord by using SkyCoord.from_pixel because it's not a normal celestial coordinate. So maybe this means the work to be done is add a convenience function that pares down this WCS to have just the 2 "regular" (meaning the RA---TAN and DEC--TAN dimensions) axes? Then it would be much easier to just create a WCS and use it with the x/y parts of the table.

[mhvk]

We already have a SkyCoord.from_table method - it would seem to me that either that method or a similar one could look through the table and see if it could find columns with the relevant WCS information, and then use those (i.e., trigger not (just) on column names but rather/also on column attributes)

[hamogu]

I ran into this today, so 👍 to doing something about it. As a related note, all the keywords mentioned above disappear when reading with Table.read(abc.fits). Presumably that's because io.fits treats them special and removes them from the meta data before passing the structure to table.Table.

I think we would want a astropy.wcs.wcs_from_table(coordsys) where coordsys tells you which coordinate system you want if there is more than one defined in the table.

[hamogu]

I found a solution to this (and it's actually pretty easy, because the WCSLIB handles it by itself without python code needed!) about a year ago, then I forgot about it. When I discovered the same problem again, I googled and landed on this issue (which I apparently commented on over a year ago).
So, here is what you need to do:

wcs = WCS(hdus[1].header, keysel=['pixel'], naxis=['celestial'])

This works at least for Chandra and XMM event data.

This definitely should be an example in the documentation!

[hamogu]

Since, at least in my view, that's simply a matter of documenting existing capability, not writing any new code, I'll change the labels on this issue.