$\textbf{From CRU homepage}$
ASCII data: The 360-lat x 720-long grid is presented exactly as that, with 720 columns, and 360 rows per timestep. The first row in each grid is the southernmost (centred on 89.75S). The first column is the westernmost (centred on 179.75W). There are scaling factors in use in the data files (see Table below). One gets the whole global grid for the first time step, then the whole grid for the second, and so on. So the first 360 rows show the data for Jan 1901, next 360 rows the data for Feb 1901, next 360 rows for March 1901 and so on.
Factor: x10
In [68]:
# Load data
first_year = 1
last_year = 10
#degWest = [120, 90]
#columns = degWest[0]+arange(2*(degWest[0]-degWest[1]))
columns = 340 + arange(90)
c = genfromtxt('rawData/cru_ts3.21.1901.2012.pre.dat.gz', dtype=int16, usecols=columns)
#c1 = genfromtxt('rawData/cru_ts3.21.1911.1920.pre.dat.gz', dtype=int16, skip_header=360*12*(first_year-1),skip_footer=360*12*(10-last_year), usecols=columns)
#c2 = genfromtxt('rawData/cru_ts3.21.1921.1930.pre.dat.gz', dtype=int16, skip_header=360*12*(first_year-1),skip_footer=360*12*(10-last_year), usecols=columns)
#c3 = genfromtxt('rawData/cru_ts3.21.1931.1940.pre.dat.gz', dtype=int16, skip_header=360*12*(first_year-1),skip_footer=360*12*(10-last_year), usecols=columns)
#c4 = genfromtxt('rawData/cru_ts3.21.1961.1970.pre.dat.gz', dtype=int16, skip_header=360*12*(first_year-1),skip_footer=360*12*(10-last_year), usecols=columns)
print c.shape
#c.reshape((last_year-first_year+1,360,len(columns)))
#c = genfromtxt('data/cru_ts3.21.1921.1930.pre.dat.gz', dtype=int16, skip_header=360*12*(year-1),skip_footer=360*12*(10-year))
In [82]:
# Full data
T = 112
y = zeros((T,360,len(columns)))
degNorth = [6,30]
rows = [2*(90-degNorth[1]), 2*(90-degNorth[0])]
for t in range(T):
for i in range(12):
temp = c[360*i+12*360*t:12*360*t+360*(i+1),:]
temp = temp[::-1,:]
if np.sum(temp[rows[0]:rows[1],:-2] < 0) > 0:
print 't: ', t, ' i: ',i
y[t,:,:] += c[360*i+12*360*t:12*360*t+360*(i+1),:]
y[t,:,:] /= 12.
y[t,:,:] = y[t,::-1,:]
imshow(y[0,rows[0]:rows[1],:-2])
colorbar()
Out[82]:
In [87]:
Y = zeros((T,rows[1]-rows[0]+1,len(columns)-2))
Y = y[:,rows[0]:rows[1],:-2]/10.
Yhalf = zeros((T,(rows[1]-rows[0]+1)/2,(len(columns)-2)/2))
for i in range((rows[1]-rows[0]+1)/2):
for j in range((len(columns)-2)/2):
Yhalf[:,i,j] = 0.25*(Y[:,2*i,2*j] + Y[:,2*i,2*j+1] + Y[:,2*i+1,2*j] + Y[:,2*i+1,2*j+1])
In [3]:
T = last_year - first_year + 1
y = zeros((3*T,360,len(columns)))
degNorth = [35, 55]
rows = [2*(90-degNorth[1]), 2*(90-degNorth[0])]
for t in range(T):
for i in range(12):
temp = c1[360*i+12*360*t:12*360*t+360*(i+1),:]
temp = temp[::-1,:]
if np.sum(temp[rows[0]:rows[1],:] < 0) > 0:
print 't: ', t, ' i: ',i
y[t,:,:] += c1[360*i+12*360*t:12*360*t+360*(i+1),:]
y[t,:,:] /= 12.
y[t,:,:] = y[t,::-1,:]
for t in range(T):
for i in range(12):
temp = c2[360*i+12*360*t:12*360*t+360*(i+1),:]
temp = temp[::-1,:]
if np.sum(temp[rows[0]:rows[1],:] < 0) > 0:
print 't: ', t, ' i: ',i
y[t+T,:,:] += c2[360*i+12*360*t:12*360*t+360*(i+1),:]
y[t+T,:,:] /= 12.
y[t+T,:,:] = y[t+T,::-1,:]
for t in range(T):
for i in range(12):
temp = c3[360*i+12*360*t:12*360*t+360*(i+1),:]
temp = temp[::-1,:]
if np.sum(temp[rows[0]:rows[1],:] < 0) > 0:
print 't: ', t, ' i: ',i
y[t+2*T,:,:] += c2[360*i+12*360*t:12*360*t+360*(i+1),:]
y[t+2*T,:,:] /= 12.
y[t+2*T,:,:] = y[t+2*T,::-1,:]
In [4]:
Y = zeros((3*T,rows[1]-rows[0]+1,len(columns)))
Y = y[:,rows[0]:rows[1],:]/10.
Yhalf = zeros((3*T,(rows[1]-rows[0]+1)/2,len(columns)/2))
for i in range((rows[1]-rows[0]+1)/2):
for j in range(len(columns)/2):
Yhalf[:,i,j] = 0.25*(Y[:,2*i,2*j] + Y[:,2*i,2*j+1] + Y[:,2*i+1,2*j] + Y[:,2*i+1,2*j+1])
In [90]:
print np.sum(Y < 0)
#for t in range(3*T):
# filename = 'processedData/dustBowlYt'+str(1910+t+1) +'_N35-55_W90-120_full.csv'
# np.savetxt(filename, Y[t,:,:], delimiter=',')
for t in range(T):
filename = 'processedData/sahelYt'+str(1900+t+1) +'_N35-55_W90-120_downsampled.csv'
np.savetxt(filename, Yhalf[t,:,:], delimiter=',')
In [116]:
mu = mean(Yhalf,axis=0)
sigma2 = mean((mu-Yhalf)**2,axis=0)
print Yhalf[0,:,:].shape
mu_ab = zeros(Yhalf[0,:,:].shape)
mu_norm = zeros(Yhalf[0,:,:].shape)
for i in range(Yhalf.shape[1]):
for j in range(Yhalf.shape[2]):
if np.sum(Yhalf[:,i,j]) == 0.:
mu_ab[i,j] = 0.
mu_norm[i,j] = 0.
else:
perc = percentile(Yhalf[:,i,j],q=15)
mu_ab[i,j] = mean(Yhalf[Yhalf[:,i,j]<=perc,i,j])
if isnan(mu_ab[i,j]):
print perc
print Yhalf[:,i,j]
if np.sum(Yhalf[:,i,j]>perc) == 0.:
mu_norm[i,j] = mu_ab[i,j]
else:
mu_norm[i,j] = mean(Yhalf[Yhalf[:,i,j]>perc,i,j])
if isnan(mu_norm[i,j]):
print Yhalf[:,i,j]>perc
print Yhalf[:,i,j]
In [130]:
filename = 'parameters/sahelSigma2_N35-55_W90-120_downsampled.csv'
np.savetxt(filename, sigma2, delimiter=',')
filename = 'parameters/sahelMuAb_N35-55_W90-120_downsampled.csv'
np.savetxt(filename, mu_ab, delimiter=',')
filename = 'parameters/sahelMuNorm_N35-55_W90-120_downsampled.csv'
np.savetxt(filename, mu_norm, delimiter=',')
In [122]:
imshow(mu_norm)
colorbar()
np.max(sigma2)
Out[122]:
In [25]:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
# setup Lambert Conformal basemap.
m = Basemap(width=12000000,height=9000000,projection='lcc',
resolution='c',lat_1=45.,lat_2=55,lat_0=50,lon_0=-107.)
# draw coastlines.
m.drawcoastlines()
# draw a boundary around the map, fill the background.
# this background will end up being the ocean color, since
# the continents will be drawn on top.
m.drawmapboundary(fill_color='aqua')
# fill continents, set lake color same as ocean color.
m.fillcontinents(color='coral',lake_color='aqua')
plt.show()
In [22]:
test = zeros(10,dtype=bool)
2.5*(0==test.astype('float'))
Out[22]:
In [89]:
imshow(Yhalf[0,:,:])
Out[89]:
In [56]:
year = 1910
imshow(-0.5*(Yhalf[year-1901,:,:]-mu_norm)**2/sigma2)
colorbar()
Out[56]:
In [59]:
year = 1926
nrNorm = (-0.5*(Yhalf[year-1901,:,:]-mu_norm)**2/sigma2) > (-0.5*(Yhalf[year-1901,:,:]-mu_ab)**2/sigma2)
print np.sum(nrNorm)
imshow(nrNorm)
colorbar()
Out[59]:
In [149]:
vecNorm = zeros(T)
for t in range(T):
nrNorm = (-0.5*(Yhalf[t,:,:]-mu_norm)**2/sigma2) > (-0.5*(Yhalf[t,:,:]-mu_ab)**2/sigma2)
vecNorm[t] = np.sum(nrNorm)
In [150]:
startYear = 1960
endYear = 2000
ind = arange(startYear-1901,endYear-1900)
plot(arange(startYear,endYear+1),vecNorm[ind])
Out[150]:
In [124]:
sum(isnan(mu_ab))
Out[124]:
In [127]:
mS2 = mean(sigma2)
In [131]:
24*44
Out[131]:
In [129]:
sum(sigma2==0.)
Out[129]:
In [152]:
test = exp(3*arange(2))
test /= sum(test)
print test
In [ ]: