Context Navigation

planetoplot.py @ 360

Last change on this file since 360 was 359, checked in by acolaitis, 14 years ago

PYTHON.
Added resolution and output filename choices.

--res X to choose X instead of 200 for png outputs
-O or --output to choose filename of the output

Added possibility to do operations between 2 files.
The plot consists of:

a main plot, given by the regular options corresponding to the choice of -var and -f.
a reference plot, given by the option --fref, in a similar way as the first plot (same var, same min/max, etc...)
a plot on which the operation between file 1 and ref file is done.

--operation to choose the operation. For now, there is + and -. It is really simple to add more.
--mope --Mope to choose min and maximum of the colorbar for the operation. Without specifying it, it is computed.
--titleref to change the title of the reference plot. To use in combination with --title for the title of the first plot.

Future work:

Add possibility to compare multiple files to the same reference file.
Add mesoscale support. (adding the options in meso.py is straightforward but some modif might have to be done in planetoplot.)

Property svn:executable set to *

File size: 18.5 KB

Line
1	#######################
2	##### PLANETOPLOT #####
3	#######################
4
5	### A. Spiga -- LMD -- 06~09/2011 -- General building and mapping capabilities
6	### T. Navarro -- LMD -- 10~11/2011 -- Improved use for GCM and added sections + 1Dplot capabilities
7	### A. Colaitis -- LMD -- -- Mostly minor improvements and inter-plot operation capabilities
8
9	def planetoplot (namefiles,\
10	level=0,\
11	vertmode=0,\
12	proj=None,\
13	back=None,\
14	target=None,
15	stride=3,\
16	numplot=None,\
17	var=None,\
18	colorb="def",\
19	winds=True,\
20	addchar=None,\
21	interv=[0,1],\
22	vmin=None,\
23	vmax=None,\
24	tile=False,\
25	zoom=None,\
26	display=True,\
27	itstep=None,\
28	hole=False,\
29	save="gui",\
30	anomaly=False,\
31	var2=None,\
32	ndiv=10,\
33	first=1,\
34	mult=1.,\
35	zetitle="fill",\
36	slon=None,\
37	slat=None,\
38	svert=None,\
39	stime=None,\
40	outputname=None,\
41	resolution=200,\
42	ope=None,\
43	fileref=None,\
44	minop=0.,\
45	maxop=0.,\
46	titleref="fill"):
47
48
49	####################################################################################################################
50	### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
51
52	#################################
53	### Load librairies and functions
54	from netCDF4 import Dataset
55	from myplot import getcoord2d,define_proj,makeplotres,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy,\
56	fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\
57	zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow,\
58	getname,localtime,polarinterv,getsindex,define_axis,determineplot
59	from mymath import deg,max,min,mean
60	from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel
61	from matplotlib.cm import get_cmap
62	import numpy as np
63	from numpy.core.defchararray import find
64
65	################################
66	### Which plot needs to be done?
67	nlon, nlat, nvert, ntime, mapmode, nslices = determineplot(slon, slat, svert, stime)
68	if mapmode == 0: winds=False
69	if not isinstance(namefiles, np.ndarray): namefiles = [namefiles]
70	zelen = len(namefiles)
71	if numplot == None: numplot = zelen*nslices
72	print "len(namefiles), nslices, numplot: ", zelen, nslices, numplot
73	if fileref is not None:
74	all_var = [[]]*(zelen+2)
75	else:
76	all_var = [[]]*zelen
77	all_var2 = [[]]*zelen
78	all_title = [[]]*zelen
79
80	#########################
81	### Loop over the files initially separated by comas to be plotted on the same figure
82	k = 0
83	firstfile = True
84	for namefile in namefiles:
85
86	######################
87	### Load NETCDF object
88	nc = Dataset(namefile)
89
90	##################################
91	### Initial checks and definitions
92	### ... TYPEFILE
93	typefile = whatkindfile(nc)
94	if firstfile:
95	typefile0 = typefile
96	elif typefile != typefile0:
97	errormess("Not the same kind of files !", [typefile0, typefile])
98	### ... VAR
99	if var not in nc.variables: var = False
100	### ... WINDS
101	if winds:
102	[uchar,vchar,metwind] = getwinddef(nc)
103	if uchar == 'not found': winds = False
104	if not var and not winds: errormess("please set at least winds or var",printvar=nc.variables)
105	### ... COORDINATES, could be moved below
106	[lon2d,lat2d] = getcoorddef(nc)
107	### ... PROJECTION
108	if proj == None: proj = getproj(nc)
109
110	##########################################################
111	if typefile == "gcm":
112	lat = nc.variables["latitude"][:]
113	lon = nc.variables["longitude"][:]
114	time = nc.variables["Time"][:]
115	vert = nc.variables["altitude"][:]
116	elif typefile in ['meso','mesoapi']:
117	if mapmode == 0:
118	if var in ['PHTOT','W']: vertdim='BOTTOM-TOP_PATCH_END_STAG'
119	else: vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
120	if var in ['V']: latdim='SOUTH-NORTH_PATCH_END_STAG'
121	else: latdim='SOUTH-NORTH_PATCH_END_UNSTAG'
122	if var in ['U']: londim='WEST-EAST_PATCH_END_STAG'
123	else: londim='WEST-EAST_PATCH_END_UNSTAG'
124	lon = np.arange(0,getattr(nc,londim),1)
125	lat = np.arange(0,getattr(nc,latdim),1)
126	if vertmode is None: vertmode=0
127	if vertmode == 0: vert = np.arange(0,getattr(nc,vertdim),1)
128	else: vert = nc.variables["vert"][:]
129	time = np.arange(0,len(nc.variables["Times"]),1)
130	#if firstfile:
131	# lat0 = lat
132	#elif len(lat0) != len(lat):
133	# errormess("Not the same latitude lengths !", [len(lat0), len(lat)])
134	#elif sum((lat == lat0) == False) != 0:
135	# errormess("Not the same latitudes !", [lat,lat0])
136	## Faire d'autre checks sur les compatibilites entre fichiers!!
137	##########################################################
138
139	if firstfile:
140	##########################
141	### Define plot boundaries
142	### todo: possible areas in latinterv in argument (ex: "Far_South_Pole")
143	if proj in ["npstere","spstere"]: [wlon,wlat] = polarinterv(lon2d,lat2d)
144	elif proj in ["lcc","laea"]: [wlon,wlat] = wrfinterv(lon2d,lat2d)
145	else: [wlon,wlat] = simplinterv(lon2d,lat2d)
146	if zoom: [wlon,wlat] = zoomset(wlon,wlat,zoom)
147
148	#########################################
149	### Name for title and graphics save file
150	basename = getname(var=var,winds=winds,anomaly=anomaly)
151	basename = basename + getstralt(nc,level) ## can be moved elsewhere for a more generic routine
152
153	print "var, var2: ", var, var2
154	if var: all_var[k] = getfield(nc,var)
155	if var2: all_var2[k] = getfield(nc,var2)
156
157	print "k", k
158	print "all_var[k].shape", all_var[k].shape
159	k += 1
160	firstfile = False
161	#### End of for namefile in namefiles
162
163	##################################
164	### Load reference file if there is one
165	if fileref is not None:
166	ncref = Dataset(fileref)
167	### ... VAR
168	if var not in nc.variables: var = False
169	if var:
170	all_var[k] = getfield(ncref,var)
171	if ope is "-":
172	all_var[k+1]= all_var[k-1] - all_var[k]
173	elif ope is "+":
174	all_var[k+1]= all_var[k-1] + all_var[k]
175	elif ope is not None:
176	print "non-implemented operation. Available op. are + and -."
177	numplot = numplot+2
178
179	##################################
180	### Open a figure and set subplots
181	fig = figure()
182	subv,subh = definesubplot( numplot, fig )
183
184	#################################
185	### Time loop for plotting device
186	found_lct = False
187	nplot = 1
188	itime = first
189	error = False
190	if itstep is None and numplot > 0: itstep = int(24./numplot)
191	elif numplot <= 0: itstep = 1
192
193	#for nplot in range(numplot):
194	while error is False:
195	print "nplot", nplot
196	print error
197
198	### Which local time ?
199	ltst = localtime ( interv[0]+itimeinterv[1], 0.5(wlon[0]+wlon[1]) )
200
201	### General plot settings
202	#print itime, int(ltst), numplot, nplot
203	if numplot >= 1:
204	if nplot > numplot: break
205	if numplot > 1:
206	if typefile not in ['geo']: subplot(subv,subh,nplot)
207	found_lct = True
208	### If only one local time is requested (numplot < 0)
209	elif numplot <= 0:
210	if int(ltst) + numplot != 0:
211	itime += 1
212	if found_lct is True: break ## because it means LT was found at previous iteration
213	else: continue ## continue to iterate to find the correct LT
214	else:
215	found_lct = True
216
217	### Map projection
218	if mapmode == 1:
219	m = define_proj(proj,wlon,wlat,back=back)
220	x, y = m(lon2d, lat2d)
221
222	####################################################################
223	if typefile in ['meso','mesoapi'] and mapmode == 1:
224	indextime = itime
225	indexlon = None
226	indexlat = None
227	indexvert = level ## ou svert ???
228	nlon = 1
229	nlat = 1
230	nvert = 1
231	ntime = 1
232	else:
233	## get all indexes to be taken into account for this subplot and then reduce field
234	## 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
235	indexlon = getsindex(slon,(nplot-1)%nlon,lon)
236	indexlat = getsindex(slat,((nplot-1)//nlon)%nlat,lat)
237	indexvert = getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert)
238	indextime = getsindex(stime,((nplot-1)//(nlonnlatnvert))%ntime,time)
239
240	if fileref is not None:
241	index_f = ((nplot-1)//(nlonnlatnvert*ntime))%(len(namefiles)+2)
242	else:
243	index_f = ((nplot-1)//(nlonnlatnvert*ntime))%len(namefiles)
244
245	print nlon, nlat, nvert, ntime
246	print slon, slat, svert, stime
247	print index_f
248	print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime)
249	####################################################################
250
251	#### Contour plot
252	if var2:
253	what_I_plot, error = reducefield(all_var2[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert )
254	#what_I_plot = what_I_plot*mult
255	if not error:
256	if mapmode == 1:
257	if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(what_I_plot,6)
258	zevmin, zevmax = calculate_bounds(what_I_plot)
259	zelevels = np.linspace(zevmin,zevmax,num=ndiv+1) #20)
260	if var2 == 'HGT': zelevels = np.arange(-10000.,30000.,2000.)
261	if mapmode == 0:
262	x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
263	indextime,what_I_plot.shape, len(all_var2[index_f].shape),vertmode)
264	### If we plot a 2-D field
265	if len(what_I_plot.shape) is 2:
266	cs = contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #0.33 colors='w' )# , alpha=0.5)
267	if typefile in ['gcm']: clabel(cs,fmt = '%1.2e')
268	### If we plot a 1-D field
269	elif len(what_I_plot.shape) is 1:
270	plot(what_I_plot,x)
271	else:
272	errormess("There is an error in reducing field !")
273
274	#### Shaded plot
275	if var:
276	what_I_plot, error = reducefield(all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert )
277	what_I_plot = what_I_plot*mult
278
279	if not error:
280	fvar = var
281	###
282	if anomaly:
283	what_I_plot = 100. * ((what_I_plot / smooth(what_I_plot,12)) - 1.)
284	fvar = 'anomaly'
285	#if mult != 1:
286	# fvar = str(mult) + "*" + var
287	###
288	if mapmode == 1:
289	if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(what_I_plot,6)
290	elif mapmode == 0:
291	what_I_plot, x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
292	indextime,what_I_plot, len(all_var[index_f].shape),vertmode)
293	if (fileref is not None) and (index_f is numplot-1):
294	zevmin, zevmax = calculate_bounds(what_I_plot,vmin=minop,vmax=maxop)
295	else:
296	zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
297	if colorb in ["def","nobar"]: palette = get_cmap(name=defcolorb(fvar.upper()))
298	else: palette = get_cmap(name=colorb)
299	##### 2D field
300	if len(what_I_plot.shape) is 2:
301	if not tile:
302	if not hole: what_I_plot = bounds(what_I_plot,zevmin,zevmax)
303	#zelevels = np.linspace(zevmin(1. + 1.e-7),zevmax(1. - 1.e-7)) #,num=20)
304	zelevels = np.linspace(zevmin,zevmax)#,num=ndiv+1)
305	#contourf(what_I_plot, zelevels, cmap = palette )
306	contourf( x, y, what_I_plot, zelevels, cmap = palette )
307	else:
308	if hole: what_I_plot = nolow(what_I_plot)
309	pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax )
310	if colorb != 'nobar' and var != 'HGT' :
311	if (fileref is not None) and (index_f is numplot-1):
312	colorbar(fraction=0.05,pad=0.03,format="%.2f",\
313	ticks=np.linspace(zevmin,zevmax,min([ndiv+1,20])),\
314	extend='neither',spacing='proportional')
315	else:
316	colorbar(fraction=0.05,pad=0.03,format=fmtvar(fvar.upper()),\
317	ticks=np.linspace(zevmin,zevmax,min([ndiv+1,20])),\
318	extend='neither',spacing='proportional')
319
320	# both min max neither
321	##### 1D field
322	elif len(what_I_plot.shape) is 1:
323	plot(x,what_I_plot)
324	#### Other cases: (maybe plot 3-D field one day or movie ??)
325	else:
326	print "WARNING!!! ",len(what_I_plot.shape),"-D PLOT NOT SUPPORTED !!!"
327	print "field dimensions: ", what_I_plot.shape
328	exit()
329
330	else:
331	errormess("There is an error in reducing field !")
332
333	### Vector plot --- a adapter
334	if winds:
335	vecx, error = reducefield( getfield(nc,uchar), d4=indextime, d3=indexvert )
336	vecy, error = reducefield( getfield(nc,vchar), d4=indextime, d3=indexvert )
337	#what_I_plot, error = reducefield(all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert )
338	if not error:
339	if typefile in ['mesoapi','meso']:
340	[vecx,vecy] = [dumpbdy(vecx,6,stag=uchar), dumpbdy(vecy,6,stag=vchar)]
341	key = True
342	elif typefile in ['gcm']:
343	key = False
344	if metwind: [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d)
345	if var == False: colorvec = definecolorvec(back)
346	else: colorvec = definecolorvec(colorb)
347	vectorfield(vecx, vecy,\
348	x, y, stride=stride, csmooth=2,\
349	#scale=15., factor=300., color=colorvec, key=key)
350	scale=20., factor=250., color=colorvec, key=key)
351	#200. ## or csmooth=stride
352
353	### Next subplot
354	if typefile in ['mesoapi','meso']:
355	plottitle = basename
356	if addchar: plottitle = plottitle + addchar + "_LT"+str(ltst)
357	else: plottitle = plottitle + "_LT"+str(ltst)
358	else:
359	if fileref is not None:
360	if index_f is numplot-1:
361	plottitle = basename+' '+"fig(1) "+ope+" fig(2)"
362	elif index_f is numplot-2:
363	plottitle = basename+' '+fileref
364	else:
365	plottitle = basename+' '+namefiles[index_f]
366	else:
367	plottitle = basename+' '+namefiles[index_f]
368	if mult != 1: plottitle = str(mult) + "*" + plottitle
369	if zetitle != "fill":
370	plottitle = zetitle
371	if titleref is "fill":
372	titleref=zetitle
373	if fileref is not None:
374	if index_f is numplot-2:
375	plottitle = titleref
376	if index_f is numplot-1:
377	plottitle = "fig(1) "+ope+" fig(2)"
378	# if indexlon is not None:
379	# plottitle = plottitle + " lon: " + str(min(lon[indexlon])) +" "+ str(max(lon[indexlon]))
380	# if indexlat is not None:
381	# plottitle = plottitle + " lat: " + str(min(lat[indexlat])) +" "+ str(max(lat[indexlat]))
382	# if indexvert is not None:
383	# plottitle = plottitle + " vert: " + str(min(vert[indexvert])) +" "+ str(max(vert[indexvert]))
384	# if indextime is not None:
385	# plottitle = plottitle + " time: " + str(min(time[indextime])) +" "+ str(max(time[indextime]))
386	ptitle( plottitle )
387	itime += itstep
388	if nplot >= numplot:
389	error = True
390	nplot += 1
391
392
393
394
395
396
397
398
399	##########################################################################
400	### Save the figure in a file in the data folder or an user-defined folder
401	if outputname is None:
402	if typefile in ['meso','mesoapi']: prefix = getprefix(nc)
403	elif typefile in ['gcm']: prefix = 'LMD_GCM_'
404	else: prefix = ''
405	###
406	zeplot = prefix + basename
407	if addchar: zeplot = zeplot + addchar
408	if numplot <= 0: zeplot = zeplot + "_LT"+str(abs(numplot))
409	###
410	if not target: zeplot = namefile[0:find(namefile,'wrfout')] + zeplot
411	else: zeplot = target + "/" + zeplot
412	###
413	else:
414	zeplot=outputname
415
416	if found_lct:
417	pad_inches_value = 0.35
418	print "save", save
419	if save == 'png':
420	if display: makeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value) #,erase=True) ## a miniature
421	makeplotres(zeplot,res=resolution,pad_inches_value=pad_inches_value,disp=False)
422	elif save in ['eps','svg','pdf']:
423	makeplotres(zeplot, pad_inches_value=pad_inches_value,disp=False,ext=save)
424	elif save == 'gui':
425	show()
426	else:
427	print "save mode not supported. using gui instead."
428	show()
429	else: print "Local time not found"
430
431	###############
432	### Now the end
433	return zeplot

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

Context Navigation

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

Download in other formats: