Reading a FITS table without a FITS library


When reading fits files, you usually want to use a FITS library like CFITSIO. Sometimes, though, you want to write a stand-alone program without any external dependencies. FITS is a simple format, and reading one directly is not that hard.

Note that this example is extremely bare-bones; I do not show any error handling, verification that the data is in FITS format, or tolerence (or even detection) of valid differences in how the FITS file might legally be written.

Writing a FITS file is somewhat more challenging, because the standard has many requirements that the reader can safely ignore but the writer cannot. (This example can give you some idea what is involved.)

  1. Start by declaring a structure that mirrors the columns of the FITS table. Each element in the structure must have the same number of bytes as that column in the fits table. For example, for a FITS header with column definitions like this:

    
TFORM1  = '1D      '
    TTYPE1  = 'ra0     '
    TFORM2  = '1D      '
    TTYPE2  = 'dec0    '
    TFORM3  = '1J      '
    TTYPE3  = 'run     '
    TFORM4  = '1J      '
    TTYPE4  = 'rerun   '
    TFORM5  = '1J      '
    TTYPE5  = 'camcol  '
    TFORM6  = '1J      '
    TTYPE6  = 'field   '
    TFORM7  = '1J      '
    TTYPE7  = 'primary '
    
    a data structure might look like this:
    
/** @file sdssfield.h
     * @brief The fits structure mirroring the FITS table row content
     * @author Eric H. Neilsen, Jr.
     */
    
    typedef struct {
      double ra;
      double dec;
      long run;
      long rerun;
      long camcol;
      long field;
      long primary;
    } sdssfield;
    
    
    
    Note that these declarations assume 32-bit long and 64-bit doubles. If you are going to use a compiler that supports the int32_t type, you may consider using that instead of a long type. To get the number of bits you need for a given FITS keyword, see table 8.5 of the FITS standard, v2.0. Derive the corresponding C type you need from your limits.h file.

  2. Next, if we want our program to run on a little-endian machine, like most Intel based computers, you need a function to swap bytes from the big-endian FITS format:

    
/** @file swap_bytes.c
     * @brief Swap bytes to switch between big and little endian
     * @author Eric H. Neilsen, Jr.
     */
    
    #include <string.h>
    
    int
    swap_bytes(void *in,size_t s) 
    {
      char *swapped; /*< Buffer for swapped bytes */
      char *outbyte; /*< Pointer to next byte to be written*/
      char *inbyte;  /*< Pointer to next byte to be read */
    
      swapped=malloc(s);
      inbyte = (char *) in + s;
      for (outbyte=swapped;outbyte<swapped+s;outbyte++) {
        inbyte--;
        *outbyte=*inbyte;
      }
      memcpy(in,swapped,s);
      free(swapped);
      return s;
    }
    
    

  3. Finally, read the file:

    
/** @file simple_fits_read.c
     * @brief Read a fits table of known structure without a fits library
     * @author Eric H. Neilsen, Jr.
     */
    
    #include <stdio.h>
    #include <sysexits.h>
    #include <string.h>
    #include <assert.h>
    #include <sdssfield.h>
    
    /**
     *  @brief The main function, where the work is done
     */
    int
    main (int argc, const char **argv)
    {
      char *filename;
      FILE *fp;      
      char hcard[81];
      char keyword[9];
      char assignment[3];
      char comment[71];
      int card_index;  
      int hpassed=0;   
      long nrows;      
      long row_index;  
      sdssfield row;   
    
      long endian_test = 1;
      int little_endian;  
      little_endian = (int) ( *(char *)&endian_test );        (1)
    
      filename=(char *) argv[1];                              (2)
      fp=fopen(filename,"r");
    
      hcard[80]='\0';                                         (3)
    
      while(hpassed < 2) {                                    (4)
        for (card_index=0; card_index<36; card_index++) {     (5)
          fread(hcard,1,80,fp);
          if( strncmp(hcard,"NAXIS2  ",8) == 0 ) {            (6)
    	sscanf(hcard,"%8s%1s%i%s",
    	       keyword,assignment,&nrows,comment);
          }
          if( strncmp(hcard,"END     ",8) == 0 ) {
    	hpassed++;
          }
        }   
      }
    
      for (row_index=0; row_index<nrows; row_index++) {
        /* It is safer to read each field individually, rather than
           the whole structure, due to memory alignment issues */
        fread(&(row.ra),sizeof(row.ra),1,fp);
        fread(&(row.dec),sizeof(row.dec),1,fp);
        fread(&(row.run),sizeof(row.run),1,fp);
        fread(&(row.rerun),sizeof(row.rerun),1,fp);
        fread(&(row.camcol),sizeof(row.camcol),1,fp);
        fread(&(row.field),sizeof(row.field),1,fp);
        fread(&(row.primary),sizeof(row.primary),1,fp);
        if ( little_endian ) {
          swap_bytes(&(row.ra),sizeof(row.ra));
          swap_bytes(&(row.dec),sizeof(row.dec));
          swap_bytes(&(row.run),sizeof(row.run));
          swap_bytes(&(row.rerun),sizeof(row.rerun));
          swap_bytes(&(row.camcol),sizeof(row.camcol));
          swap_bytes(&(row.field),sizeof(row.field));
          swap_bytes(&(row.primary),sizeof(row.primary));
        }
        printf("%g\t%g\t%d\t%d\t%d\t%d\n",row.ra,row.dec,row.run,row.rerun,row.camcol,row.field);
      }
    
      close(fp);
    
      return EX_OK;
    }
    
    
    (1)
    Check the first byte in a long of 1 to see if we are on a little endian machine.
    (3)
    Header cards, as stored in the FITS file, are 80 unterminated characters. We therefore allocate 81 bytes and add the terminator ourselves.
    (2)
    In a real program, you would add error checking here.
    (4)
    The primary HDU cannot be a table, but we are assuming that the only content of the file is our desired fits table. Therefore, we are looking for the first HDU after the primary.
    (5)
    Header data cards (line of 80 characters in a FITS header) come in groups of 36. If a block isn't full, the remainder will be blanks.
    (6)
    You can read other header keywords similarly.