##// END OF EJS Templates
schain-jc
RSS
Changes of voltageProc 2to3
George Yong -
r1183:cb5edef78583
parent child
Show More
Side by Side Unified
Context file:
r1183:cb5edef78583 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,951 +1,952
1 1 import itertools
2 2
3 3 import numpy
4 4
5 5 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
6 6 from schainpy.model.data.jrodata import Spectra
7 7 from schainpy.model.data.jrodata import hildebrand_sekhon
8 8 from schainpy.utils import log
9 9
10 10 @MPDecorator
11 11 class SpectraProc(ProcessingUnit):
12 12
13 13
14 14 def __init__(self):
15 15
16 16 ProcessingUnit.__init__(self)
17 17
18 18 self.buffer = None
19 19 self.firstdatatime = None
20 20 self.profIndex = 0
21 21 self.dataOut = Spectra()
22 22 self.id_min = None
23 23 self.id_max = None
24 24 self.setupReq = False #Agregar a todas las unidades de proc
25 25
26 26 def __updateSpecFromVoltage(self):
27 27
28 28 self.dataOut.timeZone = self.dataIn.timeZone
29 29 self.dataOut.dstFlag = self.dataIn.dstFlag
30 30 self.dataOut.errorCount = self.dataIn.errorCount
31 31 self.dataOut.useLocalTime = self.dataIn.useLocalTime
32 32 try:
33 33 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
34 34 except:
35 35 pass
36 36 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
37 37 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
38 38 self.dataOut.channelList = self.dataIn.channelList
39 39 self.dataOut.heightList = self.dataIn.heightList
40 40 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
41 41
42 42 self.dataOut.nBaud = self.dataIn.nBaud
43 43 self.dataOut.nCode = self.dataIn.nCode
44 44 self.dataOut.code = self.dataIn.code
45 45 self.dataOut.nProfiles = self.dataOut.nFFTPoints
46 46
47 47 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
48 48 self.dataOut.utctime = self.firstdatatime
49 49 # asumo q la data esta decodificada
50 50 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
51 51 # asumo q la data esta sin flip
52 52 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
53 53 self.dataOut.flagShiftFFT = False
54 54
55 55 self.dataOut.nCohInt = self.dataIn.nCohInt
56 56 self.dataOut.nIncohInt = 1
57 57
58 58 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
59 59
60 60 self.dataOut.frequency = self.dataIn.frequency
61 61 self.dataOut.realtime = self.dataIn.realtime
62 62
63 63 self.dataOut.azimuth = self.dataIn.azimuth
64 64 self.dataOut.zenith = self.dataIn.zenith
65 65
66 66 self.dataOut.beam.codeList = self.dataIn.beam.codeList
67 67 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
68 68 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
69 69
70 70 def __getFft(self):
71 71 """
72 72 Convierte valores de Voltaje a Spectra
73 73
74 74 Affected:
75 75 self.dataOut.data_spc
76 76 self.dataOut.data_cspc
77 77 self.dataOut.data_dc
78 78 self.dataOut.heightList
79 79 self.profIndex
80 80 self.buffer
81 81 self.dataOut.flagNoData
82 82 """
83 83 fft_volt = numpy.fft.fft(
84 84 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
85 85 fft_volt = fft_volt.astype(numpy.dtype('complex'))
86 86 dc = fft_volt[:, 0, :]
87 87
88 88 # calculo de self-spectra
89 89 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
90 90 spc = fft_volt * numpy.conjugate(fft_volt)
91 91 spc = spc.real
92 92
93 93 blocksize = 0
94 94 blocksize += dc.size
95 95 blocksize += spc.size
96 96
97 97 cspc = None
98 98 pairIndex = 0
99 99 if self.dataOut.pairsList != None:
100 100 # calculo de cross-spectra
101 101 cspc = numpy.zeros(
102 102 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
103 103 for pair in self.dataOut.pairsList:
104 104 if pair[0] not in self.dataOut.channelList:
105 105 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
106 106 str(pair), str(self.dataOut.channelList)))
107 107 if pair[1] not in self.dataOut.channelList:
108 108 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
109 109 str(pair), str(self.dataOut.channelList)))
110 110
111 111 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
112 112 numpy.conjugate(fft_volt[pair[1], :, :])
113 113 pairIndex += 1
114 114 blocksize += cspc.size
115 115
116 116 self.dataOut.data_spc = spc
117 117 self.dataOut.data_cspc = cspc
118 118 self.dataOut.data_dc = dc
119 119 self.dataOut.blockSize = blocksize
120 120 self.dataOut.flagShiftFFT = True
121 121
122 122 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], ippFactor=None, shift_fft=False):
123 123
124 124 if self.dataIn.type == "Spectra":
125 125 self.dataOut.copy(self.dataIn)
126 126 if shift_fft:
127 127 #desplaza a la derecha en el eje 2 determinadas posiciones
128 128 shift = int(self.dataOut.nFFTPoints/2)
129 129 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
130 130
131 131 if self.dataOut.data_cspc is not None:
132 132 #desplaza a la derecha en el eje 2 determinadas posiciones
133 133 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
134 134
135 135 return True
136 136
137 137 if self.dataIn.type == "Voltage":
138 138
139 self.dataOut.flagNoData = True
140
139 141 if nFFTPoints == None:
140 142 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
141 143
142 144 if nProfiles == None:
143 145 nProfiles = nFFTPoints
144 146
145 147 if ippFactor == None:
146 148 ippFactor = 1
147 149
148 150 self.dataOut.ippFactor = ippFactor
149 151
150 152 self.dataOut.nFFTPoints = nFFTPoints
151 153 self.dataOut.pairsList = pairsList
152 154
153 155 if self.buffer is None:
154 156 self.buffer = numpy.zeros((self.dataIn.nChannels,
155 157 nProfiles,
156 158 self.dataIn.nHeights),
157 159 dtype='complex')
158 160
159 161 if self.dataIn.flagDataAsBlock:
160 162 # data dimension: [nChannels, nProfiles, nSamples]
161 163 nVoltProfiles = self.dataIn.data.shape[1]
162 164 # nVoltProfiles = self.dataIn.nProfiles
163 165
164 166 if nVoltProfiles == nProfiles:
165 167 self.buffer = self.dataIn.data.copy()
166 168 self.profIndex = nVoltProfiles
167 169
168 170 elif nVoltProfiles < nProfiles:
169 171
170 172 if self.profIndex == 0:
171 173 self.id_min = 0
172 174 self.id_max = nVoltProfiles
173 175
174 176 self.buffer[:, self.id_min:self.id_max,
175 177 :] = self.dataIn.data
176 178 self.profIndex += nVoltProfiles
177 179 self.id_min += nVoltProfiles
178 180 self.id_max += nVoltProfiles
179 181 else:
180 182 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
181 183 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
182 184 self.dataOut.flagNoData = True
183 185 return 0
184 186 else:
185 187 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
186 188 self.profIndex += 1
187 189
188 190 if self.firstdatatime == None:
189 191 self.firstdatatime = self.dataIn.utctime
190 192
191 193 if self.profIndex == nProfiles:
192 194 self.__updateSpecFromVoltage()
193 195 self.__getFft()
194 196
195 197 self.dataOut.flagNoData = False
196 198 self.firstdatatime = None
197 199 self.profIndex = 0
198 200
199 201 return True
200 202
201 203 raise ValueError("The type of input object '%s' is not valid" % (
202 204 self.dataIn.type))
203 205
204 206 def __selectPairs(self, pairsList):
205 207
206 208 if not pairsList:
207 209 return
208 210
209 211 pairs = []
210 212 pairsIndex = []
211 213
212 214 for pair in pairsList:
213 215 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
214 216 continue
215 217 pairs.append(pair)
216 218 pairsIndex.append(pairs.index(pair))
217 219
218 220 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
219 221 self.dataOut.pairsList = pairs
220 222
221 223 return
222 224
223 225 def __selectPairsByChannel(self, channelList=None):
224 226
225 227 if channelList == None:
226 228 return
227 229
228 230 pairsIndexListSelected = []
229 231 for pairIndex in self.dataOut.pairsIndexList:
230 232 # First pair
231 233 if self.dataOut.pairsList[pairIndex][0] not in channelList:
232 234 continue
233 235 # Second pair
234 236 if self.dataOut.pairsList[pairIndex][1] not in channelList:
235 237 continue
236 238
237 239 pairsIndexListSelected.append(pairIndex)
238 240
239 241 if not pairsIndexListSelected:
240 242 self.dataOut.data_cspc = None
241 243 self.dataOut.pairsList = []
242 244 return
243 245
244 246 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
245 247 self.dataOut.pairsList = [self.dataOut.pairsList[i]
246 248 for i in pairsIndexListSelected]
247 249
248 250 return
249 251
250 252 def selectChannels(self, channelList):
251 253
252 254 channelIndexList = []
253 255
254 256 for channel in channelList:
255 257 if channel not in self.dataOut.channelList:
256 258 raise ValueError("Error selecting channels, Channel %d is not valid.\nAvailable channels = %s" % (
257 259 channel, str(self.dataOut.channelList)))
258 260
259 261 index = self.dataOut.channelList.index(channel)
260 262 channelIndexList.append(index)
261 263
262 264 self.selectChannelsByIndex(channelIndexList)
263 265
264 266 def selectChannelsByIndex(self, channelIndexList):
265 267 """
266 268 Selecciona un bloque de datos en base a canales segun el channelIndexList
267 269
268 270 Input:
269 271 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
270 272
271 273 Affected:
272 274 self.dataOut.data_spc
273 275 self.dataOut.channelIndexList
274 276 self.dataOut.nChannels
275 277
276 278 Return:
277 279 None
278 280 """
279 281
280 282 for channelIndex in channelIndexList:
281 283 if channelIndex not in self.dataOut.channelIndexList:
282 284 raise ValueError("Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " % (
283 285 channelIndex, self.dataOut.channelIndexList))
284 286
285 287 # nChannels = len(channelIndexList)
286 288
287 289 data_spc = self.dataOut.data_spc[channelIndexList, :]
288 290 data_dc = self.dataOut.data_dc[channelIndexList, :]
289 291
290 292 self.dataOut.data_spc = data_spc
291 293 self.dataOut.data_dc = data_dc
292 294
293 295 self.dataOut.channelList = [
294 296 self.dataOut.channelList[i] for i in channelIndexList]
295 297 # self.dataOut.nChannels = nChannels
296 298
297 299 self.__selectPairsByChannel(self.dataOut.channelList)
298 300
299 301 return 1
300 302
301 303 def selectHeights(self, minHei, maxHei):
302 304 """
303 305 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
304 306 minHei <= height <= maxHei
305 307
306 308 Input:
307 309 minHei : valor minimo de altura a considerar
308 310 maxHei : valor maximo de altura a considerar
309 311
310 312 Affected:
311 313 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
312 314
313 315 Return:
314 316 1 si el metodo se ejecuto con exito caso contrario devuelve 0
315 317 """
316 318
317 319 if (minHei > maxHei):
318 320 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (
319 321 minHei, maxHei))
320 322
321 323 if (minHei < self.dataOut.heightList[0]):
322 324 minHei = self.dataOut.heightList[0]
323 325
324 326 if (maxHei > self.dataOut.heightList[-1]):
325 327 maxHei = self.dataOut.heightList[-1]
326 328
327 329 minIndex = 0
328 330 maxIndex = 0
329 331 heights = self.dataOut.heightList
330 332
331 333 inda = numpy.where(heights >= minHei)
332 334 indb = numpy.where(heights <= maxHei)
333 335
334 336 try:
335 337 minIndex = inda[0][0]
336 338 except:
337 339 minIndex = 0
338 340
339 341 try:
340 342 maxIndex = indb[0][-1]
341 343 except:
342 344 maxIndex = len(heights)
343 345
344 346 self.selectHeightsByIndex(minIndex, maxIndex)
345 347
346 348 return 1
347 349
348 350 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
349 351 newheis = numpy.where(
350 352 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
351 353
352 354 if hei_ref != None:
353 355 newheis = numpy.where(self.dataOut.heightList > hei_ref)
354 356
355 357 minIndex = min(newheis[0])
356 358 maxIndex = max(newheis[0])
357 359 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
358 360 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
359 361
360 362 # determina indices
361 363 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
362 364 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
363 365 avg_dB = 10 * \
364 366 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
365 367 beacon_dB = numpy.sort(avg_dB)[-nheis:]
366 368 beacon_heiIndexList = []
367 369 for val in avg_dB.tolist():
368 370 if val >= beacon_dB[0]:
369 371 beacon_heiIndexList.append(avg_dB.tolist().index(val))
370 372
371 373 #data_spc = data_spc[:,:,beacon_heiIndexList]
372 374 data_cspc = None
373 375 if self.dataOut.data_cspc is not None:
374 376 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
375 377 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
376 378
377 379 data_dc = None
378 380 if self.dataOut.data_dc is not None:
379 381 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
380 382 #data_dc = data_dc[:,beacon_heiIndexList]
381 383
382 384 self.dataOut.data_spc = data_spc
383 385 self.dataOut.data_cspc = data_cspc
384 386 self.dataOut.data_dc = data_dc
385 387 self.dataOut.heightList = heightList
386 388 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
387 389
388 390 return 1
389 391
390 392 def selectHeightsByIndex(self, minIndex, maxIndex):
391 393 """
392 394 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
393 395 minIndex <= index <= maxIndex
394 396
395 397 Input:
396 398 minIndex : valor de indice minimo de altura a considerar
397 399 maxIndex : valor de indice maximo de altura a considerar
398 400
399 401 Affected:
400 402 self.dataOut.data_spc
401 403 self.dataOut.data_cspc
402 404 self.dataOut.data_dc
403 405 self.dataOut.heightList
404 406
405 407 Return:
406 408 1 si el metodo se ejecuto con exito caso contrario devuelve 0
407 409 """
408 410
409 411 if (minIndex < 0) or (minIndex > maxIndex):
410 412 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
411 413 minIndex, maxIndex))
412 414
413 415 if (maxIndex >= self.dataOut.nHeights):
414 416 maxIndex = self.dataOut.nHeights - 1
415 417
416 418 # Spectra
417 419 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
418 420
419 421 data_cspc = None
420 422 if self.dataOut.data_cspc is not None:
421 423 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
422 424
423 425 data_dc = None
424 426 if self.dataOut.data_dc is not None:
425 427 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
426 428
427 429 self.dataOut.data_spc = data_spc
428 430 self.dataOut.data_cspc = data_cspc
429 431 self.dataOut.data_dc = data_dc
430 432
431 433 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
432 434
433 435 return 1
434 436
435 437 def removeDC(self, mode=2):
436 438 jspectra = self.dataOut.data_spc
437 439 jcspectra = self.dataOut.data_cspc
438 440
439 441 num_chan = jspectra.shape[0]
440 442 num_hei = jspectra.shape[2]
441 443
442 444 if jcspectra is not None:
443 445 jcspectraExist = True
444 446 num_pairs = jcspectra.shape[0]
445 447 else:
446 448 jcspectraExist = False
447 449
448 450 freq_dc = int(jspectra.shape[1] / 2)
449 451 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
450 452 ind_vel = ind_vel.astype(int)
451 453
452 454 if ind_vel[0] < 0:
453 455 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
454 456
455 457 if mode == 1:
456 458 jspectra[:, freq_dc, :] = (
457 459 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
458 460
459 461 if jcspectraExist:
460 462 jcspectra[:, freq_dc, :] = (
461 463 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
462 464
463 465 if mode == 2:
464 466
465 467 vel = numpy.array([-2, -1, 1, 2])
466 468 xx = numpy.zeros([4, 4])
467 469
468 470 for fil in range(4):
469 471 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
470 472
471 473 xx_inv = numpy.linalg.inv(xx)
472 474 xx_aux = xx_inv[0, :]
473 475
474 476 for ich in range(num_chan):
475 477 yy = jspectra[ich, ind_vel, :]
476 478 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
477 479
478 480 junkid = jspectra[ich, freq_dc, :] <= 0
479 481 cjunkid = sum(junkid)
480 482
481 483 if cjunkid.any():
482 484 jspectra[ich, freq_dc, junkid.nonzero()] = (
483 485 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
484 486
485 487 if jcspectraExist:
486 488 for ip in range(num_pairs):
487 489 yy = jcspectra[ip, ind_vel, :]
488 490 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
489 491
490 492 self.dataOut.data_spc = jspectra
491 493 self.dataOut.data_cspc = jcspectra
492 494
493 495 return 1
494 496
495 497 def removeInterference(self, interf=2, hei_interf=None, nhei_interf=None, offhei_interf=None):
496 498
497 499 jspectra = self.dataOut.data_spc
498 500 jcspectra = self.dataOut.data_cspc
499 501 jnoise = self.dataOut.getNoise()
500 502 num_incoh = self.dataOut.nIncohInt
501 503
502 504 num_channel = jspectra.shape[0]
503 505 num_prof = jspectra.shape[1]
504 506 num_hei = jspectra.shape[2]
505 507
506 508 # hei_interf
507 509 if hei_interf is None:
508 510 count_hei = num_hei / 2 # Como es entero no importa
509 511 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
510 512 hei_interf = numpy.asarray(hei_interf)[0]
511 513 # nhei_interf
512 514 if (nhei_interf == None):
513 515 nhei_interf = 5
514 516 if (nhei_interf < 1):
515 517 nhei_interf = 1
516 518 if (nhei_interf > count_hei):
517 519 nhei_interf = count_hei
518 520 if (offhei_interf == None):
519 521 offhei_interf = 0
520 522
521 523 ind_hei = list(range(num_hei))
522 524 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
523 525 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
524 526 mask_prof = numpy.asarray(list(range(num_prof)))
525 527 num_mask_prof = mask_prof.size
526 528 comp_mask_prof = [0, num_prof / 2]
527 529
528 530 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
529 531 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
530 532 jnoise = numpy.nan
531 533 noise_exist = jnoise[0] < numpy.Inf
532 534
533 535 # Subrutina de Remocion de la Interferencia
534 536 for ich in range(num_channel):
535 537 # Se ordena los espectros segun su potencia (menor a mayor)
536 538 power = jspectra[ich, mask_prof, :]
537 539 power = power[:, hei_interf]
538 540 power = power.sum(axis=0)
539 541 psort = power.ravel().argsort()
540 542
541 543 # Se estima la interferencia promedio en los Espectros de Potencia empleando
542 544 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
543 545 offhei_interf, nhei_interf + offhei_interf))]]]
544 546
545 547 if noise_exist:
546 548 # tmp_noise = jnoise[ich] / num_prof
547 549 tmp_noise = jnoise[ich]
548 550 junkspc_interf = junkspc_interf - tmp_noise
549 551 #junkspc_interf[:,comp_mask_prof] = 0
550 552
551 553 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
552 554 jspc_interf = jspc_interf.transpose()
553 555 # Calculando el espectro de interferencia promedio
554 556 noiseid = numpy.where(
555 557 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
556 558 noiseid = noiseid[0]
557 559 cnoiseid = noiseid.size
558 560 interfid = numpy.where(
559 561 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
560 562 interfid = interfid[0]
561 563 cinterfid = interfid.size
562 564
563 565 if (cnoiseid > 0):
564 566 jspc_interf[noiseid] = 0
565 567
566 568 # Expandiendo los perfiles a limpiar
567 569 if (cinterfid > 0):
568 570 new_interfid = (
569 571 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
570 572 new_interfid = numpy.asarray(new_interfid)
571 573 new_interfid = {x for x in new_interfid}
572 574 new_interfid = numpy.array(list(new_interfid))
573 575 new_cinterfid = new_interfid.size
574 576 else:
575 577 new_cinterfid = 0
576 578
577 579 for ip in range(new_cinterfid):
578 580 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
579 581 jspc_interf[new_interfid[ip]
580 582 ] = junkspc_interf[ind[nhei_interf / 2], new_interfid[ip]]
581 583
582 584 jspectra[ich, :, ind_hei] = jspectra[ich, :,
583 585 ind_hei] - jspc_interf # Corregir indices
584 586
585 587 # Removiendo la interferencia del punto de mayor interferencia
586 588 ListAux = jspc_interf[mask_prof].tolist()
587 589 maxid = ListAux.index(max(ListAux))
588 590
589 591 if cinterfid > 0:
590 592 for ip in range(cinterfid * (interf == 2) - 1):
591 593 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
592 594 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
593 595 cind = len(ind)
594 596
595 597 if (cind > 0):
596 598 jspectra[ich, interfid[ip], ind] = tmp_noise * \
597 599 (1 + (numpy.random.uniform(cind) - 0.5) /
598 600 numpy.sqrt(num_incoh))
599 601
600 602 ind = numpy.array([-2, -1, 1, 2])
601 603 xx = numpy.zeros([4, 4])
602 604
603 605 for id1 in range(4):
604 606 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
605 607
606 608 xx_inv = numpy.linalg.inv(xx)
607 609 xx = xx_inv[:, 0]
608 610 ind = (ind + maxid + num_mask_prof) % num_mask_prof
609 611 yy = jspectra[ich, mask_prof[ind], :]
610 612 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
611 613 yy.transpose(), xx)
612 614
613 615 indAux = (jspectra[ich, :, :] < tmp_noise *
614 616 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
615 617 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
616 618 (1 - 1 / numpy.sqrt(num_incoh))
617 619
618 620 # Remocion de Interferencia en el Cross Spectra
619 621 if jcspectra is None:
620 622 return jspectra, jcspectra
621 623 num_pairs = jcspectra.size / (num_prof * num_hei)
622 624 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
623 625
624 626 for ip in range(num_pairs):
625 627
626 628 #-------------------------------------------
627 629
628 630 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
629 631 cspower = cspower[:, hei_interf]
630 632 cspower = cspower.sum(axis=0)
631 633
632 634 cspsort = cspower.ravel().argsort()
633 635 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
634 636 offhei_interf, nhei_interf + offhei_interf))]]]
635 637 junkcspc_interf = junkcspc_interf.transpose()
636 638 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
637 639
638 640 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
639 641
640 642 median_real = numpy.median(numpy.real(
641 643 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof / 4))]], :]))
642 644 median_imag = numpy.median(numpy.imag(
643 645 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof / 4))]], :]))
644 646 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
645 647 median_real, median_imag)
646 648
647 649 for iprof in range(num_prof):
648 650 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
649 651 jcspc_interf[iprof] = junkcspc_interf[iprof,
650 652 ind[nhei_interf / 2]]
651 653
652 654 # Removiendo la Interferencia
653 655 jcspectra[ip, :, ind_hei] = jcspectra[ip,
654 656 :, ind_hei] - jcspc_interf
655 657
656 658 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
657 659 maxid = ListAux.index(max(ListAux))
658 660
659 661 ind = numpy.array([-2, -1, 1, 2])
660 662 xx = numpy.zeros([4, 4])
661 663
662 664 for id1 in range(4):
663 665 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
664 666
665 667 xx_inv = numpy.linalg.inv(xx)
666 668 xx = xx_inv[:, 0]
667 669
668 670 ind = (ind + maxid + num_mask_prof) % num_mask_prof
669 671 yy = jcspectra[ip, mask_prof[ind], :]
670 672 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
671 673
672 674 # Guardar Resultados
673 675 self.dataOut.data_spc = jspectra
674 676 self.dataOut.data_cspc = jcspectra
675 677
676 678 return 1
677 679
678 680 def setRadarFrequency(self, frequency=None):
679 681
680 682 if frequency != None:
681 683 self.dataOut.frequency = frequency
682 684
683 685 return 1
684 686
685 687 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
686 688 # validacion de rango
687 689 if minHei == None:
688 690 minHei = self.dataOut.heightList[0]
689 691
690 692 if maxHei == None:
691 693 maxHei = self.dataOut.heightList[-1]
692 694
693 695 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
694 696 print('minHei: %.2f is out of the heights range' % (minHei))
695 697 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
696 698 minHei = self.dataOut.heightList[0]
697 699
698 700 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
699 701 print('maxHei: %.2f is out of the heights range' % (maxHei))
700 702 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
701 703 maxHei = self.dataOut.heightList[-1]
702 704
703 705 # validacion de velocidades
704 706 velrange = self.dataOut.getVelRange(1)
705 707
706 708 if minVel == None:
707 709 minVel = velrange[0]
708 710
709 711 if maxVel == None:
710 712 maxVel = velrange[-1]
711 713
712 714 if (minVel < velrange[0]) or (minVel > maxVel):
713 715 print('minVel: %.2f is out of the velocity range' % (minVel))
714 716 print('minVel is setting to %.2f' % (velrange[0]))
715 717 minVel = velrange[0]
716 718
717 719 if (maxVel > velrange[-1]) or (maxVel < minVel):
718 720 print('maxVel: %.2f is out of the velocity range' % (maxVel))
719 721 print('maxVel is setting to %.2f' % (velrange[-1]))
720 722 maxVel = velrange[-1]
721 723
722 724 # seleccion de indices para rango
723 725 minIndex = 0
724 726 maxIndex = 0
725 727 heights = self.dataOut.heightList
726 728
727 729 inda = numpy.where(heights >= minHei)
728 730 indb = numpy.where(heights <= maxHei)
729 731
730 732 try:
731 733 minIndex = inda[0][0]
732 734 except:
733 735 minIndex = 0
734 736
735 737 try:
736 738 maxIndex = indb[0][-1]
737 739 except:
738 740 maxIndex = len(heights)
739 741
740 742 if (minIndex < 0) or (minIndex > maxIndex):
741 743 raise ValueError("some value in (%d,%d) is not valid" % (
742 744 minIndex, maxIndex))
743 745
744 746 if (maxIndex >= self.dataOut.nHeights):
745 747 maxIndex = self.dataOut.nHeights - 1
746 748
747 749 # seleccion de indices para velocidades
748 750 indminvel = numpy.where(velrange >= minVel)
749 751 indmaxvel = numpy.where(velrange <= maxVel)
750 752 try:
751 753 minIndexVel = indminvel[0][0]
752 754 except:
753 755 minIndexVel = 0
754 756
755 757 try:
756 758 maxIndexVel = indmaxvel[0][-1]
757 759 except:
758 760 maxIndexVel = len(velrange)
759 761
760 762 # seleccion del espectro
761 763 data_spc = self.dataOut.data_spc[:,
762 764 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
763 765 # estimacion de ruido
764 766 noise = numpy.zeros(self.dataOut.nChannels)
765 767
766 768 for channel in range(self.dataOut.nChannels):
767 769 daux = data_spc[channel, :, :]
768 770 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
769 771
770 772 self.dataOut.noise_estimation = noise.copy()
771 773
772 774 return 1
773 775
774 776
775 777 class IncohInt(Operation):
776 778
777 779 __profIndex = 0
778 780 __withOverapping = False
779 781
780 782 __byTime = False
781 783 __initime = None
782 784 __lastdatatime = None
783 785 __integrationtime = None
784 786
785 787 __buffer_spc = None
786 788 __buffer_cspc = None
787 789 __buffer_dc = None
788 790
789 791 __dataReady = False
790 792
791 793 __timeInterval = None
792 794
793 795 n = None
794 796
795 797 def __init__(self):
796 798
797 799 Operation.__init__(self)
798 800
799 801 def setup(self, n=None, timeInterval=None, overlapping=False):
800 802 """
801 803 Set the parameters of the integration class.
802 804
803 805 Inputs:
804 806
805 807 n : Number of coherent integrations
806 808 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
807 809 overlapping :
808 810
809 811 """
810 812
811 813 self.__initime = None
812 814 self.__lastdatatime = 0
813 815
814 816 self.__buffer_spc = 0
815 817 self.__buffer_cspc = 0
816 818 self.__buffer_dc = 0
817 819
818 820 self.__profIndex = 0
819 821 self.__dataReady = False
820 822 self.__byTime = False
821 823
822 824 if n is None and timeInterval is None:
823 825 raise ValueError("n or timeInterval should be specified ...")
824 826
825 827 if n is not None:
826 828 self.n = int(n)
827 829 else:
828 830
829 831 self.__integrationtime = int(timeInterval)
830 832 self.n = None
831 833 self.__byTime = True
832 834
833 835 def putData(self, data_spc, data_cspc, data_dc):
834 836 """
835 837 Add a profile to the __buffer_spc and increase in one the __profileIndex
836 838
837 839 """
838 840
839 841 self.__buffer_spc += data_spc
840 842
841 843 if data_cspc is None:
842 844 self.__buffer_cspc = None
843 845 else:
844 846 self.__buffer_cspc += data_cspc
845 847
846 848 if data_dc is None:
847 849 self.__buffer_dc = None
848 850 else:
849 851 self.__buffer_dc += data_dc
850 852
851 853 self.__profIndex += 1
852 854
853 855 return
854 856
855 857 def pushData(self):
856 858 """
857 859 Return the sum of the last profiles and the profiles used in the sum.
858 860
859 861 Affected:
860 862
861 863 self.__profileIndex
862 864
863 865 """
864 866
865 867 data_spc = self.__buffer_spc
866 868 data_cspc = self.__buffer_cspc
867 869 data_dc = self.__buffer_dc
868 870 n = self.__profIndex
869 871
870 872 self.__buffer_spc = 0
871 873 self.__buffer_cspc = 0
872 874 self.__buffer_dc = 0
873 875 self.__profIndex = 0
874 876
875 877 return data_spc, data_cspc, data_dc, n
876 878
877 879 def byProfiles(self, *args):
878 880
879 881 self.__dataReady = False
880 882 avgdata_spc = None
881 883 avgdata_cspc = None
882 884 avgdata_dc = None
883 885
884 886 self.putData(*args)
885 887
886 888 if self.__profIndex == self.n:
887 889
888 890 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
889 891 self.n = n
890 892 self.__dataReady = True
891 893
892 894 return avgdata_spc, avgdata_cspc, avgdata_dc
893 895
894 896 def byTime(self, datatime, *args):
895 897
896 898 self.__dataReady = False
897 899 avgdata_spc = None
898 900 avgdata_cspc = None
899 901 avgdata_dc = None
900 902
901 903 self.putData(*args)
902 904
903 905 if (datatime - self.__initime) >= self.__integrationtime:
904 906 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
905 907 self.n = n
906 908 self.__dataReady = True
907 909
908 910 return avgdata_spc, avgdata_cspc, avgdata_dc
909 911
910 912 def integrate(self, datatime, *args):
911 913
912 914 if self.__profIndex == 0:
913 915 self.__initime = datatime
914 916
915 917 if self.__byTime:
916 918 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
917 919 datatime, *args)
918 920 else:
919 921 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
920 922
921 923 if not self.__dataReady:
922 924 return None, None, None, None
923 925
924 926 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
925 927
926 928 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
927 929 if n == 1:
928 930 return
929 931
930 932 dataOut.flagNoData = True
931 933
932 934 if not self.isConfig:
933 935 self.setup(n, timeInterval, overlapping)
934 936 self.isConfig = True
935 937
936 938 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
937 939 dataOut.data_spc,
938 940 dataOut.data_cspc,
939 941 dataOut.data_dc)
940 942
941 943 if self.__dataReady:
942 944
943 945 dataOut.data_spc = avgdata_spc
944 946 dataOut.data_cspc = avgdata_cspc
945 947 dataOut.data_dc = avgdata_dc
946
947 948 dataOut.nIncohInt *= self.n
948 949 dataOut.utctime = avgdatatime
949 950 dataOut.flagNoData = False
950 951
951 952 return dataOut No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,1330 +1,1329
1 1 import sys
2 2 import numpy
3 3 from scipy import interpolate
4 4 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
5 5 from schainpy.model.data.jrodata import Voltage
6 6 from schainpy.utils import log
7 7 from time import time
8 8
9 9
10 10 @MPDecorator
11 11 class VoltageProc(ProcessingUnit):
12 12
13 13 def __init__(self):
14 14
15 15 ProcessingUnit.__init__(self)
16 16
17 17 self.dataOut = Voltage()
18 18 self.flip = 1
19 19 self.setupReq = False
20 20
21 21 def run(self):
22 22
23 23 if self.dataIn.type == 'AMISR':
24 24 self.__updateObjFromAmisrInput()
25 25
26 26 if self.dataIn.type == 'Voltage':
27 27 self.dataOut.copy(self.dataIn)
28 28
29 29 # self.dataOut.copy(self.dataIn)
30 30
31 31 def __updateObjFromAmisrInput(self):
32 32
33 33 self.dataOut.timeZone = self.dataIn.timeZone
34 34 self.dataOut.dstFlag = self.dataIn.dstFlag
35 35 self.dataOut.errorCount = self.dataIn.errorCount
36 36 self.dataOut.useLocalTime = self.dataIn.useLocalTime
37 37
38 38 self.dataOut.flagNoData = self.dataIn.flagNoData
39 39 self.dataOut.data = self.dataIn.data
40 40 self.dataOut.utctime = self.dataIn.utctime
41 41 self.dataOut.channelList = self.dataIn.channelList
42 42 #self.dataOut.timeInterval = self.dataIn.timeInterval
43 43 self.dataOut.heightList = self.dataIn.heightList
44 44 self.dataOut.nProfiles = self.dataIn.nProfiles
45 45
46 46 self.dataOut.nCohInt = self.dataIn.nCohInt
47 47 self.dataOut.ippSeconds = self.dataIn.ippSeconds
48 48 self.dataOut.frequency = self.dataIn.frequency
49 49
50 50 self.dataOut.azimuth = self.dataIn.azimuth
51 51 self.dataOut.zenith = self.dataIn.zenith
52 52
53 53 self.dataOut.beam.codeList = self.dataIn.beam.codeList
54 54 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
55 55 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
56 56 #
57 57 # pass#
58 58 #
59 59 # def init(self):
60 60 #
61 61 #
62 62 # if self.dataIn.type == 'AMISR':
63 63 # self.__updateObjFromAmisrInput()
64 64 #
65 65 # if self.dataIn.type == 'Voltage':
66 66 # self.dataOut.copy(self.dataIn)
67 67 # # No necesita copiar en cada init() los atributos de dataIn
68 68 # # la copia deberia hacerse por cada nuevo bloque de datos
69 69
70 70 def selectChannels(self, channelList):
71 71
72 72 channelIndexList = []
73 73
74 74 for channel in channelList:
75 75 if channel not in self.dataOut.channelList:
76 76 raise ValueError("Channel %d is not in %s" %(channel, str(self.dataOut.channelList)))
77 77
78 78 index = self.dataOut.channelList.index(channel)
79 79 channelIndexList.append(index)
80 80
81 81 self.selectChannelsByIndex(channelIndexList)
82 82
83 83 def selectChannelsByIndex(self, channelIndexList):
84 84 """
85 85 Selecciona un bloque de datos en base a canales segun el channelIndexList
86 86
87 87 Input:
88 88 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
89 89
90 90 Affected:
91 91 self.dataOut.data
92 92 self.dataOut.channelIndexList
93 93 self.dataOut.nChannels
94 94 self.dataOut.m_ProcessingHeader.totalSpectra
95 95 self.dataOut.systemHeaderObj.numChannels
96 96 self.dataOut.m_ProcessingHeader.blockSize
97 97
98 98 Return:
99 99 None
100 100 """
101 101
102 102 for channelIndex in channelIndexList:
103 103 if channelIndex not in self.dataOut.channelIndexList:
104 104 print(channelIndexList)
105 105 raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
106 106
107 107 if self.dataOut.flagDataAsBlock:
108 108 """
109 109 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
110 110 """
111 111 data = self.dataOut.data[channelIndexList,:,:]
112 112 else:
113 113 data = self.dataOut.data[channelIndexList,:]
114 114
115 115 self.dataOut.data = data
116 116 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
117 117 # self.dataOut.nChannels = nChannels
118 118
119 119 return 1
120 120
121 121 def selectHeights(self, minHei=None, maxHei=None):
122 122 """
123 123 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
124 124 minHei <= height <= maxHei
125 125
126 126 Input:
127 127 minHei : valor minimo de altura a considerar
128 128 maxHei : valor maximo de altura a considerar
129 129
130 130 Affected:
131 131 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
132 132
133 133 Return:
134 134 1 si el metodo se ejecuto con exito caso contrario devuelve 0
135 135 """
136 136
137 137 if minHei == None:
138 138 minHei = self.dataOut.heightList[0]
139 139
140 140 if maxHei == None:
141 141 maxHei = self.dataOut.heightList[-1]
142 142
143 143 if (minHei < self.dataOut.heightList[0]):
144 144 minHei = self.dataOut.heightList[0]
145 145
146 146 if (maxHei > self.dataOut.heightList[-1]):
147 147 maxHei = self.dataOut.heightList[-1]
148 148
149 149 minIndex = 0
150 150 maxIndex = 0
151 151 heights = self.dataOut.heightList
152 152
153 153 inda = numpy.where(heights >= minHei)
154 154 indb = numpy.where(heights <= maxHei)
155 155
156 156 try:
157 157 minIndex = inda[0][0]
158 158 except:
159 159 minIndex = 0
160 160
161 161 try:
162 162 maxIndex = indb[0][-1]
163 163 except:
164 164 maxIndex = len(heights)
165 165
166 166 self.selectHeightsByIndex(minIndex, maxIndex)
167 167
168 168 return 1
169 169
170 170
171 171 def selectHeightsByIndex(self, minIndex, maxIndex):
172 172 """
173 173 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
174 174 minIndex <= index <= maxIndex
175 175
176 176 Input:
177 177 minIndex : valor de indice minimo de altura a considerar
178 178 maxIndex : valor de indice maximo de altura a considerar
179 179
180 180 Affected:
181 181 self.dataOut.data
182 182 self.dataOut.heightList
183 183
184 184 Return:
185 185 1 si el metodo se ejecuto con exito caso contrario devuelve 0
186 186 """
187 187
188 188 if (minIndex < 0) or (minIndex > maxIndex):
189 189 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
190 190
191 191 if (maxIndex >= self.dataOut.nHeights):
192 192 maxIndex = self.dataOut.nHeights
193 193
194 194 #voltage
195 195 if self.dataOut.flagDataAsBlock:
196 196 """
197 197 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
198 198 """
199 199 data = self.dataOut.data[:,:, minIndex:maxIndex]
200 200 else:
201 201 data = self.dataOut.data[:, minIndex:maxIndex]
202 202
203 203 # firstHeight = self.dataOut.heightList[minIndex]
204 204
205 205 self.dataOut.data = data
206 206 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
207 207
208 208 if self.dataOut.nHeights <= 1:
209 209 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
210 210
211 211 return 1
212 212
213 213
214 214 def filterByHeights(self, window):
215 215
216 216 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
217 217
218 218 if window == None:
219 219 window = (self.dataOut.radarControllerHeaderObj.txA/self.dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
220 220
221 221 newdelta = deltaHeight * window
222 222 r = self.dataOut.nHeights % window
223 223 newheights = (self.dataOut.nHeights-r)/window
224 224
225 225 if newheights <= 1:
226 226 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(self.dataOut.nHeights, window))
227 227
228 228 if self.dataOut.flagDataAsBlock:
229 229 """
230 230 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
231 231 """
232 buffer = self.dataOut.data[:, :, 0:self.dataOut.nHeights-r]
232 buffer = self.dataOut.data[:, :, 0:int(self.dataOut.nHeights-r)]
233 233 buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nProfiles,self.dataOut.nHeights/window,window)
234 234 buffer = numpy.sum(buffer,3)
235 235
236 236 else:
237 buffer = self.dataOut.data[:,0:self.dataOut.nHeights-r]
238 buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nHeights/window,window)
237 buffer = self.dataOut.data[:,0:int(self.dataOut.nHeights-r)]
238 buffer = buffer.reshape(self.dataOut.nChannels,int(self.dataOut.nHeights/window),int(window))
239 239 buffer = numpy.sum(buffer,2)
240 240
241 241 self.dataOut.data = buffer
242 242 self.dataOut.heightList = self.dataOut.heightList[0] + numpy.arange( newheights )*newdelta
243 243 self.dataOut.windowOfFilter = window
244 244
245 245 def setH0(self, h0, deltaHeight = None):
246 246
247 247 if not deltaHeight:
248 248 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
249 249
250 250 nHeights = self.dataOut.nHeights
251 251
252 252 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
253 253
254 254 self.dataOut.heightList = newHeiRange
255 255
256 256 def deFlip(self, channelList = []):
257 257
258 258 data = self.dataOut.data.copy()
259 259
260 260 if self.dataOut.flagDataAsBlock:
261 261 flip = self.flip
262 262 profileList = list(range(self.dataOut.nProfiles))
263 263
264 264 if not channelList:
265 265 for thisProfile in profileList:
266 266 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
267 267 flip *= -1.0
268 268 else:
269 269 for thisChannel in channelList:
270 270 if thisChannel not in self.dataOut.channelList:
271 271 continue
272 272
273 273 for thisProfile in profileList:
274 274 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
275 275 flip *= -1.0
276 276
277 277 self.flip = flip
278 278
279 279 else:
280 280 if not channelList:
281 281 data[:,:] = data[:,:]*self.flip
282 282 else:
283 283 for thisChannel in channelList:
284 284 if thisChannel not in self.dataOut.channelList:
285 285 continue
286 286
287 287 data[thisChannel,:] = data[thisChannel,:]*self.flip
288 288
289 289 self.flip *= -1.
290 290
291 291 self.dataOut.data = data
292 292
293 293 def setRadarFrequency(self, frequency=None):
294 294
295 295 if frequency != None:
296 296 self.dataOut.frequency = frequency
297 297
298 298 return 1
299 299
300 300 def interpolateHeights(self, topLim, botLim):
301 301 #69 al 72 para julia
302 302 #82-84 para meteoros
303 303 if len(numpy.shape(self.dataOut.data))==2:
304 304 sampInterp = (self.dataOut.data[:,botLim-1] + self.dataOut.data[:,topLim+1])/2
305 305 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
306 306 #self.dataOut.data[:,botLim:limSup+1] = sampInterp
307 307 self.dataOut.data[:,botLim:topLim+1] = sampInterp
308 308 else:
309 309 nHeights = self.dataOut.data.shape[2]
310 310 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
311 311 y = self.dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
312 312 f = interpolate.interp1d(x, y, axis = 2)
313 313 xnew = numpy.arange(botLim,topLim+1)
314 314 ynew = f(xnew)
315 315
316 316 self.dataOut.data[:,:,botLim:topLim+1] = ynew
317 317
318 318 # import collections
319 319
320 320 class CohInt(Operation):
321 321
322 322 isConfig = False
323 323 __profIndex = 0
324 324 __byTime = False
325 325 __initime = None
326 326 __lastdatatime = None
327 327 __integrationtime = None
328 328 __buffer = None
329 329 __bufferStride = []
330 330 __dataReady = False
331 331 __profIndexStride = 0
332 332 __dataToPutStride = False
333 333 n = None
334 334
335 335 def __init__(self, **kwargs):
336 336
337 337 Operation.__init__(self, **kwargs)
338 338
339 339 # self.isConfig = False
340 340
341 341 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
342 342 """
343 343 Set the parameters of the integration class.
344 344
345 345 Inputs:
346 346
347 347 n : Number of coherent integrations
348 348 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
349 349 overlapping :
350 350 """
351 351
352 352 self.__initime = None
353 353 self.__lastdatatime = 0
354 354 self.__buffer = None
355 355 self.__dataReady = False
356 356 self.byblock = byblock
357 357 self.stride = stride
358 358
359 359 if n == None and timeInterval == None:
360 360 raise ValueError("n or timeInterval should be specified ...")
361 361
362 362 if n != None:
363 363 self.n = n
364 364 self.__byTime = False
365 365 else:
366 366 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
367 367 self.n = 9999
368 368 self.__byTime = True
369 369
370 370 if overlapping:
371 371 self.__withOverlapping = True
372 372 self.__buffer = None
373 373 else:
374 374 self.__withOverlapping = False
375 375 self.__buffer = 0
376 376
377 377 self.__profIndex = 0
378 378
379 379 def putData(self, data):
380 380
381 381 """
382 382 Add a profile to the __buffer and increase in one the __profileIndex
383 383
384 384 """
385 385
386 386 if not self.__withOverlapping:
387 387 self.__buffer += data.copy()
388 388 self.__profIndex += 1
389 389 return
390 390
391 391 #Overlapping data
392 392 nChannels, nHeis = data.shape
393 393 data = numpy.reshape(data, (1, nChannels, nHeis))
394 394
395 395 #If the buffer is empty then it takes the data value
396 396 if self.__buffer is None:
397 397 self.__buffer = data
398 398 self.__profIndex += 1
399 399 return
400 400
401 401 #If the buffer length is lower than n then stakcing the data value
402 402 if self.__profIndex < self.n:
403 403 self.__buffer = numpy.vstack((self.__buffer, data))
404 404 self.__profIndex += 1
405 405 return
406 406
407 407 #If the buffer length is equal to n then replacing the last buffer value with the data value
408 408 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
409 409 self.__buffer[self.n-1] = data
410 410 self.__profIndex = self.n
411 411 return
412 412
413 413
414 414 def pushData(self):
415 415 """
416 416 Return the sum of the last profiles and the profiles used in the sum.
417 417
418 418 Affected:
419 419
420 420 self.__profileIndex
421 421
422 422 """
423 423
424 424 if not self.__withOverlapping:
425 425 data = self.__buffer
426 426 n = self.__profIndex
427 427
428 428 self.__buffer = 0
429 429 self.__profIndex = 0
430 430
431 431 return data, n
432 432
433 433 #Integration with Overlapping
434 434 data = numpy.sum(self.__buffer, axis=0)
435 435 # print data
436 436 # raise
437 437 n = self.__profIndex
438 438
439 439 return data, n
440 440
441 441 def byProfiles(self, data):
442 442
443 443 self.__dataReady = False
444 444 avgdata = None
445 445 # n = None
446 446 # print data
447 447 # raise
448 448 self.putData(data)
449 449
450 450 if self.__profIndex == self.n:
451 451 avgdata, n = self.pushData()
452 452 self.__dataReady = True
453 453
454 454 return avgdata
455 455
456 456 def byTime(self, data, datatime):
457 457
458 458 self.__dataReady = False
459 459 avgdata = None
460 460 n = None
461 461
462 462 self.putData(data)
463 463
464 464 if (datatime - self.__initime) >= self.__integrationtime:
465 465 avgdata, n = self.pushData()
466 466 self.n = n
467 467 self.__dataReady = True
468 468
469 469 return avgdata
470 470
471 471 def integrateByStride(self, data, datatime):
472 472 # print data
473 473 if self.__profIndex == 0:
474 474 self.__buffer = [[data.copy(), datatime]]
475 475 else:
476 476 self.__buffer.append([data.copy(),datatime])
477 477 self.__profIndex += 1
478 478 self.__dataReady = False
479 479
480 480 if self.__profIndex == self.n * self.stride :
481 481 self.__dataToPutStride = True
482 482 self.__profIndexStride = 0
483 483 self.__profIndex = 0
484 484 self.__bufferStride = []
485 485 for i in range(self.stride):
486 486 current = self.__buffer[i::self.stride]
487 487 data = numpy.sum([t[0] for t in current], axis=0)
488 488 avgdatatime = numpy.average([t[1] for t in current])
489 489 # print data
490 490 self.__bufferStride.append((data, avgdatatime))
491 491
492 492 if self.__dataToPutStride:
493 493 self.__dataReady = True
494 494 self.__profIndexStride += 1
495 495 if self.__profIndexStride == self.stride:
496 496 self.__dataToPutStride = False
497 497 # print self.__bufferStride[self.__profIndexStride - 1]
498 498 # raise
499 499 return self.__bufferStride[self.__profIndexStride - 1]
500 500
501 501
502 502 return None, None
503 503
504 504 def integrate(self, data, datatime=None):
505 505
506 506 if self.__initime == None:
507 507 self.__initime = datatime
508 508
509 509 if self.__byTime:
510 510 avgdata = self.byTime(data, datatime)
511 511 else:
512 512 avgdata = self.byProfiles(data)
513 513
514 514
515 515 self.__lastdatatime = datatime
516 516
517 517 if avgdata is None:
518 518 return None, None
519 519
520 520 avgdatatime = self.__initime
521 521
522 522 deltatime = datatime - self.__lastdatatime
523 523
524 524 if not self.__withOverlapping:
525 525 self.__initime = datatime
526 526 else:
527 527 self.__initime += deltatime
528 528
529 529 return avgdata, avgdatatime
530 530
531 531 def integrateByBlock(self, dataOut):
532 532
533 533 times = int(dataOut.data.shape[1]/self.n)
534 534 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
535 535
536 536 id_min = 0
537 537 id_max = self.n
538 538
539 539 for i in range(times):
540 540 junk = dataOut.data[:,id_min:id_max,:]
541 541 avgdata[:,i,:] = junk.sum(axis=1)
542 542 id_min += self.n
543 543 id_max += self.n
544 544
545 545 timeInterval = dataOut.ippSeconds*self.n
546 546 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
547 547 self.__dataReady = True
548 548 return avgdata, avgdatatime
549 549
550 550 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
551 551
552 552 if not self.isConfig:
553 553 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
554 554 self.isConfig = True
555 555
556 556 if dataOut.flagDataAsBlock:
557 557 """
558 558 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
559 559 """
560 560 avgdata, avgdatatime = self.integrateByBlock(dataOut)
561 561 dataOut.nProfiles /= self.n
562 562 else:
563 563 if stride is None:
564 564 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
565 565 else:
566 566 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
567 567
568 568
569 569 # dataOut.timeInterval *= n
570 570 dataOut.flagNoData = True
571 571
572 572 if self.__dataReady:
573 573 dataOut.data = avgdata
574 574 dataOut.nCohInt *= self.n
575 575 dataOut.utctime = avgdatatime
576 576 # print avgdata, avgdatatime
577 577 # raise
578 578 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
579 579 dataOut.flagNoData = False
580 580 return dataOut
581 581
582 582 class Decoder(Operation):
583 583
584 584 isConfig = False
585 585 __profIndex = 0
586 586
587 587 code = None
588 588
589 589 nCode = None
590 590 nBaud = None
591 591
592 592 def __init__(self, **kwargs):
593 593
594 594 Operation.__init__(self, **kwargs)
595 595
596 596 self.times = None
597 597 self.osamp = None
598 598 # self.__setValues = False
599 599 self.isConfig = False
600 600 self.setupReq = False
601 601 def setup(self, code, osamp, dataOut):
602 602
603 603 self.__profIndex = 0
604 604
605 605 self.code = code
606 606
607 607 self.nCode = len(code)
608 608 self.nBaud = len(code[0])
609 609
610 610 if (osamp != None) and (osamp >1):
611 611 self.osamp = osamp
612 612 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
613 613 self.nBaud = self.nBaud*self.osamp
614 614
615 615 self.__nChannels = dataOut.nChannels
616 616 self.__nProfiles = dataOut.nProfiles
617 617 self.__nHeis = dataOut.nHeights
618 618
619 619 if self.__nHeis < self.nBaud:
620 620 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
621 621
622 622 #Frequency
623 623 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
624 624
625 625 __codeBuffer[:,0:self.nBaud] = self.code
626 626
627 627 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
628 628
629 629 if dataOut.flagDataAsBlock:
630 630
631 631 self.ndatadec = self.__nHeis #- self.nBaud + 1
632 632
633 633 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
634 634
635 635 else:
636 636
637 637 #Time
638 638 self.ndatadec = self.__nHeis #- self.nBaud + 1
639 639
640 640 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
641 641
642 642 def __convolutionInFreq(self, data):
643 643
644 644 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
645 645
646 646 fft_data = numpy.fft.fft(data, axis=1)
647 647
648 648 conv = fft_data*fft_code
649 649
650 650 data = numpy.fft.ifft(conv,axis=1)
651 651
652 652 return data
653 653
654 654 def __convolutionInFreqOpt(self, data):
655 655
656 656 raise NotImplementedError
657 657
658 658 def __convolutionInTime(self, data):
659 659
660 660 code = self.code[self.__profIndex]
661 661 for i in range(self.__nChannels):
662 662 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
663 663
664 664 return self.datadecTime
665 665
666 666 def __convolutionByBlockInTime(self, data):
667 667
668 668 repetitions = self.__nProfiles / self.nCode
669 669
670 670 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
671 671 junk = junk.flatten()
672 672 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
673 673 profilesList = range(self.__nProfiles)
674 674
675 675 for i in range(self.__nChannels):
676 676 for j in profilesList:
677 677 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
678 678 return self.datadecTime
679 679
680 680 def __convolutionByBlockInFreq(self, data):
681 681
682 682 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
683 683
684 684
685 685 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
686 686
687 687 fft_data = numpy.fft.fft(data, axis=2)
688 688
689 689 conv = fft_data*fft_code
690 690
691 691 data = numpy.fft.ifft(conv,axis=2)
692 692
693 693 return data
694 694
695 695
696 696 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
697 697
698 698 if dataOut.flagDecodeData:
699 699 print("This data is already decoded, recoding again ...")
700 700
701 701 if not self.isConfig:
702 702
703 703 if code is None:
704 704 if dataOut.code is None:
705 705 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
706 706
707 707 code = dataOut.code
708 708 else:
709 709 code = numpy.array(code).reshape(nCode,nBaud)
710 710 self.setup(code, osamp, dataOut)
711 711
712 712 self.isConfig = True
713 713
714 714 if mode == 3:
715 715 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
716 716
717 717 if times != None:
718 718 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
719 719
720 720 if self.code is None:
721 721 print("Fail decoding: Code is not defined.")
722 722 return
723 723
724 724 self.__nProfiles = dataOut.nProfiles
725 725 datadec = None
726 726
727 727 if mode == 3:
728 728 mode = 0
729 729
730 730 if dataOut.flagDataAsBlock:
731 731 """
732 732 Decoding when data have been read as block,
733 733 """
734 734
735 735 if mode == 0:
736 736 datadec = self.__convolutionByBlockInTime(dataOut.data)
737 737 if mode == 1:
738 738 datadec = self.__convolutionByBlockInFreq(dataOut.data)
739 739 else:
740 740 """
741 741 Decoding when data have been read profile by profile
742 742 """
743 743 if mode == 0:
744 744 datadec = self.__convolutionInTime(dataOut.data)
745 745
746 746 if mode == 1:
747 747 datadec = self.__convolutionInFreq(dataOut.data)
748 748
749 749 if mode == 2:
750 750 datadec = self.__convolutionInFreqOpt(dataOut.data)
751 751
752 752 if datadec is None:
753 753 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
754 754
755 755 dataOut.code = self.code
756 756 dataOut.nCode = self.nCode
757 757 dataOut.nBaud = self.nBaud
758 758
759 759 dataOut.data = datadec
760 760
761 761 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
762 762
763 763 dataOut.flagDecodeData = True #asumo q la data esta decodificada
764 764
765 765 if self.__profIndex == self.nCode-1:
766 766 self.__profIndex = 0
767 767 return dataOut
768 768
769 769 self.__profIndex += 1
770 770
771 771 return dataOut
772 772 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
773 773
774 774
775 775 class ProfileConcat(Operation):
776 776
777 777 isConfig = False
778 778 buffer = None
779 779
780 780 def __init__(self, **kwargs):
781 781
782 782 Operation.__init__(self, **kwargs)
783 783 self.profileIndex = 0
784 784
785 785 def reset(self):
786 786 self.buffer = numpy.zeros_like(self.buffer)
787 787 self.start_index = 0
788 788 self.times = 1
789 789
790 790 def setup(self, data, m, n=1):
791 791 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
792 792 self.nHeights = data.shape[1]#.nHeights
793 793 self.start_index = 0
794 794 self.times = 1
795 795
796 796 def concat(self, data):
797 797
798 798 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
799 799 self.start_index = self.start_index + self.nHeights
800 800
801 801 def run(self, dataOut, m):
802
803 802 dataOut.flagNoData = True
804 803
805 804 if not self.isConfig:
806 805 self.setup(dataOut.data, m, 1)
807 806 self.isConfig = True
808 807
809 808 if dataOut.flagDataAsBlock:
810 809 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
811 810
812 811 else:
813 812 self.concat(dataOut.data)
814 813 self.times += 1
815 814 if self.times > m:
816 815 dataOut.data = self.buffer
817 816 self.reset()
818 817 dataOut.flagNoData = False
819 818 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
820 819 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
821 820 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
822 821 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
823 822 dataOut.ippSeconds *= m
824 823 return dataOut
825 824
826 825 class ProfileSelector(Operation):
827 826
828 827 profileIndex = None
829 828 # Tamanho total de los perfiles
830 829 nProfiles = None
831 830
832 831 def __init__(self, **kwargs):
833 832
834 833 Operation.__init__(self, **kwargs)
835 834 self.profileIndex = 0
836 835
837 836 def incProfileIndex(self):
838 837
839 838 self.profileIndex += 1
840 839
841 840 if self.profileIndex >= self.nProfiles:
842 841 self.profileIndex = 0
843 842
844 843 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
845 844
846 845 if profileIndex < minIndex:
847 846 return False
848 847
849 848 if profileIndex > maxIndex:
850 849 return False
851 850
852 851 return True
853 852
854 853 def isThisProfileInList(self, profileIndex, profileList):
855 854
856 855 if profileIndex not in profileList:
857 856 return False
858 857
859 858 return True
860 859
861 860 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
862 861
863 862 """
864 863 ProfileSelector:
865 864
866 865 Inputs:
867 866 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
868 867
869 868 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
870 869
871 870 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
872 871
873 872 """
874 873
875 874 if rangeList is not None:
876 875 if type(rangeList[0]) not in (tuple, list):
877 876 rangeList = [rangeList]
878 877
879 878 dataOut.flagNoData = True
880 879
881 880 if dataOut.flagDataAsBlock:
882 881 """
883 882 data dimension = [nChannels, nProfiles, nHeis]
884 883 """
885 884 if profileList != None:
886 885 dataOut.data = dataOut.data[:,profileList,:]
887 886
888 887 if profileRangeList != None:
889 888 minIndex = profileRangeList[0]
890 889 maxIndex = profileRangeList[1]
891 890 profileList = list(range(minIndex, maxIndex+1))
892 891
893 892 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
894 893
895 894 if rangeList != None:
896 895
897 896 profileList = []
898 897
899 898 for thisRange in rangeList:
900 899 minIndex = thisRange[0]
901 900 maxIndex = thisRange[1]
902 901
903 902 profileList.extend(list(range(minIndex, maxIndex+1)))
904 903
905 904 dataOut.data = dataOut.data[:,profileList,:]
906 905
907 906 dataOut.nProfiles = len(profileList)
908 907 dataOut.profileIndex = dataOut.nProfiles - 1
909 908 dataOut.flagNoData = False
910 909
911 return True
910 return dataOut
912 911
913 912 """
914 913 data dimension = [nChannels, nHeis]
915 914 """
916 915
917 916 if profileList != None:
918 917
919 918 if self.isThisProfileInList(dataOut.profileIndex, profileList):
920 919
921 920 self.nProfiles = len(profileList)
922 921 dataOut.nProfiles = self.nProfiles
923 922 dataOut.profileIndex = self.profileIndex
924 923 dataOut.flagNoData = False
925 924
926 925 self.incProfileIndex()
927 return True
926 return dataOut
928 927
929 928 if profileRangeList != None:
930 929
931 930 minIndex = profileRangeList[0]
932 931 maxIndex = profileRangeList[1]
933 932
934 933 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
935 934
936 935 self.nProfiles = maxIndex - minIndex + 1
937 936 dataOut.nProfiles = self.nProfiles
938 937 dataOut.profileIndex = self.profileIndex
939 938 dataOut.flagNoData = False
940 939
941 940 self.incProfileIndex()
942 return True
941 return dataOut
943 942
944 943 if rangeList != None:
945 944
946 945 nProfiles = 0
947 946
948 947 for thisRange in rangeList:
949 948 minIndex = thisRange[0]
950 949 maxIndex = thisRange[1]
951 950
952 951 nProfiles += maxIndex - minIndex + 1
953 952
954 953 for thisRange in rangeList:
955 954
956 955 minIndex = thisRange[0]
957 956 maxIndex = thisRange[1]
958 957
959 958 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
960 959
961 960 self.nProfiles = nProfiles
962 961 dataOut.nProfiles = self.nProfiles
963 962 dataOut.profileIndex = self.profileIndex
964 963 dataOut.flagNoData = False
965 964
966 965 self.incProfileIndex()
967 966
968 967 break
969 968
970 return True
969 return dataOut
971 970
972 971
973 972 if beam != None: #beam is only for AMISR data
974 973 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
975 974 dataOut.flagNoData = False
976 975 dataOut.profileIndex = self.profileIndex
977 976
978 977 self.incProfileIndex()
979 978
980 return True
979 return dataOut
981 980
982 981 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
983 982
984 983 #return False
985 984 return dataOut
986 985
987 986 class Reshaper(Operation):
988 987
989 988 def __init__(self, **kwargs):
990 989
991 990 Operation.__init__(self, **kwargs)
992 991
993 992 self.__buffer = None
994 993 self.__nitems = 0
995 994
996 995 def __appendProfile(self, dataOut, nTxs):
997 996
998 997 if self.__buffer is None:
999 998 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1000 999 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1001 1000
1002 1001 ini = dataOut.nHeights * self.__nitems
1003 1002 end = ini + dataOut.nHeights
1004 1003
1005 1004 self.__buffer[:, ini:end] = dataOut.data
1006 1005
1007 1006 self.__nitems += 1
1008 1007
1009 1008 return int(self.__nitems*nTxs)
1010 1009
1011 1010 def __getBuffer(self):
1012 1011
1013 1012 if self.__nitems == int(1./self.__nTxs):
1014 1013
1015 1014 self.__nitems = 0
1016 1015
1017 1016 return self.__buffer.copy()
1018 1017
1019 1018 return None
1020 1019
1021 1020 def __checkInputs(self, dataOut, shape, nTxs):
1022 1021
1023 1022 if shape is None and nTxs is None:
1024 1023 raise ValueError("Reshaper: shape of factor should be defined")
1025 1024
1026 1025 if nTxs:
1027 1026 if nTxs < 0:
1028 1027 raise ValueError("nTxs should be greater than 0")
1029 1028
1030 1029 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1031 1030 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1032 1031
1033 1032 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1034 1033
1035 1034 return shape, nTxs
1036 1035
1037 1036 if len(shape) != 2 and len(shape) != 3:
1038 1037 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))
1039 1038
1040 1039 if len(shape) == 2:
1041 1040 shape_tuple = [dataOut.nChannels]
1042 1041 shape_tuple.extend(shape)
1043 1042 else:
1044 1043 shape_tuple = list(shape)
1045 1044
1046 1045 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1047 1046
1048 1047 return shape_tuple, nTxs
1049 1048
1050 1049 def run(self, dataOut, shape=None, nTxs=None):
1051 1050
1052 1051 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1053 1052
1054 1053 dataOut.flagNoData = True
1055 1054 profileIndex = None
1056 1055
1057 1056 if dataOut.flagDataAsBlock:
1058 1057
1059 1058 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1060 1059 dataOut.flagNoData = False
1061 1060
1062 1061 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1063 1062
1064 1063 else:
1065 1064
1066 1065 if self.__nTxs < 1:
1067 1066
1068 1067 self.__appendProfile(dataOut, self.__nTxs)
1069 1068 new_data = self.__getBuffer()
1070 1069
1071 1070 if new_data is not None:
1072 1071 dataOut.data = new_data
1073 1072 dataOut.flagNoData = False
1074 1073
1075 1074 profileIndex = dataOut.profileIndex*nTxs
1076 1075
1077 1076 else:
1078 1077 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1079 1078
1080 1079 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1081 1080
1082 1081 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1083 1082
1084 1083 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1085 1084
1086 1085 dataOut.profileIndex = profileIndex
1087 1086
1088 1087 dataOut.ippSeconds /= self.__nTxs
1089 1088
1090 1089 return dataOut
1091 1090
1092 1091 class SplitProfiles(Operation):
1093 1092
1094 1093 def __init__(self, **kwargs):
1095 1094
1096 1095 Operation.__init__(self, **kwargs)
1097 1096
1098 1097 def run(self, dataOut, n):
1099 1098
1100 1099 dataOut.flagNoData = True
1101 1100 profileIndex = None
1102 1101
1103 1102 if dataOut.flagDataAsBlock:
1104 1103
1105 1104 #nchannels, nprofiles, nsamples
1106 1105 shape = dataOut.data.shape
1107 1106
1108 1107 if shape[2] % n != 0:
1109 1108 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1110 1109
1111 1110 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1112 1111
1113 1112 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1114 1113 dataOut.flagNoData = False
1115 1114
1116 1115 profileIndex = int(dataOut.nProfiles/n) - 1
1117 1116
1118 1117 else:
1119 1118
1120 1119 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1121 1120
1122 1121 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1123 1122
1124 1123 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1125 1124
1126 1125 dataOut.nProfiles = int(dataOut.nProfiles*n)
1127 1126
1128 1127 dataOut.profileIndex = profileIndex
1129 1128
1130 1129 dataOut.ippSeconds /= n
1131 1130
1132 1131 return dataOut
1133 1132
1134 1133 class CombineProfiles(Operation):
1135 1134 def __init__(self, **kwargs):
1136 1135
1137 1136 Operation.__init__(self, **kwargs)
1138 1137
1139 1138 self.__remData = None
1140 1139 self.__profileIndex = 0
1141 1140
1142 1141 def run(self, dataOut, n):
1143 1142
1144 1143 dataOut.flagNoData = True
1145 1144 profileIndex = None
1146 1145
1147 1146 if dataOut.flagDataAsBlock:
1148 1147
1149 1148 #nchannels, nprofiles, nsamples
1150 1149 shape = dataOut.data.shape
1151 1150 new_shape = shape[0], shape[1]/n, shape[2]*n
1152 1151
1153 1152 if shape[1] % n != 0:
1154 1153 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1155 1154
1156 1155 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1157 1156 dataOut.flagNoData = False
1158 1157
1159 1158 profileIndex = int(dataOut.nProfiles*n) - 1
1160 1159
1161 1160 else:
1162 1161
1163 1162 #nchannels, nsamples
1164 1163 if self.__remData is None:
1165 1164 newData = dataOut.data
1166 1165 else:
1167 1166 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1168 1167
1169 1168 self.__profileIndex += 1
1170 1169
1171 1170 if self.__profileIndex < n:
1172 1171 self.__remData = newData
1173 1172 #continue
1174 1173 return
1175 1174
1176 1175 self.__profileIndex = 0
1177 1176 self.__remData = None
1178 1177
1179 1178 dataOut.data = newData
1180 1179 dataOut.flagNoData = False
1181 1180
1182 1181 profileIndex = dataOut.profileIndex/n
1183 1182
1184 1183
1185 1184 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1186 1185
1187 1186 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1188 1187
1189 1188 dataOut.nProfiles = int(dataOut.nProfiles/n)
1190 1189
1191 1190 dataOut.profileIndex = profileIndex
1192 1191
1193 1192 dataOut.ippSeconds *= n
1194 1193
1195 1194 return dataOut
1196 1195 # import collections
1197 1196 # from scipy.stats import mode
1198 1197 #
1199 1198 # class Synchronize(Operation):
1200 1199 #
1201 1200 # isConfig = False
1202 1201 # __profIndex = 0
1203 1202 #
1204 1203 # def __init__(self, **kwargs):
1205 1204 #
1206 1205 # Operation.__init__(self, **kwargs)
1207 1206 # # self.isConfig = False
1208 1207 # self.__powBuffer = None
1209 1208 # self.__startIndex = 0
1210 1209 # self.__pulseFound = False
1211 1210 #
1212 1211 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1213 1212 #
1214 1213 # #Read data
1215 1214 #
1216 1215 # powerdB = dataOut.getPower(channel = channel)
1217 1216 # noisedB = dataOut.getNoise(channel = channel)[0]
1218 1217 #
1219 1218 # self.__powBuffer.extend(powerdB.flatten())
1220 1219 #
1221 1220 # dataArray = numpy.array(self.__powBuffer)
1222 1221 #
1223 1222 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1224 1223 #
1225 1224 # maxValue = numpy.nanmax(filteredPower)
1226 1225 #
1227 1226 # if maxValue < noisedB + 10:
1228 1227 # #No se encuentra ningun pulso de transmision
1229 1228 # return None
1230 1229 #
1231 1230 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1232 1231 #
1233 1232 # if len(maxValuesIndex) < 2:
1234 1233 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1235 1234 # return None
1236 1235 #
1237 1236 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1238 1237 #
1239 1238 # #Seleccionar solo valores con un espaciamiento de nSamples
1240 1239 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1241 1240 #
1242 1241 # if len(pulseIndex) < 2:
1243 1242 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1244 1243 # return None
1245 1244 #
1246 1245 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1247 1246 #
1248 1247 # #remover senales que se distancien menos de 10 unidades o muestras
1249 1248 # #(No deberian existir IPP menor a 10 unidades)
1250 1249 #
1251 1250 # realIndex = numpy.where(spacing > 10 )[0]
1252 1251 #
1253 1252 # if len(realIndex) < 2:
1254 1253 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1255 1254 # return None
1256 1255 #
1257 1256 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1258 1257 # realPulseIndex = pulseIndex[realIndex]
1259 1258 #
1260 1259 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1261 1260 #
1262 1261 # print "IPP = %d samples" %period
1263 1262 #
1264 1263 # self.__newNSamples = dataOut.nHeights #int(period)
1265 1264 # self.__startIndex = int(realPulseIndex[0])
1266 1265 #
1267 1266 # return 1
1268 1267 #
1269 1268 #
1270 1269 # def setup(self, nSamples, nChannels, buffer_size = 4):
1271 1270 #
1272 1271 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1273 1272 # maxlen = buffer_size*nSamples)
1274 1273 #
1275 1274 # bufferList = []
1276 1275 #
1277 1276 # for i in range(nChannels):
1278 1277 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1279 1278 # maxlen = buffer_size*nSamples)
1280 1279 #
1281 1280 # bufferList.append(bufferByChannel)
1282 1281 #
1283 1282 # self.__nSamples = nSamples
1284 1283 # self.__nChannels = nChannels
1285 1284 # self.__bufferList = bufferList
1286 1285 #
1287 1286 # def run(self, dataOut, channel = 0):
1288 1287 #
1289 1288 # if not self.isConfig:
1290 1289 # nSamples = dataOut.nHeights
1291 1290 # nChannels = dataOut.nChannels
1292 1291 # self.setup(nSamples, nChannels)
1293 1292 # self.isConfig = True
1294 1293 #
1295 1294 # #Append new data to internal buffer
1296 1295 # for thisChannel in range(self.__nChannels):
1297 1296 # bufferByChannel = self.__bufferList[thisChannel]
1298 1297 # bufferByChannel.extend(dataOut.data[thisChannel])
1299 1298 #
1300 1299 # if self.__pulseFound:
1301 1300 # self.__startIndex -= self.__nSamples
1302 1301 #
1303 1302 # #Finding Tx Pulse
1304 1303 # if not self.__pulseFound:
1305 1304 # indexFound = self.__findTxPulse(dataOut, channel)
1306 1305 #
1307 1306 # if indexFound == None:
1308 1307 # dataOut.flagNoData = True
1309 1308 # return
1310 1309 #
1311 1310 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1312 1311 # self.__pulseFound = True
1313 1312 # self.__startIndex = indexFound
1314 1313 #
1315 1314 # #If pulse was found ...
1316 1315 # for thisChannel in range(self.__nChannels):
1317 1316 # bufferByChannel = self.__bufferList[thisChannel]
1318 1317 # #print self.__startIndex
1319 1318 # x = numpy.array(bufferByChannel)
1320 1319 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1321 1320 #
1322 1321 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1323 1322 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1324 1323 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1325 1324 #
1326 1325 # dataOut.data = self.__arrayBuffer
1327 1326 #
1328 1327 # self.__startIndex += self.__newNSamples
1329 1328 #
1330 1329 # return
General Comments 0
You need to be logged in to leave comments. Login now