Context Navigation

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

Last change on this file since 207 was 207, checked in by aslmd, 13 years ago
MESOSCALE: A GENERAL CLEAN-UP FOLLOWING UPDATING THE USER MANUAL. EVERYTHING ESSENTIAL IS IN MESOSCALE (much lighter than before). EVERYTHING FOR DEVELOPPERS OR EXPERTS IS IN MESOSCALE_DEV.
Property svn:executable set to ``*
File size: 16.5 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=4,\
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	else:
55	print "typefile not supported."
56	print nc.variables
57	exit()
58
59	##############################################################
60	### Try to guess the projection from wrfout if not set by user
61	if typefile in ['mesoapi','meso']:
62	if proj == None: proj = getproj(nc)
63	### (il faudrait passer CEN_LON dans la projection ?)
64	elif typefile in ['gcm']:
65	if proj == None: proj = "cyl"
66	## pb avec les autres (de trace derriere la sphere ?)
67
68	############################################
69	#### Choose underlying topography by default
70	if not back:
71	if not var: back = "mola" ## if no var: draw mola
72	elif typefile in ['mesoapi','meso'] \
73	and proj not in ['merc','lcc','nsper','laea']: back = "molabw" ## if var but meso: draw molabw
74	else: pass ## else: draw None
75
76	####################################################
77	### Get geographical coordinates and plot boundaries
78	if typefile in ['mesoapi','meso']:
79	[lon2d,lat2d] = getcoord2d(nc)
80	lon2d = dumpbdy(lon2d)
81	lat2d = dumpbdy(lat2d)
82	elif typefile in ['gcm']:
83	[lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True)
84	if proj == "npstere": [wlon,wlat] = latinterv("North_Pole")
85	elif proj in ["lcc","laea"]: [wlon,wlat] = wrfinterv(lon2d,lat2d)
86	else: [wlon,wlat] = simplinterv(lon2d,lat2d)
87	if zoom:
88	dlon = abs(wlon[1]-wlon[0])/2.
89	dlat = abs(wlat[1]-wlat[0])/2.
90	[wlon,wlat] = [ [wlon[0]+zoomdlon/100.,wlon[1]-zoomdlon/100.],\
91	[wlat[0]+zoomdlat/100.,wlat[1]-zoomdlat/100.] ]
92	print "zoom %",zoom,wlon,wlat
93
94	##############################################################################
95	### Get winds and know if those are meteorological winds (ie. zon, mer) or not
96	if winds:
97	if typefile is 'mesoapi':
98	[u,v] = getwinds(nc)
99	metwind = True ## meteorological (zon/mer)
100	elif typefile is 'gcm':
101	[u,v] = getwinds(nc,charu='u',charv='v')
102	metwind = True ## meteorological (zon/mer)
103	elif typefile is 'meso':
104	[u,v] = getwinds(nc,charu='U',charv='V')
105	metwind = False ## geometrical (wrt grid)
106	print "Beware ! Not using meteorological winds. You trust numerical grid as being (x,y)."
107
108	#####################################################
109	### Load the chosen variables, whether it is 2D or 3D
110	if var:
111	if var not in nc.variables:
112	print "not found in file:",var
113	exit()
114	else:
115	dimension = np.array(nc.variables[var]).ndim
116	if dimension == 2: field = nc.variables[var][:,:]
117	elif dimension == 3: field = nc.variables[var][:,:,:]
118	elif dimension == 4: field = nc.variables[var][:,nvert,:,:]
119	dev = np.std(field)*3.0
120	if vmin is None: zevmin = mean(field) - dev
121	else: zevmin = vmin
122	if vmax is None: zevmax = mean(field) + dev
123	else: zevmax = vmax
124	print "bounds ", zevmin, zevmax
125	### some already defined colormaps
126	if colorb is True: colorb = defcolorb(var)
127	else:
128	dimension = 0
129
130	###########################
131	### Get length of time axis
132	if winds: nt = len(u[:,0,0,0])
133	elif var:
134	if dimension == 2: nt = 1
135	else : nt = len(field[:,0,0])
136
137	#########################################
138	### Name for title and graphics save file
139	if winds: basename = "UV_"
140	else: basename = ""
141	if var: basename = basename + var
142	###
143	if dimension == 4:
144	if typefile is 'meso': stralt = "_lvl" + str(nvert)
145	elif typefile is 'mesoapi':
146	zelevel = int(nc.variables['vert'][nvert])
147	if 'altitude' in nc.dimensions: stralt = "_"+str(zelevel)+"m-AMR"
148	elif 'altitude_abg' in nc.dimensions: stralt = "_"+str(zelevel)+"m-ALS"
149	elif 'bottom_top' in nc.dimensions: stralt = "_"+str(zelevel)+"m"
150	elif 'pressure' in nc.dimensions: stralt = "_"+str(zelevel)+"Pa"
151	else: stralt = ""
152	else:
153	stralt = ""
154	###
155	basename = basename + stralt
156
157	##################################
158	### Open a figure and set subplots
159	fig = figure()
160	if numplot > 0:
161	if numplot == 4:
162	sub = 221
163	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
164	rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
165	elif numplot == 2:
166	sub = 121
167	fig.subplots_adjust(wspace = 0.3)
168	rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
169	elif numplot == 3:
170	sub = 131
171	fig.subplots_adjust(wspace = 0.5)
172	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
173	elif numplot == 6:
174	sub = 231
175	fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
176	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
177	elif numplot == 8:
178	sub = 331 #241
179	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
180	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
181	elif numplot == 9:
182	sub = 331
183	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
184	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
185	elif numplot == 1:
186	sub = 99999
187	else:
188	print "supported: 1,2,3,4,6,8,9"
189	exit()
190	### Prepare time loop
191	if nt <= numplot or numplot == 1:
192	tabrange = [0]
193	numplot = 1
194	else:
195	tabrange = range(0,nt,int(nt/numplot)) #nt-1
196	tabrange = tabrange[0:numplot]
197	else:
198	tabrange = range(0,nt,1)
199	sub = 99999
200	print tabrange
201
202	#################################
203	### Time loop for plotting device
204	found_lct = False
205	for i in tabrange:
206
207	### Which local time ?
208	ltst = ( interv[0] + 0.5(wlon[0]+wlon[1])/15.) + iinterv[1]
209	ltst = int (ltst * 10) / 10.
210	ltst = ltst % 24
211
212	### General plot settings
213	if numplot > 1:
214	subplot(sub)
215	found_lct = True
216	elif numplot == 1:
217	found_lct = True
218	### If only one local time is requested (numplot < 0)
219	elif numplot <= 0:
220	if int(ltst) + numplot != 0: continue
221	else: found_lct = True
222
223	### Map projection
224	m = define_proj(proj,wlon,wlat,back=back)
225	x, y = m(lon2d, lat2d)
226
227	#### Contour plot
228	if var:
229	if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(field[i,:,:])
230	elif typefile in ['gcm']:
231	if dimension == 2: what_I_plot = field[:,:]
232	elif dimension == 3: what_I_plot = field[i,:,:]
233	palette = get_cmap(name=colorb)
234	#palette.set_over('b', 1.0)
235	if not tile:
236	zelevels = np.linspace(zevmin,zevmax)
237	contourf( x, y, what_I_plot, 10, cmap = palette, levels = zelevels )
238	else:
239	pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax )
240	if var in ['HGT']: pass
241	elif colorb:
242	ndiv = 10
243	colorbar(fraction=0.05,pad=0.1,format=fmtvar(var),\
244	ticks=np.linspace(zevmin,zevmax,ndiv+1),\
245	extend='max',spacing='proportional')
246	# both min max neither
247	### Vector plot
248	if winds:
249	if typefile in ['mesoapi','meso']:
250	[vecx,vecy] = [dumpbdy(u[i,nvert,:,:]), dumpbdy(v[i,nvert,:,:])]
251	key = True
252	elif typefile in ['gcm']:
253	[vecx,vecy] = [ u[i,nvert,:,:] , v[i,nvert,:,:] ]
254	key = False
255	if metwind: [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d)
256	if var == None: colorvec = definecolorvec(back)
257	else: colorvec = definecolorvec(colorb)
258	vectorfield(vecx, vecy,\
259	x, y, stride=stride, csmooth=stride,\
260	scale=15., factor=300., color=colorvec, key=key)
261	#200. ## or csmooth=2
262
263	### Next subplot
264	plottitle = basename
265	if addchar: plottitle = plottitle + addchar + "_LT"+str(ltst)
266	else: plottitle = plottitle + "_LT"+str(ltst)
267	ptitle( plottitle )
268	sub += 1
269
270	##########################################################################
271	### Save the figure in a file in the data folder or an user-defined folder
272	if typefile in ['meso','mesoapi']: prefix = getprefix(nc)
273	elif typefile in ['gcm']: prefix = 'LMD_GCM_'
274	else: prefix = ''
275	###
276	zeplot = prefix + basename
277	if addchar: zeplot = zeplot + addchar
278	if numplot <= 0: zeplot = zeplot + "_LT"+str(abs(numplot))
279	###
280	if not target: zeplot = namefile[0:find(namefile,'wrfout')] + zeplot
281	else: zeplot = target + "/" + zeplot
282	###
283	if found_lct: makeplotpng(zeplot,pad_inches_value=0.35)
284	else: print "Local time not found"
285
286	###############
287	### Now the end
288	return zeplot
289
290	###########################################################################################
291	###########################################################################################
292	### What is below relate to running the file as a command line executable (very convenient)
293	if __name__ == "__main__":
294	import sys
295	from optparse import OptionParser ### to be replaced by argparse
296	from api_wrapper import api_onelevel
297	from netCDF4 import Dataset
298	from myplot import getlschar
299
300	#############################
301	### Get options and variables
302	parser = OptionParser()
303	parser.add_option('-f', action='store', dest='namefile', type="string", default=None, help='[NEEDED] name of WRF file')
304	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)')
305	parser.add_option('-p', action='store', dest='proj', type="string", default=None, help='projection')
306	parser.add_option('-b', action='store', dest='back', type="string", default=None, help='background')
307	parser.add_option('-t', action='store', dest='target', type="string", default=None, help='destination folder')
308	parser.add_option('-s', action='store', dest='stride', type="int", default=3, help='stride vectors (def=3)')
309	parser.add_option('-v', action='store', dest='var', type="string", default=None, help='variable contoured')
310	parser.add_option('-n', action='store', dest='numplot', type="int", default=4, help='number of plots (def=1)(<0: 1 plot of LT -numplot)')
311	parser.add_option('-i', action='store', dest='interp', type="int", default=None, help='interpolation method (2: press, 3: z-amr, 4:z-als)')
312	parser.add_option('-c', action='store', dest='colorb', type="string", default=True, help='change colormap')
313	parser.add_option('-x', action='store_false', dest='winds', default=True, help='no wind vectors')
314	parser.add_option('-m', action='store', dest='vmin', type="float", default=None, help='bounding minimum value for color plot')
315	parser.add_option('-M', action='store', dest='vmax', type="float", default=None, help='bounding maximum value for color plot')
316	parser.add_option('-T', action='store_true', dest='tile', default=False, help='draw a tiled plot (no blank zone)')
317	parser.add_option('-z', action='store', dest='zoom', type="float", default=None, help='zoom factor in %')
318	parser.add_option('-N', action='store_true', dest='nocall', default=False, help='do not recreate api file')
319	#parser.add_option('-V', action='store', dest='comb', type="float", default=None, help='a defined combination of variables')
320	(opt,args) = parser.parse_args()
321	if opt.namefile is None:
322	print "I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h"
323	exit()
324	print "Options:", opt
325	zefile = opt.namefile
326	zelevel = opt.nvert
327	stralt = None
328	[lschar,zehour,zehourin] = getlschar ( zefile ) ## getlschar from wrfout (or simply return "" if another file)
329
330	#####################################################
331	### Call Fortran routines for vertical interpolations
332	if opt.interp is not None:
333	if opt.nvert is 0 and opt.interp is 4: zelevel = 0.010
334	### winds or no winds
335	if opt.winds : zefields = 'uvmet'
336	else : zefields = ''
337	### var or no var
338	if opt.var is None : pass
339	elif zefields == '' : zefields = opt.var
340	else : zefields = zefields + "," + opt.var
341	print zefields
342	zefile = api_onelevel ( path_to_input = '', \
343	input_name = zefile, \
344	path_to_output = opt.target, \
345	fields = zefields, \
346	interp_method = opt.interp, \
347	onelevel = zelevel, \
348	nocall = opt.nocall )
349	print zefile
350	zelevel = 0 ## so that zelevel could play again the role of nvert
351
352	#############
353	### Main call
354	name = winds (zefile,int(zelevel),\
355	proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var,numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\
356	addchar=lschar,interv=[zehour,zehourin],vmin=opt.vmin,vmax=opt.vmax,tile=opt.tile,zoom=opt.zoom)
357	print 'Done: '+name
358
359	#########################################################
360	### Generate a .sh file with the used command saved in it
361	command = ""
362	for arg in sys.argv: command = command + arg + ' '
363	f = open(name+'.sh', 'w')
364	f.write(command)

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

Download in other formats: