Merge e4fb9bf into 1ac364a

.travis.yml

@@ -29,7 +29,7 @@ install:
 
   # create environment and install dependencies
   - conda install numpy scipy astropy nose pip h5py six
-  - conda install -c conda-forge healpy pyephem python-casacore pycodestyle coveralls
+  - conda install -c conda-forge healpy python-casacore pycodestyle coveralls
   - conda list
 script:
   - python setup.py build_ext --force --inplace

README.md

@@ -51,7 +51,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
@@ -103,14 +102,13 @@ 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, scipy, and optionally h5py, and
-conda-forge for pyephem and optionally casacore-python and healpy (e.g. ```conda install -c conda-forge pyephem```).
+conda-forge for optionally installing python-casacore and healpy (e.g. ```conda install -c conda-forge python-casacore```).
 
 * numpy >= 1.10
 * scipy
 * astropy >= 1.2
-* pyephem
 * h5py (optional: for reading and writing uvh5 format)
-* casacore-python (optional: for CASA measurement set reading functionality)
+* python-casacore (optional: for CASA measurement set reading functionality)
 * healpy (optional: working with beams in HEALPix formats)
 
 ### For CASA measurement set functionality, install python-casacore

docs/tutorial.rst

@@ -255,23 +255,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, SkyCoord
 import os
 import shutil
 import numpy as np
@@ -99,7 +100,7 @@ def read_miriad(self, filepath, antenna_nums=None, ant_str=None, bls=None,
                 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:
@@ -165,7 +166,13 @@ def read_miriad(self, filepath, antenna_nums=None, ant_str=None, bls=None,
             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:
@@ -335,7 +342,9 @@ def read_miriad(self, filepath, antenna_nums=None, ant_str=None, bls=None,
             # 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)
 
@@ -468,12 +477,28 @@ def read_miriad(self, filepath, antenna_nums=None, ant_str=None, bls=None,
             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
+
+            # use skycoord to simplify calculating sky separations.
+            # Note, this should be done in the TEE frame, which isn't supported by astropy
+            # Frame doesn't really matter, though, because this is just geometrical, so use icrs
+            pointing_coords = SkyCoord(ra=ra_list, dec=dec_list, unit='radian', frame='icrs')
+            zenith_coord = SkyCoord(ra=self.lst_array,
+                                    dec=self.telescope_location_lat_lon_alt[0],
+                                    unit='radian', frame='icrs')
+
+            separation_angles = pointing_coords.separation(zenith_coord)
+            acceptable_offset = Angle('1d')
+            if (np.max(separation_angles.rad) > acceptable_offset.rad):
+                warnings.warn('drift RA, Dec is off from lst, 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()
@@ -504,6 +529,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)
@@ -622,6 +653,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')
@@ -726,18 +760,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 '
@@ -852,7 +886,7 @@ def _load_miriad_variables(self, uv):
                                     '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