Commits
1 1 | import os |
2 2 | import shutil |
3 3 | import numpy as np |
4 4 | |
5 5 | from casatools import ms, table |
6 6 | from casatasks import casalog |
7 7 | from .mstools import write_history |
8 8 | |
9 - | def importfitsidi(fitsidifile,vis,constobsid=None,scanreindexgap_s=None,specframe=None): |
9 + | def importfitsidi(fitsidifile,vis,constobsid=None,scanreindexgap_s=None,specframe=None,coordframe=None): |
10 10 | """Convert FITS-IDI visibility file into a CASA visibility file (MS). |
11 11 | |
12 12 | Keyword arguments: |
13 13 | fitsidifile -- Name(s) of input FITS IDI file(s) |
14 14 | default: None; example='3C273XC1.IDI' or ['3C273XC1.IDI1', '3C273XC1.IDI2'] |
15 15 | vis -- Name of output visibility file (MS) |
16 16 | default: None; example: vis='3C273XC1.ms' |
17 17 | |
18 18 | constobsid -- If True a constant obs id == 0 of is given to all input files |
19 19 | default = False (new obs id for each input file) |
20 20 | |
21 21 | scanreindexgap_s -- if > 0., a new scan is started whenever the gap between two |
22 22 | integrations is > the given value (seconds) or when a new field starts |
23 23 | or when the ARRAY_ID changes. |
24 24 | default = 0. (no reindexing) |
25 25 | |
26 26 | specframe -- this frame will be used to set the spectral reference frame |
27 27 | for all spectral windows in the output MS |
28 28 | default = GEO (geocentric), other options: TOPO, LSRK, BARY |
29 29 | NOTE: if specframe is set to TOPO, the reference location will be taken from |
30 30 | the Observatories table in the CASA data repository for the given name of |
31 - | the observatory. You can edit that table and add new rows. |
32 - | |
31 + | the observatory. You can edit that table and add new rows. |
32 + | |
33 + | coordframe -- this frame will be used to set the reference frame |
34 + | for all source positions in the output MS |
35 + | default = '' (set based on opeoch/equinox), |
36 + | other options: ICRS, B1950, J2000 |
37 + | NOTE: if the epoch/equinox is set to J2000 in the FITS-IDI file, |
38 + | the reference frame will be set to ICRS. To use FK5/J2000 as the |
39 + | reference frame, set coordframe to J2000. |
33 40 | """ |
34 41 | |
35 42 | #Python script |
36 43 | |
37 44 | try: |
38 45 | casalog.origin('importfitsidi') |
39 46 | casalog.post("") |
40 47 | myms = ms() |
41 48 | mytb = table() |
42 49 | |
43 50 | if type(specframe)==str and not specframe=='': |
44 51 | myspecframe=specframe.upper() |
45 52 | else: |
46 53 | myspecframe='GEO' |
47 54 | |
48 55 | refframes = {'REST': 0, 'LSRK': 1, 'LSRD': 2, 'BARY': 3, 'GEO': 4, 'TOPO': 5} |
49 56 | if myspecframe not in refframes: |
50 57 | raise ValueError('Value '+myspecframe+' of parameter specframe invalid. Possible values are REST, LSRK, LSRD, BARY, GEO, TOPO') |
51 58 | |
59 + | if type(coordframe)==str and not coordframe=='': |
60 + | mycoordframe=coordframe.upper() |
61 + | else: |
62 + | mycoordframe='' |
63 + | |
64 + | if mycoordframe not in ['', 'ICRS', 'B1950', 'J2000']: |
65 + | raise ValueError('Value '+mycoordframe+' of parameter coordframe invalid. Possible values are ICRS, B1950, J2000') |
66 + | |
52 67 | if(type(fitsidifile)==str): |
53 68 | casalog.post('### Reading file '+fitsidifile, 'INFO') |
54 69 | myms.fromfitsidi(vis,fitsidifile) |
55 70 | myms.close() |
56 71 | elif(type(fitsidifile)==list): |
57 72 | clist = fitsidifile |
58 73 | casalog.post('### Reading file '+clist[0], 'INFO') |
59 74 | myms.fromfitsidi(vis,clist[0]) |
60 75 | myms.close() |
61 76 | clist.pop(0) |
145 160 | |
146 161 | if myspecframe in refframes: |
147 162 | casalog.post('Setting reference frame for all spectral windows to '+myspecframe, 'INFO') |
148 163 | if myspecframe == 'TOPO': |
149 164 | casalog.post('NOTE: reference position for TOPO frame will be the observatory location', 'WARN') |
150 165 | mytb.open(vis+'/SPECTRAL_WINDOW', nomodify=False) |
151 166 | refcol = mytb.getcol('MEAS_FREQ_REF') |
152 167 | refcol = [refframes[myspecframe]]*len(refcol) |
153 168 | mytb.putcol('MEAS_FREQ_REF', refcol) |
154 169 | mytb.close() |
155 - | |
170 + | |
171 + | mytb.open(vis+'/FIELD', nomodify=False) |
172 + | for colname in ['DELAY_DIR', 'PHASE_DIR', 'REFERENCE_DIR']: |
173 + | kwds = mytb.getcolkeywords(colname) |
174 + | if not mycoordframe=='': |
175 + | kwds['MEASINFO']['Ref'] = mycoordframe |
176 + | # until casacore gets fixed to use ICRS as its default |
177 + | elif mycoordframe=='' and kwds['MEASINFO']['Ref'] == 'J2000': |
178 + | kwds['MEASINFO']['Ref'] = 'ICRS' |
179 + | mytb.putcolkeywords(colname, kwds) |
180 + | mytb.close() |
181 + | |
156 182 | # write history |
157 183 | try: |
158 184 | param_names = importfitsidi.__code__.co_varnames[:importfitsidi.__code__.co_argcount] |
159 185 | vars = locals( ) |
160 186 | param_vals = [vars[p] for p in param_names] |
161 187 | write_history(myms, vis, 'importfitsidi', param_names, param_vals, casalog) |
162 188 | |
163 189 | except Exception as instance: |
164 190 | casalog.post("*** Error \'%s\' updating HISTORY" % instance, 'WARN') |
165 191 | |