Merge 6a32d45 into a111e56

hera_pspec/__init__.py

@@ -2,7 +2,7 @@
 __init__.py file for hera_pspec
 """
 
-from hera_pspec import version, conversions, grouping, pspecbeam, plot
+from hera_pspec import version, conversions, grouping, pspecbeam, plot, pstokes
 from hera_pspec import uvpspec_utils as uvputils
 
 from hera_pspec.uvpspec import UVPSpec

hera_pspec/data/zen.all.yy.LST.1.06964.uvA/history

@@ -0,0 +1,2 @@
+sig_clip=True, sigma=5.0, min_N=4, rephase=True Input files: zen.2458042.53563.yy.HH.uvXR-zen.2458042.54309.yy.HH.uvXR-zen.2458042.55054.yy.HH.uvXR-zen.2458043.53563.yy.HH.uvXR-zen.2458043.54309.yy.HH.uvXR-zen.2458043.55054.yy.HH.uvXR-zen.2458044.53562.yy.HH.uvXR-zen.2458044.54308.yy.HH.uvXR-zen.2458044.55054.yy.HH.uvXR-zen.2458045.52817.yy.HH.uvXR-zen.2458045.53562.yy.HH.uvXR-zen.2458045.54308.yy.HH.uvXR-zen.2458046.52817.yy.HH.uvXR-zen.2458046.53563.yy.HH.uvXR-zen.2458046.54308.yy.HH.uvXR-zen.2458047.52817.yy.HH.uvXR-zen.2458047.53562.yy.HH.uvXR-zen.2458048.52071.yy.HH.uvXR-zen.2458048.52816.yy.HH.uvXR-zen.2458048.53562.yy.HH.uvXR-zen.2458049.52071.yy.HH.uvXR-zen.2458049.52817.yy.HH.uvXR-zen.2458049.53562.yy.HH.uvXR-zen.2458050.51326.yy.HH.uvXR-zen.2458050.52072.yy.HH.uvXR-zen.2458050.52817.yy.HH.uvXR-zen.2458051.51326.yy.HH.uvXR-zen.2458051.52072.yy.HH.uvXR-zen.2458051.52817.yy.HH.uvXR-zen.2458052.51325.yy.HH.uvXR-zen.2458052.52071.yy.HH.uvXR-zen.2458052.52816.yy.HH.uvXR-zen.2458054.50580.yy.HH.uvXR-zen.2458054.51325.yy.HH.uvXR-zen.2458054.52071.yy.HH.uvXR-zen.2458055.50580.yy.HH.uvXR-zen.2458055.51326.yy.HH.uvXR-zen.2458055.52071.yy.HH.uvXR-zen.2458056.49835.yy.HH.uvXR-zen.2458056.50580.yy.HH.uvXR-zen.2458056.51326.yy.HH.uvXR
+  Read/written with pyuvdata version: 1.2.1.  Git origin: https://github.com/HERA-Team/pyuvdata.  Git hash: 1764956382453fc9582477f99cc73132e3b9c0d0.  Git branch: renumber_ants_miriad.  Git description: v1.2-139-g1764956.  Downselected to specific antennas, times using pyuvdata.

hera_pspec/data/zen.all.yy.LST.1.06964.uvA/vartable

@@ -0,0 +1,32 @@
+r corr
+i nchan
+i npol
+i nspect
+d inttime
+d sdf
+a source
+a telescop
+d latitud
+d longitu
+d antdiam
+i nschan
+i ischan
+i nants
+d antpos
+d sfreq
+i ntimes
+i nbls
+i nblts
+a visunits
+a instrume
+d altitude
+d antnums
+a antnames
+d lst
+d ra
+d dec
+i pol
+d cnt
+d coord
+d time
+r baseline

hera_pspec/pstokes.py

@@ -0,0 +1,271 @@
+"""
+Module to construct pseudo-Stokes (I,Q,U,V) visibilities from miriad files or UVData objects
+"""
+import numpy as np, os
+import pyuvdata
+import copy
+from hera_pspec import version
+from collections import OrderedDict as odict
+
+
+# weights used in forming Stokes visibilities.
+# see pyuvdata.utils.polstr2num for conversion between polarization string
+# and polarization integer. Ex. {'XX': -5, ...}
+pol_weights = {
+    1: odict([(-5, 1.), (-6, 1.)]),
+    2: odict([(-5, 1.), (-6, -1.)]),
+    3: odict([(-7, 1.), (-8, 1.)]),
+    4: odict([(-7, -1.j), (-8, 1.j)])
+}
+
+def miriad2pyuvdata(dset, antenna_nums=None, ant_pairs_nums=None, polarizations=None,
+                    ant_str=None, time_range=None):
+    """
+    Reads-in a Miriad filepath to a UVData object
+
+    Parameters
+    ----------
+    dset : str
+        Miriad file to convert to UVData object containing visibilities and corresponding metadata
+
+    antenna_nums: integer list
+        The antennas numbers to read into the object.
+
+    ant_pairs_nums: list of tuples
+        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: 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: integer or string list
+        List of polarization integers or strings to read-in.
+        Ex: ['xx', 'yy', ...]
+
+    time_range: float list
+        len-2 list containing min and max range of times (Julian Date) to read-in.
+        Ex: [2458115.20, 2458115.40]
+
+    Returns
+    -------
+    uvd : pyuvdata.UVData object
+    """
+    uvd = pyuvdata.UVData()
+    uvd.read_miriad(dset, antenna_nums=antenna_nums, ant_pairs_nums=ant_pairs_nums, 
+                    polarizations=polarizations, ant_str=ant_str, time_range=time_range)
+    return uvd
+
+
+def _combine_pol(uvd1, uvd2, pol1, pol2, pstokes='pI'):
+    """
+    Combines UVData visibilities to form the desired pseudo-stokes visibilities. 
+    It returns UVData object containing the pseudo-stokes visibilities
+
+    Parameters
+    ----------
+    uvd1 : UVData object
+        First UVData object containing data that is used to
+        form Stokes visibilities
+
+    uvd2 : UVData oject
+        Second UVData objects containing data that is used to
+        form Stokes visibilities
+
+    pol1 : Polarization, type: str
+        Polarization of the first UVData object to use in constructing 
+        pStokes visibility.
+
+    pol2 : Polarization, type: str
+        Polarization of the second UVData object to use in constructing 
+        pStokes visibility.
+
+    pstokes: Pseudo-stokes polarization to form, type: str
+        Pseudo stokes polarization to form, can be 'pI' or 'pQ' or 'pU' or 'pV'. 
+        Default: pI
+
+    Returns
+    -------
+    uvdS : UVData object
+    """
+    assert isinstance(uvd1, pyuvdata.UVData), "uvd1 must be a pyuvdata.UVData instance"
+    assert isinstance(uvd2, pyuvdata.UVData), "uvd2 must be a pyuvdata.UVData instance"
+
+    # convert pol1 and/or pol2 to integer if fed as a string
+    if isinstance(pol1, (str, np.str)):
+        pol1 = pyuvdata.utils.polstr2num(pol1)
+    if isinstance(pol2, (str, np.str)):
+        pol2 = pyuvdata.utils.polstr2num(pol2)
+
+    # extracting data array from the UVData objects
+    data1 = uvd1.data_array
+    data2 = uvd2.data_array
+
+    # extracting flag array from the UVdata objects
+    flag1 = uvd1.flag_array
+    flag2 = uvd2.flag_array
+
+    # constructing flags (boolean)
+    flag = np.logical_or(flag1, flag2)
+
+    # convert pStokes to polarization integer if a string
+    if isinstance(pstokes, (str, np.str)):
+        pstokes = pyuvdata.utils.polstr2num(pstokes)
+
+    # get string form of polarizations
+    pol1_str = pyuvdata.utils.polnum2str(pol1)
+    pol2_str = pyuvdata.utils.polnum2str(pol2)
+    pstokes_str = pyuvdata.utils.polnum2str(pstokes)
+
+    # assert pstokes in pol_weights, and pol1 and pol2 in pol_weights[pstokes]
+    assert pstokes in pol_weights, "unrecognized pstokes parameter {}".format(pstokes_str)
+    assert pol1 in pol_weights[pstokes], "pol1 {} not used in constructing pstokes {}".format(pol1_str, pstokes_str)
+    assert pol2 in pol_weights[pstokes], "pol2 {} not used in constructing pstokes {}".format(pol2_str, pstokes_str)
+
+    # constructing Stokes visibilities
+    stdata = 0.5 * (pol_weights[pstokes][pol1]*data1 + pol_weights[pstokes][pol2]*data2)
+
+    # assigning and writing data, flags and metadata to UVData object
+    uvdS = copy.deepcopy(uvd1)
+    uvdS.data_array = stdata  # pseudo-stokes data
+    uvdS.flag_array = flag  # flag array
+    uvdS.polarization_array = np.array([pstokes], dtype=np.int)  # polarization number
+    uvdS.nsample_array = uvd1.nsample_array + uvd2.nsample_array  # nsamples
+    uvdS.history = "Merged into pseudo-stokes vis with hera_pspec version {} Git hash {}\n{}" \
+                    "{}{}{}{}\n".format(version.version, version.git_hash, "-"*20+'\n',
+                    'dset1 history:\n', uvd1.history, '\n'+'-'*20+'\ndset2 history:\n',
+                    uvd2.history)
+
+    return uvdS
+
+
+def construct_pstokes(dset1, dset2, pstokes='pI', run_check=True, antenna_nums=None, 
+                      ant_pairs_nums=None, polarizations=None, ant_str=None, time_range=None,
+                      history=''):
+    """
+    Validates datasets required to construct desired visibilities and 
+    constructs desired pseudo-Stokes visibilities. These are formed
+    via the following expression
+
+        ( V_pI )            ( 1  0  0  1 )   ( V_XX )
+        | V_pQ |            | 1  0  0 -1 |   | V_XY |
+        | V_pU |    = 0.5 * | 0  1  1  0 | * | V_YX |
+        ( V_pV )            ( 0 -i  i  0 )   ( V_YY )
+
+    In constructing a given pseudo-Stokes visibilities, the XX or XY polarization is 
+    taken from dset1, and the YX or YY pol is taken from dset2.
+
+    Parameters
+    ----------
+    dset1 : UVData object or Miriad file
+        First UVData object or Miriad file containing data that is used to
+        form Stokes visibilities
+
+    dset2 : UVData oject or Miriad file
+        Second UVData object or Miriad file containing data that is used to
+        form Stokes visibilities
+
+    pstokes: Stokes polarization, type: str
+        Pseudo stokes polarization to form, can be 'pI' or 'pQ' or 'pU' or 'pV'. 
+        Default: pI
+
+    run_check: boolean
+        Option to check for the existence and proper shapes of
+        parameters after downselecting data on this object. Default is True.
+
+    antenna_nums: integer list
+        The antennas numbers to read into the object.
+
+    ant_pairs_nums: list of tuples
+        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: 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: integer or string list
+        List of polarization integers or strings to read-in.
+        Ex: ['xx', 'yy', ...]
+
+    time_range: float list
+        len-2 list containing min and max range of times (Julian Date) to read-in.
+        Ex: [2458115.20, 2458115.40]
+
+    history : str
+        Extra history string to add to concatenated pseudo-Stokes visibility.
+
+    Returns
+    -------
+    uvdS : UVData object with pseudo-Stokes visibility
+    """
+    # convert dset1 and dset2 to UVData objects if they are miriad files
+    if isinstance(dset1, pyuvdata.UVData) == False:
+        assert isinstance(dset1, (str, np.str)), "dset1 must be fed as a string or UVData object"
+        uvd1 = miriad2pyuvdata(dset1, antenna_nums=antenna_nums, ant_pairs_nums=ant_pairs_nums,
+                               polarizations=polarizations, ant_str=ant_str, time_range=time_range)
+    else:
+        uvd1 = dset1
+    if isinstance(dset2, pyuvdata.UVData) == False:
+        assert isinstance(dset2, (str, np.str)), "dset2 must be fed as a string or UVData object"
+        uvd2 = miriad2pyuvdata(dset2, antenna_nums=antenna_nums, ant_pairs_nums=ant_pairs_nums,
+                               polarizations=polarizations, ant_str=ant_str, time_range=time_range)
+    else:
+        uvd2 = dset2
+
+    # convert pstokes to integer if fed as a string
+    if isinstance(pstokes, (str, np.str)):
+        pstokes = pyuvdata.utils.polstr2num(pstokes)
+
+    # check if dset1 and dset2 habe the same spectral window
+    spw1 = uvd1.spw_array
+    spw2 = uvd2.spw_array
+    assert (spw1 == spw2), "dset1 and dset2 must have the same spectral windows."
+
+    # check if dset1 and dset2 have the same frequencies
+    freqs1 = uvd1.freq_array
+    freqs2 = uvd2.freq_array
+    if np.array_equal(freqs1, freqs2) == False:
+        raise ValueError("dset1 and dset2 must have the same frequencies.")
+
+    # check if dset1 and dset2 have the same timestamps
+    times1 = uvd1.time_array
+    times2 = uvd2.time_array
+    if np.array_equal(times1, times2) == False:
+        raise ValueError("dset1 and dset2 must have the same timestamps.")
+
+    # check if dset1 and dset2 have the same baselines
+    bls1 = uvd1.baseline_array
+    bls2 = uvd2.baseline_array
+    if np.array_equal(bls1, bls2) == False:
+        raise ValueError("dset1 and dset2 must have the same baselines")
+
+    # makes the Npol length==1 so that the UVData carries data for the required polarization only
+    st_keys = pol_weights[pstokes].keys()
+    req_pol1 = st_keys[0]
+    req_pol2 = st_keys[1]
+
+    # check polarizations of UVData objects are consistent with the required polarization 
+    # to form the desired pseudo Stokes visibilities. If multiple exist, downselect on polarization.
+    assert req_pol1 in uvd1.polarization_array, "Polarization {} not found in dset1 object".format(req_pol1)
+    if uvd1.Npols > 1:
+        uvd1 = uvd1.select(polarizations=req_pol1, inplace=False)
+
+    assert req_pol2 in uvd2.polarization_array, "Polarization {} not found in dset2 object".format(req_pol2)
+    if uvd2.Npols > 1:
+        uvd2 = uvd2.select(polarizations=req_pol1, inplace=False)
+
+    # combining visibilities to form the desired Stokes visibilties
+    uvdS = _combine_pol(uvd1=uvd1, uvd2=uvd2, pol1=req_pol1, pol2=req_pol2, pstokes=pstokes)
+    uvdS.history += history
+
+    if run_check:
+        uvdS.check()
+
+    return uvdS

hera_pspec/tests/test_plot.py

@@ -46,7 +46,7 @@ def setUp(self):
                         bls, exclude_permutations=False, exclude_auto_bls=True)
 
         # Calculate the power spectrum
-        self.uvp = self.ds.pspec(self.bls1, self.bls2, (0, 1), 
+        self.uvp = self.ds.pspec(self.bls1, self.bls2, (0, 1), ('xx','xx'),  
                                  spw_ranges=[(300, 400), (600,721)], 
                                  input_data_weight='identity', norm='I', 
                                  taper='blackman-harris', verbose=False)