diff --git a/examples/tab2spfits.py b/examples/tab2spfits.py index e21230b..75b3f44 100755 --- a/examples/tab2spfits.py +++ b/examples/tab2spfits.py @@ -5,25 +5,27 @@ # tab2spfits.py J17293-CO10.tab # -_version = "24-jan-2023" +_version = "30-jan-2023" _help = """Usage: sp2sdfits [options] TABLE_FILE --h --help Give this help --d --debug Add more debugging info --v --version Report version +-r --restfreq RF Give a (new) RESTFREQ in GHz. Optional. +-h --help Give this help. +-d --debug Add more debugging info. +-v --version Report version. tab2spfits.py converts an ascii table to a simple 1D HDU spectrum fits file -""" - +A restfreq 115.2712018 is common for RSR, but RSR by default does not keep it. +""" import os import sys from astropy.io import fits from astropy.table import Table +from astropy.wcs import WCS import numpy as np import matplotlib.pyplot as plt import warnings @@ -34,12 +36,34 @@ def help(cmd): print("no help here") sys.exit(0) + +def dms2deg(dms): + """ convert a string d:m:s to degrees + """ + sign = 1 + if dms[0] == '-': + sign = -1 + dms = dms[1:] + if dms[0] == '+': + dms = dms[1:] + w = dms.split(':') + deg = float(w[0]) + if len(w) > 1: + deg = deg + float(w[1])/60.0 + if len(w) > 2: + deg = deg + float(w[2])/3600.0 + return sign*deg + if __name__ == "__main__": av = docopt(_help, options_first=True, version='sp2sdfits. %s' % _version) tab = av['TABLE_FILE'] Qdebug = av['--debug'] + if av['--restfreq']: + restfreq = float(av['--restfreq']) * 1e9 # assumed in GHz e.g. 115.2712018 + else: + restfreq = None if not os.path.exists(tab): print("Table %s does not exist" % tab) @@ -55,6 +79,13 @@ def help(cmd): crval1 = freq[naxis1//2]*scale cdelt1 = (freq[1] - freq[0])*scale + wcs = WCS(naxis=1) + wcs.wcs.ctype = ['FREQ'] + wcs.wcs.crpix = [crpix1] + wcs.wcs.crval = [crval1] + wcs.wcs.cdelt = [cdelt1] + print(wcs) + # create a HDU for fits file spname = tab + '.fits' hdu = fits.PrimaryHDU() @@ -64,9 +95,17 @@ def help(cmd): hdu.header['CRPIX1'] = crpix1 hdu.header['CRVAL1'] = crval1 hdu.header['CDELT1'] = cdelt1 + hdu.header['CUNIT1'] = 'Hz' hdu.header['ORIGIN'] = 'LMTOY' - - # reading header (not all tables have a header) + if restfreq != None: + hdu.header['RESTFREQ'] = restfreq + else: + hdu.header['RESTFREQ'] = crval1 + hdu.header['BUNIT'] = 'K' + + # reading header + # *.driver.sum.txt has a header + # *.blanking.sum.txt has no header fp = open(tab) lines = fp.readlines() for line in lines: @@ -82,12 +121,14 @@ def help(cmd): hdu.header['EQUINOX'] = 2000 continue if w[2] == 'RA:': - hdu.header['CRVAL2'] = w[3] # @todo convert to degrees + hdu.header['RA-DMS'] = w[3] + hdu.header['CRVAL2'] = dms2deg(w[3])*15.0 hdu.header['CRPIX2'] = 0.0 hdu.header['CTYPE2'] = 'RA---GLS' continue if w[2] == 'DEC:': - hdu.header['CRVAL3'] = w[3] # @todo convert to degrees + hdu.header['DEC-DMS'] = w[3] + hdu.header['CRVAL3'] = dms2deg(w[3]) hdu.header['CRPIX3'] = 0.0 hdu.header['CTYPE3'] = 'DEC--GLS' continue