-
Notifications
You must be signed in to change notification settings - Fork 1
/
miriad.py
1083 lines (988 loc) · 54.8 KB
/
miriad.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
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
"""Class for reading and writing Miriad files."""
from astropy import constants as const
import os
import shutil
import numpy as np
import copy
import warnings
import aipy
from uvdata import UVData
import telescopes as uvtel
import utils as uvutils
import itertools
class Miriad(UVData):
"""
Defines a Miriad-specific subclass of UVData for reading and writing Miriad files.
This class should not be interacted with directly, instead use the read_miriad
and write_miriad methods on the UVData class.
"""
def _pol_to_ind(self, pol):
if self.polarization_array is None:
raise(ValueError, "Can't index polarization {p} because "
"polarization_array is not set".format(p=pol))
pol_ind = np.argwhere(self.polarization_array == pol)
if len(pol_ind) != 1:
raise(ValueError, "multiple matches for pol={pol} in "
"polarization_array".format(pol=pol))
return pol_ind
def read_miriad(self, filepath, correct_lat_lon=True, run_check=True,
check_extra=True, run_check_acceptability=True, phase_type=None,
antenna_nums=None, ant_str=None, ant_pairs_nums=None,
polarizations=None, time_range=None):
"""
Read in data from a miriad file.
Args:
filepath: The miriad file directory to read from.
correct_lat_lon: flag -- that only matters if altitude is missing --
to update the latitude and longitude from the known_telescopes list
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.
antenna_nums: The antennas numbers to only read into the object.
ant_pairs_nums: A list of antenna number tuples (e.g. [(0,1), (3,2)])
specifying baselines to read into the object. Ordering of the
numbers within the tuple does not matter. A single antenna iterable
e.g. (1,) is interpreted as all visibilities with that antenna.
ant_str: A string containing information about what kinds of visibility data
to read-in. Can be 'auto', 'cross', 'all'. Cannot provide ant_str if
antenna_nums and/or ant_pairs_nums is not None.
polarizations: List of polarization integers or strings to read-in.
Ex: ['xx', 'yy', ...]
time_range: len-2 list containing min and max range of times (Julian Date) to read-in.
Ex: [2458115.20, 2458115.40]
"""
if not os.path.exists(filepath):
raise(IOError, filepath + ' not found')
uv = aipy.miriad.UV(filepath)
# list of miriad variables always read
# NB: this includes variables in try/except (i.e. not all variables are
# necessarily present in the miriad file)
default_miriad_variables = ['nchan', 'npol', 'inttime', 'sdf',
'source', 'telescop', 'latitud', 'longitu',
'altitude', 'history', 'visunits',
'instrume', 'dut1', 'gst0', 'rdate',
'timesys', 'xorient', 'cnt', 'ra', 'dec',
'lst', 'pol', 'nants', 'antnames', 'nblts',
'ntimes', 'nbls', 'sfreq', 'epoch',
'antpos', 'antnums', 'degpdy', 'antdiam',
]
# list of miriad variables not read, but also not interesting
# NB: nspect (I think) is number of spectral windows, will want one day
# NB: xyphase & xyamp are "On-line X Y phase/amplitude measurements" which we may want in
# a calibration object some day
# NB: systemp, xtsys & ytsys are "System temperatures of the antenna/X/Y feed"
# which we may want in a calibration object some day
# NB: freqs, leakage and bandpass may be part of a calibration object some day
other_miriad_variables = ['nspect', 'obsdec', 'vsource', 'ischan',
'restfreq', 'nschan', 'corr', 'freq',
'freqs', 'leakage', 'bandpass',
'tscale', 'coord', 'veldop', 'time', 'obsra',
'operator', 'version', 'axismax', 'axisrms',
'xyphase', 'xyamp', 'systemp', 'xtsys', 'ytsys',
'baseline']
extra_miriad_variables = []
for variable in uv.vars():
if (variable not in default_miriad_variables
and variable not in other_miriad_variables):
extra_miriad_variables.append(variable)
miriad_header_data = {'Nfreqs': 'nchan',
'Npols': 'npol',
'integration_time': 'inttime',
'channel_width': 'sdf', # in Ghz!
'object_name': 'source',
'telescope_name': 'telescop'
}
for item in miriad_header_data:
if isinstance(uv[miriad_header_data[item]], str):
header_value = uv[miriad_header_data[item]].replace('\x00', '')
else:
header_value = uv[miriad_header_data[item]]
setattr(self, item, header_value)
latitude = uv['latitud'] # in units of radians
longitude = uv['longitu']
try:
altitude = uv['altitude']
self.telescope_location_lat_lon_alt = (latitude, longitude, altitude)
except(KeyError):
# get info from known telescopes. Check to make sure the lat/lon values match reasonably well
telescope_obj = uvtel.get_telescope(self.telescope_name)
if telescope_obj is not False:
tol = 2 * np.pi * 1e-3 / (60.0 * 60.0 * 24.0) # 1mas in radians
lat_close = np.isclose(telescope_obj.telescope_location_lat_lon_alt[0],
latitude, rtol=0, atol=tol)
lon_close = np.isclose(telescope_obj.telescope_location_lat_lon_alt[1],
longitude, rtol=0, atol=tol)
if correct_lat_lon:
self.telescope_location_lat_lon_alt = telescope_obj.telescope_location_lat_lon_alt
else:
self.telescope_location_lat_lon_alt = (latitude, longitude, telescope_obj.telescope_location_lat_lon_alt[2])
if lat_close and lon_close:
if correct_lat_lon:
warnings.warn('Altitude is not present in Miriad file, '
'using known location values for '
'{telescope_name}.'.format(telescope_name=telescope_obj.telescope_name))
else:
warnings.warn('Altitude is not present in Miriad file, '
'using known location altitude value '
'for {telescope_name} and lat/lon from '
'file.'.format(telescope_name=telescope_obj.telescope_name))
else:
warn_string = ('Altitude is not present in file ')
if not lat_close and not lon_close:
warn_string = warn_string + 'and latitude and longitude values do not match values '
else:
if not lat_close:
warn_string = warn_string + 'and latitude value does not match value '
else:
warn_string = warn_string + 'and longitude value does not match value '
if correct_lat_lon:
warn_string = (warn_string + 'for {telescope_name} in known '
'telescopes. Using values from known '
'telescopes.'.format(telescope_name=telescope_obj.telescope_name))
warnings.warn(warn_string)
else:
warn_string = (warn_string + 'for {telescope_name} in known '
'telescopes. Using altitude value from known '
'telescopes and lat/lon from '
'file.'.format(telescope_name=telescope_obj.telescope_name))
warnings.warn(warn_string)
else:
warnings.warn('Altitude is not present in Miriad file, and '
'telescope {telescope_name} is not in '
'known_telescopes. Telescope location will be '
'set using antenna positions.'
.format(telescope_name=self.telescope_name))
self.history = uv['history']
if not uvutils.check_history_version(self.history, self.pyuvdata_version_str):
self.history += self.pyuvdata_version_str
self.channel_width *= 1e9 # change from GHz to Hz
# check for pyuvdata variables that are not recognized miriad variables
if 'visunits' in uv.vartable.keys():
self.vis_units = uv['visunits'].replace('\x00', '')
else:
self.vis_units = 'UNCALIB' # assume no calibration
if 'instrume' in uv.vartable.keys():
self.instrument = uv['instrume'].replace('\x00', '')
else:
self.instrument = self.telescope_name # set instrument = telescope
if 'dut1' in uv.vartable.keys():
self.dut1 = uv['dut1']
if 'degpdy' in uv.vartable.keys():
self.earth_omega = uv['degpdy']
if 'gst0' in uv.vartable.keys():
self.gst0 = uv['gst0']
if 'rdate' in uv.vartable.keys():
self.rdate = uv['rdate'].replace('\x00', '')
if 'timesys' in uv.vartable.keys():
self.timesys = uv['timesys'].replace('\x00', '')
if 'xorient' in uv.vartable.keys():
self.x_orientation = uv['xorient'].replace('\x00', '')
# read through the file and get the data
_source = uv['source'] # check source of initial visibility
# dict of extra variables to see if they change
check_variables = {}
for extra_variable in extra_miriad_variables:
check_variables[extra_variable] = uv[extra_variable]
history_update_string = ' Downselected to specific '
n_selects = 0
# select on ant_str if provided
if ant_str is not None:
# type check
assert isinstance(ant_str, (str, np.str)), "ant_str must be fed as a string"
assert antenna_nums is None and ant_pairs_nums is None, "ant_str must be None if antenna_nums or ant_pairs_nums is not None"
aipy.scripting.uv_selector(uv, ant_str)
if ant_str != 'all':
history_update_string += 'antenna pairs'
n_selects += 1
# select on antenna_nums and/or ant_pairs_nums using aipy.scripting.uv_selector
if antenna_nums is not None or ant_pairs_nums is not None:
antpair_str = ''
if antenna_nums is not None:
# type check
err_msg = "antenna_nums must be fed as a list of antenna number integers"
assert isinstance(antenna_nums, (np.ndarray, list)), err_msg
assert isinstance(antenna_nums[0], (int, np.int, np.int32)), err_msg
# get all possible combinations
antpairs = list(itertools.combinations_with_replacement(antenna_nums, 2))
# convert antenna numbers to string form required by aipy.scripting.uv_selector
antpair_str += ','.join(map(lambda ap: '_'.join(map(lambda a: str(a), ap)), antpairs))
history_update_string += 'antennas'
n_selects += 1
if ant_pairs_nums is not None:
# type check
err_msg = "ant_pairs_nums must be a list of antnum integer tuples, Ex: [(0, 1), ...]"
assert isinstance(ant_pairs_nums, list), err_msg
assert np.array(map(lambda ap: isinstance(ap, tuple), ant_pairs_nums)).all(), err_msg
assert np.array(map(lambda ap: map(lambda a: isinstance(a, (int, np.int, np.int32)), ap), ant_pairs_nums)).all(), err_msg
# convert ant-pair tuples to string form required by aipy.scripting.uv_selector
if len(antpair_str) > 0:
antpair_str += ','
antpair_str += ','.join(map(lambda ap: '_'.join(map(lambda a: str(a), ap)), ant_pairs_nums))
if n_selects > 0:
history_update_string += ', antenna pairs'
else:
history_update_string += 'antenna pairs'
n_selects += 1
aipy.scripting.uv_selector(uv, antpair_str)
# select on time range
if time_range is not None:
# type check
err_msg = "time_range must be a len-2 list of Julian Date floats, Ex: [2458115.2, 2458115.6]"
assert isinstance(time_range, (list, np.ndarray)), err_msg
assert len(time_range) == 2, err_msg
assert np.array(map(lambda t: isinstance(t, (float, np.float, np.float64)), time_range)).all(), err_msg
uv.select('time', time_range[0], time_range[1], include=True)
if n_selects > 0:
history_update_string += ', times'
else:
history_update_string += 'times'
n_selects += 1
# select on polarizations
if polarizations is not None:
# type check
err_msg = "pols must be a list of polarization strings or ints, Ex: ['xx', ...] or [-5, ...]"
assert isinstance(polarizations, (list, np.ndarray)), err_msg
assert np.array(map(lambda p: isinstance(p, (str, np.str, int, np.int, np.int32)), polarizations)).all(), err_msg
# convert to pol integer if string
polarizations = [p if isinstance(p, (int, np.int, np.int32)) else uvutils.polstr2num(p) for p in polarizations]
# iterate through all possible pols and reject if not in pols
pol_list = []
for p in np.arange(-8, 5):
if p not in polarizations:
uv.select('polarization', p, p, include=False)
else:
pol_list.append(p)
# assert not empty
assert len(pol_list) > 0, "No polarizations in data matched {}".format(polarizations)
if n_selects > 0:
history_update_string += ', polarizations'
else:
history_update_string += 'polarizations'
n_selects += 1
history_update_string += ' using pyuvdata.'
if n_selects > 0:
self.history += history_update_string
data_accumulator = {}
pol_list = []
for (uvw, t, (i, j)), d, f in uv.all(raw=True):
# control for the case of only a single spw not showing up in
# the dimension
# Note that the (i, j) tuple is calculated from a baseline number in
# aipy (see miriad_wrap.h). The i, j values are also adjusted by aipy
# to start at 0 rather than 1.
if len(d.shape) == 1:
d.shape = (1,) + d.shape
self.Nspws = d.shape[0]
self.spw_array = np.arange(self.Nspws)
else:
raise(ValueError, """Sorry. Files with more than one spectral
window (spw) are not yet supported. A great
project for the interested student!""")
try:
cnt = uv['cnt']
except(KeyError):
cnt = np.ones(d.shape, dtype=np.float)
ra = uv['ra']
dec = uv['dec']
lst = uv['lst']
source = uv['source']
if source != _source:
raise(ValueError, 'This appears to be a multi source file, which is not supported.')
else:
_source = source
# check extra variables for changes compared with initial value
for extra_variable in check_variables.keys():
if type(check_variables[extra_variable]) == str:
if uv[extra_variable] != check_variables[extra_variable]:
check_variables.pop(extra_variable)
else:
if not np.allclose(uv[extra_variable],
check_variables[extra_variable]):
check_variables.pop(extra_variable)
try:
data_accumulator[uv['pol']].append([uvw, t, i, j, d, f, cnt,
ra, dec])
except(KeyError):
data_accumulator[uv['pol']] = [[uvw, t, i, j, d, f, cnt,
ra, dec]]
pol_list.append(uv['pol'])
# NB: flag types in miriad are usually ints
if len(data_accumulator.keys()) == 0:
raise ValueError('No data is present, probably as a result of '
'select on read that excludes all the data')
# keep all single valued extra_variables as extra_keywords
for key in check_variables.keys():
if type(check_variables[key]) == str:
value = check_variables[key].replace('\x00', '')
# check for booleans encoded as strings
if value == 'True':
value = True
elif value == 'False':
value = False
self.extra_keywords[key] = value
else:
self.extra_keywords[key] = check_variables[key]
# Check for items in itemtable to put into extra_keywords
# These will end up as variables in written files, but is internally consistent.
for key in uv.items():
# A few items that are not needed, we read elsewhere, or is not supported
# when downselecting, so we don't read here.
if key not in ['vartable', 'history', 'obstype'] and key not in other_miriad_variables:
if type(uv[key]) == str:
value = uv[key].replace('\x00', '')
value = uv[key].replace('\x01', '')
if value == 'True':
value = True
elif value == 'False':
value = False
self.extra_keywords[key] = value
else:
self.extra_keywords[key] = uv[key]
for pol, data in data_accumulator.iteritems():
data_accumulator[pol] = np.array(data)
self.polarization_array = np.array(pol_list)
if polarizations is None:
# A select on read would make the header npols not match the pols in the data
if len(self.polarization_array) != self.Npols:
warnings.warn('npols={npols} but found {n} pols in data file'.format(
npols=self.Npols, n=len(self.polarization_array)))
self.Npols = len(pol_list)
# makes a data_array (and flag_array) of zeroes to be filled in by
# data values
# any missing data will have zeros
# use set to get the unique list of all times ever listed in the file
# iterate over polarizations and all spectra (bls and times) in two
# nested loops, then flatten into a single vector, then set
# then list again.
times = list(set(
np.concatenate([[k[1] for k in d] for d in data_accumulator.values()])))
times = np.sort(times)
ant_i_unique = list(set(
np.concatenate([[k[2] for k in d] for d in data_accumulator.values()])))
ant_j_unique = list(set(
np.concatenate([[k[3] for k in d] for d in data_accumulator.values()])))
sorted_unique_ants = sorted(list(set(ant_i_unique + ant_j_unique)))
ant_i_unique = np.array(ant_i_unique)
ant_j_unique = np.array(ant_j_unique)
# Determine maximum digits needed to distinguish different values
if sorted_unique_ants[-1] > 0:
ndig_ant = np.ceil(np.log10(sorted_unique_ants[-1])).astype(int) + 1
else:
ndig_ant = 1
# Be excessive in precision because we use the floating point values as dictionary keys later
prec_t = - 2 * np.floor(np.log10(self._time_array.tols[-1])).astype(int)
ndig_t = (np.ceil(np.log10(times[-1])).astype(int) + prec_t + 2)
blts = []
for d in data_accumulator.values():
for k in d:
blt = ["{1:.{0}f}".format(prec_t, k[1]).zfill(ndig_t),
str(k[2]).zfill(ndig_ant), str(k[3]).zfill(ndig_ant)]
blt = "_".join(blt)
blts.append(blt)
unique_blts = np.unique(np.array(blts))
reverse_inds = dict(zip(unique_blts, range(len(unique_blts))))
self.Nants_data = len(sorted_unique_ants)
# Miriad has no way to keep track of antenna numbers, so the antenna
# numbers are simply the index for each antenna in any array that
# describes antenna attributes (e.g. antpos for the antenna_postions).
# Therefore on write, nants (which gives the size of the antpos array)
# needs to be increased to be the max value of antenna_numbers+1 and the
# antpos array needs to be inflated with zeros at locations where we
# don't have antenna information. These inflations need to be undone at
# read. If the file was written by pyuvdata, then the variable antnums
# will be present and we can use it, otherwise we need to test for zeros
# in the antpos array and/or antennas with no visibilities.
try:
# The antnums variable will only exist if the file was written with pyuvdata.
# For some reason Miriad doesn't handle an array of integers properly,
# so we convert to floats on write and back here
self.antenna_numbers = uv['antnums'].astype(int)
self.Nants_telescope = len(self.antenna_numbers)
except(KeyError):
self.antenna_numbers = None
self.Nants_telescope = None
nants = uv['nants']
try:
# Miriad stores antpos values in units of ns, pyuvdata uses meters.
antpos = uv['antpos'].reshape(3, nants).T * const.c.to('m/ns').value
# first figure out what are good antenna positions so we can only
# use the non-zero ones to evaluate position information
antpos_length = np.sqrt(np.sum(np.abs(antpos)**2, axis=1))
good_antpos = np.where(antpos_length > 0)[0]
mean_antpos_length = np.mean(antpos_length[good_antpos])
if mean_antpos_length > 6.35e6 and mean_antpos_length < 6.39e6:
absolute_positions = True
else:
absolute_positions = False
# Miriad stores antpos values in a rotated ECEF coordinate system
# where the x-axis goes through the local meridan. Need to convert
# these positions back to standard ECEF and if they are absolute positions,
# subtract off the telescope position to make them relative to the array center.
ecef_antpos = uvutils.ECEF_from_rotECEF(antpos, longitude)
if self.telescope_location is not None:
if absolute_positions:
rel_ecef_antpos = ecef_antpos - self.telescope_location
# maintain zeros because they mark missing data
rel_ecef_antpos[np.where(antpos_length == 0)[0]] = ecef_antpos[np.where(antpos_length == 0)[0]]
else:
rel_ecef_antpos = ecef_antpos
else:
self.telescope_location = np.mean(ecef_antpos[good_antpos, :], axis=0)
valid_location = self._telescope_location.check_acceptability()[0]
# check to see if this could be a valid telescope_location
if valid_location:
mean_lat, mean_lon, mean_alt = self.telescope_location_lat_lon_alt
tol = 2 * np.pi / (60.0 * 60.0 * 24.0) # 1 arcsecond in radians
mean_lat_close = np.isclose(mean_lat, latitude, rtol=0, atol=tol)
mean_lon_close = np.isclose(mean_lon, longitude, rtol=0, atol=tol)
if mean_lat_close and mean_lon_close:
# this looks like a valid telescope_location, and the
# mean antenna lat & lon values are close. Set the
# telescope_location using the file lat/lons and the mean alt.
# Then subtract it off of the antenna positions
warnings.warn('Telescope location is not set, but antenna '
'positions are present. Mean antenna latitude and '
'longitude values match file values, so '
'telescope_position will be set using the '
'mean of the antenna altitudes')
self.telescope_location_lat_lon_alt = (latitude, longitude, mean_alt)
rel_ecef_antpos = ecef_antpos - self.telescope_location
else:
# this looks like a valid telescope_location, but the
# mean antenna lat & lon values are not close. Set the
# telescope_location using the file lat/lons at sea level.
# Then subtract it off of the antenna positions
self.telescope_location_lat_lon_alt = (latitude, longitude, 0)
warn_string = ('Telescope location is set at sealevel at '
'the file lat/lon coordinates. Antenna '
'positions are present, but the mean '
'antenna ')
rel_ecef_antpos = ecef_antpos - self.telescope_location
if not mean_lat_close and not mean_lon_close:
warn_string += ('latitude and longitude values do not '
'match file values so they are not used '
'for altiude.')
elif not mean_lat_close:
warn_string += ('latitude value does not '
'match file values so they are not used '
'for altiude.')
else:
warn_string += ('longitude value does not '
'match file values so they are not used '
'for altiude.')
warnings.warn(warn_string)
else:
# This does not give a valid telescope_location. Instead
# calculate it from the file lat/lon and sea level for altiude
self.telescope_location_lat_lon_alt = (latitude, longitude, 0)
warn_string = ('Telescope location is set at sealevel at '
'the file lat/lon coordinates. Antenna '
'positions are present, but the mean '
'antenna ')
warn_string += ('position does not give a '
'telescope_location on the surface of the '
'earth.')
if absolute_positions:
rel_ecef_antpos = ecef_antpos - self.telescope_location
else:
warn_string += (' Antenna positions do not appear to be '
'on the surface of the earth and will be treated '
'as relative.')
rel_ecef_antpos = ecef_antpos
warnings.warn(warn_string)
if self.Nants_telescope is not None:
# in this case there is an antnums variable
# (meaning that the file was written with pyuvdata), so we'll use it
if nants == self.Nants_telescope:
# no inflation, so just use the positions
self.antenna_positions = rel_ecef_antpos
else:
# there is some inflation, just use the ones that appear in antnums
self.antenna_positions = np.zeros((self.Nants_telescope, 3), dtype=antpos.dtype)
for ai, num in enumerate(self.antenna_numbers):
self.antenna_positions[ai, :] = rel_ecef_antpos[num, :]
else:
# there is no antnums variable (meaning that this file was not
# written by pyuvdata), so we test for antennas with non-zero
# positions and/or that appear in the visibility data
# (meaning that they have entries in ant_1_array or ant_2_array)
antpos_length = np.sqrt(np.sum(np.abs(antpos)**2, axis=1))
good_antpos = np.where(antpos_length > 0)[0]
# take the union of the antennas with good positions (good_antpos)
# and the antennas that have visisbilities (sorted_unique_ants)
# if there are antennas with visibilities but zeroed positions we issue a warning below
ants_use = set(good_antpos).union(sorted_unique_ants)
# ants_use are the antennas we'll keep track of in the UVData
# object, so they dictate Nants_telescope
self.Nants_telescope = len(ants_use)
self.antenna_numbers = np.array(list(ants_use))
self.antenna_positions = np.zeros((self.Nants_telescope, 3), dtype=rel_ecef_antpos.dtype)
for ai, num in enumerate(self.antenna_numbers):
if antpos_length[num] == 0:
warnings.warn('antenna number {n} has visibilities '
'associated with it, but it has a position'
' of (0,0,0)'.format(n=num))
else:
# leave bad locations as zeros to make them obvious
self.antenna_positions[ai, :] = rel_ecef_antpos[num, :]
except(KeyError):
# there is no antpos variable
warnings.warn('Antenna positions are not present in the file.')
self.antenna_positions = None
if self.antenna_numbers is None:
# there are no antenna_numbers or antenna_positions, so just use
# the antennas present in the visibilities
# (Nants_data will therefore match Nants_telescope)
self.antenna_numbers = np.array(sorted_unique_ants)
self.Nants_telescope = len(self.antenna_numbers)
# antenna names is a foreign concept in miriad but required in other formats.
try:
# Here we deal with the way pyuvdata tacks it on to keep the
# name information if we have it:
# make it into one long comma-separated string
ant_name_var = uv['antnames']
if isinstance(ant_name_var, str):
ant_name_str = ant_name_var.replace('\x00', '')
ant_name_list = ant_name_str[1:-1].split(', ')
self.antenna_names = ant_name_list
else:
# Backwards compatibility for old way of storing antenna_names.
# This is a horrible hack to save & recover antenna_names array.
# Miriad can't handle arrays of strings and AIPY use to not handle
# long enough single strings to put them all into one string
# so we convert them into hex values and then into floats on
# write and convert back to strings here
warnings.warn('{file} was written with an old version of '
'pyuvdata, which has been deprecated. Rewrite this '
'file with write_miriad to ensure future '
'compatibility'.format(file=filepath))
ant_name_flt = uv['antnames']
ant_name_list = [('%x' % elem.astype(np.int64)).decode('hex') for elem in ant_name_flt]
self.antenna_names = ant_name_list
except(KeyError):
self.antenna_names = self.antenna_numbers.astype(str).tolist()
# check for antenna diameters
try:
self.antenna_diameters = uv['antdiam']
except(KeyError):
# backwards compatibility for when keyword was 'diameter'
try:
self.antenna_diameters = uv['diameter']
# if we find it, we need to remove it from extra_keywords to keep from writing it out
self.extra_keywords.pop('diameter')
except(KeyError):
pass
if self.antenna_diameters is not None:
self.antenna_diameters = (self.antenna_diameters
* np.ones(self.Nants_telescope, dtype=np.float))
# form up a grid which indexes time and baselines along the 'long'
# axis of the visdata array
tij_grid = np.array(map(lambda x: map(float, x.split("_")), unique_blts))
t_grid, ant_i_grid, ant_j_grid = tij_grid.T
# set the data sizes
if antenna_nums is None and ant_pairs_nums is None and ant_str is None and time_range is None:
try:
self.Nblts = uv['nblts']
if self.Nblts != len(t_grid):
warnings.warn('Nblts does not match the number of unique blts in the data')
self.Nblts = len(t_grid)
except(KeyError):
self.Nblts = len(t_grid)
else:
# The select on read will make the header nblts not match the number of unique blts
self.Nblts = len(t_grid)
if time_range is None:
try:
self.Ntimes = uv['ntimes']
if self.Ntimes != len(times):
warnings.warn('Ntimes does not match the number of unique times in the data')
self.Ntimes = len(times)
except(KeyError):
self.Ntimes = len(times)
else:
# The select on read will make the header ntimes not match the number of unique times
self.Ntimes = len(times)
self.time_array = t_grid
self.ant_1_array = ant_i_grid.astype(int)
self.ant_2_array = ant_j_grid.astype(int)
self.baseline_array = self.antnums_to_baseline(ant_i_grid.astype(int),
ant_j_grid.astype(int))
if antenna_nums is None and ant_pairs_nums is None and ant_str is None:
try:
self.Nbls = uv['nbls']
if self.Nbls != len(np.unique(self.baseline_array)):
warnings.warn('Nbls does not match the number of unique baselines in the data')
self.Nbls = len(np.unique(self.baseline_array))
except(KeyError):
self.Nbls = len(np.unique(self.baseline_array))
else:
# The select on read will make the header nbls not match the number of unique bls
self.Nbls = len(np.unique(self.baseline_array))
# slot the data into a grid
self.data_array = np.zeros((self.Nblts, self.Nspws, self.Nfreqs,
self.Npols), dtype=np.complex64)
self.flag_array = np.ones(self.data_array.shape, dtype=np.bool)
self.uvw_array = np.zeros((self.Nblts, 3))
# NOTE: Using our lst calculator, which uses astropy,
# instead of aipy values which come from pyephem.
# The differences are of order 5 seconds.
if self.telescope_location is not None:
self.set_lsts_from_time_array()
self.nsample_array = np.ones(self.data_array.shape, dtype=np.float)
self.freq_array = (np.arange(self.Nfreqs) * self.channel_width
+ uv['sfreq'] * 1e9)
# Tile freq_array to shape (Nspws, Nfreqs).
# Currently does not actually support Nspws>1!
self.freq_array = np.tile(self.freq_array, (self.Nspws, 1))
# Temporary arrays to hold polarization axis, which will be collapsed
ra_pol_list = np.zeros((self.Nblts, self.Npols))
dec_pol_list = np.zeros((self.Nblts, self.Npols))
uvw_pol_list = np.zeros((self.Nblts, 3, self.Npols))
c_ns = const.c.to('m/ns').value
for pol, data in data_accumulator.iteritems():
pol_ind = self._pol_to_ind(pol)
for ind, d in enumerate(data):
blt = ["{1:.{0}f}".format(prec_t, d[1]).zfill(ndig_t),
str(d[2]).zfill(ndig_ant), str(d[3]).zfill(ndig_ant)]
blt = "_".join(blt)
blt_index = reverse_inds[blt]
self.data_array[blt_index, :, :, pol_ind] = d[4]
self.flag_array[blt_index, :, :, pol_ind] = d[5]
self.nsample_array[blt_index, :, :, pol_ind] = d[6]
# because there are uvws/ra/dec for each pol, and one pol may not
# have that visibility, we collapse along the polarization
# axis but avoid any missing visbilities
uvw = d[0] * c_ns
uvw.shape = (1, 3)
uvw_pol_list[blt_index, :, pol_ind] = uvw
ra_pol_list[blt_index, pol_ind] = d[7]
dec_pol_list[blt_index, pol_ind] = d[8]
# Collapse pol axis for ra_list, dec_list, and uvw_list
ra_list = np.zeros(self.Nblts)
dec_list = np.zeros(self.Nblts)
for blt_index in xrange(self.Nblts):
test = ~np.all(self.flag_array[blt_index, :, :, :], axis=(0, 1))
good_pol = np.where(test)[0]
if len(good_pol) == 1:
# Only one good pol, use it
self.uvw_array[blt_index, :] = uvw_pol_list[blt_index, :, good_pol]
ra_list[blt_index] = ra_pol_list[blt_index, good_pol]
dec_list[blt_index] = dec_pol_list[blt_index, good_pol]
elif len(good_pol) > 1:
# Multiple good pols, check for consistency. pyuvdata does not
# support pol-dependent uvw, ra, or dec.
if np.any(np.diff(uvw_pol_list[blt_index, :, good_pol], axis=0)):
raise ValueError('uvw values are different by polarization.')
else:
self.uvw_array[blt_index, :] = uvw_pol_list[blt_index, :, good_pol[0]]
if np.any(np.diff(ra_pol_list[blt_index, good_pol])):
raise ValueError('ra values are different by polarization.')
else:
ra_list[blt_index] = ra_pol_list[blt_index, good_pol[0]]
if np.any(np.diff(dec_pol_list[blt_index, good_pol])):
raise ValueError('dec values are different by polarization.')
else:
dec_list[blt_index] = dec_pol_list[blt_index, good_pol[0]]
else:
# No good pols for this blt. Fill with first one.
self.uvw_array[blt_index, :] = uvw_pol_list[blt_index, :, 0]
ra_list[blt_index] = ra_pol_list[blt_index, 0]
dec_list[blt_index] = dec_pol_list[blt_index, 0]
# get unflagged blts
blt_good = np.where(~np.all(self.flag_array, axis=(1, 2, 3)))
single_ra = np.isclose(np.mean(np.diff(ra_list[blt_good])), 0.)
single_time = np.isclose(np.mean(np.diff(self.time_array[blt_good])), 0.)
# first check to see if the phase_type was specified.
if phase_type is not None:
if phase_type is 'phased':
self.set_phased()
elif phase_type is 'drift':
self.set_drift()
else:
raise ValueError('The phase_type was not recognized. '
'Set the phase_type to "drift" or "phased" to '
'reflect the phasing status of the data')
else:
# check if ra is constant throughout file; if it is,
# file is tracking if not, file is drift scanning
# check if there's only one unflagged time
if not single_time:
if single_ra:
self.set_phased()
else:
self.set_drift()
else:
# if there's only one time, checking for consistent RAs doesn't work.
# instead check for the presence of an epoch variable, which isn't
# really a good option, but at least it prevents crashes.
if 'epoch' in uv.vartable.keys():
self.set_phased()
else:
self.set_drift()
if self.phase_type == 'phased':
# check that the RA values do not vary
if not single_ra:
raise(ValueError, 'phase_type is "phased" but the RA values are varying.')
self.phase_center_ra = float(ra_list[0])
self.phase_center_dec = float(dec_list[0])
self.phase_center_epoch = uv['epoch']
else:
# check that the RA values are not constant (if more than one time present)
if (single_ra and not single_time):
raise(ValueError, 'phase_type is "drift" but the RA values are constant.')
self.zenith_ra = ra_list
self.zenith_dec = dec_list
try:
self.set_telescope_params()
except ValueError, ve:
warnings.warn(str(ve))
# check if object has all required uv_properties set
if run_check:
self.check(check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
def write_miriad(self, filepath, run_check=True, check_extra=True,
run_check_acceptability=True,
clobber=False, no_antnums=False):
"""
Write the data to a miriad file.
Args:
filename: The miriad file directory to write to.
run_check: Option to check for the existence and proper shapes of
parameters before writing 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 before writing the file. Default is True.
clobber: Option to overwrite the filename if the file already exists.
Default is False.
no_antnums: Option to not write the antnums variable to the file.
Should only be used for testing purposes.
"""
if run_check:
self.check(check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
# check for multiple spws
if self.data_array.shape[1] > 1:
raise ValueError('write_miriad currently only handles single spw files.')
if os.path.exists(filepath):
if clobber:
print 'File exists: clobbering'
shutil.rmtree(filepath)
else:
raise ValueError('File exists: skipping')
if self.Nfreqs > 1:
freq_spacing = self.freq_array[0, 1:] - self.freq_array[0, :-1]
if not np.isclose(np.min(freq_spacing), np.max(freq_spacing),
rtol=self._freq_array.tols[0], atol=self._freq_array.tols[1]):
raise ValueError('The frequencies are not evenly spaced (probably '
'because of a select operation). The miriad format '
'does not support unevenly spaced frequencies.')
if not np.isclose(np.max(freq_spacing), self.channel_width,
rtol=self._freq_array.tols[0], atol=self._freq_array.tols[1]):
raise ValueError('The frequencies are separated by more than their '
'channel width (probably because of a select operation). '
'The miriad format does not support frequencies '
'that are spaced by more than their channel width.')
uv = aipy.miriad.UV(filepath, status='new')
# initialize header variables
uv._wrhd('obstype', 'mixed-auto-cross')
# avoid inserting extra \n.
if not self.history[-1] == '\n':
self.history += '\n'
uv._wrhd('history', self.history)
# recognized miriad variables
uv.add_var('nchan', 'i')
uv['nchan'] = self.Nfreqs
uv.add_var('npol', 'i')
uv['npol'] = self.Npols
uv.add_var('nspect', 'i')
uv['nspect'] = self.Nspws
uv.add_var('inttime', 'd')
uv['inttime'] = self.integration_time
uv.add_var('sdf', 'd')
uv['sdf'] = self.channel_width / 1e9 # in GHz
uv.add_var('source', 'a')
uv['source'] = self.object_name
uv.add_var('telescop', 'a')
uv['telescop'] = self.telescope_name
uv.add_var('latitud', 'd')
uv['latitud'] = self.telescope_location_lat_lon_alt[0]
uv.add_var('longitu', 'd')
uv['longitu'] = self.telescope_location_lat_lon_alt[1]
uv.add_var('nants', 'i')
if self.x_orientation is not None:
uv.add_var('xorient', 'a')
uv['xorient'] = self.x_orientation
if self.antenna_diameters is not None:
if not np.allclose(self.antenna_diameters, self.antenna_diameters[0]):
warnings.warn('Antenna diameters are not uniform, but miriad only'
'supports a single diameter. Skipping.')
else:
uv.add_var('antdiam', 'd')
uv['antdiam'] = float(self.antenna_diameters[0])
# These are added to make files written by pyuvdata more "miriad correct", and
# should be changed when support for more than one spectral window is added.
# 'nschan' is the number of channels per spectral window, and 'ischan' is the
# starting channel for each spectral window. Both should be arrays of size Nspws.
# Also note that indexing in Miriad is 1-based
uv.add_var('nschan', 'i')
uv['nschan'] = self.Nfreqs
uv.add_var('ischan', 'i')
uv['ischan'] = 1
# Miriad has no way to keep track of antenna numbers, so the antenna
# numbers are simply the index for each antenna in any array that
# describes antenna attributes (e.g. antpos for the antenna_postions).
# Therefore on write, nants (which gives the size of the antpos array)
# needs to be increased to be the max value of antenna_numbers+1 and the
# antpos array needs to be inflated with zeros at locations where we
# don't have antenna information. These inflations need to be undone at
# read. If the file was written by pyuvdata, then the variable antnums
# will be present and we can use it, otherwise we need to test for zeros
# in the antpos array and/or antennas with no visibilities.
nants = np.max(self.antenna_numbers) + 1
uv['nants'] = nants
if self.antenna_positions is not None:
# Miriad wants antenna_positions to be in absolute coordinates
# (not relative to array center) in a rotated ECEF frame where the
# x-axis goes through the local meridian.
rel_ecef_antpos = np.zeros((nants, 3), dtype=self.antenna_positions.dtype)
for ai, num in enumerate(self.antenna_numbers):
rel_ecef_antpos[num, :] = self.antenna_positions[ai, :]
# find zeros so antpos can be zeroed there too
antpos_length = np.sqrt(np.sum(np.abs(rel_ecef_antpos)**2, axis=1))
ecef_antpos = rel_ecef_antpos + self.telescope_location
longitude = self.telescope_location_lat_lon_alt[1]
antpos = uvutils.rotECEF_from_ECEF(ecef_antpos, longitude)
# zero out bad locations (these are checked on read)
antpos[np.where(antpos_length == 0), :] = [0, 0, 0]
uv.add_var('antpos', 'd')
# Miriad stores antpos values in units of ns, pyuvdata uses meters.
uv['antpos'] = antpos.T.flatten() / const.c.to('m/ns').value
uv.add_var('sfreq', 'd')
uv['sfreq'] = self.freq_array[0, 0] / 1e9 # first spw; in GHz
if self.phase_type == 'phased':
uv.add_var('epoch', 'r')
uv['epoch'] = self.phase_center_epoch
# required pyuvdata variables that are not recognized miriad variables
uv.add_var('ntimes', 'i')
uv['ntimes'] = self.Ntimes
uv.add_var('nbls', 'i')
uv['nbls'] = self.Nbls
uv.add_var('nblts', 'i')
uv['nblts'] = self.Nblts
uv.add_var('visunits', 'a')
uv['visunits'] = self.vis_units
uv.add_var('instrume', 'a')
uv['instrume'] = self.instrument
uv.add_var('altitude', 'd')
uv['altitude'] = self.telescope_location_lat_lon_alt[2]
# optional pyuvdata variables that are not recognized miriad variables
if self.dut1 is not None:
uv.add_var('dut1', 'd')
uv['dut1'] = self.dut1
if self.earth_omega is not None:
uv.add_var('degpdy', 'd')
uv['degpdy'] = self.earth_omega
if self.gst0 is not None:
uv.add_var('gst0', 'd')
uv['gst0'] = self.gst0
if self.rdate is not None:
uv.add_var('rdate', 'a')
uv['rdate'] = self.rdate
if self.timesys is not None:
uv.add_var('timesys', 'a')
uv['timesys'] = self.timesys
# other extra keywords
# set up dictionaries to map common python types to miriad types
# NB: arrays/lists/dicts could potentially be written as strings or 1D
# vectors. This is not supported at present!
# NB: complex numbers *should* be supportable, but are not currently
# supported due to unexplained errors in aipy.miriad and/or its underlying libraries
numpy_types = {np.int8: int,
np.int16: int,
np.int32: int,
np.int64: int,
np.uint8: int,
np.uint16: int,
np.uint32: int,
np.uint64: int,
np.float16: float,
np.float32: float,
np.float64: float,
np.float128: float,