import numpy as np
import matplotlib.pyplot as pyp
import pandas as pnd
from IPython.display import Image
import scipy.sparse.linalg as sp
import h5py
from IPython.display import Math

a = np.sin(np.arange(1,100,0.1)/np.pi*3.5)+np.sin(np.arange(1,100,0.1)/2.5)*np.random.rand(990)*1.7
b = np.sin(np.arange(1,100,0.1)/np.pi*3.5)+np.sin(np.arange(1,100,0.1)/7)*np.random.rand(990)*1.5

plot(a,'r',label='a')
plot(b,'k',label='b')
legend()
title('Two timeseries')
xlabel('Time')

a = pnd.Series(a) #Here I turn the timeseries into a "Series" to use the pandas library
b = pnd.Series(b)
correlation = a.corr(b)

print "The correlation between a and b is: ", correlation

k = pyp.xcorr(a, b, maxlags = 200)
title('xcorr of timeseries a and b')
xlabel('Lag')
ylabel('Normalised correlation strength')

ind = np.where(np.abs(k[:][0][np.where(np.abs(k[:][1]) < 0.03)]) == np.abs(k[:][0][np.where(np.abs(k[:][1]) < 0.03)]).min())
df = k[:][0][np.where(np.abs(k[:][1]) < 0.03)][ind]

print "The degrees of freedom are: ", df

Image(filename='/noc/users/mjsp106/SignifOfCorrelations.png')

f = h5py.File('/noc/msm/scratch/valor2/maike/beringDump/concatenateData/N1/ORCA1_S_T_SSH_SST_SSS.nc')
SSH=f['sossheig']
ssh = squeeze(SSH[:])
SST=f['sosstsst']
sst = squeeze(SST[:])
#Some land data, just por plotting. Oceanographers tend to politely ignore the land.
land = h5py.File("/noc/msm/scratch/valor2/maike/beringDump/land_N1_global.h5", 'r')
land1=land['land']

cclim = np.array([-2,2])
imshow(flipud(ssh[0,...]), cmap=cm.coolwarm)
clim(cclim[:])
colorbar()
imshow(flipud(land1[:]))
title('Sea Surface Height (m) January 1978')

Math(r'\mathrm{M} = \mathrm{U} \Sigma \mathrm{V}^{*}')

u_ssh,s_ssh,v_ssh=sp.svds(np.transpose(ssh.reshape((360, 292*362))), k=300)

eof1issh = u_ssh[:, -1]
eof1ssh=eof1issh.reshape(292, 362)
eof2issh = u_ssh[:, -2]
eof2ssh=eof2issh.reshape(292, 362)
eof3issh = u_ssh[:, -3]
eof3ssh=eof3issh.reshape(292, 362)

cclim = np.array([-0.01,0.01])
imshow(flipud(eof1ssh), cmap=cm.coolwarm)
clim(cclim[:])
colorbar()
imshow(flipud(land1[:]))
title('Sea Surface Height EOF1')
figure()
imshow(flipud(eof2ssh), cmap=cm.coolwarm)
clim(cclim[:])
colorbar()
imshow(flipud(land1[:]))
title('Sea Surface Height EOF2')
figure()
imshow(flipud(eof3ssh), cmap=cm.coolwarm)
clim(cclim[:])
colorbar()
imshow(flipud(land1[:]))
title('Sea Surface Height EOF3')

s2_ssh = s_ssh[:-1]**2
percent_ssh = (100*s2_ssh)/sum(s2_ssh)
sumvar_ssh = sum(percent_ssh[:])

plot(flipud(percent_ssh[259:]))
plot(flipud(percent_ssh[259:]), 'k*')
title('The first 40 EOFs SSH')
xlabel('EOF number')
ylabel('Percentage of variance (%)')
print sumvar_ssh

cclim = np.array([-1,31])
imshow(flipud(sst[0,...]), cmap=cm.coolwarm)
clim(cclim[:])
colorbar()
imshow(flipud(land1[:]))
title('Sea Surface Temperature (m) January 1978')

u_sst,s_sst,v_sst=sp.svds(np.transpose(sst.reshape((360, 292*362))), k=300)

eof1isst = u_sst[:, -1]
eof1sst=eof1isst.reshape(292, 362)
eof2isst = u_sst[:, -2]
eof2sst=eof2isst.reshape(292, 362)
eof3isst = u_sst[:, -3]
eof3sst=eof3isst.reshape(292, 362)

cclim = np.array([-0.01,0.01])
imshow(flipud(eof1sst), cmap=cm.coolwarm)
clim(cclim[:])
colorbar()
imshow(flipud(land1[:]))
title('Sea Surface Temperature EOF1')
figure()
imshow(flipud(eof2sst), cmap=cm.coolwarm)
clim(cclim[:])
colorbar()
imshow(flipud(land1[:]))
title('Sea Surface Temperature EOF2')
figure()
imshow(flipud(eof3sst), cmap=cm.coolwarm)
clim(cclim[:])
colorbar()
imshow(flipud(land1[:]))
title('Sea Surface Temperature EOF3')

s2_sst = s_sst[:-1]**2
percent_sst = (100*s2_sst)/sum(s2_sst)
sumvar_sst = sum(percent_sst[:])

plot(flipud(percent_sst[259:]))
plot(flipud(percent_sst[259:]), 'k*')
title('The first 40 EOFs SST')
xlabel('EOF number')
ylabel('Percentage of variance (%)')
print sumvar_sst

plot(flipud(percent_sst[259:]),label='SST')
plot(flipud(percent_sst[259:]), 'k*')
plot(flipud(percent_ssh[259:]),label='SSH')
plot(flipud(percent_ssh[259:]), 'k*')
legend()
title('The first 40 EOFs SST and SSH')
xlabel('EOF number')
ylabel('Percentage of variance (%)')



Image(filename='/noc/users/mjsp106/svd_memory.png')


u,s,v=sp.svds(np.transpose(ssh.reshape((360, 292*362))), k=359)
reconstr = np.dot(np.dot(u, np.diag(s)), v)

reconstrTest=reconstr[:, 0].reshape(292, 362)

cclimssh = np.array([-1,1])
imshow(flipud(reconstrTest), cmap=cm.coolwarm)
clim(cclimssh[:])
colorbar()
imshow(flipud(land1[:]))
title('SSH reconstructed from 359 EOFs')
figure()
imshow(flipud(ssh[0,...]), cmap=cm.coolwarm)
clim(cclimssh[:])
colorbar()
title('SSH original')
imshow(flipud(land1[:]))
figure()
imshow(flipud(reconstrTest-ssh[0,...]), cmap=cm.coolwarm)
clim(cclim[:])
colorbar()
imshow(flipud(land1[:]))
title('SSH reconstructed - SSH original')