File: analysis.rst

package info (click to toggle)
astropy 5.2.1-2
  • links: PTS, VCS
  • area: main
  • in suites: bookworm
  • size: 41,972 kB
  • sloc: python: 219,331; ansic: 147,297; javascript: 13,556; lex: 8,496; sh: 3,319; xml: 1,622; makefile: 185
file content (301 lines) | stat: -rw-r--r-- 9,940 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
.. _timeseries-analysis:

Manipulation and Analysis of Time Series
****************************************

Combining Time Series
=====================

The :func:`~astropy.table.vstack` and :func:`~astropy.table.hstack` functions
from the :mod:`astropy.table` module can be used to stack time series in
different ways.

Examples
--------

.. EXAMPLE START: Stacking Time Series Row-Wise Using table.vstack

Time series can be stacked "vertically" or row-wise using the
:func:`~astropy.table.vstack` function (although note that sampled time
series cannot be combined with binned time series and vice versa)::

    >>> from astropy.table import vstack
    >>> from astropy import units as u
    >>> from astropy.timeseries import TimeSeries
    >>> ts_a = TimeSeries(time_start='2016-03-22T12:30:31',
    ...                   time_delta=3 * u.s,
    ...                   data={'flux': [1, 4, 5, 3, 2] * u.mJy})
    >>> ts_b = TimeSeries(time_start='2016-03-22T12:50:31',
    ...                   time_delta=3 * u.s,
    ...                   data={'flux': [4, 3, 1, 2, 3] * u.mJy})
    >>> ts_ab = vstack([ts_a, ts_b])
    >>> ts_ab
    <TimeSeries length=10>
              time            flux
                              mJy
              Time          float64
    ----------------------- -------
    2016-03-22T12:30:31.000     1.0
    2016-03-22T12:30:34.000     4.0
    2016-03-22T12:30:37.000     5.0
    2016-03-22T12:30:40.000     3.0
    2016-03-22T12:30:43.000     2.0
    2016-03-22T12:50:31.000     4.0
    2016-03-22T12:50:34.000     3.0
    2016-03-22T12:50:37.000     1.0
    2016-03-22T12:50:40.000     2.0
    2016-03-22T12:50:43.000     3.0

Note that :func:`~astropy.table.vstack` does not automatically sort, nor get rid
of duplicates — this is something you would need to do explicitly afterwards.

.. EXAMPLE END

.. EXAMPLE START: Stacking Time Series Column-Wise Using table.vstack

Time series can also be combined "horizontally" or column-wise with other tables
using the :func:`~astropy.table.hstack` function, though these should not be
time series (as having multiple time columns would be confusing)::

    >>> from astropy.table import Table, hstack
    >>> data = Table(data={'temperature': [40., 41., 40., 39., 30.] * u.K})
    >>> ts_a_data = hstack([ts_a, data])
    >>> ts_a_data
    <TimeSeries length=5>
              time            flux  temperature
                              mJy          K
              Time          float64    float64
    ----------------------- ------- -----------
    2016-03-22T12:30:31.000     1.0        40.0
    2016-03-22T12:30:34.000     4.0        41.0
    2016-03-22T12:30:37.000     5.0        40.0
    2016-03-22T12:30:40.000     3.0        39.0
    2016-03-22T12:30:43.000     2.0        30.0

.. EXAMPLE END

Sorting Time Series
===================

.. EXAMPLE START: Sorting Time Series

Sorting time series in place can be done using the
:meth:`~astropy.table.Table.sort` method, as for |Table|::

    >>> ts = TimeSeries(time_start='2016-03-22T12:30:31',
    ...                 time_delta=3 * u.s,
    ...                 data={'flux': [1., 4., 5., 3., 2.]})
    >>> ts
    <TimeSeries length=5>
              time            flux
              Time          float64
    ----------------------- -------
    2016-03-22T12:30:31.000     1.0
    2016-03-22T12:30:34.000     4.0
    2016-03-22T12:30:37.000     5.0
    2016-03-22T12:30:40.000     3.0
    2016-03-22T12:30:43.000     2.0
    >>> ts.sort('flux')
    >>> ts
    <TimeSeries length=5>
              time            flux
              Time          float64
    ----------------------- -------
    2016-03-22T12:30:31.000     1.0
    2016-03-22T12:30:43.000     2.0
    2016-03-22T12:30:40.000     3.0
    2016-03-22T12:30:34.000     4.0
    2016-03-22T12:30:37.000     5.0

.. EXAMPLE END

Resampling
==========

We provide a :func:`~astropy.timeseries.aggregate_downsample` function
that can be used to bin values from a time series into equal-size or uneven bins,
and contiguous and non-contiguous bins, using a custom function (mean, median, etc.).
This operation returns a |BinnedTimeSeries|. Note that this is a basic function in
the sense that it does not, for example, know how to treat columns with uncertainties
differently from other values, and it will blindly apply the custom function
specified to all columns.

Example
-------

.. EXAMPLE START: Creating a BinnedTimeSeries with even contiguous bins

The following example shows how to use
:func:`~astropy.timeseries.aggregate_downsample` to bin a light curve from the
Kepler mission into 20 minute contiguous bins using a median function. First,
we read in the data using:

.. plot::
   :include-source:
   :context: reset
   :nofigs:

    from astropy.timeseries import TimeSeries
    from astropy.utils.data import get_pkg_data_filename
    example_data = get_pkg_data_filename('timeseries/kplr010666592-2009131110544_slc.fits')
    kepler = TimeSeries.read(example_data, format='kepler.fits')

(See :ref:`timeseries-io` for more details about reading in data). We can then
downsample using:

.. plot::
   :context:
   :nofigs:

   import warnings
   warnings.filterwarnings('ignore', message='All-NaN slice encountered')

.. plot::
   :include-source:
   :context:
   :nofigs:

    import numpy as np
    from astropy import units as u
    from astropy.timeseries import aggregate_downsample
    kepler_binned = aggregate_downsample(kepler, time_bin_size=20 * u.min, aggregate_func=np.nanmedian)

We can take a look at the results:

.. plot::
   :include-source:
   :context:

    import matplotlib.pyplot as plt
    plt.plot(kepler.time.jd, kepler['sap_flux'], 'k.', markersize=1)
    plt.plot(kepler_binned.time_bin_start.jd, kepler_binned['sap_flux'], 'r-', drawstyle='steps-pre')
    plt.xlabel('Julian Date')
    plt.ylabel('SAP Flux (e-/s)')

.. EXAMPLE END

.. EXAMPLE START: Creating a BinnedTimeSeries with uneven contiguous bins

The :func:`~astropy.timeseries.aggregate_downsample` can also be used
to bin the light curve into custom bins. The following example shows
the case of uneven-size contiguous bins:

.. plot::
   :context: reset
   :nofigs:

    import numpy as np
    from astropy import units as u
    import matplotlib.pyplot as plt
    from astropy.timeseries import TimeSeries
    from astropy.timeseries import aggregate_downsample
    from astropy.utils.data import get_pkg_data_filename

    example_data = get_pkg_data_filename('timeseries/kplr010666592-2009131110544_slc.fits')
    kepler = TimeSeries.read(example_data, format='kepler.fits')

    import warnings
    warnings.filterwarnings('ignore', message='All-NaN slice encountered')

.. plot::
   :include-source:
   :context:

    kepler_binned = aggregate_downsample(kepler, time_bin_size=[1000, 125, 80, 25, 150, 210, 273] * u.min,
                                         aggregate_func=np.nanmedian)

    plt.plot(kepler.time.jd, kepler['sap_flux'], 'k.', markersize=1)
    plt.plot(kepler_binned.time_bin_start.jd, kepler_binned['sap_flux'], 'r-', drawstyle='steps-pre')
    plt.xlabel('Julian Date')
    plt.ylabel('SAP Flux (e-/s)')

To learn more about the custom binning functionality in
:func:`~astropy.timeseries.aggregate_downsample`, see
:ref:`timeseries-binned-initializing`.

Folding
=======

.. EXAMPLE START: Phase Folding a Time Series

The |TimeSeries| class has a
:meth:`~astropy.timeseries.TimeSeries.fold` method that can be used to
return a new time series with a relative and folded time axis. This method
takes the period as a :class:`~astropy.units.Quantity`, and optionally takes
an epoch as a :class:`~astropy.time.Time`, which defines a zero time offset:

.. plot::
   :context: reset
   :nofigs:

   import numpy as np
   from astropy import units as u
   import matplotlib.pyplot as plt
   from astropy.timeseries import TimeSeries
   from astropy.utils.data import get_pkg_data_filename

   example_data = get_pkg_data_filename('timeseries/kplr010666592-2009131110544_slc.fits')
   kepler = TimeSeries.read(example_data, format='kepler.fits')

.. plot::
   :include-source:
   :context:

    kepler_folded = kepler.fold(period=2.2 * u.day, epoch_time='2009-05-02T20:53:40')

    plt.plot(kepler_folded.time.jd, kepler_folded['sap_flux'], 'k.', markersize=1)
    plt.xlabel('Time from midpoint epoch (days)')
    plt.ylabel('SAP Flux (e-/s)')

Note that in this example we happened to know the period and midpoint from a
previous periodogram analysis. See the example in :doc:`index` for how you
might do this.

.. EXAMPLE END

Arithmetic
==========

.. EXAMPLE START: Arithmetic with Time Series

Since |TimeSeries| objects are subclasses of |Table|, they naturally support
arithmetic on any of the data columns. As an example, we can take the folded
Kepler time series we have seen in previous examples, and normalize it to the
sigma-clipped median value.

.. plot::
   :context: reset
   :nofigs:

   import numpy as np
   from astropy import units as u
   import matplotlib.pyplot as plt
   from astropy.timeseries import TimeSeries
   from astropy.utils.data import get_pkg_data_filename

   example_data = get_pkg_data_filename('timeseries/kplr010666592-2009131110544_slc.fits')
   kepler = TimeSeries.read(example_data, format='kepler.fits')
   kepler_folded = kepler.fold(period=2.2 * u.day, epoch_time='2009-05-02T20:53:40')

.. plot::
   :context:
   :nofigs:

   import warnings
   warnings.filterwarnings('ignore', message='Input data contains invalid values')

.. plot::
   :include-source:
   :context:

    from astropy.stats import sigma_clipped_stats

    mean, median, stddev = sigma_clipped_stats(kepler_folded['sap_flux'])

    kepler_folded['sap_flux_norm'] = kepler_folded['sap_flux'] / median

    plt.plot(kepler_folded.time.jd, kepler_folded['sap_flux_norm'], 'k.', markersize=1)
    plt.xlabel('Time from midpoint epoch (days)')
    plt.ylabel('Normalized flux')

.. EXAMPLE END