Context Navigation

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

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

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

Download in other formats: