diff --git a/__init__.py b/__init__.py index 327665d..ff5a373 100644 --- a/__init__.py +++ b/__init__.py @@ -16,6 +16,7 @@ from clExport import clExport from ihm import HoverInfo from display_lib import display1d +from display_lib import display2d from display_lib import display3d from display_lib import display_colormaps diff --git a/ceflib_mods.py b/ceflib_mods.py index be52264..65194d9 100644 --- a/ceflib_mods.py +++ b/ceflib_mods.py @@ -1,67 +1,149 @@ #! /bin/env python # coding: utf8 +import numpy as np +import datetime class cefwrite(object) : - def __init__(self, data = None, filename = 'test.cef'): + def __init__(self, data = None, filename = 'test.cef', mode = 'w', checkVariableOrder = False): """ - Data of the form data['field1'] = numpy array, and data['epoch'] = numpy time array + Write data into a cef file format + data must be of the form data['field1'] = numpy array, and data['epoch'] = numpy time array + provide filename (default is test.cef) + mode : 'w' start from fresh empty file + 'a' start from existing file and add only lines (variables must have same names) + default is 'w' """ - self.data = data + self.data = {key:np.copy(data[key]) for key in data} + self.fichier = filename + self.mode = mode + + self.variables = sorted(list(set(self.data.keys()) - set(['epoch']))) + + self.deleteAllNan() + self.changeNanToFillVal() + if checkVariableOrder : print self.variables + + if mode == 'w': + self.makeHeader() + self.writeData() + elif mode == 'a': + self.writeData() + + def changeNanToFillVal(self): + """ + For CL : change NAN to -10E31 + seems to work fine :) + """ + for key in self.variables : self.data[key][np.isnan(self.data[key])] = -10**31 + + def deleteAnyNan(self): + """ + Delete all lines containing at least one np.nan + """ + for var in self.variables : + inds = np.where(np.isnan(self.data[var])) + for key in self.data : + self.data[key] = np.delete(self.data[key], inds) + + def deleteAllNan(self): + """ + Delete all lines where all data exept time are np.nan + """ + inds0 = set(np.where(np.isnan(self.data[self.variables[0]]))[0]) + for var in self.variables[1:] : + inds1 = set(np.where(np.isnan(self.data[var]))[0]).intersection(inds0) + inds0 = inds1 + + for key in self.data : + self.data[key] = np.delete(self.data[key], list(inds1)) + + def makeHeader(self): self.initializeHeader() - variables = sorted(list(set(self.data.keys()) - set(['epoch']))) - for var in variables : self.addvariableHeader(var) + for var in self.variables : self.addvariableHeader(var) self.closeHeader() - + + def writeData(self): out = open(self.fichier, 'a') for i,t in enumerate(self.data['epoch']): - - out.write(t.isoformat()) - #print >> out, t.isoformat() - for var in variables : - out.write(' ') - #print >> out, str(self.data[key][i]) + out.write(t.isoformat()+'Z') + for var in self.variables : + out.write(', ') out.write('{:f}'.format(self.data[var][i])) - out.write('\n') - + out.write(';\n') out.close() def initializeHeader(self): - #set time variable as first variable out = open(self.fichier, 'w') out.write( + "!------------------------------------------------------------------------|\n" + "! |\n" + "! Generated by the mms_tools_lib, developped by Yoann Vernisse, IRAP |\n" + "! ASCII Format |\n" + "!------------------------------------------------------------------------|\n" + "\n" + 'FILE_FORMAT_VERSION = "CEF-2.0"\n' + "\n" + "\n" + "!------------------------------------------------------------------------|\n" + "! Global Metadata |\n" + "!------------------------------------------------------------------------|\n" + "\n" + "START_META = Generation_date\n" + " VALUE_TYPE = ISO_TIME\n" + " ENTRY = "+ datetime.datetime.now().isoformat() +"\n" + "END_META = Generation_date\n" + "\n" + "START_META = Generated_by\n" + ' ENTRY = "mms_tools_lib"\n' + "END_META = Generated_by\n" + "\n" + "START_META = Mission_Group\n" + ' ENTRY = "MMS"\n' + "END_META = Mission_Group\n" + "\n" + 'END_OF_RECORD_MARKER = ";"\n' + "\n" + "!------------------------------------------------------------------------|\n" + "! Variables |\n" + "!------------------------------------------------------------------------|\n" + "\n" 'START_VARIABLE = epoch \n' - ' PARAMETER_TYPE = "Support_Data" \n' - ' VALUE_TYPE = ISO_TIME \n' - ' UNITS = "s" \n' - ' DELTA_PLUS = 0 \n' - ' DELTA_MINUS = 0 \n' - 'END_VARIABLE = epoch \n' - ' \n' + ' PARAMETER_TYPE = "Support_Data"\n' + ' VALUE_TYPE = ISO_TIME\n' + ' UNITS = "s"\n' + ' DELTA_PLUS = 0\n' + ' DELTA_MINUS = 0\n' + 'END_VARIABLE = epoch\n' + '\n' ) def addvariableHeader(self, variable): out = open(self.fichier, 'a') out.write( - 'START_VARIABLE = '+variable +'\n' - ' PARAMETER_TYPE = "Data" \n' - ' VALUE_TYPE = FLOAT \n' - ' FILLVAL = -1.000E+31 \n' - ' DEPEND_0 = epoch \n' - ' UNITS = "" \n' - ' FIELDNAM = "" \n' + 'START_VARIABLE = '+variable +'\n' + ' PARAMETER_TYPE = "Data"\n' + ' VALUE_TYPE = FLOAT\n' + ' FILLVAL = -1.000E+31\n' + ' DEPEND_0 = epoch\n' + ' UNITS = ""\n' + ' FIELDNAM = ""\n' 'END_VARIABLE = '+variable+ '\n' - ' \n' + '\n' ) def closeHeader(self): out = open(self.fichier, 'a') out.write( - 'DATA_UNTIL = EOF \n' - ) + "!------------------------------------------------------------------------|\n" + "! Data |\n" + "!------------------------------------------------------------------------|\n" + '\n' + 'DATA_UNTIL = EOF\n' + ) \ No newline at end of file diff --git a/display_lib.py b/display_lib.py index ae1fd50..ae306d2 100644 --- a/display_lib.py +++ b/display_lib.py @@ -4,8 +4,41 @@ import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D +import os -class display1d(object): + +class display_general(object): + """ + Object with general functions for plotting objects + namely : save figure, show; others? + """ + + def __init__(self): + + self.setTitle('Figure_from_displaylib') + + def setTitle(self, title) : + self.title = title + + def getTitle(self) : + return self.title + + def savefig(self, f = None) : + if not f : + i = 0 + while i < 999 : + f = self.title + str(int(i)) + '.png' + if os.path.exists(os.path.join(os.curdir, f)) : i+=1 + else : break + print 'Saving figure in ', f + self.fig.savefig(f, bbox_inches='tight') + plt.close() + + def show(self): + plt.show() + + +class display1d(display_general): """ Functions to easily add options to 1d objects plots initialize as display1d(3,2) for 3 columns and 2 rows @@ -15,8 +48,9 @@ class display1d(object): add colorscale to a scatter plot NOTE : Not made for timeseries!! """ - def __init__(self, ncols = 1, nrows = 1) : + def __init__(self, ncols = 1, nrows = 1, title = 'display1d') : + self.setTitle(title) self.fig = plt.figure(figsize = (16, 10), dpi = 96) self.ncols = int(ncols) self.nrows = int(nrows) @@ -27,7 +61,7 @@ class display1d(object): for i in range(nrows*ncols-1) : self.plots.append({}) self.plots[-1]['label'] = {'data' : None, 'xaxis' : None, 'yaxis' : None, 'cbar' : None, 'title' : None} - + #self.set_colorscale() #self.set_errorbar(x = np.arange(len(self.y))/5.) #self.set_logscale() @@ -84,7 +118,9 @@ class display1d(object): def set_colorscale(self, col_arr = [0], minmax = 'auto', cmap = 'plasma') : i = self.get_actualPlotNumber() self.plots[i]['col_arr'] = np.array(col_arr) - if self.plots[i]['col_arr'].shape != self.plots[i]['data'].shape : self.plots[i]['col_arr'] = self.plots[i]['data'] + if len(self.plots[i]['col_arr']) != len(self.plots[i]['data']) : + print 'WRONG COLOR ARRAY SIZE, SETTING DATA IN COLOR INSTEAD' + self.plots[i]['col_arr'] = self.plots[i]['data'] #self.new_val = (self.col_arr - self.col_arr.min())/(self.col_arr.max() - self.col_arr.min()) #self.carray = plt.get_cmap(cmap)(self.new_val) @@ -126,7 +162,100 @@ class display1d(object): -class display3d(object) : +class display2d(display_general): + """ + Grid used to make 2D plots. + This create ONLY linear grid plot. For non linear 2D plot, another routine should be made. + input : 2D numpy array + NOTE : meshgrid can be tricky regarding x/y indexing : check option indexing = 'ij' or 'xy' + """ + def __init__(self, z = None, title = '', Zinit = True): + #np.ndarray.__init__(self) + + self.x = {} + self.y = {} + self.z = {} + self.minmax = {} + self.title = title + self.layers = {} + + if z is None : z = np.zeros([10,10]) + + xsize = float(z.shape[0]) + ysize = float(z.shape[1]) + + self.initialize(xsize, ysize) + self.meshinit() + self.zmesh, tmp = np.meshgrid(np.zeros(z.shape[0]), np.zeros(z.shape[1])) + self.zmesh = self.zmesh[:-1, :-1] + if Zinit : self.zmesh = z.T[:-1, :-1] + + def initialize(self, xsize, ysize) : + self.minmax['xmin'] = 0. + self.minmax['xmax'] = xsize + self.x['res'] = xsize + + self.minmax['ymin'] = 0. + self.minmax['ymax'] = ysize + self.y['res'] = ysize + + def set_layers(self, layers = ['Name']): + for layer in layers : self.layers[layer] = np.zeros(self.zmesh.shape) + + def meshinit(self) : + self.x['points'] = np.linspace(self.minmax['xmin'], self.minmax['xmax'], num = self.x['res'], dtype = 'double', endpoint = True) + self.y['points'] = np.linspace(self.minmax['ymin'], self.minmax['ymax'], num = self.y['res'], dtype = 'double', endpoint = True) + self.xmesh, self.ymesh = np.meshgrid(self.x['points'], self.y['points']) + + def set_minmax(self, cle, val) : #cle = xmax, xmin, ymax, ymin, zmin, zmax + self.minmax[cle] = float(val) + self.meshinit() + + def set_layers_val(self, x, y, layers_val = {'Name' : 0}): + xind = self.postoindicex(x) + yind = self.postoindicey(y) + if xind is not None and yind is not None : + for layer in self.layers : self.layers[layer][yind, xind] = layers_val[layer] + + def set_val(self, x, y, val): #x et y en position, PAS EN INDICE!! + xind = self.postoindicex(x) + yind = self.postoindicey(y) + if xind is not None and yind is not None : self.zmesh[yind, xind] = val + + def postoindicex(self, x): + if self.minmax['xmin'] < x < self.minmax['xmax'] : return (abs(self.x['points'] - x)).argmin() + #else : print 'X OUT OF BOUNDS', x, self.x['min'], self.x['max'] + def postoindicey(self, x): + if self.minmax['ymin'] < x < self.minmax['ymax'] : return (abs(self.y['points'] - x)).argmin() + #else : print 'Y OUT OF BOUNDS', x, self.y['min'], self.y['max'] + + def add_plot(self, x, y) : + plt.plot(x, y, 'H', color = 'white') + + def plot(self, layer = None, autominmax = True, cmap = 'seismic') : + if layer : z = self.layers[layer] + else : z = self.zmesh + + if autominmax or 'zmin' not in self.minmax or 'zmax' not in self.minmax: + vmin = z.min() + vmax = z.max() + else : + vmin = self.minmax['zmin'] + vmax = self.minmax['zmax'] + + self.fig = plt.figure() + plt.pcolormesh(self.xmesh, self.ymesh, z, vmin = vmin, vmax = vmax, cmap = cmap) + plt.title(self.title) + plt.colorbar() + + def xlineplot(self, x) : + simumods.lineplot(self.zmesh[self.postoindicex(x), :], x = self.xmesh[self.postoindicex(0), :]) + def ylineplot(self, x) : + simumods.lineplot(self.zmesh[:, self.postoindicey(x)], x = self.ymesh[:, self.postoindicey(0)]) + + + +class display3d(display_general) : """ 3D plotting object works fine! init 3D coordinates as numpy arrays (x,y,z) @@ -135,14 +264,16 @@ class display3d(object) : Thus this object impose an invisible cube at first use set_autoscatter for a default scatter plot of x,y,z (NOT generated by default, must be called!!) """ - def __init__(self, x, y, z) : + def __init__(self, x, y, z, title = 'display3d') : self.x = x self.y = y self.z = z self.planet = False + self.trisurf = False self._initialize() self.plots = {} + self.setTitle(title) def _initialize(self): @@ -152,10 +283,22 @@ class display3d(object) : ### Create a cube to force aspect ratio 1:1 self.ax.set_aspect('equal') - MAX = max(np.concatenate([abs(self.x),abs(self.y),abs(self.z)])) + self._max = max(np.concatenate([abs(self.x),abs(self.y),abs(self.z)])) for direction in (-1, 1): - for point in np.diag(direction * MAX * np.array([1,1,1])): + for point in np.diag(direction * self._max * np.array([1,1,1])): self.ax.plot([point[0]], [point[1]], [point[2]], 'w') + + self.set_xyzlabels() + + def set_xyzlabels(self, x = None, y = None, z = None): + + if x is None : x = 'X-axis' + if y is None : y = 'Y-axis' + if z is None : z = 'Z-axis' + + self.ax.set_xlabel(x) + self.ax.set_ylabel(y) + self.ax.set_zlabel(z) def forceEqualRatio(self): @@ -164,12 +307,11 @@ class display3d(object) : for key in self.plots : mlist.append(abs(np.concatenate(self.plots[key]['data'])).max()) mlist.append(abs(np.concatenate([self.x,self.y,self.z])).max()) #add any other data to the list - MAX = max(mlist) + self._max = max(mlist) for direction in (-1, 1): - for point in np.diag(direction * MAX * np.array([1,1,1])): + for point in np.diag(direction * self._max * np.array([1,1,1])): self.ax.plot([point[0]], [point[1]], [point[2]], 'w') - def redraw(self): self._initialize() if self.trisurf : self.set_trisurf(tripos = self.trisurf['tripos'], trival = self.trisurf['trival'], minmax = self.trisurf['minmax']) @@ -201,10 +343,16 @@ class display3d(object) : def set_scatter(self, scatter, c = None) : self.scatter = self.ax.scatter(scatter[:,0], scatter[:,1], scatter[:,2], c = c, cmap = 'seismic') - def set_autoscatter(self, c = None, vmin = None, vmax = None): - self.autoscatter = self.ax.scatter(self.x, self.y, self.z, vmin = vmin, vmax = vmax, c = c, cmap = 'seismic') - cbar = self.fig.colorbar(self.autoscatter) - + def set_autoscatter(self, c = None, vmin = None, vmax = None, **kw): + self.autoscatter = self.ax.scatter(self.x, self.y, self.z, vmin = vmin, vmax = vmax, c = c, cmap = 'seismic', **kw) + if c is not None : cbar = self.fig.colorbar(self.autoscatter) + + + def set_centralAxis(self, linewidth = .5): + self.ax.plot([-self._max, self._max], [0,0], [0,0], '-k', linewidth = linewidth) + self.ax.plot([0,0], [-self._max, self._max], [0,0], '-k', linewidth = linewidth) + self.ax.plot([0,0], [0,0], [-self._max, self._max], '-k', linewidth = linewidth) + def set_planet(self) : self.planet = True u, v = np.mgrid[0:2*np.pi:20j, 0:np.pi:10j] @@ -223,10 +371,6 @@ class display3d(object) : if save : self.plots[len(self.plots)] = {'data' : plot, 'color' : color} self.ax.plot(plot[:,0], plot[:,1], plot[:,2], linewidth = 3., color = color) - def show(self) : - plt.show() - - diff --git a/mva.py b/mva.py index cdf40bc..e196e19 100644 --- a/mva.py +++ b/mva.py @@ -23,7 +23,7 @@ class mva(object): self.results = {} self.calc(noprint = noprint) #self.noprint = noprint - if field2 != None : + if field2 is not None : self.field_cross = field2 # en entré : deux vecteurs de 3 composantes déjà moyenné self.calc_cross(noprint = noprint) diff --git a/read_data.py b/read_data.py index d257a14..a74c233 100644 --- a/read_data.py +++ b/read_data.py @@ -7,6 +7,7 @@ from datetime import timedelta import numpy as np import os import matplotlib.pyplot as plt +import pickle class donnee(object) : @@ -14,9 +15,15 @@ class donnee(object) : """ Donnee object : open data from cef file This object is useful to manage timeseries with implemented calculation on timedeltas - get metadata from ceffile and store it in self._metadata - to load data use self._get_data() (THIS MUST BE IMPROVED) + all data are stored in self.donnee + REMAINING TO DO : Harmonize all methods conventional name (use either _ of capital) + Rename self.donnee into self.data + Rename the class into a more appropriate international word + Add Nan & Fillvall treatment to all methods + Add a write to cef file method + Add a get_as_tensor method (similar to get_as_field) + Vectorize all methods """ ceflib.verbosity(0) self._fichier = fichier @@ -28,11 +35,14 @@ class donnee(object) : ceflib.close() self.donnee = {} + self.varattributes = {} self._get_time() - + def _get_time(self) : ceflib.read(self._fichier) self.donnee['epoch'] = np.copy(np.array([datetime.utcfromtimestamp(ceflib.milli_to_timestamp(i)) for i in ceflib.var('epoch')])) + self.varattributes['epoch'] = {} + for attribute in ceflib.vattributes('epoch') : self.varattributes['epoch'][attribute] = ceflib.vattr('epoch', attribute) ceflib.close() def get_time(self) : @@ -42,14 +52,29 @@ class donnee(object) : def get_varnames(self) : return self._varnames - def _get_data(self, var = 'all') : - #self.donnee = {} - ceflib.read(self._fichier) - if var == 'all' : - for nom in self._varnames: - if nom not in ['epoch', 'K001', 'K002'] : self.donnee[nom] = np.copy(ceflib.var(nom)) - else : self.donnee[var] = np.copy(ceflib.var(var)) - ceflib.close() + #def _get_data(self, var = 'all') : + ##self.donnee = {} + #ceflib.read(self._fichier) + #if var == 'all' : + #for nom in self._varnames: + #if nom not in ['epoch', 'K001', 'K002'] : self.donnee[nom] = np.copy(ceflib.var(nom)) + #else : self.donnee[var] = np.copy(ceflib.var(var)) + #ceflib.close() + + def _modify_fillval_to_nan(self, var): + self.donnee[var][self.donnee[var] == float(self.varattributes[var]['FILLVAL'])] = np.nan + return self.donnee[var] + + def _load_var(self, var) : + if var not in self.donnee : + if var in self._varnames : + ceflib.read(self._fichier) + self.donnee[var] = np.copy(ceflib.var(var)) + self.varattributes[var] = {} + for attribute in ceflib.vattributes(var) : self.varattributes[var][attribute] = ceflib.vattr(var, attribute) + ceflib.close() + if 'FILLVAL' in self.varattributes[var] : self._modify_fillval_to_nan(var) + else : print "!!!Input Variable Not Found!!!" def get_data(self, var) : self._load_var(var) @@ -99,11 +124,16 @@ class donnee(object) : return out2 def get_mean(self, t1, t2, var) : + """ + # NOTE: a.mean() when a posseses NaN return NaN -> treatment required + get mean value of an interval. NaN treatment implemented + """ self._load_var(var) i1 = abs(self.donnee['epoch'] - t1).argmin() i2 = abs(self.donnee['epoch'] - t2).argmin() - return self.donnee[var][i1:i2].mean() - + return self.donnee[var][i1:i2][~np.isnan(self.donnee[var][i1:i2])].mean() + #return self.donnee[var][i1:i2].mean() + def getIndex(self, t = None): if t : return abs(self.donnee['epoch'] - t).argmin() else : return None @@ -111,26 +141,26 @@ class donnee(object) : def get_projected(self, varz, matrice, new_coord = ['l', 'm', 'n']) : """ !!! Only works with 3D array at the moment - user must provide a transport matrix, with L,M,N in each colums of the matrix, so that : - new_vect = old_vect([x,y,z]) * matrix, or new_vect = inv(matrix) * old_vect[[x], [y], [z]] + user must provide a transport matrix, with L,M,N in each colums of the matrix, so that: + mat[:,0] = l, mat[:,1] = m, mat[:,2] = n, with mat being a numpy.mat object + and new_vect = old_vect([x,y,z]) * matrix, or new_vect = inv(matrix) * old_vect[[x], [y], [z]] + Example : l = numpy.array([ 0.16631492, -0.08297703, 0.98257527]) + m = numpy.array([-0.4464368 , -0.89481517, 0. ]) + n = numpy.array([ 0.87922325, -0.43865776, -0.18586511]) + matrice = numpy.mat([l,m,n]).transpose() """ field = self.get_data_as_field(varz) f = np.array([vect * matrice for vect in field])[:,0] self.donnee[new_coord[0]] = f[:,0] self.donnee[new_coord[1]] = f[:,1] self.donnee[new_coord[2]] = f[:,2] - - def _load_var(self, var) : - if var not in self.donnee : - ceflib.read(self._fichier) - self.donnee[var] = np.copy(ceflib.var(var)) - ceflib.close() - + return {new_coord[i] : self.donnee[new_coord[i]] for i in range(3)} def addField(self, field, abscisse = None, name = 'New_interpolated_field') : """ Add new field in the donnee object, use old abscisse to interpolate on current absisse If no abscisse is given... do something? + TO DO : modify self._varnames ???? """ def interp(xnew, xold, yold) : def tofloat(d) : @@ -143,7 +173,6 @@ class donnee(object) : self.donnee[name] = interp(self.donnee['epoch'], abscisse, field) return self.donnee[name] - def get_interval(self, var, t1 = None, t2 = None, include = False) : """ get part of the data set, between start time t1 and stop time t2 @@ -175,7 +204,6 @@ class donnee(object) : if show : plt.show() - def time_interpolation(x1, x2, y2) : """ Time interpolation for 2D data : @@ -194,8 +222,35 @@ def time_interpolation(x1, x2, y2) : return y3 +def pickleSave(filename = 'pickleSave', data = ''): + """ + Input filename WITHOUT any extension + a .save extension will be added to the name + This routine DO NOT allow overwriting!! + """ + if os.path.exists(filename) : + i = 0 + while os.path.exists(filename) : + filename = filename + str(i) + i+=1 + filename = filename + '.save' + print 'Saving data into ', filename + print 'FYI current working directory is ', os.get_cwd() + with open(filename, 'w') as sv : pickle.dump(data, sv) + return data - +def pickleLoad(filename = 'pickleSave', filename_has_extension = False): + + if not filename_has_extension : + print 'Assume file extension is not provided, adding .save to the filename' + filename = filename + '.save' + + if os.path.exists(filename): + print 'returning data' + with open(savename, 'rb') as sv: return pickle.load(sv) + else : + print 'file not found, please check filename and extension ', filename + #class field(np.ndarray) : -- libgit2 0.21.2