Python astropy.time.Time() Examples

def cent(self, currentTime):
        """Finds time in Julian centuries since J2000 epoch
        
        This quantity is needed for many algorithms from Vallado 2013.
        
        Args:
            currentTime (astropy Time array):
                Current absolute mission time in MJD
            
        Returns:
            float ndarray:
                time in Julian centuries since the J2000 epoch 
        
        """
        
        j2000 = Time(2000., format='jyear',scale='tai')
        TDB = (currentTime.jd - j2000.jd)/36525.
        
        return TDB 

def eclip2equat(self, r_eclip, currentTime):
        """Rotates heliocentric coordinates from ecliptic to equatorial frame.
        
        Args:
            r_eclip (astropy Quantity nx3 array):
                Positions vector in heliocentric ecliptic frame in units of AU
            currentTime (astropy Time array):
                Current absolute mission time in MJD
        
        Returns:
            r_equat (astropy Quantity nx3 array):
                Positions vector in heliocentric equatorial frame in units of AU
        
        """
        
        r_equat = self.equat2eclip(r_eclip, currentTime, rotsign=-1)
        
        return r_equat 

def set_time(self, time):
        """Convert Time object to BCD timestamp elements.

        Parameters
        ----------
        time : `~astropy.time.Time`
            The time to use for this header.
        """
        old_precision = time.precision
        try:
            time.precision = 5
            yday = time.yday.split(':')
        finally:
            time.precision = old_precision
        # Set fraction first since that checks precision.
        self.fraction = float(yday[4]) % 1
        self.decade = int(yday[0][:3]) * 10
        self['bcd_unit_year'] = int(yday[0][3], base=16)
        self['bcd_day'] = int(yday[1], base=16)
        self['bcd_hour'] = int(yday[2], base=16)
        self['bcd_minute'] = int(yday[3], base=16)
        self['bcd_second'] = int(yday[4][:2], base=16) 

def test_orbit(self):
        """
        Test orbit method to ensure proper output
        """

        t_ref = Time(2000.0, format='jyear')
        for mod in self.allmods:
            with RedirectStreams(stdout=self.dev_null):
                if 'SotoStarshade' in mod.__name__:
                    obj = mod(f_nStars=4, **copy.deepcopy(self.spec))
                else:
                    obj = mod(**copy.deepcopy(self.spec))
            r_sc = obj.orbit(t_ref)
            # the r_sc attribute is set and is a 3-tuple of astropy Quantity's
            self.assertEqual(type(r_sc), type(1.0 * u.km))
            self.assertEqual(r_sc.shape, (1, 3)) 

def test_stream(self, tmpdir):
        with vdif.open(SAMPLE_VDIF, 'rs') as fr:
            vh = fr.header0
            data = fr.read(20000)  # enough to fill two Mark 5B frames.

        fl = str(tmpdir.join('test.m5b'))
        with mark5b.open(fl, 'ws', sample_rate=vh.sample_rate,
                         nchan=data.shape[1], bps=vh.bps, time=vh.time) as fw:
            fw.write(data)

        with vdif.open(SAMPLE_VDIF, 'rs') as fv, mark5b.open(
                fl, 'rs', sample_rate=32.*u.MHz,
                ref_time=Time(57000, format='mjd'), nchan=8, bps=2) as fm:
            assert fv.header0.time == fm.header0.time
            dv = fv.read(20000)
            dm = fm.read(20000)
            assert np.all(dm == dv)
            assert fm.offset == fv.offset
            assert fm.tell(unit='time') == fv.tell(unit='time') 

def test_rawdump_header_seek_offset(self):
        fh = open(SAMPLE_RAWDUMP_HEADER, 'rt')

        header = gsb.GSBHeader.fromfile(fh, verify=True)
        # Includes 1 trailing blank space, one line separator
        header_nbytes = header.nbytes
        assert (header_nbytes
                == len(' '.join(header.words) + ' ' + os.linesep))
        assert header.seek_offset(1) == header_nbytes
        assert header.seek_offset(12) == 12 * header_nbytes

        # Note that text pointers can't seek from current position.
        # Seek 2nd header.
        fh.seek(header.seek_offset(1))
        header1 = gsb.GSBHeader.fromfile(fh, verify=True)
        assert abs(header1.time
                   - Time('2015-04-27T13:15:00.251658480')) < 1.*u.ns

        fh.seek(header.seek_offset(9))
        header2 = gsb.GSBHeader.fromfile(fh, verify=True)
        assert abs(header2.time
                   - Time('2015-04-27T13:15:02.264924400')) < 1.*u.ns

        fh.close() 

def test_rawdump_header(self):
        with open(SAMPLE_RAWDUMP_HEADER, 'rt') as fh:
            header = gsb.GSBHeader.fromfile(fh, verify=True)
        assert header.mode == 'rawdump'
        assert header['gps'] == '2015 04 27 18 45 00 0.000000240'
        # Includes UTC offset.
        assert abs(header.time
                   - Time('2015-04-27T13:15:00.000000240')) < 1.*u.ns
        header2 = gsb.GSBHeader.fromkeys(**header)
        assert header2 == header
        header3 = gsb.GSBHeader.fromvalues(mode='rawdump', **header2)
        assert header3 == header2
        with pytest.raises(TypeError):
            gsb.GSBHeader.fromvalues(**header)
        with pytest.raises(TypeError):
            gsb.GSBHeader(None)
        # Check that recovering with the actual header type works as well.
        header4 = type(header).fromkeys(**header)
        assert header4 == header
        # Check that trying to initialize a phased header doesn't work.
        with pytest.raises(KeyError):
            gsb.header.GSBPhasedHeader.fromkeys(**header)
        header5 = header.copy()
        assert header5 == header 

def test_start_at_last_frame(tmpdir):
    """Regression test for files where first frame is last in a second.

    Previously, this caused the frame rate and thus the sample rate to
    be calculated wrongly.  See gh-340.
    """
    fl = str(tmpdir.join('test.m4'))
    sample_rate = 32*u.MHz
    frame_rate = sample_rate / 80000
    start_time = Time('2012-01-02') - 1/frame_rate
    with mark4.open(fl, 'ws', sample_rate=sample_rate, time=start_time,
                    ntrack=32, sample_shape=(4,), fanout=4, bps=2) as fw:
        fw.write(np.ones((80000*2, 4)))

    with mark4.open(fl, 'rs', decade=2010) as fr:
        assert fr.sample_rate == sample_rate 

def time(self, time):
        """Set header time.

        If ``start_time`` is not already set, this sets it using the time and
        ``offset``.  Otherwise, this sets ``offset`` using the time and
        ``start_time``.

        Parameters
        ----------
        time : `~astropy.time.Time`
            Time for the first sample associated with this header.
        """
        if 'MJD_START' not in self.keys():
            self.start_time = time - self.offset
        else:
            self.offset = time - self.start_time 

def set_jds(self, val1, val2):
        """Parse the time strings contained in val1 and set jd1, jd2"""
        iterator = np.nditer([val1, None, None, None, None, None, None],
                             op_dtypes=([val1.dtype] + 5 * [np.intc]
                                        + [np.double]))
        try:
            for val, iy, im, id, ihr, imin, dsec in iterator:
                timestr = val.item()
                components = timestr.split()
                iy[...], im[...], id[...], ihr[...], imin[...], sec = (
                    int(component) for component in components[:-1])
                dsec[...] = sec + float(components[-1])
        except Exception:
            raise ValueError('Time {0} does not match {1} format'
                             .format(timestr, self.name))

        self.jd1, self.jd2 = erfa.dtf2d(
            self.scale.upper().encode('utf8'), *iterator.operands[1:]) 

def make_cphase_df(obs,band='unknown',polarization='unknown',mode='all',count='max',round_s=0.1,snrcut=0.,uv_min=False):

    """generate DataFrame of closure phases

    Args: 
        obs: ObsData object
        round_s: accuracy of datetime object in seconds

    Returns:
        df: closure phase data in DataFrame format
    """

    data=obs.c_phases(mode=mode,count=count,snrcut=snrcut,uv_min=uv_min)
    sour=obs.source
    df = pd.DataFrame(data=data).copy()
    df['fmjd'] = df['time']/24.
    df['mjd'] = obs.mjd + df['fmjd']
    df['triangle'] = list(map(lambda x: x[0]+'-'+x[1]+'-'+x[2],zip(df['t1'],df['t2'],df['t3'])))
    df['datetime'] = Time(df['mjd'], format='mjd').datetime
    df['datetime'] =list(map(lambda x: round_time(x,round_s=round_s),df['datetime']))
    df['jd'] = Time(df['mjd'], format='mjd').jd
    df['polarization'] = polarization
    df['band'] = band
    df['source'] = sour
    return df 

def time(self, time):
        """Set header time.

        If ``start_time`` is not already set, this sets it using the time and
        ``offset``.  Otherwise, this sets ``offset`` using the time and
        ``start_time``.

        Parameters
        ----------
        time : `~astropy.time.Time`
            Time for the first sample associated with this header.
        """
        if 'STT_IMJD' not in self.keys():
            self.start_time = time - self.offset
        else:
            self.offset = time - self.start_time 

def test_moon_earth(self):
        r"""Test moon_earth method.

        Approach: Reference to pre-computed result from Matlab.
        """

        print('moon_earth()')
        obs = self.fixture
        # TODO: add other times besides this one
        # century = 0
        t_ref_string = '2000-01-01T12:00:00.0'
        t_ref = Time(t_ref_string, format='isot', scale='utc')
        moon = obs.moon_earth(t_ref).flatten().to(u.km)
        # print moon
        r_earth = 6378.137 # earth radius [km], as used in Vallado's code
        moon_ref = [-45.74169421, -41.80825511, -11.88954996] # pre-computed from Matlab
        for coord in range(3):
            self.assertAlmostEqual(moon[coord].value, moon_ref[coord] * r_earth, places=1) 

def test_cent(self):
        r"""Test cent method.

        Approach: Probes for a range of inputs.
        """

        print('cent()')
        obs = self.fixture
        # origin at 12:00 on 2000.Jan.01
        t_ref_string = '2000-01-01T12:00:00.0'
        t_ref = Time(t_ref_string, format='isot', scale='utc')
        self.assertEqual(obs.cent(t_ref), 0.0)

        # even-julian-year probes
        t_probe = np.linspace(1950.0, 2050.0, 101)
        for t_ref in t_probe:
            # get the Time object, required by the cent() method
            t_probe_2 = Time(t_ref, format='jyear')
            # exosims century (offset from 2000)
            t_exo = obs.cent(t_probe_2)
            # reference century
            t_ref_cent = (t_ref - 2000.0) / 100.0
            # they are not exactly equal
            self.assertAlmostEqual(t_exo, t_ref_cent, places=10) 

def start_time(self, start_time):
        start_time = Time(start_time, scale='utc', format='isot', precision=9)
        self['STT_IMJD'] = int(start_time.mjd)
        self['STT_SMJD'] = int(np.around(
            (start_time - Time(self['STT_IMJD'], format='mjd',
                               scale=start_time.scale)).sec)) 

def start_time(self, start_time):
        start_time = Time(start_time, scale='utc', format='isot', precision=9)
        self['UTC_START'] = (start_time.isot.replace('T', '-')
                             .replace('.000000000', ''))
        mjd_int = int(start_time.mjd)
        mjd_frac = (start_time - Time(mjd_int, format='mjd',
                                      scale=start_time.scale)).jd
        if mjd_frac < 0:
            mjd_int -= 1
            mjd_frac += 1.
        self['MJD_START'] = ('{0:05d}'.format(mjd_int)
                             + ('{0:17.15f}'.format(mjd_frac))[1:]) 

def test_phased_write_one_file(self, tmpdir):
        # With a single file, one does have to pass in header_mode
        # to get phased.
        with gsb.open(str(tmpdir.join('test.timstamp')), 'ws',
                      raw=str(tmpdir.join('test.raw')),
                      header_mode='phased',
                      time=Time('2010-10-10')) as fh_right:
            assert fh_right.header0.mode == 'phased'

        with gsb.open(str(tmpdir.join('test.timstamp')), 'ws',
                      raw=str(tmpdir.join('test.raw')),
                      time=Time('2010-10-10')) as fh_wrong:
            assert fh_wrong.header0.mode == 'rawdump' 

def start_time(self):
        """Start time of the observation."""
        mjd_int, frac = self['MJD_START'].split('.')
        mjd_int = int(mjd_int)
        frac = float('.' + frac)
        # replace '-' between date and time with a 'T' and convert to Time
        return Time(mjd_int, frac, scale='utc', format='mjd') 

def test_mwa_vdif_with_sample_rate():
    info1 = file_info(SAMPLE_MWA)
    assert info1.format == 'vdif'
    assert 'frame_rate' in info1.errors
    info2 = file_info(SAMPLE_MWA, sample_rate=1.28*u.MHz)
    assert info2.format == 'vdif'
    assert info2.start_time == Time('2015-10-03T20:49:45.000')
    info3 = file_info(SAMPLE_MWA, sample_rate='bla')
    assert info3.format == 'vdif'
    assert 'stream' in info3.errors 

def test_file_info(sample, format_, used, consistent, inconsistent):
    # Pass on extra arguments needed to get Mark4 and Mark5B to pass.
    # For GSB, we also need raw files, so we test that below.
    extra_args = {'ref_time': Time('2014-01-01'),
                  'nchan': 8}
    info = file_info(sample, **extra_args)
    assert info.format == format_
    info_dict = info()
    for attr in info.attr_names:
        info_value = getattr(info, attr)
        assert info_value is not None
        assert attr in info_dict or info_value == {}
    assert set(info.used_kwargs.keys()) == set(used)
    assert set(info.consistent_kwargs.keys()) == set(consistent)
    assert set(info.inconsistent_kwargs.keys()) == set(inconsistent)
    assert set(info.irrelevant_kwargs.keys()) == set()
    # Check that extraneous arguments get classified correctly.
    info2 = file_info(sample, life=42, **extra_args)
    assert info2.used_kwargs == info.used_kwargs
    assert info2.consistent_kwargs == info.consistent_kwargs
    assert info2.inconsistent_kwargs == info.inconsistent_kwargs
    assert info2.irrelevant_kwargs == {'life': 42}
    # Check we can indeed open a file with the extra arguments.
    module = importlib.import_module('.' + info.format, package='baseband')
    with module.open(sample, mode='rs', **info.used_kwargs) as fh:
        info3 = fh.info
    assert info3() == info_dict
    # Check that things properly do *not* work on a closed file.
    with module.open(sample, mode='rs', **info.used_kwargs) as fh:
        pass
    info4 = fh.info
    assert not info4
    assert 'File closed' in repr(info4)
    assert 'errors' in info4()
    assert any(isinstance(v, ValueError) for v in info4.errors.values()) 

def test_mwa_vdif():
    """Test phased VDIF format (uses EDV=0)"""
    with vdif.open(SAMPLE_MWA, 'rs', sample_rate=1.28*u.MHz) as fh:
        assert fh.samples_per_frame == 128
        assert fh.sample_rate == 1.28*u.MHz
        assert fh.time == Time('2015-10-03T20:49:45.000')
        assert fh.header0.edv == 0 

def get_time(self):
        """Convert BCD time code to Time object.

        Calculate time using bcd-encoded 'bcd_unit_year', 'bcd_day',
        'bcd_hour', 'bcd_minute', 'bcd_second' header items, as well as
        the ``fraction`` property (inferred from 'bcd_fraction') and
        ``decade`` from the initialisation.  See
        See https://www.haystack.mit.edu/tech/vlbi/mark5/docs/230.3.pdf
        """
        return Time('{decade:03d}{uy:1x}:{d:03x}:{h:02x}:{m:02x}:{s:08.5f}'
                    .format(decade=self.decade//10, uy=self['bcd_unit_year'],
                            d=self['bcd_day'], h=self['bcd_hour'],
                            m=self['bcd_minute'],
                            s=bcd_decode(self['bcd_second']) + self.fraction),
                    format='yday', scale='utc', precision=5) 

def infer_decade(self, ref_time):
        """Uses a reference time to set a header's ``decade``.

        Parameters
        ----------
        ref_time : `~astropy.time.Time`
            Reference time within 5 years of the observation time.
        """
        self.decade = np.around(ref_time.decimalyear - self['bcd_unit_year'],
                                decimals=-1).astype(int) 

def test_frame(self):
        """Check a whole Mark 5B frame can be translated to VDIF."""
        with mark5b.open(SAMPLE_M5B, 'rb', ref_time=Time(57000, format='mjd'),
                         nchan=8, bps=2) as fh:
            # pick second frame just to be different from header checks above.
            fh.seek(10016)
            m5f = fh.read_frame()

        assert m5f['frame_nr'] == 1
        frame = vdif.VDIFFrame.from_mark5b_frame(m5f)
        assert frame.nbytes == 10032
        assert frame.shape == (5000, 8)
        assert np.all(frame.data == m5f.data)
        assert frame.time == m5f.time 

def update(self, *, time=None, frame_rate=None, crc=None, verify=True,
               **kwargs):
        """Update the header by setting keywords or properties.

        Here, any keywords matching header keys are applied first, and any
        remaining ones are used to set header properties, in the order set
        by the class (in ``_properties``).

        Parameters
        ----------
        time : `~astropy.time.Time`, optional
            A possible new time.  This is updated last.
        frame_rate : `~astropy.units.Quantity`, optional
            The frame rate to use in calculating the frame number from the
            time.  Needed for times at non-integer seconds.
        crc : int or None, optional
            If `None` (default), recalculate the CRC after updating.
        verify : bool, optional
            If `True` (default), verify integrity after updating.
        **kwargs
            Arguments used to set keywords and properties.
        """
        super().update(verify=False, **kwargs)

        if time is not None:
            self.set_time(time, frame_rate=frame_rate)

        if crc is None:
            # Do not use words 2 & 3 directly, so that this works also if part
            # of a VDIF header, where the time information is in words 7 & 8.
            stream = ((((self['bcd_jday'] << 20)
                        + self['bcd_seconds']) << 16)
                      + self['bcd_fraction'])
            crc = crc16(stream)

        self['crc'] = crc
        if verify:
            self.verify() 

def test_sequentialfile_vdif_stream(tmpdir):
    vdif_sequencer = Sequencer(str(tmpdir.join('{:07d}.vdif')))
    # Try writing a very simple file, using edv=0.
    data = np.ones((16, 16, 2, 2))
    for i, dat in enumerate(data):
        dat[i, 0, 0] = -1.
        dat[i, 1, 1] = -1.
    data.shape = -1, 2, 2
    # construct first header
    header = vdif.VDIFHeader.fromvalues(
        edv=0, time=Time('2010-01-01'), nchan=2, bps=2,
        complex_data=False, frame_nr=0, thread_id=0, samples_per_frame=16,
        station='me')
    with sequentialfile.open(vdif_sequencer, 'w+b',
                             file_size=4*header.frame_nbytes) as sfh, \
            vdif.open(sfh, 'ws', header0=header, nthread=2,
                      sample_rate=256*u.Hz) as fw:
        fw.write(data)
    # check that this wrote 8 frames
    files = [vdif_sequencer[i] for i in range(8)]
    for file_ in files:
        assert os.path.isfile(file_)
    assert not os.path.isfile(vdif_sequencer[8])

    with sequentialfile.open(vdif_sequencer, 'rb') as sfh, vdif.open(
            sfh, 'rs', sample_rate=256*u.Hz) as fr:
        record1 = fr.read(21)
        assert np.all(record1 == data[:21])
        fr.seek(7*16)
        record2 = fr.read(61)
        assert np.all(record2 == data[7*16:7*16+61])
        assert fr.tell() == 7*16+61
    # Might as well check file list too.
    with sequentialfile.open(files, 'rb') as sfh, vdif.open(
            sfh, 'rs', sample_rate=256*u.Hz) as fr:
        record1 = fr.read()
        assert np.all(record1 == data) 

def test_header_non_gmrt(self):
        # Assume file is non-GMRT.
        with open(SAMPLE_PHASED_HEADER, 'rt') as fh:
            header = gsb.GSBHeader.fromfile(fh, verify=True,
                                            utc_offset=0.*u.hr)
        assert abs(header.pc_time
                   - Time('2013-07-28T02:53:55.517535')) < 1.*u.ns
        assert header.gps_time == header.time
        assert abs(header.time
                   - Time('2013-07-28T02:53:55.3241088')) < 1.*u.ns 

def pc_time(self):
        return Time(self['pc'], format='gsb',
                    precision=self._pc_time_precision) - self.utc_offset 

def gps_time(self):
        return Time(self['gps'], format='gsb',
                    precision=self._gps_time_precision) - self.utc_offset 

def setup_class(cls):
        cls.header0 = vdif.VDIFHeader.fromvalues(
            edv=1, time=Time('2010-11-12T13:14:15'), nchan=2, bps=2,
            complex_data=False, thread_id=0, samples_per_frame=16,
            station='me', sample_rate=2*u.kHz)
        cls.nthread = 2
        cls.data = np.array([[[-1, 1],
                              [-3, 3]]]*16)
        cls.frameset_nbytes = cls.header0.frame_nbytes * cls.nthread