-
Notifications
You must be signed in to change notification settings - Fork 1
/
fhd.py
462 lines (401 loc) · 21.3 KB
/
fhd.py
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
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
# -*- mode: python; coding: utf-8 -*-
# Copyright (c) 2018 Radio Astronomy Software Group
# Licensed under the 2-clause BSD License
"""Class for reading FHD save files.
"""
from __future__ import absolute_import, division, print_function
from itertools import islice
import numpy as np
import warnings
from scipy.io.idl import readsav
from astropy import constants as const
from . import UVData
from . import utils as uvutils
from . import telescopes as uvtel
def get_fhd_history(settings_file, return_user=False):
"""
Small function to get the important history from an FHD settings text file.
Includes information about the command line call, the user, machine name and date
"""
with open(settings_file, 'r') as f:
settings_lines = f.readlines()
main_loc = None
command_loc = None
obs_loc = None
user_line = None
for ind, line in enumerate(settings_lines):
if line.startswith('##MAIN'):
main_loc = ind
if line.startswith('##COMMAND_LINE'):
command_loc = ind
if line.startswith('##OBS'):
obs_loc = ind
if line.startswith('User'):
user_line = ind
if (main_loc is not None and command_loc is not None
and obs_loc is not None and user_line is not None):
break
main_lines = settings_lines[main_loc + 1:command_loc]
command_lines = settings_lines[command_loc + 1:obs_loc]
history_lines = ['FHD history\n'] + main_lines + command_lines
for ind, line in enumerate(history_lines):
history_lines[ind] = line.rstrip().replace('\t', ' ')
history = '\n'.join(history_lines)
user = settings_lines[user_line].split()[1]
if return_user:
return history, user
else:
return history
class FHD(UVData):
"""
Defines a FHD-specific subclass of UVData for reading FHD save files.
This class should not be interacted with directly, instead use the read_fhd
method on the UVData class.
"""
def _latlonalt_close(self, latlonalt1, latlonalt2):
radian_tols = self._phase_center_ra.tols
loc_tols = self._telescope_location.tols
latlon_close = np.allclose(np.array(latlonalt1[0:2]),
np.array(latlonalt2[0:2]),
rtol=radian_tols[0], atol=radian_tols[1])
alt_close = np.isclose(latlonalt1[2], latlonalt2[2],
rtol=loc_tols[0], atol=loc_tols[1])
if latlon_close and alt_close:
return True
else:
return False
def _xyz_close(self, xyz1, xyz2):
loc_tols = self._telescope_location.tols
return np.allclose(xyz1, xyz2, rtol=loc_tols[0], atol=loc_tols[1])
def read_fhd(self, filelist, use_model=False, run_check=True, check_extra=True,
run_check_acceptability=True):
"""
Read in data from a list of FHD files.
Args:
filelist: The list of FHD save files to read from. Must include at
least one polarization file, a params file and a flag file.
use_model: Option to read in the model visibilities rather than the
dirty visibilities. Default is False.
run_check: Option to check for the existence and proper shapes of
parameters after reading in the file. Default is True.
check_extra: Option to check optional parameters as well as required
ones. Default is True.
run_check_acceptability: Option to check acceptable range of the values of
parameters after reading in the file. Default is True.
"""
datafiles = {}
params_file = None
flags_file = None
layout_file = None
settings_file = None
if use_model:
data_name = '_vis_model_'
else:
data_name = '_vis_'
for file in filelist:
if file.lower().endswith(data_name + 'xx.sav'):
if 'xx' in list(datafiles.keys()):
raise ValueError('multiple xx datafiles in filelist')
datafiles['xx'] = file
elif file.lower().endswith(data_name + 'yy.sav'):
if 'yy' in list(datafiles.keys()):
raise ValueError('multiple yy datafiles in filelist')
datafiles['yy'] = file
elif file.lower().endswith(data_name + 'xy.sav'):
if 'xy' in list(datafiles.keys()):
raise ValueError('multiple xy datafiles in filelist')
datafiles['xy'] = file
elif file.lower().endswith(data_name + 'yx.sav'):
if 'yx' in list(datafiles.keys()):
raise ValueError('multiple yx datafiles in filelist')
datafiles['yx'] = file
elif file.lower().endswith('_params.sav'):
if params_file is not None:
raise ValueError('multiple params files in filelist')
params_file = file
elif file.lower().endswith('_flags.sav'):
if flags_file is not None:
raise ValueError('multiple flags files in filelist')
flags_file = file
elif file.lower().endswith('_layout.sav'):
if layout_file is not None:
raise ValueError('multiple layout files in filelist')
layout_file = file
elif file.lower().endswith('_settings.txt'):
if settings_file is not None:
raise ValueError('multiple settings files in filelist')
settings_file = file
else:
# this is reached in tests but marked as uncovered because
# CPython's peephole optimizer replaces a jump to a continue
# with a jump to the top of the loop
continue # pragma: no cover
if len(datafiles) < 1:
raise Exception('No data files included in file list')
if params_file is None:
raise Exception('No params file included in file list')
if flags_file is None:
raise Exception('No flags file included in file list')
if layout_file is None:
warnings.warn('No layout file included in file list. '
'Support for FHD data without layout files will be '
'deprecated in a future version.', PendingDeprecationWarning)
if settings_file is None:
warnings.warn('No settings file included in file list')
# TODO: add checking to make sure params, flags and datafiles are
# consistent with each other
vis_data = {}
for pol, file in datafiles.items():
this_dict = readsav(file, python_dict=True)
if use_model:
vis_data[pol] = this_dict['vis_model_ptr']
else:
vis_data[pol] = this_dict['vis_ptr']
this_obs = this_dict['obs']
data_shape = vis_data[pol].shape
obs = this_obs
bl_info = obs['BASELINE_INFO'][0]
meta_data = obs['META_DATA'][0]
astrometry = obs['ASTR'][0]
fhd_pol_list = []
for pol in obs['POL_NAMES'][0]:
fhd_pol_list.append(uvutils._bytes_to_str(pol).lower())
params_dict = readsav(params_file, python_dict=True)
params = params_dict['params']
flag_file_dict = readsav(flags_file, python_dict=True)
# The name for this variable changed recently (July 2016). Test for both.
vis_weights_data = {}
if 'flag_arr' in flag_file_dict:
weights_key = 'flag_arr'
elif 'vis_weights' in flag_file_dict:
weights_key = 'vis_weights'
else:
raise ValueError('No recognized key for visibility weights in flags_file.')
for index, w in enumerate(flag_file_dict[weights_key]):
vis_weights_data[fhd_pol_list[index]] = w
self.Ntimes = int(obs['N_TIME'][0])
self.Nbls = int(obs['NBASELINES'][0])
self.Nblts = data_shape[0]
self.Nfreqs = int(obs['N_FREQ'][0])
self.Npols = len(list(vis_data.keys()))
self.Nspws = 1
self.spw_array = np.array([0])
self.vis_units = 'JY'
lin_pol_order = ['xx', 'yy', 'xy', 'yx']
linear_pol_dict = dict(zip(lin_pol_order, np.arange(5, 9) * -1))
pol_list = []
for pol in lin_pol_order:
if pol in vis_data:
pol_list.append(linear_pol_dict[pol])
self.polarization_array = np.asarray(pol_list)
self.data_array = np.zeros((self.Nblts, self.Nspws, self.Nfreqs,
self.Npols), dtype=np.complex_)
self.nsample_array = np.zeros((self.Nblts, self.Nspws, self.Nfreqs,
self.Npols), dtype=np.float_)
self.flag_array = np.zeros((self.Nblts, self.Nspws, self.Nfreqs,
self.Npols), dtype=np.bool_)
for pol, vis in vis_data.items():
pol_i = pol_list.index(linear_pol_dict[pol])
self.data_array[:, 0, :, pol_i] = vis
self.flag_array[:, 0, :, pol_i] = vis_weights_data[pol] <= 0
self.nsample_array[:, 0, :, pol_i] = np.abs(vis_weights_data[pol])
# In FHD, uvws are in seconds not meters!
self.uvw_array = np.zeros((self.Nblts, 3))
self.uvw_array[:, 0] = params['UU'][0] * const.c.to('m/s').value
self.uvw_array[:, 1] = params['VV'][0] * const.c.to('m/s').value
self.uvw_array[:, 2] = params['WW'][0] * const.c.to('m/s').value
# bl_info.JDATE (a vector of length Ntimes) is the only safe date/time
# to use in FHD files.
# (obs.JD0 (float) and params.TIME (vector of length Nblts) are
# context dependent and are not safe
# because they depend on the phasing of the visibilities)
# the values in bl_info.JDATE are the JD for each integration.
# We need to expand up to Nblts.
int_times = bl_info['JDATE'][0]
bin_offset = bl_info['BIN_OFFSET'][0]
if self.Ntimes != len(int_times):
warnings.warn('Ntimes does not match the number of unique times in the data')
self.time_array = np.zeros(self.Nblts)
if self.Ntimes == 1:
self.time_array.fill(int_times[0])
else:
for ii in range(0, len(int_times)):
if ii < (len(int_times) - 1):
self.time_array[bin_offset[ii]:bin_offset[ii + 1]] = int_times[ii]
else:
self.time_array[bin_offset[ii]:] = int_times[ii]
# Note that FHD antenna arrays are 1-indexed so we subtract 1
# to get 0-indexed arrays
self.ant_1_array = bl_info['TILE_A'][0] - 1
self.ant_2_array = bl_info['TILE_B'][0] - 1
self.Nants_data = int(len(np.unique(self.ant_1_array.tolist() + self.ant_2_array.tolist())))
self.baseline_array = \
self.antnums_to_baseline(self.ant_1_array,
self.ant_2_array)
if self.Nbls != len(np.unique(self.baseline_array)):
warnings.warn('Nbls does not match the number of unique baselines in the data')
self.freq_array = np.zeros((self.Nspws, len(bl_info['FREQ'][0])), dtype=np.float_)
self.freq_array[0, :] = bl_info['FREQ'][0]
if not np.isclose(obs['OBSRA'][0], obs['PHASERA'][0]) or \
not np.isclose(obs['OBSDEC'][0], obs['PHASEDEC'][0]):
warnings.warn('These visibilities may have been phased '
'improperly -- without changing the uvw locations')
self.set_phased()
self.phase_center_ra_degrees = np.float(obs['OBSRA'][0])
self.phase_center_dec_degrees = np.float(obs['OBSDEC'][0])
# this is generated in FHD by subtracting the JD of neighboring
# integrations. This can have limited accuracy, so it can be slightly
# off the actual value.
# (e.g. 1.999426... rather than 2)
time_res = obs['TIME_RES']
# time_res is constrained to be a scalar currently
self.integration_time = (np.ones_like(self.time_array, dtype=np.float64)
* time_res[0])
self.channel_width = float(obs['FREQ_RES'][0])
# # --- observation information ---
self.telescope_name = uvutils._bytes_to_str(obs['INSTRUMENT'][0])
# This is a bit of a kludge because nothing like object_name exists
# in FHD files.
# At least for the MWA, obs.ORIG_PHASERA and obs.ORIG_PHASEDEC specify
# the field the telescope was nominally pointing at
# (May need to be revisited, but probably isn't too important)
self.object_name = 'Field RA(deg): ' + str(obs['ORIG_PHASERA'][0]) + \
', Dec:' + str(obs['ORIG_PHASEDEC'][0])
# For the MWA, this can sometimes be converted to EoR fields
if self.telescope_name.lower() == 'mwa':
if np.isclose(obs['ORIG_PHASERA'][0], 0) and \
np.isclose(obs['ORIG_PHASEDEC'][0], -27):
object_name = 'EoR 0 Field'
self.instrument = self.telescope_name
latitude = np.deg2rad(float(obs['LAT'][0]))
longitude = np.deg2rad(float(obs['LON'][0]))
altitude = float(obs['ALT'][0])
# history: add the first few lines from the settings file
if settings_file is not None:
self.history = get_fhd_history(settings_file)
else:
self.history = ''
if not uvutils._check_history_version(self.history, self.pyuvdata_version_str):
self.history += self.pyuvdata_version_str
self.phase_center_epoch = astrometry['EQUINOX'][0]
# get the stuff FHD read from the antenna table (in layout file)
if layout_file is not None:
layout_dict = readsav(layout_file, python_dict=True)
layout = layout_dict['layout']
layout_fields = [name.lower() for name in layout.dtype.names]
# Try to get the telescope location from the layout file &
# compare it to the position from the obs structure.
arr_center = layout['array_center'][0]
layout_fields.remove('array_center')
xyz_telescope_frame = uvutils._bytes_to_str(layout['coordinate_frame'][0]).lower()
layout_fields.remove('coordinate_frame')
if xyz_telescope_frame == 'itrf':
# compare to lat/lon/alt
location_latlonalt = uvutils.XYZ_from_LatLonAlt(latitude, longitude, altitude)
latlonalt_arr_center = uvutils.LatLonAlt_from_XYZ(arr_center)
# check both lat/lon/alt and xyz because of subtle differences in tolerances
if (self._xyz_close(location_latlonalt, arr_center)
or self._latlonalt_close((latitude, longitude, altitude), latlonalt_arr_center)):
self.telescope_location = arr_center
else:
# values do not agree with each other to within the tolerances.
# this is a known issue with FHD runs on cotter uvfits files for the MWA
# compare with the known_telescopes values
telescope_obj = uvtel.get_telescope(self.telescope_name)
# start warning message
message = ('Telescope location derived from obs lat/lon/alt '
'values does not match the location in the layout file.')
if telescope_obj is not False:
if self._latlonalt_close((latitude, longitude, altitude),
telescope_obj.telescope_location_lat_lon_alt):
# obs lat/lon/alt matches known_telescopes
message += (' Value from obs lat/lon/alt matches the '
'known_telescopes values, using them.')
self.telescope_location = location_latlonalt
elif self._xyz_close(arr_center, telescope_obj.telescope_location):
# layout xyz matches known_telescopes
message += (' Value from the layout file matches the '
'known_telescopes values, using them.')
self.telescope_location = arr_center
else:
# None of the values match each other. Defaulting to known_telescopes value.
message += (' Neither location matches the values '
'in known_telescopes. Defaulting to '
'using the known_telescopes values.')
self.telescope_location = telescope_obj.telescope_location
else:
message += (' Telescope is not in known_telescopes. '
'Defaulting to using the obs derived values.')
self.telescope_location = location_latlonalt
# issue warning
warnings.warn(message)
else:
self.telescope_location_lat_lon_alt = (latitude, longitude, altitude)
self.antenna_positions = layout['antenna_coords'][0]
layout_fields.remove('antenna_coords')
self.antenna_names = [uvutils._bytes_to_str(ant).strip() for ant in layout['antenna_names'][0].tolist()]
layout_fields.remove('antenna_names')
# make these 0-indexed (rather than one indexed)
self.antenna_numbers = layout['antenna_numbers'][0] - 1
layout_fields.remove('antenna_numbers')
self.Nants_telescope = int(layout['n_antenna'][0])
layout_fields.remove('n_antenna')
if self.telescope_name.lower() == 'mwa':
# check that obs.baseline_info.tile_names match the antenna names
# this only applies for MWA because the tile_names come from metafits files
obs_tile_names = [uvutils._bytes_to_str(ant).strip() for ant in bl_info['TILE_NAMES'][0].tolist()]
obs_tile_names = ['Tile' + '0' * (3 - len(ant)) + ant for ant in obs_tile_names]
# tile_names are assumed to be ordered: so their index gives the antenna number
# make an comparison array from self.antenna_names ordered this way.
ant_names = np.zeros((np.max(self.antenna_numbers) + 1), str).tolist()
for index, number in enumerate(self.antenna_numbers):
ant_names[number] = self.antenna_names[index]
if obs_tile_names != ant_names:
warnings.warn('tile_names from obs structure does not match antenna_names from layout')
self.gst0 = float(layout['gst0'][0])
layout_fields.remove('gst0')
self.rdate = uvutils._bytes_to_str(layout['ref_date'][0]).lower()
layout_fields.remove('ref_date')
self.earth_omega = float(layout['earth_degpd'][0])
layout_fields.remove('earth_degpd')
self.dut1 = float(layout['dut1'][0])
layout_fields.remove('dut1')
self.timesys = uvutils._bytes_to_str(layout['time_system'][0]).upper().strip()
layout_fields.remove('time_system')
if 'diameters' in layout_fields:
self.timesys = uvutils._bytes_to_str(layout['time_system'][0]).upper().strip()
layout_fields.remove('diameters')
# ignore some fields, put everything else in extra_keywords
layout_fields_ignore = ['diff_utc', 'pol_type', 'n_pol_cal_params',
'mount_type', 'axis_offset',
'pola', 'pola_orientation', 'pola_cal_params',
'polb', 'polb_orientation', 'polb_cal_params',
'beam_fwhm']
for field in layout_fields_ignore:
if field in layout_fields:
layout_fields.remove(field)
for field in layout_fields:
keyword = field
if len(keyword) > 8:
keyword = field.replace('_', '')
value = layout[field][0]
if isinstance(value, bytes):
value = uvutils._bytes_to_str(value)
self.extra_keywords[keyword.upper()] = value
else:
self.telescope_location_lat_lon_alt = (latitude, longitude, altitude)
self.antenna_names = [uvutils._bytes_to_str(ant).strip() for ant in bl_info['TILE_NAMES'][0].tolist()]
if self.telescope_name.lower() == 'mwa':
self.antenna_names = ['Tile' + '0' * (3 - len(ant)) + ant for ant in self.antenna_names]
self.Nants_telescope = len(self.antenna_names)
self.antenna_numbers = np.arange(self.Nants_telescope)
try:
self.set_telescope_params()
except ValueError as ve:
warnings.warn(str(ve))
# need to make sure telescope location is defined properly before this call
self.set_lsts_from_time_array()
# check if object has all required uv_properties set
if run_check:
self.check(check_extra=check_extra,
run_check_acceptability=run_check_acceptability)