File: readuvfits.py

package info (click to toggle)
purify 5.0.1-2
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 186,836 kB
  • sloc: cpp: 17,731; python: 510; xml: 182; makefile: 7; sh: 6
file content (44 lines) | stat: -rw-r--r-- 1,914 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
import astropy.io.fits as fits
import numpy as np

speed_of_light = 299792458. #m/s

def readData(filename, vis_name, pol):
    hdu = fits.open(filename)
    #reading in uvdata
    print "Header..."
    print hdu[0].header
    print hdu[0].data['DATA'].shape

    no_chan = hdu[0].data['DATA'].shape[3]
    no_pol =  hdu[0].data['DATA'].shape[4]

    re = np.array([])
    im = np.array([])
    sigma = np.array([])
    u = np.array([])
    v = np.array([])
    for chan in range(no_chan):
        freq = hdu[0].header["CRVAL4"] + chan * hdu[0].header["CDELT4"]
        
        if (np.count_nonzero(hdu[0].data['DATA'][:, 0, 0, chan, pol, 0]**2 + hdu[0].data['DATA'][:, 0, 0, chan, pol, 1]**2) > 0):
            print "loading frequency ", freq, chan
            flags = hdu[0].data['DATA'][:, 0, 0, chan, pol, 2]/np.abs(hdu[0].data['DATA'][:, 0, 0, chan, pol, 2]) #I think this gives the flags....
            sigma = np.concatenate((sigma, 1./np.sqrt(hdu[0].data['DATA'][:, 0, 0, chan, pol, 2][flags > 0])))
            re = np.concatenate((re, hdu[0].data['DATA'][:, 0, 0, chan, pol, 0][flags > 0]))
            im = np.concatenate((im, hdu[0].data['DATA'][:, 0, 0, chan, pol, 1][flags > 0]))
            #reading in uv-coordinates
            u = np.concatenate((u, hdu[0].data['UU'][flags > 0] * freq))
            v = np.concatenate((v, hdu[0].data['VV'][flags > 0] * freq))
    print "Total visibilities... ", u.shape, v.shape, re.shape, im.shape, sigma.shape
    #reading in weights
    table = np.column_stack((u, v, re, im, sigma))
    print table[0,:], u[0], v[0], re[0], im[0], sigma[0]
    np.savetxt(vis_name, table, delimiter = " ")

names = ["0332-391", "0114-476", "1637-77", "sumss0515"]
for i in range(len(names)):
    uv_fits = "/Users/luke/Radio_Data/" + names[i] +".uvfits"
    output_vis = names[i] + ".vis"
    readData(uv_fits, output_vis, 0)