Context Navigation

source: trunk/UTIL/PYTHON/planetoplot_v2/ppclass.py @ 1211

Last change on this file since 1211 was 1201, checked in by aslmd, 11 years ago
planetoplot. fix a bug introduced with getfd stuff.
File size: 101.0 KB

Line
1	###############################################
2	## PLANETOPLOT ##
3	## --> PPCLASS ##
4	## A generic and versatile Python module ##
5	## ... to read netCDF files and plot ##
6	###############################################
7	## Author: Aymeric Spiga. 02-03/2013 ##
8	###############################################
9	# python built-in librairies
10	import os
11	import time as timelib
12	import pickle
13	# added librairies
14	import numpy as np
15	import netCDF4
16	import matplotlib.pyplot as mpl
17	# personal librairies
18	import ppplot
19	import ppcompute
20	###############################################
21
22	###################################
23	#### HEADER ##
24	#### ... executed when imported ##
25	###################################
26	# where settings files are located...
27	# ... this can be hardcoded here
28	whereset = None
29	whereset = ppcompute.findset(whereset)
30	# ... we load user-defined automatic settings from set_ppclass.txt
31	zefile = "set_ppclass.txt"
32	glob_listx = [] ; glob_listy = [] ; glob_listz = [] ; glob_listt = []
33	glob_listarea = []
34	try:
35	f = open(whereset+zefile, 'r') ; lines = f.readlines()
36	for stuff in lines[5].strip().split(';'): glob_listx.append(stuff)
37	for stuff in lines[8].strip().split(';'): glob_listy.append(stuff)
38	for stuff in lines[11].strip().split(';'): glob_listz.append(stuff)
39	for stuff in lines[14].strip().split(';'): glob_listt.append(stuff)
40	for stuff in lines[17].strip().split(';'): glob_listarea.append(stuff)
41	except IOError:
42	print "PPCLASS warning: "+zefile+" not in "+whereset+" ; no presets."
43
44	##################################
45	#### USEFUL GENERIC FUNCTIONS ####
46	##################################
47
48	# inspect variables and dimensions in a netCDF file
49	def inspect(filename):
50	print "---------------------------------------------------------------"
51	print "**** INSPECT FILE",filename
52	test = netCDF4.Dataset(filename)
53	print "---------------------------------------------------------------"
54	print "**** DIMENSIONS:"
55	for obj in test.dimensions.values(): print obj
56	print "---------------------------------------------------------------"
57	print "**** VARIABLES: ",test.variables.keys()
58	print "---------------------------------------------------------------"
59	findinlist(test,glob_listx,"**** FOUND X-AXIS ---> ")
60	findinlist(test,glob_listy,"**** FOUND Y-AXIS ---> ")
61	findinlist(test,glob_listz,"**** FOUND Z-AXIS ---> ")
62	findinlist(test,glob_listt,"**** FOUND T-AXIS ---> ")
63	print "**** ( according to settings in "+whereset+zefile+" )"
64	print "---------------------------------------------------------------"
65	# -- function defined for the above function inspect
66	def findinlist(netcdfobj,extlist,message):
67	found = False
68	for c in extlist:
69	if c in netcdfobj.variables.keys():
70	found = True
71	output = message+str(c)
72	if c in netcdfobj.dimensions.keys():
73	output = output+" ---- has "+str(len(netcdfobj.dimensions[c]))+" values"
74	try: output = output + " ---- from "+str(netcdfobj.variables[c][0])+" to "+str(netcdfobj.variables[c][-1])
75	except: pass
76	print output ; output = ""
77	if not found:
78	print message+"not found. will simply use index (check out dimensions)."
79
80	# request a given attribute (and test if it is here)
81	def ncattr(filename,char):
82	nc = netCDF4.Dataset(filename)
83	if hasattr(nc,char): ncattr=getattr(nc,char)
84	return ncattr
85
86	# check a tab and exit if wrong. if just one string make it a list.
87	# (if allownumber, convert this into a string).
88	def checktab(tab,mess="",allownone=False,allownumber=False):
89	if tab is None:
90	if not allownone: print "pp.define: no "+mess ; exit()
91	else: pass
92	else:
93	if not isinstance(tab, list):
94	if isinstance(tab, str):
95	tab = [tab]
96	elif (isinstance(tab, int) or isinstance(tab, float)) and allownumber:
97	tab = [str(tab)]
98	else:
99	print "pp.define: "+mess+" should be either a string or a list of strings!" ; exit()
100	elif isinstance(tab, list):
101	if isinstance(tab[0],str):
102	pass
103	elif (isinstance(tab[0], int) or isinstance(tab[0], float)) and allownumber:
104	for iii in range(len(tab)): tab[iii] = str(tab[iii])
105	else:
106	print "pp.define: "+mess+" should be either a string or a list of strings!" ; exit()
107	return tab
108
109	# determine which method is to be applied to a given dimension
110	def findmethod(tab):
111	if tab is None: output = "free"
112	elif tab[0,0] != tab[0,1]: output = "comp"
113	else: output = "fixed"
114	return output
115
116	# read what is given by the user (version of T. Navarro)
117	def readslices(saxis):
118	if saxis == None:
119	zesaxis = None
120	else:
121	zesaxis = np.empty((len(saxis),2))
122	for i in range(len(saxis)):
123	a = separatenames(saxis[i])
124	if len(a) == 1:
125	zesaxis[i,:] = float(a[0])
126	else:
127	zesaxis[i,0] = float(a[0])
128	zesaxis[i,1] = float(a[1])
129	return zesaxis
130
131	# look for comas in the input name to separate different names (files, variables,etc ..)
132	# (needed by readslices)
133	def separatenames (name):
134	if name is None: names = None
135	else:
136	names = [] ; stop = 0 ; currentname = name
137	while stop == 0:
138	indexvir = currentname.find(',')
139	if indexvir == -1: stop = 1 ; name1 = currentname
140	else: name1 = currentname[0:indexvir]
141	names = np.concatenate((names,[name1]))
142	currentname = currentname[indexvir+1:len(currentname)]
143	return names
144
145	#######################
146	### THE MAIN OBJECT ###
147	#######################
148	class pp():
149
150	# print out a help string when help is invoked on the object
151	def __repr__(self):
152	whatprint = 'pp object. \"help(pp)\" for more information\n'
153	return whatprint
154
155	# default settings
156	# -- user can define settings by two methods.
157	# -- 1. yeah = pp(file="file.nc")
158	# -- 2. yeah = pp() ; yeah.file = "file.nc"
159	def __init__(self,file=None,var="notset",\
160	filegoal=None,vargoal=None,\
161	x=None,y=None,z=None,t=None,\
162	sx=1,sy=1,\
163	sz=1,st=1,\
164	svx=1,\
165	svy=1,\
166	compute="mean",\
167	kind3d="tyx",\
168	verbose=False,\
169	quiet=False,\
170	noproj=False,\
171	superpose=False,\
172	plotin=None,\
173	forcedimplot=-1,\
174	out="gui",\
175	filename="myplot",\
176	folder="./",\
177	includedate=True,\
178	res=150.,\
179	xlabel=None,ylabel=None,\
180	xcoeff=None,ycoeff=None,\
181	xmin=None,ymin=None,\
182	xmax=None,ymax=None,\
183	nxticks=10,nyticks=10,\
184	proj=None,\
185	vmin=None,vmax=None,\
186	div=None,\
187	colorbar=None,\
188	linestyle=None,\
189	marker=None,\
190	color=None,\
191	legend=None,\
192	changetime=None,\
193	units=None,\
194	savtxt=False,\
195	modx=None,\
196	fmt=None,\
197	xp=16,yp=8,\
198	missing=1.e25,\
199	trans=None,back=None,\
200	showcb=None,\
201	logy=None,\
202	title=None):
203	self.request = None
204	self.nrequest = 0
205	self.nfin = 0 ; self.nvin = 0
206	self.nplotx = None ; self.nploty = None
207	self.nplotz = None ; self.nplott = None
208	self.status = "init"
209	self.fig = None ; self.subv = None ; self.subh = None
210	self.n = 0 ; self.howmanyplots = 0
211	self.nplot = 0
212	self.p = None
213	self.customplot = False
214	self.f = None
215	self.l = None
216	self.xx = None
217	self.yy = None
218	self.zz = None
219	self.tt = None
220	## what could be defined by the user
221	self.file = file
222	self.var = var
223	self.filegoal = filegoal
224	self.vargoal = vargoal
225	self.x = x ; self.y = y ## if None, free dimension
226	self.z = z ; self.t = t ## if None, free dimension
227	self.sx = sx ; self.sy = sy
228	self.sz = sz ; self.st = st
229	self.svx = svx
230	self.svy = svy
231	self.compute = compute
232	self.kind3d = kind3d
233	self.verbose = verbose
234	self.quiet = quiet
235	self.noproj = noproj
236	self.plotin = plotin
237	self.superpose = superpose
238	self.forcedimplot = forcedimplot
239	self.out = out
240	self.filename = filename
241	self.res = res
242	self.folder = folder
243	self.includedate = includedate
244	self.changetime = changetime
245	self.savtxt = savtxt
246	self.modx = modx
247	self.missing = missing
248	## here are user-defined plot settings
249	## -- if not None, valid on all plots in the pp() objects
250	self.xlabel = xlabel ; self.xcoeff = xcoeff
251	self.ylabel = ylabel ; self.ycoeff = ycoeff
252	self.xmin = xmin ; self.xmax = xmax
253	self.ymin = ymin ; self.ymax = ymax
254	self.proj = proj
255	self.vmin = vmin ; self.vmax = vmax
256	self.div = div
257	self.colorbar = colorbar
258	self.linestyle = linestyle
259	self.marker = marker
260	self.color = color
261	self.legend = legend
262	self.units = units
263	self.title = title
264	self.xp = xp ; self.yp = yp
265	self.nxticks = nxticks ; self.nyticks = nyticks
266	self.fmt = fmt
267	self.trans = trans ; self.back = back
268	self.showcb = showcb
269	self.logy = logy
270
271	# print status
272	def printstatus(self):
273	if not self.quiet:
274	if self.filename == "THIS_IS_A_CLONE":
275	pass
276	else:
277	print "**** PPCLASS. Done step: " + self.status
278
279	# print attributes
280	def printme(self):
281	for k, v in vars(self).items():
282	print k,v
283
284	#####################################################
285	# EMULATE OPERATORS + - * / ** << FOR PP() OBJECTS #
286	#####################################################
287
288	# define the operation <<
289	# ... e.g. obj2 << obj1
290	# ... means: get init for pp object obj2 from another pp object obj1
291	# ... (this should solve the affectation trap obj2 = obj1)
292	def __lshift__(self,other):
293	if other.__class__.__name__ == "pp":
294	self.file = other.file
295	self.var = other.var
296	self.filegoal = other.filegoal
297	self.vargoal = other.vargoal
298	self.x = other.x ; self.y = other.y ## if None, free dimension
299	self.z = other.z ; self.t = other.t ## if None, free dimension
300	self.sx = other.sx ; self.sy = other.sy
301	self.sz = other.sz ; self.st = other.st
302	self.compute = other.compute
303	self.kind3d = other.kind3d
304	self.verbose = other.verbose
305	self.noproj = other.noproj
306	self.plotin = other.plotin
307	self.superpose = other.superpose
308	self.forcedimplot = other.forcedimplot
309	self.out = other.out
310	self.filename = other.filename
311	self.folder = other.folder
312	self.xlabel = other.xlabel ; self.xcoeff = other.xcoeff
313	self.ylabel = other.ylabel ; self.ycoeff = other.ycoeff
314	self.proj = other.proj
315	self.vmin = other.vmin ; self.vmax = other.vmax
316	self.xmin = other.xmin ; self.xmax = other.xmax
317	self.ymin = other.ymin ; self.ymax = other.ymax
318	self.div = other.div
319	self.colorbar = other.colorbar
320	self.linestyle = other.linestyle
321	self.marker = other.marker
322	self.color = other.color
323	self.legend = other.legend
324	self.units = other.units
325	self.title = other.title
326	self.includedate = other.includedate
327	self.changetime = other.changetime
328	self.savtxt = other.savtxt
329	self.modx = other.modx
330	self.xp = other.xp ; self.yp = other.yp
331	self.missing = other.missing
332	self.nxticks = other.nxticks ; self.nyticks = other.nyticks
333	self.fmt = other.fmt
334	self.trans = other.trans ; self.back = other.back
335	self.showcb = other.showcb
336	self.logy = other.logy
337	else:
338	print "!! ERROR !! argument must be a pp object." ; exit()
339
340	# check the compatibility of two objects for operations
341	# --> if other is a pp class, test sizes and return isnum = False
342	# --> if other is an int or a float, return isnum = True
343	# --> otherwise, just print an error and exit
344	def checktwo(self,other):
345	if other.__class__.__name__ == "pp":
346	isnum = False
347	if self.status in ["init","defined"] or other.status in ["init","define"]:
348	print "!! ERROR !! Please use .retrieve to get fields for plots with one of your pp operands." ; exit()
349	if self.nfin != other.nfin or \
350	self.nvin != other.nvin or \
351	self.nplott != other.nplott or \
352	self.nplotz != other.nploty or \
353	self.nploty != other.nploty or \
354	self.nplotx != other.nplotx :
355	print "!! ERROR !! The two operands do not have the same number of files, variables, t z y x requests."
356	exit()
357	elif isinstance(other,int) or isinstance(other,float):
358	isnum = True
359	else:
360	print "!! ERROR !! The operand is neither a pp class nor an integer or a float." ; exit()
361	return isnum
362
363	# define a selective copy of a pp() object for operations
364	# ... copy.copy() is not conservative (still acts like a pointer)
365	# ... copy.deepcopy() does not work with netCDF objects
366	# so what is done here is a copy of everything except
367	# (to avoid sharing with self and therefore modifying self through operations)
368	# - request attribute of pp() object
369	# - field attribute of the onerequest() objects
370	def selective_copy(self):
371	if self.status in ["init","defined"]:
372	print "!! ERROR !! Please use .retrieve to get fields for the object you want to copy from." ; exit()
373	the_clone = pp()
374	for k, v in vars(self).items():
375	if k != "request":
376	setattr(the_clone,k,v)
377	the_clone.verbose = False
378	the_clone.filename = "THIS_IS_A_CLONE" # trick to avoid additional outputs
379	the_clone.define()
380	for i in range(self.nfin):
381	for j in range(self.nvin):
382	for t in range(self.nplott):
383	for z in range(self.nplotz):
384	for y in range(self.nploty):
385	for x in range(self.nplotx):
386	obj_ref = self.request[i][j][t][z][y][x]
387	obj = the_clone.request[i][j][t][z][y][x]
388	for k, v in vars(obj_ref).items():
389	if k != "field":
390	setattr(obj,k,v)
391	the_clone.status = "retrieved"
392	the_clone.filename = self.filename
393	return the_clone
394
395	# define the operation + on two objects. or with an int/float.
396	# ... with selective_copy the self object is not modified.
397	def __add__(self,other):
398	isnum = self.checktwo(other)
399	the_clone = self.selective_copy()
400	for i in range(self.nfin):
401	for j in range(self.nvin):
402	for t in range(self.nplott):
403	for z in range(self.nplotz):
404	for y in range(self.nploty):
405	for x in range(self.nplotx):
406	obj = the_clone.request[i][j][t][z][y][x]
407	obj_ref = self.request[i][j][t][z][y][x]
408	if not isnum:
409	ope = other.request[i][j][t][z][y][x].field
410	if ope.ndim == 0:
411	ope = float(ope) # if no dimension then this means that ope is a single value (not to be kept as an array)
412	elif obj_ref.field.shape != ope.shape:
413	print "!! ERROR !! The two fields for operation do not have the same shape.",obj_ref.field.shape,ope.shape
414	exit()
415	else:
416	ope = other
417	goal = self.vargoal[j] + self.filegoal[i]
418	if ("vector" in goal) or ("contour" in goal):
419	if self.verbose: print "!! WARNING !! No operation was made on contours and vectors. This can be debatted actually."
420	obj.field = obj_ref.field
421	else:
422	obj.field = obj_ref.field + ope
423	return the_clone
424
425	# define the operation - on two objects. or with an int/float.
426	# ... with selective_copy the self object is not modified.
427	def __sub__(self,other):
428	isnum = self.checktwo(other)
429	the_clone = self.selective_copy()
430	for i in range(self.nfin):
431	for j in range(self.nvin):
432	for t in range(self.nplott):
433	for z in range(self.nplotz):
434	for y in range(self.nploty):
435	for x in range(self.nplotx):
436	obj = the_clone.request[i][j][t][z][y][x]
437	obj_ref = self.request[i][j][t][z][y][x]
438	if not isnum:
439	ope = other.request[i][j][t][z][y][x].field
440	if ope.ndim == 0:
441	ope = float(ope) # if no dimension then this means that ope is a single value (not to be kept as an array)
442	elif obj_ref.field.shape != ope.shape:
443	print "!! ERROR !! The two fields for operation do not have the same shape.",obj_ref.field.shape,ope.shape
444	exit()
445	else:
446	ope = other
447	goal = self.vargoal[j] + self.filegoal[i]
448	if ("vector" in goal) or ("contour" in goal):
449	if self.verbose: print "!! WARNING !! No operation was made on contours and vectors. This can be debatted actually."
450	obj.field = obj_ref.field
451	else:
452	obj.field = obj_ref.field - ope
453	return the_clone
454
455	# define the operation * on two objects. or with an int/float.
456	# ... with selective_copy the self object is not modified.
457	def __mul__(self,other):
458	isnum = self.checktwo(other)
459	the_clone = self.selective_copy()
460	for i in range(self.nfin):
461	for j in range(self.nvin):
462	for t in range(self.nplott):
463	for z in range(self.nplotz):
464	for y in range(self.nploty):
465	for x in range(self.nplotx):
466	obj = the_clone.request[i][j][t][z][y][x]
467	obj_ref = self.request[i][j][t][z][y][x]
468	if not isnum:
469	ope = other.request[i][j][t][z][y][x].field
470	if ope.ndim == 0:
471	ope = float(ope) # if no dimension then this means that ope is a single value (not to be kept as an array)
472	elif obj_ref.field.shape != ope.shape:
473	print "!! ERROR !! The two fields for operation do not have the same shape.",obj_ref.field.shape,ope.shape
474	exit()
475	else:
476	ope = other
477	goal = self.vargoal[j] + self.filegoal[i]
478	if ("vector" in goal) or ("contour" in goal):
479	if self.verbose: print "!! WARNING !! No operation was made on contours and vectors. This can be debatted actually."
480	obj.field = obj_ref.field
481	else:
482	obj.field = obj_ref.field * ope
483	return the_clone
484
485	# define the operation / on two objects. or with an int/float.
486	# ... with selective_copy the self object is not modified.
487	def __div__(self,other):
488	isnum = self.checktwo(other)
489	the_clone = self.selective_copy()
490	for i in range(self.nfin):
491	for j in range(self.nvin):
492	for t in range(self.nplott):
493	for z in range(self.nplotz):
494	for y in range(self.nploty):
495	for x in range(self.nplotx):
496	obj = the_clone.request[i][j][t][z][y][x]
497	obj_ref = self.request[i][j][t][z][y][x]
498	if not isnum:
499	ope = other.request[i][j][t][z][y][x].field
500	if ope.ndim == 0:
501	ope = float(ope) # if no dimension then this means that ope is a single value (not to be kept as an array)
502	elif obj_ref.field.shape != ope.shape:
503	print "!! ERROR !! The two fields for operation do not have the same shape.",obj_ref.field.shape,ope.shape
504	exit()
505	else:
506	ope = other
507	goal = self.vargoal[j] + self.filegoal[i]
508	if ("vector" in goal) or ("contour" in goal):
509	if self.verbose: print "!! WARNING !! No operation was made on contours and vectors. This can be debatted actually."
510	obj.field = obj_ref.field
511	else:
512	obj.field = obj_ref.field / ope
513	return the_clone
514
515	# define the reverse operation float/int + object
516	def __radd__(self,other):
517	isnum = self.checktwo(other)
518	if not isnum: print "!! ERROR !! Operand should be a number" ; exit()
519	return self.__add__(other)
520
521	# define the reverse operation float/int - object
522	def __rsub__(self,other):
523	isnum = self.checktwo(other)
524	if not isnum: print "!! ERROR !! Operand should be a number" ; exit()
525	return self.__sub__(other)
526
527	# define the reverse operation float/int * object
528	def __rmul__(self,other):
529	isnum = self.checktwo(other)
530	if not isnum: print "!! ERROR !! Operand should be a number" ; exit()
531	return self.__mul__(other)
532
533	# define the reverse operation float/int / object
534	def __rdiv__(self,other):
535	isnum = self.checktwo(other)
536	if not isnum: print "!! ERROR !! Operand should be a number" ; exit()
537	return self.__div__(other)
538
539	# define the operation ** on one object.
540	# ... with selective_copy the self object is not modified.
541	def __pow__(self,num):
542	the_clone = self.selective_copy()
543	if isinstance(num,int) or isinstance(num,float):
544	for i in range(self.nfin):
545	for j in range(self.nvin):
546	for t in range(self.nplott):
547	for z in range(self.nplotz):
548	for y in range(self.nploty):
549	for x in range(self.nplotx):
550	obj = the_clone.request[i][j][t][z][y][x]
551	obj_ref = self.request[i][j][t][z][y][x]
552	goal = self.vargoal[j] + self.filegoal[i]
553	if ("vector" in goal) or ("contour" in goal):
554	if self.verbose: print "!! WARNING !! No operation was made on contours and vectors. This can be debatted actually."
555	obj.field = obj_ref.field
556	else:
557	obj.field = obj_ref.field ** num
558	else:
559	print "!! ERROR !! To define a power, either an int or a float is needed." ; exit()
560	return the_clone
561
562	### TBD: reverse power? for exponentials?
563
564	##############################################################################################
565	# define method
566	# ---------
567	# ... (file and var are either one string or a vector of strings)
568	# ... the goal of define is to define a 2D array of onerequest() objects (see class below)
569	# given the number of file, var, x, y, z, t asked by the user
570	# ... objectives for file or var are given through filegoal and vargoal
571	# --> possible values: main contour vector
572	# ---------
573	# ... then onerequest() objects are being defined more precisely
574	# by getting index_x index_y index_z index_t
575	# and setting method_x method_y method_z method_t to either
576	# - "free" for free dimensions (plot dimensions)
577	# - "comp" for averages, max, min
578	# - "fixed" for fixed dimensions (possibly several i.e. multislice)
579	##############################################################################################
580	def define(self):
581	self.printstatus()
582	# initial check and get dimensions
583	self.file = checktab(self.file,mess="file")
584	self.nfin = len(self.file)
585	if self.verbose:
586	for i in range(self.nfin): inspect(self.file[i])
587	self.var = checktab(self.var,mess="var")
588	self.nvin = len(self.var)
589	# check goal tabs for files and variables
590	# ... default is to plot everything
591	if self.filegoal is None: self.filegoal = ["main"]*self.nfin
592	if self.vargoal is None: self.vargoal = ["main"]*self.nvin
593	self.filegoal = checktab(self.filegoal, mess="filegoal")
594	self.vargoal = checktab(self.vargoal, mess="vargoal")
595	if len(self.filegoal) != self.nfin: print "!! ERROR !! filegoal must be the same size as file." ; exit()
596	if len(self.vargoal) != self.nvin: print "!! ERROR !! vargoal must be the same size as var." ; exit()
597	# variables: initial check
598	self.x = checktab(self.x,mess="x",allownone=True,allownumber=True)
599	self.y = checktab(self.y,mess="y",allownone=True,allownumber=True)
600	self.z = checktab(self.z,mess="z",allownone=True,allownumber=True)
601	self.t = checktab(self.t,mess="t",allownone=True,allownumber=True)
602	# for the moment not var- nor file- dependent.
603	# but this could be the case.
604	sx = readslices(self.x) ; sy = readslices(self.y)
605	sz = readslices(self.z) ; st = readslices(self.t)
606	# get methods
607	mx = findmethod(sx) ; my = findmethod(sy)
608	mz = findmethod(sz) ; mt = findmethod(st)
609	# get number of plots to be done
610	if mx in ["fixed","comp"]: self.nplotx = sx.size/2
611	else: self.nplotx = 1
612	if my in ["fixed","comp"]: self.nploty = sy.size/2
613	else: self.nploty = 1
614	if mz in ["fixed","comp"]: self.nplotz = sz.size/2
615	else: self.nplotz = 1
616	if mt in ["fixed","comp"]: self.nplott = st.size/2
617	else: self.nplott = 1
618	if self.verbose: print "**** OK. Plots over x,y,z,t -->",self.nplotx,self.nploty,self.nplotz,self.nplott
619	# create the list of onerequest() objects
620	self.request = [[[[[[ \
621	onerequest() \
622	for x in range(self.nplotx)] for y in range(self.nploty)] \
623	for z in range(self.nplotz)] for t in range(self.nplott)] \
624	for j in range(self.nvin)] for i in range(self.nfin)]
625	# store how many onerequest() objects are in self.request
626	self.nrequest = self.nfinself.nvinself.nplotxself.nplotyself.nplotz*self.nplott
627	# loop on onerequest() objects
628	for i in range(self.nfin):
629	for j in range(self.nvin):
630	for t in range(self.nplott):
631	for z in range(self.nplotz):
632	for y in range(self.nploty):
633	for x in range(self.nplotx):
634	obj = self.request[i][j][t][z][y][x]
635	# fill in names for files and variables
636	obj.verbose = self.verbose
637	obj.file = self.file[i]
638	obj.var = self.var[j]
639	# get methods
640	obj.method_x = mx ; obj.method_y = my
641	obj.method_z = mz ; obj.method_t = mt
642	# indicate the computation method
643	obj.compute = self.compute
644	# indicate the kind of 3D fields
645	obj.kind3d = self.kind3d
646	# open the files (the same file might be opened several times but this is cheap)
647	obj.openfile()
648	### get x,y,z,t dimensions from file
649	obj.getdim()
650	### possible time axis change
651	obj.changetime = self.changetime
652	obj.performtimechange()
653	# get strides
654	obj.sx = self.sx ; obj.sy = self.sy
655	obj.sz = self.sz ; obj.st = self.st
656	### get index
657	obj.getindextime(dalist=st,ind=t)
658	obj.getindexvert(dalist=sz,ind=z)
659	obj.getindexhori(dalistx=sx,dalisty=sy,indx=x,indy=y)
660	# missing value
661	obj.missing = self.missing
662	# change status
663	self.status = "defined"
664	return self
665
666	##############################################################################################
667	# retrieve method
668	# --> for each element onerequest() in the array, get field .var from .f file
669	# --> see below the onerequest() class:
670	# - only get what is needed for computing and plotting
671	# - averages etc... are computed here
672	# --> RESULT: each onerequest() object has now its attribute .field filled
673	# --> if one wants to perform operations on fields, this should be done after retrieve()
674	##############################################################################################
675	def retrieve(self):
676	self.printstatus()
677	# check if things were done OK before
678	if self.status != "defined": print "!! ERROR !! Please use .define() to define your pp object." ; exit()
679	## create the list of f() and l() objects
680	## --> so that the user can easily access values (and labels for easy exploration)
681	## --> see example easy_get_field
682	self.f = [ [] for iii in range(self.nrequest) ]
683	self.l = [ [] for iii in range(self.nrequest) ]
684	self.xx = [ [] for iii in range(self.nrequest) ]
685	self.yy = [ [] for iii in range(self.nrequest) ]
686	self.zz = [ [] for iii in range(self.nrequest) ]
687	self.tt = [ [] for iii in range(self.nrequest) ]
688	count = 0
689	## first get fields
690	## ... only what is needed is extracted from the files
691	## ... and computations are performed
692	for i in range(self.nfin):
693	for j in range(self.nvin):
694	for t in range(self.nplott):
695	for z in range(self.nplotz):
696	for y in range(self.nploty):
697	for x in range(self.nplotx):
698	obj = self.request[i][j][t][z][y][x]
699	obj.getfield()
700	if self.compute != "nothing":
701	obj.computations()
702	# save fields in self.f for the user
703	self.f[count] = obj.field
704	self.xx[count] = obj.field_x
705	self.yy[count] = obj.field_y
706	self.zz[count] = obj.field_z
707	self.tt[count] = obj.field_t
708	# save a legend in self.l for the user
709	self.l[count] = "_"
710	if self.nfin > 1: self.l[count] = self.l[count] + "f=#"+str(int(i+1))+'_'
711	if self.nvin > 1: self.l[count] = self.l[count] + "v="+obj.var+'_'
712	if self.nplotx > 1: self.l[count] = self.l[count] + "x="+str(self.x[x])+'_'
713	if self.nploty > 1: self.l[count] = self.l[count] + "y="+str(self.y[y])+'_'
714	if self.nplotz > 1: self.l[count] = self.l[count] + "z="+str(self.z[z])+'_'
715	if self.nplott > 1: self.l[count] = self.l[count] + "t="+str(self.t[t])+'_'
716	# close the file
717	obj.closefile()
718	count = count + 1
719	## make it simple: self.f is simply the data array if self.nrequest=1
720	if self.nrequest == 1:
721	self.f = self.f[0]
722	self.xx = self.xx[0]
723	self.yy = self.yy[0]
724	self.zz = self.zz[0]
725	self.tt = self.tt[0]
726	# change status
727	self.status = "retrieved"
728	return self
729
730	##########################################################
731	# get: a shortcut method for the define + retrieve chain #
732	##########################################################
733	def get(self):
734	self.define()
735	self.retrieve()
736	return self
737
738	###########################################################
739	# getf: a shortcut method for the define + retrieve chain #
740	# ... in which the output is self.f #
741	# ... and the ppclass is kept quiet #
742	###########################################################
743	def getf(self):
744	self.quiet = True
745	self.get()
746	return self.f
747
748	############################################################
749	# getfd: a shortcut method for the define + retrieve chain #
750	# ... in which the output is self.f,x,y,z,t #
751	# ... and the ppclass is kept quiet #
752	############################################################
753	def getfd(self):
754	self.quiet = True
755	self.get()
756	return self.f,self.xx,self.yy,self.zz,self.tt
757
758	############################################################
759	# getfl: a shortcut method for the define + retrieve chain #
760	# ... in which the output is self.f, self.l #
761	# ... and the ppclass is kept quiet #
762	############################################################
763	def getfl(self):
764	self.quiet = True
765	self.get()
766	return self.f,self.l
767
768	########################################
769	# smooth: smooth the field in 1D or 2D #
770	########################################
771	## TBD: smooth not OK with masked array in the end of retrieve()
772	def smooth(self,window):
773	if self.verbose:
774	print "!! WARNING !! Performing a smoothing with a window size",window
775	print "!! WARNING !! To come back to unsmoothed file, use .get() again"
776	for i in range(self.nfin):
777	for j in range(self.nvin):
778	for t in range(self.nplott):
779	for z in range(self.nplotz):
780	for y in range(self.nploty):
781	for x in range(self.nplotx):
782	obj = self.request[i][j][t][z][y][x]
783	if obj.field.ndim == 1:
784	print "!! ERROR !! 1D smoothing not supported yet because reduces array sizes."
785	exit()
786	# TBD TBD TBD
787	#obj.field = ppcompute.smooth1d(obj.field,window=window)
788	elif obj.field.ndim == 2:
789	obj.field = ppcompute.smooth2d(obj.field,window=window)
790	#obj.field = ppcompute.smooth2diter(obj.field,n=window)
791
792	##############################################################################################
793	# defineplot method
794	# --> defineplot first defines arrays of plot objects and set each of them
795	# ... simple looping except cases where goal is not main (e.g. contour or vector)
796	# --> principle: each onerequest() object gives birth to a subplot
797	# --> defineplot vs. makeplot: defining plot and actually plotting it are clearly separated
798	# --> THE KEY OUPUT OF defineplot IS AN ARRAY self.p OF PLOT OBJECTS
799	# optional arguments
800	# --> extraplot: to indicate a number of plots to be added afterwards (use self.plotin)
801	# --> loadfile: to use self.p from a previously saved file
802	##############################################################################################
803	def defineplot(self,extraplot=0,loadfile=None):
804	# -----------------------------------------------------
805	# LOAD MODE: load a self.p object. count plots from it.
806	# -----------------------------------------------------
807	if loadfile is not None:
808	try: filehandler = open(loadfile, 'r') ; self.p = pickle.load(filehandler)
809	except IOError: print "!! ERROR !! Cannot find object file to load." ; exit()
810	self.status = "definedplot" ; self.plotin = None
811	self.nplot = len(self.p) ; self.howmanyplots = self.nplot
812	## [BUG FIX: apparently info about missing values is not saved correctly]
813	for count in range(self.nplot):
814	pl = self.p[count]
815	masked = np.ma.masked_where(np.abs(pl.f) > self.missing,pl.f)
816	pl.f = masked ; pl.f[pl.f.mask] = np.NaN
817	return #self?
818	# -----------------------------------------------------
819	# REGULAR MODE
820	# -----------------------------------------------------
821	self.printstatus()
822	# check if things were done OK before
823	if self.status in ["init","defined"]:
824	print "!! ERROR !! Please use .retrieve() to get fields for plots with your pp object." ; exit()
825	# check self.plotin (an existing fig on which to add plots afterwards)
826	if self.plotin.__class__.__name__ == "pp":
827	if self.plotin.fig is None:
828	self.plotin = None # this is an additional security in case
829	# a pp object is given without figure opened yet.
830	elif self.plotin is not None:
831	print "!! ERROR !! plotin argument must be a pp object." ; exit()
832	# initialize the array of subplot objects
833	# either something new or attributes coming from plotin object
834	if self.plotin is None: self.p = [ ]
835	else: self.p = self.plotin.p
836	# create an array of subplot objects
837	# ... in theory the order of looping can be changed without any harm
838	# ... the only important thing is to keep i,j,t,z,y,x resp. for file,var,t,z,y,x
839	count = 0
840	for i in range(self.nfin):
841	if self.filegoal[i] == "main":
842	for j in range(self.nvin):
843	if self.vargoal[j] == "main":
844	for t in range(self.nplott):
845	for z in range(self.nplotz):
846	for y in range(self.nploty):
847	for x in range(self.nplotx):
848	# look at dimension and append the right plot object
849	obj = self.request[i][j][t][z][y][x]
850	dp = obj.dimplot
851	if dp == 1 or self.forcedimplot == 1: plobj = ppplot.plot1d()
852	elif dp == 2 or self.forcedimplot == 2: plobj = ppplot.plot2d()
853	elif dp == 0: print "**** OK. VALUES VALUES VALUES",obj.field
854	else: print "!! ERROR !! 3D or 4D plots not supported" ; exit()
855	# load abscissa and ordinate in obj
856	obj.definecoord()
857	# we start to define things here before appending
858	# (convenient: could be overridden by the user before makeplot)
859	# ... the if loop is necessary so that we can loop above on the dp=0 case
860	if dp in [1,2]:
861	# and define what to do in plobj
862	plobj.invert = obj.invert_axes
863	plobj.swap = obj.swap_axes
864	# axis labels
865	plobj.xlabel = obj.absclab ; plobj.ylabel = obj.ordilab
866	# superpose or not (this is mostly for saving purpose)
867	plobj.superpose = self.superpose
868	# get title, colormaps, labels, etc.. from var
869	plobj.var = obj.var
870	plobj.define_from_var()
871	# generic 1D/2D: load field and coord in plot object
872	plobj.f = obj.field # field to be plotted
873	plobj.x = obj.absc # abscissa (or longitude)
874	plobj.y = obj.ordi # ordinate (or latitude)
875	# -- useless in 1D but not used anyway
876	# specific 1D plot stuff
877	if dp == 1:
878	# -- a default legend
879	plobj.legend = ""
880	if self.nfin > 1: plobj.legend = plobj.legend + " file #"+str(i+1)
881	if self.nvin > 1: plobj.legend = plobj.legend + " var "+plobj.var
882	if self.nplott > 1: plobj.legend = plobj.legend + " t="+str(self.t[t])
883	if self.nplotz > 1: plobj.legend = plobj.legend + " z="+str(self.z[z])
884	if self.nploty > 1: plobj.legend = plobj.legend + " y="+str(self.y[y])
885	if self.nplotx > 1: plobj.legend = plobj.legend + " x="+str(self.x[x])
886	# specific 2d plot stuff
887	if dp == 2:
888	# -- light grey background for missing values
889	if type(plobj.f).__name__ in 'MaskedArray': plobj.axisbg = '0.75'
890	# -- define if it is a map or a plot
891	plobj.mapmode = ( obj.method_x+obj.method_y == "freefree" \
892	and "grid points" not in obj.name_x \
893	and not self.noproj )
894	# possible user-defined plot settings shared by all plots
895	if self.div is not None: plobj.div = self.div
896	if self.xlabel is not None: plobj.xlabel = self.xlabel
897	if self.xcoeff is not None: plobj.xcoeff = self.xcoeff
898	if self.ylabel is not None: plobj.ylabel = self.ylabel
899	if self.ycoeff is not None: plobj.ycoeff = self.ycoeff
900	if self.title is not None: plobj.title = self.title
901	if self.units is not None: plobj.units = self.units
902	if self.colorbar is not None: plobj.colorbar = self.colorbar
903	if self.modx is not None: plobj.modx = self.modx
904	if self.nxticks is not None: plobj.nxticks = self.nxticks
905	if self.nyticks is not None: plobj.nyticks = self.nyticks
906	if self.fmt is not None: plobj.fmt = self.fmt
907	if self.showcb is not None: plobj.showcb = self.showcb
908	if self.logy is not None: plobj.logy = self.logy
909	if self.xmin is not None: plobj.xmin = self.xmin
910	if self.xmax is not None: plobj.xmax = self.xmax
911	if self.ymin is not None: plobj.ymin = self.ymin
912	if self.ymax is not None: plobj.ymax = self.ymax
913	# -- 1D specific
914	if dp == 1:
915	if self.linestyle is not None: plobj.linestyle = self.linestyle
916	if self.marker is not None: plobj.marker = self.marker
917	if self.color is not None: plobj.color = self.color
918	if self.legend is not None: plobj.legend = self.legend
919	# -- 2D specific
920	elif dp == 2:
921	if self.proj is not None and not self.noproj: plobj.proj = self.proj
922	if self.vmin is not None: plobj.vmin = self.vmin
923	if self.vmax is not None: plobj.vmax = self.vmax
924	if self.trans is not None: plobj.trans = self.trans
925	if self.back is not None: plobj.back = self.back
926	plobj.svx = self.svx
927	plobj.svy = self.svy
928	# finally append plot object
929	self.p.append(plobj)
930	count = count + 1
931	# self.nplot is number of plot to be defined in this call to defineplot()
932	# (because of self.plotin this might less than length of self.p)
933	self.nplot = count
934	# --- superimposed contours and vectors ---
935	# we have to start another loop because we need forward information
936	# TBD: there is probably a more flexible way to do that
937	count = 0
938	for i in range(self.nfin):
939	for j in range(self.nvin):
940	for t in range(self.nplott):
941	for z in range(self.nplotz):
942	for y in range(self.nploty):
943	for x in range(self.nplotx):
944	goal = self.vargoal[j] + self.filegoal[i]
945	obj = self.request[i][j][t][z][y][x]
946	if "mainmain" in goal and obj.dimplot == 2:
947	# the plot object we consider in the loop
948	plobj = self.p[count]
949	# -- see if there is a contour requested...
950	# (we use try because we might be at the end of the list)
951	found = 0
952	try: condvar = self.vargoal[j+1]
953	except: condvar = "itisok"
954	try: condfile = self.filegoal[i+1]
955	except: condfile = "itisok"
956	# ... get contour
957	##########################################
958	# NB: contour is expected to be right after main otherwise it is not displayed
959	##########################################
960	if condvar == "contour":
961	plobj.c = self.request[i][j+1][t][z][y][x].field ; found += 1
962	if condfile == "contour":
963	plobj.c = self.request[i+1][j][t][z][y][x].field ; found += 1
964	# see if there is a vector requested...
965	# (we use try because we might be at the end of the list)
966	try: condvar = self.vargoal[j+found+1]+self.vargoal[j+found+2]
967	except: condvar = "itisok"
968	try: condfile = self.filegoal[i+found+1]+self.filegoal[i+found+2]
969	except: condfile = "itisok"
970	# ... get vector and go directly to the next iteration
971	# (in some cases we would do this twice but this is cheap)
972	if "vector" in condvar:
973	plobj.vx = self.request[i][j+found+1][t][z][y][x].field
974	plobj.vy = self.request[i][j+found+2][t][z][y][x].field
975	if "vector" in condfile:
976	plobj.vx = self.request[i+found+1][j][t][z][y][x].field
977	plobj.vy = self.request[i+found+2][j][t][z][y][x].field
978	count = count + 1
979	# COUNT PLOTS. if 0 just exit.
980	# self.howmanyplots is self.nplot + possible extraplots
981	self.howmanyplots = self.nplot + extraplot
982	if self.howmanyplots > 0:
983	if self.verbose: print "**** OK. expect %i plots" % (self.howmanyplots)
984	else:
985	pass # because this means that we only had 0D values !
986	# final status
987	self.status = "definedplot"
988	return self
989
990	##############################################################################################
991	# makeplot method
992	# --> after defineplot and before makeplot, user-defined plot settings can be easily given
993	# simply by modifying the attributes of each elements of self.p
994	# --> to change only one plot setting, no need to call defineplot again before makeplot
995	# --> in the end, the array self.p of plot objects is saved for easy and convenient replotting
996	# --> see practical examples in the folder 'examples'
997	##############################################################################################
998	def makeplot(self):
999	if self.howmanyplots > 0:
1000	self.printstatus()
1001	# a few initial operations
1002	# ------------------------
1003	if "definedplot" not in self.status:
1004	print "!! ERROR !! Please use .defineplot() before .makeplot() can be used with your pp object." ; exit()
1005	# open a figure and define subplots
1006	# ---------------------------------
1007	if self.plotin is None:
1008	# start from scratch
1009	self.fig = ppplot.figuref(x=self.xp,y=self.yp)
1010	self.subv,self.subh = ppplot.definesubplot(self.howmanyplots,self.fig)
1011	self.n = 0
1012	## adapted space for labels etc
1013	## ... except for ortho because there is no label anyway
1014	self.customplot = self.p[0].f.ndim == 2 \
1015	and self.p[0].mapmode == True \
1016	and self.p[0].proj not in ["ortho"]
1017	if self.customplot:
1018	margin = 0.07
1019	self.fig.subplots_adjust(left=margin,right=1-margin,bottom=margin,top=1-margin)
1020	else:
1021	# start from an existing figure.
1022	# extraplot must have been set in the call to the previous figure.
1023	self.fig = self.plotin.fig
1024	self.subv,self.subh = self.plotin.subv,self.plotin.subh
1025	self.n = self.plotin.n
1026	self.howmanyplots = self.plotin.howmanyplots
1027	self.customplot = self.plotin.customplot
1028	# LOOP on all subplots
1029	# NB: cannot use 'for pl in self.p' if self.plotin not None
1030	# --------------------
1031	for count in range(self.nplot):
1032	# the plot object we consider in the loop
1033	pl = self.p[self.n]
1034	# before making the plot, create a subplot. the first one is numbered 1 not 0.
1035	# ... if pl.superpose, we use only one and only figure
1036	# ... (and we have to be careful with not doing things several times)
1037	if pl.superpose:
1038	if self.n == 0:
1039	self.fig.add_subplot(1,1,1,axisbg=pl.axisbg) # define one subplot (still needed for user-defined font sizes)
1040	sav = pl.xlabel,pl.ylabel,\
1041	pl.xcoeff,pl.ycoeff,\
1042	pl.nxticks,pl.nyticks,\
1043	pl.xmin,pl.xmax,\
1044	pl.ymin,pl.ymax,\
1045	pl.title,pl.swaplab # save titles and labels
1046	# possibility to color lines according to a color map
1047	# ... made so that all plots span the whole color map automatically.
1048	if self.colorbar is not None:
1049	if self.verbose: print "**** OK. We make a rainbow spaghetti plot with color map ",self.colorbar
1050	ppplot.rainbow(cb=self.colorbar,num=self.howmanyplots)
1051	else:
1052	pl.invert = False ; pl.linestyle = None # don't invert again axis
1053	# set saved titles and labels
1054	if self.plotin is None:
1055	pl.xlabel,pl.ylabel,\
1056	pl.xcoeff,pl.ycoeff,\
1057	pl.nxticks,pl.nyticks,\
1058	pl.xmin,pl.xmax,\
1059	pl.ymin,pl.ymax,\
1060	pl.title,pl.swaplab = sav
1061	else:
1062	prev_plot = self.plotin.p[self.n-1]
1063	pl.xlabel = prev_plot.xlabel ; pl.ylabel = prev_plot.ylabel
1064	pl.xcoeff = prev_plot.xcoeff ; pl.ycoeff = prev_plot.ycoeff
1065	pl.nxticks = prev_plot.nxticks ; pl.nyticks = prev_plot.nyticks
1066	pl.xmin = prev_plot.xmin ; pl.xmax = prev_plot.xmax
1067	pl.ymin = prev_plot.ymin ; pl.ymax = prev_plot.ymax
1068	pl.title = prev_plot.title ; pl.swaplab = prev_plot.swaplab
1069	else:
1070	self.fig.add_subplot(self.subv,self.subh,self.n+1,axisbg=pl.axisbg)
1071	if self.verbose: print "**** Done subplot %i / %i " %( self.n+1,self.howmanyplots )
1072	# finally make the plot
1073	pl.make()
1074	# possibly print results in a text file
1075	if self.savtxt:
1076	if self.verbose: print "**** Printing results in a text file"
1077	name = pl.var + "%04d" % self.n
1078	ppplot.writeascii(field=pl.f,absc=pl.x,name=name)
1079	# increment plot count (and propagate this in plotin)
1080	self.n = self.n+1
1081	if self.plotin is not None: self.plotin.n = self.n
1082	# once completed show the plot (cannot show intermediate plotin)
1083	# ... added a fix (customplot=True) for the label problem in basemap
1084	if not self.quiet: print "**** PPCLASS. Done step: makeplot"
1085	if (self.n == self.howmanyplots):
1086	ppplot.save(mode=self.out,filename=self.filename,folder=self.folder,custom=self.customplot,includedate=self.includedate,res=self.res)
1087	mpl.close()
1088	# SAVE A PICKLE FILE WITH THE self.p ARRAY OF OBJECTS
1089	if self.verbose: print "**** Saving session in "+self.filename + ".ppobj"
1090	savfile = self.folder + "/" + self.filename + ".ppobj"
1091	try:
1092	filehandler = open(savfile, 'w')
1093	pickle.dump(self.p, filehandler)
1094	except IOError:
1095	if self.verbose: print "!! WARNING !! Saved object file not written. Probably do not have permission to write here."
1096	return self
1097
1098	###########################################################
1099	# plot: a shortcut method for the defineplot + plot chain #
1100	###########################################################
1101	def plot(self,extraplot=0):
1102	self.defineplot(extraplot=extraplot)
1103	self.makeplot()
1104	return self
1105
1106	#######################################################
1107	# getplot: a shortcut method for the get + plot chain #
1108	#######################################################
1109	def getplot(self,extraplot=0):
1110	self.get()
1111	self.plot(extraplot=extraplot)
1112	return self
1113
1114	###################################################################
1115	# getdefineplot: a shortcut method for the get + defineplot chain #
1116	###################################################################
1117	def getdefineplot(self,extraplot=0):
1118	self.get()
1119	self.defineplot(extraplot=extraplot)
1120	return self
1121
1122	#################################################################
1123	# func: operation on two pp objects being on status 'definedplot'
1124	# this allows for one field being function of another one
1125	# e.g. u.func(v) means u will be displayed as a function of v
1126	# ... no need to do defineplot after u.func(v), makeplot directly
1127	#################################################################
1128	def func(self,other):
1129	# preamble: for this operation to work, defineplot() must have been done
1130	if self.status != "definedplot":
1131	if self.verbose: print "!! WARNING !! performing defineplot on operand"
1132	self.defineplot()
1133	if other.status != "definedplot":
1134	if self.verbose: print "!! WARNING !! performing defineplot on operand"
1135	other.defineplot()
1136	# check total number of plots
1137	if self.howmanyplots != other.howmanyplots:
1138	print "!! ERROR !! The two operands do not have the same number of subplots."
1139	exit()
1140	# and now operation.
1141	count = 0
1142	while count < self.howmanyplots:
1143	sobj = self.p[count] ; oobj = other.p[count]
1144	if sobj.f.ndim !=1 or oobj.f.ndim !=1:
1145	if self.verbose: print "!! WARNING !! Flattening arrays because more than one-dimensional."
1146	sobj.f = np.ravel(sobj.f)
1147	oobj.f = np.ravel(oobj.f)
1148	sobj.x = oobj.f
1149	sobj.xlabel = oobj.ylabel
1150	if sobj.x.size > sobj.f.size:
1151	if self.verbose:
1152	print "!! WARNING !! Trying to define y=f(x) with x and y not at the same size.",sobj.x.size,sobj.f.size
1153	print "!! WARNING !! Modifying x to fit y size but please check."
1154	sobj.x = sobj.x[0:sobj.f.size]
1155	count = count + 1
1156	return self
1157
1158	###########################################################
1159	# copyopt: get options from e.g. a parser
1160	# ... allow for simple scripting and user-defined settings
1161	# ... must be called between defineplot and makeplot
1162	# REQUIRED: attributes of opt must be the same as in the pp object
1163	###########################################################
1164	def getopt(self,opt):
1165	# -- if only one, or less than the number of plots --> we take the first one
1166	# -- if as many as number of plots --> OK, each plot has its own setting
1167	# (except a few cases such as trans)
1168	for iii in range(self.howmanyplots):
1169	## solve the bug about reversed labels with swaplab
1170	if opt.xlabel is None and opt.ylabel is None:
1171	self.p[iii].swaplab = True
1172	else:
1173	self.p[iii].swaplab = False
1174	##
1175	if opt.void:
1176	self.p[iii].showcb = False
1177	else:
1178	self.p[iii].showcb = True
1179	###
1180	try: self.p[iii].trans = opt.trans
1181	except: pass
1182	###
1183	try: self.p[iii].div = opt.div
1184	except: pass
1185	###
1186	try: self.p[iii].logy = opt.logy
1187	except: pass
1188	###
1189	try: self.p[iii].colorbar = opt.colorbar[iii]
1190	except:
1191	try: self.p[iii].colorbar = opt.colorbar[0] ; self.colorbar = opt.colorbar[0]
1192	except: pass
1193	###
1194	if opt.void:
1195	self.p[iii].title = ""
1196	else:
1197	try: self.p[iii].title = opt.title[iii]
1198	except:
1199	try: self.p[iii].title = opt.title[0]
1200	except: pass
1201	###
1202	if opt.void:
1203	self.p[iii].xlabel = ""
1204	else:
1205	try: self.p[iii].xlabel = opt.xlabel[iii]
1206	except:
1207	try: self.p[iii].xlabel = opt.xlabel[0]
1208	except: pass
1209	###
1210	if opt.void:
1211	self.p[iii].ylabel = ""
1212	else:
1213	try: self.p[iii].ylabel = opt.ylabel[iii]
1214	except:
1215	try: self.p[iii].ylabel = opt.ylabel[0]
1216	except: pass
1217	###
1218	try: self.p[iii].linestyle = opt.linestyle[iii]
1219	except:
1220	try: self.p[iii].linestyle = opt.linestyle[0]
1221	except: pass
1222	###
1223	try: self.p[iii].color = opt.color[iii]
1224	except:
1225	try: self.p[iii].color = opt.color[0]
1226	except: pass
1227	###
1228	try: self.p[iii].marker = opt.marker[iii]
1229	except:
1230	try: self.p[iii].marker = opt.marker[0]
1231	except: pass
1232	###
1233	try: self.p[iii].legend = opt.legend[iii]
1234	except:
1235	try: self.p[iii].legend = opt.legend[0]
1236	except: pass
1237	###
1238	try: self.p[iii].proj = opt.proj[iii]
1239	except:
1240	try: self.p[iii].proj = opt.proj[0]
1241	except: pass
1242	###
1243	try: self.p[iii].back = opt.back[iii]
1244	except:
1245	try: self.p[iii].back = opt.back[0]
1246	except: pass
1247	###
1248	try: self.p[iii].area = opt.area[iii]
1249	except:
1250	try: self.p[iii].area = opt.area[0]
1251	except: pass
1252	###
1253	try: self.p[iii].blon = opt.blon[iii]
1254	except:
1255	try: self.p[iii].blon = opt.blon[0]
1256	except: pass
1257	###
1258	try: self.p[iii].blat = opt.blat[iii]
1259	except:
1260	try: self.p[iii].blat = opt.blat[0]
1261	except: pass
1262	###
1263	try: self.p[iii].vmin = opt.vmin[iii]
1264	except:
1265	try: self.p[iii].vmin = opt.vmin[0]
1266	except: pass
1267	###
1268	try: self.p[iii].vmax = opt.vmax[iii]
1269	except:
1270	try: self.p[iii].vmax = opt.vmax[0]
1271	except: pass
1272	###
1273	try: self.p[iii].xcoeff = opt.xcoeff[iii]
1274	except:
1275	try: self.p[iii].xcoeff = opt.xcoeff[0]
1276	except: pass
1277	###
1278	try: self.p[iii].ycoeff = opt.ycoeff[iii]
1279	except:
1280	try: self.p[iii].ycoeff = opt.ycoeff[0]
1281	except: pass
1282	###
1283	try: self.p[iii].units = opt.units[iii]
1284	except:
1285	try: self.p[iii].units = opt.units[0]
1286	except: pass
1287	###
1288	try: self.p[iii].wscale = opt.wscale[iii]
1289	except:
1290	try: self.p[iii].wscale = opt.wscale[0]
1291	except: pass
1292	###
1293	try: self.p[iii].xmin = opt.xmin[iii]
1294	except:
1295	try: self.p[iii].xmin = opt.xmin[0]
1296	except: pass
1297	###
1298	try: self.p[iii].ymin = opt.ymin[iii]
1299	except:
1300	try: self.p[iii].ymin = opt.ymin[0]
1301	except: pass
1302	###
1303	try: self.p[iii].xmax = opt.xmax[iii]
1304	except:
1305	try: self.p[iii].xmax = opt.xmax[0]
1306	except: pass
1307	###
1308	try: self.p[iii].ymax = opt.ymax[iii]
1309	except:
1310	try: self.p[iii].ymax = opt.ymax[0]
1311	except: pass
1312	###
1313	try: self.p[iii].nxticks = opt.nxticks[iii]
1314	except:
1315	try: self.p[iii].nxticks = opt.nxticks[0]
1316	except: pass
1317	###
1318	try: self.p[iii].nyticks = opt.nyticks[iii]
1319	except:
1320	try: self.p[iii].nyticks = opt.nyticks[0]
1321	except: pass
1322	###
1323	try: self.p[iii].cbticks = opt.cbticks[iii]
1324	except:
1325	try: self.p[iii].cbticks = opt.cbticks[0]
1326	except: pass
1327	###
1328	try: self.p[iii].modx = opt.modx[iii]
1329	except:
1330	try: self.p[iii].modx = opt.modx[0]
1331	except: pass
1332	###
1333	try: self.p[iii].fmt = opt.fmt[iii]
1334	except:
1335	try: self.p[iii].fmt = opt.fmt[0]
1336	except: pass
1337
1338
1339	##########################################################
1340	### THE ONEREQUEST SUBOBJECT TO PP (ON WHICH IT LOOPS) ###
1341	##########################################################
1342	class onerequest():
1343
1344	# default settings. mostly initialized to diagnose problem, except dimplot, nplot, verbose, swap_axes, invert_axes
1345	# -------------------------------
1346	def __init__(self):
1347	self.file = '!! file: I am not set, damned !!'
1348	self.f = None
1349	self.dim = None
1350	self.var = '!! var: I am not set, damned !!'
1351	self.index_x = [] ; self.index_y = [] ; self.index_z = [] ; self.index_t = []
1352	self.index_x2d = [] ; self.index_y2d = []
1353	self.method_x = '!! method_x: I am not set, damned !!'
1354	self.method_y = '!! method_y: I am not set, damned !!'
1355	self.method_z = '!! method_z: I am not set, damned !!'
1356	self.method_t = '!! method_t: I am not set, damned !!'
1357	self.field = None
1358	self.name_x = None ; self.name_y = None ; self.name_z = None ; self.name_t = None
1359	self.dim_x = None ; self.dim_y = None ; self.dim_z = None ; self.dim_t = None
1360	self.field_x = None ; self.field_y = None ; self.field_z = None ; self.field_t = None
1361	self.tabtime = None
1362	self.dimplot = 0
1363	self.nplot = 1
1364	self.absc = None ; self.ordi = None ; self.absclab = None ; self.ordilab = None
1365	self.verbose = True
1366	self.swap_axes = False ; self.invert_axes = False
1367	self.compute = None
1368	self.changetime = None
1369	self.sx = 1 ; self.sy = 1 ; self.sz = 1 ; self.st = 1
1370	self.missing = '!! missing value: I am not set, damned !!'
1371	self.kind3d = '!! kind3d: I am not set, damned !!'
1372
1373	# open a file. for now it is netcdf. TBD for other formats.
1374	# check that self.var is inside.
1375	# -------------------------------
1376	def openfile(self):
1377	if not os.path.exists(self.file): print '!! ERROR !! I could not find the following file: '+self.file ; exit()
1378	if not os.path.isfile(self.file): print '!! ERROR !! This does not appear to be a file: '+self.file ; exit()
1379	self.f = netCDF4.Dataset(self.file)
1380	if self.verbose: print "**** OK. Opened file "+self.file
1381	if self.var not in self.f.variables.keys():
1382	print '!! ERROR !! File '+self.file+' does not contain variable: '+self.var
1383	print '..... try instead with ',self.f.variables.keys() ; exit()
1384
1385	# close a file
1386	# ------------
1387	def closefile(self):
1388	self.f.close()
1389
1390	# copy attributes from another existing object
1391	# --------------------------------------------
1392	def copy(self,source):
1393	for k, v in vars(source).items():
1394	setattr(self,k,v)
1395
1396	# get x,y,z,t dimensions from NETCDF file
1397	# TBD: user could request for a specific altitude dimension
1398	# TBD: staggered variables could request specific dimensions
1399	# -------------------------------
1400	def getdim(self):
1401	# GET SIZES OF EACH DIMENSION
1402	if self.verbose: print "**** OK. Found variable "+self.var
1403	shape = self.f.variables[self.var].shape
1404	self.dim = len(shape)
1405	if self.dim == 1:
1406	if self.verbose: print "**** OK. 1D field. I assume this varies with time."
1407	self.dim_x = 1 ; self.dim_y = 1 ; self.dim_z = 1 ; self.dim_t = shape[0]
1408	elif self.dim == 2:
1409	if self.verbose: print "**** OK. 2D field. I assume this is not-time-varying lat-lon map."
1410	self.dim_x = shape[1] ; self.dim_y = shape[0] ; self.dim_z = 1 ; self.dim_t = 1
1411	elif self.dim == 3:
1412	if self.verbose: print "**** OK. 3D field. I assume this is time-varying lat-lon map."
1413	## see below for comment
1414	if self.kind3d == "tyx":
1415	self.dim_x = shape[2] ; self.dim_y = shape[1] ; self.dim_z = 1 ; self.dim_t = shape[0]
1416	elif self.kind3d == "tzy":
1417	self.dim_x = 1 ; self.dim_y = shape[2] ; self.dim_z = shape[1] ; self.dim_t = shape[0]
1418	else:
1419	print "!! ERROR !! This kind of 3D field is not supported. Please send feedback."
1420	print self.kind3d
1421	exit()
1422	elif self.dim == 4:
1423	if self.verbose: print "**** OK. 4D field."
1424	self.dim_x = shape[3] ; self.dim_y = shape[2] ; self.dim_z = shape[1] ; self.dim_t = shape[0]
1425	# LONGITUDE. Try preset fields. If not present set grid points axis.
1426	self.name_x = "nothing"
1427	for c in glob_listx:
1428	if c in self.f.variables.keys():
1429	self.name_x = c
1430	if self.name_x == "nothing":
1431	self.field_x = np.array(range(self.dim_x))
1432	self.name_x = "x grid points"
1433	else:
1434	self.field_x = self.f.variables[self.name_x]
1435	# LATITUDE. Try preset fields. If not present set grid points axis.
1436	self.name_y = "nothing"
1437	for c in glob_listy:
1438	if c in self.f.variables.keys():
1439	self.name_y = c
1440	if self.name_y == "nothing":
1441	self.field_y = np.array(range(self.dim_y))
1442	self.name_y = "y grid points"
1443	else:
1444	self.field_y = self.f.variables[self.name_y]
1445	# ensure that lon and lat are 2D fields
1446	# 1. simple 1D case (not time-varying)
1447	if len(self.field_x.shape)*len(self.field_y.shape) == 1:
1448	if self.verbose: print "**** OK. recasting lon and lat as 2D fields."
1449	[self.field_x,self.field_y] = np.meshgrid(self.field_x,self.field_y)
1450	# 2. complex 3D case (time-varying, actually just copied over time axis)
1451	elif len(self.field_x.shape)*len(self.field_y.shape) == 9:
1452	if self.verbose: print "**** OK. reducing lon and lat as 2D fields. get rid of time."
1453	self.field_x = self.field_x[0,:,:]
1454	self.field_y = self.field_y[0,:,:]
1455	# if xy axis are apparently undefined, set 2D grid points axis.
1456	if "grid points" not in self.name_x:
1457	if np.all(self.field_x == self.field_x[0,0]) \
1458	or self.field_x.min() == self.field_x.max() \
1459	or self.field_y.min() == self.field_y.max():
1460	if self.verbose: print "!! WARNING !! xy axis look undefined. creating non-dummy ones."
1461	self.field_x = np.array(range(self.dim_x)) ; self.name_x = "x grid points"
1462	self.field_y = np.array(range(self.dim_y)) ; self.name_y = "y grid points"
1463	[self.field_x,self.field_y] = np.meshgrid(self.field_x,self.field_y)
1464	if self.dim_x > 1:
1465	if self.verbose: print "**** OK. x axis %4.0f values [%5.1f,%5.1f]" % (self.dim_x,self.field_x.min(),self.field_x.max())
1466	if self.dim_y > 1:
1467	if self.verbose: print "**** OK. y axis %4.0f values [%5.1f,%5.1f]" % (self.dim_y,self.field_y.min(),self.field_y.max())
1468	# ALTITUDE. Try preset fields. If not present set grid points axis.
1469	# WARNING: how do we do if several are available? the last one is chosen.
1470	self.name_z = "nothing"
1471	for c in glob_listz:
1472	if c in self.f.variables.keys():
1473	self.name_z = c
1474	if self.name_z == "nothing":
1475	self.field_z = np.array(range(self.dim_z))
1476	self.name_z = "z grid points"
1477	else:
1478	tabalt = self.f.variables[self.name_z]
1479	# (consider the case where tabtime is not dim 1) TBD: 2D and 3D cases
1480	if tabalt.ndim == 4:
1481	try:
1482	self.field_z = tabalt[1,:,0,0] # 4D case. alt is usually third dimension.
1483	# 1 for time to avoid initial 0s
1484	except:
1485	self.field_z = tabalt[0,:,0,0]
1486	if self.verbose: print "!! WARNING !! "+self.name_z+" is 4D var. We made it 1D."
1487	else:
1488	self.field_z = self.f.variables[self.name_z][:] # specify dimension
1489	# TBD: problems when self.dim_z != self.field_z.size
1490	if self.field_z.size != self.dim_z:
1491	if self.verbose: print "!! WARNING !! Cannot use this z coordinate. Not enough points. Use simple z axis."
1492	self.field_z = np.array(range(self.dim_z))
1493	self.name_z = "z grid points"
1494	if self.dim_z > 1:
1495	if self.verbose: print "**** OK. z axis %4.0f values [%5.1f,%5.1f]" % (self.dim_z,self.field_z.min(),self.field_z.max())
1496
1497	# TIME. Try preset fields.
1498	self.name_t = "nothing"
1499	for c in glob_listt:
1500	if c in self.f.variables.keys():
1501	self.name_t = c
1502	if self.verbose: print "**** OK. Found time variable: ",c
1503	try:
1504	# speed up: only get first value, last one.
1505	self.tabtime = self.f.variables[self.name_t]
1506	# (consider the case where tabtime is not dim 1)
1507	# (time is most often the first dimension)
1508	if self.tabtime.ndim == 2: self.tabtime = self.tabtime[:,0]
1509	elif self.tabtime.ndim == 3: self.tabtime = self.tabtime[:,0,0]
1510	elif self.tabtime.ndim == 4: self.tabtime = self.tabtime[:,0,0,0]
1511	# (now proceed) (the +0. just ensures this is a number)
1512	dafirst = self.tabtime[0] + 0.
1513	if self.dim_t == 1:
1514	self.field_t = np.array([dafirst])
1515	else:
1516	daint = self.tabtime[1] - dafirst
1517	dalast = dafirst + (self.dim_t-1)*daint
1518	self.field_t = np.linspace(dafirst,dalast,num=self.dim_t)
1519	if self.verbose:
1520	print "!! WARNING !! WARNING !! Time axis is supposed to be equally spaced !!"
1521	if dalast != self.tabtime[self.dim_t-1]:
1522	print "!! WARNING !! Time axis has been recast to be monotonic",dalast,self.tabtime[self.dim_t-1]
1523	except:
1524	# ... or if a problem encountered, define a simple time axis
1525	if self.verbose: print "**** OK. There is something weird. Let us go for a simple time axis."
1526	self.field_t = np.array(range(self.dim_t))
1527	self.name_t = "t grid points"
1528	if self.dim_t > 1:
1529	if self.verbose: print "**** OK. t axis %4.0f values [%5.1f,%5.1f]" % (self.dim_t,self.field_t.min(),self.field_t.max())
1530
1531	# change time axis
1532	# ... add your options here!
1533	# --------------------------
1534	def performtimechange(self):
1535	if self.changetime is not None:
1536	if self.verbose: print "**** OK. Converting time axis:",self.changetime
1537	### options added by T. Navarro
1538	if self.changetime == "mars_sol2ls":
1539	if "controle" in self.f.variables:
1540	self.field_t = self.field_t \
1541	+ self.f.variables['controle'][3]%669 \
1542	+ self.f.variables['controle'][26]
1543	self.field_t = ppcompute.mars_sol2ls(self.field_t)
1544	elif self.changetime == "mars_dayini" and "controle" in self.f.variables:
1545	self.field_t = self.field_t \
1546	+ self.f.variables['controle'][3]%669 \
1547	+ self.f.variables['controle'][26]
1548	### options added by A. Spiga
1549	elif self.changetime == "correctls":
1550	if self.tabtime is None:
1551	print "!! WARNING !! Encountered a problem with correctls. Check your time dimension. Skipping this part."
1552	else:
1553	dafirst = self.tabtime[0] + 0.
1554	if self.dim_t == 1:
1555	self.field_t = np.array([dafirst])
1556	else:
1557	daint = self.tabtime[1] - dafirst
1558	dalast = dafirst + (self.dim_t-1)*daint
1559	year = 0.
1560	add = np.linspace(dafirst,dalast,num=self.dim_t) ; add[0] = 0.
1561	for iii in range(1,self.dim_t):
1562	if self.tabtime[iii] - self.tabtime[iii-1] < 0: year = year+1.
1563	add[iii] = year*360.
1564	self.field_t = add + self.tabtime
1565	elif "mars_meso" in self.changetime:
1566	if 'Times' not in self.f.variables.keys():
1567	if self.verbose: print "!! WARNING !! Variable Times not in file. Cannot proceed to change of time axis."
1568	else:
1569	# get the array of strings describing dates
1570	dates = self.f.variables['Times']
1571	dates.set_auto_maskandscale(False) # necessary to solve the api Times bug!
1572	# get ls sol utc from those strings
1573	ls, sol, utc = ppcompute.mars_date(dates[:])
1574	# populate self.field_t with the right output from mars_date
1575	if self.changetime == "mars_meso_ls":
1576	self.field_t = ls
1577	self.name_t = "Ls"
1578	elif self.changetime == "mars_meso_sol":
1579	self.field_t = sol
1580	self.name_t = "sol"
1581	elif self.changetime == "mars_meso_utc" \
1582	and ( self.changetime == "mars_meso_lt" \
1583	and not hasattr(self.f,'CEN_LON') ):
1584	self.field_t = ppcompute.timecorrect(utc)
1585	self.name_t = "utc"
1586	if self.method_t == "fixed":
1587	self.field_t = self.field_t % 24 # so that the user is not mistaken!
1588	elif self.changetime == "mars_meso_lt":
1589	self.field_t = ppcompute.timecorrect(utc) + getattr(self.f,'CEN_LON') / 15.
1590	self.field_t = ppcompute.timecorrect(self.field_t)
1591	self.name_t = "local time (center of domain)"
1592	if self.method_t == "fixed":
1593	self.field_t = self.field_t % 24 # so that the user is not mistaken!
1594	else:
1595	if self.verbose: print "!! WARNING !! This time change is not implemented. Nothing is done."
1596	if self.verbose: print "**** OK. new t axis values [%5.1f,%5.1f]" % (self.field_t.min(),self.field_t.max())
1597
1598	# get list of index to be retrieved for time axis
1599	### TBD: il faudrait ne prendre que les indices qui correspondent a l interieur d un plot (dans all)
1600	# -------------------------------
1601	def getindextime(self,dalist=None,ind=None):
1602	if self.method_t == "free":
1603	self.index_t = np.arange(0,self.dim_t,self.st)
1604	if self.dim_t > 1:
1605	self.dimplot = self.dimplot + 1
1606	if self.verbose: print "**** OK. t values. all."
1607	else:
1608	self.method_t = "fixed"
1609	if self.verbose: print "**** OK. no t dimension."
1610	elif self.method_t == "comp":
1611	start = np.argmin( np.abs( self.field_t - dalist[ind][0] ) )
1612	stop = np.argmin( np.abs( self.field_t - dalist[ind][1] ) )
1613	self.index_t = np.arange(start,stop+1,self.st)
1614	if self.verbose: print "**** OK. t values. comp over interval ",self.field_t[start],self.field_t[stop]," nvalues=",self.index_t.size
1615	elif self.method_t == "fixed":
1616	self.index_t.append( np.argmin( np.abs( self.field_t - dalist[ind][0] ) ))
1617	if self.verbose: print "**** OK. t values",self.field_t[self.index_t]
1618	else:
1619	print "!! ERROR !! method "+self.method_t+" not supported"
1620	self.index_t = np.array(self.index_t)
1621
1622	# get list of index to be retrieved for vertical axis
1623	### TBD: il faudrait ne prendre que les indices qui correspondent a l interieur d un plot (dans all)
1624	# -------------------------------
1625	def getindexvert(self,dalist=None,ind=None):
1626	if self.method_z == "free":
1627	self.index_z = np.arange(0,self.dim_z,self.sz)
1628	if self.dim_z > 1:
1629	self.dimplot = self.dimplot + 1
1630	if self.verbose: print "**** OK. z values. all."
1631	else:
1632	self.method_z = "fixed"
1633	if self.verbose: print "**** OK. no z dimension."
1634	elif self.method_z == "comp":
1635	start = np.argmin( np.abs( self.field_z - dalist[ind][0] ) )
1636	stop = np.argmin( np.abs( self.field_z - dalist[ind][1] ) )
1637	self.index_z = np.arange(start,stop+1,self.sz)
1638	if self.verbose: print "**** OK. z values. comp over interval",self.field_z[start],self.field_z[stop]," nvalues=",self.index_z.size
1639	elif self.method_z == "fixed":
1640	self.index_z.append( np.argmin( np.abs( self.field_z - dalist[ind][0] ) ))
1641	if self.verbose: print "**** OK. z values",self.field_z[self.index_z]
1642	else:
1643	if self.verbose: print "!! ERROR !! method "+self.method_z+" not supported"
1644	self.index_z = np.array(self.index_z)
1645
1646	# get list of index to be retrieved for horizontal grid
1647	# --> index_x and index_y are slices to be retrieved from NETCDF files
1648	# --> index_x2d and index_y2d are the actual (x,y) coordinates corresponding to each relevant point
1649	# [this is slightly more complicated because 2D arrays for lat-lon projection possibly irregular]
1650	# NB: to append index we use lists (the most convenient) then we convert into a numpy.array
1651	### TBD: il faudrait ne prendre que les indices qui correspondent a l interieur d un plot (dans all)
1652	# -------------------------------
1653	def getindexhori(self,dalistx=None,dalisty=None,indx=None,indy=None):
1654	## get what is the method over x and y axis
1655	test = self.method_x+self.method_y
1656	## CASE 0, EASY CASES:
1657	## - LAT IS FREE (we do here what must be done whatever LON case is)
1658	## - LON IS FREE (we do here what must be done whatever LAT case is)
1659	## - LAT IS COMP AND LON IS FREE
1660	## - LON IS COMP AND LAT IS FREE
1661	if self.method_x == "free" or test in ["compfree","compcomp"]:
1662	self.index_x = range(0,self.dim_x,self.sx)
1663	if self.dim_x > 1:
1664	if self.method_x == "free": self.dimplot = self.dimplot + 1
1665	if self.verbose: print "**** OK. x values. all."
1666	else:
1667	self.method_x = "fixed"
1668	if self.verbose: print "**** OK. no x dimension."
1669	if self.method_y == "free" or test in ["freecomp","compcomp"]:
1670	self.index_y = range(0,self.dim_y,self.sy)
1671	if self.dim_y > 1:
1672	if self.method_y == "free": self.dimplot = self.dimplot + 1
1673	if self.verbose: print "**** OK. y values. all."
1674	else:
1675	self.method_y = "fixed"
1676	if self.verbose: print "**** OK. no y dimension."
1677	## CASE 0 above, this is just for continuity for free.
1678	## ... for comp we have to select bounds.
1679	## ... TBD: take the bool array strategy for what follows!
1680	if self.method_x in ["free","comp"] and self.method_y in ["free","comp"]:
1681	### ref1_dirty_hack
1682	### ... for the moment this is a hack. but actually this is more powerful.
1683	if self.method_x == "comp":
1684	yeah = (self.field_x >= dalistx[indx][0])*(self.field_x <= dalistx[indx][1])
1685	self.index_x = yeah[0,:]
1686	if self.method_y == "comp":
1687	yeah = (self.field_y >= dalisty[indy][0]) * (self.field_y <= dalisty[indy][1])
1688	self.index_y = yeah[:,0]
1689	self.index_x2d = self.index_x
1690	self.index_y2d = self.index_y
1691	## AND NOW THE LITTLE BIT MORE COMPLICATED CASES
1692	## CASE 1 LAT AND LON ARE FIXED
1693	elif test == "fixedfixed":
1694	idy,idx = np.unravel_index( np.argmin( ( self.field_x - dalistx[indx][0])2 + (self.field_y - dalisty[indy][0])2 ), self.field_x.shape )
1695	#TBD: pb with staggered coord
1696	if idx not in self.index_x: self.index_x.append(idx)
1697	if idy not in self.index_y: self.index_y.append(idy)
1698	self.index_x2d.append(idx)
1699	self.index_y2d.append(idy)
1700	## CASE 2 LON IS FIXED BUT NOT LAT
1701	elif test in ["fixedfree","fixedcomp"]:
1702	# find where are requested x values for each y on the free dimension
1703	# NB: this does not work for non-bijective cases e.g. polar stereographic
1704	for iy in range(self.dim_y):
1705	idx = np.argmin( np.abs( self.field_x[iy,:] - dalistx[indx][0] ) )
1706	# if comp is requested we select only indexes which yield values between requested min and max
1707	storeval = (self.method_y == "comp") and (self.field_y[iy,idx] > dalisty[indy][0]) and (self.field_y[iy,idx] < dalisty[indy][1])
1708	storeval = storeval or (self.method_y == "free")
1709	if storeval:
1710	if idx not in self.index_x: self.index_x.append(idx)
1711	if iy not in self.index_y and self.method_y == "comp": self.index_y.append(iy)
1712	if idx not in self.index_x2d or iy not in self.index_y2d:
1713	self.index_x2d.append(idx)
1714	self.index_y2d.append(iy)
1715	## CASE 3 LAT IS FIXED BUT NOT LON
1716	elif test in ["freefixed","compfixed"]:
1717	# find where are requested y values for each x on the free dimension
1718	# NB: this does not work for non-bijective cases e.g. polar stereographic
1719	for ix in range(self.dim_x):
1720	idy = np.argmin( np.abs( self.field_y[:,ix] - dalisty[indy][0] ) )
1721	# if comp is requested we select only indexes which yield values between requested min and max
1722	storeval = (self.method_x == "comp") and (self.field_x[idy,ix] > dalistx[indx][0]) and (self.field_x[idy,ix] < dalistx[indx][1])
1723	storeval = storeval or (self.method_x == "free")
1724	if storeval:
1725	if idy not in self.index_y: self.index_y.append(idy)
1726	if ix not in self.index_x and self.method_x == "comp": self.index_x.append(ix)
1727	if ix not in self.index_x2d or idy not in self.index_y2d:
1728	self.index_x2d.append(ix)
1729	self.index_y2d.append(idy)
1730	## check index tab
1731	if len(self.index_x) == 0 or len(self.index_y) == 0:
1732	print "!! ERROR !! no indices found. check prescribed latitudes or longitudes" ; exit()
1733	## ensure the array is a numpy array for getfield to work
1734	self.index_x = np.array(self.index_x)
1735	self.index_y = np.array(self.index_y)
1736	self.index_x2d = np.array(self.index_x2d)
1737	self.index_y2d = np.array(self.index_y2d)
1738	### print extrema
1739	printx = self.field_x[np.ix_(self.index_y2d, self.index_x2d)]
1740	printy = self.field_y[np.ix_(self.index_y2d, self.index_x2d)]
1741	if self.verbose:
1742	print "**** OK. x values (min,max).", printx.min(),printx.max()
1743	print "**** OK. y values (min,max).", printy.min(),printy.max()
1744
1745	# get the field from the NETCDF file and perform averages
1746	# -------------------------------
1747	def getfield(self):
1748	## first tell what is to be done
1749	if self.verbose:
1750	if self.dimplot > 2: print "**** !! WARNING !! "+str(self.dimplot)+"D plots will not be supported!"
1751	elif self.dimplot == 0 and self.verbose: print "**** OK. 0D value requested."
1752	elif self.dimplot == 1 and self.verbose: print "**** OK. 1D plot requested."
1753	elif self.verbose: print "**** OK. 2D section requested."
1754	# well, now get field from netcdf file
1755	# part below is necessary otherwise there is an index error below
1756	if self.index_x.size == 1: self.index_x = self.index_x[0]
1757	if self.index_y.size == 1: self.index_y = self.index_y[0]
1758	if self.index_z.size == 1: self.index_z = self.index_z[0]
1759	if self.index_t.size == 1: self.index_t = self.index_t[0]
1760	# then retrieve what is requested by user
1761	# each self.dim case corresponds to tests in the beginning of getdim.
1762	time0 = timelib.time()
1763	if self.verbose: print "**** OK. I am getting values from files. Please wait."
1764	if self.dim == 1:
1765	nt = self.index_t.size ; nz = 1 ; ny = 1 ; nx = 1
1766	self.field = self.f.variables[self.var][self.index_t]
1767	elif self.dim == 2:
1768	nt = 1 ; nz = 1 ; ny = self.index_y.size ; nx = self.index_x.size
1769	self.field = self.f.variables[self.var][self.index_y,self.index_x]
1770	elif self.dim == 3:
1771	## DEFAULT tyx (time-varying 2D field)
1772	if self.kind3d == "tyx":
1773	nt = self.index_t.size ; nz = 1 ; ny = self.index_y.size ; nx = self.index_x.size
1774	self.field = self.f.variables[self.var][self.index_t,self.index_y,self.index_x]
1775	## CASE tzy (e.g. time-varying zonal mean y-z field)
1776	elif self.kind3d == "tzy":
1777	nt = self.index_t.size ; nz = self.index_z.size ; ny = self.index_y.size ; nx = 1
1778	self.field = self.f.variables[self.var][self.index_t,self.index_z,self.index_y]
1779	else:
1780	print "!! ERROR !! This kind of 3D field is not supported. Please send feedback." ; exit()
1781	# this is far faster than retrieving each term with a loop
1782	elif self.dim == 4:
1783	nt = self.index_t.size ; nz = self.index_z.size ; ny = self.index_y.size ; nx = self.index_x.size
1784	self.field = self.f.variables[self.var][self.index_t,self.index_z,self.index_y,self.index_x]
1785	else:
1786	print "!! ERROR !! field would have more than four dimensions ?" ; exit()
1787	# dirty hack (AS) ref1_dirty_hack
1788	# waiting for more fundamental modifications. case when self.index_y is a bool array.
1789	# ... be careful if no point...
1790	try:
1791	if type(self.index_x[0]) == np.bool_: nx = np.sum(self.index_x) ## gives the size of the True part!
1792	if type(self.index_y[0]) == np.bool_: ny = np.sum(self.index_y) ## gives the size of the True part!
1793	except:
1794	pass
1795	# NB: ... always 4D array but possibly with "size 1" dimensions
1796	# ... if one dimension is missing because 1D 2D or 3D requests, make it appear again
1797	self.field = np.reshape(self.field,(nt,nz,ny,nx))
1798	if self.verbose: print "**** OK. I got %7.1e values. This took me %6.4f seconds" % (nxnynz*nt,timelib.time() - time0)
1799	if self.verbose: print "**** OK. I got var "+self.var+" with shape",self.field.shape
1800	# reduce coordinates to useful points
1801	# ... TBD: this should be ordered in the case of non-regular projections
1802	if self.method_x in ["free","comp"] and self.method_y in ["free","comp"]:
1803	# we need 2D coordinates (free) or we get broadcast problem (comp) so we use np.ix
1804	self.field_x = self.field_x[np.ix_(self.index_y2d, self.index_x2d)]
1805	self.field_y = self.field_y[np.ix_(self.index_y2d, self.index_x2d)]
1806	else:
1807	# we are OK with 1D coordinates
1808	self.field_x = self.field_x[self.index_y2d, self.index_x2d]
1809	self.field_y = self.field_y[self.index_y2d, self.index_x2d]
1810	# ... there are special cases with strides
1811	# ... some other fixes will be needed probably TBD
1812	if self.method_x == "free" and self.sx != 1:
1813	self.field_x = self.field_x[self.index_x]
1814	if self.method_y == "free" and self.sy != 1:
1815	self.field_y = self.field_y[self.index_y]
1816	self.field_z = self.field_z[self.index_z]
1817	self.field_t = self.field_t[self.index_t]
1818	# extract relevant horizontal points
1819	# TBD: is compfree OK with computing on irregular grid?
1820	test = self.method_x + self.method_y
1821	if test in ["fixedfixed","freefree"]:
1822	pass
1823	elif test in ["fixedfree","fixedcomp"] or test in ["freefixed","compfixed"]:
1824	if self.sx == 1 and self.sy == 1:
1825	time0 = timelib.time()
1826	# now have to obtain the new indexes which correspond to the extracted self.field
1827	# for it to work with unique index, ensure that any index_* is a numpy array
1828	if not isinstance(self.index_x, np.ndarray): self.index_x = np.array([self.index_x])
1829	if not isinstance(self.index_y, np.ndarray): self.index_y = np.array([self.index_y])
1830	if not isinstance(self.index_z, np.ndarray): self.index_z = np.array([self.index_z])
1831	if not isinstance(self.index_t, np.ndarray): self.index_t = np.array([self.index_t])
1832	for val in self.index_x: self.index_x2d[np.where(self.index_x2d == val)] = np.where(self.index_x == val)[0]
1833	for val in self.index_y: self.index_y2d[np.where(self.index_y2d == val)] = np.where(self.index_y == val)[0]
1834	##### VERY EXPENSIVE
1835	## recast self.field with 2D horizontal arrays because we might have extracted
1836	## more than what is to be actually plot or computed, in particular for comps on 2D lat,lon coordinates
1837	#self.field = self.field[np.ix_(self.index_t,self.index_z,self.index_y2d,self.index_x2d)]
1838	#(nt,nz,ny,nx) = self.field.shape
1839	# prepare the loop on all relevant horizontal points
1840	if self.method_x in ["comp","free"]:
1841	nnn = self.index_x2d.shape[0] ; what_I_am_supposed_to_do = "keepx"
1842	elif self.method_y in ["comp","free"]:
1843	nnn = self.index_y2d.shape[0] ; what_I_am_supposed_to_do = "keepy"
1844	# LOOP to extract only useful values over horizontal dimensions
1845	# only take diagonal terms, do not loop on all self.index_x2d*self.index_y2d
1846	# ... this method is fast enough, perhaps there is a faster way though
1847	# ... (for sure the method with np.diag is much slower)
1848	for iii in range(nnn):
1849	ix = self.index_x2d[iii] ; iy = self.index_y2d[iii]
1850	for iz in range(self.index_z.size):
1851	for it in range(self.index_t.size):
1852	if what_I_am_supposed_to_do == "keepx": self.field[it,iz,0,ix] = self.field[it,iz,iy,ix]
1853	elif what_I_am_supposed_to_do == "keepy": self.field[it,iz,iy,0] = self.field[it,iz,iy,ix]
1854	if self.verbose: print "**** OK. I got to pick the right values for your request. This took me %6.4f seconds" % (timelib.time() - time0)
1855	# we only keep the one value that was modified on the dimension which is not free
1856	if what_I_am_supposed_to_do == "keepx": self.field = self.field[:,:,0,:] ; ny = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
1857	elif what_I_am_supposed_to_do == "keepy": self.field = self.field[:,:,:,0] ; nx = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
1858	else:
1859	# there is a problem above if stride != 1. a fix must be found. rewrite might be necessary. TBD
1860	pass
1861	# check if 'not finite' values are present
1862	# (happens with some netCDF files where -- appears in arrays)
1863	# (but isn't it that netcdf4 returns masked arrays?)
1864	# -- we do not perform this correction for computations for which -- are handled correctly
1865	if "comp" not in self.method_t+self.method_z+self.method_y+self.method_x:
1866	w = np.where(np.isfinite(self.field) != True)
1867	self.field[w] = np.NaN
1868	## catch netCDF missing values (TBD: add a test try)
1869	#miss = self.f.variables[self.var].missing_value
1870	#if miss is not None: self.missing = miss
1871	# make a mask in case there are non-NaN missing values.
1872	# ... this is important for computations below (see ppcompute)
1873	masked = np.ma.masked_where(np.abs(self.field) >= self.missing,self.field)
1874	if masked.mask.any() == True:
1875	if self.verbose: print "!! WARNING !! Values over %5.3e are considered missing values." % self.missing
1876	self.field = masked
1877	self.field.set_fill_value([np.NaN])
1878
1879	# perform computations
1880	# -------------------------------
1881	# available: mean, max, min, meanarea
1882	# TB: integrals? for derivatives, define a function self.dx()
1883	def computations(self):
1884	nt,nz,ny,nx = self.field.shape
1885	# treat the case of mean on fields normalized with grid mesh area
1886	# ... this is done in the .area() method.
1887	# after that self.field contains field*area/totarea
1888	if "area" in self.compute:
1889	if "comp" in self.method_x+self.method_y:
1890	self.area()
1891	else:
1892	if self.verbose: print "!! WARNING !! No area accounted for (computing on t and/or z axis)."
1893	# prepare quadratic mean
1894	if "qmean" in self.compute: self.field = self.field*self.field
1895	# now ready to compute [TBD: we would like to have e.g. mean over x,y and min over t ??]
1896	if self.method_t == "comp":
1897	if self.verbose: print "**** OK. Computing over t axis."
1898	if "mean" in self.compute: self.field = ppcompute.mean(self.field,axis=0)
1899	elif self.compute == "min": self.field = ppcompute.min(self.field,axis=0)
1900	elif self.compute == "max": self.field = ppcompute.max(self.field,axis=0)
1901	else: print "!! ERROR !! operation not supported." ; exit()
1902	nt = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
1903	if self.method_z == "comp":
1904	if self.verbose: print "**** OK. Computing over z axis."
1905	if "mean" in self.compute: self.field = ppcompute.mean(self.field,axis=1)
1906	elif self.compute == "min": self.field = ppcompute.min(self.field,axis=1)
1907	elif self.compute == "max": self.field = ppcompute.max(self.field,axis=1)
1908	else: print "!! ERROR !! operation not supported." ; exit()
1909	nz = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
1910	if self.method_y == "comp":
1911	if self.verbose: print "**** OK. Computing over y axis."
1912	if self.compute == "meanarea": self.field = ppcompute.sum(self.field,axis=2)
1913	elif "mean" in self.compute: self.field = ppcompute.mean(self.field,axis=2)
1914	elif self.compute == "min": self.field = ppcompute.min(self.field,axis=2)
1915	elif self.compute == "max": self.field = ppcompute.max(self.field,axis=2)
1916	else: print "!! ERROR !! operation not supported." ; exit()
1917	ny = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
1918	if self.field_x.ndim == 2: self.field_x = self.field_x[0,:] # TBD: this is OK for regular grid but not for irregular
1919	if self.method_x == "comp":
1920	if self.verbose: print "**** OK. Computing over x axis."
1921	if self.compute == "meanarea": self.field = ppcompute.sum(self.field,axis=3)
1922	elif "mean" in self.compute: self.field = ppcompute.mean(self.field,axis=3)
1923	elif self.compute == "min": self.field = ppcompute.min(self.field,axis=3)
1924	elif self.compute == "max": self.field = ppcompute.max(self.field,axis=3)
1925	else: print "!! ERROR !! operation not supported." ; exit()
1926	nx = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
1927	if self.field_y.ndim == 2: self.field_y = self.field_y[:,0] # TBD: this is OK for regular grid but not for irregular
1928	# remove all dimensions with size 1 to prepare plot (and check the resulting dimension with dimplot)
1929	self.field = np.squeeze(self.field)
1930	# take root mean square for quadratic mean
1931	if self.compute == "qmean": self.field = np.sqrt(self.field)
1932	# error handling and verbose
1933	if self.field.ndim != self.dimplot:
1934	print "!! ERROR !! Problem: self.field is different than plot dimensions", self.field.ndim, self.dimplot ; exit()
1935	if self.verbose:
1936	print "**** OK. Final shape for "+self.var+" after averaging and squeezing",self.field.shape
1937
1938	# get areas for computations and ponderate self.field by area/totarea
1939	# -------------------------------------------------------------------
1940	def area(self):
1941	if self.verbose: print "**** OK. Get area array for computations."
1942	# create a request object for area
1943	# ... and copy known attributes from self
1944	aire = onerequest()
1945	aire.copy(self)
1946	# get area field name
1947	aire.var = "nothing"
1948	for c in glob_listarea:
1949	if c in aire.f.variables.keys():
1950	aire.var = c
1951	# do not try to calculate areas automatically
1952	if aire.var == "nothing":
1953	print "!! ERROR !! area variable not found... needs to be added in set_ppclass.txt?"
1954	exit()
1955	# define area request dimensions
1956	aire.getdim()
1957	# ensure this is a 2D horizontal request and define indexes
1958	# ... areas are not supposed to vary with time and height
1959	aire.method_x = "free" ; aire.method_y = "free"
1960	aire.getindexhori() ; aire.dimplot = 2
1961	aire.method_z = "fixed" ; aire.field_z = np.array([0]) ; aire.index_z = np.array([0])
1962	aire.method_t = "fixed" ; aire.field_t = np.array([0]) ; aire.index_t = np.array([0])
1963	# read the 2D area array in netCDF file
1964	aire.getfield()
1965	aire.field = np.squeeze(aire.field)
1966	# reduce with self horizontal indexes
1967	if "fixed" in self.method_x+self.method_y:
1968	aire.field = aire.field[self.index_y,self.index_x]
1969	# calculate total area
1970	# ... 2D comp is easy. 1D comp is a bit less easy but simple array manipulation.
1971	if "free" in self.method_x+self.method_y:
1972	if self.method_x == "free":
1973	totarea = ppcompute.sum(aire.field,axis=0)
1974	totarea = np.reshape(totarea,(1,totarea.size))
1975	totarea = np.tile(totarea,(1,self.index_x))
1976	elif self.method_y == "free":
1977	totarea = ppcompute.sum(aire.field,axis=1)
1978	totarea = np.reshape(totarea,(totarea.size,1))
1979	totarea = np.tile(totarea,(1,self.index_x.size))
1980	elif self.method_x == "comp" and self.method_y == "comp":
1981	aire.field = aire.field[np.ix_(self.index_y, self.index_x)] # reduce to requested indexes only
1982	totarea = ppcompute.sum(ppcompute.sum(aire.field,axis=1),axis=0)
1983	else:
1984	if self.verbose: print "!! WARNING !! Not account for areas. Only averaging over z and/or t axis."
1985	# normalize by total area
1986	print "**** OK. I can now normalize field by areas."
1987	aire.field = aire.field / totarea
1988	# tile area array over self t and z axis so that area field could be multiplied with self.field
1989	aire.field = np.tile(aire.field,(self.index_t.size,self.index_z.size,1,1))
1990	if self.field.shape != aire.field.shape:
1991	print "!! ERROR !! Problem in area(). Check array shapes."
1992	print "Field vs. aire:",self.field.shape,aire.field.shape ; exit()
1993	else:
1994	self.field = self.field*aire.field
1995
1996	# define coordinates for plot
1997	# -------------------------------
1998	def definecoord(self):
1999	I_got_abs = False ; I_got_ord = False
2000	# here is the thing. time is usually taken as an abscissa so we start with time.
2001	if self.method_t == "free":
2002	self.absc = self.field_t ; self.absclab = self.name_t
2003	I_got_abs = True
2004	# then we usually have x as an abscissa.
2005	if self.method_x == "free":
2006	if I_got_abs:
2007	self.ordi = self.field_x ; self.ordilab = self.name_x
2008	I_got_ord = True
2009	else:
2010	self.absc = self.field_x ; self.absclab = self.name_x
2011	I_got_abs = True
2012	# ... or we have y
2013	if self.method_y == "free":
2014	if I_got_abs:
2015	self.ordi = self.field_y ; self.ordilab = self.name_y
2016	I_got_ord = True
2017	else:
2018	self.absc = self.field_y ; self.absclab = self.name_y
2019	I_got_abs = True
2020	# ... and we end with z because it is usually not an abscissa (profiles).
2021	if self.method_z == "free":
2022	if self.field_z[0] > self.field_z[1]:
2023	self.invert_axes = True # the axis will be turned upside-down
2024	if I_got_abs:
2025	self.ordi = self.field_z ; self.ordilab = self.name_z
2026	I_got_ord = True
2027	else:
2028	self.absc = self.field_z ; self.absclab = self.name_z
2029	I_got_abs = True
2030	self.swap_axes = True # says that altitude is not supposed to remain as an abscissa
2031	if I_got_abs and self.verbose: print "**** OK. abscissa:",self.absclab, self.absc.shape
2032	if I_got_ord and self.verbose: print "**** OK. ordinate:",self.ordilab, self.ordi.shape

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

Download in other formats: