-
Notifications
You must be signed in to change notification settings - Fork 25
/
uvdata.py
2627 lines (2351 loc) · 130 KB
/
uvdata.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
# -*- coding: utf-8 -*-
"""Primary container for radio interferometer datasets.
"""
from __future__ import absolute_import, division, print_function
from astropy import constants as const
from astropy.time import Time
import os
import numpy as np
import six
import warnings
import ephem
from .uvbase import UVBase
from . import parameter as uvp
from . import telescopes as uvtel
from . import utils as uvutils
import copy
import collections
import re
class UVData(UVBase):
"""
A class for defining a radio interferometer dataset.
Currently supported file types: uvfits, miriad, fhd.
Provides phasing functions.
Attributes:
UVParameter objects: For full list see UVData Parameters
(http://pyuvdata.readthedocs.io/en/latest/uvdata_parameters.html).
Some are always required, some are required for certain phase_types
and others are always optional.
"""
def __init__(self):
"""Create a new UVData object."""
# add the UVParameters to the class
# standard angle tolerance: 10 mas in radians.
# Should perhaps be decreased to 1 mas in the future
radian_tol = 10 * 2 * np.pi * 1e-3 / (60.0 * 60.0 * 360.0)
self._Ntimes = uvp.UVParameter('Ntimes', description='Number of times',
expected_type=int)
self._Nbls = uvp.UVParameter('Nbls', description='Number of baselines',
expected_type=int)
self._Nblts = uvp.UVParameter('Nblts', description='Number of baseline-times '
'(i.e. number of spectra). Not necessarily '
'equal to Nbls * Ntimes', expected_type=int)
self._Nfreqs = uvp.UVParameter('Nfreqs', description='Number of frequency channels',
expected_type=int)
self._Npols = uvp.UVParameter('Npols', description='Number of polarizations',
expected_type=int)
desc = ('Array of the visibility data, shape: (Nblts, Nspws, Nfreqs, '
'Npols), type = complex float, in units of self.vis_units')
self._data_array = uvp.UVParameter('data_array', description=desc,
form=('Nblts', 'Nspws',
'Nfreqs', 'Npols'),
expected_type=np.complex)
desc = 'Visibility units, options are: "uncalib", "Jy" or "K str"'
self._vis_units = uvp.UVParameter('vis_units', description=desc,
form='str', expected_type=str,
acceptable_vals=["uncalib", "Jy", "K str"])
desc = ('Number of data points averaged into each data element, '
'NOT required to be an integer. type = float, same shape as data_array')
self._nsample_array = uvp.UVParameter('nsample_array', description=desc,
form=('Nblts', 'Nspws',
'Nfreqs', 'Npols'),
expected_type=(np.float))
desc = 'Boolean flag, True is flagged, same shape as data_array.'
self._flag_array = uvp.UVParameter('flag_array', description=desc,
form=('Nblts', 'Nspws',
'Nfreqs', 'Npols'),
expected_type=np.bool)
self._Nspws = uvp.UVParameter('Nspws', description='Number of spectral windows '
'(ie non-contiguous spectral chunks). '
'More than one spectral window is not '
'currently supported.', expected_type=int)
self._spw_array = uvp.UVParameter('spw_array',
description='Array of spectral window '
'Numbers, shape (Nspws)', form=('Nspws',),
expected_type=int)
desc = ('Projected baseline vectors relative to phase center, '
'shape (Nblts, 3), units meters. Convention is: uvw = xyz(ant2) - xyz(ant1).'
'Note that this is the Miriad convention but it is different '
'from the AIPS/FITS convention (where uvw = xyz(ant1) - xyz(ant2)).')
self._uvw_array = uvp.UVParameter('uvw_array', description=desc,
form=('Nblts', 3),
expected_type=np.float,
acceptable_range=(1e-3, 1e8), tols=.001)
desc = ('Array of times, center of integration, shape (Nblts), '
'units Julian Date')
self._time_array = uvp.UVParameter('time_array', description=desc,
form=('Nblts',),
expected_type=np.float,
tols=1e-3 / (60.0 * 60.0 * 24.0)) # 1 ms in days
desc = ('Array of lsts, center of integration, shape (Nblts), '
'units radians')
self._lst_array = uvp.UVParameter('lst_array', description=desc,
form=('Nblts',),
expected_type=np.float,
tols=radian_tol)
desc = ('Array of first antenna indices, shape (Nblts), '
'type = int, 0 indexed')
self._ant_1_array = uvp.UVParameter('ant_1_array', description=desc,
expected_type=int, form=('Nblts',))
desc = ('Array of second antenna indices, shape (Nblts), '
'type = int, 0 indexed')
self._ant_2_array = uvp.UVParameter('ant_2_array', description=desc,
expected_type=int, form=('Nblts',))
desc = ('Array of baseline indices, shape (Nblts), '
'type = int; baseline = 2048 * (ant1+1) + (ant2+1) + 2^16')
self._baseline_array = uvp.UVParameter('baseline_array',
description=desc,
expected_type=int, form=('Nblts',))
# this dimensionality of freq_array does not allow for different spws
# to have different dimensions
desc = 'Array of frequencies, shape (Nspws, Nfreqs), units Hz'
self._freq_array = uvp.UVParameter('freq_array', description=desc,
form=('Nspws', 'Nfreqs'),
expected_type=np.float,
tols=1e-3) # mHz
desc = ('Array of polarization integers, shape (Npols). '
'AIPS Memo 117 says: pseudo-stokes 1:4 (pI, pQ, pU, pV); '
'circular -1:-4 (RR, LL, RL, LR); linear -5:-8 (XX, YY, XY, YX). '
'NOTE: AIPS Memo 117 actually calls the pseudo-Stokes polarizations '
'"Stokes", but this is inaccurate as visibilities cannot be in '
'true Stokes polarizations for physical antennas. We adopt the '
'term pseudo-Stokes to refer to linear combinations of instrumental '
'visibility polarizations (e.g. pI = xx + yy).')
self._polarization_array = uvp.UVParameter('polarization_array',
description=desc,
expected_type=int,
acceptable_vals=list(
np.arange(-8, 0)) + list(np.arange(1, 5)),
form=('Npols',))
self._integration_time = uvp.UVParameter('integration_time',
description='Length of the integration (s)',
expected_type=np.float, tols=1e-3) # 1 ms
self._channel_width = uvp.UVParameter('channel_width',
description='Width of frequency channels (Hz)',
expected_type=np.float,
tols=1e-3) # 1 mHz
# --- observation information ---
self._object_name = uvp.UVParameter('object_name',
description='Source or field '
'observed (string)', form='str',
expected_type=str)
self._telescope_name = uvp.UVParameter('telescope_name',
description='Name of telescope '
'(string)', form='str',
expected_type=str)
self._instrument = uvp.UVParameter('instrument', description='Receiver or backend. '
'Sometimes identical to telescope_name',
form='str', expected_type=str)
desc = ('Telescope location: xyz in ITRF (earth-centered frame). '
'Can also be accessed using telescope_location_lat_lon_alt or '
'telescope_location_lat_lon_alt_degrees properties')
self._telescope_location = uvp.LocationParameter('telescope_location',
description=desc,
acceptable_range=(
6.35e6, 6.39e6),
tols=1e-3)
self._history = uvp.UVParameter('history', description='String of history, units English',
form='str', expected_type=str)
# --- phasing information ---
desc = ('String indicating phasing type. Allowed values are "drift", '
'"phased" and "unknown"')
self._phase_type = uvp.UVParameter('phase_type', form='str', expected_type=str,
description=desc, value='unknown',
acceptable_vals=['drift', 'phased', 'unknown'])
desc = ('Required if phase_type = "drift". Right ascension of zenith. '
'units: radians, shape (Nblts). Can also be accessed using zenith_ra_degrees.')
self._zenith_ra = uvp.AngleParameter('zenith_ra', required=False,
description=desc,
expected_type=np.float,
form=('Nblts',),
tols=radian_tol)
desc = ('Required if phase_type = "drift". Declination of zenith. '
'units: radians, shape (Nblts). Can also be accessed using zenith_dec_degrees.')
# in practice, dec of zenith will never change; does not need to be shape Nblts
self._zenith_dec = uvp.AngleParameter('zenith_dec', required=False,
description=desc,
expected_type=np.float,
form=('Nblts',),
tols=radian_tol)
desc = ('Required if phase_type = "phased". Epoch year of the phase '
'applied to the data (eg 2000.)')
self._phase_center_epoch = uvp.UVParameter('phase_center_epoch',
required=False,
description=desc,
expected_type=np.float)
desc = ('Required if phase_type = "phased". Right ascension of phase '
'center (see uvw_array), units radians. Can also be accessed using phase_center_ra_degrees.')
self._phase_center_ra = uvp.AngleParameter('phase_center_ra',
required=False,
description=desc,
expected_type=np.float,
tols=radian_tol)
desc = ('Required if phase_type = "phased". Declination of phase center '
'(see uvw_array), units radians. Can also be accessed using phase_center_dec_degrees.')
self._phase_center_dec = uvp.AngleParameter('phase_center_dec',
required=False,
description=desc,
expected_type=np.float,
tols=radian_tol)
# --- antenna information ----
desc = ('Number of antennas with data present (i.e. number of unique '
'entries in ant_1_array and ant_2_array). May be smaller '
'than the number of antennas in the array')
self._Nants_data = uvp.UVParameter('Nants_data', description=desc,
expected_type=int)
desc = ('Number of antennas in the array. May be larger '
'than the number of antennas with data')
self._Nants_telescope = uvp.UVParameter('Nants_telescope',
description=desc, expected_type=int)
desc = ('List of antenna names, shape (Nants_telescope), '
'with numbers given by antenna_numbers (which can be matched '
'to ant_1_array and ant_2_array). There must be one entry '
'here for each unique entry in ant_1_array and '
'ant_2_array, but there may be extras as well.')
self._antenna_names = uvp.UVParameter('antenna_names', description=desc,
form=('Nants_telescope',),
expected_type=str)
desc = ('List of integer antenna numbers corresponding to antenna_names, '
'shape (Nants_telescope). There must be one '
'entry here for each unique entry in ant_1_array and '
'ant_2_array, but there may be extras as well.')
self._antenna_numbers = uvp.UVParameter('antenna_numbers', description=desc,
form=('Nants_telescope',),
expected_type=int)
# -------- extra, non-required parameters ----------
desc = ('Orientation of the physical dipole corresponding to what is '
'labelled as the x polarization. Examples include "east" '
'(indicating east/west orientation) and "north" (indicating '
'north/south orientation)')
self._x_orientation = uvp.UVParameter('x_orientation', description=desc,
required=False, expected_type=str)
desc = ('Any user supplied extra keywords, type=dict. Keys should be '
'8 character or less strings if writing to uvfits or miriad files. '
'Use the special key "comment" for long multi-line string comments.')
self._extra_keywords = uvp.UVParameter('extra_keywords', required=False,
description=desc, value={},
spoof_val={}, expected_type=dict)
desc = ('Array giving coordinates of antennas relative to '
'telescope_location (ITRF frame), shape (Nants_telescope, 3), '
'units meters. See the tutorial page in the documentation '
'for an example of how to convert this to topocentric frame.')
self._antenna_positions = uvp.AntPositionParameter('antenna_positions',
required=False,
description=desc,
form=(
'Nants_telescope', 3),
expected_type=np.float,
tols=1e-3) # 1 mm
desc = ('Array of antenna diameters in meters. Used by CASA to '
'construct a default beam if no beam is supplied.')
self._antenna_diameters = uvp.UVParameter('antenna_diameters',
required=False,
description=desc,
form=('Nants_telescope',),
expected_type=np.float,
tols=1e-3) # 1 mm
# --- other stuff ---
# the below are copied from AIPS memo 117, but could be revised to
# merge with other sources of data.
self._gst0 = uvp.UVParameter('gst0', required=False,
description='Greenwich sidereal time at '
'midnight on reference date',
spoof_val=0.0, expected_type=np.float)
self._rdate = uvp.UVParameter('rdate', required=False,
description='Date for which the GST0 or '
'whatever... applies',
spoof_val='', form='str')
self._earth_omega = uvp.UVParameter('earth_omega', required=False,
description='Earth\'s rotation rate '
'in degrees per day',
spoof_val=360.985, expected_type=np.float)
self._dut1 = uvp.UVParameter('dut1', required=False,
description='DUT1 (google it) AIPS 117 '
'calls it UT1UTC',
spoof_val=0.0, expected_type=np.float)
self._timesys = uvp.UVParameter('timesys', required=False,
description='We only support UTC',
spoof_val='UTC', form='str')
desc = ('FHD thing we do not understand, something about the time '
'at which the phase center is normal to the chosen UV plane '
'for phasing')
self._uvplane_reference_time = uvp.UVParameter('uvplane_reference_time',
required=False,
description=desc,
spoof_val=0)
super(UVData, self).__init__()
def check(self, check_extra=True, run_check_acceptability=True):
"""
Add some extra checks on top of checks on UVBase class.
Check that required parameters exist. Check that parameters have
appropriate shapes and optionally that the values are acceptable.
Args:
check_extra: If true, check all parameters, otherwise only check
required parameters.
run_check_acceptability: Option to check if values in parameters
are acceptable. Default is True.
"""
# first run the basic check from UVBase
if self.ant_1_array is not None and np.all(self.ant_1_array == self.ant_2_array):
# Special case of only containing auto correlations, adjust uvw acceptable_range
self._uvw_array.acceptable_range = (0.0, 0.0)
# set the phase type based on object's value
if self.phase_type == 'phased':
self.set_phased()
elif self.phase_type == 'drift':
self.set_drift()
else:
self.set_unknown_phase_type()
super(UVData, self).check(check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
# Check internal consistency of numbers which don't explicitly correspond
# to the shape of another array.
nants_data_calc = int(len(np.unique(self.ant_1_array.tolist()
+ self.ant_2_array.tolist())))
if self.Nants_data != nants_data_calc:
raise ValueError('Nants_data must be equal to the number of unique '
'values in ant_1_array and ant_2_array')
if self.Nbls != len(np.unique(self.baseline_array)):
raise ValueError('Nbls must be equal to the number of unique '
'baselines in the data_array')
if self.Ntimes != len(np.unique(self.time_array)):
raise ValueError('Ntimes must be equal to the number of unique '
'times in the time_array')
# issue warning if extra_keywords keys are longer than 8 characters
for key in self.extra_keywords.keys():
if len(key) > 8:
warnings.warn('key {key} in extra_keywords is longer than 8 '
'characters. It will be truncated to 8 if written '
'to uvfits or miriad file formats.'.format(key=key))
# issue warning if extra_keywords values are lists, arrays or dicts
for key, value in self.extra_keywords.items():
if isinstance(value, (list, dict, np.ndarray)):
warnings.warn('{key} in extra_keywords is a list, array or dict, '
'which will raise an error when writing uvfits or '
'miriad file types'.format(key=key))
return True
def set_drift(self):
"""Set phase_type to 'drift' and adjust required parameters."""
self.phase_type = 'drift'
self._zenith_ra.required = True
self._zenith_dec.required = True
self._phase_center_epoch.required = False
self._phase_center_ra.required = False
self._phase_center_dec.required = False
def set_phased(self):
"""Set phase_type to 'phased' and adjust required parameters."""
self.phase_type = 'phased'
self._zenith_ra.required = False
self._zenith_dec.required = False
self._phase_center_epoch.required = True
self._phase_center_ra.required = True
self._phase_center_dec.required = True
def set_unknown_phase_type(self):
"""Set phase_type to 'unknown' and adjust required parameters."""
self.phase_type = 'unknown'
self._zenith_ra.required = False
self._zenith_dec.required = False
self._phase_center_epoch.required = False
self._phase_center_ra.required = False
self._phase_center_dec.required = False
def known_telescopes(self):
"""
Retun a list of telescopes known to pyuvdata.
This is just a shortcut to uvdata.telescopes.known_telescopes()
"""
return uvtel.known_telescopes()
def set_telescope_params(self, overwrite=False):
"""
Set telescope related parameters.
If the telescope_name is in the known_telescopes, set any missing
telescope-associated parameters (e.g. telescope location) to the value
for the known telescope.
Args:
overwrite: Option to overwrite existing telescope-associated
parameters with the values from the known telescope.
Default is False.
"""
telescope_obj = uvtel.get_telescope(self.telescope_name)
if telescope_obj is not False:
params_set = []
for p in telescope_obj:
telescope_param = getattr(telescope_obj, p)
self_param = getattr(self, p)
if telescope_param.value is not None and (overwrite is True
or self_param.value is None):
telescope_shape = telescope_param.expected_shape(telescope_obj)
self_shape = self_param.expected_shape(self)
if telescope_shape == self_shape:
params_set.append(self_param.name)
prop_name = self_param.name
setattr(self, prop_name, getattr(telescope_obj, prop_name))
else:
# expected shapes aren't equal. This can happen e.g. with diameters,
# which is a single value on the telescope object but is
# an array of length Nants_telescope on the UVData object
if telescope_shape == () and self_shape != 'str':
array_val = np.zeros(self_shape,
dtype=telescope_param.expected_type) + telescope_param.value
params_set.append(self_param.name)
prop_name = self_param.name
setattr(self, prop_name, array_val)
else:
raise ValueError('parameter {p} on the telescope '
'object does not have a compatible '
'expected shape.')
if len(params_set) > 0:
params_set_str = ', '.join(params_set)
warnings.warn('{params} is not set. Using known values '
'for {telescope_name}.'.format(params=params_set_str,
telescope_name=telescope_obj.telescope_name))
else:
raise ValueError('Telescope {telescope_name} is not in '
'known_telescopes.'.format(telescope_name=self.telescope_name))
def baseline_to_antnums(self, baseline):
"""
Get the antenna numbers corresponding to a given baseline number.
Args:
baseline: integer baseline number
Returns:
tuple with the two antenna numbers corresponding to the baseline.
"""
if self.Nants_telescope > 2048:
raise Exception('error Nants={Nants}>2048 not '
'supported'.format(Nants=self.Nants_telescope))
if np.min(baseline) > 2**16:
ant2 = (baseline - 2**16) % 2048 - 1
ant1 = (baseline - 2**16 - (ant2 + 1)) / 2048 - 1
else:
ant2 = (baseline) % 256 - 1
ant1 = (baseline - (ant2 + 1)) / 256 - 1
return np.int32(ant1), np.int32(ant2)
def antnums_to_baseline(self, ant1, ant2, attempt256=False):
"""
Get the baseline number corresponding to two given antenna numbers.
Args:
ant1: first antenna number (integer)
ant2: second antenna number (integer)
attempt256: Option to try to use the older 256 standard used in
many uvfits files (will use 2048 standard if there are more
than 256 antennas). Default is False.
Returns:
integer baseline number corresponding to the two antenna numbers.
"""
ant1, ant2 = np.int64((ant1, ant2))
if self.Nants_telescope is not None and self.Nants_telescope > 2048:
raise Exception('cannot convert ant1, ant2 to a baseline index '
'with Nants={Nants}>2048.'
.format(Nants=self.Nants_telescope))
if attempt256:
if (np.max(ant1) < 255 and np.max(ant2) < 255):
return 256 * (ant1 + 1) + (ant2 + 1)
else:
print('Max antnums are {} and {}'.format(
np.max(ant1), np.max(ant2)))
message = 'antnums_to_baseline: found > 256 antennas, using ' \
'2048 baseline indexing. Beware compatibility ' \
'with CASA etc'
warnings.warn(message)
baseline = 2048 * (ant1 + 1) + (ant2 + 1) + 2**16
if isinstance(baseline, np.ndarray):
return np.asarray(baseline, dtype=np.int64)
else:
return np.int64(baseline)
def order_pols(self, order='AIPS'):
'''
Arranges polarizations in orders corresponding to either AIPS or CASA convention.
CASA stokes types are ordered with cross-pols in between (e.g. XX,XY,YX,YY) while
AIPS orders pols with auto-pols followed by cross-pols (e.g. XX,YY,XY,YX)
Args:
order: string, either 'CASA' or 'AIPS'. Default='AIPS'
'''
if(order == 'AIPS'):
pol_inds = np.argsort(self.polarization_array)
pol_inds = pol_inds[::-1]
elif(order == 'CASA'): # sandwich
casa_order = np.array([1, 2, 3, 4, -1, -3, -4, -2, -5, -7, -8, -6])
pol_inds = []
for pol in self.polarization_array:
pol_inds.append(np.where(casa_order == pol)[0][0])
pol_inds = np.argsort(pol_inds)
else:
warnings.warn('Invalid order supplied. No sorting performed.')
pol_inds = list(range(self.Npols))
# Generate a map from original indices to new indices
if not np.array_equal(pol_inds, self.Npols):
self.reorder_pols(order=pol_inds)
def set_lsts_from_time_array(self):
"""Set the lst_array based from the time_array."""
lsts = []
self.lst_array = np.zeros(self.Nblts)
latitude, longitude, altitude = self.telescope_location_lat_lon_alt_degrees
for ind, jd in enumerate(np.unique(self.time_array)):
t = Time(jd, format='jd', location=(longitude, latitude))
self.lst_array[np.where(np.isclose(
jd, self.time_array, atol=1e-6, rtol=1e-12))] = t.sidereal_time('apparent').radian
def juldate2ephem(self, num):
"""
Convert Julian date to ephem date, measured from noon, Dec. 31, 1899.
Args:
num: Julian date
Returns:
ephem date, measured from noon, Dec. 31, 1899.
"""
return ephem.date(num - 2415020.)
def ephem2juldate(self, ephemdate):
"""
Convert ephem date to Julian date
Args:
ephemdate: ephemeral date, as ephem object or tuple (year, dd, hh, mm, ss)
Returns:
Julian date
"""
return ephem.date(ephemdate) + 2415020.
def unphase_to_drift(self):
"""Convert from a phased dataset to a drift dataset."""
if self.phase_type == 'phased':
pass
elif self.phase_type == 'drift':
raise ValueError('The data is already drift scanning; can only '
'unphase phased data.')
else:
raise ValueError('The phasing type of the data is unknown. '
'Set the phase_type to drift or phased to '
'reflect the phasing status of the data')
latitude, longitude, altitude = self.telescope_location_lat_lon_alt
obs = ephem.Observer()
# obs inits with default values for parameters -- be sure to replace them
obs.lat = latitude
obs.lon = longitude
phase_center = ephem.FixedBody()
epoch = (self.phase_center_epoch - 2000.) * 365.2422 + \
ephem.J2000 # convert years to ephemtime
phase_center._epoch = epoch
phase_center._ra = self.phase_center_ra
phase_center._dec = self.phase_center_dec
self.zenith_ra = np.zeros_like(self.time_array)
self.zenith_dec = np.zeros_like(self.time_array)
# apply -w phasor
w_lambda = (self.uvw_array[:, 2].reshape(self.Nblts, 1).astype(np.float64)
/ const.c.to('m/s').value * self.freq_array.reshape(1, self.Nfreqs))
phs = np.exp(-1j * 2 * np.pi * (-1) * w_lambda[:, None, :, None])
self.data_array *= phs
unique_times = np.unique(self.time_array)
for jd in unique_times:
inds = np.where(self.time_array == jd)[0]
obs.date, obs.epoch = self.juldate2ephem(
jd), self.juldate2ephem(jd)
phase_center.compute(obs)
phase_center_ra, phase_center_dec = phase_center.a_ra, phase_center.a_dec
zenith_ra = obs.sidereal_time()
zenith_dec = latitude
self.zenith_ra[inds] = zenith_ra
self.zenith_dec[inds] = zenith_dec
# generate rotation matrices
m0 = uvutils.top2eq_m(0., phase_center_dec)
m1 = uvutils.eq2top_m(phase_center_ra - zenith_ra, zenith_dec)
# rotate and write uvws
uvw = self.uvw_array[inds, :]
uvw = np.dot(m0, uvw.T).T
uvw = np.dot(m1, uvw.T).T
self.uvw_array[inds, :] = uvw
# remove phase center
self.phase_center_ra = None
self.phase_center_dec = None
self.phase_center_epoch = None
self.set_drift()
def phase_to_time(self, time):
"""
Phase a drift scan dataset to the ra/dec of zenith at a particular time.
Args:
time: The time to phase to.
"""
if self.phase_type == 'drift':
pass
elif self.phase_type == 'phased':
raise ValueError('The data is already phased; can only phase '
'drift scanning data.')
else:
raise ValueError('The phasing type of the data is unknown. '
'Set the phase_type to drift or phased to '
'reflect the phasing status of the data')
obs = ephem.Observer()
# obs inits with default values for parameters -- be sure to replace them
latitude, longitude, altitude = self.telescope_location_lat_lon_alt
obs.lat = latitude
obs.lon = longitude
obs.date, obs.epoch = self.juldate2ephem(
time), self.juldate2ephem(time)
ra = obs.sidereal_time()
dec = latitude
epoch = self.juldate2ephem(time)
self.phase(ra, dec, epoch)
def phase(self, ra, dec, epoch):
"""
Phase a drift scan dataset to a single ra/dec at a particular epoch.
Will not phase already phased data.
Args:
ra: The ra to phase to in radians.
dec: The dec to phase to in radians.
epoch: The epoch to use for phasing. Should be an ephem date,
measured from noon Dec. 31, 1899.
"""
if self.phase_type == 'drift':
pass
elif self.phase_type == 'phased':
raise ValueError('The data is already phased; can only phase '
'drift scan data. Use unphase_to_drift to '
'convert to a drift scan.')
else:
raise ValueError('The phasing type of the data is unknown. '
'Set the phase_type to "drift" or "phased" to '
'reflect the phasing status of the data')
obs = ephem.Observer()
# obs inits with default values for parameters -- be sure to replace them
latitude, longitude, altitude = self.telescope_location_lat_lon_alt
obs.lat = latitude
obs.lon = longitude
# create a pyephem object for the phasing position
precess_pos = ephem.FixedBody()
precess_pos._ra = ra
precess_pos._dec = dec
precess_pos._epoch = epoch
# calculate RA/DEC in J2000 and write to object
obs.date, obs.epoch = ephem.J2000, ephem.J2000
precess_pos.compute(obs)
self.phase_center_ra = precess_pos.a_ra + \
0.0 # force to be a float not ephem.Angle
self.phase_center_dec = precess_pos.a_dec + \
0.0 # force to be a float not ephem.Angle
# explicitly set epoch to J2000
self.phase_center_epoch = 2000.0
unique_times, unique_inds = np.unique(
self.time_array, return_index=True)
uvws = np.zeros(self.uvw_array.shape, dtype=np.float64)
for ind, jd in enumerate(unique_times):
inds = np.where(self.time_array == jd)[0]
lst = self.lst_array[unique_inds[ind]]
# calculate ra/dec of phase center in current epoch
obs.date, obs.epoch = self.juldate2ephem(
jd), self.juldate2ephem(jd)
precess_pos.compute(obs)
ra, dec = precess_pos.a_ra, precess_pos.a_dec
# generate rotation matrices
m0 = uvutils.top2eq_m(lst - obs.sidereal_time(), latitude)
m1 = uvutils.eq2top_m(lst - ra, dec)
# rotate and write uvws
uvw = self.uvw_array[inds, :]
uvw = np.dot(m0, uvw.T).T
uvw = np.dot(m1, uvw.T).T
self.uvw_array[inds, :] = uvw
uvws[inds, :] = uvw
# calculate data and apply phasor
w_lambda = (uvws[:, 2].reshape(self.Nblts, 1)
/ const.c.to('m/s').value * self.freq_array.reshape(1, self.Nfreqs))
phs = np.exp(-1j * 2 * np.pi * w_lambda[:, None, :, None])
self.data_array *= phs
del(obs)
self.set_phased()
def __add__(self, other, run_check=True, check_extra=True,
run_check_acceptability=True, inplace=False):
"""
Combine two UVData objects. Objects can be added along frequency,
polarization, and/or baseline-time axis.
Args:
other: Another UVData object which will be added to self.
run_check: Option to check for the existence and proper shapes of
parameters after combining objects. 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 combining objects. Default is True.
inplace: Overwrite self as we go, otherwise create a third object
as the sum of the two (default).
"""
if inplace:
this = self
else:
this = copy.deepcopy(self)
# Check that both objects are UVData and valid
this.check(check_extra=check_extra, run_check_acceptability=run_check_acceptability)
if not isinstance(other, this.__class__):
raise ValueError('Only UVData objects can be added to a UVData object')
other.check(check_extra=check_extra, run_check_acceptability=run_check_acceptability)
# Check objects are compatible
# Note zenith_ra will not necessarily be the same if times are different.
# But phase_center should be the same, even if in drift (empty parameters)
compatibility_params = ['_vis_units', '_integration_time', '_channel_width',
'_object_name', '_telescope_name', '_instrument',
'_telescope_location', '_phase_type',
'_Nants_telescope', '_antenna_names',
'_antenna_numbers', '_antenna_positions',
'_phase_center_ra', '_phase_center_dec',
'_phase_center_epoch']
for a in compatibility_params:
if getattr(this, a) != getattr(other, a):
msg = 'UVParameter ' + \
a[1:] + ' does not match. Cannot combine objects.'
raise ValueError(msg)
# Build up history string
history_update_string = ' Combined data along '
n_axes = 0
# Create blt arrays for convenience
prec_t = - 2 * \
np.floor(np.log10(this._time_array.tols[-1])).astype(int)
prec_b = 8
this_blts = np.array(["_".join(["{1:.{0}f}".format(prec_t, blt[0]),
str(blt[1]).zfill(prec_b)]) for blt in
zip(this.time_array, this.baseline_array)])
other_blts = np.array(["_".join(["{1:.{0}f}".format(prec_t, blt[0]),
str(blt[1]).zfill(prec_b)]) for blt in
zip(other.time_array, other.baseline_array)])
# Check we don't have overlapping data
both_pol = np.intersect1d(
this.polarization_array, other.polarization_array)
both_freq = np.intersect1d(
this.freq_array[0, :], other.freq_array[0, :])
both_blts = np.intersect1d(this_blts, other_blts)
if len(both_pol) > 0:
if len(both_freq) > 0:
if len(both_blts) > 0:
raise ValueError('These objects have overlapping data and'
' cannot be combined.')
temp = np.nonzero(~np.in1d(other_blts, this_blts))[0]
if len(temp) > 0:
bnew_inds = temp
new_blts = other_blts[temp]
history_update_string += 'baseline-time'
n_axes += 1
else:
bnew_inds, new_blts = ([], [])
temp = np.nonzero(
~np.in1d(other.freq_array[0, :], this.freq_array[0, :]))[0]
if len(temp) > 0:
fnew_inds = temp
new_freqs = other.freq_array[0, temp]
if n_axes > 0:
history_update_string += ', frequency'
else:
history_update_string += 'frequency'
n_axes += 1
else:
fnew_inds, new_freqs = ([], [])
temp = np.nonzero(~np.in1d(other.polarization_array,
this.polarization_array))[0]
if len(temp) > 0:
pnew_inds = temp
new_pols = other.polarization_array[temp]
if n_axes > 0:
history_update_string += ', polarization'
else:
history_update_string += 'polarization'
n_axes += 1
else:
pnew_inds, new_pols = ([], [])
# Pad out self to accommodate new data
if len(bnew_inds) > 0:
this_blts = np.concatenate((this_blts, new_blts))
blt_order = np.argsort(this_blts)
zero_pad = np.zeros(
(len(bnew_inds), this.Nspws, this.Nfreqs, this.Npols))
this.data_array = np.concatenate([this.data_array, zero_pad], axis=0)
this.nsample_array = np.concatenate([this.nsample_array, zero_pad], axis=0)
this.flag_array = np.concatenate([this.flag_array,
1 - zero_pad], axis=0).astype(np.bool)
this.uvw_array = np.concatenate([this.uvw_array,
other.uvw_array[bnew_inds, :]], axis=0)[blt_order, :]
this.time_array = np.concatenate([this.time_array,
other.time_array[bnew_inds]])[blt_order]
this.lst_array = np.concatenate(
[this.lst_array, other.lst_array[bnew_inds]])[blt_order]
this.ant_1_array = np.concatenate([this.ant_1_array,
other.ant_1_array[bnew_inds]])[blt_order]
this.ant_2_array = np.concatenate([this.ant_2_array,
other.ant_2_array[bnew_inds]])[blt_order]
this.baseline_array = np.concatenate([this.baseline_array,
other.baseline_array[bnew_inds]])[blt_order]
if this.phase_type == 'drift':
this.zenith_ra = np.concatenate([this.zenith_ra,
other.zenith_ra[bnew_inds]])[blt_order]
this.zenith_dec = np.concatenate([this.zenith_dec,
other.zenith_dec[bnew_inds]])[blt_order]
if len(fnew_inds) > 0:
zero_pad = np.zeros((this.data_array.shape[0], this.Nspws, len(fnew_inds),
this.Npols))
this.freq_array = np.concatenate([this.freq_array,
other.freq_array[:, fnew_inds]], axis=1)
f_order = np.argsort(this.freq_array[0, :])
this.data_array = np.concatenate([this.data_array, zero_pad], axis=2)
this.nsample_array = np.concatenate([this.nsample_array, zero_pad], axis=2)
this.flag_array = np.concatenate([this.flag_array, 1 - zero_pad],
axis=2).astype(np.bool)
if len(pnew_inds) > 0:
zero_pad = np.zeros((this.data_array.shape[0], this.Nspws,
this.data_array.shape[2], len(pnew_inds)))
this.polarization_array = np.concatenate([this.polarization_array,
other.polarization_array[pnew_inds]])
p_order = np.argsort(np.abs(this.polarization_array))
this.data_array = np.concatenate([this.data_array, zero_pad], axis=3)
this.nsample_array = np.concatenate([this.nsample_array, zero_pad], axis=3)
this.flag_array = np.concatenate([this.flag_array, 1 - zero_pad],
axis=3).astype(np.bool)
# Now populate the data
pol_t2o = np.nonzero(
np.in1d(this.polarization_array, other.polarization_array))[0]
freq_t2o = np.nonzero(
np.in1d(this.freq_array[0, :], other.freq_array[0, :]))[0]
blt_t2o = np.nonzero(np.in1d(this_blts, other_blts))[0]
this.data_array[np.ix_(blt_t2o, [0], freq_t2o,
pol_t2o)] = other.data_array
this.nsample_array[np.ix_(
blt_t2o, [0], freq_t2o, pol_t2o)] = other.nsample_array
this.flag_array[np.ix_(blt_t2o, [0], freq_t2o,
pol_t2o)] = other.flag_array
if len(bnew_inds) > 0:
this.data_array = this.data_array[blt_order, :, :, :]
this.nsample_array = this.nsample_array[blt_order, :, :, :]
this.flag_array = this.flag_array[blt_order, :, :, :]
if len(fnew_inds) > 0:
this.freq_array = this.freq_array[:, f_order]
this.data_array = this.data_array[:, :, f_order, :]
this.nsample_array = this.nsample_array[:, :, f_order, :]
this.flag_array = this.flag_array[:, :, f_order, :]
if len(pnew_inds) > 0:
this.polarization_array = this.polarization_array[p_order]
this.data_array = this.data_array[:, :, :, p_order]
this.nsample_array = this.nsample_array[:, :, :, p_order]
this.flag_array = this.flag_array[:, :, :, p_order]
# Update N parameters (e.g. Npols)
this.Ntimes = len(np.unique(this.time_array))
this.Nbls = len(np.unique(this.baseline_array))
this.Nblts = this.uvw_array.shape[0]
this.Nfreqs = this.freq_array.shape[1]
this.Npols = this.polarization_array.shape[0]
this.Nants_data = len(
np.unique(this.ant_1_array.tolist() + this.ant_2_array.tolist()))
# Check specific requirements
if this.Nfreqs > 1:
freq_separation = np.diff(this.freq_array[0, :])
if not np.isclose(np.min(freq_separation), np.max(freq_separation),
rtol=this._freq_array.tols[0], atol=this._freq_array.tols[1]):
warnings.warn('Combined frequencies are not evenly spaced. This will '
'make it impossible to write this data out to some file types.')
elif np.max(freq_separation) > this.channel_width:
warnings.warn('Combined frequencies are not contiguous. This will make '
'it impossible to write this data out to some file types.')
if this.Npols > 2:
pol_separation = np.diff(this.polarization_array)
if np.min(pol_separation) < np.max(pol_separation):
warnings.warn('Combined polarizations are not evenly spaced. This will '
'make it impossible to write this data out to some file types.')
if n_axes > 0:
history_update_string += ' axis using pyuvdata.'
this.history += history_update_string
this.history = uvutils.combine_histories(this.history, other.history)
# Check final object is self-consistent
if run_check:
this.check(check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
if not inplace:
return this
def __iadd__(self, other):
"""
In place add.
Args:
other: Another UVData object which will be added to self.
"""
self.__add__(other, inplace=True)
return self
def _select_preprocess(self, antenna_nums, antenna_names, ant_str, bls,
frequencies, freq_chans, times, polarizations, blt_inds):
"""
Internal function to build up blt_inds, freq_inds, pol_inds
and history_update_string for select.
"""