from mpl_toolkits.basemap import Basemap
from scipy.io.numpyio import fread
import matplotlib.pyplot as plt
import matplotlib.numerix.ma as M
import matplotlib.colors as C
from pylab import *
import numpy as np
import scipy as sc
import os, fnmatch
import array
import sys
import datetime
 
d=str(datetime.date.today())
d=d.replace('-', '')
 
path='c:/' #sys.argv[1]
for fn in os.listdir(path):
    if fnmatch.fnmatch(fn, d + '*.b*'):
        fileName=path + fn
        break
 
masking=True
treshold=0
savename=d + '_SO2.png'
dpi=120
smoothness=1
projection='mill'
titre='SO2'
legende='Delta Brightness Temperature (K)'
centre=[-90,0]
 
print('reading binary data')
Lat=[]
Lon=[]
Val=[]
 
nbreligne=long(os.stat(fileName)[6])/(8*int(fileName[-2:]))
rawfile=np.fromfile(open(fileName,'rb'),'d',-1)        
Lat=rawfile[0:nbreligne]
Lon=rawfile[nbreligne:nbreligne*2]
Val=rawfile[nbreligne*21:nbreligne*22]    
 
if(masking==True):
    i=0
    while(i<nbreligne):
        if(Val[i]>-1.1):
            Val[i]=-1.1
        i+=1
 
print('plotting data')
map=Basemap(projection='mill',lat_0=-90,lon_0=0,llcrnrlat=-90,urcrnrlat=90,\
            urcrnrlon=180,llcrnrlon=-180,resolution='i',area_thresh=30000.)
 
xi=np.linspace(-180,180,360*smoothness) #1440
yi=np.linspace(-90,90,180*smoothness)    #720
zi=griddata(Lon,Lat,Val,xi,yi)
 
topodat=map.transform_scalar(zi,xi,yi,720,1440)
map.imshow(topodat,cm.winter, interpolation='bilinear',aspect='auto',norm = C.Normalize(vmin = -5., vmax = -1.2, clip = True),extent=[-180,180,-90,90],filternorm=1)
 
cb=plt.colorbar(shrink=0.7)
cb.ax.set_ylabel(legende,fontsize=11)
for t in cb.ax.get_yticklabels():
    t.set_fontsize(7)
 
meridians = arange(-180,180,60)
parallels = arange(-90,90,30)
map.drawparallels(parallels,labels=[1,0,0,0],fontsize=7,linewidth=0.25)
map.drawmeridians(meridians,labels=[0,0,0,1],fontsize=7,linewidth=0.25)
 
title(titre)
 
print('saving map')        
map.drawcoastlines(0.25,antialiased=1)
plt.savefig(savename,dpi=dpi)
print('done')