Context Navigation

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

Last change on this file since 886 was 886, checked in by aslmd, 12 years ago
UTIL PYTHON. png thumbnails are now hidden and name of pic does not have _res anymore
Property svn:executable set to ``*
File size: 67.3 KB

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

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

Download in other formats: