source: trunk/MESOSCALE/PLOT/PYTHON/scripts/winds.py @ 186

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

MESOSCALE: python post-processing. wrapper with my Fortran interpolator (small modifications done to api.F90). added the corresponding options to winds.py which is now a pretty complete script
  • Property svn:executable set to *
File size: 10.7 KB
Line 
1#!/usr/bin/env python
2
3### A. Spiga -- LMD -- 30/05/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=4,\
17           var=None): 
18
19    #################################
20    ### Load librairies and functions
21    from netCDF4 import Dataset
22    from myplot import getcoord2d,define_proj,makeplotpng,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy
23    from matplotlib.pyplot import contourf, subplot, figure, rcParams, savefig
24    import numpy as np
25
26    ######################
27    ### Load NETCDF object
28    nc  = Dataset(namefile)
29
30    ###################################
31    ### Recognize predefined file types
32    if 'controle' in nc.variables:   typefile = 'gcm'
33    elif 'Um' in nc.variables:       typefile = 'mesoapi'
34    elif 'U' in nc.variables:        typefile = 'meso'
35    else:                           
36        print "typefile not supported."
37        exit()
38
39    ##############################################################
40    ### Try to guess the projection from wrfout if not set by user
41    if typefile in ['mesoapi','meso']:
42        if proj == None:       proj = getproj(nc)
43                                    ### (il faudrait passer CEN_LON dans la projection ?)
44    elif typefile in ['gcm']:
45        if proj == None:       proj = "cyl"   
46                                    ## pb avec les autres (de trace derriere la sphere ?)
47
48    ####################################################################
49    #### For mesoscale results plot the underlying topography by default
50    if typefile in ['mesoapi','meso']: 
51        if var == None: back="mola" #var = 'HGT'
52
53    ####################################################
54    ### Get geographical coordinates and plot boundaries
55    if typefile in ['mesoapi','meso']:
56        [lon2d,lat2d] = getcoord2d(nc)
57        lon2d = dumpbdy(lon2d)
58        lat2d = dumpbdy(lat2d)
59    elif typefile in ['gcm']:
60        [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True)
61    if proj == "npstere":    [wlon,wlat] = latinterv("North_Pole")
62    elif proj == "lcc":      [wlon,wlat] = wrfinterv(lon2d,lat2d)
63    else:                    [wlon,wlat] = simplinterv(lon2d,lat2d)
64
65    ##################
66    ### Get local time
67    if typefile in ['mesoapi','meso']:  ltst = int(getattr(nc, 'GMT') + 0.5*(wlon[0]+wlon[1])/15.)
68    elif typefile in ['gcm']:           ltst = 0
69
70    ##############################################################################
71    ### Get winds and know if those are meteorological winds (ie. zon, mer) or not
72    if typefile is 'mesoapi': 
73        [u,v] = getwinds(nc)
74        metwind = True  ## meteorological (zon/mer)
75    elif typefile is 'gcm':
76        [u,v] = getwinds(nc,charu='u',charv='v')
77        metwind = True  ## meteorological (zon/mer)
78    elif typefile is 'meso':
79        [u,v] = getwinds(nc,charu='U',charv='V')
80        metwind = False ## geometrical (wrt grid)
81        print "Beware ! Not using meteorological winds. You trust numerical grid as being (x,y)."
82
83    ####################################################
84    ### Load the chosen variables, whether it is 2D or 3D
85    if var:
86        if var not in nc.variables: 
87            print "not found in file:",var
88            exit()
89        else:   
90            dimension = np.array(nc.variables[var]).ndim
91            if dimension == 2:     field = nc.variables[var][:,:]
92            elif dimension == 3:   field = nc.variables[var][:,:,:]
93            elif dimension == 4:   field = nc.variables[var][:,nvert,:,:] 
94    nt = len(u[:,0,0,0])
95
96    #########################################
97    ### Name for title and graphics save file
98    basename = "WINDS"
99    if var:
100        basename = basename + "_" + var
101    basename = basename + "_z" + str(nvert)
102
103    ##################################
104    ### Open a figure and set subplots
105    fig = figure()
106    if   numplot > 0:   
107        if   numplot == 4: 
108            sub = 221
109            fig.subplots_adjust(wspace = 0.1, hspace = 0.3)
110            rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
111        elif numplot == 2: 
112            sub = 121
113            fig.subplots_adjust(wspace = 0.3)
114            rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
115        elif numplot == 3: 
116            sub = 131
117            fig.subplots_adjust(wspace = 0.3)
118            rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
119        elif numplot == 6: 
120            sub = 231
121            fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
122            rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
123        elif numplot == 8: 
124            sub = 331 #241
125            fig.subplots_adjust(wspace = 0.1, hspace = 0.3)
126            rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
127        elif numplot == 9:
128            sub = 331 
129            fig.subplots_adjust(wspace = 0.1, hspace = 0.3)
130            rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
131        elif numplot == 1:
132            pass
133        else:
134            print "supported: 1,2,3,4,6,8"
135            exit()
136        ### Prepare time loop
137        if nt <= numplot or numplot == 1: 
138            tabrange = [0]
139            numplot = 1
140        else:                         
141            tabrange = range(0,nt,int(nt/numplot))  #nt-1
142            tabrange = tabrange[0:numplot]
143    else: 
144        tabrange = range(0,nt,1)
145        sub = 99999
146    print tabrange
147
148    #################################
149    ### Time loop for plotting device
150    found_lct = False
151    for i in tabrange:
152
153       ### General plot settings
154       if numplot > 1: 
155           subplot(sub)
156           found_lct = True
157       elif numplot == 1:
158           found_lct = True 
159        ### If only one local time is requested (numplot < 0)
160       elif numplot <= 0: 
161           #print (ltst+i)%24, numplot, (ltst+i)%24+numplot
162           if (ltst+i)%24 + numplot != 0:   continue
163           else:                            found_lct = True
164
165       ### Map projection
166       m = define_proj(proj,wlon,wlat,back=back)
167       x, y = m(lon2d, lat2d)
168
169       #### Contour plot
170       if var:
171           if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(field[i,:,:])
172           elif typefile in ['gcm']:             
173               if dimension == 2:                what_I_plot = field[:,:]
174               elif dimension == 3:              what_I_plot = field[i,:,:]
175           contourf(x, y, what_I_plot, 30)
176
177       ### Vector plot
178       if   typefile in ['mesoapi','meso']:   
179           [vecx,vecy] = [dumpbdy(u[i,nvert,:,:]), dumpbdy(v[i,nvert,:,:])]
180           key = True
181       elif typefile in ['gcm']:               
182           [vecx,vecy] = [        u[i,nvert,:,:] ,         v[i,nvert,:,:] ]
183           key = False
184       if metwind:  [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d)
185       vectorfield(vecx, vecy,\
186                      x, y, stride=stride, csmooth=stride,\
187                      scale=15., factor=200., color='k', key=key)
188                                                   ## or csmooth=2
189       
190       ### Next subplot
191       ptitle( basename + "_LT"+str((ltst+i)%24) )
192       sub += 1
193
194    ##########################################################################
195    ### Save the figure in a file in the data folder or an user-defined folder
196    if not target:    zeplot = namefile+"_"+basename
197    else:             zeplot = target+"/"+basename
198    if numplot <= 0:  zeplot = zeplot + "_LT"+str(abs(numplot))
199    if found_lct:     makeplotpng(zeplot,pad_inches_value=0.35)   
200    else:             print "Local time not found"
201
202
203
204####################################################
205####################################################
206### A simple program to get wind vectors' components
207def getwinds (nc,charu='Um',charv='Vm'):
208    import numpy as np
209    u = nc.variables[charu]
210    v = nc.variables[charv]
211    if charu == 'U': u = u[:, :, :, 0:len(u[0,0,0,:])-1]
212    if charv == 'V': v = v[:, :, 0:len(v[0,0,:,0])-1, :]
213                     ### ou alors prendre les coordonnees speciales
214    return u,v
215
216
217
218###########################################################################################
219###########################################################################################
220### What is below relate to running the file as a command line executable (very convenient)
221if __name__ == "__main__":
222    import sys
223    from optparse import OptionParser    ### to be replaced by argparse
224    from api_wrapper import api_onelevel
225    parser = OptionParser()
226    parser.add_option('-f', action='store', dest='namefile',    type="string",  default=None,  help='name of WRF file [NEEDED]')
227    parser.add_option('-l', action='store', dest='nvert',       type="float",   default=0,     help='vertical level (def=0)')
228    parser.add_option('-p', action='store', dest='proj',        type="string",  default=None,  help='projection')
229    parser.add_option('-b', action='store', dest='back',        type="string",  default=None,  help='background')
230    parser.add_option('-t', action='store', dest='target',      type="string",  default=None,  help='destination folder')
231    parser.add_option('-s', action='store', dest='stride',      type="int",     default=3,     help='stride vectors (def=3)')
232    parser.add_option('-v', action='store', dest='var',         type="string",  default=None,  help='variable contoured')
233    parser.add_option('-n', action='store', dest='numplot',     type="int",     default=4,     help='number of plots (def=1)(if <0: 1 plot of LT -*numplot*)')
234    parser.add_option('-i', action='store', dest='interp',      type="int",     default=None,  help='interpolation method (done at level *nvert* km)')
235    (opt,args) = parser.parse_args()
236    if opt.namefile is None: 
237        print "I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h"
238        exit()
239    print "Options:", opt
240
241    zefile = opt.namefile   
242    zelevel = opt.nvert   
243    if opt.nvert is 0 and opt.interp:   zelevel = 0.020 
244    if opt.interp is not None:
245        if   opt.var is None    :  zefields = 'uvmet'
246        else                    :  zefields = 'uvmet,'+opt.var
247        zefile = api_onelevel (  path_to_input   = '', \
248                                 input_name      = zefile, \
249                                 path_to_output  = opt.target, \
250                                 fields          = zefields, \
251                                 interp_method   = opt.interp, \
252                                 onelevel        = zelevel )
253        zelevel = 0
254
255    winds (zefile,int(zelevel),proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var,numplot=opt.numplot)
Note: See TracBrowser for help on using the repository browser.