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