-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathclimstats.py
160 lines (140 loc) · 5.11 KB
/
climstats.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
import numpy as np
import scipy as sp
from scipy import stats
import inout as io
import rstats as rs
# Correlation
def cormap(txt_ts,ncarray,lag=0):
"""input: txt_file... a text file string
ncfilein... an ncfile - string
"""
col=0
tdim = ncarray.shape[0]
if ncarray.ndim==4:
map_ts=ncarray[:,0,:,:]
elif ncarray.ndim==3:
map_ts=ncarray
t,lat,lon = map_ts.shape
cor_map = np.zeros((lat,lon))
for i in xrange(lat):
for j in xrange(lon):
#slope, intercept, r_val, p_val, stderr = stats.linregress(map_ts[:,i,j],txt_ts[:,col])
r_mat = np.corrcoef(map_ts[lag:tdim,i,j],txt_ts[0:tdim-lag,col])
cor_map[i,j] = r_mat[1,0]
out_cor_map = np.nan_to_num(cor_map)
return out_cor_map
# Regression
"""Maps the regression values (or other values) between a timeseries, e.g. NINO3 and an array, e.g. Hadisst
txt_ts is a text file of the timeseries, and col is the column to use (if there are more than one, otherwise leave blank)
statvar is the statistical variable: slope=0, intercept=1, r_val=2, p_val=3, stderr=4
Lag is in MONTHS
"""
def regmap(txt_ts,data,lag=0,statvar=0):
"""input: txt_ts = 1D timeseries array
data = 3D or 4D data array
statvar=0 for regression, 2 for correlation
"""
#txt_ts = io.readtext(txt_file)
#ncarray = io.readnc(ncfilein,ncvar)[0]
### Note: program changed, originally the input was a string of a text or nc file, now it's an array
### Also removed 'col=0' from input, so txt array has to be 1D
tdim=data.shape[0]
if data.ndim==4:
map_ts=data[:,0,:,:]
elif data.ndim==3:
map_ts=data
t,lat,lon = map_ts.shape
reg_map = np.zeros((lat,lon))
for i in xrange(lat):
for j in xrange(lon):
#slope, intercept, r_val, p_val, stderr = stats.linregress(map_ts[:,i,j],txt_ts[:,col])
#stattup = stats.linregress(map_ts[lag:1706,i,j],txt_ts[0:1706-lag,col])
print('i = '+str(i))
print('j = '+str(j))
stattup = stats.linregress(txt_ts[0:tdim-lag],map_ts[lag:tdim,i,j])
#stattup = stats.linregress(txt_ts[11+lag:1695,col],map_ts[11:1695-lag,i,j])
print(stattup[statvar])
reg_map[i,j] = stattup[statvar] #.copy()
reg_II = np.nan_to_num(reg_map)
out_reg = np.ma.masked_outside(reg_II,-1e9,1e9)
return out_reg
def regmap_box(temp,var,stat=0,lag=0):
"""
Input a temperature array [time,lat,lon] and array of variable in question
e.g. pressure [time,lat,lon].
Outputs an array [lat,lon] of regression values, where each grid box is the temperature regessed with the variale
stat=0: use 0 for regression, 2 for correlatino
"""
t,lat,lon = temp.shape
reg_map = np.zeros((lat,lon))
for i in xrange(lat):
for j in xrange(lon):
#slope, intercept, r_val, p_val, stderr = stats.linregress(map_ts[:,i,j],txt_ts[:,col])
#stattup = stats.linregress(map_ts[lag:1706,i,j],txt_ts[0:1706-lag,col])
stattup = stats.linregress(temp[lag:t,i,j],var[0:t-lag,i,j])
#stattup = stats.linregress(txt_ts[11+lag:1695,col],map_ts[11:1695-lag,i,j])
reg_map[i,j] = stattup[stat]
reg_II = np.nan_to_num(reg_map)
out_reg = np.ma.masked_outside(reg_II,-100000,100000)
return out_reg
def crosscormap(txt_ts,data,maxlag=24):
"""input: txt_ts = 1D timeseries array (monthly data is best)
data = 3D or 4D data array
This program will caculate the maximum cross-correlation, and lag time of that max, between a given timeseries and each gridbox
"""
tdim=data.shape[0]
if data.ndim==4:
map_ts=data[:,0,:,:]
elif data.ndim==3:
map_ts=data
t,lat,lon = map_ts.shape
maxcor_map = np.zeros((lat,lon))
maxcorlag_map = np.zeros((lat,lon))
for i in xrange(lat):
for j in xrange(lon):
xcv, lags=rs.xcorr(txt_ts[0:tdim],ts2=map_ts[0:tdim,i,j],maxlag=maxlag)
maxcor=abs(xcv).max()
maxcorlag=np.array([abs(np.argmax(xcv)-maxlag),abs(np.argmin(xcv)-maxlag)]).min()
maxcor_map[i,j] = np.nan_to_num(maxcor.copy())
maxcorlag_map[i,j] = np.nan_to_num(maxcorlag.copy())
#out_reg = np.ma.masked_outside(reg_II,-1e-9,1e9)
return maxcor_map, maxcorlag_map
# Tl/To
def tltomap(datain,txt_ts,lag=0):
"""input: txt_file... a text file string
ncfilein... an ncfile - string
"""
col=0
tdim = datain.shape[0]
if datain.ndim==4:
map_ts=datain[:,0,:,:]
elif datain.ndim==3:
map_ts=datain
t,lat,lon = map_ts.shape
tlto_map = np.zeros((lat,lon))
c_map = np.zeros((lat,lon))
rs_map = np.zeros((lat,lon))
for i in xrange(lat):
for j in xrange(lon):
# Calculate the land/ocean temp contrast, phi
#def phi(ts1,ts2):
#Input: Two 1D timeseries, e.g. land/ocean.
#Output: Phi, c (correlation coefficient), std(ts1)/std(ts2)
c=np.corrcoef(map_ts[lag:tdim,i,j],txt_ts[0:tdim-lag])[0,1]
rs=np.std(map_ts[lag:tdim,i,j])/np.std(txt_ts[0:tdim-lag])
phi = c*rs
# return phi,c,rs
#slope, intercept, r_val, p_val, stderr = stats.linregress(map_ts[:,i,j],txt_ts[:,col])
tlto_map[i,j] = phi
c_map[i,j] = c
rs_map[i,j] = rs
#out_cor_map = np.nan_to_num(cor_map)
return tlto_map, c_map, rs_map
# Calculate the land/ocean temp contrast, phi
def phi(ts1,ts2):
""" Input: Two 1D timeseries, e.g. land/ocean.
Output: Phi, c (correlation coefficient), std(ts1)/std(ts2) """
c=np.corrcoef(ts1,ts2)[0,1]
rs=np.std(ts1)/np.std(ts2)
phi = c*rs
return phi,c,rs