Context Navigation

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

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