``````

In [1]:

%matplotlib inline
from scipy import interpolate
from scipy import special
from scipy.signal import butter, lfilter, filtfilt
import matplotlib.pyplot as plt
import numpy as np
from numpy import genfromtxt
from nitime import algorithms as alg
from nitime import utils
from scipy.stats import t
import xray
import pandas as pd
from rpy2.robjects import FloatVector
from rpy2.robjects.vectors import StrVector
import rpy2.robjects as robjects
from rpy2.robjects.packages import importr
r = robjects.r

``````

Define functions for filtering, moving averages, and normalizing data

``````

In [2]:

def butter_lowpass(cutoff, fs, order=3):
nyq = 0.5 * fs
normal_cutoff = cutoff / nyq
b, a = butter(order, normal_cutoff, btype='low', analog=False)
return b, a

def filter(x, cutoff, axis, fs=1.0, order=3):
b, a = butter_lowpass(cutoff, fs, order=order)
y = filtfilt(b, a, x, axis=axis)
return y

def movingaverage(interval, window_size):
window = np.ones(int(window_size))/float(window_size)
return np.convolve(interval, window, 'valid')

def owncorr(x,y,n):
x_ano=np.ma.anomalies(x)
x_sd=np.sum(x_ano**2,axis=0)
y_ano=np.ma.anomalies(y)
y_sd=np.sum(y_ano**2,axis=0)
nomi = np.dot(x_ano,y_ano)
corr = nomi/np.sqrt(np.dot(x_sd[None],y_sd[None]))

x_coef, x_sigma = alg.AR_est_YW (x, 1)
y_coef, y_sigma = alg.AR_est_YW (y, 1)
neff = n*(1-x_coef*y_coef)/(1+x_coef*y_coef)

if neff <3:
neff = 3

coef = []
coef.append(x_coef)
coef.append(y_coef)
tval = corr/np.sqrt(1-corr**2)*np.sqrt(neff-2)
pval = t.sf(abs(tval),neff-2)*2

return corr,pval,coef

def gaussianize(X):
n = X.shape[0]
#p = X.shape[1]

Xn = np.empty((n,))
Xn[:] = np.NAN
nz = np.logical_not(np.isnan(X))

index = np.argsort(X[nz])
rank = np.argsort(index)

CDF = 1.*(rank+1)/(1.*n) -1./(2*n)
Xn[nz] = np.sqrt(2)*special.erfinv(2*CDF -1)
return Xn

``````

Read bandwidth and rain/temperature data and normalize them

``````

In [4]:

data = genfromtxt('data/scotland.csv', delimiter=',')
bandw = data[0:115,4] # band width (1879-1993), will be correlated with T/P
bandwl = data[3:129,4] # band width (1865-1990), will be correlation with winter NAO
bandwn = gaussianize(bandw) #normalized band width
bandwln = gaussianize(bandwl) #normalized band width

rain = genfromtxt('data/Assynt_P.txt') #precipitaiton
temp = genfromtxt('data/Assynt_T.txt') #temperature
wnao = genfromtxt('data/wnao.txt') #winter NAO
wnao = wnao[::-1]

rainn = gaussianize(rain)
tempn = gaussianize(temp)

#calculate the ratio of temperature over precipitation
ratio = temp/rain
ration = gaussianize(ratio)

``````

Smoothing data (11-year running average)

``````

In [5]:

bandw_fil = movingaverage(bandw, 11)
bandwn_fil = movingaverage(bandwn, 11)
bandwl_fil = movingaverage(bandwl, 11)
rain_fil = movingaverage(rain, 11)
rainn_fil = movingaverage(rainn, 11)
ratio_fil = movingaverage(ratio, 11)
wnao_fil = movingaverage(wnao, 11)

``````

Calculate correlation and p-values with considering autocorrelation, and the autocorrelations (coef)

``````

In [6]:

corr_ratio,pval_ratio,coef = owncorr(bandw_fil,ratio_fil,115) #correlation between smoothed bandwidth and ratio
corr_nao,pval_nao,coef_nao = owncorr(bandwl_fil,wnao_fil,126) #correlation between smoothed bandwidth and winter NAO
corr_n,pval_n,coef_n = owncorr(bandwn,ration,115) #correlation between normalized bandwidth and ratio
corr_naon,pval_naon,coef_naon = owncorr(bandwln,wnao,126) #correlation between normalized bandwidtha and winter NAO

``````

Check the correlation results

``````

In [7]:

print(corr_ratio)
print(pval_ratio)
print(coef)
print(corr_nao)
print(pval_nao)
print(coef_nao)

print(corr_n)
print(pval_n)
print(coef_n)
print(corr_naon)
print(pval_naon)
print(coef_naon)

``````
``````

0.769730956482
0.440780720831
[array([ 0.99411616]), array([ 0.99106323])]
-0.794624754752
[ 0.37152434]
[array([ 0.99594871]), array([ 0.95510202])]
0.306696255799
[ 0.00285928]
[array([ 0.77346338]), array([ 0.13991603])]
-0.208749844418
[ 0.0420236]
[array([ 0.77795155]), array([ 0.17840258])]

``````
``````

In [ ]:

``````