import numpy as N
import Image,ImageFont,ImageDraw,ImageOps

pfad="./"#"/storage/sdcard0/prediction/"

jj=N.genfromtxt(pfad+"tmean.txt" ,usecols=(0),skip_header=0,dtype="I")
mm=N.genfromtxt(pfad+"tmean.txt",usecols=(1),skip_header=0,dtype="I")
tt=N.genfromtxt(pfad+"tmean.txt",usecols=(3),skip_header=0,dtype="f")	
jx=N.genfromtxt(pfad+"tmax.txt",usecols=(0),skip_header=0,dtype="f")
tx=N.genfromtxt(pfad+"tmax.txt",usecols=(3),skip_header=0,dtype="f")

ja=1993
je=2017

jo=N.arange(ja,je+1,1)
mo=N.arange(1,13,1)
nj=len(jo)
nm=len(mo)

mx=N.zeros((nj,nm),float)
sx=N.zeros((nj,nm),float)
hx=N.zeros((nj),float)

for j in range(nj):
    for m in range(nm):
  
        w=N.ma.masked_where(jj!=jo[j],tt)
        w=N.ma.masked_where(mm!=mo[m],w)
        w=N.ma.masked_where(tt<-100.,w)
        w=w.compressed()
    
        mx[j,m]=N.mean(w)
        sx[j,m]=N.std(w)
   
    w=N.ma.masked_where(jj!=jo[j],tx)
    w=N.ma.masked_where(tx<30.,w)
    w=w.compressed()
    
    hx[j]=len(w)
    
print jo
    
my=N.diff(mx,axis=0)
sy=N.diff(sx,axis=0)
hy=N.diff(hx)

print '%5s'%('Year')

for j in range(nj-1):

   print '%5i'%(jo[j+1])

print nj,hy.shape,N.shape(hy[:nj-2])

#########

per=['JUN-AUG MEAN TEMPERATURE:','JUL-AUG MEAN TEMPERATURE:','HOT DAYS:']
fac=[-25.,-25.,-5.]

z=-1

for opt in ['jja','ja','heta']:

   z=z+1

   x=N.ones((nj-1,6),float)

   x[:,1]=sy[:,0]
   x[:,2]=sy[:,1]
   x[:,3]=sy[:,2]
   x[:,4]=sy[:,3]
   x[:,5]=my[:,3]

   if(opt=="jja"):y=N.mean(my[:,5:8],1)
   if(opt=="ja"):y=N.mean(my[:,6:8],1)
   if(opt=="heta"):y=hy

   a=N.linalg.lstsq(x[:nj-2,:],y[:nj-2])[0]

   #print a

   fit=N.dot(x[:nj-2,:],a)

   b=N.cumsum(a[1:]*x[:nj-2,1:],1)
   #print b.shape

   c=N.corrcoef(y[:nj-2],fit)
   co=c[0,1]**2

   xy=[]
   xz=[]

   for j in range(nj-1):
    
      if(j<nj-2):
         
         xy.append((40+j*20,140+fac[z]*y[j]))
         xz.append((40+j*20,140+fac[z]*fit[j]))

      else:
    
         fit=N.dot(x[nj-2,:],a)
         if(z==0): az = a      
         print "%4i: obs: %3.1f mod: %3.1f r2: %3.2f"%(jo[j+1],y[j],fit,co)
         con=N.cumsum(a[1:]*x[nj-2,1:])
   
         if(opt=='jja'): jja = 240.-((fit+N.mean(mx[nj-2,5:8]))-16.)*40.; jja0 = 240.-(N.mean(mx[nj-2,5:8])-16.)*40.; jja1 = 240.-(N.mean(mx[nj-1,5:8])-16.)*40.;print fit+N.mean(mx[nj-2,5:8])
         if(opt=='ja'): ja = 240.-((fit+N.mean(mx[nj-2,6:8]))-16.)*40.; ja0 = 240.-(N.mean(mx[nj-2,6:8])-16.)*40.; ja1 = 240.-(N.mean(mx[nj-1,6:8])-16.)*40.
         if(opt=='heta'): heta = 240.-((fit+hx[nj-2])-2.)*10.; heta0 = 240.-(hx[nj-2]-2.)*10.; heta1 = 240.-(hx[nj-1]-1.9)*10.
   
   fits = []
   
   for j in range(nj-3):
         
      xdum = x[:nj-2,:]
      ydum = y[:nj-2]      
            
      xdum = N.delete(xdum,j,0)
      ydum = N.delete(ydum,j)   
            
      aa=N.linalg.lstsq(xdum,ydum)[0]

      ff=N.dot(x[nj-2,:],aa)      
      
      fits.append(ff)   

   som = y

   ############

   im = Image.new('RGB',(500,280),(255,255,255))
   dr = ImageDraw.Draw(im)
   font=ImageFont.truetype(pfad+"arial.ttf",16)
   fonk=ImageFont.truetype(pfad+"arial.ttf",12)

   dr.rectangle((40,60,480,220),fill=(240,240,240,50))
   dr.rectangle((40,80,480,200),fill=(220,220,220,50))
   dr.rectangle((40,100,480,180),fill=(200,200,200,50))
   dr.rectangle((40,120,480,160),fill=(180,180,180,50))

   dr.rectangle((20,60,40,80),fill=(240,0,0,50))
   dr.rectangle((20,80,40,100),fill=(200,0,0,50))
   dr.rectangle((20,100,40,120),fill=(160,0,0,50))
   dr.rectangle((20,120,40,140),fill=(120,0,0,50))

   dr.rectangle((20,140,40,160),fill=(0,0,120,50))
   dr.rectangle((20,160,40,180),fill=(0,0,160,50))
   dr.rectangle((20,180,40,200),fill=(0,0,200,50))
   dr.rectangle((20,200,40,220),fill=(0,0,240,50))

   for x in range(20,500,20):
      gl=[]
      for y in range(40,260,20):
         gl.append((x,y))
      dr.line(gl,"rgb(150,150,150)")   
    
   for y in range(40,260,20):
      gl=[]
      for x in range(20,500,20):
         gl.append((x,y))
      dr.line(gl,"rgb(220,220,220)")  
     
   gl=[]
   for x in range(20,500,20):
      gl.append((x,140))
   dr.line(gl,"rgb(0,0,0)")  
         
   gl=[]
   for y in range(40,260,20):
      gl.append((40,y))
   dr.line(gl,"rgb(0,0,0)")   

   dr.line(xy,"rgb(255,0,0)")
   dr.line(xz,"rgb(0,0,255)")   
   
   for ik in fits:
   
       dr.ellipse((39+(nj-2)*20,139+fac[z]*ik,41+(nj-2)*20,141+fac[z]*ik),fill='grey')
   
   dr.ellipse((39+(nj-2)*20,139+fac[z]*fit,41+(nj-2)*20,141+fac[z]*fit),fill='black')
   dr.text((40+(nj-2)*20,125+fac[z]*fit),str(N.round(fit,1)),fill=(0,0,0),font=fonk)

   dr.text((360,40),"OBS",fill=(255,0,0),font=font) 
   dr.text((400,40),"SIM",fill=(0,0,255),font=font)
   dr.text((440,40),"EST",fill=(0,0,0),font=font)         
   
   dr.text((40,20),per[z]+" POTSDAM (r2="+str(N.round(co,2))+")",fill=(0,0,0,100),font=font)
   dr.text((40,250),"(c) peterh@pik-potsdam.de",fill=(200,200,200,50),font=font)

   for j in xrange(0,nj,5):

      dr.text((22+20*j,220),str(jo[j+1]),fill=(0,0,0),font=font)

   txt=Image.new('L',(280,280)) 
   d = ImageDraw.Draw(txt) 
   d.text((80,0),"DIFFERENTIALS",font=font,fill=200) 
   w=txt.rotate(90,expand=1)

   im.paste(ImageOps.colorize(w,(0,0,0),(0,0,0)),(0,0),w)

   im.save(pfad+"pred_"+opt+".png")

   ##########

   im = Image.new('RGB',(280,280),(255,255,255))
   dr = ImageDraw.Draw(im)
   font=ImageFont.truetype(pfad+"arial.ttf",14)

   dr.rectangle((40,40,240,240),fill=(240,240,240,50))
   dr.rectangle((90,40,240,240),fill=(220,220,220,50))
   dr.rectangle((140,40,240,240),fill=(200,200,200,50))
   dr.rectangle((190,40,240,240),fill=(180,180,180,50))

   for y in range(40,260,20):
      gl=[]
      for x in range(40,290,50):
         gl.append((x,y))
      dr.line(gl,"rgb(150,150,150)")   
  
   for x in range(40,290,50):
      gl=[]
      for y in range(40,260,20):
         gl.append((x,y))
      dr.line(gl,"rgb(220,220,220)")  
       
   gl=[]
   for x in range(40,290,50):
      gl.append((x,140))
   dr.line(gl,"rgb(0,0,0)")   

   for j in range(nj-2):
      xy=[]
      for k in range(4):
         xy.append((40+50*k,140+fac[z]*b[j,k]))
      xy.append((40+50*4,140+fac[z]*som[j]))
  
      if(som[j]<-0.5):
    
         dr.line(xy,"rgb(0,0,255)")
         dr.text((245,140+fac[z]*som[j]),str(jo[j+1]),fill="rgb(0,0,255)",font=font)

      if(som[j]>0.5):
   
         dr.line(xy,"rgb(255,0,0)")
         dr.text((245,140+fac[z]*som[j]),str(jo[j+1]),fill="rgb(255,0,0)",font=font)

   for k in range(4):

      dr.ellipse((39+50*k,139.+fac[z]*con[k],41+50*k,141.+fac[z]*con[k]),fill='black')

   dr.text((50,52),"+0.5K",fill=(255,0,0,100),font=font) 
   dr.text((50,72),"-0.5K",fill=(0,0,255,100),font=font)
   dr.text((40,20),"CONTRIBUTIONS: POTSDAM",fill=(0,0,0),font=font)
   dr.text((40,250),"(c) peterh@pik-potsdam.de",fill=(200,200,200),font=font)

   dr.text((30,220),"Jan",fill=(0,0,0),font=font)
   dr.text((80,220),"Feb",fill=(0,0,0),font=font)
   dr.text((130,220),"Mar",fill=(0,0,0),font=font)
   dr.text((180,220),"Apr",fill=(0,0,0),font=font)
   dr.text((230,220),"JJA",fill=(0,0,0),font=font)

   txt=Image.new('L',(280,280)) 
   d = ImageDraw.Draw(txt) 
   d.text((80,10),"DIFFERENTIALS",font=font,fill=200) 
   w=txt.rotate(90,expand=1)

   im.paste(ImageOps.colorize(w,(0,0,0),(0,0,0)),(0,0),w)
   im.save(pfad+"cont_"+opt+".png")

###############

im = Image.new('RGB',(280,280),(255,255,255))
dr = ImageDraw.Draw(im)
font=ImageFont.truetype(pfad+"arial.ttf",12)

dr.rectangle((50,jja,80,240),fill=(251,192,45))
dr.rectangle((120,ja,150,240),fill=(230,74,25))
dr.rectangle((190,heta,220,240),fill=(223,31,162))

dr.line(((50,jja0),(80,jja0)),"rgb(251,192,45)")
dr.line(((120,ja0),(150,ja0)),"rgb(230,74,25)")
dr.line(((190,heta0),(220,heta0)),"rgb(223,31,162)")

dr.line(((50,jja1),(80,jja1)),"rgb(0,0,0)")
dr.line(((120,ja1),(150,ja1)),"rgb(0,0,0)")
dr.line(((190,heta1),(220,heta1)),"rgb(0,0,0)")

tem = ['16.0','16.5','17.0','17.5','18.0','18.5','19.0','19.5','20.0','20.5','21.0']
tem = tem[::-1]

day = ['2','4','6','8','10','12','14','16','18','20','22']
day = day[::-1]

i = -1

for y in range(40,260,20):
   
   i = i+1
   
   gl=[]
   for x in range(40,290,50):
      gl.append((x,y))
   dr.line(gl,"rgb(150,150,150)")

   if(y>40):

      dr.text((85,y-15),tem[i],fill=(0,0,0),font=font)
      dr.text((155,y-15),tem[i],fill=(0,0,0),font=font)
      dr.text((225,y-15),day[i],fill=(0,0,0),font=font)

dr.text((45,20),"Jun-Aug",fill=(0,0,0),font=font)
dr.text((115,20),"Jul-Aug",fill=(0,0,0),font=font)
dr.text((185,20),"Hot-Days",fill=(0,0,0),font=font)

dr.text((40,250),"(c) peterh@pik-potsdam.de",fill=(200,200,200),font=font)

font=ImageFont.truetype(pfad+"arial.ttf",20)

txt=Image.new('L',(280,280)) 
d = ImageDraw.Draw(txt) 
d.text((80,10),"2015 vs. 2016",font=font,fill=200) 
w=txt.rotate(90,expand=1)

im.paste(ImageOps.colorize(w,(0,0,0),(0,0,0)),(0,0),w)

im.save(pfad+"pred_all.png")

###############

f=open(pfad+'app.html','w')
f.write('<html>\n')
f.write('<body bgcolor="Khaki">\n')
#f.write('<h1><a name="L0">Seasonal Prediction</a></h1>\n')

f.write('</table>\n')
f.write('<h1>Seasonal Prediction: Summary</a></h1>\n')
f.write('<table border=1 bgcolor="#82b1ff">\n')
f.write('<td><img src="%s" height=345></td>\n'%(pfad+'pred_all.png'))
f.write('<td><img src="%s" height=345></td>\n'%(pfad+'1985-2014.png'))
f.write('</table>\n')

f.write('<h1><a name="L0">Background</a></h1>\n')
f.write('<table border=1 bgcolor="#82b1ff">\n')
f.write('<colgroup><col width="200"><col width="495"></colgroup>')
f.write('<tr align="left"><th>Method:</th><td>Univariate Linear Regression</td></tr>\n')
f.write('<tr align="left"><th>Station:</th><td>Potsdam (Germany)</td></tr>\n')
f.write('<tr align="left"><th>Parameter:</th><td>Daily Mean Temperatur (Tmit)</td></tr>\n')
f.write('<tr align="left"><th><a href="#L1">Aggregation:</a></th><td>Monthly Means (m)</td></tr>\n')
f.write('<tr align="left"><th></th><td>Monthly Standard Deviations (&sigma;)</td></tr>\n')
f.write('<tr align="left"><th><a href="#L2">Detrending:</a></th><td>Differentials</td></tr>\n')
f.write('<tr align="left"><th>Target Values:</th><td><a href="#L3">Summer Mean Temperature (Jun-Aug)</a></td></tr>\n')
f.write('<tr align="left"><th></th><td><a href="#L4">Mean Temperature (Jul-Aug)</a></td></tr>\n')
f.write('<tr align="left"><th></th><td><a href="#L5">Number of Hot Days (hoda)</a></td></tr>\n')
f.write('<tr align="left"><th>Formula:</th><td>T<sub>JJA</sub>=a&bull;&sigma;{T<sub>Jan</sub>}+b&bull;&sigma;{T<sub>Feb</sub>}+c&bull;&sigma;{T<sub>Mar</sub>}+d&bull;&sigma;{T<sub>Apr</sub>}+e&bull;m{T<sub>Apr</sub>}</td></tr>\n')
f.write('<tr align="left"><th>Coefficients:</th><td>a = %3.1f</td></tr>\n'%(az[1]))
f.write('<tr align="left"><th></th><td>b = %3.1f</td></tr>\n'%(az[2]))
f.write('<tr align="left"><th></th><td>c = %3.1f</td></tr>\n'%(az[3]))
f.write('<tr align="left"><th></th><td>d = %3.1f</td></tr>\n'%(az[4]))
f.write('<tr align="left"><th></th><td>e = %3.1f</td></tr>\n'%(az[5]))
f.write('</table>\n')

f.write('<h1><a name="L1">Absolute Values (OBS)</a></h1><a href="#L0">TOP</a>\n')
f.write('<table border=1 bgcolor="#82b1ff">\n')
f.write('<colgroup><col width="60"><col width="120"><col width="120"><col width="120"><col width="60"><col width="60"><col width="60"><col width="60"></colgroup>')
f.write('<tr><th>Potsdam</th><th colspan="3">Targets</th><th colspan="4">Predictors</th></tr>\n') 
f.write('<tr><th>Year</th><th>Jun-Aug (&deg;C)</th><th>Jul-Aug (&deg;C)</th><th>Hot-Days (d)</th><th>sFeb</th><th>sMar</th><th>sApr</th><th>mApr</th></tr>\n')

for j in range(nj):

   f.write('<tr align=right>\n')
   f.write('<td>%5i</td><td>%5.1f</td><td>%5.1f</td><td>%3i</td><td>%3.1f</td><td>%3.1f</td><td>%3.1f</td><td>%3.1f</td>\n'%(jo[j],N.mean(mx[j,5:8]),N.mean(mx[j,6:8]),hx[j],sx[j,1],sx[j,2],sx[j,3],mx[j,3]))
   f.write('</tr>\n')

f.write('</table>\n')

f.write('<h1><a name="L2">Differentials (OBS)</a></h1><a href="#L0">TOP</a>\n')
f.write('<table border=1 bgcolor="#82b1ff">\n')
f.write('<colgroup><col width="60"><col width="120"><col width="120"><col width="120"><col width="60"><col width="60"><col width="60"><col width="60"></colgroup>')
f.write('<tr><th>Potsdam</th><th colspan="3">Targets</th><th colspan="4">Predictors</th></tr>\n')  
f.write('<tr><th>Year</th><th>Jun-Aug (K)</th><th>Jul-Aug (K)</th><th>Hot-Days (d)</th><th>sFeb</th><th>sMar</th><th>sApr</th><th>mApr</th></tr>\n')

for j in range(nj-1):

   f.write('<tr align=right>\n')
   f.write('<td>%5i</td><td>%5.1f</td><td>%5.1f</td><td>%3i</td><td>%3.1f</td><td>%3.1f</td><td>%3.1f</td><td>%3.1f</td>\n'%(jo[j+1],N.mean(my[j,5:8]),N.mean(my[j,6:8]),hy[j],sy[j,1],sy[j,2],sy[j,3],my[j,3]))
   f.write('</tr>\n')

f.write('</table>\n')

f.write('<h1><a name="L3">Estimate: Jun-Aug</a></h1><a href="#L0">TOP</a>\n')
f.write('<table border=1 bgcolor="#82b1ff">\n')
f.write('<td><img src="%s"></td>\n'%(pfad+'pred_jja.png'))
f.write('<td><img src="%s"></td>\n'%(pfad+'cont_jja.png'))
f.write('</table>\n')

f.write('<h1><a name="L4">Estimate: Jul-Aug</a></h1><a href="#L0">TOP</a>\n')
f.write('<table border=1 bgcolor="#82b1ff">\n')
f.write('<td><img src="%s"></td>\n'%(pfad+'pred_ja.png'))
f.write('<td><img src="%s"></td>\n'%(pfad+'cont_ja.png'))
f.write('</table>\n')

f.write('<h1><a name="L5">Estimate: Hot-Days</a></h1><a href="#L0">TOP</a>\n')
f.write('<table border=1 bgcolor="#82b1ff">\n')
f.write('<td><img src="%s"></td>\n'%(pfad+'pred_heta.png'))
f.write('<td><img src="%s"></td>\n'%(pfad+'cont_heta.png'))
f.write('</table>\n')

f.write('</body>\n')
f.write('</html>\n')


f.close()