Context Navigation

source: trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/myplot.py @ 251

Last change on this file since 251 was 247, checked in by aslmd, 14 years ago
MESOSCALE: minor bug fixes in execution script and graphical python routines.
File size: 24.9 KB

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

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

Download in other formats: