Context Navigation

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

Last change on this file since 642 was 642, checked in by tnavarro, 13 years ago
Correction for histogram of difference. Write both variable and axis for 1D txt output - useful for gnuplot or sparse.py
Property svn:executable set to ``*
File size: 46.0 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	def planetoplot (namefiles,\
14	level=0,\
15	vertmode=0,\
16	proj=None,\
17	back=None,\
18	target=None,
19	stride=3,\
20	var=None,\
21	colorb="def",\
22	winds=False,\
23	addchar=None,\
24	interv=[0,1],\
25	vmin=None,\
26	vmax=None,\
27	tile=False,\
28	zoom=None,\
29	display=True,\
30	hole=False,\
31	save="gui",\
32	anomaly=False,\
33	var2=None,\
34	ndiv=10,\
35	mult=1.,\
36	zetitle=["fill"],\
37	slon=None,\
38	slat=None,\
39	svert=None,\
40	stime=None,\
41	outputname=None,\
42	resolution=200,\
43	ope=None,\
44	fileref=None,\
45	minop=0.,\
46	maxop=0.,\
47	titleref="fill",\
48	invert_y=False,\
49	xaxis=[None,None],\
50	yaxis=[None,None],\
51	ylog=False,\
52	yintegral=False,\
53	blat=None,\
54	blon=None,\
55	tsat=False,\
56	flagnolow=False,\
57	mrate=None,\
58	mquality=False,\
59	trans=1,\
60	zarea=None,\
61	axtime=None,\
62	redope=None,\
63	seevar=False,\
64	xlab=None,\
65	ylab=None,\
66	lbls=None,\
67	lstyle=None,\
68	cross=None):
69
70	####################################################################################################################
71	### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
72
73	#################################
74	### Load librairies and functions
75	from netCDF4 import Dataset
76	from myplot import getcoord2d,define_proj,makeplotres,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy,\
77	fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\
78	zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow,\
79	getname,localtime,check_localtime,polarinterv,getsindex,define_axis,determineplot,readslices,bidimfind,getlschar,hole_bounds,\
80	windamplitude,slopeamplitude,deltat0t1
81	from mymath import deg,max,min,mean,get_tsat,writeascii,fig2data,fig2img
82	import matplotlib as mpl
83	from matplotlib.pyplot import contour,contourf,hist, text,subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel, title, close, legend, xlabel, axis, ylabel, subplots_adjust, axes, clabel
84	from matplotlib.cm import get_cmap
85	#from mpl_toolkits.basemap import cm
86	import numpy as np
87	from numpy.core.defchararray import find
88	from videosink import VideoSink
89	from timestuff import sol2ls
90	import subprocess
91	#from singlet import singlet
92	from itertools import cycle
93
94	#########################
95	### PRELIMINARY STUFF ###
96	#########################
97	### we say hello
98	print "********************************************"
99	print "******** WELCOME TO PLANETOPLOT ********"
100	print "********************************************"
101	### we ensure namefiles and var are arrays
102	if not isinstance(namefiles, np.ndarray): namefiles = [namefiles]
103	if not isinstance(var, np.ndarray): var = [var]
104	### we initialize a few variables
105	initime=-1 ; sslon = None ; sslat = None
106	if slon is not None: sslon = np.zeros([1,2])
107	if slat is not None: sslat = np.zeros([1,2])
108	k = 0 ; firstfile = True ; count = 0
109	### we avoid specific cases not yet implemented
110	if mrate is not None and len(var) > 1: errormess("multivar not allowed in movies. should be fixed soon!")
111	### we prepare number of plot fields [zelen] and plot instances [numplot] according to user choices
112	### --> we support multifile and multivar plots : files + vars separated by commas are on the same figure
113	nlon, nlat, nvert, ntime, mapmode, nslices = determineplot(slon, slat, svert, stime)
114	zelen = len(namefiles)*len(var)
115	numplot = zelen*nslices
116	print "********** FILES, SLICES, VARS, TOTAL PLOTS: ", len(namefiles), nslices, len(var), numplot
117	print "********** MAPMODE: ", mapmode
118	### we correct number of plot fields for possible operation (substract, etc...)
119	if ope is not None:
120	if fileref is not None: zelen = zelen + 2
121	elif "var" in ope: zelen = zelen + 1
122	### we define the arrays for plot fields -- which will be placed within multiplot loops
123	all_var = [[]]zelen ; all_var2 = [[]]zelen
124	all_title = [[]]zelen ; all_varname = [[]]zelen ; all_namefile = [[]]zelen ; all_time = [[]]zelen
125	all_windu = [[]]zelen ; all_windv = [[]]zelen
126
127	#############################
128	### LOOP OVER PLOT FIELDS ###
129	#############################
130	for nnn in range(len(namefiles)):
131	for vvv in range(len(var)):
132
133	### we load each NETCDF objects in namefiles
134	namefile = namefiles[nnn]
135	nc = Dataset(namefile)
136	### we explore the variables in the file
137	varinfile = nc.variables.keys()
138	if seevar: print varinfile ; exit()
139	### we define the type of file we have (gcm, meso, etc...)
140	typefile = whatkindfile(nc)
141	if firstfile: typefile0 = typefile
142	elif typefile != typefile0: errormess("Not the same kind of files !", [typefile0, typefile])
143	### we care for input file being 1D simulations
144	if "longitude" in nc.dimensions and "latitude" in nc.dimensions: is1d = len(nc.variables["longitude"][:])*len(nc.variables["latitude"][:])
145	elif "lon" in nc.dimensions and "lat" in nc.dimensions: is1d = len(nc.variables["lon"][:])*len(nc.variables["lat"][:])
146	if typefile in ['gcm','earthgcm'] and is1d == 1: mapmode=0 ; winds=False
147	### we create default vert and time prescriptions if not here in case mapping mode is on (lat/lon maps)
148	if redope is None and mapmode == 1:
149	if svert is None: svert = readslices(str(level)) ; nvert=1
150	if stime is None and mrate is None:
151	stime = readslices(str(0)) ; ntime=1 ## this is a default choice
152	print "WELL... nothing about time axis. I took default: first time reference stored in file."
153	### we get the names of variables to be read. in case only one available, we choose this one.
154	varname=var[vvv]
155	if ((varname not in nc.variables) and not ((typefile in ['meso']) and (varname in ['UV','uv','uvmet','slopexy','SLOPEXY','deltat','DELTAT']))):
156	if len(varinfile) == 1: varname = varinfile[0]
157	else: varname = False
158	### we get the names of wind variables to be read (if any)
159	if winds:
160	[uchar,vchar,metwind] = getwinddef(nc)
161	if uchar == 'not found': winds = False
162	### we tell the user that either no var or no wind is not acceptable
163	if not varname and not winds: errormess("please set at least winds or var",printvar=nc.variables)
164	### we get the coordinates lat/lon to be used
165	[lon2d,lat2d] = getcoorddef(nc)
166	### we get an adapted map projection if none is provided by the user
167	if proj == None: proj = getproj(nc)
168	### we define plot boundaries given projection or user choices
169	if firstfile:
170	if proj in ["npstere","spstere"]: [wlon,wlat] = polarinterv(lon2d,lat2d)
171	elif proj in ["lcc","laea"]: [wlon,wlat] = wrfinterv(lon2d,lat2d)
172	else: [wlon,wlat] = simplinterv(lon2d,lat2d)
173	if zoom: [wlon,wlat] = zoomset(wlon,wlat,zoom)
174	elif zarea is not None: [wlon,wlat] = latinterv(area=zarea)
175
176	##########################################################
177	############ LOAD 4D DIMENSIONS : x, y, z, t #############
178	##########################################################
179	if typefile in ["gcm","earthgcm","ecmwf"]:
180	if slon is not None: sslon = slon
181	if slat is not None: sslat = slat
182	### LAT
183	if "latitude" in nc.variables: lat = nc.variables["latitude"][:]
184	elif "lat" in nc.variables: lat = nc.variables["lat"][:]
185	### LON
186	if "longitude" in nc.variables: lon = nc.variables["longitude"][:]
187	elif "lon" in nc.variables: lon = nc.variables["lon"][:]
188	### ALT
189	if "altitude" in nc.variables: vert = nc.variables["altitude"][:]
190	elif "Alt" in nc.variables: vert = nc.variables["Alt"][:]
191	elif "lev" in nc.variables: vert = nc.variables["lev"][:]
192	else: vert = [0.]
193	### AIRE (to weight means with the area)
194	if "aire" in nc.variables: area = nc.variables["aire"][:,:]
195	else: area = None
196	### TIME
197	if "Time" in nc.variables: time = nc.variables["Time"][:]
198	elif "time" in nc.variables: time = nc.variables["time"][:]
199	elif "time_counter" in nc.variables: time = nc.variables["time_counter"][:]/86400. #### time counter cinverstion from s-> days
200	else: time = [0.] #errormess("no time axis found.")
201	if axtime in ["ls"]:
202	print "converting to Ls ..."
203	for iii in range(len(time)):
204	time[iii] = sol2ls(time[iii]) # to use with timestuff
205	if iii > 0:
206	while abs(time[iii]-time[iii-1]) > 300:
207	time[iii] = time[iii]+360
208	# for 1D plots (no need for longitude computation):
209	if axtime in ["lt"]:
210	if initime == -1: initime=input("Please type initial local time:")
211	time = (initime+time*24)%24
212	print "LOCAL TIMES.... ", time
213	elif typefile in ['meso','geo']:
214	area = None ## not active for the moment
215	###### STUFF TO GET THE CORRECT LAT/LON FROM MESOSCALE FILES WITH 2D LAT/LON ARRAYS
216	###### principle: calculate correct indices then repopulate slon and slat
217	if slon is not None or slat is not None:
218	if firstfile and save == 'png' and typefile == 'meso' and "HGT" in varinfile and display: iwantawhereplot = nc #show a topo map with a cross on the chosen point
219	else: iwantawhereplot = None #do not show anything, just select indices
220	numlon = 1 ; numlat = 1
221	if slon is not None: numlon = slon.shape[1]
222	if slat is not None: numlat = slat.shape[1]
223	indices = np.ones([numlon,numlat,2]) ; vlon = None ; vlat = None
224	for iii in range(numlon):
225	for jjj in range(numlat):
226	if slon is not None: vlon = slon[0][iii] ### note: slon[:][0] does not work
227	if slat is not None: vlat = slat[0][jjj] ### note: slon[:][0] does not work
228	indices[iii,jjj,:] = bidimfind(lon2d,lat2d,vlon,vlat,file=iwantawhereplot)
229	lonp,latp = ( lon2d[indices[iii,jjj,0],indices[iii,jjj,1]] , lat2d[indices[iii,jjj,0],indices[iii,jjj,1]] )
230	### possible bug here if several --lat
231	for iii in range(numlon):
232	for jjj in range(numlat):
233	if slon is not None: sslon[0][iii] = indices[iii,0,1] #...this is idx
234	if slat is not None: sslat[0][jjj] = indices[0,jjj,0] #...this is idy
235	lonp,latp = ( lon2d[indices[0,0,0],indices[0,0,1]] , lat2d[indices[0,0,0],indices[0,0,1]] )
236	######
237	if typefile in ['meso'] and mapmode == 1: lon2d = dumpbdy(lon2d,6) ; lat2d = dumpbdy(lat2d,6) ### important to do that now and not before
238	######
239	if varname in ['PHTOT','W']: vertdim='BOTTOM-TOP_PATCH_END_STAG'
240	else: vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
241	if (var2 is not None and var2 not in ['PHTOT','W']): dumped_vert_stag=True ; vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
242	else: dumped_vert_stag=False
243	if varname in ['V']: latdim='SOUTH-NORTH_PATCH_END_STAG'
244	else: latdim='SOUTH-NORTH_PATCH_END_UNSTAG'
245	if varname in ['U']: londim='WEST-EAST_PATCH_END_STAG'
246	else: londim='WEST-EAST_PATCH_END_UNSTAG'
247	lon = np.arange(0,getattr(nc,londim),1) ; lat = np.arange(0,getattr(nc,latdim),1)
248	###
249	if axtime in ["ls","sol"]:
250	lstab, soltab, lttab = getlschar ( namefile, getaxis = True )
251	if axtime == "ls": time = lstab
252	elif axtime == "sol": time = soltab
253	else:
254	if "Times" in nc.dimensions: time = count + np.arange(0,len(nc.dimensions["Times"]),1)
255	elif "Time" in nc.dimensions: time = count + np.arange(0,len(nc.dimensions["Time"]),1)
256	else: time = count + np.arange(0,1,1)
257	if ope in ["cat"]:
258	if nnn > 0: count = time[-1] + 1 ## so that a cat is possible with simple subscripts
259	else: count = 0
260	else: count = 0
261	if axtime in ["lt"]:
262	ftime = np.zeros(len(time))
263	for i in range(len(time)):
264	ftime[i] = localtime ( interv[0]+time[i]interv[1], 0.5(wlon[0]+wlon[1]) )
265	time=ftime
266	time=check_localtime(time)
267	print "LOCAL TIMES.... ", time
268
269	###
270	if typefile in ['geo']: vert = [0.] ; stime = readslices(str(0))
271	else:
272	if vertmode is None: vertmode=0
273	if vertmode == 0: vert = np.arange(0,getattr(nc,vertdim),1)
274	elif vertmode == -1: vert = nc.variables["PHTOT"][0,:,0,0]/3.72 ; vert = np.array(vert[0:len(vert)-1]) #; print vert
275	elif vertmode == 1 or vertmode == 2: vert = nc.variables["vert"][:] ## pressure in Pa
276	else: vert = nc.variables["vert"][:]/1000. ## altitude in km
277	#if firstfile:
278	# lat0 = lat
279	#elif len(lat0) != len(lat):
280	# errormess("Not the same latitude lengths !", [len(lat0), len(lat)])
281	#elif sum((lat == lat0) == False) != 0:
282	# errormess("Not the same latitudes !", [lat,lat0])
283	## Faire d'autre checks sur les compatibilites entre fichiers!!
284	##########################################################
285	##########################################################
286	##########################################################
287
288	all_varname[k] = varname
289	all_namefile[k] = namefile
290	all_time[k] = time
291	if var2:
292	all_var2[k] = getfield(nc,var2)
293	if ((var2 in ['slopexy','SLOPEXY']) and (typefile in ['meso']) and (var2 not in nc.variables)): all_var2[k] = slopeamplitude(nc)
294	if winds: all_windu[k] = getfield(nc,uchar) ; all_windv[k] = getfield(nc,vchar)
295	#####################
296	##### GETFIELDS #####
297	#####################
298	##### SPECIFIC CASES
299	##### 1. saturation temperature
300	if varname in ["temp","t","T_nadir_nit","T_nadir_day","temp_day","temp_night"] and tsat:
301	tt=getfield(nc,varname) ; print "computing Tsat-T, I ASSUME Z-AXIS IS PRESSURE"
302	if type(tt).__name__=='MaskedArray': tt.set_fill_value([np.NaN]) ; tinput=tt.filled()
303	else: tinput=tt
304	all_var[k]=get_tsat(vert,tinput,zlon=lon,zlat=lat,zalt=vert,ztime=time)
305	##### 2. wind amplitude
306	elif ((varname in ['UV','uv','uvmet']) and (typefile in ['meso']) and (varname not in nc.variables)):
307	all_var[k]=windamplitude(nc)
308	elif ((varname in ['slopexy','SLOPEXY']) and (typefile in ['meso']) and (varname not in nc.variables)):
309	all_var[k]=slopeamplitude(nc)
310	elif ((varname in ['DELTAT','deltat']) and (typefile in ['meso']) and (varname not in nc.variables)):
311	all_var[k]=deltat0t1(nc)
312	else:
313	##### GENERAL STUFF HERE
314	all_var[k] = getfield(nc,varname)
315	print "**** GOT SUBDATA:",k," NAMEFILE:",namefile," VAR:",varname, var2 ; k += 1 ; firstfile = False
316	#### End of for namefile in namefiles
317
318	##################################
319	### Operation on files
320	if ope is not None:
321	print "********** OPERATION: ",ope
322	if "var" not in ope:
323	if len(var) > 1: errormess("for this operation... please set only one var !")
324	if ope in ["-","+","-%"]:
325	if fileref is not None:
326	ncref = Dataset(fileref)
327	if ((all_varname[k-1] in ['UV','uv','uvmet']) and (typefile in ['meso']) and (all_varname[k-1] not in ncref.variables)):
328	all_var[k] = windamplitude(ncref)
329	elif ((all_varname[k-1] in ['slopexy','SLOPEXY']) and (typefile in ['meso']) and (all_varname[k-1] not in ncref.variables)):
330	all_var[k] = slopeamplitude(ncref)
331	elif ((all_varname[k-1] in ['DELTAT','deltat']) and (typefile in ['meso']) and (all_varname[k-1] not in ncref.variables)):
332	all_var[k]=deltat0t1(ncref)
333	else:
334	all_var[k] = getfield(ncref,all_varname[k-1])
335	all_varname[k] = all_varname[k-1] ; all_time[k] = all_time[k-1] ; all_namefile[k] = all_namefile[k-1] ; all_var2[k] = all_var2[k-1]
336	if winds: all_windu[k] = getfield(ncref,uchar) ; all_windv[k] = getfield(ncref,vchar)
337	else: errormess("fileref is missing!")
338	if ope == "-": all_var[k+1]= all_var[k-1] - all_var[k]
339	elif ope == "+": all_var[k+1]= all_var[k-1] + all_var[k]
340	elif ope == "-%":
341	masked = np.ma.masked_where(all_var[k] == 0,all_var[k])
342	masked.set_fill_value([np.NaN])
343	all_var[k+1]= 100.*(all_var[k-1] - masked)/masked
344	all_varname[k+1] = all_varname[k] ; all_time[k+1] = all_time[k] ; all_namefile[k+1] = all_namefile[k] ; all_var2[k+1] = all_var2[k] ; numplot = numplot+2
345	if winds: all_windu[k+1] = all_windu[k-1]-all_windu[k] ; all_windv[k+1] = all_windv[k-1] - all_windv[k]
346	elif ope in ["cat"]:
347	tabtime = all_time[0];tab = all_var[0];k = 1
348	if var2: tab2 = all_var2[0]
349	while k != len(namefiles) and len(all_time[k]) != 0:
350	if var2: tab2 = np.append(tab2,all_var2[k],axis=0)
351	tabtime = np.append(tabtime,all_time[k]) ; tab = np.append(tab,all_var[k],axis=0) ; k += 1
352	all_time[0] = np.array(tabtime) ; all_var[0] = np.array(tab) ; numplot = 1
353	if var2: all_var2[0] = np.array(tab2)
354	else: errormess(ope+" : non-implemented operation. Check pp.py --help")
355	else:
356	if len(namefiles) > 1: errormess("for this operation... please set only one file !")
357	if len(var) > 2: errormess("not sure this works for more than 2 vars... please check.")
358	if ope in ["div_var"]: all_var[k] = all_var[k-2] / all_var[k-1] ; insert = '_div_'
359	elif ope in ["mul_var"]: all_var[k] = all_var[k-2] * all_var[k-1] ; insert = '_mul_'
360	elif ope in ["add_var"]: all_var[k] = all_var[k-2] + all_var[k-1] ; insert = '_add_'
361	elif ope in ["sub_var"]: all_var[k] = all_var[k-2] - all_var[k-1] ; insert = '_sub_'
362	else: errormess(ope+" : non-implemented operation. Check pp.py --help")
363	numplot = numplot + 1 ; all_time[k] = all_time[k-1] ; all_namefile[k] = all_namefile[k-1]
364	all_varname[k] = all_varname[k-2] + insert + all_varname[k-1]
365	##################################
366	### Open a figure and set subplots
367	fig = figure()
368	subv,subh = definesubplot( numplot, fig )
369	if ope in ['-','-%']: subv,subh = 2,2
370
371	#################################
372	### Time loop for plotting device
373	nplot = 1 ; error = False
374	if lstyle is not None: linecycler = cycle(lstyle)
375	else: linecycler = cycle(["-","--","-.",":"])
376	print "********************************************"
377	while error is False:
378
379	print "********** NPLOT", nplot
380	if nplot > numplot: break
381
382	####################################################################
383	## get all indexes to be taken into account for this subplot and then reduce field
384	## 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
385	indexlon = getsindex(sslon,(nplot-1)%nlon,lon)
386	indexlat = getsindex(sslat,((nplot-1)//nlon)%nlat,lat)
387	indexvert = getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert)
388	if ope is not None:
389	if fileref is not None: index_f = ((nplot-1)//(nlonnlatnvert*ntime))%(len(namefiles)+2) ## OK only 1 var, see test in the beginning
390	elif "var" in ope: index_f = ((nplot-1)//(nlonnlatnvert*ntime))%(len(var)+1) ## OK only 1 file, see test in the beginning
391	elif "cat" in ope: index_f = 0
392	else: yeah = len(namefiles)len(var) ; index_f = ((nplot-1)//(nlonnlatnvertntime))%yeah
393	time = all_time[index_f]
394	if mrate is not None: indextime = None
395	else: indextime = getsindex(stime,((nplot-1)//(nlonnlatnvert))%ntime,time)
396	ltst = None
397	if typefile in ['meso'] and indextime is not None: ltst = localtime ( interv[0]+indextimeinterv[1], 0.5(wlon[0]+wlon[1]) )
398	print "********** INDEX LON:",indexlon," LAT:",indexlat," VERT:",indexvert," TIME:",indextime
399	##var = nc.variables["phisinit"][:,:]
400	##contourf(np.transpose(var),30,cmap = get_cmap(name="Greys_r") ) ; axis('off') ; plot(indexlat,indexlon,'mx',mew=4.0,ms=20.0)
401	##show()
402	##exit()
403	#truc = True
404	#truc = False
405	#if truc: indexvert = None
406	####################################################################
407	########## REDUCE FIELDS
408	####################################################################
409	error = False
410	varname = all_varname[index_f]
411	if varname: ### what is shaded.
412	what_I_plot, error = reducefield( all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \
413	yint=yintegral, alt=vert, anomaly=anomaly, redope=redope, mesharea=area )
414	if mult != 2718.: what_I_plot = what_I_plot*mult
415	else: what_I_plot = np.log10(what_I_plot) ; print "log plot"
416	if var2: ### what is contoured.
417	what_I_plot_contour, error = reducefield( all_var2[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert, \
418	yint=yintegral, alt=vert )
419	if winds: ### what is plot as vectors.
420	vecx, error = reducefield( all_windu[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert)
421	vecy, error = reducefield( all_windv[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert)
422	if varname in [uchar,vchar]: what_I_plot = np.sqrt( np.square(vecx) + np.square(vecy) ) ; varname = "wind"
423
424	#####################################################################
425	#if truc:
426	# nx = what_I_plot.shape[2] ; ny = what_I_plot.shape[1] ; nz = what_I_plot.shape[0]
427	# for k in range(nz): print k,' over ',nz ; what_I_plot[k,:,:] = what_I_plot[k,:,:] / smooth(what_I_plot[k,:,:],12)
428	# for iii in range(nx):
429	# for jjj in range(ny):
430	# deviation = what_I_plot[:,jjj,iii] ; mx = max(deviation) ; mn = min(deviation)
431	# 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
432	# else: what_I_plot[0,jjj,iii] = 0.
433	# if np.abs(what_I_plot[0,jjj,iii]) > 1.5:
434	# print iii,jjj,what_I_plot[0,jjj,iii],int(abs(1.-mx)100.),int(abs(1.-mn)100.)
435	# plot(rel)
436	# show()
437	# anomaly = True ### pour avoir les bons reglages plots
438	# what_I_plot = what_I_plot[0,:,:]
439	#####################################################################
440
441	####################################################################
442	### General plot settings
443	changesubplot = (numplot > 1) and (len(what_I_plot.shape) != 1) ## default for 1D plots: superimposed. to be reworked for better flexibility.
444	if changesubplot: subplot(subv,subh,nplot) #; subplots_adjust(wspace=0,hspace=0)
445	####################################################################
446	if error:
447	errormess("There is an error in reducing field !")
448	else:
449	ticks = ndiv + 1
450	fvar = varname
451	if anomaly: fvar = 'anomaly'
452	###
453	if mapmode == 0: ### could this be moved inside imov loop ?
454	itime=indextime
455	if len(what_I_plot.shape) == 3: itime=[0]
456	m = None ; x = None ; y = None
457	latyeah = lat ; lonyeah = lon
458	if typefile in ['meso']:
459	# now that the section is determined we can set the real lat
460	# ... or for now, a temptative one.
461	if slon is not None: latyeah = lat2d[:,0]
462	if slat is not None: lonyeah = lon2d[0,:]
463	what_I_plot, x, y = define_axis(lonyeah,latyeah,vert,time,indexlon,indexlat,indexvert,\
464	itime,what_I_plot, len(all_var[index_f].shape),vertmode)
465	###
466	if (fileref is not None) and (index_f == numplot-1): zevmin, zevmax = calculate_bounds(what_I_plot,vmin=minop,vmax=maxop)
467	else: zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
468	if (fileref is not None) and (index_f == numplot-1): colorb = "RdBu_r"
469	if colorb in ["def","nobar","onebar"]: palette = get_cmap(name=defcolorb(fvar.upper()))
470	else: palette = get_cmap(name=colorb)
471	#palette = cm.GMT_split
472	##### 1. ELIMINATE 0D or >3D CASES
473	if len(what_I_plot.shape) == 0:
474	print "VALUE VALUE VALUE VALUE ::: ",what_I_plot
475	save = 'donothing'
476	elif len(what_I_plot.shape) >= 4:
477	print "WARNING!!! ",len(what_I_plot.shape),"-D PLOT NOT SUPPORTED !!! dimensions: ",what_I_plot.shape
478	errormess("Are you sure you did not forget to prescribe a dimension ?")
479	##### 2. HANDLE simple 1D/2D field and movies of 1D/2D fields
480	else:
481	if mrate is not None: iend=len(time)-1
482	else: iend=0
483	imov = 0
484	if len(what_I_plot.shape) == 3:
485	if var2: which = "contour" ## have to start with contours rather than shading
486	else: which = "regular"
487	if mrate is None: errormess("3D field. Use --rate RATE for movie or specify --time TIME. Exit.")
488	elif len(what_I_plot.shape) == 2:
489	if var2: which = "contour" ## have to start with contours rather than shading
490	else: which = "regular"
491	if mrate is not None: which = "unidim"
492	elif len(what_I_plot.shape) == 1:
493	which = "unidim"
494	if what_I_plot.shape[-1] == 1: print "VALUE VALUE VALUE VALUE ::: ", what_I_plot[0] ; save = 'donothing'
495	##### IMOV LOOP #### IMOV LOOP
496	while imov <= iend:
497	print "-> frame ",imov+1, which
498	if which == "regular":
499	if mrate is None: what_I_plot_frame = what_I_plot
500	else: what_I_plot_frame = what_I_plot[imov,:,:]
501	if winds:
502	if mrate is None: vecx_frame = vecx ; vecy_frame = vecy
503	else: vecx_frame = vecx[imov,:,:] ; vecy_frame = vecy[imov,:,:]
504	elif which == "contour":
505	if mrate is None or what_I_plot_contour.ndim < 3: what_I_plot_frame = what_I_plot_contour
506	else: what_I_plot_frame = what_I_plot_contour[imov,:,:]
507	elif which == "unidim":
508	if mrate is None: what_I_plot_frame = what_I_plot
509	else: what_I_plot_frame = what_I_plot[:,imov] ## because swapaxes
510	#if mrate is not None:
511	if mapmode == 1:
512	m = define_proj(proj,wlon,wlat,back=back,blat=blat,blon=blon) ## this is dirty, defined above but out of imov loop
513	x, y = m(lon2d, lat2d) ## this is dirty, defined above but out of imov loop
514	if typefile in ['meso'] and mapmode == 1: what_I_plot_frame = dumpbdy(what_I_plot_frame,6,condition=True)
515	# if typefile in ['mesoideal']: what_I_plot_frame = dumpbdy(what_I_plot_frame,0,stag='W',condition=dumped_vert_stag)
516
517	if which == "unidim":
518	if lbls is not None: lbl=lbls[nplot-1]
519	else:
520	lbl = ""
521	if indexlat is not None: lbl = lbl + " ix" + str(indexlat[0])
522	if indexlon is not None: lbl = lbl + " iy" + str(indexlon[0])
523	if indexvert is not None: lbl = lbl + " iz" + str(indexvert[0])
524	if indextime is not None: lbl = lbl + " it" + str(indextime[0])
525	if lbl == "": lbl = namefiles[index_f]
526
527	if mrate is not None: x = y ## because swapaxes...
528	#what_I_plot_frame = np.diff(what_I_plot_frame, n=1) ; x = x[1:]
529
530	zeline = next(linecycler) ## "-" for simple lines
531	if tile: zemarker = 'x'
532	else: zemarker = None
533	this_is_a_regular_plot = (indexvert is not None) or (indextime is None) or (indexlat is None) or (indexlon is None)
534	if this_is_a_regular_plot: plot(x,what_I_plot_frame,zeline,label=lbl,marker=zemarker) ## vertical profile
535	else: plot(what_I_plot_frame,x,zeline,label=lbl,marker=zemarker) ## regular plot
536	if nplot > 1: legend(loc='best')
537	if indextime is None and axtime is not None and xlab is None: xlabel(axtime.upper()) ## define the right label
538	if save == 'txt': writeascii(np.transpose([x,np.transpose(what_I_plot)]),'profile'+str(nplot*1000+imov)+'.txt')
539
540	elif which == "regular":
541	if hole: what_I_plot_frame = hole_bounds(what_I_plot_frame,zevmin,zevmax)
542	else: what_I_plot_frame = bounds(what_I_plot_frame,zevmin,zevmax)
543	if flagnolow: what_I_plot_frame = nolow(what_I_plot_frame)
544	if not tile:
545	#zelevels = np.linspace(zevmin(1. + 1.e-7),zevmax(1. - 1.e-7)) #,num=20)
546	zelevels = np.linspace(zevmin,zevmax,num=ticks)
547	#what_I_plot_frame = smooth(what_I_plot_frame,100)
548	if mapmode == 1: m.contourf( x, y, what_I_plot_frame, zelevels, cmap = palette, alpha=trans)
549	elif mapmode == 0: contourf( x, y, what_I_plot_frame, zelevels, cmap = palette, alpha=trans)
550	if cross is not None and mapmode == 1:
551	idx,idy=m(cross[0][0],cross[0][1])
552	mpl.pyplot.plot([idx],[idy],'+k',mew=0.5,ms=18)
553	else:
554	if mapmode == 1: m.pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax, alpha=trans)
555	elif mapmode == 0: pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax, alpha=trans)
556
557	if colorb not in ['nobar','onebar']:
558	if (fileref is not None) and (index_f == numplot-1): daformat = "%.3f"
559	elif mult != 1: daformat = "%.1f"
560	else: daformat = fmtvar(fvar.upper())
561	#if proj in ['moll','cyl']: zeorientation="horizontal" ; zepad = 0.05
562	if proj in ['moll']: zeorientation="horizontal" ; zepad = 0.05
563	else: zeorientation="vertical" ; zepad = 0.03
564	zecb = colorbar( fraction=0.05,pad=zepad,format=daformat,orientation=zeorientation,\
565	ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,21])),extend='neither',spacing='proportional' )
566	if zeorientation == "horizontal" and zetitle[0] != "fill": zecb.ax.set_xlabel(zetitle[index_f]) ; zetitle[index_f]=""
567	if winds:
568	if typefile in ['meso']:
569	[vecx_frame,vecy_frame] = [dumpbdy(vecx_frame,6,stag=uchar,condition=True), dumpbdy(vecy_frame,6,stag=vchar,condition=True)]
570	key = True
571	if fvar in ['UV','uv','uvmet']: key = False
572	elif typefile in ['gcm']:
573	key = False
574	if metwind and mapmode == 1: [vecx_frame,vecy_frame] = m.rotate_vector(vecx_frame, vecy_frame, lon2d, lat2d)
575	if var: colorvec = definecolorvec(colorb)
576	else: colorvec = definecolorvec(back)
577	vectorfield(vecx_frame, vecy_frame, x, y, stride=stride, csmooth=2,\
578	#scale=15., factor=300., color=colorvec, key=key)
579	scale=20., factor=250., color=colorvec, key=key)
580	#200. ## or csmooth=stride
581	if ope == '-' and nplot == numplot: # this should work as long as ope is '-' guarantees 3 plots for 4 panels without contour
582	subplot(subv,subh,nplot+1)
583	rcParams["legend.fontsize"] = 'xx-large'
584	if indexlat is None:
585	latmin = -50.; latmax = 50. # latitude range for histogram of difference
586	zeindexlat = (lat<latmax)*(lat>latmin)
587	# this follows the define_axis logic in myplot.py:
588	if indextime is None or indexlon is None: what_I_plot_frame = what_I_plot_frame[zeindexlat,:]
589	else: what_I_plot_frame = what_I_plot_frame[:,zeindexlat]
590	toplot = np.ravel(what_I_plot_frame[np.isnan(what_I_plot_frame)==False])
591	zebins = np.linspace(zevmin,zevmax,num=30)
592	hist(toplot,bins=zebins,histtype='step',linewidth=2,color='k',normed=True)
593	zestd = np.std(toplot)
594	zemean = np.mean(toplot)
595	zebins = np.linspace(zevmin,zevmax,num=300)
596	zegauss = (1./(zestd * np.sqrt(2 * np.pi)) * np.exp( - (zebins - zemean)*2 / (2 zestd**2) ) )
597	plot(zebins, zegauss, linewidth=1, color='r',label="mean: "+str(zemean)[0:5]+"\nstd: "+str(zestd)[0:5])
598	legend(loc=0,frameon=False)
599	title("Histogram fig(1) "+ope+" fig(2)")
600	subplot(subv,subh,nplot) # go back to last plot for title of contour difference
601
602	elif which == "contour":
603	zevminc, zevmaxc = calculate_bounds(what_I_plot_frame, vmin=min(what_I_plot_frame), vmax=max(what_I_plot_frame))
604	zelevels = np.linspace(zevminc,zevmaxc,ticks/2) #20)
605	if var2 == 'HGT': zelevels = np.arange(-10000.,30000.,1000.)
606	elif var2 == 'tpot': zelevels = np.arange(270,370,5)
607	elif var2 == 'tk': zelevels = np.arange(150,250,5)
608	if mapmode == 0:
609	what_I_plot_frame, x, y = define_axis( lonyeah,latyeah,vert,time,indexlon,indexlat,indexvert,\
610	itime,what_I_plot_frame, len(all_var2[index_f].shape),vertmode )
611	## is this needed? only if len(all_var2[index_f].shape) != len(all_var[index_f].shape)
612	cs = contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 0.33)#, alpha=0.5, linestyles='solid')
613	clabel(cs, inline=1, fontsize = 4.*rcParams['font.size']/5., fmt=fmtvar(var2.upper()))
614	elif mapmode == 1: cs = m.contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 0.33)#, alpha=0.5, linestyles='solid')
615	if which in ["regular","unidim"]:
616
617	if nplot > 1 and which == "unidim":
618	pass ## because we superimpose nplot instances
619	else:
620	# Axis directives for movie frames [including the first one).
621	zxmin, zxmax = xaxis ; zymin, zymax = yaxis
622	if zxmin is not None: mpl.pyplot.xlim(xmin=zxmin)
623	if zxmax is not None: mpl.pyplot.xlim(xmax=zxmax)
624	if zymin is not None: mpl.pyplot.ylim(ymin=zymin)
625	if zymax is not None: mpl.pyplot.ylim(ymax=zymax)
626	if ylog: mpl.pyplot.semilogy()
627	if invert_y: ax = mpl.pyplot.gca() ; ax.set_ylim(ax.get_ylim()[::-1])
628	if xlab is not None: xlabel(xlab)
629	if ylab is not None: ylabel(ylab)
630
631	if mrate is not None:
632	### THIS IS A MENCODER MOVIE
633	if mrate > 0:
634	figframe=mpl.pyplot.gcf()
635	if mquality: figframe.set_dpi(600.)
636	else: figframe.set_dpi(200.)
637	mframe=fig2img(figframe)
638	if imov == 0:
639	moviename='movie' ;W,H = figframe.canvas.get_width_height()
640	video = VideoSink((H,W), moviename, rate=mrate, byteorder="rgba")
641	video.run(mframe) ; close()
642	if imov == iend: video.close()
643	### THIS IS A WEBPAGE MOVIE
644	else:
645	nameframe = "image"+str(1000+imov)
646	makeplotres(nameframe,res=100.,disp=False) ; close()
647	if imov == 0: myfile = open("zepics", 'w')
648	myfile.write("modImages["+str(imov)+"] = '"+nameframe+"_100.png';"+ '\n')
649	if imov == iend:
650	myfile.write("first_image = 0;"+ '\n')
651	myfile.write("last_image = "+str(iend)+";"+ '\n')
652	myfile.close()
653	if var2 and which == "regular": which = "contour"
654	imov = imov+1
655	elif which == "contour":
656	which = "regular"
657
658	### Next subplot
659	zevarname = varname
660	if redope is not None: zevarname = zevarname + "_" + redope
661	basename = getname(var=zevarname,var2=var2,winds=winds,anomaly=anomaly)
662	if len(what_I_plot.shape) > 3:
663	basename = basename + getstralt(nc,level)
664	if mrate is not None: basename = "movie_" + basename
665	if typefile in ['meso']:
666	if sslon is not None: basename = basename + "_lon_" + str(int(round(lonp)))
667	if sslat is not None: basename = basename + "_lat_" + str(int(round(latp)))
668	plottitle = basename
669	### dans le nouveau systeme time=ls,sol,lt cette ligne pourrait ne servir a rien (ou deplacer au dessus)
670	if addchar and indextime is not None: [addchar,gogol,gogol2] = getlschar ( all_namefile[index_f] ) ; plottitle = plottitle + addchar
671	### en fait redope is None doit etre remplace par : n'est ni maxt ni mint
672	if redope is None and ltst is not None and ( (mapmode == 0) or (proj in ["lcc","laea","merc","nsper"]) ): plottitle = plottitle + "_LT" + str(ltst)
673	else:
674	if fileref is not None:
675	if index_f is numplot-1: plottitle = basename+' '+"fig(1) "+ope+" fig(2)"
676	elif index_f is numplot-2: plottitle = basename+' '+fileref
677	else: plottitle = basename+' '+namefiles[0]#index_f]
678	else: plottitle = basename+' '+namefiles[0]#index_f]
679	if mult != 1: plottitle = '{:.0e}'.format(mult) + "*" + plottitle
680	if zetitle[0] != "fill":
681	if index_f < len(zetitle): plottitle = zetitle[index_f]
682	else: plottitle = zetitle[0]
683	if titleref is "fill": titleref=zetitle[index_f]
684	if fileref is not None:
685	if index_f is numplot-2: plottitle = titleref
686	if index_f is numplot-1: plottitle = "fig(1) "+ope+" fig(2)"
687	# if indexlon is not None: plottitle = plottitle + " lon: " + str(min(lon[indexlon])) +" "+ str(max(lon[indexlon]))
688	# if indexlat is not None: plottitle = plottitle + " lat: " + str(min(lat[indexlat])) +" "+ str(max(lat[indexlat]))
689	# if indexvert is not None: plottitle = plottitle + " vert: " + str(min(vert[indexvert])) +" "+ str(max(vert[indexvert]))
690	# if indextime is not None: plottitle = plottitle + " time: " + str(min(time[indextime])) +" "+ str(max(time[indextime]))
691	if colorb != "onebar": title( plottitle )
692	if nplot >= numplot: error = True
693	nplot += 1
694
695	if colorb == "onebar":
696	cax = axes([0.1, 0.2, 0.8, 0.03]) # a ameliorer
697	zecb = colorbar(cax=cax, orientation="horizontal", format=fmtvar(fvar.upper()),\
698	ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,21])),extend='neither',spacing='proportional')
699	if zetitle[0] != "fill": zecb.ax.set_xlabel(zetitle[index_f]) ; zetitle[index_f]=""
700
701
702	##########################################################################
703	### Save the figure in a file in the data folder or an user-defined folder
704	if outputname is None:
705	if typefile in ['meso']: prefix = getprefix(nc)
706	elif typefile in ['gcm']: prefix = 'LMD_GCM_'
707	else: prefix = ''
708	###
709	zeplot = prefix + basename
710	if addchar: zeplot = zeplot + addchar
711	if numplot <= 0: zeplot = zeplot + "_LT"+str(abs(numplot))
712	###
713	if not target: zeplot = namefile[0:find(namefile,'wrfout')] + zeplot
714	else: zeplot = target + "/" + zeplot
715	###
716	else:
717	zeplot=outputname
718
719	if mrate is None:
720	pad_inches_value = 0.35
721	print "********** SAVE ", save
722	if save == 'png':
723	if display: makeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value) #,erase=True) ## a miniature
724	makeplotres(zeplot,res=resolution,pad_inches_value=pad_inches_value,disp=False)
725	elif save in ['eps','svg','pdf']: makeplotres(zeplot,pad_inches_value=pad_inches_value,disp=False,ext=save)
726	elif save == 'gui': show()
727	elif save == 'donothing': pass
728	elif save == 'txt': print "Saved results in txt file."
729	else:
730	print "INFO: save mode not supported. using gui instead."
731	show()
732
733	###################################
734	#### Getting more out of this video -- PROBLEMS WITH CREATED VIDEOS
735	#
736	#if mrate is not None:
737	# print "Re-encoding movie.. first pass"
738	# video.first_pass(filename=moviename,quality=mquality,rate=mrate)
739	# print "Re-encoding movie.. second pass"
740	# video.second_pass(filename=moviename,quality=mquality,rate=mrate)
741
742	###############
743	### Now the end
744	return zeplot

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

Download in other formats: