[AstroPy] Adding a column to a FITS binary HDU

[E. Madison Bray]

On Fri, Jan 15, 2021 at 5:55 PM Alberto <alberto.manfreda at pi.infn.it> wrote:
>
> Hi,
>
> despite I have been using astropy.io for a few years now, this is my first message to the list (which speaks highly of the quality of the module and of the documentation - great job, really!). Recently, however, I found myself stuck doing an apparently simple task. I have a FITS file containing a bunch of binary HDU tables, which I open with the:
>
> hdu_list = open_event_file('file_path')
>
> command. I need to add a new column to one of the HDU tables, then write the entire hdu_list to a different file location.
>
> What is the simplest way to do this? The methods for adding columns that I found in the online FAQs work for Table objects, but what I can get from my hdu_list is a BinTableHDU object, which doesn't support dictionary-like assignment or the add_column() method. I know I can get a Table object from the hdu_list like this:
>
> table = Table(hdu_list[1].data)
>
> but this is a copy, so the changes to the table does not apply to the original hdu_list.
>
> Any help is appreciated.

Hi Alberto,

The `data` arrays accessible on the HDU objects are basically Numpy
structured arrays, so techniques and utilities for adding a column to
a structured array apply here, though unfortunately it's not as easy
as one would like.  Fortunately, newer versions of Numpy have added a
few utility functions for this; in this case in particular
numpy.lib.recfunctions.append_fields:
https://numpy.org/doc/stable/user/basics.rec.html#numpy.lib.recfunctions.append_fields

So you should be able to add a column like:

>>> from numpy.lib.recfunctions import append_fields
>>> hdu = hdu_list[1]
>>> hdu.data = append_fields(hdu.data, 'NAME_OF_COLUMN', column_data, usemask=False)

I tested that this works (it's a bit annoying that it returns a masked
array by default so you need to give usemask=False) for this to work.

In principle you should also be able to use Table.add_column with the
Table class as you mentioned.  Then you can do:

>>> hdul_list[1].data = table.as_array()

Unfortunately this has a bug related to string columns that I just
discovered in testing this.

I'm sorry that even now this is not as easy to do as it should be.
For the most part the Table class has obviated the need for anything
like this in the "lower-level" interface, but I agree it's not obvious
how to do if you want to manipulate existing FITS files.

Madison