Context Navigation

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

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

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

Download in other formats: