Source code for irfpy.ica.data_class

## -*-coding:utf-8 -*
'''
Description:

A wrapper class handling RPC-ICA level 1 data. Centralises various attributes from the initial matlab files under a single class, as well as new attributes for a broader use (see below). The class includes methods for various purposes, and aims to ease the filtering of different modes, software versions and post-acceleration values. Filtering is however optional.
The object will contain all data available for one day. It is not meant to span over several days (because of the size of the arrays) and do not allow a shorter selection.

To create one object::

    date = '20141128'
    ICADataSet = ICAData(date, dataPath='/<bli>/<bla>/')	## (Absolute) path to the data directory, e.g. containing /level1.

A more detailed example is given in pltSpec.py


Arguments:

day    Day of interest.
dataPath (default rcDataPath)    Where data can be found. Should lead to the top level of the data tree - ex ‘/…/processed’ - optional - string
modeSel (default 'mostPresent’)   Desired mode(s) - optional - float
swSel (default 'accordingToMode’)    Desired software version(s) - optional - float
paccSel (default 'all’)    Desired pacc value(s) - optional - float
cleanCT (default True)    If cross-talk cleaning wanted - optional - boolean


Dependencies:

irfpy.ica.modules
irfpy.ica.io.readproc()

'''

import sys
import time
import datetime
from collections import Counter

import numpy as np
import matplotlib as mpl

#Qt5Agg breaks in Jenkins. Ignore error.
try:
    mpl.use(' Qt5Agg')
except:
    pass


from matplotlib.dates import date2num, num2date
from irfpy.util.irfpyrc import Rc
import irfpy.ica.io as icaio




font = {'family': 'sans serif',
        'serif': 'Helvetica',
        'size': 20}
mpl.rc('font', **font)
mpl.rcParams['pdf.fonttype'] = 42  # For exporting pdf plots with editable texts.

rcDataPath = Rc().get('ica', 'dataPath')
rcOwnPath = Rc().get('ica', 'ownPath')
rcIcaPyPath = Rc().get('ica', 'icaPyPath')
rcSPICEPath = Rc().get('ica', 'SPICEPath')
rcMAGPath = Rc().get('ica', 'MAGPath')
rcProductsPath = Rc().get('ica', 'productsPath')
rcNpyPath = Rc().get('ica', 'npyPath')
rcOwnNpyPath = Rc().get('ica', 'ownNpyPath')
rcSelec5DFile = Rc().get('ica', 'selec5DFile')
rcICAGFFile = Rc().get('ica', 'ICAGFFile')
rcSelecFile = Rc().get('ica', 'selecFile')
rcSummaryFile = Rc().get('ica', 'summaryFile')
rcCTFile = Rc().get('ica', 'CTFile')



[docs]class ICAData: def __init__(self, day, dataPath=rcDataPath, modeSel='mostPresent', swSel='accordingToMode', paccSel='all', cleanCT=True): r''' day [necessary] = 'yyyymmdd' ; dataPath [optionnal but...] = '/<bli>/<bla>/ICAPy/' ; the absolute/relative path to the ICAPy directory. ''' if dataPath[-1] != '/': dataPath = dataPath + '/' self.date = day # ICAData will load and handle all data of this precise date, only. print('\nICAData___________________________________________________________________') print('|\n| - ICAPy -') print('| - {} {} {} -\n|'.format(day[:4], day[4:6], day[6:8])) ########################################################################## ## We load the hourly mat files and declare the unfiltered ICAData attributes. self.variables = ['orig_ionspectra', 'time_instances', 'decoder_error_flag', 'E', 'elevation_rad', 'cur_pacc', 'mode', 'sw_version', 'version_list'] # All variables of interest for this class. print('| Loading level 1 data ... ', end='') tic = time.time() ## Next line: the data are actually loaded. dataDict = icaio.readproc(dataPath, from_day=day, variables=self.variables, verbose=False, flat=False, includeerrors=False, partialelevationscan=False) print('{:.2f}s\n|'.format(time.time() - tic)) if not bool(dataDict): ## Nothing came out of readproc: no data for this path/date. sys.exit('No data loaded from {}\n'.format(dataPath)) ## try: self.selec = np.load(selecFile).item()[self.date] ## Manual species selection. ## except FileNotFoundError: print('| Selection file not found.') ## The following attributes are unfiltered. They are the exact out put of loadlevel1. ## len = nb time instances; datetime format (dt). self.timeDtAll = dataDict['time_instances'] ## ICAData.ionSpectra: 16 * 96 * 32 * nb time instances , changed in a few lines. self.ionSpectraAll = dataDict['orig_ionspectra'] self.modeAll = dataDict['mode'] # Mode, for each time instance. self.swAll = dataDict['sw_version'] # Software version used, for each time instance. self.paccAll = dataDict['cur_pacc'] # Pacc value for each time instance. From 0 to 7. self.EAll = dataDict['E'] # Best estimations for the energy tables. self.elevTableAll = dataDict['elevation_rad'] # Physical elevation tables. ######################## ########################################################################## ## ICAData.ionSpectra has the time axis as FIRST axis. Yes. ## Decoder errors are deleted here, systematically. It is assumed that ## decoder errors are always over an entire 16 elevation 'scan'. ## nb time instances * 16 * 96 * 32 self.ionSpectraAll = np.transpose(self.ionSpectraAll, (3, 0, 1, 2)) ## i.e. an incomplete angular scan was introduced in the data. if self.ionSpectraAll.shape[0] % 16 != 0: sys.exit( '\n\n The file contains an incomplete angular scan\n (the number of time instances is not a multiple of 16).\n\n') ########################################################################## ## Cleaning sector cross-talk. if cleanCT: self.cleanCT() ######################## ## ICAData.ionSpectraScanAll = nbScanx16(el)x16(az)x96x32 self.ionSpectraScanAll = np.reshape(self.ionSpectraAll, (-1, 16, 16, 96, 32)) ## msp stands for mode/software/postacc. self.mspAll = self.modeAll * 1000 + self.swAll * 10 + self.paccAll ## All msp values operated during the date. self.mspList = np.array(list(Counter(self.mspAll))) self.mspCounter = Counter(self.mspAll) # All msp values, with the number of occurences. ######################## ########################################################################## # The following attributes are the filtered ones. self.ionSpectra = None # The spectra after filtering. self.timeDt = None # The selected time instances after filtering. self.mode = None # The selected mode after filtering. self.sw = None # The selected software version after filtering. self.pacc = None # The selected pacc value after filtering. self.E = None # Energy table actually used by ICAData after filtering. self.elevTable = None # Physical elevation table actually used by ICAData after filtering. self.azimTable = None # Angle of the center of each azimuth anode, or sector. RADIAN. self.msp = None # Selected combinations of modes software versions and pacc values. ######################## ########################################################################## # self.modeTable[i] gives for the ith mode: [Masses, Azimuth angles, # Energies, Elevation angles] self.modeTable = np.array([[0, 0, 0, 0], [0, 0, 0, 0], [2, 1, 32, 1], [0, 0, 0, 0], [6, 1, 96, 1], [32, 1, 96, 1], [0, 0, 0, 0], [0, 0, 0, 0], [6, 16, 96, 16], [6, 16, 96, 8], [6, 16, 96, 4], [6, 16, 96, 2], [6, 8, 96, 2], [6, 4, 96, 2], [3, 4, 96, 2], [3, 4, 96, 1], [16, 16, 96, 16], [16, 16, 96, 8], [16, 16, 96, 4], [8, 16, 96, 4], [4, 16, 96, 4], [2, 16, 96, 4], [2, 8, 96, 4], [2, 8, 96, 4], [32, 16, 96, 16], [32, 16, 96, 8], [32, 16, 96, 4], [32, 16, 96, 2], [32, 8, 96, 2], [32, 4, 96, 2], [32, 2, 96, 2], [32, 2, 96, 1]]) ######################## ########################################################################## ## # Filtering the data. ## self.filterData(modeSel=modeSel, swSel=swSel, paccSel=paccSel) ## ## ######################## ########################################################################## ## # Angles attributes. ## # Number of sector, even though this will most likely never change. Binned # data are expanded back to 16 sectors. self.nbSec = 16 self.deltaAzim = 2. * np.pi / self.nbSec # Angular width of one sector. RADIAN. # Angle of the center of each azimuth anode, or sector. RADIAN. self.azimTable = self.indToAz(np.arange(16)) self.elevLim = 45. * np.pi / \ 180. # Limit in elevation of what we want to consider. Idealy, should be the same than the one used to create the elevation table. RADIAN. self.nbElev = 17 # Number of elevation bins in the re-binned data. # Angular width of the above elevations bins. RADIAN. self.deltaElev = 2 * self.elevLim / (self.nbElev - 1) # See comment just bellow: self.elevTableNew = np.linspace(-self.elevLim, self.elevLim, self.nbElev) '''Since the given elevation table is not necessarely "alligned" (elevation/energy dependency), one cannot simply use the indices along elevation dimension and calculate a physical elevation from it (as done with azimuth anodes, because anodes are always at the same physical angle, no matter the energy!) So we define a new "artificial" elevation table, and will sort all (elevation, energy) combination according to this new table, ICAData.elevTableNew . This is in fact a re-binning.''' ## ## ######################## ########################################################################## ## # Time attributes. ## if self.timeDt.size == 0: # The filtering removed all time instances, nothing is left. self.nbTime = 0 self.nbScan = 0 print('|\n|\n| No data left after filtering.\n|\n|') print('|_________________________________________________________________________\n') else: self.timeNum = date2num(self.timeDt) # python timedate ordinal format (float). # Important value 366: matlab datenum start the 00-00-0000 whereas python # ordinal start the 01-01-0001. self.day = num2date(self.timeNum[-1]) # Day in datetime format. # Date of the first time instance as a string. self.dayString = num2date(self.timeNum[-1]).strftime('%FT%T')[:10] self.nbTime = self.timeDt.size # Number of time instances. # Data are expanded to 16 elevation bins during pre-processing, whatever the mode. self.nbScan = self.nbTime // 16 print('|\n| {} time instances after filtering.\n|'.format(self.nbTime)) print('|_________________________________________________________________________\n') ## ## ######################## ########################################################################## ## # Spectrogram, energy-masse matrix and sector-masse matrix, # direct subproducts of ICAData.ionSpectra ## # Division of the orig_ionspecrtra in full angular scans: nbScanx16(el)x16(az)x96x32 self.ionSpectraScan = self.ionSpectra.reshape(-1, 16, 16, 96, 32) # ## Now ICAData.ionSpectraScan.shape = nbScanx16(az)x16(el)x96x32 # Spectrogram, ICAData.spectrogram: nbTimex96 # self.ionSpectrogram = np.sum(np.sum(self.ionSpectra, axis=3), axis=1) # # Energy mass matrix integrated over the whole data set, ICAData.EMTot: 96x32 # self.EMTot = np.sum(np.sum(self.ionSpectra, axis=1), axis=0) # # Energy mass matrix integrated over each FOV, ICAData.EMScan: nbScanx96x32 # self.EMScan = np.sum(np.sum(self.ionSpectraScan, axis=2), axis=1) # # Sector versus Mass ring integrated over the whole data set, ICAData.SMTot: 16x32 # self.SMTot = np.sum(np.sum(self.ionSpectra, axis=2), axis=0) # # Sector versus Mass ring , ICAData.SMScan: nbScanx6x32 # self.SMScan = np.sum(np.sum(self.ionSpectraScan, axis=3), axis=1) # # Energy versus Sector integrated over the whole data set, ICAData.SMTot: 16x96 # self.ESTot = np.sum(np.sum(self.ionSpectra, axis=3), axis=0) # # Energy versus Sector, ICAData.SMScan: nbScanx16x96 # self.ESScan = np.sum(np.sum(self.ionSpectraScan, axis=4), axis=1) self.nbScan = self.ionSpectra.shape[0]/16 self.timeDtScan = self.timeDt[8::16] ## ## ######################## ########################################################################## ## # Field of View related attributes. ## # Species selection informations, defined as a selection window in EM space. self.EMWindow = {'H+': np.array([[34, 54], [26, 27]])} # The boundaries are included in the selection [[ELow,EHigh],[massLow,massHigh]]. # ICAData.FOVSpecies[s].shape = ICAData.nbScan x 16(az) x ICAData.nbElev(newElBin) self.FOVSpecies = {} self.FOVEMean = {} # ICAData.specSpecies[s].shape = ICAData.nbScan x 16(az) x 16(elBin) x 96 x 32 self.specSpecies = {} ## ## ########################
[docs] def filterData(self, modeSel='mostPresent', swSel='accordingToMode', paccSel='all'): ################################################ # What do we have in the data set? print('| Modes: ', end='') for k in Counter(self.modeAll).keys(): print('{} ({}), '.format(k, Counter(self.modeAll)[k]), end="") print('') print('| SW: ', end="") for k in Counter(self.swAll).keys(): print('{} ({}), '.format(k, Counter(self.swAll)[k]), end="") print('') print('| Pacc: ', end="") for k in Counter(self.paccAll).keys(): print('{} ({}), '.format(k, Counter(self.paccAll)[k]), end="") print('') print('|', self.mspCounter, '\n|') ######################## ################################################ # Modes. if modeSel == 'mostPresent': # Only the most present mode in the data set is considered. self.mode = max(Counter(self.modeAll), key=Counter(self.modeAll).get) indMode = (self.modeAll == self.mode) print('| Mode {} is selected, all other time instances are ignored.'.format(self.mode)) elif modeSel == 'all': indMode = np.arange(self.timeDtAll.size)[0] self.mode = 'all' print('| All modes are selected.') else: try: # Several desired mode. indMode = np.zeros(self.timeDtAll.size) for m in modeSel: indMode += (self.modeAll == m) self.mode = modeSel # Here it might be dangerous. If the user choose two different modes that have not the same # number of elevation for example. Should be enhanced according to a # precise scientific goal. print('| Modes {} are selected, all other time instances are ignored.'.format(modeSel)) except TypeError: # Only one desired mode is considered. indMode = (self.modeAll == modeSel) self.mode = modeSel print('| Mode {} is selected, all other time instances are ignored.'.format(self.mode)) ######################## ################################################ # Software version. # The software version is chosen according to the selected mode. if swSel == 'accordingToMode': c = Counter(self.swAll[indMode.astype(np.int8)]) if len(list(c)) == 1: # Only one sw for the indMode self.sw = list(c)[0] indSw = (self.swAll == self.sw) self.E = self.EAll[self.sw - 1] # Now only one energy table and # one elevation table are considered. self.elevTable = self.elevTableAll[self.sw - 1] print('| Software version {} is selected, all other time instances are ignored.'.format(self.sw)) else: # swMax = max(c, key=c.get) sys.exit('\n|\n| Multiple software versions for the selected modes: {}\n|\n|'.format(list(c))) else: # One desired software version is considered. swSel = int(swSel) self.sw = swSel indSw = (self.swAll == swSel) self.E = self.EAll[swSel - 1] # Now only one energy table and self.elevTable = self.elevTableAll[swSel - 1] # one elevation table are considered. print('| Software version {} is selected, all other time instances are ignored.'.format(swSel)) # Here, the software version is chosen according to the selected mode. # If the scientific need requires to consider different sw versions, the same thing used for # the mode can be done, but different energy and elevation tables as to be taken # into account in the analysis. Cannot be dealt with the current version of the class. ######################## ################################################ # pacc value. if not paccSel == 'all': self.pacc = paccSel indPacc = self.paccAll == paccSel if paccSel < 3: print('| Low pacc value is selected, all other time instances are ignored.'.format(paccSel)) elif paccSel > 2 and paccSel < 6: print('| Medium pacc value is selected, all other time instances are ignored.'.format(paccSel)) elif paccSel > 5: print('| High pacc value is selected, all other time instances are ignored.'.format(paccSel)) else: sys.exit('| Unvalid pacc value: {}'.format(paccSel)) elif paccSel == 'all': self.pacc = 'all' indPacc = 1 # Will be broadcasted. print('| All pacc values are selected.') ######################## ################################################ # 01-01-01. Who did, and why?.. indWhy = (self.timeDtAll == datetime.datetime(1, 1, 1)) # Quick fix, but of course self.timeDt[0] should be 01-01-01... self.timeDtAll[indWhy] = self.timeDtAll[0] ######################## # We put in indToKeep all the indices we want to keep (the intersection of # all desired indices). indToKeep = indMode * indSw * indPacc # And finally, the 'scientific' time series (ICAData.ionSpectra, # ICAData.timeDt) are reduced to these desired indices. self.ionSpectra = self.ionSpectraAll[np.where(indToKeep)] self.timeDt = self.timeDtAll[np.where(indToKeep)] # We also note the used msp. if np.sum(indToKeep) == 0: self.msp = None print('| No filtered msp.') else: self.msp = np.array(list(Counter(self.mspAll[np.where(indToKeep)]))) print('| Filtered msp: {}'.format(self.msp))
[docs] def cleanCT(self): ''' Removing electronic cross-talk between sectors, empirical approach. With E the counts in the emitter sector, R the counts in the receiver sector: R -= (a*ln(E)+b) * E ''' print('| Cleaning sector cross-talk.\n|') ct = np.nan_to_num(np.load(rcCTFile)) # .item()['LogCrosstalk'] eR = np.array([[5, 5, 6, 6, 7, 7, 7, 8, 9, 9, 10, 10, 11, 11, 12, 12, 14, 14, 15, 15], [0, 1, 0, 2, 0, 2, 3, 0, 0, 1, 0, 2, 0, 2, 4, 0, 6, 2, 6, 2]]) # Emitters, receivers. azElScan = np.nansum(self.ionSpectraAll, axis=(2, 3)) # Sum over mass and energy. for e, r in eR.T: self.ionSpectraAll[:, r] -= (ct[0, r, e] * np.log(azElScan[:, e].clip(1e-10)) + ct[1, r, e])[:, None, None] * self.ionSpectraAll[:, e] # R -= (a*ln(E)+b) * E self.ionSpectraAll = self.ionSpectraAll.clip(0)
[docs] def alignElev(self): ''' Produces a new ionSPectraScan array with aligned elevation bins: for a same elevation index, all energies correspond to one physical elevation, and this new elevation table is given in ICAData.elevTableNew ''' print('ICAData___________________________________________________________________\n| Aligne elevation ... ', end='') tic = time.time() self.ionSpectraScanAligned = np.zeros((self.nbScan, self.nbElev, 16, 96, 32)) elevT = self.elevTable # Not NaN anymore, but unreachable value (not covered by ICAData.elevTablNew) elevT[np.isnan(elevT)] = 100 # ICAData.nbElev loops. nel stands for "new elevation". for i, nel in enumerate(self.elevTableNew): # This second mask select the count with the right elevation. elM stands # for elevation mask. elM = (elevT >= (nel - self.deltaElev * .5)) & (elevT < (nel + self.deltaElev * .5)) self.ionSpectraScanAligned[:, i] = np.sum( self.ionSpectraScan * elM.T[None, :, None, :, None], axis=1) print(' {:.2f}s'.format(time.time() - tic)) print('|_________________________________________________________________________\n')
[docs] def initFOVEM(self, species=['H+']): """ WARNING: initFOVEM integrate 16 time instances to create 1 FOV for each species. SO using ICAData.initFOVEM without considering the different modes operated for the data set considered is not appropriate. ICAData.FOVSpecies can be initialized with any data set since any binned or absent azimuth/elevation is expanded again in the process. Modes 8, 16, 24 correspond to the highest angular resolution (16x16). Modes 9, 17, 25 have binned elevation (from 16 down to 8), and maximum azimuth resolution (16). See Nilsson et al., 2007, ssr, p.682-683 """ print('ICAData___________________________________________________________________\n| Initialising FOVs ... ', end='') tic = time.time() if self.mode == 'all': sys.exit('| ICAData does not initialise FOVs if all modes are selected.\n|') elevT = self.elevTable # Not NaN anymore, but unreachable value (not covered by ICAData.elevTablNew) elevT[np.isnan(elevT)] = 100 for s in species: # ~2 loops. # ionSScan.shape = nbFOVx16(az)x16(el)x96x32 ionSScan = np.transpose(self.ionSpectraScan.copy(), (0, 2, 1, 3, 4)) self.FOVSpecies[s] = np.zeros((self.nbScan, 16, self.nbElev)) # nbScan x azim x nbElev # Mask for energy. eM = (np.arange(96) >= self.EMWindow[s][0, 0]) & ( np.arange(96) <= self.EMWindow[s][0, 1]) # Mask for mass. mM stands for mass Mask. mM = (np.arange(32) >= self.EMWindow[s][1, 0]) & ( np.arange(32) <= self.EMWindow[s][1, 1]) # eMM stands for energy Mass Mask. eMM = eM[None, None, None, :, None] * mM[None, None, None, None, :] # This first 'mask' m takes care of energy range and mass selection: # counts out of this selection are put to 0. ionSScan *= eMM # ionSScan = ionSScan[:,:,:,self.EMWindow[s][0,0]:self.EMWindow[s][0,1], self.EMWindow[s][1,0]:self.EMWindow[s][1,1]] # ICAData.nbElev loops. nel stands for "new elevation". for i, nel in enumerate(self.elevTableNew): # This second mask select the count with the right elevation. elM stands # for elevation mask. elM = (elevT >= (nel - self.deltaElev * .5)) & (elevT < (nel + self.deltaElev * .5)) self.FOVSpecies[s][:, :, i] = np.sum(np.sum(ionSScan[:, :, elM.T], axis=3), axis=2) print(' {:.2f}s'.format(time.time() - tic)) print('|_________________________________________________________________________\n')
[docs] def initFOVEMScan(self, masks): """ WARNING: initFOVEMScan integrate 16 time instances to create 1 FOV for each species. SO using ICAData.initFOVEMScan without considering the different modes operated for the data set considered is not appropriate. ICAData.FOVSpecies can be initialized with any data set since any binned or absent azimuth/elevation is expanded again in the process. Modes 8, 16, 24 correspond to the highest angular resolution (16x16). Modes 9, 17, 25 have binned elevation (from 16 down to 8), and maximum azimuth resolution (16). See Nilsson et al., 2007, ssr, p.682-683 """ if self.mode == 'all': sys.exit('ICAData does not initialise FOVs if all modes are selected.') elevT = self.elevTable # Not NaN anymore, but unreachable value (not covered by ICAData.elevTablNew) elevT[np.isnan(elevT)] = 100 # ionSScan.shape = nbFOVx16(az)x16(el)x96x32 ionSScan = np.transpose(self.ionSpectraScan, (0, 2, 1, 3, 4)) for s in list(masks): # ~2 loops. self.FOVSpecies[s] = np.zeros((self.nbScan, 16, self.nbElev)) # nbScan x azim x nbElev ss = ionSScan * masks[s] # [:,None,None,:] self.specSpecies[s] = ss for i, nel in enumerate(self.elevTableNew): # 17 loops. nel stands for "new elevation". # This second mask select the count with the right elevation. elM stands # for Energy Mass mask elM = (elevT >= (nel - self.deltaElev * .5)) & (elevT < (nel + self.deltaElev * .5)) self.FOVSpecies[s][:, :, i] = np.sum(np.sum(ss[:, :, elM.T], axis=3), axis=2)
[docs] def initSpecSpecies(self, species=['H+']): """ Initialise a 5D matrix for the specified species, with elevation aligned. """ print('ICAData___________________________________________________________________\n| Initialising specSpecies ... ', end='') tic = time.time() if self.mode == 'all': sys.exit('| ICAData does not initialise FOVs if all modes are selected.\n| ') elevT = self.elevTable # Not NaN anymore, but unreachable value (not covered by ICAData.elevTablNew) elevT[np.isnan(elevT)] = 100 for s in species: # ~2 loops. # ionSScan.shape = nbFOVx16(az)x16(el)x96x32 ionSScan = np.transpose(self.ionSpectraScan.copy(), (0, 2, 1, 3, 4)) # Mask for energy. eM = (np.arange(96) >= self.EMWindow[s][0, 0]) & ( np.arange(96) <= self.EMWindow[s][0, 1]) # Mask for mass. mM stands for mass Mask. mM = (np.arange(32) >= self.EMWindow[s][1, 0]) & ( np.arange(32) <= self.EMWindow[s][1, 1]) # eMM stands for energy Mass Mask. eMM = eM[None, None, None, :, None] * mM[None, None, None, None, :] # This first 'mask' m takes care of energy range and mass selection: # counts out of this selection are put to 0. ionSScan *= eMM self.specSpecies[s] = np.zeros((self.nbScan, 16, self.nbElev, 96, 32)) # ICAData.nbElev loops. nel stands for "new elevation". for i, nel in enumerate(self.elevTableNew): # This second mask select the count with the right elevation. elM stands # for elevation mask. elM = (elevT >= (nel - self.deltaElev * .5)) & (elevT < (nel + self.deltaElev * .5)) ss = ionSScan * elM.T[None, None, :, :, None] # Extremely slow. self.specSpecies[s][:, :, i, :, :] = np.sum(ss, axis=2) print('{:.2f}s'.format(time.time() - tic)) print('|_________________________________________________________________________\n')
[docs] def initSpectraSpecies(self, species, msp=None): r''' Returns an array of same shape as ICAData.ionSpectraScan. The counts outside the manual selection for the specified species are put to 0. ''' if msp is None: msp = self.msp[0] indScan = (self.mspAll[::16] == msp) spec = np.zeros(self.ionSpectraScanAll[indScan].shape) e1, e2, m1, m2 = self.selec[msp][species].flatten() spec[:, :, :, e1:e2, m1:m2] = self.ionSpectraScanAll[indScan, :, :, e1:e2, m1:m2] return spec
########################################################################## ## # Very important methods here, allowing to translate angles, directions, indexes. ## # They sit here because they need ICAData angle attributes. Other wise, non-sens. ## # WARNING! Example: index 8 (azimuth) corresponds to the center of the 8th sector, # and the value 8.75 (real values used for the Sun and comet azimuth for instance) # would correspond to the first quarter of the 9th sector. Sketches are the key! ##
[docs] def indToAz(self, ind): '''az=0 inbetween sectors 8 and 9, along +Xica = +Xsc.''' return ((ind + self.nbSec / 2. - .5) * self.deltaAzim) % (2 * np.pi) # Here the -.5 means that an integer index corresponds to the
# center of a sector: we have to remove half a sector-angle to get the real angle. # i.e. with 16 sectors we can't get an angle equal to zero.
[docs] def indToEl(self, ind): ''' From an elevation index (as binned previously) [0,nbElev] to [-elevLim,elevLim] ''' # return (ind - self.nbElev/2.)*self.deltaElev Old way to translate, but wrong now! An even number of elevation bins would make # it possible to have elev=0°, which is impossible (symmetry), and an odd number of bins couldn't have the physical elevation 0°! # A +.5 is needed in the parenthesis, but this is useless since we use an "artificial" elevation table now. # We now just have to interpolate this new aertificial elevation table. # One possibility: np.interp(ind,np.arange(self.nbElev),self.elevTableNew) return (ind - self.nbElev / 2. + .5) * self.deltaElev
[docs] def azToInd(self, azim): ''' azim physical angle in RADIAN, in [0 , 2*pi] ''' return (azim / self.deltaAzim + self.nbSec / 2. + .5) % self.nbSec # Here the +.5 means that an integer index corresponds to the
# center of a sector: we have to add half a sector-index to get the real index. # i.e. with 16 sectors we can't get an angle equal to zero.
[docs] def elToInd(self, elev): ''' elev physical angle in RADIAN, in [-pi/4 , pi/4] ''' return elev / self.deltaElev + self.nbElev / 2.
[docs] def indToVec(self, secInd, elevInd): ''' From (secInd,elevInd) to a unit VIEWING DIRECTION vector in ICA reference frame ''' # From indexes to physical angles elev = self.indToEl(elevInd) azim = self.indToAz(secInd) # From angles to vector, in ICA reference frame. try: V = np.zeros((3, elevInd.shape[0], secInd.shape[1])) except: try: V = np.zeros((3, secInd.shape[0])) except: V = np.zeros(3) # azim in [0,2*pi], trigonometric circle in (x,y). Thus, cartesian # components computation is straightforward. V[0] = np.cos(elev) * np.cos(azim) V[1] = np.cos(elev) * np.sin(azim) V[2] = np.sin(elev) # elev in [-pi/4,pi/4], no problem. return V
[docs] def indToVecRaw(self, secInd, elevInd): ''' From (secInd,elevInd) to a unit VIEWING DIRECTION vector in ICA reference frame, for unaligned elevation. ''' # From indexes to physical angles elev = self.elevTable[70, elevInd] # 40 is an arbitrary medium energy index. azim = self.indToAz(secInd) # From angles to vector, in ICA reference frame. try: V = np.zeros((3, elevInd.shape[0], secInd.shape[1])) except: try: V = np.zeros((3, secInd.shape[0])) except: V = np.zeros(3) # azim in [0,2*pi], trigonometric circle in (x,y). Thus, cartesian # components computation is straightforward. V[0] = np.cos(elev) * np.cos(azim) V[1] = np.cos(elev) * np.sin(azim) V[2] = np.sin(elev) # elev in [-pi/4,pi/4], no problem. return V
[docs] def azElToVec(self, azim, elev): ''' From (azim,elev) in RADIAN to a unit VIEWING DIRECTION vector in ICA reference frame.''' print('To be refurbished!') sys.exit() try: V = np.zeros((3, elev.shape[0], azim.shape[1])) except: try: V = np.zeros((3, azim.shape[0])) except: V = np.zeros(3) print(V.shape) # azim in [0,2*pi], trigonometric circle in (x,y). Thus, cartesian # components computation is straightforward. V[0] = np.cos(elev) * np.cos(azim) V[1] = np.cos(elev) * np.sin(azim) V[2] = np.sin(elev) # elev in [-pi/4,pi/4], no problem. return V
## ## ##########################################################################