Merge fa18102 into 145874a

README.md

@@ -46,7 +46,6 @@ pyuvdata has three major user classes:
 ## Known Issues and Planned Improvements
 * UVData: different multiple spectral windows or multiple sources are not currently supported
 * UVData: testing against miriad package
-* UVData: replacing pyephem with astropy+NOVAS for time and phase calculations
 * UVData: add support for writing CASA measurement sets
 * UVData: phasing is tested to a part in 10^3, and assumes planar array. Improvements are tracked on Issue \#148.
 * UVCal: expand support for calibration solutions: support other formats beyond FITS
@@ -98,12 +97,11 @@ python setup.py build_ext --inplace
 The numpy and astropy versions are important, so be sure to make sure these are up to date before you install.
 
 For anaconda users, we suggest using conda to install astropy, numpy and scipy and conda-forge
-for pyephem and optionally casacore-python and healpy (e.g. ```conda install -c conda-forge pyephem```).
+for optionally installing casacore-python and healpy (e.g. ```conda install -c conda-forge python-casacore```).
 
 * numpy >= 1.10
 * scipy
 * astropy >= 1.2
-* pyephem
 * casacore-python (optional: for CASA measurement set reading functionality)
 * healpy (optional: working with beams in HEALPix formats)
 

docs/tutorial.rst

@@ -222,23 +222,23 @@ Phasing/unphasing data
 ::
 
   >>> from pyuvdata import UVData
-  >>> import ephem
+  >>> from astropy.time import Time
   >>> UV = UVData()
   >>> miriad_file = 'pyuvdata/data/new.uvA'
   >>> UV.read_miriad(miriad_file)
   >>> print(UV.phase_type)
   drift
 
   # Phase the data to the zenith at first time step
-  >>> UV.phase_to_time(UV.time_array[0])
+  >>> UV.phase_to_time(Time(UV.time_array[0], format='jd'))
   >>> print(UV.phase_type)
   phased
 
   # Undo phasing to try another phase center
   >>> UV.unphase_to_drift()
 
   # Phase to a specific ra/dec/epoch (in radians)
-  >>> UV.phase(5.23368, 0.710940, ephem.J2000)
+  >>> UV.phase(5.23368, 0.710940, epoch="J2000")
 
 UVData: Plotting
 ------------------

pyuvdata/miriad.py

@@ -6,6 +6,7 @@
 from __future__ import absolute_import, division, print_function
 
 from astropy import constants as const
+from astropy.coordinates import Angle
 import os
 import shutil
 import numpy as np
@@ -83,7 +84,7 @@ def read_miriad(self, filepath, correct_lat_lon=True, run_check=True,
                                     'lst', 'pol', 'nants', 'antnames', 'nblts',
                                     'ntimes', 'nbls', 'sfreq', 'epoch',
                                     'antpos', 'antnums', 'degpdy', 'antdiam',
-                                    ]
+                                    'phsframe']
         # list of miriad variables not read, but also not interesting
         # NB: nspect (I think) is number of spectral windows, will want one day
         # NB: xyphase & xyamp are "On-line X Y phase/amplitude measurements" which we may want in
@@ -222,7 +223,7 @@ def read_miriad(self, filepath, correct_lat_lon=True, run_check=True,
             if ant_str != 'all':
                 history_update_string += 'antenna pairs'
                 n_selects += 1
-        # select on antenna_nums and/or bls using aipy.uv_selector
+        # select on antenna_nums and/or bls using aipy_extracts.uv_selector
         if antenna_nums is not None or bls is not None:
             antpair_str = ''
             if antenna_nums is not None:
@@ -287,7 +288,13 @@ def read_miriad(self, filepath, correct_lat_lon=True, run_check=True,
             assert isinstance(time_range, (list, np.ndarray)), err_msg
             assert len(time_range) == 2, err_msg
             assert np.array([isinstance(t, (float, np.float, np.float64)) for t in time_range]).all(), err_msg
-            uv.select('time', time_range[0], time_range[1], include=True)
+
+            # UVData.time_array marks center of integration, while Miriad 'time' marks beginning
+            # assume time_range refers to the center of the integrations,
+            # so subtract 1/2 an integration before using with miriad select
+            time_range_use = np.array(time_range) - uv['inttime'] / (24 * 3600.) / 2
+
+            uv.select('time', time_range_use[0], time_range_use[1], include=True)
             if n_selects > 0:
                 history_update_string += ', times'
             else:
@@ -697,7 +704,9 @@ def read_miriad(self, filepath, correct_lat_lon=True, run_check=True,
             # The select on read will make the header ntimes not match the number of unique times
             self.Ntimes = len(times)
 
-        self.time_array = t_grid
+        # UVData.time_array marks center of integration, while Miriad 'time' marks beginning
+        self.time_array = t_grid + uv['inttime'] / (24 * 3600.) / 2
+
         self.ant_1_array = ant_i_grid.astype(int)
         self.ant_2_array = ant_j_grid.astype(int)
 
@@ -830,12 +839,23 @@ def read_miriad(self, filepath, correct_lat_lon=True, run_check=True,
             self.phase_center_ra = float(ra_list[0])
             self.phase_center_dec = float(dec_list[0])
             self.phase_center_epoch = uv['epoch']
+            if 'phsframe' in uv.vartable.keys():
+                self.phase_center_frame = uv['phsframe'].replace('\x00', '')
         else:
             # check that the RA values are not constant (if more than one time present)
             if (single_ra and not single_time):
                 raise ValueError('phase_type is "drift" but the RA values are constant.')
-            self.zenith_ra = ra_list
-            self.zenith_dec = dec_list
+
+            ra_diff = np.abs(ra_list - self.lst_array)
+            dec_diff = np.abs(dec_list - latitude)
+            acceptable_offset = Angle('1d')
+            if (np.max(ra_diff) > acceptable_offset.rad
+                    or np.max(dec_diff) > acceptable_offset.rad):
+                warnings.warn('drift RA and/or Dec is off from lst and/or '
+                              'latitude by more than {}, '
+                              'so it appears that it is not a zenith drift scan. '
+                              'Setting phase_type to "unknown"'.format(acceptable_offset))
+                self.set_unknown_phase_type()
 
         try:
             self.set_telescope_params()
@@ -866,6 +886,12 @@ def write_miriad(self, filepath, run_check=True, check_extra=True,
             no_antnums: Option to not write the antnums variable to the file.
                 Should only be used for testing purposes.
         """
+        # change time_array and lst_array to mark beginning of integration, per Miriad format
+        miriad_time_array = self.time_array - self.integration_time / (24 * 3600.) / 2
+        if self.telescope_location is not None:
+            latitude, longitude, altitude = self.telescope_location_lat_lon_alt_degrees
+            miriad_lsts = uvutils.get_lst_for_time(miriad_time_array, latitude, longitude, altitude)
+
         if run_check:
             self.check(check_extra=check_extra,
                        run_check_acceptability=run_check_acceptability)
@@ -984,6 +1010,9 @@ def write_miriad(self, filepath, run_check=True, check_extra=True,
         if self.phase_type == 'phased':
             uv.add_var('epoch', 'r')
             uv['epoch'] = self.phase_center_epoch
+            if self.phase_center_frame is not None:
+                uv.add_var('phsframe', 'a')
+                uv['phsframe'] = self.phase_center_frame
 
         # required pyuvdata variables that are not recognized miriad variables
         uv.add_var('ntimes', 'i')
@@ -1088,18 +1117,18 @@ def write_miriad(self, filepath, run_check=True, check_extra=True,
         # write data
         c_ns = const.c.to('m/ns').value
         for viscnt, blt in enumerate(self.data_array):
-            uvw = (self.uvw_array[viscnt] / c_ns).astype(np.double)  # NOTE issue 50 on conjugation
-            t = self.time_array[viscnt]
+            uvw = (self.uvw_array[viscnt] / c_ns).astype(np.double)
+            t = miriad_time_array[viscnt]
             i = self.ant_1_array[viscnt]
             j = self.ant_2_array[viscnt]
 
-            uv['lst'] = self.lst_array[viscnt]
+            uv['lst'] = miriad_lsts[viscnt]
             if self.phase_type == 'phased':
                 uv['ra'] = self.phase_center_ra
                 uv['dec'] = self.phase_center_dec
             elif self.phase_type == 'drift':
-                uv['ra'] = self.zenith_ra[viscnt]
-                uv['dec'] = self.zenith_dec[viscnt]
+                uv['ra'] = miriad_lsts[viscnt]
+                uv['dec'] = self.telescope_location_lat_lon_alt[0]
             else:
                 raise ValueError('The phasing type of the data is unknown. '
                                  'Set the phase_type to "drift" or "phased" to '

pyuvdata/tests/test_miriad.py

@@ -9,8 +9,8 @@
 import shutil
 import copy
 import numpy as np
-import ephem
 import nose.tools as nt
+from astropy.time import Time, TimeDelta
 from pyuvdata import UVData
 import pyuvdata.utils as uvutils
 import pyuvdata.tests as uvtest
@@ -72,15 +72,6 @@ def test_ReadMiriadWriteUVFits():
                          known_warning='miriad')
     miriad_uv.write_uvfits(testfile, spoof_nonessential=True, force_phase=True)
     uvfits_uv.read_uvfits(testfile)
-    # these are not equal because miriad_uv still retains zenith_dec and
-    # zenith_ra, which are not present in uvfits_uv
-    nt.assert_false(miriad_uv == uvfits_uv)
-    # they are equal if only required parameters are checked:
-    nt.assert_true(miriad_uv.__eq__(uvfits_uv, check_extra=False))
-
-    # remove zenith_ra and zenith_dec to test that the rest of the objects are equal
-    miriad_uv.zenith_ra = None
-    miriad_uv.zenith_dec = None
     nt.assert_equal(miriad_uv, uvfits_uv)
 
     # check error if phase_type is wrong and force_phase not set
@@ -153,9 +144,10 @@ def test_wronglatlon():
     lonfile = os.path.join(DATA_PATH, 'zen.2456865.60537_wronglon.xy.uvcRREAA')
     telescopefile = os.path.join(DATA_PATH, 'zen.2456865.60537_wrongtelecope.xy.uvcRREAA')
 
-    uvtest.checkWarnings(uv_in.read_miriad, [latfile], nwarnings=2,
+    uvtest.checkWarnings(uv_in.read_miriad, [latfile], nwarnings=3,
                          message=['Altitude is not present in file and latitude value does not match',
-                                  latfile + ' was written with an old version of pyuvdata'])
+                                  latfile + ' was written with an old version of pyuvdata',
+                                  'drift RA and/or Dec is off from lst and/or latitude by more than 1.0 deg'])
     uvtest.checkWarnings(uv_in.read_miriad, [lonfile], nwarnings=2,
                          message=['Altitude is not present in file and longitude value does not match',
                                   lonfile + ' was written with an old version of pyuvdata'])
@@ -240,14 +232,16 @@ def test_miriad_location_handling():
     # close file properly
     del(aipy_uv2)
 
-    uvtest.checkWarnings(uv_out.read_miriad, [testfile], nwarnings=3,
+    uvtest.checkWarnings(uv_out.read_miriad, [testfile], nwarnings=4,
                          message=['Altitude is not present in Miriad file, and '
                                   'telescope foo is not in known_telescopes. '
                                   'Telescope location will be set using antenna positions.',
                                   'Telescope location is set at sealevel at the '
                                   'file lat/lon coordinates. Antenna positions '
                                   'are present, but the mean antenna latitude '
                                   'value does not match',
+                                  'drift RA and/or Dec is off from lst and/or '
+                                  'latitude by more than 1.0 deg',
                                   'Telescope foo is not in known_telescopes.'])
 
     # Test for handling when antenna positions have a different mean longitude than the file longitude
@@ -268,14 +262,16 @@ def test_miriad_location_handling():
     # close file properly
     del(aipy_uv2)
 
-    uvtest.checkWarnings(uv_out.read_miriad, [testfile], nwarnings=3,
+    uvtest.checkWarnings(uv_out.read_miriad, [testfile], nwarnings=4,
                          message=['Altitude is not present in Miriad file, and '
                                   'telescope foo is not in known_telescopes. '
                                   'Telescope location will be set using antenna positions.',
                                   'Telescope location is set at sealevel at the '
                                   'file lat/lon coordinates. Antenna positions '
                                   'are present, but the mean antenna longitude '
                                   'value does not match',
+                                  'drift RA and/or Dec is off from lst and/or '
+                                  'latitude by more than 1.0 deg',
                                   'Telescope foo is not in known_telescopes.'])
 
     # Test for handling when antenna positions have a different mean longitude &
@@ -295,14 +291,16 @@ def test_miriad_location_handling():
     # close file properly
     del(aipy_uv2)
 
-    uvtest.checkWarnings(uv_out.read_miriad, [testfile], nwarnings=3,
+    uvtest.checkWarnings(uv_out.read_miriad, [testfile], nwarnings=4,
                          message=['Altitude is not present in Miriad file, and '
                                   'telescope foo is not in known_telescopes. '
                                   'Telescope location will be set using antenna positions.',
                                   'Telescope location is set at sealevel at the '
                                   'file lat/lon coordinates. Antenna positions '
                                   'are present, but the mean antenna latitude and '
                                   'longitude values do not match',
+                                  'drift RA and/or Dec is off from lst and/or '
+                                  'latitude by more than 1.0 deg',
                                   'Telescope foo is not in known_telescopes.'])
 
     # Test for handling when antenna positions are far enough apart to make the
@@ -522,6 +520,14 @@ def test_readWriteReadMiriad():
 
     nt.assert_equal(uv_in, uv_out)
 
+    # check that we can read & write phased data
+    uv_in2 = copy.deepcopy(uv_in)
+    uv_in2.phase_to_time(Time(np.mean(uv_in2.time_array), format='jd'))
+    uv_in2.write_miriad(write_file, clobber=True)
+    uv_out.read_miriad(write_file)
+
+    nt.assert_equal(uv_in2, uv_out)
+
     # check that trying to overwrite without clobber raises an error
     nt.assert_raises(ValueError, uv_in.write_miriad, write_file)
 
@@ -761,7 +767,7 @@ def test_ReadMiriadPhase():
     phased_uv = UVData()
     # test that phasing makes files equal
     uvtest.checkWarnings(unphased.read, [unphasedfile, 'miriad'], known_warning='miriad')
-    unphased.phase(ra=0.0, dec=0.0, epoch=ephem.J2000)
+    unphased.phase(ra=0.0, dec=0.0, epoch='J2000')
     uvtest.checkWarnings(phased.read, [phasedfile, 'miriad'], known_warning='miriad')
     nt.assert_equal(unphased, phased)
 '''
@@ -809,3 +815,37 @@ def test_antpos_units():
     aantpos = (uvutils.ECEF_from_rotECEF(aantpos, uv.telescope_location_lat_lon_alt[1])
                - uv.telescope_location)
     nt.assert_true(np.allclose(aantpos, uv.antenna_positions))
+
+
+def test_readMiriadwriteMiraid_check_time_format():
+    """
+    test time_array is converted properly from Miriad format
+    """
+    # test read-in
+    fname = os.path.join(DATA_PATH, 'zen.2457698.40355.xx.HH.uvcA')
+    uvd = UVData()
+    uvd.read_miriad(fname)
+    uvd_t = uvd.time_array.min()
+    uvd_l = uvd.lst_array.min()
+    uv = aipy_extracts.UV(fname)
+    uv_t = uv['time'] + uv['inttime'] / (24 * 3600.) / 2
+
+    lat, lon, alt = uvd.telescope_location_lat_lon_alt
+    t1 = Time(uv['time'], format='jd', location=(lon, lat))
+    dt = TimeDelta(uv['inttime'] / 2, format='sec')
+    t2 = t1 + dt
+    delta_lst = t2.sidereal_time('apparent').radian - t1.sidereal_time('apparent').radian
+    uv_l = uv['lst'] + delta_lst
+
+    # assert starting time array and lst array are shifted by half integration
+    nt.assert_almost_equal(uvd_t, uv_t)
+    nt.assert_almost_equal(uvd_l, uv_l, delta=1e-8)
+    # test write-out
+    fout = os.path.join(DATA_PATH, 'ex_miriad')
+    uvd.write_miriad(fout, clobber=True)
+    # assert equal to original miriad time
+    uv2 = aipy_extracts.UV(fout)
+    nt.assert_almost_equal(uv['time'], uv2['time'])
+    nt.assert_almost_equal(uv['lst'], uv2['lst'])
+    if os.path.exists(fout):
+        shutil.rmtree(fout)