Context Navigation

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

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