source: trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/winds.py @ 230

Last change on this file since 230 was 225, checked in by aslmd, 14 years ago

MESOSCALE: python graphics. now possible for multi-files and multi-var. very easy to make an atlas.

  • Property svn:executable set to *
File size: 19.2 KB
Line 
1#!/usr/bin/env python
2
3### A. Spiga -- LMD -- 30/06/2011 to 10/07/2011
4### Thanks to A. Colaitis for the parser trick
5
6
7####################################
8####################################
9### The main program to plot vectors
10def winds (namefile,\
11           nvert,\
12           proj=None,\
13           back=None,\
14           target=None,
15           stride=3,\
16           numplot=2,\
17           var=None,\
18           colorb=True,\
19           winds=True,\
20           addchar=None,\
21           interv=[0,1],\
22           vmin=None,\
23           vmax=None,\
24           tile=False,\
25           zoom=None):
26
27    ####################################################################################################################
28    ### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
29
30    #################################
31    ### Load librairies and functions
32    from netCDF4 import Dataset
33    from myplot import getcoord2d,define_proj,makeplotpng,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy,\
34                       fmtvar,definecolorvec,getwinds,defcolorb,getprefix
35    from mymath import deg,max,min,mean
36    from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor
37    from matplotlib.cm import get_cmap
38    import numpy as np
39    from numpy.core.defchararray import find
40
41    ###
42    #rcParams['text.usetex'] = True
43    #rcParams['cairo.format'] = 'svg'
44
45    ######################
46    ### Load NETCDF object
47    nc  = Dataset(namefile)
48
49    ###################################
50    ### Recognize predefined file types
51    if 'controle' in nc.variables:   typefile = 'gcm'
52    elif 'vert' in nc.variables:     typefile = 'mesoapi'
53    elif 'U' in nc.variables:        typefile = 'meso'
54    elif 'HGT_M' in nc.variables:    typefile = 'geo'
55    else:                           
56        print "typefile not supported."
57        print nc.variables
58        exit()
59
60    ##############################################################
61    ### Try to guess the projection from wrfout if not set by user
62    if typefile in ['mesoapi','meso','geo']:
63        if proj == None:       proj = getproj(nc)
64                                    ### (il faudrait passer CEN_LON dans la projection ?)
65    elif typefile in ['gcm']:
66        if proj == None:       proj = "cyl"   
67                                    ## pb avec les autres (de trace derriere la sphere ?)
68
69    ############################################
70    #### Choose underlying topography by default
71    if not back:
72        if not var:                                        back = "mola"    ## if no var:         draw mola
73        elif typefile in ['mesoapi','meso','geo'] \
74           and proj not in ['merc','lcc','nsper','laea']:  back = "molabw"  ## if var but meso:   draw molabw
75        else:                                              pass             ## else:              draw None
76
77    ####################################################
78    ### Get geographical coordinates and plot boundaries
79    if typefile in ['mesoapi','meso']:
80        [lon2d,lat2d] = getcoord2d(nc)
81        lon2d = dumpbdy(lon2d)
82        lat2d = dumpbdy(lat2d)
83    elif typefile in ['gcm']:
84        [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True)
85    elif typefile in ['geo']:
86        [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M')
87    if proj == "npstere":             [wlon,wlat] = latinterv("North_Pole")
88    elif proj in ["lcc","laea"]:      [wlon,wlat] = wrfinterv(lon2d,lat2d)
89    else:                             [wlon,wlat] = simplinterv(lon2d,lat2d)
90    print wlon
91    print wlat
92    if zoom: 
93        dlon = abs(wlon[1]-wlon[0])/2.
94        dlat = abs(wlat[1]-wlat[0])/2.
95        [wlon,wlat] = [ [wlon[0]+zoom*dlon/100.,wlon[1]-zoom*dlon/100.],\
96                        [wlat[0]+zoom*dlat/100.,wlat[1]-zoom*dlat/100.] ]
97        print "zoom %",zoom,wlon,wlat
98
99    ##############################################################################
100    ### Get winds and know if those are meteorological winds (ie. zon, mer) or not
101    if winds:
102        if typefile is 'mesoapi': 
103            [u,v] = getwinds(nc)
104            metwind = True  ## meteorological (zon/mer)
105        elif typefile is 'gcm':
106            [u,v] = getwinds(nc,charu='u',charv='v')
107            metwind = True  ## meteorological (zon/mer)
108        elif typefile is 'meso':
109            [u,v] = getwinds(nc,charu='U',charv='V')
110            metwind = False ## geometrical (wrt grid)
111            print "Beware ! Not using meteorological winds. You trust numerical grid as being (x,y)."
112        elif typefile is 'geo':
113            winds = None
114
115    #####################################################
116    ### Load the chosen variables, whether it is 2D or 3D
117    if var:
118        if var not in nc.variables: 
119            print "not found in file:",var
120            exit()
121        else:   
122            dimension = np.array(nc.variables[var]).ndim
123            if dimension == 2:     field = nc.variables[var][:,:]
124            elif dimension == 3:   field = nc.variables[var][:,:,:]
125            elif dimension == 4:   field = nc.variables[var][:,nvert,:,:] 
126        dev = np.std(field)*3.0
127        if vmin is None:  zevmin = mean(field) - dev
128        else:             zevmin = vmin
129        if vmax is None:  zevmax = mean(field) + dev
130        else:             zevmax = vmax
131        if vmin == vmax:   
132                          zevmin = mean(field) - dev  ### for continuity
133                          zevmax = mean(field) + dev  ### for continuity           
134        print "field ", min(field), max(field)
135        print "bounds ", zevmin, zevmax
136        ### some already defined colormaps
137        if colorb is True:    colorb = defcolorb(var)
138    else:
139        dimension = 0
140 
141    ###########################
142    ### Get length of time axis
143    if winds:               nt = len(u[:,0,0,0])
144    elif var: 
145        if dimension == 2:  nt = 1
146        else             :  nt = len(field[:,0,0])
147
148    #########################################
149    ### Name for title and graphics save file
150    if var and winds:     basename = var + '_UV'
151    elif var:             basename = var
152    elif winds:           basename = 'UV'
153    else:                 exit()
154    ###
155    if dimension == 4 or winds:
156        if typefile is 'meso':                      stralt = "_lvl" + str(nvert)
157        elif typefile is 'mesoapi': 
158            zelevel = int(nc.variables['vert'][nvert])
159            if 'altitude'       in nc.dimensions:   stralt = "_"+str(zelevel)+"m-AMR"
160            elif 'altitude_abg' in nc.dimensions:   stralt = "_"+str(zelevel)+"m-ALS"
161            elif 'bottom_top'   in nc.dimensions:   stralt = "_"+str(zelevel)+"m"
162            elif 'pressure'     in nc.dimensions:   stralt = "_"+str(zelevel)+"Pa" 
163        else:                                       stralt = ""         
164    else:
165        stralt = ""
166    ###
167    basename = basename + stralt
168
169    ##################################
170    ### Open a figure and set subplots
171    fig = figure()
172    rcParams['font.size'] = 12.
173    if   typefile in ['geo']:  numplot = 1
174    if   numplot > 0:   
175        if   numplot == 4: 
176            sub = 221
177            fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
178            rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
179        elif numplot == 2: 
180            sub = 121
181            fig.subplots_adjust(wspace = 0.3)
182            rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
183        elif numplot == 3: 
184            sub = 131
185            fig.subplots_adjust(wspace = 0.5)
186            rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
187        elif numplot == 6: 
188            sub = 231
189            fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
190            rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
191        elif numplot == 8: 
192            sub = 331 #241
193            fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
194            rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
195        elif numplot == 9:
196            sub = 331
197            fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
198            rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
199        elif numplot == 1:
200            sub = 99999
201        else:
202            print "supported: 1,2,3,4,6,8,9"
203            exit()
204        ### Prepare time loop
205        if nt <= numplot or numplot == 1: 
206            tabrange = [0]
207            numplot = 1
208        else:                         
209            tabrange = range(0,nt,int(nt/numplot))  #nt-1
210            tabrange = tabrange[0:numplot]
211    else: 
212        tabrange = range(0,nt,1)
213        sub = 99999
214    print tabrange
215
216    #################################
217    ### Time loop for plotting device
218    found_lct = False
219    for i in tabrange:
220
221       ### Which local time ?
222       ltst = ( interv[0] + 0.5*(wlon[0]+wlon[1])/15.) + i*interv[1]
223       ltst = int (ltst * 10) / 10.
224       ltst = ltst % 24
225
226       ### General plot settings
227       if numplot > 1: 
228           subplot(sub)
229           found_lct = True
230       elif numplot == 1:
231           found_lct = True 
232            ### If only one local time is requested (numplot < 0)
233       elif numplot <= 0: 
234           if int(ltst) + numplot != 0:         continue
235           else:                                found_lct = True
236
237       ### Map projection
238       m = define_proj(proj,wlon,wlat,back=back)
239       x, y = m(lon2d, lat2d)
240
241       #### Contour plot
242       if var:
243           if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(field[i,:,:])
244           elif typefile in ['geo']:             what_I_plot = field[0,:,:]
245           elif typefile in ['gcm']:             
246               if dimension == 2:                what_I_plot = field[:,:]
247               elif dimension == 3:              what_I_plot = field[i,:,:]
248           palette = get_cmap(name=colorb)
249           #palette.set_over('b', 1.0)
250           #print np.array(x).shape
251           #print np.array(y).shape
252           #print np.array(what_I_plot).shape
253           if not tile:
254               zelevels = np.linspace(zevmin,zevmax)
255               what_I_plot[ what_I_plot < zevmin ] = zevmin*(1. + 1.e-7)
256               what_I_plot[ what_I_plot > zevmax ] = zevmax*(1. - 1.e-7)
257               contourf( x, y, what_I_plot, 10, cmap = palette, levels = zelevels )
258           else:   
259               pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax )
260           if var in ['HGT']: pass
261           elif colorb:             
262                              ndiv = 10
263                              colorbar(fraction=0.05,pad=0.1,format=fmtvar(var),\
264                                       ticks=np.linspace(zevmin,zevmax,ndiv+1),\
265                                       extend='max',spacing='proportional')
266                                       # both min max neither
267       ### Vector plot
268       if winds:
269           if   typefile in ['mesoapi','meso']:   
270               [vecx,vecy] = [dumpbdy(u[i,nvert,:,:]), dumpbdy(v[i,nvert,:,:])]
271               key = True
272           elif typefile in ['gcm']:               
273               [vecx,vecy] = [        u[i,nvert,:,:] ,         v[i,nvert,:,:] ]
274               key = False
275           if metwind:  [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d)
276           if var == None:        colorvec = definecolorvec(back)
277           else:                  colorvec = definecolorvec(colorb)
278           vectorfield(vecx, vecy,\
279                      x, y, stride=stride, csmooth=2,\
280                      scale=15., factor=300., color=colorvec, key=key)
281                                        #200.         ## or csmooth=stride
282       
283       ### Next subplot
284       plottitle = basename
285       if typefile in ['mesoapi','meso']:
286           if addchar:  plottitle = plottitle + addchar + "_LT"+str(ltst)
287           else:        plottitle = plottitle + "_LT"+str(ltst)
288       ptitle( plottitle )
289       sub += 1
290
291    ##########################################################################
292    ### Save the figure in a file in the data folder or an user-defined folder
293    if typefile in ['meso','mesoapi']:   prefix = getprefix(nc)   
294    elif typefile in ['gcm']:            prefix = 'LMD_GCM_'
295    else:                                prefix = ''
296    ###
297    zeplot = prefix + basename
298    if addchar:         zeplot = zeplot + addchar
299    if numplot <= 0:    zeplot = zeplot + "_LT"+str(abs(numplot))
300    ###
301    if not target:      zeplot = namefile[0:find(namefile,'wrfout')] + zeplot
302    else:               zeplot = target + "/" + zeplot 
303    ###
304    if found_lct:     makeplotpng(zeplot,pad_inches_value=0.35)   
305    else:             print "Local time not found"
306
307    ###############
308    ### Now the end
309    return zeplot
310
311###########################################################################################
312###########################################################################################
313### What is below relate to running the file as a command line executable (very convenient)
314if __name__ == "__main__":
315    import sys
316    from optparse import OptionParser    ### to be replaced by argparse
317    from api_wrapper import api_onelevel
318    from netCDF4 import Dataset
319    from myplot import getlschar
320    from numpy import ones
321
322    #############################
323    ### Get options and variables
324    parser = OptionParser()
325    parser.add_option('-f', action='append', dest='namefile',    type="string",  default=None,  help='[NEEDED] name of WRF file (append)')
326    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)')
327    parser.add_option('-p', action='store', dest='proj',        type="string",  default=None,  help='projection')
328    parser.add_option('-b', action='store', dest='back',        type="string",  default=None,  help='background')
329    parser.add_option('-t', action='store', dest='target',      type="string",  default=None,  help='destination folder')
330    parser.add_option('-s', action='store', dest='stride',      type="int",     default=3,     help='stride vectors (def=3)')
331    parser.add_option('-v', action='append', dest='var',         type="string",  default=None,  help='variable contoured (append)')
332    parser.add_option('-n', action='store', dest='numplot',     type="int",     default=2,     help='number of plots (def=1)(<0: 1 plot of LT -*numplot*)')
333    parser.add_option('-i', action='store', dest='interp',      type="int",     default=None,  help='interpolation method (2: press, 3: z-amr, 4:z-als)')
334    parser.add_option('-c', action='store', dest='colorb',      type="string",  default=True,  help='change colormap')
335    parser.add_option('-x', action='store_false', dest='winds',                 default=True,  help='no wind vectors')
336    parser.add_option('-m', action='append', dest='vmin',        type="float",   default=None,  help='bounding minimum value (append)')   
337    parser.add_option('-M', action='append', dest='vmax',        type="float",   default=None,  help='bounding maximum value (append)') 
338    parser.add_option('-T', action='store_true', dest='tile',                   default=False, help='draw a tiled plot (no blank zone)')
339    parser.add_option('-z', action='store', dest='zoom',        type="float",   default=None,  help='zoom factor in %')
340    parser.add_option('-N', action='store_true', dest='nocall',                 default=False, help='do not recreate api file')
341    #parser.add_option('-V', action='store', dest='comb',        type="float",   default=None,  help='a defined combination of variables')
342    (opt,args) = parser.parse_args()
343    if opt.namefile is None: 
344        print "I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h"
345        exit()
346    print "Options:", opt
347   
348    if opt.var is None:
349        pass
350    else:
351        listvar = ''
352        zelen = len(opt.var)
353        if zelen == 1: listvar = opt.var[0] + ','
354        else         : 
355            for jjj in range(zelen): listvar += opt.var[jjj] + ','
356        listvar = listvar[0:len(listvar)-1]
357        if opt.vmin:
358            if zelen != len(opt.vmin):
359                print "not enough or too much vmin values... setting same values all variables"
360                vmintab = ones(zelen) * opt.vmin[0]
361            else:
362                vmintab = opt.vmin
363        else:
364            vmintab = None
365        if opt.vmax:
366            if zelen != len(opt.vmax):
367                print "not enough or too much vmax values... setting same values all variables"
368                vmaxtab = ones(zelen) * opt.vmax[0]
369            else:
370                vmaxtab = opt.vmax
371        else:
372            vmaxtab = None
373
374    for i in range(len(opt.namefile)):
375
376        zefile = opt.namefile[i]
377        print zefile   
378        zelevel = opt.nvert   
379        stralt = None
380        [lschar,zehour,zehourin] = getlschar ( zefile )  ## getlschar from wrfout (or simply return "" if another file)
381   
382        #####################################################
383        ### Call Fortran routines for vertical interpolations
384        if opt.interp is not None:
385            if opt.nvert is 0 and opt.interp is 4:  zelevel = 0.010
386            ### winds or no winds
387            if opt.winds            :  zefields = 'uvmet'
388            else                    :  zefields = ''
389            ### var or no var
390            if opt.var is None      :  pass
391            elif zefields == ''     :  zefields = listvar
392            else                    :  zefields = zefields + "," + listvar
393            print zefields
394            zefile = api_onelevel (  path_to_input   = '', \
395                                     input_name      = zefile, \
396                                     #path_to_output  = opt.target, \
397                                     fields          = zefields, \
398                                     interp_method   = opt.interp, \
399                                     onelevel        = zelevel, \
400                                     nocall          = opt.nocall )
401            print zefile
402            zelevel = 0 ## so that zelevel could play again the role of nvert
403
404        if opt.var is None: 
405            #############
406            ### Main call
407            name = winds (zefile,int(zelevel),\
408               proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var,numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\
409               addchar=lschar,interv=[zehour,zehourin],vmin=opt.vmin,vmax=opt.vmax,tile=opt.tile,zoom=opt.zoom)
410            print 'Done: '+name
411        else:       
412            for jjj in range(len(opt.var)): 
413                if vmintab: argvmin = vmintab[jjj]
414                else:       argvmin = None
415                if vmaxtab: argvmax = vmaxtab[jjj] 
416                else:       argvmax = None
417                #############
418                ### Main call
419                name = winds (zefile,int(zelevel),\
420                   proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var[jjj],numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\
421                   addchar=lschar,interv=[zehour,zehourin],vmin=argvmin,vmax=argvmax,tile=opt.tile,zoom=opt.zoom)
422                print 'Done: '+name
423   
424        #########################################################
425        ### Generate a .sh file with the used command saved in it
426        command = ""
427        for arg in sys.argv: command = command + arg + ' '
428        f = open(name+'.sh', 'w')
429        f.write(command)
Note: See TracBrowser for help on using the repository browser.