How to generate a fits file from the beginning

Question

In this post, they explain how to generate a fits file from ascii file. However, I also would like to know how to define header and data into fits file. (Converting ASCII Table to FITS image)

For example, when I call a spectral fits file with astropy (which is downloaded from a telescope), I can call data and header separately.

I.E

  In [1]:hdu = fits.open('observation.fits', memmap=True)

  In [2]:header = hdu[0].header

  In [3]:header 
   Out [3]: 
   SIMPLE  =                    T / conforms to FITS standard                      
   BITPIX  =                     8                                                 
   NAXIS   =                     1
   NAXIS1  =                     47356                      
   EXTEND  =                    T                                                  
   DATE    = 'date' / file creation date (YYYY-MM-DDThh:mm:ss UT)   
   ORIGIN  = 'XXX     '           / European Southern Observatory                  
   TELESCOP= 'XXX'         / ESO Telescope Name                             
   INSTRUME= 'Instrument'           / Instrument used.                               
   OBJECT  = 'ABC '           / Original target.                               
   RA      =           30.4993 / xx:xx:xx.x RA (J2000) pointing                 
   DEC     =            -20.0009 / xx:xx:xx.x DEC (J2000) pointing               

   CTYPE1  = 'WAVE    '           / wavelength axis in nm                          
   CRPIX1  =                   0. / Reference pixel in z                           
   CRVAL1  =     298.903594970703 / central wavelength                             
   CDELT1  =   0.0199999995529652 / nm per pixel                                   
   CUNIT1  = 'nm      '           / spectral unit                                  

  ..
   bla bla
  ..                                                                               
   END                                                                             




   In [3]:data = hdu[0].data
   In [4]:data 
   Out [4]:array([  1000,   1001,   1002, ...,
     5.18091546e-13,   4.99434453e-13,   4.91280864e-13])


   Lets assume, I have data like below

   WAVE FLUX
   1000 2.02e-12
   1001 3.03e-12
   1002 4.04e-12
  ..
   bla bla
  ..

So, I'd like to generate a spectral fits file with my own data (with its own header).

Mini question : Now lets assume, I generate spectral fits file correctly, but I realised that I forgot to take logarithm of WAVE values in X axis (1000, 1001, 1002, ....) . How can I do that without touching FLUX values of Y-axis (2.02e-12, 3.03e-13, 4.04e-13) ?

Answer 1

FITS files are organized as one or more HDUs (Header Data Units) consisting, as the name suggests, as one data object (generally, a single array for an observation, though sometimes something else like a table), and the header of metadata that goes with that data.

To create a file from scratch, especially an image, the simplest way is to directly create an ImageHDU object:

>>> from astropy.io import fits
>>> hdu = fits.ImageHDU()

Just as with an HDU read from an existing file, this HDU has a (mostly empty) header, and an empty data attribute that you can then assign to:

>>> hdu.data = np.array(<some numpy array>)
>>> hdu.header['TELESCOP'] = 'Gemini'

When you're satisfied you can write the HDU out to a file with:

>>> hdu.writeto('filename.fits')

(Note: A lot of the documentation you'll see demonstrates a more complex process of creating an HDUList object, appending the HDU to the HDU list, and then writing the full HDU list. This is only necessary if you're creating a multi-extension FITS file. For a single HDU, you can use hdu.writeto directly and the framework will handle the other structural details.)

In general you don't need to manipulate the headers that describe the format of the data itself--that is automatic and should not be touched by hand (FITS has the unfortunate misfeature of mixing information about data structure with actual metadata). You can see more examples on how to manipulate FITS data here: http://docs.astropy.org/en/stable/generated/examples/index.html#astropy-io

Your other question pertains to manipulating the WCS (World Coordinate System) of the image, and in particular for spectral data this can be non-trivial. I would ask a separate question about that with more details about what you hope to accomplish.