Context Navigation

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

Last change on this file since 232 was 232, checked in by aslmd, 14 years ago
MESOSCALE: corrections petits bugs graphiques.
Property svn:executable set to ``*
File size: 19.6 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
10	def 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,putpoints
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]+zoomdlon/100.,wlon[1]-zoomdlon/100.],\
96	[wlat[0]+zoomdlat/100.,wlat[1]-zoomdlat/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	fieldcalc = field[ field < 9e+35 ]
127	dev = np.std(fieldcalc)*2.0
128	if vmin is None: zevmin = mean(fieldcalc) - dev
129	else: zevmin = vmin
130	if vmax is None: zevmax = mean(fieldcalc) + dev
131	else: zevmax = vmax
132	if vmin == vmax:
133	zevmin = mean(fieldcalc) - dev ### for continuity
134	zevmax = mean(fieldcalc) + dev ### for continuity
135	print "field ", min(fieldcalc), max(fieldcalc)
136	print "bounds ", zevmin, zevmax
137	### some already defined colormaps
138	if colorb is True: colorb = defcolorb(var)
139	else:
140	dimension = 0
141
142	###########################
143	### Get length of time axis
144	if winds: nt = len(u[:,0,0,0])
145	elif var:
146	if dimension == 2: nt = 1
147	else : nt = len(field[:,0,0])
148
149	#########################################
150	### Name for title and graphics save file
151	if var and winds: basename = var + '_UV'
152	elif var: basename = var
153	elif winds: basename = 'UV'
154	else: exit()
155	###
156	if dimension == 4 or winds:
157	if typefile is 'meso': stralt = "_lvl" + str(nvert)
158	elif typefile is 'mesoapi':
159	zelevel = int(nc.variables['vert'][nvert])
160	if abs(zelevel) < 10000.: strheight=str(zelevel)+"m"
161	else: strheight=str(int(zelevel/1000.))+"km"
162	if 'altitude' in nc.dimensions: stralt = "_"+strheight+"-AMR"
163	elif 'altitude_abg' in nc.dimensions: stralt = "_"+strheight+"-ALS"
164	elif 'bottom_top' in nc.dimensions: stralt = "_"+strheight
165	elif 'pressure' in nc.dimensions: stralt = "_"+str(zelevel)+"Pa"
166	else: stralt = ""
167	else:
168	stralt = ""
169	###
170	basename = basename + stralt
171
172	##################################
173	### Open a figure and set subplots
174	fig = figure()
175	rcParams['font.size'] = 12.
176	if typefile in ['geo']: numplot = 1
177	if numplot > 0:
178	if numplot == 4:
179	sub = 221
180	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
181	rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
182	elif numplot == 2:
183	sub = 121
184	fig.subplots_adjust(wspace = 0.35)
185	rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
186	elif numplot == 3:
187	sub = 131
188	fig.subplots_adjust(wspace = 0.5)
189	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
190	elif numplot == 6:
191	sub = 231
192	fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
193	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
194	elif numplot == 8:
195	sub = 331 #241
196	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
197	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
198	elif numplot == 9:
199	sub = 331
200	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
201	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
202	elif numplot == 1:
203	sub = 99999
204	else:
205	print "supported: 1,2,3,4,6,8,9"
206	exit()
207	### Prepare time loop
208	if nt <= numplot or numplot == 1:
209	tabrange = [0]
210	numplot = 1
211	else:
212	tabrange = range(0,nt,int(nt/numplot)) #nt-1
213	tabrange = tabrange[0:numplot]
214	else:
215	tabrange = range(0,nt,1)
216	sub = 99999
217	print tabrange
218
219	#################################
220	### Time loop for plotting device
221	found_lct = False
222	for i in tabrange:
223
224	### Which local time ?
225	ltst = ( interv[0] + 0.5(wlon[0]+wlon[1])/15.) + iinterv[1]
226	ltst = int (ltst * 10) / 10.
227	ltst = ltst % 24
228
229	### General plot settings
230	if numplot > 1:
231	subplot(sub)
232	found_lct = True
233	elif numplot == 1:
234	found_lct = True
235	### If only one local time is requested (numplot < 0)
236	elif numplot <= 0:
237	if int(ltst) + numplot != 0: continue
238	else: found_lct = True
239
240	### Map projection
241	m = define_proj(proj,wlon,wlat,back=back)
242	x, y = m(lon2d, lat2d)
243
244	#### Contour plot
245	if var:
246	if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(field[i,:,:])
247	elif typefile in ['geo']: what_I_plot = field[0,:,:]
248	elif typefile in ['gcm']:
249	if dimension == 2: what_I_plot = field[:,:]
250	elif dimension == 3: what_I_plot = field[i,:,:]
251	palette = get_cmap(name=colorb)
252	#palette.set_over('b', 1.0)
253	#print np.array(x).shape
254	#print np.array(y).shape
255	#print np.array(what_I_plot).shape
256	if not tile:
257	zelevels = np.linspace(zevmin,zevmax)
258	hole = True
259	if not hole:
260	what_I_plot[ what_I_plot < zevmin ] = zevmin*(1. + 1.e-7)
261	what_I_plot[ what_I_plot > 9e+35 ] = -9e+35
262	what_I_plot[ what_I_plot > zevmax ] = zevmax*(1. - 1.e-7)
263	contourf( x, y, what_I_plot, 10, cmap = palette, levels = zelevels )
264	else:
265	pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax )
266	#putpoints(m,fig)
267	if var in ['HGT']: pass
268	elif colorb:
269	ndiv = 10
270	colorbar(fraction=0.05,pad=0.1,format=fmtvar(var),\
271	ticks=np.linspace(zevmin,zevmax,ndiv+1),\
272	extend='max',spacing='proportional')
273	# both min max neither
274	### Vector plot
275	if winds:
276	if typefile in ['mesoapi','meso']:
277	[vecx,vecy] = [dumpbdy(u[i,nvert,:,:]), dumpbdy(v[i,nvert,:,:])]
278	key = True
279	elif typefile in ['gcm']:
280	[vecx,vecy] = [ u[i,nvert,:,:] , v[i,nvert,:,:] ]
281	key = False
282	if metwind: [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d)
283	if var == None: colorvec = definecolorvec(back)
284	else: colorvec = definecolorvec(colorb)
285	vectorfield(vecx, vecy,\
286	x, y, stride=stride, csmooth=2,\
287	scale=15., factor=300., color=colorvec, key=key)
288	#200. ## or csmooth=stride
289
290	### Next subplot
291	plottitle = basename
292	if typefile in ['mesoapi','meso']:
293	if addchar: plottitle = plottitle + addchar + "_LT"+str(ltst)
294	else: plottitle = plottitle + "_LT"+str(ltst)
295	ptitle( plottitle )
296	sub += 1
297
298	##########################################################################
299	### Save the figure in a file in the data folder or an user-defined folder
300	if typefile in ['meso','mesoapi']: prefix = getprefix(nc)
301	elif typefile in ['gcm']: prefix = 'LMD_GCM_'
302	else: prefix = ''
303	###
304	zeplot = prefix + basename
305	if addchar: zeplot = zeplot + addchar
306	if numplot <= 0: zeplot = zeplot + "_LT"+str(abs(numplot))
307	###
308	if not target: zeplot = namefile[0:find(namefile,'wrfout')] + zeplot
309	else: zeplot = target + "/" + zeplot
310	###
311	if found_lct: makeplotpng(zeplot,pad_inches_value=0.35)
312	else: print "Local time not found"
313
314	###############
315	### Now the end
316	return zeplot
317
318	###########################################################################################
319	###########################################################################################
320	### What is below relate to running the file as a command line executable (very convenient)
321	if __name__ == "__main__":
322	import sys
323	from optparse import OptionParser ### to be replaced by argparse
324	from api_wrapper import api_onelevel
325	from netCDF4 import Dataset
326	from myplot import getlschar
327	from numpy import ones
328
329	#############################
330	### Get options and variables
331	parser = OptionParser()
332	parser.add_option('-f', action='append', dest='namefile', type="string", default=None, help='[NEEDED] name of WRF file (append)')
333	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)')
334	parser.add_option('-p', action='store', dest='proj', type="string", default=None, help='projection')
335	parser.add_option('-b', action='store', dest='back', type="string", default=None, help='background')
336	parser.add_option('-t', action='store', dest='target', type="string", default=None, help='destination folder')
337	parser.add_option('-s', action='store', dest='stride', type="int", default=3, help='stride vectors (def=3)')
338	parser.add_option('-v', action='append', dest='var', type="string", default=None, help='variable contoured (append)')
339	parser.add_option('-n', action='store', dest='numplot', type="int", default=2, help='number of plots (def=1)(<0: 1 plot of LT -numplot)')
340	parser.add_option('-i', action='store', dest='interp', type="int", default=None, help='interpolation method (2: press, 3: z-amr, 4:z-als)')
341	parser.add_option('-c', action='store', dest='colorb', type="string", default=True, help='change colormap')
342	parser.add_option('-x', action='store_false', dest='winds', default=True, help='no wind vectors')
343	parser.add_option('-m', action='append', dest='vmin', type="float", default=None, help='bounding minimum value (append)')
344	parser.add_option('-M', action='append', dest='vmax', type="float", default=None, help='bounding maximum value (append)')
345	parser.add_option('-T', action='store_true', dest='tile', default=False, help='draw a tiled plot (no blank zone)')
346	parser.add_option('-z', action='store', dest='zoom', type="float", default=None, help='zoom factor in %')
347	parser.add_option('-N', action='store_true', dest='nocall', default=False, help='do not recreate api file')
348	#parser.add_option('-V', action='store', dest='comb', type="float", default=None, help='a defined combination of variables')
349	(opt,args) = parser.parse_args()
350	if opt.namefile is None:
351	print "I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h"
352	exit()
353	print "Options:", opt
354
355	if opt.var is None:
356	pass
357	else:
358	listvar = ''
359	zelen = len(opt.var)
360	if zelen == 1: listvar = opt.var[0] + ','
361	else :
362	for jjj in range(zelen): listvar += opt.var[jjj] + ','
363	listvar = listvar[0:len(listvar)-1]
364	if opt.vmin:
365	if zelen != len(opt.vmin):
366	print "not enough or too much vmin values... setting same values all variables"
367	vmintab = ones(zelen) * opt.vmin[0]
368	else:
369	vmintab = opt.vmin
370	else:
371	vmintab = None
372	if opt.vmax:
373	if zelen != len(opt.vmax):
374	print "not enough or too much vmax values... setting same values all variables"
375	vmaxtab = ones(zelen) * opt.vmax[0]
376	else:
377	vmaxtab = opt.vmax
378	else:
379	vmaxtab = None
380
381	for i in range(len(opt.namefile)):
382
383	zefile = opt.namefile[i]
384	print zefile
385	zelevel = opt.nvert
386	stralt = None
387	[lschar,zehour,zehourin] = getlschar ( zefile ) ## getlschar from wrfout (or simply return "" if another file)
388
389	#####################################################
390	### Call Fortran routines for vertical interpolations
391	if opt.interp is not None:
392	if opt.nvert is 0 and opt.interp is 4: zelevel = 0.010
393	### winds or no winds
394	if opt.winds : zefields = 'uvmet'
395	else : zefields = ''
396	### var or no var
397	if opt.var is None : pass
398	elif zefields == '' : zefields = listvar
399	else : zefields = zefields + "," + listvar
400	print zefields
401	zefile = api_onelevel ( path_to_input = '', \
402	input_name = zefile, \
403	#path_to_output = opt.target, \
404	fields = zefields, \
405	interp_method = opt.interp, \
406	onelevel = zelevel, \
407	nocall = opt.nocall )
408	print zefile
409	zelevel = 0 ## so that zelevel could play again the role of nvert
410
411	if opt.var is None:
412	#############
413	### Main call
414	name = winds (zefile,int(zelevel),\
415	proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var,numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\
416	addchar=lschar,interv=[zehour,zehourin],vmin=opt.vmin,vmax=opt.vmax,tile=opt.tile,zoom=opt.zoom)
417	print 'Done: '+name
418	else:
419	for jjj in range(len(opt.var)):
420	if vmintab: argvmin = vmintab[jjj]
421	else: argvmin = None
422	if vmaxtab: argvmax = vmaxtab[jjj]
423	else: argvmax = None
424	#############
425	### Main call
426	name = winds (zefile,int(zelevel),\
427	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,\
428	addchar=lschar,interv=[zehour,zehourin],vmin=argvmin,vmax=argvmax,tile=opt.tile,zoom=opt.zoom)
429	print 'Done: '+name
430
431	#########################################################
432	### Generate a .sh file with the used command saved in it
433	command = ""
434	for arg in sys.argv: command = command + arg + ' '
435	f = open(name+'.sh', 'w')
436	f.write(command)

Note: See TracBrowser for help on using the repository browser.

Download in other formats: