https://github.com/RadioAstronomySoftwareGroup/pyuvdata
Tip revision: c2cd6c296ccb7aa31f3333b27193baf5f547682a authored by Nicholas Kern on 04 June 2018, 15:34:20 UTC
moved combine_uvdata to uvdata.py and updated tests
moved combine_uvdata to uvdata.py and updated tests
Tip revision: c2cd6c2
uvdata.py
"""Primary container for radio interferometer datasets."""
from astropy import constants as const
from astropy.time import Time
import os
import numpy as np
import warnings
import ephem
from uvbase import UVBase
import parameter as uvp
import telescopes as uvtel
import utils as uvutils
import copy
import collections
import re
import itertools
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.iteritems():
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 StandardError('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 > 2048:
raise StandardError('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)
return np.int64(2048 * (ant1 + 1) + (ant2 + 1) + 2**16)
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 = 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 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).
"""
# combine
combined = combine_uvdata([self, other], inplace=inplace, run_check=run_check, check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
if not inplace:
return combined
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, ant_pairs_nums,
frequencies, freq_chans, times, polarizations, blt_inds):
"""
Internal function to build up blt_inds, freq_inds, pol_inds
and history_update_string for select.
"""
# build up history string as we go
history_update_string = ' Downselected to specific '
n_selects = 0
if ant_str is not None:
if not (antenna_nums is None and antenna_names is None
and ant_pairs_nums is None and polarizations is None):
raise ValueError(
'Cannot provide ant_str with antenna_nums, antenna_names, '
'ant_pairs_nums, or polarizations.')
else:
ant_pairs_nums, polarizations = self.parse_ants(ant_str)
# Antennas, times and blt_inds all need to be combined into a set of
# blts indices to keep.
# test for blt_inds presence before adding inds from antennas & times
if blt_inds is not None:
blt_inds = uvutils.get_iterable(blt_inds)
history_update_string += 'baseline-times'
n_selects += 1
if antenna_names is not None:
if antenna_nums is not None:
raise ValueError(
'Only one of antenna_nums and antenna_names can be provided.')
antenna_names = uvutils.get_iterable(antenna_names)
antenna_nums = []
for s in antenna_names:
if s not in self.antenna_names:
raise ValueError(
'Antenna name {a} is not present in the antenna_names array'.format(a=s))
antenna_nums.append(self.antenna_numbers[np.where(
np.array(self.antenna_names) == s)[0]])
if antenna_nums is not None:
antenna_nums = uvutils.get_iterable(antenna_nums)
if n_selects > 0:
history_update_string += ', antennas'
else:
history_update_string += 'antennas'
n_selects += 1
inds1 = np.zeros(0, dtype=np.int)
inds2 = np.zeros(0, dtype=np.int)
for ant in antenna_nums:
if ant in self.ant_1_array or ant in self.ant_2_array:
wh1 = np.where(self.ant_1_array == ant)[0]
wh2 = np.where(self.ant_2_array == ant)[0]
if len(wh1) > 0:
inds1 = np.append(inds1, list(wh1))
if len(wh2) > 0:
inds2 = np.append(inds2, list(wh2))
else:
raise ValueError('Antenna number {a} is not present in the '
'ant_1_array or ant_2_array'.format(a=ant))
ant_blt_inds = np.array(
list(set(inds1).intersection(inds2)), dtype=np.int)
else:
ant_blt_inds = None
if ant_pairs_nums is not None:
if isinstance(ant_pairs_nums, tuple) and len(ant_pairs_nums) == 2:
ant_pairs_nums = [ant_pairs_nums]
if not all(isinstance(item, tuple) for item in ant_pairs_nums):
raise ValueError(
'ant_pairs_nums must be a list of tuples of antenna numbers.')
if not all([isinstance(item[0], (int, long, np.integer)) for item in ant_pairs_nums]
+ [isinstance(item[1], (int, long, np.integer)) for item in ant_pairs_nums]):
raise ValueError(
'ant_pairs_nums must be a list of tuples of antenna numbers.')
if n_selects > 0:
history_update_string += ', antenna pairs'
else:
history_update_string += 'antenna pairs'
n_selects += 1
ant_pair_blt_inds = np.zeros(0, dtype=np.int)
for pair in ant_pairs_nums:
if not (pair[0] in self.ant_1_array or pair[0] in self.ant_2_array):
raise ValueError('Antenna number {a} is not present in the '
'ant_1_array or ant_2_array'.format(a=pair[0]))
if not (pair[1] in self.ant_1_array or pair[1] in self.ant_2_array):
raise ValueError('Antenna number {a} is not present in the '
'ant_1_array or ant_2_array'.format(a=pair[1]))
wh1 = np.where(np.logical_and(
self.ant_1_array == pair[0], self.ant_2_array == pair[1]))[0]
wh2 = np.where(np.logical_and(
self.ant_1_array == pair[1], self.ant_2_array == pair[0]))[0]
if len(wh1) > 0:
ant_pair_blt_inds = np.append(ant_pair_blt_inds, list(wh1))
if len(wh2) > 0:
ant_pair_blt_inds = np.append(ant_pair_blt_inds, list(wh2))
if len(wh1) == 0 and len(wh2) == 0:
raise ValueError('Antenna pair {p} does not have any data '
'associated with it.'.format(p=pair))
if ant_blt_inds is not None:
# Use union (or) to join antenna_names/nums & ant_pairs_nums
ant_blt_inds = np.array(list(set(ant_blt_inds).union(ant_pair_blt_inds)))
else:
ant_blt_inds = ant_pair_blt_inds
if ant_blt_inds is not None:
if blt_inds is not None:
# Use intesection (and) to join antenna_names/nums/ant_pairs_nums with blt_inds
# handled differently because of the time aspect (which is anded with antennas below)
blt_inds = np.array(
list(set(blt_inds).intersection(ant_blt_inds)), dtype=np.int)
else:
blt_inds = ant_blt_inds
if times is not None:
times = uvutils.get_iterable(times)
if n_selects > 0:
history_update_string += ', times'
else:
history_update_string += 'times'
n_selects += 1
time_blt_inds = np.zeros(0, dtype=np.int)
for jd in times:
if jd in self.time_array:
time_blt_inds = np.append(
time_blt_inds, np.where(self.time_array == jd)[0])
else:
raise ValueError(
'Time {t} is not present in the time_array'.format(t=jd))
if blt_inds is not None:
# Use intesection (and) to join antenna_names/nums/ant_pairs_nums/blt_inds with times
blt_inds = np.array(
list(set(blt_inds).intersection(time_blt_inds)), dtype=np.int)
else:
blt_inds = time_blt_inds
if blt_inds is not None:
if len(blt_inds) == 0:
raise ValueError(
'No baseline-times were found that match criteria')
if max(blt_inds) >= self.Nblts:
raise ValueError(
'blt_inds contains indices that are too large')
if min(blt_inds) < 0:
raise ValueError('blt_inds contains indices that are negative')
blt_inds = list(sorted(set(list(blt_inds))))
if freq_chans is not None:
freq_chans = uvutils.get_iterable(freq_chans)
if frequencies is None:
frequencies = self.freq_array[0, freq_chans]
else:
frequencies = uvutils.get_iterable(frequencies)
frequencies = np.sort(list(set(frequencies)
| set(self.freq_array[0, freq_chans])))
if frequencies is not None:
frequencies = uvutils.get_iterable(frequencies)
if n_selects > 0:
history_update_string += ', frequencies'
else:
history_update_string += 'frequencies'
n_selects += 1
freq_inds = np.zeros(0, dtype=np.int)
# this works because we only allow one SPW. This will have to be reworked when we support more.
freq_arr_use = self.freq_array[0, :]
for f in frequencies:
if f in freq_arr_use:
freq_inds = np.append(
freq_inds, np.where(freq_arr_use == f)[0])
else:
raise ValueError(
'Frequency {f} is not present in the freq_array'.format(f=f))
if len(frequencies) > 1:
freq_ind_separation = freq_inds[1:] - freq_inds[:-1]
if np.min(freq_ind_separation) < np.max(freq_ind_separation):
warnings.warn('Selected frequencies are not evenly spaced. This '
'will make it impossible to write this data out to '
'some file types')
elif np.max(freq_ind_separation) > 1:
warnings.warn('Selected frequencies are not contiguous. This '
'will make it impossible to write this data out to '
'some file types.')
freq_inds = list(sorted(set(list(freq_inds))))
else:
freq_inds = None
if polarizations is not None:
polarizations = uvutils.get_iterable(polarizations)
if n_selects > 0:
history_update_string += ', polarizations'
else:
history_update_string += 'polarizations'
n_selects += 1
pol_inds = np.zeros(0, dtype=np.int)
for p in polarizations:
if isinstance(p, str):
p_num = uvutils.polstr2num(p)
else:
p_num = p
if p_num in self.polarization_array:
pol_inds = np.append(pol_inds, np.where(
self.polarization_array == p_num)[0])
else:
raise ValueError(
'Polarization {p} is not present in the polarization_array'.format(p=p))
if len(pol_inds) > 2:
pol_ind_separation = pol_inds[1:] - pol_inds[:-1]
if np.min(pol_ind_separation) < np.max(pol_ind_separation):
warnings.warn('Selected polarization values are not evenly spaced. This '
'will make it impossible to write this data out to '
'some file types')
pol_inds = list(sorted(set(list(pol_inds))))
else:
pol_inds = None
history_update_string += ' using pyuvdata.'
return blt_inds, freq_inds, pol_inds, history_update_string
def _select_metadata(self, blt_inds, freq_inds, pol_inds, history_update_string):
"""
Internal function to perform select on everything except the data_array,
flag_array and nsample_array to allow for re-use in uvfits reading
"""
if blt_inds is not None:
self.Nblts = len(blt_inds)
self.baseline_array = self.baseline_array[blt_inds]
self.Nbls = len(np.unique(self.baseline_array))
self.time_array = self.time_array[blt_inds]
self.lst_array = self.lst_array[blt_inds]
self.uvw_array = self.uvw_array[blt_inds, :]
self.ant_1_array = self.ant_1_array[blt_inds]
self.ant_2_array = self.ant_2_array[blt_inds]
self.Nants_data = int(
len(set(self.ant_1_array.tolist() + self.ant_2_array.tolist())))
self.Ntimes = len(np.unique(self.time_array))
if self.phase_type == 'drift':
self.zenith_ra = self.zenith_ra[blt_inds]
self.zenith_dec = self.zenith_dec[blt_inds]
if freq_inds is not None:
self.Nfreqs = len(freq_inds)
self.freq_array = self.freq_array[:, freq_inds]
if pol_inds is not None:
self.Npols = len(pol_inds)
self.polarization_array = self.polarization_array[pol_inds]
self.history = self.history + history_update_string
def select(self, antenna_nums=None, antenna_names=None, ant_str=None,
ant_pairs_nums=None, frequencies=None, freq_chans=None,
times=None, polarizations=None, blt_inds=None, run_check=True,
check_extra=True, run_check_acceptability=True, inplace=True):
"""
Select specific antennas, antenna pairs, frequencies, times and
polarizations to keep in the object while discarding others.
Also supports selecting specific baseline-time indices to keep while
discarding others, but this is not commonly used. The history attribute
on the object will be updated to identify the operations performed.
Args:
antenna_nums: The antennas numbers to keep in the object (antenna
positions and names for the removed antennas will be retained).
This cannot be provided if antenna_names is also provided.
antenna_names: The antennas names to keep in the object (antenna
positions and names for the removed antennas will be retained).
This cannot be provided if antenna_nums is also provided.
ant_pairs_nums: A list of antenna number tuples (e.g. [(0,1), (3,2)])
specifying baselines to keep in the object. Ordering of the
numbers within the tuple does not matter.
ant_str: A string containing information about what antenna numbers
and polarizations to keep in the object. Can be 'auto', 'cross', 'all',
or combinations of antenna numbers and polarizations (e.g. '1',
'1_2', '1x_2y'). See tutorial for more examples of valid strings and
the behavior of different forms for ant_str.
If '1x_2y,2y_3y' is passed, both polarizations 'xy' and 'yy' will
be kept for both baselines (1,2) and (2,3) to return a valid
pyuvdata object.
An ant_str cannot be passed in addition to any of the above antenna
args or the polarizations arg. If this occurs, select raises a ValueError.
frequencies: The frequencies to keep in the object.
freq_chans: The frequency channel numbers to keep in the object.
times: The times to keep in the object.
polarizations: The polarizations to keep in the object.
blt_inds: The baseline-time indices to keep in the object. This is
not commonly used.
run_check: Option to check for the existence and proper shapes of
parameters after downselecting data on this object. 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 downselecting data on this object. Default is True.
inplace: Option to perform the select directly on self (True, default) or return
a new UVData object, which is a subselection of self (False)
"""
if inplace:
uv_object = self
else:
uv_object = copy.deepcopy(self)
blt_inds, freq_inds, pol_inds, history_update_string = \
uv_object._select_preprocess(antenna_nums, antenna_names, ant_str, ant_pairs_nums,
frequencies, freq_chans, times, polarizations, blt_inds)
# do select operations on everything except data_array, flag_array and nsample_array
uv_object._select_metadata(blt_inds, freq_inds, pol_inds, history_update_string)
if blt_inds is not None:
uv_object.data_array = uv_object.data_array[blt_inds, :, :, :]
uv_object.flag_array = uv_object.flag_array[blt_inds, :, :, :]
uv_object.nsample_array = uv_object.nsample_array[blt_inds, :, :, :]
if freq_inds is not None:
uv_object.data_array = uv_object.data_array[:, :, freq_inds, :]
uv_object.flag_array = uv_object.flag_array[:, :, freq_inds, :]
uv_object.nsample_array = uv_object.nsample_array[:, :, freq_inds, :]
if pol_inds is not None:
uv_object.data_array = uv_object.data_array[:, :, :, pol_inds]
uv_object.flag_array = uv_object.flag_array[:, :, :, pol_inds]
uv_object.nsample_array = uv_object.nsample_array[:, :, :, pol_inds]
# check if object is uv_object-consistent
if run_check:
uv_object.check(check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
if not inplace:
return uv_object
def _convert_from_filetype(self, other):
for p in other:
param = getattr(other, p)
setattr(self, p, param)
def _convert_to_filetype(self, filetype):
if filetype is 'uvfits':
import uvfits
other_obj = uvfits.UVFITS()
elif filetype is 'fhd':
import fhd
other_obj = fhd.FHD()
elif filetype is 'miriad':
import miriad
other_obj = miriad.Miriad()
elif filetype is 'uvh5':
import uvh5
other_obj = uvh5.UVH5()
else:
raise ValueError('filetype must be uvfits, miriad, fhd, or uvh5')
for p in self:
param = getattr(self, p)
setattr(other_obj, p, param)
return other_obj
def read_uvfits(self, filename, antenna_nums=None, antenna_names=None,
ant_str=None, ant_pairs_nums=None, frequencies=None,
freq_chans=None, times=None, polarizations=None, blt_inds=None,
read_data=True, read_metadata=True, run_check=True,
check_extra=True, run_check_acceptability=True):
"""
Read in header, metadata and data from uvfits file(s).
Args:
filename: The uvfits file or list of files to read from.
antenna_nums: The antennas numbers to include when reading data into
the object (antenna positions and names for the excluded antennas
will be retained). This cannot be provided if antenna_names is
also provided. Ignored if read_data is False.
antenna_names: The antennas names to include when reading data into
the object (antenna positions and names for the excluded antennas
will be retained). This cannot be provided if antenna_nums is
also provided. Ignored if read_data is False.
ant_pairs_nums: A list of antenna number tuples (e.g. [(0,1), (3,2)])
specifying baselines to include when reading data into the object.
Ordering of the numbers within the tuple does not matter.
Ignored if read_data is False.
ant_str: A string containing information about what antenna numbers
and polarizations to include when reading data into the object.
Can be 'auto', 'cross', 'all', or combinations of antenna numbers
and polarizations (e.g. '1', '1_2', '1x_2y').
See tutorial for more examples of valid strings and
the behavior of different forms for ant_str.
If '1x_2y,2y_3y' is passed, both polarizations 'xy' and 'yy' will
be kept for both baselines (1,2) and (2,3) to return a valid
pyuvdata object.
An ant_str cannot be passed in addition to any of the above antenna
args or the polarizations arg.
Ignored if read_data is False.
frequencies: The frequencies to include when reading data into the
object. Ignored if read_data is False.
freq_chans: The frequency channel numbers to include when reading
data into the object. Ignored if read_data is False.
times: The times to include when reading data into the object.
Ignored if read_data is False.
polarizations: The polarizations to include when reading data into
the object. Ignored if read_data is False.
blt_inds: The baseline-time indices to include when reading data into
the object. This is not commonly used. Ignored if read_data is False.
read_data: Read in the visibility and flag data. If set to false,
only the basic header info and metadata (if read_metadata is True)
will be read in. Results in an incompletely defined object
(check will not pass). Default True.
read_metadata: Read in metadata (times, baselines, uvws) as well as
basic header info. Only used if read_data is False
(metadata will be read if data is read). If both read_data and
read_metadata are false, only basic header info is read in. Default True.
run_check: Option to check for the existence and proper shapes of
parameters after reading in the file. Default is True.
Ignored if read_data is False.
check_extra: Option to check optional parameters as well as required
ones. Default is True. Ignored if read_data is False.
run_check_acceptability: Option to check acceptable range of the values of
parameters after reading in the file. Default is True.
Ignored if read_data is False.
"""
import uvfits
if isinstance(filename, (list, tuple)):
if not read_data:
raise ValueError('read_data cannot be False for a list of uvfits files')
self.read_uvfits(filename[0], antenna_nums=antenna_nums,
antenna_names=antenna_names, ant_str=ant_str,
ant_pairs_nums=ant_pairs_nums, frequencies=frequencies,
freq_chans=freq_chans, times=times,
polarizations=polarizations, blt_inds=blt_inds,
run_check=run_check, check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
if len(filename) > 1:
for f in filename[1:]:
uv2 = UVData()
uv2.read_uvfits(f, antenna_nums=antenna_nums,
antenna_names=antenna_names, ant_str=ant_str,
ant_pairs_nums=ant_pairs_nums, frequencies=frequencies,
freq_chans=freq_chans, times=times,
polarizations=polarizations, blt_inds=blt_inds,
run_check=run_check, check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
self += uv2
del(uv2)
else:
uvfits_obj = uvfits.UVFITS()
uvfits_obj.read_uvfits(filename, antenna_nums=antenna_nums,
antenna_names=antenna_names, ant_str=ant_str,
ant_pairs_nums=ant_pairs_nums, frequencies=frequencies,
freq_chans=freq_chans, times=times,
polarizations=polarizations, blt_inds=blt_inds,
read_data=read_data, read_metadata=read_metadata,
run_check=run_check, check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
self._convert_from_filetype(uvfits_obj)
del(uvfits_obj)
def read_uvfits_metadata(self, filename):
"""
Read in metadata (random parameter info) but not data from a uvfits file
(useful for an object that already has the associated header info and
full visibility data isn't needed).
Args:
filename: The uvfits file to read from.
"""
import uvfits
if isinstance(filename, (list, tuple)):
raise ValueError('A list of files cannot be used when just reading metadata')
uvfits_obj = self._convert_to_filetype('uvfits')
uvfits_obj.read_uvfits_metadata(filename)
self._convert_from_filetype(uvfits_obj)
del(uvfits_obj)
def read_uvfits_data(self, filename, antenna_nums=None, antenna_names=None,
ant_str=None, ant_pairs_nums=None, frequencies=None,
freq_chans=None, times=None, polarizations=None,
blt_inds=None, run_check=True, check_extra=True,
run_check_acceptability=True):
"""
Read in data but not header info from a uvfits file
(useful for an object that already has the associated header info).
Args:
filename: The uvfits file
antenna_nums: The antennas numbers to include when reading data into
the object (antenna positions and names for the excluded antennas
will be retained). This cannot be provided if antenna_names is
also provided.
antenna_names: The antennas names to include when reading data into
the object (antenna positions and names for the excluded antennas
will be retained). This cannot be provided if antenna_nums is
also provided.
ant_pairs_nums: A list of antenna number tuples (e.g. [(0,1), (3,2)])
specifying baselines to include when reading data into the object.
Ordering of the numbers within the tuple does not matter.
ant_str: A string containing information about what antenna numbers
and polarizations to include when reading data into the object.
Can be 'auto', 'cross', 'all', or combinations of antenna numbers
and polarizations (e.g. '1', '1_2', '1x_2y').
See tutorial for more examples of valid strings and
the behavior of different forms for ant_str.
If '1x_2y,2y_3y' is passed, both polarizations 'xy' and 'yy' will
be kept for both baselines (1,2) and (2,3) to return a valid
pyuvdata object.
An ant_str cannot be passed in addition to any of the above antenna
args or the polarizations arg.
frequencies: The frequencies to include when reading data into the
object.
freq_chans: The frequency channel numbers to include when reading
data into the object.
times: The times to include when reading data into the object.
polarizations: The polarizations to include when reading data into
the object.
blt_inds: The baseline-time indices to include when reading data into
the object. This is not commonly used.
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.
"""
import uvfits
if isinstance(filename, (list, tuple)):
raise ValueError('A list of files cannot be used when just reading data')
uvfits_obj = self._convert_to_filetype('uvfits')
uvfits_obj.read_uvfits_data(filename, antenna_nums=antenna_nums,
antenna_names=antenna_names, ant_str=ant_str,
ant_pairs_nums=ant_pairs_nums, frequencies=frequencies,
freq_chans=freq_chans, times=times,
polarizations=polarizations, blt_inds=blt_inds,
run_check=run_check, check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
self._convert_from_filetype(uvfits_obj)
del(uvfits_obj)
def write_uvfits(self, filename, spoof_nonessential=False,
force_phase=False, run_check=True, check_extra=True,
run_check_acceptability=True):
"""
Write the data to a uvfits file.
Args:
filename: The uvfits file to write to.
spoof_nonessential: Option to spoof the values of optional
UVParameters that are not set but are required for uvfits files.
Default is False.
force_phase: Option to automatically phase drift scan data to
zenith of the first timestamp. Default is False.
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.
"""
uvfits_obj = self._convert_to_filetype('uvfits')
uvfits_obj.write_uvfits(filename, spoof_nonessential=spoof_nonessential,
force_phase=force_phase, run_check=run_check,
check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
del(uvfits_obj)
def read_ms(self, filepath, run_check=True, check_extra=True,
run_check_acceptability=True, data_column='DATA', pol_order='AIPS'):
"""
Read in data from a measurement set
Args:
filepath: The measurement set file directory or list of directories to read from.
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 the values of parameters
after reading in the file. Default is True.
data_column: name of CASA data column to read into data_array.
'DATA', 'MODEL', or 'CORRECTED_DATA'
pol_order: specify whether you want polarizations ordered by
'CASA' or 'AIPS' conventions.
"""
# check if casacore is installed
try:
import casacore
except(ImportError):
# only import skip if importerror raised
# this way, nose import errors aren't raised
# unless this errors is already encountered.
from nose.plugins.skip import SkipTest
raise SkipTest('casacore not detected. casacore'
' must be installed for measurement set functionality')
import ms
if isinstance(filepath, (list, tuple)):
self.read_ms(filepath[0], run_check=run_check, check_extra=check_extra,
run_check_acceptability=run_check_acceptability,
data_column=data_column, pol_order=pol_order)
if len(filepath) > 1:
for f in filepath[1:]:
uv2 = UVData()
uv2.read_ms(f, run_check=run_check, check_extra=check_extra,
run_check_acceptability=run_check_acceptability,
data_column=data_column, pol_order=pol_order)
self += uv2
del(uv2)
else:
ms_obj = ms.MS()
ms_obj.read_ms(filepath, run_check=run_check, check_extra=check_extra,
run_check_acceptability=run_check_acceptability,
data_column=data_column, pol_order=pol_order)
self._convert_from_filetype(ms_obj)
del(ms_obj)
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.
Can also be a list of lists to read multiple data sets.
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.
"""
import fhd
if isinstance(filelist[0], (list, tuple)):
self.read_fhd(filelist[0], use_model=use_model, run_check=run_check,
check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
if len(filelist) > 1:
for f in filelist[1:]:
uv2 = UVData()
uv2.read_fhd(f, use_model=use_model, run_check=run_check,
check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
self += uv2
del(uv2)
else:
fhd_obj = fhd.FHD()
fhd_obj.read_fhd(filelist, use_model=use_model, run_check=run_check,
check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
self._convert_from_filetype(fhd_obj)
del(fhd_obj)
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 or list of directories to read from.
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 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]
"""
import miriad
if isinstance(filepath, (list, tuple)):
self.read_miriad(filepath[0], correct_lat_lon=correct_lat_lon,
run_check=run_check, check_extra=check_extra,
run_check_acceptability=run_check_acceptability,
phase_type=phase_type, antenna_nums=antenna_nums,
ant_str=ant_str, ant_pairs_nums=ant_pairs_nums,
polarizations=polarizations, time_range=time_range)
if len(filepath) > 1:
for f in filepath[1:]:
uv2 = UVData()
uv2.read_miriad(f, correct_lat_lon=correct_lat_lon,
run_check=run_check, check_extra=check_extra,
run_check_acceptability=run_check_acceptability,
phase_type=phase_type, antenna_nums=antenna_nums,
ant_str=ant_str, ant_pairs_nums=ant_pairs_nums,
polarizations=polarizations, time_range=time_range)
self += uv2
del(uv2)
else:
miriad_obj = miriad.Miriad()
miriad_obj.read_miriad(filepath, correct_lat_lon=correct_lat_lon,
run_check=run_check, check_extra=check_extra,
run_check_acceptability=run_check_acceptability,
phase_type=phase_type, antenna_nums=antenna_nums,
ant_str=ant_str, ant_pairs_nums=ant_pairs_nums,
polarizations=polarizations, time_range=time_range)
self._convert_from_filetype(miriad_obj)
del(miriad_obj)
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.
"""
miriad_obj = self._convert_to_filetype('miriad')
miriad_obj.write_miriad(filepath, run_check=run_check, check_extra=check_extra,
run_check_acceptability=run_check_acceptability,
clobber=clobber, no_antnums=no_antnums)
del(miriad_obj)
def read_uvh5(self, filename, run_check=True, check_extra=True,
run_check_acceptability=True):
"""
Read a UVH5 file.
Args:
filename: The UVH5 file to read.
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.
Returns:
None
"""
import uvh5
if isinstance(filename, (list, tuple)):
self.read_uvh5(filename[0], run_check=run_check, check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
if len(filename) > 1:
for f in filename[1:]:
uv2 = UVData()
uv2.read_uvh5(f, run_check=run_check, check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
self += uv2
del(uv2)
else:
uvh5_obj = uvh5.UVH5()
uvh5_obj.read_uvh5(filename, run_check=run_check, check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
self._convert_from_filetype(uvh5_obj)
del(uvh5_obj)
def write_uvh5(self, filename, run_check=True, check_extra=True,
run_check_acceptability=True, clobber=False):
"""
Write a UVData object to a UVH5 file.
Args:
filename: The UVH5 file 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 file if it already exists. Default is False.
Returns:
None
"""
uvh5_obj = self._convert_to_filetype('uvh5')
uvh5_obj.write_uvh5(filename, run_check=run_check,
check_extra=check_extra,
run_check_acceptability=run_check_acceptability,
clobber=clobber)
del(uvh5_obj)
def reorder_pols(self, order=None, run_check=True, check_extra=True,
run_check_acceptability=True):
"""
Rearrange polarizations in the event they are not uvfits compatible.
Args:
order: Provide the order which to shuffle the data. Default will
sort by absolute value of pol values.
run_check: Option to check for the existence and proper shapes of
parameters after reordering. 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 reordering. Default is True.
"""
if order is None:
order = np.argsort(np.abs(self.polarization_array))
self.polarization_array = self.polarization_array[order]
self.data_array = self.data_array[:, :, :, order]
self.nsample_array = self.nsample_array[:, :, :, order]
self.flag_array = self.flag_array[:, :, :, order]
# check if object is self-consistent
if run_check:
self.check(check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
def get_ants(self):
"""
Returns numpy array of unique antennas in the data.
"""
return np.unique(np.append(self.ant_1_array, self.ant_2_array))
def get_ENU_antpos(self, center=True, pick_data_ants=False):
"""
Returns antenna positions in ENU (topocentric) coordinates in units of meters.
Parameters:
-----------
center : bool, if True: subtract median of array position from antpos
pick_data_ants : bool, if True: return only antennas found in data
Returns: (antpos, ants)
--------
antpos : ndarray, antenna positions in TOPO frame and units of meters, shape=(Nants, 3)
ants : ndarray, antenna numbers matching ordering of antpos, shape=(Nants,)
"""
antpos = uvutils.ENU_from_ECEF((self.antenna_positions + self.telescope_location).T, *self.telescope_location_lat_lon_alt).T
ants = self.antenna_numbers
if pick_data_ants:
data_ants = np.unique(np.concatenate([self.ant_1_array, self.ant_2_array]))
telescope_ants = self.antenna_numbers
select = map(lambda x: x in data_ants, telescope_ants)
antpos = antpos[select, :]
ants = telescope_ants[select]
if center is True:
antpos -= np.median(antpos, axis=0)
return antpos, ants
def get_baseline_nums(self):
"""
Returns numpy array of unique baseline numbers in data.
"""
return np.unique(self.baseline_array)
def get_antpairs(self):
"""
Returns list of unique antpair tuples (ant1, ant2) in data.
"""
return [self.baseline_to_antnums(bl) for bl in self.get_baseline_nums()]
def get_pols(self):
"""
Returns list of pols in string format.
"""
return uvutils.polnum2str(self.polarization_array)
def get_antpairpols(self):
"""
Returns list of unique antpair + pol tuples (ant1, ant2, pol) in data.
"""
bli = 0
pols = self.get_pols()
bls = self.get_antpairs()
return [(bl) + (pol,) for bl in bls for pol in pols]
def get_feedpols(self):
"""
Returns list of unique antenna feed polarizations (e.g. ['X', 'Y']) in data.
NOTE: Will return ValueError if any pseudo-Stokes visibilities are present.
"""
if np.any(self.polarization_array > 0):
raise ValueError('Pseudo-Stokes visibilities cannot be interpreted as feed polarizations')
else:
return list(set(''.join(self.get_pols())))
def antpair2ind(self, ant1, ant2):
"""
Get blt indices for given (ordered) antenna pair.
"""
return np.where((self.ant_1_array == ant1) & (self.ant_2_array == ant2))[0]
def _key2inds(self, key):
"""
Interpret user specified key as a combination of antenna pair and/or polarization.
Args:
key: Identifier of data. Key can be 1, 2, or 3 numbers:
if len(key) == 1:
if (key < 5) or (type(key) is str): interpreted as a
polarization number/name, return all blts for that pol.
else: interpreted as a baseline number. Return all times and
polarizations for that baseline.
if len(key) == 2: interpreted as an antenna pair. Return all
times and pols for that baseline.
if len(key) == 3: interpreted as antenna pair and pol (ant1, ant2, pol).
Return all times for that baseline, pol. pol may be a string.
Returns:
blt_ind1: numpy array with blt indices for antenna pair.
blt_ind2: numpy array with blt indices for conjugate antenna pair.
pol_ind: numpy array with polarization indices
"""
key = uvutils.get_iterable(key)
if type(key) is str:
# Single string given, assume it is polarization
pol_ind = np.where(self.polarization_array == uvutils.polstr2num(key))[0]
if len(pol_ind) == 0:
raise KeyError('Polarization {pol} not found in data.'.format(pol=key))
blt_ind1 = np.arange(self.Nblts)
blt_ind2 = np.array([], dtype=np.int64)
elif len(key) == 1:
key = key[0] # For simplicity
if isinstance(key, collections.Iterable):
# Nested tuple. Call function again.
blt_ind1, blt_ind2, pol_ind = self._key2inds(key)
elif key < 5:
# Small number, assume it is a polarization number a la AIPS memo
pol_ind = np.where(self.polarization_array == key)[0]
if len(pol_ind) == 0:
raise KeyError('Polarization {pol} not found in data.'.format(pol=key))
blt_ind1 = np.arange(self.Nblts)
blt_ind2 = np.array([], dtype=np.int64)
else:
# Larger number, assume it is a baseline number
inv_bl = self.antnums_to_baseline(self.baseline_to_antnums(key)[1],
self.baseline_to_antnums(key)[0])
blt_ind1 = np.where(self.baseline_array == key)[0]
blt_ind2 = np.where(self.baseline_array == inv_bl)[0]
if len(blt_ind1) + len(blt_ind2) == 0:
raise KeyError('Baseline {bl} not found in data.'.format(bl=key))
pol_ind = np.arange(self.Npols)
elif len(key) == 2:
# Key is an antenna pair
blt_ind1 = self.antpair2ind(key[0], key[1])
blt_ind2 = self.antpair2ind(key[1], key[0])
if len(blt_ind1) + len(blt_ind2) == 0:
raise KeyError('Antenna pair {pair} not found in data'.format(pair=key))
pol_ind = np.arange(self.Npols)
elif len(key) == 3:
# Key is an antenna pair + pol
blt_ind1 = self.antpair2ind(key[0], key[1])
blt_ind2 = self.antpair2ind(key[1], key[0])
if len(blt_ind1) + len(blt_ind2) == 0:
raise KeyError('Antenna pair {pair} not found in '
'data'.format(pair=(key[0], key[1])))
if type(key[2]) is str:
# pol is str
pol_ind = np.where(self.polarization_array == uvutils.polstr2num(key[2]))[0]
else:
# polarization number a la AIPS memo
pol_ind = np.where(self.polarization_array == key[2])[0]
if len(pol_ind) == 0:
raise KeyError('Polarization {pol} not found in data.'.format(pol=key[2]))
# Catch autos
if np.array_equal(blt_ind1, blt_ind2):
blt_ind2 = np.array([])
return (blt_ind1, blt_ind2, pol_ind)
def _smart_slicing(self, data, ind1, ind2, indp, **kwargs):
"""
Method for quickly picking out the relevant section of data for get_data or get_flags
Args:
data: 4-dimensional array in the format of self.data_array
ind1: list with blt indices for antenna pair (e.g. from self._key2inds)
ind2: list with blt indices for conjugate antenna pair. (e.g. from self._key2inds)
indp: list with polarization indices (e.g. from self._key2inds)
kwargs:
squeeze: 'default': squeeze pol and spw dimensions if possible (default)
'none': no squeezing of resulting numpy array
'full': squeeze all length 1 dimensions
force_copy: Option to explicitly make a copy of the data. Default is False.
Returns:
out: numpy array copy (or if possible, a read-only view) of relevant section of data
"""
force_copy = kwargs.pop('force_copy', False)
squeeze = kwargs.pop('squeeze', 'default')
if len(set(np.ediff1d(indp))) <= 1:
p_reg_spaced = True
p_start = indp[0]
p_stop = indp[-1] + 1
if len(indp) == 1:
dp = 1
else:
dp = indp[1] - indp[0]
else:
p_reg_spaced = False
if len(ind2) == 0:
# only unconjugated baselines
if len(set(np.ediff1d(ind1))) <= 1:
blt_start = ind1[0]
blt_stop = ind1[-1] + 1
if len(ind1) == 1:
dblt = 1
else:
dblt = ind1[1] - ind1[0]
if p_reg_spaced:
out = data[blt_start:blt_stop:dblt, :, :, p_start:p_stop:dp]
else:
out = data[blt_start:blt_stop:dblt, :, :, indp]
else:
out = data[ind1, :, :, :]
if p_reg_spaced:
out = out[:, :, :, p_start:p_stop:dp]
else:
out = out[:, :, :, indp]
elif len(ind1) == 0:
# only conjugated baselines
if len(set(np.ediff1d(ind2))) <= 1:
blt_start = ind2[0]
blt_stop = ind2[-1] + 1
if len(ind2) == 1:
dblt = 1
else:
dblt = ind2[1] - ind2[0]
if p_reg_spaced:
out = np.conj(data[blt_start:blt_stop:dblt, :, :, p_start:p_stop:dp])
else:
out = np.conj(data[blt_start:blt_stop:dblt, :, :, indp])
else:
out = data[ind2, :, :, :]
if p_reg_spaced:
out = np.conj(out[:, :, :, p_start:p_stop:dp])
else:
out = np.conj(out[:, :, :, indp])
else:
# both conjugated and unconjugated baselines
out = np.append(data[ind1, :, :, :], np.conj(data[ind2, :, :, :]), axis=0)
if p_reg_spaced:
out = out[:, :, :, p_start:p_stop:dp]
else:
out = out[:, :, :, indp]
if squeeze is 'full':
out = np.squeeze(out)
elif squeeze is 'default':
if out.shape[3] is 1:
# one polarization dimension
out = np.squeeze(out, axis=3)
if out.shape[1] is 1:
# one spw dimension
out = np.squeeze(out, axis=1)
elif squeeze is not 'none':
raise ValueError('"' + str(squeeze) + '" is not a valid option for squeeze.'
'Only "default", "none", or "full" are allowed.')
if force_copy:
out = np.array(out)
elif out.base is not None:
# if out is a view rather than a copy, make it read-only
out.flags.writeable = False
return out
def get_data(self, *args, **kwargs):
"""
Function for quick access to numpy array with data corresponding to
a baseline and/or polarization. Returns a read-only view if possible, otherwise a copy.
Args:
*args: parameters or tuple of parameters defining the key to identify
desired data. See _key2inds for formatting.
**kwargs: Keyword arguments:
squeeze: 'default': squeeze pol and spw dimensions if possible
'none': no squeezing of resulting numpy array
'full': squeeze all length 1 dimensions
force_copy: Option to explicitly make a copy of the data.
Default is False.
Returns:
Numpy array of data corresponding to key.
If data exists conjugate to requested antenna pair, it will be conjugated
before returning.
"""
ind1, ind2, indp = self._key2inds(args)
out = self._smart_slicing(self.data_array, ind1, ind2, indp, **kwargs)
return out
def get_flags(self, *args, **kwargs):
"""
Function for quick access to numpy array with flags corresponding to
a baseline and/or polarization. Returns a read-only view if possible, otherwise a copy.
Args:
*args: parameters or tuple of parameters defining the key to identify
desired data. See _key2inds for formatting.
**kwargs: Keyword arguments:
squeeze: 'default': squeeze pol and spw dimensions if possible
'none': no squeezing of resulting numpy array
'full': squeeze all length 1 dimensions
force_copy: Option to explicitly make a copy of the data.
Default is False.
Returns:
Numpy array of flags corresponding to key.
"""
ind1, ind2, indp = self._key2inds(args)
out = self._smart_slicing(self.flag_array, ind1, ind2, indp, **kwargs)
return out
def get_nsamples(self, *args, **kwargs):
"""
Function for quick access to numpy array with nsamples corresponding to
a baseline and/or polarization. Returns a read-only view if possible, otherwise a copy.
Args:
*args: parameters or tuple of parameters defining the key to identify
desired data. See _key2inds for formatting.
**kwargs: Keyword arguments:
squeeze: 'default': squeeze pol and spw dimensions if possible
'none': no squeezing of resulting numpy array
'full': squeeze all length 1 dimensions
force_copy: Option to explicitly make a copy of the data.
Default is False.
Returns:
Numpy array of nsamples corresponding to key.
"""
ind1, ind2, indp = self._key2inds(args)
out = self._smart_slicing(self.nsample_array, ind1, ind2, indp, **kwargs)
return out
def get_times(self, *args):
"""
Find the time_array entries for a given antpair or baseline number.
Meant to be used in conjunction with get_data function.
Args:
*args: parameters or tuple of parameters defining the key to identify
desired data. See _key2inds for formatting.
Returns:
Numpy array of times corresonding to key.
"""
ind1, ind2, indp = self._key2inds(args)
return np.append(self.time_array[ind1], self.time_array[ind2])
def antpairpol_iter(self, squeeze='default'):
"""
Generates numpy arrays of data for each antpair, pol combination.
Args:
squeeze: 'default': squeeze pol and spw dimensions if possible
'none': no squeezing of resulting numpy array
'full': squeeze all length 1 dimensions
Returns (for each iteration):
key: tuple with antenna1, antenna2, and polarization string
data: Numpy array with data which is the result of self[key]
"""
antpairpols = self.get_antpairpols()
for key in antpairpols:
yield (key, self.get_data(key, squeeze=squeeze))
def parse_ants(self, ant_str, print_toggle=False):
"""
Generates two lists of antenna pair tuples and polarization indices based
on parsing of the string ant_str. If no valid polarizations (pseudo-Stokes
params, or combinations of [lr] or [xy]) or antenna numbers are found in
ant_str, ant_pairs_nums and polarizations are returned as None.
Args:
ant_str: String containing antenna information to pass to select
function. Can be 'all', 'auto', cross, or combinations of antenna
numbers and polarization indicators l and r or x and y. Minus
signs can also be used in front of an antenna number or baseline
to exclude it from being output in ant_pairs_nums. If the antenna
number attached with a minus sign is present in the outputted list
ant_pairs_nums, it will be removed from the list. If ant_str
passed with a minus sign as the first character, 'all,' will be
appended to the beginning of the string. See the
tutorial for examples of valid strings and their behavior.
print_toggle: Boolean for printing parsed baselines for a visual user
check.
Output:
ant_pairs_nums: List of tuples containing the parsed pairs of
antenna numbers. If 'all' or pseudo-Stokes parameters are passed as
ant_str, returned as None.
polarizations: List of desired polarizations. If no polarizations
found in ant_str then returned as None.
"""
ant_re = r'(\(((-?\d+[lrxy]?,?)+)\)|-?\d+[lrxy]?)'
bl_re = '(^(%s_%s|%s),?)' % (ant_re, ant_re, ant_re)
str_pos = 0
ant_pairs_nums = []
polarizations = []
ants_data = self.get_ants()
ant_pairs_data = self.get_antpairs()
pols_data = self.get_pols()
warned_ants = []
warned_pols = []
warnings.simplefilter('always')
if ant_str.startswith('-'):
ant_str = 'all,' + ant_str
while str_pos < len(ant_str):
m = re.search(bl_re, ant_str[str_pos:])
if m is None:
if ant_str[str_pos:].upper().startswith('ALL'):
if len(ant_str[str_pos:].split(',')) > 1:
ant_pairs_nums = self.get_antpairs()
elif ant_str[str_pos:].upper().startswith('AUTO'):
for pair in ant_pairs_data:
if (pair[0] == pair[1]
and pair not in ant_pairs_nums):
ant_pairs_nums.append(pair)
elif ant_str[str_pos:].upper().startswith('CROSS'):
for pair in ant_pairs_data:
if not (pair[0] == pair[1]
or pair in ant_pairs_nums):
ant_pairs_nums.append(pair)
elif ant_str[str_pos:].upper().startswith('PI'):
polarizations.append(uvutils.polstr2num('pI'))
elif ant_str[str_pos:].upper().startswith('PQ'):
polarizations.append(uvutils.polstr2num('pQ'))
elif ant_str[str_pos:].upper().startswith('PU'):
polarizations.append(uvutils.polstr2num('pU'))
elif ant_str[str_pos:].upper().startswith('PV'):
polarizations.append(uvutils.polstr2num('pV'))
else:
raise ValueError('Unparsible argument {s}'.format(s=ant_str))
comma_cnt = ant_str[str_pos:].find(',')
if comma_cnt >= 0:
str_pos += comma_cnt + 1
else:
str_pos = len(ant_str)
else:
m = m.groups()
str_pos += len(m[0])
if m[2] is None:
ant_i_list = [m[8]]
ant_j_list = list(self.get_ants())
else:
if m[3] is None:
ant_i_list = [m[2]]
else:
ant_i_list = m[3].split(',')
if m[6] is None:
ant_j_list = [m[5]]
else:
ant_j_list = m[6].split(',')
for ant_i in ant_i_list:
include_i = True
if type(ant_i) == str and ant_i.startswith('-'):
ant_i = ant_i[1:] # nibble the - off the string
include_i = False
for ant_j in ant_j_list:
include_j = True
if type(ant_j) == str and ant_j.startswith('-'):
ant_j = ant_j[1:]
include_j = False
pols = None
ant_i, ant_j = str(ant_i), str(ant_j)
if not ant_i.isdigit():
ai = re.search(r'(\d+)([x,y,l,r])', ant_i).groups()
if not ant_j.isdigit():
aj = re.search(r'(\d+)([x,y,l,r])', ant_j).groups()
if ant_i.isdigit() and ant_j.isdigit():
ai = [ant_i, '']
aj = [ant_j, '']
elif ant_i.isdigit() and not ant_j.isdigit():
if ('x' in ant_j or 'y' in ant_j):
pols = ['x' + aj[1], 'y' + aj[1]]
else:
pols = ['l' + aj[1], 'r' + aj[1]]
ai = [ant_i, '']
elif not ant_i.isdigit() and ant_j.isdigit():
if ('x' in ant_i or 'y' in ant_i):
pols = [ai[1] + 'x', ai[1] + 'y']
else:
pols = [ai[1] + 'l', ai[1] + 'r']
aj = [ant_j, '']
elif not ant_i.isdigit() and not ant_j.isdigit():
pols = [ai[1] + aj[1]]
ant_tuple = tuple((abs(int(ai[0])), abs(int(aj[0]))))
# Order tuple according to order in object
if ant_tuple in ant_pairs_data:
pass
elif ant_tuple[::-1] in ant_pairs_data:
ant_tuple = ant_tuple[::-1]
else:
if not (ant_tuple[0] in ants_data
or ant_tuple[0] in warned_ants):
warned_ants.append(ant_tuple[0])
if not (ant_tuple[1] in ants_data
or ant_tuple[1] in warned_ants):
warned_ants.append(ant_tuple[1])
if pols is not None:
for pol in pols:
if not (pol.upper() in pols_data
or pol in warned_pols):
warned_pols.append(pol)
continue
if include_i and include_j:
if ant_tuple not in ant_pairs_nums:
ant_pairs_nums.append(ant_tuple)
if pols is not None:
for pol in pols:
if (pol.upper() in pols_data
and uvutils.polstr2num(pol) not in polarizations):
polarizations.append(uvutils.polstr2num(pol))
elif not (pol.upper() in pols_data
or pol in warned_pols):
warned_pols.append(pol)
else:
if pols is not None:
for pol in pols:
if pol.upper() in pols_data:
if (self.Npols == 1
and [pol.upper()] == pols_data):
ant_pairs_nums.remove(ant_tuple)
if uvutils.polstr2num(pol) in polarizations:
polarizations.remove(uvutils.polstr2num(pol))
elif not (pol.upper() in pols_data
or pol in warned_pols):
warned_pols.append(pol)
elif ant_tuple in ant_pairs_nums:
ant_pairs_nums.remove(ant_tuple)
if ant_str.upper() == 'ALL':
ant_pairs_nums = None
elif len(ant_pairs_nums) == 0:
if (not ant_str.upper() in ['AUTO', 'CROSS']):
ant_pairs_nums = None
if len(polarizations) == 0:
polarizations = None
else:
polarizations.sort(reverse=True)
if print_toggle:
print '\nParsed antenna pairs:'
if ant_pairs_nums is not None:
for pair in ant_pairs_nums:
print pair
print '\nParsed polarizations:'
if polarizations is not None:
for pol in polarizations:
print uvutils.polnum2str(pol)
if len(warned_ants) > 0:
warnings.warn('Warning: Antenna number {a} passed, but not present '
'in the ant_1_array or ant_2_array'
.format(a=(',').join(map(str, warned_ants))))
if len(warned_pols) > 0:
warnings.warn('Warning: Polarization {p} is not present in '
'the polarization_array'
.format(p=(',').join(warned_pols).upper()))
return ant_pairs_nums, polarizations
def combine_uvdata(uvds, inplace=True, run_check=True, check_extra=True,
run_check_acceptability=True):
"""
Combine multiple, non-overlapping UVData objects into a single UVData.
Args:
uvds: list of UVData objects with non-overlapping data
inplace: boolean, if True, create the output UVData object from
the first element in uvds in-place, otherwise make a deepcopy.
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).
"""
# Check that all objects are UVData and valid
for uvd in uvds:
assert isinstance(uvd, UVData), 'Only UVData objects can be added to a UVData object'
uvd.check(check_extra=check_extra, run_check_acceptability=run_check_acceptability)
# assert at least two UVData objects present
assert len(uvds) > 1, "Input uvds must have at least two UVData objects"
# 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']
if inplace:
this = uvds[0]
else:
this = copy.deepcopy(uvds[0])
for uvd in uvds[1:]:
for a in compatibility_params:
if getattr(this, a) != getattr(uvd, a):
msg = 'UVParameter ' + \
a[1:] + ' does not match. Cannot combine objects.'
raise(ValueError(msg))
# Create blt arrays for each uvd for convenience
prec_t = - 2 * \
np.floor(np.log10(this._time_array.tols[-1])).astype(int)
prec_b = 8
uvd_blts = []
for uvd in uvds:
uvd_blts.append(np.array(["_".join(["{1:.{0}f}".format(prec_t, blt[0]),
str(blt[1]).zfill(prec_b)]) for blt in
zip(uvd.time_array, uvd.baseline_array)]))
# get freq and pol arrays for each uvd
uvd_freqs = [uvd.freq_array.squeeze(0) for uvd in uvds]
uvd_pols = [uvd.polarization_array for uvd in uvds]
# iterate over each uvd pair and check for overlapping data
join_axes = {"baseline-time": False, "frequency": False, "polarization": False}
for uvd_pair in list(itertools.combinations(range(len(uvds)), 2)):
uvd1 = uvds[uvd_pair[0]]
uvd2 = uvds[uvd_pair[1]]
# Check we don't have overlapping data
both_pol = np.intersect1d(
uvd1.polarization_array, uvd2.polarization_array)
both_freq = np.intersect1d(
uvd1.freq_array[0, :], uvd2.freq_array[0, :])
both_blts = np.intersect1d(uvd_blts[uvd_pair[0]], uvd_blts[uvd_pair[1]])
if len(both_pol) > 0:
if len(both_freq) > 0:
if len(both_blts) > 0:
raise(ValueError('Objects {} and {} have overlapping data and'
' cannot be combined.'.format(uvd1, uvd2)))
# Set join_axes to True
if len(both_blts) == 0:
join_axes['baseline-time'] = True
if len(both_freq) == 0:
join_axes['frequency'] = True
if len(both_pol) == 0:
join_axes['polarization'] = True
# Get new data axes
new_blts = np.array(sorted(set(np.concatenate(uvd_blts))))
Nblts = len(new_blts)
new_freqs = np.array(sorted(set(np.concatenate(uvd_freqs))))
Nfreqs = len(new_freqs)
new_pols = np.array(sorted(set(np.concatenate(uvd_pols))))
new_pols = new_pols[np.argsort(np.abs(new_pols))]
Npols = len(new_pols)
# Pad out new arrays
data_array = np.zeros((Nblts, 1, Nfreqs, Npols), dtype=np.complex128)
flag_array = np.ones((Nblts, 1, Nfreqs, Npols), dtype=np.bool)
nsample_array = np.zeros((Nblts, 1, Nfreqs, Npols), dtype=np.float64)
uvw_array = np.empty((Nblts, 3), dtype=np.float64)
lst_array = np.empty((Nblts), dtype=np.float64)
ant_1_array = np.empty((Nblts), dtype=np.int64)
ant_2_array = np.empty((Nblts), dtype=np.int64)
zenith_ra = np.empty((Nblts), dtype=np.float64)
zenith_dec = np.empty((Nblts), dtype=np.float64)
time_array, baseline_array = np.split(np.array([blt.split('_') for blt in new_blts]).astype(np.float64), 2, axis=1)
time_array = time_array[:, 0]
baseline_array = baseline_array[:, 0].astype(np.int64)
freq_array = new_freqs[np.newaxis, :]
polarization_array = new_pols.astype(np.int64)
# iterate over uvds and insert data
for i, uvd in enumerate(uvds):
# get indexing arrays
blts_inds = np.nonzero(np.in1d(new_blts, uvd_blts[i]))[0][np.argsort(uvd_blts[i])][:, None, None]
freq_inds = np.nonzero(np.in1d(new_freqs, uvd_freqs[i]))[0][np.argsort(uvd_freqs[i])][None, :, None]
pol_inds = np.nonzero(np.in1d(new_pols, uvd_pols[i]))[0][np.argsort(np.abs(uvd_pols[i]))][None, None, :]
# insert data
data_array[blts_inds, :, freq_inds, pol_inds] = np.moveaxis(uvds[i].data_array, 1, 3)
flag_array[blts_inds, :, freq_inds, pol_inds] = np.moveaxis(uvds[i].flag_array, 1, 3)
nsample_array[blts_inds, :, freq_inds, pol_inds] = np.moveaxis(uvds[i].nsample_array, 1, 3)
uvw_array[blts_inds.squeeze()] = uvds[i].uvw_array
lst_array[blts_inds.squeeze()] = uvds[i].lst_array
ant_1_array[blts_inds.squeeze()] = uvds[i].ant_1_array
ant_2_array[blts_inds.squeeze()] = uvds[i].ant_2_array
zenith_ra[blts_inds.squeeze()] = uvds[i].zenith_ra
zenith_dec[blts_inds.squeeze()] = uvds[i].zenith_dec
# assign to object
this.data_array = data_array
this.flag_array = flag_array
this.nsample_array = nsample_array
this.uvw_array = uvw_array
this.lst_array = lst_array
this.ant_1_array = ant_1_array
this.ant_2_array = ant_2_array
this.time_array = time_array
this.lst_array = lst_array
this.baseline_array = baseline_array
this.freq_array = freq_array
this.polarization_array = polarization_array
this.Ntimes = len(np.unique(this.time_array))
this.Nbls = len(np.unique(this.baseline_array))
this.Nblts = Nblts
this.Nfreqs = Nfreqs
this.Npols = Npols
this.Nants_data = len(
np.unique(this.ant_1_array.tolist() + this.ant_2_array.tolist()))
# only assign zenith arrays if not None in input uvds
if this.zenith_ra is not None:
this.zenith_ra = zenith_ra
if this.zenith_dec is not None:
this.zenith_dec = zenith_dec
# 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.')
# Build up history string
this.history += ' Combined data along {} axis using pyuvdata.'.format(', '.join([k for k in join_axes if join_axes[k]]))
for uvd in uvds[1:]:
this.history = uvutils.combine_histories(this.history, uvd.history)
# Check final object is self-consistent
if run_check:
this.check(check_extra=check_extra,
run_check_acceptability=run_check_acceptability)
return this
