Context Navigation

source: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py @ 294

Last change on this file since 294 was 287, checked in by aslmd, 14 years ago
MESOSCALE: minor updates to parameter files, plotting python module, launching and post-processing
File size: 25.3 KB

Rev	Line
[252]	1	def errormess(text,printvar=None):
[233]	2	print text
[252]	3	if printvar: print printvar
[233]	4	exit()
	5	return
	6
[252]	7	def getname(var=False,winds=False,anomaly=False):
	8	if var and winds: basename = var + '_UV'
	9	elif var: basename = var
	10	elif winds: basename = 'UV'
	11	else: errormess("please set at least winds or var",printvar=nc.variables)
	12	if anomaly: basename = 'd' + basename
	13	return basename
	14
	15	def localtime(utc,lon):
	16	ltst = utc + lon / 15.
	17	ltst = int (ltst * 10) / 10.
	18	ltst = ltst % 24
	19	return ltst
	20
[233]	21	def whatkindfile (nc):
	22	if 'controle' in nc.variables: typefile = 'gcm'
[254]	23	elif 'phisinit' in nc.variables: typefile = 'gcmex'
[233]	24	elif 'vert' in nc.variables: typefile = 'mesoapi'
	25	elif 'U' in nc.variables: typefile = 'meso'
	26	elif 'HGT_M' in nc.variables: typefile = 'geo'
	27	else: errormess("whatkindfile: typefile not supported.")
	28	return typefile
	29
	30	def getfield (nc,var):
	31	## this allows to get much faster (than simply referring to nc.variables[var])
	32	dimension = len(nc.variables[var].dimensions)
	33	if dimension == 2: field = nc.variables[var][:,:]
	34	elif dimension == 3: field = nc.variables[var][:,:,:]
	35	elif dimension == 4: field = nc.variables[var][:,:,:,:]
	36	return field
	37
	38	def reducefield (input,d4=None,d3=None,d2=None,d1=None):
[252]	39	### we do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x"
[233]	40	### it would be actually better to name d4 d3 d2 d1 as t z y x
	41	import numpy as np
	42	dimension = np.array(input).ndim
	43	shape = np.array(input).shape
	44	print 'dim,shape: ',dimension,shape
	45	output = input
	46	error = False
	47	if dimension == 2:
	48	if d2 >= shape[0]: error = True
	49	elif d1 >= shape[1]: error = True
	50	elif d1 is not None and d2 is not None: output = input[d2,d1]
	51	elif d1 is not None: output = input[:,d1]
	52	elif d2 is not None: output = input[d2,:]
	53	elif dimension == 3:
	54	if d4 >= shape[0]: error = True
	55	elif d2 >= shape[1]: error = True
	56	elif d1 >= shape[2]: error = True
	57	elif d4 is not None and d2 is not None and d1 is not None: output = input[d4,d2,d1]
	58	elif d4 is not None and d2 is not None: output = input[d4,d2,:]
	59	elif d4 is not None and d1 is not None: output = input[d4,:,d1]
	60	elif d2 is not None and d1 is not None: output = input[:,d2,d1]
	61	elif d1 is not None: output = input[:,:,d1]
	62	elif d2 is not None: output = input[:,d2,:]
	63	elif d4 is not None: output = input[d4,:,:]
	64	elif dimension == 4:
	65	if d4 >= shape[0]: error = True
	66	elif d3 >= shape[1]: error = True
	67	elif d2 >= shape[2]: error = True
	68	elif d1 >= shape[3]: error = True
	69	elif d4 is not None and d3 is not None and d2 is not None and d1 is not None: output = input[d4,d3,d2,d1]
	70	elif d4 is not None and d3 is not None and d2 is not None: output = input[d4,d3,d2,:]
	71	elif d4 is not None and d3 is not None and d1 is not None: output = input[d4,d3,:,d1]
	72	elif d4 is not None and d2 is not None and d1 is not None: output = input[d4,:,d2,d1]
	73	elif d3 is not None and d2 is not None and d1 is not None: output = input[:,d3,d2,d1]
	74	elif d4 is not None and d3 is not None: output = input[d4,d3,:,:]
	75	elif d4 is not None and d2 is not None: output = input[d4,:,d2,:]
	76	elif d4 is not None and d1 is not None: output = input[d4,:,:,d1]
	77	elif d3 is not None and d2 is not None: output = input[:,d3,d2,:]
	78	elif d3 is not None and d1 is not None: output = input[:,d3,:,d1]
	79	elif d2 is not None and d1 is not None: output = input[:,:,d2,d1]
	80	elif d1 is not None: output = input[:,:,:,d1]
	81	elif d2 is not None: output = input[:,:,d2,:]
	82	elif d3 is not None: output = input[:,d3,:,:]
	83	elif d4 is not None: output = input[d4,:,:,:]
	84	dimension = np.array(output).ndim
	85	shape = np.array(output).shape
	86	print 'dim,shape: ',dimension,shape
	87	return output, error
	88
	89	def definesubplot ( numplot, fig ):
	90	from matplotlib.pyplot import rcParams
	91	rcParams['font.size'] = 12. ## default (important for multiple calls)
	92	if numplot == 4:
	93	sub = 221
	94	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
	95	rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
	96	elif numplot == 2:
	97	sub = 121
	98	fig.subplots_adjust(wspace = 0.35)
	99	rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
	100	elif numplot == 3:
	101	sub = 131
	102	fig.subplots_adjust(wspace = 0.5)
	103	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
	104	elif numplot == 6:
	105	sub = 231
	106	fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
	107	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
	108	elif numplot == 8:
	109	sub = 331 #241
	110	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
	111	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
	112	elif numplot == 9:
	113	sub = 331
	114	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
	115	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
	116	elif numplot == 1:
	117	sub = 99999
[241]	118	elif numplot <= 0:
[233]	119	sub = 99999
	120	else:
	121	print "supported: 1,2,3,4,6,8,9"
	122	exit()
	123	return sub
	124
	125	def getstralt(nc,nvert):
	126	typefile = whatkindfile(nc)
	127	if typefile is 'meso':
	128	stralt = "_lvl" + str(nvert)
	129	elif typefile is 'mesoapi':
	130	zelevel = int(nc.variables['vert'][nvert])
	131	if abs(zelevel) < 10000.: strheight=str(zelevel)+"m"
	132	else: strheight=str(int(zelevel/1000.))+"km"
	133	if 'altitude' in nc.dimensions: stralt = "_"+strheight+"-AMR"
	134	elif 'altitude_abg' in nc.dimensions: stralt = "_"+strheight+"-ALS"
	135	elif 'bottom_top' in nc.dimensions: stralt = "_"+strheight
	136	elif 'pressure' in nc.dimensions: stralt = "_"+str(zelevel)+"Pa"
	137	else: stralt = ""
	138	else:
	139	stralt = ""
	140	return stralt
	141
[195]	142	def getlschar ( namefile ):
	143	from netCDF4 import Dataset
	144	from timestuff import sol2ls
[233]	145	from numpy import array
[195]	146	nc = Dataset(namefile)
[237]	147	zetime = None
[233]	148	if 'Times' in nc.variables:
	149	zetime = nc.variables['Times'][0]
	150	shape = array(nc.variables['Times']).shape
	151	if shape[0] < 2: zetime = None
	152	if zetime is not None \
[225]	153	and 'vert' not in nc.variables:
[233]	154	#### strangely enough this does not work for api or ncrcat results!
[195]	155	zetimestart = getattr(nc, 'START_DATE')
	156	zeday = int(zetime[8]+zetime[9]) - int(zetimestart[8]+zetimestart[9])
	157	if zeday < 0: lschar="" ## might have crossed a month... fix soon
[241]	158	else: lschar="_Ls"+str( int( 10. * sol2ls ( getattr( nc, 'JULDAY' ) + zeday ) ) / 10. )
[197]	159	###
	160	zetime2 = nc.variables['Times'][1]
	161	one = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37.
	162	next = int(zetime2[11]+zetime2[12]) + int(zetime2[14]+zetime2[15])/37.
	163	zehour = one
	164	zehourin = abs ( next - one )
[195]	165	else:
	166	lschar=""
[197]	167	zehour = 0
	168	zehourin = 1
	169	return lschar, zehour, zehourin
[195]	170
[202]	171	def getprefix (nc):
	172	prefix = 'LMD_MMM_'
	173	prefix = prefix + 'd'+str(getattr(nc,'GRID_ID'))+'_'
	174	prefix = prefix + str(int(getattr(nc,'DX')/1000.))+'km_'
	175	return prefix
	176
[184]	177	def getproj (nc):
[233]	178	typefile = whatkindfile(nc)
	179	if typefile in ['mesoapi','meso','geo']:
	180	### (il faudrait passer CEN_LON dans la projection ?)
	181	map_proj = getattr(nc, 'MAP_PROJ')
	182	cen_lat = getattr(nc, 'CEN_LAT')
	183	if map_proj == 2:
	184	if cen_lat > 10.:
	185	proj="npstere"
	186	print "NP stereographic polar domain"
	187	else:
	188	proj="spstere"
	189	print "SP stereographic polar domain"
	190	elif map_proj == 1:
	191	print "lambert projection domain"
	192	proj="lcc"
	193	elif map_proj == 3:
	194	print "mercator projection"
	195	proj="merc"
	196	else:
	197	proj="merc"
[252]	198	elif typefile in ['gcm']: proj="cyl" ## pb avec les autres (de trace derriere la sphere ?)
[233]	199	else: proj="ortho"
[184]	200	return proj
	201
[180]	202	def ptitle (name):
	203	from matplotlib.pyplot import title
	204	title(name)
	205	print name
	206
[252]	207	def polarinterv (lon2d,lat2d):
	208	import numpy as np
	209	wlon = [np.min(lon2d),np.max(lon2d)]
	210	ind = np.array(lat2d).shape[0] / 2 ## to get a good boundlat and to get the pole
	211	wlat = [np.min(lat2d[ind,:]),np.max(lat2d[ind,:])]
	212	return [wlon,wlat]
	213
[180]	214	def simplinterv (lon2d,lat2d):
	215	import numpy as np
	216	return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]]
	217
[184]	218	def wrfinterv (lon2d,lat2d):
	219	nx = len(lon2d[0,:])-1
	220	ny = len(lon2d[:,0])-1
[225]	221	lon1 = lon2d[0,0]
	222	lon2 = lon2d[nx,ny]
	223	lat1 = lat2d[0,0]
	224	lat2 = lat2d[nx,ny]
[233]	225	if abs(0.5(lat1+lat2)) > 60.: wider = 0.5 (abs(lon1)+abs(lon2)) * 0.1
	226	else: wider = 0.
	227	if lon1 < lon2: wlon = [lon1, lon2 + wider]
[225]	228	else: wlon = [lon2, lon1 + wider]
	229	if lat1 < lat2: wlat = [lat1, lat2]
	230	else: wlat = [lat2, lat1]
	231	return [wlon,wlat]
[184]	232
[240]	233	def makeplotres (filename,res=None,pad_inches_value=0.25,folder='',disp=True,ext='png',erase=False):
[180]	234	import matplotlib.pyplot as plt
[240]	235	from os import system
	236	addstr = ""
	237	if res is not None:
	238	res = int(res)
	239	addstr = "_"+str(res)
	240	name = filename+addstr+"."+ext
[186]	241	if folder != '': name = folder+'/'+name
[180]	242	plt.savefig(name,dpi=res,bbox_inches='tight',pad_inches=pad_inches_value)
[240]	243	if disp: display(name)
	244	if ext in ['eps','ps','svg']: system("tar czvf "+name+".tar.gz "+name+" ; rm -f "+name)
	245	if erase: system("mv "+name+" to_be_erased")
[180]	246	return
	247
[240]	248	def dumpbdy (field,n,stag=None):
[184]	249	nx = len(field[0,:])-1
	250	ny = len(field[:,0])-1
[233]	251	if stag == 'U': nx = nx-1
	252	if stag == 'V': ny = ny-1
[240]	253	return field[n:ny-n,n:nx-n]
[180]	254
[233]	255	def getcoorddef ( nc ):
	256	## getcoord2d for predefined types
	257	typefile = whatkindfile(nc)
	258	if typefile in ['mesoapi','meso']:
	259	[lon2d,lat2d] = getcoord2d(nc)
[240]	260	lon2d = dumpbdy(lon2d,6)
	261	lat2d = dumpbdy(lat2d,6)
[254]	262	elif typefile in ['gcm','gcmex']:
[233]	263	[lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True)
	264	elif typefile in ['geo']:
	265	[lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M')
	266	return lon2d,lat2d
	267
[184]	268	def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False):
	269	import numpy as np
	270	if is1d:
	271	lat = nc.variables[nlat][:]
	272	lon = nc.variables[nlon][:]
	273	[lon2d,lat2d] = np.meshgrid(lon,lat)
	274	else:
	275	lat = nc.variables[nlat][0,:,:]
	276	lon = nc.variables[nlon][0,:,:]
	277	[lon2d,lat2d] = [lon,lat]
	278	return lon2d,lat2d
	279
[180]	280	def smooth (field, coeff):
	281	## actually blur_image could work with different coeff on x and y
	282	if coeff > 1: result = blur_image(field,int(coeff))
	283	else: result = field
	284	return result
	285
	286	def gauss_kern(size, sizey=None):
	287	import numpy as np
	288	## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
	289	# Returns a normalized 2D gauss kernel array for convolutions
	290	size = int(size)
	291	if not sizey:
	292	sizey = size
	293	else:
	294	sizey = int(sizey)
	295	x, y = np.mgrid[-size:size+1, -sizey:sizey+1]
	296	g = np.exp(-(x2/float(size)+y2/float(sizey)))
	297	return g / g.sum()
	298
	299	def blur_image(im, n, ny=None) :
	300	from scipy.signal import convolve
	301	## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
	302	# blurs the image by convolving with a gaussian kernel of typical size n.
	303	# The optional keyword argument ny allows for a different size in the y direction.
	304	g = gauss_kern(n, sizey=ny)
	305	improc = convolve(im, g, mode='same')
	306	return improc
	307
[233]	308	def getwinddef (nc):
	309	## getwinds for predefined types
	310	typefile = whatkindfile(nc)
	311	###
	312	if typefile is 'mesoapi': [uchar,vchar] = ['Um','Vm']
	313	elif typefile is 'gcm': [uchar,vchar] = ['u','v']
	314	elif typefile is 'meso': [uchar,vchar] = ['U','V']
	315	else: [uchar,vchar] = ['not found','not found']
	316	###
	317	if typefile in ['meso']: metwind = False ## geometrical (wrt grid)
	318	else: metwind = True ## meteorological (zon/mer)
	319	if metwind is False: print "Not using meteorological winds. You trust numerical grid as being (x,y)"
	320	###
	321	return uchar,vchar,metwind
[202]	322
[184]	323	def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True):
	324	## scale regle la reference du vecteur
	325	## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs.
	326	import matplotlib.pyplot as plt
	327	import numpy as np
[187]	328	posx = np.min(x) - np.std(x) / 10.
	329	posy = np.min(y) - np.std(y) / 10.
[184]	330	u = smooth(u,csmooth)
	331	v = smooth(v,csmooth)
[188]	332	widthvec = 0.003 #0.005 #0.003
[184]	333	q = plt.quiver( x[::stride,::stride],\
	334	y[::stride,::stride],\
	335	u[::stride,::stride],\
	336	v[::stride,::stride],\
[228]	337	angles='xy',color=color,pivot='middle',\
[184]	338	scale=factor,width=widthvec )
[202]	339	if color in ['white','yellow']: kcolor='black'
	340	else: kcolor=color
[184]	341	if key: p = plt.quiverkey(q,posx,posy,scale,\
[194]	342	str(int(scale)),coordinates='data',color=kcolor,labelpos='S',labelsep = 0.03)
[184]	343	return
[180]	344
	345	def display (name):
[184]	346	from os import system
	347	system("display "+name+" > /dev/null 2> /dev/null &")
	348	return name
[180]	349
	350	def findstep (wlon):
[184]	351	steplon = int((wlon[1]-wlon[0])/4.) #3
	352	step = 120.
	353	while step > steplon and step > 15. : step = step / 2.
	354	if step <= 15.:
	355	while step > steplon and step > 5. : step = step - 5.
	356	if step <= 5.:
	357	while step > steplon and step > 1. : step = step - 1.
	358	if step <= 1.:
	359	step = 1.
[180]	360	return step
	361
[233]	362	def define_proj (char,wlon,wlat,back=None):
[180]	363	from mpl_toolkits.basemap import Basemap
	364	import numpy as np
	365	import matplotlib as mpl
[240]	366	from mymath import max
[180]	367	meanlon = 0.5*(wlon[0]+wlon[1])
	368	meanlat = 0.5*(wlat[0]+wlat[1])
[184]	369	if wlat[0] >= 80.: blat = 40.
[247]	370	elif wlat[1] <= -80.: blat = -40.
	371	elif wlat[1] >= 0.: blat = wlat[0]
	372	elif wlat[0] <= 0.: blat = wlat[1]
[252]	373	print "blat ", blat
[207]	374	h = 50. ## en km
[202]	375	radius = 3397200.
[184]	376	if char == "cyl": m = Basemap(rsphere=radius,projection='cyl',\
[180]	377	llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
[184]	378	elif char == "moll": m = Basemap(rsphere=radius,projection='moll',lon_0=meanlon)
	379	elif char == "ortho": m = Basemap(rsphere=radius,projection='ortho',lon_0=meanlon,lat_0=meanlat)
	380	elif char == "lcc": m = Basemap(rsphere=radius,projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\
	381	llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
	382	elif char == "npstere": m = Basemap(rsphere=radius,projection='npstere', boundinglat=blat, lon_0=0.)
	383	elif char == "spstere": m = Basemap(rsphere=radius,projection='spstere', boundinglat=blat, lon_0=0.)
[207]	384	elif char == "nplaea": m = Basemap(rsphere=radius,projection='nplaea', boundinglat=wlat[0], lon_0=meanlon)
	385	elif char == "laea": m = Basemap(rsphere=radius,projection='laea',lon_0=meanlon,lat_0=meanlat,lat_ts=meanlat,\
	386	llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
[184]	387	elif char == "nsper": m = Basemap(rsphere=radius,projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=h*1000.)
	388	elif char == "merc": m = Basemap(rsphere=radius,projection='merc',lat_ts=0.,\
	389	llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
	390	fontsizemer = int(mpl.rcParams['font.size']*3./4.)
[207]	391	if char in ["cyl","lcc","merc","nsper","laea"]: step = findstep(wlon)
	392	else: step = 10.
[238]	393	steplon = step*2.
	394	#if back in ["geolocal"]:
	395	# step = np.min([5.,step])
	396	# steplon = step
[225]	397	print step
[238]	398	m.drawmeridians(np.r_[-180.:180.:steplon], labels=[0,0,0,1], color='grey', fontsize=fontsizemer)
[180]	399	m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color='grey', fontsize=fontsizemer)
[233]	400	if back: m.warpimage(marsmap(back),scale=0.75)
	401	#if not back:
	402	# if not var: back = "mola" ## if no var: draw mola
	403	# elif typefile in ['mesoapi','meso','geo'] \
	404	# and proj not in ['merc','lcc','nsper','laea']: back = "molabw" ## if var but meso: draw molabw
	405	# else: pass ## else: draw None
[180]	406	return m
	407
[232]	408	#### test temporaire
	409	def putpoints (map,plot):
	410	#### from http://www.scipy.org/Cookbook/Matplotlib/Maps
	411	# lat/lon coordinates of five cities.
	412	lats = [18.4]
	413	lons = [-134.0]
	414	points=['Olympus Mons']
	415	# compute the native map projection coordinates for cities.
	416	x,y = map(lons,lats)
	417	# plot filled circles at the locations of the cities.
	418	map.plot(x,y,'bo')
	419	# plot the names of those five cities.
	420	wherept = 0 #1000 #50000
	421	for name,xpt,ypt in zip(points,x,y):
	422	plot.text(xpt+wherept,ypt+wherept,name)
	423	## le nom ne s'affiche pas...
	424	return
	425
[233]	426	def calculate_bounds(field,vmin=None,vmax=None):
	427	import numpy as np
	428	from mymath import max,min,mean
	429	ind = np.where(field < 9e+35)
	430	fieldcalc = field[ ind ] # la syntaxe compacte ne marche si field est un tuple
	431	###
	432	dev = np.std(fieldcalc)*3.0
	433	###
	434	if vmin is None:
	435	zevmin = mean(fieldcalc) - dev
	436	else: zevmin = vmin
	437	###
	438	if vmax is None: zevmax = mean(fieldcalc) + dev
	439	else: zevmax = vmax
	440	if vmin == vmax:
	441	zevmin = mean(fieldcalc) - dev ### for continuity
	442	zevmax = mean(fieldcalc) + dev ### for continuity
	443	###
	444	if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0.
	445	print "field ", min(fieldcalc), max(fieldcalc)
	446	print "bounds ", zevmin, zevmax
	447	return zevmin, zevmax
[232]	448
[233]	449	def bounds(what_I_plot,zevmin,zevmax):
[247]	450	from mymath import max,min,mean
[233]	451	### might be convenient to add the missing value in arguments
	452	what_I_plot[ what_I_plot < zevmin ] = zevmin*(1. + 1.e-7)
[247]	453	print "new min ", min(what_I_plot)
[233]	454	what_I_plot[ what_I_plot > 9e+35 ] = -9e+35
	455	what_I_plot[ what_I_plot > zevmax ] = zevmax*(1. - 1.e-7)
[247]	456	print "new max ", max(what_I_plot)
[233]	457	return what_I_plot
	458
[241]	459	def nolow(what_I_plot):
	460	from mymath import max,min
	461	lim = 0.150.5(abs(max(what_I_plot))+abs(min(what_I_plot)))
	462	print "on vire en dessous de ", lim
	463	what_I_plot [ abs(what_I_plot) < lim ] = 1.e40
	464	return what_I_plot
	465
[233]	466	def zoomset (wlon,wlat,zoom):
	467	dlon = abs(wlon[1]-wlon[0])/2.
	468	dlat = abs(wlat[1]-wlat[0])/2.
	469	[wlon,wlat] = [ [wlon[0]+zoomdlon/100.,wlon[1]-zoomdlon/100.],\
	470	[wlat[0]+zoomdlat/100.,wlat[1]-zoomdlat/100.] ]
	471	print "zoom %",zoom,wlon,wlat
	472	return wlon,wlat
	473
[201]	474	def fmtvar (whichvar="def"):
[204]	475	fmtvar = { \
	476	"tk": "%.0f",\
	477	"tpot": "%.0f",\
	478	"def": "%.1e",\
	479	"PTOT": "%.0f",\
	480	"HGT": "%.1e",\
	481	"USTM": "%.2f",\
[225]	482	"HFX": "%.0f",\
[232]	483	"ICETOT": "%.1e",\
[237]	484	"TAU_ICE": "%.2f",\
[252]	485	"VMR_ICE": "%.1e",\
	486	"MTOT": "%.0f",\
[240]	487	"anomaly": "%.1f",\
[241]	488	"W": "%.1f",\
[287]	489	"WMAX_TH": "%.1f",\
	490	"QSURFICE": "%.0f",\
[204]	491	}
	492	if whichvar not in fmtvar:
	493	whichvar = "def"
	494	return fmtvar[whichvar]
[201]	495
[233]	496	####################################################################################################################
	497	### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
[202]	498	def defcolorb (whichone="def"):
[204]	499	whichcolorb = { \
	500	"def": "spectral",\
	501	"HGT": "spectral",\
	502	"tk": "gist_heat",\
	503	"QH2O": "PuBu",\
	504	"USTM": "YlOrRd",\
[225]	505	"HFX": "RdYlBu",\
[232]	506	"ICETOT": "YlGnBu",\
[252]	507	"MTOT": "PuBu",\
[238]	508	"TAU_ICE": "Blues",\
[252]	509	"VMR_ICE": "Blues",\
[241]	510	"W": "jet",\
[287]	511	"WMAX_TH": "spectral",\
[240]	512	"anomaly": "RdBu_r",\
[287]	513	"QSURFICE": "hot_r",\
[204]	514	}
[241]	515	#W --> spectral ou jet
[240]	516	#spectral BrBG RdBu_r
[241]	517	print "predefined colorbars"
[204]	518	if whichone not in whichcolorb:
	519	whichone = "def"
	520	return whichcolorb[whichone]
[202]	521
	522	def definecolorvec (whichone="def"):
	523	whichcolor = { \
	524	"def": "black",\
	525	"vis": "yellow",\
	526	"vishires": "yellow",\
	527	"molabw": "yellow",\
	528	"mola": "black",\
	529	"gist_heat": "white",\
	530	"hot": "tk",\
	531	"gist_rainbow": "black",\
	532	"spectral": "black",\
	533	"gray": "red",\
	534	"PuBu": "black",\
	535	}
	536	if whichone not in whichcolor:
	537	whichone = "def"
	538	return whichcolor[whichone]
	539
[180]	540	def marsmap (whichone="vishires"):
[233]	541	from os import uname
	542	mymachine = uname()[1]
	543	### not sure about speed-up with this method... looks the same
	544	if "lmd.jussieu.fr" in mymachine: domain = "/u/aslmd/WWW/maps/"
	545	else: domain = "http://www.lmd.jussieu.fr/~aslmd/maps/"
[180]	546	whichlink = { \
[233]	547	#"vis": "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\
	548	#"vishires": "http://www.lmd.jussieu.fr/~aslmd/maps/MarsMap_2500x1250.jpg",\
	549	#"geolocal": "http://dl.dropbox.com/u/11078310/geolocal.jpg",\
	550	#"mola": "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\
	551	#"molabw": "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\
	552	"vis": domain+"mar0kuu2.jpg",\
	553	"vishires": domain+"MarsMap_2500x1250.jpg",\
	554	"geolocal": domain+"geolocal.jpg",\
	555	"mola": domain+"mars-mola-2k.jpg",\
	556	"molabw": domain+"MarsElevation_2500x1250.jpg",\
[238]	557	"clouds": "http://www.johnstonsarchive.net/spaceart/marswcloudmap.jpg",\
	558	"jupiter": "http://www.mmedia.is/~bjj/data/jupiter_css/jupiter_css.jpg",\
	559	"jupiter_voy": "http://www.mmedia.is/~bjj/data/jupiter/jupiter_vgr2.jpg",\
[273]	560	"bw": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg",\
	561	"contrast": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg",\
	562	"nice": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg",\
	563	"blue": "http://eoimages.gsfc.nasa.gov/ve/2430/land_ocean_ice_2048.jpg",\
[180]	564	}
[238]	565	### see http://www.mmedia.is/~bjj/planetary_maps.html
[180]	566	if whichone not in whichlink:
	567	print "marsmap: choice not defined... you'll get the default one... "
	568	whichone = "vishires"
	569	return whichlink[whichone]
	570
[273]	571	#def earthmap (whichone):
	572	# if whichone == "contrast": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg"
	573	# elif whichone == "bw": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg"
	574	# elif whichone == "nice": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg"
	575	# return whichlink
[180]	576
[241]	577	def latinterv (area="Whole"):
	578	list = { \
	579	"Europe": [[ 20., 80.],[- 50., 50.]],\
	580	"Central_America": [[-10., 40.],[ 230., 300.]],\
	581	"Africa": [[-20., 50.],[- 50., 50.]],\
[273]	582	"Whole": [[-90., 90.],[-180., 180.]],\
	583	"Southern_Hemisphere": [[-90., 60.],[-180., 180.]],\
	584	"Northern_Hemisphere": [[-60., 90.],[-180., 180.]],\
[241]	585	"Tharsis": [[-30., 60.],[-170.,- 10.]],\
	586	"Whole_No_High": [[-60., 60.],[-180., 180.]],\
	587	"Chryse": [[-60., 60.],[- 60., 60.]],\
	588	"North_Pole": [[ 50., 90.],[-180., 180.]],\
	589	"Close_North_Pole": [[ 75., 90.],[-180., 180.]],\
	590	"Far_South_Pole": [[-90.,-40.],[-180., 180.]],\
	591	"South_Pole": [[-90.,-50.],[-180., 180.]],\
	592	"Close_South_Pole": [[-90.,-75.],[-180., 180.]],\
	593	}
	594	if area not in list: area = "Whole"
	595	[olat,olon] = list[area]
	596	return olon,olat
	597

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

Download in other formats: