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