Sunsplots

sunspots_periodogram.py

# -*- coding: utf-8 -*-
"""
Created on Sat Feb 20 2016
"""

from __future__ import division, print_function
import Quandl as q
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from numpy.fft import rfft, rfftfreq
from scipy.signal import blackmanharris

file_name = 'sunspots.pkl'

try:
    sunspots = pd.read_pickle(file_name)
except IOError:
    sunspots = q.get("SIDC/SUNSPOTS_M")
    sunspots.to_pickle(file_name)

# Checked with plot(diff(sunspots.index)) that there are no gaps or jumps

column = 'Monthly Mean Total Sunspot Number'

ss_time = sunspots.index.year + sunspots.index.month/12
ss_dt = np.diff(ss_time).mean()
ss = np.array(sunspots[column])

t = ss_time
sig = ss - ss.mean()

N = len(sig)

fs = 1/(t[1] - t[0])

w = blackmanharris(N)

N_interp = N * 51

ampl = 1/N * np.absolute(rfft(w * sig, N_interp))

freqs = 1/rfftfreq(N_interp, 1/fs)

plt.figure('Spectrum')
plt.plot(freqs, 20*np.log10(ampl))
plt.xscale('log')
plt.xlim(1, 200)
plt.ylim(-40, 30)

peak = freqs[np.argmax(ampl)]
print(peak, 'years')
for cycle in (peak, 2*peak, 4*peak):
    plt.axvline(cycle, color='black', ls='--', alpha=0.5)

plt.title('Sunspot periodogram ({:.1f} year cycle)'.format(peak))
plt.xlabel('Period [yr]')
plt.ylabel('Amplitude [dB]')
plt.margins(0.01, 0.1)
plt.grid(True, color='0.7', linestyle='-', which='major')
plt.grid(True, color='0.9', linestyle='-', which='minor')
plt.show()

sunspots_scalogram.py

# -*- coding: utf-8 -*-
"""
Created on Sat Feb 20 2016

Copy of plot in
http://wattsupwiththat.com/2008/09/22/new-cycle-24-sunspot/
"""

from __future__ import division
import Quandl as q
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from wavelets import WaveletAnalysis, Ricker, Morlet, DOG  # arren
from bipolar import bipolar

file_name = 'sunspots.pkl'

try:
    sunspots = pd.read_pickle(file_name)
except IOError:
    sunspots = q.get("SIDC/SUNSPOTS_M")
    sunspots.to_pickle(file_name)

# Checked with plot(diff(sunspots.index)) that there are no gaps or jumps

column = 'Monthly Mean Total Sunspot Number'

# TODO: WHY DOESN'T SHAREX WORK
fig, (ax_t, ax_w) = plt.subplots(2, 1) #, sharex=True)
ax_t.plot(sunspots.index, sunspots[column])
# 12 = start at 1750 instead of 1749 TODO automate this
ax_t.set_xticks(sunspots.index[12::12*50])
ax_t.set_title(column)
ax_t.set_ylabel('sunspot number')


# TODO: Would be nice if it could handle pandas series type things directly
ss_time = sunspots.index.year + sunspots.index.month/12
ss_dt = np.diff(ss_time).mean()
ss = np.array(sunspots[column])

# shouldn't dt be calculated from time internally?
wa = WaveletAnalysis(ss, time=ss_time, dt=ss_dt, dj = .125/4, wavelet=Morlet(6))

t, s = wa.time, wa.scales

# plot the wavelet power
T, S = np.meshgrid(t, s)

## which is better?
ax_w.contourf(T, S, wa.wavelet_power, 256)
#ax_w.pcolormesh(T, S, wa.wavelet_power)

# cmap = 'RdBu'
cmap = bipolar(neutral=0)

#data = wa.wavelet_transform.real
#ax_w.pcolormesh(T, S, data, cmap=cmap, vmax=np.amax(abs(data)),
#                vmin=-np.amax(abs(data)))

ax_w.set_xlabel('year')
ax_w.set_yscale('log')
ax_w.grid(True)
ax_w.set_ylabel('fourier period (years)')
ax_w.set_ylim(wa.scales.min(), wa.scales.max())

# shade the region between the edge and coi
C, S = wa.coi
ax_w.fill_between(x=C, y1=S, y2=s.max(), color='gray', alpha=0.3)
ax_w.set_xlim(t.min(), t.max())
for cycle in (11.1, 22.2, 44.4):
    plt.axhline(cycle, color='white', alpha=0.3, ls='--')

plt.tight_layout()

Output:

Real part with hotcold bipolar colormap:

Power with viridis colormap:

hey! Is there a script where in I can find rms of large amount of data and plot it on matplotlib?