Context Navigation

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

Last change on this file since 920 was 920, checked in by aslmd, 12 years ago
UTIL PYTHON planetoplot_v2. Added options to pp.py. Added help to pp.py. Set warnings in case files cannot be written. Set warnings for bounds too tight. Corrected one possible bug when several plots in a row.
File size: 69.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	print "************************************"
	27	print "** WELCOME TO PLANETOPLOT 2.0 **"
	28	print "************************************"
	29	# where settings files are located...
	30	whereset = None
	31	whereset = ppcompute.findset(whereset)
	32	# ... we load user-defined automatic settings from set_ppclass.txt
	33	zefile = "set_ppclass.txt"
	34	glob_listx = [] ; glob_listy = [] ; glob_listz = [] ; glob_listt = []
	35	try:
	36	f = open(whereset+zefile, 'r') ; lines = f.readlines()
	37	for stuff in lines[5].strip().split(';'): glob_listx.append(stuff)
	38	for stuff in lines[8].strip().split(';'): glob_listy.append(stuff)
	39	for stuff in lines[11].strip().split(';'): glob_listz.append(stuff)
	40	for stuff in lines[14].strip().split(';'): glob_listt.append(stuff)
	41	except IOError:
	42	print "warning: "+zefile+" not in "+whereset+" ; no presets."
	43
	44	##################################
	45	#### USEFUL GENERIC FUNCTIONS ####
	46	##################################
	47
	48	# inspect variables and dimensions in a netCDF file
	49	def inspect(filename):
	50	print "**** INSPECT FILE",filename
	51	test = netCDF4.Dataset(filename)
	52	print "**** VARIABLES: ",test.variables.keys()
	53	for dim in test.dimensions.keys():
	54	output = "**** DIMENSION: "+str(dim)+" "+str(len(test.dimensions[dim]))
	55	try: output = output + " ----> "+str(test.variables[dim][0])+" "+str(test.variables[dim][-1])
	56	except: pass
	57	print output ; output = ""
	58
	59	# 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
	103	# (needed by readslices)
	104	# look for comas in the input name to separate different names (files, variables,etc ..)
	105	def separatenames (name):
	106	if name is None: names = None
	107	else:
	108	names = [] ; stop = 0 ; currentname = name
	109	while stop == 0:
	110	indexvir = currentname.find(',')
	111	if indexvir == -1: stop = 1 ; name1 = currentname
	112	else: name1 = currentname[0:indexvir]
	113	names = np.concatenate((names,[name1]))
	114	currentname = currentname[indexvir+1:len(currentname)]
	115	return names
	116
	117	#######################
	118	### THE MAIN OBJECT ###
	119	#######################
	120	class pp():
	121
	122	# print out a help string when help is invoked on the object
	123	def __repr__(self):
	124	whatprint = 'pp object. \"help(pp)\" for more information\n'
	125	return whatprint
	126
	127	# default settings
	128	# -- user can define settings by two methods.
	129	# -- 1. yeah = pp(file="file.nc")
	130	# -- 2. yeah = pp() ; yeah.file = "file.nc"
	131	def __init__(self,file=None,var="notset",\
	132	filegoal=None,vargoal=None,\
	133	x=None,y=None,z=None,t=None,\
	134	stridex=1,stridey=1,\
	135	stridez=1,stridet=1,\
	136	compute="mean",\
	137	verbose=False,noproj=False,\
	138	superpose=False,\
	139	plotin=None,\
	140	forcedimplot=-1,\
	141	out="gui",\
	142	filename="myplot",\
	143	folder="./"):
	144	self.request = None
	145	self.nfin = 0 ; self.nvin = 0
	146	self.nplotx = None ; self.nploty = None
	147	self.nplotz = None ; self.nplott = None
	148	self.status = "init"
	149	self.fig = None ; self.subv = None ; self.subh = None
	150	self.n = 0 ; self.howmanyplots = 0
	151	self.nplot = 0
	152	self.p = None
	153	## what could be defined by the user
	154	self.file = file
	155	self.var = var
	156	self.filegoal = filegoal
	157	self.vargoal = vargoal
	158	self.x = x ; self.y = y ## if None, free dimension
	159	self.z = z ; self.t = t ## if None, free dimension
	160	self.stridex = stridex ; self.stridey = stridey
	161	self.stridez = stridez ; self.stridet = stridet
	162	self.compute = compute
	163	self.verbose = verbose
	164	self.noproj = noproj
	165	self.plotin = plotin
	166	self.superpose = superpose
	167	self.forcedimplot = forcedimplot
	168	self.out = out
	169	self.filename = filename
	170	self.folder = folder
	171
	172	# print status
	173	def printstatus(self):
[920]	174	if self.filename == "THIS_IS_A_CLONE":
	175	pass
	176	else:
	177	print "**** Done step: " + self.status
[910]	178
	179	#####################################################
	180	# EMULATE OPERATORS + - * / ** << FOR PP() OBJECTS #
	181	#####################################################
	182
	183	# check the compatibility of two objects for operations
	184	# --> if other is a pp class, test sizes and return isnum = False
	185	# --> if other is an int or a float, return isnum = True
	186	# --> otherwise, just print an error and exit
	187	def checktwo(self,other):
	188	if other.__class__.__name__ == "pp":
	189	isnum = False
	190	if self.status in ["init","defined"] or other.status in ["init","define"]:
	191	print "!! ERROR !! Please use .retrieve to get fields for plots with one of your pp operands." ; exit()
	192	if self.nfin != other.nfin or \
	193	self.nvin != other.nvin or \
	194	self.nplott != other.nplott or \
	195	self.nplotz != other.nploty or \
	196	self.nploty != other.nploty or \
	197	self.nplotx != other.nplotx :
	198	print "!! ERROR !! The two operands do not have the same number of files, variables, t z y x requests."
	199	exit()
	200	elif isinstance(other,int) or isinstance(other,float):
	201	isnum = True
	202	else:
	203	print "!! ERROR !! The operand is neither a pp class nor an integer or a float." ; exit()
	204	return isnum
	205
[914]	206	# define a selective copy of a pp() object for operations
	207	# ... copy.copy() is not conservative (still acts like a pointer)
	208	# ... copy.deepcopy() does not work with netCDF objects
	209	# so what is done here is a copy of everything except
	210	# (to avoid sharing with self and therefore modifying self through operations)
	211	# - request attribute of pp() object
	212	# - field attribute of the onerequest() objects
	213	def selective_copy(self):
	214	if self.status in ["init","defined"]:
	215	print "!! ERROR !! Please use .retrieve to get fields for the object you want to copy from." ; exit()
	216	the_clone = pp()
	217	for k, v in vars(self).items():
	218	if k != "request":
	219	setattr(the_clone,k,v)
	220	the_clone.verbose = False
[920]	221	the_clone.filename = "THIS_IS_A_CLONE" # trick to avoid additional outputs
[914]	222	the_clone.define()
	223	for i in range(self.nfin):
	224	for j in range(self.nvin):
	225	for t in range(self.nplott):
	226	for z in range(self.nplotz):
	227	for y in range(self.nploty):
	228	for x in range(self.nplotx):
[915]	229	obj_ref = self.request[i][j][t][z][y][x]
	230	obj = the_clone.request[i][j][t][z][y][x]
	231	for k, v in vars(obj_ref).items():
[914]	232	if k != "field":
[915]	233	setattr(obj,k,v)
[914]	234	the_clone.status = "retrieved"
	235	return the_clone
	236
[910]	237	# define the operation + on two objects. or with an int/float.
[914]	238	# ... with selective_copy the self object is not modified.
[910]	239	def __add__(self,other):
	240	isnum = self.checktwo(other)
[914]	241	the_clone = self.selective_copy()
[910]	242	for i in range(self.nfin):
	243	for j in range(self.nvin):
	244	for t in range(self.nplott):
	245	for z in range(self.nplotz):
	246	for y in range(self.nploty):
	247	for x in range(self.nplotx):
[914]	248	obj = the_clone.request[i][j][t][z][y][x]
	249	obj_ref = self.request[i][j][t][z][y][x]
	250	if not isnum:
	251	ope = other.request[i][j][t][z][y][x].field
	252	if obj_ref.field.shape != ope.shape:
	253	print "!! ERROR !! The two fields for operation do not have the same shape.",self.field.shape,other.field.shape
	254	exit()
	255	else:
	256	ope = other
[910]	257	if "vector" in self.vargoal[j] + self.filegoal[i]:
[914]	258	print "!! ERROR !! we do not operate on vectors yet."
	259	exit()
[910]	260	else:
[914]	261	obj.field = obj_ref.field + ope
	262	return the_clone
[910]	263
	264	# define the operation - on two objects. or with an int/float.
[914]	265	# ... with selective_copy the self object is not modified.
[910]	266	def __sub__(self,other):
	267	isnum = self.checktwo(other)
[914]	268	the_clone = self.selective_copy()
[910]	269	for i in range(self.nfin):
	270	for j in range(self.nvin):
	271	for t in range(self.nplott):
	272	for z in range(self.nplotz):
	273	for y in range(self.nploty):
	274	for x in range(self.nplotx):
[914]	275	obj = the_clone.request[i][j][t][z][y][x]
	276	obj_ref = self.request[i][j][t][z][y][x]
	277	if not isnum:
	278	ope = other.request[i][j][t][z][y][x].field
	279	if obj_ref.field.shape != ope.shape:
	280	print "!! ERROR !! The two fields for operation do not have the same shape.",self.field.shape,other.field.shape
	281	exit()
	282	else:
	283	ope = other
[910]	284	if "vector" in self.vargoal[j] + self.filegoal[i]:
[914]	285	print "!! ERROR !! we do not operate on vectors yet."
	286	exit()
[910]	287	else:
[914]	288	obj.field = obj_ref.field - ope
	289	return the_clone
[910]	290
	291	# define the operation * on two objects. or with an int/float.
[914]	292	# ... with selective_copy the self object is not modified.
[910]	293	def __mul__(self,other):
	294	isnum = self.checktwo(other)
[914]	295	the_clone = self.selective_copy()
[910]	296	for i in range(self.nfin):
	297	for j in range(self.nvin):
	298	for t in range(self.nplott):
	299	for z in range(self.nplotz):
	300	for y in range(self.nploty):
	301	for x in range(self.nplotx):
[914]	302	obj = the_clone.request[i][j][t][z][y][x]
	303	obj_ref = self.request[i][j][t][z][y][x]
	304	if not isnum:
	305	ope = other.request[i][j][t][z][y][x].field
	306	if obj_ref.field.shape != ope.shape:
	307	print "!! ERROR !! The two fields for operation do not have the same shape.",self.field.shape,other.field.shape
	308	exit()
	309	else:
	310	ope = other
[910]	311	if "vector" in self.vargoal[j] + self.filegoal[i]:
[914]	312	print "!! ERROR !! we do not operate on vectors yet."
	313	exit()
[910]	314	else:
[914]	315	obj.field = obj_ref.field * ope
	316	return the_clone
[910]	317
	318	# define the operation / on two objects. or with an int/float.
[914]	319	# ... with selective_copy the self object is not modified.
[910]	320	def __div__(self,other):
	321	isnum = self.checktwo(other)
[914]	322	the_clone = self.selective_copy()
[910]	323	for i in range(self.nfin):
	324	for j in range(self.nvin):
	325	for t in range(self.nplott):
	326	for z in range(self.nplotz):
	327	for y in range(self.nploty):
	328	for x in range(self.nplotx):
[914]	329	obj = the_clone.request[i][j][t][z][y][x]
	330	obj_ref = self.request[i][j][t][z][y][x]
	331	if not isnum:
	332	ope = other.request[i][j][t][z][y][x].field
	333	if obj_ref.field.shape != ope.shape:
	334	print "!! ERROR !! The two fields for operation do not have the same shape.",self.field.shape,other.field.shape
	335	exit()
	336	else:
	337	ope = other
[910]	338	if "vector" in self.vargoal[j] + self.filegoal[i]:
[914]	339	print "!! ERROR !! we do not operate on vectors yet."
	340	exit()
[910]	341	else:
[914]	342	obj.field = obj_ref.field / ope
	343	return the_clone
[910]	344
[914]	345	# define the reverse operation float/int + object
	346	def __radd__(self,other):
	347	isnum = self.checktwo(other)
	348	if not isnum: print "!! ERROR !! Operand should be a number" ; exit()
	349	return self.__add__(other)
	350
	351	# define the reverse operation float/int - object
	352	def __rsub__(self,other):
	353	isnum = self.checktwo(other)
	354	if not isnum: print "!! ERROR !! Operand should be a number" ; exit()
	355	return self.__sub__(other)
	356
	357	# define the reverse operation float/int * object
	358	def __rmul__(self,other):
	359	isnum = self.checktwo(other)
	360	if not isnum: print "!! ERROR !! Operand should be a number" ; exit()
	361	return self.__mul__(other)
	362
	363	# define the reverse operation float/int / object
	364	def __rdiv__(self,other):
	365	isnum = self.checktwo(other)
	366	if not isnum: print "!! ERROR !! Operand should be a number" ; exit()
	367	return self.__div__(other)
	368
[910]	369	# define the operation ** on one object.
[914]	370	# ... with selective_copy the self object is not modified.
[910]	371	def __pow__(self,num):
[914]	372	the_clone = self.selective_copy()
[910]	373	if isinstance(num,int) or isinstance(num,float):
	374	for i in range(self.nfin):
	375	for j in range(self.nvin):
	376	for t in range(self.nplott):
	377	for z in range(self.nplotz):
	378	for y in range(self.nploty):
	379	for x in range(self.nplotx):
[914]	380	obj = the_clone.request[i][j][t][z][y][x]
	381	obj_ref = self.request[i][j][t][z][y][x]
[910]	382	if "vector" in self.vargoal[j] + self.filegoal[i]:
[914]	383	print "!! ERROR !! we do not operate on vectors yet."
	384	exit()
[910]	385	else:
[914]	386	obj.field = obj_ref.field ** num
[910]	387	else:
	388	print "!! ERROR !! To define a power, either an int or a float is needed." ; exit()
[914]	389	return the_clone
[910]	390
	391	# define the operation <<
	392	# ... e.g. obj2 << obj1
	393	# ... means: get init for pp object obj2 from another pp object obj1
[914]	394	# ... (this should solve the affectation trap obj2 = obj1)
[910]	395	def __lshift__(self,other):
	396	if other.__class__.__name__ == "pp":
	397	self.file = other.file
	398	self.var = other.var
	399	self.filegoal = other.filegoal
	400	self.vargoal = other.vargoal
	401	self.x = other.x ; self.y = other.y ## if None, free dimension
	402	self.z = other.z ; self.t = other.t ## if None, free dimension
	403	self.stridex = other.stridex ; self.stridey = other.stridey
	404	self.stridez = other.stridez ; self.stridet = other.stridet
	405	self.verbose = other.verbose
	406	self.noproj = other.noproj
	407	self.plotin = other.plotin
	408	self.superpose = other.superpose
	409	self.forcedimplot = other.forcedimplot
	410	self.out = other.out
	411	self.filename = other.filename
	412	self.folder = other.folder
	413	else:
	414	print "!! ERROR !! argument must be a pp object." ; exit()
	415
	416	##############################################################################################
	417	# define method
	418	# ---------
	419	# ... (file and var are either one string or a vector of strings)
	420	# ... the goal of define is to define a 2D array of onerequest() objects (see class below)
	421	# given the number of file, var, x, y, z, t asked by the user
	422	# ... objectives for file or var are given through filegoal and vargoal
	423	# --> possible values: main contour vector
	424	# ---------
	425	# ... then onerequest() objects are being defined more precisely
	426	# by getting index_x index_y index_z index_t
	427	# and setting method_x method_y method_z method_t to either
	428	# - "free" for free dimensions (plot dimensions)
	429	# - "comp" for averages, max, min
	430	# - "fixed" for fixed dimensions (possibly several i.e. multislice)
	431	##############################################################################################
	432	def define(self):
	433	self.printstatus()
	434	# initial check and get dimensions
	435	self.file = checktab(self.file,mess="file")
	436	self.nfin = len(self.file)
	437	if self.verbose:
	438	for i in range(self.nfin): inspect(self.file[i])
	439	self.var = checktab(self.var,mess="var")
	440	self.nvin = len(self.var)
	441	# check goal tabs for files and variables
	442	# ... default is to plot everything
	443	if self.filegoal is None: self.filegoal = ["main"]*self.nfin
	444	if self.vargoal is None: self.vargoal = ["main"]*self.nvin
	445	self.filegoal = checktab(self.filegoal, mess="filegoal")
	446	self.vargoal = checktab(self.vargoal, mess="vargoal")
	447	if len(self.filegoal) != self.nfin: print "!! ERROR !! filegoal must be the same size as file." ; exit()
	448	if len(self.vargoal) != self.nvin: print "!! ERROR !! vargoal must be the same size as var." ; exit()
	449	# variables: initial check
	450	self.x = checktab(self.x,mess="x",allownone=True,allownumber=True)
	451	self.y = checktab(self.y,mess="y",allownone=True,allownumber=True)
	452	self.z = checktab(self.z,mess="z",allownone=True,allownumber=True)
	453	self.t = checktab(self.t,mess="t",allownone=True,allownumber=True)
	454	# for the moment not var- nor file- dependent.
	455	# but this could be the case.
	456	sx = readslices(self.x) ; sy = readslices(self.y)
	457	sz = readslices(self.z) ; st = readslices(self.t)
	458	# get methods
	459	mx = findmethod(sx) ; my = findmethod(sy)
	460	mz = findmethod(sz) ; mt = findmethod(st)
	461	# get number of plots to be done
	462	if mx == "fixed": self.nplotx = sx.size/2
	463	else: self.nplotx = 1
	464	if my == "fixed": self.nploty = sy.size/2
	465	else: self.nploty = 1
	466	if mz == "fixed": self.nplotz = sz.size/2
	467	else: self.nplotz = 1
	468	if mt == "fixed": self.nplott = st.size/2
	469	else: self.nplott = 1
	470	if self.verbose: print "**** OK. Plots over x,y,z,t -->",self.nplotx,self.nploty,self.nplotz,self.nplott
	471	# create the list of onerequest() objects
	472	self.request = [[[[[[ \
	473	onerequest() \
	474	for x in range(self.nplotx)] for y in range(self.nploty)] \
	475	for z in range(self.nplotz)] for t in range(self.nplott)] \
	476	for j in range(self.nvin)] for i in range(self.nfin)]
	477	# loop on onerequest() objects
	478	for i in range(self.nfin):
	479	for j in range(self.nvin):
	480	for t in range(self.nplott):
	481	for z in range(self.nplotz):
	482	for y in range(self.nploty):
	483	for x in range(self.nplotx):
	484	obj = self.request[i][j][t][z][y][x]
	485	# fill in names for files and variables
	486	obj.verbose = self.verbose
	487	obj.file = self.file[i]
	488	obj.var = self.var[j]
	489	# indicate the computation method
	490	obj.compute = self.compute
	491	# open the files (the same file might be opened several times but this is cheap)
	492	obj.openfile()
	493	### get x,y,z,t dimensions from file
	494	obj.getdim()
	495	### get methods
	496	obj.method_x = mx ; obj.method_y = my
	497	obj.method_z = mz ; obj.method_t = mt
	498	### get index
	499	obj.getindextime(st,t,self.stridet)
	500	obj.getindexvert(sz,z,self.stridez)
	501	obj.getindexhori(sx,sy,x,y,self.stridex,self.stridey)
	502	# change status
	503	self.status = "defined"
	504
	505	##############################################################################################
	506	# retrieve method
	507	# --> for each element onerequest() in the array, get field .var from .f file
	508	# --> see below the onerequest() class:
	509	# - only get what is needed for computing and plotting
	510	# - averages etc... are computed here
	511	# --> RESULT: each onerequest() object has now its attribute .field filled
	512	# --> if one wants to perform operations on fields, this should be done after retrieve()
	513	##############################################################################################
	514	def retrieve(self):
	515	self.printstatus()
	516	# check if things were done OK before
[920]	517	if self.status != "defined": print "!! ERROR !! Please use .define() to define your pp object." ; exit()
[910]	518	## first get fields
	519	## ... only what is needed is extracted from the files
	520	## ... and averages are computed
	521	for i in range(self.nfin):
	522	for j in range(self.nvin):
	523	for t in range(self.nplott):
	524	for z in range(self.nplotz):
	525	for y in range(self.nploty):
	526	for x in range(self.nplotx):
	527	obj = self.request[i][j][t][z][y][x]
	528	obj.getfield()
	529	# change status
	530	self.status = "retrieved"
	531
	532	##########################################################
	533	# get: a shortcut method for the define + retrieve chain #
	534	##########################################################
	535	def get(self):
	536	self.define()
	537	self.retrieve()
	538
	539	########################################
	540	# smooth: smooth the field in 1D or 2D #
	541	########################################
	542	## TBD: smooth not OK with masked array in the end of retrieve()
	543	def smooth(self,window):
	544	if self.verbose:
	545	print "!! WARNING !! Performing a smoothing with a window size",window
	546	print "!! WARNING !! To come back to unsmoothed file, use .get() again"
	547	for i in range(self.nfin):
	548	for j in range(self.nvin):
	549	for t in range(self.nplott):
	550	for z in range(self.nplotz):
	551	for y in range(self.nploty):
	552	for x in range(self.nplotx):
	553	obj = self.request[i][j][t][z][y][x]
	554	if obj.field.ndim == 1:
	555	print "!! ERROR !! 1D smoothing not supported yet because reduces array sizes."
	556	exit()
	557	# TBD TBD TBD
	558	#obj.field = ppcompute.smooth1d(obj.field,window=window)
	559	elif obj.field.ndim == 2:
	560	obj.field = ppcompute.smooth2d(obj.field,window=window)
	561
	562	##############################################################################################
	563	# defineplot method
	564	# --> defineplot first defines arrays of plot objects and set each of them
	565	# ... simple looping except cases where goal is not main (e.g. contour or vector)
	566	# --> principle: each onerequest() object gives birth to a subplot
	567	# --> defineplot vs. makeplot: defining plot and actually plotting it are clearly separated
	568	# --> THE KEY OUPUT OF defineplot IS AN ARRAY self.p OF PLOT OBJECTS
	569	# optional arguments
	570	# --> extraplot: to indicate a number of plots to be added afterwards (use self.plotin)
	571	# --> loadfile: to use self.p from a previously saved file
	572	##############################################################################################
	573	def defineplot(self,extraplot=0,loadfile=None):
	574	# -----------------------------------------------------
	575	# LOAD MODE: load a self.p object. count plots from it.
	576	# -----------------------------------------------------
	577	if loadfile is not None:
	578	try: filehandler = open(loadfile, 'r') ; self.p = pickle.load(filehandler)
	579	except IOError: print "!! ERROR !! Cannot find object file to load." ; exit()
	580	self.status = "definedplot" ; self.plotin = None
	581	self.nplot = len(self.p) ; self.howmanyplots = self.nplot
	582	return
	583	# -----------------------------------------------------
	584	# REGULAR MODE
	585	# -----------------------------------------------------
	586	self.printstatus()
	587	# check if things were done OK before
	588	if self.status in ["init","defined"]:
	589	print "!! ERROR !! Please use .retrieve() to get fields for plots with your pp object." ; exit()
	590	# check self.plotin (an existing fig on which to add plots afterwards)
	591	if self.plotin.__class__.__name__ == "pp":
	592	if self.plotin.fig is None:
	593	self.plotin = None # this is an additional security in case
	594	# a pp object is given without figure opened yet.
	595	elif self.plotin is not None:
	596	print "!! ERROR !! plotin argument must be a pp object." ; exit()
	597	# initialize the array of subplot objects
	598	# either something new or attributes coming from plotin object
	599	if self.plotin is None: self.p = [ ]
	600	else: self.p = self.plotin.p
	601	# create an array of subplot objects
	602	# ... in theory the order of looping can be changed without any harm
	603	# ... the only important thing is to keep i,j,t,z,y,x resp. for file,var,t,z,y,x
	604	count = 0
	605	for i in range(self.nfin):
	606	if self.filegoal[i] == "main":
	607	for j in range(self.nvin):
	608	if self.vargoal[j] == "main":
	609	for t in range(self.nplott):
	610	for z in range(self.nplotz):
	611	for y in range(self.nploty):
	612	for x in range(self.nplotx):
	613	# look at dimension and append the right plot object
	614	obj = self.request[i][j][t][z][y][x]
	615	dp = obj.dimplot
	616	if dp == 1 or self.forcedimplot == 1: plobj = ppplot.plot1d()
	617	elif dp == 2 or self.forcedimplot == 2: plobj = ppplot.plot2d()
	618	elif dp == 0: print "**** OK. VALUES VALUES VALUES",obj.field
	619	else: print "!! ERROR !! 3D or 4D plots not supported" ; exit()
	620	# load abscissa and ordinate in obj
	621	obj.definecoord()
	622	# we start to define things here before appending
	623	# (convenient: could be overridden by the user before makeplot)
	624	# ... the if loop is necessary so that we can loop above on the dp=0 case
	625	if dp in [1,2]:
	626	# and define what to do in plobj
	627	plobj.invert = obj.invert_axes
	628	plobj.swap = obj.swap_axes
	629	# axis labels
	630	plobj.xlabel = obj.absclab ; plobj.ylabel = obj.ordilab
	631	# superpose or not (this is mostly for saving purpose)
	632	plobj.superpose = self.superpose
	633	# get title, colormaps, labels, etc.. from var
	634	plobj.var = obj.var
	635	plobj.define_from_var()
	636	# generic 1D/2D: load field and coord in plot object
	637	plobj.field = obj.field # field to be plotted
	638	plobj.absc = obj.absc # abscissa (or longitude)
	639	plobj.ordi = obj.ordi # ordinate (or latitude)
	640	# -- useless in 1D but not used anyway
	641	if dp == 2:
	642	# specific 2d plot stuff
	643	# -- light grey background for missing values
	644	if type(plobj.field).__name__ in 'MaskedArray': plobj.axisbg = '0.75'
	645	# -- define if it is a map or a plot
	646	plobj.mapmode = ( obj.method_x+obj.method_y == "freefree" \
	647	and "grid points" not in obj.name_x \
	648	and not self.noproj )
	649	# finally append
	650	self.p.append(plobj)
	651	count = count + 1
	652	# self.nplot is number of plot to be defined in this call to defineplot()
	653	# (because of self.plotin this might less than length of self.p)
	654	self.nplot = count
	655	# --- superimposed contours and vectors ---
	656	# we have to start another loop because we need forward information
	657	# TBD: there is probably a more flexible way to do that
	658	count = 0
	659	for i in range(self.nfin):
	660	for j in range(self.nvin):
	661	for t in range(self.nplott):
	662	for z in range(self.nplotz):
	663	for y in range(self.nploty):
	664	for x in range(self.nplotx):
	665	goal = self.vargoal[j] + self.filegoal[i]
	666	obj = self.request[i][j][t][z][y][x]
	667	if "mainmain" in goal and obj.dimplot == 2:
	668	# the plot object we consider in the loop
	669	pl = self.p[count]
	670	# -- see if there is a contour requested...
	671	# (we use try because we might be at the end of the list)
	672	found = 0
	673	try: condvar = self.vargoal[j+1]
	674	except: condvar = "itisok"
	675	try: condfile = self.filegoal[i+1]
	676	except: condfile = "itisok"
	677	# ... get contour
	678	##########################################
	679	# NB: contour is expected to be right after main otherwise it is not displayed
	680	##########################################
	681	if condvar == "contour":
	682	plobj.addcontour = self.request[i][j+1][t][z][y][x].field ; found += 1
	683	if condfile == "contour":
	684	plobj.addcontour = self.request[i+1][j][t][z][y][x].field ; found += 1
	685	# see if there is a vector requested...
	686	# (we use try because we might be at the end of the list)
	687	try: condvar = self.vargoal[j+found+1]+self.vargoal[j+found+2]
	688	except: condvar = "itisok"
	689	try: condfile = self.filegoal[i+found+1]+self.filegoal[i+found+2]
	690	except: condfile = "itisok"
	691	# ... get vector and go directly to the next iteration
	692	# (in some cases we would do this twice but this is cheap)
	693	if "vector" in condvar:
	694	plobj.addvecx = self.request[i][j+found+1][t][z][y][x].field
	695	plobj.addvecy = self.request[i][j+found+2][t][z][y][x].field
	696	if "vector" in condfile:
	697	plobj.addvecx = self.request[i+found+1][j][t][z][y][x].field
	698	plobj.addvecy = self.request[i+found+2][j][t][z][y][x].field
	699	count = count + 1
	700	# COUNT PLOTS. if 0 just exit.
	701	# self.howmanyplots is self.nplot + possible extraplots
	702	self.howmanyplots = self.nplot + extraplot
	703	if self.howmanyplots > 0:
	704	if self.verbose: print "**** OK. expect %i plots" % (self.howmanyplots)
	705	else:
	706	exit() # because this means that we only had 0D values !
	707	# final status
	708	self.status = "definedplot"
	709
	710	##############################################################################################
	711	# makeplot method
	712	# --> after defineplot and before makeplot, user-defined plot settings can be easily given
	713	# simply by modifying the attributes of each elements of self.p
	714	# --> to change only one plot setting, no need to call defineplot again before makeplot
	715	# --> in the end, the array self.p of plot objects is saved for easy and convenient replotting
	716	# --> see practical examples in the folder 'examples'
	717	##############################################################################################
	718	def makeplot(self):
	719	self.printstatus()
	720	# a few initial operations
	721	# ------------------------
	722	if "definedplot" not in self.status:
	723	print "!! ERROR !! Please use .defineplot() before .makeplot() can be used with your pp object." ; exit()
	724	# open a figure and define subplots
	725	# ---------------------------------
	726	if self.plotin is None:
	727	# start from scratch
	728	self.fig = mpl.figure(figsize=(16,8))
	729	self.subv,self.subh = ppplot.definesubplot(self.howmanyplots,self.fig)
	730	self.n = 0
[917]	731	## adapted space for labels etc
	732	## ... except for ortho because there is no label anyway
	733	customplot = self.p[0].field.ndim == 2 \
	734	and self.p[0].mapmode == True \
	735	and self.p[0].proj not in ["ortho"]
	736	if customplot:
	737	margin = 0.07
	738	self.fig.subplots_adjust(left=margin,right=1-margin,bottom=margin,top=1-margin)
[910]	739	else:
	740	# start from an existing figure.
	741	# extraplot must have been set in the call to the previous figure.
	742	self.fig = self.plotin.fig
	743	self.subv,self.subh = self.plotin.subv,self.plotin.subh
	744	self.n = self.plotin.n
	745	self.howmanyplots = self.plotin.howmanyplots
	746	# LOOP on all subplots
	747	# NB: cannot use 'for pl in self.p' if self.plotin not None
	748	# --------------------
	749	for count in range(self.nplot):
	750	# the plot object we consider in the loop
	751	pl = self.p[self.n]
	752	# before making the plot, create a subplot. the first one is numbered 1 not 0.
	753	# ... if pl.superpose, we use only one and only figure
	754	# ... (and we have to be careful with not doing things several times)
	755	if pl.superpose:
	756	if self.n == 0:
	757	self.fig.add_subplot(1,1,1,axisbg=pl.axisbg) # define one subplot (still needed for user-defined font sizes)
[914]	758	sav = pl.xlabel,pl.ylabel,pl.xcoeff,pl.ycoeff,pl.title,pl.swaplab # save titles and labels
[910]	759	else:
	760	pl.invert = False ; pl.lstyle = None # don't invert again axis
[914]	761	# set saved titles and labels
	762	if self.plotin is None:
	763	pl.xlabel,pl.ylabel,pl.xcoeff,pl.ycoeff,pl.title,pl.swaplab = sav
	764	else:
	765	prev_plot = self.plotin.p[self.n-1]
	766	pl.xlabel = prev_plot.xlabel
	767	pl.ylabel = prev_plot.ylabel
	768	pl.xcoeff = prev_plot.xcoeff
	769	pl.ycoeff = prev_plot.ycoeff
	770	pl.title = prev_plot.title
	771	pl.swaplab = prev_plot.swaplab
[910]	772	else:
	773	self.fig.add_subplot(self.subv,self.subh,self.n+1,axisbg=pl.axisbg)
	774	if self.verbose: print "**** Done subplot %i / %i " %( self.n+1,self.howmanyplots )
	775	# finally make the plot
	776	pl.make()
	777	self.n = self.n+1
	778	# once completed show the plot (cannot show intermediate plotin)
[917]	779	# ... added a fix (customplot=True) for the label problem in basemap
[910]	780	print "**** Done step: makeplot"
[917]	781	if (self.n == self.howmanyplots):
	782	ppplot.save(mode=self.out,filename=self.filename,folder=self.folder,custom=customplot)
[920]	783	mpl.close()
[910]	784	# SAVE A PICKLE FILE WITH THE self.p ARRAY OF OBJECTS
	785	if self.verbose: print "**** Saving session in "+self.filename + ".ppobj"
	786	savfile = self.folder + "/" + self.filename + ".ppobj"
[920]	787	try:
	788	filehandler = open(savfile, 'w')
	789	pickle.dump(self.p, filehandler)
	790	except IOError:
	791	print "!! WARNING !! Saved object file not written. Probably do not have permission to write here."
[910]	792
	793	###########################################################
	794	# plot: a shortcut method for the defineplot + plot chain #
	795	###########################################################
	796	def plot(self,extraplot=0):
	797	self.defineplot(extraplot=extraplot)
	798	self.makeplot()
	799
	800	#######################################################
	801	# getplot: a shortcut method for the get + plot chain #
	802	#######################################################
	803	def getplot(self,extraplot=0):
	804	self.get()
	805	self.plot(extraplot=extraplot)
	806
	807	###################################################################
	808	# getdefineplot: a shortcut method for the get + defineplot chain #
	809	###################################################################
	810	def getdefineplot(self,extraplot=0):
	811	self.get()
	812	self.defineplot(extraplot=extraplot)
	813
	814	##############################################################
	815	# f: operation on two pp objects being on status 'definedplot'
	816	# this allows for one field being function of another one
	817	# e.g. u.f(v) means u will be displayed as a function of v
	818	# ... no need to do defineplot after u.f(v), makeplot directly
	819	##############################################################
	820	def f(self,other):
	821	# preamble: for this operation to work, defineplot() must have been done
	822	if self.status != "definedplot":
	823	if self.verbose: print "!! WARNING !! performing defineplot on operand"
	824	self.defineplot()
	825	if other.status != "definedplot":
	826	if self.verbose: print "!! WARNING !! performing defineplot on operand"
	827	other.defineplot()
	828	# check total number of plots
	829	if self.howmanyplots != other.howmanyplots:
	830	print "!! ERROR !! The two operands do not have the same number of subplots."
	831	exit()
	832	# and now operation.
	833	count = 0
	834	while count < self.howmanyplots:
	835	sobj = self.p[count] ; oobj = other.p[count]
	836	if sobj.field.ndim !=1 or oobj.field.ndim !=1:
	837	if self.verbose: print "!! WARNING !! Flattening arrays because more than one-dimensional."
	838	sobj.field = np.ravel(sobj.field)
	839	oobj.field = np.ravel(oobj.field)
	840	sobj.absc = oobj.field
	841	sobj.xlabel = oobj.ylabel
	842	if sobj.absc.size > sobj.field.size:
	843	if self.verbose:
	844	print "!! WARNING !! Trying to define y=f(x) with x and y not at the same size.",sobj.absc.size,sobj.field.size
	845	print "!! WARNING !! Modifying x to fit y size but please check."
	846	sobj.absc = sobj.absc[0:sobj.field.size]
	847	count = count + 1
	848	return self
	849
	850	###########################################################
	851	# copyopt: get options from e.g. a parser
	852	# ... allow for simple scripting and user-defined settings
	853	# ... must be called between defineplot and makeplot
	854	# REQUIRED: attributes of opt must be the same as in the pp object
	855	###########################################################
	856	def getopt(self,opt):
	857	# -- if only one, or less than the number of plots --> we take the first one
	858	# -- if as many as number of plots --> OK, each plot has its own setting
	859	# (except a few cases such as trans)
	860	for iii in range(self.howmanyplots):
	861	###
	862	try: self.p[iii].trans = opt.trans
	863	except: pass
	864	###
	865	try: self.p[iii].div = opt.div
	866	except: pass
	867	###
	868	try: self.p[iii].logy = opt.logy
	869	except: pass
	870	###
	871	try: self.p[iii].colorb = opt.colorb[iii]
	872	except:
	873	try: self.p[iii].colorb = opt.colorb[0]
	874	except: pass
	875	###
	876	try: self.p[iii].title = opt.title[iii]
	877	except:
	878	try: self.p[iii].title = opt.title[0]
	879	except: pass
	880	###
	881	try: self.p[iii].xlabel = opt.xlabel[iii]
	882	except:
	883	try: self.p[iii].xlabel = opt.xlabel[0]
	884	except: pass
	885	###
	886	try: self.p[iii].ylabel = opt.ylabel[iii]
	887	except:
	888	try: self.p[iii].ylabel = opt.ylabel[0]
	889	except: pass
	890	###
	891	try: self.p[iii].lstyle = opt.lstyle[iii]
	892	except:
	893	try: self.p[iii].lstyle = opt.lstyle[0]
	894	except: pass
	895	###
[920]	896	try: self.p[iii].color = opt.color[iii]
	897	except:
	898	try: self.p[iii].color = opt.color[0]
	899	except: pass
	900	###
	901	try: self.p[iii].marker = opt.marker[iii]
	902	except:
	903	try: self.p[iii].marker = opt.marker[0]
	904	except: pass
	905	###
[910]	906	try: self.p[iii].proj = opt.proj[iii]
	907	except:
	908	try: self.p[iii].proj = opt.proj[0]
	909	except: pass
	910	###
	911	try: self.p[iii].back = opt.back[iii]
	912	except:
	913	try: self.p[iii].back = opt.back[0]
	914	except: pass
	915	###
	916	try: self.p[iii].area = opt.area[iii]
	917	except:
	918	try: self.p[iii].area = opt.area[0]
	919	except: pass
	920	###
	921	try: self.p[iii].blon = opt.blon[iii]
	922	except:
	923	try: self.p[iii].blon = opt.blon[0]
	924	except: pass
	925	###
	926	try: self.p[iii].blat = opt.blat[iii]
	927	except:
	928	try: self.p[iii].blat = opt.blat[0]
	929	except: pass
	930	###
	931	try: self.p[iii].vmin = opt.vmin[iii]
	932	except:
	933	try: self.p[iii].vmin = opt.vmin[0]
	934	except: pass
	935	###
	936	try: self.p[iii].vmax = opt.vmax[iii]
	937	except:
	938	try: self.p[iii].vmax = opt.vmax[0]
	939	except: pass
	940
	941	##########################################################
	942	### THE ONEREQUEST SUBOBJECT TO PP (ON WHICH IT LOOPS) ###
	943	##########################################################
	944	class onerequest():
	945
	946	# default settings. mostly initialized to diagnose problem, except dimplot, nplot, verbose, swap_axes, invert_axes
	947	# -------------------------------
	948	def __init__(self):
	949	self.file = '!! file: I am not set, damned !!'
	950	self.f = None
	951	self.dim = None
	952	self.var = '!! var: I am not set, damned !!'
	953	self.index_x = [] ; self.index_y = [] ; self.index_z = [] ; self.index_t = []
	954	self.index_x2d = [] ; self.index_y2d = []
	955	self.method_x = '!! method_x: I am not set, damned !!'
	956	self.method_y = '!! method_y: I am not set, damned !!'
	957	self.method_z = '!! method_z: I am not set, damned !!'
	958	self.method_t = '!! method_t: I am not set, damned !!'
	959	self.field = None
	960	self.name_x = None ; self.name_y = None ; self.name_z = None ; self.name_t = None
	961	self.dim_x = None ; self.dim_y = None ; self.dim_z = None ; self.dim_t = None
	962	self.field_x = None ; self.field_y = None ; self.field_z = None ; self.field_t = None
	963	self.dimplot = 0
	964	self.nplot = 1
	965	self.absc = None ; self.ordi = None ; self.absclab = None ; self.ordilab = None
	966	self.verbose = True
	967	self.swap_axes = False ; self.invert_axes = False
	968	self.compute = None
	969
	970	# open a file. for now it is netcdf. TBD for other formats.
	971	# check that self.var is inside.
	972	# -------------------------------
	973	def openfile(self):
	974	if not os.path.exists(self.file): print '!! ERROR !! I could not find the following file: '+self.file ; exit()
	975	if not os.path.isfile(self.file): print '!! ERROR !! This does not appear to be a file: '+self.file ; exit()
	976	self.f = netCDF4.Dataset(self.file)
	977	if self.verbose: print "**** OK. Opened file "+self.file
	978	if self.var not in self.f.variables.keys():
	979	print '!! ERROR !! File '+self.file+' does not contain variable: '+self.var
	980	print '..... try instead with ',self.f.variables.keys() ; exit()
	981
	982	# get x,y,z,t dimensions from NETCDF file
	983	# TBD: user could request for a specific altitude dimension
	984	# TBD: staggered variables could request specific dimensions
	985	# -------------------------------
	986	def getdim(self):
	987	# GET SIZES OF EACH DIMENSION
	988	if self.verbose: print "**** OK. Found variable "+self.var
	989	shape = self.f.variables[self.var].shape
	990	self.dim = len(shape)
	991	if self.dim == 1:
	992	if self.verbose: print "**** OK. 1D field. I assume this varies with time."
	993	self.dim_x = 1 ; self.dim_y = 1 ; self.dim_z = 1 ; self.dim_t = shape[0]
	994	elif self.dim == 2:
	995	if self.verbose: print "**** OK. 2D field. I assume this is not-time-varying lat-lon map."
	996	self.dim_x = shape[1] ; self.dim_y = shape[0] ; self.dim_z = 1 ; self.dim_t = 1
	997	elif self.dim == 3:
	998	if self.verbose: print "**** OK. 3D field. I assume this is time-varying lat-lon map."
	999	self.dim_x = shape[2] ; self.dim_y = shape[1] ; self.dim_z = 1 ; self.dim_t = shape[0]
	1000	elif self.dim == 4:
	1001	if self.verbose: print "**** OK. 4D field."
	1002	self.dim_x = shape[3] ; self.dim_y = shape[2] ; self.dim_z = shape[1] ; self.dim_t = shape[0]
	1003	# LONGITUDE. Try preset fields. If not present set grid points axis.
	1004	self.name_x = "nothing"
	1005	for c in glob_listx:
	1006	if c in self.f.variables.keys():
	1007	self.name_x = c
	1008	if self.name_x == "nothing":
	1009	self.field_x = np.array(range(self.dim_x))
	1010	self.name_x = "x grid points"
	1011	else:
	1012	self.field_x = self.f.variables[self.name_x]
	1013	# LATITUDE. Try preset fields. If not present set grid points axis.
	1014	self.name_y = "nothing"
	1015	for c in glob_listy:
	1016	if c in self.f.variables.keys():
	1017	self.name_y = c
	1018	if self.name_y == "nothing":
	1019	self.field_y = np.array(range(self.dim_y))
	1020	self.name_y = "y grid points"
	1021	else:
	1022	self.field_y = self.f.variables[self.name_y]
	1023	# ensure that lon and lat are 2D fields
	1024	# 1. simple 1D case (not time-varying)
	1025	if len(self.field_x.shape)*len(self.field_y.shape) == 1:
	1026	if self.verbose: print "**** OK. recasting lon and lat as 2D fields."
	1027	[self.field_x,self.field_y] = np.meshgrid(self.field_x,self.field_y)
	1028	# 2. complex 3D case (time-varying, actually just copied over time axis)
	1029	elif len(self.field_x.shape)*len(self.field_y.shape) == 9:
	1030	if self.verbose: print "**** OK. reducing lon and lat as 2D fields. get rid of time."
	1031	self.field_x = self.field_x[0,:,:]
	1032	self.field_y = self.field_y[0,:,:]
	1033	# if xy axis are apparently undefined, set 2D grid points axis.
	1034	if "grid points" not in self.name_x:
	1035	if self.field_x.all() == self.field_x[0,0]:
	1036	print "!! WARNING !! xy axis look undefined. creating a non-dummy ones."
	1037	self.field_x = np.array(range(self.dim_x)) ; self.name_x = "x grid points"
	1038	self.field_y = np.array(range(self.dim_y)) ; self.name_y = "y grid points"
	1039	[self.field_x,self.field_y] = np.meshgrid(self.field_x,self.field_y)
	1040	if self.dim_x > 1:
	1041	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())
	1042	if self.dim_y > 1:
	1043	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())
	1044	# ALTITUDE. Try preset fields. If not present set grid points axis.
	1045	# WARNING: how do we do if several are available?
	1046	self.name_z = "nothing"
	1047	for c in glob_listz:
	1048	if c in self.f.variables.keys():
	1049	self.name_z = c
	1050	if self.name_z == "nothing":
	1051	self.field_z = np.array(range(self.dim_z))
	1052	self.name_z = "z grid points"
	1053	else:
	1054	self.field_z = self.f.variables[self.name_z][:] # specify dimension
	1055	# TBD: have to check that this is not a 3D field
	1056	if self.dim_z > 1:
	1057	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())
	1058	# TIME. Try preset fields.
	1059	self.name_t = "nothing"
	1060	for c in glob_listt:
	1061	if c in self.f.dimensions.keys():
	1062	self.name_t = c
	1063	try:
[911]	1064	# speed up: only get first value, last one.
[910]	1065	dafirst = self.f.variables[self.name_t][0]
	1066	dalast = self.f.variables[self.name_t][self.dim_t-1]
[911]	1067	self.field_t = np.linspace(dafirst,dalast,num=self.dim_t)
[910]	1068	if dafirst == dalast: self.field_t = np.array([dafirst])
	1069	except:
	1070	# ... or if a problem encountered, define a simple time axis
	1071	self.field_t = np.array(range(self.dim_t))
	1072	self.name_t = "t grid points"
	1073	if self.dim_t > 1:
	1074	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())
	1075
	1076	# get list of index to be retrieved for time axis
	1077	### TBD: il faudrait ne prendre que les indices qui correspondent a l interieur d un plot (dans all)
	1078	# -------------------------------
	1079	def getindextime(self,dalist,ind,stride):
	1080	if self.method_t == "free":
	1081	self.index_t = np.arange(0,self.dim_t,stride)
	1082	if self.dim_t > 1:
	1083	self.dimplot = self.dimplot + 1
	1084	if self.verbose: print "**** OK. t values. all."
	1085	else:
	1086	self.method_t = "fixed"
	1087	if self.verbose: print "**** OK. no t dimension."
	1088	elif self.method_t == "comp":
	1089	start = np.argmin( np.abs( self.field_t - dalist[ind][0] ) )
	1090	stop = np.argmin( np.abs( self.field_t - dalist[ind][1] ) )
	1091	self.index_t = np.arange(start,stop,stride)
	1092	if self.verbose: print "**** OK. t values. comp over interval ",self.field_t[start],self.field_t[stop]," nvalues=",self.index_t.size
	1093	elif self.method_t == "fixed":
	1094	self.index_t.append( np.argmin( np.abs( self.field_t - dalist[ind][0] ) ))
	1095	if self.verbose: print "**** OK. t values",self.field_t[self.index_t]
	1096	else:
	1097	print "!! ERROR !! method "+self.method_t+" not supported"
	1098	self.index_t = np.array(self.index_t)
	1099
	1100	# get list of index to be retrieved for vertical axis
	1101	### TBD: il faudrait ne prendre que les indices qui correspondent a l interieur d un plot (dans all)
	1102	# -------------------------------
	1103	def getindexvert(self,dalist,ind,stride):
	1104	if self.method_z == "free":
	1105	self.index_z = np.arange(0,self.dim_z,stride)
	1106	if self.dim_z > 1:
	1107	self.dimplot = self.dimplot + 1
	1108	if self.verbose: print "**** OK. z values. all."
	1109	else:
	1110	self.method_z = "fixed"
	1111	if self.verbose: print "**** OK. no z dimension."
	1112	elif self.method_z == "comp":
	1113	start = np.argmin( np.abs( self.field_z - dalist[ind][0] ) )
	1114	stop = np.argmin( np.abs( self.field_z - dalist[ind][1] ) )
	1115	self.index_z = np.arange(start,stop,stride)
	1116	if self.verbose: print "**** OK. z values. comp over interval",self.field_z[start],self.field_z[stop]," nvalues=",self.index_z.size
	1117	elif self.method_z == "fixed":
	1118	self.index_z.append( np.argmin( np.abs( self.field_z - dalist[ind][0] ) ))
	1119	if self.verbose: print "**** OK. z values",self.field_z[self.index_z]
	1120	else:
	1121	if self.verbose: print "!! ERROR !! method "+self.method_z+" not supported"
	1122	self.index_z = np.array(self.index_z)
	1123
	1124	# get list of index to be retrieved for horizontal grid
	1125	# --> index_x and index_y are slices to be retrieved from NETCDF files
	1126	# --> index_x2d and index_y2d are the actual (x,y) coordinates corresponding to each relevant point
	1127	# [this is slightly more complicated because 2D arrays for lat-lon projection possibly irregular]
	1128	# NB: to append index we use lists (the most convenient) then we convert into a numpy.array
	1129	### TBD: il faudrait ne prendre que les indices qui correspondent a l interieur d un plot (dans all)
	1130	# -------------------------------
	1131	def getindexhori(self,dalistx,dalisty,indx,indy,stridex,stridey):
	1132	## get what is the method over x and y axis
	1133	test = self.method_x+self.method_y
	1134	## CASE 0, EASY CASES:
	1135	## - LAT IS FREE (we do here what must be done whatever LON case is)
	1136	## - LON IS FREE (we do here what must be done whatever LAT case is)
	1137	## - LAT IS COMP AND LON IS FREE
	1138	## - LON IS COMP AND LAT IS FREE
	1139	if self.method_x == "free" or test in ["compfree","compcomp"]:
	1140	self.index_x = range(0,self.dim_x,stridex)
	1141	if self.dim_x > 1:
	1142	if self.method_x == "free": self.dimplot = self.dimplot + 1
	1143	if self.verbose: print "**** OK. x values. all."
	1144	else:
	1145	self.method_x = "fixed"
	1146	if self.verbose: print "**** OK. no x dimension."
[916]	1147	if self.method_y == "free" or test in ["freecomp","compcomp"]:
[910]	1148	self.index_y = range(0,self.dim_y,stridey)
	1149	if self.dim_y > 1:
	1150	if self.method_y == "free": self.dimplot = self.dimplot + 1
	1151	if self.verbose: print "**** OK. y values. all."
	1152	else:
	1153	self.method_y = "fixed"
	1154	if self.verbose: print "**** OK. no y dimension."
	1155	## CASE 0 above, this is just for continuity.
	1156	if self.method_x in ["free","comp"] and self.method_y in ["free","comp"]:
	1157	self.index_x2d = self.index_x
	1158	self.index_y2d = self.index_y
	1159	## AND NOW THE LITTLE BIT MORE COMPLICATED CASES
	1160	## CASE 1 LAT AND LON ARE FIXED
	1161	elif test == "fixedfixed":
	1162	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 )
	1163	#TBD: pb with staggered coord
	1164	if idx not in self.index_x: self.index_x.append(idx)
	1165	if idy not in self.index_y: self.index_y.append(idy)
	1166	self.index_x2d.append(idx)
	1167	self.index_y2d.append(idy)
	1168	## CASE 2 LON IS FIXED BUT NOT LAT
	1169	elif test in ["fixedfree","fixedcomp"]:
	1170	# find where are requested x values for each y on the free dimension
	1171	# NB: this does not work for non-bijective cases e.g. polar stereographic
	1172	for iy in range(self.dim_y):
	1173	idx = np.argmin( np.abs( self.field_x[iy,:] - dalistx[indx][0] ) )
	1174	# if comp is requested we select only indexes which yield values between requested min and max
	1175	storeval = (self.method_y == "comp") and (self.field_y[iy,idx] > dalisty[indy][0]) and (self.field_y[iy,idx] < dalisty[indy][1])
	1176	storeval = storeval or (self.method_y == "free")
	1177	if storeval:
	1178	if idx not in self.index_x: self.index_x.append(idx)
	1179	if iy not in self.index_y and self.method_y == "comp": self.index_y.append(iy)
	1180	if idx not in self.index_x2d or iy not in self.index_y2d:
	1181	self.index_x2d.append(idx)
	1182	self.index_y2d.append(iy)
	1183	## CASE 3 LAT IS FIXED BUT NOT LON
	1184	elif test in ["freefixed","compfixed"]:
	1185	# find where are requested y values for each x on the free dimension
	1186	# NB: this does not work for non-bijective cases e.g. polar stereographic
	1187	for ix in range(self.dim_x):
	1188	idy = np.argmin( np.abs( self.field_y[:,ix] - dalisty[indy][0] ) )
	1189	# if comp is requested we select only indexes which yield values between requested min and max
	1190	storeval = (self.method_x == "comp") and (self.field_x[idy,ix] > dalistx[indx][0]) and (self.field_x[idy,ix] < dalistx[indx][1])
	1191	storeval = storeval or (self.method_x == "free")
	1192	if storeval:
	1193	if idy not in self.index_y: self.index_y.append(idy)
	1194	if ix not in self.index_x and self.method_x == "comp": self.index_x.append(ix)
	1195	if ix not in self.index_x2d or idy not in self.index_y2d:
	1196	self.index_x2d.append(ix)
	1197	self.index_y2d.append(idy)
	1198	## check index tab
	1199	if len(self.index_x) == 0 or len(self.index_y) == 0:
	1200	print "!! ERROR !! no indices found. check prescribed latitudes or longitudes" ; exit()
	1201	## ensure the array is a numpy array for getfield to work
	1202	self.index_x = np.array(self.index_x)
	1203	self.index_y = np.array(self.index_y)
	1204	self.index_x2d = np.array(self.index_x2d)
	1205	self.index_y2d = np.array(self.index_y2d)
	1206	### print extrema
	1207	printx = self.field_x[np.ix_(self.index_y2d, self.index_x2d)]
	1208	printy = self.field_y[np.ix_(self.index_y2d, self.index_x2d)]
	1209	if self.verbose:
	1210	print "**** OK. x values (min,max).", printx.min(),printx.max()
	1211	print "**** OK. y values (min,max).", printy.min(),printy.max()
	1212
	1213	# get the field from the NETCDF file and perform averages
	1214	# -------------------------------
	1215	def getfield(self):
	1216	## first tell what is to be done
	1217	if self.dimplot > 2: print "**** !! ERROR !! "+str(self.dimplot)+"D plots not supported!" ; exit()
	1218	elif self.dimplot == 0 and self.verbose: print "**** OK. 0D value requested."
	1219	elif self.dimplot == 1 and self.verbose: print "**** OK. 1D plot requested."
	1220	elif self.verbose: print "**** OK. 2D section requested."
	1221	# well, now get field from netcdf file
	1222	# part below is necessary otherwise there is an index error below
	1223	if self.index_x.size == 1: self.index_x = self.index_x[0]
	1224	if self.index_y.size == 1: self.index_y = self.index_y[0]
	1225	if self.index_z.size == 1: self.index_z = self.index_z[0]
	1226	if self.index_t.size == 1: self.index_t = self.index_t[0]
	1227	# then retrieve what is requested by user
	1228	# each self.dim case corresponds to tests in the beginning of getdim.
	1229	time0 = timelib.time()
	1230	if self.verbose: print "**** OK. I am getting values from files. Please wait."
	1231	if self.dim == 1:
	1232	nt = self.index_t.size ; nz = 1 ; ny = 1 ; nx = 1
	1233	self.field = self.f.variables[self.var][self.index_t]
	1234	elif self.dim == 2:
	1235	nt = 1 ; nz = 1 ; ny = self.index_y.size ; nx = self.index_x.size
	1236	self.field = self.f.variables[self.var][self.index_y,self.index_x]
	1237	elif self.dim == 3:
	1238	nt = self.index_t.size ; nz = 1 ; ny = self.index_y.size ; nx = self.index_x.size
	1239	self.field = self.f.variables[self.var][self.index_t,self.index_y,self.index_x]
	1240	# this is far faster than retrieving each term with a loop
	1241	elif self.dim == 4:
	1242	nt = self.index_t.size ; nz = self.index_z.size ; ny = self.index_y.size ; nx = self.index_x.size
	1243	self.field = self.f.variables[self.var][self.index_t,self.index_z,self.index_y,self.index_x]
	1244	else:
	1245	print "!! ERROR !! field would have more than four dimensions ?" ; exit()
	1246	# NB: ... always 4D array but possibly with "size 1" dimensions
	1247	# ... if one dimension is missing because 1D 2D or 3D requests, make it appear again
	1248	self.field = np.reshape(self.field,(nt,nz,ny,nx))
	1249	if self.verbose: print "**** OK. I got %7.1e values. This took me %6.4f seconds" % (nxnynz*nt,timelib.time() - time0)
	1250	if self.verbose: print "**** OK. I got var "+self.var+" with shape",self.field.shape
	1251	# reduce coordinates to useful points
	1252	# ... TBD: this should be ordered in the case of non-regular projections
	1253	if self.method_x in ["free","comp"] and self.method_y in ["free","comp"]:
	1254	# we need 2D coordinates (free) or we get broadcast problem (comp) so we use np.ix
	1255	self.field_x = self.field_x[np.ix_(self.index_y2d, self.index_x2d)]
	1256	self.field_y = self.field_y[np.ix_(self.index_y2d, self.index_x2d)]
	1257	else:
	1258	# we are OK with 1D coordinates
	1259	self.field_x = self.field_x[self.index_y2d, self.index_x2d]
	1260	self.field_y = self.field_y[self.index_y2d, self.index_x2d]
	1261	self.field_z = self.field_z[self.index_z]
	1262	self.field_t = self.field_t[self.index_t]
	1263	# now have to obtain the new indexes which correspond to the extracted self.field
	1264	# for it to work with unique index, ensure that any index_* is a numpy array
	1265	if not isinstance(self.index_x, np.ndarray): self.index_x = np.array([self.index_x])
	1266	if not isinstance(self.index_y, np.ndarray): self.index_y = np.array([self.index_y])
	1267	if not isinstance(self.index_z, np.ndarray): self.index_z = np.array([self.index_z])
	1268	if not isinstance(self.index_t, np.ndarray): self.index_t = np.array([self.index_t])
	1269	for val in self.index_x: self.index_x2d[np.where(self.index_x2d == val)] = np.where(self.index_x == val)[0]
	1270	for val in self.index_y: self.index_y2d[np.where(self.index_y2d == val)] = np.where(self.index_y == val)[0]
	1271	for val in self.index_z: self.index_z [np.where(self.index_z == val)] = np.where(self.index_z == val)[0]
	1272	for val in self.index_t: self.index_t [np.where(self.index_t == val)] = np.where(self.index_t == val)[0]
	1273	##### VERY EXPENSIVE
	1274	## recast self.field with 2D horizontal arrays because we might have extracted
	1275	## more than what is to be actually plot or computed, in particular for comps on 2D lat,lon coordinates
	1276	#self.field = self.field[np.ix_(self.index_t,self.index_z,self.index_y2d,self.index_x2d)]
	1277	#(nt,nz,ny,nx) = self.field.shape
	1278	# extract relevant horizontal points
	1279	# TBD: is compfree OK with computing on irregular grid?
	1280	test = self.method_x + self.method_y
	1281	if test in ["fixedfixed","freefree"]:
	1282	pass
	1283	elif test in ["fixedfree","fixedcomp"] or test in ["freefixed","compfixed"]:
	1284	time0 = timelib.time()
	1285	# prepare the loop on all relevant horizontal points
	1286	if self.method_x in ["comp","free"]:
	1287	nnn = self.index_x2d.shape[0] ; what_I_am_supposed_to_do = "keepx"
	1288	elif self.method_y in ["comp","free"]:
	1289	nnn = self.index_y2d.shape[0] ; what_I_am_supposed_to_do = "keepy"
	1290	# LOOP to extract only useful values over horizontal dimensions
	1291	# only take diagonal terms, do not loop on all self.index_x2d*self.index_y2d
	1292	# ... this method is fast enough, perhaps there is a faster way though
	1293	# ... (for sure the method with np.diag is much slower)
	1294	for iii in range(nnn):
	1295	ix = self.index_x2d[iii] ; iy = self.index_y2d[iii]
	1296	for iz in self.index_z:
	1297	for it in self.index_t:
	1298	if what_I_am_supposed_to_do == "keepx": self.field[it,iz,0,ix] = self.field[it,iz,iy,ix]
	1299	elif what_I_am_supposed_to_do == "keepy": self.field[it,iz,iy,0] = self.field[it,iz,iy,ix]
	1300	if self.verbose: print "**** OK. I got to pick the right values for your request. This took me %6.4f seconds" % (timelib.time() - time0)
	1301	# we only keep the one value that was modified on the dimension which is not free
	1302	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))
	1303	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))
	1304	# make a mask in case there are non-NaN missing values. (what about NaN missing values?)
	1305	# ... this is important for computations below (see ppcompute)
	1306	masked = np.ma.masked_where(np.abs(self.field) > 1e25,self.field)
	1307	if masked.mask.any() == True:
	1308	if self.verbose: print "!! WARNING !! Values over +-1e25 are considered missing values."
	1309	self.field = masked
	1310	self.field.set_fill_value([np.NaN])
	1311	# now ready to compute [TBD?? we would like to have e.g. mean over x,y and min over t]
	1312	if self.method_t == "comp":
	1313	if self.verbose: print "**** OK. Computing over t axis."
	1314	if self.compute == "mean": self.field = ppcompute.mean(self.field,axis=0)
	1315	elif self.compute == "min": self.field = ppcompute.min(self.field,axis=0)
	1316	elif self.compute == "max": self.field = ppcompute.max(self.field,axis=0)
	1317	else: print "!! ERROR !! operation not supported." ; exit()
	1318	nt = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
	1319	if self.method_z == "comp":
	1320	if self.verbose: print "**** OK. Computing over z axis."
	1321	if self.compute == "mean": self.field = ppcompute.mean(self.field,axis=1)
	1322	elif self.compute == "min": self.field = ppcompute.min(self.field,axis=1)
	1323	elif self.compute == "max": self.field = ppcompute.max(self.field,axis=1)
	1324	nz = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
	1325	if self.method_y == "comp":
	1326	if self.verbose: print "**** OK. Computing over y axis."
	1327	if self.compute == "mean": self.field = ppcompute.mean(self.field,axis=2)
	1328	elif self.compute == "min": self.field = ppcompute.min(self.field,axis=2)
	1329	elif self.compute == "max": self.field = ppcompute.max(self.field,axis=2)
	1330	ny = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
	1331	if self.field_x.ndim == 2: self.field_x = self.field_x[0,:] # TBD: this is OK for regular grid but not for irregular
	1332	if self.method_x == "comp":
	1333	if self.verbose: print "**** OK. Computing over x axis."
	1334	if self.compute == "mean": self.field = ppcompute.mean(self.field,axis=3)
	1335	elif self.compute == "min": self.field = ppcompute.min(self.field,axis=3)
	1336	elif self.compute == "max": self.field = ppcompute.max(self.field,axis=3)
	1337	nx = 1 ; self.field = np.reshape(self.field,(nt,nz,ny,nx))
	1338	if self.field_y.ndim == 2: self.field_y = self.field_y[:,0] # TBD: this is OK for regular grid but not for irregular
	1339	# remove all dimensions with size 1 to prepare plot (and check the resulting dimension with dimplot)
	1340	self.field = np.squeeze(self.field)
	1341	if self.field.ndim != self.dimplot:
	1342	print "!! ERROR !! Problem: self.field is different than plot dimensions", self.field.ndim, self.dimplot ; exit()
	1343	if self.verbose:
	1344	print "**** OK. Final shape for "+self.var+" after averaging and squeezing",self.field.shape
	1345
	1346	# define coordinates for plot
	1347	# -------------------------------
	1348	def definecoord(self):
	1349	I_got_abs = False ; I_got_ord = False
	1350	# here is the thing. time is usually taken as an abscissa so we start with time.
	1351	if self.method_t == "free":
	1352	self.absc = self.field_t ; self.absclab = self.name_t
	1353	I_got_abs = True
	1354	# then we usually have x as an abscissa.
	1355	if self.method_x == "free":
	1356	if I_got_abs:
	1357	self.ordi = self.field_x ; self.ordilab = self.name_x
	1358	I_got_ord = True
	1359	else:
	1360	self.absc = self.field_x ; self.absclab = self.name_x
	1361	I_got_abs = True
	1362	# ... or we have y
	1363	if self.method_y == "free":
	1364	if I_got_abs:
	1365	self.ordi = self.field_y ; self.ordilab = self.name_y
	1366	I_got_ord = True
	1367	else:
	1368	self.absc = self.field_y ; self.absclab = self.name_y
	1369	I_got_abs = True
	1370	# ... and we end with z because it is usually not an abscissa (profiles).
	1371	if self.method_z == "free":
	1372	if self.field_z[0] > self.field_z[1]:
	1373	self.invert_axes = True # the axis will be turned upside-down
	1374	if I_got_abs:
	1375	self.ordi = self.field_z ; self.ordilab = self.name_z
	1376	I_got_ord = True
	1377	else:
	1378	self.absc = self.field_z ; self.absclab = self.name_z
	1379	I_got_abs = True
	1380	self.swap_axes = True # says that altitude is not supposed to remain as an abscissa
	1381	if I_got_abs and self.verbose: print "**** OK. abscissa:",self.absclab, self.absc.shape
	1382	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: