Context Navigation

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

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

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

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

Download in other formats: