#!/usr/bin/env python

### A. Spiga -- LMD -- 30/06/2011 to 10/07/2011
### Thanks to A. Colaitis for the parser trick


####################################
####################################
### The main program to plot vectors
def winds (namefile,\
           nvert,\
           proj=None,\
           back=None,\
           target=None,
           stride=3,\
           numplot=4,\
           var=None,\
           colorb=True,\
           winds=True,\
           addchar=None,\
           interv=[0,1],\
           vmin=None,\
           vmax=None,\
           tile=False,\
           zoom=None):

    ####################################################################################################################
    ### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png

    #################################
    ### Load librairies and functions
    from netCDF4 import Dataset
    from myplot import getcoord2d,define_proj,makeplotpng,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy,\
                       fmtvar,definecolorvec,getwinds,defcolorb,getprefix
    from mymath import deg,max,min,mean
    from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor
    from matplotlib.cm import get_cmap
    import numpy as np
    from numpy.core.defchararray import find

    ###
    #rcParams['text.usetex'] = True
    #rcParams['cairo.format'] = 'svg'

    ######################
    ### Load NETCDF object
    nc  = Dataset(namefile)

    ###################################
    ### Recognize predefined file types
    if 'controle' in nc.variables:   typefile = 'gcm'
    elif 'vert' in nc.variables:     typefile = 'mesoapi'
    elif 'U' in nc.variables:        typefile = 'meso'
    else:                            
        print "typefile not supported."
        print nc.variables
        exit()

    ##############################################################
    ### Try to guess the projection from wrfout if not set by user 
    if typefile in ['mesoapi','meso']:
        if proj == None:       proj = getproj(nc)
                                    ### (il faudrait passer CEN_LON dans la projection ?)
    elif typefile in ['gcm']:
        if proj == None:       proj = "cyl"   
                                    ## pb avec les autres (de trace derriere la sphere ?)

    ############################################
    #### Choose underlying topography by default 
    if not back:
        if not var:                                        back = "mola"    ## if no var:         draw mola
        elif typefile in ['mesoapi','meso'] \
             and proj not in ['merc','lcc','nsper','laea']:       back = "molabw"  ## if var but meso:   draw molabw
        else:                                              pass             ## else:              draw None

    ####################################################
    ### Get geographical coordinates and plot boundaries
    if typefile in ['mesoapi','meso']:
        [lon2d,lat2d] = getcoord2d(nc)
        lon2d = dumpbdy(lon2d)
        lat2d = dumpbdy(lat2d)
    elif typefile in ['gcm']:
        [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True)
    if proj == "npstere":             [wlon,wlat] = latinterv("North_Pole")
    elif proj in ["lcc","laea"]:      [wlon,wlat] = wrfinterv(lon2d,lat2d)
    else:                             [wlon,wlat] = simplinterv(lon2d,lat2d)
    if zoom:  
        dlon = abs(wlon[1]-wlon[0])/2.
        dlat = abs(wlat[1]-wlat[0])/2.
        [wlon,wlat] = [ [wlon[0]+zoom*dlon/100.,wlon[1]-zoom*dlon/100.],\
                        [wlat[0]+zoom*dlat/100.,wlat[1]-zoom*dlat/100.] ]
        print "zoom %",zoom,wlon,wlat

    ##############################################################################
    ### Get winds and know if those are meteorological winds (ie. zon, mer) or not
    if winds:
        if typefile is 'mesoapi': 
            [u,v] = getwinds(nc)
            metwind = True  ## meteorological (zon/mer)
        elif typefile is 'gcm':
            [u,v] = getwinds(nc,charu='u',charv='v')
            metwind = True  ## meteorological (zon/mer)
        elif typefile is 'meso':
            [u,v] = getwinds(nc,charu='U',charv='V')
            metwind = False ## geometrical (wrt grid)
            print "Beware ! Not using meteorological winds. You trust numerical grid as being (x,y)."

    #####################################################
    ### Load the chosen variables, whether it is 2D or 3D
    if var:
        if var not in nc.variables: 
            print "not found in file:",var
            exit()
        else:   
            dimension = np.array(nc.variables[var]).ndim
            if dimension == 2:     field = nc.variables[var][:,:]
            elif dimension == 3:   field = nc.variables[var][:,:,:]
            elif dimension == 4:   field = nc.variables[var][:,nvert,:,:]  
        dev = np.std(field)*3.0
        if vmin is None:  zevmin = mean(field) - dev
        else:             zevmin = vmin
        if vmax is None:  zevmax = mean(field) + dev
        else:             zevmax = vmax
        print "bounds ", zevmin, zevmax
        ### some already defined colormaps
        if colorb is True:    colorb = defcolorb(var)
    else:
        dimension = 0
 
    ###########################
    ### Get length of time axis
    if winds:               nt = len(u[:,0,0,0])
    elif var:  
        if dimension == 2:  nt = 1
        else             :  nt = len(field[:,0,0])

    #########################################
    ### Name for title and graphics save file
    if winds:   basename = "UV_" 
    else:       basename = ""
    if var:     basename = basename + var
    ###
    if dimension == 4:
        if typefile is 'meso':                      stralt = "_lvl" + str(nvert)
        elif typefile is 'mesoapi':  
            zelevel = int(nc.variables['vert'][nvert])
            if 'altitude'       in nc.dimensions:   stralt = "_"+str(zelevel)+"m-AMR"
            elif 'altitude_abg' in nc.dimensions:   stralt = "_"+str(zelevel)+"m-ALS"
            elif 'bottom_top'   in nc.dimensions:   stralt = "_"+str(zelevel)+"m"
            elif 'pressure'     in nc.dimensions:   stralt = "_"+str(zelevel)+"Pa" 
        else:                                       stralt = ""          
    else:
        stralt = ""
    ###
    basename = basename + stralt

    ##################################
    ### Open a figure and set subplots
    fig = figure()
    if   numplot > 0:    
        if   numplot == 4: 
            sub = 221
            fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
            rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
        elif numplot == 2: 
            sub = 121
            fig.subplots_adjust(wspace = 0.3)
            rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
        elif numplot == 3: 
            sub = 131
            fig.subplots_adjust(wspace = 0.5)
            rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
        elif numplot == 6: 
            sub = 231
            fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
            rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
        elif numplot == 8: 
            sub = 331 #241
            fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
            rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
        elif numplot == 9:
            sub = 331
            fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
            rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
        elif numplot == 1:
            sub = 99999
        else:
            print "supported: 1,2,3,4,6,8,9"
            exit()
        ### Prepare time loop
        if nt <= numplot or numplot == 1:  
            tabrange = [0]
            numplot = 1
        else:                          
            tabrange = range(0,nt,int(nt/numplot))  #nt-1
            tabrange = tabrange[0:numplot]
    else: 
        tabrange = range(0,nt,1)
        sub = 99999
    print tabrange

    #################################
    ### Time loop for plotting device
    found_lct = False
    for i in tabrange:

       ### Which local time ?
       ltst = ( interv[0] + 0.5*(wlon[0]+wlon[1])/15.) + i*interv[1]
       ltst = int (ltst * 10) / 10.
       ltst = ltst % 24

       ### General plot settings
       if numplot > 1: 
           subplot(sub)
           found_lct = True
       elif numplot == 1:
           found_lct = True 
            ### If only one local time is requested (numplot < 0)
       elif numplot <= 0: 
           if int(ltst) + numplot != 0:		continue
           else:                  		found_lct = True

       ### Map projection
       m = define_proj(proj,wlon,wlat,back=back)
       x, y = m(lon2d, lat2d)

       #### Contour plot
       if var:
           if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(field[i,:,:])
           elif typefile in ['gcm']:             
               if dimension == 2:                what_I_plot = field[:,:]
               elif dimension == 3:              what_I_plot = field[i,:,:]
           palette = get_cmap(name=colorb)
           #palette.set_over('b', 1.0)
           if not tile:
               zelevels = np.linspace(zevmin,zevmax)
               contourf( x, y, what_I_plot, 10, cmap = palette, levels = zelevels )
           else:    
               pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax )
           if var in ['HGT']:        pass
           elif colorb:              
                                     ndiv = 10
                                     colorbar(fraction=0.05,pad=0.1,format=fmtvar(var),\
                                              ticks=np.linspace(zevmin,zevmax,ndiv+1),\
                                              extend='max',spacing='proportional')
                                                     # both min max neither
       ### Vector plot
       if winds:
           if   typefile in ['mesoapi','meso']:    
               [vecx,vecy] = [dumpbdy(u[i,nvert,:,:]), dumpbdy(v[i,nvert,:,:])]
               key = True
           elif typefile in ['gcm']:               
               [vecx,vecy] = [        u[i,nvert,:,:] ,         v[i,nvert,:,:] ]
               key = False
           if metwind:  [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d)
           if var == None:        colorvec = definecolorvec(back)
           else:                  colorvec = definecolorvec(colorb)
           vectorfield(vecx, vecy,\
                      x, y, stride=stride, csmooth=stride,\
                      scale=15., factor=300., color=colorvec, key=key)
                                        #200.         ## or csmooth=2
       
       ### Next subplot
       plottitle = basename
       if addchar:  plottitle = plottitle + addchar + "_LT"+str(ltst)
       else:        plottitle = plottitle + "_LT"+str(ltst)
       ptitle( plottitle )
       sub += 1

    ##########################################################################
    ### Save the figure in a file in the data folder or an user-defined folder
    if typefile in ['meso','mesoapi']:   prefix = getprefix(nc)   
    elif typefile in ['gcm']:            prefix = 'LMD_GCM_'
    else:                                prefix = ''
    ###
    zeplot = prefix + basename 
    if addchar:         zeplot = zeplot + addchar
    if numplot <= 0:    zeplot = zeplot + "_LT"+str(abs(numplot))
    ###
    if not target:      zeplot = namefile[0:find(namefile,'wrfout')] + zeplot
    else:               zeplot = target + "/" + zeplot  
    ###
    if found_lct:     makeplotpng(zeplot,pad_inches_value=0.35)   
    else:             print "Local time not found"

    ###############
    ### Now the end
    return zeplot

###########################################################################################
###########################################################################################
### What is below relate to running the file as a command line executable (very convenient)
if __name__ == "__main__":
    import sys
    from optparse import OptionParser    ### to be replaced by argparse
    from api_wrapper import api_onelevel
    from netCDF4 import Dataset
    from myplot import getlschar

    #############################
    ### Get options and variables
    parser = OptionParser()
    parser.add_option('-f', action='store', dest='namefile',    type="string",  default=None,  help='[NEEDED] name of WRF file')
    parser.add_option('-l', action='store', dest='nvert',       type="float",   default=0,     help='vertical level (def=0)(-i 2: p,mbar)(-i 3,4: z,km)')
    parser.add_option('-p', action='store', dest='proj',        type="string",  default=None,  help='projection')
    parser.add_option('-b', action='store', dest='back',        type="string",  default=None,  help='background')
    parser.add_option('-t', action='store', dest='target',      type="string",  default=None,  help='destination folder')
    parser.add_option('-s', action='store', dest='stride',      type="int",     default=3,     help='stride vectors (def=3)')
    parser.add_option('-v', action='store', dest='var',         type="string",  default=None,  help='variable contoured')
    parser.add_option('-n', action='store', dest='numplot',     type="int",     default=4,     help='number of plots (def=1)(<0: 1 plot of LT -*numplot*)')
    parser.add_option('-i', action='store', dest='interp',      type="int",     default=None,  help='interpolation method (2: press, 3: z-amr, 4:z-als)')
    parser.add_option('-c', action='store', dest='colorb',      type="string",  default=True,  help='change colormap')
    parser.add_option('-x', action='store_false', dest='winds',                 default=True,  help='no wind vectors')
    parser.add_option('-m', action='store', dest='vmin',        type="float",   default=None,  help='bounding minimum value for color plot')    
    parser.add_option('-M', action='store', dest='vmax',        type="float",   default=None,  help='bounding maximum value for color plot') 
    parser.add_option('-T', action='store_true', dest='tile',                   default=False, help='draw a tiled plot (no blank zone)')
    parser.add_option('-z', action='store', dest='zoom',        type="float",   default=None,  help='zoom factor in %')
    parser.add_option('-N', action='store_true', dest='nocall',                 default=False, help='do not recreate api file')
    #parser.add_option('-V', action='store', dest='comb',        type="float",   default=None,  help='a defined combination of variables')
    (opt,args) = parser.parse_args()
    if opt.namefile is None: 
        print "I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h"
        exit()
    print "Options:", opt
    zefile = opt.namefile    
    zelevel = opt.nvert   
    stralt = None
    [lschar,zehour,zehourin] = getlschar ( zefile )  ## getlschar from wrfout (or simply return "" if another file)

    #####################################################
    ### Call Fortran routines for vertical interpolations
    if opt.interp is not None:
        if opt.nvert is 0 and opt.interp is 4:  zelevel = 0.010
        ### winds or no winds
        if opt.winds            :  zefields = 'uvmet'
        else                    :  zefields = ''
        ### var or no var
        if opt.var is None      :  pass
        elif zefields == ''     :  zefields = opt.var
        else                    :  zefields = zefields + "," + opt.var
        print zefields
        zefile = api_onelevel (  path_to_input   = '', \
                                 input_name      = zefile, \
                                 path_to_output  = opt.target, \
                                 fields          = zefields, \
                                 interp_method   = opt.interp, \
                                 onelevel        = zelevel, \
                                 nocall          = opt.nocall )
        print zefile
        zelevel = 0 ## so that zelevel could play again the role of nvert

    #############
    ### Main call
    name = winds (zefile,int(zelevel),\
           proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var,numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\
           addchar=lschar,interv=[zehour,zehourin],vmin=opt.vmin,vmax=opt.vmax,tile=opt.tile,zoom=opt.zoom)
    print 'Done: '+name

    #########################################################
    ### Generate a .sh file with the used command saved in it
    command = ""
    for arg in sys.argv: command = command + arg + ' '
    f = open(name+'.sh', 'w')
    f.write(command)