Context Navigation

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

Last change on this file since 723 was 721, checked in by aslmd, 13 years ago
UTIL PYTHON. Added an option to set magnifying factor for winds. added support for GCM v5 files. added a handful of pretty cool background maps (titan, venus, triton, cobe, pluto). added the possibility to use div_var_only etc... instead of div_var to get only the operation plot. corrected bugs : if transparency is zero color for vectors must be set according to background not colorbar; added a line which avoid a bug with histograms and mesoscale files.
Property svn:executable set to ``*
File size: 49.3 KB

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

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

Download in other formats: