Context Navigation

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

Last change on this file since 1242 was 1219, checked in by aslmd, 11 years ago
PLANETOPLOT. improve flexibility of code for computations. added computation of perturbations and diff
File size: 100.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 "---------------------------------------------------------------"
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	# if longitudes are crossing 180 deg limit, then modify longitude vector
1454	self.field_x = self.field_x[0,:,:]
1455	for j in range(self.field_x[0,:].size-1):
1456	if (self.field_x[0,j+1]-self.field_x[0,j] < 0.):
1457	self.field_x[:,j+1]=self.field_x[:,j+1]+360.
1458	self.field_y = self.field_y[0,:,:]
1459	# if xy axis are apparently undefined, set 2D grid points axis.
1460	if "grid points" not in self.name_x:
1461	if np.all(self.field_x == self.field_x[0,0]) \
1462	or self.field_x.min() == self.field_x.max() \
1463	or self.field_y.min() == self.field_y.max():
1464	if self.verbose: print "!! WARNING !! xy axis look undefined. creating non-dummy ones."
1465	self.field_x = np.array(range(self.dim_x)) ; self.name_x = "x grid points"
1466	self.field_y = np.array(range(self.dim_y)) ; self.name_y = "y grid points"
1467	[self.field_x,self.field_y] = np.meshgrid(self.field_x,self.field_y)
1468	if self.dim_x > 1:
1469	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())
1470	if self.dim_y > 1:
1471	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())
1472	# ALTITUDE. Try preset fields. If not present set grid points axis.
1473	# WARNING: how do we do if several are available? the last one is chosen.
1474	self.name_z = "nothing"
1475	for c in glob_listz:
1476	if c in self.f.variables.keys():
1477	self.name_z = c
1478	if self.name_z == "nothing":
1479	self.field_z = np.array(range(self.dim_z))
1480	self.name_z = "z grid points"
1481	else:
1482	tabalt = self.f.variables[self.name_z]
1483	# (consider the case where tabtime is not dim 1) TBD: 2D and 3D cases
1484	if tabalt.ndim == 4:
1485	try:
1486	self.field_z = tabalt[1,:,0,0] # 4D case. alt is usually third dimension.
1487	# 1 for time to avoid initial 0s
1488	except:
1489	self.field_z = tabalt[0,:,0,0]
1490	if self.verbose: print "!! WARNING !! "+self.name_z+" is 4D var. We made it 1D."
1491	else:
1492	self.field_z = self.f.variables[self.name_z][:] # specify dimension
1493	# TBD: problems when self.dim_z != self.field_z.size
1494	if self.field_z.size != self.dim_z:
1495	if self.verbose: print "!! WARNING !! Cannot use this z coordinate. Not enough points. Use simple z axis."
1496	self.field_z = np.array(range(self.dim_z))
1497	self.name_z = "z grid points"
1498	if self.dim_z > 1:
1499	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())
1500
1501	# TIME. Try preset fields.
1502	self.name_t = "nothing"
1503	for c in glob_listt:
1504	if c in self.f.variables.keys():
1505	self.name_t = c
1506	if self.verbose: print "**** OK. Found time variable: ",c
1507	try:
1508	# speed up: only get first value, last one.
1509	self.tabtime = self.f.variables[self.name_t]
1510	# (consider the case where tabtime is not dim 1)
1511	# (time is most often the first dimension)
1512	if self.tabtime.ndim == 2: self.tabtime = self.tabtime[:,0]
1513	elif self.tabtime.ndim == 3: self.tabtime = self.tabtime[:,0,0]
1514	elif self.tabtime.ndim == 4: self.tabtime = self.tabtime[:,0,0,0]
1515	# (now proceed) (the +0. just ensures this is a number)
1516	dafirst = self.tabtime[0] + 0.
1517	if self.dim_t == 1:
1518	self.field_t = np.array([dafirst])
1519	else:
1520	daint = self.tabtime[1] - dafirst
1521	dalast = dafirst + (self.dim_t-1)*daint
1522	self.field_t = np.linspace(dafirst,dalast,num=self.dim_t)
1523	if self.verbose:
1524	print "!! WARNING !! WARNING !! Time axis is supposed to be equally spaced !!"
1525	if dalast != self.tabtime[self.dim_t-1]:
1526	print "!! WARNING !! Time axis has been recast to be monotonic",dalast,self.tabtime[self.dim_t-1]
1527	except:
1528	# ... or if a problem encountered, define a simple time axis
1529	if self.verbose: print "**** OK. There is something weird. Let us go for a simple time axis."
1530	self.field_t = np.array(range(self.dim_t))
1531	self.name_t = "t grid points"
1532	if self.dim_t > 1:
1533	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())
1534
1535	# change time axis
1536	# ... add your options here!
1537	# --------------------------
1538	def performtimechange(self):
1539	if self.changetime is not None:
1540	if self.verbose: print "**** OK. Converting time axis:",self.changetime
1541	### options added by T. Navarro
1542	if self.changetime == "mars_sol2ls":
1543	if "controle" in self.f.variables:
1544	self.field_t = self.field_t \
1545	+ self.f.variables['controle'][3]%669 \
1546	+ self.f.variables['controle'][26]
1547	self.field_t = ppcompute.mars_sol2ls(self.field_t)
1548	elif self.changetime == "mars_dayini" and "controle" in self.f.variables:
1549	self.field_t = self.field_t \
1550	+ self.f.variables['controle'][3]%669 \
1551	+ self.f.variables['controle'][26]
1552	### options added by A. Spiga
1553	elif self.changetime == "correctls":
1554	if self.tabtime is None:
1555	print "!! WARNING !! Encountered a problem with correctls. Check your time dimension. Skipping this part."
1556	else:
1557	dafirst = self.tabtime[0] + 0.
1558	if self.dim_t == 1:
1559	self.field_t = np.array([dafirst])
1560	else:
1561	daint = self.tabtime[1] - dafirst
1562	dalast = dafirst + (self.dim_t-1)*daint
1563	year = 0.
1564	add = np.linspace(dafirst,dalast,num=self.dim_t) ; add[0] = 0.
1565	for iii in range(1,self.dim_t):
1566	if self.tabtime[iii] - self.tabtime[iii-1] < 0: year = year+1.
1567	add[iii] = year*360.
1568	self.field_t = add + self.tabtime
1569	elif "mars_meso" in self.changetime:
1570	if 'Times' not in self.f.variables.keys():
1571	if self.verbose: print "!! WARNING !! Variable Times not in file. Cannot proceed to change of time axis."
1572	else:
1573	# get the array of strings describing dates
1574	dates = self.f.variables['Times']
1575	dates.set_auto_maskandscale(False) # necessary to solve the api Times bug!
1576	# get ls sol utc from those strings
1577	ls, sol, utc = ppcompute.mars_date(dates[:])
1578	# populate self.field_t with the right output from mars_date
1579	if self.changetime == "mars_meso_ls":
1580	self.field_t = ls
1581	self.name_t = "Ls"
1582	elif self.changetime == "mars_meso_sol":
1583	self.field_t = sol
1584	self.name_t = "sol"
1585	elif self.changetime == "mars_meso_utc" \
1586	and ( self.changetime == "mars_meso_lt" \
1587	and not hasattr(self.f,'CEN_LON') ):
1588	self.field_t = ppcompute.timecorrect(utc)
1589	self.name_t = "utc"
1590	if self.method_t == "fixed":
1591	self.field_t = self.field_t % 24 # so that the user is not mistaken!
1592	elif self.changetime == "mars_meso_lt":
1593	self.field_t = ppcompute.timecorrect(utc) + getattr(self.f,'CEN_LON') / 15.
1594	self.field_t = ppcompute.timecorrect(self.field_t)
1595	self.name_t = "local time (center of domain)"
1596	if self.method_t == "fixed":
1597	self.field_t = self.field_t % 24 # so that the user is not mistaken!
1598	else:
1599	if self.verbose: print "!! WARNING !! This time change is not implemented. Nothing is done."
1600	if self.verbose: print "**** OK. new t axis values [%5.1f,%5.1f]" % (self.field_t.min(),self.field_t.max())
1601
1602	# get list of index to be retrieved for time axis
1603	### TBD: il faudrait ne prendre que les indices qui correspondent a l interieur d un plot (dans all)
1604	# -------------------------------
1605	def getindextime(self,dalist=None,ind=None):
1606	if self.method_t == "free":
1607	self.index_t = np.arange(0,self.dim_t,self.st)
1608	if self.dim_t > 1:
1609	self.dimplot = self.dimplot + 1
1610	if self.verbose: print "**** OK. t values. all."
1611	else:
1612	self.method_t = "fixed"
1613	if self.verbose: print "**** OK. no t dimension."
1614	elif self.method_t == "comp":
1615	start = np.argmin( np.abs( self.field_t - dalist[ind][0] ) )
1616	stop = np.argmin( np.abs( self.field_t - dalist[ind][1] ) )
1617	self.index_t = np.arange(start,stop+1,self.st)
1618	if self.verbose: print "**** OK. t values. comp over interval ",self.field_t[start],self.field_t[stop]," nvalues=",self.index_t.size
1619	elif self.method_t == "fixed":
1620	self.index_t.append( np.argmin( np.abs( self.field_t - dalist[ind][0] ) ))
1621	if self.verbose: print "**** OK. t values",self.field_t[self.index_t]
1622	else:
1623	print "!! ERROR !! method "+self.method_t+" not supported"
1624	self.index_t = np.array(self.index_t)
1625
1626	# get list of index to be retrieved for vertical axis
1627	### TBD: il faudrait ne prendre que les indices qui correspondent a l interieur d un plot (dans all)
1628	# -------------------------------
1629	def getindexvert(self,dalist=None,ind=None):
1630	if self.method_z == "free":
1631	self.index_z = np.arange(0,self.dim_z,self.sz)
1632	if self.dim_z > 1:
1633	self.dimplot = self.dimplot + 1
1634	if self.verbose: print "**** OK. z values. all."
1635	else:
1636	self.method_z = "fixed"
1637	if self.verbose: print "**** OK. no z dimension."
1638	elif self.method_z == "comp":
1639	start = np.argmin( np.abs( self.field_z - dalist[ind][0] ) )
1640	stop = np.argmin( np.abs( self.field_z - dalist[ind][1] ) )
1641	self.index_z = np.arange(start,stop+1,self.sz)
1642	if self.verbose: print "**** OK. z values. comp over interval",self.field_z[start],self.field_z[stop]," nvalues=",self.index_z.size
1643	elif self.method_z == "fixed":
1644	self.index_z.append( np.argmin( np.abs( self.field_z - dalist[ind][0] ) ))
1645	if self.verbose: print "**** OK. z values",self.field_z[self.index_z]
1646	else:
1647	if self.verbose: print "!! ERROR !! method "+self.method_z+" not supported"
1648	self.index_z = np.array(self.index_z)
1649
1650	# get list of index to be retrieved for horizontal grid
1651	# --> index_x and index_y are slices to be retrieved from NETCDF files
1652	# --> index_x2d and index_y2d are the actual (x,y) coordinates corresponding to each relevant point
1653	# [this is slightly more complicated because 2D arrays for lat-lon projection possibly irregular]
1654	# NB: to append index we use lists (the most convenient) then we convert into a numpy.array
1655	### TBD: il faudrait ne prendre que les indices qui correspondent a l interieur d un plot (dans all)
1656	# -------------------------------
1657	def getindexhori(self,dalistx=None,dalisty=None,indx=None,indy=None):
1658	## get what is the method over x and y axis
1659	test = self.method_x+self.method_y
1660	## CASE 0, EASY CASES:
1661	## - LAT IS FREE (we do here what must be done whatever LON case is)
1662	## - LON IS FREE (we do here what must be done whatever LAT case is)
1663	## - LAT IS COMP AND LON IS FREE
1664	## - LON IS COMP AND LAT IS FREE
1665	if self.method_x == "free" or test in ["compfree","compcomp"]:
1666	self.index_x = range(0,self.dim_x,self.sx)
1667	if self.dim_x > 1:
1668	if self.method_x == "free": self.dimplot = self.dimplot + 1
1669	if self.verbose: print "**** OK. x values. all."
1670	else:
1671	self.method_x = "fixed"
1672	if self.verbose: print "**** OK. no x dimension."
1673	if self.method_y == "free" or test in ["freecomp","compcomp"]:
1674	self.index_y = range(0,self.dim_y,self.sy)
1675	if self.dim_y > 1:
1676	if self.method_y == "free": self.dimplot = self.dimplot + 1
1677	if self.verbose: print "**** OK. y values. all."
1678	else:
1679	self.method_y = "fixed"
1680	if self.verbose: print "**** OK. no y dimension."
1681	## CASE 0 above, this is just for continuity for free.
1682	## ... for comp we have to select bounds.
1683	## ... TBD: take the bool array strategy for what follows!
1684	if self.method_x in ["free","comp"] and self.method_y in ["free","comp"]:
1685	### ref1_dirty_hack
1686	### ... for the moment this is a hack. but actually this is more powerful.
1687	if self.method_x == "comp":
1688	yeah = (self.field_x >= dalistx[indx][0])*(self.field_x <= dalistx[indx][1])
1689	self.index_x = yeah[0,:]
1690	if self.method_y == "comp":
1691	yeah = (self.field_y >= dalisty[indy][0])*(self.field_y <= dalisty[indy][1])
1692	self.index_y = yeah[:,0]
1693	self.index_x2d = self.index_x
1694	self.index_y2d = self.index_y
1695	## AND NOW THE LITTLE BIT MORE COMPLICATED CASES
1696	## CASE 1 LAT AND LON ARE FIXED
1697	elif test == "fixedfixed":
1698	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 )
1699	#TBD: pb with staggered coord
1700	if idx not in self.index_x: self.index_x.append(idx)
1701	if idy not in self.index_y: self.index_y.append(idy)
1702	self.index_x2d.append(idx)
1703	self.index_y2d.append(idy)
1704	## CASE 2 LON IS FIXED BUT NOT LAT
1705	elif test in ["fixedfree","fixedcomp"]:
1706	# find where are requested x values for each y on the free dimension
1707	# NB: this does not work for non-bijective cases e.g. polar stereographic
1708	for iy in range(self.dim_y):
1709	idx = np.argmin( np.abs( self.field_x[iy,:] - dalistx[indx][0] ) )
1710	# if comp is requested we select only indexes which yield values between requested min and max
1711	storeval = (self.method_y == "comp") and (self.field_y[iy,idx] > dalisty[indy][0]) and (self.field_y[iy,idx] < dalisty[indy][1])
1712	storeval = storeval or (self.method_y == "free")
1713	if storeval:
1714	if idx not in self.index_x: self.index_x.append(idx)
1715	if iy not in self.index_y and self.method_y == "comp": self.index_y.append(iy)
1716	if idx not in self.index_x2d or iy not in self.index_y2d:
1717	self.index_x2d.append(idx)
1718	self.index_y2d.append(iy)
1719	## CASE 3 LAT IS FIXED BUT NOT LON
1720	elif test in ["freefixed","compfixed"]:
1721	# find where are requested y values for each x on the free dimension
1722	# NB: this does not work for non-bijective cases e.g. polar stereographic
1723	for ix in range(self.dim_x):
1724	idy = np.argmin( np.abs( self.field_y[:,ix] - dalisty[indy][0] ) )
1725	# if comp is requested we select only indexes which yield values between requested min and max
1726	storeval = (self.method_x == "comp") and (self.field_x[idy,ix] > dalistx[indx][0]) and (self.field_x[idy,ix] < dalistx[indx][1])
1727	storeval = storeval or (self.method_x == "free")
1728	if storeval:
1729	if idy not in self.index_y: self.index_y.append(idy)
1730	if ix not in self.index_x and self.method_x == "comp": self.index_x.append(ix)
1731	if ix not in self.index_x2d or idy not in self.index_y2d:
1732	self.index_x2d.append(ix)
1733	self.index_y2d.append(idy)
1734	## check index tab
1735	if len(self.index_x) == 0 or len(self.index_y) == 0:
1736	print "!! ERROR !! no indices found. check prescribed latitudes or longitudes" ; exit()
1737	## ensure the array is a numpy array for getfield to work
1738	self.index_x = np.array(self.index_x)
1739	self.index_y = np.array(self.index_y)
1740	self.index_x2d = np.array(self.index_x2d)
1741	self.index_y2d = np.array(self.index_y2d)
1742	### print extrema
1743	printx = self.field_x[np.ix_(self.index_y2d, self.index_x2d)]
1744	printy = self.field_y[np.ix_(self.index_y2d, self.index_x2d)]
1745	if self.verbose:
1746	print "**** OK. x values (min,max).", printx.min(),printx.max()
1747	print "**** OK. y values (min,max).", printy.min(),printy.max()
1748
1749	# get the field from the NETCDF file and perform averages
1750	# -------------------------------
1751	def getfield(self):
1752	## first tell what is to be done
1753	if self.verbose:
1754	if self.dimplot > 2: print "**** !! WARNING !! "+str(self.dimplot)+"D plots will not be supported!"
1755	elif self.dimplot == 0 and self.verbose: print "**** OK. 0D value requested."
1756	elif self.dimplot == 1 and self.verbose: print "**** OK. 1D plot requested."
1757	elif self.verbose: print "**** OK. 2D section requested."
1758	# well, now get field from netcdf file
1759	# part below is necessary otherwise there is an index error below
1760	if self.index_x.size == 1: self.index_x = self.index_x[0]
1761	if self.index_y.size == 1: self.index_y = self.index_y[0]
1762	if self.index_z.size == 1: self.index_z = self.index_z[0]
1763	if self.index_t.size == 1: self.index_t = self.index_t[0]
1764	# set everything about dimensions and slices
1765	# -- each self.dim case corresponds to tests in the beginning of getdim.
1766	if self.dim == 1:
1767	nt = self.index_t.size ; nz = 1 ; ny = 1 ; nx = 1
1768	tupledim = self.index_t
1769	elif self.dim == 2:
1770	nt = 1 ; nz = 1 ; ny = self.index_y.size ; nx = self.index_x.size
1771	tupledim = self.index_y,self.index_x
1772	elif self.dim == 3:
1773	## DEFAULT tyx (time-varying 2D field)
1774	if self.kind3d == "tyx":
1775	nt = self.index_t.size ; nz = 1 ; ny = self.index_y.size ; nx = self.index_x.size
1776	tupledim = self.index_t,self.index_y,self.index_x
1777	## CASE tzy (e.g. time-varying zonal mean y-z field)
1778	elif self.kind3d == "tzy":
1779	nt = self.index_t.size ; nz = self.index_z.size ; ny = self.index_y.size ; nx = 1
1780	tupledim = self.index_t,self.index_z,self.index_y
1781	else:
1782	print "!! ERROR !! This kind of 3D field is not supported. Please send feedback." ; exit()
1783	# this is far faster than retrieving each term with a loop
1784	elif self.dim == 4:
1785	nt = self.index_t.size ; nz = self.index_z.size ; ny = self.index_y.size ; nx = self.index_x.size
1786	tupledim = self.index_t,self.index_z,self.index_y,self.index_x
1787	else:
1788	print "!! ERROR !! field would have more than four dimensions ?" ; exit()
1789	# then retrieve what is requested by user!
1790	time0 = timelib.time()
1791	if self.verbose: print "**** OK. I am getting values from files. Please wait."
1792	self.field = self.f.variables[self.var][tupledim]
1793	# dirty hack (AS) ref1_dirty_hack
1794	# waiting for more fundamental modifications. case when self.index_y is a bool array.
1795	# ... be careful if no point...
1796	try:
1797	if type(self.index_x[0]) == np.bool_: nx = np.sum(self.index_x) ## gives the size of the True part!
1798	if type(self.index_y[0]) == np.bool_: ny = np.sum(self.index_y) ## gives the size of the True part!
1799	except:
1800	pass
1801	# NB: ... always 4D array but possibly with "size 1" dimensions
1802	# ... if one dimension is missing because 1D 2D or 3D requests, make it appear again
1803	self.field = np.reshape(self.field,(nt,nz,ny,nx))
1804	if self.verbose: print "**** OK. I got %7.1e values. This took me %6.4f seconds" % (nxnynz*nt,timelib.time() - time0)
1805	if self.verbose: print "**** OK. I got var "+self.var+" with shape",self.field.shape
1806	# reduce coordinates to useful points
1807	# ... TBD: this should be ordered in the case of non-regular projections
1808	if self.method_x in ["free","comp"] and self.method_y in ["free","comp"]:
1809	# we need 2D coordinates (free) or we get broadcast problem (comp) so we use np.ix
1810	self.field_x = self.field_x[np.ix_(self.index_y2d, self.index_x2d)]
1811	self.field_y = self.field_y[np.ix_(self.index_y2d, self.index_x2d)]
1812	else:
1813	# we are OK with 1D coordinates
1814	self.field_x = self.field_x[self.index_y2d, self.index_x2d]
1815	self.field_y = self.field_y[self.index_y2d, self.index_x2d]
1816	# ... there are special cases with strides
1817	# ... some other fixes will be needed probably TBD
1818	if self.method_x == "free" and self.sx != 1:
1819	self.field_x = self.field_x[self.index_x]
1820	if self.method_y == "free" and self.sy != 1:
1821	self.field_y = self.field_y[self.index_y]
1822	self.field_z = self.field_z[self.index_z]
1823	self.field_t = self.field_t[self.index_t]
1824	# extract relevant horizontal points
1825	# TBD: is compfree OK with computing on irregular grid?
1826	test = self.method_x + self.method_y
1827	if test in ["fixedfixed","freefree"]:
1828	pass
1829	elif test in ["fixedfree","fixedcomp"] or test in ["freefixed","compfixed"]:
1830	if self.sx == 1 and self.sy == 1:
1831	time0 = timelib.time()
1832	# now have to obtain the new indexes which correspond to the extracted self.field
1833	# for it to work with unique index, ensure that any index_* is a numpy array
1834	if not isinstance(self.index_x, np.ndarray): self.index_x = np.array([self.index_x])
1835	if not isinstance(self.index_y, np.ndarray): self.index_y = np.array([self.index_y])
1836	if not isinstance(self.index_z, np.ndarray): self.index_z = np.array([self.index_z])
1837	if not isinstance(self.index_t, np.ndarray): self.index_t = np.array([self.index_t])
1838	for val in self.index_x: self.index_x2d[np.where(self.index_x2d == val)] = np.where(self.index_x == val)[0]
1839	for val in self.index_y: self.index_y2d[np.where(self.index_y2d == val)] = np.where(self.index_y == val)[0]
1840	##### VERY EXPENSIVE
1841	## recast self.field with 2D horizontal arrays because we might have extracted
1842	## more than what is to be actually plot or computed, in particular for comps on 2D lat,lon coordinates
1843	#self.field = self.field[np.ix_(self.index_t,self.index_z,self.index_y2d,self.index_x2d)]
1844	#(nt,nz,ny,nx) = self.field.shape
1845	# prepare the loop on all relevant horizontal points
1846	if self.method_x in ["comp","free"]:
1847	nnn = self.index_x2d.shape[0] ; what_I_am_supposed_to_do = "keepx"
1848	elif self.method_y in ["comp","free"]:
1849	nnn = self.index_y2d.shape[0] ; what_I_am_supposed_to_do = "keepy"
1850	# LOOP to extract only useful values over horizontal dimensions
1851	# only take diagonal terms, do not loop on all self.index_x2d*self.index_y2d
1852	# ... this method is fast enough, perhaps there is a faster way though
1853	# ... (for sure the method with np.diag is much slower)
1854	for iii in range(nnn):
1855	ix = self.index_x2d[iii] ; iy = self.index_y2d[iii]
1856	for iz in range(self.index_z.size):
1857	for it in range(self.index_t.size):
1858	if what_I_am_supposed_to_do == "keepx": self.field[it,iz,0,ix] = self.field[it,iz,iy,ix]
1859	elif what_I_am_supposed_to_do == "keepy": self.field[it,iz,iy,0] = self.field[it,iz,iy,ix]
1860	if self.verbose: print "**** OK. I got to pick the right values for your request. This took me %6.4f seconds" % (timelib.time() - time0)
1861	# we only keep the one value that was modified on the dimension which is not free
1862	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))
1863	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))
1864	else:
1865	# there is a problem above if stride != 1. a fix must be found. rewrite might be necessary. TBD
1866	pass
1867	# check if 'not finite' values are present
1868	# (happens with some netCDF files where -- appears in arrays)
1869	# (but isn't it that netcdf4 returns masked arrays?)
1870	# -- we do not perform this correction for computations for which -- are handled correctly
1871	if "comp" not in self.method_t+self.method_z+self.method_y+self.method_x:
1872	w = np.where(np.isfinite(self.field) != True)
1873	self.field[w] = np.NaN
1874	## catch netCDF missing values (TBD: add a test try)
1875	#miss = self.f.variables[self.var].missing_value
1876	#if miss is not None: self.missing = miss
1877	# make a mask in case there are non-NaN missing values.
1878	# ... this is important for computations below (see ppcompute)
1879	masked = np.ma.masked_where(np.abs(self.field) >= self.missing,self.field)
1880	if masked.mask.any() == True:
1881	if self.verbose: print "!! WARNING !! Values over %5.3e are considered missing values." % self.missing
1882	self.field = masked
1883	self.field.set_fill_value([np.NaN])
1884
1885	# perform computations
1886	# -------------------------------
1887	# available: mean, max, min, meanarea
1888	# TB: integrals? for derivatives, define a function self.dx()
1889	def computations(self):
1890	nt,nz,ny,nx = self.field.shape
1891	# treat the case of mean on fields normalized with grid mesh area
1892	# ... this is done in the .area() method.
1893	# after that self.field contains field*area/totarea
1894	if "area" in self.compute:
1895	if "comp" in self.method_x+self.method_y:
1896	self.area()
1897	else:
1898	if self.verbose: print "!! WARNING !! No area accounted for (computing on t and/or z axis)."
1899	# prepare quadratic mean
1900	if "qmean" in self.compute: self.field = self.field*self.field
1901	# prepare and execute (possibly sequential) means
1902	roll = [self.method_t, self.method_z, self.method_y, self.method_x]
1903	reshapedim = [nt,nz,ny,nx]
1904	for nr in range(len(roll)):
1905	rrr = roll[nr]
1906	if "comp" in rrr:
1907	# a. computing
1908	if self.verbose: print "**** OK. Computing over axis number ",zeaxis
1909	if self.compute == "meanarea": self.field = ppcompute.sum (self.field,axis=nr)
1910	elif "mean" in self.compute: self.field = ppcompute.mean (self.field,axis=nr)
1911	elif self.compute == "min": self.field = ppcompute.min (self.field,axis=nr)
1912	elif self.compute == "max": self.field = ppcompute.max (self.field,axis=nr)
1913	else: print "!! ERROR !! operation not supported." ; exit()
1914	# b. reshaping
1915	reshapedim[nr] = 1
1916	self.field = np.reshape(self.field,reshapedim)
1917	# c. particular cases for 2D coordinates
1918	if nr == 2 and 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 nr == 3 and self.field_y.ndim == 2: self.field_y = self.field_y[:,0] # TBD: this is OK for regular grid but not for irregular
1920	# computations for which only setting compute is needed
1921	if "_x" in self.compute: zeaxis = 3
1922	elif "_y" in self.compute: zeaxis = 2
1923	elif "_z" in self.compute: zeaxis = 1
1924	elif "_t" in self.compute: zeaxis = 0
1925	if "pert" in self.compute:
1926	self.field = ppcompute.perturbation(self.field,axis=zeaxis)
1927	elif "diff" in self.compute:
1928	dadiff = np.diff(self.field,axis=zeaxis)
1929	# (diff ouput has one value less in impacted dimension. fix this.)
1930	reshapedim[zeaxis] = reshapedim[zeaxis]-1
1931	self.field[:,:,:,:] = np.nan
1932	self.field[0:reshapedim[0],0:reshapedim[1],0:reshapedim[2],0:reshapedim[3]] = dadiff
1933	# remove all dimensions with size 1 to prepare plot (and check the resulting dimension with dimplot)
1934	self.field = np.squeeze(self.field)
1935	# take root mean square for quadratic mean
1936	if self.compute == "qmean": self.field = np.sqrt(self.field)
1937	# error handling and verbose
1938	if self.field.ndim != self.dimplot:
1939	print "!! ERROR !! Problem: self.field is different than plot dimensions", self.field.ndim, self.dimplot ; exit()
1940	if self.verbose:
1941	print "**** OK. Final shape for "+self.var+" after averaging and squeezing",self.field.shape
1942
1943	# get areas for computations and ponderate self.field by area/totarea
1944	# -------------------------------------------------------------------
1945	def area(self):
1946	if self.verbose: print "**** OK. Get area array for computations."
1947	# create a request object for area
1948	# ... and copy known attributes from self
1949	aire = onerequest()
1950	aire.copy(self)
1951	# get area field name
1952	aire.var = "nothing"
1953	for c in glob_listarea:
1954	if c in aire.f.variables.keys():
1955	aire.var = c
1956	# do not try to calculate areas automatically
1957	if aire.var == "nothing":
1958	print "!! ERROR !! area variable not found... needs to be added in set_ppclass.txt?"
1959	exit()
1960	# define area request dimensions
1961	aire.getdim()
1962	# ensure this is a 2D horizontal request and define indexes
1963	# ... areas are not supposed to vary with time and height
1964	aire.method_x = "free" ; aire.method_y = "free"
1965	aire.getindexhori() ; aire.dimplot = 2
1966	aire.method_z = "fixed" ; aire.field_z = np.array([0]) ; aire.index_z = np.array([0])
1967	aire.method_t = "fixed" ; aire.field_t = np.array([0]) ; aire.index_t = np.array([0])
1968	# read the 2D area array in netCDF file
1969	aire.getfield()
1970	aire.field = np.squeeze(aire.field)
1971	# reduce with self horizontal indexes
1972	if "fixed" in self.method_x+self.method_y:
1973	aire.field = aire.field[self.index_y,self.index_x]
1974	# calculate total area
1975	# ... 2D comp is easy. 1D comp is a bit less easy but simple array manipulation.
1976	if "free" in self.method_x+self.method_y:
1977	if self.method_x == "free":
1978	totarea = ppcompute.sum(aire.field,axis=0)
1979	totarea = np.reshape(totarea,(1,totarea.size))
1980	totarea = np.tile(totarea,(1,self.index_x))
1981	elif self.method_y == "free":
1982	totarea = ppcompute.sum(aire.field,axis=1)
1983	totarea = np.reshape(totarea,(totarea.size,1))
1984	totarea = np.tile(totarea,(1,self.index_x.size))
1985	elif self.method_x == "comp" and self.method_y == "comp":
1986	aire.field = aire.field[np.ix_(self.index_y, self.index_x)] # reduce to requested indexes only
1987	totarea = ppcompute.sum(ppcompute.sum(aire.field,axis=1),axis=0)
1988	else:
1989	if self.verbose: print "!! WARNING !! Not account for areas. Only averaging over z and/or t axis."
1990	# normalize by total area
1991	print "**** OK. I can now normalize field by areas."
1992	aire.field = aire.field / totarea
1993	# tile area array over self t and z axis so that area field could be multiplied with self.field
1994	aire.field = np.tile(aire.field,(self.index_t.size,self.index_z.size,1,1))
1995	if self.field.shape != aire.field.shape:
1996	print "!! ERROR !! Problem in area(). Check array shapes."
1997	print "Field vs. aire:",self.field.shape,aire.field.shape ; exit()
1998	else:
1999	self.field = self.field*aire.field
2000
2001	# define coordinates for plot
2002	# -------------------------------
2003	def definecoord(self):
2004	I_got_abs = False ; I_got_ord = False
2005	# here is the thing. time is usually taken as an abscissa so we start with time.
2006	if self.method_t == "free":
2007	self.absc = self.field_t ; self.absclab = self.name_t
2008	I_got_abs = True
2009	# then we usually have x as an abscissa.
2010	if self.method_x == "free":
2011	if I_got_abs:
2012	self.ordi = self.field_x ; self.ordilab = self.name_x
2013	I_got_ord = True
2014	else:
2015	self.absc = self.field_x ; self.absclab = self.name_x
2016	I_got_abs = True
2017	# ... or we have y
2018	if self.method_y == "free":
2019	if I_got_abs:
2020	self.ordi = self.field_y ; self.ordilab = self.name_y
2021	I_got_ord = True
2022	else:
2023	self.absc = self.field_y ; self.absclab = self.name_y
2024	I_got_abs = True
2025	# ... and we end with z because it is usually not an abscissa (profiles).
2026	if self.method_z == "free":
2027	if self.field_z[0] > self.field_z[1]:
2028	self.invert_axes = True # the axis will be turned upside-down
2029	if I_got_abs:
2030	self.ordi = self.field_z ; self.ordilab = self.name_z
2031	I_got_ord = True
2032	else:
2033	self.absc = self.field_z ; self.absclab = self.name_z
2034	I_got_abs = True
2035	self.swap_axes = True # says that altitude is not supposed to remain as an abscissa
2036	if I_got_abs and self.verbose: print "**** OK. abscissa:",self.absclab, self.absc.shape
2037	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: