# Install needed packages
# pip install numpy
#     Etc etc

import os
import sys
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import scipy.io as sio
from scipy import signal

# To Download files try
# wget http://iodlabs.ucsd.edu/sccrosby/RBR_therm_01.mat
# Therm is a dictionary, and has the following keys
#     Type Therm.keys()
#     time, name, temp, hob (height above bottom)

Therm1 = sio.loadmat('RBR_therm_01.mat')
Therm2 = sio.loadmat('RBR_therm_02.mat')
print Therm1.keys()

# Setup for plot time series
label_1= Therm1['hob'][0]
y1 = Therm1['temp_LR'][0]
x1 = Therm1['year_days_LR'][0]

label_2= Therm2['hob'][0]
y2 = Therm2['temp_LR'][0]
x2 = Therm2['year_days_LR'][0]

plt.figure(num=1,figsize=(12,6),dpi=150)

# Plot temp and time
plt.plot(x1, y1, '-', label=label_1)
plt.plot(x2, y2, '-', label=label_2)
plt.legend()
plt.xlabel('Year Day')
plt.ylabel('Temperature [C]')
plt.show()
plt.savefig('TimeSeries.png')
plt.close()

# Now estimate power spectra
# Time step for LR is 5min
Ts = 5./(24*60) # 5 min sampling period in days
fs = 1./Ts
f1, Txx1 = signal.welch(y1, fs, nperseg=256)
f2, Txx2 = signal.welch(y2, fs, nperseg=256)

# Now plot power spectra
plt.figure(num=1,figsize=(12,6),dpi=150)
plt.loglog(f1, Txx1, '-', label=label_1)
plt.loglog(f2, Txx2, '-', label=label_2)
plt.legend()
plt.xlabel('Cycle-per-day (cpd)')
plt.ylabel('Power Spectra Density C^2/Hz')
plt.show()
plt.savefig('PowerSpec.png')
plt.close()


# Try plotting histograms in violin plots
plt.figure(num=1,figsize=(12,6),dpi=150)
sns.violinplot(scale="width",data=[y1,y2])
plt.ylabel('Temperature [C]')
plt.show()
plt.savefig('Histograms.png')
plt.close()