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()