Adds Power Spectrum example for a timeseries #2496

yashkgp · 2018-03-02T19:13:30Z

Issue #2483
Here is an example :

import sunpy.timeseries
import sunpy.data.sample
import matplotlib.pyplot as plt
import pandas as pd 
from scipy import signal
ts = sunpy.timeseries.TimeSeries(sunpy.data.sample.GOES_XRS_TIMESERIES)
k,l=ts.power_spectra()
plt.semilogy(k[0], l[0])
plt.show()

Apart from adding the unit tests,is anything else needed for this method?

power_spectra

pep8speaks · 2018-03-02T19:13:32Z

Hello @yashkgp! Thanks for updating the PR.

In the file examples/power_spectra_example.py, following are the PEP8 issues :

Line 12:32: W291 trailing whitespace
Line 21:79: W291 trailing whitespace
Line 22:11: W291 trailing whitespace
Line 27:65: W291 trailing whitespace
Line 29:79: W291 trailing whitespace

In the file sunpy/timeseries/timeseriesbase.py, following are the PEP8 issues :

Line 138:9: E265 block comment should start with '# '
Line 139:9: E265 block comment should start with '# '
Line 184:13: E221 multiple spaces before operator
Line 216:13: E221 multiple spaces before operator
Line 217:13: E221 multiple spaces before operator
Line 246:101: E501 line too long (138 > 100 characters)
Line 277:16: E221 multiple spaces before operator
Line 281:16: E221 multiple spaces before operator
Line 318:101: E501 line too long (107 > 100 characters)
Line 324:101: E501 line too long (127 > 100 characters)
Line 435:101: E501 line too long (113 > 100 characters)
Line 507:33: E231 missing whitespace after ':'
Line 604:101: E501 line too long (110 > 100 characters)

Comment last updated on April 12, 2018 at 11:42 Hours UTC

dstansby

What happens if the timeseries data isn't on continuously spaced timestamps?

dpshelio · 2018-03-07T17:56:31Z

sunpy/timeseries/timeseriesbase.py

+            data_ = self.data[col]
+            freq, spectra = signal.periodogram(data_, fs)
+            Freqs.append(freq)
+            Spectra.append(spectra)


from functools import partial periodogram_fs = partial(signal.periodogram, fs=fs) Freqs, Spectra = zip(*map(periodogram_fs, self.data))

This will remove all the loop, it will produce tuples instead of lists though.

dpshelio · 2018-03-07T17:56:52Z

sunpy/timeseries/timeseriesbase.py

+
+        Returns
+        -------
+        Freqs : array of arrays


You are not returning array but a list of arrays.

dpshelio · 2018-03-07T17:57:29Z

sunpy/timeseries/timeseriesbase.py

+
+        Parameter
+        ----------
+        fs : `float, optional`


optional should be outside the backtick.

nabobalis · 2018-03-07T18:24:37Z

Hey @yashkgp, sorry for the late reply.

We aren't sure we want to add methods on to time series right now. It sparked quite a discussion about methods in various parts of SunPy. That shake the foundation of SunPy itself.

What we want from this now is if you could turn it into an example we could add to our Gallery? If you have a look at the other examples here for how we format them!

yashkgp · 2018-03-07T18:32:36Z

Sure @nabobalis Not a Problem .

…ectra

yashkgp · 2018-03-12T16:18:28Z

Hey @nabobalis Could you give it a review?

nabobalis

I think there needs to be more explanations and warnings for people who read the example.

nabobalis · 2018-03-12T16:23:32Z

examples/power_spectra_example.py

@@ -0,0 +1,43 @@
+"""
+===============================
+POWER SPECTRUM OF A TIMESERIES


Not all caps, just normal case.

nabobalis · 2018-03-12T16:24:10Z

examples/power_spectra_example.py

+POWER SPECTRUM OF A TIMESERIES
+===============================
+
+An example showing how to estimate power spectrum of a TimeSeries


I think a the should between estimate and power. Full stop at the end.

nabobalis · 2018-03-12T16:25:22Z

examples/power_spectra_example.py

+###############################################################################
+# Let's first load a TimeSeries from sample data
+
+ts = sunpy.timeseries.TimeSeries(GOES_XRS_TIMESERIES)


You might want to describe the data you will working on to give whoever reads this an idea of what the input is.

I have added a brief description about the input data now . Is it okay ?

nabobalis · 2018-03-12T16:26:24Z

examples/power_spectra_example.py

+ts = sunpy.timeseries.TimeSeries(GOES_XRS_TIMESERIES)
+
+###############################################################################
+# We now use scipy's periodogram to estimate the power spectra of the first 


You want to explain why you use the periodogram instead of say the FFT. Also we probably want to mention the limitations or requirements of using these methods. We would like people not to assume they can use these methods on any signal.

I think there is one in Astropy as well now. But I am unsure of what version.

Yeah , just discovered that . Astropy's Periodogram has some amazing functionalities . We should probably switch over to Astropy's Periodogram .

nabobalis · 2018-03-12T16:27:53Z

examples/power_spectra_example.py

+# column of the Timeseries
+
+cols=ts.columns
+freq, spectra = signal.periodogram(ts.data[cols[0]], fs=1)


Is fs=1 correct for the data?

Also I wonder if we could use the false alarm probability in this example as well.

This feature is covered in Astropy's periodogram . Should we switch over to Astropy's periodogram ?

If it works out regarding the version of Astropy it is supported in.

Yeah , it's there in Astropy 3.0

Well we still have to support Astropy 2 until our next release.

nabobalis · 2018-03-12T16:28:09Z

examples/power_spectra_example.py

+
+
+
+


All these extra lines need to go.

nabobalis · 2018-03-12T16:32:53Z

examples/power_spectra_example.py

+
+plt.semilogy(freq, spectra)
+plt.title('Power Spectrum of {}'.format(cols[0]))
+plt.ylabel('Power Spectral Density')


Missing the unit

yashkgp · 2018-03-12T17:47:54Z

Hey @nabobalis , Just discovered that Astropy's Periodogram is much more powerful than scipy's Periodogram . They have full fledged tutorial dedicated to Periodogram http://docs.astropy.org/en/stable/stats/lombscargle.html#.It has various functionalities like frequency grid spacing, false alarm probabilities, normalization etc.We should definitely switch over to Astrpoy's peiodogram . What do you say @nabobalis ?

dpshelio · 2018-03-12T20:01:23Z

examples/power_spectra_example.py

+# We now use scipy's periodogram to estimate the power spectra of the first 
+# column of the Timeseries
+
+cols=ts.columns


Since probably these examples need to be somehow meaningful I would do the following instead of using a variable called cols:

xrays_short = ts.columns[0]

and change cols[0] for xrays_short everywhere. That makes the code more readable.

…ectra

yashkgp · 2018-03-14T10:57:24Z

Hey @nabobalis @dpshelio , I tried switching to Astropy's version of periodogram (yeah, its supported in Astropy 2 ) but its throwing a ValueError

TypeError: ufunc true_divide cannot use operands with types dtype('float64') and dtype('<m8[ns]')

It seems like Astropy's periodogram does accept Pandas Timeseries object . One possible solution might to use Timedelta and convert dtype('<m8[ns]') into float . Is there any easy way to solve this problem ?

nabobalis · 2018-03-15T08:59:10Z

If you can't get the astropy one to work, It would just be worth noting that its available and linking to it in the example.

nabobalis · 2018-03-15T09:01:00Z

examples/power_spectra_example.py

+ts = sunpy.timeseries.TimeSeries(GOES_XRS_TIMESERIES)
+# Data is taken from GOES X-Ray Sensor(XRS), which observes full-disk-integrated
+# solar flux in two broadband channels: 1--8 angstrom (long); and
+# 0.5--4 angstrom (short)


This should be moved above the python code.

But what I wanted here was not what is GOES, it was can we use signal analysis on this data. Is it evenly spaced in time, whats its time step, are there any problems that might mean that signal analysis would be tricky etc. We don't want to give the impression you can just singal analysis any signal.

Data is not completely evenly sampled,it has 2 outliers .Should I mention this?

nabobalis · 2018-03-15T09:01:24Z

examples/power_spectra_example.py

+
+###############################################################################
+# We now use scipy's periodogram to estimate the power spectra of the first 
+# column of the Timeseries


Make a link to the astropy version here, saying it's available instead.

dpshelio · 2018-03-15T12:02:16Z

examples/power_spectra_example.py

-cols=ts.columns
-freq, spectra = signal.periodogram(ts.data[cols[0]], fs=1)
+xray_short = ts.columns[0]
+freq, spectra = signal.periodogram(ts.data[xray_short], fs=1)


I've looked in detail into the data and I don't believe this is correct. fs=1 means the cadence of the data in Hz, the data seems to not be completely evenly sampled!!

>>> sample = np.array([(x - y).total_seconds() for x, y in zip(ts.index[1:], ts.index)]) >>> sample.min(), sample.max(), sample.mean(), sample.std() (2.0430000000000001, 14.337000000000002, 2.0485031771623676, 0.077912702228747382)

I don't know how there's a 14.33 seconds gap!!! If you plot that you see it happens twice! This means that this analysis would be wrong all together!

In case the cadence would be 2.043, then I think the fs should be 1/2.043 s = 0.48 Hz

Should I clean the data by removing those 2 outliers ?

You cannot clean that by removing them... you would have bigger gaps! The only way to "clean" it is to resample the data.

Should I just change fs to 0.48 and notify the users about the irregularities in data or should I change the input data ?

dpshelio · 2018-03-15T12:05:35Z

examples/power_spectra_example.py

-plt.title('Power Spectrum of {}'.format(cols[0]))
-plt.ylabel('Power Spectral Density')
+plt.title('Power Spectrum of {}'.format(xray_short))
+plt.ylabel('Power Spectral Density(({})**2/Hz)'.format(ts.units[xray_short]))


You could use units here as:

import astropy.units as u plt.ylabel('Power Spectral Density ({:LaTeX})'.format(ts.units[xray_short] ** 2 / u.Hz))

Also, put a space between the Density and the ()

I believe square brackets for the units are more common: [Hz] instead of (Hz)

Yeah, this is much cleaner.

dpshelio · 2018-03-15T12:09:46Z

examples/power_spectra_example.py

-plt.title('Power Spectrum of {}'.format(cols[0]))
-plt.ylabel('Power Spectral Density')
+plt.title('Power Spectrum of {}'.format(xray_short))
+plt.ylabel('Power Spectral Density(({})**2/Hz)'.format(ts.units[xray_short]))
 plt.xlabel('Frequency(Hz)')


Also space between the title and the units.

…ectra

yashkgp · 2018-03-19T10:38:35Z

Hey @nabobalis @dpshelio , I have changed the input data, as the previous one was not evenly sampled . I have switched to RHESSI_TIMESERIES , which is a evenly sampled timeseries with a time step of 4 sec. Could you please give it a review?

nabobalis · 2018-03-19T10:51:30Z

examples/power_spectra_example.py

+import matplotlib.pyplot as plt
+import pandas as pd 
+from scipy import signal
+


I would remove this blank line.

nabobalis · 2018-03-19T10:51:44Z

examples/power_spectra_example.py

+
+###############################################################################
+# Let's first load a TimeSeries from sample data.
+


Remove this line.

nabobalis · 2018-03-19T10:52:33Z

examples/power_spectra_example.py

+ts = sunpy.timeseries.TimeSeries(RHESSI_TIMESERIES)
+
+###############################################################################
+# We now use scipy's periodogram to estimate the power spectra of the first 


SciPy, also what is the first column?

nabobalis · 2018-03-19T10:52:52Z

examples/power_spectra_example.py

+###############################################################################
+# We now use scipy's periodogram to estimate the power spectra of the first 
+# column of the Timeseries. Astropy's Lomb-Scargle Periodograms can be used
+# instead of scipy's periodogram.


nabobalis · 2018-03-19T10:53:14Z

examples/power_spectra_example.py

+
+###############################################################################
+# We now use scipy's periodogram to estimate the power spectra of the first 
+# column of the Timeseries. Astropy's Lomb-Scargle Periodograms can be used


Link to Astropy's documentation for it.

nabobalis · 2018-03-19T10:53:30Z

examples/power_spectra_example.py

+x_ray = ts.columns[0]
+# The suitable value for fs would be 0.25 Hz as the time step is 4 s.
+freq, spectra = signal.periodogram(ts.data[xray_short], fs=0.25)
+###############################################################################


Blank line here above the ####...

nabobalis · 2018-03-19T10:55:04Z

examples/power_spectra_example.py

+
+x_ray = ts.columns[0]
+# The suitable value for fs would be 0.25 Hz as the time step is 4 s.
+freq, spectra = signal.periodogram(ts.data[xray_short], fs=0.25)


Sorry @nabobalis ,had updated in my local branch , but somehow forgot to push it.

nabobalis · 2018-03-19T10:55:11Z

examples/power_spectra_example.py

+# Plot the power spectrum
+
+plt.semilogy(freq, spectra)
+plt.title('Power Spectrum of {}'.format(xray_short))


nabobalis · 2018-03-19T12:25:45Z

examples/power_spectra_example.py

+# Let's first load a TimeSeries from sample data.
+# Input data contains 9 columns, which are evenly sampled with a time
+# step of 4 seconds. So, the data is suitable for signal analysis.
+ts = sunpy.timeseries.TimeSeries(RHESSI_TIMESERIES)


Line break above here just to be consistent.

nabobalis · 2018-03-19T12:27:40Z

examples/power_spectra_example.py

+from sunpy.data.sample import RHESSI_TIMESERIES
+
+###############################################################################
+# Let's first load a TimeSeries from sample data.


Lets first load a RHESSI TimeSeries from SunPy's sample data

nabobalis · 2018-03-19T12:30:41Z

examples/power_spectra_example.py

+
+###############################################################################
+# Let's first load a TimeSeries from sample data.
+# Input data contains 9 columns, which are evenly sampled with a time


This data contains

nabobalis · 2018-03-19T12:35:49Z

examples/power_spectra_example.py

+###############################################################################
+# Let's first load a TimeSeries from sample data.
+# Input data contains 9 columns, which are evenly sampled with a time
+# step of 4 seconds. So, the data is suitable for signal analysis.


Remove the last sentence.

nabobalis · 2018-03-19T12:37:18Z

examples/power_spectra_example.py

+ts = sunpy.timeseries.TimeSeries(RHESSI_TIMESERIES)
+
+###############################################################################
+# We now use SciPy's periodogram to estimate the power spectra of the first 


I think that ~scipy.signal.lombscargle should work as a link to the documentation

Ohh! Thats quite useful. Thanks a lot @nabobalis

nabobalis · 2018-03-19T12:38:50Z

examples/power_spectra_example.py

+# We now use SciPy's periodogram to estimate the power spectra of the first 
+# column of the Timeseries. The first column contains X-Ray emmisions in the range 
+# of 3-6 keV. Astropy's Lomb-Scargle Periodograms 
+# (http://docs.astropy.org/en/v2.0.5/stats/lombscargle.html) can be used


Try ~astropy.stats.LombScargle instead. Sphinx should link to it since we have that setup.

I would rephrase the last sentence as well. An alterantive version is .... or something.

dpshelio · 2018-03-19T21:52:23Z

examples/power_spectra_example.py

+# Start by importing the necessary modules.
+
+import matplotlib.pyplot as plt
+import pandas as pd 


pandas is not used

dpshelio · 2018-03-19T21:53:45Z

examples/power_spectra_example.py

+ts = sunpy.timeseries.TimeSeries(RHESSI_TIMESERIES)
+
+###############################################################################
+# We now use SciPy's periodogram(~scipy.signal.lombscargle) to estimate the power 


Wait... signal.periodogram is the same than signal.lombscargle? If so, then you could have used the previous series that was not evenly spaced.

@dpshelio Its not the same , I'll just change it in the next commit .

Sorry @dpshelio , I am not that well acquainted with lomb scargle periodogram . Is there any advantage of using signal.lombscargle over signal.periodogram?

They are not the same, that should point to scipy.signal.periodogram instead.

…ectra

nabobalis · 2018-03-23T08:42:27Z

Odd how the sphinx interlink hasn't worked. I wonder what is wrong.

nabobalis · 2018-04-02T12:07:00Z

examples/power_spectra_example.py

+ts = sunpy.timeseries.TimeSeries(RHESSI_TIMESERIES)
+
+###############################################################################
+# We now use SciPy's periodogram(~scipy.signal.periodogram) to estimate the power 


There needs to be ` on each side for these links to work.

nabobalis · 2018-04-02T12:07:05Z

examples/power_spectra_example.py

+# We now use SciPy's periodogram(~scipy.signal.periodogram) to estimate the power 
+# spectra of the first column of the Timeseries. The first column contains X-Ray 
+# emmisions in the range of 3-6 keV. An alterantive version is Astropy's 
+# Lomb-Scargle Periodograms(~astropy.stats.LombScargle).


There needs to be ` on each side for these links to work.

nabobalis

Other than the backticks, it looks good to me.

yashkgp · 2018-04-10T21:03:09Z

Sorry for the late update , had my exams going on.

…ectra

nabobalis · 2018-04-11T08:49:30Z

examples/power_spectra_example.py

+ts = sunpy.timeseries.TimeSeries(RHESSI_TIMESERIES)
+
+###############################################################################
+# We now use SciPy's periodogram(`~scipy.signal.periodogram`) to estimate the 


Looking at the final result, i think we want this to be SciPy's ~scipy.signal.periodogram

nabobalis · 2018-04-11T08:50:43Z

examples/power_spectra_example.py

+# We now use SciPy's periodogram(`~scipy.signal.periodogram`) to estimate the 
+# power spectra of the first column of the Timeseries. The first column contains
+# X-Ray emmisions in the range of 3-6 keV. An alterantive version is Astropy's 
+# Lomb-Scargle Periodograms(`~astropy.stats.LombScargle`).


Similar here. Astropy's ~astropy.stats.LombScargle.

nabobalis

Need two minor changes

…ectra

nabobalis · 2018-04-13T07:35:24Z

Thanks @yashkgp!

yashkgp added 2 commits March 3, 2018 00:30

Implements power spectra for TimeSeries

6349469

Merge branch 'master' of https://github.com/sunpy/sunpy into

a488175

power_spectra

PEP-8 errors

ce8a3d0

dstansby reviewed Mar 7, 2018

View reviewed changes

dpshelio reviewed Mar 7, 2018

View reviewed changes

yashkgp added 4 commits March 11, 2018 00:07

Merge branch 'master' of https://github.com/sunpy/sunpy into power_sp…

9d26e08

…ectra

Adds Power spectra example

4453d37

Removes power spectra from sunpy core

02b9da8

Minor typos

46605b3

yashkgp changed the title ~~Implements Power Spectrum for a TimeSeries Object~~ Adds Power Spectrum example for a timeseries Mar 10, 2018

nabobalis reviewed Mar 12, 2018

View reviewed changes

dpshelio reviewed Mar 12, 2018

View reviewed changes

yashkgp added 3 commits March 14, 2018 16:16

Added units for graph

b63a8e6

Merge branch 'master' of https://github.com/sunpy/sunpy into power_sp…

c9480f9

…ectra

Removed extra lines

ae100e5

nabobalis reviewed Mar 15, 2018

View reviewed changes

dpshelio requested changes Mar 15, 2018

View reviewed changes

nabobalis added this to the 0.9 milestone Mar 16, 2018

yashkgp added 4 commits March 17, 2018 23:36

Added small description about data

f34ba74

Merge branch 'master' of https://github.com/sunpy/sunpy into power_sp…

65242cb

…ectra

Changed the input data

fc58fc0

minor changes

17db1aa

nabobalis reviewed Mar 19, 2018

View reviewed changes

yashkgp added 2 commits March 19, 2018 17:23

Documentation

15dce34

Documentation

0fa1376

nabobalis reviewed Mar 19, 2018

View reviewed changes

minor changes

2d4c012

dpshelio requested changes Mar 19, 2018

View reviewed changes

yashkgp added 2 commits March 23, 2018 03:07

Removed Pandas

ea781d2

Merge branch 'master' of https://github.com/sunpy/sunpy into power_sp…

5a6dcf5

…ectra

nabobalis reviewed Apr 2, 2018

View reviewed changes

nabobalis requested changes Apr 2, 2018

View reviewed changes

yashkgp added 2 commits April 11, 2018 06:56

Added links

3045553

Merge branch 'master' of https://github.com/sunpy/sunpy into power_sp…

e40ee2c

…ectra

nabobalis approved these changes Apr 11, 2018

View reviewed changes

nabobalis reviewed Apr 11, 2018

View reviewed changes

dpshelio approved these changes Apr 11, 2018

View reviewed changes

nabobalis requested changes Apr 12, 2018

View reviewed changes

yashkgp added 2 commits April 12, 2018 21:37

improves Documentation

9e673fd

Merge branch 'master' of https://github.com/sunpy/sunpy into power_sp…

5f5abec

…ectra

nabobalis approved these changes Apr 13, 2018

View reviewed changes

nabobalis merged commit ae70bda into sunpy:master Apr 13, 2018

yashkgp deleted the power_spectra branch April 16, 2018 10:38


		###############################################################################
		# Let's first load a TimeSeries from sample data.

Adds Power Spectrum example for a timeseries #2496

Adds Power Spectrum example for a timeseries #2496

Conversation

yashkgp commented Mar 2, 2018 • edited

pep8speaks commented Mar 2, 2018 • edited

Comment last updated on April 12, 2018 at 11:42 Hours UTC

dstansby left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nabobalis commented Mar 7, 2018

yashkgp commented Mar 7, 2018

yashkgp commented Mar 12, 2018 • edited

nabobalis left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

yashkgp Mar 12, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

yashkgp commented Mar 12, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

yashkgp commented Mar 14, 2018 • edited

nabobalis commented Mar 15, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

yashkgp Mar 17, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

yashkgp commented Mar 19, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

yashkgp Mar 20, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nabobalis commented Mar 23, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nabobalis left a comment

Choose a reason for hiding this comment

yashkgp commented Apr 10, 2018

Choose a reason for hiding this comment

yashkgp commented Mar 2, 2018 •

edited

pep8speaks commented Mar 2, 2018 •

edited

yashkgp commented Mar 12, 2018 •

edited

yashkgp Mar 12, 2018 •

edited

yashkgp commented Mar 14, 2018 •

edited

yashkgp Mar 17, 2018 •

edited

yashkgp Mar 20, 2018 •

edited