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

Rev	Line
[910]	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...
[960]	27	# ... this can be hardcoded here
[910]	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 = []
[921]	33	glob_listarea = []
[910]	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)
[921]	40	for stuff in lines[17].strip().split(';'): glob_listarea.append(stuff)
[910]	41	except IOError:
[960]	42	print "PPCLASS warning: "+zefile+" not in "+whereset+" ; no presets."
[910]	43
	44	##################################
	45	#### USEFUL GENERIC FUNCTIONS ####
	46	##################################
	47
	48	# inspect variables and dimensions in a netCDF file
	49	def inspect(filename):
[1083]	50	print "---------------------------------------------------------------"
[910]	51	print "**** INSPECT FILE",filename
	52	test = netCDF4.Dataset(filename)
[1083]	53	print "---------------------------------------------------------------"
	54	print "**** DIMENSIONS:"
	55	for obj in test.dimensions.values(): print obj
	56	print "---------------------------------------------------------------"
[910]	57	print "**** VARIABLES: ",test.variables.keys()
[1083]	58	print "---------------------------------------------------------------"
[1082]	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+" )"
[1083]	64	print "---------------------------------------------------------------"
[1082]	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:
[1083]	78	print message+"not found. will simply use index (check out dimensions)."
[910]	79
[1054]	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
[910]	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
[930]	130
	131	# look for comas in the input name to separate different names (files, variables,etc ..)
[910]	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,\
[1050]	162	sx=1,sy=1,\
	163	sz=1,st=1,\
	164	svx=1,\
	165	svy=1,\
[910]	166	compute="mean",\
[1108]	167	kind3d="tyx",\
[960]	168	verbose=False,\
	169	quiet=False,\
	170	noproj=False,\
[910]	171	superpose=False,\
	172	plotin=None,\
	173	forcedimplot=-1,\
	174	out="gui",\
	175	filename="myplot",\
[923]	176	folder="./",\
[930]	177	includedate=True,\
[964]	178	res=150.,\
[923]	179	xlabel=None,ylabel=None,\
	180	xcoeff=None,ycoeff=None,\
[1158]	181	xmin=None,ymin=None,\
	182	xmax=None,ymax=None,\
[1029]	183	nxticks=10,nyticks=10,\
[923]	184	proj=None,\
	185	vmin=None,vmax=None,\
	186	div=None,\
[1050]	187	colorbar=None,\
	188	linestyle=None,\
[923]	189	marker=None,\
	190	color=None,\
[1050]	191	legend=None,\
[934]	192	changetime=None,\
[960]	193	units=None,\
[990]	194	savtxt=False,\
[991]	195	modx=None,\
[1054]	196	fmt=None,\
[1029]	197	xp=16,yp=8,\
	198	missing=1.e25,\
[1062]	199	trans=None,back=None,\
	200	showcb=None,\
[1067]	201	logy=None,\
[923]	202	title=None):
[910]	203	self.request = None
[960]	204	self.nrequest = 0
[910]	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
[921]	213	self.customplot = False
[960]	214	self.f = None
	215	self.l = None
[1201]	216	self.xx = None
	217	self.yy = None
	218	self.zz = None
	219	self.tt = None
[910]	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
[1050]	227	self.sx = sx ; self.sy = sy
	228	self.sz = sz ; self.st = st
	229	self.svx = svx
	230	self.svy = svy
[910]	231	self.compute = compute
[1108]	232	self.kind3d = kind3d
[910]	233	self.verbose = verbose
[960]	234	self.quiet = quiet
[910]	235	self.noproj = noproj
	236	self.plotin = plotin
	237	self.superpose = superpose
	238	self.forcedimplot = forcedimplot
	239	self.out = out
	240	self.filename = filename
[964]	241	self.res = res
[910]	242	self.folder = folder
[930]	243	self.includedate = includedate
[934]	244	self.changetime = changetime
[990]	245	self.savtxt = savtxt
[991]	246	self.modx = modx
[1029]	247	self.missing = missing
[923]	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
[1158]	252	self.xmin = xmin ; self.xmax = xmax
	253	self.ymin = ymin ; self.ymax = ymax
[923]	254	self.proj = proj
	255	self.vmin = vmin ; self.vmax = vmax
	256	self.div = div
[1050]	257	self.colorbar = colorbar
	258	self.linestyle = linestyle
[923]	259	self.marker = marker
	260	self.color = color
[1050]	261	self.legend = legend
[960]	262	self.units = units
[923]	263	self.title = title
[1029]	264	self.xp = xp ; self.yp = yp
	265	self.nxticks = nxticks ; self.nyticks = nyticks
[1054]	266	self.fmt = fmt
[1062]	267	self.trans = trans ; self.back = back
	268	self.showcb = showcb
[1067]	269	self.logy = logy
[910]	270
	271	# print status
	272	def printstatus(self):
[960]	273	if not self.quiet:
[920]	274	if self.filename == "THIS_IS_A_CLONE":
	275	pass
	276	else:
[923]	277	print "**** PPCLASS. Done step: " + self.status
[910]	278
[1029]	279	# print attributes
	280	def printme(self):
	281	for k, v in vars(self).items():
	282	print k,v
	283
[910]	284	#####################################################
	285	# EMULATE OPERATORS + - * / ** << FOR PP() OBJECTS #
	286	#####################################################
	287
[923]	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
[1050]	300	self.sx = other.sx ; self.sy = other.sy
	301	self.sz = other.sz ; self.st = other.st
[1108]	302	self.compute = other.compute
	303	self.kind3d = other.kind3d
[923]	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
[1158]	316	self.xmin = other.xmin ; self.xmax = other.xmax
	317	self.ymin = other.ymin ; self.ymax = other.ymax
[923]	318	self.div = other.div
[1050]	319	self.colorbar = other.colorbar
	320	self.linestyle = other.linestyle
[923]	321	self.marker = other.marker
	322	self.color = other.color
[1050]	323	self.legend = other.legend
[960]	324	self.units = other.units
[923]	325	self.title = other.title
[930]	326	self.includedate = other.includedate
[934]	327	self.changetime = other.changetime
[990]	328	self.savtxt = other.savtxt
[991]	329	self.modx = other.modx
[1029]	330	self.xp = other.xp ; self.yp = other.yp
	331	self.missing = other.missing
	332	self.nxticks = other.nxticks ; self.nyticks = other.nyticks
[1054]	333	self.fmt = other.fmt
[1062]	334	self.trans = other.trans ; self.back = other.back
	335	self.showcb = other.showcb
[1067]	336	self.logy = other.logy
[923]	337	else:
	338	print "!! ERROR !! argument must be a pp object." ; exit()
	339
[910]	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
[914]	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
[920]	378	the_clone.filename = "THIS_IS_A_CLONE" # trick to avoid additional outputs
[914]	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):
[915]	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():
[914]	389	if k != "field":
[915]	390	setattr(obj,k,v)
[914]	391	the_clone.status = "retrieved"
[923]	392	the_clone.filename = self.filename
[914]	393	return the_clone
	394
[910]	395	# define the operation + on two objects. or with an int/float.
[914]	396	# ... with selective_copy the self object is not modified.
[910]	397	def __add__(self,other):
	398	isnum = self.checktwo(other)
[914]	399	the_clone = self.selective_copy()
[910]	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):
[914]	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
[960]	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
[914]	414	exit()
	415	else:
	416	ope = other
[923]	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
[910]	421	else:
[914]	422	obj.field = obj_ref.field + ope
	423	return the_clone
[910]	424
	425	# define the operation - on two objects. or with an int/float.
[914]	426	# ... with selective_copy the self object is not modified.
[910]	427	def __sub__(self,other):
	428	isnum = self.checktwo(other)
[914]	429	the_clone = self.selective_copy()
[910]	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):
[914]	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
[960]	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
[914]	444	exit()
	445	else:
	446	ope = other
[923]	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
[910]	451	else:
[914]	452	obj.field = obj_ref.field - ope
	453	return the_clone
[910]	454
	455	# define the operation * on two objects. or with an int/float.
[914]	456	# ... with selective_copy the self object is not modified.
[910]	457	def __mul__(self,other):
	458	isnum = self.checktwo(other)
[914]	459	the_clone = self.selective_copy()
[910]	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):
[914]	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
[960]	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
[914]	474	exit()
	475	else:
	476	ope = other
[923]	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
[910]	481	else:
[914]	482	obj.field = obj_ref.field * ope
	483	return the_clone
[910]	484
	485	# define the operation / on two objects. or with an int/float.
[914]	486	# ... with selective_copy the self object is not modified.
[910]	487	def __div__(self,other):
	488	isnum = self.checktwo(other)
[914]	489	the_clone = self.selective_copy()
[910]	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):
[914]	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
[960]	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
[914]	504	exit()
	505	else:
	506	ope = other
[923]	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
[910]	511	else:
[914]	512	obj.field = obj_ref.field / ope
	513	return the_clone
[910]	514
[914]	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
[910]	539	# define the operation ** on one object.
[914]	540	# ... with selective_copy the self object is not modified.
[910]	541	def __pow__(self,num):
[914]	542	the_clone = self.selective_copy()
[910]	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):
[914]	550	obj = the_clone.request[i][j][t][z][y][x]
	551	obj_ref = self.request[i][j][t][z][y][x]
[923]	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
[910]	556	else:
[914]	557	obj.field = obj_ref.field ** num
[910]	558	else:
	559	print "!! ERROR !! To define a power, either an int or a float is needed." ; exit()
[914]	560	return the_clone
[910]	561
[923]	562	### TBD: reverse power? for exponentials?
[910]	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
[932]	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
[910]	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)]
[960]	625	# store how many onerequest() objects are in self.request
	626	self.nrequest = self.nfinself.nvinself.nplotxself.nplotyself.nplotz*self.nplott
[910]	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]
[942]	639	# get methods
	640	obj.method_x = mx ; obj.method_y = my
	641	obj.method_z = mz ; obj.method_t = mt
[910]	642	# indicate the computation method
	643	obj.compute = self.compute
[1108]	644	# indicate the kind of 3D fields
	645	obj.kind3d = self.kind3d
[910]	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()
[934]	650	### possible time axis change
	651	obj.changetime = self.changetime
	652	obj.performtimechange()
[981]	653	# get strides
[1050]	654	obj.sx = self.sx ; obj.sy = self.sy
	655	obj.sz = self.sz ; obj.st = self.st
[910]	656	### get index
[981]	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)
[1029]	660	# missing value
	661	obj.missing = self.missing
[910]	662	# change status
	663	self.status = "defined"
[923]	664	return self
[910]	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
[920]	678	if self.status != "defined": print "!! ERROR !! Please use .define() to define your pp object." ; exit()
[960]	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) ]
[1201]	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) ]
[960]	688	count = 0
[910]	689	## first get fields
	690	## ... only what is needed is extracted from the files
[921]	691	## ... and computations are performed
[910]	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()
[1002]	700	if self.compute != "nothing":
	701	obj.computations()
[960]	702	# save fields in self.f for the user
	703	self.f[count] = obj.field
[1201]	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
[1050]	708	# save a legend in self.l for the user
[960]	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])+'_'
[983]	716	# close the file
	717	obj.closefile()
[960]	718	count = count + 1
	719	## make it simple: self.f is simply the data array if self.nrequest=1
[1192]	720	if self.nrequest == 1:
	721	self.f = self.f[0]
[1201]	722	self.xx = self.xx[0]
	723	self.yy = self.yy[0]
	724	self.zz = self.zz[0]
	725	self.tt = self.tt[0]
[910]	726	# change status
	727	self.status = "retrieved"
[923]	728	return self
[910]	729
	730	##########################################################
	731	# get: a shortcut method for the define + retrieve chain #
	732	##########################################################
	733	def get(self):
	734	self.define()
	735	self.retrieve()
[923]	736	return self
[910]	737
[960]	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	############################################################
[1192]	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()
[1201]	756	return self.f,self.xx,self.yy,self.zz,self.tt
[1192]	757
	758	############################################################
[960]	759	# getfl: a shortcut method for the define + retrieve chain #
	760	# ... in which the output is self.f, self.l #
[1192]	761	# ... and the ppclass is kept quiet #
[960]	762	############################################################
	763	def getfl(self):
	764	self.quiet = True
	765	self.get()
	766	return self.f,self.l
	767
[910]	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)
[1108]	790	#obj.field = ppcompute.smooth2diter(obj.field,n=window)
[910]	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
[1029]	812	## [BUG FIX: apparently info about missing values is not saved correctly]
	813	for count in range(self.nplot):
	814	pl = self.p[count]
[1050]	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
[1029]	817	return #self?
[910]	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
[1050]	872	plobj.f = obj.field # field to be plotted
[1219]	873	plobj.x = obj.absc # abscissa (or longitude)
	874	plobj.y = obj.ordi # ordinate (or latitude)
	875	# -- useless in 1D but not used anyway
[923]	876	# specific 1D plot stuff
	877	if dp == 1:
[1050]	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])
[923]	886	# specific 2d plot stuff
[910]	887	if dp == 2:
	888	# -- light grey background for missing values
[1050]	889	if type(plobj.f).__name__ in 'MaskedArray': plobj.axisbg = '0.75'
[910]	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 )
[923]	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
[960]	901	if self.units is not None: plobj.units = self.units
[1050]	902	if self.colorbar is not None: plobj.colorbar = self.colorbar
[991]	903	if self.modx is not None: plobj.modx = self.modx
[1029]	904	if self.nxticks is not None: plobj.nxticks = self.nxticks
	905	if self.nyticks is not None: plobj.nyticks = self.nyticks
[1054]	906	if self.fmt is not None: plobj.fmt = self.fmt
[1062]	907	if self.showcb is not None: plobj.showcb = self.showcb
[1067]	908	if self.logy is not None: plobj.logy = self.logy
[1158]	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
[923]	913	# -- 1D specific
	914	if dp == 1:
[1050]	915	if self.linestyle is not None: plobj.linestyle = self.linestyle
[923]	916	if self.marker is not None: plobj.marker = self.marker
	917	if self.color is not None: plobj.color = self.color
[1050]	918	if self.legend is not None: plobj.legend = self.legend
[923]	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
[1062]	924	if self.trans is not None: plobj.trans = self.trans
	925	if self.back is not None: plobj.back = self.back
[1050]	926	plobj.svx = self.svx
	927	plobj.svy = self.svy
[923]	928	# finally append plot object
[910]	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
[1116]	948	plobj = self.p[count]
[910]	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":
[1050]	961	plobj.c = self.request[i][j+1][t][z][y][x].field ; found += 1
[910]	962	if condfile == "contour":
[1050]	963	plobj.c = self.request[i+1][j][t][z][y][x].field ; found += 1
[910]	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:
[1050]	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
[910]	975	if "vector" in condfile:
[1050]	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
[910]	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:
[960]	985	pass # because this means that we only had 0D values !
[910]	986	# final status
	987	self.status = "definedplot"
[923]	988	return self
[910]	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):
[960]	999	if self.howmanyplots > 0:
[910]	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
[1029]	1009	self.fig = ppplot.figuref(x=self.xp,y=self.yp)
[910]	1010	self.subv,self.subh = ppplot.definesubplot(self.howmanyplots,self.fig)
	1011	self.n = 0
[917]	1012	## adapted space for labels etc
	1013	## ... except for ortho because there is no label anyway
[1050]	1014	self.customplot = self.p[0].f.ndim == 2 \
[917]	1015	and self.p[0].mapmode == True \
	1016	and self.p[0].proj not in ["ortho"]
[921]	1017	if self.customplot:
[917]	1018	margin = 0.07
	1019	self.fig.subplots_adjust(left=margin,right=1-margin,bottom=margin,top=1-margin)
[910]	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
[921]	1027	self.customplot = self.plotin.customplot
[910]	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)
[1170]	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
[977]	1046	# possibility to color lines according to a color map
	1047	# ... made so that all plots span the whole color map automatically.
[1050]	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)
[910]	1051	else:
[1050]	1052	pl.invert = False ; pl.linestyle = None # don't invert again axis
[914]	1053	# set saved titles and labels
	1054	if self.plotin is None:
[1158]	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
[914]	1061	else:
	1062	prev_plot = self.plotin.p[self.n-1]
[1158]	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
[910]	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()
[990]	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
[1050]	1078	ppplot.writeascii(field=pl.f,absc=pl.x,name=name)
[923]	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
[910]	1082	# once completed show the plot (cannot show intermediate plotin)
[917]	1083	# ... added a fix (customplot=True) for the label problem in basemap
[977]	1084	if not self.quiet: print "**** PPCLASS. Done step: makeplot"
[923]	1085	if (self.n == self.howmanyplots):
[964]	1086	ppplot.save(mode=self.out,filename=self.filename,folder=self.folder,custom=self.customplot,includedate=self.includedate,res=self.res)
[920]	1087	mpl.close()
[910]	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"
[920]	1091	try:
	1092	filehandler = open(savfile, 'w')
	1093	pickle.dump(self.p, filehandler)
	1094	except IOError:
[960]	1095	if self.verbose: print "!! WARNING !! Saved object file not written. Probably do not have permission to write here."
[923]	1096	return self
[910]	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()
[923]	1104	return self
[910]	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)
[923]	1112	return self
[910]	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)
[923]	1120	return self
[910]	1121
[960]	1122	#################################################################
	1123	# func: operation on two pp objects being on status 'definedplot'
[910]	1124	# this allows for one field being function of another one
[960]	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):
[910]	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]
[1050]	1144	if sobj.f.ndim !=1 or oobj.f.ndim !=1:
[910]	1145	if self.verbose: print "!! WARNING !! Flattening arrays because more than one-dimensional."
[1050]	1146	sobj.f = np.ravel(sobj.f)
	1147	oobj.f = np.ravel(oobj.f)
	1148	sobj.x = oobj.f
[910]	1149	sobj.xlabel = oobj.ylabel
[1050]	1150	if sobj.x.size > sobj.f.size:
[910]	1151	if self.verbose:
[1050]	1152	print "!! WARNING !! Trying to define y=f(x) with x and y not at the same size.",sobj.x.size,sobj.f.size
[910]	1153	print "!! WARNING !! Modifying x to fit y size but please check."
[1050]	1154	sobj.x = sobj.x[0:sobj.f.size]
[910]	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):
[1092]	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	##
[972]	1175	if opt.void:
	1176	self.p[iii].showcb = False
	1177	else:
	1178	self.p[iii].showcb = True
[910]	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	###
[1050]	1189	try: self.p[iii].colorbar = opt.colorbar[iii]
[910]	1190	except:
[1050]	1191	try: self.p[iii].colorbar = opt.colorbar[0] ; self.colorbar = opt.colorbar[0]
[910]	1192	except: pass
	1193	###
[972]	1194	if opt.void:
	1195	self.p[iii].title = ""
	1196	else:
	1197	try: self.p[iii].title = opt.title[iii]
	1198	except:
[910]	1199	try: self.p[iii].title = opt.title[0]
	1200	except: pass
	1201	###
[972]	1202	if opt.void:
	1203	self.p[iii].xlabel = ""
	1204	else:
	1205	try: self.p[iii].xlabel = opt.xlabel[iii]
	1206	except:
[910]	1207	try: self.p[iii].xlabel = opt.xlabel[0]
	1208	except: pass
	1209	###
[972]	1210	if opt.void:
	1211	self.p[iii].ylabel = ""
	1212	else:
	1213	try: self.p[iii].ylabel = opt.ylabel[iii]
	1214	except:
[910]	1215	try: self.p[iii].ylabel = opt.ylabel[0]
	1216	except: pass
	1217	###
[1050]	1218	try: self.p[iii].linestyle = opt.linestyle[iii]
[910]	1219	except:
[1050]	1220	try: self.p[iii].linestyle = opt.linestyle[0]
[910]	1221	except: pass
	1222	###
[920]	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	###
[1050]	1233	try: self.p[iii].legend = opt.legend[iii]
[923]	1234	except:
[1050]	1235	try: self.p[iii].legend = opt.legend[0]
[923]	1236	except: pass
	1237	###
[910]	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
[963]	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
[972]	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
[977]	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
[991]	1312	###
[1084]	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	###
[991]	1328	try: self.p[iii].modx = opt.modx[iii]
	1329	except:
	1330	try: self.p[iii].modx = opt.modx[0]
	1331	except: pass
[1054]	1332	###
	1333	try: self.p[iii].fmt = opt.fmt[iii]
	1334	except:
	1335	try: self.p[iii].fmt = opt.fmt[0]
	1336	except: pass
[910]	1337
[963]	1338
[910]	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
[1068]	1361	self.tabtime = None
[910]	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
[934]	1368	self.changetime = None
[1050]	1369	self.sx = 1 ; self.sy = 1 ; self.sz = 1 ; self.st = 1
[1029]	1370	self.missing = '!! missing value: I am not set, damned !!'
[1108]	1371	self.kind3d = '!! kind3d: I am not set, damned !!'
[910]	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
[983]	1385	# close a file
	1386	# ------------
	1387	def closefile(self):
	1388	self.f.close()
	1389
[921]	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
[910]	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."
[1108]	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()
[910]	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."
[1215]	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.
[910]	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:
[1090]	1461	if np.all(self.field_x == self.field_x[0,0]) \
[1002]	1462	or self.field_x.min() == self.field_x.max() \
	1463	or self.field_y.min() == self.field_y.max():
[960]	1464	if self.verbose: print "!! WARNING !! xy axis look undefined. creating non-dummy ones."
[910]	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.
[1029]	1473	# WARNING: how do we do if several are available? the last one is chosen.
[910]	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:
[990]	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
[1029]	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"
[910]	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())
[1029]	1500
[910]	1501	# TIME. Try preset fields.
	1502	self.name_t = "nothing"
	1503	for c in glob_listt:
[936]	1504	if c in self.f.variables.keys():
[910]	1505	self.name_t = c
[936]	1506	if self.verbose: print "**** OK. Found time variable: ",c
[910]	1507	try:
[911]	1508	# speed up: only get first value, last one.
[1068]	1509	self.tabtime = self.f.variables[self.name_t]
[936]	1510	# (consider the case where tabtime is not dim 1)
	1511	# (time is most often the first dimension)
[1068]	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]
[972]	1515	# (now proceed) (the +0. just ensures this is a number)
[1068]	1516	dafirst = self.tabtime[0] + 0.
[936]	1517	if self.dim_t == 1:
	1518	self.field_t = np.array([dafirst])
	1519	else:
[1068]	1520	daint = self.tabtime[1] - dafirst
[1029]	1521	dalast = dafirst + (self.dim_t-1)*daint
[936]	1522	self.field_t = np.linspace(dafirst,dalast,num=self.dim_t)
[1029]	1523	if self.verbose:
[1031]	1524	print "!! WARNING !! WARNING !! Time axis is supposed to be equally spaced !!"
[1068]	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]
[910]	1527	except:
	1528	# ... or if a problem encountered, define a simple time axis
[972]	1529	if self.verbose: print "**** OK. There is something weird. Let us go for a simple time axis."
[910]	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
[934]	1535	# change time axis
	1536	# ... add your options here!
	1537	# --------------------------
	1538	def performtimechange(self):
[942]	1539	if self.changetime is not None:
[934]	1540	if self.verbose: print "**** OK. Converting time axis:",self.changetime
[1029]	1541	### options added by T. Navarro
[1015]	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]
[934]	1547	self.field_t = ppcompute.mars_sol2ls(self.field_t)
[1015]	1548	elif self.changetime == "mars_dayini" and "controle" in self.f.variables:
[1003]	1549	self.field_t = self.field_t \
	1550	+ self.f.variables['controle'][3]%669 \
	1551	+ self.f.variables['controle'][26]
[942]	1552	### options added by A. Spiga
[1029]	1553	elif self.changetime == "correctls":
[1163]	1554	if self.tabtime is None:
	1555	print "!! WARNING !! Encountered a problem with correctls. Check your time dimension. Skipping this part."
	1556	else:
[1069]	1557	dafirst = self.tabtime[0] + 0.
	1558	if self.dim_t == 1:
	1559	self.field_t = np.array([dafirst])
	1560	else:
[1068]	1561	daint = self.tabtime[1] - dafirst
[1029]	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):
[1068]	1566	if self.tabtime[iii] - self.tabtime[iii-1] < 0: year = year+1.
[1029]	1567	add[iii] = year*360.
[1068]	1568	self.field_t = add + self.tabtime
[942]	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!
[934]	1598	else:
[960]	1599	if self.verbose: print "!! WARNING !! This time change is not implemented. Nothing is done."
[942]	1600	if self.verbose: print "**** OK. new t axis values [%5.1f,%5.1f]" % (self.field_t.min(),self.field_t.max())
[934]	1601
[910]	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	# -------------------------------
[981]	1605	def getindextime(self,dalist=None,ind=None):
[910]	1606	if self.method_t == "free":
[1050]	1607	self.index_t = np.arange(0,self.dim_t,self.st)
[910]	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] ) )
[1178]	1617	self.index_t = np.arange(start,stop+1,self.st)
[910]	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	# -------------------------------
[981]	1629	def getindexvert(self,dalist=None,ind=None):
[910]	1630	if self.method_z == "free":
[1050]	1631	self.index_z = np.arange(0,self.dim_z,self.sz)
[910]	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] ) )
[1178]	1641	self.index_z = np.arange(start,stop+1,self.sz)
[910]	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	# -------------------------------
[981]	1657	def getindexhori(self,dalistx=None,dalisty=None,indx=None,indy=None):
[910]	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"]:
[1050]	1666	self.index_x = range(0,self.dim_x,self.sx)
[910]	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."
[916]	1673	if self.method_y == "free" or test in ["freecomp","compcomp"]:
[1050]	1674	self.index_y = range(0,self.dim_y,self.sy)
[910]	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."
[926]	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!
[910]	1684	if self.method_x in ["free","comp"] and self.method_y in ["free","comp"]:
[926]	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":
[1219]	1691	yeah = (self.field_y >= dalisty[indy][0])*(self.field_y <= dalisty[indy][1])
[926]	1692	self.index_y = yeah[:,0]
[910]	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
[960]	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."
[1219]	1758	# well, now get field from netcdf file
[910]	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]
[1219]	1764	# set everything about dimensions and slices
	1765	# -- each self.dim case corresponds to tests in the beginning of getdim.
[910]	1766	if self.dim == 1:
	1767	nt = self.index_t.size ; nz = 1 ; ny = 1 ; nx = 1
[1219]	1768	tupledim = self.index_t
[910]	1769	elif self.dim == 2:
	1770	nt = 1 ; nz = 1 ; ny = self.index_y.size ; nx = self.index_x.size
[1219]	1771	tupledim = self.index_y,self.index_x
[910]	1772	elif self.dim == 3:
[1108]	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
[1219]	1776	tupledim = self.index_t,self.index_y,self.index_x
[1108]	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
[1219]	1780	tupledim = self.index_t,self.index_z,self.index_y
[1108]	1781	else:
	1782	print "!! ERROR !! This kind of 3D field is not supported. Please send feedback." ; exit()
[910]	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
[1219]	1786	tupledim = self.index_t,self.index_z,self.index_y,self.index_x
[910]	1787	else:
	1788	print "!! ERROR !! field would have more than four dimensions ?" ; exit()
[1219]	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]
[932]	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...
[933]	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
[910]	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
[932]	1803	self.field = np.reshape(self.field,(nt,nz,ny,nx))
[910]	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]
[981]	1816	# ... there are special cases with strides
	1817	# ... some other fixes will be needed probably TBD
[1050]	1818	if self.method_x == "free" and self.sx != 1:
[981]	1819	self.field_x = self.field_x[self.index_x]
[1050]	1820	if self.method_y == "free" and self.sy != 1:
[981]	1821	self.field_y = self.field_y[self.index_y]
[910]	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"]:
[1050]	1830	if self.sx == 1 and self.sy == 1:
[930]	1831	time0 = timelib.time()
[925]	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
[910]	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]
[925]	1856	for iz in range(self.index_z.size):
	1857	for it in range(self.index_t.size):
[910]	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))
[981]	1864	else:
	1865	# there is a problem above if stride != 1. a fix must be found. rewrite might be necessary. TBD
	1866	pass
[1089]	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?)
[1094]	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
[1089]	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.
[910]	1878	# ... this is important for computations below (see ppcompute)
[1089]	1879	masked = np.ma.masked_where(np.abs(self.field) >= self.missing,self.field)
[910]	1880	if masked.mask.any() == True:
[1029]	1881	if self.verbose: print "!! WARNING !! Values over %5.3e are considered missing values." % self.missing
[910]	1882	self.field = masked
	1883	self.field.set_fill_value([np.NaN])
[921]	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
[923]	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)."
[1011]	1899	# prepare quadratic mean
	1900	if "qmean" in self.compute: self.field = self.field*self.field
[1219]	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
[910]	1933	# remove all dimensions with size 1 to prepare plot (and check the resulting dimension with dimplot)
	1934	self.field = np.squeeze(self.field)
[1011]	1935	# take root mean square for quadratic mean
	1936	if self.compute == "qmean": self.field = np.sqrt(self.field)
	1937	# error handling and verbose
[910]	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
[921]	1942
[923]	1943	# get areas for computations and ponderate self.field by area/totarea
	1944	# -------------------------------------------------------------------
[921]	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]
[923]	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":
[930]	1986	aire.field = aire.field[np.ix_(self.index_y, self.index_x)] # reduce to requested indexes only
[923]	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
[921]	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:
[930]	1996	print "!! ERROR !! Problem in area(). Check array shapes."
	1997	print "Field vs. aire:",self.field.shape,aire.field.shape ; exit()
[921]	1998	else:
	1999	self.field = self.field*aire.field
[910]	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: