Context Navigation

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

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

Rev	Line
[345]	1	## Author: AS
[252]	2	def errormess(text,printvar=None):
[233]	3	print text
[399]	4	if printvar is not None: print printvar
[233]	5	exit()
	6	return
	7
[345]	8	## Author: AS
[349]	9	def adjust_length (tab, zelen):
	10	from numpy import ones
	11	if tab is None:
	12	outtab = ones(zelen) * -999999
	13	else:
	14	if zelen != len(tab):
	15	print "not enough or too much values... setting same values all variables"
	16	outtab = ones(zelen) * tab[0]
	17	else:
	18	outtab = tab
	19	return outtab
	20
	21	## Author: AS
[468]	22	def getname(var=False,var2=False,winds=False,anomaly=False):
[252]	23	if var and winds: basename = var + '_UV'
	24	elif var: basename = var
	25	elif winds: basename = 'UV'
	26	else: errormess("please set at least winds or var",printvar=nc.variables)
	27	if anomaly: basename = 'd' + basename
[468]	28	if var2: basename = basename + '_' + var2
[252]	29	return basename
	30
[763]	31	## Author: AS + AC
[782]	32	def localtime(time,lon,namefile): # lon is the mean longitude of the plot, not of the domain. central lon of domain is taken from cen_lon
[763]	33	import numpy as np
[782]	34	from netCDF4 import Dataset
	35	## THIS IS FOR MESOSCALE
	36	nc = Dataset(namefile)
	37	## get start date and intervals
	38	dt_hour=1. ; start=0.
	39	if hasattr(nc,'TITLE'):
	40	title=getattr(nc, 'TITLE')
	41	if hasattr(nc,'DT') and hasattr(nc,'START_DATE') and 'MRAMS' in title:
	42	## we must adapt what is done in getlschar to MRAMS (outputs from ic.py)
	43	dt_hour=getattr(nc, 'DT')/60.
	44	start_date=getattr(nc, 'START_DATE')
	45	start_hour=np.float(start_date[11:13])
	46	start_minute=np.float(start_date[14:16])/60.
	47	start=start_hour+start_minute # start is the local time of simu at longitude 0
	48	#LMD MMM is 1 output/hour (and not 1 output/timestep)
	49	#MRAMS is 1 output/timestep, unless stride is added in ic.py
	50	elif 'WRF' in title:
	51	[dummy,start,dt_hour] = getlschar ( namefile ) # get start hour and interval hour
	52	## get longitude
[763]	53	if lon is not None:
	54	if lon[0,1]!=lon[0,0]: mean_lon_plot = 0.5*(lon[0,1]-lon[0,0])
	55	else: mean_lon_plot=lon[0,0]
	56	elif hasattr(nc, 'CEN_LON'): mean_lon_plot=getattr(nc, 'CEN_LON')
	57	else: mean_lon_plot=0.
[782]	58	## calculate local time
	59	ltst = start + time*dt_hour + mean_lon_plot / 15.
[252]	60	ltst = int (ltst * 10) / 10.
	61	ltst = ltst % 24
	62	return ltst
	63
[569]	64	## Author: AC
	65	def check_localtime(time):
	66	a=-1
[763]	67	print time
[569]	68	for i in range(len(time)-1):
[583]	69	if (time[i] > time[i+1]): a=i
	70	if a >= 0 and a < (len(time)-1)/2.:
[569]	71	print "Sorry, time axis is not regular."
	72	print "Contourf needs regular axis... recasting"
	73	for i in range(a+1):
	74	time[i]=time[i]-24.
[583]	75	if a >= 0 and a >= (len(time)-1)/2.:
	76	print "Sorry, time axis is not regular."
	77	print "Contourf needs regular axis... recasting"
	78	for i in range((len(time)-1) - a):
	79	time[a+1+i]=time[a+1+i]+24.
[569]	80	return time
	81
[525]	82	## Author: AS, AC, JL
[233]	83	def whatkindfile (nc):
[647]	84	typefile = 'gcm' # default
[429]	85	if 'controle' in nc.variables: typefile = 'gcm'
	86	elif 'phisinit' in nc.variables: typefile = 'gcm'
[721]	87	elif 'phis' in nc.variables: typefile = 'gcm'
[525]	88	elif 'time_counter' in nc.variables: typefile = 'earthgcm'
[548]	89	elif hasattr(nc,'START_DATE'): typefile = 'meso'
[429]	90	elif 'HGT_M' in nc.variables: typefile = 'geo'
[558]	91	elif hasattr(nc,'institution'):
	92	if "European Centre" in getattr(nc,'institution'): typefile = 'ecmwf'
[233]	93	return typefile
	94
[345]	95	## Author: AS
[233]	96	def getfield (nc,var):
	97	## this allows to get much faster (than simply referring to nc.variables[var])
[395]	98	import numpy as np
[233]	99	dimension = len(nc.variables[var].dimensions)
[392]	100	#print " Opening variable with", dimension, "dimensions ..."
[233]	101	if dimension == 2: field = nc.variables[var][:,:]
	102	elif dimension == 3: field = nc.variables[var][:,:,:]
	103	elif dimension == 4: field = nc.variables[var][:,:,:,:]
[494]	104	elif dimension == 1: field = nc.variables[var][:]
[395]	105	# if there are NaNs in the ncdf, they should be loaded as a masked array which will be
	106	# recasted as a regular array later in reducefield
	107	if (np.isnan(np.sum(field)) and (type(field).__name__ not in 'MaskedArray')):
	108	print "Warning: netcdf as nan values but is not loaded as a Masked Array."
	109	print "recasting array type"
	110	out=np.ma.masked_invalid(field)
	111	out.set_fill_value([np.NaN])
	112	else:
[464]	113	# missing values from zrecast or hrecast are -1e-33
[469]	114	masked=np.ma.masked_where(field < -1e30,field)
[578]	115	masked2=np.ma.masked_where(field > 1e35,field)
	116	masked.set_fill_value([np.NaN]) ; masked2.set_fill_value([np.NaN])
	117	mask = np.ma.getmask(masked) ; mask2 = np.ma.getmask(masked2)
	118	if (True in np.array(mask)):
	119	out=masked
	120	print "Masked array... Missing value is NaN"
	121	elif (True in np.array(mask2)):
	122	out=masked2
	123	print "Masked array... Missing value is NaN"
	124	# else:
	125	# # missing values from api are 1e36
	126	# masked=np.ma.masked_where(field > 1e35,field)
	127	# masked.set_fill_value([np.NaN])
	128	# mask = np.ma.getmask(masked)
	129	# if (True in np.array(mask)):out=masked
	130	# else:out=field
[763]	131	else:
	132	# # missing values from MRAMS files are 0.100E+32
	133	masked=np.ma.masked_where(field > 1e30,field)
	134	masked.set_fill_value([np.NaN])
	135	mask = np.ma.getmask(masked)
	136	if (True in np.array(mask)):out=masked
	137	else:out=field
	138	# else:out=field
[395]	139	return out
[233]	140
[571]	141	## Author: AC
[763]	142	# Compute the norm of the winds or return an hodograph
[612]	143	# The corresponding variable to call is UV or uvmet (to use api)
[763]	144	def windamplitude (nc,mode):
[581]	145	import numpy as np
[571]	146	varinfile = nc.variables.keys()
	147	if "U" in varinfile: zu=getfield(nc,'U')
	148	elif "Um" in varinfile: zu=getfield(nc,'Um')
[763]	149	else: errormess("you need slopex or U or Um in your file.")
[571]	150	if "V" in varinfile: zv=getfield(nc,'V')
	151	elif "Vm" in varinfile: zv=getfield(nc,'Vm')
[763]	152	else: errormess("you need V or Vm in your file.")
[571]	153	znt,znz,zny,znx = np.array(zu).shape
[763]	154	if hasattr(nc,'WEST-EAST_PATCH_END_UNSTAG'):znx=getattr(nc, 'WEST-EAST_PATCH_END_UNSTAG')
[571]	155	zuint = np.zeros([znt,znz,zny,znx])
	156	zvint = np.zeros([znt,znz,zny,znx])
	157	if "U" in varinfile:
[763]	158	if hasattr(nc,'SOUTH-NORTH_PATCH_END_STAG'): zny_stag=getattr(nc, 'SOUTH-NORTH_PATCH_END_STAG')
	159	if hasattr(nc,'WEST-EAST_PATCH_END_STAG'): znx_stag=getattr(nc, 'WEST-EAST_PATCH_END_STAG')
	160	if zny_stag == zny: zvint=zv
	161	else:
	162	for yy in np.arange(zny): zvint[:,:,yy,:] = (zv[:,:,yy,:] + zv[:,:,yy+1,:])/2.
	163	if znx_stag == znx: zuint=zu
	164	else:
	165	for xx in np.arange(znx): zuint[:,:,:,xx] = (zu[:,:,:,xx] + zu[:,:,:,xx+1])/2.
[571]	166	else:
	167	zuint=zu
	168	zvint=zv
[763]	169	if mode=='amplitude': return np.sqrt(zuint2 + zvint2)
	170	if mode=='hodograph': return zuint,zvint
	171	if mode=='hodograph_2': return None, 360.np.arctan(zvint/zuint)/(2.np.pi)
[571]	172
[612]	173	## Author: AC
[701]	174	# Compute the enrichment factor of non condensible gases
	175	# The corresponding variable to call is enfact
[753]	176	# enrichment factor is computed as in Yuan Lian et al. 2012
	177	# i.e. you need to have VL2 site at LS 135 in your data
	178	# this only requires co2col so that you can concat.nc at low cost
	179	def enrichment_factor(nc,lon,lat,time):
[701]	180	import numpy as np
[753]	181	from myplot import reducefield
[701]	182	varinfile = nc.variables.keys()
[753]	183	if "co2col" in varinfile: co2col=getfield(nc,'co2col')
	184	else: print "error, you need co2col var in your file"
[701]	185	if "ps" in varinfile: ps=getfield(nc,'ps')
	186	else: print "error, you need ps var in your file"
[753]	187	dimension = len(nc.variables['co2col'].dimensions)
	188	if dimension == 2:
	189	zny,znx = np.array(co2col).shape
[701]	190	znt=1
[753]	191	elif dimension == 3: znt,zny,znx = np.array(co2col).shape
[701]	192	mmrarcol = np.zeros([znt,zny,znx])
	193	enfact = np.zeros([znt,zny,znx])
	194	grav=3.72
[753]	195	mmrarcol[:,:,:] = 1. - grav*co2col[:,:,:]/ps[:,:,:]
	196	# Computation with reference argon mmr at VL2 Ls 135 (as in Yuan Lian et al 2012)
	197	lonvl2=np.zeros([1,2])
	198	latvl2=np.zeros([1,2])
	199	timevl2=np.zeros([1,2])
	200	lonvl2[0,0]=-180
	201	lonvl2[0,1]=180
	202	latvl2[:,:]=48.16
	203	timevl2[:,:]=135.
	204	indexlon = getsindex(lonvl2,0,lon)
	205	indexlat = getsindex(latvl2,0,lat)
	206	indextime = getsindex(timevl2,0,time)
	207	mmrvl2, error = reducefield( mmrarcol, d4=indextime, d1=indexlon, d2=indexlat)
	208	print "VL2 Ls 135 mmr arcol:", mmrvl2
	209	enfact[:,:,:] = mmrarcol[:,:,:]/mmrvl2
[701]	210	return enfact
	211
	212	## Author: AC
[612]	213	# Compute the norm of the slope angles
	214	# The corresponding variable to call is SLOPEXY
	215	def slopeamplitude (nc):
	216	import numpy as np
	217	varinfile = nc.variables.keys()
	218	if "slopex" in varinfile: zu=getfield(nc,'slopex')
	219	elif "SLOPEX" in varinfile: zu=getfield(nc,'SLOPEX')
[754]	220	else: errormess("you need slopex or SLOPEX in your file.")
[612]	221	if "slopey" in varinfile: zv=getfield(nc,'slopey')
	222	elif "SLOPEY" in varinfile: zv=getfield(nc,'SLOPEY')
[754]	223	else: errormess("you need slopey or SLOPEY in your file.")
[612]	224	znt,zny,znx = np.array(zu).shape
	225	zuint = np.zeros([znt,zny,znx])
	226	zvint = np.zeros([znt,zny,znx])
	227	zuint=zu
	228	zvint=zv
	229	return np.sqrt(zuint2 + zvint2)
	230
	231	## Author: AC
	232	# Compute the temperature difference between surface and first level.
	233	# API is automatically called to get TSURF and TK.
	234	# The corresponding variable to call is DELTAT
	235	def deltat0t1 (nc):
	236	import numpy as np
	237	varinfile = nc.variables.keys()
	238	if "tsurf" in varinfile: zu=getfield(nc,'tsurf')
	239	elif "TSURF" in varinfile: zu=getfield(nc,'TSURF')
[754]	240	else: errormess("You need tsurf or TSURF in your file")
[612]	241	if "tk" in varinfile: zv=getfield(nc,'tk')
	242	elif "TK" in varinfile: zv=getfield(nc,'TK')
[754]	243	else: errormess("You need tk or TK in your file. (might need to use API. try to add -i 4 -l XXX)")
[612]	244	znt,zny,znx = np.array(zu).shape
	245	zuint = np.zeros([znt,zny,znx])
	246	zuint=zu - zv[:,0,:,:]
	247	return zuint
	248
[382]	249	## Author: AS + TN + AC
[717]	250	def reducefield (input,d4=None,d3=None,d2=None,d1=None,yint=False,alt=None,anomaly=False,redope=None,mesharea=None,unidim=999):
[252]	251	### we do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x"
[233]	252	### it would be actually better to name d4 d3 d2 d1 as t z y x
[405]	253	### ... note, anomaly is only computed over d1 and d2 for the moment
[233]	254	import numpy as np
[647]	255	from mymath import max,mean,min,sum,getmask
[422]	256	csmooth = 12 ## a fair amount of grid points (too high results in high computation time)
[483]	257	if redope is not None:
	258	if redope == "mint": input = min(input,axis=0) ; d1 = None
	259	elif redope == "maxt": input = max(input,axis=0) ; d1 = None
[763]	260	elif redope == "edge_y1": input = input[:,:,0,:] ; d2 = None
	261	elif redope == "edge_y2": input = input[:,:,-1,:] ; d2 = None
	262	elif redope == "edge_x1": input = input[:,:,:,0] ; d1 = None
	263	elif redope == "edge_x2": input = input[:,:,:,-1] ; d1 = None
[483]	264	else: errormess("not supported. but try lines in reducefield beforehand.")
	265	#elif redope == "minz": input = min(input,axis=1) ; d2 = None
	266	#elif redope == "maxz": input = max(input,axis=1) ; d2 = None
	267	#elif redope == "miny": input = min(input,axis=2) ; d3 = None
	268	#elif redope == "maxy": input = max(input,axis=2) ; d3 = None
	269	#elif redope == "minx": input = min(input,axis=3) ; d4 = None
	270	#elif redope == "maxx": input = max(input,axis=3) ; d4 = None
[233]	271	dimension = np.array(input).ndim
[525]	272	shape = np.array(np.array(input).shape)
[349]	273	#print 'd1,d2,d3,d4: ',d1,d2,d3,d4
[405]	274	if anomaly: print 'ANOMALY ANOMALY'
[233]	275	output = input
	276	error = False
[350]	277	#### this is needed to cope the case where d4,d3,d2,d1 are single integers and not arrays
[345]	278	if d4 is not None and not isinstance(d4, np.ndarray): d4=[d4]
	279	if d3 is not None and not isinstance(d3, np.ndarray): d3=[d3]
	280	if d2 is not None and not isinstance(d2, np.ndarray): d2=[d2]
	281	if d1 is not None and not isinstance(d1, np.ndarray): d1=[d1]
	282	### now the main part
[233]	283	if dimension == 2:
[717]	284	#### this is needed for 1d-type files (where dim=2 but axes are time-vert and not lat-lon)
[753]	285	if unidim==1: d2=d4 ; d1=d3 ; d4 = None ; d3 = None
[525]	286	if mesharea is None: mesharea=np.ones(shape)
	287	if max(d2) >= shape[0]: error = True
	288	elif max(d1) >= shape[1]: error = True
	289	elif d1 is not None and d2 is not None:
[687]	290	try:
	291	totalarea = np.ma.masked_where(getmask(output),mesharea)
	292	totalarea = mean(totalarea[d2,:],axis=0);totalarea = mean(totalarea[d1])
	293	except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
[647]	294	output = output*mesharea; output = mean(output[d2,:],axis=0); output = mean(output[d1])/totalarea
[350]	295	elif d1 is not None: output = mean(input[:,d1],axis=1)
[647]	296	elif d2 is not None:
[687]	297	try:
	298	totalarea = np.ma.masked_where(getmask(output),mesharea)
	299	totalarea = mean(totalarea[d2,:],axis=0)
	300	except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
[647]	301	output = output*mesharea; output = mean(output[d2,:],axis=0)/totalarea
[233]	302	elif dimension == 3:
[525]	303	if mesharea is None: mesharea=np.ones(shape[[1,2]])
[345]	304	if max(d4) >= shape[0]: error = True
	305	elif max(d2) >= shape[1]: error = True
	306	elif max(d1) >= shape[2]: error = True
[647]	307	elif d4 is not None and d2 is not None and d1 is not None:
[525]	308	output = mean(input[d4,:,:],axis=0)
[687]	309	try:
	310	totalarea = np.ma.masked_where(getmask(output),mesharea)
	311	totalarea = mean(totalarea[d2,:],axis=0);totalarea = mean(totalarea[d1])
	312	except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
[647]	313	output = output*mesharea; output = mean(output[d2,:],axis=0); output = mean(output[d1])/totalarea
[525]	314	elif d4 is not None and d2 is not None:
	315	output = mean(input[d4,:,:],axis=0)
[687]	316	try:
	317	totalarea = np.ma.masked_where(getmask(output),mesharea)
	318	totalarea = mean(totalarea[d2,:],axis=0)
	319	except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
[647]	320	output = output*mesharea; output = mean(output[d2,:],axis=0)/totalarea
[349]	321	elif d4 is not None and d1 is not None: output = mean(input[d4,:,:],axis=0); output=mean(output[:,d1],axis=1)
[525]	322	elif d2 is not None and d1 is not None:
[687]	323	try:
	324	totalarea = np.tile(mesharea,(output.shape[0],1,1))
	325	totalarea = np.ma.masked_where(getmask(output),totalarea)
	326	totalarea = mean(totalarea[:,d2,:],axis=1);totalarea = mean(totalarea[:,d1],axis=1)
	327	except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
[647]	328	output = output*mesharea; output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)/totalarea
[525]	329	elif d1 is not None: output = mean(input[:,:,d1],axis=2)
[647]	330	elif d2 is not None:
[687]	331	try:
	332	totalarea = np.tile(mesharea,(output.shape[0],1,1))
	333	totalarea = np.ma.masked_where(getmask(output),totalarea)
	334	totalarea = mean(totalarea[:,d2,:],axis=1)
	335	except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
[647]	336	output = output*mesharea; output = mean(output[:,d2,:],axis=1)/totalarea
[525]	337	elif d4 is not None: output = mean(input[d4,:,:],axis=0)
[233]	338	elif dimension == 4:
[647]	339	if mesharea is None: mesharea=np.ones(shape[[2,3]]) # mesharea=np.random.random_sample(shape[[2,3]])*5. + 2. # pour tester
[345]	340	if max(d4) >= shape[0]: error = True
	341	elif max(d3) >= shape[1]: error = True
	342	elif max(d2) >= shape[2]: error = True
	343	elif max(d1) >= shape[3]: error = True
[382]	344	elif d4 is not None and d3 is not None and d2 is not None and d1 is not None:
[392]	345	output = mean(input[d4,:,:,:],axis=0)
	346	output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
[427]	347	if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
[687]	348	try:
	349	totalarea = np.ma.masked_where(np.isnan(output),mesharea)
	350	totalarea = mean(totalarea[d2,:],axis=0); totalarea = mean(totalarea[d1])
	351	except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
[647]	352	output = output*mesharea; output = mean(output[d2,:],axis=0); output = mean(output[d1])/totalarea
[350]	353	elif d4 is not None and d3 is not None and d2 is not None:
[392]	354	output = mean(input[d4,:,:,:],axis=0)
	355	output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
[525]	356	if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
[687]	357	try:
	358	totalarea = np.ma.masked_where(np.isnan(output),mesharea)
	359	totalarea = mean(totalarea[d2,:],axis=0)
	360	except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
[647]	361	output = output*mesharea; output = mean(output[d2,:],axis=0)/totalarea
[350]	362	elif d4 is not None and d3 is not None and d1 is not None:
[392]	363	output = mean(input[d4,:,:,:],axis=0)
	364	output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
[405]	365	if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
[392]	366	output = mean(output[:,d1],axis=1)
[350]	367	elif d4 is not None and d2 is not None and d1 is not None:
[392]	368	output = mean(input[d4,:,:,:],axis=0)
[405]	369	if anomaly:
	370	for k in range(output.shape[0]): output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.)
[687]	371	try:
	372	totalarea = np.tile(mesharea,(output.shape[0],1,1))
	373	totalarea = np.ma.masked_where(getmask(output),mesharea)
	374	totalarea = mean(totalarea[:,d2,:],axis=1); totalarea = mean(totalarea[:,d1],axis=1)
	375	except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
[647]	376	output = output*mesharea; output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)/totalarea
[405]	377	#noperturb = smooth1d(output,window_len=7)
	378	#lenlen = len(output) ; output = output[1:lenlen-7] ; yeye = noperturb[4:lenlen-4]
	379	#plot(output) ; plot(yeye) ; show() ; plot(output-yeye) ; show()
[647]	380	elif d3 is not None and d2 is not None and d1 is not None:
	381	output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
[405]	382	if anomaly:
	383	for k in range(output.shape[0]): output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.)
[687]	384	try:
	385	totalarea = np.tile(mesharea,(output.shape[0],1,1))
	386	totalarea = np.ma.masked_where(getmask(output),totalarea)
	387	totalarea = mean(totalarea[:,d2,:],axis=1); totalarea = mean(totalarea[:,d1],axis=1)
	388	except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
[647]	389	output = output*mesharea; output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)/totalarea
[392]	390	elif d4 is not None and d3 is not None:
	391	output = mean(input[d4,:,:,:],axis=0)
	392	output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
[405]	393	if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
[392]	394	elif d4 is not None and d2 is not None:
[647]	395	output = mean(input[d4,:,:,:],axis=0)
[687]	396	try:
	397	totalarea = np.tile(mesharea,(output.shape[0],1,1))
	398	totalarea = np.ma.masked_where(getmask(output),mesharea)
	399	totalarea = mean(totalarea[:,d2,:],axis=1)
	400	except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
[647]	401	output = output*mesharea; output = mean(output[:,d2,:],axis=1)/totalarea
[392]	402	elif d4 is not None and d1 is not None:
	403	output = mean(input[d4,:,:,:],axis=0)
	404	output = mean(output[:,:,d1],axis=2)
[647]	405	elif d3 is not None and d2 is not None:
[392]	406	output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
[687]	407	try:
	408	totalarea = np.tile(mesharea,(output.shape[0],1,1))
	409	totalarea = np.ma.masked_where(getmask(output),mesharea)
	410	totalarea = mean(totalarea[:,d2,:],axis=1)
	411	except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
[647]	412	output = output*mesharea; output = mean(output[:,d2,:],axis=1)/totalarea
[392]	413	elif d3 is not None and d1 is not None:
	414	output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
[448]	415	output = mean(output[:,:,d1],axis=2)
[647]	416	elif d2 is not None and d1 is not None:
[687]	417	try:
	418	totalarea = np.tile(mesharea,(output.shape[0],output.shape[1],1,1))
	419	totalarea = np.ma.masked_where(getmask(output),totalarea)
	420	totalarea = mean(totalarea[:,:,d2,:],axis=2); totalarea = mean(totalarea[:,:,d1],axis=1)
	421	except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
[763]	422	output = output*mesharea; output = mean(output[:,:,d2,:],axis=2); output = mean(output[:,:,d1],axis=2)/totalarea
[392]	423	elif d1 is not None: output = mean(input[:,:,:,d1],axis=3)
[647]	424	elif d2 is not None:
[687]	425	try:
	426	totalarea = np.tile(mesharea,(output.shape[0],output.shape[1],1,output.shape[3]))
	427	totalarea = np.ma.masked_where(getmask(output),totalarea)
	428	totalarea = mean(totalarea[:,:,d2,:],axis=2)
	429	except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
[647]	430	output = output*mesharea; output = mean(output[:,:,d2,:],axis=2)/totalarea
[437]	431	elif d3 is not None: output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
[392]	432	elif d4 is not None: output = mean(input[d4,:,:,:],axis=0)
[468]	433	dimension2 = np.array(output).ndim
	434	shape2 = np.array(output).shape
	435	print 'REDUCEFIELD dim,shape: ',dimension,shape,' >>> ',dimension2,shape2
[233]	436	return output, error
	437
[392]	438	## Author: AC + AS
	439	def reduce_zaxis (input,ax=None,yint=False,vert=None,indice=None):
[382]	440	from mymath import max,mean
	441	from scipy import integrate
[637]	442	import numpy as np
[392]	443	if yint and vert is not None and indice is not None:
[391]	444	if type(input).__name__=='MaskedArray':
	445	input.set_fill_value([np.NaN])
[392]	446	output = integrate.trapz(input.filled(),x=vert[indice],axis=ax)
[391]	447	else:
[396]	448	output = integrate.trapz(input,x=vert[indice],axis=ax)
[382]	449	else:
	450	output = mean(input,axis=ax)
	451	return output
	452
[345]	453	## Author: AS + TN
[233]	454	def definesubplot ( numplot, fig ):
	455	from matplotlib.pyplot import rcParams
	456	rcParams['font.size'] = 12. ## default (important for multiple calls)
[345]	457	if numplot <= 0:
	458	subv = 99999
	459	subh = 99999
	460	elif numplot == 1:
[568]	461	subv = 1
	462	subh = 1
[233]	463	elif numplot == 2:
[483]	464	subv = 1 #2
	465	subh = 2 #1
[233]	466	fig.subplots_adjust(wspace = 0.35)
	467	rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
	468	elif numplot == 3:
[453]	469	subv = 3
	470	subh = 1
[613]	471	fig.subplots_adjust(hspace = 0.5)
[233]	472	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
[637]	473	elif numplot == 4:
[345]	474	subv = 2
	475	subh = 2
[781]	476	#fig.subplots_adjust(wspace = 0.4, hspace = 0.6)
[610]	477	fig.subplots_adjust(wspace = 0.4, hspace = 0.3)
[345]	478	rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
	479	elif numplot <= 6:
	480	subv = 2
	481	subh = 3
[638]	482	#fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
	483	fig.subplots_adjust(wspace = 0.5, hspace = 0.3)
[233]	484	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
[345]	485	elif numplot <= 8:
	486	subv = 2
	487	subh = 4
[233]	488	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
	489	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
[345]	490	elif numplot <= 9:
	491	subv = 3
	492	subh = 3
[233]	493	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
	494	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
[345]	495	elif numplot <= 12:
	496	subv = 3
	497	subh = 4
[453]	498	fig.subplots_adjust(wspace = 0, hspace = 0.1)
[345]	499	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
	500	elif numplot <= 16:
	501	subv = 4
	502	subh = 4
	503	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
	504	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
[233]	505	else:
[345]	506	print "number of plot supported: 1 to 16"
[233]	507	exit()
[345]	508	return subv,subh
[233]	509
[345]	510	## Author: AS
[233]	511	def getstralt(nc,nvert):
[548]	512	varinfile = nc.variables.keys()
	513	if 'vert' not in varinfile:
[233]	514	stralt = "_lvl" + str(nvert)
[548]	515	else:
[233]	516	zelevel = int(nc.variables['vert'][nvert])
	517	if abs(zelevel) < 10000.: strheight=str(zelevel)+"m"
	518	else: strheight=str(int(zelevel/1000.))+"km"
	519	if 'altitude' in nc.dimensions: stralt = "_"+strheight+"-AMR"
	520	elif 'altitude_abg' in nc.dimensions: stralt = "_"+strheight+"-ALS"
	521	elif 'bottom_top' in nc.dimensions: stralt = "_"+strheight
	522	elif 'pressure' in nc.dimensions: stralt = "_"+str(zelevel)+"Pa"
	523	else: stralt = ""
	524	return stralt
	525
[345]	526	## Author: AS
[468]	527	def getlschar ( namefile, getaxis=False ):
[195]	528	from netCDF4 import Dataset
	529	from timestuff import sol2ls
[233]	530	from numpy import array
[400]	531	from string import rstrip
[687]	532	import os as daos
	533	namefiletest = rstrip( rstrip( rstrip( namefile, chars="_z"), chars="_zabg"), chars="_p")
	534	testexist = daos.path.isfile(namefiletest)
[237]	535	zetime = None
[687]	536	if testexist:
	537	namefile = namefiletest
	538	#### we assume that wrfout is next to wrfout_z and wrfout_zabg
	539	nc = Dataset(namefile)
	540	zetime = None
	541	days_in_month = [61, 66, 66, 65, 60, 54, 50, 46, 47, 47, 51, 56]
	542	plus_in_month = [ 0, 61,127,193,258,318,372,422,468,515,562,613]
	543	if 'Times' in nc.variables:
[233]	544	zetime = nc.variables['Times'][0]
	545	shape = array(nc.variables['Times']).shape
	546	if shape[0] < 2: zetime = None
	547	if zetime is not None \
[225]	548	and 'vert' not in nc.variables:
[489]	549	##### strangely enough this does not work for api or ncrcat results!
	550	zesol = plus_in_month[ int(zetime[5]+zetime[6])-1 ] + int(zetime[8]+zetime[9]) - 1 ##les sols GCM commencent a 0
	551	dals = int( 10. * sol2ls ( zesol ) ) / 10.
[197]	552	###
	553	zetime2 = nc.variables['Times'][1]
	554	one = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37.
	555	next = int(zetime2[11]+zetime2[12]) + int(zetime2[14]+zetime2[15])/37.
	556	zehour = one
	557	zehourin = abs ( next - one )
[489]	558	if not getaxis:
	559	lschar = "_Ls"+str(dals)
	560	else:
[468]	561	zelen = len(nc.variables['Times'][:])
	562	yeye = range(zelen) ; lsaxis = range(zelen) ; solaxis = range(zelen) ; ltaxis = range(zelen)
	563	for iii in yeye:
[489]	564	zetime = nc.variables['Times'][iii]
[468]	565	ltaxis[iii] = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37.
[489]	566	solaxis[iii] = ltaxis[iii] / 24. + plus_in_month[ int(zetime[5]+zetime[6])-1 ] + int(zetime[8]+zetime[9]) - 1 ##les sols GCM commencent a 0
[468]	567	lsaxis[iii] = sol2ls ( solaxis[iii] )
	568	if ltaxis[iii] < ltaxis[iii-1]: ltaxis[iii] = ltaxis[iii] + 24.
[489]	569	#print ltaxis[iii], solaxis[iii], lsaxis[iii], getattr( nc, 'JULDAY' )
[468]	570	lschar = lsaxis ; zehour = solaxis ; zehourin = ltaxis
[195]	571	else:
	572	lschar=""
[197]	573	zehour = 0
	574	zehourin = 1
	575	return lschar, zehour, zehourin
[195]	576
[345]	577	## Author: AS
[202]	578	def getprefix (nc):
	579	prefix = 'LMD_MMM_'
	580	prefix = prefix + 'd'+str(getattr(nc,'GRID_ID'))+'_'
	581	prefix = prefix + str(int(getattr(nc,'DX')/1000.))+'km_'
	582	return prefix
	583
[345]	584	## Author: AS
[184]	585	def getproj (nc):
[233]	586	typefile = whatkindfile(nc)
[548]	587	if typefile in ['meso','geo']:
[233]	588	### (il faudrait passer CEN_LON dans la projection ?)
	589	map_proj = getattr(nc, 'MAP_PROJ')
	590	cen_lat = getattr(nc, 'CEN_LAT')
	591	if map_proj == 2:
	592	if cen_lat > 10.:
	593	proj="npstere"
[392]	594	#print "NP stereographic polar domain"
[233]	595	else:
	596	proj="spstere"
[392]	597	#print "SP stereographic polar domain"
[233]	598	elif map_proj == 1:
[392]	599	#print "lambert projection domain"
[233]	600	proj="lcc"
	601	elif map_proj == 3:
[392]	602	#print "mercator projection"
[233]	603	proj="merc"
	604	else:
	605	proj="merc"
[548]	606	elif typefile in ['gcm']: proj="cyl" ## pb avec les autres (de trace derriere la sphere ?)
	607	else: proj="ortho"
[184]	608	return proj
	609
[345]	610	## Author: AS
[180]	611	def ptitle (name):
	612	from matplotlib.pyplot import title
	613	title(name)
	614	print name
	615
[345]	616	## Author: AS
[252]	617	def polarinterv (lon2d,lat2d):
	618	import numpy as np
	619	wlon = [np.min(lon2d),np.max(lon2d)]
	620	ind = np.array(lat2d).shape[0] / 2 ## to get a good boundlat and to get the pole
	621	wlat = [np.min(lat2d[ind,:]),np.max(lat2d[ind,:])]
	622	return [wlon,wlat]
	623
[345]	624	## Author: AS
[180]	625	def simplinterv (lon2d,lat2d):
	626	import numpy as np
	627	return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]]
	628
[345]	629	## Author: AS
[184]	630	def wrfinterv (lon2d,lat2d):
	631	nx = len(lon2d[0,:])-1
	632	ny = len(lon2d[:,0])-1
[225]	633	lon1 = lon2d[0,0]
	634	lon2 = lon2d[nx,ny]
	635	lat1 = lat2d[0,0]
	636	lat2 = lat2d[nx,ny]
[233]	637	if abs(0.5(lat1+lat2)) > 60.: wider = 0.5 (abs(lon1)+abs(lon2)) * 0.1
	638	else: wider = 0.
	639	if lon1 < lon2: wlon = [lon1, lon2 + wider]
[225]	640	else: wlon = [lon2, lon1 + wider]
	641	if lat1 < lat2: wlat = [lat1, lat2]
	642	else: wlat = [lat2, lat1]
	643	return [wlon,wlat]
[184]	644
[345]	645	## Author: AS
[240]	646	def makeplotres (filename,res=None,pad_inches_value=0.25,folder='',disp=True,ext='png',erase=False):
[180]	647	import matplotlib.pyplot as plt
[240]	648	from os import system
	649	addstr = ""
	650	if res is not None:
	651	res = int(res)
	652	addstr = "_"+str(res)
	653	name = filename+addstr+"."+ext
[186]	654	if folder != '': name = folder+'/'+name
[180]	655	plt.savefig(name,dpi=res,bbox_inches='tight',pad_inches=pad_inches_value)
[240]	656	if disp: display(name)
	657	if ext in ['eps','ps','svg']: system("tar czvf "+name+".tar.gz "+name+" ; rm -f "+name)
	658	if erase: system("mv "+name+" to_be_erased")
[180]	659	return
	660
[430]	661	## Author: AS + AC
[451]	662	def dumpbdy (field,n,stag=None,condition=False,onlyx=False,onlyy=False):
[447]	663	nx = len(field[0,:])-1
	664	ny = len(field[:,0])-1
[444]	665	if condition:
	666	if stag == 'U': nx = nx-1
	667	if stag == 'V': ny = ny-1
	668	if stag == 'W': nx = nx+1 #special les case when we dump stag on W
[451]	669	if onlyx: result = field[:,n:nx-n]
	670	elif onlyy: result = field[n:ny-n,:]
	671	else: result = field[n:ny-n,n:nx-n]
	672	return result
[180]	673
[444]	674	## Author: AS + AC
[233]	675	def getcoorddef ( nc ):
[317]	676	import numpy as np
[233]	677	## getcoord2d for predefined types
	678	typefile = whatkindfile(nc)
[548]	679	if typefile in ['meso']:
	680	if '9999' not in getattr(nc,'START_DATE') :
[753]	681	## regular mesoscale
[548]	682	[lon2d,lat2d] = getcoord2d(nc)
	683	else:
	684	## idealized mesoscale
	685	nx=getattr(nc,'WEST-EAST_GRID_DIMENSION')
	686	ny=getattr(nc,'SOUTH-NORTH_GRID_DIMENSION')
	687	dlat=getattr(nc,'DX')
	688	## this is dirty because Martian-specific
	689	# ... but this just intended to get "lat-lon" like info
	690	falselon = np.arange(-dlat(nx-1)/2.,dlat(nx-1)/2.,dlat)/60000.
	691	falselat = np.arange(-dlat(ny-1)/2.,dlat(ny-1)/2.,dlat)/60000.
	692	[lon2d,lat2d] = np.meshgrid(falselon,falselat) ## dummy coordinates
	693	print "WARNING: domain plot artificially centered on lat,lon 0,0"
[637]	694	elif typefile in ['gcm','earthgcm','ecmwf']:
[724]	695	#### n est ce pas nc.variables ?
[637]	696	if "longitude" in nc.dimensions: dalon = "longitude"
	697	elif "lon" in nc.dimensions: dalon = "lon"
[724]	698	else: dalon = "nothing"
[637]	699	if "latitude" in nc.dimensions: dalat = "latitude"
	700	elif "lat" in nc.dimensions: dalat = "lat"
[724]	701	else: dalat = "nothing"
[637]	702	[lon2d,lat2d] = getcoord2d(nc,nlat=dalat,nlon=dalon,is1d=True)
[233]	703	elif typefile in ['geo']:
	704	[lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M')
	705	return lon2d,lat2d
	706
[345]	707	## Author: AS
[184]	708	def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False):
	709	import numpy as np
[724]	710	if nlon == "nothing" or nlat == "nothing":
	711	print "NO LAT LON FIELDS. I AM TRYING MY BEST. I ASSUME GLOBAL FIELD."
	712	lon = np.linspace(-180.,180.,getdimfromvar(nc)[-1])
	713	lat = np.linspace(-90.,90.,getdimfromvar(nc)[-2])
[184]	714	[lon2d,lat2d] = np.meshgrid(lon,lat)
	715	else:
[724]	716	if is1d:
	717	lat = nc.variables[nlat][:]
	718	lon = nc.variables[nlon][:]
	719	[lon2d,lat2d] = np.meshgrid(lon,lat)
	720	else:
	721	lat = nc.variables[nlat][0,:,:]
	722	lon = nc.variables[nlon][0,:,:]
	723	[lon2d,lat2d] = [lon,lat]
[184]	724	return lon2d,lat2d
	725
[724]	726	## Author: AS
	727	def getdimfromvar (nc):
	728	varinfile = nc.variables.keys()
	729	dim = nc.variables[varinfile[-1]].shape ## usually the last variable is 4D or 3D
	730	return dim
	731
[405]	732	## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
	733	def smooth1d(x,window_len=11,window='hanning'):
	734	import numpy
	735	"""smooth the data using a window with requested size.
	736	This method is based on the convolution of a scaled window with the signal.
	737	The signal is prepared by introducing reflected copies of the signal
	738	(with the window size) in both ends so that transient parts are minimized
	739	in the begining and end part of the output signal.
	740	input:
	741	x: the input signal
	742	window_len: the dimension of the smoothing window; should be an odd integer
	743	window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
	744	flat window will produce a moving average smoothing.
	745	output:
	746	the smoothed signal
	747	example:
	748	t=linspace(-2,2,0.1)
	749	x=sin(t)+randn(len(t))*0.1
	750	y=smooth(x)
	751	see also:
	752	numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
	753	scipy.signal.lfilter
	754	TODO: the window parameter could be the window itself if an array instead of a string
	755	"""
	756	x = numpy.array(x)
	757	if x.ndim != 1:
	758	raise ValueError, "smooth only accepts 1 dimension arrays."
	759	if x.size < window_len:
	760	raise ValueError, "Input vector needs to be bigger than window size."
	761	if window_len<3:
	762	return x
	763	if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
	764	raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
	765	s=numpy.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]
	766	#print(len(s))
	767	if window == 'flat': #moving average
	768	w=numpy.ones(window_len,'d')
	769	else:
	770	w=eval('numpy.'+window+'(window_len)')
	771	y=numpy.convolve(w/w.sum(),s,mode='valid')
	772	return y
	773
[345]	774	## Author: AS
[180]	775	def smooth (field, coeff):
	776	## actually blur_image could work with different coeff on x and y
	777	if coeff > 1: result = blur_image(field,int(coeff))
	778	else: result = field
	779	return result
	780
[345]	781	## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
[180]	782	def gauss_kern(size, sizey=None):
	783	import numpy as np
	784	# Returns a normalized 2D gauss kernel array for convolutions
	785	size = int(size)
	786	if not sizey:
	787	sizey = size
	788	else:
	789	sizey = int(sizey)
	790	x, y = np.mgrid[-size:size+1, -sizey:sizey+1]
	791	g = np.exp(-(x2/float(size)+y2/float(sizey)))
	792	return g / g.sum()
	793
[345]	794	## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
[180]	795	def blur_image(im, n, ny=None) :
	796	from scipy.signal import convolve
	797	# blurs the image by convolving with a gaussian kernel of typical size n.
	798	# The optional keyword argument ny allows for a different size in the y direction.
	799	g = gauss_kern(n, sizey=ny)
	800	improc = convolve(im, g, mode='same')
	801	return improc
	802
[345]	803	## Author: AS
[233]	804	def getwinddef (nc):
	805	###
[548]	806	varinfile = nc.variables.keys()
	807	if 'Um' in varinfile: [uchar,vchar] = ['Um','Vm'] #; print "this is API meso file"
	808	elif 'U' in varinfile: [uchar,vchar] = ['U','V'] #; print "this is RAW meso file"
	809	elif 'u' in varinfile: [uchar,vchar] = ['u','v'] #; print "this is GCM file"
[721]	810	elif 'vitu' in varinfile: [uchar,vchar] = ['vitu','vitv'] #; print "this is GCM v5 file"
[548]	811	### you can add choices here !
	812	else: [uchar,vchar] = ['not found','not found']
[233]	813	###
[548]	814	if uchar in ['U']: metwind = False ## geometrical (wrt grid)
	815	else: metwind = True ## meteorological (zon/mer)
	816	if metwind is False: print "Not using meteorological winds. You trust numerical grid as being (x,y)"
[233]	817	###
	818	return uchar,vchar,metwind
[202]	819
[345]	820	## Author: AS
[184]	821	def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True):
	822	## scale regle la reference du vecteur
	823	## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs.
	824	import matplotlib.pyplot as plt
	825	import numpy as np
[638]	826	#posx = np.min(x) - np.std(x) / 10.
	827	#posy = np.min(y) - np.std(y) / 10.
	828	posx = np.min(x)
	829	posy = np.min(y) - 4.*np.std(y) / 10.
[184]	830	u = smooth(u,csmooth)
	831	v = smooth(v,csmooth)
[188]	832	widthvec = 0.003 #0.005 #0.003
[184]	833	q = plt.quiver( x[::stride,::stride],\
	834	y[::stride,::stride],\
	835	u[::stride,::stride],\
	836	v[::stride,::stride],\
[228]	837	angles='xy',color=color,pivot='middle',\
[184]	838	scale=factor,width=widthvec )
[202]	839	if color in ['white','yellow']: kcolor='black'
	840	else: kcolor=color
[184]	841	if key: p = plt.quiverkey(q,posx,posy,scale,\
[194]	842	str(int(scale)),coordinates='data',color=kcolor,labelpos='S',labelsep = 0.03)
[184]	843	return
[180]	844
[345]	845	## Author: AS
[180]	846	def display (name):
[184]	847	from os import system
	848	system("display "+name+" > /dev/null 2> /dev/null &")
	849	return name
[180]	850
[345]	851	## Author: AS
[180]	852	def findstep (wlon):
[184]	853	steplon = int((wlon[1]-wlon[0])/4.) #3
	854	step = 120.
	855	while step > steplon and step > 15. : step = step / 2.
	856	if step <= 15.:
	857	while step > steplon and step > 5. : step = step - 5.
	858	if step <= 5.:
	859	while step > steplon and step > 1. : step = step - 1.
	860	if step <= 1.:
	861	step = 1.
[180]	862	return step
	863
[345]	864	## Author: AS
[451]	865	def define_proj (char,wlon,wlat,back=None,blat=None,blon=None):
[180]	866	from mpl_toolkits.basemap import Basemap
	867	import numpy as np
	868	import matplotlib as mpl
[240]	869	from mymath import max
[180]	870	meanlon = 0.5*(wlon[0]+wlon[1])
	871	meanlat = 0.5*(wlat[0]+wlat[1])
[637]	872	zewidth = np.abs(wlon[0]-wlon[1])60000.np.cos(3.14*meanlat/180.)
	873	zeheight = np.abs(wlat[0]-wlat[1])*60000.
[385]	874	if blat is None:
[398]	875	ortholat=meanlat
[453]	876	if wlat[0] >= 80.: blat = -40.
[345]	877	elif wlat[1] <= -80.: blat = -40.
	878	elif wlat[1] >= 0.: blat = wlat[0]
	879	elif wlat[0] <= 0.: blat = wlat[1]
[398]	880	else: ortholat=blat
[451]	881	if blon is None: ortholon=meanlon
	882	else: ortholon=blon
[207]	883	h = 50. ## en km
[202]	884	radius = 3397200.
[637]	885	#print meanlat, meanlon
[184]	886	if char == "cyl": m = Basemap(rsphere=radius,projection='cyl',\
[180]	887	llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
[184]	888	elif char == "moll": m = Basemap(rsphere=radius,projection='moll',lon_0=meanlon)
[451]	889	elif char == "ortho": m = Basemap(rsphere=radius,projection='ortho',lon_0=ortholon,lat_0=ortholat)
[184]	890	elif char == "lcc": m = Basemap(rsphere=radius,projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\
[637]	891	width=zewidth,height=zeheight)
	892	#llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
[184]	893	elif char == "npstere": m = Basemap(rsphere=radius,projection='npstere', boundinglat=blat, lon_0=0.)
[395]	894	elif char == "spstere": m = Basemap(rsphere=radius,projection='spstere', boundinglat=blat, lon_0=180.)
[207]	895	elif char == "nplaea": m = Basemap(rsphere=radius,projection='nplaea', boundinglat=wlat[0], lon_0=meanlon)
	896	elif char == "laea": m = Basemap(rsphere=radius,projection='laea',lon_0=meanlon,lat_0=meanlat,lat_ts=meanlat,\
[637]	897	width=zewidth,height=zeheight)
	898	#llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
[184]	899	elif char == "nsper": m = Basemap(rsphere=radius,projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=h*1000.)
	900	elif char == "merc": m = Basemap(rsphere=radius,projection='merc',lat_ts=0.,\
	901	llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
	902	fontsizemer = int(mpl.rcParams['font.size']*3./4.)
[207]	903	if char in ["cyl","lcc","merc","nsper","laea"]: step = findstep(wlon)
	904	else: step = 10.
[238]	905	steplon = step*2.
[453]	906	zecolor ='grey'
	907	zelinewidth = 1
	908	zelatmax = 80
[637]	909	if meanlat > 75.: zelatmax = 90. ; step = step/2.
[453]	910	# to show gcm grid:
	911	#zecolor = 'r'
	912	#zelinewidth = 1
[647]	913	#step = 180./48.
	914	#steplon = 360./64.
	915	#zelatmax = 90. - step/3
[516]	916	if char not in ["moll"]:
[760]	917	if wlon[1]-wlon[0] < 2.: ## LOCAL MODE
	918	m.drawmeridians(np.r_[-1.:1.:0.05], labels=[0,0,0,1], color=zecolor, linewidth=zelinewidth, fontsize=fontsizemer, fmt='%5.2f')
	919	m.drawparallels(np.r_[-1.:1.:0.05], labels=[1,0,0,0], color=zecolor, linewidth=zelinewidth, fontsize=fontsizemer, fmt='%5.2f')
	920	else: ## GLOBAL OR REGIONAL MODE
	921	m.drawmeridians(np.r_[-180.:180.:steplon], labels=[0,0,0,1], color=zecolor, linewidth=zelinewidth, fontsize=fontsizemer, latmax=zelatmax)
	922	m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color=zecolor, linewidth=zelinewidth, fontsize=fontsizemer, latmax=zelatmax)
[233]	923	if back: m.warpimage(marsmap(back),scale=0.75)
	924	#if not back:
	925	# if not var: back = "mola" ## if no var: draw mola
	926	# elif typefile in ['mesoapi','meso','geo'] \
	927	# and proj not in ['merc','lcc','nsper','laea']: back = "molabw" ## if var but meso: draw molabw
	928	# else: pass ## else: draw None
[180]	929	return m
	930
[345]	931	## Author: AS
[232]	932	#### test temporaire
	933	def putpoints (map,plot):
	934	#### from http://www.scipy.org/Cookbook/Matplotlib/Maps
	935	# lat/lon coordinates of five cities.
	936	lats = [18.4]
	937	lons = [-134.0]
	938	points=['Olympus Mons']
	939	# compute the native map projection coordinates for cities.
	940	x,y = map(lons,lats)
	941	# plot filled circles at the locations of the cities.
	942	map.plot(x,y,'bo')
	943	# plot the names of those five cities.
	944	wherept = 0 #1000 #50000
	945	for name,xpt,ypt in zip(points,x,y):
	946	plot.text(xpt+wherept,ypt+wherept,name)
	947	## le nom ne s'affiche pas...
	948	return
	949
[345]	950	## Author: AS
[233]	951	def calculate_bounds(field,vmin=None,vmax=None):
	952	import numpy as np
	953	from mymath import max,min,mean
	954	ind = np.where(field < 9e+35)
	955	fieldcalc = field[ ind ] # la syntaxe compacte ne marche si field est un tuple
	956	###
	957	dev = np.std(fieldcalc)*3.0
	958	###
[562]	959	if vmin is None: zevmin = mean(fieldcalc) - dev
[233]	960	else: zevmin = vmin
	961	###
	962	if vmax is None: zevmax = mean(fieldcalc) + dev
	963	else: zevmax = vmax
	964	if vmin == vmax:
	965	zevmin = mean(fieldcalc) - dev ### for continuity
	966	zevmax = mean(fieldcalc) + dev ### for continuity
	967	###
	968	if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0.
[468]	969	print "BOUNDS field ", min(fieldcalc), max(fieldcalc), " //// adopted", zevmin, zevmax
[233]	970	return zevmin, zevmax
[232]	971
[345]	972	## Author: AS
[233]	973	def bounds(what_I_plot,zevmin,zevmax):
[247]	974	from mymath import max,min,mean
[233]	975	### might be convenient to add the missing value in arguments
[310]	976	#what_I_plot[ what_I_plot < zevmin ] = zevmin#*(1. + 1.e-7)
	977	if zevmin < 0: what_I_plot[ what_I_plot < zevmin(1. - 1.e-7) ] = zevmin(1. - 1.e-7)
	978	else: what_I_plot[ what_I_plot < zevmin(1. + 1.e-7) ] = zevmin(1. + 1.e-7)
[451]	979	#print "NEW MIN ", min(what_I_plot)
[233]	980	what_I_plot[ what_I_plot > 9e+35 ] = -9e+35
[587]	981	what_I_plot[ what_I_plot > zevmax ] = zevmax*(1. - 1.e-7)
[451]	982	#print "NEW MAX ", max(what_I_plot)
[233]	983	return what_I_plot
	984
[345]	985	## Author: AS
[241]	986	def nolow(what_I_plot):
	987	from mymath import max,min
	988	lim = 0.150.5(abs(max(what_I_plot))+abs(min(what_I_plot)))
[392]	989	print "NO PLOT BELOW VALUE ", lim
[241]	990	what_I_plot [ abs(what_I_plot) < lim ] = 1.e40
	991	return what_I_plot
	992
[418]	993
	994	## Author : AC
	995	def hole_bounds(what_I_plot,zevmin,zevmax):
	996	import numpy as np
	997	zi=0
	998	for i in what_I_plot:
	999	zj=0
	1000	for j in i:
	1001	if ((j < zevmin) or (j > zevmax)):what_I_plot[zi,zj]=np.NaN
	1002	zj=zj+1
	1003	zi=zi+1
	1004
	1005	return what_I_plot
	1006
[345]	1007	## Author: AS
[233]	1008	def zoomset (wlon,wlat,zoom):
	1009	dlon = abs(wlon[1]-wlon[0])/2.
	1010	dlat = abs(wlat[1]-wlat[0])/2.
	1011	[wlon,wlat] = [ [wlon[0]+zoomdlon/100.,wlon[1]-zoomdlon/100.],\
	1012	[wlat[0]+zoomdlat/100.,wlat[1]-zoomdlat/100.] ]
[392]	1013	print "ZOOM %",zoom,wlon,wlat
[233]	1014	return wlon,wlat
	1015
[345]	1016	## Author: AS
[201]	1017	def fmtvar (whichvar="def"):
[204]	1018	fmtvar = { \
[502]	1019	"MIXED": "%.0f",\
	1020	"UPDRAFT": "%.0f",\
	1021	"DOWNDRAFT": "%.0f",\
[405]	1022	"TK": "%.0f",\
[637]	1023	"T": "%.0f",\
[516]	1024	#"ZMAX_TH": "%.0f",\
	1025	#"WSTAR": "%.0f",\
[425]	1026	# Variables from TES ncdf format
[363]	1027	"T_NADIR_DAY": "%.0f",\
[376]	1028	"T_NADIR_NIT": "%.0f",\
[425]	1029	# Variables from tes.py ncdf format
[398]	1030	"TEMP_DAY": "%.0f",\
	1031	"TEMP_NIGHT": "%.0f",\
[425]	1032	# Variables from MCS and mcs.py ncdf format
[427]	1033	"DTEMP": "%.0f",\
	1034	"NTEMP": "%.0f",\
	1035	"DNUMBINTEMP": "%.0f",\
	1036	"NNUMBINTEMP": "%.0f",\
[425]	1037	# other stuff
[405]	1038	"TPOT": "%.0f",\
[295]	1039	"TSURF": "%.0f",\
[612]	1040	"U_OUT1": "%.0f",\
	1041	"T_OUT1": "%.0f",\
[204]	1042	"def": "%.1e",\
	1043	"PTOT": "%.0f",\
[760]	1044	"PSFC": "%.1f",\
[204]	1045	"HGT": "%.1e",\
	1046	"USTM": "%.2f",\
[225]	1047	"HFX": "%.0f",\
[232]	1048	"ICETOT": "%.1e",\
[237]	1049	"TAU_ICE": "%.2f",\
[451]	1050	"TAUICE": "%.2f",\
[252]	1051	"VMR_ICE": "%.1e",\
[345]	1052	"MTOT": "%.1f",\
[405]	1053	"ANOMALY": "%.1f",\
[771]	1054	"W": "%.2f",\
[287]	1055	"WMAX_TH": "%.1f",\
[562]	1056	"WSTAR": "%.1f",\
[287]	1057	"QSURFICE": "%.0f",\
[405]	1058	"UM": "%.0f",\
[490]	1059	"WIND": "%.0f",\
[612]	1060	"UVMET": "%.0f",\
	1061	"UV": "%.0f",\
[295]	1062	"ALBBARE": "%.2f",\
[317]	1063	"TAU": "%.1f",\
[382]	1064	"CO2": "%.2f",\
[701]	1065	"ENFACT": "%.1f",\
[771]	1066	"QDUST": "%.6f",\
[345]	1067	#### T.N.
	1068	"TEMP": "%.0f",\
	1069	"VMR_H2OICE": "%.0f",\
	1070	"VMR_H2OVAP": "%.0f",\
	1071	"TAUTES": "%.2f",\
	1072	"TAUTESAP": "%.2f",\
	1073
[204]	1074	}
[518]	1075	if "TSURF" in whichvar: whichvar = "TSURF"
[204]	1076	if whichvar not in fmtvar:
	1077	whichvar = "def"
	1078	return fmtvar[whichvar]
[201]	1079
[345]	1080	## Author: AS
[233]	1081	####################################################################################################################
	1082	### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
[202]	1083	def defcolorb (whichone="def"):
[204]	1084	whichcolorb = { \
	1085	"def": "spectral",\
	1086	"HGT": "spectral",\
[426]	1087	"HGT_M": "spectral",\
[405]	1088	"TK": "gist_heat",\
[425]	1089	"TPOT": "Paired",\
[295]	1090	"TSURF": "RdBu_r",\
[204]	1091	"QH2O": "PuBu",\
	1092	"USTM": "YlOrRd",\
[490]	1093	"WIND": "YlOrRd",\
[363]	1094	#"T_nadir_nit": "RdBu_r",\
	1095	#"T_nadir_day": "RdBu_r",\
[225]	1096	"HFX": "RdYlBu",\
[310]	1097	"ICETOT": "YlGnBu_r",\
[345]	1098	#"MTOT": "PuBu",\
	1099	"CCNQ": "YlOrBr",\
	1100	"CCNN": "YlOrBr",\
	1101	"TEMP": "Jet",\
[238]	1102	"TAU_ICE": "Blues",\
[451]	1103	"TAUICE": "Blues",\
[252]	1104	"VMR_ICE": "Blues",\
[241]	1105	"W": "jet",\
[287]	1106	"WMAX_TH": "spectral",\
[405]	1107	"ANOMALY": "RdBu_r",\
[287]	1108	"QSURFICE": "hot_r",\
[295]	1109	"ALBBARE": "spectral",\
[317]	1110	"TAU": "YlOrBr_r",\
[382]	1111	"CO2": "YlOrBr_r",\
[753]	1112	"MIXED": "GnBu",\
[345]	1113	#### T.N.
[647]	1114	"MTOT": "spectral",\
[345]	1115	"H2O_ICE_S": "RdBu",\
	1116	"VMR_H2OICE": "PuBu",\
	1117	"VMR_H2OVAP": "PuBu",\
[453]	1118	"WATERCAPTAG": "Blues",\
[204]	1119	}
[241]	1120	#W --> spectral ou jet
[240]	1121	#spectral BrBG RdBu_r
[392]	1122	#print "predefined colorbars"
[518]	1123	if "TSURF" in whichone: whichone = "TSURF"
[204]	1124	if whichone not in whichcolorb:
	1125	whichone = "def"
	1126	return whichcolorb[whichone]
[202]	1127
[345]	1128	## Author: AS
[202]	1129	def definecolorvec (whichone="def"):
	1130	whichcolor = { \
	1131	"def": "black",\
	1132	"vis": "yellow",\
[781]	1133	"vishires": "green",\
[202]	1134	"molabw": "yellow",\
	1135	"mola": "black",\
	1136	"gist_heat": "white",\
	1137	"hot": "tk",\
	1138	"gist_rainbow": "black",\
	1139	"spectral": "black",\
	1140	"gray": "red",\
	1141	"PuBu": "black",\
[721]	1142	"titan": "red",\
[202]	1143	}
	1144	if whichone not in whichcolor:
	1145	whichone = "def"
	1146	return whichcolor[whichone]
	1147
[345]	1148	## Author: AS
[180]	1149	def marsmap (whichone="vishires"):
[233]	1150	from os import uname
	1151	mymachine = uname()[1]
	1152	### not sure about speed-up with this method... looks the same
[511]	1153	if "lmd.jussieu.fr" in mymachine: domain = "/u/aslmd/WWW/maps/"
	1154	elif "aymeric-laptop" in mymachine: domain = "/home/aymeric/Dropbox/Public/"
	1155	else: domain = "http://www.lmd.jussieu.fr/~aslmd/maps/"
[180]	1156	whichlink = { \
[233]	1157	#"vis": "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\
	1158	#"vishires": "http://www.lmd.jussieu.fr/~aslmd/maps/MarsMap_2500x1250.jpg",\
	1159	#"geolocal": "http://dl.dropbox.com/u/11078310/geolocal.jpg",\
	1160	#"mola": "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\
	1161	#"molabw": "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\
[453]	1162	"thermalday": domain+"thermalday.jpg",\
	1163	"thermalnight": domain+"thermalnight.jpg",\
	1164	"tesalbedo": domain+"tesalbedo.jpg",\
[233]	1165	"vis": domain+"mar0kuu2.jpg",\
	1166	"vishires": domain+"MarsMap_2500x1250.jpg",\
	1167	"geolocal": domain+"geolocal.jpg",\
	1168	"mola": domain+"mars-mola-2k.jpg",\
	1169	"molabw": domain+"MarsElevation_2500x1250.jpg",\
[238]	1170	"clouds": "http://www.johnstonsarchive.net/spaceart/marswcloudmap.jpg",\
	1171	"jupiter": "http://www.mmedia.is/~bjj/data/jupiter_css/jupiter_css.jpg",\
	1172	"jupiter_voy": "http://www.mmedia.is/~bjj/data/jupiter/jupiter_vgr2.jpg",\
[558]	1173	#"bw": domain+"EarthElevation_2500x1250.jpg",\
[273]	1174	"bw": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg",\
	1175	"contrast": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg",\
	1176	"nice": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg",\
	1177	"blue": "http://eoimages.gsfc.nasa.gov/ve/2430/land_ocean_ice_2048.jpg",\
[296]	1178	"blueclouds": "http://eoimages.gsfc.nasa.gov/ve/2431/land_ocean_ice_cloud_2048.jpg",\
	1179	"justclouds": "http://eoimages.gsfc.nasa.gov/ve/2432/cloud_combined_2048.jpg",\
[721]	1180	"pluto": "http://www.boulder.swri.edu/~buie/pluto/pluto_all.png",\
	1181	"triton": "http://laps.noaa.gov/albers/sos/neptune/triton/triton_rgb_cyl_www.jpg",\
	1182	"titan": "http://laps.noaa.gov/albers/sos/saturn/titan/titan_rgb_cyl_www.jpg",\
	1183	#"titan": "http://laps.noaa.gov/albers/sos/celestia/titan_50.jpg",\
	1184	"titanuni": "http://maps.jpl.nasa.gov/pix/sat6fss1.jpg",\
	1185	"venus": "http://laps.noaa.gov/albers/sos/venus/venus4/venus4_rgb_cyl_www.jpg",\
	1186	"cosmic": "http://laps.noaa.gov/albers/sos/universe/wmap/wmap_rgb_cyl_www.jpg",\
[180]	1187	}
[238]	1188	### see http://www.mmedia.is/~bjj/planetary_maps.html
[180]	1189	if whichone not in whichlink:
	1190	print "marsmap: choice not defined... you'll get the default one... "
	1191	whichone = "vishires"
	1192	return whichlink[whichone]
	1193
[273]	1194	#def earthmap (whichone):
	1195	# if whichone == "contrast": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg"
	1196	# elif whichone == "bw": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg"
	1197	# elif whichone == "nice": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg"
	1198	# return whichlink
[180]	1199
[345]	1200	## Author: AS
[241]	1201	def latinterv (area="Whole"):
	1202	list = { \
	1203	"Europe": [[ 20., 80.],[- 50., 50.]],\
	1204	"Central_America": [[-10., 40.],[ 230., 300.]],\
	1205	"Africa": [[-20., 50.],[- 50., 50.]],\
[273]	1206	"Whole": [[-90., 90.],[-180., 180.]],\
	1207	"Southern_Hemisphere": [[-90., 60.],[-180., 180.]],\
	1208	"Northern_Hemisphere": [[-60., 90.],[-180., 180.]],\
[241]	1209	"Tharsis": [[-30., 60.],[-170.,- 10.]],\
	1210	"Whole_No_High": [[-60., 60.],[-180., 180.]],\
	1211	"Chryse": [[-60., 60.],[- 60., 60.]],\
	1212	"North_Pole": [[ 50., 90.],[-180., 180.]],\
	1213	"Close_North_Pole": [[ 75., 90.],[-180., 180.]],\
	1214	"Far_South_Pole": [[-90.,-40.],[-180., 180.]],\
	1215	"South_Pole": [[-90.,-50.],[-180., 180.]],\
	1216	"Close_South_Pole": [[-90.,-75.],[-180., 180.]],\
[637]	1217	"Sirenum_Crater_large": [[-46.,-34.],[-166.,-151.]],\
	1218	"Sirenum_Crater_small": [[-36.,-26.],[-168.,-156.]],\
	1219	"Rupes": [[ 72., 90.],[-120.,- 20.]],\
[721]	1220	"Xanadu": [[-40., 20.],[ 40., 120.]],\
[781]	1221	"Hyperboreae": [[ 80., 87.],[- 70.,- 10.]],\
[241]	1222	}
	1223	if area not in list: area = "Whole"
	1224	[olat,olon] = list[area]
	1225	return olon,olat
	1226
[345]	1227	## Author: TN
	1228	def separatenames (name):
	1229	from numpy import concatenate
	1230	# look for comas in the input name to separate different names (files, variables,etc ..)
	1231	if name is None:
	1232	names = None
	1233	else:
	1234	names = []
	1235	stop = 0
	1236	currentname = name
	1237	while stop == 0:
	1238	indexvir = currentname.find(',')
	1239	if indexvir == -1:
	1240	stop = 1
	1241	name1 = currentname
	1242	else:
	1243	name1 = currentname[0:indexvir]
	1244	names = concatenate((names,[name1]))
	1245	currentname = currentname[indexvir+1:len(currentname)]
	1246	return names
	1247
	1248
	1249	## Author: TN
	1250	def readslices(saxis):
	1251	from numpy import empty
	1252	if saxis == None:
	1253	zesaxis = None
	1254	else:
	1255	zesaxis = empty((len(saxis),2))
	1256	for i in range(len(saxis)):
	1257	a = separatenames(saxis[i])
	1258	if len(a) == 1:
	1259	zesaxis[i,:] = float(a[0])
	1260	else:
	1261	zesaxis[i,0] = float(a[0])
	1262	zesaxis[i,1] = float(a[1])
	1263
	1264	return zesaxis
	1265
[568]	1266	## Author: TN
	1267	def readdata(data,datatype,coord1,coord2):
	1268	## Read sparse data
	1269	from numpy import empty
[572]	1270	if datatype == 'txt':
[568]	1271	if len(data[coord1].shape) == 1:
	1272	return data[coord1][:]
	1273	elif len(data[coord1].shape) == 2:
	1274	return data[coord1][:,int(coord2)-1]
	1275	else:
	1276	errormess('error in readdata')
[572]	1277	elif datatype == 'sav':
[568]	1278	return data[coord1][coord2]
	1279	else:
	1280	errormess(datatype+' type is not supported!')
	1281
	1282
[399]	1283	## Author: AS
[475]	1284	def bidimfind(lon2d,lat2d,vlon,vlat,file=None):
[399]	1285	import numpy as np
[475]	1286	import matplotlib.pyplot as mpl
[399]	1287	if vlat is None: array = (lon2d - vlon)**2
	1288	elif vlon is None: array = (lat2d - vlat)**2
	1289	else: array = (lon2d - vlon)2 + (lat2d - vlat)2
	1290	idy,idx = np.unravel_index( np.argmin(array), lon2d.shape )
	1291	if vlon is not None:
[475]	1292	if (np.abs(lon2d[idy,idx]-vlon)) > 5: errormess("longitude not found ",printvar=lon2d)
[399]	1293	if vlat is not None:
[475]	1294	if (np.abs(lat2d[idy,idx]-vlat)) > 5: errormess("latitude not found ",printvar=lat2d)
	1295	if file is not None:
	1296	print idx,idy,lon2d[idy,idx],vlon
	1297	print idx,idy,lat2d[idy,idx],vlat
	1298	var = file.variables["HGT"][:,:,:]
[489]	1299	mpl.contourf(var[0,:,:],30,cmap = mpl.get_cmap(name="Greys_r") ) ; mpl.axis('off') ; mpl.plot(idx,idy,'mx',mew=4.0,ms=20.0)
[475]	1300	mpl.show()
	1301	return idy,idx
[399]	1302
[345]	1303	## Author: TN
[399]	1304	def getsindex(saxis,index,axis):
[345]	1305	# input : all the desired slices and the good index
	1306	# output : all indexes to be taken into account for reducing field
	1307	import numpy as np
[349]	1308	if ( np.array(axis).ndim == 2):
	1309	axis = axis[:,0]
[345]	1310	if saxis is None:
	1311	zeindex = None
	1312	else:
	1313	aaa = int(np.argmin(abs(saxis[index,0] - axis)))
	1314	bbb = int(np.argmin(abs(saxis[index,1] - axis)))
	1315	[imin,imax] = np.sort(np.array([aaa,bbb]))
	1316	zeindex = np.array(range(imax-imin+1))+imin
	1317	# because -180 and 180 are the same point in longitude,
	1318	# we get rid of one for averaging purposes.
	1319	if axis[imin] == -180 and axis[imax] == 180:
	1320	zeindex = zeindex[0:len(zeindex)-1]
[392]	1321	print "INFO: whole longitude averaging asked, so last point is not taken into account."
[345]	1322	return zeindex
	1323
	1324	## Author: TN
[763]	1325	def define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,indextime,what_I_plot,dim0,vertmode,redope):
[345]	1326	# Purpose of define_axis is to find x and y axis scales in a smart way
	1327	# x axis priority: 1/time 2/lon 3/lat 4/vertical
	1328	# To be improved !!!...
	1329	from numpy import array,swapaxes
	1330	x = None
	1331	y = None
	1332	count = 0
	1333	what_I_plot = array(what_I_plot)
	1334	shape = what_I_plot.shape
[477]	1335	if indextime is None and len(time) > 1:
[392]	1336	print "AXIS is time"
[345]	1337	x = time
	1338	count = count+1
[763]	1339	if indexlon is None and len(lon) > 1 and redope not in ['edge_x1','edge_x2']:
[392]	1340	print "AXIS is lon"
[345]	1341	if count == 0: x = lon
	1342	else: y = lon
	1343	count = count+1
[763]	1344	if indexlat is None and len(lat) > 1 and redope not in ['edge_y1','edge_y2']:
[392]	1345	print "AXIS is lat"
[345]	1346	if count == 0: x = lat
	1347	else: y = lat
	1348	count = count+1
[579]	1349	if indexvert is None and ((dim0 == 4) or (y is None)):
[392]	1350	print "AXIS is vert"
[345]	1351	if vertmode == 0: # vertical axis is as is (GCM grid)
	1352	if count == 0: x=range(len(vert))
	1353	else: y=range(len(vert))
	1354	count = count+1
	1355	else: # vertical axis is in kms
	1356	if count == 0: x = vert
	1357	else: y = vert
	1358	count = count+1
	1359	x = array(x)
	1360	y = array(y)
[468]	1361	print "CHECK SHAPE: what_I_plot, x, y", what_I_plot.shape, x.shape, y.shape
[345]	1362	if len(shape) == 1:
[562]	1363	if shape[0] != len(x): print "WARNING: shape[0] != len(x). Correcting." ; what_I_plot = what_I_plot[0:len(x)]
[579]	1364	if len(y.shape) > 0: y = ()
[345]	1365	elif len(shape) == 2:
[562]	1366	if shape[1] == len(y) and shape[0] == len(x) and shape[0] != shape[1]:
	1367	print "INFO: swapaxes: ",what_I_plot.shape,shape ; what_I_plot = swapaxes(what_I_plot,0,1)
	1368	else:
	1369	if shape[0] != len(y): print "WARNING: shape[0] != len(y). Correcting." ; what_I_plot = what_I_plot[0:len(y),:]
	1370	elif shape[1] != len(x): print "WARNING: shape[1] != len(x). Correcting." ; what_I_plot = what_I_plot[:,0:len(x)]
	1371	elif len(shape) == 3:
	1372	if vertmode < 0: print "not supported. must check array dimensions at some point. not difficult to implement though."
[345]	1373	return what_I_plot,x,y
[349]	1374
[763]	1375	# Author: TN + AS + AC
	1376	def determineplot(slon, slat, svert, stime, redope):
[349]	1377	nlon = 1 # number of longitudinal slices -- 1 is None
	1378	nlat = 1
	1379	nvert = 1
	1380	ntime = 1
	1381	nslices = 1
	1382	if slon is not None:
[770]	1383	length=len(slon[:,0])
[763]	1384	nslices = nslices*length
[349]	1385	nlon = len(slon)
	1386	if slat is not None:
[770]	1387	length=len(slat[:,0])
[763]	1388	nslices = nslices*length
[349]	1389	nlat = len(slat)
	1390	if svert is not None:
[771]	1391	length=len(svert[:,0])
[763]	1392	nslices = nslices*length
[349]	1393	nvert = len(svert)
	1394	if stime is not None:
[770]	1395	length=len(stime[:,0])
[763]	1396	nslices = nslices*length
[349]	1397	ntime = len(stime)
	1398	#else:
	1399	# nslices = 2
	1400	mapmode = 0
[763]	1401	if slon is None and slat is None and redope not in ['edge_x1','edge_x2','edge_y1','edge_y2']:
[349]	1402	mapmode = 1 # in this case we plot a map, with the given projection
	1403	return nlon, nlat, nvert, ntime, mapmode, nslices
[440]	1404
[638]	1405	## Author : AS
	1406	def maplatlon( lon,lat,field,\
	1407	proj="cyl",colorb="jet",ndiv=10,zeback="molabw",trans=0.6,title="",\
	1408	vecx=None,vecy=None,stride=2 ):
	1409	### an easy way to map a field over lat/lon grid
	1410	import numpy as np
	1411	import matplotlib.pyplot as mpl
	1412	from matplotlib.cm import get_cmap
	1413	## get lon and lat in 2D version. get lat/lon intervals
	1414	numdim = len(np.array(lon).shape)
	1415	if numdim == 2: [lon2d,lat2d] = [lon,lat]
	1416	elif numdim == 1: [lon2d,lat2d] = np.meshgrid(lon,lat)
	1417	else: errormess("lon and lat arrays must be 1D or 2D")
	1418	[wlon,wlat] = latinterv()
	1419	## define projection and background. define x and y given the projection
	1420	m = define_proj(proj,wlon,wlat,back=zeback,blat=None,blon=None)
	1421	x, y = m(lon2d, lat2d)
	1422	## define field. bound field.
	1423	what_I_plot = np.transpose(field)
	1424	zevmin, zevmax = calculate_bounds(what_I_plot) ## vmin=min(what_I_plot_frame), vmax=max(what_I_plot_frame))
	1425	what_I_plot = bounds(what_I_plot,zevmin,zevmax)
	1426	## define contour field levels. define color palette
	1427	ticks = ndiv + 1
	1428	zelevels = np.linspace(zevmin,zevmax,ticks)
	1429	palette = get_cmap(name=colorb)
	1430	## contour field
	1431	m.contourf( x, y, what_I_plot, zelevels, cmap = palette, alpha = trans )
	1432	## draw colorbar
	1433	if proj in ['moll','cyl']: zeorientation="horizontal" ; zepad = 0.07
	1434	else: zeorientation="vertical" ; zepad = 0.03
	1435	#daformat = fmtvar(fvar.upper())
	1436	daformat = "%.0f"
	1437	zecb = mpl.colorbar( fraction=0.05,pad=zepad,format=daformat,orientation=zeorientation,\
	1438	ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,21])),extend='neither',spacing='proportional' )
	1439	## give a title
	1440	if zeorientation == "horizontal": zecb.ax.set_xlabel(title)
	1441	else: ptitle(title)
	1442	## draw vector
	1443	if vecx is not None and vecy is not None:
	1444	[vecx_frame,vecy_frame] = m.rotate_vector( np.transpose(vecx), np.transpose(vecy), lon2d, lat2d ) ## for metwinds
	1445	vectorfield(vecx_frame, vecy_frame, x, y, stride=stride, csmooth=2,\
	1446	scale=30., factor=500., color=definecolorvec(colorb), key=True)
	1447	## scale regle la reference du vecteur. factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs.
	1448	return
[754]	1449	## Author : AC
	1450	## Handles calls to specific computations (e.g. wind norm, enrichment factor...)
[763]	1451	def select_getfield(zvarname=None,znc=None,ztypefile=None,mode=None,ztsat=None,ylon=None,ylat=None,yalt=None,ytime=None,analysis=None):
[754]	1452	from mymath import get_tsat
	1453
	1454	## Specific variables are described here:
	1455	# for the mesoscale:
[763]	1456	specificname_meso = ['UV','uv','uvmet','slopexy','SLOPEXY','deltat','DELTAT','hodograph','tk','hodograph_2']
[754]	1457	# for the gcm:
	1458	specificname_gcm = ['enfact']
	1459
	1460	## Check for variable in file:
	1461	if mode == 'check':
	1462	varname = zvarname
	1463	varinfile=znc.variables.keys()
	1464	logical_novarname = zvarname not in znc.variables
	1465	logical_nospecificname_meso = not ((ztypefile in ['meso']) and (zvarname in specificname_meso))
	1466	logical_nospecificname_gcm = not ((ztypefile in ['gcm']) and (zvarname in specificname_gcm))
	1467	if ( logical_novarname and logical_nospecificname_meso and logical_nospecificname_gcm ):
	1468	if len(varinfile) == 1: varname = varinfile[0]
	1469	else: varname = False
	1470	## Return the variable name:
	1471	return varname
	1472
	1473	## Get the corresponding variable:
	1474	if mode == 'getvar':
[763]	1475	plot_x = None ; plot_y = None ;
[754]	1476	### ----------- 1. saturation temperature
	1477	if zvarname in ["temp","t","T_nadir_nit","T_nadir_day","temp_day","temp_night"] and ztsat:
	1478	tt=getfield(znc,zvarname) ; print "computing Tsat-T, I ASSUME Z-AXIS IS PRESSURE"
	1479	if type(tt).__name__=='MaskedArray': tt.set_fill_value([np.NaN]) ; tinput=tt.filled()
	1480	else: tinput=tt
	1481	all_var=get_tsat(yalt,tinput,zlon=ylon,zlat=ylat,zalt=yalt,ztime=ytime)
	1482	### ----------- 2. wind amplitude
	1483	elif ((zvarname in ['UV','uv','uvmet']) and (ztypefile in ['meso']) and (zvarname not in znc.variables)):
[763]	1484	all_var=windamplitude(znc,'amplitude')
	1485	elif ((zvarname in ['hodograph','hodograph_2']) and (ztypefile in ['meso']) and (zvarname not in znc.variables)):
	1486	plot_x, plot_y = windamplitude(znc,zvarname)
	1487	if plot_x is not None: all_var=plot_x # dummy
	1488	else: all_var=plot_y ; plot_x = None ; plot_y = None # Hodograph type 2 is not 'xy' mode
[754]	1489	elif ((zvarname in ['slopexy','SLOPEXY']) and (ztypefile in ['meso']) and (zvarname not in znc.variables)):
	1490	all_var=slopeamplitude(znc)
	1491	### ------------ 3. Near surface instability
	1492	elif ((zvarname in ['DELTAT','deltat']) and (ztypefile in ['meso']) and (zvarname not in znc.variables)):
	1493	all_var=deltat0t1(znc)
	1494	### ------------ 4. Enrichment factor
	1495	elif ((ztypefile in ['gcm']) and (zvarname in ['enfact'])):
	1496	all_var=enrichment_factor(znc,ylon,ylat,ytime)
[763]	1497	### ------------ 5. teta -> temp
	1498	elif ((ztypefile in ['meso']) and (zvarname in ['tk']) and ('tk' not in znc.variables.keys())):
	1499	all_var=teta_to_tk(znc)
[754]	1500	else:
	1501	### ----------- 999. Normal case
	1502	all_var = getfield(znc,zvarname)
[763]	1503	if analysis is not None:
	1504	if analysis in ['histo','density','histodensity']: plot_y=all_var ; plot_x = plot_y
	1505	elif analysis == 'fft': plot_y, plot_x = spectrum(all_var,ytime,yalt,ylat,ylon) ; all_var = plot_y
	1506	return all_var, plot_x, plot_y
	1507
	1508	# Author : A.C
	1509	# FFT is computed before reducefield voluntarily, because we dont want to compute
	1510	# ffts on averaged fields (which would kill all waves). Instead, we take the fft everywhere
	1511	# (which is not efficient but it is still ok) and then, make the average (if the user wants to)
	1512	def spectrum(var,time,vert,lat,lon):
	1513	import numpy as np
	1514	fft=np.fft.fft(var,axis=1)
	1515	N=len(vert)
	1516	step=(vert[1]-vert[0])*1000.
	1517	print "step is: ",step
	1518	fftfreq=np.fft.fftfreq(N,d=step)
	1519	fftfreq=np.fft.fftshift(fftfreq) # spatial FFT => this is the wavenumber
	1520	fft=np.fft.fftshift(fft)
	1521	fftfreq = 1./fftfreq # => wavelength (div by 0 expected, don't panic)
	1522	fft=np.abs(fft) # => amplitude spectrum
	1523	# fft=np.abs(fft)**2 # => power spectrum
	1524	return fft,fftfreq
	1525
	1526	# Author : A.C.
	1527	# Computes temperature from potential temperature for mesoscale files, without the need to use API, i.e. using natural vertical grid
	1528	def teta_to_tk(nc):
	1529	import numpy as np
	1530	varinfile = nc.variables.keys()
	1531	p0=610.
	1532	t0=220.
	1533	r_cp=1./3.89419
	1534	if "T" in varinfile: zteta=getfield(nc,'T')
	1535	else: errormess("you need T in your file.")
	1536	if "PTOT" in varinfile: zptot=getfield(nc,'PTOT')
	1537	else: errormess("you need PTOT in your file.")
	1538	zt=(zteta+220.)(zptot/p0)*(r_cp)
	1539	return zt

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

Download in other formats: