##// END OF EJS Templates
schain
RSS
Se adiciona script Julia_EEJ.py para test y la operaciones GetSNR
Alexander Valdez -
r1684:1f582a059484
parent child
Show More
Side by Side Unified
Context file:
r1684:1f582a059484 Loading diff...:
Show file before | Show file after | Hide comments
@@ -0,0 +1,291
1 #!/usr/bin/env python
2 '''
3 Created on Jul 7, 2014
4
5 @author: roj-idl71
6 '''
7 import os, sys, json
8
9 #path = os.path.dirname(os.getcwd())
10 #path = os.path.join(path, 'source')
11 #sys.path.insert(0, path)
12
13 from schainpy.controller import Project
14
15 if __name__ == '__main__':
16
17 desc = "JULIA raw experiment "
18 filename = "schain.xml"
19
20 dpath = '/home/roberto/puma/JULIA_NEW/JULIA_EW/D2022/'
21 dpath = '/home/roberto/puma/JULIA_NEW/JULIA_EW/D2021/'
22 dpath = '/home/roberto/puma/JULIA_NEW/JULIA_EW/D2017/'
23 dpath = '/home/roberto/puma/JULIA_NEW/JULIA_EW/D2021/'
24 dpath = '/home/roberto/puma/JULIA_NEW/JULIA_EW/D2018/'
25 #dpath = '/home/roberto/puma/JULIA_NEW/JULIA_EW/D2017/'
26 #dpath = '/home/roberto/Folder_aux/D2019/'
27
28 ##dpath = '/home/roberto/puma/JULIA_EW_IMAGING/JULIA_EW/D2017/'
29 dpath = '/home/soporte/PUMA/JULIA_EW_IMAGING/JULIA_EW/D2017/'
30
31 figpath = '/home/soporte/DATA/Pictures/JULIA/EEJ/Skew_but_dop_is_shift'+'/'+dpath[-5:-1]
32 ppath = "/home/soporte/DATA/MLT/Oblique/2022_03/data_reshape"
33 spcdata = '/home/soporte/DATA/JULIA/EEJ/SPC'
34 path_parameters = '/home/soporte/DATA/JULIA/EEJ/Params'
35 db_range=['25','35']
36 db_range=['10','20']
37 #db_range=['14','20']
38 db_range=['10','23']
39 db_range=['13','20']
40 db_range=['21','30']
41 db_range=['15','30']
42 db_range=['13','28']
43 tiempo=['7','18']
44 altura1=[2,20]
45 altura1=[90,220]
46 velocity=['-80','80']
47 period=60
48 # PROJECT 1
49
50 show_spc = 0
51 save_spc = 0
52 fitting = 1
53 save_params = 1
54 plot_params = 0
55
56 controllerObj = Project()
57 controllerObj.setup(id = '191', name='altura1', description=desc)
58
59 readUnitConfObj1 = controllerObj.addReadUnit(datatype='SpectraReader',
60 path=dpath,
61 startDate='2017/09/01', #Check 21-29 Jun 2021, 18 May 2022
62 endDate='2017/09/30', #05,06,07-01-18
63 #startTime='06:00:00',
64 #endTime='18:00:00',
65 startTime='07:00:00',
66 #startTime='10:00:00',
67 #startTime='08:38:01',
68 #startTime='11:16:00',
69 #startTime='08:13:00',
70 endTime='17:59:59',
71 #startTime='16:30:00',
72 #endTime='17:30:59',
73 online=0,
74 walk=1,
75 expLabel='150EEJ',
76 getByBlock=1,
77 delay=20)
78
79 # opObj00 = readUnitConfObj.addOperation(name='printInfo')
80 # opObj00 = readUnitConfObj.addOperation(name='printNumberOfBlock')
81
82 procUnitConfObj1 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=readUnitConfObj1.getId())
83
84 opObj11 = procUnitConfObj1.addOperation(name='selectChannels')
85 opObj11.addParameter(name='channelList', value='2,', format='intlist')
86
87 '''
88 opObj11 = procUnitConfObj1SPC.addOperation(name='removeDC')
89 opObj11.addParameter(name='mode', value='2', format='int')
90 '''
91
92 #opObj11 = procUnitConfObj1.addOperation(name='dopplerFlip') #It fixes the Doppler
93 #opObj11.addParameter(name='chann', value='0', format='int')
94
95 opObj11 = procUnitConfObj1.addOperation(name='removeInterference')
96
97 opObj11 = procUnitConfObj1.addOperation(name='IncohInt', optype='other')
98 #opObj11.addParameter(name='n', value='20', format='int')
99 opObj11.addParameter(name='n', value='1', format='int')
100
101 opObj11 = procUnitConfObj1.addOperation(name='GetSNR', optype='other')
102
103 if save_spc:
104 dataList=['data_spc',
105 'utctime']
106 metadataList=['nFFTPoints','VelRange','normFactor',
107 'heightList','timeZone']
108
109 op221 = procUnitConfObj1.addOperation(name='HDFWriter', optype='external')
110 op221.addParameter(name='path', value=spcdata)
111 #op221.addParameter(name='mode', value=1, format='int')
112 op221.addParameter(name='dataList', value=dataList)
113 op221.addParameter(name='metadataList', value=metadataList)
114 #op221.addParameter(name='blocksPerFile', value=500)
115 op221.addParameter(name='blocksPerFile', value=2000)
116
117 if show_spc:
118 #'''
119 # opObj11 = procUnitConfObj1SPC.addOperation(name='removeInterference')
120
121 opObj11 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other')
122 opObj11.addParameter(name='id', value='1', format='int')
123 opObj11.addParameter(name='wintitle', value='Oblique', format='str')
124 opObj11.addParameter(name='zmin', value=db_range[0], format='int')
125 opObj11.addParameter(name='zmax', value=db_range[1], format='int')
126 opObj11.addParameter(name='xaxis', value='velocity', format='str')
127 opObj11.addParameter(name='ymin', value=altura1[0], format='int')
128 opObj11.addParameter(name='ymax', value=altura1[1], format='int')
129 # opObj11.addParameter(name='xmin', value=velocity[0], format='int')
130 # opObj11.addParameter(name='xmax', value=velocity[1], format='int')
131 opObj11.addParameter(name='showprofile', value='1', format='int')
132 opObj11.addParameter(name='save', value=figpath, format='str')
133 #'''
134 #'''
135 opObj11 = procUnitConfObj1.addOperation(name='RTIPlot', optype='other')
136 opObj11.addParameter(name='id', value='10', format='int')
137 opObj11.addParameter(name='wintitle', value='JULIA EEJ RTI', format='str')
138 opObj11.addParameter(name='xmin', value=tiempo[0], format='float')
139 opObj11.addParameter(name='xmax', value=tiempo[1], format='float')
140 opObj11.addParameter(name='ymin', value=altura1[0], format='int')
141 opObj11.addParameter(name='ymax', value=altura1[1], format='int')
142 opObj11.addParameter(name='zmin', value=db_range[0], format='int')
143 opObj11.addParameter(name='zmax', value=db_range[1], format='int')
144 opObj11.addParameter(name='showprofile', value='1', format='int')
145 #opObj11.addParameter(name='save_period', value=40, format='str')
146 opObj11.addParameter(name='save', value=figpath, format='str')
147 #opObj11.addParameter(name='throttle', value=1, format='str')
148 #'''
149 '''
150 opObj11 = procUnitConfObj1.addOperation(name='SnrPlot', optype='other')
151 opObj11.addParameter(name='id', value='10', format='int')
152 opObj11.addParameter(name='wintitle', value='JULIA EEJ SNR', format='str')
153 opObj11.addParameter(name='xmin', value=tiempo[0], format='float')
154 opObj11.addParameter(name='xmax', value=tiempo[1], format='float')
155 opObj11.addParameter(name='ymin', value=altura1[0], format='int')
156 opObj11.addParameter(name='ymax', value=altura1[1], format='int')
157 opObj11.addParameter(name='zmin', value=0.1, format='int')
158 opObj11.addParameter(name='zmax', value=50, format='int')
159 #opObj11.addParameter(name='showprofile', value='1', format='int')
160 opObj11.addParameter(name='save', value=figpath, format='str')
161 '''
162 if fitting:
163 Dop = 'Max' #Plot and Save the Pos[Max_val] as the Doppler Shift
164 Dop = 'Shift' #Plot and Save the Skew Gaussian Shift as the Doppler Shift
165 opObj11 = procUnitConfObj1.addOperation(name='Oblique_Gauss_Fit', optype='other')
166 opObj11.addParameter(name='mode', value=9, format='int') #Skew
167 opObj11.addParameter(name='Dop', value=Dop)
168 #opObj11.addParameter(name='mode', value=11, format='int')
169
170 if save_params:
171 '''
172 dataList=['powerdB', 'Oblique_params', 'Oblique_param_errors', 'dplr_2_u', 'data_snr',
173 'utctime']
174 metadataList=['VelRange',
175 'heightList','timeZone']
176
177 op221 = procUnitConfObj1.addOperation(name='HDFWriter', optype='external')
178 op221.addParameter(name='path', value=path_parameters)
179 #op221.addParameter(name='mode', value=1, format='int')
180 op221.addParameter(name='dataList', value=dataList)
181 op221.addParameter(name='metadataList', value=metadataList)
182 #op221.addParameter(name='blocksPerFile', value=500)
183 op221.addParameter(name='blocksPerFile', value=2000)
184 '''
185
186 one = {'gdlatr': 'lat', 'gdlonr': 'lon', 'inttms': 'paramInterval'} #reader gdlatr-->lat only 1D
187
188 two = {
189 'snl': 'snl', #DeberΓ­a salir como el original pero mΓ‘s limpio
190 'RANGE': 'heightList', #<----- nmonics
191 'DOPP_T1_EEJ': ('Dop_EEJ_T1', (0)),
192 'DDOPP_T1_EEJ': ('Err_Dop_EEJ_T1', (0)),
193 'SPEC_W_T1_EEJ': ('Spec_W_T1', (0)),
194 'DSPEC_W_T1_EEJ': ('Err_Spec_W_T1', (0)),
195 'DOPP_T2_EEJ': ('Dop_EEJ_T2', (0)),
196 'DDOPP_T2_EEJ': ('Err_Dop_EEJ_T2', (0)),
197 'SPEC_W_T2_EEJ': ('Spec_W_T2', (0)),
198 'DSPEC_W_T2_EEJ': ('Err_Spec_W_T2', (0)),
199 } #writer
200
201
202 ind = ['range']
203
204 meta = {
205 'kinst': 840, #instrumnet code, 840 for JULIA, 14 for JULIA MP CSR
206 'kindat': 1962, #type of data #Este es el nuevo e igual para JULIA y JULIA CSR
207 'catalog': {
208 'principleInvestigator': 'Danny ScipiΓ³n',
209 'expPurpose': 'Equatorial Electrojet Parameters',
210 },
211 'header': {
212 'analyst': 'D. Hysell'
213 }
214 }
215
216
217 op_writer = procUnitConfObj1.addOperation(name='MADWriter', optype='external')
218 #op_writer.addParameter(name='path', value='/home/roberto/DATA/hdf5_outputs/Madrigal/EEJ')
219 #op_writer.addParameter(name='path', value='/home/roberto/DATA/hdf5_outputs/Madrigal/EEJ/Dop_Max_Val')
220 #op_writer.addParameter(name='path', value='/home/roberto/DATA/hdf5_outputs/Madrigal/EEJ/'+Dop+'/'+dpath[-5:-1])
221 #op_writer.addParameter(name='path', value='/home/roberto/DATA/hdf5_outputs/Madrigal/EEJ/CorrectFiles/no_snl_Test/01/'+Dop+'/'+dpath[-5:-1])
222 op_writer.addParameter(name='path', value='/home/soporte/DATA/hdf5_outputs/Madrigal/EEJ/FinalFiles/'+Dop+'/'+dpath[-5:-1])
223 op_writer.addParameter(name='format', value='hdf5', format='str')
224 op_writer.addParameter(name='oneDDict', value=json.dumps(one), format='str')
225 op_writer.addParameter(name='twoDDict', value=json.dumps(two), format='str')
226 op_writer.addParameter(name='ind2DList', value=json.dumps(ind), format='str')
227 op_writer.addParameter(name='metadata', value=json.dumps(meta), format='str')
228
229
230 if plot_params:
231 '''
232 opObj11 = procUnitConfObj1.addOperation(name='DopplerEEJPlot', optype='other')
233 opObj11.addParameter(name='id', value='10', format='int')
234 opObj11.addParameter(name='wintitle', value='Doppler EEJ', format='str')
235 opObj11.addParameter(name='xmin', value=tiempo[0], format='float')
236 opObj11.addParameter(name='xmax', value=tiempo[1], format='float')
237 opObj11.addParameter(name='ymin', value=altura1[0], format='int')
238 opObj11.addParameter(name='ymax', value=altura1[1], format='int')
239 #opObj11.addParameter(name='zmin', value=-250, format='int')
240 #opObj11.addParameter(name='zmax', value=250, format='int')
241 opObj11.addParameter(name='zlimits', value='(-400,400),(-250,250)')
242 #opObj11.addParameter(name='showprofile', value='1', format='int')
243 opObj11.addParameter(name='save', value=figpath, format='str')
244 #opObj11.addParameter(name='EEJtype', value=1, format='int')
245
246 opObj11 = procUnitConfObj1.addOperation(name='SpcWidthEEJPlot', optype='other')
247 opObj11.addParameter(name='id', value='10', format='int')
248 opObj11.addParameter(name='wintitle', value='Spectral Width EEJ', format='str')
249 opObj11.addParameter(name='xmin', value=tiempo[0], format='float')
250 opObj11.addParameter(name='xmax', value=tiempo[1], format='float')
251 opObj11.addParameter(name='ymin', value=altura1[0], format='int')
252 opObj11.addParameter(name='ymax', value=altura1[1], format='int')
253 #opObj11.addParameter(name='zmin', value=0., format='int')
254 #opObj11.addParameter(name='zmax', value=250, format='int')
255 opObj11.addParameter(name='zlimits', value='(0.1,100),(0.1,250)')
256 #opObj11.addParameter(name='showprofile', value='1', format='int')
257 opObj11.addParameter(name='save', value=figpath, format='str')
258 #opObj11.addParameter(name='EEJtype', value=1, format='int')
259 '''
260 opObj11 = procUnitConfObj1.addOperation(name='SpectraObliquePlot', optype='other')
261 opObj11.addParameter(name='id', value='1', format='int')
262 opObj11.addParameter(name='wintitle', value='Oblique', format='str')
263 opObj11.addParameter(name='zmin', value=db_range[0], format='int')
264 opObj11.addParameter(name='zmax', value=db_range[1], format='int')
265 opObj11.addParameter(name='xaxis', value='velocity', format='str')
266 opObj11.addParameter(name='ymin', value=altura1[0], format='int')
267 opObj11.addParameter(name='ymax', value=altura1[1], format='int')
268 # opObj11.addParameter(name='xmin', value=velocity[0], format='int')
269 # opObj11.addParameter(name='xmax', value=velocity[1], format='int')
270 opObj11.addParameter(name='showprofile', value='1', format='int')
271 opObj11.addParameter(name='save', value=figpath+'/400', format='str')
272
273 '''
274 opObj11 = procUnitConfObj1.addOperation(name='DopplerEEJPlot', optype='other')
275 opObj11.addParameter(name='id', value='10', format='int')
276 opObj11.addParameter(name='wintitle', value='Doppler EEJ Type II', format='str')
277 opObj11.addParameter(name='xmin', value=tiempo[0], format='float')
278 opObj11.addParameter(name='xmax', value=tiempo[1], format='float')
279 # opObj11.addParameter(name='ymin', value=altura1[0], format='int')
280 # opObj11.addParameter(name='ymax', value=altura1[1], format='int')
281 opObj11.addParameter(name='zmin', value=-250, format='int')
282 opObj11.addParameter(name='zmax', value=250, format='int')
283 #opObj11.addParameter(name='showprofile', value='1', format='int')
284 opObj11.addParameter(name='save', value=figpath, format='str')
285 opObj11.addParameter(name='EEJtype', value=2, format='int')
286 '''
287
288
289
290
291 controllerObj.start() No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,937 +1,967
1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 # All rights reserved.
2 # All rights reserved.
3 #
3 #
4 # Distributed under the terms of the BSD 3-clause license.
4 # Distributed under the terms of the BSD 3-clause license.
5 """Spectra processing Unit and operations
5 """Spectra processing Unit and operations
6
6
7 Here you will find the processing unit `SpectraProc` and several operations
7 Here you will find the processing unit `SpectraProc` and several operations
8 to work with Spectra data type
8 to work with Spectra data type
9 """
9 """
10
10
11 import time
11 import time
12 import itertools
12 import itertools
13 import numpy
13 import numpy
14 # repositorio
14 # repositorio
15 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
15 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
16 from schainpy.model.data.jrodata import Spectra
16 from schainpy.model.data.jrodata import Spectra
17 from schainpy.model.data.jrodata import hildebrand_sekhon
17 from schainpy.model.data.jrodata import hildebrand_sekhon
18 from schainpy.utils import log
18 from schainpy.utils import log
19
19
20
20
21 class SpectraProc(ProcessingUnit):
21 class SpectraProc(ProcessingUnit):
22
22
23 def __init__(self):
23 def __init__(self):
24
24
25 ProcessingUnit.__init__(self)
25 ProcessingUnit.__init__(self)
26
26
27 self.buffer = None
27 self.buffer = None
28 self.firstdatatime = None
28 self.firstdatatime = None
29 self.profIndex = 0
29 self.profIndex = 0
30 self.dataOut = Spectra()
30 self.dataOut = Spectra()
31 self.id_min = None
31 self.id_min = None
32 self.id_max = None
32 self.id_max = None
33 self.setupReq = False #Agregar a todas las unidades de proc
33 self.setupReq = False #Agregar a todas las unidades de proc
34
34
35 def __updateSpecFromVoltage(self):
35 def __updateSpecFromVoltage(self):
36
36
37 self.dataOut.timeZone = self.dataIn.timeZone
37 self.dataOut.timeZone = self.dataIn.timeZone
38 self.dataOut.dstFlag = self.dataIn.dstFlag
38 self.dataOut.dstFlag = self.dataIn.dstFlag
39 self.dataOut.errorCount = self.dataIn.errorCount
39 self.dataOut.errorCount = self.dataIn.errorCount
40 self.dataOut.useLocalTime = self.dataIn.useLocalTime
40 self.dataOut.useLocalTime = self.dataIn.useLocalTime
41 try:
41 try:
42 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
42 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
43 except:
43 except:
44 pass
44 pass
45 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
45 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
46 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
46 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
47 self.dataOut.channelList = self.dataIn.channelList
47 self.dataOut.channelList = self.dataIn.channelList
48 self.dataOut.heightList = self.dataIn.heightList
48 self.dataOut.heightList = self.dataIn.heightList
49 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
49 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
50 self.dataOut.nProfiles = self.dataOut.nFFTPoints
50 self.dataOut.nProfiles = self.dataOut.nFFTPoints
51 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
51 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
52 self.dataOut.utctime = self.firstdatatime
52 self.dataOut.utctime = self.firstdatatime
53 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
53 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
54 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
54 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
55 self.dataOut.flagShiftFFT = False
55 self.dataOut.flagShiftFFT = False
56 self.dataOut.nCohInt = self.dataIn.nCohInt
56 self.dataOut.nCohInt = self.dataIn.nCohInt
57 self.dataOut.nIncohInt = 1
57 self.dataOut.nIncohInt = 1
58 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
58 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
59 self.dataOut.frequency = self.dataIn.frequency
59 self.dataOut.frequency = self.dataIn.frequency
60 self.dataOut.realtime = self.dataIn.realtime
60 self.dataOut.realtime = self.dataIn.realtime
61 self.dataOut.azimuth = self.dataIn.azimuth
61 self.dataOut.azimuth = self.dataIn.azimuth
62 self.dataOut.zenith = self.dataIn.zenith
62 self.dataOut.zenith = self.dataIn.zenith
63 self.dataOut.beam.codeList = self.dataIn.beam.codeList
63 self.dataOut.beam.codeList = self.dataIn.beam.codeList
64 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
64 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
65 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
65 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
66 self.dataOut.runNextUnit = self.dataIn.runNextUnit
66 self.dataOut.runNextUnit = self.dataIn.runNextUnit
67 try:
67 try:
68 self.dataOut.step = self.dataIn.step
68 self.dataOut.step = self.dataIn.step
69 except:
69 except:
70 pass
70 pass
71
71
72 def __getFft(self):
72 def __getFft(self):
73 """
73 """
74 Convierte valores de Voltaje a Spectra
74 Convierte valores de Voltaje a Spectra
75
75
76 Affected:
76 Affected:
77 self.dataOut.data_spc
77 self.dataOut.data_spc
78 self.dataOut.data_cspc
78 self.dataOut.data_cspc
79 self.dataOut.data_dc
79 self.dataOut.data_dc
80 self.dataOut.heightList
80 self.dataOut.heightList
81 self.profIndex
81 self.profIndex
82 self.buffer
82 self.buffer
83 self.dataOut.flagNoData
83 self.dataOut.flagNoData
84 """
84 """
85 fft_volt = numpy.fft.fft(
85 fft_volt = numpy.fft.fft(
86 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
86 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
87 fft_volt = fft_volt.astype(numpy.dtype('complex'))
87 fft_volt = fft_volt.astype(numpy.dtype('complex'))
88 dc = fft_volt[:, 0, :]
88 dc = fft_volt[:, 0, :]
89
89
90 # calculo de self-spectra
90 # calculo de self-spectra
91 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
91 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
92 spc = fft_volt * numpy.conjugate(fft_volt)
92 spc = fft_volt * numpy.conjugate(fft_volt)
93 spc = spc.real
93 spc = spc.real
94
94
95 blocksize = 0
95 blocksize = 0
96 blocksize += dc.size
96 blocksize += dc.size
97 blocksize += spc.size
97 blocksize += spc.size
98
98
99 cspc = None
99 cspc = None
100 pairIndex = 0
100 pairIndex = 0
101 if self.dataOut.pairsList != None:
101 if self.dataOut.pairsList != None:
102 # calculo de cross-spectra
102 # calculo de cross-spectra
103 cspc = numpy.zeros(
103 cspc = numpy.zeros(
104 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
104 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
105 for pair in self.dataOut.pairsList:
105 for pair in self.dataOut.pairsList:
106 if pair[0] not in self.dataOut.channelList:
106 if pair[0] not in self.dataOut.channelList:
107 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
107 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
108 str(pair), str(self.dataOut.channelList)))
108 str(pair), str(self.dataOut.channelList)))
109 if pair[1] not in self.dataOut.channelList:
109 if pair[1] not in self.dataOut.channelList:
110 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
110 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
111 str(pair), str(self.dataOut.channelList)))
111 str(pair), str(self.dataOut.channelList)))
112
112
113 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
113 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
114 numpy.conjugate(fft_volt[pair[1], :, :])
114 numpy.conjugate(fft_volt[pair[1], :, :])
115 pairIndex += 1
115 pairIndex += 1
116 blocksize += cspc.size
116 blocksize += cspc.size
117
117
118 self.dataOut.data_spc = spc
118 self.dataOut.data_spc = spc
119 self.dataOut.data_cspc = cspc
119 self.dataOut.data_cspc = cspc
120 self.dataOut.data_dc = dc
120 self.dataOut.data_dc = dc
121 self.dataOut.blockSize = blocksize
121 self.dataOut.blockSize = blocksize
122 self.dataOut.flagShiftFFT = False
122 self.dataOut.flagShiftFFT = False
123
123
124 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False, runNextUnit = 0):
124 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False, runNextUnit = 0):
125
125
126 self.dataIn.runNextUnit = runNextUnit
126 self.dataIn.runNextUnit = runNextUnit
127 if self.dataIn.type == "Spectra":
127 if self.dataIn.type == "Spectra":
128 self.dataOut.copy(self.dataIn)
128 self.dataOut.copy(self.dataIn)
129 if shift_fft:
129 if shift_fft:
130 #desplaza a la derecha en el eje 2 determinadas posiciones
130 #desplaza a la derecha en el eje 2 determinadas posiciones
131 shift = int(self.dataOut.nFFTPoints/2)
131 shift = int(self.dataOut.nFFTPoints/2)
132 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
132 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
133
133
134 if self.dataOut.data_cspc is not None:
134 if self.dataOut.data_cspc is not None:
135 #desplaza a la derecha en el eje 2 determinadas posiciones
135 #desplaza a la derecha en el eje 2 determinadas posiciones
136 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
136 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
137 if pairsList:
137 if pairsList:
138 self.__selectPairs(pairsList)
138 self.__selectPairs(pairsList)
139
139
140 elif self.dataIn.type == "Voltage":
140 elif self.dataIn.type == "Voltage":
141
141
142 self.dataOut.flagNoData = True
142 self.dataOut.flagNoData = True
143
143
144 if nFFTPoints == None:
144 if nFFTPoints == None:
145 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
145 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
146
146
147 if nProfiles == None:
147 if nProfiles == None:
148 nProfiles = nFFTPoints
148 nProfiles = nFFTPoints
149
149
150 if ippFactor == None:
150 if ippFactor == None:
151 self.dataOut.ippFactor = 1
151 self.dataOut.ippFactor = 1
152
152
153 self.dataOut.nFFTPoints = nFFTPoints
153 self.dataOut.nFFTPoints = nFFTPoints
154
154
155 if self.buffer is None:
155 if self.buffer is None:
156 self.buffer = numpy.zeros((self.dataIn.nChannels,
156 self.buffer = numpy.zeros((self.dataIn.nChannels,
157 nProfiles,
157 nProfiles,
158 self.dataIn.nHeights),
158 self.dataIn.nHeights),
159 dtype='complex')
159 dtype='complex')
160
160
161 if self.dataIn.flagDataAsBlock:
161 if self.dataIn.flagDataAsBlock:
162 nVoltProfiles = self.dataIn.data.shape[1]
162 nVoltProfiles = self.dataIn.data.shape[1]
163 if nVoltProfiles == nProfiles:
163 if nVoltProfiles == nProfiles:
164 self.buffer = self.dataIn.data.copy()
164 self.buffer = self.dataIn.data.copy()
165 self.profIndex = nVoltProfiles
165 self.profIndex = nVoltProfiles
166
166
167 elif nVoltProfiles < nProfiles:
167 elif nVoltProfiles < nProfiles:
168
168
169 if self.profIndex == 0:
169 if self.profIndex == 0:
170 self.id_min = 0
170 self.id_min = 0
171 self.id_max = nVoltProfiles
171 self.id_max = nVoltProfiles
172
172
173 self.buffer[:, self.id_min:self.id_max,
173 self.buffer[:, self.id_min:self.id_max,
174 :] = self.dataIn.data
174 :] = self.dataIn.data
175 self.profIndex += nVoltProfiles
175 self.profIndex += nVoltProfiles
176 self.id_min += nVoltProfiles
176 self.id_min += nVoltProfiles
177 self.id_max += nVoltProfiles
177 self.id_max += nVoltProfiles
178 elif nVoltProfiles > nProfiles:
178 elif nVoltProfiles > nProfiles:
179 self.reader.bypass = True
179 self.reader.bypass = True
180 if self.profIndex == 0:
180 if self.profIndex == 0:
181 self.id_min = 0
181 self.id_min = 0
182 self.id_max = nProfiles
182 self.id_max = nProfiles
183
183
184 self.buffer = self.dataIn.data[:, self.id_min:self.id_max,:]
184 self.buffer = self.dataIn.data[:, self.id_min:self.id_max,:]
185 self.profIndex += nProfiles
185 self.profIndex += nProfiles
186 self.id_min += nProfiles
186 self.id_min += nProfiles
187 self.id_max += nProfiles
187 self.id_max += nProfiles
188 if self.id_max == nVoltProfiles:
188 if self.id_max == nVoltProfiles:
189 self.reader.bypass = False
189 self.reader.bypass = False
190
190
191 else:
191 else:
192 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
192 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
193 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
193 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
194 self.dataOut.flagNoData = True
194 self.dataOut.flagNoData = True
195 else:
195 else:
196 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
196 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
197 self.profIndex += 1
197 self.profIndex += 1
198
198
199 if self.firstdatatime == None:
199 if self.firstdatatime == None:
200 self.firstdatatime = self.dataIn.utctime
200 self.firstdatatime = self.dataIn.utctime
201
201
202 if self.profIndex == nProfiles:
202 if self.profIndex == nProfiles:
203 self.__updateSpecFromVoltage()
203 self.__updateSpecFromVoltage()
204 if pairsList == None:
204 if pairsList == None:
205 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
205 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
206 else:
206 else:
207 self.dataOut.pairsList = pairsList
207 self.dataOut.pairsList = pairsList
208 self.__getFft()
208 self.__getFft()
209 self.dataOut.flagNoData = False
209 self.dataOut.flagNoData = False
210 self.firstdatatime = None
210 self.firstdatatime = None
211 #if not self.reader.bypass:
211 #if not self.reader.bypass:
212 self.profIndex = 0
212 self.profIndex = 0
213 else:
213 else:
214 raise ValueError("The type of input object '%s' is not valid".format(
214 raise ValueError("The type of input object '%s' is not valid".format(
215 self.dataIn.type))
215 self.dataIn.type))
216
216
217 def __selectPairs(self, pairsList):
217 def __selectPairs(self, pairsList):
218
218
219 if not pairsList:
219 if not pairsList:
220 return
220 return
221
221
222 pairs = []
222 pairs = []
223 pairsIndex = []
223 pairsIndex = []
224
224
225 for pair in pairsList:
225 for pair in pairsList:
226 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
226 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
227 continue
227 continue
228 pairs.append(pair)
228 pairs.append(pair)
229 pairsIndex.append(pairs.index(pair))
229 pairsIndex.append(pairs.index(pair))
230
230
231 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
231 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
232 self.dataOut.pairsList = pairs
232 self.dataOut.pairsList = pairs
233
233
234 return
234 return
235
235
236 def selectFFTs(self, minFFT, maxFFT ):
236 def selectFFTs(self, minFFT, maxFFT ):
237 """
237 """
238 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
238 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
239 minFFT<= FFT <= maxFFT
239 minFFT<= FFT <= maxFFT
240 """
240 """
241
241
242 if (minFFT > maxFFT):
242 if (minFFT > maxFFT):
243 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
243 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
244
244
245 if (minFFT < self.dataOut.getFreqRange()[0]):
245 if (minFFT < self.dataOut.getFreqRange()[0]):
246 minFFT = self.dataOut.getFreqRange()[0]
246 minFFT = self.dataOut.getFreqRange()[0]
247
247
248 if (maxFFT > self.dataOut.getFreqRange()[-1]):
248 if (maxFFT > self.dataOut.getFreqRange()[-1]):
249 maxFFT = self.dataOut.getFreqRange()[-1]
249 maxFFT = self.dataOut.getFreqRange()[-1]
250
250
251 minIndex = 0
251 minIndex = 0
252 maxIndex = 0
252 maxIndex = 0
253 FFTs = self.dataOut.getFreqRange()
253 FFTs = self.dataOut.getFreqRange()
254
254
255 inda = numpy.where(FFTs >= minFFT)
255 inda = numpy.where(FFTs >= minFFT)
256 indb = numpy.where(FFTs <= maxFFT)
256 indb = numpy.where(FFTs <= maxFFT)
257
257
258 try:
258 try:
259 minIndex = inda[0][0]
259 minIndex = inda[0][0]
260 except:
260 except:
261 minIndex = 0
261 minIndex = 0
262
262
263 try:
263 try:
264 maxIndex = indb[0][-1]
264 maxIndex = indb[0][-1]
265 except:
265 except:
266 maxIndex = len(FFTs)
266 maxIndex = len(FFTs)
267
267
268 self.selectFFTsByIndex(minIndex, maxIndex)
268 self.selectFFTsByIndex(minIndex, maxIndex)
269
269
270 return 1
270 return 1
271
271
272 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
272 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
273 newheis = numpy.where(
273 newheis = numpy.where(
274 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
274 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
275
275
276 if hei_ref != None:
276 if hei_ref != None:
277 newheis = numpy.where(self.dataOut.heightList > hei_ref)
277 newheis = numpy.where(self.dataOut.heightList > hei_ref)
278
278
279 minIndex = min(newheis[0])
279 minIndex = min(newheis[0])
280 maxIndex = max(newheis[0])
280 maxIndex = max(newheis[0])
281 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
281 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
282 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
282 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
283
283
284 # determina indices
284 # determina indices
285 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
285 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
286 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
286 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
287 avg_dB = 10 * \
287 avg_dB = 10 * \
288 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
288 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
289 beacon_dB = numpy.sort(avg_dB)[-nheis:]
289 beacon_dB = numpy.sort(avg_dB)[-nheis:]
290 beacon_heiIndexList = []
290 beacon_heiIndexList = []
291 for val in avg_dB.tolist():
291 for val in avg_dB.tolist():
292 if val >= beacon_dB[0]:
292 if val >= beacon_dB[0]:
293 beacon_heiIndexList.append(avg_dB.tolist().index(val))
293 beacon_heiIndexList.append(avg_dB.tolist().index(val))
294
294
295 data_cspc = None
295 data_cspc = None
296 if self.dataOut.data_cspc is not None:
296 if self.dataOut.data_cspc is not None:
297 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
297 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
298
298
299 data_dc = None
299 data_dc = None
300 if self.dataOut.data_dc is not None:
300 if self.dataOut.data_dc is not None:
301 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
301 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
302
302
303 self.dataOut.data_spc = data_spc
303 self.dataOut.data_spc = data_spc
304 self.dataOut.data_cspc = data_cspc
304 self.dataOut.data_cspc = data_cspc
305 self.dataOut.data_dc = data_dc
305 self.dataOut.data_dc = data_dc
306 self.dataOut.heightList = heightList
306 self.dataOut.heightList = heightList
307 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
307 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
308
308
309 return 1
309 return 1
310
310
311 def selectFFTsByIndex(self, minIndex, maxIndex):
311 def selectFFTsByIndex(self, minIndex, maxIndex):
312 """
312 """
313
313
314 """
314 """
315
315
316 if (minIndex < 0) or (minIndex > maxIndex):
316 if (minIndex < 0) or (minIndex > maxIndex):
317 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
317 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
318
318
319 if (maxIndex >= self.dataOut.nProfiles):
319 if (maxIndex >= self.dataOut.nProfiles):
320 maxIndex = self.dataOut.nProfiles-1
320 maxIndex = self.dataOut.nProfiles-1
321
321
322 #Spectra
322 #Spectra
323 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
323 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
324
324
325 data_cspc = None
325 data_cspc = None
326 if self.dataOut.data_cspc is not None:
326 if self.dataOut.data_cspc is not None:
327 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
327 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
328
328
329 data_dc = None
329 data_dc = None
330 if self.dataOut.data_dc is not None:
330 if self.dataOut.data_dc is not None:
331 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
331 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
332
332
333 self.dataOut.data_spc = data_spc
333 self.dataOut.data_spc = data_spc
334 self.dataOut.data_cspc = data_cspc
334 self.dataOut.data_cspc = data_cspc
335 self.dataOut.data_dc = data_dc
335 self.dataOut.data_dc = data_dc
336
336
337 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
337 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
338 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
338 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
339 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
339 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
340
340
341 return 1
341 return 1
342
342
343 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
343 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
344 # validacion de rango
344 # validacion de rango
345 if minHei == None:
345 if minHei == None:
346 minHei = self.dataOut.heightList[0]
346 minHei = self.dataOut.heightList[0]
347
347
348 if maxHei == None:
348 if maxHei == None:
349 maxHei = self.dataOut.heightList[-1]
349 maxHei = self.dataOut.heightList[-1]
350
350
351 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
351 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
352 print('minHei: %.2f is out of the heights range' % (minHei))
352 print('minHei: %.2f is out of the heights range' % (minHei))
353 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
353 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
354 minHei = self.dataOut.heightList[0]
354 minHei = self.dataOut.heightList[0]
355
355
356 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
356 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
357 print('maxHei: %.2f is out of the heights range' % (maxHei))
357 print('maxHei: %.2f is out of the heights range' % (maxHei))
358 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
358 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
359 maxHei = self.dataOut.heightList[-1]
359 maxHei = self.dataOut.heightList[-1]
360
360
361 # validacion de velocidades
361 # validacion de velocidades
362 velrange = self.dataOut.getVelRange(1)
362 velrange = self.dataOut.getVelRange(1)
363
363
364 if minVel == None:
364 if minVel == None:
365 minVel = velrange[0]
365 minVel = velrange[0]
366
366
367 if maxVel == None:
367 if maxVel == None:
368 maxVel = velrange[-1]
368 maxVel = velrange[-1]
369
369
370 if (minVel < velrange[0]) or (minVel > maxVel):
370 if (minVel < velrange[0]) or (minVel > maxVel):
371 print('minVel: %.2f is out of the velocity range' % (minVel))
371 print('minVel: %.2f is out of the velocity range' % (minVel))
372 print('minVel is setting to %.2f' % (velrange[0]))
372 print('minVel is setting to %.2f' % (velrange[0]))
373 minVel = velrange[0]
373 minVel = velrange[0]
374
374
375 if (maxVel > velrange[-1]) or (maxVel < minVel):
375 if (maxVel > velrange[-1]) or (maxVel < minVel):
376 print('maxVel: %.2f is out of the velocity range' % (maxVel))
376 print('maxVel: %.2f is out of the velocity range' % (maxVel))
377 print('maxVel is setting to %.2f' % (velrange[-1]))
377 print('maxVel is setting to %.2f' % (velrange[-1]))
378 maxVel = velrange[-1]
378 maxVel = velrange[-1]
379
379
380 # seleccion de indices para rango
380 # seleccion de indices para rango
381 minIndex = 0
381 minIndex = 0
382 maxIndex = 0
382 maxIndex = 0
383 heights = self.dataOut.heightList
383 heights = self.dataOut.heightList
384
384
385 inda = numpy.where(heights >= minHei)
385 inda = numpy.where(heights >= minHei)
386 indb = numpy.where(heights <= maxHei)
386 indb = numpy.where(heights <= maxHei)
387
387
388 try:
388 try:
389 minIndex = inda[0][0]
389 minIndex = inda[0][0]
390 except:
390 except:
391 minIndex = 0
391 minIndex = 0
392
392
393 try:
393 try:
394 maxIndex = indb[0][-1]
394 maxIndex = indb[0][-1]
395 except:
395 except:
396 maxIndex = len(heights)
396 maxIndex = len(heights)
397
397
398 if (minIndex < 0) or (minIndex > maxIndex):
398 if (minIndex < 0) or (minIndex > maxIndex):
399 raise ValueError("some value in (%d,%d) is not valid" % (
399 raise ValueError("some value in (%d,%d) is not valid" % (
400 minIndex, maxIndex))
400 minIndex, maxIndex))
401
401
402 if (maxIndex >= self.dataOut.nHeights):
402 if (maxIndex >= self.dataOut.nHeights):
403 maxIndex = self.dataOut.nHeights - 1
403 maxIndex = self.dataOut.nHeights - 1
404
404
405 # seleccion de indices para velocidades
405 # seleccion de indices para velocidades
406 indminvel = numpy.where(velrange >= minVel)
406 indminvel = numpy.where(velrange >= minVel)
407 indmaxvel = numpy.where(velrange <= maxVel)
407 indmaxvel = numpy.where(velrange <= maxVel)
408 try:
408 try:
409 minIndexVel = indminvel[0][0]
409 minIndexVel = indminvel[0][0]
410 except:
410 except:
411 minIndexVel = 0
411 minIndexVel = 0
412
412
413 try:
413 try:
414 maxIndexVel = indmaxvel[0][-1]
414 maxIndexVel = indmaxvel[0][-1]
415 except:
415 except:
416 maxIndexVel = len(velrange)
416 maxIndexVel = len(velrange)
417
417
418 # seleccion del espectro
418 # seleccion del espectro
419 data_spc = self.dataOut.data_spc[:,
419 data_spc = self.dataOut.data_spc[:,
420 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
420 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
421 # estimacion de ruido
421 # estimacion de ruido
422 noise = numpy.zeros(self.dataOut.nChannels)
422 noise = numpy.zeros(self.dataOut.nChannels)
423
423
424 for channel in range(self.dataOut.nChannels):
424 for channel in range(self.dataOut.nChannels):
425 daux = data_spc[channel, :, :]
425 daux = data_spc[channel, :, :]
426 sortdata = numpy.sort(daux, axis=None)
426 sortdata = numpy.sort(daux, axis=None)
427 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
427 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
428
428
429 self.dataOut.noise_estimation = noise.copy()
429 self.dataOut.noise_estimation = noise.copy()
430
430
431 return 1
431 return 1
432
432
433
434 class GetSNR(Operation):
435 '''
436 Written by R. Flores
437 '''
438 """Operation to get SNR.
439
440 Parameters:
441 -----------
442
443 Example
444 --------
445
446 op = proc_unit.addOperation(name='GetSNR', optype='other')
447
448 """
449
450 def __init__(self, **kwargs):
451
452 Operation.__init__(self, **kwargs)
453
454 def run(self,dataOut):
455
456 noise = dataOut.getNoise(ymin_index=-10) #RegiΓ³n superior donde solo deberΓ­a de haber ruido
457 dataOut.data_snr = (dataOut.data_spc.sum(axis=1)-noise[:,None]*dataOut.nFFTPoints)/(noise[:,None]*dataOut.nFFTPoints) #It works apparently
458 dataOut.snl = numpy.log10(dataOut.data_snr)
459 dataOut.snl = numpy.where(dataOut.data_snr<.01, numpy.nan, dataOut.snl)
460
461 return dataOut
462
433 class removeDC(Operation):
463 class removeDC(Operation):
434
464
435 def run(self, dataOut, mode=2):
465 def run(self, dataOut, mode=2):
436 self.dataOut = dataOut
466 self.dataOut = dataOut
437 jspectra = self.dataOut.data_spc
467 jspectra = self.dataOut.data_spc
438 jcspectra = self.dataOut.data_cspc
468 jcspectra = self.dataOut.data_cspc
439
469
440 num_chan = jspectra.shape[0]
470 num_chan = jspectra.shape[0]
441 num_hei = jspectra.shape[2]
471 num_hei = jspectra.shape[2]
442
472
443 if jcspectra is not None:
473 if jcspectra is not None:
444 jcspectraExist = True
474 jcspectraExist = True
445 num_pairs = jcspectra.shape[0]
475 num_pairs = jcspectra.shape[0]
446 else:
476 else:
447 jcspectraExist = False
477 jcspectraExist = False
448
478
449 freq_dc = int(jspectra.shape[1] / 2)
479 freq_dc = int(jspectra.shape[1] / 2)
450 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
480 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
451 ind_vel = ind_vel.astype(int)
481 ind_vel = ind_vel.astype(int)
452
482
453 if ind_vel[0] < 0:
483 if ind_vel[0] < 0:
454 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
484 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
455
485
456 if mode == 1:
486 if mode == 1:
457 jspectra[:, freq_dc, :] = (
487 jspectra[:, freq_dc, :] = (
458 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
488 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
459
489
460 if jcspectraExist:
490 if jcspectraExist:
461 jcspectra[:, freq_dc, :] = (
491 jcspectra[:, freq_dc, :] = (
462 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
492 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
463
493
464 if mode == 2:
494 if mode == 2:
465
495
466 vel = numpy.array([-2, -1, 1, 2])
496 vel = numpy.array([-2, -1, 1, 2])
467 xx = numpy.zeros([4, 4])
497 xx = numpy.zeros([4, 4])
468
498
469 for fil in range(4):
499 for fil in range(4):
470 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
500 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
471
501
472 xx_inv = numpy.linalg.inv(xx)
502 xx_inv = numpy.linalg.inv(xx)
473 xx_aux = xx_inv[0, :]
503 xx_aux = xx_inv[0, :]
474
504
475 for ich in range(num_chan):
505 for ich in range(num_chan):
476 yy = jspectra[ich, ind_vel, :]
506 yy = jspectra[ich, ind_vel, :]
477 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
507 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
478
508
479 junkid = jspectra[ich, freq_dc, :] <= 0
509 junkid = jspectra[ich, freq_dc, :] <= 0
480 cjunkid = sum(junkid)
510 cjunkid = sum(junkid)
481
511
482 if cjunkid.any():
512 if cjunkid.any():
483 jspectra[ich, freq_dc, junkid.nonzero()] = (
513 jspectra[ich, freq_dc, junkid.nonzero()] = (
484 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
514 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
485
515
486 if jcspectraExist:
516 if jcspectraExist:
487 for ip in range(num_pairs):
517 for ip in range(num_pairs):
488 yy = jcspectra[ip, ind_vel, :]
518 yy = jcspectra[ip, ind_vel, :]
489 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
519 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
490
520
491 self.dataOut.data_spc = jspectra
521 self.dataOut.data_spc = jspectra
492 self.dataOut.data_cspc = jcspectra
522 self.dataOut.data_cspc = jcspectra
493
523
494 return self.dataOut
524 return self.dataOut
495
525
496 class removeInterference(Operation):
526 class removeInterference(Operation):
497
527
498 def removeInterference2(self):
528 def removeInterference2(self):
499
529
500 cspc = self.dataOut.data_cspc
530 cspc = self.dataOut.data_cspc
501 spc = self.dataOut.data_spc
531 spc = self.dataOut.data_spc
502 Heights = numpy.arange(cspc.shape[2])
532 Heights = numpy.arange(cspc.shape[2])
503 realCspc = numpy.abs(cspc)
533 realCspc = numpy.abs(cspc)
504
534
505 for i in range(cspc.shape[0]):
535 for i in range(cspc.shape[0]):
506 LinePower= numpy.sum(realCspc[i], axis=0)
536 LinePower= numpy.sum(realCspc[i], axis=0)
507 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
537 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
508 SelectedHeights = Heights[ numpy.where(LinePower < Threshold) ]
538 SelectedHeights = Heights[ numpy.where(LinePower < Threshold) ]
509 InterferenceSum = numpy.sum(realCspc[i,:,SelectedHeights],axis=0)
539 InterferenceSum = numpy.sum(realCspc[i,:,SelectedHeights],axis=0)
510 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
540 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
511 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
541 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
512
542
513
543
514 InterferenceRange = numpy.where(([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
544 InterferenceRange = numpy.where(([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
515 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
545 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
516 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
546 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
517 cspc[i,InterferenceRange,:] = numpy.NaN
547 cspc[i,InterferenceRange,:] = numpy.NaN
518
548
519 self.dataOut.data_cspc = cspc
549 self.dataOut.data_cspc = cspc
520
550
521 def removeInterference(self, interf=2, hei_interf=None, nhei_interf=None, offhei_interf=None):
551 def removeInterference(self, interf=2, hei_interf=None, nhei_interf=None, offhei_interf=None):
522
552
523 jspectra = self.dataOut.data_spc
553 jspectra = self.dataOut.data_spc
524 jcspectra = self.dataOut.data_cspc
554 jcspectra = self.dataOut.data_cspc
525 jnoise = self.dataOut.getNoise()
555 jnoise = self.dataOut.getNoise()
526 num_incoh = self.dataOut.nIncohInt
556 num_incoh = self.dataOut.nIncohInt
527
557
528 num_channel = jspectra.shape[0]
558 num_channel = jspectra.shape[0]
529 num_prof = jspectra.shape[1]
559 num_prof = jspectra.shape[1]
530 num_hei = jspectra.shape[2]
560 num_hei = jspectra.shape[2]
531
561
532 # hei_interf
562 # hei_interf
533 if hei_interf is None:
563 if hei_interf is None:
534 count_hei = int(num_hei / 2)
564 count_hei = int(num_hei / 2)
535 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
565 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
536 hei_interf = numpy.asarray(hei_interf)[0]
566 hei_interf = numpy.asarray(hei_interf)[0]
537 # nhei_interf
567 # nhei_interf
538 if (nhei_interf == None):
568 if (nhei_interf == None):
539 nhei_interf = 5
569 nhei_interf = 5
540 if (nhei_interf < 1):
570 if (nhei_interf < 1):
541 nhei_interf = 1
571 nhei_interf = 1
542 if (nhei_interf > count_hei):
572 if (nhei_interf > count_hei):
543 nhei_interf = count_hei
573 nhei_interf = count_hei
544 if (offhei_interf == None):
574 if (offhei_interf == None):
545 offhei_interf = 0
575 offhei_interf = 0
546
576
547 ind_hei = list(range(num_hei))
577 ind_hei = list(range(num_hei))
548 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
578 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
549 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
579 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
550 mask_prof = numpy.asarray(list(range(num_prof)))
580 mask_prof = numpy.asarray(list(range(num_prof)))
551 num_mask_prof = mask_prof.size
581 num_mask_prof = mask_prof.size
552 comp_mask_prof = [0, num_prof / 2]
582 comp_mask_prof = [0, num_prof / 2]
553
583
554 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
584 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
555 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
585 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
556 jnoise = numpy.nan
586 jnoise = numpy.nan
557 noise_exist = jnoise[0] < numpy.Inf
587 noise_exist = jnoise[0] < numpy.Inf
558
588
559 # Subrutina de Remocion de la Interferencia
589 # Subrutina de Remocion de la Interferencia
560 for ich in range(num_channel):
590 for ich in range(num_channel):
561 # Se ordena los espectros segun su potencia (menor a mayor)
591 # Se ordena los espectros segun su potencia (menor a mayor)
562 power = jspectra[ich, mask_prof, :]
592 power = jspectra[ich, mask_prof, :]
563 power = power[:, hei_interf]
593 power = power[:, hei_interf]
564 power = power.sum(axis=0)
594 power = power.sum(axis=0)
565 psort = power.ravel().argsort()
595 psort = power.ravel().argsort()
566
596
567 # Se estima la interferencia promedio en los Espectros de Potencia empleando
597 # Se estima la interferencia promedio en los Espectros de Potencia empleando
568 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
598 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
569 offhei_interf, nhei_interf + offhei_interf))]]]
599 offhei_interf, nhei_interf + offhei_interf))]]]
570
600
571 if noise_exist:
601 if noise_exist:
572 # tmp_noise = jnoise[ich] / num_prof
602 # tmp_noise = jnoise[ich] / num_prof
573 tmp_noise = jnoise[ich]
603 tmp_noise = jnoise[ich]
574 junkspc_interf = junkspc_interf - tmp_noise
604 junkspc_interf = junkspc_interf - tmp_noise
575 #junkspc_interf[:,comp_mask_prof] = 0
605 #junkspc_interf[:,comp_mask_prof] = 0
576
606
577 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
607 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
578 jspc_interf = jspc_interf.transpose()
608 jspc_interf = jspc_interf.transpose()
579 # Calculando el espectro de interferencia promedio
609 # Calculando el espectro de interferencia promedio
580 noiseid = numpy.where(
610 noiseid = numpy.where(
581 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
611 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
582 noiseid = noiseid[0]
612 noiseid = noiseid[0]
583 cnoiseid = noiseid.size
613 cnoiseid = noiseid.size
584 interfid = numpy.where(
614 interfid = numpy.where(
585 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
615 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
586 interfid = interfid[0]
616 interfid = interfid[0]
587 cinterfid = interfid.size
617 cinterfid = interfid.size
588
618
589 if (cnoiseid > 0):
619 if (cnoiseid > 0):
590 jspc_interf[noiseid] = 0
620 jspc_interf[noiseid] = 0
591
621
592 # Expandiendo los perfiles a limpiar
622 # Expandiendo los perfiles a limpiar
593 if (cinterfid > 0):
623 if (cinterfid > 0):
594 new_interfid = (
624 new_interfid = (
595 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
625 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
596 new_interfid = numpy.asarray(new_interfid)
626 new_interfid = numpy.asarray(new_interfid)
597 new_interfid = {x for x in new_interfid}
627 new_interfid = {x for x in new_interfid}
598 new_interfid = numpy.array(list(new_interfid))
628 new_interfid = numpy.array(list(new_interfid))
599 new_cinterfid = new_interfid.size
629 new_cinterfid = new_interfid.size
600 else:
630 else:
601 new_cinterfid = 0
631 new_cinterfid = 0
602
632
603 for ip in range(new_cinterfid):
633 for ip in range(new_cinterfid):
604 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
634 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
605 jspc_interf[new_interfid[ip]
635 jspc_interf[new_interfid[ip]
606 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
636 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
607
637
608 jspectra[ich, :, ind_hei] = jspectra[ich, :,
638 jspectra[ich, :, ind_hei] = jspectra[ich, :,
609 ind_hei] - jspc_interf # Corregir indices
639 ind_hei] - jspc_interf # Corregir indices
610
640
611 # Removiendo la interferencia del punto de mayor interferencia
641 # Removiendo la interferencia del punto de mayor interferencia
612 ListAux = jspc_interf[mask_prof].tolist()
642 ListAux = jspc_interf[mask_prof].tolist()
613 maxid = ListAux.index(max(ListAux))
643 maxid = ListAux.index(max(ListAux))
614
644
615 if cinterfid > 0:
645 if cinterfid > 0:
616 for ip in range(cinterfid * (interf == 2) - 1):
646 for ip in range(cinterfid * (interf == 2) - 1):
617 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
647 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
618 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
648 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
619 cind = len(ind)
649 cind = len(ind)
620
650
621 if (cind > 0):
651 if (cind > 0):
622 jspectra[ich, interfid[ip], ind] = tmp_noise * \
652 jspectra[ich, interfid[ip], ind] = tmp_noise * \
623 (1 + (numpy.random.uniform(cind) - 0.5) /
653 (1 + (numpy.random.uniform(cind) - 0.5) /
624 numpy.sqrt(num_incoh))
654 numpy.sqrt(num_incoh))
625
655
626 ind = numpy.array([-2, -1, 1, 2])
656 ind = numpy.array([-2, -1, 1, 2])
627 xx = numpy.zeros([4, 4])
657 xx = numpy.zeros([4, 4])
628
658
629 for id1 in range(4):
659 for id1 in range(4):
630 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
660 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
631
661
632 xx_inv = numpy.linalg.inv(xx)
662 xx_inv = numpy.linalg.inv(xx)
633 xx = xx_inv[:, 0]
663 xx = xx_inv[:, 0]
634 ind = (ind + maxid + num_mask_prof) % num_mask_prof
664 ind = (ind + maxid + num_mask_prof) % num_mask_prof
635 yy = jspectra[ich, mask_prof[ind], :]
665 yy = jspectra[ich, mask_prof[ind], :]
636 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
666 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
637 yy.transpose(), xx)
667 yy.transpose(), xx)
638
668
639 indAux = (jspectra[ich, :, :] < tmp_noise *
669 indAux = (jspectra[ich, :, :] < tmp_noise *
640 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
670 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
641 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
671 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
642 (1 - 1 / numpy.sqrt(num_incoh))
672 (1 - 1 / numpy.sqrt(num_incoh))
643
673
644 # Remocion de Interferencia en el Cross Spectra
674 # Remocion de Interferencia en el Cross Spectra
645 if jcspectra is None:
675 if jcspectra is None:
646 return jspectra, jcspectra
676 return jspectra, jcspectra
647 num_pairs = int(jcspectra.size / (num_prof * num_hei))
677 num_pairs = int(jcspectra.size / (num_prof * num_hei))
648 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
678 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
649
679
650 for ip in range(num_pairs):
680 for ip in range(num_pairs):
651
681
652 #-------------------------------------------
682 #-------------------------------------------
653
683
654 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
684 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
655 cspower = cspower[:, hei_interf]
685 cspower = cspower[:, hei_interf]
656 cspower = cspower.sum(axis=0)
686 cspower = cspower.sum(axis=0)
657
687
658 cspsort = cspower.ravel().argsort()
688 cspsort = cspower.ravel().argsort()
659 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
689 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
660 offhei_interf, nhei_interf + offhei_interf))]]]
690 offhei_interf, nhei_interf + offhei_interf))]]]
661 junkcspc_interf = junkcspc_interf.transpose()
691 junkcspc_interf = junkcspc_interf.transpose()
662 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
692 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
663
693
664 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
694 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
665
695
666 median_real = int(numpy.median(numpy.real(
696 median_real = int(numpy.median(numpy.real(
667 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
697 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
668 median_imag = int(numpy.median(numpy.imag(
698 median_imag = int(numpy.median(numpy.imag(
669 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
699 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
670 comp_mask_prof = [int(e) for e in comp_mask_prof]
700 comp_mask_prof = [int(e) for e in comp_mask_prof]
671 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
701 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
672 median_real, median_imag)
702 median_real, median_imag)
673
703
674 for iprof in range(num_prof):
704 for iprof in range(num_prof):
675 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
705 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
676 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
706 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
677
707
678 # Removiendo la Interferencia
708 # Removiendo la Interferencia
679 jcspectra[ip, :, ind_hei] = jcspectra[ip,
709 jcspectra[ip, :, ind_hei] = jcspectra[ip,
680 :, ind_hei] - jcspc_interf
710 :, ind_hei] - jcspc_interf
681
711
682 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
712 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
683 maxid = ListAux.index(max(ListAux))
713 maxid = ListAux.index(max(ListAux))
684
714
685 ind = numpy.array([-2, -1, 1, 2])
715 ind = numpy.array([-2, -1, 1, 2])
686 xx = numpy.zeros([4, 4])
716 xx = numpy.zeros([4, 4])
687
717
688 for id1 in range(4):
718 for id1 in range(4):
689 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
719 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
690
720
691 xx_inv = numpy.linalg.inv(xx)
721 xx_inv = numpy.linalg.inv(xx)
692 xx = xx_inv[:, 0]
722 xx = xx_inv[:, 0]
693
723
694 ind = (ind + maxid + num_mask_prof) % num_mask_prof
724 ind = (ind + maxid + num_mask_prof) % num_mask_prof
695 yy = jcspectra[ip, mask_prof[ind], :]
725 yy = jcspectra[ip, mask_prof[ind], :]
696 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
726 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
697
727
698 # Guardar Resultados
728 # Guardar Resultados
699 self.dataOut.data_spc = jspectra
729 self.dataOut.data_spc = jspectra
700 self.dataOut.data_cspc = jcspectra
730 self.dataOut.data_cspc = jcspectra
701
731
702 return 1
732 return 1
703
733
704 def run(self, dataOut, interf=2,hei_interf=None, nhei_interf=None, offhei_interf=None, mode=1):
734 def run(self, dataOut, interf=2,hei_interf=None, nhei_interf=None, offhei_interf=None, mode=1):
705
735
706 self.dataOut = dataOut
736 self.dataOut = dataOut
707
737
708 if mode == 1:
738 if mode == 1:
709 self.removeInterference(interf=2,hei_interf=None, nhei_interf=None, offhei_interf=None)
739 self.removeInterference(interf=2,hei_interf=None, nhei_interf=None, offhei_interf=None)
710 elif mode == 2:
740 elif mode == 2:
711 self.removeInterference2()
741 self.removeInterference2()
712
742
713 return self.dataOut
743 return self.dataOut
714
744
715
745
716 class deflip(Operation):
746 class deflip(Operation):
717
747
718 def run(self, dataOut):
748 def run(self, dataOut):
719 # arreglo 1: (num_chan, num_profiles, num_heights)
749 # arreglo 1: (num_chan, num_profiles, num_heights)
720 self.dataOut = dataOut
750 self.dataOut = dataOut
721
751
722 # JULIA-oblicua, indice 2
752 # JULIA-oblicua, indice 2
723 # arreglo 2: (num_profiles, num_heights)
753 # arreglo 2: (num_profiles, num_heights)
724 jspectra = self.dataOut.data_spc[2]
754 jspectra = self.dataOut.data_spc[2]
725 jspectra_tmp=numpy.zeros(jspectra.shape)
755 jspectra_tmp=numpy.zeros(jspectra.shape)
726 num_profiles=jspectra.shape[0]
756 num_profiles=jspectra.shape[0]
727 freq_dc = int(num_profiles / 2)
757 freq_dc = int(num_profiles / 2)
728 # Flip con for
758 # Flip con for
729 for j in range(num_profiles):
759 for j in range(num_profiles):
730 jspectra_tmp[num_profiles-j-1]= jspectra[j]
760 jspectra_tmp[num_profiles-j-1]= jspectra[j]
731 # Intercambio perfil de DC con perfil inmediato anterior
761 # Intercambio perfil de DC con perfil inmediato anterior
732 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
762 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
733 jspectra_tmp[freq_dc]= jspectra[freq_dc]
763 jspectra_tmp[freq_dc]= jspectra[freq_dc]
734 # canal modificado es re-escrito en el arreglo de canales
764 # canal modificado es re-escrito en el arreglo de canales
735 self.dataOut.data_spc[2] = jspectra_tmp
765 self.dataOut.data_spc[2] = jspectra_tmp
736
766
737 return self.dataOut
767 return self.dataOut
738
768
739
769
740 class IncohInt(Operation):
770 class IncohInt(Operation):
741
771
742 __profIndex = 0
772 __profIndex = 0
743 __withOverapping = False
773 __withOverapping = False
744
774
745 __byTime = False
775 __byTime = False
746 __initime = None
776 __initime = None
747 __lastdatatime = None
777 __lastdatatime = None
748 __integrationtime = None
778 __integrationtime = None
749
779
750 __buffer_spc = None
780 __buffer_spc = None
751 __buffer_cspc = None
781 __buffer_cspc = None
752 __buffer_dc = None
782 __buffer_dc = None
753
783
754 __dataReady = False
784 __dataReady = False
755
785
756 __timeInterval = None
786 __timeInterval = None
757
787
758 n = None
788 n = None
759
789
760 def __init__(self):
790 def __init__(self):
761
791
762 Operation.__init__(self)
792 Operation.__init__(self)
763
793
764 def setup(self, n=None, timeInterval=None, overlapping=False):
794 def setup(self, n=None, timeInterval=None, overlapping=False):
765 """
795 """
766 Set the parameters of the integration class.
796 Set the parameters of the integration class.
767
797
768 Inputs:
798 Inputs:
769
799
770 n : Number of coherent integrations
800 n : Number of coherent integrations
771 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
801 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
772 overlapping :
802 overlapping :
773
803
774 """
804 """
775
805
776 self.__initime = None
806 self.__initime = None
777 self.__lastdatatime = 0
807 self.__lastdatatime = 0
778
808
779 self.__buffer_spc = 0
809 self.__buffer_spc = 0
780 self.__buffer_cspc = 0
810 self.__buffer_cspc = 0
781 self.__buffer_dc = 0
811 self.__buffer_dc = 0
782
812
783 self.__profIndex = 0
813 self.__profIndex = 0
784 self.__dataReady = False
814 self.__dataReady = False
785 self.__byTime = False
815 self.__byTime = False
786
816
787 if n is None and timeInterval is None:
817 if n is None and timeInterval is None:
788 raise ValueError("n or timeInterval should be specified ...")
818 raise ValueError("n or timeInterval should be specified ...")
789
819
790 if n is not None:
820 if n is not None:
791 self.n = int(n)
821 self.n = int(n)
792 else:
822 else:
793
823
794 self.__integrationtime = int(timeInterval)
824 self.__integrationtime = int(timeInterval)
795 self.n = None
825 self.n = None
796 self.__byTime = True
826 self.__byTime = True
797
827
798 def putData(self, data_spc, data_cspc, data_dc):
828 def putData(self, data_spc, data_cspc, data_dc):
799 """
829 """
800 Add a profile to the __buffer_spc and increase in one the __profileIndex
830 Add a profile to the __buffer_spc and increase in one the __profileIndex
801
831
802 """
832 """
803
833
804 self.__buffer_spc += data_spc
834 self.__buffer_spc += data_spc
805
835
806 if data_cspc is None:
836 if data_cspc is None:
807 self.__buffer_cspc = None
837 self.__buffer_cspc = None
808 else:
838 else:
809 self.__buffer_cspc += data_cspc
839 self.__buffer_cspc += data_cspc
810
840
811 if data_dc is None:
841 if data_dc is None:
812 self.__buffer_dc = None
842 self.__buffer_dc = None
813 else:
843 else:
814 self.__buffer_dc += data_dc
844 self.__buffer_dc += data_dc
815
845
816 self.__profIndex += 1
846 self.__profIndex += 1
817
847
818 return
848 return
819
849
820 def pushData(self):
850 def pushData(self):
821 """
851 """
822 Return the sum of the last profiles and the profiles used in the sum.
852 Return the sum of the last profiles and the profiles used in the sum.
823
853
824 Affected:
854 Affected:
825
855
826 self.__profileIndex
856 self.__profileIndex
827
857
828 """
858 """
829
859
830 data_spc = self.__buffer_spc
860 data_spc = self.__buffer_spc
831 data_cspc = self.__buffer_cspc
861 data_cspc = self.__buffer_cspc
832 data_dc = self.__buffer_dc
862 data_dc = self.__buffer_dc
833 n = self.__profIndex
863 n = self.__profIndex
834
864
835 self.__buffer_spc = 0
865 self.__buffer_spc = 0
836 self.__buffer_cspc = 0
866 self.__buffer_cspc = 0
837 self.__buffer_dc = 0
867 self.__buffer_dc = 0
838 self.__profIndex = 0
868 self.__profIndex = 0
839
869
840 return data_spc, data_cspc, data_dc, n
870 return data_spc, data_cspc, data_dc, n
841
871
842 def byProfiles(self, *args):
872 def byProfiles(self, *args):
843
873
844 self.__dataReady = False
874 self.__dataReady = False
845 avgdata_spc = None
875 avgdata_spc = None
846 avgdata_cspc = None
876 avgdata_cspc = None
847 avgdata_dc = None
877 avgdata_dc = None
848
878
849 self.putData(*args)
879 self.putData(*args)
850
880
851 if self.__profIndex == self.n:
881 if self.__profIndex == self.n:
852
882
853 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
883 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
854 self.n = n
884 self.n = n
855 self.__dataReady = True
885 self.__dataReady = True
856
886
857 return avgdata_spc, avgdata_cspc, avgdata_dc
887 return avgdata_spc, avgdata_cspc, avgdata_dc
858
888
859 def byTime(self, datatime, *args):
889 def byTime(self, datatime, *args):
860
890
861 self.__dataReady = False
891 self.__dataReady = False
862 avgdata_spc = None
892 avgdata_spc = None
863 avgdata_cspc = None
893 avgdata_cspc = None
864 avgdata_dc = None
894 avgdata_dc = None
865
895
866 self.putData(*args)
896 self.putData(*args)
867
897
868 if (datatime - self.__initime) >= self.__integrationtime:
898 if (datatime - self.__initime) >= self.__integrationtime:
869 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
899 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
870 self.n = n
900 self.n = n
871 self.__dataReady = True
901 self.__dataReady = True
872
902
873 return avgdata_spc, avgdata_cspc, avgdata_dc
903 return avgdata_spc, avgdata_cspc, avgdata_dc
874
904
875 def integrate(self, datatime, *args):
905 def integrate(self, datatime, *args):
876
906
877 if self.__profIndex == 0:
907 if self.__profIndex == 0:
878 self.__initime = datatime
908 self.__initime = datatime
879
909
880 if self.__byTime:
910 if self.__byTime:
881 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
911 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
882 datatime, *args)
912 datatime, *args)
883 else:
913 else:
884 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
914 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
885
915
886 if not self.__dataReady:
916 if not self.__dataReady:
887 return None, None, None, None
917 return None, None, None, None
888
918
889 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
919 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
890
920
891 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
921 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
892 if n == 1:
922 if n == 1:
893 return dataOut
923 return dataOut
894
924
895 dataOut.flagNoData = True
925 dataOut.flagNoData = True
896
926
897 if not self.isConfig:
927 if not self.isConfig:
898 self.setup(n, timeInterval, overlapping)
928 self.setup(n, timeInterval, overlapping)
899 self.isConfig = True
929 self.isConfig = True
900
930
901 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
931 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
902 dataOut.data_spc,
932 dataOut.data_spc,
903 dataOut.data_cspc,
933 dataOut.data_cspc,
904 dataOut.data_dc)
934 dataOut.data_dc)
905
935
906 if self.__dataReady:
936 if self.__dataReady:
907
937
908 dataOut.data_spc = avgdata_spc
938 dataOut.data_spc = avgdata_spc
909 dataOut.data_cspc = avgdata_cspc
939 dataOut.data_cspc = avgdata_cspc
910 dataOut.data_dc = avgdata_dc
940 dataOut.data_dc = avgdata_dc
911 dataOut.nIncohInt *= self.n
941 dataOut.nIncohInt *= self.n
912 dataOut.utctime = avgdatatime
942 dataOut.utctime = avgdatatime
913 dataOut.flagNoData = False
943 dataOut.flagNoData = False
914
944
915 return dataOut
945 return dataOut
916
946
917 class dopplerFlip(Operation):
947 class dopplerFlip(Operation):
918
948
919 def run(self, dataOut):
949 def run(self, dataOut):
920 # arreglo 1: (num_chan, num_profiles, num_heights)
950 # arreglo 1: (num_chan, num_profiles, num_heights)
921 self.dataOut = dataOut
951 self.dataOut = dataOut
922 # JULIA-oblicua, indice 2
952 # JULIA-oblicua, indice 2
923 # arreglo 2: (num_profiles, num_heights)
953 # arreglo 2: (num_profiles, num_heights)
924 jspectra = self.dataOut.data_spc[2]
954 jspectra = self.dataOut.data_spc[2]
925 jspectra_tmp = numpy.zeros(jspectra.shape)
955 jspectra_tmp = numpy.zeros(jspectra.shape)
926 num_profiles = jspectra.shape[0]
956 num_profiles = jspectra.shape[0]
927 freq_dc = int(num_profiles / 2)
957 freq_dc = int(num_profiles / 2)
928 # Flip con for
958 # Flip con for
929 for j in range(num_profiles):
959 for j in range(num_profiles):
930 jspectra_tmp[num_profiles-j-1]= jspectra[j]
960 jspectra_tmp[num_profiles-j-1]= jspectra[j]
931 # Intercambio perfil de DC con perfil inmediato anterior
961 # Intercambio perfil de DC con perfil inmediato anterior
932 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
962 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
933 jspectra_tmp[freq_dc]= jspectra[freq_dc]
963 jspectra_tmp[freq_dc]= jspectra[freq_dc]
934 # canal modificado es re-escrito en el arreglo de canales
964 # canal modificado es re-escrito en el arreglo de canales
935 self.dataOut.data_spc[2] = jspectra_tmp
965 self.dataOut.data_spc[2] = jspectra_tmp
936
966
937 return self.dataOut No newline at end of file
967 return self.dataOut
General Comments 0
You need to be logged in to leave comments. Login now