Ignore:
Timestamp:
Nov 5, 2011, 6:54:37 PM (13 years ago)
Author:
aslmd
Message:

PYTHON GRAPHICS: A REFERENCE VERSION FOR BOTH GCM AND MESOSCALE USE!

Merged recent changes by TN for better GCM handling + 1D/2D plots with AS mesoscale version [winds.py]

Now

  • there is a common function named planetoplot in planetoplot.py which supports both GCM and mesoscale plots

    this will be a common reference for future development

  • mesoscale users might want to use the script meso.py
  • GCM users might want to use the script gcm.py

Todo

  • a few capabilities that were in winds.py are broken in meso.py, this will be fixed soon. the 'old' reference script winds.py is still here anyway!
  • there is still a bug with 'ortho' projection with GCM fields. still searching...

Misc

  • A few modifications were done in myplot.py and mymath.py
  • Added a simple stupid little.py script illustrating very simple use of Python graphics
  • Added the script doing statistical analysis for dust devils
Location:
trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON
Files:
3 added
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/examples/scripts/find_devils.py

    r310 r349  
    2727   
    2828    ### LOAD NETCDF DATA
    29     filename = "psfc.nc"
     29    #filename = "/home/aymeric/Big_Data/psfc.nc"
    3030    nc = Dataset(filename)
    3131    psfc = nc.variables["PSFC"]
    … …  
    3636    allsizesx = []
    3737    allsizesy = []
     38    depression = []
    3839    stride = 1 #5
    3940    for i in range(0,shape[0],stride):
    … …  
    5253        fac = 3.5
    5354        #fac = 2.5
     55        fac = 3.
    5456        lim = ave - fac*std
    5557        where = np.where(psfc2d < lim)
    5658        ############### END CRITERION
     59
     60        depression = np.append(depression,np.ravel(psfc2d[where])-ave)
    5761   
    5862        lab = np.zeros(np.array(psfc2d).shape) ## points to be treated by the minimum_position routine
    … …  
    115119    allsizesy.sort()
    116120   
    117     return allsizesx, allsizesy
     121    return allsizesx, allsizesy, depression
    118122
    119123#########################################################################
    … …  
    125129import matplotlib.mlab as mlab
    126130import mymath as mym
     131import myplot as myp
     132
     133import plfit
     134import plplot
     135import randht
     136import plpva
    127137
    128138save = True
    129139save = False
     140pression = False
     141pression = True
    130142
    131143if save:
    132     allsizesx, allsizesy = getsize("psfc.nc")
     144    allsizesx, allsizesy, depression = getsize("/home/aymeric/Big_Data/psfc.nc")
    133145    ### sauvegarde texte pour inspection
    134146    mym.writeascii(allsizesx,'allsizex.txt')
    135147    mym.writeascii(allsizesy,'allsizey.txt')
     148    mym.writeascii(depression,'alldepression.txt')
    136149    ### sauvegarde binaire pour utilisation python
    137150    myfile = open('allsizex.bin', 'wb')
    … …  
    141154    pickle.dump(allsizesy, myfile)
    142155    myfile.close()
     156    myfile = open('alldepression.bin', 'wb')
     157    pickle.dump(depression, myfile)
     158    myfile.close()
    143159
    144160### load files
    … …  
    147163myfile = open('allsizey.bin', 'r')
    148164allsizesy = pickle.load(myfile)
    149 
    150 ### append sizes
    151 plothist = np.append(allsizesx,allsizesy)
     165myfile = open('alldepression.bin', 'r')
     166depression = pickle.load(myfile)
     167depression = np.array(abs(depression))#*1000.
     168
     169### sizes
     170#plothist = np.append(allsizesx,allsizesy)
    152171plothist = allsizesx
    153 #plothist = allsizesy
     172if pression: plothist = depression
    154173plothist.sort()
    155 
    156 ### make bins
     174print 'mean ', np.mean(plothist,dtype=np.float64) 
     175print 'std ', np.std(plothist,dtype=np.float64)
     176print 'max ', np.max(plothist)
     177print 'min ', np.min(plothist)
     178print 'len ', len(plothist)
     179
     180
     181### MAKE BINS
    157182nbins = 100
    158183zebins = [2.0]
    … …  
    165190zebins = [2./np.sqrt(2.)]  ## ne pas tomber sur une dizaine ronde
    166191nbins = 200
    167 zebins = [0.2/np.sqrt(2.)]  ## ne pas tomber sur une dizaine ronde
     192
     193if pression: zebins = [0.3]
    168194
    169195for i in range(0,nbins):  zebins.append(zebins[i]*np.sqrt(2))
    170196zebins = np.array(zebins)
    171 ### select reasonable bins for DD
    172 zebins = zebins [ zebins > 10. ]
    173 zebins = zebins [ zebins < 1000. ]
    174 #print 'bins ',zebins
     197#### select reasonable bins for DD
     198if not pression:
     199    zebins = zebins [ zebins > 15. ]
     200    zebins = zebins [ zebins < 1000. ]
     201else:
     202    zebins = zebins [ zebins < 10. ]
    175203print 'corrected bins ',zebins
    176 print 'max value ',mym.max(plothist),mym.max(allsizesx),mym.max(allsizesy)
    177 
    178 #plothist = [20.,30.,40.,50.,70.,100.,130.,190.,260.,370.,520.]  ## pour calibrer
    179204
    180205#### HISTOGRAM
     206plt.figure(1)
    181207plt.hist(    plothist,\
    182208             log=True,\
    183209             bins=zebins,\
    184210#             cumulative=-1,\
     211             normed=True,\
    185212        )
    186213plt.xscale('log')
    187 plt.show()
    188 
     214if pression:    plt.xlabel('Pressure (Pa)')
     215else:           plt.xlabel('Diameter (m)')
     216plt.ylabel('Population (normalized)')
     217if pression: prefix="p"
     218else:        prefix=""
     219myp.makeplotres(prefix+"histogram",res=200,disp=False)
     220plt.close(1)
     221
     222### COMPARED HISTOGRAMS
     223### --- FIT WITH POWER LAW
     224if pression: [alpha, xmin, L] = plfit.plfit(plothist,'xmin',0.3)
     225else:        [alpha, xmin, L] = plfit.plfit(plothist,'limit',20.)
     226print alpha,xmin
     227
     228#a = plpva.plpva(plothist,0.75,'xmin',0.75)
     229#print a
     230
     231#### DEUXIEME ROUTINE
     232####IL FAUT UTILISER LE DISCRET POUR LA TAILLE !!!
     233#if pression:   myplfit = plfit.plfit(plothist,verbose=True,xmin=0.75)
     234#else:          myplfit = plfit.plfit(plothist,verbose=True,xmin=20.)
     235#myplfit.plotppf()
     236#plt.show()
     237#exit()
     238
     239
     240#plt.figure(1)
     241#h = plplot.plplot(plothist,xmin,alpha)
     242#myp.makeplotres(prefix+"fit",res=200,disp=False)
     243plt.figure(2)
     244### --- POWER LAW (factor does not really matter)
     245power = (xmin/2.2)*np.array(randht.randht(10000,'powerlaw',alpha))
     246#power = (xmin/2.2)*np.array(randht.randht(10000,'cutoff',alpha,10.)) ##marche pas si trop grand
     247print 'mean ', np.mean(power,dtype=np.float64)
     248### --- EXPONENTIAL LAW
     249expo = randht.randht(10000,'exponential',1./(np.mean(power,dtype=np.float64)*1.00))
     250print 'mean ', np.mean(expo,dtype=np.float64)
     251### --- PLOT
     252plt.hist(    [plothist,power,expo],\
     253             label=['LES vortices','Power law '+'{:.1f}'.format(alpha),'Exponential law'],\
     254             log=True,\
     255             bins=zebins,\
     256#            cumulative=-1,\
     257             normed=True,\
     258        )
     259plt.legend()
     260plt.xscale('log')
     261if pression:    plt.xlabel('Pressure (Pa)')
     262else:           plt.xlabel('Diameter (m)')
     263plt.ylabel('Population (normalized)')
     264myp.makeplotres(prefix+"comparison",res=200,disp=False)
     265plt.close(2)
     266
     267########################
     268########################
     269zebins = [30.,42.,60.,84.,120.,170.,240.,340.]
     270plothist = []
     271plothist = np.append(plothist,30 *np.ones(306))
     272plothist = np.append(plothist,42 *np.ones(58) )
     273plothist = np.append(plothist,60 *np.ones(66) )
     274plothist = np.append(plothist,84 *np.ones(41) )
     275plothist = np.append(plothist,120*np.ones(19) )
     276plothist = np.append(plothist,170*np.ones(9)  )
     277plothist = np.append(plothist,240*np.ones(2)  )
     278plothist = np.append(plothist,340*np.ones(1)  )
     279
     280#zebins = [50.,71.,100.,141.,200.,282.,400.]
     281#plothist = []
     282#plothist = np.append(plothist,50. *np.ones(36))
     283#plothist = np.append(plothist,71. *np.ones(18) )
     284#plothist = np.append(plothist,100. *np.ones(12) )
     285#plothist = np.append(plothist,141. *np.ones(6) )
     286#plothist = np.append(plothist,200.*np.ones(4) )
     287#plothist = np.append(plothist,282.*np.ones(1)  )
     288#plothist = np.append(plothist,400.*np.ones(2)  )
     289
     290plt.figure(3)
     291[alpha, xmin, L] = plfit.plfit(plothist,'xmin',30)#50.)
     292print alpha,xmin
     293#a = plpva.plpva(plothist,30,'xmin',30)
     294h = plplot.plplot(plothist,xmin,alpha)
     295plt.loglog(h[0], h[1], 'k--',linewidth=2)
     296plt.hist(    plothist,\
     297             log=True,\
     298             bins=zebins,\
     299#             cumulative=-1,\
     300             normed=True,\
     301        )
     302plt.xscale('log')
     303plt.xlabel('Pressure (micro Pa)')
     304plt.ylabel('Population (normalized)')
     305myp.makeplotres("data",res=200,disp=False)
     306
     307#plt.figure(4)
     308#[alpha, xmin, L] = plfit.plfit(plothist,'xmin',50.) #,'xmin',30.)
     309#print alpha,xmin
     310#h = plplot.plplot(plothist,xmin,alpha)
     311#myp.makeplotres("datafit",res=200,disp=False)

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/mymath.py

    r345 r349  
    11def min (field,axis=None):
    22        import numpy as np
    3         return np.array(field).min(axis=axis)
     3        if field is None: return None
     4        else: return np.array(field).min(axis=axis)
    45
    56def max (field,axis=None):
    … …  
    1011def mean (field,axis=None):
    1112        import numpy as np
    12         return np.array(field).mean(axis=axis)
     13        if field is None: return None
     14        else: return np.array(field).mean(axis=axis)
    1315
    1416def deg ():

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py

    r345 r349  
    55    exit()
    66    return
     7
     8## Author: AS
     9def adjust_length (tab, zelen):
     10    from numpy import ones
     11    if tab is None:
     12        outtab = ones(zelen) * -999999
     13    else:
     14        if zelen != len(tab):
     15            print "not enough or too much values... setting same values all variables"
     16            outtab = ones(zelen) * tab[0]
     17        else:
     18            outtab = tab
     19    return outtab
    720
    821## Author: AS
    … …  
    4861    ### it would be actually better to name d4 d3 d2 d1 as t z y x
    4962    import numpy as np
    50     from mymath import max
     63    from mymath import max,mean
    5164    dimension = np.array(input).ndim
    5265    shape = np.array(input).shape
     66    #print 'd1,d2,d3,d4: ',d1,d2,d3,d4
    5367    print 'dim,shape: ',dimension,shape
    5468    output = input
    … …  
    7185        elif max(d1) >= shape[2]: error = True
    7286        elif d4 is not None and d2 is not None and d1 is not None:  output = input[d4,d2,d1]
    73         elif d4 is not None and d2 is not None:    output=np.mean(input[d4,:,:],axis=0); output=np.mean(output[d2,:],axis=0)
    74         elif d4 is not None and d1 is not None:    output=np.mean(input[d4,:,:],0); output=np.mean(output[:,d1],1)
    75         elif d2 is not None and d1 is not None:    output=np.mean(input[:,d2,:],axis=1); output=np.mean(output[:,d1],axis=1)
    76         elif d1 is not None:                output = np.mean(input[:,:,d1],axis=2)
    77         elif d2 is not None:                output = np.mean(input[:,d2,:],axis=1)
    78         elif d4 is not None:                output = np.mean(input[d4,:,:],axis=0)
     87        elif d4 is not None and d2 is not None:    output = mean(input[d4,:,:],axis=0); output=mean(output[d2,:],axis=0)
     88        elif d4 is not None and d1 is not None:    output = mean(input[d4,:,:],axis=0); output=mean(output[:,d1],axis=1)
     89        elif d2 is not None and d1 is not None:    output = mean(input[:,d2,:],axis=1); output=mean(output[:,d1],axis=1)
     90        elif d1 is not None:                       output = mean(input[:,:,d1],axis=2)
     91        elif d2 is not None:                       output = mean(input[:,d2,:],axis=1)
     92        elif d4 is not None:                       output = mean(input[d4,:,:],axis=0)
    7993    elif dimension == 4:
    8094        if   max(d4) >= shape[0]: error = True
    … …  
    87101        elif d4 is not None and d2 is not None and d1 is not None:         output = input[d4,:,d2,d1]
    88102        elif d3 is not None and d2 is not None and d1 is not None:         output = input[:,d3,d2,d1]
    89         elif d4 is not None and d3 is not None:  output = np.mean(input[d4,:,:,:],0); output = np.mean(output[d3,:,:],0)
    90         elif d4 is not None and d2 is not None:  output = np.mean(input[d4,:,:,:],0); output = np.mean(output[:,d2,:],1)
    91         elif d4 is not None and d1 is not None:  output = np.mean(input[d4,:,:,:],0); output = np.mean(output[:,:,d1],2)
    92         elif d3 is not None and d2 is not None:  output = np.mean(input[:,d3,:,:],1); output = np.mean(output[:,d2,:],1)
    93         elif d3 is not None and d1 is not None:  output = np.mean(input[:,d3,:,:],1); output = np.mean(output[:,:,d1],0)
    94         elif d2 is not None and d1 is not None:  output = np.mean(input[:,:,d2,:],2); output = np.mean(output[:,:,d1],2)
     103        elif d4 is not None and d3 is not None:  output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0)
     104        elif d4 is not None and d2 is not None:  output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,d2,:],axis=1)
     105        elif d4 is not None and d1 is not None:  output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,:,d1],axis=2)
     106        elif d3 is not None and d2 is not None:  output = mean(input[:,d3,:,:],axis=1); output = mean(output[:,d2,:],axis=1)
     107        elif d3 is not None and d1 is not None:  output = mean(input[:,d3,:,:],axis=1); output = mean(output[:,:,d1],axis=0)
     108        elif d2 is not None and d1 is not None:  output = mean(input[:,:,d2,:],axis=2); output = mean(output[:,:,d1],axis=2)
    95109        elif d1 is not None:                       output = input[:,:,:,d1]
    96110        elif d2 is not None:                       output = input[:,:,d2,:]
    … …  
    759773# output : all indexes to be taken into account for reducing field
    760774  import numpy as np
     775  if ( np.array(axis).ndim == 2):
     776      axis = axis[:,0]
    761777  if saxis is None:
    762778      zeindex = None
    … …  
    820836   #print "define_axis", x, y
    821837   return what_I_plot,x,y
     838
     839# Author: TN + AS
     840def determineplot(slon, slat, svert, stime):
     841    nlon = 1 # number of longitudinal slices -- 1 is None
     842    nlat = 1
     843    nvert = 1
     844    ntime = 1
     845    nslices = 1
     846    if slon is not None:
     847        nslices = nslices*len(slon)
     848        nlon = len(slon)
     849    if slat is not None:
     850        nslices = nslices*len(slat)
     851        nlat = len(slat)
     852    if svert is not None:
     853        nslices = nslices*len(svert)
     854        nvert = len(svert)
     855    if stime is not None:
     856        nslices = nslices*len(stime)
     857        ntime = len(stime)
     858    #else:
     859    #    nslices = 2 
     860
     861    mapmode = 0
     862    if slon is None and slat is None:
     863       mapmode = 1 # in this case we plot a map, with the given projection
     864    #elif proj is not None:
     865    #   print "WARNING: you specified a", proj,\
     866    #   "projection but asked for slices", slon,"in longitude and",slat,"in latitude"
     867    print "mapmode: ", mapmode
     868
     869    return nlon, nlat, nvert, ntime, mapmode, nslices

  • trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/planetoplot.py

    r345 r349  
    1 #!/usr/bin/env python
    2 
    3 ### A. Spiga  -- LMD -- 30/06/2011 to 10/07/2011
    4 ### T.Navarro -- LMD -- 10/2011
    5 ### Thanks to A. Colaitis for the parser trick
    6 
    7 
    8 ####################################
    9 ####################################
    10 ### The main program to plot vectors
     1#######################
     2##### PLANETOPLOT #####
     3#######################
     4
     5### A. Spiga     -- LMD -- 06~09/2011 -- General building and mapping capabilities
     6### T. Navarro   -- LMD -- 10~11/2011 -- Improved use for GCM and added sections + 1Dplot capabilities
     7
    118def planetoplot (namefiles,\
    129           vertmode,\
    … …  
    4946                       fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\
    5047                       zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow,\
    51                        getname,localtime,polarinterv,getsindex,define_axis
     48                       getname,localtime,polarinterv,getsindex,define_axis,determineplot
    5249    from mymath import deg,max,min,mean
    5350    from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel
    … …  
    5653    from numpy.core.defchararray import find
    5754       
    58     #whichmode = ""
    59     nlon = 1 # number of longitudinal slices -- 1 is None
    60     nlat = 1
    61     nvert = 1
    62     ntime = 1
    63     nslices = 1
    64     if slon is not None:
    65         nslices = nslices*len(slon)
    66         nlon = len(slon)
    67     if slat is not None:
    68         nslices = nslices*len(slat)
    69         nlat = len(slat)
    70     if svert is not None:
    71         nslices = nslices*len(svert)
    72         nvert = len(svert)
    73     if stime is not None:
    74         nslices = nslices*len(stime)
    75         ntime = len(stime)
    76     #else:
    77     #    nslices = 2
    78     numplot = len(namefiles)*nslices
     55    ################################
     56    ### Which plot needs to be done?
     57    nlon, nlat, nvert, ntime, mapmode, nslices = determineplot(slon, slat, svert, stime)
     58    if mapmode == 0: winds=False
     59    if not isinstance(namefiles, np.ndarray): namefiles = [namefiles]
     60    zelen = len(namefiles)
     61    numplot = zelen*nslices
     62    print "len(namefiles), nslices, numplot: ", zelen, nslices, numplot
     63    all_var   = [[]]*zelen
     64    all_var2  = [[]]*zelen
     65    all_title = [[]]*zelen
    7966   
    80     mapmode = 0
    81     if slon is None and slat is None:
    82        mapmode = 1 # in this case we plot a map, with the given projection
    83     elif proj is not None:
    84        print "WARNING: you specified a", proj,\
    85        "projection but asked for slices", slon,"in longitude and",slat,"in latitude"
    86     print "mapmode", mapmode
    87        
    88    
    89     all_var   = [[]]*len(namefiles)
    90     all_var2  = [[]]*len(namefiles)
    91     all_title = [[]]*len(namefiles)
    92    
    93     print "len(namefiles), nslices", len(namefiles), nslices
    94     print "numplot", numplot
    95 
    9667    #########################
    9768    ### Loop over the files initially separated by comas to be plotted on the same figure
    … …  
    10677      ##################################
    10778      ### Initial checks and definitions
    108       typefile = whatkindfile(nc)                                  ## TYPEFILE
     79      ### ... TYPEFILE
     80      typefile = whatkindfile(nc)                                 
    10981      if firstfile:
    11082         typefile0 = typefile
    11183      elif typefile != typefile0:
    112          errormess("Not the same files !", [typefile0, typefile])
    113       if var not in nc.variables: var = False                      ## VAR
    114       if winds:                                                    ## WINDS
     84         errormess("Not the same kind of files !", [typefile0, typefile])
     85      ### ... VAR
     86      if var not in nc.variables: var = False
     87      ### ... WINDS
     88      if winds:                                                   
    11589         [uchar,vchar,metwind] = getwinddef(nc)             
    11690         if uchar == 'not found': winds = False
    11791      if not var and not winds: errormess("please set at least winds or var",printvar=nc.variables)
    118       [lon2d,lat2d] = getcoorddef(nc)                              ## COORDINATES, could be moved below
    119       if proj == None:   proj = getproj(nc)                        ## PROJECTION
    120 
    121       lat = nc.variables["latitude"][:]
    122       lon = nc.variables["longitude"][:]
    123       time = nc.variables["Time"][:]
    124       vert = nc.variables["altitude"][:]
    125       if firstfile:
    126          lat0 = lat
    127       elif len(lat0) != len(lat):
    128          errormess("Not the same latitude lengths !", [len(lat0), len(lat)])
    129       elif sum((lat == lat0) == False) != 0:
    130          errormess("Not the same latitudes !", [lat,lat0])
    131       # Faire d'autre checks sur les compatibilites entre fichiers!!
     92      ### ... COORDINATES, could be moved below
     93      [lon2d,lat2d] = getcoorddef(nc)                       
     94      ### ... PROJECTION
     95      if proj == None:   proj = getproj(nc)                 
     96
     97##########################################################
     98      if typefile == "gcm":
     99          lat = nc.variables["latitude"][:]
     100          lon = nc.variables["longitude"][:]
     101          time = nc.variables["Time"][:]
     102          vert = nc.variables["altitude"][:]
     103          #if firstfile:
     104          #   lat0 = lat
     105          #elif len(lat0) != len(lat):
     106          #   errormess("Not the same latitude lengths !", [len(lat0), len(lat)])
     107          #elif sum((lat == lat0) == False) != 0:
     108          #   errormess("Not the same latitudes !", [lat,lat0])
     109          ## Faire d'autre checks sur les compatibilites entre fichiers!!
     110##########################################################
    132111
    133112      if firstfile:
    … …  
    143122         ### Name for title and graphics save file
    144123         basename = getname(var=var,winds=winds,anomaly=anomaly)
    145          basename = basename + getstralt(nc,nvert)  ## can be moved elsewhere for a more generic routine
     124         basename = basename + getstralt(nc,vertmode)  ## can be moved elsewhere for a more generic routine
    146125     
    147       print "var", var
    148       print "var2", var2
    149       #print all_var
    150       #print nc
     126      print "var, var2: ", var, var2
    151127      if var: all_var[k] = getfield(nc,var)
    152128      if var2: all_var2[k] = getfield(nc,var2)
    … …  
    173149    elif numplot <= 0:                 itstep = 1
    174150   
    175     ### Map projection
    176     if mapmode == 1:
    177         m = define_proj(proj,wlon,wlat,back=back)
    178         x, y = m(lon2d, lat2d)
    179        
    180151    #for nplot in range(numplot):
    181152    while error is False:
    … …  
    201172           else:
    202173                 found_lct = True
    203                  
    204                  
    205        ## get all indexes to be taken into account for this subplot and then reduce field
    206        ## We plot 1) all lon slices 2) all lat slices 3) all vert slices 4) all time slices and then go to the next slice
    207        #print "(nplot-1)%nlon", (nplot-1)%nlon
    208        indexlon  = getsindex(slon,(nplot-1)%nlon,lon)
    209        #print "indexlon", indexlon
    210        indexlat  = getsindex(slat,((nplot-1)//nlon)%nlat,lat)
    211        indexvert = getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert)
    212        indextime = getsindex(stime,((nplot-1)//(nlon*nlat*nvert))%ntime,time)
    213        index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%len(namefiles)
    214        #print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime)
     174
     175       ### Map projection                   
     176       if mapmode == 1:
     177           m = define_proj(proj,wlon,wlat,back=back)
     178           x, y = m(lon2d, lat2d)
     179
     180####################################################################
     181       if typefile in ['mesoapi','meso']:
     182           indexlon = None
     183           indexlat = None
     184           indexvert = vertmode 
     185           indextime = itime
     186           nlon = 1
     187           nlat = 1
     188           nvert = 1
     189           ntime = 1 
     190           index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%len(namefiles)
     191           print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime)
     192       elif typefile in ['gcm']:
     193           ## get all indexes to be taken into account for this subplot and then reduce field
     194           ## We plot 1) all lon slices 2) all lat slices 3) all vert slices 4) all time slices and then go to the next slice
     195           #print "(nplot-1)%nlon", (nplot-1)%nlon
     196           indexlon  = getsindex(slon,(nplot-1)%nlon,lon)
     197           print "indexlon", indexlon
     198           indexlat  = getsindex(slat,((nplot-1)//nlon)%nlat,lat)
     199           indexvert = getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert) 
     200           indextime = getsindex(stime,((nplot-1)//(nlon*nlat*nvert))%ntime,time) 
     201           index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%len(namefiles)
     202           #print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime)
     203####################################################################
    215204
    216205       #### Contour plot
    217206       if var2:
    218207           what_I_plot, error = reducefield(all_var2[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert )
    219            what_I_plot = what_I_plot*mult
     208           #what_I_plot = what_I_plot*mult
    220209           if not error:
    221               if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_contour,6)
     210              if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot,6)
    222211              zevmin, zevmax = calculate_bounds(what_I_plot)
    223               zelevels = np.linspace(zevmin,zevmax,num=ndiv)
     212              zelevels = np.linspace(zevmin,zevmax,num=ndiv+1) #20)
     213              if var2 == 'HGT':  zelevels = np.arange(-10000.,30000.,2000.)
    224214              if mapmode == 0:
    225215                  x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
    … …  
    227217              ### If we plot a 2-D field
    228218              if len(what_I_plot.shape) is 2:
    229                   cs = contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #colors='w' )# , alpha=0.5)
    230                   clabel(cs,fmt = '%1.2e')
     219                  cs = contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #0.33 colors='w' )# , alpha=0.5)
     220                  if typefile in ['gcm']: clabel(cs,fmt = '%1.2e')
    231221              ### If we plot a 1-D field
    232222              elif len(what_I_plot.shape) is 1:
    … …  
    256246               if len(what_I_plot.shape) is 2:
    257247                 if mapmode == 0:
    258                      what_I_plot,x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
     248                     what_I_plot, x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
    259249                            indextime,what_I_plot, len(all_var[index_f].shape),vertmode)
    260250                 if not tile:
    … …  
    263253                     print "what_I_plot.shape", what_I_plot.shape
    264254                     print "x.shape, y.shape", x.shape, y.shape
    265                      zelevels = np.linspace(zevmin,zevmax,num=ndiv)
     255                     zelevels = np.linspace(zevmin,zevmax,num=ndiv+1)
    266256                     #contourf(what_I_plot, zelevels, cmap = palette )
    267257                     contourf(x,y,what_I_plot, zelevels, cmap = palette )
    … …  
    274264                         #colorbar()         
    275265                         colorbar(fraction=0.05,pad=0.03,format=fmtvar(fvar.upper()),\
    276                                            ticks=np.linspace(zevmin,zevmax,min([ndiv,20])),\
     266                                           ticks=np.linspace(zevmin,zevmax,min([ndiv+1,20])),\
    277267                                           extend='neither',spacing='proportional')
    278268                                           # both min max neither
    … …  
    295285           else:
    296286               errormess("There is an error in reducing field !")
    297  
     287
     288       ### Vector plot --- a adapter
     289       if winds:
     290           vecx, error = reducefield( getfield(nc,uchar), d4=indextime, d3=indexvert )
     291           vecy, error = reducefield( getfield(nc,vchar), d4=indextime, d3=indexvert )
     292           #what_I_plot, error = reducefield(all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert )
     293           if not error:
     294               if typefile in ['mesoapi','meso']:   
     295                   [vecx,vecy] = [dumpbdy(vecx,6,stag=uchar), dumpbdy(vecy,6,stag=vchar)]
     296                   key = True
     297               elif typefile in ['gcm']:           
     298                   key = False
     299               if metwind:  [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d)
     300               if var == False:       colorvec = definecolorvec(back)
     301               else:                  colorvec = definecolorvec(colorb)
     302               vectorfield(vecx, vecy,\
     303                          x, y, stride=stride, csmooth=2,\
     304                          #scale=15., factor=300., color=colorvec, key=key)
     305                          scale=20., factor=250., color=colorvec, key=key)
     306                                            #200.         ## or csmooth=stride
    298307   
    299308       ### Next subplot
    300        plottitle = basename+' '+namefiles[index_f]
    301309       if typefile in ['mesoapi','meso']:
     310            plottitle = basename
    302311            if addchar:  plottitle = plottitle + addchar + "_LT"+str(ltst)
    303312            else:        plottitle = plottitle + "_LT"+str(ltst)
     313       else:
     314            plottitle = basename+' '+namefiles[index_f]
    304315       if mult != 1:           plottitle = str(mult) + "*" + plottitle
    305316       if zetitle != "fill":   plottitle = zetitle
    … …  
    314325       ptitle( plottitle )
    315326       itime += itstep
    316        if nplot == numplot:
     327       if nplot >= numplot:
    317328          error = True
    318329       nplot += 1
    … …  
    356367    ### Now the end
    357368    return zeplot
    358 
    359 ##############################
    360 ### A specific stuff for below
    361 def adjust_length (tab, zelen):
    362     from numpy import ones
    363     if tab is None:
    364         outtab = ones(zelen) * -999999
    365     else:
    366         if zelen != len(tab):
    367             print "not enough or too much values... setting same values all variables"
    368             outtab = ones(zelen) * tab[0]
    369         else:
    370             outtab = tab
    371     return outtab
    372 
    373 ###########################################################################################
    374 ###########################################################################################
    375 ### What is below relate to running the file as a command line executable (very convenient)
    376 if __name__ == "__main__":
    377     import sys
    378     from optparse import OptionParser    ### to be replaced by argparse
    379     from api_wrapper import api_onelevel
    380     from netCDF4 import Dataset
    381     from myplot import getlschar, separatenames, readslices
    382     from os import system
    383     import numpy as np
    384 
    385     #############################
    386     ### Get options and variables
    387     parser = OptionParser()
    388     parser.add_option('-f', '--file',   action='append',dest='namefile', type="string",  default=None,  help='[NEEDED] name of WRF file (append). Plot files separated by comas in the same figure')
    389     parser.add_option('-l', '--level',  action='store',dest='nvert',     type="float",   default=0,     help='level (def=0)(-i 2: p,mbar)(-i 3,4: z,km)')
    390     parser.add_option('-p', '--proj',   action='store',dest='proj',      type="string",  default=None,  help='projection')
    391     parser.add_option('-b', '--back',   action='store',dest='back',      type="string",  default=None,  help='background image (def: None)')
    392     parser.add_option('-t', '--target', action='store',dest='target',    type="string",  default=None,  help='destination folder')
    393     parser.add_option('-s', '--stride', action='store',dest='stride',    type="int",     default=3,     help='stride vectors (def=3)')
    394     parser.add_option('-v', '--var',    action='append',dest='var',      type="string",  default=None,  help='variable color-shaded (append)')
    395     parser.add_option('-n', '--num',    action='store',dest='numplot',   type="int",     default=2,     help='plot number (def=2)(<0: plot LT -*numplot*)')
    396     parser.add_option('-i', '--interp', action='store',dest='interp',    type="int",     default=None,  help='interpolation (2: p, 3: z-amr, 4:z-als)')
    397     parser.add_option('-c', '--color',  action='store',dest='colorb',    type="string",  default="def", help='change colormap (nobar: no colorbar)')
    398     parser.add_option('-x', '--no-vect',action='store_false',dest='winds',               default=True,  help='no wind vectors')
    399     parser.add_option('-m', '--min',    action='append',dest='vmin',     type="float",   default=None,  help='bounding minimum value (append)')   
    400     parser.add_option('-M', '--max',    action='append',dest='vmax',     type="float",   default=None,  help='bounding maximum value (append)')
    401     parser.add_option('-T', '--tiled',  action='store_true',dest='tile',                 default=False, help='draw a tiled plot (no blank zone)')
    402     parser.add_option('-z', '--zoom',   action='store',dest='zoom',      type="float",   default=None,  help='zoom factor in %')
    403     parser.add_option('-N', '--no-api', action='store_true',dest='nocall',               default=False, help='do not recreate api file')
    404     parser.add_option('-d', '--display',action='store_false',dest='display',             default=True,  help='do not pop up created images')
    405     parser.add_option('-e', '--itstep', action='store',dest='itstep',    type="int",     default=None,  help='stride time (def=4)')
    406     parser.add_option('-H', '--hole',   action='store_true',dest='hole',                 default=False, help='holes above max and below min')
    407     parser.add_option('-S', '--save',   action='store',dest='save',      type="string",  default="gui", help='save mode (png,eps,svg,pdf or gui)(def=gui)')
    408     parser.add_option('-a', '--anomaly',action='store_true',dest='anomaly',              default=False, help='compute and plot relative anomaly in %')
    409     parser.add_option('-w', '--with',   action='store',dest='var2',      type="string",  default=None,  help='variable contoured')
    410     parser.add_option('--div',          action='store',dest='ndiv',      type="int",     default=10,    help='number of divisions in colorbar (def: 10)')
    411     parser.add_option('-F', '--first',  action='store',dest='first',     type="int",     default=1,     help='first subscript to plot (def: 1)')
    412     parser.add_option('--mult',         action='store',dest='mult',      type="float",   default=1.,    help='a multiplicative factor to plotted field')
    413     parser.add_option('--title',        action='store',dest='zetitle',   type="string",  default="fill",help='customize the whole title')
    414     #parser.add_option('-V', action='store', dest='comb',        type="float",   default=None,  help='a defined combination of variables')
    415 
    416     ############# T.N. changes
    417     #parser.add_option('-o','--operation',action='store',dest='operation',type="string",   default=None ,help='matrix of operations between files (for now see code, sorry)')
    418     parser.add_option('--lat',          action='append',dest='slat',type="string",   default=None, help='slices along latitude. One value, or two values separated by comas for averaging')
    419     parser.add_option('--lon',          action='append',dest='slon', type="string",   default=None, help='slices along longitude. One value, or two values separated by comas for averaging')
    420     parser.add_option('--vert',         action='append',dest='svert',type="string",   default=None, help='slices along vertical axis. One value, or two values separated by comas for averaging') # quelles coordonnees ?
    421     parser.add_option('--time',         action='append',dest='stime',type="string",   default=None, help='slices along time. One value, or two values separated by comas for averaging') # quelles coordonnees ?
    422 
    423     (opt,args) = parser.parse_args()
    424     if opt.namefile is None:
    425         print "I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h"
    426         exit()
    427     if opt.var is None and opt.anomaly is True:
    428         print "Cannot ask to compute anomaly if no variable is set"
    429         exit()
    430     print "Options:", opt
    431    
    432     #listvar = ''
    433     #if opt.var is None:
    434     #    zerange = [-999999]
    435     #else:
    436     #    zelen = len(opt.var)
    437     #    zerange = range(zelen)
    438     #    #if zelen == 1: listvar = opt.var[0] + ','
    439     #    #else         :
    440     #    for jjj in zerange: listvar += opt.var[jjj] + ','
    441     #    listvar = listvar[0:len(listvar)-1]
    442     #    vmintab = adjust_length (opt.vmin, zelen) 
    443     #    vmaxtab = adjust_length (opt.vmax, zelen)
    444        
    445     print "namefile, length", opt.namefile, len(opt.namefile)
    446 
    447     zeslat  = readslices(opt.slat)
    448     zeslon  = readslices(opt.slon)
    449     zesvert = readslices(opt.svert)
    450     zestime = readslices(opt.stime)
    451     print "slat,zeslat", opt.slat, zeslat
    452     print "slon,zeslon", opt.slon, zeslon
    453     print "svert,zesvert", opt.svert, zesvert
    454     print "stime,zestime", opt.stime, zestime
    455 
    456      
    457     for i in range(len(opt.namefile)):
    458       for j in range(len(opt.var)):
    459 
    460         zenamefiles = separatenames(opt.namefile[i])
    461         print "zenamefiles", zenamefiles
    462        
    463         if opt.vmin is not None : zevmin  = opt.vmin[min(i,len(opt.vmin)-1)]
    464         else: zevmin = None
    465         if opt.vmax is not None : zevmax  = opt.vmax[min(i,len(opt.vmax)-1)]
    466         else: zevmax = None
    467         print "vmin, zevmin", opt.vmin, zevmin
    468         print "vmax, zevmax", opt.vmax, zevmax
    469        
    470         zevar = separatenames(opt.var[j])
    471         zevar = zevar[0]
    472         print "var, zevar", opt.var, zevar
    473        
    474         #checkcoherence(len(zenamefiles),len(opt.slat),len(opt.slon),len(opt.stime))
    475        
    476         zefile = zenamefiles[0]
    477              
    478         zelevel = opt.nvert   
    479         stralt = None
    480         [lschar,zehour,zehourin] = getlschar ( zefile )  ## getlschar from wrfout (or simply return "" if another file)
    481         #print "lschar ",lschar
    482         #print "zehour ",zehour
    483         #print "zehourin ",zehourin
    484    
    485         #####################################################
    486         ### Call Fortran routines for vertical interpolations
    487         if opt.interp is not None:
    488             if opt.nvert is 0 and opt.interp is 4:  zelevel = 0.010
    489             ### winds or no winds
    490             if opt.winds            :  zefields = 'uvmet'
    491             else                    :  zefields = ''
    492             ### var or no var
    493             #if opt.var is None      :  pass
    494             if zefields == ''       :  zefields = listvar
    495             else                    :  zefields = zefields + "," + listvar
    496             if opt.var2 is not None : zefields = zefields + "," + opt.var2 
    497             print zefields
    498             zefile = api_onelevel (  path_to_input   = '', \
    499                                      input_name      = zefile, \
    500                                      fields          = zefields, \
    501                                      interp_method   = opt.interp, \
    502                                      onelevel        = zelevel, \
    503                                      nocall          = opt.nocall )
    504             print zefile
    505             zelevel = 0 ## so that zelevel could play again the role of nvert
    506 
    507 
    508         #############
    509         ### Main call
    510         name = planetoplot (zenamefiles,int(zelevel),\
    511                 proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=zevar,\
    512                 numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\
    513                 addchar=lschar,interv=[zehour,zehourin],vmin=zevmin,vmax=zevmax,\
    514                 tile=opt.tile,zoom=opt.zoom,display=opt.display,\
    515                 itstep=opt.itstep,hole=opt.hole,save=opt.save,\
    516                 anomaly=opt.anomaly,var2=opt.var2,ndiv=opt.ndiv,first=opt.first,\
    517                 mult=opt.mult,zetitle=opt.zetitle,\
    518                 slon=zeslon,slat=zeslat,svert=zesvert,stime=zestime)
    519         print 'Done: '+name
    520         system("rm -f to_be_erased")
    521  
    522     #########################################################
    523     ### Generate a .sh file with the used command saved in it
    524     command = ""
    525     for arg in sys.argv: command = command + arg + ' '
    526     name = 'pycommand'
    527     f = open(name+'.sh', 'w')
    528     f.write(command)
Note: See TracChangeset for help on using the changeset viewer.