##// END OF EJS Templates
schain
RSS
atributo chann en class dopplerFlip
Alexander Valdez -
r1697:40634503f931
parent child
Show More
Side by Side Unified
Context file:
r1697:40634503f931 Loading diff...:

The requested changes are too big and content was truncated. Show full diff

Show file before | Show file after | Hide comments
@@ -0,0 +1,223
1 '''
2 Base clases to create Processing units and operations, the MPDecorator
3 must be used in plotting and writing operations to allow to run as an
4 external process.
5 '''
6 import os
7 import inspect
8 import zmq
9 import time
10 import pickle
11 import traceback
12 from threading import Thread
13 from multiprocessing import Process, Queue
14 from schainpy.utils import log
15 #isr-jro_proc_base.py
16 import copy
17 QUEUE_SIZE = int(os.environ.get('QUEUE_MAX_SIZE', '10'))
18
19 class ProcessingUnit(object):
20 '''
21 Base class to create Signal Chain Units
22 '''
23
24 proc_type = 'processing'
25
26 def __init__(self):
27
28 self.dataIn = None
29 self.dataOut = None
30 self.isConfig = False
31 self.operations = []
32 self.name = 'Test'
33 self.inputs = []
34
35 def setInput(self, unit):
36
37 attr = 'dataIn'
38 for i, u in enumerate(unit):
39 if i==0:
40 self.dataIn = u.dataOut#.copy()
41 self.inputs.append('dataIn')
42 else:
43 setattr(self, 'dataIn{}'.format(i), u.dataOut)#.copy())
44 self.inputs.append('dataIn{}'.format(i))
45
46 def getAllowedArgs(self):
47 if hasattr(self, '__attrs__'):
48 return self.__attrs__
49 else:
50 return inspect.getargspec(self.run).args
51
52 def addOperation(self, conf, operation):
53 '''
54 '''
55
56 self.operations.append((operation, conf.type, conf.getKwargs()))
57
58 def getOperationObj(self, objId):
59
60 if objId not in list(self.operations.keys()):
61 return None
62
63 return self.operations[objId]
64
65 def call(self, **kwargs):
66 '''
67 '''
68
69 try:
70 if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
71 if self.dataIn.runNextUnit:
72 return not self.dataIn.isReady()
73 else:
74 return self.dataIn.isReady()
75 elif self.dataIn is None or not self.dataIn.error:
76 self.run(**kwargs)
77 elif self.dataIn.error:
78 self.dataOut.error = self.dataIn.error
79 self.dataOut.flagNoData = True
80 except:
81 err = traceback.format_exc()
82 if 'SchainWarning' in err:
83 log.warning(err.split('SchainWarning:')[-1].split('\n')[0].strip(), self.name)
84 elif 'SchainError' in err:
85 log.error(err.split('SchainError:')[-1].split('\n')[0].strip(), self.name)
86 else:
87 log.error(err, self.name)
88 self.dataOut.error = True
89 for op, optype, opkwargs in self.operations:
90 aux = self.dataOut.copy()
91 if optype == 'other' and not self.dataOut.flagNoData:
92 self.dataOut = op.run(self.dataOut, **opkwargs)
93 elif optype == 'external' and not self.dataOut.flagNoData:
94 op.queue.put(aux)
95 elif optype == 'external' and self.dataOut.error:
96 op.queue.put(aux)
97 try:
98 if self.dataOut.runNextUnit:
99 runNextUnit = self.dataOut.runNextUnit
100 else:
101 runNextUnit = self.dataOut.isReady()
102 except:
103 runNextUnit = self.dataOut.isReady()
104 return 'Error' if self.dataOut.error else runNextUnit# self.dataOut.isReady()
105
106 def setup(self):
107
108 raise NotImplementedError
109
110 def run(self):
111
112 raise NotImplementedError
113
114 def close(self):
115
116 return
117
118
119 class Operation(object):
120
121 '''
122 '''
123
124 proc_type = 'operation'
125
126 def __init__(self):
127
128 self.id = None
129 self.isConfig = False
130
131 if not hasattr(self, 'name'):
132 self.name = self.__class__.__name__
133
134 def getAllowedArgs(self):
135 if hasattr(self, '__attrs__'):
136 return self.__attrs__
137 else:
138 return inspect.getargspec(self.run).args
139
140 def setup(self):
141
142 self.isConfig = True
143
144 raise NotImplementedError
145
146 def run(self, dataIn, **kwargs):
147 """
148 Realiza las operaciones necesarias sobre la dataIn.data y actualiza los
149 atributos del objeto dataIn.
150
151 Input:
152
153 dataIn : objeto del tipo JROData
154
155 Return:
156
157 None
158
159 Affected:
160 __buffer : buffer de recepcion de datos.
161
162 """
163 if not self.isConfig:
164 self.setup(**kwargs)
165
166 raise NotImplementedError
167
168 def close(self):
169
170 return
171
172
173 def MPDecorator(BaseClass):
174 """
175 Multiprocessing class decorator
176
177 This function add multiprocessing features to a BaseClass.
178 """
179
180 class MPClass(BaseClass, Process):
181
182 def __init__(self, *args, **kwargs):
183 super(MPClass, self).__init__()
184 Process.__init__(self)
185
186 self.args = args
187 self.kwargs = kwargs
188 self.t = time.time()
189 self.op_type = 'external'
190 self.name = BaseClass.__name__
191 self.__doc__ = BaseClass.__doc__
192
193 if 'plot' in self.name.lower() and not self.name.endswith('_'):
194 self.name = '{}{}'.format(self.CODE.upper(), 'Plot')
195
196 self.start_time = time.time()
197 self.err_queue = args[3]
198 self.queue = Queue(maxsize=QUEUE_SIZE)
199 self.myrun = BaseClass.run
200
201 def run(self):
202
203 while True:
204
205 dataOut = self.queue.get()
206
207 if not dataOut.error:
208 try:
209 BaseClass.run(self, dataOut, **self.kwargs)
210 except:
211 err = traceback.format_exc()
212 log.error(err, self.name)
213 else:
214 break
215
216 self.close()
217
218 def close(self):
219
220 BaseClass.close(self)
221 log.success('Done...(Time:{:4.2f} secs)'.format(time.time() - self.start_time), self.name)
222
223 return MPClass No newline at end of file
Show file before | Show file after | Hide comments
1 NO CONTENT: new file 100644
NO CONTENT: new file 100644
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (982 lines changed) Show them Hide them
@@ -0,0 +1,982
1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 # All rights reserved.
3 #
4 # Distributed under the terms of the BSD 3-clause license.
5 """Spectra processing Unit and operations
6
7 Here you will find the processing unit `SpectraProc` and several operations
8 to work with Spectra data type
9 """
10
11 import time
12 import itertools
13
14 import numpy
15
16 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
17 from schainpy.model.data.jrodata import Spectra
18 from schainpy.model.data.jrodata import hildebrand_sekhon
19 from schainpy.utils import log
20
21
22 class SpectraProc(ProcessingUnit):
23
24 def __init__(self):
25
26 ProcessingUnit.__init__(self)
27
28 self.buffer = None
29 self.firstdatatime = None
30 self.profIndex = 0
31 self.dataOut = Spectra()
32 self.id_min = None
33 self.id_max = None
34 self.setupReq = False #Agregar a todas las unidades de proc
35
36 def __updateSpecFromVoltage(self):
37
38 self.dataOut.timeZone = self.dataIn.timeZone
39 self.dataOut.dstFlag = self.dataIn.dstFlag
40 self.dataOut.errorCount = self.dataIn.errorCount
41 self.dataOut.useLocalTime = self.dataIn.useLocalTime
42 try:
43 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
44 except:
45 pass
46 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
47 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
48 self.dataOut.channelList = self.dataIn.channelList
49 self.dataOut.heightList = self.dataIn.heightList
50 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
51 self.dataOut.nProfiles = self.dataOut.nFFTPoints
52 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
53 self.dataOut.utctime = self.firstdatatime
54 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
55 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
56 self.dataOut.flagShiftFFT = False
57 self.dataOut.nCohInt = self.dataIn.nCohInt
58 self.dataOut.nIncohInt = 1
59 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
60 self.dataOut.frequency = self.dataIn.frequency
61 self.dataOut.realtime = self.dataIn.realtime
62 self.dataOut.azimuth = self.dataIn.azimuth
63 self.dataOut.zenith = self.dataIn.zenith
64 self.dataOut.beam.codeList = self.dataIn.beam.codeList
65 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
66 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
67 self.dataOut.runNextUnit = self.dataIn.runNextUnit
68 try:
69 self.dataOut.step = self.dataIn.step
70 except:
71 pass
72
73 def __getFft(self):
74 """
75 Convierte valores de Voltaje a Spectra
76
77 Affected:
78 self.dataOut.data_spc
79 self.dataOut.data_cspc
80 self.dataOut.data_dc
81 self.dataOut.heightList
82 self.profIndex
83 self.buffer
84 self.dataOut.flagNoData
85 """
86 fft_volt = numpy.fft.fft(
87 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
88 fft_volt = fft_volt.astype(numpy.dtype('complex'))
89 dc = fft_volt[:, 0, :]
90
91 # calculo de self-spectra
92 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
93 spc = fft_volt * numpy.conjugate(fft_volt)
94 spc = spc.real
95
96 blocksize = 0
97 blocksize += dc.size
98 blocksize += spc.size
99
100 cspc = None
101 pairIndex = 0
102 if self.dataOut.pairsList != None:
103 # calculo de cross-spectra
104 cspc = numpy.zeros(
105 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
106 for pair in self.dataOut.pairsList:
107 if pair[0] not in self.dataOut.channelList:
108 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
109 str(pair), str(self.dataOut.channelList)))
110 if pair[1] not in self.dataOut.channelList:
111 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
112 str(pair), str(self.dataOut.channelList)))
113
114 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
115 numpy.conjugate(fft_volt[pair[1], :, :])
116 pairIndex += 1
117 blocksize += cspc.size
118
119 self.dataOut.data_spc = spc
120 self.dataOut.data_cspc = cspc
121 self.dataOut.data_dc = dc
122 self.dataOut.blockSize = blocksize
123 self.dataOut.flagShiftFFT = False
124
125 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False, runNextUnit = 0):
126
127 self.dataIn.runNextUnit = runNextUnit
128 if self.dataIn.type == "Spectra":
129
130 self.dataOut.copy(self.dataIn)
131 if shift_fft:
132 #desplaza a la derecha en el eje 2 determinadas posiciones
133 shift = int(self.dataOut.nFFTPoints/2)
134 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
135
136 if self.dataOut.data_cspc is not None:
137 #desplaza a la derecha en el eje 2 determinadas posiciones
138 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
139 if pairsList:
140 self.__selectPairs(pairsList)
141
142 elif self.dataIn.type == "Voltage":
143
144 self.dataOut.flagNoData = True
145
146 if nFFTPoints == None:
147 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
148
149 if nProfiles == None:
150 nProfiles = nFFTPoints
151 #print(self.dataOut.ipp)
152 #exit(1)
153 if ippFactor == None:
154 self.dataOut.ippFactor = 1
155 #if ippFactor is not None:
156 #self.dataOut.ippFactor = ippFactor
157 #print(ippFactor)
158 #print(self.dataOut.ippFactor)
159 #exit(1)
160
161 self.dataOut.nFFTPoints = nFFTPoints
162
163 if self.buffer is None:
164 self.buffer = numpy.zeros((self.dataIn.nChannels,
165 nProfiles,
166 self.dataIn.nHeights),
167 dtype='complex')
168
169 if self.dataIn.flagDataAsBlock:
170 nVoltProfiles = self.dataIn.data.shape[1]
171
172 if nVoltProfiles == nProfiles:
173 self.buffer = self.dataIn.data.copy()
174 self.profIndex = nVoltProfiles
175
176 elif nVoltProfiles < nProfiles:
177
178 if self.profIndex == 0:
179 self.id_min = 0
180 self.id_max = nVoltProfiles
181 #print(self.id_min)
182 #print(self.id_max)
183 #print(numpy.shape(self.buffer))
184 self.buffer[:, self.id_min:self.id_max,
185 :] = self.dataIn.data
186 self.profIndex += nVoltProfiles
187 self.id_min += nVoltProfiles
188 self.id_max += nVoltProfiles
189 else:
190 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
191 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
192 self.dataOut.flagNoData = True
193 else:
194 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
195 self.profIndex += 1
196
197 if self.firstdatatime == None:
198 self.firstdatatime = self.dataIn.utctime
199
200 if self.profIndex == nProfiles:
201 self.__updateSpecFromVoltage()
202 if pairsList == None:
203 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
204 else:
205 self.dataOut.pairsList = pairsList
206 self.__getFft()
207 self.dataOut.flagNoData = False
208 self.firstdatatime = None
209 self.profIndex = 0
210 else:
211 raise ValueError("The type of input object '%s' is not valid".format(
212 self.dataIn.type))
213
214
215 def __selectPairs(self, pairsList):
216
217 if not pairsList:
218 return
219
220 pairs = []
221 pairsIndex = []
222
223 for pair in pairsList:
224 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
225 continue
226 pairs.append(pair)
227 pairsIndex.append(pairs.index(pair))
228
229 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
230 self.dataOut.pairsList = pairs
231
232 return
233
234 def selectFFTs(self, minFFT, maxFFT ):
235 """
236 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
237 minFFT<= FFT <= maxFFT
238 """
239
240 if (minFFT > maxFFT):
241 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
242
243 if (minFFT < self.dataOut.getFreqRange()[0]):
244 minFFT = self.dataOut.getFreqRange()[0]
245
246 if (maxFFT > self.dataOut.getFreqRange()[-1]):
247 maxFFT = self.dataOut.getFreqRange()[-1]
248
249 minIndex = 0
250 maxIndex = 0
251 FFTs = self.dataOut.getFreqRange()
252
253 inda = numpy.where(FFTs >= minFFT)
254 indb = numpy.where(FFTs <= maxFFT)
255
256 try:
257 minIndex = inda[0][0]
258 except:
259 minIndex = 0
260
261 try:
262 maxIndex = indb[0][-1]
263 except:
264 maxIndex = len(FFTs)
265
266 self.selectFFTsByIndex(minIndex, maxIndex)
267
268 return 1
269
270 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
271 newheis = numpy.where(
272 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
273
274 if hei_ref != None:
275 newheis = numpy.where(self.dataOut.heightList > hei_ref)
276
277 minIndex = min(newheis[0])
278 maxIndex = max(newheis[0])
279 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
280 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
281
282 # determina indices
283 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
284 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
285 avg_dB = 10 * \
286 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
287 beacon_dB = numpy.sort(avg_dB)[-nheis:]
288 beacon_heiIndexList = []
289 for val in avg_dB.tolist():
290 if val >= beacon_dB[0]:
291 beacon_heiIndexList.append(avg_dB.tolist().index(val))
292
293 #data_spc = data_spc[:,:,beacon_heiIndexList]
294 data_cspc = None
295 if self.dataOut.data_cspc is not None:
296 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
297 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
298
299 data_dc = None
300 if self.dataOut.data_dc is not None:
301 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
302 #data_dc = data_dc[:,beacon_heiIndexList]
303
304 self.dataOut.data_spc = data_spc
305 self.dataOut.data_cspc = data_cspc
306 self.dataOut.data_dc = data_dc
307 self.dataOut.heightList = heightList
308 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
309
310 return 1
311
312 def selectFFTsByIndex(self, minIndex, maxIndex):
313 """
314
315 """
316
317 if (minIndex < 0) or (minIndex > maxIndex):
318 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
319
320 if (maxIndex >= self.dataOut.nProfiles):
321 maxIndex = self.dataOut.nProfiles-1
322
323 #Spectra
324 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
325
326 data_cspc = None
327 if self.dataOut.data_cspc is not None:
328 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
329
330 data_dc = None
331 if self.dataOut.data_dc is not None:
332 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
333
334 self.dataOut.data_spc = data_spc
335 self.dataOut.data_cspc = data_cspc
336 self.dataOut.data_dc = data_dc
337
338 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
339 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
340 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
341
342 return 1
343
344 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
345 # validacion de rango
346 print("NOISeeee")
347 if minHei == None:
348 minHei = self.dataOut.heightList[0]
349
350 if maxHei == None:
351 maxHei = self.dataOut.heightList[-1]
352
353 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
354 print('minHei: %.2f is out of the heights range' % (minHei))
355 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
356 minHei = self.dataOut.heightList[0]
357
358 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
359 print('maxHei: %.2f is out of the heights range' % (maxHei))
360 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
361 maxHei = self.dataOut.heightList[-1]
362
363 # validacion de velocidades
364 velrange = self.dataOut.getVelRange(1)
365
366 if minVel == None:
367 minVel = velrange[0]
368
369 if maxVel == None:
370 maxVel = velrange[-1]
371
372 if (minVel < velrange[0]) or (minVel > maxVel):
373 print('minVel: %.2f is out of the velocity range' % (minVel))
374 print('minVel is setting to %.2f' % (velrange[0]))
375 minVel = velrange[0]
376
377 if (maxVel > velrange[-1]) or (maxVel < minVel):
378 print('maxVel: %.2f is out of the velocity range' % (maxVel))
379 print('maxVel is setting to %.2f' % (velrange[-1]))
380 maxVel = velrange[-1]
381
382 # seleccion de indices para rango
383 minIndex = 0
384 maxIndex = 0
385 heights = self.dataOut.heightList
386
387 inda = numpy.where(heights >= minHei)
388 indb = numpy.where(heights <= maxHei)
389
390 try:
391 minIndex = inda[0][0]
392 except:
393 minIndex = 0
394
395 try:
396 maxIndex = indb[0][-1]
397 except:
398 maxIndex = len(heights)
399
400 if (minIndex < 0) or (minIndex > maxIndex):
401 raise ValueError("some value in (%d,%d) is not valid" % (
402 minIndex, maxIndex))
403
404 if (maxIndex >= self.dataOut.nHeights):
405 maxIndex = self.dataOut.nHeights - 1
406
407 # seleccion de indices para velocidades
408 indminvel = numpy.where(velrange >= minVel)
409 indmaxvel = numpy.where(velrange <= maxVel)
410 try:
411 minIndexVel = indminvel[0][0]
412 except:
413 minIndexVel = 0
414
415 try:
416 maxIndexVel = indmaxvel[0][-1]
417 except:
418 maxIndexVel = len(velrange)
419
420 # seleccion del espectro
421 data_spc = self.dataOut.data_spc[:,
422 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
423 # estimacion de ruido
424 noise = numpy.zeros(self.dataOut.nChannels)
425
426 for channel in range(self.dataOut.nChannels):
427 daux = data_spc[channel, :, :]
428 sortdata = numpy.sort(daux, axis=None)
429 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
430
431 self.dataOut.noise_estimation = noise.copy()
432
433 return 1
434
435 class GetSNR(Operation):
436 '''
437 Written by R. Flores
438 '''
439 """Operation to get SNR.
440
441 Parameters:
442 -----------
443
444 Example
445 --------
446
447 op = proc_unit.addOperation(name='GetSNR', optype='other')
448
449 """
450
451 def __init__(self, **kwargs):
452
453 Operation.__init__(self, **kwargs)
454
455
456 def run(self,dataOut):
457
458 #noise = dataOut.getNoise()
459 noise = dataOut.getNoise(ymin_index=-10) #RegiΓ³n superior donde solo deberΓ­a de haber ruido
460 #print("Noise: ", noise)
461 #print("Noise_dB: ", 10*numpy.log10(noise/dataOut.normFactor))
462 #print("Heights: ", dataOut.heightList)
463 #dataOut.data_snr = (dataOut.data_spc.sum(axis=1))/(noise[:,None]*dataOut.normFactor)
464 ################dataOut.data_snr = (dataOut.data_spc.sum(axis=1))/(noise[:,None]*dataOut.nFFTPoints) #Before 12Jan2023
465 #dataOut.data_snr = (dataOut.data_spc.sum(axis=1)-noise[:,None])/(noise[:,None])
466 dataOut.data_snr = (dataOut.data_spc.sum(axis=1)-noise[:,None]*dataOut.nFFTPoints)/(noise[:,None]*dataOut.nFFTPoints) #It works apparently
467 dataOut.snl = numpy.log10(dataOut.data_snr)
468 #print("snl: ", dataOut.snl)
469 #exit(1)
470 #print(dataOut.heightList[-11])
471 #print(numpy.shape(dataOut.heightList))
472 #print(dataOut.data_snr)
473 #print(dataOut.data_snr[0,-11])
474 #exit(1)
475 #dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.5, numpy.nan, dataOut.data_snr)
476 #dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.1, numpy.nan, dataOut.data_snr)
477 #dataOut.data_snr = numpy.where(10*numpy.log10(dataOut.data_snr)<.0, numpy.nan, dataOut.data_snr)
478 #dataOut.data_snr = numpy.where(dataOut.data_snr<.05, numpy.nan, dataOut.data_snr)
479 #dataOut.snl = numpy.where(dataOut.data_snr<.01, numpy.nan, dataOut.snl)
480 dataOut.snl = numpy.where(dataOut.snl<-1, numpy.nan, dataOut.snl)
481 '''
482 import matplotlib.pyplot as plt
483 #plt.plot(10*numpy.log10(dataOut.data_snr[0]),dataOut.heightList)
484 plt.plot(dataOut.data_snr[0],dataOut.heightList)#,marker='*')
485 plt.xlim(-1,10)
486 plt.axvline(1,color='k')
487 plt.axvline(.1,color='k',linestyle='--')
488 plt.grid()
489 plt.show()
490 '''
491 #dataOut.data_snr = 10*numpy.log10(dataOut.data_snr)
492 #dataOut.data_snr = numpy.expand_dims(dataOut.data_snr,axis=0)
493 #print(dataOut.data_snr.shape)
494 #exit(1)
495 #print("Before: ", dataOut.data_snr[0])
496
497
498 return dataOut
499
500 class removeDC(Operation):
501
502 def run(self, dataOut, mode=2):
503 self.dataOut = dataOut
504 jspectra = self.dataOut.data_spc
505 jcspectra = self.dataOut.data_cspc
506
507 num_chan = jspectra.shape[0]
508 num_hei = jspectra.shape[2]
509
510 if jcspectra is not None:
511 jcspectraExist = True
512 num_pairs = jcspectra.shape[0]
513 else:
514 jcspectraExist = False
515
516 freq_dc = int(jspectra.shape[1] / 2)
517 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
518 ind_vel = ind_vel.astype(int)
519
520 if ind_vel[0] < 0:
521 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
522
523 if mode == 1:
524 jspectra[:, freq_dc, :] = (
525 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
526
527 if jcspectraExist:
528 jcspectra[:, freq_dc, :] = (
529 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
530
531 if mode == 2:
532
533 vel = numpy.array([-2, -1, 1, 2])
534 xx = numpy.zeros([4, 4])
535
536 for fil in range(4):
537 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
538
539 xx_inv = numpy.linalg.inv(xx)
540 xx_aux = xx_inv[0, :]
541
542 for ich in range(num_chan):
543 yy = jspectra[ich, ind_vel, :]
544 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
545
546 junkid = jspectra[ich, freq_dc, :] <= 0
547 cjunkid = sum(junkid)
548
549 if cjunkid.any():
550 jspectra[ich, freq_dc, junkid.nonzero()] = (
551 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
552
553 if jcspectraExist:
554 for ip in range(num_pairs):
555 yy = jcspectra[ip, ind_vel, :]
556 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
557
558 self.dataOut.data_spc = jspectra
559 self.dataOut.data_cspc = jcspectra
560
561 return self.dataOut
562
563 class removeInterference(Operation):
564
565 def removeInterference2(self):
566
567 cspc = self.dataOut.data_cspc
568 spc = self.dataOut.data_spc
569 Heights = numpy.arange(cspc.shape[2])
570 realCspc = numpy.abs(cspc)
571
572 for i in range(cspc.shape[0]):
573 LinePower= numpy.sum(realCspc[i], axis=0)
574 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
575 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
576 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
577 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
578 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
579
580
581 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
582 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
583 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
584 cspc[i,InterferenceRange,:] = numpy.NaN
585
586 self.dataOut.data_cspc = cspc
587
588 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
589
590 jspectra = self.dataOut.data_spc
591 jcspectra = self.dataOut.data_cspc
592 jnoise = self.dataOut.getNoise()
593 num_incoh = self.dataOut.nIncohInt
594
595 num_channel = jspectra.shape[0]
596 num_prof = jspectra.shape[1]
597 num_hei = jspectra.shape[2]
598
599 # hei_interf
600 if hei_interf is None:
601 count_hei = int(num_hei / 2)
602 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
603 hei_interf = numpy.asarray(hei_interf)[0]
604 # nhei_interf
605 if (nhei_interf == None):
606 nhei_interf = 5
607 if (nhei_interf < 1):
608 nhei_interf = 1
609 if (nhei_interf > count_hei):
610 nhei_interf = count_hei
611 if (offhei_interf == None):
612 offhei_interf = 0
613
614 ind_hei = list(range(num_hei))
615 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
616 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
617 mask_prof = numpy.asarray(list(range(num_prof)))
618 num_mask_prof = mask_prof.size
619 comp_mask_prof = [0, num_prof / 2]
620
621 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
622 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
623 jnoise = numpy.nan
624 noise_exist = jnoise[0] < numpy.Inf
625
626 # Subrutina de Remocion de la Interferencia
627 for ich in range(num_channel):
628 # Se ordena los espectros segun su potencia (menor a mayor)
629 power = jspectra[ich, mask_prof, :]
630 power = power[:, hei_interf]
631 power = power.sum(axis=0)
632 psort = power.ravel().argsort()
633
634 # Se estima la interferencia promedio en los Espectros de Potencia empleando
635 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
636 offhei_interf, nhei_interf + offhei_interf))]]]
637
638 if noise_exist:
639 # tmp_noise = jnoise[ich] / num_prof
640 tmp_noise = jnoise[ich]
641 junkspc_interf = junkspc_interf - tmp_noise
642 #junkspc_interf[:,comp_mask_prof] = 0
643
644 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
645 jspc_interf = jspc_interf.transpose()
646 # Calculando el espectro de interferencia promedio
647 noiseid = numpy.where(
648 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
649 noiseid = noiseid[0]
650 cnoiseid = noiseid.size
651 interfid = numpy.where(
652 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
653 interfid = interfid[0]
654 cinterfid = interfid.size
655
656 if (cnoiseid > 0):
657 jspc_interf[noiseid] = 0
658
659 # Expandiendo los perfiles a limpiar
660 if (cinterfid > 0):
661 new_interfid = (
662 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
663 new_interfid = numpy.asarray(new_interfid)
664 new_interfid = {x for x in new_interfid}
665 new_interfid = numpy.array(list(new_interfid))
666 new_cinterfid = new_interfid.size
667 else:
668 new_cinterfid = 0
669
670 for ip in range(new_cinterfid):
671 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
672 jspc_interf[new_interfid[ip]
673 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
674
675 jspectra[ich, :, ind_hei] = jspectra[ich, :,
676 ind_hei] - jspc_interf # Corregir indices
677
678 # Removiendo la interferencia del punto de mayor interferencia
679 ListAux = jspc_interf[mask_prof].tolist()
680 maxid = ListAux.index(max(ListAux))
681
682 if cinterfid > 0:
683 for ip in range(cinterfid * (interf == 2) - 1):
684 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
685 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
686 cind = len(ind)
687
688 if (cind > 0):
689 jspectra[ich, interfid[ip], ind] = tmp_noise * \
690 (1 + (numpy.random.uniform(cind) - 0.5) /
691 numpy.sqrt(num_incoh))
692
693 ind = numpy.array([-2, -1, 1, 2])
694 xx = numpy.zeros([4, 4])
695
696 for id1 in range(4):
697 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
698
699 xx_inv = numpy.linalg.inv(xx)
700 xx = xx_inv[:, 0]
701 ind = (ind + maxid + num_mask_prof) % num_mask_prof
702 yy = jspectra[ich, mask_prof[ind], :]
703 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
704 yy.transpose(), xx)
705
706 indAux = (jspectra[ich, :, :] < tmp_noise *
707 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
708 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
709 (1 - 1 / numpy.sqrt(num_incoh))
710
711 # Remocion de Interferencia en el Cross Spectra
712 if jcspectra is None:
713 return jspectra, jcspectra
714 num_pairs = int(jcspectra.size / (num_prof * num_hei))
715 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
716
717 for ip in range(num_pairs):
718
719 #-------------------------------------------
720
721 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
722 cspower = cspower[:, hei_interf]
723 cspower = cspower.sum(axis=0)
724
725 cspsort = cspower.ravel().argsort()
726 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
727 offhei_interf, nhei_interf + offhei_interf))]]]
728 junkcspc_interf = junkcspc_interf.transpose()
729 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
730
731 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
732
733 median_real = int(numpy.median(numpy.real(
734 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
735 median_imag = int(numpy.median(numpy.imag(
736 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
737 comp_mask_prof = [int(e) for e in comp_mask_prof]
738 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
739 median_real, median_imag)
740
741 for iprof in range(num_prof):
742 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
743 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
744
745 # Removiendo la Interferencia
746 jcspectra[ip, :, ind_hei] = jcspectra[ip,
747 :, ind_hei] - jcspc_interf
748
749 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
750 maxid = ListAux.index(max(ListAux))
751
752 ind = numpy.array([-2, -1, 1, 2])
753 xx = numpy.zeros([4, 4])
754
755 for id1 in range(4):
756 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
757
758 xx_inv = numpy.linalg.inv(xx)
759 xx = xx_inv[:, 0]
760
761 ind = (ind + maxid + num_mask_prof) % num_mask_prof
762 yy = jcspectra[ip, mask_prof[ind], :]
763 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
764
765 # Guardar Resultados
766 self.dataOut.data_spc = jspectra
767 self.dataOut.data_cspc = jcspectra
768
769 return 1
770
771 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
772
773 self.dataOut = dataOut
774
775 if mode == 1:
776 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
777 elif mode == 2:
778 self.removeInterference2()
779
780 return self.dataOut
781
782
783 class IncohInt(Operation):
784
785 __profIndex = 0
786 __withOverapping = False
787
788 __byTime = False
789 __initime = None
790 __lastdatatime = None
791 __integrationtime = None
792
793 __buffer_spc = None
794 __buffer_cspc = None
795 __buffer_dc = None
796
797 __dataReady = False
798
799 __timeInterval = None
800
801 n = None
802
803 def __init__(self):
804
805 Operation.__init__(self)
806
807 def setup(self, n=None, timeInterval=None, overlapping=False):
808 """
809 Set the parameters of the integration class.
810
811 Inputs:
812
813 n : Number of coherent integrations
814 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
815 overlapping :
816
817 """
818
819 self.__initime = None
820 self.__lastdatatime = 0
821
822 self.__buffer_spc = 0
823 self.__buffer_cspc = 0
824 self.__buffer_dc = 0
825
826 self.__profIndex = 0
827 self.__dataReady = False
828 self.__byTime = False
829
830 if n is None and timeInterval is None:
831 raise ValueError("n or timeInterval should be specified ...")
832
833 if n is not None:
834 self.n = int(n)
835 else:
836
837 self.__integrationtime = int(timeInterval)
838 self.n = None
839 self.__byTime = True
840
841 def putData(self, data_spc, data_cspc, data_dc):
842 """
843 Add a profile to the __buffer_spc and increase in one the __profileIndex
844
845 """
846
847 self.__buffer_spc += data_spc
848
849 if data_cspc is None:
850 self.__buffer_cspc = None
851 else:
852 self.__buffer_cspc += data_cspc
853
854 if data_dc is None:
855 self.__buffer_dc = None
856 else:
857 self.__buffer_dc += data_dc
858
859 self.__profIndex += 1
860
861 return
862
863 def pushData(self):
864 """
865 Return the sum of the last profiles and the profiles used in the sum.
866
867 Affected:
868
869 self.__profileIndex
870
871 """
872
873 data_spc = self.__buffer_spc
874 data_cspc = self.__buffer_cspc
875 data_dc = self.__buffer_dc
876 n = self.__profIndex
877
878 self.__buffer_spc = 0
879 self.__buffer_cspc = 0
880 self.__buffer_dc = 0
881 self.__profIndex = 0
882
883 return data_spc, data_cspc, data_dc, n
884
885 def byProfiles(self, *args):
886
887 self.__dataReady = False
888 avgdata_spc = None
889 avgdata_cspc = None
890 avgdata_dc = None
891
892 self.putData(*args)
893
894 if self.__profIndex == self.n:
895
896 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
897 self.n = n
898 self.__dataReady = True
899
900 return avgdata_spc, avgdata_cspc, avgdata_dc
901
902 def byTime(self, datatime, *args):
903
904 self.__dataReady = False
905 avgdata_spc = None
906 avgdata_cspc = None
907 avgdata_dc = None
908
909 self.putData(*args)
910
911 if (datatime - self.__initime) >= self.__integrationtime:
912 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
913 self.n = n
914 self.__dataReady = True
915
916 return avgdata_spc, avgdata_cspc, avgdata_dc
917
918 def integrate(self, datatime, *args):
919
920 if self.__profIndex == 0:
921 self.__initime = datatime
922
923 if self.__byTime:
924 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
925 datatime, *args)
926 else:
927 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
928
929 if not self.__dataReady:
930 return None, None, None, None
931
932 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
933
934 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
935 if n == 1:
936 return dataOut
937 print("JERE")
938 dataOut.flagNoData = True
939
940 if not self.isConfig:
941 self.setup(n, timeInterval, overlapping)
942 self.isConfig = True
943
944 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
945 dataOut.data_spc,
946 dataOut.data_cspc,
947 dataOut.data_dc)
948
949 if self.__dataReady:
950
951 dataOut.data_spc = avgdata_spc
952 print(numpy.sum(dataOut.data_spc))
953 exit(1)
954 dataOut.data_cspc = avgdata_cspc
955 dataOut.data_dc = avgdata_dc
956 dataOut.nIncohInt *= self.n
957 dataOut.utctime = avgdatatime
958 dataOut.flagNoData = False
959
960 return dataOut
961
962 class dopplerFlip(Operation):
963
964 def run(self, dataOut, chann = None):
965 # arreglo 1: (num_chan, num_profiles, num_heights)
966 self.dataOut = dataOut
967 # JULIA-oblicua, indice 2
968 # arreglo 2: (num_profiles, num_heights)
969 jspectra = self.dataOut.data_spc[chann]
970 jspectra_tmp = numpy.zeros(jspectra.shape)
971 num_profiles = jspectra.shape[0]
972 freq_dc = int(num_profiles / 2)
973 # Flip con for
974 for j in range(num_profiles):
975 jspectra_tmp[num_profiles-j-1]= jspectra[j]
976 # Intercambio perfil de DC con perfil inmediato anterior
977 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
978 jspectra_tmp[freq_dc]= jspectra[freq_dc]
979 # canal modificado es re-escrito en el arreglo de canales
980 self.dataOut.data_spc[chann] = jspectra_tmp
981
982 return self.dataOut
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (1644 lines changed) Show them Hide them
@@ -0,0 +1,1644
1
2 import os
3 import sys
4 import numpy, math
5 from scipy import interpolate
6 from scipy.optimize import nnls
7 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
8 from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
9 from schainpy.utils import log
10 from time import time, mktime, strptime, gmtime, ctime
11 from scipy.optimize import least_squares
12 import datetime
13 import collections.abc
14
15 try:
16 from schainpy.model.proc import fitacf_guess
17 from schainpy.model.proc import fitacf_fit_short
18 from schainpy.model.proc import fitacf_acf2
19 from schainpy.model.proc import full_profile_profile
20 except:
21 log.warning('Missing Faraday fortran libs')
22
23 class VoltageProc(ProcessingUnit):
24
25 def __init__(self):
26
27 ProcessingUnit.__init__(self)
28
29 self.dataOut = Voltage()
30 self.flip = 1
31 self.setupReq = False
32
33 def run(self, runNextUnit = 0):
34
35 if self.dataIn.type == 'AMISR':
36 self.__updateObjFromAmisrInput()
37
38 if self.dataIn.type == 'Voltage':
39 self.dataOut.copy(self.dataIn)
40 self.dataOut.runNextUnit = runNextUnit
41
42 def __updateObjFromAmisrInput(self):
43
44 self.dataOut.timeZone = self.dataIn.timeZone
45 self.dataOut.dstFlag = self.dataIn.dstFlag
46 self.dataOut.errorCount = self.dataIn.errorCount
47 self.dataOut.useLocalTime = self.dataIn.useLocalTime
48
49 self.dataOut.flagNoData = self.dataIn.flagNoData
50 self.dataOut.data = self.dataIn.data
51 self.dataOut.utctime = self.dataIn.utctime
52 self.dataOut.channelList = self.dataIn.channelList
53 #self.dataOut.timeInterval = self.dataIn.timeInterval
54 self.dataOut.heightList = self.dataIn.heightList
55 self.dataOut.nProfiles = self.dataIn.nProfiles
56
57 self.dataOut.nCohInt = self.dataIn.nCohInt
58 self.dataOut.ippSeconds = self.dataIn.ippSeconds
59 self.dataOut.frequency = self.dataIn.frequency
60
61 self.dataOut.azimuth = self.dataIn.azimuth
62 self.dataOut.zenith = self.dataIn.zenith
63
64 self.dataOut.beam.codeList = self.dataIn.beam.codeList
65 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
66 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
67
68
69 class selectChannels(Operation):
70
71 def run(self, dataOut, channelList):
72
73 channelIndexList = []
74 self.dataOut = dataOut
75 for channel in channelList:
76 if channel not in self.dataOut.channelList:
77 raise ValueError("Channel %d is not in %s" %(channel, str(self.dataOut.channelList)))
78
79 index = self.dataOut.channelList.index(channel)
80 channelIndexList.append(index)
81 self.selectChannelsByIndex(channelIndexList)
82 return self.dataOut
83
84 def selectChannelsByIndex(self, channelIndexList):
85 """
86 Selecciona un bloque de datos en base a canales segun el channelIndexList
87
88 Input:
89 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
90
91 Affected:
92 self.dataOut.data
93 self.dataOut.channelIndexList
94 self.dataOut.nChannels
95 self.dataOut.m_ProcessingHeader.totalSpectra
96 self.dataOut.systemHeaderObj.numChannels
97 self.dataOut.m_ProcessingHeader.blockSize
98
99 Return:
100 None
101 """
102
103 for channelIndex in channelIndexList:
104 if channelIndex not in self.dataOut.channelIndexList:
105 raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
106
107 if self.dataOut.type == 'Voltage':
108 if self.dataOut.flagDataAsBlock:
109 """
110 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
111 """
112 data = self.dataOut.data[channelIndexList,:,:]
113 else:
114 data = self.dataOut.data[channelIndexList,:]
115
116 self.dataOut.data = data
117 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
118 self.dataOut.channelList = range(len(channelIndexList))
119
120 elif self.dataOut.type == 'Spectra':
121 data_spc = self.dataOut.data_spc[channelIndexList, :]
122 data_dc = self.dataOut.data_dc[channelIndexList, :]
123
124 self.dataOut.data_spc = data_spc
125 self.dataOut.data_dc = data_dc
126
127 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
128 self.dataOut.channelList = range(len(channelIndexList))
129 self.__selectPairsByChannel(channelIndexList)
130
131 return 1
132
133 def __selectPairsByChannel(self, channelList=None):
134
135 if channelList == None:
136 return
137
138 pairsIndexListSelected = []
139 for pairIndex in self.dataOut.pairsIndexList:
140 # First pair
141 if self.dataOut.pairsList[pairIndex][0] not in channelList:
142 continue
143 # Second pair
144 if self.dataOut.pairsList[pairIndex][1] not in channelList:
145 continue
146
147 pairsIndexListSelected.append(pairIndex)
148
149 if not pairsIndexListSelected:
150 self.dataOut.data_cspc = None
151 self.dataOut.pairsList = []
152 return
153
154 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
155 self.dataOut.pairsList = [self.dataOut.pairsList[i]
156 for i in pairsIndexListSelected]
157
158 return
159
160 class selectHeights(Operation):
161
162 def run(self, dataOut, minHei=None, maxHei=None, minIndex=None, maxIndex=None):
163 """
164 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
165 minHei <= height <= maxHei
166
167 Input:
168 minHei : valor minimo de altura a considerar
169 maxHei : valor maximo de altura a considerar
170
171 Affected:
172 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
173
174 Return:
175 1 si el metodo se ejecuto con exito caso contrario devuelve 0
176 """
177
178 self.dataOut = dataOut
179
180 if minHei and maxHei:
181
182 if (minHei < self.dataOut.heightList[0]):
183 minHei = self.dataOut.heightList[0]
184
185 if (maxHei > self.dataOut.heightList[-1]):
186 maxHei = self.dataOut.heightList[-1]
187
188 minIndex = 0
189 maxIndex = 0
190 heights = self.dataOut.heightList
191
192 inda = numpy.where(heights >= minHei)
193 indb = numpy.where(heights <= maxHei)
194
195 try:
196 minIndex = inda[0][0]
197 except:
198 minIndex = 0
199
200 try:
201 maxIndex = indb[0][-1]
202 except:
203 maxIndex = len(heights)
204
205 self.selectHeightsByIndex(minIndex, maxIndex)
206
207 return self.dataOut
208
209 def selectHeightsByIndex(self, minIndex, maxIndex):
210 """
211 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
212 minIndex <= index <= maxIndex
213
214 Input:
215 minIndex : valor de indice minimo de altura a considerar
216 maxIndex : valor de indice maximo de altura a considerar
217
218 Affected:
219 self.dataOut.data
220 self.dataOut.heightList
221
222 Return:
223 1 si el metodo se ejecuto con exito caso contrario devuelve 0
224 """
225
226 if self.dataOut.type == 'Voltage':
227 if (minIndex < 0) or (minIndex > maxIndex):
228 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
229
230 if (maxIndex >= self.dataOut.nHeights):
231 maxIndex = self.dataOut.nHeights
232
233 #voltage
234 if self.dataOut.flagDataAsBlock:
235 """
236 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
237 """
238 data = self.dataOut.data[:,:, minIndex:maxIndex]
239 else:
240 data = self.dataOut.data[:, minIndex:maxIndex]
241
242 # firstHeight = self.dataOut.heightList[minIndex]
243
244 self.dataOut.data = data
245 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
246
247 if self.dataOut.nHeights <= 1:
248 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
249 elif self.dataOut.type == 'Spectra':
250 if (minIndex < 0) or (minIndex > maxIndex):
251 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
252 minIndex, maxIndex))
253
254 if (maxIndex >= self.dataOut.nHeights):
255 maxIndex = self.dataOut.nHeights - 1
256
257 # Spectra
258 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
259
260 data_cspc = None
261 if self.dataOut.data_cspc is not None:
262 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
263
264 data_dc = None
265 if self.dataOut.data_dc is not None:
266 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
267
268 self.dataOut.data_spc = data_spc
269 self.dataOut.data_cspc = data_cspc
270 self.dataOut.data_dc = data_dc
271
272 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
273
274 return 1
275
276
277 class filterByHeights(Operation):
278
279 def run(self, dataOut, window):
280
281 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
282
283 if window == None:
284 window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
285
286 newdelta = deltaHeight * window
287 r = dataOut.nHeights % window
288 newheights = (dataOut.nHeights-r)/window
289
290 if newheights <= 1:
291 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
292
293 if dataOut.flagDataAsBlock:
294 """
295 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
296 """
297 buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)]
298 buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window)
299 buffer = numpy.sum(buffer,3)
300
301 else:
302 buffer = dataOut.data[:,0:int(dataOut.nHeights-r)]
303 buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
304 buffer = numpy.sum(buffer,2)
305
306 dataOut.data = buffer
307 dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta
308 dataOut.windowOfFilter = window
309
310 return dataOut
311
312
313 class setH0(Operation):
314
315 def run(self, dataOut, h0, deltaHeight = None):
316
317 if not deltaHeight:
318 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
319
320 nHeights = dataOut.nHeights
321
322 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
323
324 dataOut.heightList = newHeiRange
325
326 return dataOut
327
328
329 class deFlip(Operation):
330 def __init__(self):
331
332 self.flip = 1
333
334 def run(self, dataOut, channelList = []):
335
336 data = dataOut.data.copy()
337
338 if channelList==1: #PARCHE #Lista de un solo canal produce error
339 channelList=[1]
340
341 dataOut.FlipChannels=channelList
342 if dataOut.flagDataAsBlock:
343 flip = self.flip
344 profileList = list(range(dataOut.nProfiles))
345
346 if not channelList:
347 for thisProfile in profileList:
348 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
349 flip *= -1.0
350 else:
351 for thisChannel in channelList:
352 if thisChannel not in dataOut.channelList:
353 continue
354
355 for thisProfile in profileList:
356 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
357 flip *= -1.0
358
359 self.flip = flip
360
361 else:
362 if not channelList:
363 data[:,:] = data[:,:]*self.flip
364 else:
365 for thisChannel in channelList:
366 if thisChannel not in dataOut.channelList:
367 continue
368
369 data[thisChannel,:] = data[thisChannel,:]*self.flip
370
371 self.flip *= -1.
372
373 dataOut.data = data
374
375 return dataOut
376
377
378 class setAttribute(Operation):
379 '''
380 Set an arbitrary attribute(s) to dataOut
381 '''
382
383 def __init__(self):
384
385 Operation.__init__(self)
386 self._ready = False
387
388 def run(self, dataOut, **kwargs):
389
390 for key, value in kwargs.items():
391 setattr(dataOut, key, value)
392
393 return dataOut
394
395
396 @MPDecorator
397 class printAttribute(Operation):
398 '''
399 Print an arbitrary attribute of dataOut
400 '''
401
402 def __init__(self):
403
404 Operation.__init__(self)
405
406 def run(self, dataOut, attributes):
407
408 if isinstance(attributes, str):
409 attributes = [attributes]
410 for attr in attributes:
411 if hasattr(dataOut, attr):
412 log.log(getattr(dataOut, attr), attr)
413
414
415 class interpolateHeights(Operation):
416
417 def run(self, dataOut, topLim, botLim):
418 #69 al 72 para julia
419 #82-84 para meteoros
420 if len(numpy.shape(dataOut.data))==2:
421 sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
422 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
423 #dataOut.data[:,botLim:limSup+1] = sampInterp
424 dataOut.data[:,botLim:topLim+1] = sampInterp
425 else:
426 nHeights = dataOut.data.shape[2]
427 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
428 y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
429 f = interpolate.interp1d(x, y, axis = 2)
430 xnew = numpy.arange(botLim,topLim+1)
431 ynew = f(xnew)
432 dataOut.data[:,:,botLim:topLim+1] = ynew
433
434 return dataOut
435
436
437 class CohInt(Operation):
438
439 isConfig = False
440 __profIndex = 0
441 __byTime = False
442 __initime = None
443 __lastdatatime = None
444 __integrationtime = None
445 __buffer = None
446 __bufferStride = []
447 __dataReady = False
448 __profIndexStride = 0
449 __dataToPutStride = False
450 n = None
451
452 def __init__(self, **kwargs):
453
454 Operation.__init__(self, **kwargs)
455
456 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
457 """
458 Set the parameters of the integration class.
459
460 Inputs:
461
462 n : Number of coherent integrations
463 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
464 overlapping :
465 """
466
467 self.__initime = None
468 self.__lastdatatime = 0
469 self.__buffer = None
470 self.__dataReady = False
471 self.byblock = byblock
472 self.stride = stride
473
474 if n == None and timeInterval == None:
475 raise ValueError("n or timeInterval should be specified ...")
476
477 if n != None:
478 self.n = n
479 self.__byTime = False
480 else:
481 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
482 self.n = 9999
483 self.__byTime = True
484
485 if overlapping:
486 self.__withOverlapping = True
487 self.__buffer = None
488 else:
489 self.__withOverlapping = False
490 self.__buffer = 0
491
492 self.__profIndex = 0
493
494 def putData(self, data):
495
496 """
497 Add a profile to the __buffer and increase in one the __profileIndex
498
499 """
500
501 if not self.__withOverlapping:
502 self.__buffer += data.copy()
503 self.__profIndex += 1
504 return
505
506 #Overlapping data
507 nChannels, nHeis = data.shape
508 data = numpy.reshape(data, (1, nChannels, nHeis))
509
510 #If the buffer is empty then it takes the data value
511 if self.__buffer is None:
512 self.__buffer = data
513 self.__profIndex += 1
514 return
515
516 #If the buffer length is lower than n then stakcing the data value
517 if self.__profIndex < self.n:
518 self.__buffer = numpy.vstack((self.__buffer, data))
519 self.__profIndex += 1
520 return
521
522 #If the buffer length is equal to n then replacing the last buffer value with the data value
523 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
524 self.__buffer[self.n-1] = data
525 self.__profIndex = self.n
526 return
527
528
529 def pushData(self):
530 """
531 Return the sum of the last profiles and the profiles used in the sum.
532
533 Affected:
534
535 self.__profileIndex
536
537 """
538
539 if not self.__withOverlapping:
540 data = self.__buffer
541 n = self.__profIndex
542
543 self.__buffer = 0
544 self.__profIndex = 0
545
546 return data, n
547
548 #Integration with Overlapping
549 data = numpy.sum(self.__buffer, axis=0)
550 # print data
551 # raise
552 n = self.__profIndex
553
554 return data, n
555
556 def byProfiles(self, data):
557
558 self.__dataReady = False
559 avgdata = None
560 # n = None
561 # print data
562 # raise
563 self.putData(data)
564
565 if self.__profIndex == self.n:
566 avgdata, n = self.pushData()
567 self.__dataReady = True
568
569 return avgdata
570
571 def byTime(self, data, datatime):
572
573 self.__dataReady = False
574 avgdata = None
575 n = None
576
577 self.putData(data)
578
579 if (datatime - self.__initime) >= self.__integrationtime:
580 avgdata, n = self.pushData()
581 self.n = n
582 self.__dataReady = True
583
584 return avgdata
585
586 def integrateByStride(self, data, datatime):
587 # print data
588 if self.__profIndex == 0:
589 self.__buffer = [[data.copy(), datatime]]
590 else:
591 self.__buffer.append([data.copy(),datatime])
592 self.__profIndex += 1
593 self.__dataReady = False
594
595 if self.__profIndex == self.n * self.stride :
596 self.__dataToPutStride = True
597 self.__profIndexStride = 0
598 self.__profIndex = 0
599 self.__bufferStride = []
600 for i in range(self.stride):
601 current = self.__buffer[i::self.stride]
602 data = numpy.sum([t[0] for t in current], axis=0)
603 avgdatatime = numpy.average([t[1] for t in current])
604 # print data
605 self.__bufferStride.append((data, avgdatatime))
606
607 if self.__dataToPutStride:
608 self.__dataReady = True
609 self.__profIndexStride += 1
610 if self.__profIndexStride == self.stride:
611 self.__dataToPutStride = False
612 # print self.__bufferStride[self.__profIndexStride - 1]
613 # raise
614 return self.__bufferStride[self.__profIndexStride - 1]
615
616
617 return None, None
618
619 def integrate(self, data, datatime=None):
620
621 if self.__initime == None:
622 self.__initime = datatime
623
624 if self.__byTime:
625 avgdata = self.byTime(data, datatime)
626 else:
627 avgdata = self.byProfiles(data)
628
629
630 self.__lastdatatime = datatime
631
632 if avgdata is None:
633 return None, None
634
635 avgdatatime = self.__initime
636
637 deltatime = datatime - self.__lastdatatime
638
639 if not self.__withOverlapping:
640 self.__initime = datatime
641 else:
642 self.__initime += deltatime
643
644 return avgdata, avgdatatime
645
646 def integrateByBlock(self, dataOut):
647
648 times = int(dataOut.data.shape[1]/self.n)
649 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
650
651 id_min = 0
652 id_max = self.n
653
654 for i in range(times):
655 junk = dataOut.data[:,id_min:id_max,:]
656 avgdata[:,i,:] = junk.sum(axis=1)
657 id_min += self.n
658 id_max += self.n
659
660 timeInterval = dataOut.ippSeconds*self.n
661 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
662 self.__dataReady = True
663 return avgdata, avgdatatime
664
665 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
666
667 if not self.isConfig:
668 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
669 self.isConfig = True
670 if dataOut.flagDataAsBlock:
671 """
672 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
673 """
674 avgdata, avgdatatime = self.integrateByBlock(dataOut)
675 dataOut.nProfiles /= self.n
676 else:
677 if stride is None:
678 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
679 else:
680 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
681
682
683 # dataOut.timeInterval *= n
684 dataOut.flagNoData = True
685
686 if self.__dataReady:
687 dataOut.data = avgdata
688 if not dataOut.flagCohInt:
689 dataOut.nCohInt *= self.n
690 dataOut.flagCohInt = True
691 dataOut.utctime = avgdatatime
692 # print avgdata, avgdatatime
693 # raise
694 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
695 dataOut.flagNoData = False
696 return dataOut
697
698 class Decoder(Operation):
699
700 isConfig = False
701 __profIndex = 0
702
703 code = None
704
705 nCode = None
706 nBaud = None
707
708 def __init__(self, **kwargs):
709
710 Operation.__init__(self, **kwargs)
711
712 self.times = None
713 self.osamp = None
714 # self.__setValues = False
715 self.isConfig = False
716 self.setupReq = False
717 def setup(self, code, osamp, dataOut):
718
719 self.__profIndex = 0
720
721 self.code = code
722
723 self.nCode = len(code)
724 self.nBaud = len(code[0])
725
726 if (osamp != None) and (osamp >1):
727 self.osamp = osamp
728 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
729 self.nBaud = self.nBaud*self.osamp
730
731 self.__nChannels = dataOut.nChannels
732 self.__nProfiles = dataOut.nProfiles
733 self.__nHeis = dataOut.nHeights
734
735 if self.__nHeis < self.nBaud:
736 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
737
738 #Frequency
739 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
740
741 __codeBuffer[:,0:self.nBaud] = self.code
742
743 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
744
745 if dataOut.flagDataAsBlock:
746
747 self.ndatadec = self.__nHeis #- self.nBaud + 1
748
749 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
750
751 else:
752
753 #Time
754 self.ndatadec = self.__nHeis #- self.nBaud + 1
755
756 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
757
758 def __convolutionInFreq(self, data):
759
760 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
761
762 fft_data = numpy.fft.fft(data, axis=1)
763
764 conv = fft_data*fft_code
765
766 data = numpy.fft.ifft(conv,axis=1)
767
768 return data
769
770 def __convolutionInFreqOpt(self, data):
771
772 raise NotImplementedError
773
774 def __convolutionInTime(self, data):
775
776 code = self.code[self.__profIndex]
777 for i in range(self.__nChannels):
778 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
779
780 return self.datadecTime
781
782 def __convolutionByBlockInTime(self, data):
783
784 repetitions = int(self.__nProfiles / self.nCode)
785 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
786 junk = junk.flatten()
787 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
788 profilesList = range(self.__nProfiles)
789
790 for i in range(self.__nChannels):
791 for j in profilesList:
792 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
793 return self.datadecTime
794
795 def __convolutionByBlockInFreq(self, data):
796
797 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
798
799
800 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
801
802 fft_data = numpy.fft.fft(data, axis=2)
803
804 conv = fft_data*fft_code
805
806 data = numpy.fft.ifft(conv,axis=2)
807
808 return data
809
810
811 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
812
813 if dataOut.flagDecodeData:
814 print("This data is already decoded, recoding again ...")
815
816 if not self.isConfig:
817
818 if code is None:
819 if dataOut.code is None:
820 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
821
822 code = dataOut.code
823 else:
824 code = numpy.array(code).reshape(nCode,nBaud)
825 self.setup(code, osamp, dataOut)
826
827 self.isConfig = True
828
829 if mode == 3:
830 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
831
832 if times != None:
833 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
834
835 if self.code is None:
836 print("Fail decoding: Code is not defined.")
837 return
838
839 self.__nProfiles = dataOut.nProfiles
840 datadec = None
841
842 if mode == 3:
843 mode = 0
844
845 if dataOut.flagDataAsBlock:
846 """
847 Decoding when data have been read as block,
848 """
849
850 if mode == 0:
851 datadec = self.__convolutionByBlockInTime(dataOut.data)
852 if mode == 1:
853 datadec = self.__convolutionByBlockInFreq(dataOut.data)
854 else:
855 """
856 Decoding when data have been read profile by profile
857 """
858 if mode == 0:
859 datadec = self.__convolutionInTime(dataOut.data)
860
861 if mode == 1:
862 datadec = self.__convolutionInFreq(dataOut.data)
863
864 if mode == 2:
865 datadec = self.__convolutionInFreqOpt(dataOut.data)
866
867 if datadec is None:
868 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
869
870 dataOut.code = self.code
871 dataOut.nCode = self.nCode
872 dataOut.nBaud = self.nBaud
873
874 dataOut.data = datadec
875 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
876
877 dataOut.flagDecodeData = True #asumo q la data esta decodificada
878
879 if self.__profIndex == self.nCode-1:
880 self.__profIndex = 0
881 return dataOut
882
883 self.__profIndex += 1
884
885 return dataOut
886
887
888 class ProfileConcat(Operation):
889
890 isConfig = False
891 buffer = None
892
893 def __init__(self, **kwargs):
894
895 Operation.__init__(self, **kwargs)
896 self.profileIndex = 0
897
898 def reset(self):
899 self.buffer = numpy.zeros_like(self.buffer)
900 self.start_index = 0
901 self.times = 1
902
903 def setup(self, data, m, n=1):
904 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
905 self.nHeights = data.shape[1]#.nHeights
906 self.start_index = 0
907 self.times = 1
908
909 def concat(self, data):
910
911 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
912 self.start_index = self.start_index + self.nHeights
913
914 def run(self, dataOut, m):
915 dataOut.flagNoData = True
916
917 if not self.isConfig:
918 self.setup(dataOut.data, m, 1)
919 self.isConfig = True
920
921 if dataOut.flagDataAsBlock:
922 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
923
924 else:
925 self.concat(dataOut.data)
926 self.times += 1
927 if self.times > m:
928 dataOut.data = self.buffer
929 self.reset()
930 dataOut.flagNoData = False
931 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
932 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
933 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
934 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
935 dataOut.ippSeconds *= m
936 return dataOut
937
938 class ProfileSelector(Operation):
939
940 profileIndex = None
941 # Tamanho total de los perfiles
942 nProfiles = None
943
944 def __init__(self, **kwargs):
945
946 Operation.__init__(self, **kwargs)
947 self.profileIndex = 0
948
949 def incProfileIndex(self):
950
951 self.profileIndex += 1
952
953 if self.profileIndex >= self.nProfiles:
954 self.profileIndex = 0
955
956 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
957
958 if profileIndex < minIndex:
959 return False
960
961 if profileIndex > maxIndex:
962 return False
963
964 return True
965
966 def isThisProfileInList(self, profileIndex, profileList):
967
968 if profileIndex not in profileList:
969 return False
970
971 return True
972
973 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
974
975 """
976 ProfileSelector:
977
978 Inputs:
979 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
980
981 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
982
983 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
984
985 """
986
987 if rangeList is not None:
988 if type(rangeList[0]) not in (tuple, list):
989 rangeList = [rangeList]
990
991 dataOut.flagNoData = True
992
993 if dataOut.flagDataAsBlock:
994 """
995 data dimension = [nChannels, nProfiles, nHeis]
996 """
997 if profileList != None:
998 dataOut.data = dataOut.data[:,profileList,:]
999
1000 if profileRangeList != None:
1001 minIndex = profileRangeList[0]
1002 maxIndex = profileRangeList[1]
1003 profileList = list(range(minIndex, maxIndex+1))
1004
1005 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
1006
1007 if rangeList != None:
1008
1009 profileList = []
1010
1011 for thisRange in rangeList:
1012 minIndex = thisRange[0]
1013 maxIndex = thisRange[1]
1014
1015 profileList.extend(list(range(minIndex, maxIndex+1)))
1016
1017 dataOut.data = dataOut.data[:,profileList,:]
1018
1019 dataOut.nProfiles = len(profileList)
1020 dataOut.profileIndex = dataOut.nProfiles - 1
1021 dataOut.flagNoData = False
1022
1023 return dataOut
1024
1025 """
1026 data dimension = [nChannels, nHeis]
1027 """
1028
1029 if profileList != None:
1030
1031 if self.isThisProfileInList(dataOut.profileIndex, profileList):
1032
1033 self.nProfiles = len(profileList)
1034 dataOut.nProfiles = self.nProfiles
1035 dataOut.profileIndex = self.profileIndex
1036 dataOut.flagNoData = False
1037
1038 self.incProfileIndex()
1039 return dataOut
1040
1041 if profileRangeList != None:
1042
1043 minIndex = profileRangeList[0]
1044 maxIndex = profileRangeList[1]
1045
1046 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1047
1048 self.nProfiles = maxIndex - minIndex + 1
1049 dataOut.nProfiles = self.nProfiles
1050 dataOut.profileIndex = self.profileIndex
1051 dataOut.flagNoData = False
1052
1053 self.incProfileIndex()
1054 return dataOut
1055
1056 if rangeList != None:
1057
1058 nProfiles = 0
1059
1060 for thisRange in rangeList:
1061 minIndex = thisRange[0]
1062 maxIndex = thisRange[1]
1063
1064 nProfiles += maxIndex - minIndex + 1
1065
1066 for thisRange in rangeList:
1067
1068 minIndex = thisRange[0]
1069 maxIndex = thisRange[1]
1070
1071 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1072
1073 self.nProfiles = nProfiles
1074 dataOut.nProfiles = self.nProfiles
1075 dataOut.profileIndex = self.profileIndex
1076 dataOut.flagNoData = False
1077
1078 self.incProfileIndex()
1079
1080 break
1081
1082 return dataOut
1083
1084
1085 if beam != None: #beam is only for AMISR data
1086 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
1087 dataOut.flagNoData = False
1088 dataOut.profileIndex = self.profileIndex
1089
1090 self.incProfileIndex()
1091
1092 return dataOut
1093
1094 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
1095
1096 #return False
1097 return dataOut
1098
1099 class Reshaper(Operation):
1100
1101 def __init__(self, **kwargs):
1102
1103 Operation.__init__(self, **kwargs)
1104
1105 self.__buffer = None
1106 self.__nitems = 0
1107
1108 def __appendProfile(self, dataOut, nTxs):
1109
1110 if self.__buffer is None:
1111 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1112 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1113
1114 ini = dataOut.nHeights * self.__nitems
1115 end = ini + dataOut.nHeights
1116
1117 self.__buffer[:, ini:end] = dataOut.data
1118
1119 self.__nitems += 1
1120
1121 return int(self.__nitems*nTxs)
1122
1123 def __getBuffer(self):
1124
1125 if self.__nitems == int(1./self.__nTxs):
1126
1127 self.__nitems = 0
1128
1129 return self.__buffer.copy()
1130
1131 return None
1132
1133 def __checkInputs(self, dataOut, shape, nTxs):
1134
1135 if shape is None and nTxs is None:
1136 raise ValueError("Reshaper: shape of factor should be defined")
1137
1138 if nTxs:
1139 if nTxs < 0:
1140 raise ValueError("nTxs should be greater than 0")
1141
1142 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1143 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1144
1145 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1146
1147 return shape, nTxs
1148
1149 if len(shape) != 2 and len(shape) != 3:
1150 raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights))
1151
1152 if len(shape) == 2:
1153 shape_tuple = [dataOut.nChannels]
1154 shape_tuple.extend(shape)
1155 else:
1156 shape_tuple = list(shape)
1157
1158 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1159
1160 return shape_tuple, nTxs
1161
1162 def run(self, dataOut, shape=None, nTxs=None):
1163
1164 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1165
1166 dataOut.flagNoData = True
1167 profileIndex = None
1168
1169 if dataOut.flagDataAsBlock:
1170
1171 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1172 dataOut.flagNoData = False
1173
1174 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1175
1176 else:
1177
1178 if self.__nTxs < 1:
1179
1180 self.__appendProfile(dataOut, self.__nTxs)
1181 new_data = self.__getBuffer()
1182
1183 if new_data is not None:
1184 dataOut.data = new_data
1185 dataOut.flagNoData = False
1186
1187 profileIndex = dataOut.profileIndex*nTxs
1188
1189 else:
1190 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1191
1192 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1193
1194 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1195
1196 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1197
1198 dataOut.profileIndex = profileIndex
1199
1200 dataOut.ippSeconds /= self.__nTxs
1201
1202 return dataOut
1203
1204 class SplitProfiles(Operation):
1205
1206 def __init__(self, **kwargs):
1207
1208 Operation.__init__(self, **kwargs)
1209
1210 def run(self, dataOut, n):
1211
1212 dataOut.flagNoData = True
1213 profileIndex = None
1214
1215 if dataOut.flagDataAsBlock:
1216
1217 #nchannels, nprofiles, nsamples
1218 shape = dataOut.data.shape
1219
1220 if shape[2] % n != 0:
1221 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1222
1223 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1224
1225 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1226 dataOut.flagNoData = False
1227
1228 profileIndex = int(dataOut.nProfiles/n) - 1
1229
1230 else:
1231
1232 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1233
1234 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1235
1236 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1237
1238 dataOut.nProfiles = int(dataOut.nProfiles*n)
1239
1240 dataOut.profileIndex = profileIndex
1241
1242 dataOut.ippSeconds /= n
1243
1244 return dataOut
1245
1246 class CombineProfiles(Operation):
1247 def __init__(self, **kwargs):
1248
1249 Operation.__init__(self, **kwargs)
1250
1251 self.__remData = None
1252 self.__profileIndex = 0
1253
1254 def run(self, dataOut, n):
1255
1256 dataOut.flagNoData = True
1257 profileIndex = None
1258
1259 if dataOut.flagDataAsBlock:
1260
1261 #nchannels, nprofiles, nsamples
1262 shape = dataOut.data.shape
1263 new_shape = shape[0], shape[1]/n, shape[2]*n
1264
1265 if shape[1] % n != 0:
1266 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1267
1268 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1269 dataOut.flagNoData = False
1270
1271 profileIndex = int(dataOut.nProfiles*n) - 1
1272
1273 else:
1274
1275 #nchannels, nsamples
1276 if self.__remData is None:
1277 newData = dataOut.data
1278 else:
1279 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1280
1281 self.__profileIndex += 1
1282
1283 if self.__profileIndex < n:
1284 self.__remData = newData
1285 #continue
1286 return
1287
1288 self.__profileIndex = 0
1289 self.__remData = None
1290
1291 dataOut.data = newData
1292 dataOut.flagNoData = False
1293
1294 profileIndex = dataOut.profileIndex/n
1295
1296
1297 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1298
1299 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1300
1301 dataOut.nProfiles = int(dataOut.nProfiles/n)
1302
1303 dataOut.profileIndex = profileIndex
1304
1305 dataOut.ippSeconds *= n
1306
1307 return dataOut
1308
1309 class PulsePairVoltage(Operation):
1310 '''
1311 Function PulsePair(Signal Power, Velocity)
1312 The real component of Lag[0] provides Intensity Information
1313 The imag component of Lag[1] Phase provides Velocity Information
1314
1315 Configuration Parameters:
1316 nPRF = Number of Several PRF
1317 theta = Degree Azimuth angel Boundaries
1318
1319 Input:
1320 self.dataOut
1321 lag[N]
1322 Affected:
1323 self.dataOut.spc
1324 '''
1325 isConfig = False
1326 __profIndex = 0
1327 __initime = None
1328 __lastdatatime = None
1329 __buffer = None
1330 noise = None
1331 __dataReady = False
1332 n = None
1333 __nch = 0
1334 __nHeis = 0
1335 removeDC = False
1336 ipp = None
1337 lambda_ = 0
1338
1339 def __init__(self,**kwargs):
1340 Operation.__init__(self,**kwargs)
1341
1342 def setup(self, dataOut, n = None, removeDC=False):
1343 '''
1344 n= Numero de PRF's de entrada
1345 '''
1346 self.__initime = None
1347 self.__lastdatatime = 0
1348 self.__dataReady = False
1349 self.__buffer = 0
1350 self.__profIndex = 0
1351 self.noise = None
1352 self.__nch = dataOut.nChannels
1353 self.__nHeis = dataOut.nHeights
1354 self.removeDC = removeDC
1355 self.lambda_ = 3.0e8/(9345.0e6)
1356 self.ippSec = dataOut.ippSeconds
1357 self.nCohInt = dataOut.nCohInt
1358 print("IPPseconds",dataOut.ippSeconds)
1359
1360 print("ELVALOR DE n es:", n)
1361 if n == None:
1362 raise ValueError("n should be specified.")
1363
1364 if n != None:
1365 if n<2:
1366 raise ValueError("n should be greater than 2")
1367
1368 self.n = n
1369 self.__nProf = n
1370
1371 self.__buffer = numpy.zeros((dataOut.nChannels,
1372 n,
1373 dataOut.nHeights),
1374 dtype='complex')
1375
1376 def putData(self,data):
1377 '''
1378 Add a profile to he __buffer and increase in one the __profiel Index
1379 '''
1380 self.__buffer[:,self.__profIndex,:]= data
1381 self.__profIndex += 1
1382 return
1383
1384 def pushData(self,dataOut):
1385 '''
1386 Return the PULSEPAIR and the profiles used in the operation
1387 Affected : self.__profileIndex
1388 '''
1389 #----------------- Remove DC-----------------------------------
1390 if self.removeDC==True:
1391 mean = numpy.mean(self.__buffer,1)
1392 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1393 dc= numpy.tile(tmp,[1,self.__nProf,1])
1394 self.__buffer = self.__buffer - dc
1395 #------------------Calculo de Potencia ------------------------
1396 pair0 = self.__buffer*numpy.conj(self.__buffer)
1397 pair0 = pair0.real
1398 lag_0 = numpy.sum(pair0,1)
1399 #------------------Calculo de Ruido x canal--------------------
1400 self.noise = numpy.zeros(self.__nch)
1401 for i in range(self.__nch):
1402 daux = numpy.sort(pair0[i,:,:],axis= None)
1403 self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt)
1404
1405 self.noise = self.noise.reshape(self.__nch,1)
1406 self.noise = numpy.tile(self.noise,[1,self.__nHeis])
1407 noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis)
1408 noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1])
1409 #------------------ Potencia recibida= P , Potencia senal = S , Ruido= N--
1410 #------------------ P= S+N ,P=lag_0/N ---------------------------------
1411 #-------------------- Power --------------------------------------------------
1412 data_power = lag_0/(self.n*self.nCohInt)
1413 #------------------ Senal ---------------------------------------------------
1414 data_intensity = pair0 - noise_buffer
1415 data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt)
1416 #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt)
1417 for i in range(self.__nch):
1418 for j in range(self.__nHeis):
1419 if data_intensity[i][j] < 0:
1420 data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j]))
1421
1422 #----------------- Calculo de Frecuencia y Velocidad doppler--------
1423 pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
1424 lag_1 = numpy.sum(pair1,1)
1425 data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
1426 data_velocity = (self.lambda_/2.0)*data_freq
1427
1428 #---------------- Potencia promedio estimada de la Senal-----------
1429 lag_0 = lag_0/self.n
1430 S = lag_0-self.noise
1431
1432 #---------------- Frecuencia Doppler promedio ---------------------
1433 lag_1 = lag_1/(self.n-1)
1434 R1 = numpy.abs(lag_1)
1435
1436 #---------------- Calculo del SNR----------------------------------
1437 data_snrPP = S/self.noise
1438 for i in range(self.__nch):
1439 for j in range(self.__nHeis):
1440 if data_snrPP[i][j] < 1.e-20:
1441 data_snrPP[i][j] = 1.e-20
1442
1443 #----------------- Calculo del ancho espectral ----------------------
1444 L = S/R1
1445 L = numpy.where(L<0,1,L)
1446 L = numpy.log(L)
1447 tmp = numpy.sqrt(numpy.absolute(L))
1448 data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L)
1449 n = self.__profIndex
1450
1451 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1452 self.__profIndex = 0
1453 return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,n
1454
1455
1456 def pulsePairbyProfiles(self,dataOut):
1457
1458 self.__dataReady = False
1459 data_power = None
1460 data_intensity = None
1461 data_velocity = None
1462 data_specwidth = None
1463 data_snrPP = None
1464 self.putData(data=dataOut.data)
1465 if self.__profIndex == self.n:
1466 data_power,data_intensity, data_velocity,data_snrPP,data_specwidth, n = self.pushData(dataOut=dataOut)
1467 self.__dataReady = True
1468
1469 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth
1470
1471
1472 def pulsePairOp(self, dataOut, datatime= None):
1473
1474 if self.__initime == None:
1475 self.__initime = datatime
1476 data_power, data_intensity, data_velocity, data_snrPP, data_specwidth = self.pulsePairbyProfiles(dataOut)
1477 self.__lastdatatime = datatime
1478
1479 if data_power is None:
1480 return None, None, None,None,None,None
1481
1482 avgdatatime = self.__initime
1483 deltatime = datatime - self.__lastdatatime
1484 self.__initime = datatime
1485
1486 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth, avgdatatime
1487
1488 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1489
1490 if not self.isConfig:
1491 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1492 self.isConfig = True
1493 data_power, data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
1494 dataOut.flagNoData = True
1495
1496 if self.__dataReady:
1497 dataOut.nCohInt *= self.n
1498 dataOut.dataPP_POW = data_intensity # S
1499 dataOut.dataPP_POWER = data_power # P
1500 dataOut.dataPP_DOP = data_velocity
1501 dataOut.dataPP_SNR = data_snrPP
1502 dataOut.dataPP_WIDTH = data_specwidth
1503 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1504 dataOut.utctime = avgdatatime
1505 dataOut.flagNoData = False
1506 return dataOut
1507
1508
1509
1510 # import collections
1511 # from scipy.stats import mode
1512 #
1513 # class Synchronize(Operation):
1514 #
1515 # isConfig = False
1516 # __profIndex = 0
1517 #
1518 # def __init__(self, **kwargs):
1519 #
1520 # Operation.__init__(self, **kwargs)
1521 # # self.isConfig = False
1522 # self.__powBuffer = None
1523 # self.__startIndex = 0
1524 # self.__pulseFound = False
1525 #
1526 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1527 #
1528 # #Read data
1529 #
1530 # powerdB = dataOut.getPower(channel = channel)
1531 # noisedB = dataOut.getNoise(channel = channel)[0]
1532 #
1533 # self.__powBuffer.extend(powerdB.flatten())
1534 #
1535 # dataArray = numpy.array(self.__powBuffer)
1536 #
1537 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1538 #
1539 # maxValue = numpy.nanmax(filteredPower)
1540 #
1541 # if maxValue < noisedB + 10:
1542 # #No se encuentra ningun pulso de transmision
1543 # return None
1544 #
1545 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1546 #
1547 # if len(maxValuesIndex) < 2:
1548 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1549 # return None
1550 #
1551 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1552 #
1553 # #Seleccionar solo valores con un espaciamiento de nSamples
1554 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1555 #
1556 # if len(pulseIndex) < 2:
1557 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1558 # return None
1559 #
1560 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1561 #
1562 # #remover senales que se distancien menos de 10 unidades o muestras
1563 # #(No deberian existir IPP menor a 10 unidades)
1564 #
1565 # realIndex = numpy.where(spacing > 10 )[0]
1566 #
1567 # if len(realIndex) < 2:
1568 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1569 # return None
1570 #
1571 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1572 # realPulseIndex = pulseIndex[realIndex]
1573 #
1574 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1575 #
1576 # print "IPP = %d samples" %period
1577 #
1578 # self.__newNSamples = dataOut.nHeights #int(period)
1579 # self.__startIndex = int(realPulseIndex[0])
1580 #
1581 # return 1
1582 #
1583 #
1584 # def setup(self, nSamples, nChannels, buffer_size = 4):
1585 #
1586 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1587 # maxlen = buffer_size*nSamples)
1588 #
1589 # bufferList = []
1590 #
1591 # for i in range(nChannels):
1592 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1593 # maxlen = buffer_size*nSamples)
1594 #
1595 # bufferList.append(bufferByChannel)
1596 #
1597 # self.__nSamples = nSamples
1598 # self.__nChannels = nChannels
1599 # self.__bufferList = bufferList
1600 #
1601 # def run(self, dataOut, channel = 0):
1602 #
1603 # if not self.isConfig:
1604 # nSamples = dataOut.nHeights
1605 # nChannels = dataOut.nChannels
1606 # self.setup(nSamples, nChannels)
1607 # self.isConfig = True
1608 #
1609 # #Append new data to internal buffer
1610 # for thisChannel in range(self.__nChannels):
1611 # bufferByChannel = self.__bufferList[thisChannel]
1612 # bufferByChannel.extend(dataOut.data[thisChannel])
1613 #
1614 # if self.__pulseFound:
1615 # self.__startIndex -= self.__nSamples
1616 #
1617 # #Finding Tx Pulse
1618 # if not self.__pulseFound:
1619 # indexFound = self.__findTxPulse(dataOut, channel)
1620 #
1621 # if indexFound == None:
1622 # dataOut.flagNoData = True
1623 # return
1624 #
1625 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1626 # self.__pulseFound = True
1627 # self.__startIndex = indexFound
1628 #
1629 # #If pulse was found ...
1630 # for thisChannel in range(self.__nChannels):
1631 # bufferByChannel = self.__bufferList[thisChannel]
1632 # #print self.__startIndex
1633 # x = numpy.array(bufferByChannel)
1634 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1635 #
1636 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1637 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1638 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1639 #
1640 # dataOut.data = self.__arrayBuffer
1641 #
1642 # self.__startIndex += self.__newNSamples
1643 #
1644 # return No newline at end of file
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (833 lines changed) Show them Hide them
@@ -0,0 +1,833
1 class SpectralFitting(Operation):
2 '''
3 Function GetMoments()
4
5 Input:
6 Output:
7 Variables modified:
8 '''
9 isConfig = False
10 __dataReady = False
11 bloques = None
12 bloque0 = None
13 index = 0
14 fint = 0
15 buffer = 0
16 buffer2 = 0
17 buffer3 = 0
18
19 def __init__(self):
20 Operation.__init__(self)
21 self.i=0
22 self.isConfig = False
23
24
25 def setup(self,nChan,nProf,nHei,nBlocks):
26 self.__dataReady = False
27 self.bloques = numpy.zeros([2, nProf, nHei,nBlocks], dtype= complex)
28 self.bloque0 = numpy.zeros([nChan, nProf, nHei, nBlocks])
29
30 def __calculateMoments(self,oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
31
32 if (nicoh is None): nicoh = 1
33 if (graph is None): graph = 0
34 if (smooth is None): smooth = 0
35 elif (self.smooth < 3): smooth = 0
36
37 if (type1 is None): type1 = 0
38 if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
39 if (snrth is None): snrth = -3
40 if (dc is None): dc = 0
41 if (aliasing is None): aliasing = 0
42 if (oldfd is None): oldfd = 0
43 if (wwauto is None): wwauto = 0
44
45 if (n0 < 1.e-20): n0 = 1.e-20
46
47 freq = oldfreq
48 vec_power = numpy.zeros(oldspec.shape[1])
49 vec_fd = numpy.zeros(oldspec.shape[1])
50 vec_w = numpy.zeros(oldspec.shape[1])
51 vec_snr = numpy.zeros(oldspec.shape[1])
52
53 oldspec = numpy.ma.masked_invalid(oldspec)
54
55 for ind in range(oldspec.shape[1]):
56
57 spec = oldspec[:,ind]
58 aux = spec*fwindow
59 max_spec = aux.max()
60 m = list(aux).index(max_spec)
61
62 #Smooth
63 if (smooth == 0): spec2 = spec
64 else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
65
66 # Calculo de Momentos
67 bb = spec2[list(range(m,spec2.size))]
68 bb = (bb<n0).nonzero()
69 bb = bb[0]
70
71 ss = spec2[list(range(0,m + 1))]
72 ss = (ss<n0).nonzero()
73 ss = ss[0]
74
75 if (bb.size == 0):
76 bb0 = spec.size - 1 - m
77 else:
78 bb0 = bb[0] - 1
79 if (bb0 < 0):
80 bb0 = 0
81
82 if (ss.size == 0): ss1 = 1
83 else: ss1 = max(ss) + 1
84
85 if (ss1 > m): ss1 = m
86
87 valid = numpy.asarray(list(range(int(m + bb0 - ss1 + 1)))) + ss1
88 power = ((spec2[valid] - n0)*fwindow[valid]).sum()
89 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
90 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
91 snr = (spec2.mean()-n0)/n0
92
93 if (snr < 1.e-20) :
94 snr = 1.e-20
95
96 vec_power[ind] = power
97 vec_fd[ind] = fd
98 vec_w[ind] = w
99 vec_snr[ind] = snr
100
101 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
102 return moments
103
104 def __DiffCoherent(self, spectra, cspectra, dataOut, noise, snrth, coh_th, hei_th):
105
106 nProf = dataOut.nProfiles
107 heights = dataOut.heightList
108 nHei = len(heights)
109 channels = dataOut.channelList
110 nChan = len(channels)
111 crosspairs = dataOut.groupList
112 nPairs = len(crosspairs)
113 #Separar espectros incoherentes de coherentes snr > 20 dB'
114 snr_th = 10**(snrth/10.0)
115 my_incoh_spectra = numpy.zeros([nChan, nProf,nHei], dtype='float')
116 my_incoh_cspectra = numpy.zeros([nPairs,nProf, nHei], dtype='complex')
117 my_incoh_aver = numpy.zeros([nChan, nHei])
118 my_coh_aver = numpy.zeros([nChan, nHei])
119
120 coh_spectra = numpy.zeros([nChan, nProf, nHei], dtype='float')
121 coh_cspectra = numpy.zeros([nPairs, nProf, nHei], dtype='complex')
122 coh_aver = numpy.zeros([nChan, nHei])
123
124 incoh_spectra = numpy.zeros([nChan, nProf, nHei], dtype='float')
125 incoh_cspectra = numpy.zeros([nPairs, nProf, nHei], dtype='complex')
126 incoh_aver = numpy.zeros([nChan, nHei])
127 power = numpy.sum(spectra, axis=1)
128
129 if coh_th == None : coh_th = numpy.array([0.75,0.65,0.15]) # 0.65
130 if hei_th == None : hei_th = numpy.array([60,300,650])
131 for ic in range(nPairs):
132 pair = crosspairs[ic]
133 #si el SNR es mayor que el SNR threshold los datos se toman coherentes
134 s_n0 = power[pair[0],:]/noise[pair[0]]
135 s_n1 = power[pair[1],:]/noise[pair[1]]
136 valid1 =(s_n0>=snr_th).nonzero()
137 valid2 = (s_n1>=snr_th).nonzero()
138 valid1 = numpy.array(valid1[0])
139 valid2 = numpy.array(valid2[0])
140 valid = valid1
141 for iv in range(len(valid2)):
142 indv = numpy.array((valid1 == valid2[iv]).nonzero())
143 if len(indv[0]) == 0 :
144 valid = numpy.concatenate((valid,valid2[iv]), axis=None)
145 if len(valid)>0:
146 my_coh_aver[pair[0],valid]=1
147 my_coh_aver[pair[1],valid]=1
148 # si la coherencia es mayor a la coherencia threshold los datos se toman
149 coh = numpy.squeeze(numpy.nansum(cspectra[ic,:,:], axis=0)/numpy.sqrt(numpy.nansum(spectra[pair[0],:,:], axis=0)*numpy.nansum(spectra[pair[1],:,:], axis=0)))
150 for ih in range(len(hei_th)):
151 hvalid = (heights>hei_th[ih]).nonzero()
152 hvalid = hvalid[0]
153 if len(hvalid)>0:
154 valid = (numpy.absolute(coh[hvalid])>coh_th[ih]).nonzero()
155 valid = valid[0]
156 if len(valid)>0:
157 my_coh_aver[pair[0],hvalid[valid]] =1
158 my_coh_aver[pair[1],hvalid[valid]] =1
159
160 coh_echoes = (my_coh_aver[pair[0],:] == 1).nonzero()
161 incoh_echoes = (my_coh_aver[pair[0],:] != 1).nonzero()
162 incoh_echoes = incoh_echoes[0]
163 if len(incoh_echoes) > 0:
164 my_incoh_spectra[pair[0],:,incoh_echoes] = spectra[pair[0],:,incoh_echoes]
165 my_incoh_spectra[pair[1],:,incoh_echoes] = spectra[pair[1],:,incoh_echoes]
166 my_incoh_cspectra[ic,:,incoh_echoes] = cspectra[ic,:,incoh_echoes]
167 my_incoh_aver[pair[0],incoh_echoes] = 1
168 my_incoh_aver[pair[1],incoh_echoes] = 1
169
170
171 for ic in range(nPairs):
172 pair = crosspairs[ic]
173
174 valid1 =(my_coh_aver[pair[0],:]==1 ).nonzero()
175 valid2 = (my_coh_aver[pair[1],:]==1).nonzero()
176 valid1 = numpy.array(valid1[0])
177 valid2 = numpy.array(valid2[0])
178 valid = valid1
179
180 for iv in range(len(valid2)):
181
182 indv = numpy.array((valid1 == valid2[iv]).nonzero())
183 if len(indv[0]) == 0 :
184 valid = numpy.concatenate((valid,valid2[iv]), axis=None)
185 valid1 =(my_coh_aver[pair[0],:] !=1 ).nonzero()
186 valid2 = (my_coh_aver[pair[1],:] !=1).nonzero()
187 valid1 = numpy.array(valid1[0])
188 valid2 = numpy.array(valid2[0])
189 incoh_echoes = valid1
190 for iv in range(len(valid2)):
191
192 indv = numpy.array((valid1 == valid2[iv]).nonzero())
193 if len(indv[0]) == 0 :
194 incoh_echoes = numpy.concatenate(( incoh_echoes,valid2[iv]), axis=None)
195
196 if len(valid)>0:
197 coh_spectra[pair[0],:,valid] = spectra[pair[0],:,valid]
198 coh_spectra[pair[1],:,valid] = spectra[pair[1],:,valid]
199 coh_cspectra[ic,:,valid] = cspectra[ic,:,valid]
200 coh_aver[pair[0],valid]=1
201 coh_aver[pair[1],valid]=1
202 if len(incoh_echoes)>0:
203 incoh_spectra[pair[0],:,incoh_echoes] = spectra[pair[0],:,incoh_echoes]
204 incoh_spectra[pair[1],:,incoh_echoes] = spectra[pair[1],:,incoh_echoes]
205 incoh_cspectra[ic,:,incoh_echoes] = cspectra[ic,:,incoh_echoes]
206 incoh_aver[pair[0],incoh_echoes]=1
207 incoh_aver[pair[1],incoh_echoes]=1
208 return my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver
209
210
211 def __CleanCoherent(self,snrth, spectra, cspectra, coh_aver,dataOut, noise,clean_coh_echoes,index):
212
213 nProf = dataOut.nProfiles
214 heights = dataOut.heightList
215 nHei = len(heights)
216 channels = dataOut.channelList
217 nChan = len(channels)
218 crosspairs = dataOut.groupList
219 nPairs = len(crosspairs)
220
221 absc = dataOut.abscissaList[:-1]
222 data_param = numpy.zeros((nChan, 4, spectra.shape[2]))
223 clean_coh_spectra = spectra.copy()
224 clean_coh_cspectra = cspectra.copy()
225 clean_coh_aver = coh_aver.copy()
226
227 spwd_th=[10,6] #spwd_th[0] --> For satellites ; spwd_th[1] --> For special events like SUN.
228 coh_th = 0.75
229
230 rtime0 = [6,18] # periodo sin ESF
231 rtime1 = [10.5,13.5] # periodo con alta coherencia y alto ancho espectral (esperado): SOL.
232
233 time = index*5./60 # en base a 5 min de proceso
234 if clean_coh_echoes == 1 :
235 for ind in range(nChan):
236 data_param[ind,:,:] = self.__calculateMoments( spectra[ind,:,:] , absc , noise[ind] )
237 spwd = data_param[:,3]
238 # SPECB_JULIA,header=anal_header,jspectra=spectra,vel=velocities,hei=heights, num_aver=1, mode_fit=0,smoothing=smoothing,jvelr=velr,jspwd=spwd,jsnr=snr,jnoise=noise,jstdvnoise=stdvnoise
239 # para obtener spwd
240 for ic in range(nPairs):
241 pair = crosspairs[ic]
242 coh = numpy.squeeze(numpy.sum(cspectra[ic,:,:], axis=1)/numpy.sqrt(numpy.sum(spectra[pair[0],:,:], axis=1)*numpy.sum(spectra[pair[1],:,:], axis=1)))
243 for ih in range(nHei) :
244 # Considering heights higher than 200km in order to avoid removing phenomena like EEJ.
245 if heights[ih] >= 200 and coh_aver[pair[0],ih] == 1 and coh_aver[pair[1],ih] == 1 :
246 # Checking coherence
247 if (numpy.abs(coh[ih]) <= coh_th) or (time >= rtime0[0] and time <= rtime0[1]) :
248 # Checking spectral widths
249 if (spwd[pair[0],ih] > spwd_th[0]) or (spwd[pair[1],ih] > spwd_th[0]) :
250 # satelite
251 clean_coh_spectra[pair,ih,:] = 0.0
252 clean_coh_cspectra[ic,ih,:] = 0.0
253 clean_coh_aver[pair,ih] = 0
254 else :
255 if ((spwd[pair[0],ih] < spwd_th[1]) or (spwd[pair[1],ih] < spwd_th[1])) :
256 # Especial event like sun.
257 clean_coh_spectra[pair,ih,:] = 0.0
258 clean_coh_cspectra[ic,ih,:] = 0.0
259 clean_coh_aver[pair,ih] = 0
260
261 return clean_coh_spectra, clean_coh_cspectra, clean_coh_aver
262
263 def CleanRayleigh(self,dataOut,spectra,cspectra,save_drifts):
264
265 rfunc = cspectra.copy()
266 n_funct = len(rfunc[0,:,0,0])
267 val_spc = spectra*0.0
268 val_cspc = cspectra*0.0
269 in_sat_spectra = spectra.copy()
270 in_sat_cspectra = cspectra.copy()
271
272 min_hei = 200
273 nProf = dataOut.nProfiles
274 heights = dataOut.heightList
275 nHei = len(heights)
276 channels = dataOut.channelList
277 nChan = len(channels)
278 crosspairs = dataOut.groupList
279 nPairs = len(crosspairs)
280 hval=(heights >= min_hei).nonzero()
281 ih=hval[0]
282 for ih in range(hval[0][0],nHei):
283 for ifreq in range(nProf):
284 for ii in range(n_funct):
285
286 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih]))
287 val = (numpy.isfinite(func2clean)==True).nonzero()
288 if len(val)>0:
289 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
290 if min_val <= -40 : min_val = -40
291 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
292 if max_val >= 200 : max_val = 200
293 step = 1
294 #Getting bins and the histogram
295 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
296 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
297 mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
298 sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
299 parg = [numpy.amax(y_dist),mean,sigma]
300 try :
301 gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
302 mode = gauss_fit[1]
303 stdv = gauss_fit[2]
304 except:
305 mode = mean
306 stdv = sigma
307
308 #Removing echoes greater than mode + 3*stdv
309 factor_stdv = 2.5
310 noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
311
312 if len(noval[0]) > 0:
313 novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
314 cross_pairs = crosspairs[ii]
315 #Getting coherent echoes which are removed.
316 if len(novall[0]) > 0:
317 val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
318 val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
319 val_cspc[novall[0],ii,ifreq,ih] = 1
320 #Removing coherent from ISR data
321 spectra[noval,cross_pairs[0],ifreq,ih] = numpy.nan
322 spectra[noval,cross_pairs[1],ifreq,ih] = numpy.nan
323 cspectra[noval,ii,ifreq,ih] = numpy.nan
324
325 #Getting average of the spectra and cross-spectra from incoherent echoes.
326 out_spectra = numpy.zeros([nChan,nProf,nHei], dtype=float) #+numpy.nan
327 out_cspectra = numpy.zeros([nPairs,nProf,nHei], dtype=complex) #+numpy.nan
328 for ih in range(nHei):
329 for ifreq in range(nProf):
330 for ich in range(nChan):
331 tmp = spectra[:,ich,ifreq,ih]
332 valid = (numpy.isfinite(tmp[:])==True).nonzero()
333 if len(valid[0]) >0 :
334 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
335 for icr in range(nPairs):
336 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
337 valid = (numpy.isfinite(tmp)==True).nonzero()
338 if len(valid[0]) > 0:
339 out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
340 #Removing fake coherent echoes (at least 4 points around the point)
341 val_spectra = numpy.sum(val_spc,0)
342 val_cspectra = numpy.sum(val_cspc,0)
343
344 val_spectra = self.REM_ISOLATED_POINTS(val_spectra,4)
345 val_cspectra = self.REM_ISOLATED_POINTS(val_cspectra,4)
346
347 for i in range(nChan):
348 for j in range(nProf):
349 for k in range(nHei):
350 if numpy.isfinite(val_spectra[i,j,k]) and val_spectra[i,j,k] < 1 :
351 val_spc[:,i,j,k] = 0.0
352 for i in range(nPairs):
353 for j in range(nProf):
354 for k in range(nHei):
355 if numpy.isfinite(val_cspectra[i,j,k]) and val_cspectra[i,j,k] < 1 :
356 val_cspc[:,i,j,k] = 0.0
357
358 tmp_sat_spectra = spectra.copy()
359 tmp_sat_spectra = tmp_sat_spectra*numpy.nan
360 tmp_sat_cspectra = cspectra.copy()
361 tmp_sat_cspectra = tmp_sat_cspectra*numpy.nan
362 val = (val_spc > 0).nonzero()
363 if len(val[0]) > 0:
364 tmp_sat_spectra[val] = in_sat_spectra[val]
365
366 val = (val_cspc > 0).nonzero()
367 if len(val[0]) > 0:
368 tmp_sat_cspectra[val] = in_sat_cspectra[val]
369
370 #Getting average of the spectra and cross-spectra from incoherent echoes.
371 sat_spectra = numpy.zeros((nChan,nProf,nHei), dtype=float)
372 sat_cspectra = numpy.zeros((nPairs,nProf,nHei), dtype=complex)
373 for ih in range(nHei):
374 for ifreq in range(nProf):
375 for ich in range(nChan):
376 tmp = numpy.squeeze(tmp_sat_spectra[:,ich,ifreq,ih])
377 valid = (numpy.isfinite(tmp)).nonzero()
378 if len(valid[0]) > 0:
379 sat_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
380
381 for icr in range(nPairs):
382 tmp = numpy.squeeze(tmp_sat_cspectra[:,icr,ifreq,ih])
383 valid = (numpy.isfinite(tmp)).nonzero()
384 if len(valid[0]) > 0:
385 sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
386 return out_spectra, out_cspectra,sat_spectra,sat_cspectra
387 def REM_ISOLATED_POINTS(self,array,rth):
388 if rth == None : rth = 4
389 num_prof = len(array[0,:,0])
390 num_hei = len(array[0,0,:])
391 n2d = len(array[:,0,0])
392
393 for ii in range(n2d) :
394 tmp = array[ii,:,:]
395 tmp = numpy.reshape(tmp,num_prof*num_hei)
396 indxs1 = (numpy.isfinite(tmp)==True).nonzero()
397 indxs2 = (tmp > 0).nonzero()
398 indxs1 = (indxs1[0])
399 indxs2 = indxs2[0]
400 indxs = None
401 for iv in range(len(indxs2)):
402 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
403 if len(indv[0]) > 0 :
404 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
405 indxs = indxs[1:]
406 if len(indxs) < 4 :
407 array[ii,:,:] = 0.
408 return
409
410 xpos = numpy.mod(indxs ,num_hei)
411 ypos = (indxs / num_hei)
412 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
413 xpos = xpos[sx]
414 ypos = ypos[sx]
415 # *********************************** Cleaning isolated points **********************************
416 ic = 0
417 while True :
418 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
419 no_coh1 = (numpy.isfinite(r)==True).nonzero()
420 no_coh2 = (r <= rth).nonzero()
421 no_coh1 = numpy.array(no_coh1[0])
422 no_coh2 = numpy.array(no_coh2[0])
423 no_coh = None
424 for iv in range(len(no_coh2)):
425 indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
426 if len(indv[0]) > 0 :
427 no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
428 no_coh = no_coh[1:]
429 if len(no_coh) < 4 :
430 xpos[ic] = numpy.nan
431 ypos[ic] = numpy.nan
432
433 ic = ic + 1
434 if (ic == len(indxs)) :
435 break
436 indxs = (numpy.isfinite(list(xpos))==True).nonzero()
437 if len(indxs[0]) < 4 :
438 array[ii,:,:] = 0.
439 return
440
441 xpos = xpos[indxs[0]]
442 ypos = ypos[indxs[0]]
443 for i in range(0,len(ypos)):
444 ypos[i]=int(ypos[i])
445 junk = tmp
446 tmp = junk*0.0
447
448 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
449 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
450 return array
451
452 def moments(self,doppler,yarray,npoints):
453 ytemp = yarray
454 val = (ytemp > 0).nonzero()
455 val = val[0]
456 if len(val) == 0 : val = range(npoints-1)
457
458 ynew = 0.5*(ytemp[val[0]]+ytemp[val[len(val)-1]])
459 ytemp[len(ytemp):] = [ynew]
460
461 index = 0
462 index = numpy.argmax(ytemp)
463 ytemp = numpy.roll(ytemp,int(npoints/2)-1-index)
464 ytemp = ytemp[0:npoints-1]
465
466 fmom = numpy.sum(doppler*ytemp)/numpy.sum(ytemp)+(index-(npoints/2-1))*numpy.abs(doppler[1]-doppler[0])
467 smom = numpy.sum(doppler*doppler*ytemp)/numpy.sum(ytemp)
468 return [fmom,numpy.sqrt(smom)]
469
470
471
472
473
474 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None, filec=None,coh_th=None, hei_th=None,taver=None,proc=None,nhei=None,nprofs=None,ipp=None,channelList=None):
475 if not numpy.any(proc):
476 nChannels = dataOut.nChannels
477 nHeights= dataOut.heightList.size
478 nProf = dataOut.nProfiles
479 if numpy.any(taver): taver=int(taver)
480 else : taver = 5
481 tini=time.localtime(dataOut.utctime)
482 if (tini.tm_min % taver) == 0 and (tini.tm_sec < 5 and self.fint==0):
483 self.index = 0
484 jspc = self.buffer
485 jcspc = self.buffer2
486 jnoise = self.buffer3
487 self.buffer = dataOut.data_spc
488 self.buffer2 = dataOut.data_cspc
489 self.buffer3 = dataOut.noise
490 self.fint = 1
491 if numpy.any(jspc) :
492 jspc= numpy.reshape(jspc,(int(len(jspc)/nChannels),nChannels,nProf,nHeights))
493 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/int(nChannels/2)),int(nChannels/2),nProf,nHeights))
494 jnoise= numpy.reshape(jnoise,(int(len(jnoise)/nChannels),nChannels))
495 else:
496 dataOut.flagNoData = True
497 return dataOut
498 else :
499 if (tini.tm_min % taver) == 0 : self.fint = 1
500 else : self.fint = 0
501 self.index += 1
502 if numpy.any(self.buffer):
503 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
504 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
505 self.buffer3 = numpy.concatenate((self.buffer3,dataOut.noise), axis=0)
506 else:
507 self.buffer = dataOut.data_spc
508 self.buffer2 = dataOut.data_cspc
509 self.buffer3 = dataOut.noise
510 dataOut.flagNoData = True
511 return dataOut
512 if path != None:
513 sys.path.append(path)
514 self.library = importlib.import_module(file)
515 if filec != None:
516 self.weightf = importlib.import_module(filec)
517
518 #To be inserted as a parameter
519 groupArray = numpy.array(groupList)
520 #groupArray = numpy.array([[0,1],[2,3]])
521 dataOut.groupList = groupArray
522 nGroups = groupArray.shape[0]
523 nChannels = dataOut.nChannels
524 nHeights = dataOut.heightList.size
525
526 #Parameters Array
527 dataOut.data_param = None
528 dataOut.data_paramC = None
529 dataOut.clean_num_aver = None
530 dataOut.coh_num_aver = None
531 dataOut.tmp_spectra_i = None
532 dataOut.tmp_cspectra_i = None
533 dataOut.tmp_spectra_c = None
534 dataOut.tmp_cspectra_c = None
535 dataOut.index = None
536
537 #Set constants
538 constants = self.library.setConstants(dataOut)
539 dataOut.constants = constants
540 M = dataOut.normFactor
541 N = dataOut.nFFTPoints
542 ippSeconds = dataOut.ippSeconds
543 K = dataOut.nIncohInt
544 pairsArray = numpy.array(dataOut.pairsList)
545 snrth= 20
546 spectra = dataOut.data_spc
547 cspectra = dataOut.data_cspc
548 nProf = dataOut.nProfiles
549 heights = dataOut.heightList
550 nHei = len(heights)
551 channels = dataOut.channelList
552 nChan = len(channels)
553 nIncohInt = dataOut.nIncohInt
554 crosspairs = dataOut.groupList
555 noise = dataOut.noise
556 jnoise = jnoise/N
557 noise = numpy.nansum(jnoise,axis=0)#/len(jnoise)
558 power = numpy.sum(spectra, axis=1)
559 nPairs = len(crosspairs)
560 absc = dataOut.abscissaList[:-1]
561
562 if not self.isConfig:
563 self.isConfig = True
564
565 index = tini.tm_hour*12+tini.tm_min/taver
566 dataOut.index= index
567 jspc = jspc/N/N
568 jcspc = jcspc/N/N
569 tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.CleanRayleigh(dataOut,jspc,jcspc,2)
570 jspectra = tmp_spectra*len(jspc[:,0,0,0])
571 jcspectra = tmp_cspectra*len(jspc[:,0,0,0])
572 my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver = self.__DiffCoherent(jspectra, jcspectra, dataOut, noise, snrth,coh_th, hei_th)
573 clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index)
574 dataOut.data_spc = incoh_spectra
575 dataOut.data_cspc = incoh_cspectra
576 clean_num_aver = incoh_aver*len(jspc[:,0,0,0])
577 coh_num_aver = clean_coh_aver*len(jspc[:,0,0,0])
578 dataOut.clean_num_aver = clean_num_aver
579 dataOut.coh_num_aver = coh_num_aver
580 dataOut.tmp_spectra_i = incoh_spectra
581 dataOut.tmp_cspectra_i = incoh_cspectra
582 dataOut.tmp_spectra_c = clean_coh_spectra
583 dataOut.tmp_cspectra_c = clean_coh_cspectra
584 #List of possible combinations
585 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
586 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
587 if getSNR:
588 listChannels = groupArray.reshape((groupArray.size))
589 listChannels.sort()
590 dataOut.data_SNR = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise[listChannels])
591 else:
592 clean_num_aver = dataOut.clean_num_aver
593 coh_num_aver = dataOut.coh_num_aver
594 dataOut.data_spc = dataOut.tmp_spectra_i
595 dataOut.data_cspc = dataOut.tmp_cspectra_i
596 clean_coh_spectra = dataOut.tmp_spectra_c
597 clean_coh_cspectra = dataOut.tmp_cspectra_c
598 jspectra = dataOut.data_spc+clean_coh_spectra
599 nHeights = len(dataOut.heightList) # nhei
600 nProf = int(dataOut.nProfiles)
601 dataOut.nProfiles = nProf
602 dataOut.data_param = None
603 dataOut.data_paramC = None
604 dataOut.code = numpy.array([[-1.,-1.,1.],[1.,1.,-1.]])
605 #M=600
606 #N=200
607 dataOut.flagDecodeData=True
608 M = int(dataOut.normFactor)
609 N = int(dataOut.nFFTPoints)
610 dataOut.nFFTPoints = N
611 dataOut.nIncohInt= int(dataOut.nIncohInt)
612 dataOut.nProfiles = int(dataOut.nProfiles)
613 dataOut.nCohInt = int(dataOut.nCohInt)
614 print('sale',dataOut.nProfiles,dataOut.nHeights)
615 #dataOut.nFFTPoints=nprofs
616 #dataOut.normFactor = nprofs
617 dataOut.channelList = channelList
618 #dataOut.ippFactor=1
619 #ipp = ipp/150*1.e-3
620 vmax = (300000000/49920000.0/2) / (dataOut.ippSeconds)
621 #dataOut.ippSeconds=ipp
622 absc = vmax*( numpy.arange(nProf,dtype='float')-nProf/2.)/nProf
623 print('sale 2',dataOut.ippSeconds,M,N)
624 print('Empieza procesamiento offline')
625 if path != None:
626 sys.path.append(path)
627 self.library = importlib.import_module(file)
628 constants = self.library.setConstants(dataOut)
629 constants['M'] = M
630 dataOut.constants = constants
631
632 groupArray = numpy.array(groupList)
633 dataOut.groupList = groupArray
634 nGroups = groupArray.shape[0]
635 #List of possible combinations
636 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
637 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
638 if dataOut.data_paramC is None:
639 dataOut.data_paramC = numpy.zeros((nGroups*4, nHeights,2))*numpy.nan
640 for i in range(nGroups):
641 coord = groupArray[i,:]
642 #Input data array
643 data = dataOut.data_spc[coord,:,:]/(M*N)
644 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
645
646 #Cross Spectra data array for Covariance Matrixes
647 ind = 0
648 for pairs in listComb:
649 pairsSel = numpy.array([coord[x],coord[y]])
650 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
651 ind += 1
652 dataCross = dataOut.data_cspc[indCross,:,:]/(M*N)
653 dataCross = dataCross**2
654 nhei = nHeights
655 poweri = numpy.sum(dataOut.data_spc[:,1:nProf-0,:],axis=1)/clean_num_aver[:,:]
656 if i == 0 : my_noises = numpy.zeros(4,dtype=float)
657 n0i = numpy.nanmin(poweri[0+i*2,0:nhei-0])/(nProf-1)
658 n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(nProf-1)
659 n0 = n0i
660 n1= n1i
661 my_noises[2*i+0] = n0
662 my_noises[2*i+1] = n1
663 snrth = -15.0 # -4 -16 -25
664 snrth = 10**(snrth/10.0)
665 jvelr = numpy.zeros(nHeights, dtype = 'float')
666 hvalid = [0]
667 coh2 = abs(dataOut.data_cspc[i,1:nProf,:])**2/(dataOut.data_spc[0+i*2,1:nProf-0,:]*dataOut.data_spc[1+i*2,1:nProf-0,:])
668 for h in range(nHeights):
669 smooth = clean_num_aver[i+1,h]
670 signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth
671 signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth
672 signal0 = signalpn0-n0
673 signal1 = signalpn1-n1
674 snr0 = numpy.sum(signal0/n0)/(nProf-1)
675 snr1 = numpy.sum(signal1/n1)/(nProf-1)
676 gamma = coh2[:,h]
677 indxs = (numpy.isfinite(list(gamma))==True).nonzero()
678 if len(indxs) >0:
679 if numpy.nanmean(gamma) > 0.07:
680 maxp0 = numpy.argmax(signal0*gamma)
681 maxp1 = numpy.argmax(signal1*gamma)
682 #print('usa gamma',numpy.nanmean(gamma))
683 else:
684 maxp0 = numpy.argmax(signal0)
685 maxp1 = numpy.argmax(signal1)
686 jvelr[h] = (absc[maxp0]+absc[maxp1])/2.
687 else: jvelr[h] = absc[0]
688 if snr0 > 0.1 and snr1 > 0.1: hvalid = numpy.concatenate((hvalid,h), axis=None)
689 #print(maxp0,absc[maxp0],snr0,jvelr[h])
690
691 if len(hvalid)> 1: fd0 = numpy.median(jvelr[hvalid[1:]])*-1
692 else: fd0 = numpy.nan
693 for h in range(nHeights):
694 d = data[:,h]
695 smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:]
696 signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth
697 signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth
698 signal0 = signalpn0-n0
699 signal1 = signalpn1-n1
700 snr0 = numpy.sum(signal0/n0)/(nProf-1)
701 snr1 = numpy.sum(signal1/n1)/(nProf-1)
702 if snr0 > snrth and snr1 > snrth and clean_num_aver[i+1,h] > 0 :
703 #Covariance Matrix
704 D = numpy.diag(d**2)
705 ind = 0
706 for pairs in listComb:
707 #Coordinates in Covariance Matrix
708 x = pairs[0]
709 y = pairs[1]
710 #Channel Index
711 S12 = dataCross[ind,:,h]
712 D12 = numpy.diag(S12)
713 #Completing Covariance Matrix with Cross Spectras
714 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
715 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
716 ind += 1
717 diagD = numpy.zeros(256)
718
719 try:
720 Dinv=numpy.linalg.inv(D)
721 L=numpy.linalg.cholesky(Dinv)
722 except:
723 Dinv = D*numpy.nan
724 L= D*numpy.nan
725 LT=L.T
726
727 dp = numpy.dot(LT,d)
728 #Initial values
729 data_spc = dataOut.data_spc[coord,:,h]
730 w = data_spc/data_spc
731 if filec != None:
732 w = self.weightf.weightfit(w,tini.tm_year,tini.tm_yday,index,h,i)
733 if (h>6)and(error1[3]<25):
734 p0 = dataOut.data_param[i,:,h-1]
735 else:
736 p0 = numpy.array(self.library.initialValuesFunction(data_spc*w, constants))# sin el i(data_spc, constants, i)
737 p0[3] = fd0
738 if filec != None:
739 p0 = self.weightf.Vrfit(p0,tini.tm_year,tini.tm_yday,index,h,i)
740 try:
741 #Least Squares
742 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
743 #minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
744 #Chi square error
745 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
746 #Error with Jacobian
747 error1 = self.library.errorFunction(minp,constants,LT)
748
749 except:
750 minp = p0*numpy.nan
751 error0 = numpy.nan
752 error1 = p0*numpy.nan
753 else :
754 data_spc = dataOut.data_spc[coord,:,h]
755 p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants))
756 minp = p0*numpy.nan
757 error0 = numpy.nan
758 error1 = p0*numpy.nan
759 if dataOut.data_param is None:
760 dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
761 dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
762
763 dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
764 dataOut.data_param[i,:,h] = minp
765 for ht in range(nHeights-1) :
766 smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam]
767 dataOut.data_paramC[4*i,ht,1] = smooth
768 signalpn0 = (clean_coh_spectra[i*2 ,1:(nProf-0),ht])/smooth #coh_spectra
769 signalpn1 = (clean_coh_spectra[i*2+1,1:(nProf-0),ht])/smooth
770 val0 = (signalpn0 > 0).nonzero()
771 val0 = val0[0]
772 if len(val0) == 0 : val0_npoints = nProf
773 else : val0_npoints = len(val0)
774
775 val1 = (signalpn1 > 0).nonzero()
776 val1 = val1[0]
777 if len(val1) == 0 : val1_npoints = nProf
778 else : val1_npoints = len(val1)
779
780 dataOut.data_paramC[0+4*i,ht,0] = numpy.sum((signalpn0/val0_npoints))/n0
781 dataOut.data_paramC[1+4*i,ht,0] = numpy.sum((signalpn1/val1_npoints))/n1
782
783 signal0 = (signalpn0-n0)
784 vali = (signal0 < 0).nonzero()
785 vali = vali[0]
786 if len(vali) > 0 : signal0[vali] = 0
787 signal1 = (signalpn1-n1)
788 vali = (signal1 < 0).nonzero()
789 vali = vali[0]
790 if len(vali) > 0 : signal1[vali] = 0
791 snr0 = numpy.sum(signal0/n0)/(nProf-1)
792 snr1 = numpy.sum(signal1/n1)/(nProf-1)
793 doppler = absc[1:]
794 if snr0 >= snrth and snr1 >= snrth and smooth :
795 signalpn0_n0 = signalpn0
796 signalpn0_n0[val0] = signalpn0[val0] - n0
797 mom0 = self.moments(doppler,signalpn0-n0,nProf)
798 signalpn1_n1 = signalpn1
799 signalpn1_n1[val1] = signalpn1[val1] - n1
800 mom1 = self.moments(doppler,signalpn1_n1,nProf)
801 dataOut.data_paramC[2+4*i,ht,0] = (mom0[0]+mom1[0])/2.
802 dataOut.data_paramC[3+4*i,ht,0] = (mom0[1]+mom1[1])/2.
803
804 dataOut.data_spc = jspectra
805 dataOut.spc_noise = my_noises*nProf*M
806 if numpy.any(proc): dataOut.spc_noise = my_noises*nProf*M
807 if getSNR:
808 listChannels = groupArray.reshape((groupArray.size))
809 listChannels.sort()
810
811 dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], my_noises[listChannels])
812 return dataOut
813
814 def __residFunction(self, p, dp, LT, constants):
815
816 fm = self.library.modelFunction(p, constants)
817 fmp=numpy.dot(LT,fm)
818 return dp-fmp
819
820 def __getSNR(self, z, noise):
821
822 avg = numpy.average(z, axis=1)
823 SNR = (avg.T-noise)/noise
824 SNR = SNR.T
825 return SNR
826
827 def __chisq(self, p, chindex, hindex):
828 #similar to Resid but calculates CHI**2
829 [LT,d,fm]=setupLTdfm(p,chindex,hindex)
830 dp=numpy.dot(LT,d)
831 fmp=numpy.dot(LT,fm)
832 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
833 return chisq
Show file before | Show file after | Hide comments
@@ -1,966 +1,966
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 class GetSNR(Operation):
433 class GetSNR(Operation):
434 '''
434 '''
435 Written by R. Flores
435 Written by R. Flores
436 '''
436 '''
437 """Operation to get SNR.
437 """Operation to get SNR.
438
438
439 Parameters:
439 Parameters:
440 -----------
440 -----------
441
441
442 Example
442 Example
443 --------
443 --------
444
444
445 op = proc_unit.addOperation(name='GetSNR', optype='other')
445 op = proc_unit.addOperation(name='GetSNR', optype='other')
446
446
447 """
447 """
448
448
449 def __init__(self, **kwargs):
449 def __init__(self, **kwargs):
450
450
451 Operation.__init__(self, **kwargs)
451 Operation.__init__(self, **kwargs)
452
452
453 def run(self,dataOut):
453 def run(self,dataOut):
454
454
455 noise = dataOut.getNoise(ymin_index=-10) #RegiΓ³n superior donde solo deberΓ­a de haber ruido
455 noise = dataOut.getNoise(ymin_index=-10) #RegiΓ³n superior donde solo deberΓ­a de haber ruido
456 dataOut.data_snr = (dataOut.data_spc.sum(axis=1)-noise[:,None]*dataOut.nFFTPoints)/(noise[:,None]*dataOut.nFFTPoints) #It works apparently
456 dataOut.data_snr = (dataOut.data_spc.sum(axis=1)-noise[:,None]*dataOut.nFFTPoints)/(noise[:,None]*dataOut.nFFTPoints) #It works apparently
457 dataOut.snl = numpy.log10(dataOut.data_snr)
457 dataOut.snl = numpy.log10(dataOut.data_snr)
458 dataOut.snl = numpy.where(dataOut.data_snr<.01, numpy.nan, dataOut.snl)
458 dataOut.snl = numpy.where(dataOut.data_snr<.01, numpy.nan, dataOut.snl)
459
459
460 return dataOut
460 return dataOut
461
461
462 class removeDC(Operation):
462 class removeDC(Operation):
463
463
464 def run(self, dataOut, mode=2):
464 def run(self, dataOut, mode=2):
465 self.dataOut = dataOut
465 self.dataOut = dataOut
466 jspectra = self.dataOut.data_spc
466 jspectra = self.dataOut.data_spc
467 jcspectra = self.dataOut.data_cspc
467 jcspectra = self.dataOut.data_cspc
468
468
469 num_chan = jspectra.shape[0]
469 num_chan = jspectra.shape[0]
470 num_hei = jspectra.shape[2]
470 num_hei = jspectra.shape[2]
471
471
472 if jcspectra is not None:
472 if jcspectra is not None:
473 jcspectraExist = True
473 jcspectraExist = True
474 num_pairs = jcspectra.shape[0]
474 num_pairs = jcspectra.shape[0]
475 else:
475 else:
476 jcspectraExist = False
476 jcspectraExist = False
477
477
478 freq_dc = int(jspectra.shape[1] / 2)
478 freq_dc = int(jspectra.shape[1] / 2)
479 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
479 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
480 ind_vel = ind_vel.astype(int)
480 ind_vel = ind_vel.astype(int)
481
481
482 if ind_vel[0] < 0:
482 if ind_vel[0] < 0:
483 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
483 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
484
484
485 if mode == 1:
485 if mode == 1:
486 jspectra[:, freq_dc, :] = (
486 jspectra[:, freq_dc, :] = (
487 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
487 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
488
488
489 if jcspectraExist:
489 if jcspectraExist:
490 jcspectra[:, freq_dc, :] = (
490 jcspectra[:, freq_dc, :] = (
491 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
491 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
492
492
493 if mode == 2:
493 if mode == 2:
494
494
495 vel = numpy.array([-2, -1, 1, 2])
495 vel = numpy.array([-2, -1, 1, 2])
496 xx = numpy.zeros([4, 4])
496 xx = numpy.zeros([4, 4])
497
497
498 for fil in range(4):
498 for fil in range(4):
499 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
499 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
500
500
501 xx_inv = numpy.linalg.inv(xx)
501 xx_inv = numpy.linalg.inv(xx)
502 xx_aux = xx_inv[0, :]
502 xx_aux = xx_inv[0, :]
503
503
504 for ich in range(num_chan):
504 for ich in range(num_chan):
505 yy = jspectra[ich, ind_vel, :]
505 yy = jspectra[ich, ind_vel, :]
506 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
506 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
507
507
508 junkid = jspectra[ich, freq_dc, :] <= 0
508 junkid = jspectra[ich, freq_dc, :] <= 0
509 cjunkid = sum(junkid)
509 cjunkid = sum(junkid)
510
510
511 if cjunkid.any():
511 if cjunkid.any():
512 jspectra[ich, freq_dc, junkid.nonzero()] = (
512 jspectra[ich, freq_dc, junkid.nonzero()] = (
513 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
513 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
514
514
515 if jcspectraExist:
515 if jcspectraExist:
516 for ip in range(num_pairs):
516 for ip in range(num_pairs):
517 yy = jcspectra[ip, ind_vel, :]
517 yy = jcspectra[ip, ind_vel, :]
518 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
518 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
519
519
520 self.dataOut.data_spc = jspectra
520 self.dataOut.data_spc = jspectra
521 self.dataOut.data_cspc = jcspectra
521 self.dataOut.data_cspc = jcspectra
522
522
523 return self.dataOut
523 return self.dataOut
524
524
525 class removeInterference(Operation):
525 class removeInterference(Operation):
526
526
527 def removeInterference2(self):
527 def removeInterference2(self):
528
528
529 cspc = self.dataOut.data_cspc
529 cspc = self.dataOut.data_cspc
530 spc = self.dataOut.data_spc
530 spc = self.dataOut.data_spc
531 Heights = numpy.arange(cspc.shape[2])
531 Heights = numpy.arange(cspc.shape[2])
532 realCspc = numpy.abs(cspc)
532 realCspc = numpy.abs(cspc)
533
533
534 for i in range(cspc.shape[0]):
534 for i in range(cspc.shape[0]):
535 LinePower= numpy.sum(realCspc[i], axis=0)
535 LinePower= numpy.sum(realCspc[i], axis=0)
536 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
536 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
537 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
537 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
538 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
538 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
539 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
539 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
540 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
540 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
541
541
542
542
543 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
543 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
544 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
544 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
545 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
545 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
546 cspc[i,InterferenceRange,:] = numpy.NaN
546 cspc[i,InterferenceRange,:] = numpy.NaN
547
547
548 self.dataOut.data_cspc = cspc
548 self.dataOut.data_cspc = cspc
549
549
550 def removeInterference(self, interf=2, hei_interf=None, nhei_interf=None, offhei_interf=None):
550 def removeInterference(self, interf=2, hei_interf=None, nhei_interf=None, offhei_interf=None):
551
551
552 jspectra = self.dataOut.data_spc
552 jspectra = self.dataOut.data_spc
553 jcspectra = self.dataOut.data_cspc
553 jcspectra = self.dataOut.data_cspc
554 jnoise = self.dataOut.getNoise()
554 jnoise = self.dataOut.getNoise()
555 num_incoh = self.dataOut.nIncohInt
555 num_incoh = self.dataOut.nIncohInt
556
556
557 num_channel = jspectra.shape[0]
557 num_channel = jspectra.shape[0]
558 num_prof = jspectra.shape[1]
558 num_prof = jspectra.shape[1]
559 num_hei = jspectra.shape[2]
559 num_hei = jspectra.shape[2]
560
560
561 # hei_interf
561 # hei_interf
562 if hei_interf is None:
562 if hei_interf is None:
563 count_hei = int(num_hei / 2)
563 count_hei = int(num_hei / 2)
564 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
564 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
565 hei_interf = numpy.asarray(hei_interf)[0]
565 hei_interf = numpy.asarray(hei_interf)[0]
566 # nhei_interf
566 # nhei_interf
567 if (nhei_interf == None):
567 if (nhei_interf == None):
568 nhei_interf = 5
568 nhei_interf = 5
569 if (nhei_interf < 1):
569 if (nhei_interf < 1):
570 nhei_interf = 1
570 nhei_interf = 1
571 if (nhei_interf > count_hei):
571 if (nhei_interf > count_hei):
572 nhei_interf = count_hei
572 nhei_interf = count_hei
573 if (offhei_interf == None):
573 if (offhei_interf == None):
574 offhei_interf = 0
574 offhei_interf = 0
575
575
576 ind_hei = list(range(num_hei))
576 ind_hei = list(range(num_hei))
577 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
577 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
578 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
578 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
579 mask_prof = numpy.asarray(list(range(num_prof)))
579 mask_prof = numpy.asarray(list(range(num_prof)))
580 num_mask_prof = mask_prof.size
580 num_mask_prof = mask_prof.size
581 comp_mask_prof = [0, num_prof / 2]
581 comp_mask_prof = [0, num_prof / 2]
582
582
583 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
583 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
584 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
584 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
585 jnoise = numpy.nan
585 jnoise = numpy.nan
586 noise_exist = jnoise[0] < numpy.Inf
586 noise_exist = jnoise[0] < numpy.Inf
587
587
588 # Subrutina de Remocion de la Interferencia
588 # Subrutina de Remocion de la Interferencia
589 for ich in range(num_channel):
589 for ich in range(num_channel):
590 # Se ordena los espectros segun su potencia (menor a mayor)
590 # Se ordena los espectros segun su potencia (menor a mayor)
591 power = jspectra[ich, mask_prof, :]
591 power = jspectra[ich, mask_prof, :]
592 power = power[:, hei_interf]
592 power = power[:, hei_interf]
593 power = power.sum(axis=0)
593 power = power.sum(axis=0)
594 psort = power.ravel().argsort()
594 psort = power.ravel().argsort()
595
595
596 # Se estima la interferencia promedio en los Espectros de Potencia empleando
596 # Se estima la interferencia promedio en los Espectros de Potencia empleando
597 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
597 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
598 offhei_interf, nhei_interf + offhei_interf))]]]
598 offhei_interf, nhei_interf + offhei_interf))]]]
599
599
600 if noise_exist:
600 if noise_exist:
601 # tmp_noise = jnoise[ich] / num_prof
601 # tmp_noise = jnoise[ich] / num_prof
602 tmp_noise = jnoise[ich]
602 tmp_noise = jnoise[ich]
603 junkspc_interf = junkspc_interf - tmp_noise
603 junkspc_interf = junkspc_interf - tmp_noise
604 #junkspc_interf[:,comp_mask_prof] = 0
604 #junkspc_interf[:,comp_mask_prof] = 0
605
605
606 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
606 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
607 jspc_interf = jspc_interf.transpose()
607 jspc_interf = jspc_interf.transpose()
608 # Calculando el espectro de interferencia promedio
608 # Calculando el espectro de interferencia promedio
609 noiseid = numpy.where(
609 noiseid = numpy.where(
610 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
610 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
611 noiseid = noiseid[0]
611 noiseid = noiseid[0]
612 cnoiseid = noiseid.size
612 cnoiseid = noiseid.size
613 interfid = numpy.where(
613 interfid = numpy.where(
614 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
614 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
615 interfid = interfid[0]
615 interfid = interfid[0]
616 cinterfid = interfid.size
616 cinterfid = interfid.size
617
617
618 if (cnoiseid > 0):
618 if (cnoiseid > 0):
619 jspc_interf[noiseid] = 0
619 jspc_interf[noiseid] = 0
620
620
621 # Expandiendo los perfiles a limpiar
621 # Expandiendo los perfiles a limpiar
622 if (cinterfid > 0):
622 if (cinterfid > 0):
623 new_interfid = (
623 new_interfid = (
624 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
624 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
625 new_interfid = numpy.asarray(new_interfid)
625 new_interfid = numpy.asarray(new_interfid)
626 new_interfid = {x for x in new_interfid}
626 new_interfid = {x for x in new_interfid}
627 new_interfid = numpy.array(list(new_interfid))
627 new_interfid = numpy.array(list(new_interfid))
628 new_cinterfid = new_interfid.size
628 new_cinterfid = new_interfid.size
629 else:
629 else:
630 new_cinterfid = 0
630 new_cinterfid = 0
631
631
632 for ip in range(new_cinterfid):
632 for ip in range(new_cinterfid):
633 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
633 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
634 jspc_interf[new_interfid[ip]
634 jspc_interf[new_interfid[ip]
635 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
635 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
636
636
637 jspectra[ich, :, ind_hei] = jspectra[ich, :,
637 jspectra[ich, :, ind_hei] = jspectra[ich, :,
638 ind_hei] - jspc_interf # Corregir indices
638 ind_hei] - jspc_interf # Corregir indices
639
639
640 # Removiendo la interferencia del punto de mayor interferencia
640 # Removiendo la interferencia del punto de mayor interferencia
641 ListAux = jspc_interf[mask_prof].tolist()
641 ListAux = jspc_interf[mask_prof].tolist()
642 maxid = ListAux.index(max(ListAux))
642 maxid = ListAux.index(max(ListAux))
643
643
644 if cinterfid > 0:
644 if cinterfid > 0:
645 for ip in range(cinterfid * (interf == 2) - 1):
645 for ip in range(cinterfid * (interf == 2) - 1):
646 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
646 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
647 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
647 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
648 cind = len(ind)
648 cind = len(ind)
649
649
650 if (cind > 0):
650 if (cind > 0):
651 jspectra[ich, interfid[ip], ind] = tmp_noise * \
651 jspectra[ich, interfid[ip], ind] = tmp_noise * \
652 (1 + (numpy.random.uniform(cind) - 0.5) /
652 (1 + (numpy.random.uniform(cind) - 0.5) /
653 numpy.sqrt(num_incoh))
653 numpy.sqrt(num_incoh))
654
654
655 ind = numpy.array([-2, -1, 1, 2])
655 ind = numpy.array([-2, -1, 1, 2])
656 xx = numpy.zeros([4, 4])
656 xx = numpy.zeros([4, 4])
657
657
658 for id1 in range(4):
658 for id1 in range(4):
659 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
659 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
660
660
661 xx_inv = numpy.linalg.inv(xx)
661 xx_inv = numpy.linalg.inv(xx)
662 xx = xx_inv[:, 0]
662 xx = xx_inv[:, 0]
663 ind = (ind + maxid + num_mask_prof) % num_mask_prof
663 ind = (ind + maxid + num_mask_prof) % num_mask_prof
664 yy = jspectra[ich, mask_prof[ind], :]
664 yy = jspectra[ich, mask_prof[ind], :]
665 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
665 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
666 yy.transpose(), xx)
666 yy.transpose(), xx)
667
667
668 indAux = (jspectra[ich, :, :] < tmp_noise *
668 indAux = (jspectra[ich, :, :] < tmp_noise *
669 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
669 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
670 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
670 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
671 (1 - 1 / numpy.sqrt(num_incoh))
671 (1 - 1 / numpy.sqrt(num_incoh))
672
672
673 # Remocion de Interferencia en el Cross Spectra
673 # Remocion de Interferencia en el Cross Spectra
674 if jcspectra is None:
674 if jcspectra is None:
675 return jspectra, jcspectra
675 return jspectra, jcspectra
676 num_pairs = int(jcspectra.size / (num_prof * num_hei))
676 num_pairs = int(jcspectra.size / (num_prof * num_hei))
677 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
677 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
678
678
679 for ip in range(num_pairs):
679 for ip in range(num_pairs):
680
680
681 #-------------------------------------------
681 #-------------------------------------------
682
682
683 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
683 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
684 cspower = cspower[:, hei_interf]
684 cspower = cspower[:, hei_interf]
685 cspower = cspower.sum(axis=0)
685 cspower = cspower.sum(axis=0)
686
686
687 cspsort = cspower.ravel().argsort()
687 cspsort = cspower.ravel().argsort()
688 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
688 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
689 offhei_interf, nhei_interf + offhei_interf))]]]
689 offhei_interf, nhei_interf + offhei_interf))]]]
690 junkcspc_interf = junkcspc_interf.transpose()
690 junkcspc_interf = junkcspc_interf.transpose()
691 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
691 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
692
692
693 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
693 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
694
694
695 median_real = int(numpy.median(numpy.real(
695 median_real = int(numpy.median(numpy.real(
696 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
696 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
697 median_imag = int(numpy.median(numpy.imag(
697 median_imag = int(numpy.median(numpy.imag(
698 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
698 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
699 comp_mask_prof = [int(e) for e in comp_mask_prof]
699 comp_mask_prof = [int(e) for e in comp_mask_prof]
700 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
700 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
701 median_real, median_imag)
701 median_real, median_imag)
702
702
703 for iprof in range(num_prof):
703 for iprof in range(num_prof):
704 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
704 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
705 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
705 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
706
706
707 # Removiendo la Interferencia
707 # Removiendo la Interferencia
708 jcspectra[ip, :, ind_hei] = jcspectra[ip,
708 jcspectra[ip, :, ind_hei] = jcspectra[ip,
709 :, ind_hei] - jcspc_interf
709 :, ind_hei] - jcspc_interf
710
710
711 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
711 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
712 maxid = ListAux.index(max(ListAux))
712 maxid = ListAux.index(max(ListAux))
713
713
714 ind = numpy.array([-2, -1, 1, 2])
714 ind = numpy.array([-2, -1, 1, 2])
715 xx = numpy.zeros([4, 4])
715 xx = numpy.zeros([4, 4])
716
716
717 for id1 in range(4):
717 for id1 in range(4):
718 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
718 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
719
719
720 xx_inv = numpy.linalg.inv(xx)
720 xx_inv = numpy.linalg.inv(xx)
721 xx = xx_inv[:, 0]
721 xx = xx_inv[:, 0]
722
722
723 ind = (ind + maxid + num_mask_prof) % num_mask_prof
723 ind = (ind + maxid + num_mask_prof) % num_mask_prof
724 yy = jcspectra[ip, mask_prof[ind], :]
724 yy = jcspectra[ip, mask_prof[ind], :]
725 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
725 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
726
726
727 # Guardar Resultados
727 # Guardar Resultados
728 self.dataOut.data_spc = jspectra
728 self.dataOut.data_spc = jspectra
729 self.dataOut.data_cspc = jcspectra
729 self.dataOut.data_cspc = jcspectra
730
730
731 return 1
731 return 1
732
732
733 def run(self, dataOut, interf=2,hei_interf=None, nhei_interf=None, offhei_interf=None, mode=1):
733 def run(self, dataOut, interf=2,hei_interf=None, nhei_interf=None, offhei_interf=None, mode=1):
734
734
735 self.dataOut = dataOut
735 self.dataOut = dataOut
736
736
737 if mode == 1:
737 if mode == 1:
738 self.removeInterference(interf=2,hei_interf=None, nhei_interf=None, offhei_interf=None)
738 self.removeInterference(interf=2,hei_interf=None, nhei_interf=None, offhei_interf=None)
739 elif mode == 2:
739 elif mode == 2:
740 self.removeInterference2()
740 self.removeInterference2()
741
741
742 return self.dataOut
742 return self.dataOut
743
743
744
744
745 class deflip(Operation):
745 class deflip(Operation):
746
746
747 def run(self, dataOut):
747 def run(self, dataOut):
748 # arreglo 1: (num_chan, num_profiles, num_heights)
748 # arreglo 1: (num_chan, num_profiles, num_heights)
749 self.dataOut = dataOut
749 self.dataOut = dataOut
750
750
751 # JULIA-oblicua, indice 2
751 # JULIA-oblicua, indice 2
752 # arreglo 2: (num_profiles, num_heights)
752 # arreglo 2: (num_profiles, num_heights)
753 jspectra = self.dataOut.data_spc[2]
753 jspectra = self.dataOut.data_spc[2]
754 jspectra_tmp=numpy.zeros(jspectra.shape)
754 jspectra_tmp=numpy.zeros(jspectra.shape)
755 num_profiles=jspectra.shape[0]
755 num_profiles=jspectra.shape[0]
756 freq_dc = int(num_profiles / 2)
756 freq_dc = int(num_profiles / 2)
757 # Flip con for
757 # Flip con for
758 for j in range(num_profiles):
758 for j in range(num_profiles):
759 jspectra_tmp[num_profiles-j-1]= jspectra[j]
759 jspectra_tmp[num_profiles-j-1]= jspectra[j]
760 # Intercambio perfil de DC con perfil inmediato anterior
760 # Intercambio perfil de DC con perfil inmediato anterior
761 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
761 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
762 jspectra_tmp[freq_dc]= jspectra[freq_dc]
762 jspectra_tmp[freq_dc]= jspectra[freq_dc]
763 # canal modificado es re-escrito en el arreglo de canales
763 # canal modificado es re-escrito en el arreglo de canales
764 self.dataOut.data_spc[2] = jspectra_tmp
764 self.dataOut.data_spc[2] = jspectra_tmp
765
765
766 return self.dataOut
766 return self.dataOut
767
767
768
768
769 class IncohInt(Operation):
769 class IncohInt(Operation):
770
770
771 __profIndex = 0
771 __profIndex = 0
772 __withOverapping = False
772 __withOverapping = False
773
773
774 __byTime = False
774 __byTime = False
775 __initime = None
775 __initime = None
776 __lastdatatime = None
776 __lastdatatime = None
777 __integrationtime = None
777 __integrationtime = None
778
778
779 __buffer_spc = None
779 __buffer_spc = None
780 __buffer_cspc = None
780 __buffer_cspc = None
781 __buffer_dc = None
781 __buffer_dc = None
782
782
783 __dataReady = False
783 __dataReady = False
784
784
785 __timeInterval = None
785 __timeInterval = None
786
786
787 n = None
787 n = None
788
788
789 def __init__(self):
789 def __init__(self):
790
790
791 Operation.__init__(self)
791 Operation.__init__(self)
792
792
793 def setup(self, n=None, timeInterval=None, overlapping=False):
793 def setup(self, n=None, timeInterval=None, overlapping=False):
794 """
794 """
795 Set the parameters of the integration class.
795 Set the parameters of the integration class.
796
796
797 Inputs:
797 Inputs:
798
798
799 n : Number of coherent integrations
799 n : Number of coherent integrations
800 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
800 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
801 overlapping :
801 overlapping :
802
802
803 """
803 """
804
804
805 self.__initime = None
805 self.__initime = None
806 self.__lastdatatime = 0
806 self.__lastdatatime = 0
807
807
808 self.__buffer_spc = 0
808 self.__buffer_spc = 0
809 self.__buffer_cspc = 0
809 self.__buffer_cspc = 0
810 self.__buffer_dc = 0
810 self.__buffer_dc = 0
811
811
812 self.__profIndex = 0
812 self.__profIndex = 0
813 self.__dataReady = False
813 self.__dataReady = False
814 self.__byTime = False
814 self.__byTime = False
815
815
816 if n is None and timeInterval is None:
816 if n is None and timeInterval is None:
817 raise ValueError("n or timeInterval should be specified ...")
817 raise ValueError("n or timeInterval should be specified ...")
818
818
819 if n is not None:
819 if n is not None:
820 self.n = int(n)
820 self.n = int(n)
821 else:
821 else:
822
822
823 self.__integrationtime = int(timeInterval)
823 self.__integrationtime = int(timeInterval)
824 self.n = None
824 self.n = None
825 self.__byTime = True
825 self.__byTime = True
826
826
827 def putData(self, data_spc, data_cspc, data_dc):
827 def putData(self, data_spc, data_cspc, data_dc):
828 """
828 """
829 Add a profile to the __buffer_spc and increase in one the __profileIndex
829 Add a profile to the __buffer_spc and increase in one the __profileIndex
830
830
831 """
831 """
832
832
833 self.__buffer_spc += data_spc
833 self.__buffer_spc += data_spc
834
834
835 if data_cspc is None:
835 if data_cspc is None:
836 self.__buffer_cspc = None
836 self.__buffer_cspc = None
837 else:
837 else:
838 self.__buffer_cspc += data_cspc
838 self.__buffer_cspc += data_cspc
839
839
840 if data_dc is None:
840 if data_dc is None:
841 self.__buffer_dc = None
841 self.__buffer_dc = None
842 else:
842 else:
843 self.__buffer_dc += data_dc
843 self.__buffer_dc += data_dc
844
844
845 self.__profIndex += 1
845 self.__profIndex += 1
846
846
847 return
847 return
848
848
849 def pushData(self):
849 def pushData(self):
850 """
850 """
851 Return the sum of the last profiles and the profiles used in the sum.
851 Return the sum of the last profiles and the profiles used in the sum.
852
852
853 Affected:
853 Affected:
854
854
855 self.__profileIndex
855 self.__profileIndex
856
856
857 """
857 """
858
858
859 data_spc = self.__buffer_spc
859 data_spc = self.__buffer_spc
860 data_cspc = self.__buffer_cspc
860 data_cspc = self.__buffer_cspc
861 data_dc = self.__buffer_dc
861 data_dc = self.__buffer_dc
862 n = self.__profIndex
862 n = self.__profIndex
863
863
864 self.__buffer_spc = 0
864 self.__buffer_spc = 0
865 self.__buffer_cspc = 0
865 self.__buffer_cspc = 0
866 self.__buffer_dc = 0
866 self.__buffer_dc = 0
867 self.__profIndex = 0
867 self.__profIndex = 0
868
868
869 return data_spc, data_cspc, data_dc, n
869 return data_spc, data_cspc, data_dc, n
870
870
871 def byProfiles(self, *args):
871 def byProfiles(self, *args):
872
872
873 self.__dataReady = False
873 self.__dataReady = False
874 avgdata_spc = None
874 avgdata_spc = None
875 avgdata_cspc = None
875 avgdata_cspc = None
876 avgdata_dc = None
876 avgdata_dc = None
877
877
878 self.putData(*args)
878 self.putData(*args)
879
879
880 if self.__profIndex == self.n:
880 if self.__profIndex == self.n:
881
881
882 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
882 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
883 self.n = n
883 self.n = n
884 self.__dataReady = True
884 self.__dataReady = True
885
885
886 return avgdata_spc, avgdata_cspc, avgdata_dc
886 return avgdata_spc, avgdata_cspc, avgdata_dc
887
887
888 def byTime(self, datatime, *args):
888 def byTime(self, datatime, *args):
889
889
890 self.__dataReady = False
890 self.__dataReady = False
891 avgdata_spc = None
891 avgdata_spc = None
892 avgdata_cspc = None
892 avgdata_cspc = None
893 avgdata_dc = None
893 avgdata_dc = None
894
894
895 self.putData(*args)
895 self.putData(*args)
896
896
897 if (datatime - self.__initime) >= self.__integrationtime:
897 if (datatime - self.__initime) >= self.__integrationtime:
898 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
898 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
899 self.n = n
899 self.n = n
900 self.__dataReady = True
900 self.__dataReady = True
901
901
902 return avgdata_spc, avgdata_cspc, avgdata_dc
902 return avgdata_spc, avgdata_cspc, avgdata_dc
903
903
904 def integrate(self, datatime, *args):
904 def integrate(self, datatime, *args):
905
905
906 if self.__profIndex == 0:
906 if self.__profIndex == 0:
907 self.__initime = datatime
907 self.__initime = datatime
908
908
909 if self.__byTime:
909 if self.__byTime:
910 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
910 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
911 datatime, *args)
911 datatime, *args)
912 else:
912 else:
913 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
913 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
914
914
915 if not self.__dataReady:
915 if not self.__dataReady:
916 return None, None, None, None
916 return None, None, None, None
917
917
918 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
918 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
919
919
920 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
920 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
921 if n == 1:
921 if n == 1:
922 return dataOut
922 return dataOut
923
923
924 dataOut.flagNoData = True
924 dataOut.flagNoData = True
925
925
926 if not self.isConfig:
926 if not self.isConfig:
927 self.setup(n, timeInterval, overlapping)
927 self.setup(n, timeInterval, overlapping)
928 self.isConfig = True
928 self.isConfig = True
929
929
930 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
930 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
931 dataOut.data_spc,
931 dataOut.data_spc,
932 dataOut.data_cspc,
932 dataOut.data_cspc,
933 dataOut.data_dc)
933 dataOut.data_dc)
934
934
935 if self.__dataReady:
935 if self.__dataReady:
936
936
937 dataOut.data_spc = avgdata_spc
937 dataOut.data_spc = avgdata_spc
938 dataOut.data_cspc = avgdata_cspc
938 dataOut.data_cspc = avgdata_cspc
939 dataOut.data_dc = avgdata_dc
939 dataOut.data_dc = avgdata_dc
940 dataOut.nIncohInt *= self.n
940 dataOut.nIncohInt *= self.n
941 dataOut.utctime = avgdatatime
941 dataOut.utctime = avgdatatime
942 dataOut.flagNoData = False
942 dataOut.flagNoData = False
943
943
944 return dataOut
944 return dataOut
945
945
946 class dopplerFlip(Operation):
946 class dopplerFlip(Operation):
947
947
948 def run(self, dataOut):
948 def run(self, dataOut, chann = None):
949 # arreglo 1: (num_chan, num_profiles, num_heights)
949 # arreglo 1: (num_chan, num_profiles, num_heights)
950 self.dataOut = dataOut
950 self.dataOut = dataOut
951 # JULIA-oblicua, indice 2
951 # JULIA-oblicua, indice 2
952 # arreglo 2: (num_profiles, num_heights)
952 # arreglo 2: (num_profiles, num_heights)
953 jspectra = self.dataOut.data_spc[2]
953 jspectra = self.dataOut.data_spc[chann]
954 jspectra_tmp = numpy.zeros(jspectra.shape)
954 jspectra_tmp = numpy.zeros(jspectra.shape)
955 num_profiles = jspectra.shape[0]
955 num_profiles = jspectra.shape[0]
956 freq_dc = int(num_profiles / 2)
956 freq_dc = int(num_profiles / 2)
957 # Flip con for
957 # Flip con for
958 for j in range(num_profiles):
958 for j in range(num_profiles):
959 jspectra_tmp[num_profiles-j-1]= jspectra[j]
959 jspectra_tmp[num_profiles-j-1]= jspectra[j]
960 # Intercambio perfil de DC con perfil inmediato anterior
960 # Intercambio perfil de DC con perfil inmediato anterior
961 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
961 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
962 jspectra_tmp[freq_dc]= jspectra[freq_dc]
962 jspectra_tmp[freq_dc]= jspectra[freq_dc]
963 # canal modificado es re-escrito en el arreglo de canales
963 # canal modificado es re-escrito en el arreglo de canales
964 self.dataOut.data_spc[2] = jspectra_tmp
964 self.dataOut.data_spc[chann] = jspectra_tmp
965
965
966 return self.dataOut No newline at end of file
966 return self.dataOut
General Comments 0
You need to be logged in to leave comments. Login now