Context Navigation

source: trunk/UTIL/PYTHON/planetoplot.py @ 782

Last change on this file since 782 was 782, checked in by acolaitis, 12 years ago
UTIL PYTHON fixed mesoscale local time bugs introduced at revision 763. now working fine for LMD and MRAMS-generated wrfout files. plus polar domain now OK
Property svn:executable set to ``*
File size: 59.4 KB

Line
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	### A. Colaitis -- LMD -- 11/2011 -- Mostly minor improvements and inter-plot operation capabilities + zrecast interpolation for gcm
8	### A. Spiga -- LMD -- 11~12/2011 -- Extended multivar subplot capabilities + cosmetic changes + general cleaning and tests
9	### A. Colaitis -- LMD -- 12/2011 -- Added movie capability [mencoder must be installed]
10	### A. Spiga -- LMD -- 12/2011 -- Added HTML animated page capability + general tests of consistency [winds, etc...] + consistent generic movie loop
11	### J. Leconte -- LMD -- 02/2012 -- Added area weighted averaging. Compatibility with terrestrial gcm.
12	### A. Spiga -- LMD -- 03/2012 -- Cleaning and improved comments
13	### T. Navarro -- LMD -- 04/2012 -- Added capabilities (e.g. histograms for difference maps)
14	### A. Colaitis -- LMD -- 05~06/2012 -- Added capabilities for analysis of mesoscale files (wind speed, etc...)
15	### A. Spiga -- LMD -- 04~07/2012 -- Added larger support of files + planets. Corrected a few bugs. Cleaning and improved comments
16	### A. Colaitis -- LMD -- 08/2012 -- Added functionalities for further analysis: spectra, hodographs, etc... plus improved flexibility (e.g. grid)
17
18	def planetoplot (namefiles,\
19	level=0,\
20	vertmode=0,\
21	proj=None,\
22	back=None,\
23	target=None,
24	stride=3,\
25	var=None,\
26	clb=None,\
27	winds=False,\
28	addchar=None,\
29	vmin=None,\
30	vmax=None,\
31	tile=False,\
32	zoom=None,\
33	display=True,\
34	hole=False,\
35	save="gui",\
36	anomaly=False,\
37	var2=None,\
38	ndiv=10,\
39	mult=1.,\
40	zetitle=["fill"],\
41	slon=None,\
42	slat=None,\
43	svert=None,\
44	stime=None,\
45	outputname=None,\
46	resolution=200,\
47	ope=None,\
48	fileref=None,\
49	minop=0.,\
50	maxop=0.,\
51	titleref="fill",\
52	invert_y=False,\
53	xaxis=[None,None],\
54	yaxis=[None,None],\
55	ylog=False,\
56	xlog=False,\
57	yintegral=False,\
58	blat=None,\
59	blon=None,\
60	tsat=False,\
61	flagnolow=False,\
62	mrate=None,\
63	mquality=False,\
64	trans=1,\
65	zarea=None,\
66	axtime=None,\
67	redope=None,\
68	seevar=False,\
69	xlab=None,\
70	ylab=None,\
71	lbls=None,\
72	lstyle=None,\
73	cross=None,\
74	facwind=1.,\
75	trycol=False,\
76	streamflag=False,\
77	analysis=None):
78
79	####################################################################################################################
80	### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
81
82	#################################
83	### Load librairies and functions
84	from netCDF4 import Dataset
85	from myplot import getcoord2d,define_proj,makeplotres,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy,\
86	fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\
87	zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow,\
88	getname,localtime,check_localtime,polarinterv,getsindex,define_axis,determineplot,readslices,bidimfind,getlschar,hole_bounds,\
89	getdimfromvar,select_getfield
90	from mymath import deg,max,min,mean,writeascii,fig2data,fig2img
91	import matplotlib as mpl
92	from matplotlib.pyplot import contour,contourf,hist, text,subplot, figure, rcParams, savefig, colorbar, \
93	pcolor, show, plot, clabel, title, close, legend, xlabel, axis, ylabel, \
94	subplots_adjust, axes, clabel
95	from matplotlib.cm import get_cmap
96	#from mpl_toolkits.basemap import cm
97	import numpy as np
98	from numpy.core.defchararray import find
99	from scipy.stats import gaussian_kde,describe
100	from videosink import VideoSink
101	from times import sol2ls
102	import subprocess
103	#from singlet import singlet
104	from itertools import cycle
105	import os
106
107	#########################
108	### PRELIMINARY STUFF ###
109	#########################
110	### we say hello
111	print "********************************************"
112	print "******** WELCOME TO PLANETOPLOT ********"
113	print "********************************************"
114	### we ensure namefiles and var are arrays
115	if not isinstance(namefiles, np.ndarray): namefiles = [namefiles]
116	if not isinstance(var, np.ndarray): var = [var]
117	### we initialize a few variables
118	initime=-1 ; sslon = None ; sslat = None
119	k = 0 ; firstfile = True ; count = 0
120	### we perform sanity checks and correct for insufficient information
121	if slon is not None: sslon = np.zeros([1,2])
122	if slat is not None: sslat = np.zeros([1,2])
123	if clb is None: clb = ["def"]*len(var)
124	elif len(clb) < len(var): clb = [clb[0]]*len(var) ; print "WARNING: less color than vars! setting all to 1st value."
125	### we set option trycol i.e. the user wants to try a set of colorbars
126	if trycol: clb = ["Greys","Blues","YlOrRd","jet","spectral","hot","RdBu","RdYlBu","Paired"] ; zetitle = clb ; var = [var[0]]*9
127	### we avoid specific cases not yet implemented
128	if mrate is not None and len(var) > 1: errormess("multivar not allowed in movies. should be fixed soon!")
129	### we prepare number of plot fields [zelen] and number of plot instances [numplot] according to user choices
130	### --> we support multifile and multivar plots : files + vars separated by commas are on the same figure
131	nlon, nlat, nvert, ntime, mapmode, nslices = determineplot(slon, slat, svert, stime, redope)
132	zelen = len(namefiles)*len(var)
133	### we correct number of plot fields for possible operation (substract, etc...)
134	if ope is not None:
135	if fileref is not None: zelen = 3len(var)len(namefiles)
136	elif "var" in ope: zelen = zelen + 1
137	numplot = zelen*nslices
138	print "********** FILES, SLICES, VARS, TOTAL PLOTS: ", len(namefiles), nslices, len(var), numplot
139	print "********** MAPMODE: ", mapmode
140	### we define the arrays for plot fields -- which will be placed within multiplot loops
141	all_var = [[]]zelen ; all_var2 = [[]]zelen
142	all_title = [[]]zelen ; all_varname = [[]]zelen ; all_namefile = [[]]zelen ; all_colorb = [[]]zelen
143	all_time = [[]]zelen ; all_vert = [[]]zelen ; all_lat = [[]]zelen ; all_lon = [[]]zelen
144	all_windu = [[]]zelen ; all_windv = [[]]zelen
145	plot_x = [[]]zelen ; plot_y = [[]]zelen ; multiplot = [[]]*zelen
146	### tool (should be move to mymath)
147	getVar = lambda searchList, ind: [searchList[i] for i in ind]
148	#############################
149	### LOOP OVER PLOT FIELDS ###
150	#############################
151
152	for nnn in range(len(namefiles)):
153	for vvv in range(len(var)):
154
155	### we load each NETCDF objects in namefiles
156	namefile = namefiles[nnn]
157	nc = Dataset(namefile)
158	### we explore the variables in the file
159	varinfile = nc.variables.keys()
160	if seevar: print varinfile ; exit()
161	### we define the type of file we have (gcm, meso, etc...)
162	typefile = whatkindfile(nc)
163	if firstfile: typefile0 = typefile
164	elif typefile != typefile0: errormess("Not the same kind of files !", [typefile0, typefile])
165	### we care for input file being 1D simulations
166	is1d=999
167	if "longitude" in nc.dimensions and "latitude" in nc.dimensions: is1d = len(nc.variables["longitude"][:])*len(nc.variables["latitude"][:])
168	elif "lon" in nc.dimensions and "lat" in nc.dimensions: is1d = len(nc.variables["lon"][:])*len(nc.variables["lat"][:])
169	if typefile in ['gcm','earthgcm'] and is1d == 1: mapmode=0 ; winds=False
170	### we create default vert and time prescriptions if not here in case mapping mode is on (lat/lon maps)
171	if redope is None and mapmode == 1:
172	if svert is None: svert = readslices(str(level)) ; nvert=1
173	if stime is None and mrate is None:
174	stime = readslices(str(0)) ; ntime=1 ## this is a default choice
175	print "WELL... nothing about time axis. I took default: first time reference stored in file."
176	### we get the names of variables to be read. in case only one available, we choose this one.
177	### (we take out of the test specific names e.g. UV is not in the file but used to ask a wind speed computation)
178	varname = select_getfield(zvarname=var[vvv],znc=nc,ztypefile=typefile,mode='check')
179	### we get the names of wind variables to be read (if any)
180	if winds:
181	[uchar,vchar,metwind] = getwinddef(nc)
182	if uchar == 'not found': winds = False
183	### we tell the user that either no var or no wind is not acceptable
184	if not varname and not winds: errormess("please set at least winds or var",printvar=nc.variables)
185	### we get the coordinates lat/lon to be used
186	[lon2d,lat2d] = getcoorddef(nc)
187	### we get an adapted map projection if none is provided by the user
188	if proj == None: proj = getproj(nc)
189	### we define plot boundaries given projection or user choices
190	if firstfile:
191	if proj in ["npstere","spstere"]: [wlon,wlat] = polarinterv(lon2d,lat2d)
192	elif proj in ["lcc","laea"]: [wlon,wlat] = wrfinterv(lon2d,lat2d)
193	else: [wlon,wlat] = simplinterv(lon2d,lat2d)
194	if zoom: [wlon,wlat] = zoomset(wlon,wlat,zoom)
195	elif zarea is not None: [wlon,wlat] = latinterv(area=zarea)
196
197	#############################################################
198	############ WE LOAD 4D DIMENSIONS : x, y, z, t #############
199	#############################################################
200
201	###should be available for GCM or simple files ?
202	##if ope in ["cat"] and nnn > 0: count = time[-1] + 1 ## so that a cat is possible with simple subscripts
203
204	### TYPE 1 : GCM files or simple files
205	if typefile in ["gcm","earthgcm","ecmwf"]:
206	### this is needed for continuity
207	if slon is not None: sslon = slon
208	if slat is not None: sslat = slat
209	### we define lat/lon vectors. we get what was done in getcoorddef.
210	lat = lat2d[:,0] ; lon = lon2d[0,:]
211	### we define areas. this is needed for calculate means and weight with area. this is not compulsory (see reduce_field).
212	if "aire" in nc.variables: area = nc.variables["aire"][:,:]
213	### --> add a line here if your reference is not present
214	else: area = None
215	### we define altitude vector. either it is referenced or it is guessed based on last variable's dimensions.
216	if "altitude" in nc.variables: vert = nc.variables["altitude"][:]
217	elif "Alt" in nc.variables: vert = nc.variables["Alt"][:]
218	elif "lev" in nc.variables: vert = nc.variables["lev"][:]
219	elif "presnivs" in nc.variables: vert = nc.variables["presnivs"][:]
220	### --> add a line here if your reference is not present
221	else:
222	print "No altitude found. Try to build a simple axis." ; dadim = getdimfromvar(nc)
223	if len(dadim) == 4: print "-- 4D field. Assume z is dim 2." ; vert = np.arange(dadim[-3])
224	elif len(dadim) == 3: print "-- 3D field. Assume z is dim 1." ; vert = [0.]
225	else: vert = [0.]
226	### we define time vector. either it is referenced or it is guessed based on last variable's dimensions.
227	if "Time" in nc.variables: time = nc.variables["Time"][:]
228	elif "time_counter" in nc.variables: time = nc.variables["time_counter"][:]/86400. #### convert from s to days
229	elif "time" in nc.variables: time = nc.variables["time"][:]
230	### --> add a line here if your reference is not present
231	else:
232	print "No time found. Try to build a simple axis. Assume t is dim 1."
233	dadim = getdimfromvar(nc)
234	if len(dadim) == 4: time = np.arange(dadim[-4])
235	elif len(dadim) == 3: time = np.arange(dadim[-3])
236	else: time = [0.] #errormess("no time axis found.")
237	### (SPECIFIC. convert to Ls for Martian GCM files.)
238	if axtime in ["ls"]:
239	print "converting to Ls ..."
240	for iii in range(len(time)):
241	time[iii] = sol2ls(time[iii])
242	if iii > 0:
243	while abs(time[iii]-time[iii-1]) > 300: time[iii] = time[iii]+360
244	### (a case where the user would like to set local times. e.g. : 1D plot with no longitude reference)
245	if axtime in ["lt"]:
246	if initime == -1: initime=input("Please type initial local time:")
247	time = (initime+time*24)%24 ; print "LOCAL TIMES.... ", time
248
249	### TYPE 2 : MESOSCALE FILES
250	elif typefile in ['meso','geo']:
251	### area not active with mesoscale files
252	area = None
253	### HACK TO GET THE CORRECT LAT/LON FROM MESOSCALE FILES WITH 2D LAT/LON ARRAYS
254	### principle: calculate correct indices then repopulate slon and slat
255	if slon is not None or slat is not None:
256	show_topo_map_user = ((save == 'png') and (typefile == 'meso') and ("HGT" in varinfile) and display)
257	if firstfile and show_topo_map_user: iwantawhereplot = nc #show a topo map with a cross on the chosen point
258	else: iwantawhereplot = None #do not show anything, just select indices
259	numlon = 1 ; numlat = 1
260	if slon is not None: numlon = slon.shape[1]
261	if slat is not None: numlat = slat.shape[1]
262	indices = np.ones([numlon,numlat,2]) ; vlon = None ; vlat = None
263	for iii in range(numlon):
264	for jjj in range(numlat):
265	if slon is not None: vlon = slon[0][iii] ### note: slon[:][0] does not work
266	if slat is not None: vlat = slat[0][jjj] ### note: slon[:][0] does not work
267	indices[iii,jjj,:] = bidimfind(lon2d,lat2d,vlon,vlat,file=iwantawhereplot)
268	lonp,latp = ( lon2d[indices[iii,jjj,0],indices[iii,jjj,1]] , lat2d[indices[iii,jjj,0],indices[iii,jjj,1]] )
269	### possible bug here if several --lat
270	for iii in range(numlon):
271	for jjj in range(numlat):
272	if slon is not None: sslon[0][iii] = indices[iii,0,1] #...this is idx
273	if slat is not None: sslat[0][jjj] = indices[0,jjj,0] #...this is idy
274	lonp,latp = ( lon2d[indices[0,0,0],indices[0,0,1]] , lat2d[indices[0,0,0],indices[0,0,1]] )
275	### we get rid of boundary relaxation zone for plots. important to do that now and not before.
276	if typefile in ['meso'] and mapmode == 1: lon2d = dumpbdy(lon2d,6) ; lat2d = dumpbdy(lat2d,6)
277	### we read the keyword for vertical dimension. we take care for vertical staggering.
278	if varname in ['PHTOT','W']: vertdim='BOTTOM-TOP_PATCH_END_STAG'
279	else: vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
280	if (var2 is not None and var2 not in ['PHTOT','W']): dumped_vert_stag=True ; vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
281	else: dumped_vert_stag=False
282	### we read the keyword for horizontal dimensions. we take care for horizontal staggering.
283	if varname in ['V']: latdim='SOUTH-NORTH_PATCH_END_STAG'
284	else: latdim='SOUTH-NORTH_PATCH_END_UNSTAG'
285	if varname in ['U']: londim='WEST-EAST_PATCH_END_STAG'
286	else: londim='WEST-EAST_PATCH_END_UNSTAG'
287	lon = np.arange(0,getattr(nc,londim),1) ; lat = np.arange(0,getattr(nc,latdim),1)
288	### we define the time axis and take care of various specificities (lt, ls, sol) or issues (concatenation)
289	if axtime in ["ls","sol"]:
290	lstab, soltab, lttab = getlschar ( namefile, getaxis = True )
291	if axtime == "ls": time = lstab
292	elif axtime == "sol": time = soltab
293	else:
294	if ope in ["cat"] and nnn > 0: count = time[-1] + 1 ## so that a cat is possible with simple subscripts
295	else: count = 0
296	if "Times" in nc.dimensions: time = count + np.arange(0,len(nc.dimensions["Times"]),1)
297	elif "Time" in nc.dimensions: time = count + np.arange(0,len(nc.dimensions["Time"]),1)
298	else: time = count + np.arange(0,1,1)
299	if axtime in ["lt"]:
300	ftime = np.zeros(len(time))
301	for i in range(len(time)): ftime[i] = localtime ( time[i], slon , namefile )
302	time=ftime ; time=check_localtime(time)
303	print "LOCAL TIMES.... ", time
304	### we define the vertical axis (or lack thereof) and cover possibilities for it to be altitude, pressure, geopotential. quite SPECIFIC.
305	if typefile in ['geo']: vert = [0.] ; stime = readslices(str(0))
306	else:
307	if vertmode is None: vertmode=0
308	if vertmode == 0:
309	if "vert" in nc.variables: vert = nc.variables["vert"][:]/1000. ; vertmode = 1
310	else: vert = np.arange(0,getattr(nc,vertdim),1)
311	elif vertmode == -1: vert = nc.variables["PHTOT"][0,:,0,0]/3.72 ; vert = np.array(vert[0:len(vert)-1]) #; print vert
312	elif vertmode == 1 or vertmode == 2: vert = nc.variables["vert"][:] ## pressure in Pa
313	else: vert = nc.variables["vert"][:]/1000. ## altitude in km
314	####################################################################
315	############ END of WE LOAD 4D DIMENSIONS : x, y, z, t #############
316	####################################################################
317
318	### we fill the arrays of varname, namefile, time, colorbar at the current step considered (NB: why use both k and nnn ?)
319	all_varname[k] = varname
320	all_namefile[k] = namefile
321	all_time[k] = time
322	all_vert[k] = vert
323	all_lat[k] = lat
324	all_lon[k] = lon
325	all_colorb[k] = clb[vvv]
326	if var2: all_var2[k], plot_x[k], plot_y[k] = select_getfield(zvarname=var2,znc=nc,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=all_lon[k],ylat=all_lat[k],yalt=all_vert[k],ytime=all_time[k],analysis=analysis)
327	if winds: all_windu[k] = getfield(nc,uchar) ; all_windv[k] = getfield(nc,vchar)
328	### we fill the arrays of fields to be plotted at the current step considered
329	all_var[k], plot_x[k], plot_y[k] = select_getfield(zvarname=all_varname[k],znc=nc,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=all_lon[k],ylat=all_lat[k],yalt=all_vert[k],ytime=all_time[k],analysis=analysis)
330
331	# last line of for namefile in namefiles
332	print "**** GOT SUBDATA:",k," NAMEFILE:",namefile," VAR:",varname, var2 ; k += 1 ; firstfile = False
333
334	### note to self : stopped comment rewriting here.
335
336	##################################
337	### Operation on files (I) with _var
338	if ope is not None and "var" in ope:
339	print "********** OPERATION: ",ope
340	if len(namefiles) > 1: errormess("for this operation... please set only one file !")
341	if len(var) > 2: errormess("not sure this works for more than 2 vars... please check.")
342	if "div_var" in ope: all_var[k] = all_var[k-2] / all_var[k-1] ; insert = '_div_'
343	elif "mul_var" in ope: all_var[k] = all_var[k-2] * all_var[k-1] ; insert = '_mul_'
344	elif "add_var" in ope: all_var[k] = all_var[k-2] + all_var[k-1] ; insert = '_add_'
345	elif "sub_var" in ope: all_var[k] = all_var[k-2] - all_var[k-1] ; insert = '_sub_'
346	else: errormess(ope+" : non-implemented operation. Check pp.py --help")
347	plot_x[k] = None ; plot_y[k] = None
348	all_time[k] = all_time[k-1] ; all_vert[k] = all_vert[k-1] ; all_lat[k] = all_lat[k-1] ; all_lon[k] = all_lon[k-1] ; all_namefile[k] = all_namefile[k-1]
349	all_varname[k] = all_varname[k-2] + insert + all_varname[k-1]
350	if len(clb) >= zelen: all_colorb[k] = clb[-1] # last additional user-defined color is for operation plot
351	else: all_colorb[k] = all_colorb[k-1] # if no additional user-defined color... set same as var
352	### only the operation plot. do not mention colorb so that it is user-defined?
353	if "only" in ope:
354	numplot = 1 ; all_var[0] = all_var[k]
355	all_time[0] = all_time[k] ; all_vert[0] = all_vert[k] ; all_lat[0] = all_lat[k] ; all_lon[0] = all_lon[k] ; all_namefile[0] = all_namefile[k]
356	all_varname[0] = all_varname[k-2] + insert + all_varname[k-1]
357
358
359	##################################
360	### Operation on files (II) without _var
361	# we re-iterate on the plots to set operation subplots to make it compatible with multifile (using the same ref file)
362	# (k+1)%3==0 is the index of operation plots
363	# (k+2)%3==0 is the index of reference plots
364	# (k+3)%3==0 is the index of first plots
365	opefirstpass=True
366	if ope is not None and "var" not in ope:
367	print "********** OPERATION: ",ope
368	for k in np.arange(zelen):
369	if len(var) > 1: errormess("for this operation... please set only one var !")
370	if ope in ["-","+","-%","-_only","+_only","-%_only","-_histo"]:
371	if fileref is None: errormess("fileref is missing!")
372	else:ncref = Dataset(fileref)
373
374	if opefirstpass: ## first plots
375	for ll in np.arange(len(namefiles)):
376	print "SETTING FIRST PLOT"
377	all_varname[3ll] = all_varname[ll] ; all_time[3ll] = all_time[ll] ; all_vert[3ll] = all_vert[ll] ; all_lat[3ll] = all_lat[ll] ; all_lon[3ll] = all_lon[ll] ; all_namefile[3ll] = all_namefile[ll] ; all_var2[3ll] = all_var2[ll] ; all_colorb[3ll] = all_colorb[ll] ; all_var[3*ll] = all_var[ll]
378	if plot_y[ll] is not None: plot_y[3ll] = plot_y[ll] ; plot_x[3ll] = plot_x[ll]
379	else: plot_y[3ll] = None ; plot_x[3ll] = None
380	opefirstpass=False
381
382	if (k+2)%3==0: ## reference plots
383	print "SETTING REFERENCE PLOT"
384	all_varname[k] = all_varname[k-1] ; all_time[k] = all_time[k-1] ; all_vert[k] = all_vert[k-1] ; all_lat[k] = all_lat[k-1] ; all_lon[k] = all_lon[k-1] ; all_namefile[k] = all_namefile[k-1] ; all_var2[k] = all_var2[k-1] ; all_colorb[k] = all_colorb[k-1]
385	all_var[k], plot_x[k], plot_y[k] = select_getfield(zvarname=all_varname[k-1],znc=ncref,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=all_lon[k],ylat=all_lat[k],yalt=all_vert[k],ytime=all_time[k],analysis=analysis)
386	if winds: all_windu[k] = getfield(ncref,uchar) ; all_windv[k] = getfield(ncref,vchar)
387
388	if (k+1)%3==0: ## operation plots
389	print "SETTING OPERATION PLOT"
390	all_varname[k] = all_varname[k-1] ; all_time[k] = all_time[k-1] ; all_vert[k] = all_vert[k-1] ; all_lat[k] = all_lat[k-1] ; all_lon[k] = all_lon[k-1] ; all_namefile[k] = all_namefile[k-1] ; all_var2[k] = all_var2[k-1]
391	if ope in ["-","-_only","-_histo"]:
392	all_var[k]= all_var[k-2] - all_var[k-1]
393	if plot_y[k-1] is not None and plot_y[k-2] is not None: plot_y[k] = plot_y[k-2] - plot_y[k-1]
394	if plot_y[k-2] is None: plot_y[k] = None; plot_x[k] = None
395	elif ope in ["+","+_only"]:
396	all_var[k]= all_var[k-2] + all_var[k-1]
397	if plot_y[k-1] is not None and plot_y[k-2] is not None: plot_y[k] = plot_y[k-2] + plot_y[k-1]
398	if plot_y[k-2] is None: plot_y[k] = None; plot_x[k] = None
399	elif ope in ["-%","-%_only"]:
400	masked = np.ma.masked_where(all_var[k-1] == 0,all_var[k-1])
401	masked.set_fill_value([np.NaN])
402	all_var[k]= 100.*(all_var[k-2] - masked)/masked
403	if plot_y[k-1] is not None and plot_y[k-2] is not None:
404	masked = np.ma.masked_where(plot_y[k-1] == 0,plot_y[k-1])
405	masked.set_fill_value([np.NaN])
406	plot_y[k]= 100.*(plot_y[k-2] - masked)/masked
407	if plot_y[k-2] is None: plot_y[k] = None; plot_x[k] = None
408	if len(clb) >= zelen: all_colorb[k] = clb[-1]
409	else: all_colorb[k] = "RdBu_r" # if no additional user-defined color... set a good default one
410	if winds: all_windu[k] = all_windu[k-2]-all_windu[k-1] ; all_windv[k] = all_windv[k-2] - all_windv[k-1]
411
412	elif ope in ["cat"]:
413	tabtime = all_time[0];tab = all_var[0];k = 1
414	if var2: tab2 = all_var2[0]
415	while k != len(namefiles) and len(all_time[k]) != 0:
416	if var2: tab2 = np.append(tab2,all_var2[k],axis=0)
417	tabtime = np.append(tabtime,all_time[k]) ; tab = np.append(tab,all_var[k],axis=0) ; k += 1
418	all_time[0] = np.array(tabtime) ; all_var[0] = np.array(tab) ; numplot = 1
419	if var2: all_var2[0] = np.array(tab2)
420	else: errormess(ope+" : non-implemented operation. Check pp.py --help")
421	if "only" in ope:
422	numplot = 1 ; all_var[0] = all_var[k]
423	all_time[0] = all_time[k] ; all_vert[0] = all_vert[k] ; all_lat[0] = all_lat[k] ; all_lon[0] = all_lon[k] ; all_namefile[0] = all_namefile[k] ; plot_x[0]=plot_x[k] ; plot_y[0]=plot_y[k]
424	all_varname[0] = all_varname[k]
425
426	##################################
427	### Open a figure and set subplots
428	fig = figure()
429	subv,subh = definesubplot( numplot, fig )
430	if ope in ['-','-%','-_histo'] and len(namefiles) ==1 : subv,subh = 2,2
431	elif ope in ['-','-%'] and len(namefiles)>1 : subv, subh = len(namefiles), 3
432
433	#################################
434	### Time loop for plotting device
435	nplot = 1 ; error = False ; firstpass = True
436	if lstyle is not None: linecycler = cycle(lstyle)
437	else: linecycler = cycle(["-","--","-.",":"])
438	print "********************************************"
439	while error is False:
440
441	print "********** PLOT", nplot, " OF ",numplot
442	if nplot > numplot: break
443
444	####################################################################
445	## get all indexes to be taken into account for this subplot and then reduce field
446	## 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
447	if ope is not None:
448	if fileref is not None: index_f = ((nplot-1)//(nlonnlatnvertntime))%(3len(namefiles)) ## OK only 1 var, see test in the beginning
449	elif "var" in ope: index_f = ((nplot-1)//(nlonnlatnvert*ntime))%(len(var)+1) ## OK only 1 file, see test in the beginning
450	elif "cat" in ope: index_f = 0
451	elif not firstpass:
452	if len(namefiles) > 1 and len(var) == 1 and which == "unidim": pass
453	else: yeah = len(namefiles)len(var) ; index_f = ((nplot-1)//(nlonnlatnvertntime))%yeah
454	else: yeah = len(namefiles)len(var) ; index_f = ((nplot-1)//(nlonnlatnvertntime))%yeah
455	time = all_time[index_f] ; vert = all_vert[index_f] ; lat = all_lat[index_f] ; lon = all_lon[index_f]
456	indexlon = getsindex(sslon,(nplot-1)%nlon,lon)
457	indexlat = getsindex(sslat,((nplot-1)//nlon)%nlat,lat)
458	indexvert = getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert)
459	plotx=plot_x[index_f] ; ploty=plot_y[index_f]
460	if mrate is not None: indextime = None
461	else: indextime = getsindex(stime,((nplot-1)//(nlonnlatnvert))%ntime,time)
462	ltst = None
463	if typefile in ['meso'] and indextime is not None and len(indextime) < 2: ltst = localtime ( indextime, slon , all_namefile[index_f])
464	print "********** INDEX LON:",indexlon," LAT:",indexlat," VERT:",indexvert," TIME:",indextime
465	##var = nc.variables["phisinit"][:,:]
466	##contourf(np.transpose(var),30,cmap = get_cmap(name="Greys_r") ) ; axis('off') ; plot(indexlat,indexlon,'mx',mew=4.0,ms=20.0)
467	##show()
468	##exit()
469	#truc = True
470	#truc = False
471	#if truc: indexvert = None
472	####################################################################
473	########## REDUCE FIELDS
474	####################################################################
475	error = False
476	varname = all_varname[index_f]
477	which=''
478	if varname: ### what is shaded.
479	what_I_plot, error = reducefield( all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \
480	yint=yintegral, alt=vert, anomaly=anomaly, redope=redope, mesharea=area, unidim=is1d)
481	if mult != 2718.: what_I_plot = what_I_plot*mult
482	else: what_I_plot = np.log10(what_I_plot) ; print "log plot"
483
484	if var2: ### what is contoured.
485	what_I_plot_contour, error = reducefield( all_var2[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert, \
486	yint=yintegral, alt=vert )
487	if winds: ### what is plot as vectors.
488	vecx, error = reducefield( all_windu[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert)
489	vecy, error = reducefield( all_windv[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert)
490	if varname in [uchar,vchar]: what_I_plot = np.sqrt( np.square(vecx) + np.square(vecy) ) ; varname = "wind"
491
492	if plotx is not None:
493	plotx, error = reducefield( plotx, d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \
494	yint=yintegral, alt=vert, anomaly=anomaly, redope=redope, mesharea=area, unidim=is1d)
495	ploty, error = reducefield( ploty, d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \
496	yint=yintegral, alt=vert, anomaly=anomaly, redope=redope, mesharea=area, unidim=is1d)
497	which='xy'
498	#####################################################################
499	#if truc:
500	# nx = what_I_plot.shape[2] ; ny = what_I_plot.shape[1] ; nz = what_I_plot.shape[0]
501	# for k in range(nz): print k,' over ',nz ; what_I_plot[k,:,:] = what_I_plot[k,:,:] / smooth(what_I_plot[k,:,:],12)
502	# for iii in range(nx):
503	# for jjj in range(ny):
504	# deviation = what_I_plot[:,jjj,iii] ; mx = max(deviation) ; mn = min(deviation)
505	# if iii > 6 and iii < nx-6 and jjj > 6 and jjj < ny-6: what_I_plot[0,jjj,iii],rel = singlet(deviation,vert/1000.) ### z must be in km
506	# else: what_I_plot[0,jjj,iii] = 0.
507	# if np.abs(what_I_plot[0,jjj,iii]) > 1.5:
508	# print iii,jjj,what_I_plot[0,jjj,iii],int(abs(1.-mx)100.),int(abs(1.-mn)100.)
509	# plot(rel)
510	# show()
511	# anomaly = True ### pour avoir les bons reglages plots
512	# what_I_plot = what_I_plot[0,:,:]
513	#####################################################################
514
515	####################################################################
516	### General plot settings
517	changesubplot = (numplot > 1) and (len(what_I_plot.shape) != 1) and (which != "xy") ## default for 1D plots: superimposed. to be reworked for better flexibility.
518	if changesubplot: subplot(subv,subh,nplot) #; subplots_adjust(wspace=0,hspace=0)
519	colorb = all_colorb[index_f]
520	####################################################################
521	if error:
522	errormess("There is an error in reducing field !")
523	else:
524	ticks = ndiv + 1
525	fvar = varname
526	if anomaly: fvar = 'anomaly'
527	###
528	if mapmode == 0: ### could this be moved inside imov loop ?
529	itime=indextime
530	if len(what_I_plot.shape) == 3: itime=[0]
531	m = None ; x = None ; y = None
532	latyeah = lat ; lonyeah = lon
533	if typefile in ['meso']:
534	# now that the section is determined we can set the real lat
535	# ... or for now, a temptative one.
536	if slon is not None: latyeah = lat2d[:,0]
537	if slat is not None: lonyeah = lon2d[0,:]
538	what_I_plot, x, y = define_axis(lonyeah,latyeah,vert,time,indexlon,indexlat,indexvert,\
539	itime,what_I_plot, len(all_var[index_f].shape),vertmode,redope)
540	###
541	if (fileref is not None) and ((index_f+1)%3 == 0): zevmin, zevmax = calculate_bounds(what_I_plot,vmin=minop,vmax=maxop)
542	else: zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
543	#if (fileref is not None) and (index_f == numplot-1): colorb = "RdBu_r"
544	if colorb in ["def","nobar","onebar"]: palette = get_cmap(name=defcolorb(fvar.upper()))
545	else: palette = get_cmap(name=colorb)
546	#palette = cm.GMT_split
547	##### 1. ELIMINATE 0D or >3D CASES
548	if len(what_I_plot.shape) == 0:
549	print "VALUE VALUE VALUE VALUE ::: ",what_I_plot
550	save = 'donothing'
551	elif len(what_I_plot.shape) >= 4:
552	print "WARNING!!! ",len(what_I_plot.shape),"-D PLOT NOT SUPPORTED !!! dimensions: ",what_I_plot.shape
553	errormess("Are you sure you did not forget to prescribe a dimension ?")
554	##### 2. HANDLE simple 1D/2D field and movies of 1D/2D fields
555	else:
556	if mrate is not None: iend=len(time)-1
557	else: iend=0
558	imov = 0
559	if analysis in ['density','histo','fft','histodensity']: which="xy"
560	elif len(what_I_plot.shape) == 3:
561	if var2 and which == '': which = "contour" ## have to start with contours rather than shading
562	elif which == '': which = "regular"
563	if mrate is None: errormess("3D field. Use --rate RATE for movie or specify --time TIME. Exit.")
564	elif len(what_I_plot.shape) == 2:
565	if var2 and which == '': which = "contour" ## have to start with contours rather than shading
566	elif which == '': which = "regular"
567	if mrate is not None and which == '': which = "unidim"
568	elif len(what_I_plot.shape) == 1 and which == '' :
569	which = "unidim"
570	if what_I_plot.shape[-1] == 1: print "VALUE VALUE VALUE VALUE ::: ", what_I_plot[0] ; save = 'donothing'
571	## if which == "unidim" and len(namefiles) > 1: numplot = 1 # this case is similar to several vars from one file
572	##### IMOV LOOP #### IMOV LOOP
573	while imov <= iend:
574	print "-> frame ",imov+1, which
575	if which == "regular":
576	if mrate is None: what_I_plot_frame = what_I_plot
577	else: what_I_plot_frame = what_I_plot[imov,:,:]
578	if winds:
579	if mrate is None: vecx_frame = vecx ; vecy_frame = vecy
580	else: vecx_frame = vecx[imov,:,:] ; vecy_frame = vecy[imov,:,:]
581	elif which == "contour":
582	if mrate is None or what_I_plot_contour.ndim < 3: what_I_plot_frame = what_I_plot_contour
583	else: what_I_plot_frame = what_I_plot_contour[imov,:,:]
584	elif which == "unidim":
585	if mrate is None: what_I_plot_frame = what_I_plot
586	else: what_I_plot_frame = what_I_plot[:,imov] ## because swapaxes
587	#if mrate is not None:
588	if mapmode == 1:
589	m = define_proj(proj,wlon,wlat,back=back,blat=blat,blon=blon) ## this is dirty, defined above but out of imov loop
590	x, y = m(lon2d, lat2d) ## this is dirty, defined above but out of imov loop
591	if typefile in ['meso'] and mapmode == 1: what_I_plot_frame = dumpbdy(what_I_plot_frame,6,condition=True)
592	# if typefile in ['mesoideal']: what_I_plot_frame = dumpbdy(what_I_plot_frame,0,stag='W',condition=dumped_vert_stag)
593
594	if which == "unidim":
595	# if lbls is not None: lbl=lbls[nplot-1]
596	if lbls is not None: lbl=lbls[index_f]
597	else:
598	lbl = ""
599	if indexlat is not None: lbl = lbl + " ix" + str(indexlat[0])
600	if indexlon is not None: lbl = lbl + " iy" + str(indexlon[0])
601	if indexvert is not None: lbl = lbl + " iz" + str(indexvert[0])
602	if indextime is not None: lbl = lbl + " it" + str(indextime[0])
603	if lbl == "": lbl = all_namefile[index_f]
604
605	if mrate is not None: x = y ## because swapaxes...
606	#what_I_plot_frame = np.diff(what_I_plot_frame, n=1) ; x = x[1:]
607
608	zeline = next(linecycler) ## "-" for simple lines
609	if tile: zemarker = 'x'
610	else: zemarker = None
611	this_is_a_regular_plot = (indexvert is not None) or (indextime is None) or (indexlat is None) or (indexlon is None)
612	if this_is_a_regular_plot: plot(x,what_I_plot_frame,zeline,label=lbl,marker=zemarker) ## vertical profile
613	else: plot(what_I_plot_frame,x,zeline,label=lbl,marker=zemarker) ## regular plot
614	if nplot > 1: legend(loc='best')
615	if indextime is None and axtime is not None and xlab is None: xlabel(axtime.upper()) ## define the right label
616	if save == 'txt': writeascii(np.transpose([x,np.transpose(what_I_plot)]),'profile'+str(nplot*1000+imov)+'.txt')
617
618	elif which == "xy":
619	if lbls is not None: lbl=lbls[index_f]
620	else: lbl=None
621	if analysis is not None:
622	if analysis == 'histo':
623	if zelen == 1:
624	mpl.pyplot.hist(ploty.flatten(),bins=ndiv,normed=True, alpha=0.67, facecolor = 'green', label = lbls)
625	mpl.pyplot.legend()
626	mpl.pyplot.grid(True)
627	mpl.pyplot.title(zetitle)
628	else:
629	multiplot[index_f]=ploty.flatten()
630	if index_f == zelen-1:
631	if ope is not None: multiplot = getVar(multiplot,3*np.arange((index_f+1)/3)+2) ## we only compute histograms for the operation plots
632	mpl.pyplot.hist(multiplot,bins=ndiv,normed=True, alpha=0.75, label = lbls)
633	mpl.pyplot.legend()
634	mpl.pyplot.grid(True)
635	mpl.pyplot.title(zetitle)
636	elif analysis in ['density','histodensity']:
637	if ope is not None and (index_f+1)%3 !=0: pass
638	else:
639	plotx = np.linspace(min(ploty.flatten()),max(ploty.flatten()),1000)
640	density = gaussian_kde(ploty.flatten())
641	# density.covariance_factor = lambda : .25 # adjust the covariance factor to change the bandwidth if needed
642	# density._compute_covariance()
643	# display the mean and variance of the kde:
644	sample = density.resample(size=20000)
645	n, (smin, smax), sm, sv, ss, sk = describe(sample[0])
646	mpl.pyplot.plot(plotx,density(plotx), label = lbl+'\nmean: '+str(sm)[0:5]+' std: '+str(np.sqrt(sv))[0:5]+'\nskewness: '+str(ss)[0:5]+' kurtosis: '+str(sk)[0:5])
647	if analysis == 'histodensity': # plot both histo and density (to assess the rightness of the kernel density estimate for exemple) and display the estimated variance
648	mpl.pyplot.hist(ploty.flatten(),bins=ndiv,normed=True, alpha=0.30, label = lbl)
649	if index_f == zelen-1: mpl.pyplot.legend() ; mpl.pyplot.title(zetitle)
650	else:
651	plot(plotx,ploty,label = lbl)
652	if index_f == zelen-1: mpl.pyplot.legend() ; mpl.pyplot.title(zetitle)
653	mpl.pyplot.grid(True)
654	else:
655	plot(plotx,ploty,label = lbl)
656	if index_f == zelen-1: mpl.pyplot.legend() ; mpl.pyplot.title(zetitle)
657	if varname == 'hodograph':
658	a=0
659	for ii in np.arange(len(time)):
660	if a%6 == 0: mpl.pyplot.text(plotx[ii],ploty[ii],time[ii])
661	a=a+1
662	mpl.pyplot.grid(True)
663
664	elif which == "regular":
665
666	# plot stream lines if there is a stream file and a vert/lat slice. Might not work with movies ??
667	if streamflag and sslat is None and svert is None:
668	streamfile = all_namefile[index_f].replace('.nc','_stream.nc')
669	teststream = os.path.exists(streamfile)
670	if teststream:
671	print 'INFO: Using stream file',streamfile, 'for stream lines'
672	ncstream = Dataset(streamfile)
673	psi = getfield(ncstream,'psi')
674	psi = psi[0,:,:,0] # time and longitude are dummy dimensions
675	if psi.shape[1] != len(x) or psi.shape[0] != len(y):
676	errormess('stream file does not fit! Dimensions: '+str(psi.shape)+' '+str(x.shape)+' '+str(y.shape))
677	zelevels = np.arange(-1.e10,1.e10,1.e9)
678	zemin = np.min(abs(zelevels))
679	zemax = np.max(abs(zelevels))
680	zewidth = (abs(zelevels)-zemin)*(5.- 0.5)/(zemax - zemin) + 0.5 # linewidth ranges from 5 to 0.5
681	cs = contour( x,y,psi, zelevels, colors='k', linewidths = zewidth)
682	clabel(cs, inline=True, fontsize = 4.*rcParams['font.size']/5., fmt="%1.1e")
683	else:
684	print 'WARNING: STREAM FILE',streamfile, 'DOES NOT EXIST !'
685
686	if hole: what_I_plot_frame = hole_bounds(what_I_plot_frame,zevmin,zevmax)
687	else: what_I_plot_frame = bounds(what_I_plot_frame,zevmin,zevmax)
688	if flagnolow: what_I_plot_frame = nolow(what_I_plot_frame)
689	if not tile:
690	#zelevels = np.linspace(zevmin(1. + 1.e-7),zevmax(1. - 1.e-7)) #,num=20)
691	zelevels = np.linspace(zevmin,zevmax,num=ticks)
692	#what_I_plot_frame = smooth(what_I_plot_frame,100)
693	if mapmode == 1: m.contourf( x, y, what_I_plot_frame, zelevels, cmap = palette, alpha=trans)
694	elif mapmode == 0: contourf( x, y, what_I_plot_frame, zelevels, cmap = palette, alpha=trans)
695	if cross is not None and mapmode == 1:
696	idx,idy=m(cross[0][0],cross[0][1])
697	mpl.pyplot.plot([idx],[idy],'+k',mew=0.5,ms=18)
698	else:
699	if mapmode == 1: m.pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax, alpha=trans)
700	elif mapmode == 0: pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax, alpha=trans)
701
702
703	if colorb not in ['nobar','onebar']:
704	if (fileref is not None) and ((index_f+1)%3 == 0): daformat = "%.3f"
705	elif mult != 1: daformat = "%.1f"
706	else: daformat = fmtvar(fvar.upper())
707	#if proj in ['moll','cyl']: zeorientation="horizontal" ; zepad = 0.05
708	if proj in ['moll']: zeorientation="horizontal" ; zepad = 0.05
709	else: zeorientation="vertical" ; zepad = 0.03
710	zecb = colorbar( fraction=0.05,pad=zepad,format=daformat,orientation=zeorientation,\
711	ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,21])),extend='neither',spacing='proportional' )
712	if zeorientation == "horizontal" and zetitle[0] != "fill": zecb.ax.set_xlabel(zetitle[index_f]) ; zetitle[index_f]=""
713	if winds:
714	if typefile in ['meso']:
715	[vecx_frame,vecy_frame] = [dumpbdy(vecx_frame,6,stag=uchar,condition=True), dumpbdy(vecy_frame,6,stag=vchar,condition=True)]
716	key = True
717	if fvar in ['UV','uv','uvmet']: key = False
718	elif typefile in ['gcm']:
719	key = False
720	if metwind and mapmode == 1: [vecx_frame,vecy_frame] = m.rotate_vector(vecx_frame, vecy_frame, lon2d, lat2d)
721	if trans != 0.0: colorvec = definecolorvec(colorb)
722	else: colorvec = definecolorvec(back)
723	vectorfield(vecx_frame, vecy_frame, x, y, stride=stride, csmooth=2,\
724	#scale=15., factor=300., color=colorvec, key=key)
725	scale=20., factor=250./facwind, color=colorvec, key=key)
726	#200. ## or csmooth=stride
727	### THIS IS A QUITE SPECIFIC PIECE (does not work for mesoscale files)
728	if ope == '-_histo' and nplot == numplot: # this should work as long as ope is '-' guarantees 3 plots for 4 panels without contour
729	subplot(subv,subh,nplot+1)
730	rcParams["legend.fontsize"] = 'xx-large'
731	if indexlat is None:
732	latmin = -50.; latmax = 50. # latitude range for histogram of difference
733	zeindexlat = (lat<latmax)*(lat>latmin)
734	if typefile in ['meso']: zeindexlat = 10
735	# this follows the define_axis logic in myplot.py:
736	if indextime is None or indexlon is None: what_I_plot_frame = what_I_plot_frame[zeindexlat,:]
737	else: what_I_plot_frame = what_I_plot_frame[:,zeindexlat]
738	toplot = np.ravel(what_I_plot_frame[np.isnan(what_I_plot_frame)==False])
739	zebins = np.linspace(zevmin,zevmax,num=30)
740	hist(toplot,bins=zebins,histtype='step',linewidth=2,color='k',normed=True)
741	zestd = np.std(toplot);zemean = np.mean(toplot)
742	zebins = np.linspace(zevmin,zevmax,num=300)
743	zegauss = (1./(zestd * np.sqrt(2 * np.pi)) * np.exp( - (zebins - zemean)*2 / (2 zestd**2) ) )
744	plot(zebins, zegauss, linewidth=1, color='r',label="mean: "+str(zemean)[0:5]+"\nstd: "+str(zestd)[0:5])
745	legend(loc=0,frameon=False)
746	subplot(subv,subh,nplot) # go back to last plot for title of contour difference
747	if ope is not None and "only" not in ope: title("fig(1) "+ope+" fig(2)")
748	elif ope is not None and "only" in ope: title("fig(1) "+ope[0]+" fig(2)")
749
750	elif which == "contour":
751	zevminc, zevmaxc = calculate_bounds(what_I_plot_frame, vmin=min(what_I_plot_frame), vmax=max(what_I_plot_frame))
752	zelevels = np.linspace(zevminc,zevmaxc,ticks/2) #20)
753	### another dirty specific stuff in the wall
754	if var2 == 'HGT': zelevels = np.arange(-10000.,30000.,500.) #1000.)
755	elif var2 == 'tpot': zelevels = np.arange(270,370,5)
756	elif var2 == 'tk': zelevels = np.arange(150,250,5)
757	elif var2 == 'wstar': zelevels = np.arange(0,10,1.0)
758	elif var2 == 'zmax_th': zelevels = np.arange(0,10,2.0) ; what_I_plot_frame = what_I_plot_frame / 1000.
759	###
760	if mapmode == 0:
761	what_I_plot_frame, x, y = define_axis( lonyeah,latyeah,vert,time,indexlon,indexlat,indexvert,\
762	itime,what_I_plot_frame, len(all_var2[index_f].shape),vertmode,redope)
763	## is this needed? only if len(all_var2[index_f].shape) != len(all_var[index_f].shape)
764	cs = contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 0.33)#, alpha=0.5, linestyles='solid')
765	#cs = contour( x,y,what_I_plot_frame, zelevels, colors='w', linewidths = 0.5)#, alpha=0.5, linestyles='solid')
766	clabel(cs, inline=1, fontsize = 4.*rcParams['font.size']/5., fmt=fmtvar(var2.upper()))
767	elif mapmode == 1:
768	cs = m.contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 0.33)#, alpha=0.5, linestyles='solid')
769	#clabel(cs, inline=0, fontsize = rcParams['font.size'], fmt="%.0f") #fmtvar(var2.upper()))
770	if which in ["regular","unidim","xy"]:
771
772	if nplot > 1 and which in ["unidim","xy"]:
773	pass ## because we superimpose nplot instances
774	else:
775	# Axis directives for movie frames [including the first one).
776	zxmin, zxmax = xaxis ; zymin, zymax = yaxis
777	if zxmin is not None: mpl.pyplot.xlim(xmin=zxmin)
778	if zxmax is not None: mpl.pyplot.xlim(xmax=zxmax)
779	if zymin is not None: mpl.pyplot.ylim(ymin=zymin)
780	if zymax is not None: mpl.pyplot.ylim(ymax=zymax)
781	if ylog and not xlog: mpl.pyplot.semilogy()
782	if xlog and not ylog: mpl.pyplot.semilogx()
783	if xlog and ylog:
784	mpl.pyplot.xscale('log')
785	mpl.pyplot.yscale('log')
786	if invert_y: ax = mpl.pyplot.gca() ; ax.set_ylim(ax.get_ylim()[::-1])
787	if xlab is not None: xlabel(xlab)
788	if ylab is not None: ylabel(ylab)
789
790	if mrate is not None:
791	### THIS IS A MENCODER MOVIE
792	if mrate > 0:
793	figframe=mpl.pyplot.gcf()
794	if mquality: figframe.set_dpi(600.)
795	else: figframe.set_dpi(200.)
796	mframe=fig2img(figframe)
797	if imov == 0:
798	moviename='movie' ;W,H = figframe.canvas.get_width_height()
799	video = VideoSink((H,W), moviename, rate=mrate, byteorder="rgba")
800	video.run(mframe) ; close()
801	if imov == iend: video.close()
802	### THIS IS A WEBPAGE MOVIE
803	else:
804	nameframe = "image"+str(1000+imov)
805	makeplotres(nameframe,res=100.,disp=False) ; close()
806	if imov == 0: myfile = open("zepics", 'w')
807	myfile.write("modImages["+str(imov)+"] = '"+nameframe+"_100.png';"+ '\n')
808	if imov == iend:
809	myfile.write("first_image = 0;"+ '\n')
810	myfile.write("last_image = "+str(iend)+";"+ '\n')
811	myfile.close()
812	if var2 and which == "regular": which = "contour"
813	imov = imov+1
814	elif which == "contour":
815	which = "regular"
816
817	### Next subplot
818	zevarname = varname
819	if redope is not None: zevarname = zevarname + "_" + redope
820	basename = getname(var=zevarname,var2=var2,winds=winds,anomaly=anomaly)
821	if len(what_I_plot.shape) > 3:
822	basename = basename + getstralt(nc,level)
823	if mrate is not None: basename = "movie_" + basename
824	if typefile in ['meso']:
825	if sslon is not None: basename = basename + "_lon_" + str(int(round(lonp)))
826	if sslat is not None: basename = basename + "_lat_" + str(int(round(latp)))
827	plottitle = basename
828	### dans le nouveau systeme time=ls,sol,lt cette ligne pourrait ne servir a rien (ou deplacer au dessus)
829	if addchar and indextime is not None: [addchar,gogol,gogol2] = getlschar ( all_namefile[index_f] ) ; plottitle = plottitle + addchar
830	### en fait redope is None doit etre remplace par : n'est ni maxt ni mint
831	if redope is None and ltst is not None and ( (mapmode == 0) or (proj in ["lcc","laea","merc","nsper"]) ): plottitle = plottitle + "_LT" + str(ltst)
832	else:
833	if fileref is not None:
834	if (index_f+1)%3 == 0: plottitle = basename+' '+"fig(1) "+ope+" fig(2)"
835	elif (index_f+2)%3 == 0: plottitle = basename+' '+fileref
836	else: plottitle = basename+' '+all_namefile[index_f]
837	else: plottitle = basename+' '+all_namefile[index_f]
838	if mult != 1: plottitle = '{:.0e}'.format(mult) + "*" + plottitle
839	if zetitle[0] != "fill":
840	if index_f < len(zetitle): plottitle = zetitle[index_f]
841	else: plottitle = zetitle[0]
842	if titleref is "fill": titleref=zetitle[index_f]
843	if fileref is not None and which != "xy":
844	if (index_f+2)%3 == 0: plottitle = titleref
845	if (index_f+1)%3 == 0: plottitle = "fig(1) "+ope+" fig(2)"
846	# if indexlon is not None: plottitle = plottitle + " lon: " + str(min(lon[indexlon])) +" "+ str(max(lon[indexlon]))
847	# if indexlat is not None: plottitle = plottitle + " lat: " + str(min(lat[indexlat])) +" "+ str(max(lat[indexlat]))
848	# if indexvert is not None: plottitle = plottitle + " vert: " + str(min(vert[indexvert])) +" "+ str(max(vert[indexvert]))
849	# if indextime is not None: plottitle = plottitle + " time: " + str(min(time[indextime])) +" "+ str(max(time[indextime]))
850	if colorb != "onebar": title( plottitle )
851	if nplot >= numplot: error = True
852	nplot += 1
853
854	if len(namefiles) > 1 and len(var) == 1 and which == "unidim": index_f=index_f+1
855	firstpass=False
856
857	if colorb == "onebar":
858	cax = axes([0.1, 0.2, 0.8, 0.03]) # a ameliorer
859	zecb = colorbar(cax=cax, orientation="horizontal", format=fmtvar(fvar.upper()),\
860	ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,21])),extend='neither',spacing='proportional')
861	if zetitle[0] != "fill": zecb.ax.set_xlabel(zetitle[index_f]) ; zetitle[index_f]=""
862
863
864	##########################################################################
865	### Save the figure in a file in the data folder or an user-defined folder
866	if outputname is None:
867	if typefile in ['meso']: prefix = getprefix(nc)
868	elif typefile in ['gcm']: prefix = 'LMD_GCM_'
869	else: prefix = ''
870	###
871	zeplot = prefix + basename
872	if zoom: zeplot = zeplot + "zoom"+str(abs(zoom))
873	if addchar: zeplot = zeplot + addchar
874	if numplot <= 0: zeplot = zeplot + "_LT"+str(abs(numplot))
875	###
876	if not target: zeplot = namefile[0:find(namefile,'wrfout')] + zeplot
877	else: zeplot = target + "/" + zeplot
878	###
879	else:
880	zeplot=outputname
881
882	if mrate is None:
883	pad_inches_value = 0.35
884	if wlon[1]-wlon[0] < 2.: pad_inches_value = 0.5 # LOCAL MODE (small values)
885	print "********** SAVE ", save
886	if save == 'png':
887	if display: makeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value) #,erase=True) ## a miniature
888	makeplotres(zeplot,res=resolution,pad_inches_value=pad_inches_value,disp=False)
889	elif save in ['eps','svg','pdf']: makeplotres(zeplot,pad_inches_value=pad_inches_value,disp=False,ext=save)
890	elif save == 'gui': show()
891	elif save == 'donothing': pass
892	elif save == 'txt': print "Saved results in txt file."
893	else:
894	print "INFO: save mode not supported. using gui instead."
895	show()
896
897	###################################
898	#### Getting more out of this video -- PROBLEMS WITH CREATED VIDEOS
899	#
900	#if mrate is not None:
901	# print "Re-encoding movie.. first pass"
902	# video.first_pass(filename=moviename,quality=mquality,rate=mrate)
903	# print "Re-encoding movie.. second pass"
904	# video.second_pass(filename=moviename,quality=mquality,rate=mrate)
905
906	###############
907	### Now the end
908	return zeplot

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

Download in other formats: