Raw File
Tip revision: 4b5dfe6ed6e3067d14f2e762d768a82229681b18 authored by EXTERNAL-Ewall-Wice on 02 November 2018, 05:08 UTC
now output separate file per polarization.
Tip revision: 4b5dfe6
# -*- coding: utf-8 -*-
# Copyright 2018 the HERA Project
# Licensed under the 2-clause BSD License

"""This module extracts some Python code from AIPY used in our MIRIAD I/O
routines. It was copied from AIPY commit
6cb5a70876f33dccdd68d4063b076f8d42d9edae, then reformatted. The only items
used by pyuvdata are ``uv_selector`` and ``UV``.

from __future__ import absolute_import, division, print_function

__all__ = [

import numpy as np
import re

    from . import _miriad
except ImportError:
    # See for an explanation. The workaround here is kind of annoying
    # since the code below needs the name `_miriad.UV` to be a class in order
    # to import successfully.
    import os

    if os.environ.get('PYUVDATA_IGNORE_EXTMOD_IMPORT_FAIL', '') != '1':

    class UV(object):

    _miriad = UV()  # eww gross but it works
    _miriad.UV = UV

str2pol = {
    'I': 1,  # Stokes Paremeters
    'Q': 2,
    'U': 3,
    'V': 4,
    'rr': -1,  # Circular Polarizations
    'll': -2,
    'rl': -3,
    'lr': -4,
    'xx': -5,  # Linear Polarizations
    'yy': -6,
    'xy': -7,
    'yx': -8,

def bl2ij(bl):
    bl = int(bl)

    if bl > 65536:
        bl -= 65536
        mant = 2048
        mant = 256

    return (bl // mant - 1, bl % mant - 1)

def ij2bl(i, j):
    if i > j:
        i, j = j, i

    if j + 1 < 256:
        return 256 * (i + 1) + (j + 1)

    return 2048 * (i + 1) + (j + 1) + 65536

ant_re = r'(\(((-?\d+[xy]?,?)+)\)|-?\d+[xy]?)'
bl_re = '(^(%s_%s|%s),?)' % (ant_re, ant_re, ant_re)

def parse_ants(ant_str, nants):
    """Generate list of (baseline, include, pol) tuples based on parsing of the
    string associated with the 'ants' command-line option.

    rv, cnt = [], 0

    while cnt < len(ant_str):
        m =, ant_str[cnt:])

        if m is None:
            if ant_str[cnt:].startswith('all'):
                rv = []
            elif ant_str[cnt:].startswith('auto'):
                rv.append(('auto', 1, -1))
            elif ant_str[cnt:].startswith('cross'):
                rv.append(('auto', 0, -1))
                raise ValueError('Unparsable ant argument "%s"' % ant_str)

            c = ant_str[cnt:].find(',')

            if c >= 0:
                cnt += c + 1
                cnt = len(ant_str)
            m = m.groups()
            cnt += len(m[0])

            if m[2] is None:
                ais = [m[8]]
                ajs = list(range(nants))
                if m[3] is None:
                    ais = [m[2]]
                    ais = m[3].split(',')

                if m[6] is None:
                    ajs = [m[5]]
                    ajs = m[6].split(',')

            for i in ais:
                if type(i) == str and i.startswith('-'):
                    i = i[1:]  # nibble the - off the string
                    include_i = 0
                    include_i = 1

                for j in ajs:
                    include = None

                    if type(j) == str and j.startswith('-'):
                        j = j[1:]
                        include_j = 0
                        include_j = 1

                    include = int(include_i and include_j)
                    pol = None
                    i, j = str(i), str(j)

                    if not i.isdigit():
                        ai ='(\d+)([x,y])', i).groups()
                    if not j.isdigit():
                        aj ='(\d+)([x,y])', j).groups()

                    if i.isdigit() and not j.isdigit():
                        pol = ['x' + aj[1], 'y' + aj[1]]
                        ai = [i, '']
                    elif not i.isdigit() and j.isdigit():
                        pol = [ai[1] + 'x', ai[1] + 'y']
                        aj = [j, '']
                    elif not i.isdigit() and not j.isdigit():
                        pol = [ai[1] + aj[1]]

                    if pol is not None:
                        bl = ij2bl(abs(int(ai[0])), abs(int(aj[0])))
                        for p in pol:
                            rv.append((bl, include, p))
                        bl = ij2bl(abs(int(i)), abs(int(j)))
                        rv.append((bl, include, -1))
    return rv

def uv_selector(uv, ants=-1, pol_str=-1):
    """Call with appropriate options based on string argument for
    antennas (can be 'all', 'auto', 'cross', '0,1,2', or '0_1,0_2') and string
    for polarization ('xx','yy','xy','yx').

    if ants != -1:
        if type(ants) == str:
            ants = parse_ants(ants, uv['nants'])

        for cnt, (bl, include, pol) in enumerate(ants):
            if cnt > 0:
                if include:
          'or', -1, -1)
          'and', -1, -1)

            if pol == -1:
                pol = pol_str  # default to explicit pol parameter

            if bl == 'auto':
      'auto', 0, 0, include=include)
                i, j = bl2ij(bl)
      'antennae', i, j, include=include)

            if pol != -1:
                for p in pol.split(','):
                    polopt = str2pol[p]
          'polarization', polopt, 0)
    elif pol_str != -1:
        for p in pol_str.split(','):
            polopt = str2pol[p]
  'polarization', polopt, 0)

itemtable = {
    'obstype': 'a',
    'history': 'a',
    'vartable': 'a',
    'ngains': 'i',
    'nfeeds': 'i',
    'ntau': 'i',
    'nsols': 'i',
    'interval': 'd',
    'leakage': 'c',
    'freq0': 'd',
    'freqs': '?',
    'bandpass': 'c',
    'nspect0': 'i',
    'nchan0': 'i',
    'stopt': 'd',
    'duration': 'd',

def _uv_pipe_default_action(uv, p, d):
    return p, d

class UV(_miriad.UV):
    """Top-level interface to a Miriad UV data set.


    def __init__(self, filename, status='old', corrmode='r'):
        """Open a miriad file. status can be ('old','new','append'). corrmode can be
        'r' (float32 data storage) or 'j' (int16 with shared exponent).
        Default is 'r'.

        assert status in ['old', 'new', 'append']
        assert corrmode in ['r', 'j']

        _miriad.UV.__init__(self, filename, status, corrmode)

        self.status = status
        self.nchan = _miriad.MAXCHAN

        if status == 'old':
            self.vartable = self._gen_vartable()
            self.rewind()  # Update variables for the user
                self.nchan = self['nchan']
            except KeyError:
            self.vartable = {'corr': corrmode}

    def _gen_vartable(self):
        """Generate table of variables and types from the vartable header.

        vartable = {}
        for line in self._rdhd('vartable').split('\n'):
                type, name = line.split()
                vartable[name] = type
            except ValueError:

        return vartable

    def vars(self):
        """Return a list of available variables.

        return list(self.vartable.keys())

    def items(self):
        """Return a list of available header items.

        items = []

        for i in itemtable:
                _miriad.hdaccess(self.haccess(i, 'read'))
            except IOError:
        return items

    def _rdhd(self, name):
        """Provide read access to header items via low-level calls.

        itype = itemtable[name]

        if itype == '?':
            return self._rdhd_special(name)

        h = self.haccess(name, 'read')
        rv = []

        if len(itype) == 1:
            if itype == 'a':
                offset = 0
                t, offset = _miriad.hread_init(h)
                assert itype == t

            while True:
                    c, o = _miriad.hread(h, offset, itype)
                except IOError:

                if itype == 'a':
                        c = str(c[:o], 'utf-8')
                    except TypeError:
                        c = c[:o]

                offset += o

            if itype == 'a':
                rv = ''.join(rv)
            t, offset = _miriad.hread_init(h)
            assert t == 'b'

            for t in itype:
                v, o = _miread.hread(h, offset, t)
                offset += o


        if len(rv) == 1:
            return rv[0]
        elif type(rv) == str:
            return rv
            return np.array(rv)

    def _wrhd(self, name, val):
        """Provide write access to header items via low-level calls.

        type = itemtable[name]

        if type == '?':
            return self._wrhd_special(name, val)

        h = self.haccess(name, 'write')

        if len(type) == 1:
            except TypeError:
                val = [val]

            if type == 'a':
                offset = 0
                offset = _miriad.hwrite_init(h, type)

            for v in val:
                offset += _miriad.hwrite(h, offset, v, type)
            offset = _miriad.hwrite_init(h, 'b')
            for v, t in zip(val, type):
                offset += _miriad.hwrite(h, offset, v, t)


    def _rdhd_special(self, name):
        """Provide read access to special header items of type '?' to _rdhd.

        if name == 'freqs':
            h = self.haccess(name, 'read')
            c, o = _miriad.hread(h, 0, 'i')
            rv = [c]
            offset = 8

            while True:
                    c, o = _miriad.hread(h, offset, 'i')
                    offset += 8

                    c, o = _miriad.hread(h, offset, 'd')
                    offset += 8

                    c, o = _miriad.hread(h, offset, 'd')
                    offset += 8
                except IOError:

            return rv
            raise ValueError('Unknown special header: ' + name)

    def _wrhd_special(self, name, val):
        """Provide write access to special header items of type '?' to _wrhd

        if name == 'freqs':
            h = self.haccess(name, 'write')
            o = _miriad.hwrite(h, 0, val[0], 'i')
            offset = 8

            for i, v in enumerate(val[1:]):
                if i % 3 == 0:
                    o = _miriad.hwrite(h, offset, v, 'i')
                    o = _miriad.hwrite(h, offset, v, 'd')
                offset += 8

            raise ValueError('Unknown special header: ' + name)

    def __getitem__(self, name):
        """Allow access to variables and header items via ``uv[name]``.

            type = self.vartable[name]
            return self._rdvr(name, type)
        except KeyError:
            type = itemtable[name]
            return self._rdhd(name)

    def __setitem__(self, name, val):
        """Allow setting variables and header items via ``uv[name] = val``.

            type = self.vartable[name]
            self._wrvr(name, type, val)
        except KeyError:
            self._wrhd(name, val)

    def select(self, name, n1, n2, include=True):
        """Choose which data are returned by read().

            name    This can be: 'decimate','time','antennae','visibility',
                    'ddec','uvnrange','increment','ra','dec','and', 'clear',
            n1,n2   Generally this is the range of values to select. For
                    'antennae', this is the two antennae pair to select
                    (indexed from 0); a -1 indicates 'all antennae'.
                    For 'decimate', n1 is every Nth integration to use, and
                    n2 is which integration within a block of N to use.
                    For 'shadow', a zero indicates use 'antdiam' variable.
                    For 'on','window','polarization','increment','shadow' only
                    p1 is used.
                    For 'and','or','clear','auto' p1 and p2 are ignored.
            include If true, the data is selected. If false, the data is
                    discarded. Ignored for 'and','or','clear'.

        if name == 'antennae':
            n1 += 1
            n2 += 1
        self._select(name, float(n1), float(n2), int(include))

    def read(self, raw=False):
        """Return the next data record. Calling this function causes vars to change to
        reflect the record which this function returns. 'raw' causes data and
        flags to be returned seperately.

        preamble, data, flags, nread = self.raw_read(self.nchan)

        if nread == 0:
            raise IOError("No data read")

        flags = np.logical_not(flags)

        if raw:
            return preamble, data, flags
        return preamble,, mask=flags)

    def all(self, raw=False):
        """Provide an iterator over preamble, data. Allows constructs like: ``for
        preamble, data in uv.all(): ...``

        while True:
            except IOError:

    def write(self, preamble, data, flags=None):
        """Write the next data record. data must be a complex, masked array. preamble
        must be (uvw, t, (i,j)), where uvw is an array of u,v,w, t is the
        Julian date, and (i,j) is an antenna pair.

        if data is None:

        if flags is not None:
            flags = np.logical_not(flags)
        elif len(data.mask.shape) == 0:
            flags = np.ones(data.shape)
            data = data.unmask()
            flags = np.logical_not(data.mask)
            data =

        self.raw_write(preamble, data.astype(np.complex64), flags.astype(np.int32))

    def init_from_uv(self, uv, override={}, exclude=[]):
        """Initialize header items and variables from another UV. Those in override
        will be overwritten by override[k], and tracking will be turned off
        (meaning they will not be updated in pipe()). Those in exclude are
        omitted completely.

        for k in uv.items():
            if k in exclude:
            elif k in override:
                self._wrhd(k, override[k])
                self._wrhd(k, uv[k])

        self.vartable = {}

        for k in uv.vars():
            if k in exclude:
            elif k == 'corr':
                # I don't understand why reading 'corr' segfaults miriad,
                # but it does.  This is a cludgy work-around.
            elif k in override:
                self.vartable[k] = uv.vartable[k]
                self._wrvr(k, uv.vartable[k], override[k])
                self.vartable[k] = uv.vartable[k]
                self._wrvr(k, uv.vartable[k], uv[k])
                uv.trackvr(k, 'c')  # Set to copy when copyvr() called

    def pipe(self, uv, mfunc=_uv_pipe_default_action, append2hist='', raw=False):
        """Pipe in data from another UV through the function ``mfunc(uv, preamble,
        data)``, which should return ``(preamble, data)``. If mfunc is not
        provided, the dataset will just be cloned, and if the returned data is
        None, it will be omitted. The string 'append2hist' will be appended to

        self._wrhd('history', self['history'] + append2hist)

        if raw:
            for p, d, f in uv.all(raw=raw):
                np, nd, nf = mfunc(uv, p, d, f)
                self.write(np, nd, nf)
            for p, d in uv.all():
                np, nd = mfunc(uv, p, d)
                self.write(np, nd)

    def add_var(self, name, type):
        """Add a variable of the specified type to a UV file.

        self.vartable[name] = type
back to top