##// END OF EJS Templates
schain
RSS
Update jroproc_base.py Queue size y jroproc_spectra.py
Alexander Valdez -
r1668:8880a814f9ea
parent child
Show More
Side by Side Unified
Context file:
r1668:8880a814f9ea Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,203 +1,203
1 1 '''
2 2 Base clases to create Processing units and operations, the MPDecorator
3 3 must be used in plotting and writing operations to allow to run as an
4 4 external process.
5 5 '''
6
6 # repositorio master
7 7 import inspect
8 8 import zmq
9 9 import time
10 10 import pickle
11 11 import traceback
12 12 from threading import Thread
13 13 from multiprocessing import Process, Queue
14 14 from schainpy.utils import log
15 15
16 16
17 17 class ProcessingUnit(object):
18 18 '''
19 19 Base class to create Signal Chain Units
20 20 '''
21 21
22 22 proc_type = 'processing'
23 23
24 24 def __init__(self):
25 25
26 26 self.dataIn = None
27 27 self.dataOut = None
28 28 self.isConfig = False
29 29 self.operations = []
30 30
31 31 def setInput(self, unit):
32 32
33 33 self.dataIn = unit.dataOut
34 34
35 35 def getAllowedArgs(self):
36 36 if hasattr(self, '__attrs__'):
37 37 return self.__attrs__
38 38 else:
39 39 return inspect.getargspec(self.run).args
40 40
41 41 def addOperation(self, conf, operation):
42 42 '''
43 43 '''
44 44
45 45 self.operations.append((operation, conf.type, conf.getKwargs()))
46 46
47 47 def getOperationObj(self, objId):
48 48
49 49 if objId not in list(self.operations.keys()):
50 50 return None
51 51
52 52 return self.operations[objId]
53 53
54 54 def call(self, **kwargs):
55 55 '''
56 56 '''
57 57
58 58 try:
59 59 if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
60 60 return self.dataIn.isReady()
61 61 elif self.dataIn is None or not self.dataIn.error:
62 62 self.run(**kwargs)
63 63 elif self.dataIn.error:
64 64 self.dataOut.error = self.dataIn.error
65 65 self.dataOut.flagNoData = True
66 66 except:
67 67 err = traceback.format_exc()
68 68 if 'SchainWarning' in err:
69 69 log.warning(err.split('SchainWarning:')[-1].split('\n')[0].strip(), self.name)
70 70 elif 'SchainError' in err:
71 71 log.error(err.split('SchainError:')[-1].split('\n')[0].strip(), self.name)
72 72 else:
73 73 log.error(err, self.name)
74 74 self.dataOut.error = True
75 75
76 76 for op, optype, opkwargs in self.operations:
77 77 if optype == 'other' and not self.dataOut.flagNoData:
78 78 self.dataOut = op.run(self.dataOut, **opkwargs)
79 79 elif optype == 'external' and not self.dataOut.flagNoData:
80 80 op.queue.put(self.dataOut)
81 81 elif optype == 'external' and self.dataOut.error:
82 82 op.queue.put(self.dataOut)
83 83
84 84 return 'Error' if self.dataOut.error else self.dataOut.isReady()
85 85
86 86 def setup(self):
87 87
88 88 raise NotImplementedError
89 89
90 90 def run(self):
91 91
92 92 raise NotImplementedError
93 93
94 94 def close(self):
95 95
96 96 return
97 97
98 98
99 99 class Operation(object):
100 100
101 101 '''
102 102 '''
103 103
104 104 proc_type = 'operation'
105 105
106 106 def __init__(self):
107 107
108 108 self.id = None
109 109 self.isConfig = False
110 110
111 111 if not hasattr(self, 'name'):
112 112 self.name = self.__class__.__name__
113 113
114 114 def getAllowedArgs(self):
115 115 if hasattr(self, '__attrs__'):
116 116 return self.__attrs__
117 117 else:
118 118 return inspect.getargspec(self.run).args
119 119
120 120 def setup(self):
121 121
122 122 self.isConfig = True
123 123
124 124 raise NotImplementedError
125 125
126 126 def run(self, dataIn, **kwargs):
127 127 """
128 128 Realiza las operaciones necesarias sobre la dataIn.data y actualiza los
129 129 atributos del objeto dataIn.
130 130
131 131 Input:
132 132
133 133 dataIn : objeto del tipo JROData
134 134
135 135 Return:
136 136
137 137 None
138 138
139 139 Affected:
140 140 __buffer : buffer de recepcion de datos.
141 141
142 142 """
143 143 if not self.isConfig:
144 144 self.setup(**kwargs)
145 145
146 146 raise NotImplementedError
147 147
148 148 def close(self):
149 149
150 150 return
151 151
152 152
153 153 def MPDecorator(BaseClass):
154 154 """
155 155 Multiprocessing class decorator
156 156
157 157 This function add multiprocessing features to a BaseClass.
158 158 """
159 159
160 160 class MPClass(BaseClass, Process):
161 161
162 162 def __init__(self, *args, **kwargs):
163 163 super(MPClass, self).__init__()
164 164 Process.__init__(self)
165 165
166 166 self.args = args
167 167 self.kwargs = kwargs
168 168 self.t = time.time()
169 169 self.op_type = 'external'
170 170 self.name = BaseClass.__name__
171 171 self.__doc__ = BaseClass.__doc__
172 172
173 173 if 'plot' in self.name.lower() and not self.name.endswith('_'):
174 174 self.name = '{}{}'.format(self.CODE.upper(), 'Plot')
175 175
176 176 self.start_time = time.time()
177 177 self.err_queue = args[3]
178 178 self.queue = Queue(maxsize=1)
179 179 self.myrun = BaseClass.run
180 180
181 181 def run(self):
182 182
183 183 while True:
184 184
185 185 dataOut = self.queue.get()
186 186
187 187 if not dataOut.error:
188 188 try:
189 189 BaseClass.run(self, dataOut, **self.kwargs)
190 190 except:
191 191 err = traceback.format_exc()
192 192 log.error(err, self.name)
193 193 else:
194 194 break
195 195
196 196 self.close()
197 197
198 198 def close(self):
199 199
200 200 BaseClass.close(self)
201 201 log.success('Done...(Time:{:4.2f} secs)'.format(time.time()-self.start_time), self.name)
202 202
203 203 return MPClass
Show file before | Show file after | Hide comments
@@ -1,898 +1,931
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Spectra processing Unit and operations
6 6
7 7 Here you will find the processing unit `SpectraProc` and several operations
8 8 to work with Spectra data type
9 9 """
10 10
11 11 import time
12 12 import itertools
13
14 13 import numpy
15 14 # repositorio
16 15 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
17 16 from schainpy.model.data.jrodata import Spectra
18 17 from schainpy.model.data.jrodata import hildebrand_sekhon
19 18 from schainpy.utils import log
20 19
21 20
22 21 class SpectraProc(ProcessingUnit):
23 22
24 23 def __init__(self):
25 24
26 25 ProcessingUnit.__init__(self)
27 26
28 27 self.buffer = None
29 28 self.firstdatatime = None
30 29 self.profIndex = 0
31 30 self.dataOut = Spectra()
32 31 self.id_min = None
33 32 self.id_max = None
34 33 self.setupReq = False #Agregar a todas las unidades de proc
35 34
36 35 def __updateSpecFromVoltage(self):
37 36
38 37 self.dataOut.timeZone = self.dataIn.timeZone
39 38 self.dataOut.dstFlag = self.dataIn.dstFlag
40 39 self.dataOut.errorCount = self.dataIn.errorCount
41 40 self.dataOut.useLocalTime = self.dataIn.useLocalTime
42 41 try:
43 42 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
44 43 except:
45 44 pass
46 45 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
47 46 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
48 47 self.dataOut.channelList = self.dataIn.channelList
49 48 self.dataOut.heightList = self.dataIn.heightList
50 49 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
51 50 self.dataOut.nProfiles = self.dataOut.nFFTPoints
52 51 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
53 52 self.dataOut.utctime = self.firstdatatime
54 53 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
55 54 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
56 55 self.dataOut.flagShiftFFT = False
57 56 self.dataOut.nCohInt = self.dataIn.nCohInt
58 57 self.dataOut.nIncohInt = 1
59 58 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
60 59 self.dataOut.frequency = self.dataIn.frequency
61 60 self.dataOut.realtime = self.dataIn.realtime
62 61 self.dataOut.azimuth = self.dataIn.azimuth
63 62 self.dataOut.zenith = self.dataIn.zenith
64 63 self.dataOut.beam.codeList = self.dataIn.beam.codeList
65 64 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
66 65 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
67 66
68 67 def __getFft(self):
69 68 """
70 69 Convierte valores de Voltaje a Spectra
71 70
72 71 Affected:
73 72 self.dataOut.data_spc
74 73 self.dataOut.data_cspc
75 74 self.dataOut.data_dc
76 75 self.dataOut.heightList
77 76 self.profIndex
78 77 self.buffer
79 78 self.dataOut.flagNoData
80 79 """
81 80 fft_volt = numpy.fft.fft(
82 81 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
83 82 fft_volt = fft_volt.astype(numpy.dtype('complex'))
84 83 dc = fft_volt[:, 0, :]
85 84
86 85 # calculo de self-spectra
87 86 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
88 87 spc = fft_volt * numpy.conjugate(fft_volt)
89 88 spc = spc.real
90 89
91 90 blocksize = 0
92 91 blocksize += dc.size
93 92 blocksize += spc.size
94 93
95 94 cspc = None
96 95 pairIndex = 0
97 96 if self.dataOut.pairsList != None:
98 97 # calculo de cross-spectra
99 98 cspc = numpy.zeros(
100 99 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
101 100 for pair in self.dataOut.pairsList:
102 101 if pair[0] not in self.dataOut.channelList:
103 102 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
104 103 str(pair), str(self.dataOut.channelList)))
105 104 if pair[1] not in self.dataOut.channelList:
106 105 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
107 106 str(pair), str(self.dataOut.channelList)))
108 107
109 108 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
110 109 numpy.conjugate(fft_volt[pair[1], :, :])
111 110 pairIndex += 1
112 111 blocksize += cspc.size
113 112
114 113 self.dataOut.data_spc = spc
115 114 self.dataOut.data_cspc = cspc
116 115 self.dataOut.data_dc = dc
117 116 self.dataOut.blockSize = blocksize
118 117 self.dataOut.flagShiftFFT = False
119 118
120 119 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False):
121 120
122 121 if self.dataIn.type == "Spectra":
123 122 self.dataOut.copy(self.dataIn)
124 123 if shift_fft:
125 124 #desplaza a la derecha en el eje 2 determinadas posiciones
126 125 shift = int(self.dataOut.nFFTPoints/2)
127 126 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
128 127
129 128 if self.dataOut.data_cspc is not None:
130 129 #desplaza a la derecha en el eje 2 determinadas posiciones
131 130 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
132 131 if pairsList:
133 132 self.__selectPairs(pairsList)
134 133
135 134 elif self.dataIn.type == "Voltage":
136 135
137 136 self.dataOut.flagNoData = True
138 137
139 138 if nFFTPoints == None:
140 139 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
141 140
142 141 if nProfiles == None:
143 142 nProfiles = nFFTPoints
144 143
145 144 if ippFactor == None:
146 145 self.dataOut.ippFactor = 1
147 146
148 147 self.dataOut.nFFTPoints = nFFTPoints
149 148
150 149 if self.buffer is None:
151 150 self.buffer = numpy.zeros((self.dataIn.nChannels,
152 151 nProfiles,
153 152 self.dataIn.nHeights),
154 153 dtype='complex')
155 154
156 155 if self.dataIn.flagDataAsBlock:
157 156 nVoltProfiles = self.dataIn.data.shape[1]
158
159 157 if nVoltProfiles == nProfiles:
160 158 self.buffer = self.dataIn.data.copy()
161 159 self.profIndex = nVoltProfiles
162 160
163 161 elif nVoltProfiles < nProfiles:
164 162
165 163 if self.profIndex == 0:
166 164 self.id_min = 0
167 165 self.id_max = nVoltProfiles
168 166
169 167 self.buffer[:, self.id_min:self.id_max,
170 168 :] = self.dataIn.data
171 169 self.profIndex += nVoltProfiles
172 170 self.id_min += nVoltProfiles
173 171 self.id_max += nVoltProfiles
172 elif nVoltProfiles > nProfiles:
173 self.reader.bypass = True
174 if self.profIndex == 0:
175 self.id_min = 0
176 self.id_max = nProfiles
177
178 self.buffer = self.dataIn.data[:, self.id_min:self.id_max,:]
179 self.profIndex += nProfiles
180 self.id_min += nProfiles
181 self.id_max += nProfiles
182 if self.id_max == nVoltProfiles:
183 self.reader.bypass = False
184
174 185 else:
175 186 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
176 187 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
177 188 self.dataOut.flagNoData = True
178 189 else:
179 190 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
180 191 self.profIndex += 1
181 192
182 193 if self.firstdatatime == None:
183 194 self.firstdatatime = self.dataIn.utctime
184 195
185 196 if self.profIndex == nProfiles:
186 197 self.__updateSpecFromVoltage()
187 198 if pairsList == None:
188 199 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
189 200 else:
190 201 self.dataOut.pairsList = pairsList
191 202 self.__getFft()
192 203 self.dataOut.flagNoData = False
193 204 self.firstdatatime = None
205 if not self.reader.bypass:
194 206 self.profIndex = 0
195 207 else:
196 208 raise ValueError("The type of input object '%s' is not valid".format(
197 209 self.dataIn.type))
198 210
199 211 def __selectPairs(self, pairsList):
200 212
201 213 if not pairsList:
202 214 return
203 215
204 216 pairs = []
205 217 pairsIndex = []
206 218
207 219 for pair in pairsList:
208 220 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
209 221 continue
210 222 pairs.append(pair)
211 223 pairsIndex.append(pairs.index(pair))
212 224
213 225 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
214 226 self.dataOut.pairsList = pairs
215 227
216 228 return
217 229
218 230 def selectFFTs(self, minFFT, maxFFT ):
219 231 """
220 232 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
221 233 minFFT<= FFT <= maxFFT
222 234 """
223 235
224 236 if (minFFT > maxFFT):
225 237 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
226 238
227 239 if (minFFT < self.dataOut.getFreqRange()[0]):
228 240 minFFT = self.dataOut.getFreqRange()[0]
229 241
230 242 if (maxFFT > self.dataOut.getFreqRange()[-1]):
231 243 maxFFT = self.dataOut.getFreqRange()[-1]
232 244
233 245 minIndex = 0
234 246 maxIndex = 0
235 247 FFTs = self.dataOut.getFreqRange()
236 248
237 249 inda = numpy.where(FFTs >= minFFT)
238 250 indb = numpy.where(FFTs <= maxFFT)
239 251
240 252 try:
241 253 minIndex = inda[0][0]
242 254 except:
243 255 minIndex = 0
244 256
245 257 try:
246 258 maxIndex = indb[0][-1]
247 259 except:
248 260 maxIndex = len(FFTs)
249 261
250 262 self.selectFFTsByIndex(minIndex, maxIndex)
251 263
252 264 return 1
253 265
254 266 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
255 267 newheis = numpy.where(
256 268 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
257 269
258 270 if hei_ref != None:
259 271 newheis = numpy.where(self.dataOut.heightList > hei_ref)
260 272
261 273 minIndex = min(newheis[0])
262 274 maxIndex = max(newheis[0])
263 275 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
264 276 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
265 277
266 278 # determina indices
267 279 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
268 280 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
269 281 avg_dB = 10 * \
270 282 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
271 283 beacon_dB = numpy.sort(avg_dB)[-nheis:]
272 284 beacon_heiIndexList = []
273 285 for val in avg_dB.tolist():
274 286 if val >= beacon_dB[0]:
275 287 beacon_heiIndexList.append(avg_dB.tolist().index(val))
276 288
277 #data_spc = data_spc[:,:,beacon_heiIndexList]
278 289 data_cspc = None
279 290 if self.dataOut.data_cspc is not None:
280 291 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
281 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
282 292
283 293 data_dc = None
284 294 if self.dataOut.data_dc is not None:
285 295 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
286 #data_dc = data_dc[:,beacon_heiIndexList]
287 296
288 297 self.dataOut.data_spc = data_spc
289 298 self.dataOut.data_cspc = data_cspc
290 299 self.dataOut.data_dc = data_dc
291 300 self.dataOut.heightList = heightList
292 301 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
293 302
294 303 return 1
295 304
296 305 def selectFFTsByIndex(self, minIndex, maxIndex):
297 306 """
298 307
299 308 """
300 309
301 310 if (minIndex < 0) or (minIndex > maxIndex):
302 311 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
303 312
304 313 if (maxIndex >= self.dataOut.nProfiles):
305 314 maxIndex = self.dataOut.nProfiles-1
306 315
307 316 #Spectra
308 317 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
309 318
310 319 data_cspc = None
311 320 if self.dataOut.data_cspc is not None:
312 321 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
313 322
314 323 data_dc = None
315 324 if self.dataOut.data_dc is not None:
316 325 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
317 326
318 327 self.dataOut.data_spc = data_spc
319 328 self.dataOut.data_cspc = data_cspc
320 329 self.dataOut.data_dc = data_dc
321 330
322 331 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
323 332 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
324 333 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
325 334
326 335 return 1
327 336
328 337 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
329 338 # validacion de rango
330 339 if minHei == None:
331 340 minHei = self.dataOut.heightList[0]
332 341
333 342 if maxHei == None:
334 343 maxHei = self.dataOut.heightList[-1]
335 344
336 345 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
337 346 print('minHei: %.2f is out of the heights range' % (minHei))
338 347 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
339 348 minHei = self.dataOut.heightList[0]
340 349
341 350 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
342 351 print('maxHei: %.2f is out of the heights range' % (maxHei))
343 352 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
344 353 maxHei = self.dataOut.heightList[-1]
345 354
346 355 # validacion de velocidades
347 356 velrange = self.dataOut.getVelRange(1)
348 357
349 358 if minVel == None:
350 359 minVel = velrange[0]
351 360
352 361 if maxVel == None:
353 362 maxVel = velrange[-1]
354 363
355 364 if (minVel < velrange[0]) or (minVel > maxVel):
356 365 print('minVel: %.2f is out of the velocity range' % (minVel))
357 366 print('minVel is setting to %.2f' % (velrange[0]))
358 367 minVel = velrange[0]
359 368
360 369 if (maxVel > velrange[-1]) or (maxVel < minVel):
361 370 print('maxVel: %.2f is out of the velocity range' % (maxVel))
362 371 print('maxVel is setting to %.2f' % (velrange[-1]))
363 372 maxVel = velrange[-1]
364 373
365 374 # seleccion de indices para rango
366 375 minIndex = 0
367 376 maxIndex = 0
368 377 heights = self.dataOut.heightList
369 378
370 379 inda = numpy.where(heights >= minHei)
371 380 indb = numpy.where(heights <= maxHei)
372 381
373 382 try:
374 383 minIndex = inda[0][0]
375 384 except:
376 385 minIndex = 0
377 386
378 387 try:
379 388 maxIndex = indb[0][-1]
380 389 except:
381 390 maxIndex = len(heights)
382 391
383 392 if (minIndex < 0) or (minIndex > maxIndex):
384 393 raise ValueError("some value in (%d,%d) is not valid" % (
385 394 minIndex, maxIndex))
386 395
387 396 if (maxIndex >= self.dataOut.nHeights):
388 397 maxIndex = self.dataOut.nHeights - 1
389 398
390 399 # seleccion de indices para velocidades
391 400 indminvel = numpy.where(velrange >= minVel)
392 401 indmaxvel = numpy.where(velrange <= maxVel)
393 402 try:
394 403 minIndexVel = indminvel[0][0]
395 404 except:
396 405 minIndexVel = 0
397 406
398 407 try:
399 408 maxIndexVel = indmaxvel[0][-1]
400 409 except:
401 410 maxIndexVel = len(velrange)
402 411
403 412 # seleccion del espectro
404 413 data_spc = self.dataOut.data_spc[:,
405 414 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
406 415 # estimacion de ruido
407 416 noise = numpy.zeros(self.dataOut.nChannels)
408 417
409 418 for channel in range(self.dataOut.nChannels):
410 419 daux = data_spc[channel, :, :]
411 420 sortdata = numpy.sort(daux, axis=None)
412 421 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
413 422
414 423 self.dataOut.noise_estimation = noise.copy()
415 424
416 425 return 1
417 426
418 427 class removeDC(Operation):
419 428
420 429 def run(self, dataOut, mode=2):
421 430 self.dataOut = dataOut
422 431 jspectra = self.dataOut.data_spc
423 432 jcspectra = self.dataOut.data_cspc
424 433
425 434 num_chan = jspectra.shape[0]
426 435 num_hei = jspectra.shape[2]
427 436
428 437 if jcspectra is not None:
429 438 jcspectraExist = True
430 439 num_pairs = jcspectra.shape[0]
431 440 else:
432 441 jcspectraExist = False
433 442
434 443 freq_dc = int(jspectra.shape[1] / 2)
435 444 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
436 445 ind_vel = ind_vel.astype(int)
437 446
438 447 if ind_vel[0] < 0:
439 448 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
440 449
441 450 if mode == 1:
442 451 jspectra[:, freq_dc, :] = (
443 452 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
444 453
445 454 if jcspectraExist:
446 455 jcspectra[:, freq_dc, :] = (
447 456 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
448 457
449 458 if mode == 2:
450 459
451 460 vel = numpy.array([-2, -1, 1, 2])
452 461 xx = numpy.zeros([4, 4])
453 462
454 463 for fil in range(4):
455 464 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
456 465
457 466 xx_inv = numpy.linalg.inv(xx)
458 467 xx_aux = xx_inv[0, :]
459 468
460 469 for ich in range(num_chan):
461 470 yy = jspectra[ich, ind_vel, :]
462 471 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
463 472
464 473 junkid = jspectra[ich, freq_dc, :] <= 0
465 474 cjunkid = sum(junkid)
466 475
467 476 if cjunkid.any():
468 477 jspectra[ich, freq_dc, junkid.nonzero()] = (
469 478 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
470 479
471 480 if jcspectraExist:
472 481 for ip in range(num_pairs):
473 482 yy = jcspectra[ip, ind_vel, :]
474 483 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
475 484
476 485 self.dataOut.data_spc = jspectra
477 486 self.dataOut.data_cspc = jcspectra
478 487
479 488 return self.dataOut
480 489
481 490 class removeInterference(Operation):
482 491
483 492 def removeInterference2(self):
484 493
485 494 cspc = self.dataOut.data_cspc
486 495 spc = self.dataOut.data_spc
487 496 Heights = numpy.arange(cspc.shape[2])
488 497 realCspc = numpy.abs(cspc)
489 498
490 499 for i in range(cspc.shape[0]):
491 500 LinePower= numpy.sum(realCspc[i], axis=0)
492 501 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
493 502 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
494 503 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
495 504 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
496 505 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
497 506
498 507
499 508 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
500 509 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
501 510 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
502 511 cspc[i,InterferenceRange,:] = numpy.NaN
503 512
504 513 self.dataOut.data_cspc = cspc
505 514
506 515 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
507 516
508 517 jspectra = self.dataOut.data_spc
509 518 jcspectra = self.dataOut.data_cspc
510 519 jnoise = self.dataOut.getNoise()
511 520 num_incoh = self.dataOut.nIncohInt
512 521
513 522 num_channel = jspectra.shape[0]
514 523 num_prof = jspectra.shape[1]
515 524 num_hei = jspectra.shape[2]
516 525
517 526 # hei_interf
518 527 if hei_interf is None:
519 528 count_hei = int(num_hei / 2)
520 529 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
521 530 hei_interf = numpy.asarray(hei_interf)[0]
522 531 # nhei_interf
523 532 if (nhei_interf == None):
524 533 nhei_interf = 5
525 534 if (nhei_interf < 1):
526 535 nhei_interf = 1
527 536 if (nhei_interf > count_hei):
528 537 nhei_interf = count_hei
529 538 if (offhei_interf == None):
530 539 offhei_interf = 0
531 540
532 541 ind_hei = list(range(num_hei))
533 542 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
534 543 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
535 544 mask_prof = numpy.asarray(list(range(num_prof)))
536 545 num_mask_prof = mask_prof.size
537 546 comp_mask_prof = [0, num_prof / 2]
538 547
539 548 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
540 549 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
541 550 jnoise = numpy.nan
542 551 noise_exist = jnoise[0] < numpy.Inf
543 552
544 553 # Subrutina de Remocion de la Interferencia
545 554 for ich in range(num_channel):
546 555 # Se ordena los espectros segun su potencia (menor a mayor)
547 556 power = jspectra[ich, mask_prof, :]
548 557 power = power[:, hei_interf]
549 558 power = power.sum(axis=0)
550 559 psort = power.ravel().argsort()
551 560
552 561 # Se estima la interferencia promedio en los Espectros de Potencia empleando
553 562 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
554 563 offhei_interf, nhei_interf + offhei_interf))]]]
555 564
556 565 if noise_exist:
557 566 # tmp_noise = jnoise[ich] / num_prof
558 567 tmp_noise = jnoise[ich]
559 568 junkspc_interf = junkspc_interf - tmp_noise
560 569 #junkspc_interf[:,comp_mask_prof] = 0
561 570
562 571 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
563 572 jspc_interf = jspc_interf.transpose()
564 573 # Calculando el espectro de interferencia promedio
565 574 noiseid = numpy.where(
566 575 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
567 576 noiseid = noiseid[0]
568 577 cnoiseid = noiseid.size
569 578 interfid = numpy.where(
570 579 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
571 580 interfid = interfid[0]
572 581 cinterfid = interfid.size
573 582
574 583 if (cnoiseid > 0):
575 584 jspc_interf[noiseid] = 0
576 585
577 586 # Expandiendo los perfiles a limpiar
578 587 if (cinterfid > 0):
579 588 new_interfid = (
580 589 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
581 590 new_interfid = numpy.asarray(new_interfid)
582 591 new_interfid = {x for x in new_interfid}
583 592 new_interfid = numpy.array(list(new_interfid))
584 593 new_cinterfid = new_interfid.size
585 594 else:
586 595 new_cinterfid = 0
587 596
588 597 for ip in range(new_cinterfid):
589 598 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
590 599 jspc_interf[new_interfid[ip]
591 600 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
592 601
593 602 jspectra[ich, :, ind_hei] = jspectra[ich, :,
594 603 ind_hei] - jspc_interf # Corregir indices
595 604
596 605 # Removiendo la interferencia del punto de mayor interferencia
597 606 ListAux = jspc_interf[mask_prof].tolist()
598 607 maxid = ListAux.index(max(ListAux))
599 608
600 609 if cinterfid > 0:
601 610 for ip in range(cinterfid * (interf == 2) - 1):
602 611 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
603 612 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
604 613 cind = len(ind)
605 614
606 615 if (cind > 0):
607 616 jspectra[ich, interfid[ip], ind] = tmp_noise * \
608 617 (1 + (numpy.random.uniform(cind) - 0.5) /
609 618 numpy.sqrt(num_incoh))
610 619
611 620 ind = numpy.array([-2, -1, 1, 2])
612 621 xx = numpy.zeros([4, 4])
613 622
614 623 for id1 in range(4):
615 624 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
616 625
617 626 xx_inv = numpy.linalg.inv(xx)
618 627 xx = xx_inv[:, 0]
619 628 ind = (ind + maxid + num_mask_prof) % num_mask_prof
620 629 yy = jspectra[ich, mask_prof[ind], :]
621 630 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
622 631 yy.transpose(), xx)
623 632
624 633 indAux = (jspectra[ich, :, :] < tmp_noise *
625 634 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
626 635 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
627 636 (1 - 1 / numpy.sqrt(num_incoh))
628 637
629 638 # Remocion de Interferencia en el Cross Spectra
630 639 if jcspectra is None:
631 640 return jspectra, jcspectra
632 641 num_pairs = int(jcspectra.size / (num_prof * num_hei))
633 642 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
634 643
635 644 for ip in range(num_pairs):
636 645
637 646 #-------------------------------------------
638 647
639 648 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
640 649 cspower = cspower[:, hei_interf]
641 650 cspower = cspower.sum(axis=0)
642 651
643 652 cspsort = cspower.ravel().argsort()
644 653 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
645 654 offhei_interf, nhei_interf + offhei_interf))]]]
646 655 junkcspc_interf = junkcspc_interf.transpose()
647 656 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
648 657
649 658 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
650 659
651 660 median_real = int(numpy.median(numpy.real(
652 661 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
653 662 median_imag = int(numpy.median(numpy.imag(
654 663 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
655 664 comp_mask_prof = [int(e) for e in comp_mask_prof]
656 665 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
657 666 median_real, median_imag)
658 667
659 668 for iprof in range(num_prof):
660 669 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
661 670 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
662 671
663 672 # Removiendo la Interferencia
664 673 jcspectra[ip, :, ind_hei] = jcspectra[ip,
665 674 :, ind_hei] - jcspc_interf
666 675
667 676 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
668 677 maxid = ListAux.index(max(ListAux))
669 678
670 679 ind = numpy.array([-2, -1, 1, 2])
671 680 xx = numpy.zeros([4, 4])
672 681
673 682 for id1 in range(4):
674 683 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
675 684
676 685 xx_inv = numpy.linalg.inv(xx)
677 686 xx = xx_inv[:, 0]
678 687
679 688 ind = (ind + maxid + num_mask_prof) % num_mask_prof
680 689 yy = jcspectra[ip, mask_prof[ind], :]
681 690 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
682 691
683 692 # Guardar Resultados
684 693 self.dataOut.data_spc = jspectra
685 694 self.dataOut.data_cspc = jcspectra
686 695
687 696 return 1
688 697
689 698 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
690 699
691 700 self.dataOut = dataOut
692 701
693 702 if mode == 1:
694 703 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
695 704 elif mode == 2:
696 705 self.removeInterference2()
697 706
698 707 return self.dataOut
699 708
700 709
710 class deflip(Operation):
711
712 def run(self, dataOut):
713 # arreglo 1: (num_chan, num_profiles, num_heights)
714 self.dataOut = dataOut
715
716 # JULIA-oblicua, indice 2
717 # arreglo 2: (num_profiles, num_heights)
718 jspectra = self.dataOut.data_spc[2]
719 jspectra_tmp=numpy.zeros(jspectra.shape)
720 num_profiles=jspectra.shape[0]
721 freq_dc = int(num_profiles / 2)
722 # Flip con for
723 for j in range(num_profiles):
724 jspectra_tmp[num_profiles-j-1]= jspectra[j]
725 # Intercambio perfil de DC con perfil inmediato anterior
726 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
727 jspectra_tmp[freq_dc]= jspectra[freq_dc]
728 # canal modificado es re-escrito en el arreglo de canales
729 self.dataOut.data_spc[2] = jspectra_tmp
730
731 return self.dataOut
732
733
701 734 class IncohInt(Operation):
702 735
703 736 __profIndex = 0
704 737 __withOverapping = False
705 738
706 739 __byTime = False
707 740 __initime = None
708 741 __lastdatatime = None
709 742 __integrationtime = None
710 743
711 744 __buffer_spc = None
712 745 __buffer_cspc = None
713 746 __buffer_dc = None
714 747
715 748 __dataReady = False
716 749
717 750 __timeInterval = None
718 751
719 752 n = None
720 753
721 754 def __init__(self):
722 755
723 756 Operation.__init__(self)
724 757
725 758 def setup(self, n=None, timeInterval=None, overlapping=False):
726 759 """
727 760 Set the parameters of the integration class.
728 761
729 762 Inputs:
730 763
731 764 n : Number of coherent integrations
732 765 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
733 766 overlapping :
734 767
735 768 """
736 769
737 770 self.__initime = None
738 771 self.__lastdatatime = 0
739 772
740 773 self.__buffer_spc = 0
741 774 self.__buffer_cspc = 0
742 775 self.__buffer_dc = 0
743 776
744 777 self.__profIndex = 0
745 778 self.__dataReady = False
746 779 self.__byTime = False
747 780
748 781 if n is None and timeInterval is None:
749 782 raise ValueError("n or timeInterval should be specified ...")
750 783
751 784 if n is not None:
752 785 self.n = int(n)
753 786 else:
754 787
755 788 self.__integrationtime = int(timeInterval)
756 789 self.n = None
757 790 self.__byTime = True
758 791
759 792 def putData(self, data_spc, data_cspc, data_dc):
760 793 """
761 794 Add a profile to the __buffer_spc and increase in one the __profileIndex
762 795
763 796 """
764 797
765 798 self.__buffer_spc += data_spc
766 799
767 800 if data_cspc is None:
768 801 self.__buffer_cspc = None
769 802 else:
770 803 self.__buffer_cspc += data_cspc
771 804
772 805 if data_dc is None:
773 806 self.__buffer_dc = None
774 807 else:
775 808 self.__buffer_dc += data_dc
776 809
777 810 self.__profIndex += 1
778 811
779 812 return
780 813
781 814 def pushData(self):
782 815 """
783 816 Return the sum of the last profiles and the profiles used in the sum.
784 817
785 818 Affected:
786 819
787 820 self.__profileIndex
788 821
789 822 """
790 823
791 824 data_spc = self.__buffer_spc
792 825 data_cspc = self.__buffer_cspc
793 826 data_dc = self.__buffer_dc
794 827 n = self.__profIndex
795 828
796 829 self.__buffer_spc = 0
797 830 self.__buffer_cspc = 0
798 831 self.__buffer_dc = 0
799 832 self.__profIndex = 0
800 833
801 834 return data_spc, data_cspc, data_dc, n
802 835
803 836 def byProfiles(self, *args):
804 837
805 838 self.__dataReady = False
806 839 avgdata_spc = None
807 840 avgdata_cspc = None
808 841 avgdata_dc = None
809 842
810 843 self.putData(*args)
811 844
812 845 if self.__profIndex == self.n:
813 846
814 847 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
815 848 self.n = n
816 849 self.__dataReady = True
817 850
818 851 return avgdata_spc, avgdata_cspc, avgdata_dc
819 852
820 853 def byTime(self, datatime, *args):
821 854
822 855 self.__dataReady = False
823 856 avgdata_spc = None
824 857 avgdata_cspc = None
825 858 avgdata_dc = None
826 859
827 860 self.putData(*args)
828 861
829 862 if (datatime - self.__initime) >= self.__integrationtime:
830 863 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
831 864 self.n = n
832 865 self.__dataReady = True
833 866
834 867 return avgdata_spc, avgdata_cspc, avgdata_dc
835 868
836 869 def integrate(self, datatime, *args):
837 870
838 871 if self.__profIndex == 0:
839 872 self.__initime = datatime
840 873
841 874 if self.__byTime:
842 875 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
843 876 datatime, *args)
844 877 else:
845 878 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
846 879
847 880 if not self.__dataReady:
848 881 return None, None, None, None
849 882
850 883 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
851 884
852 885 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
853 886 if n == 1:
854 887 return dataOut
855 888
856 889 dataOut.flagNoData = True
857 890
858 891 if not self.isConfig:
859 892 self.setup(n, timeInterval, overlapping)
860 893 self.isConfig = True
861 894
862 895 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
863 896 dataOut.data_spc,
864 897 dataOut.data_cspc,
865 898 dataOut.data_dc)
866 899
867 900 if self.__dataReady:
868 901
869 902 dataOut.data_spc = avgdata_spc
870 903 dataOut.data_cspc = avgdata_cspc
871 904 dataOut.data_dc = avgdata_dc
872 905 dataOut.nIncohInt *= self.n
873 906 dataOut.utctime = avgdatatime
874 907 dataOut.flagNoData = False
875 908
876 909 return dataOut
877 910
878 911 class dopplerFlip(Operation):
879 912
880 913 def run(self, dataOut):
881 914 # arreglo 1: (num_chan, num_profiles, num_heights)
882 915 self.dataOut = dataOut
883 916 # JULIA-oblicua, indice 2
884 917 # arreglo 2: (num_profiles, num_heights)
885 918 jspectra = self.dataOut.data_spc[2]
886 919 jspectra_tmp = numpy.zeros(jspectra.shape)
887 920 num_profiles = jspectra.shape[0]
888 921 freq_dc = int(num_profiles / 2)
889 922 # Flip con for
890 923 for j in range(num_profiles):
891 924 jspectra_tmp[num_profiles-j-1]= jspectra[j]
892 925 # Intercambio perfil de DC con perfil inmediato anterior
893 926 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
894 927 jspectra_tmp[freq_dc]= jspectra[freq_dc]
895 928 # canal modificado es re-escrito en el arreglo de canales
896 929 self.dataOut.data_spc[2] = jspectra_tmp
897 930
898 931 return self.dataOut No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now