# This Python module contains two functions for retrieving or displaying the # full contents of a FITS binary table file. get_fits returns a dictionary # containing all of the header keyword values labeled by the keyword names and # all of the column arrays labled by column name. show_fits displays the HDU, # keyword, and column tree with keyword values. Only column formats are shown # since their arrays can be quite long. A few utilities for retrieving and # parsing the keyword and column information are also defined at the end of # this file. The myfits_get_colval is just a copy of the one in the scan_info # module. # To run from the Python prompt: # >>> path = '/home/gbtdata/AGBT02A_069_01/SpectralProcessor/' # >>> import get_fits # >>> ft = get_fits.get_fits(path + '2003_06_17_09:22:34.fits') # >>> get_fits.show_fits(path + '2003_06_17_09:22:34.fits') # or # >>> from get_fits import * # >>> ft = get_fits(path + '2003_06_17_09:22:34.fits') # >>> show_fits(path + '2003_06_17_09:22:34.fits') import os, string from fitsio import cfitsio import Numeric # This routine returns a dictionary containing all of the header keyword values # labeled by the keyword names and all of the column arrays labled by column # name. Contents of this dictionary are retrieved by HDU extension name and # keyword or column name, e.g., # ft['RECEIVER']['UTCSTART'] # ft['DATA']['Column']['DATA'][0,0,1,0:10] # The first or main HDU is called 'main_header' in the dictionary. The # dictionary method keys() can be used to explore the name space, e.g., # ft.keys() # ft['RECEIVER'].keys() # ft['DATA']['Column'].keys() # (Note that the dictionary names are usually not printed in the order of HDU # or keyword number.) or the show_fits function in this module below can # display the whole name tree in the order that it appears in the FITS file. def get_fits ( file_path ) : # open the specified FITS file for access through cfitsio function calls status, fits = cfitsio.fits_open_file(file_path, cfitsio.READONLY) if status : print 'Cannot open file', file_path return # get the number of HDU's in the file. status, num_hdus = cfitsio.fits_get_num_hdus(fits) if status : print 'Error reading number of HDU\'s' return 0 # create an empty dictionary ft = {} # extract info from each HDU in turn for hdu_num in range(1, num_hdus + 1) : # move cfitsio's attention to this HDU qq, qqq = cfitsio.fits_movabs_hdu(fits, hdu_num) # get the number of keywords in the header status, num_kw, numfree = cfitsio.fits_get_hdrspace(fits) if status : print 'Error reading number of keywords in HDU', hdu_num return 0 # the first HDU doesn't have a name so call it 'main_header' if hdu_num == 1 : hdu_name = 'main_header' else : # for all other HDU's use the extension name status, hdu_name, comment = \ cfitsio.fits_read_key_str(fits, 'EXTNAME') if status : print 'Error reading EXTNAME of HDU', hdu_num return 0 # create an empty subdictionary to contain this HDU's info ft[hdu_name] = {} # since Python loses track of the order of creation of its dictionary # entries, tag the HDU with its serial number in the FITS file ft[hdu_name]['hdu_num'] = hdu_num # extract the keyword name and value from each 80-character FITS # header card and put them in the HDU subdictionary. Don't save # comment nor column information cards. for kw_num in range(1, num_kw + 1) : status, line_str = cfitsio.fits_read_record(fits, kw_num) status, kw, val, comment = parse_card(line_str) if (not is_col_kw(kw)) and (kw != 'COMMENT') : ft[hdu_name][kw] = val # create a special subdictionary for binary table column data ft[hdu_name]['Column'] = {} if hdu_num > 1 : # get the number of binary table columns in this HDU status, num_cols = cfitsio.fits_get_num_cols(fits) if status : print 'Error reading number of columns in HDU', hdu_num return 0 # get the number of rows in each binary table column status, num_rows = cfitsio.fits_get_num_rows(fits) if status : print 'Error reading number of rows in HDU', hdu_num return 0 # read, infer, or otherwise finagle the properties of each # column from the TFORM, TUNIT, and TDIM HDU header keywords # This function assumes that cfitsio is pointed to the current HDU. col_prop = get_col_properties(fits) for col in range(1,num_cols + 1) : col_name = col_prop['name'][col - 1] col_dim = col_prop['dim'][col - 1] repeat = col_prop['repeat'][col - 1] col_form = col_prop['format'][col - 1] # create an empty list for string array columns, othrwise # create an empty Numeric array of the right data type with # enough space for all rows ad the number of data elements in # each row if is_A_format(col_form) : col_arry = [] else : if repeat == 1 : col_arry = Numeric.zeros([num_rows], \ col_prop['ntype'][col - 1]) else : col_arry = Numeric.zeros([num_rows, repeat], \ col_prop['ntype'][col - 1]) # fill the array row by rowusing the appropriate assignment # syntax depending on whether we are filling a list, or a 1-D # or 2-D array for row_num in range(1,num_rows + 1) : x = myfits_get_colval(fits, \ hdu_num, col_name, row_num, 1, repeat) if is_A_format(col_form) : col_arry.append(x) else : if repeat == 1 : col_arry[row_num - 1] = Numeric.array(x) else: col_arry[row_num - 1, :] = Numeric.array(x) # if the data in each row is a multi-dimensional array, reshape # the array using the TDIMn keyword string if (col_dim != '') and (not is_A_format(col_form)) : dim = [num_rows] dim_set = parse_dimension(col_dim) for i in dim_set : dim.append(i) col_arry.shape = dim ft[hdu_name]['Column'][col_name] = col_arry # repeat for each column and HDU # close FITS file access and return the filled dictionary cfitsio.fits_close_file(fits) return ft # show_fits displays the HDU, keyword, and column tree with keyword values. # Only column formats are shown since their arrays can be quite long. def show_fits ( file_path ) : # open the specified FITS file for access through cfitsio function calls status, fits = cfitsio.fits_open_file(file_path, cfitsio.READONLY) if status : print 'Cannot open file', file_path return # get the number of HDU's in the file. status, num_hdus = cfitsio.fits_get_num_hdus(fits) if status : print 'Error reading number of HDU\'s' return 0 # extract info from each HDU in turn for hdu_num in range(1, num_hdus + 1) : # move cfitsio's attention to this HDU qq, qqq = cfitsio.fits_movabs_hdu(fits, hdu_num) # get the number of keywords in the header status, num_kw, numfree = cfitsio.fits_get_hdrspace(fits) if status : print 'Error reading number of keywords in HDU', hdu_num return 0 # the first HDU doesn't have a name so call it 'main_header' if hdu_num == 1 : hdu_name = 'main_header' else : # for all other HDU's use the extension name status, hdu_name, comment = \ cfitsio.fits_read_key_str(fits, 'EXTNAME') if status : print 'Error reading EXTNAME of HDU', hdu_num return 0 # start printing the tree print hdu_name, ':' for kw_num in range(1, num_kw + 1) : # extract the keyword name and value from each 80-character FITS # header card and print it. Don't print comment cards. status, line_str = cfitsio.fits_read_record(fits, kw_num) status, kw, val, comment = parse_card(line_str) if (not is_col_kw(kw)) and (kw != 'COMMENT') : if len(kw) < 6 : print ' ' + kw + '\t\t:', val else: print ' ' + kw + '\t:', val # The first HDU never has a binary table so don't look for columns if hdu_num > 1 : # get the number of binary table columns in this HDU status, num_cols = cfitsio.fits_get_num_cols(fits) if status : print 'Error reading number of columns in HDU', hdu_num return 0 # get the number of rows in each binary table column status, num_rows = cfitsio.fits_get_num_rows(fits) if status : print 'Error reading number of rows in HDU', hdu_num return 0 print ' ', num_cols, 'Columns,', num_rows, 'Rows :' # read, infer, or otherwise finagle the properties of each # column from the TFORM, TUNIT, and TDIM HDU header keywords # This function assumes that cfitsio is pointed to the current HDU. col_prop = get_col_properties(fits) for col in range(1,num_cols + 1) : col_name = col_prop['name'][col - 1] col_dim = col_prop['dim'][col - 1] repeat = col_prop['repeat'][col - 1] col_form = col_prop['format'][col - 1] print ' ' + col_name + '\t(' + col_form + ') ' + col_dim # close the FITS file cfitsio.fits_close_file(fits) return # I couldn't find a cfitsio function that would return the data type of a # of a keyword value (maybe there is one) so I wrote my own card parser. # This isn't guaranteed to work on all FITS cards, but it works on all # encountered so far in the GBT FITS files. This function returns the # keyword name, the value in the appropriate Python data type, and the # comment string. def parse_card ( card_str ) : # is it a COMENT card? If so, just return the COMMENT keyword and the # comment string - no value com_chk = string.split(card_str) if com_chk[0] == 'COMMENT' : return (0, com_chk[0], 0, string.join(com_chk[1:], ' ')) # split the card into two strings, before and after the '=' character # before is the keyword name, after is the value and comment string s1 = string.split(card_str, '=') # remove blank spaces from the keyword name kw = string.strip(s1[0]) # divide the remaining string into keyword value and comment strings and # remove unwanted whitespace at the ends s2 = string.split(s1[1], '/') val_str = string.rstrip(string.lstrip(s2[0])) comment = string.rstrip(string.lstrip(s2[1])) # if the keyword value is enclosed in single quotes, it is a string # data type if string.find(val_str, "'") >= 0 : # remove the quotes and whitespace on the ends temp = string.split(val_str, "'") val = string.rstrip(string.lstrip(temp[1])) # otherwise, if the value contains an exponential character, it is a # floating point number elif (string.find(val_str, 'E') >= 0) or (string.find(val_str, 'e') >= 0) : val = string.atof(val_str) # otherwise, if the value is a single T or F character, it is a # logical True or False elif (val_str == 'T') : val = 1 elif (val_str == 'F') : val = 0 # if all else fails, it must be an integer else : val = string.atoi(val_str) return (0, kw, val, comment) # test the keyword name to see if it is one of the reserved names for binary # table interpretation. Return True or False. def is_col_kw ( kw_str ) : if (kw_str[:5] == 'TTYPE') or (kw_str[:5] == 'TFORM') or \ (kw_str[:5] == 'TUNIT') or (kw_str[:4] == 'TDIM') : return 1 return 0 # this function determines the properties of data in a binary table from the # FITS keywords reserved for this purpose. It returns a dictionary of lists # (one list item per column) with thedictionary members being column name, # units of the data (K, MHz...), column format in FITS designations (1J, 4A...), # column dimensions if it is not a 1-dimensional array, the data type in # cfitsio codes (cfitsio.TDOUBLE, cfitsio.TSTRING...), the Numeric module data # type (Numeric.Int32, Numeric.Float32...), and the number of items in the # column entry. This function assume that cfitsio is pointed at the desired # HDU def get_col_properties ( fits ) : # get the number of keyword in this HDU status, num_kw, numfree = cfitsio.fits_get_hdrspace(fits) if status : print 'Error reading number of keywords in HDU' return 0 # create a set of empty lists that will be return by this function col_name = [] col_unit = [] col_format = [] col_dim = [] col_cfits_type = [] col_numeric_type = [] col_repeat = [] # fill these lists with null values to begin status, num_cols = cfitsio.fits_get_num_cols(fits) for col_num in range(num_cols) : col_name.append('') col_unit.append('') col_format.append('') col_dim.append('') col_cfits_type.append(0) col_numeric_type.append(0) col_repeat.append(0) # scan all keywords in the header for table info and extract when found for kw_num in range(1, num_kw + 1) : # read and parse each keyword card status, line_str = cfitsio.fits_read_record(fits, kw_num) status, kw, val, comment = parse_card(line_str) # if it is a column type keyword, if kw[:5] == 'TTYPE' : # extract the column number col_num = string.atoi(kw[5:]) - 1 # save the column name col_name[col_num] = val # get and save the cfitsio data type status, typecode, repeat, width = \ cfitsio.fits_get_coltype(fits, col_num + 1, 0) if status : print 'Error reading column', col_num + 1, 'type' return 0 col_cfits_type[col_num] = typecode # convert and save this as a Numeric module data type col_numeric_type[col_num] = cfts_to_numeric_type(typecode) # save the number of items in the column entry col_repeat[col_num] = repeat # if it is a column format keyword, elif kw[:5] == 'TFORM' : # extract the column number, and save the format string col_num = string.atoi(kw[5:]) - 1 col_format[col_num] = val # if it is a column units keyword, elif kw[:5] == 'TUNIT' : # extract the column number, and save the units string col_num = string.atoi(kw[5:]) - 1 col_unit[col_num] = val # if it is a column dimension keyword, elif kw[:4] == 'TDIM' : # extract the column number, and save the dimension string, # e.g., (1024,2,2). Only multidimensional arrays have dimension # values. Otherwise, it's the same as the number of repeats. col_num = string.atoi(kw[4:]) - 1 col_dim[col_num] = val # assemble and return the dictionary of lists return {'name' : col_name, 'unit' : col_unit, 'format' : col_format, 'dim' : col_dim, 'typecode' : col_cfits_type, 'ntype' : col_numeric_type, 'repeat' : col_repeat} # This function returns a list value from the specified HDU binary table column # designated by the column name. The arguments are the cfitsio fits file # handle, the HDU number (1 to ...), the column name string, the table row # number, the index (1 to ...) of the first value to be retrieved from the # table row/column entry array, and the number of values to be retrieved. If # num_vals == 0 or larger than the number available, all values from # first_index to the end of the array are retrieved. def myfits_get_colval ( fits, hdu_num, col_name, row_num, first_index, num_vals ) : # get the number of HDU's in this fits file and check to be sure that # the requested HDU number is within range status, nhdus = cfitsio.fits_get_num_hdus(fits) if status : print 'Error reading number of HDU\'s' return 0 if (hdu_num < 1 or hdu_num > nhdus) : if print_err: print 'HDU number:', hdu_num, 'out of range: 1 to', nhdus return 0 # switch cfitsio's attention to the specified HDU q, qq = cfitsio.fits_movabs_hdu(fits, hdu_num) # tell cfitsio to use case sensitivity in finding column by name casesens = 1 # get the column number from the column name status, col_num = cfitsio.fits_get_colnum(fits, casesens, col_name) if status : print 'Cannot find', col_name, 'column in HDU number', hdu_num return 0 # get the column's data typeand the number of elements in the array (repeat) status, typecode, repeat, width = cfitsio.fits_get_coltype(fits, col_num, 0) if status : print 'Error reading', col_name, 'column type' return 0 if num_vals < 1 : num_vals = repeat if repeat < num_vals : print 'Fewer values (', repeat, ') in', col_name, 'column than the', \ num_vals, 'requested' num_vals = repeat # use the column's data type to control which retrieval function to use if typecode == cfitsio.TDOUBLE : status, val = cfitsio.fits_read_col_dbl(fits, col_num, row_num, \ first_index, num_vals, 0) elif typecode == cfitsio.TFLOAT : status, val = cfitsio.fits_read_col_flt(fits, col_num, row_num, \ first_index, num_vals, 0) elif typecode == cfitsio.TINT : status, val = cfitsio.fits_read_col_int(fits, col_num, row_num, \ first_index, num_vals, 0) elif typecode == cfitsio.TUINT : status, val = cfitsio.fits_read_col_uint(fits, col_num, row_num, \ first_index, num_vals, 0) elif typecode == cfitsio.TLONG : status, val = cfitsio.fits_read_col_lng(fits, col_num, row_num, \ first_index, num_vals, 0) elif typecode == cfitsio.TULONG : status, val = cfitsio.fits_read_col_ulng(fits, col_num, row_num, \ first_index, num_vals, 0) elif typecode == cfitsio.TSHORT : status, val = cfitsio.fits_read_col_sht(fits, col_num, row_num, \ first_index, num_vals, 0) elif typecode == cfitsio.TUSHORT : status, val = cfitsio.fits_read_col_usht(fits, col_num, row_num, \ first_index, num_vals, 0) elif typecode == cfitsio.TSTRING : num_vals = 1 status, val = cfitsio.fits_read_col_str(fits, col_num, row_num, \ first_index, num_vals, 0) elif typecode == cfitsio.TBYTE : status, val = cfitsio.fits_read_col_byt(fits, col_num, row_num, \ first_index, num_vals, 0) else : print 'Unrecognized', col_name, 'column type:', typecode return 0 if status : print 'Error reading', col_name, 'column in HDU number', hdu_num print ' typecode:', typecode, 'repeat:', repeat, 'width:', width return 0 return val # this routine does a straightforward conversion from cfitsio data type codes # to Numeric module type codes def cfts_to_numeric_type ( typecode ) : if typecode == cfitsio.TDOUBLE : return Numeric.Float64 elif typecode == cfitsio.TFLOAT : return Numeric.Float32 elif typecode == cfitsio.TINT : return Numeric.Int elif typecode == cfitsio.TUINT : return Numeric.UInt elif typecode == cfitsio.TLONG : return Numeric.Int elif typecode == cfitsio.TULONG : return Numeric.UnsignedInt32 elif typecode == cfitsio.TSHORT : return Numeric.Int16 elif typecode == cfitsio.TUSHORT : return Numeric.UnsignedInt16 elif typecode == cfitsio.TSTRING : return Numeric.Character elif typecode == cfitsio.TBYTE : return Numeric.UnsignedInt8 else : print 'Unrecognized', col_name, 'column type:', typecode return 0 # this routine parses a TDIMn dimension string, e.g., '(1024,2,2)' into # an integer list of dimensions def parse_dimension ( col_dim ) : dlist = string.split(string.replace(string.replace(col_dim, '(', ''), \ ')', ''), ',') # the GBT list of dimensions is reversed from the way we'd like to use them # in a Numeric array dlist.reverse() dim = [] for dstr in dlist : dim.append(string.atoi(dstr)) return dim # this function tests a column format value string for the 'A' FITS # column format def is_A_format( col_form ) : if string.find(col_form, 'A') >= 0 : return 1 return 0 # To test, uncover one of the calls below and from the Python prompt # execute either: # import get_fits # or # reload(get_fits) # ft = get_fits('/home/gbtdata/AGBT02A_069_01/SpectralProcessor/2003_06_17_09:22:34.fits') # print ft['DATA']['Column']['DATA'].shape #print ft.keys(), '\n' #for k in ft.keys() : # print k, ':', ft[k].keys() # if ft[k].has_key('Column') : # print '----' # print ft[k]['Column'].keys() # print '\n' #show_fits('/home/gbtdata/AGBT02A_069_01/LO1A/2003_06_17_09:22:34.fits')