import numpy as N
from scipy import stats as S
import commands
import string
def gsl(dat):
nd = len(dat)
num = 0
for d in range(5,nd,1):
if((dat[d-4]>5.)&(dat[d-3]>5.)&(dat[d-2]>5.)&(dat[d-1]>5.)&(dat[d]>5.)):
num = num + 1
return num
#lastyear = commands.getstatusoutput('date +"%Y"')[1]
#lastyear = string.atoi(lastyear)
file = '03987.dat'
dd = N.genfromtxt(file,usecols=(0),unpack=True,dtype=(int))
mm = N.genfromtxt(file,usecols=(1),unpack=True,dtype=(int))
jj = N.genfromtxt(file,usecols=(2),unpack=True,dtype=(int))
tn = N.genfromtxt(file,usecols=(3),unpack=True,dtype=(float))
tg = N.genfromtxt(file,usecols=(4),unpack=True,dtype=(float))
tx = N.genfromtxt(file,usecols=(5),unpack=True,dtype=(float))
rr = N.genfromtxt(file,usecols=(13),unpack=True,dtype=(float))
ws = N.genfromtxt(file,usecols=(7),unpack=True,dtype=(float))
so = N.genfromtxt(file,usecols=(8),unpack=True,dtype=(float))
lastyear = jj.max()#-1
jo = N.arange(1961,lastyear+1,1);nj = len(jo)
#ni = lastyear-1991
#######################################################
tab = N.genfromtxt('potsdam2seas.tab',names=True,delimiter=",",dtype=None,comments='#')
tabs = {'tmit':'jz','nied':'jz','wind':'jz','heta':'ja','eist':'ja','rr20':'ja','tx99':'ja','rr99':'ja','cumt':'jz','taon':'ja','trni':'ja','cumr':'jz','tx01':'ja','sonn':'jz','grsl':'ja','eisr':'ja'}
np = len(tabs)
seas = ['ja','fr','so','he','wi']
nk = len(seas)
para = ['tmit','nied','wind','heta','eist','rr20','tx99','rr99','cumt','taon','trni','cumr','tx01','sonn','grsl','eisr']
print (np,nk)
f = open('potsdam2seas.csv','w')
f.write('%5s'%('YEAR'))
for tab in para:
if(tabs[tab]=='ja'):
f.write(',%5s%2s'%(tab,'ja'))
else:
for jz in seas:
f.write(',%5s%2s'%(tab,jz))
f.write('\n')
cumtja = 0
cumtfr = 0
cumtso = 0
cumthe = 0
cumtwi = 0
cumrja = 0
cumrfr = 0
cumrso = 0
cumrhe = 0
cumrwi = 0
nj = nj-1
for j in jo[:nj]:
f.write('%5i'%(j))
for tab in para:
if(tab=='tmit'): data = tg
if(tab=='nied'): data = rr
if(tab=='wind'): data = ws
if(tab=='heta'): data = tx
if(tab=='eist'): data = tx
if(tab=='rr20'): data = rr
if(tab=='tx99'): data = tx
if(tab=='rr99'): data = rr
if(tab=='cumt'): data = tg
if(tab=='taon'): data = rr
if(tab=='trni'): data = tn
if(tab=='cumr'): data = rr
if(tab=='tx01'): data = tx
if(tab=='sonn'): data = so
if(tab=='grsl'): data = tg
if(tab=='eisr'): data = rr#;data[tx>0]=0
if(tabs[tab]=='ja'):
#ind = N.where((jj==j)&(data>-999.))[0];jnd = N.where(((jj>=1961)&(jj<=1990))&(data>-999.))[0]
if(tab=='heta'): ind = N.where((jj==j)&(data>30.))[0]; tmp = len(data[ind])
if(tab=='eist'): ind = N.where((jj==j)&(data<0.))[0]; tmp = len(data[ind])
if(tab=='rr20'): ind = N.where((jj==j)&(data>20.))[0]; tmp = len(data[ind])
if(tab=='tx99'): ind = N.where((jj==j))[0]; tmp = data[ind]; tmp = N.mean(tmp[tmp>=N.percentile(data[ind],99)])
if(tab=='rr99'): ind = N.where((jj==j))[0]; tmp = data[ind]; tmp = N.sum(tmp[tmp>=N.percentile(data[ind],99)])
if(tab=='taon'): ind = N.where((jj==j)&(data<1.))[0]; tmp = len(data[ind])
if(tab=='trni'): ind = N.where((jj==j)&(data>20.))[0]; tmp = len(data[ind])
if(tab=='tx01'): ind = N.where((jj==j))[0]; tmp = N.percentile(data[ind],1)
if(tab=='grsl'): ind = N.where((jj==j))[0]; tmp = gsl(data[ind])
if(tab=='eisr'): ind = N.where((jj==j)&(tx<0.))[0]; tmp = N.sum(data[ind])
#if(tab=='cumt'): cumtja = cumtja+(N.mean(data[ind])-N.mean(data[jnd]));tmp = cumtja/30.
#if(tab=='cumt'): cumtfr = cumtfr+(N.mean(data[ind])-N.mean(data[jnd]));tmp = cumtfr/30.
#if(tab=='cumt'): cumtso = cumtso+(N.mean(data[ind])-N.mean(data[jnd]));tmp = cumtso/30.
#if(tab=='cumt'): cumthe = cumthe+(N.mean(data[ind])-N.mean(data[jnd]));tmp = cumthe/30.
#if(tab=='cumt'): cumtwi = cumtwi+(N.mean(data[ind])-N.mean(data[jnd]));tmp = cumtwi/30.
#
#if(tab=='cumr'): cumrja = cumrja+(N.sum(data[ind])-N.sum(data[jnd])/30.);tmp = cumrja/30.
#if(tab=='cumr'): cumrfr = cumrfr+(N.sum(data[ind])-N.sum(data[jnd])/30.);tmp = cumrfr/30.
#if(tab=='cumr'): cumrso = cumrso+(N.sum(data[ind])-N.sum(data[jnd])/30.);tmp = cumrso/30.
#if(tab=='cumr'): cumrhe = cumrhe+(N.sum(data[ind])-N.sum(data[jnd])/30.);tmp = cumrhe/30.
#if(tab=='cumr'): cumrwi = cumrwi+(N.sum(data[ind])-N.sum(data[jnd])/30.);tmp = cumrwi/30.
f.write(',%7.2f'%(tmp))
else:
for jz in seas:
if(jz=='ja'): ind = N.where((jj==j)&(data>-999.))[0];jnd = N.where(((jj>=1961)&(jj<=1990))&(data>-999.))[0]
if(jz=='fr'): ind = N.where((jj==j)&(mm>=3)&(mm<=5)&(data>-999.))[0];jnd = N.where(((jj>=1961)&(jj<=1990))&(mm>=3)&(mm<=5)&(data>-999.))[0]
if(jz=='so'): ind = N.where((jj==j)&(mm>=6)&(mm<=8)&(data>-999.))[0];jnd = N.where(((jj>=1961)&(jj<=1990))&(mm>=6)&(mm<=8)&(data>-999.))[0]
if(jz=='he'): ind = N.where((jj==j)&(mm>=9)&(mm<=11)&(data>-999.))[0];jnd = N.where(((jj>=1961)&(jj<=1990))&(mm>=9)&(mm<=11)&(data>-999.))[0]
if(jz=='wi'): ind = N.where((jj==j)&((mm==1)|(mm==2)|(mm==12))&(data>-999.))[0];jnd = N.where(((jj>=1961)&(jj<=1990))&((mm==1)|(mm==2)|(mm==12))&(data>-999.))[0]
if(tab=='nied'):
tmp = N.sum(data[ind])
else:
tmp = N.mean(data[ind])
if((tab=='cumt')&(jz=='ja')): cumtja = cumtja+(N.mean(data[ind])-N.mean(data[jnd]));tmp = cumtja/30.
if((tab=='cumt')&(jz=='fr')): cumtfr = cumtfr+(N.mean(data[ind])-N.mean(data[jnd]));tmp = cumtfr/30.
if((tab=='cumt')&(jz=='so')): cumtso = cumtso+(N.mean(data[ind])-N.mean(data[jnd]));tmp = cumtso/30.
if((tab=='cumt')&(jz=='he')): cumthe = cumthe+(N.mean(data[ind])-N.mean(data[jnd]));tmp = cumthe/30.
if((tab=='cumt')&(jz=='wi')): cumtwi = cumtwi+(N.mean(data[ind])-N.mean(data[jnd]));tmp = cumtwi/30.
if((tab=='cumr')&(jz=='ja')): cumrja = cumrja+(N.sum(data[ind])-N.sum(data[jnd])/30.);tmp = cumrja/30.
if((tab=='cumr')&(jz=='fr')): cumrfr = cumrfr+(N.sum(data[ind])-N.sum(data[jnd])/30.);tmp = cumrfr/30.
if((tab=='cumr')&(jz=='so')): cumrso = cumrso+(N.sum(data[ind])-N.sum(data[jnd])/30.);tmp = cumrso/30.
if((tab=='cumr')&(jz=='he')): cumrhe = cumrhe+(N.sum(data[ind])-N.sum(data[jnd])/30.);tmp = cumrhe/30.
if((tab=='cumr')&(jz=='wi')): cumrwi = cumrwi+(N.sum(data[ind])-N.sum(data[jnd])/30.);tmp = cumrwi/30.
f.write(',%7.2f'%(tmp))
f.write('\n')
f.close()
###########################
dat = N.genfromtxt('potsdam2seas.csv',names=True,delimiter=",",dtype=None)
names = dat.dtype.names
f = open('potsdam2seas.csv','a')
f.write('%5s'%('NaN'))
for name in names[1:]:
f.write(',%7s'%('NaN'))
f.write('\n')
f.write('%5i'%(jo.min()))
for name in names[1:]:
slope,intercept, r_value, p_value, std_err = S.linregress(N.arange(nj),dat[name])
if((str(name[0:4])=='cumt')or(str(name[0:4])=='cumr')):
f.write(',%7s'%('NaN'))
else:
f.write(',%7.2f'%(intercept))
f.write('\n')
f.write('%5i'%(jo.max()+15))
for name in names[1:]:
slope,intercept, r_value, p_value, std_err = S.linregress(N.arange(nj),dat[name])
if((str(name[0:4])=='cumt')or(str(name[0:4])=='cumr')):
f.write(',%7s'%('NaN'))
else:
f.write(',%7.2f'%(intercept+slope*(nj+15)))
f.write('\n')
############################################
f.write('%5s'%('NaN'))
for name in names[1:]:
f.write(',%7s'%('NaN'))
f.write('\n')
f.write('%5i'%(jo.min()+10))
for name in names[1:]:
slope,intercept, r_value, p_value, std_err = S.linregress(N.arange(nj-10),dat[name][10:])
if(str(name[0:4])!='tmit'):
f.write(',%7s'%('NaN'))
else:
f.write(',%7.2f'%(intercept))
f.write('\n')
f.write('%5i'%(jo.max()+15))
for name in names[1:]:
slope,intercept, r_value, p_value, std_err = S.linregress(N.arange(nj-10),dat[name][10:])
if(str(name[0:4])!='tmit'):
f.write(',%7s'%('NaN'))
else:
f.write(',%7.2f'%(intercept+slope*(nj+15)))
f.write('\n')
################################
f.write('%5s'%('NaN'))
for name in names[1:]:
f.write(',%7s'%('NaN'))
f.write('\n')
f.write('%5i'%(jo.min()+20))
for name in names[1:]:
slope,intercept, r_value, p_value, std_err = S.linregress(N.arange(nj-20),dat[name][20:])
if(str(name[0:4])!='tmit'):
f.write(',%7s'%('NaN'))
else:
f.write(',%7.2f'%(intercept))
f.write('\n')
f.write('%5i'%(jo.max()+15))
for name in names[1:]:
slope,intercept, r_value, p_value, std_err = S.linregress(N.arange(nj-20),dat[name][20:])
if(str(name[0:4])!='tmit'):
f.write(',%7s'%('NaN'))
else:
f.write(',%7.2f'%(intercept+slope*(nj+15)))
f.write('\n')
f.close()
######################################################
nj = nj+1
f = open('potsdam2year.csv','w')
f.write('%5s,%7s,%7s,%7s,%7s,%7s\n'%('YEAR','TMIT','NIED','HETA','EIST','WIND'))
tmit = 0
nied = 0
nd = 365
tmits = N.zeros((nj,nd),float); tmits[:,:] = N.nan
tmita = N.zeros((nj,nd),float); tmita[:,:] = N.nan
nieds = N.zeros((nj,nd),float); nieds[:,:] = N.nan
nieda = N.zeros((nj,nd),float); nieda[:,:] = N.nan
tmitja = N.zeros(nj,float)
hetaja = N.zeros(nj,float)
eistja = N.zeros(nj,float)
niedja = N.zeros(nj,float)
windja = N.zeros(nj,float)
for j in range(nj):
ind = N.where((jj==jo[j])&(tg>-900))[0]
tmit = N.mean(tg[ind])
nied = N.sum(rr[ind])
tmp = ws[ind]
wind = N.mean(tmp[tmp>-100])
tmp = tx[ind]
heta = len(tmp[tmp>30.])
eist = len(tmp[tmp<0.])
f.write('%5i,%7.2f,%7.2f,%7.2f,%7.2f,%7.2f\n'%(jo[j],tmit,nied,heta,eist,wind))
tmitja[j] = tmit
niedja[j] = nied
hetaja[j] = heta
eistja[j] = eist
windja[j] = wind
if(len(ind)<365):
nn = len(ind)
tmits[j,:nn] = N.cumsum(tg[ind][:nn])/float(nd)
tmita[j,:nn] = tg[ind][:nn]
nieds[j,:nn] = N.cumsum(rr[ind][:nn])
nieda[j,:nn] = rr[ind][:nn]
else:
tmits[j,:] = N.cumsum(tg[ind][0:365])/float(nd)
tmita[j,:] = tg[ind][0:365]
nieds[j,:] = N.cumsum(rr[ind][0:365])
nieda[j,:] = rr[ind][0:365]
f.write('%5s,%7s,%7s,%7s,%7s,%7s\n'%('NaN','NaN','NaN','NaN','NaN','NaN'))
slope,intercept, r_value, p_value, std_err = S.linregress(N.arange(nj),tmitja); slope,intercept; tmitja0 = intercept; tmitja1 = intercept+slope*nj
slope,intercept, r_value, p_value, std_err = S.linregress(N.arange(nj),niedja); niedja0 = intercept; niedja1 = intercept+slope*nj
slope,intercept, r_value, p_value, std_err = S.linregress(N.arange(nj),hetaja); hetaja0 = intercept; hetaja1 = intercept+slope*nj
slope,intercept, r_value, p_value, std_err = S.linregress(N.arange(nj),eistja); eistja0 = intercept; eistja1 = intercept+slope*nj
slope,intercept, r_value, p_value, std_err = S.linregress(N.arange(nj),windja); windja0 = intercept; windja1 = intercept+slope*nj
f.write('%5i,%7.2f,%7.2f,%7.2f,%7.2f,%7.2f\n'%(1961,tmitja0,niedja0,hetaja0,eistja0,windja0))
f.write('%5i,%7.2f,%7.2f,%7.2f,%7.2f,%7.2f\n'%(lastyear,tmitja1,niedja1,hetaja1,eistja1,windja1))
f.close()
########################################################
do = N.arange(0,nd+1,1)
tmits = N.insert(tmits,0,jo,axis=1)
tmits = N.insert(tmits,0,do,axis=0)
N.savetxt('potsdam2cum2tmit.csv',tmits.T, delimiter=",",fmt='%8.2f')
tmita = N.insert(tmita,0,jo,axis=1)
tmita = N.insert(tmita,0,do,axis=0)
N.savetxt('potsdam2abs2tmit.csv',tmita.T, delimiter=",",fmt='%8.2f')
nieds = N.insert(nieds,0,jo,axis=1)
nieds = N.insert(nieds,0,do,axis=0)
N.savetxt('potsdam2cum2nied.csv',nieds.T, delimiter=",",fmt='%8.2f')
nieda = N.insert(nieda,0,jo,axis=1)
nieda = N.insert(nieda,0,do,axis=0)
N.savetxt('potsdam2abs2nied.csv',nieda.T, delimiter=",",fmt='%8.2f')
#########################################################
xnew = N.arange(-20,41,1);nx = len(xnew)
ind = N.where((jj>=1961)&(jj<=1990))[0]
density0 = S.kde.gaussian_kde(tx[ind])
ynew0 = 365*density0(xnew)
ind = N.where((jj>=1989)&(jj<=(lastyear)))[0]
density1 = S.kde.gaussian_kde(tx[ind])
ynew1 = 365*density1(xnew)
f = open('potsdam2dist.csv','w')
f.write('%10s,%10s,%10s\n'%('TMAX','PAST','RECENT'))
for x in range(nx):
f.write('%10i,%10.2f,%10.2f\n'%(x,ynew0[x],ynew1[x]))
f.close()