Context Navigation

source: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/planetoplot.py @ 345

Last change on this file since 345 was 345, checked in by aslmd, 14 years ago
PYTHON: updates from TN. planetoplot.py is a version of winds.py which works for GCM files and adds section, 1Dplot capabilities. new functions in myplot. checked compatibility with winds.py, OK.
Property svn:executable set to ``*
File size: 24.8 KB

Line
1	#!/usr/bin/env python
2
3	### A. Spiga -- LMD -- 30/06/2011 to 10/07/2011
4	### T.Navarro -- LMD -- 10/2011
5	### Thanks to A. Colaitis for the parser trick
6
7
8	####################################
9	####################################
10	### The main program to plot vectors
11	def planetoplot (namefiles,\
12	vertmode,\
13	proj=None,\
14	back=None,\
15	target=None,
16	stride=3,\
17	numplot=2,\
18	var=None,\
19	colorb="def",\
20	winds=True,\
21	addchar=None,\
22	interv=[0,1],\
23	vmin=None,\
24	vmax=None,\
25	tile=False,\
26	zoom=None,\
27	display=True,\
28	itstep=None,\
29	hole=False,\
30	save="gui",\
31	anomaly=False,\
32	var2=None,\
33	ndiv=10,\
34	first=1,\
35	mult=1.,\
36	zetitle="fill",\
37	slon=None,\
38	slat=None,\
39	svert=None,\
40	stime=None):
41
42	####################################################################################################################
43	### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
44
45	#################################
46	### Load librairies and functions
47	from netCDF4 import Dataset
48	from myplot import getcoord2d,define_proj,makeplotres,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy,\
49	fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\
50	zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow,\
51	getname,localtime,polarinterv,getsindex,define_axis
52	from mymath import deg,max,min,mean
53	from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel
54	from matplotlib.cm import get_cmap
55	import numpy as np
56	from numpy.core.defchararray import find
57
58	#whichmode = ""
59	nlon = 1 # number of longitudinal slices -- 1 is None
60	nlat = 1
61	nvert = 1
62	ntime = 1
63	nslices = 1
64	if slon is not None:
65	nslices = nslices*len(slon)
66	nlon = len(slon)
67	if slat is not None:
68	nslices = nslices*len(slat)
69	nlat = len(slat)
70	if svert is not None:
71	nslices = nslices*len(svert)
72	nvert = len(svert)
73	if stime is not None:
74	nslices = nslices*len(stime)
75	ntime = len(stime)
76	#else:
77	# nslices = 2
78	numplot = len(namefiles)*nslices
79
80	mapmode = 0
81	if slon is None and slat is None:
82	mapmode = 1 # in this case we plot a map, with the given projection
83	elif proj is not None:
84	print "WARNING: you specified a", proj,\
85	"projection but asked for slices", slon,"in longitude and",slat,"in latitude"
86	print "mapmode", mapmode
87
88
89	all_var = [[]]*len(namefiles)
90	all_var2 = [[]]*len(namefiles)
91	all_title = [[]]*len(namefiles)
92
93	print "len(namefiles), nslices", len(namefiles), nslices
94	print "numplot", numplot
95
96	#########################
97	### Loop over the files initially separated by comas to be plotted on the same figure
98	k = 0
99	firstfile = True
100	for namefile in namefiles:
101
102	######################
103	### Load NETCDF object
104	nc = Dataset(namefile)
105
106	##################################
107	### Initial checks and definitions
108	typefile = whatkindfile(nc) ## TYPEFILE
109	if firstfile:
110	typefile0 = typefile
111	elif typefile != typefile0:
112	errormess("Not the same files !", [typefile0, typefile])
113	if var not in nc.variables: var = False ## VAR
114	if winds: ## WINDS
115	[uchar,vchar,metwind] = getwinddef(nc)
116	if uchar == 'not found': winds = False
117	if not var and not winds: errormess("please set at least winds or var",printvar=nc.variables)
118	[lon2d,lat2d] = getcoorddef(nc) ## COORDINATES, could be moved below
119	if proj == None: proj = getproj(nc) ## PROJECTION
120
121	lat = nc.variables["latitude"][:]
122	lon = nc.variables["longitude"][:]
123	time = nc.variables["Time"][:]
124	vert = nc.variables["altitude"][:]
125	if firstfile:
126	lat0 = lat
127	elif len(lat0) != len(lat):
128	errormess("Not the same latitude lengths !", [len(lat0), len(lat)])
129	elif sum((lat == lat0) == False) != 0:
130	errormess("Not the same latitudes !", [lat,lat0])
131	# Faire d'autre checks sur les compatibilites entre fichiers!!
132
133	if firstfile:
134	##########################
135	### Define plot boundaries
136	### todo: possible areas in latinterv in argument (ex: "Far_South_Pole")
137	if proj in ["npstere","spstere"]: [wlon,wlat] = polarinterv(lon2d,lat2d)
138	elif proj in ["lcc","laea"]: [wlon,wlat] = wrfinterv(lon2d,lat2d)
139	else: [wlon,wlat] = simplinterv(lon2d,lat2d)
140	if zoom: [wlon,wlat] = zoomset(wlon,wlat,zoom)
141
142	#########################################
143	### Name for title and graphics save file
144	basename = getname(var=var,winds=winds,anomaly=anomaly)
145	basename = basename + getstralt(nc,nvert) ## can be moved elsewhere for a more generic routine
146
147	print "var", var
148	print "var2", var2
149	#print all_var
150	#print nc
151	if var: all_var[k] = getfield(nc,var)
152	if var2: all_var2[k] = getfield(nc,var2)
153
154	print "k", k
155	print "all_var[k].shape", all_var[k].shape
156	k += 1
157	firstfile = False
158	#### End of for namefile in namefiles
159
160
161	##################################
162	### Open a figure and set subplots
163	fig = figure()
164	subv,subh = definesubplot( numplot, fig )
165
166	#################################
167	### Time loop for plotting device
168	found_lct = False
169	nplot = 1
170	itime = first
171	error = False
172	if itstep is None and numplot > 0: itstep = int(24./numplot)
173	elif numplot <= 0: itstep = 1
174
175	### Map projection
176	if mapmode == 1:
177	m = define_proj(proj,wlon,wlat,back=back)
178	x, y = m(lon2d, lat2d)
179
180	#for nplot in range(numplot):
181	while error is False:
182	print "nplot", nplot
183	print error
184
185	### Which local time ?
186	ltst = localtime ( interv[0]+itimeinterv[1], 0.5(wlon[0]+wlon[1]) )
187
188	### General plot settings
189	#print itime, int(ltst), numplot, nplot
190	if numplot >= 1:
191	if nplot > numplot: break
192	if numplot > 1:
193	if typefile not in ['geo']: subplot(subv,subh,nplot)
194	found_lct = True
195	### If only one local time is requested (numplot < 0)
196	elif numplot <= 0:
197	if int(ltst) + numplot != 0:
198	itime += 1
199	if found_lct is True: break ## because it means LT was found at previous iteration
200	else: continue ## continue to iterate to find the correct LT
201	else:
202	found_lct = True
203
204
205	## get all indexes to be taken into account for this subplot and then reduce field
206	## We plot 1) all lon slices 2) all lat slices 3) all vert slices 4) all time slices and then go to the next slice
207	#print "(nplot-1)%nlon", (nplot-1)%nlon
208	indexlon = getsindex(slon,(nplot-1)%nlon,lon)
209	#print "indexlon", indexlon
210	indexlat = getsindex(slat,((nplot-1)//nlon)%nlat,lat)
211	indexvert = getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert)
212	indextime = getsindex(stime,((nplot-1)//(nlonnlatnvert))%ntime,time)
213	index_f = ((nplot-1)//(nlonnlatnvert*ntime))%len(namefiles)
214	#print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime)
215
216	#### Contour plot
217	if var2:
218	what_I_plot, error = reducefield(all_var2[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert )
219	what_I_plot = what_I_plot*mult
220	if not error:
221	if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(what_I_contour,6)
222	zevmin, zevmax = calculate_bounds(what_I_plot)
223	zelevels = np.linspace(zevmin,zevmax,num=ndiv)
224	if mapmode == 0:
225	x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
226	indextime,what_I_plot.shape, len(all_var[index_f].shape),vertmode)
227	### If we plot a 2-D field
228	if len(what_I_plot.shape) is 2:
229	cs = contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #colors='w' )# , alpha=0.5)
230	clabel(cs,fmt = '%1.2e')
231	### If we plot a 1-D field
232	elif len(what_I_plot.shape) is 1:
233	plot(x,y,what_I_plot, zelevels, colors='g', linewidths = 0.33 ) #colors='w' )# , alpha=0.5)
234	else:
235	errormess("There is an error in reducing field !")
236
237	#### Shaded plot
238	if var:
239	what_I_plot, error = reducefield(all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert )
240	what_I_plot = what_I_plot*mult
241
242	if not error:
243	fvar = var
244	###
245	if anomaly:
246	what_I_plot = 100. * ((what_I_plot / smooth(what_I_plot,12)) - 1.)
247	fvar = 'anomaly'
248	#if mult != 1:
249	# fvar = str(mult) + "*" + var
250	###
251	if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(what_I_plot,6)
252	zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
253	if colorb in ["def","nobar"]: palette = get_cmap(name=defcolorb(fvar.upper()))
254	else: palette = get_cmap(name=colorb)
255	### If we plot a 2-D field
256	if len(what_I_plot.shape) is 2:
257	if mapmode == 0:
258	what_I_plot,x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
259	indextime,what_I_plot, len(all_var[index_f].shape),vertmode)
260	if not tile:
261	if not hole: what_I_plot = bounds(what_I_plot,zevmin,zevmax)
262	#zelevels = np.linspace(zevmin(1. + 1.e-7),zevmax(1. - 1.e-7)) #,num=20)
263	print "what_I_plot.shape", what_I_plot.shape
264	print "x.shape, y.shape", x.shape, y.shape
265	zelevels = np.linspace(zevmin,zevmax,num=ndiv)
266	#contourf(what_I_plot, zelevels, cmap = palette )
267	contourf(x,y,what_I_plot, zelevels, cmap = palette )
268	else:
269	if hole: what_I_plot = nolow(what_I_plot)
270	pcolor(x,y,what_I_plot, cmap = palette, \
271	vmin=zevmin, vmax=zevmax )
272	if colorb != 'nobar' and var != 'HGT' :
273	#subplot(111)
274	#colorbar()
275	colorbar(fraction=0.05,pad=0.03,format=fmtvar(fvar.upper()),\
276	ticks=np.linspace(zevmin,zevmax,min([ndiv,20])),\
277	extend='neither',spacing='proportional')
278	# both min max neither
279	### If we plot a 1-D field
280	elif len(what_I_plot.shape) is 1:
281	what_I_plot,x,y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
282	indextime,what_I_plot,len(all_var[index_f].shape),vertmode)
283	# x = np.array(x)
284	# y = np.array(y)
285	print "what_I_plot.shape", what_I_plot.shape
286	print "x.shape", x.shape
287	plot(x,what_I_plot)
288	### If we plot something that is not a 2-D or 1-D field
289	### (maybe plot 3-D field one day or movie ??)
290	else:
291	print "WARNING!!! ",len(what_I_plot.shape),"-D PLOT NOT SUPPORTED !!!"
292	print "field dimensions: ", what_I_plot.shape
293	exit()
294
295	else:
296	errormess("There is an error in reducing field !")
297
298
299	### Next subplot
300	plottitle = basename+' '+namefiles[index_f]
301	if typefile in ['mesoapi','meso']:
302	if addchar: plottitle = plottitle + addchar + "_LT"+str(ltst)
303	else: plottitle = plottitle + "_LT"+str(ltst)
304	if mult != 1: plottitle = str(mult) + "*" + plottitle
305	if zetitle != "fill": plottitle = zetitle
306	# if indexlon is not None:
307	# plottitle = plottitle + " lon: " + str(min(lon[indexlon])) +" "+ str(max(lon[indexlon]))
308	# if indexlat is not None:
309	# plottitle = plottitle + " lat: " + str(min(lat[indexlat])) +" "+ str(max(lat[indexlat]))
310	# if indexvert is not None:
311	# plottitle = plottitle + " vert: " + str(min(vert[indexvert])) +" "+ str(max(vert[indexvert]))
312	# if indextime is not None:
313	# plottitle = plottitle + " time: " + str(min(time[indextime])) +" "+ str(max(time[indextime]))
314	ptitle( plottitle )
315	itime += itstep
316	if nplot == numplot:
317	error = True
318	nplot += 1
319
320
321
322
323
324
325
326
327	##########################################################################
328	### Save the figure in a file in the data folder or an user-defined folder
329	if typefile in ['meso','mesoapi']: prefix = getprefix(nc)
330	elif typefile in ['gcm']: prefix = 'LMD_GCM_'
331	else: prefix = ''
332	###
333	zeplot = prefix + basename
334	if addchar: zeplot = zeplot + addchar
335	if numplot <= 0: zeplot = zeplot + "_LT"+str(abs(numplot))
336	###
337	if not target: zeplot = namefile[0:find(namefile,'wrfout')] + zeplot
338	else: zeplot = target + "/" + zeplot
339	###
340	if found_lct:
341	pad_inches_value = 0.35
342	print "save", save
343	if save == 'png':
344	if display: makeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value) #,erase=True) ## a miniature
345	makeplotres(zeplot,res=200.,pad_inches_value=pad_inches_value,disp=False)
346	elif save in ['eps','svg','pdf']:
347	makeplotres(zeplot, pad_inches_value=pad_inches_value,disp=False,ext=save)
348	elif save == 'gui':
349	show()
350	else:
351	print "save mode not supported. using gui instead."
352	show()
353	else: print "Local time not found"
354
355	###############
356	### Now the end
357	return zeplot
358
359	##############################
360	### A specific stuff for below
361	def adjust_length (tab, zelen):
362	from numpy import ones
363	if tab is None:
364	outtab = ones(zelen) * -999999
365	else:
366	if zelen != len(tab):
367	print "not enough or too much values... setting same values all variables"
368	outtab = ones(zelen) * tab[0]
369	else:
370	outtab = tab
371	return outtab
372
373	###########################################################################################
374	###########################################################################################
375	### What is below relate to running the file as a command line executable (very convenient)
376	if __name__ == "__main__":
377	import sys
378	from optparse import OptionParser ### to be replaced by argparse
379	from api_wrapper import api_onelevel
380	from netCDF4 import Dataset
381	from myplot import getlschar, separatenames, readslices
382	from os import system
383	import numpy as np
384
385	#############################
386	### Get options and variables
387	parser = OptionParser()
388	parser.add_option('-f', '--file', action='append',dest='namefile', type="string", default=None, help='[NEEDED] name of WRF file (append). Plot files separated by comas in the same figure')
389	parser.add_option('-l', '--level', action='store',dest='nvert', type="float", default=0, help='level (def=0)(-i 2: p,mbar)(-i 3,4: z,km)')
390	parser.add_option('-p', '--proj', action='store',dest='proj', type="string", default=None, help='projection')
391	parser.add_option('-b', '--back', action='store',dest='back', type="string", default=None, help='background image (def: None)')
392	parser.add_option('-t', '--target', action='store',dest='target', type="string", default=None, help='destination folder')
393	parser.add_option('-s', '--stride', action='store',dest='stride', type="int", default=3, help='stride vectors (def=3)')
394	parser.add_option('-v', '--var', action='append',dest='var', type="string", default=None, help='variable color-shaded (append)')
395	parser.add_option('-n', '--num', action='store',dest='numplot', type="int", default=2, help='plot number (def=2)(<0: plot LT -numplot)')
396	parser.add_option('-i', '--interp', action='store',dest='interp', type="int", default=None, help='interpolation (2: p, 3: z-amr, 4:z-als)')
397	parser.add_option('-c', '--color', action='store',dest='colorb', type="string", default="def", help='change colormap (nobar: no colorbar)')
398	parser.add_option('-x', '--no-vect',action='store_false',dest='winds', default=True, help='no wind vectors')
399	parser.add_option('-m', '--min', action='append',dest='vmin', type="float", default=None, help='bounding minimum value (append)')
400	parser.add_option('-M', '--max', action='append',dest='vmax', type="float", default=None, help='bounding maximum value (append)')
401	parser.add_option('-T', '--tiled', action='store_true',dest='tile', default=False, help='draw a tiled plot (no blank zone)')
402	parser.add_option('-z', '--zoom', action='store',dest='zoom', type="float", default=None, help='zoom factor in %')
403	parser.add_option('-N', '--no-api', action='store_true',dest='nocall', default=False, help='do not recreate api file')
404	parser.add_option('-d', '--display',action='store_false',dest='display', default=True, help='do not pop up created images')
405	parser.add_option('-e', '--itstep', action='store',dest='itstep', type="int", default=None, help='stride time (def=4)')
406	parser.add_option('-H', '--hole', action='store_true',dest='hole', default=False, help='holes above max and below min')
407	parser.add_option('-S', '--save', action='store',dest='save', type="string", default="gui", help='save mode (png,eps,svg,pdf or gui)(def=gui)')
408	parser.add_option('-a', '--anomaly',action='store_true',dest='anomaly', default=False, help='compute and plot relative anomaly in %')
409	parser.add_option('-w', '--with', action='store',dest='var2', type="string", default=None, help='variable contoured')
410	parser.add_option('--div', action='store',dest='ndiv', type="int", default=10, help='number of divisions in colorbar (def: 10)')
411	parser.add_option('-F', '--first', action='store',dest='first', type="int", default=1, help='first subscript to plot (def: 1)')
412	parser.add_option('--mult', action='store',dest='mult', type="float", default=1., help='a multiplicative factor to plotted field')
413	parser.add_option('--title', action='store',dest='zetitle', type="string", default="fill",help='customize the whole title')
414	#parser.add_option('-V', action='store', dest='comb', type="float", default=None, help='a defined combination of variables')
415
416	############# T.N. changes
417	#parser.add_option('-o','--operation',action='store',dest='operation',type="string", default=None ,help='matrix of operations between files (for now see code, sorry)')
418	parser.add_option('--lat', action='append',dest='slat',type="string", default=None, help='slices along latitude. One value, or two values separated by comas for averaging')
419	parser.add_option('--lon', action='append',dest='slon', type="string", default=None, help='slices along longitude. One value, or two values separated by comas for averaging')
420	parser.add_option('--vert', action='append',dest='svert',type="string", default=None, help='slices along vertical axis. One value, or two values separated by comas for averaging') # quelles coordonnees ?
421	parser.add_option('--time', action='append',dest='stime',type="string", default=None, help='slices along time. One value, or two values separated by comas for averaging') # quelles coordonnees ?
422
423	(opt,args) = parser.parse_args()
424	if opt.namefile is None:
425	print "I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h"
426	exit()
427	if opt.var is None and opt.anomaly is True:
428	print "Cannot ask to compute anomaly if no variable is set"
429	exit()
430	print "Options:", opt
431
432	#listvar = ''
433	#if opt.var is None:
434	# zerange = [-999999]
435	#else:
436	# zelen = len(opt.var)
437	# zerange = range(zelen)
438	# #if zelen == 1: listvar = opt.var[0] + ','
439	# #else :
440	# for jjj in zerange: listvar += opt.var[jjj] + ','
441	# listvar = listvar[0:len(listvar)-1]
442	# vmintab = adjust_length (opt.vmin, zelen)
443	# vmaxtab = adjust_length (opt.vmax, zelen)
444
445	print "namefile, length", opt.namefile, len(opt.namefile)
446
447	zeslat = readslices(opt.slat)
448	zeslon = readslices(opt.slon)
449	zesvert = readslices(opt.svert)
450	zestime = readslices(opt.stime)
451	print "slat,zeslat", opt.slat, zeslat
452	print "slon,zeslon", opt.slon, zeslon
453	print "svert,zesvert", opt.svert, zesvert
454	print "stime,zestime", opt.stime, zestime
455
456
457	for i in range(len(opt.namefile)):
458	for j in range(len(opt.var)):
459
460	zenamefiles = separatenames(opt.namefile[i])
461	print "zenamefiles", zenamefiles
462
463	if opt.vmin is not None : zevmin = opt.vmin[min(i,len(opt.vmin)-1)]
464	else: zevmin = None
465	if opt.vmax is not None : zevmax = opt.vmax[min(i,len(opt.vmax)-1)]
466	else: zevmax = None
467	print "vmin, zevmin", opt.vmin, zevmin
468	print "vmax, zevmax", opt.vmax, zevmax
469
470	zevar = separatenames(opt.var[j])
471	zevar = zevar[0]
472	print "var, zevar", opt.var, zevar
473
474	#checkcoherence(len(zenamefiles),len(opt.slat),len(opt.slon),len(opt.stime))
475
476	zefile = zenamefiles[0]
477
478	zelevel = opt.nvert
479	stralt = None
480	[lschar,zehour,zehourin] = getlschar ( zefile ) ## getlschar from wrfout (or simply return "" if another file)
481	#print "lschar ",lschar
482	#print "zehour ",zehour
483	#print "zehourin ",zehourin
484
485	#####################################################
486	### Call Fortran routines for vertical interpolations
487	if opt.interp is not None:
488	if opt.nvert is 0 and opt.interp is 4: zelevel = 0.010
489	### winds or no winds
490	if opt.winds : zefields = 'uvmet'
491	else : zefields = ''
492	### var or no var
493	#if opt.var is None : pass
494	if zefields == '' : zefields = listvar
495	else : zefields = zefields + "," + listvar
496	if opt.var2 is not None : zefields = zefields + "," + opt.var2
497	print zefields
498	zefile = api_onelevel ( path_to_input = '', \
499	input_name = zefile, \
500	fields = zefields, \
501	interp_method = opt.interp, \
502	onelevel = zelevel, \
503	nocall = opt.nocall )
504	print zefile
505	zelevel = 0 ## so that zelevel could play again the role of nvert
506
507
508	#############
509	### Main call
510	name = planetoplot (zenamefiles,int(zelevel),\
511	proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=zevar,\
512	numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\
513	addchar=lschar,interv=[zehour,zehourin],vmin=zevmin,vmax=zevmax,\
514	tile=opt.tile,zoom=opt.zoom,display=opt.display,\
515	itstep=opt.itstep,hole=opt.hole,save=opt.save,\
516	anomaly=opt.anomaly,var2=opt.var2,ndiv=opt.ndiv,first=opt.first,\
517	mult=opt.mult,zetitle=opt.zetitle,\
518	slon=zeslon,slat=zeslat,svert=zesvert,stime=zestime)
519	print 'Done: '+name
520	system("rm -f to_be_erased")
521
522	#########################################################
523	### Generate a .sh file with the used command saved in it
524	command = ""
525	for arg in sys.argv: command = command + arg + ' '
526	name = 'pycommand'
527	f = open(name+'.sh', 'w')
528	f.write(command)

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

Download in other formats: