Line_Profiler_full.out
This file contains the line by line timing results as run on the Oscar cluster at Brown. The task involved reading a 3.1G MIRIAD file, writing it out to miriad, phasing, writing to uvfits, then reading in the uvfits file and unphasing it.
Timer unit: 1e-06 s
Total time: 286.273 s
File: /users/alanman/optimize_pyuvdata/pyuvdata/pyuvdata/miriad.py
Function: read_miriad at line 28
Line # Hits Time Per Hit % Time Line Contents
==============================================================
28 @profile
29 def read_miriad(self, filepath, correct_lat_lon=True, run_check=True, run_check_acceptability=True):
30 """
31 Read in data from a miriad file.
32
33 Args:
34 filepath: The miriad file directory to read from.
35 correct_lat_lon: flag -- that only matters if altitude is missing --
36 to update the latitude and longitude from the known_telescopes list
37 run_check: Option to check for the existence and proper shapes of
38 required parameters after reading in the file. Default is True.
39 run_check_acceptability: Option to check acceptable range of the values of
40 required parameters after reading in the file. Default is True.
41 """
42 1 11715 11715.0 0.0 if not os.path.exists(filepath):
43 raise(IOError, filepath + ' not found')
44 1 60578 60578.0 0.0 uv = a.miriad.UV(filepath)
45
46 1 5 5.0 0.0 miriad_header_data = {'Nfreqs': 'nchan',
47 1 4 4.0 0.0 'Npols': 'npol',
48 1 4 4.0 0.0 'integration_time': 'inttime',
49 1 3 3.0 0.0 'channel_width': 'sdf', # in Ghz!
50 1 4 4.0 0.0 'object_name': 'source',
51 # NB: telescope_name and instrument are treated
52 # as the same
53 1 4 4.0 0.0 'telescope_name': 'telescop',
54 1 4 4.0 0.0 'instrument': 'telescop'
55 }
56 8 32 4.0 0.0 for item in miriad_header_data:
57 7 48 6.9 0.0 if isinstance(uv[miriad_header_data[item]], str):
58 3 22 7.3 0.0 header_value = uv[miriad_header_data[item]].replace('\x00', '')
59 else:
60 4 24 6.0 0.0 header_value = uv[miriad_header_data[item]]
61 7 45 6.4 0.0 setattr(self, item, header_value)
62
63 1 6 6.0 0.0 latitude = uv['latitud'] # in units of radians
64 1 6 6.0 0.0 longitude = uv['longitu']
65 1 4 4.0 0.0 try:
66 1 5 5.0 0.0 altitude = uv['altitude']
67 1 58 58.0 0.0 self.telescope_location_lat_lon_alt = (latitude, longitude, altitude)
68 except(KeyError):
69 # get info from known telescopes. Check to make sure the lat/lon values match reasonably well
70 telescope_obj = uvtel.get_telescope(self.telescope_name)
71 if telescope_obj is not False:
72
73 tol = 2 * np.pi * 1e-3 / (60.0 * 60.0 * 24.0) # 1mas in radians
74 lat_close = np.isclose(telescope_obj.telescope_location_lat_lon_alt[0],
75 latitude, rtol=0, atol=tol)
76 lon_close = np.isclose(telescope_obj.telescope_location_lat_lon_alt[1],
77 longitude, rtol=0, atol=tol)
78 if correct_lat_lon:
79 self.telescope_location_lat_lon_alt = telescope_obj.telescope_location_lat_lon_alt
80 else:
81 self.telescope_location_lat_lon_alt = (latitude, longitude, telescope_obj.telescope_location_lat_lon_alt[2])
82 if lat_close and lon_close:
83 if correct_lat_lon:
84 warnings.warn('Altitude is not present in Miriad file, '
85 'using known location values for '
86 '{telescope_name}.'.format(telescope_name=telescope_obj.telescope_name))
87 else:
88 warnings.warn('Altitude is not present in Miriad file, '
89 'using known location altitude value '
90 'for {telescope_name} and lat/lon from '
91 'file.'.format(telescope_name=telescope_obj.telescope_name))
92 else:
93 warn_string = ('Altitude is not present in file ')
94 if not lat_close and not lon_close:
95 warn_string = warn_string + 'and latitude and longitude values do not match values '
96 else:
97 if not lat_close:
98 warn_string = warn_string + 'and latitude value does not match value '
99 else:
100 warn_string = warn_string + 'and longitude value does not match value '
101 if correct_lat_lon:
102 warn_string = (warn_string + 'for {telescope_name} in known '
103 'telescopes. Using values from known '
104 'telescopes.'.format(telescope_name=telescope_obj.telescope_name))
105 warnings.warn(warn_string)
106 else:
107 warn_string = (warn_string + 'for {telescope_name} in known '
108 'telescopes. Using altitude value from known '
109 'telescopes and lat/lon from '
110 'file.'.format(telescope_name=telescope_obj.telescope_name))
111 warnings.warn(warn_string)
112 else:
113 warnings.warn('Altitude is not present in Miriad file, and ' +
114 'telescope {telescope_name} is not in ' +
115 'known_telescopes. Telescope location not '
116 'set.'.format(telescope_name=self.telescope_name))
117
118 1 1061 1061.0 0.0 self.history = uv['history']
119 1 10 10.0 0.0 if self.pyuvdata_version_str not in self.history.replace('\n', ''):
120 1 8 8.0 0.0 self.history += self.pyuvdata_version_str
121 1 7 7.0 0.0 self.channel_width *= 1e9 # change from GHz to Hz
122
123 # read through the file and get the data
124 1 8 8.0 0.0 _source = uv['source'] # check source of initial visibility
125 1 5 5.0 0.0 data_accumulator = {}
126 1 4 4.0 0.0 pol_list = []
127 924673 18606339 20.1 6.5 for (uvw, t, (i, j)), d, f in uv.all(raw=True):
128 # control for the case of only a single spw not showing up in
129 # the dimension
130 # Note that the (i, j) tuple is calculated from a baseline number in
131 # aipy (see miriad_wrap.h). The i, j values are also adjusted by aipy
132 # to start at 0 rather than 1.
133 924672 4337237 4.7 1.5 if len(d.shape) == 1:
134 924672 5044929 5.5 1.8 d.shape = (1,) + d.shape
135 924672 6079467 6.6 2.1 self.Nspws = d.shape[0]
136 924672 8114680 8.8 2.8 self.spw_array = np.arange(self.Nspws)
137 else:
138 raise(ValueError, """Sorry. Files with more than one spectral
139 window (spw) are not yet supported. A great
140 project for the interested student!""")
141 924672 3594801 3.9 1.3 try:
142 924672 8281805 9.0 2.9 cnt = uv['cnt']
143 except(KeyError):
144 cnt = np.ones(d.shape, dtype=np.float)
145 924672 5535689 6.0 1.9 ra = uv['ra']
146 924672 5230173 5.7 1.8 dec = uv['dec']
147 924672 5241237 5.7 1.8 lst = uv['lst']
148 924672 5383193 5.8 1.9 source = uv['source']
149 924672 3729196 4.0 1.3 if source != _source:
150 raise(ValueError, 'This appears to be a multi source file, which is not supported.')
151 else:
152 924672 3573790 3.9 1.2 _source = source
153
154 924672 3526011 3.8 1.2 try:
155 924672 5529045 6.0 1.9 data_accumulator[uv['pol']].append([uvw, t, i, j, d, f, cnt,
156 924670 5690368 6.2 2.0 ra, dec])
157 2 10 5.0 0.0 except(KeyError):
158 2 8 4.0 0.0 data_accumulator[uv['pol']] = [[uvw, t, i, j, d, f, cnt,
159 2 25 12.5 0.0 ra, dec]]
160 2 12 6.0 0.0 pol_list.append(uv['pol'])
161 # NB: flag types in miriad are usually ints
162
163 3 116059 38686.3 0.0 for pol,data in data_accumulator.iteritems():
164 2 960005 480002.5 0.3 data_accumulator[pol] = np.array(data)
165
166 1 44 44.0 0.0 self.polarization_array = np.array(pol_list)
167 1 11 11.0 0.0 if len(self.polarization_array) != self.Npols:
168 warnings.warn('npols={npols} but found {l} pols in data file'.format(
169 npols=self.Npols, l=len(self.polarization_array)))
170
171 # makes a data_array (and flag_array) of zeroes to be filled in by
172 # data values
173 # any missing data will have zeros
174
175 # use set to get the unique list of all times ever listed in the file
176 # iterate over polarizations and all spectra (bls and times) in two
177 # nested loops, then flatten into a single vector, then set
178 # then list again.
179
180 1 4 4.0 0.0 times = list(set(
181 924675 3855604 4.2 1.3 np.concatenate([[k[1] for k in d] for d in data_accumulator.values()])))
182 1 4357 4357.0 0.0 times = np.sort(times)
183
184 1 4 4.0 0.0 ant_i_unique = list(set(
185 924675 3710282 4.0 1.3 np.concatenate([[k[2] for k in d] for d in data_accumulator.values()])))
186 1 8 8.0 0.0 ant_j_unique = list(set(
187 924675 3692915 4.0 1.3 np.concatenate([[k[3] for k in d] for d in data_accumulator.values()])))
188
189 1 40 40.0 0.0 sorted_unique_ants = sorted(list(set(ant_i_unique + ant_j_unique)))
190 1 29 29.0 0.0 ant_i_unique = np.array(ant_i_unique)
191 1 16 16.0 0.0 ant_j_unique = np.array(ant_j_unique)
192
193 # Determine maximum digits needed to distinguish different values
194 1 86 86.0 0.0 ndig_ant = np.ceil(np.log10(sorted_unique_ants[-1])).astype(int) + 1
195 # Be excessive in precision because we use the floating point values as dictionary keys later
196 1 39 39.0 0.0 prec_t = - 2 * np.floor(np.log10(self._time_array.tols[-1])).astype(int)
197 1 15 15.0 0.0 ndig_t = (np.ceil(np.log10(times[-1])).astype(int) + prec_t + 2)
198 1 4 4.0 0.0 blts = []
199 3 16 5.3 0.0 for d in data_accumulator.values():
200 924674 3863501 4.2 1.3 for k in d:
201 924672 6986079 7.6 2.4 blt = ["{1:.{0}f}".format(prec_t, k[1]).zfill(ndig_t),
202 924672 4985296 5.4 1.7 str(k[2]).zfill(ndig_ant), str(k[3]).zfill(ndig_ant)]
203 924672 3945795 4.3 1.4 blt = "_".join(blt)
204 924672 3696188 4.0 1.3 blts.append(blt)
205 1 346633 346633.0 0.1 reverse_inds= np.argsort(np.unique(np.array(blts),return_index=True)[1])
206 1 1577533 1577533.0 0.6 unique_blts = np.sort(np.unique(np.array(blts)))
207
208 1 1632683 1632683.0 0.6 reverse_inds = dict(zip(np.unique(np.array(blts)),reverse_inds))
209 1 29 29.0 0.0 self.Nants_data = len(sorted_unique_ants)
210
211 # Miriad has no way to keep track of antenna numbers, so the antenna
212 # numbers are simply the index for each antenna in any array that
213 # describes antenna attributes (e.g. antpos for the antenna_postions).
214 # Therefore on write, nants (which gives the size of the antpos array)
215 # needs to be increased to be the max value of antenna_numbers+1 and the
216 # antpos array needs to be inflated with zeros at locations where we
217 # don't have antenna information. These inflations need to be undone at
218 # read. If the file was written by pyuvdata, then the variable antnums
219 # will be present and we can use it, otherwise we need to test for zeros
220 # in the antpos array and/or antennas with no visibilities.
221 1 5 5.0 0.0 try:
222 # The antnums variable will only exist if the file was written with pyuvdata.
223 # For some reason Miriad doesn't handle an array of integers properly,
224 # so we convert to floats on write and back here
225 1 43 43.0 0.0 self.antenna_numbers = uv['antnums'].astype(int)
226 1 10 10.0 0.0 self.Nants_telescope = len(self.antenna_numbers)
227 except(KeyError):
228 self.antenna_numbers = None
229 self.Nants_telescope = None
230
231 1 7 7.0 0.0 nants = uv['nants']
232 1 4 4.0 0.0 try:
233 1 24 24.0 0.0 antpos = uv['antpos'].reshape(3, nants).T
234 if self.Nants_telescope is not None:
235 # in this case there is an antnums variable
236 # (meaning that the file was written with pyuvdata), so we'll use it
237 if nants == self.Nants_telescope:
238 # no inflation, so just use the positions
239 self.antenna_positions = antpos
240 else:
241 # there is some inflation, just use the ones that appear in antnums
242 self.antenna_positions = np.zeros((self.Nants_telescope, 3), dtype=antpos.dtype)
243 for ai, num in enumerate(self.antenna_numbers):
244 self.antenna_positions[ai, :] = antpos[num, :]
245 else:
246 # there is no antnums variable (meaning that this file was not
247 # written by pyuvdata), so we test for antennas with non-zero
248 # positions and/or that appear in the visibility data
249 # (meaning that they have entries in ant_1_array or ant_2_array)
250 antpos_length = np.sqrt(np.sum(np.abs(antpos)**2, axis=1))
251 good_antpos = np.where(antpos_length > 0)[0]
252 # take the union of the antennas with good positions (good_antpos)
253 # and the antennas that have visisbilities (sorted_unique_ants)
254 # if there are antennas with visibilities but zeroed positions we issue a warning below
255 ants_use = set(good_antpos).union(sorted_unique_ants)
256 # ants_use are the antennas we'll keep track of in the UVData
257 # object, so they dictate Nants_telescope
258 self.Nants_telescope = len(ants_use)
259 self.antenna_numbers = np.array(list(ants_use))
260 self.antenna_positions = np.zeros((self.Nants_telescope, 3), dtype=antpos.dtype)
261 for ai, num in enumerate(self.antenna_numbers):
262 self.antenna_positions[ai, :] = antpos[num, :]
263 if antpos_length[num] == 0:
264 warnings.warn('antenna number {n} has visibilities '
265 'associated with it, but it has a position'
266 ' of (0,0,0)'.format(n=num))
267 1 4 4.0 0.0 except(KeyError):
268 # there is no antpos variable
269 1 7 7.0 0.0 self.antenna_positions = None
270
271 1 5 5.0 0.0 if self.antenna_numbers is None:
272 # there are no antenna_numbers or antenna_positions, so just use
273 # the antennas present in the visibilities
274 # (Nants_data will therefore match Nants_telescope)
275 self.antenna_numbers = np.array(sorted_unique_ants)
276 self.Nants_telescope = len(self.antenna_numbers)
277
278 # antenna names is a foreign concept in miriad but required in other formats.
279 1 4 4.0 0.0 try:
280 # Here we deal with the way pyuvdata tacks it on to keep the name information if we have it:
281 # This is a horrible hack to save & recover antenna_names array. Miriad can't handle arrays
282 # of strings or a long enough single string to put them all into one string
283 # so we convert them into hex values and then into floats on write and convert
284 # back to strings here
285 1 9 9.0 0.0 ant_name_flt = uv['antnames']
286 129 1309 10.1 0.0 ant_name_list = [('%x' % elem.astype(np.int64)).decode('hex') for elem in ant_name_flt]
287 1 8 8.0 0.0 self.antenna_names = ant_name_list
288 except(KeyError):
289 self.antenna_names = self.antenna_numbers.astype(str).tolist()
290
291 # form up a grid which indexes time and baselines along the 'long'
292 # axis of the visdata array
293
294 1 2008998 2008998.0 0.7 tij_grid = np.array(map(lambda x: map(float, x.split("_")), unique_blts))
295 1 39 39.0 0.0 t_grid, ant_i_grid, ant_j_grid = tij_grid.T
296 # set the data sizes
297 1 4 4.0 0.0 try:
298 1 37 37.0 0.0 self.Nblts = uv['nblts']
299 1 10 10.0 0.0 if self.Nblts != len(t_grid):
300 warnings.warn('Nblts does not match the number of unique blts in the data')
301 except(KeyError):
302 self.Nblts = len(t_grid)
303 1 4 4.0 0.0 try:
304 1 12 12.0 0.0 self.Ntimes = uv['ntimes']
305 1 6 6.0 0.0 if self.Ntimes != len(times):
306 warnings.warn('Ntimes does not match the number of unique times in the data')
307 except(KeyError):
308 self.Ntimes = len(times)
309
310 1 6 6.0 0.0 self.time_array = t_grid
311 1 1172 1172.0 0.0 self.ant_1_array = ant_i_grid.astype(int)
312 1 2441 2441.0 0.0 self.ant_2_array = ant_j_grid.astype(int)
313
314 1 3470 3470.0 0.0 self.baseline_array = self.antnums_to_baseline(ant_i_grid.astype(int),
315 1 16148 16148.0 0.0 ant_j_grid.astype(int))
316 1 5 5.0 0.0 try:
317 1 18 18.0 0.0 self.Nbls = uv['nbls']
318 1 21750 21750.0 0.0 if self.Nbls != len(np.unique(self.baseline_array)):
319 warnings.warn('Nbls does not match the number of unique baselines in the data')
320 except(KeyError):
321 self.Nbls = len(np.unique(self.baseline_array))
322
323 # slot the data into a grid
324 1 12 12.0 0.0 self.data_array = np.zeros((self.Nblts, self.Nspws, self.Nfreqs,
325 1 31 31.0 0.0 self.Npols), dtype=np.complex64)
326 1 107702 107702.0 0.0 self.flag_array = np.ones(self.data_array.shape, dtype=np.bool)
327 1 1679 1679.0 0.0 self.uvw_array = np.zeros((self.Nblts, 3))
328 # NOTE: Using our lst calculator, which uses astropy,
329 # instead of aipy values which come from pyephem.
330 # The differences are of order 5 seconds.
331 1 3212863 3212863.0 1.1 self.set_lsts_from_time_array()
332 1 660518 660518.0 0.2 self.nsample_array = np.ones(self.data_array.shape, dtype=np.float)
333 1 61 61.0 0.0 self.freq_array = (np.arange(self.Nfreqs) * self.channel_width +
334 1 31 31.0 0.0 uv['sfreq'] * 1e9)
335 # Tile freq_array to shape (Nspws, Nfreqs).
336 # Currently does not actually support Nspws>1!
337 1 60 60.0 0.0 self.freq_array = np.tile(self.freq_array, (self.Nspws, 1))
338
339 # Temporary arrays to hold polarization axis, which will be collapsed
340 1 1289 1289.0 0.0 ra_pol_list = np.zeros((self.Nblts, self.Npols))
341 1 800 800.0 0.0 dec_pol_list = np.zeros((self.Nblts, self.Npols))
342 1 453 453.0 0.0 uvw_pol_list = np.zeros((self.Nblts, 3, self.Npols))
343 1 1562 1562.0 0.0 c_ns = const.c.to('m/ns').value
344 3 119371 39790.3 0.0 for pol, data in data_accumulator.iteritems():
345 2 217 108.5 0.0 pol_ind = self._pol_to_ind(pol)
346 924674 4384762 4.7 1.5 for ind, d in enumerate(data):
347 924672 9734240 10.5 3.4 blt = ["{1:.{0}f}".format(prec_t, d[1]).zfill(ndig_t),
348 924672 5654097 6.1 2.0 str(d[2]).zfill(ndig_ant), str(d[3]).zfill(ndig_ant)]
349 924672 4244202 4.6 1.5 blt = "_".join(blt)
350 924672 4044832 4.4 1.4 blt_index = reverse_inds[blt]
351
352 924672 12337336 13.3 4.3 self.data_array[blt_index, :, :, pol_ind] = d[4]
353 924672 9781916 10.6 3.4 self.flag_array[blt_index, :, :, pol_ind] = d[5]
354 924672 9661016 10.4 3.4 self.nsample_array[blt_index, :, :, pol_ind] = d[6]
355
356 # because there are uvws/ra/dec for each pol, and one pol may not
357 # have that visibility, we collapse along the polarization
358 # axis but avoid any missing visbilities
359 924672 8283971 9.0 2.9 uvw = d[0] * c_ns
360 924672 5304313 5.7 1.9 uvw.shape = (1, 3)
361 924672 7964836 8.6 2.8 uvw_pol_list[blt_index, :, pol_ind] = uvw
362 924672 6968192 7.5 2.4 ra_pol_list[blt_index, pol_ind] = d[7]
363 924672 6485602 7.0 2.3 dec_pol_list[blt_index, pol_ind] = d[8]
364
365 # Collapse pol axis for ra_list, dec_list, and uvw_list
366 1 655 655.0 0.0 ra_list = np.zeros(self.Nblts)
367 1 648 648.0 0.0 dec_list = np.zeros(self.Nblts)
368 462337 1846826 4.0 0.6 for blt_index in xrange(self.Nblts):
369 462336 11539119 25.0 4.0 test = ~np.all(self.flag_array[blt_index, :, :, :], axis=(0, 1))
370 462336 2873731 6.2 1.0 good_pol = np.where(test)[0]
371 462336 2003156 4.3 0.7 if len(good_pol) == 1:
372 # Only one good pol, use it
373 self.uvw_array[blt_index, :] = uvw_pol_list[blt_index, :, good_pol]
374 ra_list[blt_index] = ra_pol_list[blt_index, good_pol]
375 dec_list[blt_index] = dec_pol_list[blt_index, good_pol]
376 462336 1887218 4.1 0.7 elif len(good_pol) > 1:
377 # Multiple good pols, check for consistency. pyuvdata does not
378 # support pol-dependent uvw, ra, or dec.
379 216991 6986427 32.2 2.4 if np.any(np.diff(uvw_pol_list[blt_index, :, good_pol], axis=0)):
380 raise ValueError('uvw values are different by polarization.')
381 else:
382 216991 2026127 9.3 0.7 self.uvw_array[blt_index, :] = uvw_pol_list[blt_index, :, good_pol[0]]
383 216991 6339086 29.2 2.2 if np.any(np.diff(ra_pol_list[blt_index, good_pol])):
384 raise ValueError('ra values are different by polarization.')
385 else:
386 216991 1191884 5.5 0.4 ra_list[blt_index] = ra_pol_list[blt_index, good_pol[0]]
387 216991 6089454 28.1 2.1 if np.any(np.diff(dec_pol_list[blt_index, good_pol])):
388 raise ValueError('dec values are different by polarization.')
389 else:
390 216991 1162411 5.4 0.4 dec_list[blt_index] = dec_pol_list[blt_index, good_pol[0]]
391 else:
392 # No good pols for this blt. Fill with first one.
393 245345 1995663 8.1 0.7 self.uvw_array[blt_index, :] = uvw_pol_list[blt_index, :, 0]
394 245345 1160717 4.7 0.4 ra_list[blt_index] = ra_pol_list[blt_index, 0]
395 245345 1047888 4.3 0.4 dec_list[blt_index] = dec_pol_list[blt_index, 0]
396
397 # check if ra is constant throughout file; if it is,
398 # file is tracking if not, file is drift scanning
399 1 34598 34598.0 0.0 blt_good = np.where(~np.all(self.flag_array, axis=(1, 2, 3)))
400 1 2293 2293.0 0.0 if np.isclose(np.mean(np.diff(ra_list[blt_good])), 0.):
401 self.set_phased()
402 self.phase_center_ra = ra_list[0]
403 self.phase_center_dec = dec_list[0]
404 self.phase_center_epoch = uv['epoch']
405 else:
406 1 21 21.0 0.0 self.set_drift()
407 1 7 7.0 0.0 self.zenith_ra = ra_list
408 1 6 6.0 0.0 self.zenith_dec = dec_list
409 1 6 6.0 0.0 self.vis_units = 'UNCALIB' # assume no calibration
410
411 1 4 4.0 0.0 try:
412 1 1022 1022.0 0.0 self.set_telescope_params()
413 except ValueError, ve:
414 warnings.warn(str(ve))
415
416 # check if object has all required uv_properties set
417 1 4 4.0 0.0 if run_check:
418 1 127303 127303.0 0.0 self.check(run_check_acceptability=run_check_acceptability)
Total time: 76.8959 s
File: /users/alanman/optimize_pyuvdata/pyuvdata/pyuvdata/miriad.py
Function: write_miriad at line 419
Line # Hits Time Per Hit % Time Line Contents
==============================================================
419 @profile
420 def write_miriad(self, filepath, run_check=True, run_check_acceptability=True,
421 clobber=False, no_antnums=False):
422 """
423 Write the data to a miriad file.
424
425 Args:
426 filename: The miriad file directory to write to.
427 run_check: Option to check for the existence and proper shapes of
428 required parameters before writing the file. Default is True.
429 run_check_acceptability: Option to check acceptable range of the values of
430 required parameters before writing the file. Default is True.
431 clobber: Option to overwrite the filename if the file already exists.
432 Default is False.
433 no_antnums: Option to not write the antnums variable to the file.
434 Should only be used for testing purposes.
435 """
436 # check for multiple spws
437 1 8 8.0 0.0 if self.data_array.shape[1] > 1:
438 raise ValueError('write_miriad currently only handles single spw files.')
439
440 1 85 85.0 0.0 if os.path.exists(filepath):
441 if clobber:
442 print 'File exists: clobbering'
443 shutil.rmtree(filepath)
444 else:
445 raise ValueError('File exists: skipping')
446
447 1 16 16.0 0.0 freq_spacing = self.freq_array[0, 1:] - self.freq_array[0, :-1]
448 1 45 45.0 0.0 if not np.isclose(np.min(freq_spacing), np.max(freq_spacing),
449 1 114 114.0 0.0 rtol=self._freq_array.tols[0], atol=self._freq_array.tols[1]):
450 raise ValueError('The frequencies are not evenly spaced (probably '
451 'because of a select operation). The miriad format '
452 'does not support unevenly spaced frequencies.')
453 1 16 16.0 0.0 if not np.isclose(np.max(freq_spacing), self.channel_width,
454 1 81 81.0 0.0 rtol=self._freq_array.tols[0], atol=self._freq_array.tols[1]):
455 raise ValueError('The frequencies are separated by more than their '
456 'channel width (probably because of a select operation). '
457 'The miriad format does not support frequencies '
458 'that are spaced by more than their channel width.')
459
460 1 20189 20189.0 0.0 uv = a.miriad.UV(filepath, status='new')
461
462 # initialize header variables
463 1 31 31.0 0.0 uv._wrhd('obstype', 'mixed-auto-cross')
464 1 1498 1498.0 0.0 uv._wrhd('history', self.history + '\n')
465
466 # recognized miriad variables
467 1 8 8.0 0.0 uv.add_var('nchan', 'i')
468 1 13 13.0 0.0 uv['nchan'] = self.Nfreqs
469 1 3 3.0 0.0 uv.add_var('npol', 'i')
470 1 6 6.0 0.0 uv['npol'] = self.Npols
471 1 2 2.0 0.0 uv.add_var('nspect', 'i')
472 1 6 6.0 0.0 uv['nspect'] = self.Nspws
473 1 3 3.0 0.0 uv.add_var('inttime', 'd')
474 1 1219 1219.0 0.0 uv['inttime'] = self.integration_time
475 1 17 17.0 0.0 uv.add_var('sdf', 'd')
476 1 11 11.0 0.0 uv['sdf'] = self.channel_width / 1e9 # in GHz
477 1 3 3.0 0.0 uv.add_var('source', 'a')
478 1 6 6.0 0.0 uv['source'] = self.object_name
479 1 3 3.0 0.0 uv.add_var('telescop', 'a')
480 1 6 6.0 0.0 uv['telescop'] = self.telescope_name
481 1 3 3.0 0.0 uv.add_var('latitud', 'd')
482 1 89 89.0 0.0 uv['latitud'] = self.telescope_location_lat_lon_alt[0]
483 1 3 3.0 0.0 uv.add_var('longitu', 'd')
484 1 71 71.0 0.0 uv['longitu'] = self.telescope_location_lat_lon_alt[1]
485 1 3 3.0 0.0 uv.add_var('nants', 'i')
486
487 # Miriad has no way to keep track of antenna numbers, so the antenna
488 # numbers are simply the index for each antenna in any array that
489 # describes antenna attributes (e.g. antpos for the antenna_postions).
490 # Therefore on write, nants (which gives the size of the antpos array)
491 # needs to be increased to be the max value of antenna_numbers+1 and the
492 # antpos array needs to be inflated with zeros at locations where we
493 # don't have antenna information. These inflations need to be undone at
494 # read. If the file was written by pyuvdata, then the variable antnums
495 # will be present and we can use it, otherwise we need to test for zeros
496 # in the antpos array and/or antennas with no visibilities.
497 1 22 22.0 0.0 nants = np.max(self.antenna_numbers) + 1
498 1 6 6.0 0.0 uv['nants'] = nants
499 1 4 4.0 0.0 if self.antenna_positions is not None:
500 antpos = np.zeros((nants, 3), dtype=self.antenna_positions.dtype)
501 for ai, num in enumerate(self.antenna_numbers):
502 antpos[num, :] = self.antenna_positions[ai, :]
503 uv.add_var('antpos', 'd')
504 uv['antpos'] = antpos.T.flatten()
505
506 1 3 3.0 0.0 uv.add_var('sfreq', 'd')
507 1 7 7.0 0.0 uv['sfreq'] = self.freq_array[0, 0] / 1e9 # first spw; in GHz
508 1 3 3.0 0.0 if self.phase_type == 'phased':
509 uv.add_var('epoch', 'r')
510 uv['epoch'] = self.phase_center_epoch
511
512 # required pyuvdata variables that are not recognized miriad variables
513 1 3 3.0 0.0 uv.add_var('ntimes', 'i')
514 1 5 5.0 0.0 uv['ntimes'] = self.Ntimes
515 1 3 3.0 0.0 uv.add_var('nbls', 'i')
516 1 5 5.0 0.0 uv['nbls'] = self.Nbls
517 1 3 3.0 0.0 uv.add_var('nblts', 'i')
518 1 6 6.0 0.0 uv['nblts'] = self.Nblts
519 1 3 3.0 0.0 uv.add_var('visunits', 'a')
520 1 6 6.0 0.0 uv['visunits'] = self.vis_units
521 1 3 3.0 0.0 uv.add_var('instrume', 'a')
522 1 5 5.0 0.0 uv['instrume'] = self.instrument
523 1 3 3.0 0.0 uv.add_var('altitude', 'd')
524 1 70 70.0 0.0 uv['altitude'] = self.telescope_location_lat_lon_alt[2]
525
526 1 2 2.0 0.0 if not no_antnums:
527 # Add in the antenna_numbers so we have them if we read this file back in
528 # for some reason Miriad doesn't handle an array of integers properly,
529 # so convert to floats here and integers on read
530 1 3 3.0 0.0 uv.add_var('antnums', 'd')
531 1 12 12.0 0.0 uv['antnums'] = self.antenna_numbers.astype(np.float64)
532
533 # antenna names is a foreign concept in miriad but required in other formats.
534 # This is a horrible hack to save antenna_names array. Miriad can't handle arrays
535 # of strings or a long enough single string to put them all into one string
536 # so we convert them into hex values and then into floats and convert
537 # back to strings on read
538 129 622 4.8 0.0 ant_name_flt = np.array([int(elem.encode("hex"), 16) for elem in self.antenna_names]).astype(np.float64)
539 1 3 3.0 0.0 uv.add_var('antnames', 'd')
540 1 8 8.0 0.0 uv['antnames'] = ant_name_flt
541
542 # variables that can get updated with every visibility
543 1 3 3.0 0.0 uv.add_var('pol', 'i')
544 1 4 4.0 0.0 uv.add_var('lst', 'd')
545 1 3 3.0 0.0 uv.add_var('cnt', 'd')
546 1 2 2.0 0.0 uv.add_var('ra', 'd')
547 1 3 3.0 0.0 uv.add_var('dec', 'd')
548
549 # write data
550 1 858 858.0 0.0 c_ns = const.c.to('m/ns').value
551 462337 1317473 2.8 1.7 for viscnt, blt in enumerate(self.data_array):
552 462336 5310432 11.5 6.9 uvw = (self.uvw_array[viscnt] / c_ns).astype(np.double) # NOTE issue 50 on conjugation
553 462336 1944312 4.2 2.5 t = self.time_array[viscnt]
554 462336 1568615 3.4 2.0 i = self.ant_1_array[viscnt]
555 462336 1490438 3.2 1.9 j = self.ant_2_array[viscnt]
556
557 462336 2884012 6.2 3.8 uv['lst'] = self.lst_array[viscnt]
558 462336 1600097 3.5 2.1 if self.phase_type == 'phased':
559 uv['ra'] = self.phase_center_ra
560 uv['dec'] = self.phase_center_dec
561 462336 1443938 3.1 1.9 elif self.phase_type == 'drift':
562 462336 2515366 5.4 3.3 uv['ra'] = self.zenith_ra[viscnt]
563 462336 2590451 5.6 3.4 uv['dec'] = self.zenith_dec[viscnt]
564 else:
565 raise ValueError('The phasing type of the data is unknown. '
566 'Set the phase_type to "drift" or "phased" to '
567 'reflect the phasing status of the data')
568
569 # NOTE only writing spw 0, not supporting multiple spws for write
570 1387008 6176944 4.5 8.0 for polcnt, pol in enumerate(self.polarization_array):
571 924672 7376416 8.0 9.6 uv['pol'] = pol.astype(np.int)
572 924672 11173991 12.1 14.5 uv['cnt'] = self.nsample_array[viscnt, 0, :, polcnt].astype(np.double)
573
574 924672 4099620 4.4 5.3 data = self.data_array[viscnt, 0, :, polcnt]
575 924672 3612908 3.9 4.7 flags = self.flag_array[viscnt, 0, :, polcnt]
576 924672 2280073 2.5 3.0 if i > j:
577 i, j, data = j, i, np.conjugate(data)
578 924672 2187587 2.4 2.8 preamble = (uvw, t, (i, j))
579
580 924672 17172875 18.6 22.3 uv.write(preamble, data, flags)
581 1 2 2.0 0.0 if run_check:
582 """Check for acceptable units."""
583 1 125145 125145.0 0.2 self.check(run_check_acceptability=run_check_acceptability)
Total time: 3.82707 s
File: /users/alanman/optimize_pyuvdata/pyuvdata/pyuvdata/uvdata.py
Function: set_lsts_from_time_array at line 457
Line # Hits Time Per Hit % Time Line Contents
==============================================================
457 @profile
458 def set_lsts_from_time_array(self):
459 """Set the lst_array based from the time_array."""
460 2 3 1.5 0.0 lsts = []
461 2 9 4.5 0.0 curtime = self.time_array[0]
462 2 1093 546.5 0.0 self.lst_array = np.zeros(self.Nblts)
463 2 213 106.5 0.0 latitude, longitude, altitude = self.telescope_location_lat_lon_alt_degrees
464 114 18347 160.9 0.5 for ind, jd in enumerate(np.unique(self.time_array)):
465 112 151277 1350.7 4.0 t = Time(jd, format='jd', location=(longitude, latitude))
466 112 3656128 32644.0 95.5 self.lst_array[np.where(np.isclose(jd, self.time_array, atol=1e-6, rtol=1e-12))] = t.sidereal_time('apparent').radian
Total time: 0.00046 s
File: /users/alanman/optimize_pyuvdata/pyuvdata/pyuvdata/uvdata.py
Function: juldate2ephem at line 468
Line # Hits Time Per Hit % Time Line Contents
==============================================================
468 @profile
469 def juldate2ephem(self, num):
470 """
471 Convert Julian date to ephem date, measured from noon, Dec. 31, 1899.
472
473 Args:
474 num: Julian date
475
476 Returns:
477 ephem date, measured from noon, Dec. 31, 1899.
478 """
479 227 460 2.0 100.0 return ephem.date(num - 2415020.)
Total time: 16.0526 s
File: /users/alanman/optimize_pyuvdata/pyuvdata/pyuvdata/uvdata.py
Function: unphase_to_drift at line 481
Line # Hits Time Per Hit % Time Line Contents
==============================================================
481 @profile
482 def unphase_to_drift(self):
483 """Convert from a phased dataset to a drift dataset."""
484 1 4 4.0 0.0 if self.phase_type == 'phased':
485 1 2 2.0 0.0 pass
486 elif self.phase_type == 'drift':
487 raise ValueError('The data is already drift scanning; can only ' +
488 'unphase phased data.')
489 else:
490 raise ValueError('The phasing type of the data is unknown. '
491 'Set the phase_type to drift or phased to '
492 'reflect the phasing status of the data')
493
494 1 69 69.0 0.0 latitude, longitude, altitude = self.telescope_location_lat_lon_alt
495
496 1 5 5.0 0.0 obs = ephem.Observer()
497 # obs inits with default values for parameters -- be sure to replace them
498 1 4 4.0 0.0 obs.lat = latitude
499 1 1 1.0 0.0 obs.lon = longitude
500
501 1 3 3.0 0.0 phase_center = ephem.FixedBody()
502 1 5 5.0 0.0 epoch = (self.phase_center_epoch - 2000.) * 365.2422 + ephem.J2000 # convert years to ephemtime
503 1 2 2.0 0.0 phase_center._epoch = epoch
504 1 3 3.0 0.0 phase_center._ra = self.phase_center_ra
505 1 3 3.0 0.0 phase_center._dec = self.phase_center_dec
506
507 1 190 190.0 0.0 self.zenith_ra = np.zeros_like(self.time_array)
508 1 332 332.0 0.0 self.zenith_dec = np.zeros_like(self.time_array)
509
510 # apply -w phasor
511 w_lambda = (self.uvw_array[:, 2].reshape(self.Nblts, 1) /
512 1 570278 570278.0 3.6 const.c.to('m/s').value * self.freq_array.reshape(1, self.Nfreqs))
513 1 11725581 11725581.0 73.0 phs = np.exp(-1j * 2 * np.pi * (-1) * w_lambda[:, None, :, None])
514 1 3695210 3695210.0 23.0 self.data_array *= phs
515
516 1 8690 8690.0 0.1 unique_times = np.unique(self.time_array)
517 57 144 2.5 0.0 for jd in unique_times:
518 56 17743 316.8 0.1 inds = np.where(self.time_array == jd)[0]
519 56 829 14.8 0.0 obs.date, obs.epoch = self.juldate2ephem(jd), self.juldate2ephem(jd)
520 56 120 2.1 0.0 phase_center.compute(obs)
521 56 2698 48.2 0.0 phase_center_ra, phase_center_dec = phase_center.a_ra, phase_center.a_dec
522 56 106 1.9 0.0 zenith_ra = obs.sidereal_time()
523 56 79 1.4 0.0 zenith_dec = latitude
524 56 1388 24.8 0.0 self.zenith_ra[inds] = zenith_ra
525 56 1278 22.8 0.0 self.zenith_dec[inds] = zenith_dec
526
527 # generate rotation matrices
528 56 1780 31.8 0.0 m0 = uvutils.top2eq_m(0., phase_center_dec)
529 56 1357 24.2 0.0 m1 = uvutils.eq2top_m(phase_center_ra - zenith_ra, zenith_dec)
530
531 # rotate and write uvws
532 56 8072 144.1 0.1 uvw = self.uvw_array[inds, :]
533 56 4558 81.4 0.0 uvw = np.dot(m0, uvw.T).T
534 56 3378 60.3 0.0 uvw = np.dot(m1, uvw.T).T
535 56 8619 153.9 0.1 self.uvw_array[inds, :] = uvw
536
537 # remove phase center
538 1 6 6.0 0.0 self.phase_center_ra = None
539 1 4 4.0 0.0 self.phase_center_dec = None
540 1 3 3.0 0.0 self.phase_center_epoch = None
541 1 9 9.0 0.0 self.set_drift()
Total time: 15.9646 s
File: /users/alanman/optimize_pyuvdata/pyuvdata/pyuvdata/uvdata.py
Function: phase_to_time at line 543
Line # Hits Time Per Hit % Time Line Contents
==============================================================
543 @profile
544 def phase_to_time(self, time):
545 """
546 Phase a drift scan dataset to the ra/dec of zenith at a particular time.
547
548 Args:
549 time: The time to phase to.
550 """
551 1 2 2.0 0.0 if self.phase_type == 'drift':
552 1 1 1.0 0.0 pass
553 elif self.phase_type == 'phased':
554 raise ValueError('The data is already phased; can only phase ' +
555 'drift scanning data.')
556 else:
557 raise ValueError('The phasing type of the data is unknown. '
558 'Set the phase_type to drift or phased to '
559 'reflect the phasing status of the data')
560
561 1 60 60.0 0.0 obs = ephem.Observer()
562 # obs inits with default values for parameters -- be sure to replace them
563 1 119 119.0 0.0 latitude, longitude, altitude = self.telescope_location_lat_lon_alt
564 1 4 4.0 0.0 obs.lat = latitude
565 1 1 1.0 0.0 obs.lon = longitude
566
567 1 35 35.0 0.0 obs.date, obs.epoch = self.juldate2ephem(time), self.juldate2ephem(time)
568
569 1 85 85.0 0.0 ra = obs.sidereal_time()
570 1 1 1.0 0.0 dec = latitude
571 1 7 7.0 0.0 epoch = self.juldate2ephem(time)
572 1 15964312 15964312.0 100.0 self.phase(ra, dec, epoch)
Total time: 15.9553 s
File: /users/alanman/optimize_pyuvdata/pyuvdata/pyuvdata/uvdata.py
Function: phase at line 574
Line # Hits Time Per Hit % Time Line Contents
==============================================================
574 @profile
575 def phase(self, ra, dec, epoch):
576 """
577 Phase a drift scan dataset to a single ra/dec at a particular epoch.
578
579 Will not phase already phased data.
580
581 Args:
582 ra: The ra to phase to in radians.
583 dec: The dec to phase to in radians.
584 epoch: The epoch to use for phasing. Should be an ephem date,
585 measured from noon Dec. 31, 1899.
586 """
587 1 3 3.0 0.0 if self.phase_type == 'drift':
588 1 2 2.0 0.0 pass
589 elif self.phase_type == 'phased':
590 raise ValueError('The data is already phased; can only phase '
591 'drift scan data. Use unphase_to_drift to '
592 'convert to a drift scan.')
593 else:
594 raise ValueError('The phasing type of the data is unknown. '
595 'Set the phase_type to "drift" or "phased" to '
596 'reflect the phasing status of the data')
597
598 1 2 2.0 0.0 obs = ephem.Observer()
599 # obs inits with default values for parameters -- be sure to replace them
600 1 68 68.0 0.0 latitude, longitude, altitude = self.telescope_location_lat_lon_alt
601 1 2 2.0 0.0 obs.lat = latitude
602 1 1 1.0 0.0 obs.lon = longitude
603
604 # create a pyephem object for the phasing position
605 1 2 2.0 0.0 precess_pos = ephem.FixedBody()
606 1 2 2.0 0.0 precess_pos._ra = ra
607 1 1 1.0 0.0 precess_pos._dec = dec
608 1 1 1.0 0.0 precess_pos._epoch = epoch
609
610 # calculate RA/DEC in J2000 and write to object
611 1 2 2.0 0.0 obs.date, obs.epoch = ephem.J2000, ephem.J2000
612 1 3 3.0 0.0 precess_pos.compute(obs)
613
614 1 142 142.0 0.0 self.phase_center_ra = precess_pos.a_ra + 0.0 # force to be a float not ephem.Angle
615 1 5 5.0 0.0 self.phase_center_dec = precess_pos.a_dec + 0.0 # force to be a float not ephem.Angle
616 # explicitly set epoch to J2000
617 1 4 4.0 0.0 self.phase_center_epoch = 2000.0
618
619 1 11379 11379.0 0.1 unique_times, unique_inds = np.unique(self.time_array, return_index=True)
620 57 157 2.8 0.0 for ind, jd in enumerate(unique_times):
621 56 45881 819.3 0.3 inds = np.where(self.time_array == jd)[0]
622 56 308 5.5 0.0 lst = self.lst_array[unique_inds[ind]]
623 # calculate ra/dec of phase center in current epoch
624 56 784 14.0 0.0 obs.date, obs.epoch = self.juldate2ephem(jd), self.juldate2ephem(jd)
625 56 117 2.1 0.0 precess_pos.compute(obs)
626 56 2746 49.0 0.0 ra, dec = precess_pos.a_ra, precess_pos.a_dec
627
628 # generate rotation matrices
629 56 1842 32.9 0.0 m0 = uvutils.top2eq_m(lst - obs.sidereal_time(), latitude)
630 56 1439 25.7 0.0 m1 = uvutils.eq2top_m(lst - ra, dec)
631
632 # rotate and write uvws
633 56 10616 189.6 0.1 uvw = self.uvw_array[inds, :]
634 56 3756 67.1 0.0 uvw = np.dot(m0, uvw.T).T
635 56 3291 58.8 0.0 uvw = np.dot(m1, uvw.T).T
636 56 8710 155.5 0.1 self.uvw_array[inds, :] = uvw
637
638 # calculate data and apply phasor
639 w_lambda = (self.uvw_array[:, 2].reshape(self.Nblts, 1) /
640 1 491238 491238.0 3.1 const.c.to('m/s').value * self.freq_array.reshape(1, self.Nfreqs))
641 1 11694982 11694982.0 73.3 phs = np.exp(-1j * 2 * np.pi * w_lambda[:, None, :, None])
642 1 3677817 3677817.0 23.1 self.data_array *= phs
643
644 1 4 4.0 0.0 del(obs)
645 1 14 14.0 0.0 self.set_phased()
Total time: 14.4262 s
File: /users/alanman/optimize_pyuvdata/pyuvdata/pyuvdata/uvfits.py
Function: read_uvfits at line 45
Line # Hits Time Per Hit % Time Line Contents
==============================================================
45 @profile
46 def read_uvfits(self, filename, run_check=True, run_check_acceptability=True):
47 """
48 Read in data from a uvfits file.
49
50 Args:
51 filename: The uvfits file to read from.
52 run_check: Option to check for the existence and proper shapes of
53 required parameters after reading in the file. Default is True.
54 run_check_acceptability: Option to check acceptable range of the values of
55 required parameters after reading in the file. Default is True.
56 """
57 1 12587 12587.0 0.1 F = fits.open(filename)
58 1 10 10.0 0.0 D = F[0] # assumes the visibilities are in the primary hdu
59 1 2860 2860.0 0.0 hdr = D.header.copy()
60 1 191 191.0 0.0 hdunames = self._indexhdus(F) # find the rest of the tables
61
62 # astropy.io fits reader scales date according to relevant PZER0 (?)
63 1 3456621 3456621.0 24.0 time0_array = D.data['DATE']
64 1 4 4.0 0.0 try:
65 # uvfits standard is to have 2 DATE parameters, both floats:
66 # DATE (full day) and _DATE (fractional day)
67 1 165 165.0 0.0 time1_array = D.data['_DATE']
68 self.time_array = (time0_array.astype(np.double) +
69 time1_array.astype(np.double))
70 1 4 4.0 0.0 except(KeyError):
71 # cotter uvfits files have one DATE that is a double
72 1 11 11.0 0.0 self.time_array = time0_array
73 1 29 29.0 0.0 if np.finfo(time0_array[0]).precision < 5:
74 raise ValueError('JDs in this file are not precise to '
75 'better than a second.')
76 1 7 7.0 0.0 if (np.finfo(time0_array[0]).precision > 5 and
77 1 7 7.0 0.0 np.finfo(time0_array[0]).precision < 8):
78 warnings.warn('The JDs in this file have sub-second '
79 'precision, but not sub-millisecond. '
80 'Use with caution.')
81
82 1 9420 9420.0 0.1 self.Ntimes = len(np.unique(self.time_array))
83
84 # if antenna arrays are present, use them. otherwise use baseline array
85 1 3 3.0 0.0 try:
86 # Note: uvfits antennas are 1 indexed,
87 # need to subtract one to get to 0-indexed
88 1 127985 127985.0 0.9 self.ant_1_array = np.int32(D.data.field('ANTENNA1')) - 1
89 1 27542 27542.0 0.2 self.ant_2_array = np.int32(D.data.field('ANTENNA2')) - 1
90 1 26001 26001.0 0.2 subarray = np.int32(D.data.field('SUBARRAY')) - 1
91 # error on files with multiple subarrays
92 1 40825 40825.0 0.3 if len(set(subarray)) > 1:
93 raise ValueError('This file appears to have multiple subarray '
94 'values; only files with one subarray are '
95 'supported.')
96 except(KeyError):
97 # cannot set this to be the baseline array because it uses the
98 # 256 convention, not our 2048 convention
99 bl_input_array = np.int64(D.data.field('BASELINE'))
100
101 # get antenna arrays based on uvfits baseline array
102 self.ant_1_array, self.ant_2_array = \
103 self.baseline_to_antnums(bl_input_array)
104 # check for multi source files
105 1 3 3.0 0.0 try:
106 1 170 170.0 0.0 source = D.data.field('SOURCE')
107 if len(set(source)) > 1:
108 raise ValueError('This file has multiple sources. Only single '
109 'source observations are supported.')
110 1 3 3.0 0.0 except:
111 1 2 2.0 0.0 pass
112
113 # get self.baseline_array using our convention
114 self.baseline_array = \
115 1 9 9.0 0.0 self.antnums_to_baseline(self.ant_1_array,
116 1 3254 3254.0 0.0 self.ant_2_array)
117 1 23001 23001.0 0.2 self.Nbls = len(np.unique(self.baseline_array))
118
119 # initialize internal variables based on the antenna table
120 1 85620 85620.0 0.6 self.Nants_data = int(len(np.unique(self.ant_1_array.tolist() + self.ant_2_array.tolist())))
121
122 1 16 16.0 0.0 self.set_phased()
123 # check if we have an spw dimension
124 1 340 340.0 0.0 if hdr.pop('NAXIS') == 7:
125 1 64 64.0 0.0 if hdr['NAXIS5'] > 1:
126 raise ValueError('Sorry. Files with more than one spectral' +
127 'window (spw) are not yet supported. A ' +
128 'great project for the interested student!')
129 1 338 338.0 0.0 self.data_array = (D.data.field('DATA')[:, 0, 0, :, :, :, 0] +
130 1 8472373 8472373.0 58.7 1j * D.data.field('DATA')[:, 0, 0, :, :, :, 1])
131 1 546361 546361.0 3.8 self.flag_array = (D.data.field('DATA')[:, 0, 0, :, :, :, 2] <= 0)
132 1 781385 781385.0 5.4 self.nsample_array = np.abs(D.data.field('DATA')[:, 0, 0, :, :, :, 2])
133 1 345 345.0 0.0 self.Nspws = hdr.pop('NAXIS5')
134 1 9 9.0 0.0 assert(self.Nspws == self.data_array.shape[1])
135
136 # the axis number for phase center depends on if the spw exists
137 # subtract 1 to be zero-indexed
138 1 614 614.0 0.0 self.spw_array = np.int32(self._gethduaxis(D, 5)) - 1
139
140 1 294 294.0 0.0 self.phase_center_ra_degrees = np.float(hdr.pop('CRVAL6'))
141 1 282 282.0 0.0 self.phase_center_dec_degrees = np.float(hdr.pop('CRVAL7'))
142 else:
143 # in many uvfits files the spw axis is left out,
144 # here we put it back in so the dimensionality stays the same
145 self.data_array = (D.data.field('DATA')[:, 0, 0, :, :, 0] +
146 1j * D.data.field('DATA')[:, 0, 0, :, :, 1])
147 self.data_array = self.data_array[:, np.newaxis, :, :]
148 self.flag_array = (D.data.field('DATA')[:, 0, 0, :, :, 2] <= 0)
149 self.flag_array = self.flag_array[:, np.newaxis, :, :]
150 self.nsample_array = np.abs(D.data.field('DATA')[:, 0, 0, :, :, 2])
151 self.nsample_array = (self.nsample_array[:, np.newaxis, :, :])
152
153 # the axis number for phase center depends on if the spw exists
154 self.Nspws = 1
155 self.spw_array = np.array([0])
156
157 self.phase_center_ra_degrees = np.float(hdr.pop('CRVAL5'))
158 self.phase_center_dec_degrees = np.float(hdr.pop('CRVAL6'))
159
160 # get shapes
161 1 191 191.0 0.0 self.Nfreqs = hdr.pop('NAXIS4')
162 1 6 6.0 0.0 if self.data_array.shape[2] != self.Nfreqs:
163 warnings.warn('Nfreqs does not match the number of frequencies in the data')
164 1 188 188.0 0.0 self.Npols = hdr.pop('NAXIS3')
165 1 6 6.0 0.0 if self.data_array.shape[3] != self.Npols:
166 warnings.warn('npols={npols} but found {l} pols in data file'.format(
167 npols=self.Npols, l=len(self.polarization_array)))
168 1 184 184.0 0.0 self.Nblts = hdr.pop('GCOUNT')
169 1 5 5.0 0.0 if self.data_array.shape[0] != self.Nblts:
170 warnings.warn('Nblts does not match the number of unique blts in the data')
171
172 # read baseline vectors in units of seconds, return in meters
173 1 304 304.0 0.0 self.uvw_array = (np.array(np.stack((D.data.field('UU'),
174 1 235 235.0 0.0 D.data.field('VV'),
175 1 34285 34285.0 0.2 D.data.field('WW')))) *
176 1 10224 10224.0 0.1 const.c.to('m/s').value).T
177
178 1 618 618.0 0.0 self.freq_array = self._gethduaxis(D, 4)
179 1 292 292.0 0.0 self.channel_width = hdr.pop('CDELT4')
180
181 1 3 3.0 0.0 try:
182 1 308 308.0 0.0 self.integration_time = float(D.data.field('INTTIM')[0])
183 except:
184 if self.Ntimes > 1:
185 self.integration_time = \
186 float(np.diff(np.sort(list(set(self.time_array))))[0]) * 86400
187 else:
188 raise ValueError('integration time not specified and only '
189 'one time present')
190
191 1 11 11.0 0.0 self.freq_array.shape = (self.Nspws,) + self.freq_array.shape
192
193 1 498 498.0 0.0 self.polarization_array = np.int32(self._gethduaxis(D, 3))
194
195 # other info -- not required but frequently used
196 1 270 270.0 0.0 self.object_name = hdr.pop('OBJECT', None)
197 1 256 256.0 0.0 self.telescope_name = hdr.pop('TELESCOP', None)
198 1 245 245.0 0.0 self.instrument = hdr.pop('INSTRUME', None)
199 1 261 261.0 0.0 latitude_degrees = hdr.pop('LAT', None)
200 1 258 258.0 0.0 longitude_degrees = hdr.pop('LON', None)
201 1 252 252.0 0.0 altitude = hdr.pop('ALT', None)
202 1 319 319.0 0.0 self.history = str(hdr.get('HISTORY', ''))
203 1 8 8.0 0.0 if self.pyuvdata_version_str not in self.history.replace('\n', ''):
204 self.history += self.pyuvdata_version_str
205 6 525 87.5 0.0 while 'HISTORY' in hdr.keys():
206 5 473 94.6 0.0 hdr.remove('HISTORY')
207 1 245 245.0 0.0 self.vis_units = hdr.pop('BUNIT', 'UNCALIB')
208 1 254 254.0 0.0 self.phase_center_epoch = hdr.pop('EPOCH', None)
209
210 # remove standard FITS header items that are still around
211 1 3 3.0 0.0 std_fits_substrings = ['SIMPLE', 'BITPIX', 'EXTEND', 'BLOCKED',
212 1 3 3.0 0.0 'GROUPS', 'PCOUNT', 'BSCALE', 'BZERO', 'NAXIS',
213 1 3 3.0 0.0 'PTYPE', 'PSCAL', 'PZERO', 'CTYPE', 'CRVAL',
214 1 3 3.0 0.0 'CRPIX', 'CDELT', 'CROTA']
215 57 230 4.0 0.0 for key in hdr.keys():
216 1008 2689 2.7 0.0 for sub in std_fits_substrings:
217 952 2834 3.0 0.0 if key.find(sub) > -1:
218 55 3645 66.3 0.0 hdr.remove(key)
219
220 # find all the remaining header items and keep them as extra_keywords
221 2 10 5.0 0.0 for key in hdr:
222 1 3 3.0 0.0 if key == '':
223 continue
224 1 3 3.0 0.0 if key == 'COMMENT':
225 self.extra_keywords[key] = str(hdr.get(key))
226 else:
227 1 147 147.0 0.0 self.extra_keywords[key] = hdr.get(key)
228
229 # READ the antenna table
230 1 12 12.0 0.0 ant_hdu = F[hdunames['AIPS AN']]
231
232 # stuff in columns
233 1 15194 15194.0 0.1 ant_names = ant_hdu.data.field('ANNAME').tolist()
234 1 7 7.0 0.0 self.antenna_names = []
235 129 350 2.7 0.0 for name in ant_names:
236 128 509 4.0 0.0 self.antenna_names.append(name.replace('\x00!', ''))
237
238 # subtract one to get to 0-indexed values rather than 1-indexed values
239 1 278 278.0 0.0 self.antenna_numbers = ant_hdu.data.field('NOSTA') - 1
240
241 1 7 7.0 0.0 self.Nants_telescope = len(self.antenna_numbers)
242
243 # stuff in the header
244 1 4 4.0 0.0 if self.telescope_name is None:
245 self.telescope_name = ant_hdu.header['ARRNAM']
246
247 1 3 3.0 0.0 try:
248 1 69 69.0 0.0 xyz_telescope_frame = ant_hdu.header['FRAME']
249 except:
250 warnings.warn('Required Antenna frame keyword not set, '
251 'setting to ????')
252 xyz_telescope_frame = '????'
253
254 # VLA incorrectly sets ARRAYX/ARRAYY/ARRAYZ to 0, and puts array center
255 # in the antenna positions themselves
256
257 1 168 168.0 0.0 if (np.isclose(ant_hdu.header['ARRAYX'], 0) and
258 np.isclose(ant_hdu.header['ARRAYY'], 0) and
259 np.isclose(ant_hdu.header['ARRAYZ'], 0)):
260 x_telescope = np.mean(ant_hdu.data['STABXYZ'][:, 0])
261 y_telescope = np.mean(ant_hdu.data['STABXYZ'][:, 1])
262 z_telescope = np.mean(ant_hdu.data['STABXYZ'][:, 2])
263 self.antenna_positions = (ant_hdu.data.field('STABXYZ') -
264 np.array([x_telescope,
265 y_telescope,
266 z_telescope]))
267
268 else:
269 1 61 61.0 0.0 x_telescope = ant_hdu.header['ARRAYX']
270 1 55 55.0 0.0 y_telescope = ant_hdu.header['ARRAYY']
271 1 53 53.0 0.0 z_telescope = ant_hdu.header['ARRAYZ']
272 1 248 248.0 0.0 self.antenna_positions = ant_hdu.data.field('STABXYZ')
273
274 1 3 3.0 0.0 if xyz_telescope_frame == 'ITRF':
275 1 11 11.0 0.0 self.telescope_location = np.array([x_telescope, y_telescope, z_telescope])
276 else:
277 if latitude_degrees is not None and longitude_degrees is not None and altitude is not None:
278 self.telescope_location_lat_lon_alt_degrees = (latitude_degrees, longitude_degrees, altitude)
279
280 1 63 63.0 0.0 self.gst0 = ant_hdu.header['GSTIA0']
281 1 58 58.0 0.0 self.rdate = ant_hdu.header['RDATE']
282 1 56 56.0 0.0 self.earth_omega = ant_hdu.header['DEGPDY']
283 1 56 56.0 0.0 self.dut1 = ant_hdu.header['UT1UTC']
284 1 3 3.0 0.0 try:
285 1 24 24.0 0.0 self.timesys = ant_hdu.header['TIMESYS']
286 1 3 3.0 0.0 except(KeyError):
287 # CASA misspells this one
288 1 58 58.0 0.0 self.timesys = ant_hdu.header['TIMSYS']
289
290 1 3 3.0 0.0 del(D)
291
292 1 3 3.0 0.0 try:
293 1 942 942.0 0.0 self.set_telescope_params()
294 except ValueError, ve:
295 warnings.warn(str(ve))
296
297 1 614854 614854.0 4.3 self.set_lsts_from_time_array()
298
299 # check if object has all required UVParameters set
300 1 3 3.0 0.0 if run_check:
301 1 113947 113947.0 0.8 self.check(run_check_acceptability=run_check_acceptability)
Total time: 32.6351 s
File: /users/alanman/optimize_pyuvdata/pyuvdata/pyuvdata/uvfits.py
Function: write_uvfits at line 303
Line # Hits Time Per Hit % Time Line Contents
==============================================================
303 @profile
304 def write_uvfits(self, filename, spoof_nonessential=False,
305 force_phase=False, run_check=True, run_check_acceptability=True):
306 """
307 Write the data to a uvfits file.
308
309 Args:
310 filename: The uvfits file to write to.
311 spoof_nonessential: Option to spoof the values of optional
312 UVParameters that are not set but are required for uvfits files.
313 Default is False.
314 force_phase: Option to automatically phase drift scan data to
315 zenith of the first timestamp. Default is False.
316 run_check: Option to check for the existence and proper shapes of
317 required parameters before writing the file. Default is True.
318 run_check_acceptability: Option to check acceptable range of the values of
319 required parameters before writing the file. Default is True.
320 """
321 1 5 5.0 0.0 if run_check:
322 1 124998 124998.0 0.4 self.check(run_check_acceptability=run_check_acceptability)
323
324 1 14 14.0 0.0 if self.phase_type == 'phased':
325 pass
326 1 5 5.0 0.0 elif self.phase_type == 'drift':
327 1 3 3.0 0.0 if force_phase:
328 1 53 53.0 0.0 print('The data are in drift mode and do not have a '
329 'defined phase center. Phasing to zenith of the first '
330 'timestamp.')
331 1 15964681 15964681.0 48.9 self.phase_to_time(self.time_array[0])
332 else:
333 raise ValueError('The data are in drift mode. ' +
334 'Set force_phase to true to phase the data ' +
335 'to zenith of the first timestamp before ' +
336 'writing a uvfits file.')
337 else:
338 raise ValueError('The phasing type of the data is unknown. '
339 'Set the phase_type to drift or phased to '
340 'reflect the phasing status of the data')
341
342 1 62 62.0 0.0 freq_spacing = self.freq_array[0, 1:] - self.freq_array[0, :-1]
343 1 64 64.0 0.0 if not np.isclose(np.min(freq_spacing), np.max(freq_spacing),
344 1 150 150.0 0.0 rtol=self._freq_array.tols[0], atol=self._freq_array.tols[1]):
345 raise ValueError('The frequencies are not evenly spaced (probably '
346 'because of a select operation). The uvfits format '
347 'does not support unevenly spaced frequencies.')
348 1 18 18.0 0.0 if not np.isclose(np.max(freq_spacing), self.channel_width,
349 1 83 83.0 0.0 rtol=self._freq_array.tols[0], atol=self._freq_array.tols[1]):
350 raise ValueError('The frequencies are separated by more than their '
351 'channel width (probably because of a select operation). '
352 'The uvfits format does not support frequencies '
353 'that are spaced by more than their channel width.')
354 1 12 12.0 0.0 if self.Npols > 2:
355 pol_spacing = self.polarization_array[1:] - self.polarization_array[:-1]
356 if np.min(pol_spacing) < np.max(pol_spacing):
357 raise ValueError('The polarization values are not evenly spaced (probably '
358 'because of a select operation). The uvfits format '
359 'does not support unevenly spaced polarizations.')
360
361 11 8640 785.5 0.0 for p in self.extra():
362 10 36 3.6 0.0 param = getattr(self, p)
363 10 40 4.0 0.0 if param.name in self.uvfits_required_extra:
364 6 23 3.8 0.0 if param.value is None:
365 6 19 3.2 0.0 if spoof_nonessential:
366 # spoof extra keywords required for uvfits
367 6 25 4.2 0.0 if isinstance(param, uvp.AntPositionParameter):
368 1 19 19.0 0.0 param.apply_spoof(self, 'Nants_telescope')
369 else:
370 5 28 5.6 0.0 param.apply_spoof()
371 6 25 4.2 0.0 setattr(self, p, param)
372 else:
373 raise ValueError('Required attribute {attribute} '
374 'for uvfits not defined. Define or '
375 'set spoof_nonessential to True to '
376 'spoof this attribute.'
377 .format(attribute=p))
378
379 1 5 5.0 0.0 weights_array = self.nsample_array * \
380 1 2142990 2142990.0 6.6 np.where(self.flag_array, -1, 1)
381 1 39 39.0 0.0 data_array = self.data_array[:, np.newaxis, np.newaxis, :, :, :, np.newaxis]
382 1 4 4.0 0.0 weights_array = weights_array[:, np.newaxis, np.newaxis, :, :, :,
383 1 5 5.0 0.0 np.newaxis]
384 # uvfits_array_data shape will be (Nblts,1,1,[Nspws],Nfreqs,Npols,3)
385 1 6 6.0 0.0 uvfits_array_data = np.concatenate([data_array.real,
386 1 4 4.0 0.0 data_array.imag,
387 1 4178199 4178199.0 12.8 weights_array], axis=6)
388
389 1 5501 5501.0 0.0 uvw_array_sec = self.uvw_array / const.c.to('m/s').value
390 # jd_midnight = np.floor(self.time_array[0] - 0.5) + 0.5
391 1 39 39.0 0.0 tzero = np.float32(self.time_array[0])
392
393 # uvfits convention is that time_array + relevant PZERO = actual JD
394 # We are setting PZERO4 = float32(first time of observation)
395 1 1438 1438.0 0.0 time_array = np.float32(self.time_array - np.float64(tzero))
396
397 1 198 198.0 0.0 int_time_array = (np.zeros_like((time_array), dtype=np.float) +
398 1 456 456.0 0.0 self.integration_time)
399
400 1 10 10.0 0.0 baselines_use = self.antnums_to_baseline(self.ant_1_array,
401 1 5 5.0 0.0 self.ant_2_array,
402 1 3807 3807.0 0.0 attempt256=True)
403 # Set up dictionaries for populating hdu
404 # Note that uvfits antenna arrays are 1-indexed so we add 1
405 # to our 0-indexed arrays
406 1 12 12.0 0.0 group_parameter_dict = {'UU ': uvw_array_sec[:, 0],
407 1 4 4.0 0.0 'VV ': uvw_array_sec[:, 1],
408 1 4 4.0 0.0 'WW ': uvw_array_sec[:, 2],
409 1 4 4.0 0.0 'DATE ': time_array,
410 1 4 4.0 0.0 'BASELINE': baselines_use,
411 1 461 461.0 0.0 'ANTENNA1': self.ant_1_array + 1,
412 1 515 515.0 0.0 'ANTENNA2': self.ant_2_array + 1,
413 1 409 409.0 0.0 'SUBARRAY': np.ones_like(self.ant_1_array),
414 1 4 4.0 0.0 'INTTIM': int_time_array}
415 1 5 5.0 0.0 pscal_dict = {'UU ': 1.0, 'VV ': 1.0, 'WW ': 1.0,
416 1 4 4.0 0.0 'DATE ': 1.0, 'BASELINE': 1.0, 'ANTENNA1': 1.0,
417 1 3 3.0 0.0 'ANTENNA2': 1.0, 'SUBARRAY': 1.0, 'INTTIM': 1.0}
418 1 5 5.0 0.0 pzero_dict = {'UU ': 0.0, 'VV ': 0.0, 'WW ': 0.0,
419 1 4 4.0 0.0 'DATE ': tzero, 'BASELINE': 0.0, 'ANTENNA1': 0.0,
420 1 4 4.0 0.0 'ANTENNA2': 0.0, 'SUBARRAY': 0.0, 'INTTIM': 0.0}
421
422 # list contains arrays of [u,v,w,date,baseline];
423 # each array has shape (Nblts)
424 1 515 515.0 0.0 if (np.max(self.ant_1_array) < 255 and
425 1 458 458.0 0.0 np.max(self.ant_2_array) < 255):
426 # if the number of antennas is less than 256 then include both the
427 # baseline array and the antenna arrays in the group parameters.
428 # Otherwise just use the antenna arrays
429 1 4 4.0 0.0 parnames_use = ['UU ', 'VV ', 'WW ',
430 1 4 4.0 0.0 'DATE ', 'BASELINE', 'ANTENNA1',
431 1 5 5.0 0.0 'ANTENNA2', 'SUBARRAY', 'INTTIM']
432 else:
433 parnames_use = ['UU ', 'VV ', 'WW ', 'DATE ',
434 'ANTENNA1', 'ANTENNA2', 'SUBARRAY', 'INTTIM']
435
436 1 4 4.0 0.0 group_parameter_list = [group_parameter_dict[parname] for
437 10 42 4.2 0.0 parname in parnames_use]
438 1 5 5.0 0.0 hdu = fits.GroupData(uvfits_array_data, parnames=parnames_use,
439 1 1963159 1963159.0 6.0 pardata=group_parameter_list, bitpix=-32)
440 1 9105 9105.0 0.0 hdu = fits.GroupsHDU(hdu)
441
442 10 42 4.2 0.0 for i, key in enumerate(parnames_use):
443 9 1908 212.0 0.0 hdu.header['PSCAL' + str(i + 1) + ' '] = pscal_dict[key]
444 9 1969 218.8 0.0 hdu.header['PZERO' + str(i + 1) + ' '] = pzero_dict[key]
445
446 # ISO string of first time in self.time_array
447 1 13 13.0 0.0 hdu.header['DATE-OBS'] = Time(self.time_array[0], scale='utc',
448 1 759 759.0 0.0 format='jd').isot
449
450 1 225 225.0 0.0 hdu.header['CTYPE2 '] = 'COMPLEX '
451 1 226 226.0 0.0 hdu.header['CRVAL2 '] = 1.0
452 1 224 224.0 0.0 hdu.header['CRPIX2 '] = 1.0
453 1 222 222.0 0.0 hdu.header['CDELT2 '] = 1.0
454
455 1 221 221.0 0.0 hdu.header['CTYPE3 '] = 'STOKES '
456 1 236 236.0 0.0 hdu.header['CRVAL3 '] = self.polarization_array[0]
457 1 228 228.0 0.0 hdu.header['CRPIX3 '] = 1.0
458 1 4 4.0 0.0 try:
459 1 252 252.0 0.0 hdu.header['CDELT3 '] = np.diff(self.polarization_array)[0]
460 except(IndexError):
461 hdu.header['CDELT3 '] = 1.0
462
463 1 226 226.0 0.0 hdu.header['CTYPE4 '] = 'FREQ '
464 1 247 247.0 0.0 hdu.header['CRVAL4 '] = self.freq_array[0, 0]
465 1 232 232.0 0.0 hdu.header['CRPIX4 '] = 1.0
466 1 262 262.0 0.0 hdu.header['CDELT4 '] = np.diff(self.freq_array[0])[0]
467
468 1 230 230.0 0.0 hdu.header['CTYPE5 '] = 'IF '
469 1 238 238.0 0.0 hdu.header['CRVAL5 '] = 1.0
470 1 234 234.0 0.0 hdu.header['CRPIX5 '] = 1.0
471 1 235 235.0 0.0 hdu.header['CDELT5 '] = 1.0
472
473 1 232 232.0 0.0 hdu.header['CTYPE6 '] = 'RA'
474 1 250 250.0 0.0 hdu.header['CRVAL6 '] = self.phase_center_ra_degrees
475
476 1 235 235.0 0.0 hdu.header['CTYPE7 '] = 'DEC'
477 1 243 243.0 0.0 hdu.header['CRVAL7 '] = self.phase_center_dec_degrees
478
479 1 237 237.0 0.0 hdu.header['BUNIT '] = self.vis_units
480 1 247 247.0 0.0 hdu.header['BSCALE '] = 1.0
481 1 243 243.0 0.0 hdu.header['BZERO '] = 0.0
482
483 1 246 246.0 0.0 hdu.header['OBJECT '] = self.object_name
484 1 241 241.0 0.0 hdu.header['TELESCOP'] = self.telescope_name
485 1 362 362.0 0.0 hdu.header['LAT '] = self.telescope_location_lat_lon_alt_degrees[0]
486 1 329 329.0 0.0 hdu.header['LON '] = self.telescope_location_lat_lon_alt_degrees[1]
487 1 329 329.0 0.0 hdu.header['ALT '] = self.telescope_location_lat_lon_alt[2]
488 1 258 258.0 0.0 hdu.header['INSTRUME'] = self.instrument
489 1 255 255.0 0.0 hdu.header['EPOCH '] = float(self.phase_center_epoch)
490
491 3 16 5.3 0.0 for line in self.history.splitlines():
492 2 1092 546.0 0.0 hdu.header.add_history(line)
493
494 # end standard keywords; begin user-defined keywords
495 1 8 8.0 0.0 for key, value in self.extra_keywords.iteritems():
496 # header keywords have to be 8 characters or less
497 keyword = key[:8].upper()
498 # print "keyword=-value-", keyword+'='+'-'+str(value)+'-'
499 if keyword == 'COMMENT':
500 for line in value.splitlines():
501 hdu.header.add_comment(line)
502 else:
503 hdu.header[keyword] = value
504
505 # ADD the ANTENNA table
506 1 10 10.0 0.0 staxof = np.zeros(self.Nants_telescope)
507
508 # 0 specifies alt-az, 6 would specify a phased array
509 1 7 7.0 0.0 mntsta = np.zeros(self.Nants_telescope)
510
511 # beware, X can mean just about anything
512 1 53 53.0 0.0 poltya = np.full((self.Nants_telescope), 'X', dtype=np.object_)
513 1 25 25.0 0.0 polaa = [90.0] + np.zeros(self.Nants_telescope)
514 1 18 18.0 0.0 poltyb = np.full((self.Nants_telescope), 'Y', dtype=np.object_)
515 1 17 17.0 0.0 polab = [0.0] + np.zeros(self.Nants_telescope)
516
517 1 4 4.0 0.0 col1 = fits.Column(name='ANNAME', format='8A',
518 1 492 492.0 0.0 array=self.antenna_names)
519 1 4 4.0 0.0 col2 = fits.Column(name='STABXYZ', format='3D',
520 1 484 484.0 0.0 array=self.antenna_positions)
521 # convert to 1-indexed from 0-indexed indicies
522 1 5 5.0 0.0 col3 = fits.Column(name='NOSTA', format='1J',
523 1 422 422.0 0.0 array=self.antenna_numbers + 1)
524 1 398 398.0 0.0 col4 = fits.Column(name='MNTSTA', format='1J', array=mntsta)
525 1 393 393.0 0.0 col5 = fits.Column(name='STAXOF', format='1E', array=staxof)
526 1 425 425.0 0.0 col6 = fits.Column(name='POLTYA', format='1A', array=poltya)
527 1 419 419.0 0.0 col7 = fits.Column(name='POLAA', format='1E', array=polaa)
528 # col8 = fits.Column(name='POLCALA', format='3E', array=polcala)
529 1 393 393.0 0.0 col9 = fits.Column(name='POLTYB', format='1A', array=poltyb)
530 1 392 392.0 0.0 col10 = fits.Column(name='POLAB', format='1E', array=polab)
531 # col11 = fits.Column(name='POLCALB', format='3E', array=polcalb)
532 # note ORBPARM is technically required, but we didn't put it in
533
534 1 4 4.0 0.0 cols = fits.ColDefs([col1, col2, col3, col4, col5, col6, col7, col9,
535 1 1626 1626.0 0.0 col10])
536
537 1 24300 24300.0 0.1 ant_hdu = fits.BinTableHDU.from_columns(cols)
538
539 1 208 208.0 0.0 ant_hdu.header['EXTNAME'] = 'AIPS AN'
540 1 198 198.0 0.0 ant_hdu.header['EXTVER'] = 1
541
542 # write XYZ coordinates if not already defined
543 1 227 227.0 0.0 ant_hdu.header['ARRAYX'] = self.telescope_location[0]
544 1 218 218.0 0.0 ant_hdu.header['ARRAYY'] = self.telescope_location[1]
545 1 216 216.0 0.0 ant_hdu.header['ARRAYZ'] = self.telescope_location[2]
546 1 205 205.0 0.0 ant_hdu.header['FRAME'] = 'ITRF'
547 1 216 216.0 0.0 ant_hdu.header['GSTIA0'] = self.gst0
548 1 220 220.0 0.0 ant_hdu.header['FREQ'] = self.freq_array[0, 0]
549 1 205 205.0 0.0 ant_hdu.header['RDATE'] = self.rdate
550 1 214 214.0 0.0 ant_hdu.header['UT1UTC'] = self.dut1
551
552 1 211 211.0 0.0 ant_hdu.header['TIMSYS'] = self.timesys
553 1 5 5.0 0.0 if self.timesys == 'IAT':
554 warnings.warn('This file has an "IAT" time system. Files of '
555 'this type are not properly supported')
556 1 215 215.0 0.0 ant_hdu.header['ARRNAM'] = self.telescope_name
557 1 207 207.0 0.0 ant_hdu.header['NO_IF'] = self.Nspws
558 1 218 218.0 0.0 ant_hdu.header['DEGPDY'] = self.earth_omega
559 # ant_hdu.header['IATUTC'] = 35.
560
561 # set mandatory parameters which are not supported by this object
562 # (or that we just don't understand)
563 1 207 207.0 0.0 ant_hdu.header['NUMORB'] = 0
564
565 # note: Bart had this set to 3. We've set it 0 after aips 117. -jph
566 1 211 211.0 0.0 ant_hdu.header['NOPCAL'] = 0
567
568 1 214 214.0 0.0 ant_hdu.header['POLTYPE'] = 'X-Y LIN'
569
570 # note: we do not support the concept of "frequency setups"
571 # -- lists of spws given in a SU table.
572 1 209 209.0 0.0 ant_hdu.header['FREQID'] = -1
573
574 # if there are offsets in images, this could be the culprit
575 1 284 284.0 0.0 ant_hdu.header['POLARX'] = 0.0
576 1 228 228.0 0.0 ant_hdu.header['POLARY'] = 0.0
577
578 1 216 216.0 0.0 ant_hdu.header['DATUTC'] = 0 # ONLY UTC SUPPORTED
579
580 # we always output right handed coordinates
581 1 219 219.0 0.0 ant_hdu.header['XYZHAND'] = 'RIGHT'
582
583 # ADD the FQ table
584 # skipping for now and limiting to a single spw
585
586 # write the file
587 1 92 92.0 0.0 hdulist = fits.HDUList(hdus=[hdu, ant_hdu])
588 1 9 9.0 0.0 if float(astropy.__version__[0:3]) < 1.3:
589 1 8180674 8180674.0 25.1 hdulist.writeto(filename, clobber=True)
590 else:
591 hdulist.writeto(filename, overwrite=True)