##// END OF EJS Templates
fix meteors processing bugs
José Chávez -
r992:7a6f35bfb333
parent child
Show More
@@ -1,851 +1,851
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROHeaderIO.py 151 2012-10-31 19:00:51Z murco $
5 5 '''
6 6 import sys
7 7 import numpy
8 8 import copy
9 9 import datetime
10 10
11 11 SPEED_OF_LIGHT = 299792458
12 12 SPEED_OF_LIGHT = 3e8
13 13
14 14 BASIC_STRUCTURE = numpy.dtype([
15 15 ('nSize','<u4'),
16 16 ('nVersion','<u2'),
17 17 ('nDataBlockId','<u4'),
18 18 ('nUtime','<u4'),
19 19 ('nMilsec','<u2'),
20 20 ('nTimezone','<i2'),
21 21 ('nDstflag','<i2'),
22 22 ('nErrorCount','<u4')
23 23 ])
24 24
25 25 SYSTEM_STRUCTURE = numpy.dtype([
26 26 ('nSize','<u4'),
27 27 ('nNumSamples','<u4'),
28 28 ('nNumProfiles','<u4'),
29 29 ('nNumChannels','<u4'),
30 30 ('nADCResolution','<u4'),
31 31 ('nPCDIOBusWidth','<u4'),
32 32 ])
33 33
34 34 RADAR_STRUCTURE = numpy.dtype([
35 35 ('nSize','<u4'),
36 36 ('nExpType','<u4'),
37 37 ('nNTx','<u4'),
38 38 ('fIpp','<f4'),
39 39 ('fTxA','<f4'),
40 40 ('fTxB','<f4'),
41 41 ('nNumWindows','<u4'),
42 42 ('nNumTaus','<u4'),
43 43 ('nCodeType','<u4'),
44 44 ('nLine6Function','<u4'),
45 45 ('nLine5Function','<u4'),
46 46 ('fClock','<f4'),
47 47 ('nPrePulseBefore','<u4'),
48 48 ('nPrePulseAfter','<u4'),
49 49 ('sRangeIPP','<a20'),
50 50 ('sRangeTxA','<a20'),
51 51 ('sRangeTxB','<a20'),
52 52 ])
53 53
54 54 SAMPLING_STRUCTURE = numpy.dtype([('h0','<f4'),('dh','<f4'),('nsa','<u4')])
55 55
56 56
57 57 PROCESSING_STRUCTURE = numpy.dtype([
58 58 ('nSize','<u4'),
59 59 ('nDataType','<u4'),
60 60 ('nSizeOfDataBlock','<u4'),
61 61 ('nProfilesperBlock','<u4'),
62 62 ('nDataBlocksperFile','<u4'),
63 63 ('nNumWindows','<u4'),
64 64 ('nProcessFlags','<u4'),
65 65 ('nCoherentIntegrations','<u4'),
66 66 ('nIncoherentIntegrations','<u4'),
67 67 ('nTotalSpectra','<u4')
68 68 ])
69 69
70 70 class Header(object):
71 71
72 72 def __init__(self):
73 73 raise NotImplementedError
74 74
75 75 def copy(self):
76 76 return copy.deepcopy(self)
77 77
78 78 def read(self):
79 79
80 80 raise NotImplementedError
81 81
82 82 def write(self):
83 83
84 84 raise NotImplementedError
85 85
86 86 def printInfo(self):
87 87
88 88 message = "#"*50 + "\n"
89 89 message += self.__class__.__name__.upper() + "\n"
90 90 message += "#"*50 + "\n"
91 91
92 92 keyList = self.__dict__.keys()
93 93 keyList.sort()
94 94
95 95 for key in keyList:
96 96 message += "%s = %s" %(key, self.__dict__[key]) + "\n"
97 97
98 98 if "size" not in keyList:
99 99 attr = getattr(self, "size")
100 100
101 101 if attr:
102 102 message += "%s = %s" %("size", attr) + "\n"
103 103
104 104 print message
105 105
106 106 class BasicHeader(Header):
107 107
108 108 size = None
109 109 version = None
110 110 dataBlock = None
111 111 utc = None
112 112 ltc = None
113 113 miliSecond = None
114 114 timeZone = None
115 115 dstFlag = None
116 116 errorCount = None
117 117 datatime = None
118 118 __LOCALTIME = None
119 119
120 120 def __init__(self, useLocalTime=True):
121 121
122 122 self.size = 24
123 123 self.version = 0
124 124 self.dataBlock = 0
125 125 self.utc = 0
126 126 self.miliSecond = 0
127 127 self.timeZone = 0
128 128 self.dstFlag = 0
129 129 self.errorCount = 0
130 130
131 131 self.useLocalTime = useLocalTime
132 132
133 133 def read(self, fp):
134 134
135 135 self.length = 0
136 136 try:
137 137 if hasattr(fp, 'read'):
138 138 header = numpy.fromfile(fp, BASIC_STRUCTURE,1)
139 139 else:
140 140 header = numpy.fromstring(fp, BASIC_STRUCTURE,1)
141 141 except Exception, e:
142 142 print "BasicHeader: "
143 143 print e
144 144 return 0
145 145
146 146 self.size = int(header['nSize'][0])
147 147 self.version = int(header['nVersion'][0])
148 148 self.dataBlock = int(header['nDataBlockId'][0])
149 149 self.utc = int(header['nUtime'][0])
150 150 self.miliSecond = int(header['nMilsec'][0])
151 151 self.timeZone = int(header['nTimezone'][0])
152 152 self.dstFlag = int(header['nDstflag'][0])
153 153 self.errorCount = int(header['nErrorCount'][0])
154 154
155 155 if self.size < 24:
156 156 return 0
157 157
158 158 self.length = header.nbytes
159 159 return 1
160 160
161 161 def write(self, fp):
162 162
163 163 headerTuple = (self.size,self.version,self.dataBlock,self.utc,self.miliSecond,self.timeZone,self.dstFlag,self.errorCount)
164 164 header = numpy.array(headerTuple, BASIC_STRUCTURE)
165 165 header.tofile(fp)
166 166
167 167 return 1
168 168
169 169 def get_ltc(self):
170 170
171 171 return self.utc - self.timeZone*60
172 172
173 173 def set_ltc(self, value):
174 174
175 175 self.utc = value + self.timeZone*60
176 176
177 177 def get_datatime(self):
178 178
179 179 return datetime.datetime.utcfromtimestamp(self.ltc)
180 180
181 181 ltc = property(get_ltc, set_ltc)
182 182 datatime = property(get_datatime)
183 183
184 184 class SystemHeader(Header):
185 185
186 186 size = None
187 187 nSamples = None
188 188 nProfiles = None
189 189 nChannels = None
190 190 adcResolution = None
191 191 pciDioBusWidth = None
192 192
193 193 def __init__(self, nSamples=0, nProfiles=0, nChannels=0, adcResolution=14, pciDioBusWith=0):
194 194
195 195 self.size = 24
196 196 self.nSamples = nSamples
197 197 self.nProfiles = nProfiles
198 198 self.nChannels = nChannels
199 199 self.adcResolution = adcResolution
200 200 self.pciDioBusWidth = pciDioBusWith
201 201
202 202 def read(self, fp):
203 203 self.length = 0
204 204 try:
205 205 startFp = fp.tell()
206 206 except Exception, e:
207 207 startFp = None
208 208 pass
209 209
210 210 try:
211 211 if hasattr(fp, 'read'):
212 212 header = numpy.fromfile(fp, SYSTEM_STRUCTURE,1)
213 213 else:
214 214 header = numpy.fromstring(fp, SYSTEM_STRUCTURE,1)
215 215 except Exception, e:
216 216 print "System Header: " + str(e)
217 217 return 0
218 218
219 219 self.size = header['nSize'][0]
220 220 self.nSamples = header['nNumSamples'][0]
221 221 self.nProfiles = header['nNumProfiles'][0]
222 222 self.nChannels = header['nNumChannels'][0]
223 223 self.adcResolution = header['nADCResolution'][0]
224 224 self.pciDioBusWidth = header['nPCDIOBusWidth'][0]
225 225
226 226
227 227 if startFp is not None:
228 228 endFp = self.size + startFp
229 229
230 230 if fp.tell() > endFp:
231 231 sys.stderr.write("Warning %s: Size value read from System Header is lower than it has to be\n" %fp.name)
232 232 return 0
233 233
234 234 if fp.tell() < endFp:
235 235 sys.stderr.write("Warning %s: Size value read from System Header size is greater than it has to be\n" %fp.name)
236 236 return 0
237 237
238 238 self.length = header.nbytes
239 239 return 1
240 240
241 241 def write(self, fp):
242 242
243 243 headerTuple = (self.size,self.nSamples,self.nProfiles,self.nChannels,self.adcResolution,self.pciDioBusWidth)
244 244 header = numpy.array(headerTuple,SYSTEM_STRUCTURE)
245 245 header.tofile(fp)
246 246
247 247 return 1
248 248
249 249 class RadarControllerHeader(Header):
250 250
251 251 expType = None
252 252 nTx = None
253 253 ipp = None
254 254 txA = None
255 255 txB = None
256 256 nWindows = None
257 257 numTaus = None
258 258 codeType = None
259 259 line6Function = None
260 260 line5Function = None
261 261 fClock = None
262 262 prePulseBefore = None
263 263 prePulserAfter = None
264 264 rangeIpp = None
265 265 rangeTxA = None
266 266 rangeTxB = None
267 267
268 268 __size = None
269 269
270 270 def __init__(self, expType=2, nTx=1,
271 271 ippKm=None, txA=0, txB=0,
272 272 nWindows=None, nHeights=None, firstHeight=None, deltaHeight=None,
273 273 numTaus=0, line6Function=0, line5Function=0, fClock=None,
274 274 prePulseBefore=0, prePulseAfter=0,
275 275 codeType=0, nCode=0, nBaud=0, code=None,
276 276 flip1=0, flip2=0):
277 277
278 278 # self.size = 116
279 279 self.expType = expType
280 280 self.nTx = nTx
281 281 self.ipp = ippKm
282 282 self.txA = txA
283 283 self.txB = txB
284 284 self.rangeIpp = ippKm
285 285 self.rangeTxA = txA
286 286 self.rangeTxB = txB
287 287
288 288 self.nWindows = nWindows
289 289 self.numTaus = numTaus
290 290 self.codeType = codeType
291 291 self.line6Function = line6Function
292 292 self.line5Function = line5Function
293 293 self.fClock = fClock
294 294 self.prePulseBefore = prePulseBefore
295 295 self.prePulserAfter = prePulseAfter
296 296
297 297 self.nHeights = nHeights
298 298 self.firstHeight = firstHeight
299 299 self.deltaHeight = deltaHeight
300 300 self.samplesWin = nHeights
301 301
302 302 self.nCode = nCode
303 303 self.nBaud = nBaud
304 304 self.code = code
305 305 self.flip1 = flip1
306 306 self.flip2 = flip2
307 307
308 308 self.code_size = int(numpy.ceil(self.nBaud/32.))*self.nCode*4
309 309 # self.dynamic = numpy.array([],numpy.dtype('byte'))
310 310
311 311 if self.fClock is None and self.deltaHeight is not None:
312 312 self.fClock = 0.15/(deltaHeight*1e-6) #0.15Km / (height * 1u)
313 313
314 314 def read(self, fp):
315 315 self.length = 0
316 316 try:
317 317 startFp = fp.tell()
318 318 except Exception, e:
319 319 startFp = None
320 320 pass
321 321
322 322 try:
323 323 if hasattr(fp, 'read'):
324 324 header = numpy.fromfile(fp, RADAR_STRUCTURE,1)
325 325 else:
326 326 header = numpy.fromstring(fp, RADAR_STRUCTURE,1)
327 327 self.length += header.nbytes
328 328 except Exception, e:
329 329 print "RadarControllerHeader: " + str(e)
330 330 return 0
331 331
332 332 size = int(header['nSize'][0])
333 333 self.expType = int(header['nExpType'][0])
334 334 self.nTx = int(header['nNTx'][0])
335 335 self.ipp = float(header['fIpp'][0])
336 336 self.txA = float(header['fTxA'][0])
337 337 self.txB = float(header['fTxB'][0])
338 338 self.nWindows = int(header['nNumWindows'][0])
339 339 self.numTaus = int(header['nNumTaus'][0])
340 340 self.codeType = int(header['nCodeType'][0])
341 341 self.line6Function = int(header['nLine6Function'][0])
342 342 self.line5Function = int(header['nLine5Function'][0])
343 343 self.fClock = float(header['fClock'][0])
344 344 self.prePulseBefore = int(header['nPrePulseBefore'][0])
345 345 self.prePulserAfter = int(header['nPrePulseAfter'][0])
346 346 self.rangeIpp = header['sRangeIPP'][0]
347 347 self.rangeTxA = header['sRangeTxA'][0]
348 348 self.rangeTxB = header['sRangeTxB'][0]
349 349
350 350 try:
351 351 if hasattr(fp, 'read'):
352 352 samplingWindow = numpy.fromfile(fp, SAMPLING_STRUCTURE, self.nWindows)
353 353 else:
354 354 samplingWindow = numpy.fromstring(fp[self.length:], SAMPLING_STRUCTURE, self.nWindows)
355 355 self.length += samplingWindow.nbytes
356 356 except Exception, e:
357 357 print "RadarControllerHeader: " + str(e)
358 358 return 0
359 359 self.nHeights = int(numpy.sum(samplingWindow['nsa']))
360 360 self.firstHeight = samplingWindow['h0']
361 361 self.deltaHeight = samplingWindow['dh']
362 362 self.samplesWin = samplingWindow['nsa']
363 363
364 364
365 365
366 366 try:
367 367 if hasattr(fp, 'read'):
368 368 self.Taus = numpy.fromfile(fp, '<f4', self.numTaus)
369 369 else:
370 370 self.Taus = numpy.fromstring(fp[self.length:], '<f4', self.numTaus)
371 371 self.length += self.Taus.nbytes
372 372 except Exception, e:
373 373 print "RadarControllerHeader: " + str(e)
374 374 return 0
375 375
376 376
377 377
378 378 self.code_size = 0
379 379 if self.codeType != 0:
380 380
381 381 try:
382 382 if hasattr(fp, 'read'):
383 self.nCode = numpy.fromfile(fp, '<u4', 1)
383 self.nCode = numpy.fromfile(fp, '<u4', 1)[0]
384 384 self.length += self.nCode.nbytes
385 self.nBaud = numpy.fromfile(fp, '<u4', 1)
385 self.nBaud = numpy.fromfile(fp, '<u4', 1)[0]
386 386 self.length += self.nBaud.nbytes
387 387 else:
388 388 self.nCode = numpy.fromstring(fp[self.length:], '<u4', 1)[0]
389 389 self.length += self.nCode.nbytes
390 390 self.nBaud = numpy.fromstring(fp[self.length:], '<u4', 1)[0]
391 391 self.length += self.nBaud.nbytes
392 392 except Exception, e:
393 393 print "RadarControllerHeader: " + str(e)
394 394 return 0
395 395 code = numpy.empty([self.nCode,self.nBaud],dtype='i1')
396 396
397 397 for ic in range(self.nCode):
398 398 try:
399 399 if hasattr(fp, 'read'):
400 400 temp = numpy.fromfile(fp,'u4', int(numpy.ceil(self.nBaud/32.)))
401 401 else:
402 402 temp = numpy.fromstring(fp,'u4', int(numpy.ceil(self.nBaud/32.)))
403 403 self.length += temp.nbytes
404 404 except Exception, e:
405 405 print "RadarControllerHeader: " + str(e)
406 406 return 0
407 407
408 408 for ib in range(self.nBaud-1,-1,-1):
409 409 code[ic,ib] = temp[ib/32]%2
410 410 temp[ib/32] = temp[ib/32]/2
411 411
412 412 self.code = 2.0*code - 1.0
413 413 self.code_size = int(numpy.ceil(self.nBaud/32.))*self.nCode*4
414 414
415 415 # if self.line5Function == RCfunction.FLIP:
416 416 # self.flip1 = numpy.fromfile(fp,'<u4',1)
417 417 #
418 418 # if self.line6Function == RCfunction.FLIP:
419 419 # self.flip2 = numpy.fromfile(fp,'<u4',1)
420 420 if startFp is not None:
421 421 endFp = size + startFp
422 422
423 423 if fp.tell() != endFp:
424 424 # fp.seek(endFp)
425 425 print "%s: Radar Controller Header size is not consistent: from data [%d] != from header field [%d]" %(fp.name, fp.tell()-startFp, size)
426 426 # return 0
427 427
428 428 if fp.tell() > endFp:
429 429 sys.stderr.write("Warning %s: Size value read from Radar Controller header is lower than it has to be\n" %fp.name)
430 430 # return 0
431 431
432 432 if fp.tell() < endFp:
433 433 sys.stderr.write("Warning %s: Size value read from Radar Controller header is greater than it has to be\n" %fp.name)
434 434
435 435
436 436 return 1
437 437
438 438 def write(self, fp):
439 439
440 440 headerTuple = (self.size,
441 441 self.expType,
442 442 self.nTx,
443 443 self.ipp,
444 444 self.txA,
445 445 self.txB,
446 446 self.nWindows,
447 447 self.numTaus,
448 448 self.codeType,
449 449 self.line6Function,
450 450 self.line5Function,
451 451 self.fClock,
452 452 self.prePulseBefore,
453 453 self.prePulserAfter,
454 454 self.rangeIpp,
455 455 self.rangeTxA,
456 456 self.rangeTxB)
457 457
458 458 header = numpy.array(headerTuple,RADAR_STRUCTURE)
459 459 header.tofile(fp)
460 460
461 461 sampleWindowTuple = (self.firstHeight,self.deltaHeight,self.samplesWin)
462 462 samplingWindow = numpy.array(sampleWindowTuple,SAMPLING_STRUCTURE)
463 463 samplingWindow.tofile(fp)
464 464
465 465 if self.numTaus > 0:
466 466 self.Taus.tofile(fp)
467 467
468 468 if self.codeType !=0:
469 469 nCode = numpy.array(self.nCode, '<u4')
470 470 nCode.tofile(fp)
471 471 nBaud = numpy.array(self.nBaud, '<u4')
472 472 nBaud.tofile(fp)
473 473 code1 = (self.code + 1.0)/2.
474 474
475 475 for ic in range(self.nCode):
476 476 tempx = numpy.zeros(numpy.ceil(self.nBaud/32.))
477 477 start = 0
478 478 end = 32
479 479 for i in range(len(tempx)):
480 480 code_selected = code1[ic,start:end]
481 481 for j in range(len(code_selected)-1,-1,-1):
482 482 if code_selected[j] == 1:
483 483 tempx[i] = tempx[i] + 2**(len(code_selected)-1-j)
484 484 start = start + 32
485 485 end = end + 32
486 486
487 487 tempx = tempx.astype('u4')
488 488 tempx.tofile(fp)
489 489
490 490 # if self.line5Function == RCfunction.FLIP:
491 491 # self.flip1.tofile(fp)
492 492 #
493 493 # if self.line6Function == RCfunction.FLIP:
494 494 # self.flip2.tofile(fp)
495 495
496 496 return 1
497 497
498 498 def get_ippSeconds(self):
499 499 '''
500 500 '''
501 501 ippSeconds = 2.0 * 1000 * self.ipp / SPEED_OF_LIGHT
502 502
503 503 return ippSeconds
504 504
505 505 def set_ippSeconds(self, ippSeconds):
506 506 '''
507 507 '''
508 508
509 509 self.ipp = ippSeconds * SPEED_OF_LIGHT / (2.0*1000)
510 510
511 511 return
512 512
513 513 def get_size(self):
514 514
515 515 self.__size = 116 + 12*self.nWindows + 4*self.numTaus
516 516
517 517 if self.codeType != 0:
518 518 self.__size += 4 + 4 + 4*self.nCode*numpy.ceil(self.nBaud/32.)
519 519
520 520 return self.__size
521 521
522 522 def set_size(self, value):
523 523
524 524 raise IOError, "size is a property and it cannot be set, just read"
525 525
526 526 return
527 527
528 528 ippSeconds = property(get_ippSeconds, set_ippSeconds)
529 529 size = property(get_size, set_size)
530 530
531 531 class ProcessingHeader(Header):
532 532
533 533 # size = None
534 534 dtype = None
535 535 blockSize = None
536 536 profilesPerBlock = None
537 537 dataBlocksPerFile = None
538 538 nWindows = None
539 539 processFlags = None
540 540 nCohInt = None
541 541 nIncohInt = None
542 542 totalSpectra = None
543 543
544 544 flag_dc = None
545 545 flag_cspc = None
546 546
547 547 def __init__(self):
548 548
549 549 # self.size = 0
550 550 self.dtype = 0
551 551 self.blockSize = 0
552 552 self.profilesPerBlock = 0
553 553 self.dataBlocksPerFile = 0
554 554 self.nWindows = 0
555 555 self.processFlags = 0
556 556 self.nCohInt = 0
557 557 self.nIncohInt = 0
558 558 self.totalSpectra = 0
559 559
560 560 self.nHeights = 0
561 561 self.firstHeight = 0
562 562 self.deltaHeight = 0
563 563 self.samplesWin = 0
564 564 self.spectraComb = 0
565 565 self.nCode = None
566 566 self.code = None
567 567 self.nBaud = None
568 568
569 569 self.shif_fft = False
570 570 self.flag_dc = False
571 571 self.flag_cspc = False
572 572 self.flag_decode = False
573 573 self.flag_deflip = False
574 574 self.length = 0
575 575 def read(self, fp):
576 576 self.length = 0
577 577 try:
578 578 startFp = fp.tell()
579 579 except Exception, e:
580 580 startFp = None
581 581 pass
582 582
583 583 try:
584 584 if hasattr(fp, 'read'):
585 585 header = numpy.fromfile(fp, PROCESSING_STRUCTURE, 1)
586 586 else:
587 587 header = numpy.fromstring(fp, PROCESSING_STRUCTURE, 1)
588 588 self.length += header.nbytes
589 589 except Exception, e:
590 590 print "ProcessingHeader: " + str(e)
591 591 return 0
592 592
593 593 size = int(header['nSize'][0])
594 594 self.dtype = int(header['nDataType'][0])
595 595 self.blockSize = int(header['nSizeOfDataBlock'][0])
596 596 self.profilesPerBlock = int(header['nProfilesperBlock'][0])
597 597 self.dataBlocksPerFile = int(header['nDataBlocksperFile'][0])
598 598 self.nWindows = int(header['nNumWindows'][0])
599 599 self.processFlags = header['nProcessFlags']
600 600 self.nCohInt = int(header['nCoherentIntegrations'][0])
601 601 self.nIncohInt = int(header['nIncoherentIntegrations'][0])
602 602 self.totalSpectra = int(header['nTotalSpectra'][0])
603 603
604 604 try:
605 605 if hasattr(fp, 'read'):
606 606 samplingWindow = numpy.fromfile(fp, SAMPLING_STRUCTURE, self.nWindows)
607 607 else:
608 608 samplingWindow = numpy.fromstring(fp[self.length:], SAMPLING_STRUCTURE, self.nWindows)
609 609 self.length += samplingWindow.nbytes
610 610 except Exception, e:
611 611 print "ProcessingHeader: " + str(e)
612 612 return 0
613 613
614 614 self.nHeights = int(numpy.sum(samplingWindow['nsa']))
615 615 self.firstHeight = float(samplingWindow['h0'][0])
616 616 self.deltaHeight = float(samplingWindow['dh'][0])
617 617 self.samplesWin = samplingWindow['nsa'][0]
618 618
619 619
620 620 try:
621 621 if hasattr(fp, 'read'):
622 622 self.spectraComb = numpy.fromfile(fp, 'u1', 2*self.totalSpectra)
623 623 else:
624 624 self.spectraComb = numpy.fromstring(fp[self.length:], 'u1', 2*self.totalSpectra)
625 625 self.length += self.spectraComb.nbytes
626 626 except Exception, e:
627 627 print "ProcessingHeader: " + str(e)
628 628 return 0
629 629
630 630 if ((self.processFlags & PROCFLAG.DEFINE_PROCESS_CODE) == PROCFLAG.DEFINE_PROCESS_CODE):
631 631 self.nCode = int(numpy.fromfile(fp,'<u4',1))
632 632 self.nBaud = int(numpy.fromfile(fp,'<u4',1))
633 633 self.code = numpy.fromfile(fp,'<f4',self.nCode*self.nBaud).reshape(self.nCode,self.nBaud)
634 634
635 635 if ((self.processFlags & PROCFLAG.EXP_NAME_ESP) == PROCFLAG.EXP_NAME_ESP):
636 636 exp_name_len = int(numpy.fromfile(fp,'<u4',1))
637 637 exp_name = numpy.fromfile(fp,'u1',exp_name_len+1)
638 638
639 639 if ((self.processFlags & PROCFLAG.SHIFT_FFT_DATA) == PROCFLAG.SHIFT_FFT_DATA):
640 640 self.shif_fft = True
641 641 else:
642 642 self.shif_fft = False
643 643
644 644 if ((self.processFlags & PROCFLAG.SAVE_CHANNELS_DC) == PROCFLAG.SAVE_CHANNELS_DC):
645 645 self.flag_dc = True
646 646 else:
647 647 self.flag_dc = False
648 648
649 649 if ((self.processFlags & PROCFLAG.DECODE_DATA) == PROCFLAG.DECODE_DATA):
650 650 self.flag_decode = True
651 651 else:
652 652 self.flag_decode = False
653 653
654 654 if ((self.processFlags & PROCFLAG.DEFLIP_DATA) == PROCFLAG.DEFLIP_DATA):
655 655 self.flag_deflip = True
656 656 else:
657 657 self.flag_deflip = False
658 658
659 659 nChannels = 0
660 660 nPairs = 0
661 661 pairList = []
662 662
663 663 for i in range( 0, self.totalSpectra*2, 2 ):
664 664 if self.spectraComb[i] == self.spectraComb[i+1]:
665 665 nChannels = nChannels + 1 #par de canales iguales
666 666 else:
667 667 nPairs = nPairs + 1 #par de canales diferentes
668 668 pairList.append( (self.spectraComb[i], self.spectraComb[i+1]) )
669 669
670 670 self.flag_cspc = False
671 671 if nPairs > 0:
672 672 self.flag_cspc = True
673 673
674 674
675 675
676 676 if startFp is not None:
677 677 endFp = size + startFp
678 678 if fp.tell() > endFp:
679 679 sys.stderr.write("Warning: Processing header size is lower than it has to be")
680 680 return 0
681 681
682 682 if fp.tell() < endFp:
683 683 sys.stderr.write("Warning: Processing header size is greater than it is considered")
684 684
685 685 return 1
686 686
687 687 def write(self, fp):
688 688 #Clear DEFINE_PROCESS_CODE
689 689 self.processFlags = self.processFlags & (~PROCFLAG.DEFINE_PROCESS_CODE)
690 690
691 691 headerTuple = (self.size,
692 692 self.dtype,
693 693 self.blockSize,
694 694 self.profilesPerBlock,
695 695 self.dataBlocksPerFile,
696 696 self.nWindows,
697 697 self.processFlags,
698 698 self.nCohInt,
699 699 self.nIncohInt,
700 700 self.totalSpectra)
701 701
702 702 header = numpy.array(headerTuple,PROCESSING_STRUCTURE)
703 703 header.tofile(fp)
704 704
705 705 if self.nWindows != 0:
706 706 sampleWindowTuple = (self.firstHeight,self.deltaHeight,self.samplesWin)
707 707 samplingWindow = numpy.array(sampleWindowTuple,SAMPLING_STRUCTURE)
708 708 samplingWindow.tofile(fp)
709 709
710 710 if self.totalSpectra != 0:
711 711 # spectraComb = numpy.array([],numpy.dtype('u1'))
712 712 spectraComb = self.spectraComb
713 713 spectraComb.tofile(fp)
714 714
715 715 # if self.processFlags & PROCFLAG.DEFINE_PROCESS_CODE == PROCFLAG.DEFINE_PROCESS_CODE:
716 716 # nCode = numpy.array([self.nCode], numpy.dtype('u4')) #Probar con un dato que almacene codigo, hasta el momento no se hizo la prueba
717 717 # nCode.tofile(fp)
718 718 #
719 719 # nBaud = numpy.array([self.nBaud], numpy.dtype('u4'))
720 720 # nBaud.tofile(fp)
721 721 #
722 722 # code = self.code.reshape(self.nCode*self.nBaud)
723 723 # code = code.astype(numpy.dtype('<f4'))
724 724 # code.tofile(fp)
725 725
726 726 return 1
727 727
728 728 def get_size(self):
729 729
730 730 self.__size = 40 + 12*self.nWindows + 2*self.totalSpectra
731 731
732 732 # if self.processFlags & PROCFLAG.DEFINE_PROCESS_CODE == PROCFLAG.DEFINE_PROCESS_CODE:
733 733 # self.__size += 4 + 4 + 4*self.nCode*numpy.ceil(self.nBaud/32.)
734 734 # self.__size += 4 + 4 + 4 * self.nCode * self.nBaud
735 735
736 736 return self.__size
737 737
738 738 def set_size(self, value):
739 739
740 740 raise IOError, "size is a property and it cannot be set, just read"
741 741
742 742 return
743 743
744 744 size = property(get_size, set_size)
745 745
746 746 class RCfunction:
747 747 NONE=0
748 748 FLIP=1
749 749 CODE=2
750 750 SAMPLING=3
751 751 LIN6DIV256=4
752 752 SYNCHRO=5
753 753
754 754 class nCodeType:
755 755 NONE=0
756 756 USERDEFINE=1
757 757 BARKER2=2
758 758 BARKER3=3
759 759 BARKER4=4
760 760 BARKER5=5
761 761 BARKER7=6
762 762 BARKER11=7
763 763 BARKER13=8
764 764 AC128=9
765 765 COMPLEMENTARYCODE2=10
766 766 COMPLEMENTARYCODE4=11
767 767 COMPLEMENTARYCODE8=12
768 768 COMPLEMENTARYCODE16=13
769 769 COMPLEMENTARYCODE32=14
770 770 COMPLEMENTARYCODE64=15
771 771 COMPLEMENTARYCODE128=16
772 772 CODE_BINARY28=17
773 773
774 774 class PROCFLAG:
775 775
776 776 COHERENT_INTEGRATION = numpy.uint32(0x00000001)
777 777 DECODE_DATA = numpy.uint32(0x00000002)
778 778 SPECTRA_CALC = numpy.uint32(0x00000004)
779 779 INCOHERENT_INTEGRATION = numpy.uint32(0x00000008)
780 780 POST_COHERENT_INTEGRATION = numpy.uint32(0x00000010)
781 781 SHIFT_FFT_DATA = numpy.uint32(0x00000020)
782 782
783 783 DATATYPE_CHAR = numpy.uint32(0x00000040)
784 784 DATATYPE_SHORT = numpy.uint32(0x00000080)
785 785 DATATYPE_LONG = numpy.uint32(0x00000100)
786 786 DATATYPE_INT64 = numpy.uint32(0x00000200)
787 787 DATATYPE_FLOAT = numpy.uint32(0x00000400)
788 788 DATATYPE_DOUBLE = numpy.uint32(0x00000800)
789 789
790 790 DATAARRANGE_CONTIGUOUS_CH = numpy.uint32(0x00001000)
791 791 DATAARRANGE_CONTIGUOUS_H = numpy.uint32(0x00002000)
792 792 DATAARRANGE_CONTIGUOUS_P = numpy.uint32(0x00004000)
793 793
794 794 SAVE_CHANNELS_DC = numpy.uint32(0x00008000)
795 795 DEFLIP_DATA = numpy.uint32(0x00010000)
796 796 DEFINE_PROCESS_CODE = numpy.uint32(0x00020000)
797 797
798 798 ACQ_SYS_NATALIA = numpy.uint32(0x00040000)
799 799 ACQ_SYS_ECHOTEK = numpy.uint32(0x00080000)
800 800 ACQ_SYS_ADRXD = numpy.uint32(0x000C0000)
801 801 ACQ_SYS_JULIA = numpy.uint32(0x00100000)
802 802 ACQ_SYS_XXXXXX = numpy.uint32(0x00140000)
803 803
804 804 EXP_NAME_ESP = numpy.uint32(0x00200000)
805 805 CHANNEL_NAMES_ESP = numpy.uint32(0x00400000)
806 806
807 807 OPERATION_MASK = numpy.uint32(0x0000003F)
808 808 DATATYPE_MASK = numpy.uint32(0x00000FC0)
809 809 DATAARRANGE_MASK = numpy.uint32(0x00007000)
810 810 ACQ_SYS_MASK = numpy.uint32(0x001C0000)
811 811
812 812 dtype0 = numpy.dtype([('real','<i1'),('imag','<i1')])
813 813 dtype1 = numpy.dtype([('real','<i2'),('imag','<i2')])
814 814 dtype2 = numpy.dtype([('real','<i4'),('imag','<i4')])
815 815 dtype3 = numpy.dtype([('real','<i8'),('imag','<i8')])
816 816 dtype4 = numpy.dtype([('real','<f4'),('imag','<f4')])
817 817 dtype5 = numpy.dtype([('real','<f8'),('imag','<f8')])
818 818
819 819 NUMPY_DTYPE_LIST = [dtype0, dtype1, dtype2, dtype3, dtype4, dtype5]
820 820
821 821 PROCFLAG_DTYPE_LIST = [PROCFLAG.DATATYPE_CHAR,
822 822 PROCFLAG.DATATYPE_SHORT,
823 823 PROCFLAG.DATATYPE_LONG,
824 824 PROCFLAG.DATATYPE_INT64,
825 825 PROCFLAG.DATATYPE_FLOAT,
826 826 PROCFLAG.DATATYPE_DOUBLE]
827 827
828 828 DTYPE_WIDTH = [1, 2, 4, 8, 4, 8]
829 829
830 830 def get_dtype_index(numpy_dtype):
831 831
832 832 index = None
833 833
834 834 for i in range(len(NUMPY_DTYPE_LIST)):
835 835 if numpy_dtype == NUMPY_DTYPE_LIST[i]:
836 836 index = i
837 837 break
838 838
839 839 return index
840 840
841 841 def get_numpy_dtype(index):
842 842
843 843 return NUMPY_DTYPE_LIST[index]
844 844
845 845 def get_procflag_dtype(index):
846 846
847 847 return PROCFLAG_DTYPE_LIST[index]
848 848
849 849 def get_dtype_width(index):
850 850
851 851 return DTYPE_WIDTH[index] No newline at end of file
@@ -1,2805 +1,2805
1 1 import numpy
2 2 import math
3 3 from scipy import optimize, interpolate, signal, stats, ndimage
4 4 import re
5 5 import datetime
6 6 import copy
7 7 import sys
8 8 import importlib
9 9 import itertools
10 10
11 11 from jroproc_base import ProcessingUnit, Operation
12 12 from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
13 13
14 14
15 15 class ParametersProc(ProcessingUnit):
16 16
17 17 nSeconds = None
18 18
19 19 def __init__(self):
20 20 ProcessingUnit.__init__(self)
21 21
22 22 # self.objectDict = {}
23 23 self.buffer = None
24 24 self.firstdatatime = None
25 25 self.profIndex = 0
26 26 self.dataOut = Parameters()
27 27
28 28 def __updateObjFromInput(self):
29 29
30 30 self.dataOut.inputUnit = self.dataIn.type
31 31
32 32 self.dataOut.timeZone = self.dataIn.timeZone
33 33 self.dataOut.dstFlag = self.dataIn.dstFlag
34 34 self.dataOut.errorCount = self.dataIn.errorCount
35 35 self.dataOut.useLocalTime = self.dataIn.useLocalTime
36 36
37 37 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
38 38 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
39 39 self.dataOut.channelList = self.dataIn.channelList
40 40 self.dataOut.heightList = self.dataIn.heightList
41 41 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
42 42 # self.dataOut.nHeights = self.dataIn.nHeights
43 43 # self.dataOut.nChannels = self.dataIn.nChannels
44 44 self.dataOut.nBaud = self.dataIn.nBaud
45 45 self.dataOut.nCode = self.dataIn.nCode
46 46 self.dataOut.code = self.dataIn.code
47 47 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
48 48 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
49 49 # self.dataOut.utctime = self.firstdatatime
50 50 self.dataOut.utctime = self.dataIn.utctime
51 51 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
52 52 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
53 53 self.dataOut.nCohInt = self.dataIn.nCohInt
54 54 # self.dataOut.nIncohInt = 1
55 55 self.dataOut.ippSeconds = self.dataIn.ippSeconds
56 56 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
57 57 self.dataOut.timeInterval1 = self.dataIn.timeInterval
58 58 self.dataOut.heightList = self.dataIn.getHeiRange()
59 59 self.dataOut.frequency = self.dataIn.frequency
60 60 #self.dataOut.noise = self.dataIn.noise
61 61
62 62 def run(self):
63 63
64 64 #---------------------- Voltage Data ---------------------------
65 65
66 66 if self.dataIn.type == "Voltage":
67 67
68 68 self.__updateObjFromInput()
69 69 self.dataOut.data_pre = self.dataIn.data.copy()
70 70 self.dataOut.flagNoData = False
71 71 self.dataOut.utctimeInit = self.dataIn.utctime
72 72 self.dataOut.paramInterval = self.dataIn.nProfiles*self.dataIn.nCohInt*self.dataIn.ippSeconds
73 73 return
74 74
75 75 #---------------------- Spectra Data ---------------------------
76 76
77 77 if self.dataIn.type == "Spectra":
78 78
79 79 self.dataOut.data_pre = (self.dataIn.data_spc, self.dataIn.data_cspc)
80 80 self.dataOut.data_spc = self.dataIn.data_spc
81 81 self.dataOut.data_cspc = self.dataIn.data_cspc
82 82 self.dataOut.nProfiles = self.dataIn.nProfiles
83 83 self.dataOut.nIncohInt = self.dataIn.nIncohInt
84 84 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
85 85 self.dataOut.ippFactor = self.dataIn.ippFactor
86 86 #self.dataOut.normFactor = self.dataIn.getNormFactor()
87 87 self.dataOut.pairsList = self.dataIn.pairsList
88 88 self.dataOut.groupList = self.dataIn.pairsList
89 89 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
90 90 self.dataOut.flagNoData = False
91 91
92 92 #---------------------- Correlation Data ---------------------------
93 93
94 94 if self.dataIn.type == "Correlation":
95 95 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.dataIn.splitFunctions()
96 96
97 97 self.dataOut.data_pre = (self.dataIn.data_cf[acf_ind,:], self.dataIn.data_cf[ccf_ind,:,:])
98 98 self.dataOut.normFactor = (self.dataIn.normFactor[acf_ind,:], self.dataIn.normFactor[ccf_ind,:])
99 99 self.dataOut.groupList = (acf_pairs, ccf_pairs)
100 100
101 101 self.dataOut.abscissaList = self.dataIn.lagRange
102 102 self.dataOut.noise = self.dataIn.noise
103 103 self.dataOut.data_SNR = self.dataIn.SNR
104 104 self.dataOut.flagNoData = False
105 105 self.dataOut.nAvg = self.dataIn.nAvg
106 106
107 107 #---------------------- Parameters Data ---------------------------
108 108
109 109 if self.dataIn.type == "Parameters":
110 110 self.dataOut.copy(self.dataIn)
111 111 self.dataOut.utctimeInit = self.dataIn.utctime
112 112 self.dataOut.flagNoData = False
113 113
114 114 return True
115 115
116 116 self.__updateObjFromInput()
117 117 self.dataOut.utctimeInit = self.dataIn.utctime
118 118 self.dataOut.paramInterval = self.dataIn.timeInterval
119 119
120 120 return
121 121
122 122 class SpectralMoments(Operation):
123 123
124 124 '''
125 125 Function SpectralMoments()
126 126
127 127 Calculates moments (power, mean, standard deviation) and SNR of the signal
128 128
129 129 Type of dataIn: Spectra
130 130
131 131 Configuration Parameters:
132 132
133 133 dirCosx : Cosine director in X axis
134 134 dirCosy : Cosine director in Y axis
135 135
136 136 elevation :
137 137 azimuth :
138 138
139 139 Input:
140 140 channelList : simple channel list to select e.g. [2,3,7]
141 141 self.dataOut.data_pre : Spectral data
142 142 self.dataOut.abscissaList : List of frequencies
143 143 self.dataOut.noise : Noise level per channel
144 144
145 145 Affected:
146 146 self.dataOut.data_param : Parameters per channel
147 147 self.dataOut.data_SNR : SNR per channel
148 148
149 149 '''
150 150
151 151 def run(self, dataOut):
152 152
153 153 #dataOut.data_pre = dataOut.data_pre[0]
154 154 data = dataOut.data_pre[0]
155 155 absc = dataOut.abscissaList[:-1]
156 156 noise = dataOut.noise
157 157 nChannel = data.shape[0]
158 158 data_param = numpy.zeros((nChannel, 4, data.shape[2]))
159 159
160 160 for ind in range(nChannel):
161 161 data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
162 162
163 163 dataOut.data_param = data_param[:,1:,:]
164 164 dataOut.data_SNR = data_param[:,0]
165 165 dataOut.data_DOP = data_param[:,1]
166 166 dataOut.data_MEAN = data_param[:,2]
167 167 dataOut.data_STD = data_param[:,3]
168 168 return
169 169
170 170 def __calculateMoments(self, oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
171 171
172 172 if (nicoh is None): nicoh = 1
173 173 if (graph is None): graph = 0
174 174 if (smooth is None): smooth = 0
175 175 elif (self.smooth < 3): smooth = 0
176 176
177 177 if (type1 is None): type1 = 0
178 178 if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
179 179 if (snrth is None): snrth = -3
180 180 if (dc is None): dc = 0
181 181 if (aliasing is None): aliasing = 0
182 182 if (oldfd is None): oldfd = 0
183 183 if (wwauto is None): wwauto = 0
184 184
185 185 if (n0 < 1.e-20): n0 = 1.e-20
186 186
187 187 freq = oldfreq
188 188 vec_power = numpy.zeros(oldspec.shape[1])
189 189 vec_fd = numpy.zeros(oldspec.shape[1])
190 190 vec_w = numpy.zeros(oldspec.shape[1])
191 191 vec_snr = numpy.zeros(oldspec.shape[1])
192 192
193 193 for ind in range(oldspec.shape[1]):
194 194
195 195 spec = oldspec[:,ind]
196 196 aux = spec*fwindow
197 197 max_spec = aux.max()
198 198 m = list(aux).index(max_spec)
199 199
200 200 #Smooth
201 201 if (smooth == 0): spec2 = spec
202 202 else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
203 203
204 204 # Calculo de Momentos
205 205 bb = spec2[range(m,spec2.size)]
206 206 bb = (bb<n0).nonzero()
207 207 bb = bb[0]
208 208
209 209 ss = spec2[range(0,m + 1)]
210 210 ss = (ss<n0).nonzero()
211 211 ss = ss[0]
212 212
213 213 if (bb.size == 0):
214 214 bb0 = spec.size - 1 - m
215 215 else:
216 216 bb0 = bb[0] - 1
217 217 if (bb0 < 0):
218 218 bb0 = 0
219 219
220 220 if (ss.size == 0): ss1 = 1
221 221 else: ss1 = max(ss) + 1
222 222
223 223 if (ss1 > m): ss1 = m
224 224
225 225 valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
226 226 power = ((spec2[valid] - n0)*fwindow[valid]).sum()
227 227 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
228 228 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
229 229 snr = (spec2.mean()-n0)/n0
230 230
231 231 if (snr < 1.e-20) :
232 232 snr = 1.e-20
233 233
234 234 vec_power[ind] = power
235 235 vec_fd[ind] = fd
236 236 vec_w[ind] = w
237 237 vec_snr[ind] = snr
238 238
239 239 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
240 240 return moments
241 241
242 242 #------------------ Get SA Parameters --------------------------
243 243
244 244 def GetSAParameters(self):
245 245 #SA en frecuencia
246 246 pairslist = self.dataOut.groupList
247 247 num_pairs = len(pairslist)
248 248
249 249 vel = self.dataOut.abscissaList
250 250 spectra = self.dataOut.data_pre[0]
251 251 cspectra = self.dataOut.data_pre[1]
252 252 delta_v = vel[1] - vel[0]
253 253
254 254 #Calculating the power spectrum
255 255 spc_pow = numpy.sum(spectra, 3)*delta_v
256 256 #Normalizing Spectra
257 257 norm_spectra = spectra/spc_pow
258 258 #Calculating the norm_spectra at peak
259 259 max_spectra = numpy.max(norm_spectra, 3)
260 260
261 261 #Normalizing Cross Spectra
262 262 norm_cspectra = numpy.zeros(cspectra.shape)
263 263
264 264 for i in range(num_chan):
265 265 norm_cspectra[i,:,:] = cspectra[i,:,:]/numpy.sqrt(spc_pow[pairslist[i][0],:]*spc_pow[pairslist[i][1],:])
266 266
267 267 max_cspectra = numpy.max(norm_cspectra,2)
268 268 max_cspectra_index = numpy.argmax(norm_cspectra, 2)
269 269
270 270 for i in range(num_pairs):
271 271 cspc_par[i,:,:] = __calculateMoments(norm_cspectra)
272 272 #------------------- Get Lags ----------------------------------
273 273
274 274 class SALags(Operation):
275 275 '''
276 276 Function GetMoments()
277 277
278 278 Input:
279 279 self.dataOut.data_pre
280 280 self.dataOut.abscissaList
281 281 self.dataOut.noise
282 282 self.dataOut.normFactor
283 283 self.dataOut.data_SNR
284 284 self.dataOut.groupList
285 285 self.dataOut.nChannels
286 286
287 287 Affected:
288 288 self.dataOut.data_param
289 289
290 290 '''
291 291 def run(self, dataOut):
292 292 data_acf = dataOut.data_pre[0]
293 293 data_ccf = dataOut.data_pre[1]
294 294 normFactor_acf = dataOut.normFactor[0]
295 295 normFactor_ccf = dataOut.normFactor[1]
296 296 pairs_acf = dataOut.groupList[0]
297 297 pairs_ccf = dataOut.groupList[1]
298 298
299 299 nHeights = dataOut.nHeights
300 300 absc = dataOut.abscissaList
301 301 noise = dataOut.noise
302 302 SNR = dataOut.data_SNR
303 303 nChannels = dataOut.nChannels
304 304 # pairsList = dataOut.groupList
305 305 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
306 306
307 307 for l in range(len(pairs_acf)):
308 308 data_acf[l,:,:] = data_acf[l,:,:]/normFactor_acf[l,:]
309 309
310 310 for l in range(len(pairs_ccf)):
311 311 data_ccf[l,:,:] = data_ccf[l,:,:]/normFactor_ccf[l,:]
312 312
313 313 dataOut.data_param = numpy.zeros((len(pairs_ccf)*2 + 1, nHeights))
314 314 dataOut.data_param[:-1,:] = self.__calculateTaus(data_acf, data_ccf, absc)
315 315 dataOut.data_param[-1,:] = self.__calculateLag1Phase(data_acf, absc)
316 316 return
317 317
318 318 # def __getPairsAutoCorr(self, pairsList, nChannels):
319 319 #
320 320 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
321 321 #
322 322 # for l in range(len(pairsList)):
323 323 # firstChannel = pairsList[l][0]
324 324 # secondChannel = pairsList[l][1]
325 325 #
326 326 # #Obteniendo pares de Autocorrelacion
327 327 # if firstChannel == secondChannel:
328 328 # pairsAutoCorr[firstChannel] = int(l)
329 329 #
330 330 # pairsAutoCorr = pairsAutoCorr.astype(int)
331 331 #
332 332 # pairsCrossCorr = range(len(pairsList))
333 333 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
334 334 #
335 335 # return pairsAutoCorr, pairsCrossCorr
336 336
337 337 def __calculateTaus(self, data_acf, data_ccf, lagRange):
338 338
339 339 lag0 = data_acf.shape[1]/2
340 340 #Funcion de Autocorrelacion
341 341 mean_acf = stats.nanmean(data_acf, axis = 0)
342 342
343 343 #Obtencion Indice de TauCross
344 344 ind_ccf = data_ccf.argmax(axis = 1)
345 345 #Obtencion Indice de TauAuto
346 346 ind_acf = numpy.zeros(ind_ccf.shape,dtype = 'int')
347 347 ccf_lag0 = data_ccf[:,lag0,:]
348 348
349 349 for i in range(ccf_lag0.shape[0]):
350 350 ind_acf[i,:] = numpy.abs(mean_acf - ccf_lag0[i,:]).argmin(axis = 0)
351 351
352 352 #Obtencion de TauCross y TauAuto
353 353 tau_ccf = lagRange[ind_ccf]
354 354 tau_acf = lagRange[ind_acf]
355 355
356 356 Nan1, Nan2 = numpy.where(tau_ccf == lagRange[0])
357 357
358 358 tau_ccf[Nan1,Nan2] = numpy.nan
359 359 tau_acf[Nan1,Nan2] = numpy.nan
360 360 tau = numpy.vstack((tau_ccf,tau_acf))
361 361
362 362 return tau
363 363
364 364 def __calculateLag1Phase(self, data, lagTRange):
365 365 data1 = stats.nanmean(data, axis = 0)
366 366 lag1 = numpy.where(lagTRange == 0)[0][0] + 1
367 367
368 368 phase = numpy.angle(data1[lag1,:])
369 369
370 370 return phase
371 371
372 372 class SpectralFitting(Operation):
373 373 '''
374 374 Function GetMoments()
375 375
376 376 Input:
377 377 Output:
378 378 Variables modified:
379 379 '''
380 380
381 381 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):
382 382
383 383
384 384 if path != None:
385 385 sys.path.append(path)
386 386 self.dataOut.library = importlib.import_module(file)
387 387
388 388 #To be inserted as a parameter
389 389 groupArray = numpy.array(groupList)
390 390 # groupArray = numpy.array([[0,1],[2,3]])
391 391 self.dataOut.groupList = groupArray
392 392
393 393 nGroups = groupArray.shape[0]
394 394 nChannels = self.dataIn.nChannels
395 395 nHeights=self.dataIn.heightList.size
396 396
397 397 #Parameters Array
398 398 self.dataOut.data_param = None
399 399
400 400 #Set constants
401 401 constants = self.dataOut.library.setConstants(self.dataIn)
402 402 self.dataOut.constants = constants
403 403 M = self.dataIn.normFactor
404 404 N = self.dataIn.nFFTPoints
405 405 ippSeconds = self.dataIn.ippSeconds
406 406 K = self.dataIn.nIncohInt
407 407 pairsArray = numpy.array(self.dataIn.pairsList)
408 408
409 409 #List of possible combinations
410 410 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
411 411 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
412 412
413 413 if getSNR:
414 414 listChannels = groupArray.reshape((groupArray.size))
415 415 listChannels.sort()
416 416 noise = self.dataIn.getNoise()
417 417 self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
418 418
419 419 for i in range(nGroups):
420 420 coord = groupArray[i,:]
421 421
422 422 #Input data array
423 423 data = self.dataIn.data_spc[coord,:,:]/(M*N)
424 424 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
425 425
426 426 #Cross Spectra data array for Covariance Matrixes
427 427 ind = 0
428 428 for pairs in listComb:
429 429 pairsSel = numpy.array([coord[x],coord[y]])
430 430 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
431 431 ind += 1
432 432 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
433 433 dataCross = dataCross**2/K
434 434
435 435 for h in range(nHeights):
436 436 # print self.dataOut.heightList[h]
437 437
438 438 #Input
439 439 d = data[:,h]
440 440
441 441 #Covariance Matrix
442 442 D = numpy.diag(d**2/K)
443 443 ind = 0
444 444 for pairs in listComb:
445 445 #Coordinates in Covariance Matrix
446 446 x = pairs[0]
447 447 y = pairs[1]
448 448 #Channel Index
449 449 S12 = dataCross[ind,:,h]
450 450 D12 = numpy.diag(S12)
451 451 #Completing Covariance Matrix with Cross Spectras
452 452 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
453 453 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
454 454 ind += 1
455 455 Dinv=numpy.linalg.inv(D)
456 456 L=numpy.linalg.cholesky(Dinv)
457 457 LT=L.T
458 458
459 459 dp = numpy.dot(LT,d)
460 460
461 461 #Initial values
462 462 data_spc = self.dataIn.data_spc[coord,:,h]
463 463
464 464 if (h>0)and(error1[3]<5):
465 465 p0 = self.dataOut.data_param[i,:,h-1]
466 466 else:
467 467 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
468 468
469 469 try:
470 470 #Least Squares
471 471 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
472 472 # minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
473 473 #Chi square error
474 474 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
475 475 #Error with Jacobian
476 476 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
477 477 except:
478 478 minp = p0*numpy.nan
479 479 error0 = numpy.nan
480 480 error1 = p0*numpy.nan
481 481
482 482 #Save
483 483 if self.dataOut.data_param is None:
484 484 self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
485 485 self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
486 486
487 487 self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
488 488 self.dataOut.data_param[i,:,h] = minp
489 489 return
490 490
491 491 def __residFunction(self, p, dp, LT, constants):
492 492
493 493 fm = self.dataOut.library.modelFunction(p, constants)
494 494 fmp=numpy.dot(LT,fm)
495 495
496 496 return dp-fmp
497 497
498 498 def __getSNR(self, z, noise):
499 499
500 500 avg = numpy.average(z, axis=1)
501 501 SNR = (avg.T-noise)/noise
502 502 SNR = SNR.T
503 503 return SNR
504 504
505 505 def __chisq(p,chindex,hindex):
506 506 #similar to Resid but calculates CHI**2
507 507 [LT,d,fm]=setupLTdfm(p,chindex,hindex)
508 508 dp=numpy.dot(LT,d)
509 509 fmp=numpy.dot(LT,fm)
510 510 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
511 511 return chisq
512 512
513 513 class WindProfiler(Operation):
514 514
515 515 __isConfig = False
516 516
517 517 __initime = None
518 518 __lastdatatime = None
519 519 __integrationtime = None
520 520
521 521 __buffer = None
522 522
523 523 __dataReady = False
524 524
525 525 __firstdata = None
526 526
527 527 n = None
528 528
529 529 def __calculateCosDir(self, elev, azim):
530 530 zen = (90 - elev)*numpy.pi/180
531 531 azim = azim*numpy.pi/180
532 532 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
533 533 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
534 534
535 535 signX = numpy.sign(numpy.cos(azim))
536 536 signY = numpy.sign(numpy.sin(azim))
537 537
538 538 cosDirX = numpy.copysign(cosDirX, signX)
539 539 cosDirY = numpy.copysign(cosDirY, signY)
540 540 return cosDirX, cosDirY
541 541
542 542 def __calculateAngles(self, theta_x, theta_y, azimuth):
543 543
544 544 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
545 545 zenith_arr = numpy.arccos(dir_cosw)
546 546 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
547 547
548 548 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
549 549 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
550 550
551 551 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
552 552
553 553 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
554 554
555 555 #
556 556 if horOnly:
557 557 A = numpy.c_[dir_cosu,dir_cosv]
558 558 else:
559 559 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
560 560 A = numpy.asmatrix(A)
561 561 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
562 562
563 563 return A1
564 564
565 565 def __correctValues(self, heiRang, phi, velRadial, SNR):
566 566 listPhi = phi.tolist()
567 567 maxid = listPhi.index(max(listPhi))
568 568 minid = listPhi.index(min(listPhi))
569 569
570 570 rango = range(len(phi))
571 571 # rango = numpy.delete(rango,maxid)
572 572
573 573 heiRang1 = heiRang*math.cos(phi[maxid])
574 574 heiRangAux = heiRang*math.cos(phi[minid])
575 575 indOut = (heiRang1 < heiRangAux[0]).nonzero()
576 576 heiRang1 = numpy.delete(heiRang1,indOut)
577 577
578 578 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
579 579 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
580 580
581 581 for i in rango:
582 582 x = heiRang*math.cos(phi[i])
583 583 y1 = velRadial[i,:]
584 584 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
585 585
586 586 x1 = heiRang1
587 587 y11 = f1(x1)
588 588
589 589 y2 = SNR[i,:]
590 590 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
591 591 y21 = f2(x1)
592 592
593 593 velRadial1[i,:] = y11
594 594 SNR1[i,:] = y21
595 595
596 596 return heiRang1, velRadial1, SNR1
597 597
598 598 def __calculateVelUVW(self, A, velRadial):
599 599
600 600 #Operacion Matricial
601 601 # velUVW = numpy.zeros((velRadial.shape[1],3))
602 602 # for ind in range(velRadial.shape[1]):
603 603 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
604 604 # velUVW = velUVW.transpose()
605 605 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
606 606 velUVW[:,:] = numpy.dot(A,velRadial)
607 607
608 608
609 609 return velUVW
610 610
611 611 # def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
612 612
613 613 def techniqueDBS(self, kwargs):
614 614 """
615 615 Function that implements Doppler Beam Swinging (DBS) technique.
616 616
617 617 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
618 618 Direction correction (if necessary), Ranges and SNR
619 619
620 620 Output: Winds estimation (Zonal, Meridional and Vertical)
621 621
622 622 Parameters affected: Winds, height range, SNR
623 623 """
624 624 velRadial0 = kwargs['velRadial']
625 625 heiRang = kwargs['heightList']
626 626 SNR0 = kwargs['SNR']
627 627
628 628 if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
629 629 theta_x = numpy.array(kwargs['dirCosx'])
630 630 theta_y = numpy.array(kwargs['dirCosy'])
631 631 else:
632 632 elev = numpy.array(kwargs['elevation'])
633 633 azim = numpy.array(kwargs['azimuth'])
634 634 theta_x, theta_y = self.__calculateCosDir(elev, azim)
635 635 azimuth = kwargs['correctAzimuth']
636 636 if kwargs.has_key('horizontalOnly'):
637 637 horizontalOnly = kwargs['horizontalOnly']
638 638 else: horizontalOnly = False
639 639 if kwargs.has_key('correctFactor'):
640 640 correctFactor = kwargs['correctFactor']
641 641 else: correctFactor = 1
642 642 if kwargs.has_key('channelList'):
643 643 channelList = kwargs['channelList']
644 644 if len(channelList) == 2:
645 645 horizontalOnly = True
646 646 arrayChannel = numpy.array(channelList)
647 647 param = param[arrayChannel,:,:]
648 648 theta_x = theta_x[arrayChannel]
649 649 theta_y = theta_y[arrayChannel]
650 650
651 651 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
652 652 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0)
653 653 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
654 654
655 655 #Calculo de Componentes de la velocidad con DBS
656 656 winds = self.__calculateVelUVW(A,velRadial1)
657 657
658 658 return winds, heiRang1, SNR1
659 659
660 660 def __calculateDistance(self, posx, posy, pairs_ccf, azimuth = None):
661 661
662 662 nPairs = len(pairs_ccf)
663 663 posx = numpy.asarray(posx)
664 664 posy = numpy.asarray(posy)
665 665
666 666 #Rotacion Inversa para alinear con el azimuth
667 667 if azimuth!= None:
668 668 azimuth = azimuth*math.pi/180
669 669 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
670 670 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
671 671 else:
672 672 posx1 = posx
673 673 posy1 = posy
674 674
675 675 #Calculo de Distancias
676 676 distx = numpy.zeros(nPairs)
677 677 disty = numpy.zeros(nPairs)
678 678 dist = numpy.zeros(nPairs)
679 679 ang = numpy.zeros(nPairs)
680 680
681 681 for i in range(nPairs):
682 682 distx[i] = posx1[pairs_ccf[i][1]] - posx1[pairs_ccf[i][0]]
683 683 disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]]
684 684 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
685 685 ang[i] = numpy.arctan2(disty[i],distx[i])
686 686
687 687 return distx, disty, dist, ang
688 688 #Calculo de Matrices
689 689 # nPairs = len(pairs)
690 690 # ang1 = numpy.zeros((nPairs, 2, 1))
691 691 # dist1 = numpy.zeros((nPairs, 2, 1))
692 692 #
693 693 # for j in range(nPairs):
694 694 # dist1[j,0,0] = dist[pairs[j][0]]
695 695 # dist1[j,1,0] = dist[pairs[j][1]]
696 696 # ang1[j,0,0] = ang[pairs[j][0]]
697 697 # ang1[j,1,0] = ang[pairs[j][1]]
698 698 #
699 699 # return distx,disty, dist1,ang1
700 700
701 701
702 702 def __calculateVelVer(self, phase, lagTRange, _lambda):
703 703
704 704 Ts = lagTRange[1] - lagTRange[0]
705 705 velW = -_lambda*phase/(4*math.pi*Ts)
706 706
707 707 return velW
708 708
709 709 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
710 710 nPairs = tau1.shape[0]
711 711 nHeights = tau1.shape[1]
712 712 vel = numpy.zeros((nPairs,3,nHeights))
713 713 dist1 = numpy.reshape(dist, (dist.size,1))
714 714
715 715 angCos = numpy.cos(ang)
716 716 angSin = numpy.sin(ang)
717 717
718 718 vel0 = dist1*tau1/(2*tau2**2)
719 719 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
720 720 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
721 721
722 722 ind = numpy.where(numpy.isinf(vel))
723 723 vel[ind] = numpy.nan
724 724
725 725 return vel
726 726
727 727 # def __getPairsAutoCorr(self, pairsList, nChannels):
728 728 #
729 729 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
730 730 #
731 731 # for l in range(len(pairsList)):
732 732 # firstChannel = pairsList[l][0]
733 733 # secondChannel = pairsList[l][1]
734 734 #
735 735 # #Obteniendo pares de Autocorrelacion
736 736 # if firstChannel == secondChannel:
737 737 # pairsAutoCorr[firstChannel] = int(l)
738 738 #
739 739 # pairsAutoCorr = pairsAutoCorr.astype(int)
740 740 #
741 741 # pairsCrossCorr = range(len(pairsList))
742 742 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
743 743 #
744 744 # return pairsAutoCorr, pairsCrossCorr
745 745
746 746 # def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
747 747 def techniqueSA(self, kwargs):
748 748
749 749 """
750 750 Function that implements Spaced Antenna (SA) technique.
751 751
752 752 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
753 753 Direction correction (if necessary), Ranges and SNR
754 754
755 755 Output: Winds estimation (Zonal, Meridional and Vertical)
756 756
757 757 Parameters affected: Winds
758 758 """
759 759 position_x = kwargs['positionX']
760 760 position_y = kwargs['positionY']
761 761 azimuth = kwargs['azimuth']
762 762
763 763 if kwargs.has_key('correctFactor'):
764 764 correctFactor = kwargs['correctFactor']
765 765 else:
766 766 correctFactor = 1
767 767
768 768 groupList = kwargs['groupList']
769 769 pairs_ccf = groupList[1]
770 770 tau = kwargs['tau']
771 771 _lambda = kwargs['_lambda']
772 772
773 773 #Cross Correlation pairs obtained
774 774 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairssList, nChannels)
775 775 # pairsArray = numpy.array(pairsList)[pairsCrossCorr]
776 776 # pairsSelArray = numpy.array(pairsSelected)
777 777 # pairs = []
778 778 #
779 779 # #Wind estimation pairs obtained
780 780 # for i in range(pairsSelArray.shape[0]/2):
781 781 # ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
782 782 # ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
783 783 # pairs.append((ind1,ind2))
784 784
785 785 indtau = tau.shape[0]/2
786 786 tau1 = tau[:indtau,:]
787 787 tau2 = tau[indtau:-1,:]
788 788 # tau1 = tau1[pairs,:]
789 789 # tau2 = tau2[pairs,:]
790 790 phase1 = tau[-1,:]
791 791
792 792 #---------------------------------------------------------------------
793 793 #Metodo Directo
794 794 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairs_ccf,azimuth)
795 795 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
796 796 winds = stats.nanmean(winds, axis=0)
797 797 #---------------------------------------------------------------------
798 798 #Metodo General
799 799 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
800 800 # #Calculo Coeficientes de Funcion de Correlacion
801 801 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
802 802 # #Calculo de Velocidades
803 803 # winds = self.calculateVelUV(F,G,A,B,H)
804 804
805 805 #---------------------------------------------------------------------
806 806 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
807 807 winds = correctFactor*winds
808 808 return winds
809 809
810 810 def __checkTime(self, currentTime, paramInterval, outputInterval):
811 811
812 812 dataTime = currentTime + paramInterval
813 813 deltaTime = dataTime - self.__initime
814 814
815 815 if deltaTime >= outputInterval or deltaTime < 0:
816 816 self.__dataReady = True
817 817 return
818 818
819 819 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax, binkm=2):
820 820 '''
821 821 Function that implements winds estimation technique with detected meteors.
822 822
823 823 Input: Detected meteors, Minimum meteor quantity to wind estimation
824 824
825 825 Output: Winds estimation (Zonal and Meridional)
826 826
827 827 Parameters affected: Winds
828 828 '''
829 829 # print arrayMeteor.shape
830 830 #Settings
831 831 nInt = (heightMax - heightMin)/binkm
832 832 # print nInt
833 833 nInt = int(nInt)
834 834 # print nInt
835 835 winds = numpy.zeros((2,nInt))*numpy.nan
836 836
837 837 #Filter errors
838 838 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
839 839 finalMeteor = arrayMeteor[error,:]
840 840
841 841 #Meteor Histogram
842 842 finalHeights = finalMeteor[:,2]
843 843 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
844 844 nMeteorsPerI = hist[0]
845 845 heightPerI = hist[1]
846 846
847 847 #Sort of meteors
848 848 indSort = finalHeights.argsort()
849 849 finalMeteor2 = finalMeteor[indSort,:]
850 850
851 851 # Calculating winds
852 852 ind1 = 0
853 853 ind2 = 0
854 854
855 855 for i in range(nInt):
856 856 nMet = nMeteorsPerI[i]
857 857 ind1 = ind2
858 858 ind2 = ind1 + nMet
859 859
860 860 meteorAux = finalMeteor2[ind1:ind2,:]
861 861
862 862 if meteorAux.shape[0] >= meteorThresh:
863 863 vel = meteorAux[:, 6]
864 864 zen = meteorAux[:, 4]*numpy.pi/180
865 865 azim = meteorAux[:, 3]*numpy.pi/180
866 866
867 867 n = numpy.cos(zen)
868 868 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
869 869 # l = m*numpy.tan(azim)
870 870 l = numpy.sin(zen)*numpy.sin(azim)
871 871 m = numpy.sin(zen)*numpy.cos(azim)
872 872
873 873 A = numpy.vstack((l, m)).transpose()
874 874 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
875 875 windsAux = numpy.dot(A1, vel)
876 876
877 877 winds[0,i] = windsAux[0]
878 878 winds[1,i] = windsAux[1]
879 879
880 880 return winds, heightPerI[:-1]
881 881
882 882 def techniqueNSM_SA(self, **kwargs):
883 883 metArray = kwargs['metArray']
884 884 heightList = kwargs['heightList']
885 885 timeList = kwargs['timeList']
886 886
887 887 rx_location = kwargs['rx_location']
888 888 groupList = kwargs['groupList']
889 889 azimuth = kwargs['azimuth']
890 890 dfactor = kwargs['dfactor']
891 891 k = kwargs['k']
892 892
893 893 azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)
894 894 d = dist*dfactor
895 895 #Phase calculation
896 896 metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)
897 897
898 898 metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities
899 899
900 900 velEst = numpy.zeros((heightList.size,2))*numpy.nan
901 901 azimuth1 = azimuth1*numpy.pi/180
902 902
903 903 for i in range(heightList.size):
904 904 h = heightList[i]
905 905 indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]
906 906 metHeight = metArray1[indH,:]
907 907 if metHeight.shape[0] >= 2:
908 908 velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities
909 909 iazim = metHeight[:,1].astype(int)
910 910 azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths
911 911 A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))
912 912 A = numpy.asmatrix(A)
913 913 A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose()
914 914 velHor = numpy.dot(A1,velAux)
915 915
916 916 velEst[i,:] = numpy.squeeze(velHor)
917 917 return velEst
918 918
919 919 def __getPhaseSlope(self, metArray, heightList, timeList):
920 920 meteorList = []
921 921 #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2
922 922 #Putting back together the meteor matrix
923 923 utctime = metArray[:,0]
924 924 uniqueTime = numpy.unique(utctime)
925 925
926 926 phaseDerThresh = 0.5
927 927 ippSeconds = timeList[1] - timeList[0]
928 928 sec = numpy.where(timeList>1)[0][0]
929 929 nPairs = metArray.shape[1] - 6
930 930 nHeights = len(heightList)
931 931
932 932 for t in uniqueTime:
933 933 metArray1 = metArray[utctime==t,:]
934 934 # phaseDerThresh = numpy.pi/4 #reducir Phase thresh
935 935 tmet = metArray1[:,1].astype(int)
936 936 hmet = metArray1[:,2].astype(int)
937 937
938 938 metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))
939 939 metPhase[:,:] = numpy.nan
940 940 metPhase[:,hmet,tmet] = metArray1[:,6:].T
941 941
942 942 #Delete short trails
943 943 metBool = ~numpy.isnan(metPhase[0,:,:])
944 944 heightVect = numpy.sum(metBool, axis = 1)
945 945 metBool[heightVect<sec,:] = False
946 946 metPhase[:,heightVect<sec,:] = numpy.nan
947 947
948 948 #Derivative
949 949 metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])
950 950 phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))
951 951 metPhase[phDerAux] = numpy.nan
952 952
953 953 #--------------------------METEOR DETECTION -----------------------------------------
954 954 indMet = numpy.where(numpy.any(metBool,axis=1))[0]
955 955
956 956 for p in numpy.arange(nPairs):
957 957 phase = metPhase[p,:,:]
958 958 phDer = metDer[p,:,:]
959 959
960 960 for h in indMet:
961 961 height = heightList[h]
962 962 phase1 = phase[h,:] #82
963 963 phDer1 = phDer[h,:]
964 964
965 965 phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap
966 966
967 967 indValid = numpy.where(~numpy.isnan(phase1))[0]
968 968 initMet = indValid[0]
969 969 endMet = 0
970 970
971 971 for i in range(len(indValid)-1):
972 972
973 973 #Time difference
974 974 inow = indValid[i]
975 975 inext = indValid[i+1]
976 976 idiff = inext - inow
977 977 #Phase difference
978 978 phDiff = numpy.abs(phase1[inext] - phase1[inow])
979 979
980 980 if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor
981 981 sizeTrail = inow - initMet + 1
982 982 if sizeTrail>3*sec: #Too short meteors
983 983 x = numpy.arange(initMet,inow+1)*ippSeconds
984 984 y = phase1[initMet:inow+1]
985 985 ynnan = ~numpy.isnan(y)
986 986 x = x[ynnan]
987 987 y = y[ynnan]
988 988 slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
989 989 ylin = x*slope + intercept
990 990 rsq = r_value**2
991 991 if rsq > 0.5:
992 992 vel = slope#*height*1000/(k*d)
993 993 estAux = numpy.array([utctime,p,height, vel, rsq])
994 994 meteorList.append(estAux)
995 995 initMet = inext
996 996 metArray2 = numpy.array(meteorList)
997 997
998 998 return metArray2
999 999
1000 1000 def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):
1001 1001
1002 1002 azimuth1 = numpy.zeros(len(pairslist))
1003 1003 dist = numpy.zeros(len(pairslist))
1004 1004
1005 1005 for i in range(len(rx_location)):
1006 1006 ch0 = pairslist[i][0]
1007 1007 ch1 = pairslist[i][1]
1008 1008
1009 1009 diffX = rx_location[ch0][0] - rx_location[ch1][0]
1010 1010 diffY = rx_location[ch0][1] - rx_location[ch1][1]
1011 1011 azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi
1012 1012 dist[i] = numpy.sqrt(diffX**2 + diffY**2)
1013 1013
1014 1014 azimuth1 -= azimuth0
1015 1015 return azimuth1, dist
1016 1016
1017 1017 def techniqueNSM_DBS(self, **kwargs):
1018 1018 metArray = kwargs['metArray']
1019 1019 heightList = kwargs['heightList']
1020 1020 timeList = kwargs['timeList']
1021 1021 azimuth = kwargs['azimuth']
1022 1022 theta_x = numpy.array(kwargs['theta_x'])
1023 1023 theta_y = numpy.array(kwargs['theta_y'])
1024 1024
1025 1025 utctime = metArray[:,0]
1026 1026 cmet = metArray[:,1].astype(int)
1027 1027 hmet = metArray[:,3].astype(int)
1028 1028 SNRmet = metArray[:,4]
1029 1029 vmet = metArray[:,5]
1030 1030 spcmet = metArray[:,6]
1031 1031
1032 1032 nChan = numpy.max(cmet) + 1
1033 1033 nHeights = len(heightList)
1034 1034
1035 1035 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
1036 1036 hmet = heightList[hmet]
1037 1037 h1met = hmet*numpy.cos(zenith_arr[cmet]) #Corrected heights
1038 1038
1039 1039 velEst = numpy.zeros((heightList.size,2))*numpy.nan
1040 1040
1041 1041 for i in range(nHeights - 1):
1042 1042 hmin = heightList[i]
1043 1043 hmax = heightList[i + 1]
1044 1044
1045 1045 thisH = (h1met>=hmin) & (h1met<hmax) & (cmet!=2) & (SNRmet>8) & (vmet<50) & (spcmet<10)
1046 1046 indthisH = numpy.where(thisH)
1047 1047
1048 1048 if numpy.size(indthisH) > 3:
1049 1049
1050 1050 vel_aux = vmet[thisH]
1051 1051 chan_aux = cmet[thisH]
1052 1052 cosu_aux = dir_cosu[chan_aux]
1053 1053 cosv_aux = dir_cosv[chan_aux]
1054 1054 cosw_aux = dir_cosw[chan_aux]
1055 1055
1056 1056 nch = numpy.size(numpy.unique(chan_aux))
1057 1057 if nch > 1:
1058 1058 A = self.__calculateMatA(cosu_aux, cosv_aux, cosw_aux, True)
1059 1059 velEst[i,:] = numpy.dot(A,vel_aux)
1060 1060
1061 1061 return velEst
1062 1062
1063 def run(self, dataOut, technique, **kwargs):
1063 def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs):
1064 1064
1065 1065 param = dataOut.data_param
1066 1066 if dataOut.abscissaList != None:
1067 1067 absc = dataOut.abscissaList[:-1]
1068 noise = dataOut.noise
1068 # noise = dataOut.noise
1069 1069 heightList = dataOut.heightList
1070 1070 SNR = dataOut.data_SNR
1071 1071
1072 1072 if technique == 'DBS':
1073 1073
1074 1074 kwargs['velRadial'] = param[:,1,:] #Radial velocity
1075 1075 kwargs['heightList'] = heightList
1076 1076 kwargs['SNR'] = SNR
1077 1077
1078 1078 dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(kwargs) #DBS Function
1079 1079 dataOut.utctimeInit = dataOut.utctime
1080 1080 dataOut.outputInterval = dataOut.paramInterval
1081 1081
1082 1082 elif technique == 'SA':
1083 1083
1084 1084 #Parameters
1085 1085 # position_x = kwargs['positionX']
1086 1086 # position_y = kwargs['positionY']
1087 1087 # azimuth = kwargs['azimuth']
1088 1088 #
1089 1089 # if kwargs.has_key('crosspairsList'):
1090 1090 # pairs = kwargs['crosspairsList']
1091 1091 # else:
1092 1092 # pairs = None
1093 1093 #
1094 1094 # if kwargs.has_key('correctFactor'):
1095 1095 # correctFactor = kwargs['correctFactor']
1096 1096 # else:
1097 1097 # correctFactor = 1
1098 1098
1099 1099 # tau = dataOut.data_param
1100 1100 # _lambda = dataOut.C/dataOut.frequency
1101 1101 # pairsList = dataOut.groupList
1102 1102 # nChannels = dataOut.nChannels
1103 1103
1104 1104 kwargs['groupList'] = dataOut.groupList
1105 1105 kwargs['tau'] = dataOut.data_param
1106 1106 kwargs['_lambda'] = dataOut.C/dataOut.frequency
1107 1107 # dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
1108 1108 dataOut.data_output = self.techniqueSA(kwargs)
1109 1109 dataOut.utctimeInit = dataOut.utctime
1110 1110 dataOut.outputInterval = dataOut.timeInterval
1111 1111
1112 1112 elif technique == 'Meteors':
1113 1113 dataOut.flagNoData = True
1114 1114 self.__dataReady = False
1115 1115
1116 1116 if kwargs.has_key('nHours'):
1117 1117 nHours = kwargs['nHours']
1118 1118 else:
1119 1119 nHours = 1
1120 1120
1121 1121 if kwargs.has_key('meteorsPerBin'):
1122 1122 meteorThresh = kwargs['meteorsPerBin']
1123 1123 else:
1124 1124 meteorThresh = 6
1125 1125
1126 1126 if kwargs.has_key('hmin'):
1127 1127 hmin = kwargs['hmin']
1128 1128 else: hmin = 70
1129 1129 if kwargs.has_key('hmax'):
1130 1130 hmax = kwargs['hmax']
1131 1131 else: hmax = 110
1132 1132
1133 1133 if kwargs.has_key('BinKm'):
1134 1134 binkm = kwargs['BinKm']
1135 1135 else:
1136 1136 binkm = 2
1137 1137
1138 1138 dataOut.outputInterval = nHours*3600
1139 1139
1140 1140 if self.__isConfig == False:
1141 1141 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1142 1142 #Get Initial LTC time
1143 1143 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
1144 1144 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1145 1145
1146 1146 self.__isConfig = True
1147 1147
1148 1148 if self.__buffer is None:
1149 1149 self.__buffer = dataOut.data_param
1150 1150 self.__firstdata = copy.copy(dataOut)
1151 1151
1152 1152 else:
1153 1153 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1154 1154
1155 1155 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1156 1156
1157 1157 if self.__dataReady:
1158 1158 dataOut.utctimeInit = self.__initime
1159 1159
1160 1160 self.__initime += dataOut.outputInterval #to erase time offset
1161 1161
1162 1162 dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax, binkm)
1163 1163 dataOut.flagNoData = False
1164 1164 self.__buffer = None
1165 1165
1166 1166 elif technique == 'Meteors1':
1167 1167 dataOut.flagNoData = True
1168 1168 self.__dataReady = False
1169 1169
1170 1170 if kwargs.has_key('nMins'):
1171 1171 nMins = kwargs['nMins']
1172 1172 else: nMins = 20
1173 1173 if kwargs.has_key('rx_location'):
1174 1174 rx_location = kwargs['rx_location']
1175 1175 else: rx_location = [(0,1),(1,1),(1,0)]
1176 1176 if kwargs.has_key('azimuth'):
1177 1177 azimuth = kwargs['azimuth']
1178 1178 else: azimuth = 51.06
1179 1179 if kwargs.has_key('dfactor'):
1180 1180 dfactor = kwargs['dfactor']
1181 1181 if kwargs.has_key('mode'):
1182 1182 mode = kwargs['mode']
1183 1183 if kwargs.has_key('theta_x'):
1184 1184 theta_x = kwargs['theta_x']
1185 1185 if kwargs.has_key('theta_y'):
1186 1186 theta_y = kwargs['theta_y']
1187 1187 else: mode = 'SA'
1188 1188
1189 1189 #Borrar luego esto
1190 1190 if dataOut.groupList is None:
1191 1191 dataOut.groupList = [(0,1),(0,2),(1,2)]
1192 1192 groupList = dataOut.groupList
1193 1193 C = 3e8
1194 1194 freq = 50e6
1195 1195 lamb = C/freq
1196 1196 k = 2*numpy.pi/lamb
1197 1197
1198 1198 timeList = dataOut.abscissaList
1199 1199 heightList = dataOut.heightList
1200 1200
1201 1201 if self.__isConfig == False:
1202 1202 dataOut.outputInterval = nMins*60
1203 1203 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1204 1204 #Get Initial LTC time
1205 1205 initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
1206 1206 minuteAux = initime.minute
1207 1207 minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)
1208 1208 self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1209 1209
1210 1210 self.__isConfig = True
1211 1211
1212 1212 if self.__buffer is None:
1213 1213 self.__buffer = dataOut.data_param
1214 1214 self.__firstdata = copy.copy(dataOut)
1215 1215
1216 1216 else:
1217 1217 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1218 1218
1219 1219 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1220 1220
1221 1221 if self.__dataReady:
1222 1222 dataOut.utctimeInit = self.__initime
1223 1223 self.__initime += dataOut.outputInterval #to erase time offset
1224 1224
1225 1225 metArray = self.__buffer
1226 1226 if mode == 'SA':
1227 1227 dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)
1228 1228 elif mode == 'DBS':
1229 1229 dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList, azimuth=azimuth, theta_x=theta_x, theta_y=theta_y)
1230 1230 dataOut.data_output = dataOut.data_output.T
1231 1231 dataOut.flagNoData = False
1232 1232 self.__buffer = None
1233 1233
1234 1234 return
1235 1235
1236 1236 class EWDriftsEstimation(Operation):
1237 1237
1238 1238
1239 1239 def __correctValues(self, heiRang, phi, velRadial, SNR):
1240 1240 listPhi = phi.tolist()
1241 1241 maxid = listPhi.index(max(listPhi))
1242 1242 minid = listPhi.index(min(listPhi))
1243 1243
1244 1244 rango = range(len(phi))
1245 1245 # rango = numpy.delete(rango,maxid)
1246 1246
1247 1247 heiRang1 = heiRang*math.cos(phi[maxid])
1248 1248 heiRangAux = heiRang*math.cos(phi[minid])
1249 1249 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1250 1250 heiRang1 = numpy.delete(heiRang1,indOut)
1251 1251
1252 1252 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1253 1253 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1254 1254
1255 1255 for i in rango:
1256 1256 x = heiRang*math.cos(phi[i])
1257 1257 y1 = velRadial[i,:]
1258 1258 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1259 1259
1260 1260 x1 = heiRang1
1261 1261 y11 = f1(x1)
1262 1262
1263 1263 y2 = SNR[i,:]
1264 1264 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1265 1265 y21 = f2(x1)
1266 1266
1267 1267 velRadial1[i,:] = y11
1268 1268 SNR1[i,:] = y21
1269 1269
1270 1270 return heiRang1, velRadial1, SNR1
1271 1271
1272 1272 def run(self, dataOut, zenith, zenithCorrection):
1273 1273 heiRang = dataOut.heightList
1274 1274 velRadial = dataOut.data_param[:,3,:]
1275 1275 SNR = dataOut.data_SNR
1276 1276
1277 1277 zenith = numpy.array(zenith)
1278 1278 zenith -= zenithCorrection
1279 1279 zenith *= numpy.pi/180
1280 1280
1281 1281 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
1282 1282
1283 1283 alp = zenith[0]
1284 1284 bet = zenith[1]
1285 1285
1286 1286 w_w = velRadial1[0,:]
1287 1287 w_e = velRadial1[1,:]
1288 1288
1289 1289 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
1290 1290 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
1291 1291
1292 1292 winds = numpy.vstack((u,w))
1293 1293
1294 1294 dataOut.heightList = heiRang1
1295 1295 dataOut.data_output = winds
1296 1296 dataOut.data_SNR = SNR1
1297 1297
1298 1298 dataOut.utctimeInit = dataOut.utctime
1299 1299 dataOut.outputInterval = dataOut.timeInterval
1300 1300 return
1301 1301
1302 1302 #--------------- Non Specular Meteor ----------------
1303 1303
1304 1304 class NonSpecularMeteorDetection(Operation):
1305 1305
1306 1306 def run(self, dataOut, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False):
1307 1307 data_acf = dataOut.data_pre[0]
1308 1308 data_ccf = dataOut.data_pre[1]
1309 1309 pairsList = dataOut.groupList[1]
1310 1310
1311 1311 lamb = dataOut.C/dataOut.frequency
1312 1312 tSamp = dataOut.ippSeconds*dataOut.nCohInt
1313 1313 paramInterval = dataOut.paramInterval
1314 1314
1315 1315 nChannels = data_acf.shape[0]
1316 1316 nLags = data_acf.shape[1]
1317 1317 nProfiles = data_acf.shape[2]
1318 1318 nHeights = dataOut.nHeights
1319 1319 nCohInt = dataOut.nCohInt
1320 1320 sec = numpy.round(nProfiles/dataOut.paramInterval)
1321 1321 heightList = dataOut.heightList
1322 1322 ippSeconds = dataOut.ippSeconds*dataOut.nCohInt*dataOut.nAvg
1323 1323 utctime = dataOut.utctime
1324 1324
1325 1325 dataOut.abscissaList = numpy.arange(0,paramInterval+ippSeconds,ippSeconds)
1326 1326
1327 1327 #------------------------ SNR --------------------------------------
1328 1328 power = data_acf[:,0,:,:].real
1329 1329 noise = numpy.zeros(nChannels)
1330 1330 SNR = numpy.zeros(power.shape)
1331 1331 for i in range(nChannels):
1332 1332 noise[i] = hildebrand_sekhon(power[i,:], nCohInt)
1333 1333 SNR[i] = (power[i]-noise[i])/noise[i]
1334 1334 SNRm = numpy.nanmean(SNR, axis = 0)
1335 1335 SNRdB = 10*numpy.log10(SNR)
1336 1336
1337 1337 if mode == 'SA':
1338 1338 dataOut.groupList = dataOut.groupList[1]
1339 1339 nPairs = data_ccf.shape[0]
1340 1340 #---------------------- Coherence and Phase --------------------------
1341 1341 phase = numpy.zeros(data_ccf[:,0,:,:].shape)
1342 1342 # phase1 = numpy.copy(phase)
1343 1343 coh1 = numpy.zeros(data_ccf[:,0,:,:].shape)
1344 1344
1345 1345 for p in range(nPairs):
1346 1346 ch0 = pairsList[p][0]
1347 1347 ch1 = pairsList[p][1]
1348 1348 ccf = data_ccf[p,0,:,:]/numpy.sqrt(data_acf[ch0,0,:,:]*data_acf[ch1,0,:,:])
1349 1349 phase[p,:,:] = ndimage.median_filter(numpy.angle(ccf), size = (5,1)) #median filter
1350 1350 # phase1[p,:,:] = numpy.angle(ccf) #median filter
1351 1351 coh1[p,:,:] = ndimage.median_filter(numpy.abs(ccf), 5) #median filter
1352 1352 # coh1[p,:,:] = numpy.abs(ccf) #median filter
1353 1353 coh = numpy.nanmax(coh1, axis = 0)
1354 1354 # struc = numpy.ones((5,1))
1355 1355 # coh = ndimage.morphology.grey_dilation(coh, size=(10,1))
1356 1356 #---------------------- Radial Velocity ----------------------------
1357 1357 phaseAux = numpy.mean(numpy.angle(data_acf[:,1,:,:]), axis = 0)
1358 1358 velRad = phaseAux*lamb/(4*numpy.pi*tSamp)
1359 1359
1360 1360 if allData:
1361 1361 boolMetFin = ~numpy.isnan(SNRm)
1362 1362 # coh[:-1,:] = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
1363 1363 else:
1364 1364 #------------------------ Meteor mask ---------------------------------
1365 1365 # #SNR mask
1366 1366 # boolMet = (SNRdB>SNRthresh)#|(~numpy.isnan(SNRdB))
1367 1367 #
1368 1368 # #Erase small objects
1369 1369 # boolMet1 = self.__erase_small(boolMet, 2*sec, 5)
1370 1370 #
1371 1371 # auxEEJ = numpy.sum(boolMet1,axis=0)
1372 1372 # indOver = auxEEJ>nProfiles*0.8 #Use this later
1373 1373 # indEEJ = numpy.where(indOver)[0]
1374 1374 # indNEEJ = numpy.where(~indOver)[0]
1375 1375 #
1376 1376 # boolMetFin = boolMet1
1377 1377 #
1378 1378 # if indEEJ.size > 0:
1379 1379 # boolMet1[:,indEEJ] = False #Erase heights with EEJ
1380 1380 #
1381 1381 # boolMet2 = coh > cohThresh
1382 1382 # boolMet2 = self.__erase_small(boolMet2, 2*sec,5)
1383 1383 #
1384 1384 # #Final Meteor mask
1385 1385 # boolMetFin = boolMet1|boolMet2
1386 1386
1387 1387 #Coherence mask
1388 1388 boolMet1 = coh > 0.75
1389 1389 struc = numpy.ones((30,1))
1390 1390 boolMet1 = ndimage.morphology.binary_dilation(boolMet1, structure=struc)
1391 1391
1392 1392 #Derivative mask
1393 1393 derPhase = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
1394 1394 boolMet2 = derPhase < 0.2
1395 1395 # boolMet2 = ndimage.morphology.binary_opening(boolMet2)
1396 1396 # boolMet2 = ndimage.morphology.binary_closing(boolMet2, structure = numpy.ones((10,1)))
1397 1397 boolMet2 = ndimage.median_filter(boolMet2,size=5)
1398 1398 boolMet2 = numpy.vstack((boolMet2,numpy.full((1,nHeights), True, dtype=bool)))
1399 1399 # #Final mask
1400 1400 # boolMetFin = boolMet2
1401 1401 boolMetFin = boolMet1&boolMet2
1402 1402 # boolMetFin = ndimage.morphology.binary_dilation(boolMetFin)
1403 1403 #Creating data_param
1404 1404 coordMet = numpy.where(boolMetFin)
1405 1405
1406 1406 tmet = coordMet[0]
1407 1407 hmet = coordMet[1]
1408 1408
1409 1409 data_param = numpy.zeros((tmet.size, 6 + nPairs))
1410 1410 data_param[:,0] = utctime
1411 1411 data_param[:,1] = tmet
1412 1412 data_param[:,2] = hmet
1413 1413 data_param[:,3] = SNRm[tmet,hmet]
1414 1414 data_param[:,4] = velRad[tmet,hmet]
1415 1415 data_param[:,5] = coh[tmet,hmet]
1416 1416 data_param[:,6:] = phase[:,tmet,hmet].T
1417 1417
1418 1418 elif mode == 'DBS':
1419 1419 dataOut.groupList = numpy.arange(nChannels)
1420 1420
1421 1421 #Radial Velocities
1422 1422 phase = numpy.angle(data_acf[:,1,:,:])
1423 1423 # phase = ndimage.median_filter(numpy.angle(data_acf[:,1,:,:]), size = (1,5,1))
1424 1424 velRad = phase*lamb/(4*numpy.pi*tSamp)
1425 1425
1426 1426 #Spectral width
1427 1427 # acf1 = ndimage.median_filter(numpy.abs(data_acf[:,1,:,:]), size = (1,5,1))
1428 1428 # acf2 = ndimage.median_filter(numpy.abs(data_acf[:,2,:,:]), size = (1,5,1))
1429 1429 acf1 = data_acf[:,1,:,:]
1430 1430 acf2 = data_acf[:,2,:,:]
1431 1431
1432 1432 spcWidth = (lamb/(2*numpy.sqrt(6)*numpy.pi*tSamp))*numpy.sqrt(numpy.log(acf1/acf2))
1433 1433 # velRad = ndimage.median_filter(velRad, size = (1,5,1))
1434 1434 if allData:
1435 1435 boolMetFin = ~numpy.isnan(SNRdB)
1436 1436 else:
1437 1437 #SNR
1438 1438 boolMet1 = (SNRdB>SNRthresh) #SNR mask
1439 1439 boolMet1 = ndimage.median_filter(boolMet1, size=(1,5,5))
1440 1440
1441 1441 #Radial velocity
1442 1442 boolMet2 = numpy.abs(velRad) < 20
1443 1443 boolMet2 = ndimage.median_filter(boolMet2, (1,5,5))
1444 1444
1445 1445 #Spectral Width
1446 1446 boolMet3 = spcWidth < 30
1447 1447 boolMet3 = ndimage.median_filter(boolMet3, (1,5,5))
1448 1448 # boolMetFin = self.__erase_small(boolMet1, 10,5)
1449 1449 boolMetFin = boolMet1&boolMet2&boolMet3
1450 1450
1451 1451 #Creating data_param
1452 1452 coordMet = numpy.where(boolMetFin)
1453 1453
1454 1454 cmet = coordMet[0]
1455 1455 tmet = coordMet[1]
1456 1456 hmet = coordMet[2]
1457 1457
1458 1458 data_param = numpy.zeros((tmet.size, 7))
1459 1459 data_param[:,0] = utctime
1460 1460 data_param[:,1] = cmet
1461 1461 data_param[:,2] = tmet
1462 1462 data_param[:,3] = hmet
1463 1463 data_param[:,4] = SNR[cmet,tmet,hmet].T
1464 1464 data_param[:,5] = velRad[cmet,tmet,hmet].T
1465 1465 data_param[:,6] = spcWidth[cmet,tmet,hmet].T
1466 1466
1467 1467 # self.dataOut.data_param = data_int
1468 1468 if len(data_param) == 0:
1469 1469 dataOut.flagNoData = True
1470 1470 else:
1471 1471 dataOut.data_param = data_param
1472 1472
1473 1473 def __erase_small(self, binArray, threshX, threshY):
1474 1474 labarray, numfeat = ndimage.measurements.label(binArray)
1475 1475 binArray1 = numpy.copy(binArray)
1476 1476
1477 1477 for i in range(1,numfeat + 1):
1478 1478 auxBin = (labarray==i)
1479 1479 auxSize = auxBin.sum()
1480 1480
1481 1481 x,y = numpy.where(auxBin)
1482 1482 widthX = x.max() - x.min()
1483 1483 widthY = y.max() - y.min()
1484 1484
1485 1485 #width X: 3 seg -> 12.5*3
1486 1486 #width Y:
1487 1487
1488 1488 if (auxSize < 50) or (widthX < threshX) or (widthY < threshY):
1489 1489 binArray1[auxBin] = False
1490 1490
1491 1491 return binArray1
1492 1492
1493 1493 #--------------- Specular Meteor ----------------
1494 1494
1495 1495 class SMDetection(Operation):
1496 1496 '''
1497 1497 Function DetectMeteors()
1498 1498 Project developed with paper:
1499 1499 HOLDSWORTH ET AL. 2004
1500 1500
1501 1501 Input:
1502 1502 self.dataOut.data_pre
1503 1503
1504 1504 centerReceiverIndex: From the channels, which is the center receiver
1505 1505
1506 1506 hei_ref: Height reference for the Beacon signal extraction
1507 1507 tauindex:
1508 1508 predefinedPhaseShifts: Predefined phase offset for the voltge signals
1509 1509
1510 1510 cohDetection: Whether to user Coherent detection or not
1511 1511 cohDet_timeStep: Coherent Detection calculation time step
1512 1512 cohDet_thresh: Coherent Detection phase threshold to correct phases
1513 1513
1514 1514 noise_timeStep: Noise calculation time step
1515 1515 noise_multiple: Noise multiple to define signal threshold
1516 1516
1517 1517 multDet_timeLimit: Multiple Detection Removal time limit in seconds
1518 1518 multDet_rangeLimit: Multiple Detection Removal range limit in km
1519 1519
1520 1520 phaseThresh: Maximum phase difference between receiver to be consider a meteor
1521 1521 SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
1522 1522
1523 1523 hmin: Minimum Height of the meteor to use it in the further wind estimations
1524 1524 hmax: Maximum Height of the meteor to use it in the further wind estimations
1525 1525 azimuth: Azimuth angle correction
1526 1526
1527 1527 Affected:
1528 1528 self.dataOut.data_param
1529 1529
1530 1530 Rejection Criteria (Errors):
1531 1531 0: No error; analysis OK
1532 1532 1: SNR < SNR threshold
1533 1533 2: angle of arrival (AOA) ambiguously determined
1534 1534 3: AOA estimate not feasible
1535 1535 4: Large difference in AOAs obtained from different antenna baselines
1536 1536 5: echo at start or end of time series
1537 1537 6: echo less than 5 examples long; too short for analysis
1538 1538 7: echo rise exceeds 0.3s
1539 1539 8: echo decay time less than twice rise time
1540 1540 9: large power level before echo
1541 1541 10: large power level after echo
1542 1542 11: poor fit to amplitude for estimation of decay time
1543 1543 12: poor fit to CCF phase variation for estimation of radial drift velocity
1544 1544 13: height unresolvable echo: not valid height within 70 to 110 km
1545 1545 14: height ambiguous echo: more then one possible height within 70 to 110 km
1546 1546 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
1547 1547 16: oscilatory echo, indicating event most likely not an underdense echo
1548 1548
1549 1549 17: phase difference in meteor Reestimation
1550 1550
1551 1551 Data Storage:
1552 1552 Meteors for Wind Estimation (8):
1553 1553 Utc Time | Range Height
1554 1554 Azimuth Zenith errorCosDir
1555 1555 VelRad errorVelRad
1556 1556 Phase0 Phase1 Phase2 Phase3
1557 1557 TypeError
1558 1558
1559 1559 '''
1560 1560
1561 1561 def run(self, dataOut, hei_ref = None, tauindex = 0,
1562 1562 phaseOffsets = None,
1563 1563 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
1564 1564 noise_timeStep = 4, noise_multiple = 4,
1565 1565 multDet_timeLimit = 1, multDet_rangeLimit = 3,
1566 1566 phaseThresh = 20, SNRThresh = 5,
1567 1567 hmin = 50, hmax=150, azimuth = 0,
1568 1568 channelPositions = None) :
1569 1569
1570 1570
1571 1571 #Getting Pairslist
1572 1572 if channelPositions is None:
1573 1573 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
1574 1574 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
1575 1575 meteorOps = SMOperations()
1576 1576 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
1577 1577 heiRang = dataOut.getHeiRange()
1578 1578 #Get Beacon signal - No Beacon signal anymore
1579 1579 # newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
1580 1580 #
1581 1581 # if hei_ref != None:
1582 1582 # newheis = numpy.where(self.dataOut.heightList>hei_ref)
1583 1583 #
1584 1584
1585 1585
1586 1586 #****************REMOVING HARDWARE PHASE DIFFERENCES***************
1587 1587 # see if the user put in pre defined phase shifts
1588 1588 voltsPShift = dataOut.data_pre.copy()
1589 1589
1590 1590 # if predefinedPhaseShifts != None:
1591 1591 # hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
1592 1592 #
1593 1593 # # elif beaconPhaseShifts:
1594 1594 # # #get hardware phase shifts using beacon signal
1595 1595 # # hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
1596 1596 # # hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
1597 1597 #
1598 1598 # else:
1599 1599 # hardwarePhaseShifts = numpy.zeros(5)
1600 1600 #
1601 1601 # voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
1602 1602 # for i in range(self.dataOut.data_pre.shape[0]):
1603 1603 # voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
1604 1604
1605 1605 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
1606 1606
1607 1607 #Remove DC
1608 1608 voltsDC = numpy.mean(voltsPShift,1)
1609 1609 voltsDC = numpy.mean(voltsDC,1)
1610 1610 for i in range(voltsDC.shape[0]):
1611 1611 voltsPShift[i] = voltsPShift[i] - voltsDC[i]
1612 1612
1613 1613 #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
1614 1614 # voltsPShift = voltsPShift[:,:,:newheis[0][0]]
1615 1615
1616 1616 #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
1617 1617 #Coherent Detection
1618 1618 if cohDetection:
1619 1619 #use coherent detection to get the net power
1620 1620 cohDet_thresh = cohDet_thresh*numpy.pi/180
1621 1621 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, dataOut.timeInterval, pairslist0, cohDet_thresh)
1622 1622
1623 1623 #Non-coherent detection!
1624 1624 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
1625 1625 #********** END OF COH/NON-COH POWER CALCULATION**********************
1626 1626
1627 1627 #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
1628 1628 #Get noise
1629 1629 noise, noise1 = self.__getNoise(powerNet, noise_timeStep, dataOut.timeInterval)
1630 1630 # noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
1631 1631 #Get signal threshold
1632 1632 signalThresh = noise_multiple*noise
1633 1633 #Meteor echoes detection
1634 1634 listMeteors = self.__findMeteors(powerNet, signalThresh)
1635 1635 #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
1636 1636
1637 1637 #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
1638 1638 #Parameters
1639 1639 heiRange = dataOut.getHeiRange()
1640 1640 rangeInterval = heiRange[1] - heiRange[0]
1641 1641 rangeLimit = multDet_rangeLimit/rangeInterval
1642 1642 timeLimit = multDet_timeLimit/dataOut.timeInterval
1643 1643 #Multiple detection removals
1644 1644 listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
1645 1645 #************ END OF REMOVE MULTIPLE DETECTIONS **********************
1646 1646
1647 1647 #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************
1648 1648 #Parameters
1649 1649 phaseThresh = phaseThresh*numpy.pi/180
1650 1650 thresh = [phaseThresh, noise_multiple, SNRThresh]
1651 1651 #Meteor reestimation (Errors N 1, 6, 12, 17)
1652 1652 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist0, thresh, noise, dataOut.timeInterval, dataOut.frequency)
1653 1653 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
1654 1654 #Estimation of decay times (Errors N 7, 8, 11)
1655 1655 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, dataOut.timeInterval, dataOut.frequency)
1656 1656 #******************* END OF METEOR REESTIMATION *******************
1657 1657
1658 1658 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
1659 1659 #Calculating Radial Velocity (Error N 15)
1660 1660 radialStdThresh = 10
1661 1661 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, dataOut.timeInterval)
1662 1662
1663 1663 if len(listMeteors4) > 0:
1664 1664 #Setting New Array
1665 1665 date = dataOut.utctime
1666 1666 arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
1667 1667
1668 1668 #Correcting phase offset
1669 1669 if phaseOffsets != None:
1670 1670 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
1671 1671 arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
1672 1672
1673 1673 #Second Pairslist
1674 1674 pairsList = []
1675 1675 pairx = (0,1)
1676 1676 pairy = (2,3)
1677 1677 pairsList.append(pairx)
1678 1678 pairsList.append(pairy)
1679 1679
1680 1680 jph = numpy.array([0,0,0,0])
1681 1681 h = (hmin,hmax)
1682 1682 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
1683 1683
1684 1684 # #Calculate AOA (Error N 3, 4)
1685 1685 # #JONES ET AL. 1998
1686 1686 # error = arrayParameters[:,-1]
1687 1687 # AOAthresh = numpy.pi/8
1688 1688 # phases = -arrayParameters[:,9:13]
1689 1689 # arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)
1690 1690 #
1691 1691 # #Calculate Heights (Error N 13 and 14)
1692 1692 # error = arrayParameters[:,-1]
1693 1693 # Ranges = arrayParameters[:,2]
1694 1694 # zenith = arrayParameters[:,5]
1695 1695 # arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)
1696 1696 # error = arrayParameters[:,-1]
1697 1697 #********************* END OF PARAMETERS CALCULATION **************************
1698 1698
1699 1699 #***************************+ PASS DATA TO NEXT STEP **********************
1700 1700 # arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
1701 1701 dataOut.data_param = arrayParameters
1702 1702
1703 1703 if arrayParameters is None:
1704 1704 dataOut.flagNoData = True
1705 1705 else:
1706 1706 dataOut.flagNoData = True
1707 1707
1708 1708 return
1709 1709
1710 1710 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
1711 1711
1712 1712 minIndex = min(newheis[0])
1713 1713 maxIndex = max(newheis[0])
1714 1714
1715 1715 voltage = voltage0[:,:,minIndex:maxIndex+1]
1716 1716 nLength = voltage.shape[1]/n
1717 1717 nMin = 0
1718 1718 nMax = 0
1719 1719 phaseOffset = numpy.zeros((len(pairslist),n))
1720 1720
1721 1721 for i in range(n):
1722 1722 nMax += nLength
1723 1723 phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
1724 1724 phaseCCF = numpy.mean(phaseCCF, axis = 2)
1725 1725 phaseOffset[:,i] = phaseCCF.transpose()
1726 1726 nMin = nMax
1727 1727 # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
1728 1728
1729 1729 #Remove Outliers
1730 1730 factor = 2
1731 1731 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
1732 1732 dw = numpy.std(wt,axis = 1)
1733 1733 dw = dw.reshape((dw.size,1))
1734 1734 ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
1735 1735 phaseOffset[ind] = numpy.nan
1736 1736 phaseOffset = stats.nanmean(phaseOffset, axis=1)
1737 1737
1738 1738 return phaseOffset
1739 1739
1740 1740 def __shiftPhase(self, data, phaseShift):
1741 1741 #this will shift the phase of a complex number
1742 1742 dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
1743 1743 return dataShifted
1744 1744
1745 1745 def __estimatePhaseDifference(self, array, pairslist):
1746 1746 nChannel = array.shape[0]
1747 1747 nHeights = array.shape[2]
1748 1748 numPairs = len(pairslist)
1749 1749 # phaseCCF = numpy.zeros((nChannel, 5, nHeights))
1750 1750 phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
1751 1751
1752 1752 #Correct phases
1753 1753 derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
1754 1754 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
1755 1755
1756 1756 if indDer[0].shape[0] > 0:
1757 1757 for i in range(indDer[0].shape[0]):
1758 1758 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
1759 1759 phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
1760 1760
1761 1761 # for j in range(numSides):
1762 1762 # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
1763 1763 # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
1764 1764 #
1765 1765 #Linear
1766 1766 phaseInt = numpy.zeros((numPairs,1))
1767 1767 angAllCCF = phaseCCF[:,[0,1,3,4],0]
1768 1768 for j in range(numPairs):
1769 1769 fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
1770 1770 phaseInt[j] = fit[1]
1771 1771 #Phase Differences
1772 1772 phaseDiff = phaseInt - phaseCCF[:,2,:]
1773 1773 phaseArrival = phaseInt.reshape(phaseInt.size)
1774 1774
1775 1775 #Dealias
1776 1776 phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival))
1777 1777 # indAlias = numpy.where(phaseArrival > numpy.pi)
1778 1778 # phaseArrival[indAlias] -= 2*numpy.pi
1779 1779 # indAlias = numpy.where(phaseArrival < -numpy.pi)
1780 1780 # phaseArrival[indAlias] += 2*numpy.pi
1781 1781
1782 1782 return phaseDiff, phaseArrival
1783 1783
1784 1784 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
1785 1785 #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
1786 1786 #find the phase shifts of each channel over 1 second intervals
1787 1787 #only look at ranges below the beacon signal
1788 1788 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
1789 1789 numBlocks = int(volts.shape[1]/numProfPerBlock)
1790 1790 numHeights = volts.shape[2]
1791 1791 nChannel = volts.shape[0]
1792 1792 voltsCohDet = volts.copy()
1793 1793
1794 1794 pairsarray = numpy.array(pairslist)
1795 1795 indSides = pairsarray[:,1]
1796 1796 # indSides = numpy.array(range(nChannel))
1797 1797 # indSides = numpy.delete(indSides, indCenter)
1798 1798 #
1799 1799 # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
1800 1800 listBlocks = numpy.array_split(volts, numBlocks, 1)
1801 1801
1802 1802 startInd = 0
1803 1803 endInd = 0
1804 1804
1805 1805 for i in range(numBlocks):
1806 1806 startInd = endInd
1807 1807 endInd = endInd + listBlocks[i].shape[1]
1808 1808
1809 1809 arrayBlock = listBlocks[i]
1810 1810 # arrayBlockCenter = listCenter[i]
1811 1811
1812 1812 #Estimate the Phase Difference
1813 1813 phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
1814 1814 #Phase Difference RMS
1815 1815 arrayPhaseRMS = numpy.abs(phaseDiff)
1816 1816 phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
1817 1817 indPhase = numpy.where(phaseRMSaux==4)
1818 1818 #Shifting
1819 1819 if indPhase[0].shape[0] > 0:
1820 1820 for j in range(indSides.size):
1821 1821 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
1822 1822 voltsCohDet[:,startInd:endInd,:] = arrayBlock
1823 1823
1824 1824 return voltsCohDet
1825 1825
1826 1826 def __calculateCCF(self, volts, pairslist ,laglist):
1827 1827
1828 1828 nHeights = volts.shape[2]
1829 1829 nPoints = volts.shape[1]
1830 1830 voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
1831 1831
1832 1832 for i in range(len(pairslist)):
1833 1833 volts1 = volts[pairslist[i][0]]
1834 1834 volts2 = volts[pairslist[i][1]]
1835 1835
1836 1836 for t in range(len(laglist)):
1837 1837 idxT = laglist[t]
1838 1838 if idxT >= 0:
1839 1839 vStacked = numpy.vstack((volts2[idxT:,:],
1840 1840 numpy.zeros((idxT, nHeights),dtype='complex')))
1841 1841 else:
1842 1842 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
1843 1843 volts2[:(nPoints + idxT),:]))
1844 1844 voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
1845 1845
1846 1846 vStacked = None
1847 1847 return voltsCCF
1848 1848
1849 1849 def __getNoise(self, power, timeSegment, timeInterval):
1850 1850 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
1851 1851 numBlocks = int(power.shape[0]/numProfPerBlock)
1852 1852 numHeights = power.shape[1]
1853 1853
1854 1854 listPower = numpy.array_split(power, numBlocks, 0)
1855 1855 noise = numpy.zeros((power.shape[0], power.shape[1]))
1856 1856 noise1 = numpy.zeros((power.shape[0], power.shape[1]))
1857 1857
1858 1858 startInd = 0
1859 1859 endInd = 0
1860 1860
1861 1861 for i in range(numBlocks): #split por canal
1862 1862 startInd = endInd
1863 1863 endInd = endInd + listPower[i].shape[0]
1864 1864
1865 1865 arrayBlock = listPower[i]
1866 1866 noiseAux = numpy.mean(arrayBlock, 0)
1867 1867 # noiseAux = numpy.median(noiseAux)
1868 1868 # noiseAux = numpy.mean(arrayBlock)
1869 1869 noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
1870 1870
1871 1871 noiseAux1 = numpy.mean(arrayBlock)
1872 1872 noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
1873 1873
1874 1874 return noise, noise1
1875 1875
1876 1876 def __findMeteors(self, power, thresh):
1877 1877 nProf = power.shape[0]
1878 1878 nHeights = power.shape[1]
1879 1879 listMeteors = []
1880 1880
1881 1881 for i in range(nHeights):
1882 1882 powerAux = power[:,i]
1883 1883 threshAux = thresh[:,i]
1884 1884
1885 1885 indUPthresh = numpy.where(powerAux > threshAux)[0]
1886 1886 indDNthresh = numpy.where(powerAux <= threshAux)[0]
1887 1887
1888 1888 j = 0
1889 1889
1890 1890 while (j < indUPthresh.size - 2):
1891 1891 if (indUPthresh[j + 2] == indUPthresh[j] + 2):
1892 1892 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
1893 1893 indDNthresh = indDNthresh[indDNAux]
1894 1894
1895 1895 if (indDNthresh.size > 0):
1896 1896 indEnd = indDNthresh[0] - 1
1897 1897 indInit = indUPthresh[j] if isinstance(indUPthresh[j], (int, float)) else indUPthresh[j][0] ##CHECK!!!!
1898 1898
1899 1899 meteor = powerAux[indInit:indEnd + 1]
1900 1900 indPeak = meteor.argmax() + indInit
1901 1901 FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
1902 1902
1903 1903 listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
1904 1904 j = numpy.where(indUPthresh == indEnd)[0] + 1
1905 1905 else: j+=1
1906 1906 else: j+=1
1907 1907
1908 1908 return listMeteors
1909 1909
1910 1910 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
1911 1911
1912 1912 arrayMeteors = numpy.asarray(listMeteors)
1913 1913 listMeteors1 = []
1914 1914
1915 1915 while arrayMeteors.shape[0] > 0:
1916 1916 FLAs = arrayMeteors[:,4]
1917 1917 maxFLA = FLAs.argmax()
1918 1918 listMeteors1.append(arrayMeteors[maxFLA,:])
1919 1919
1920 1920 MeteorInitTime = arrayMeteors[maxFLA,1]
1921 1921 MeteorEndTime = arrayMeteors[maxFLA,3]
1922 1922 MeteorHeight = arrayMeteors[maxFLA,0]
1923 1923
1924 1924 #Check neighborhood
1925 1925 maxHeightIndex = MeteorHeight + rangeLimit
1926 1926 minHeightIndex = MeteorHeight - rangeLimit
1927 1927 minTimeIndex = MeteorInitTime - timeLimit
1928 1928 maxTimeIndex = MeteorEndTime + timeLimit
1929 1929
1930 1930 #Check Heights
1931 1931 indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
1932 1932 indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
1933 1933 indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
1934 1934
1935 1935 arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
1936 1936
1937 1937 return listMeteors1
1938 1938
1939 1939 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
1940 1940 numHeights = volts.shape[2]
1941 1941 nChannel = volts.shape[0]
1942 1942
1943 1943 thresholdPhase = thresh[0]
1944 1944 thresholdNoise = thresh[1]
1945 1945 thresholdDB = float(thresh[2])
1946 1946
1947 1947 thresholdDB1 = 10**(thresholdDB/10)
1948 1948 pairsarray = numpy.array(pairslist)
1949 1949 indSides = pairsarray[:,1]
1950 1950
1951 1951 pairslist1 = list(pairslist)
1952 1952 pairslist1.append((0,4))
1953 1953 pairslist1.append((1,3))
1954 1954
1955 1955 listMeteors1 = []
1956 1956 listPowerSeries = []
1957 1957 listVoltageSeries = []
1958 1958 #volts has the war data
1959 1959
1960 1960 if frequency == 30.175e6:
1961 1961 timeLag = 45*10**-3
1962 1962 else:
1963 1963 timeLag = 15*10**-3
1964 1964 lag = int(numpy.ceil(timeLag/timeInterval))
1965 1965
1966 1966 for i in range(len(listMeteors)):
1967 1967
1968 1968 ###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
1969 1969 meteorAux = numpy.zeros(16)
1970 1970
1971 1971 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
1972 1972 mHeight = int(listMeteors[i][0])
1973 1973 mStart = int(listMeteors[i][1])
1974 1974 mPeak = int(listMeteors[i][2])
1975 1975 mEnd = int(listMeteors[i][3])
1976 1976
1977 1977 #get the volt data between the start and end times of the meteor
1978 1978 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
1979 1979 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
1980 1980
1981 1981 #3.6. Phase Difference estimation
1982 1982 phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
1983 1983
1984 1984 #3.7. Phase difference removal & meteor start, peak and end times reestimated
1985 1985 #meteorVolts0.- all Channels, all Profiles
1986 1986 meteorVolts0 = volts[:,:,mHeight]
1987 1987 meteorThresh = noise[:,mHeight]*thresholdNoise
1988 1988 meteorNoise = noise[:,mHeight]
1989 1989 meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
1990 1990 powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
1991 1991
1992 1992 #Times reestimation
1993 1993 mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
1994 1994 if mStart1.size > 0:
1995 1995 mStart1 = mStart1[-1] + 1
1996 1996
1997 1997 else:
1998 1998 mStart1 = mPeak
1999 1999
2000 2000 mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
2001 2001 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
2002 2002 if mEndDecayTime1.size == 0:
2003 2003 mEndDecayTime1 = powerNet0.size
2004 2004 else:
2005 2005 mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
2006 2006 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
2007 2007
2008 2008 #meteorVolts1.- all Channels, from start to end
2009 2009 meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
2010 2010 meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
2011 2011 if meteorVolts2.shape[1] == 0:
2012 2012 meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
2013 2013 meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
2014 2014 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
2015 2015 ##################### END PARAMETERS REESTIMATION #########################
2016 2016
2017 2017 ##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
2018 2018 # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
2019 2019 if meteorVolts2.shape[1] > 0:
2020 2020 #Phase Difference re-estimation
2021 2021 phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
2022 2022 # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
2023 2023 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
2024 2024 phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
2025 2025 meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
2026 2026
2027 2027 #Phase Difference RMS
2028 2028 phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
2029 2029 powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
2030 2030 #Data from Meteor
2031 2031 mPeak1 = powerNet1.argmax() + mStart1
2032 2032 mPeakPower1 = powerNet1.max()
2033 2033 noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
2034 2034 mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
2035 2035 Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
2036 2036 Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
2037 2037 PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
2038 2038 #Vectorize
2039 2039 meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
2040 2040 meteorAux[7:11] = phaseDiffint[0:4]
2041 2041
2042 2042 #Rejection Criterions
2043 2043 if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
2044 2044 meteorAux[-1] = 17
2045 2045 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
2046 2046 meteorAux[-1] = 1
2047 2047
2048 2048
2049 2049 else:
2050 2050 meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
2051 2051 meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
2052 2052 PowerSeries = 0
2053 2053
2054 2054 listMeteors1.append(meteorAux)
2055 2055 listPowerSeries.append(PowerSeries)
2056 2056 listVoltageSeries.append(meteorVolts1)
2057 2057
2058 2058 return listMeteors1, listPowerSeries, listVoltageSeries
2059 2059
2060 2060 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
2061 2061
2062 2062 threshError = 10
2063 2063 #Depending if it is 30 or 50 MHz
2064 2064 if frequency == 30.175e6:
2065 2065 timeLag = 45*10**-3
2066 2066 else:
2067 2067 timeLag = 15*10**-3
2068 2068 lag = int(numpy.ceil(timeLag/timeInterval))
2069 2069
2070 2070 listMeteors1 = []
2071 2071
2072 2072 for i in range(len(listMeteors)):
2073 2073 meteorPower = listPower[i]
2074 2074 meteorAux = listMeteors[i]
2075 2075
2076 2076 if meteorAux[-1] == 0:
2077 2077
2078 2078 try:
2079 2079 indmax = meteorPower.argmax()
2080 2080 indlag = indmax + lag
2081 2081
2082 2082 y = meteorPower[indlag:]
2083 2083 x = numpy.arange(0, y.size)*timeLag
2084 2084
2085 2085 #first guess
2086 2086 a = y[0]
2087 2087 tau = timeLag
2088 2088 #exponential fit
2089 2089 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
2090 2090 y1 = self.__exponential_function(x, *popt)
2091 2091 #error estimation
2092 2092 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
2093 2093
2094 2094 decayTime = popt[1]
2095 2095 riseTime = indmax*timeInterval
2096 2096 meteorAux[11:13] = [decayTime, error]
2097 2097
2098 2098 #Table items 7, 8 and 11
2099 2099 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
2100 2100 meteorAux[-1] = 7
2101 2101 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
2102 2102 meteorAux[-1] = 8
2103 2103 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
2104 2104 meteorAux[-1] = 11
2105 2105
2106 2106
2107 2107 except:
2108 2108 meteorAux[-1] = 11
2109 2109
2110 2110
2111 2111 listMeteors1.append(meteorAux)
2112 2112
2113 2113 return listMeteors1
2114 2114
2115 2115 #Exponential Function
2116 2116
2117 2117 def __exponential_function(self, x, a, tau):
2118 2118 y = a*numpy.exp(-x/tau)
2119 2119 return y
2120 2120
2121 2121 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
2122 2122
2123 2123 pairslist1 = list(pairslist)
2124 2124 pairslist1.append((0,4))
2125 2125 pairslist1.append((1,3))
2126 2126 numPairs = len(pairslist1)
2127 2127 #Time Lag
2128 2128 timeLag = 45*10**-3
2129 2129 c = 3e8
2130 2130 lag = numpy.ceil(timeLag/timeInterval)
2131 2131 freq = 30.175e6
2132 2132
2133 2133 listMeteors1 = []
2134 2134
2135 2135 for i in range(len(listMeteors)):
2136 2136 meteorAux = listMeteors[i]
2137 2137 if meteorAux[-1] == 0:
2138 2138 mStart = listMeteors[i][1]
2139 2139 mPeak = listMeteors[i][2]
2140 2140 mLag = mPeak - mStart + lag
2141 2141
2142 2142 #get the volt data between the start and end times of the meteor
2143 2143 meteorVolts = listVolts[i]
2144 2144 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
2145 2145
2146 2146 #Get CCF
2147 2147 allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
2148 2148
2149 2149 #Method 2
2150 2150 slopes = numpy.zeros(numPairs)
2151 2151 time = numpy.array([-2,-1,1,2])*timeInterval
2152 2152 angAllCCF = numpy.angle(allCCFs[:,[0,4,2,3],0])
2153 2153
2154 2154 #Correct phases
2155 2155 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
2156 2156 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
2157 2157
2158 2158 if indDer[0].shape[0] > 0:
2159 2159 for i in range(indDer[0].shape[0]):
2160 2160 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
2161 2161 angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
2162 2162
2163 2163 # fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
2164 2164 for j in range(numPairs):
2165 2165 fit = stats.linregress(time, angAllCCF[j,:])
2166 2166 slopes[j] = fit[0]
2167 2167
2168 2168 #Remove Outlier
2169 2169 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
2170 2170 # slopes = numpy.delete(slopes,indOut)
2171 2171 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
2172 2172 # slopes = numpy.delete(slopes,indOut)
2173 2173
2174 2174 radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
2175 2175 radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
2176 2176 meteorAux[-2] = radialError
2177 2177 meteorAux[-3] = radialVelocity
2178 2178
2179 2179 #Setting Error
2180 2180 #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
2181 2181 if numpy.abs(radialVelocity) > 200:
2182 2182 meteorAux[-1] = 15
2183 2183 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
2184 2184 elif radialError > radialStdThresh:
2185 2185 meteorAux[-1] = 12
2186 2186
2187 2187 listMeteors1.append(meteorAux)
2188 2188 return listMeteors1
2189 2189
2190 2190 def __setNewArrays(self, listMeteors, date, heiRang):
2191 2191
2192 2192 #New arrays
2193 2193 arrayMeteors = numpy.array(listMeteors)
2194 2194 arrayParameters = numpy.zeros((len(listMeteors), 13))
2195 2195
2196 2196 #Date inclusion
2197 2197 # date = re.findall(r'\((.*?)\)', date)
2198 2198 # date = date[0].split(',')
2199 2199 # date = map(int, date)
2200 2200 #
2201 2201 # if len(date)<6:
2202 2202 # date.append(0)
2203 2203 #
2204 2204 # date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
2205 2205 # arrayDate = numpy.tile(date, (len(listMeteors), 1))
2206 2206 arrayDate = numpy.tile(date, (len(listMeteors)))
2207 2207
2208 2208 #Meteor array
2209 2209 # arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
2210 2210 # arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
2211 2211
2212 2212 #Parameters Array
2213 2213 arrayParameters[:,0] = arrayDate #Date
2214 2214 arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range
2215 2215 arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error
2216 2216 arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases
2217 2217 arrayParameters[:,-1] = arrayMeteors[:,-1] #Error
2218 2218
2219 2219
2220 2220 return arrayParameters
2221 2221
2222 2222 class CorrectSMPhases(Operation):
2223 2223
2224 2224 def run(self, dataOut, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None):
2225 2225
2226 2226 arrayParameters = dataOut.data_param
2227 2227 pairsList = []
2228 2228 pairx = (0,1)
2229 2229 pairy = (2,3)
2230 2230 pairsList.append(pairx)
2231 2231 pairsList.append(pairy)
2232 2232 jph = numpy.zeros(4)
2233 2233
2234 2234 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
2235 2235 # arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
2236 2236 arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))
2237 2237
2238 2238 meteorOps = SMOperations()
2239 2239 if channelPositions is None:
2240 2240 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
2241 2241 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
2242 2242
2243 2243 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
2244 2244 h = (hmin,hmax)
2245 2245
2246 2246 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
2247 2247
2248 2248 dataOut.data_param = arrayParameters
2249 2249 return
2250 2250
2251 2251 class SMPhaseCalibration(Operation):
2252 2252
2253 2253 __buffer = None
2254 2254
2255 2255 __initime = None
2256 2256
2257 2257 __dataReady = False
2258 2258
2259 2259 __isConfig = False
2260 2260
2261 2261 def __checkTime(self, currentTime, initTime, paramInterval, outputInterval):
2262 2262
2263 2263 dataTime = currentTime + paramInterval
2264 2264 deltaTime = dataTime - initTime
2265 2265
2266 2266 if deltaTime >= outputInterval or deltaTime < 0:
2267 2267 return True
2268 2268
2269 2269 return False
2270 2270
2271 2271 def __getGammas(self, pairs, d, phases):
2272 2272 gammas = numpy.zeros(2)
2273 2273
2274 2274 for i in range(len(pairs)):
2275 2275
2276 2276 pairi = pairs[i]
2277 2277
2278 2278 phip3 = phases[:,pairi[0]]
2279 2279 d3 = d[pairi[0]]
2280 2280 phip2 = phases[:,pairi[1]]
2281 2281 d2 = d[pairi[1]]
2282 2282 #Calculating gamma
2283 2283 # jdcos = alp1/(k*d1)
2284 2284 # jgamma = numpy.angle(numpy.exp(1j*(d0*alp1/d1 - alp0)))
2285 2285 jgamma = -phip2*d3/d2 - phip3
2286 2286 jgamma = numpy.angle(numpy.exp(1j*jgamma))
2287 2287 # jgamma[jgamma>numpy.pi] -= 2*numpy.pi
2288 2288 # jgamma[jgamma<-numpy.pi] += 2*numpy.pi
2289 2289
2290 2290 #Revised distribution
2291 2291 jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi))
2292 2292
2293 2293 #Histogram
2294 nBins = 64.0
2294 nBins = 64
2295 2295 rmin = -0.5*numpy.pi
2296 2296 rmax = 0.5*numpy.pi
2297 2297 phaseHisto = numpy.histogram(jgammaArray, bins=nBins, range=(rmin,rmax))
2298 2298
2299 2299 meteorsY = phaseHisto[0]
2300 2300 phasesX = phaseHisto[1][:-1]
2301 2301 width = phasesX[1] - phasesX[0]
2302 2302 phasesX += width/2
2303 2303
2304 2304 #Gaussian aproximation
2305 2305 bpeak = meteorsY.argmax()
2306 2306 peak = meteorsY.max()
2307 2307 jmin = bpeak - 5
2308 2308 jmax = bpeak + 5 + 1
2309 2309
2310 2310 if jmin<0:
2311 2311 jmin = 0
2312 2312 jmax = 6
2313 2313 elif jmax > meteorsY.size:
2314 2314 jmin = meteorsY.size - 6
2315 2315 jmax = meteorsY.size
2316 2316
2317 2317 x0 = numpy.array([peak,bpeak,50])
2318 2318 coeff = optimize.leastsq(self.__residualFunction, x0, args=(meteorsY[jmin:jmax], phasesX[jmin:jmax]))
2319 2319
2320 2320 #Gammas
2321 2321 gammas[i] = coeff[0][1]
2322 2322
2323 2323 return gammas
2324 2324
2325 2325 def __residualFunction(self, coeffs, y, t):
2326 2326
2327 2327 return y - self.__gauss_function(t, coeffs)
2328 2328
2329 2329 def __gauss_function(self, t, coeffs):
2330 2330
2331 2331 return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2)
2332 2332
2333 2333 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
2334 2334 meteorOps = SMOperations()
2335 2335 nchan = 4
2336 2336 pairx = pairsList[0] #x es 0
2337 2337 pairy = pairsList[1] #y es 1
2338 2338 center_xangle = 0
2339 2339 center_yangle = 0
2340 2340 range_angle = numpy.array([10*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4])
2341 2341 ntimes = len(range_angle)
2342 2342
2343 nstepsx = 20.0
2344 nstepsy = 20.0
2343 nstepsx = 20
2344 nstepsy = 20
2345 2345
2346 2346 for iz in range(ntimes):
2347 2347 min_xangle = -range_angle[iz]/2 + center_xangle
2348 2348 max_xangle = range_angle[iz]/2 + center_xangle
2349 2349 min_yangle = -range_angle[iz]/2 + center_yangle
2350 2350 max_yangle = range_angle[iz]/2 + center_yangle
2351 2351
2352 2352 inc_x = (max_xangle-min_xangle)/nstepsx
2353 2353 inc_y = (max_yangle-min_yangle)/nstepsy
2354 2354
2355 2355 alpha_y = numpy.arange(nstepsy)*inc_y + min_yangle
2356 2356 alpha_x = numpy.arange(nstepsx)*inc_x + min_xangle
2357 2357 penalty = numpy.zeros((nstepsx,nstepsy))
2358 2358 jph_array = numpy.zeros((nchan,nstepsx,nstepsy))
2359 2359 jph = numpy.zeros(nchan)
2360 2360
2361 2361 # Iterations looking for the offset
2362 2362 for iy in range(int(nstepsy)):
2363 2363 for ix in range(int(nstepsx)):
2364 2364 d3 = d[pairsList[1][0]]
2365 2365 d2 = d[pairsList[1][1]]
2366 2366 d5 = d[pairsList[0][0]]
2367 2367 d4 = d[pairsList[0][1]]
2368 2368
2369 2369 alp2 = alpha_y[iy] #gamma 1
2370 2370 alp4 = alpha_x[ix] #gamma 0
2371 2371
2372 2372 alp3 = -alp2*d3/d2 - gammas[1]
2373 2373 alp5 = -alp4*d5/d4 - gammas[0]
2374 2374 # jph[pairy[1]] = alpha_y[iy]
2375 2375 # jph[pairy[0]] = -gammas[1] - alpha_y[iy]*d[pairy[1]]/d[pairy[0]]
2376 2376
2377 2377 # jph[pairx[1]] = alpha_x[ix]
2378 2378 # jph[pairx[0]] = -gammas[0] - alpha_x[ix]*d[pairx[1]]/d[pairx[0]]
2379 2379 jph[pairsList[0][1]] = alp4
2380 2380 jph[pairsList[0][0]] = alp5
2381 2381 jph[pairsList[1][0]] = alp3
2382 2382 jph[pairsList[1][1]] = alp2
2383 2383 jph_array[:,ix,iy] = jph
2384 2384 # d = [2.0,2.5,2.5,2.0]
2385 2385 #falta chequear si va a leer bien los meteoros
2386 2386 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, d, jph)
2387 2387 error = meteorsArray1[:,-1]
2388 2388 ind1 = numpy.where(error==0)[0]
2389 2389 penalty[ix,iy] = ind1.size
2390 2390
2391 2391 i,j = numpy.unravel_index(penalty.argmax(), penalty.shape)
2392 2392 phOffset = jph_array[:,i,j]
2393 2393
2394 2394 center_xangle = phOffset[pairx[1]]
2395 2395 center_yangle = phOffset[pairy[1]]
2396 2396
2397 2397 phOffset = numpy.angle(numpy.exp(1j*jph_array[:,i,j]))
2398 2398 phOffset = phOffset*180/numpy.pi
2399 2399 return phOffset
2400 2400
2401 2401
2402 2402 def run(self, dataOut, hmin, hmax, channelPositions=None, nHours = 1):
2403 2403
2404 2404 dataOut.flagNoData = True
2405 2405 self.__dataReady = False
2406 2406 dataOut.outputInterval = nHours*3600
2407 2407
2408 2408 if self.__isConfig == False:
2409 2409 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
2410 2410 #Get Initial LTC time
2411 2411 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
2412 2412 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
2413 2413
2414 2414 self.__isConfig = True
2415 2415
2416 2416 if self.__buffer is None:
2417 2417 self.__buffer = dataOut.data_param.copy()
2418 2418
2419 2419 else:
2420 2420 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
2421 2421
2422 2422 self.__dataReady = self.__checkTime(dataOut.utctime, self.__initime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
2423 2423
2424 2424 if self.__dataReady:
2425 2425 dataOut.utctimeInit = self.__initime
2426 2426 self.__initime += dataOut.outputInterval #to erase time offset
2427 2427
2428 2428 freq = dataOut.frequency
2429 2429 c = dataOut.C #m/s
2430 2430 lamb = c/freq
2431 2431 k = 2*numpy.pi/lamb
2432 2432 azimuth = 0
2433 2433 h = (hmin, hmax)
2434 2434 # pairs = ((0,1),(2,3)) #Estrella
2435 2435 # pairs = ((1,0),(2,3)) #T
2436 2436
2437 2437 if channelPositions is None:
2438 2438 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
2439 2439 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
2440 2440 meteorOps = SMOperations()
2441 2441 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
2442 2442
2443 2443 #Checking correct order of pairs
2444 2444 pairs = []
2445 2445 if distances[1] > distances[0]:
2446 2446 pairs.append((1,0))
2447 2447 else:
2448 2448 pairs.append((0,1))
2449 2449
2450 2450 if distances[3] > distances[2]:
2451 2451 pairs.append((3,2))
2452 2452 else:
2453 2453 pairs.append((2,3))
2454 2454 # distances1 = [-distances[0]*lamb, distances[1]*lamb, -distances[2]*lamb, distances[3]*lamb]
2455 2455
2456 2456 meteorsArray = self.__buffer
2457 2457 error = meteorsArray[:,-1]
2458 2458 boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14)
2459 2459 ind1 = numpy.where(boolError)[0]
2460 2460 meteorsArray = meteorsArray[ind1,:]
2461 2461 meteorsArray[:,-1] = 0
2462 2462 phases = meteorsArray[:,8:12]
2463 2463
2464 2464 #Calculate Gammas
2465 2465 gammas = self.__getGammas(pairs, distances, phases)
2466 2466 # gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180
2467 2467 #Calculate Phases
2468 2468 phasesOff = self.__getPhases(azimuth, h, pairs, distances, gammas, meteorsArray)
2469 2469 phasesOff = phasesOff.reshape((1,phasesOff.size))
2470 2470 dataOut.data_output = -phasesOff
2471 2471 dataOut.flagNoData = False
2472 2472 dataOut.channelList = pairslist0
2473 2473 self.__buffer = None
2474 2474
2475 2475
2476 2476 return
2477 2477
2478 2478 class SMOperations():
2479 2479
2480 2480 def __init__(self):
2481 2481
2482 2482 return
2483 2483
2484 2484 def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, distances, jph):
2485 2485
2486 2486 arrayParameters = arrayParameters0.copy()
2487 2487 hmin = h[0]
2488 2488 hmax = h[1]
2489 2489
2490 2490 #Calculate AOA (Error N 3, 4)
2491 2491 #JONES ET AL. 1998
2492 2492 AOAthresh = numpy.pi/8
2493 2493 error = arrayParameters[:,-1]
2494 2494 phases = -arrayParameters[:,8:12] + jph
2495 2495 # phases = numpy.unwrap(phases)
2496 2496 arrayParameters[:,3:6], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, distances, error, AOAthresh, azimuth)
2497 2497
2498 2498 #Calculate Heights (Error N 13 and 14)
2499 2499 error = arrayParameters[:,-1]
2500 2500 Ranges = arrayParameters[:,1]
2501 2501 zenith = arrayParameters[:,4]
2502 2502 arrayParameters[:,2], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
2503 2503
2504 2504 #----------------------- Get Final data ------------------------------------
2505 2505 # error = arrayParameters[:,-1]
2506 2506 # ind1 = numpy.where(error==0)[0]
2507 2507 # arrayParameters = arrayParameters[ind1,:]
2508 2508
2509 2509 return arrayParameters
2510 2510
2511 2511 def __getAOA(self, phases, pairsList, directions, error, AOAthresh, azimuth):
2512 2512
2513 2513 arrayAOA = numpy.zeros((phases.shape[0],3))
2514 2514 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList,directions)
2515 2515
2516 2516 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
2517 2517 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
2518 2518 arrayAOA[:,2] = cosDirError
2519 2519
2520 2520 azimuthAngle = arrayAOA[:,0]
2521 2521 zenithAngle = arrayAOA[:,1]
2522 2522
2523 2523 #Setting Error
2524 2524 indError = numpy.where(numpy.logical_or(error == 3, error == 4))[0]
2525 2525 error[indError] = 0
2526 2526 #Number 3: AOA not fesible
2527 2527 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
2528 2528 error[indInvalid] = 3
2529 2529 #Number 4: Large difference in AOAs obtained from different antenna baselines
2530 2530 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
2531 2531 error[indInvalid] = 4
2532 2532 return arrayAOA, error
2533 2533
2534 2534 def __getDirectionCosines(self, arrayPhase, pairsList, distances):
2535 2535
2536 2536 #Initializing some variables
2537 2537 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
2538 2538 ang_aux = ang_aux.reshape(1,ang_aux.size)
2539 2539
2540 2540 cosdir = numpy.zeros((arrayPhase.shape[0],2))
2541 2541 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
2542 2542
2543 2543
2544 2544 for i in range(2):
2545 2545 ph0 = arrayPhase[:,pairsList[i][0]]
2546 2546 ph1 = arrayPhase[:,pairsList[i][1]]
2547 2547 d0 = distances[pairsList[i][0]]
2548 2548 d1 = distances[pairsList[i][1]]
2549 2549
2550 2550 ph0_aux = ph0 + ph1
2551 2551 ph0_aux = numpy.angle(numpy.exp(1j*ph0_aux))
2552 2552 # ph0_aux[ph0_aux > numpy.pi] -= 2*numpy.pi
2553 2553 # ph0_aux[ph0_aux < -numpy.pi] += 2*numpy.pi
2554 2554 #First Estimation
2555 2555 cosdir0[:,i] = (ph0_aux)/(2*numpy.pi*(d0 - d1))
2556 2556
2557 2557 #Most-Accurate Second Estimation
2558 2558 phi1_aux = ph0 - ph1
2559 2559 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
2560 2560 #Direction Cosine 1
2561 2561 cosdir1 = (phi1_aux + ang_aux)/(2*numpy.pi*(d0 + d1))
2562 2562
2563 2563 #Searching the correct Direction Cosine
2564 2564 cosdir0_aux = cosdir0[:,i]
2565 2565 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
2566 2566 #Minimum Distance
2567 2567 cosDiff = (cosdir1 - cosdir0_aux)**2
2568 2568 indcos = cosDiff.argmin(axis = 1)
2569 2569 #Saving Value obtained
2570 2570 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
2571 2571
2572 2572 return cosdir0, cosdir
2573 2573
2574 2574 def __calculateAOA(self, cosdir, azimuth):
2575 2575 cosdirX = cosdir[:,0]
2576 2576 cosdirY = cosdir[:,1]
2577 2577
2578 2578 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
2579 2579 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth#0 deg north, 90 deg east
2580 2580 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
2581 2581
2582 2582 return angles
2583 2583
2584 2584 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
2585 2585
2586 2586 Ramb = 375 #Ramb = c/(2*PRF)
2587 2587 Re = 6371 #Earth Radius
2588 2588 heights = numpy.zeros(Ranges.shape)
2589 2589
2590 2590 R_aux = numpy.array([0,1,2])*Ramb
2591 2591 R_aux = R_aux.reshape(1,R_aux.size)
2592 2592
2593 2593 Ranges = Ranges.reshape(Ranges.size,1)
2594 2594
2595 2595 Ri = Ranges + R_aux
2596 2596 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
2597 2597
2598 2598 #Check if there is a height between 70 and 110 km
2599 2599 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
2600 2600 ind_h = numpy.where(h_bool == 1)[0]
2601 2601
2602 2602 hCorr = hi[ind_h, :]
2603 2603 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
2604 2604
2605 hCorr = hi[ind_hCorr]
2605 hCorr = hi[ind_hCorr][:len(ind_h)]
2606 2606 heights[ind_h] = hCorr
2607 2607
2608 2608 #Setting Error
2609 2609 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
2610 2610 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
2611 2611 indError = numpy.where(numpy.logical_or(error == 13, error == 14))[0]
2612 2612 error[indError] = 0
2613 2613 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
2614 2614 error[indInvalid2] = 14
2615 2615 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
2616 2616 error[indInvalid1] = 13
2617 2617
2618 2618 return heights, error
2619 2619
2620 2620 def getPhasePairs(self, channelPositions):
2621 2621 chanPos = numpy.array(channelPositions)
2622 2622 listOper = list(itertools.combinations(range(5),2))
2623 2623
2624 2624 distances = numpy.zeros(4)
2625 2625 axisX = []
2626 2626 axisY = []
2627 2627 distX = numpy.zeros(3)
2628 2628 distY = numpy.zeros(3)
2629 2629 ix = 0
2630 2630 iy = 0
2631 2631
2632 2632 pairX = numpy.zeros((2,2))
2633 2633 pairY = numpy.zeros((2,2))
2634 2634
2635 2635 for i in range(len(listOper)):
2636 2636 pairi = listOper[i]
2637 2637
2638 2638 posDif = numpy.abs(chanPos[pairi[0],:] - chanPos[pairi[1],:])
2639 2639
2640 2640 if posDif[0] == 0:
2641 2641 axisY.append(pairi)
2642 2642 distY[iy] = posDif[1]
2643 2643 iy += 1
2644 2644 elif posDif[1] == 0:
2645 2645 axisX.append(pairi)
2646 2646 distX[ix] = posDif[0]
2647 2647 ix += 1
2648 2648
2649 2649 for i in range(2):
2650 2650 if i==0:
2651 2651 dist0 = distX
2652 2652 axis0 = axisX
2653 2653 else:
2654 2654 dist0 = distY
2655 2655 axis0 = axisY
2656 2656
2657 2657 side = numpy.argsort(dist0)[:-1]
2658 2658 axis0 = numpy.array(axis0)[side,:]
2659 2659 chanC = int(numpy.intersect1d(axis0[0,:], axis0[1,:])[0])
2660 2660 axis1 = numpy.unique(numpy.reshape(axis0,4))
2661 2661 side = axis1[axis1 != chanC]
2662 2662 diff1 = chanPos[chanC,i] - chanPos[side[0],i]
2663 2663 diff2 = chanPos[chanC,i] - chanPos[side[1],i]
2664 2664 if diff1<0:
2665 2665 chan2 = side[0]
2666 2666 d2 = numpy.abs(diff1)
2667 2667 chan1 = side[1]
2668 2668 d1 = numpy.abs(diff2)
2669 2669 else:
2670 2670 chan2 = side[1]
2671 2671 d2 = numpy.abs(diff2)
2672 2672 chan1 = side[0]
2673 2673 d1 = numpy.abs(diff1)
2674 2674
2675 2675 if i==0:
2676 2676 chanCX = chanC
2677 2677 chan1X = chan1
2678 2678 chan2X = chan2
2679 2679 distances[0:2] = numpy.array([d1,d2])
2680 2680 else:
2681 2681 chanCY = chanC
2682 2682 chan1Y = chan1
2683 2683 chan2Y = chan2
2684 2684 distances[2:4] = numpy.array([d1,d2])
2685 2685 # axisXsides = numpy.reshape(axisX[ix,:],4)
2686 2686 #
2687 2687 # channelCentX = int(numpy.intersect1d(pairX[0,:], pairX[1,:])[0])
2688 2688 # channelCentY = int(numpy.intersect1d(pairY[0,:], pairY[1,:])[0])
2689 2689 #
2690 2690 # ind25X = numpy.where(pairX[0,:] != channelCentX)[0][0]
2691 2691 # ind20X = numpy.where(pairX[1,:] != channelCentX)[0][0]
2692 2692 # channel25X = int(pairX[0,ind25X])
2693 2693 # channel20X = int(pairX[1,ind20X])
2694 2694 # ind25Y = numpy.where(pairY[0,:] != channelCentY)[0][0]
2695 2695 # ind20Y = numpy.where(pairY[1,:] != channelCentY)[0][0]
2696 2696 # channel25Y = int(pairY[0,ind25Y])
2697 2697 # channel20Y = int(pairY[1,ind20Y])
2698 2698
2699 2699 # pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)]
2700 2700 pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]
2701 2701
2702 2702 return pairslist, distances
2703 2703 # def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
2704 2704 #
2705 2705 # arrayAOA = numpy.zeros((phases.shape[0],3))
2706 2706 # cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
2707 2707 #
2708 2708 # arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
2709 2709 # cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
2710 2710 # arrayAOA[:,2] = cosDirError
2711 2711 #
2712 2712 # azimuthAngle = arrayAOA[:,0]
2713 2713 # zenithAngle = arrayAOA[:,1]
2714 2714 #
2715 2715 # #Setting Error
2716 2716 # #Number 3: AOA not fesible
2717 2717 # indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
2718 2718 # error[indInvalid] = 3
2719 2719 # #Number 4: Large difference in AOAs obtained from different antenna baselines
2720 2720 # indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
2721 2721 # error[indInvalid] = 4
2722 2722 # return arrayAOA, error
2723 2723 #
2724 2724 # def __getDirectionCosines(self, arrayPhase, pairsList):
2725 2725 #
2726 2726 # #Initializing some variables
2727 2727 # ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
2728 2728 # ang_aux = ang_aux.reshape(1,ang_aux.size)
2729 2729 #
2730 2730 # cosdir = numpy.zeros((arrayPhase.shape[0],2))
2731 2731 # cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
2732 2732 #
2733 2733 #
2734 2734 # for i in range(2):
2735 2735 # #First Estimation
2736 2736 # phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
2737 2737 # #Dealias
2738 2738 # indcsi = numpy.where(phi0_aux > numpy.pi)
2739 2739 # phi0_aux[indcsi] -= 2*numpy.pi
2740 2740 # indcsi = numpy.where(phi0_aux < -numpy.pi)
2741 2741 # phi0_aux[indcsi] += 2*numpy.pi
2742 2742 # #Direction Cosine 0
2743 2743 # cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
2744 2744 #
2745 2745 # #Most-Accurate Second Estimation
2746 2746 # phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
2747 2747 # phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
2748 2748 # #Direction Cosine 1
2749 2749 # cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
2750 2750 #
2751 2751 # #Searching the correct Direction Cosine
2752 2752 # cosdir0_aux = cosdir0[:,i]
2753 2753 # cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
2754 2754 # #Minimum Distance
2755 2755 # cosDiff = (cosdir1 - cosdir0_aux)**2
2756 2756 # indcos = cosDiff.argmin(axis = 1)
2757 2757 # #Saving Value obtained
2758 2758 # cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
2759 2759 #
2760 2760 # return cosdir0, cosdir
2761 2761 #
2762 2762 # def __calculateAOA(self, cosdir, azimuth):
2763 2763 # cosdirX = cosdir[:,0]
2764 2764 # cosdirY = cosdir[:,1]
2765 2765 #
2766 2766 # zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
2767 2767 # azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
2768 2768 # angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
2769 2769 #
2770 2770 # return angles
2771 2771 #
2772 2772 # def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
2773 2773 #
2774 2774 # Ramb = 375 #Ramb = c/(2*PRF)
2775 2775 # Re = 6371 #Earth Radius
2776 2776 # heights = numpy.zeros(Ranges.shape)
2777 2777 #
2778 2778 # R_aux = numpy.array([0,1,2])*Ramb
2779 2779 # R_aux = R_aux.reshape(1,R_aux.size)
2780 2780 #
2781 2781 # Ranges = Ranges.reshape(Ranges.size,1)
2782 2782 #
2783 2783 # Ri = Ranges + R_aux
2784 2784 # hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
2785 2785 #
2786 2786 # #Check if there is a height between 70 and 110 km
2787 2787 # h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
2788 2788 # ind_h = numpy.where(h_bool == 1)[0]
2789 2789 #
2790 2790 # hCorr = hi[ind_h, :]
2791 2791 # ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
2792 2792 #
2793 2793 # hCorr = hi[ind_hCorr]
2794 2794 # heights[ind_h] = hCorr
2795 2795 #
2796 2796 # #Setting Error
2797 2797 # #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
2798 2798 # #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
2799 2799 #
2800 2800 # indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
2801 2801 # error[indInvalid2] = 14
2802 2802 # indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
2803 2803 # error[indInvalid1] = 13
2804 2804 #
2805 2805 # return heights, error
@@ -1,1283 +1,1285
1 1 import sys
2 2 import numpy
3 3 from scipy import interpolate
4 4
5 5 from jroproc_base import ProcessingUnit, Operation
6 6 from schainpy.model.data.jrodata import Voltage
7 7
8 8 class VoltageProc(ProcessingUnit):
9 9
10 10
11 11 def __init__(self, **kwargs):
12 12
13 13 ProcessingUnit.__init__(self, **kwargs)
14 14
15 15 # self.objectDict = {}
16 16 self.dataOut = Voltage()
17 17 self.flip = 1
18 18
19 19 def run(self):
20 20 if self.dataIn.type == 'AMISR':
21 21 self.__updateObjFromAmisrInput()
22 22
23 23 if self.dataIn.type == 'Voltage':
24 24 self.dataOut.copy(self.dataIn)
25 25
26 26 # self.dataOut.copy(self.dataIn)
27 27
28 28 def __updateObjFromAmisrInput(self):
29 29
30 30 self.dataOut.timeZone = self.dataIn.timeZone
31 31 self.dataOut.dstFlag = self.dataIn.dstFlag
32 32 self.dataOut.errorCount = self.dataIn.errorCount
33 33 self.dataOut.useLocalTime = self.dataIn.useLocalTime
34 34
35 35 self.dataOut.flagNoData = self.dataIn.flagNoData
36 36 self.dataOut.data = self.dataIn.data
37 37 self.dataOut.utctime = self.dataIn.utctime
38 38 self.dataOut.channelList = self.dataIn.channelList
39 39 # self.dataOut.timeInterval = self.dataIn.timeInterval
40 40 self.dataOut.heightList = self.dataIn.heightList
41 41 self.dataOut.nProfiles = self.dataIn.nProfiles
42 42
43 43 self.dataOut.nCohInt = self.dataIn.nCohInt
44 44 self.dataOut.ippSeconds = self.dataIn.ippSeconds
45 45 self.dataOut.frequency = self.dataIn.frequency
46 46
47 47 self.dataOut.azimuth = self.dataIn.azimuth
48 48 self.dataOut.zenith = self.dataIn.zenith
49 49
50 50 self.dataOut.beam.codeList = self.dataIn.beam.codeList
51 51 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
52 52 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
53 53 #
54 54 # pass#
55 55 #
56 56 # def init(self):
57 57 #
58 58 #
59 59 # if self.dataIn.type == 'AMISR':
60 60 # self.__updateObjFromAmisrInput()
61 61 #
62 62 # if self.dataIn.type == 'Voltage':
63 63 # self.dataOut.copy(self.dataIn)
64 64 # # No necesita copiar en cada init() los atributos de dataIn
65 65 # # la copia deberia hacerse por cada nuevo bloque de datos
66 66
67 67 def selectChannels(self, channelList):
68 68
69 69 channelIndexList = []
70 70
71 71 for channel in channelList:
72 72 if channel not in self.dataOut.channelList:
73 73 raise ValueError, "Channel %d is not in %s" %(channel, str(self.dataOut.channelList))
74 74
75 75 index = self.dataOut.channelList.index(channel)
76 76 channelIndexList.append(index)
77 77
78 78 self.selectChannelsByIndex(channelIndexList)
79 79
80 80 def selectChannelsByIndex(self, channelIndexList):
81 81 """
82 82 Selecciona un bloque de datos en base a canales segun el channelIndexList
83 83
84 84 Input:
85 85 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
86 86
87 87 Affected:
88 88 self.dataOut.data
89 89 self.dataOut.channelIndexList
90 90 self.dataOut.nChannels
91 91 self.dataOut.m_ProcessingHeader.totalSpectra
92 92 self.dataOut.systemHeaderObj.numChannels
93 93 self.dataOut.m_ProcessingHeader.blockSize
94 94
95 95 Return:
96 96 None
97 97 """
98 98
99 99 for channelIndex in channelIndexList:
100 100 if channelIndex not in self.dataOut.channelIndexList:
101 101 print channelIndexList
102 102 raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
103 103
104 104 if self.dataOut.flagDataAsBlock:
105 105 """
106 106 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
107 107 """
108 108 data = self.dataOut.data[channelIndexList,:,:]
109 109 else:
110 110 data = self.dataOut.data[channelIndexList,:]
111 111
112 112 self.dataOut.data = data
113 113 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
114 114 # self.dataOut.nChannels = nChannels
115 115
116 116 return 1
117 117
118 118 def selectHeights(self, minHei=None, maxHei=None):
119 119 """
120 120 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
121 121 minHei <= height <= maxHei
122 122
123 123 Input:
124 124 minHei : valor minimo de altura a considerar
125 125 maxHei : valor maximo de altura a considerar
126 126
127 127 Affected:
128 128 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
129 129
130 130 Return:
131 131 1 si el metodo se ejecuto con exito caso contrario devuelve 0
132 132 """
133 133
134 134 if minHei == None:
135 135 minHei = self.dataOut.heightList[0]
136 136
137 137 if maxHei == None:
138 138 maxHei = self.dataOut.heightList[-1]
139 139
140 140 if (minHei < self.dataOut.heightList[0]):
141 141 minHei = self.dataOut.heightList[0]
142 142
143 143 if (maxHei > self.dataOut.heightList[-1]):
144 144 maxHei = self.dataOut.heightList[-1]
145 145
146 146 minIndex = 0
147 147 maxIndex = 0
148 148 heights = self.dataOut.heightList
149 149
150 150 inda = numpy.where(heights >= minHei)
151 151 indb = numpy.where(heights <= maxHei)
152 152
153 153 try:
154 154 minIndex = inda[0][0]
155 155 except:
156 156 minIndex = 0
157 157
158 158 try:
159 159 maxIndex = indb[0][-1]
160 160 except:
161 161 maxIndex = len(heights)
162 162
163 163 self.selectHeightsByIndex(minIndex, maxIndex)
164 164
165 165 return 1
166 166
167 167
168 168 def selectHeightsByIndex(self, minIndex, maxIndex):
169 169 """
170 170 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
171 171 minIndex <= index <= maxIndex
172 172
173 173 Input:
174 174 minIndex : valor de indice minimo de altura a considerar
175 175 maxIndex : valor de indice maximo de altura a considerar
176 176
177 177 Affected:
178 178 self.dataOut.data
179 179 self.dataOut.heightList
180 180
181 181 Return:
182 182 1 si el metodo se ejecuto con exito caso contrario devuelve 0
183 183 """
184 184
185 185 if (minIndex < 0) or (minIndex > maxIndex):
186 186 raise ValueError, "Height index range (%d,%d) is not valid" % (minIndex, maxIndex)
187 187
188 188 if (maxIndex >= self.dataOut.nHeights):
189 189 maxIndex = self.dataOut.nHeights
190 190
191 191 #voltage
192 192 if self.dataOut.flagDataAsBlock:
193 193 """
194 194 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
195 195 """
196 196 data = self.dataOut.data[:,:, minIndex:maxIndex]
197 197 else:
198 198 data = self.dataOut.data[:, minIndex:maxIndex]
199 199
200 200 # firstHeight = self.dataOut.heightList[minIndex]
201 201
202 202 self.dataOut.data = data
203 203 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
204 204
205 205 if self.dataOut.nHeights <= 1:
206 206 raise ValueError, "selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights)
207 207
208 208 return 1
209 209
210 210
211 211 def filterByHeights(self, window):
212 212
213 213 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
214 214
215 215 if window == None:
216 216 window = (self.dataOut.radarControllerHeaderObj.txA/self.dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
217 217
218 218 newdelta = deltaHeight * window
219 219 r = self.dataOut.nHeights % window
220 220 newheights = (self.dataOut.nHeights-r)/window
221 221
222 222 if newheights <= 1:
223 223 raise ValueError, "filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(self.dataOut.nHeights, window)
224 224
225 225 if self.dataOut.flagDataAsBlock:
226 226 """
227 227 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
228 228 """
229 229 buffer = self.dataOut.data[:, :, 0:self.dataOut.nHeights-r]
230 230 buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nProfiles,self.dataOut.nHeights/window,window)
231 231 buffer = numpy.sum(buffer,3)
232 232
233 233 else:
234 234 buffer = self.dataOut.data[:,0:self.dataOut.nHeights-r]
235 235 buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nHeights/window,window)
236 236 buffer = numpy.sum(buffer,2)
237 237
238 238 self.dataOut.data = buffer
239 239 self.dataOut.heightList = self.dataOut.heightList[0] + numpy.arange( newheights )*newdelta
240 240 self.dataOut.windowOfFilter = window
241 241
242 242 def setH0(self, h0, deltaHeight = None):
243 243
244 244 if not deltaHeight:
245 245 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
246 246
247 247 nHeights = self.dataOut.nHeights
248 248
249 249 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
250 250
251 251 self.dataOut.heightList = newHeiRange
252 252
253 253 def deFlip(self, channelList = []):
254 254
255 255 data = self.dataOut.data.copy()
256 256
257 257 if self.dataOut.flagDataAsBlock:
258 258 flip = self.flip
259 259 profileList = range(self.dataOut.nProfiles)
260 260
261 261 if not channelList:
262 262 for thisProfile in profileList:
263 263 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
264 264 flip *= -1.0
265 265 else:
266 266 for thisChannel in channelList:
267 267 if thisChannel not in self.dataOut.channelList:
268 268 continue
269 269
270 270 for thisProfile in profileList:
271 271 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
272 272 flip *= -1.0
273 273
274 274 self.flip = flip
275 275
276 276 else:
277 277 if not channelList:
278 278 data[:,:] = data[:,:]*self.flip
279 279 else:
280 280 for thisChannel in channelList:
281 281 if thisChannel not in self.dataOut.channelList:
282 282 continue
283 283
284 284 data[thisChannel,:] = data[thisChannel,:]*self.flip
285 285
286 286 self.flip *= -1.
287 287
288 288 self.dataOut.data = data
289 289
290 290 def setRadarFrequency(self, frequency=None):
291 291
292 292 if frequency != None:
293 293 self.dataOut.frequency = frequency
294 294
295 295 return 1
296 296
297 297 def interpolateHeights(self, topLim, botLim):
298 298 #69 al 72 para julia
299 299 #82-84 para meteoros
300 300 if len(numpy.shape(self.dataOut.data))==2:
301 301 sampInterp = (self.dataOut.data[:,botLim-1] + self.dataOut.data[:,topLim+1])/2
302 302 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
303 303 #self.dataOut.data[:,botLim:limSup+1] = sampInterp
304 304 self.dataOut.data[:,botLim:topLim+1] = sampInterp
305 305 else:
306 306 nHeights = self.dataOut.data.shape[2]
307 307 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
308 308 y = self.dataOut.data[:,:,range(botLim)+range(topLim+1,nHeights)]
309 309 f = interpolate.interp1d(x, y, axis = 2)
310 310 xnew = numpy.arange(botLim,topLim+1)
311 311 ynew = f(xnew)
312 312
313 313 self.dataOut.data[:,:,botLim:topLim+1] = ynew
314 314
315 315 # import collections
316 316
317 317 class CohInt(Operation):
318 318
319 319 isConfig = False
320 320
321 321 __profIndex = 0
322 322 __withOverapping = False
323 323
324 324 __byTime = False
325 325 __initime = None
326 326 __lastdatatime = None
327 327 __integrationtime = None
328 328
329 329 __buffer = None
330 330
331 331 __dataReady = False
332 332
333 333 n = None
334 334
335 335
336 336 def __init__(self, **kwargs):
337 337
338 338 Operation.__init__(self, **kwargs)
339 339
340 340 # self.isConfig = False
341 341
342 342 def setup(self, n=None, timeInterval=None, overlapping=False, byblock=False):
343 343 """
344 344 Set the parameters of the integration class.
345 345
346 346 Inputs:
347 347
348 348 n : Number of coherent integrations
349 349 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
350 350 overlapping :
351 351
352 352 """
353 353
354 354 self.__initime = None
355 355 self.__lastdatatime = 0
356 356 self.__buffer = None
357 357 self.__dataReady = False
358 358 self.byblock = byblock
359 359
360 360 if n == None and timeInterval == None:
361 361 raise ValueError, "n or timeInterval should be specified ..."
362 362
363 363 if n != None:
364 364 self.n = n
365 365 self.__byTime = False
366 366 else:
367 367 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
368 368 self.n = 9999
369 369 self.__byTime = True
370 370
371 371 if overlapping:
372 372 self.__withOverapping = True
373 373 self.__buffer = None
374 374 else:
375 375 self.__withOverapping = False
376 376 self.__buffer = 0
377 377
378 378 self.__profIndex = 0
379 379
380 380 def putData(self, data):
381 381
382 382 """
383 383 Add a profile to the __buffer and increase in one the __profileIndex
384 384
385 385 """
386 386
387 387 if not self.__withOverapping:
388 388 self.__buffer += data.copy()
389 389 self.__profIndex += 1
390 390 return
391 391
392 392 #Overlapping data
393 393 nChannels, nHeis = data.shape
394 394 data = numpy.reshape(data, (1, nChannels, nHeis))
395 395
396 396 #If the buffer is empty then it takes the data value
397 397 if self.__buffer is None:
398 398 self.__buffer = data
399 399 self.__profIndex += 1
400 400 return
401 401
402 402 #If the buffer length is lower than n then stakcing the data value
403 403 if self.__profIndex < self.n:
404 404 self.__buffer = numpy.vstack((self.__buffer, data))
405 405 self.__profIndex += 1
406 406 return
407 407
408 408 #If the buffer length is equal to n then replacing the last buffer value with the data value
409 409 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
410 410 self.__buffer[self.n-1] = data
411 411 self.__profIndex = self.n
412 412 return
413 413
414 414
415 415 def pushData(self):
416 416 """
417 417 Return the sum of the last profiles and the profiles used in the sum.
418 418
419 419 Affected:
420 420
421 421 self.__profileIndex
422 422
423 423 """
424 424
425 425 if not self.__withOverapping:
426 426 data = self.__buffer
427 427 n = self.__profIndex
428 428
429 429 self.__buffer = 0
430 430 self.__profIndex = 0
431 431
432 432 return data, n
433 433
434 434 #Integration with Overlapping
435 435 data = numpy.sum(self.__buffer, axis=0)
436 436 n = self.__profIndex
437 437
438 438 return data, n
439 439
440 440 def byProfiles(self, data):
441 441
442 442 self.__dataReady = False
443 443 avgdata = None
444 444 # n = None
445 445
446 446 self.putData(data)
447 447
448 448 if self.__profIndex == self.n:
449 449
450 450 avgdata, n = self.pushData()
451 451 self.__dataReady = True
452 452
453 453 return avgdata
454 454
455 455 def byTime(self, data, datatime):
456 456
457 457 self.__dataReady = False
458 458 avgdata = None
459 459 n = None
460 460
461 461 self.putData(data)
462 462
463 463 if (datatime - self.__initime) >= self.__integrationtime:
464 464 avgdata, n = self.pushData()
465 465 self.n = n
466 466 self.__dataReady = True
467 467
468 468 return avgdata
469 469
470 470 def integrate(self, data, datatime=None):
471 471
472 472 if self.__initime == None:
473 473 self.__initime = datatime
474 474
475 475 if self.__byTime:
476 476 avgdata = self.byTime(data, datatime)
477 477 else:
478 478 avgdata = self.byProfiles(data)
479 479
480 480
481 481 self.__lastdatatime = datatime
482 482
483 483 if avgdata is None:
484 484 return None, None
485 485
486 486 avgdatatime = self.__initime
487 487
488 488 deltatime = datatime -self.__lastdatatime
489 489
490 490 if not self.__withOverapping:
491 491 self.__initime = datatime
492 492 else:
493 493 self.__initime += deltatime
494 494
495 495 return avgdata, avgdatatime
496 496
497 497 def integrateByBlock(self, dataOut):
498 498
499 499 times = int(dataOut.data.shape[1]/self.n)
500 500 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
501 501
502 502 id_min = 0
503 503 id_max = self.n
504 504
505 505 for i in range(times):
506 506 junk = dataOut.data[:,id_min:id_max,:]
507 507 avgdata[:,i,:] = junk.sum(axis=1)
508 508 id_min += self.n
509 509 id_max += self.n
510 510
511 511 timeInterval = dataOut.ippSeconds*self.n
512 512 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
513 513 self.__dataReady = True
514 514 return avgdata, avgdatatime
515 515
516 516
517 517 def run(self, dataOut, n=None, timeInterval=None, overlapping=False, byblock=False, **kwargs):
518 518 if not self.isConfig:
519 519 self.setup(n=n, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
520 520 self.isConfig = True
521 521
522 522 if dataOut.flagDataAsBlock:
523 523 """
524 524 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
525 525 """
526 526 avgdata, avgdatatime = self.integrateByBlock(dataOut)
527 527 dataOut.nProfiles /= self.n
528 528 else:
529 529 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
530 530
531 531 # dataOut.timeInterval *= n
532 532 dataOut.flagNoData = True
533 533
534 534 if self.__dataReady:
535 535 dataOut.data = avgdata
536 536 dataOut.nCohInt *= self.n
537 537 dataOut.utctime = avgdatatime
538 538 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
539 539 dataOut.flagNoData = False
540 540
541 541 class Decoder(Operation):
542 542
543 543 isConfig = False
544 544 __profIndex = 0
545 545
546 546 code = None
547 547
548 548 nCode = None
549 549 nBaud = None
550 550
551 551
552 552 def __init__(self, **kwargs):
553 553
554 554 Operation.__init__(self, **kwargs)
555 555
556 556 self.times = None
557 557 self.osamp = None
558 558 # self.__setValues = False
559 559 self.isConfig = False
560 560
561 561 def setup(self, code, osamp, dataOut):
562 562
563 563 self.__profIndex = 0
564 564
565 565 self.code = code
566 566
567 567 self.nCode = len(code)
568 568 self.nBaud = len(code[0])
569 569
570 570 if (osamp != None) and (osamp >1):
571 571 self.osamp = osamp
572 572 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
573 573 self.nBaud = self.nBaud*self.osamp
574 574
575 575 self.__nChannels = dataOut.nChannels
576 576 self.__nProfiles = dataOut.nProfiles
577 577 self.__nHeis = dataOut.nHeights
578 578
579 579 if self.__nHeis < self.nBaud:
580 580 raise ValueError, 'Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud)
581 581
582 582 #Frequency
583 583 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
584 584
585 585 __codeBuffer[:,0:self.nBaud] = self.code
586 586
587 587 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
588 588
589 589 if dataOut.flagDataAsBlock:
590 590
591 591 self.ndatadec = self.__nHeis #- self.nBaud + 1
592 592
593 593 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
594 594
595 595 else:
596 596
597 597 #Time
598 598 self.ndatadec = self.__nHeis #- self.nBaud + 1
599 599
600 600 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
601 601
602 602 def __convolutionInFreq(self, data):
603 603
604 604 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
605 605
606 606 fft_data = numpy.fft.fft(data, axis=1)
607 607
608 608 conv = fft_data*fft_code
609 609
610 610 data = numpy.fft.ifft(conv,axis=1)
611 611
612 612 return data
613 613
614 614 def __convolutionInFreqOpt(self, data):
615 615
616 616 raise NotImplementedError
617 617
618 618 def __convolutionInTime(self, data):
619 619
620 620 code = self.code[self.__profIndex]
621 621
622 622 for i in range(self.__nChannels):
623 623 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
624 624
625 625 return self.datadecTime
626 626
627 627 def __convolutionByBlockInTime(self, data):
628 628
629 629 repetitions = self.__nProfiles / self.nCode
630 630
631 631 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
632 632 junk = junk.flatten()
633 633 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
634 634
635 635 for i in range(self.__nChannels):
636 636 for j in range(self.__nProfiles):
637 print self.datadecTime[i,j,:].shape
638 print numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:].shape
637 639 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
638 640
639 641 return self.datadecTime
640 642
641 643 def __convolutionByBlockInFreq(self, data):
642 644
643 645 raise NotImplementedError, "Decoder by frequency fro Blocks not implemented"
644 646
645 647
646 648 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
647 649
648 650 fft_data = numpy.fft.fft(data, axis=2)
649 651
650 652 conv = fft_data*fft_code
651 653
652 654 data = numpy.fft.ifft(conv,axis=2)
653 655
654 656 return data
655 657
656 658 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
657 659
658 660 if dataOut.flagDecodeData:
659 661 print "This data is already decoded, recoding again ..."
660 662
661 663 if not self.isConfig:
662 664
663 665 if code is None:
664 666 if dataOut.code is None:
665 667 raise ValueError, "Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type
666 668
667 669 code = dataOut.code
668 670 else:
669 671 code = numpy.array(code).reshape(nCode,nBaud)
670 672
671 673 self.setup(code, osamp, dataOut)
672 674
673 675 self.isConfig = True
674 676
675 677 if mode == 3:
676 678 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
677 679
678 680 if times != None:
679 681 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
680 682
681 683 if self.code is None:
682 684 print "Fail decoding: Code is not defined."
683 685 return
684 686
685 687 datadec = None
686 688 if mode == 3:
687 689 mode = 0
688 690
689 691 if dataOut.flagDataAsBlock:
690 692 """
691 693 Decoding when data have been read as block,
692 694 """
693 695
694 696 if mode == 0:
695 697 datadec = self.__convolutionByBlockInTime(dataOut.data)
696 698 if mode == 1:
697 699 datadec = self.__convolutionByBlockInFreq(dataOut.data)
698 700 else:
699 701 """
700 702 Decoding when data have been read profile by profile
701 703 """
702 704 if mode == 0:
703 705 datadec = self.__convolutionInTime(dataOut.data)
704 706
705 707 if mode == 1:
706 708 datadec = self.__convolutionInFreq(dataOut.data)
707 709
708 710 if mode == 2:
709 711 datadec = self.__convolutionInFreqOpt(dataOut.data)
710 712
711 713 if datadec is None:
712 714 raise ValueError, "Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode
713 715
714 716 dataOut.code = self.code
715 717 dataOut.nCode = self.nCode
716 718 dataOut.nBaud = self.nBaud
717 719
718 720 dataOut.data = datadec
719 721
720 722 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
721 723
722 724 dataOut.flagDecodeData = True #asumo q la data esta decodificada
723 725
724 726 if self.__profIndex == self.nCode-1:
725 727 self.__profIndex = 0
726 728 return 1
727 729
728 730 self.__profIndex += 1
729 731
730 732 return 1
731 733 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
732 734
733 735
734 736 class ProfileConcat(Operation):
735 737
736 738 isConfig = False
737 739 buffer = None
738 740
739 741 def __init__(self, **kwargs):
740 742
741 743 Operation.__init__(self, **kwargs)
742 744 self.profileIndex = 0
743 745
744 746 def reset(self):
745 747 self.buffer = numpy.zeros_like(self.buffer)
746 748 self.start_index = 0
747 749 self.times = 1
748 750
749 751 def setup(self, data, m, n=1):
750 752 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
751 753 self.nHeights = data.shape[1]#.nHeights
752 754 self.start_index = 0
753 755 self.times = 1
754 756
755 757 def concat(self, data):
756 758
757 759 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
758 760 self.start_index = self.start_index + self.nHeights
759 761
760 762 def run(self, dataOut, m):
761 763
762 764 dataOut.flagNoData = True
763 765
764 766 if not self.isConfig:
765 767 self.setup(dataOut.data, m, 1)
766 768 self.isConfig = True
767 769
768 770 if dataOut.flagDataAsBlock:
769 771 raise ValueError, "ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False"
770 772
771 773 else:
772 774 self.concat(dataOut.data)
773 775 self.times += 1
774 776 if self.times > m:
775 777 dataOut.data = self.buffer
776 778 self.reset()
777 779 dataOut.flagNoData = False
778 780 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
779 781 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
780 782 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
781 783 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
782 784 dataOut.ippSeconds *= m
783 785
784 786 class ProfileSelector(Operation):
785 787
786 788 profileIndex = None
787 789 # Tamanho total de los perfiles
788 790 nProfiles = None
789 791
790 792 def __init__(self, **kwargs):
791 793
792 794 Operation.__init__(self, **kwargs)
793 795 self.profileIndex = 0
794 796
795 797 def incProfileIndex(self):
796 798
797 799 self.profileIndex += 1
798 800
799 801 if self.profileIndex >= self.nProfiles:
800 802 self.profileIndex = 0
801 803
802 804 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
803 805
804 806 if profileIndex < minIndex:
805 807 return False
806 808
807 809 if profileIndex > maxIndex:
808 810 return False
809 811
810 812 return True
811 813
812 814 def isThisProfileInList(self, profileIndex, profileList):
813 815
814 816 if profileIndex not in profileList:
815 817 return False
816 818
817 819 return True
818 820
819 821 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
820 822
821 823 """
822 824 ProfileSelector:
823 825
824 826 Inputs:
825 827 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
826 828
827 829 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
828 830
829 831 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
830 832
831 833 """
832 834
833 835 if rangeList is not None:
834 836 if type(rangeList[0]) not in (tuple, list):
835 837 rangeList = [rangeList]
836 838
837 839 dataOut.flagNoData = True
838 840
839 841 if dataOut.flagDataAsBlock:
840 842 """
841 843 data dimension = [nChannels, nProfiles, nHeis]
842 844 """
843 845 if profileList != None:
844 846 dataOut.data = dataOut.data[:,profileList,:]
845 847
846 848 if profileRangeList != None:
847 849 minIndex = profileRangeList[0]
848 850 maxIndex = profileRangeList[1]
849 851 profileList = range(minIndex, maxIndex+1)
850 852
851 853 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
852 854
853 855 if rangeList != None:
854 856
855 857 profileList = []
856 858
857 859 for thisRange in rangeList:
858 860 minIndex = thisRange[0]
859 861 maxIndex = thisRange[1]
860 862
861 863 profileList.extend(range(minIndex, maxIndex+1))
862 864
863 865 dataOut.data = dataOut.data[:,profileList,:]
864 866
865 867 dataOut.nProfiles = len(profileList)
866 868 dataOut.profileIndex = dataOut.nProfiles - 1
867 869 dataOut.flagNoData = False
868 870
869 871 return True
870 872
871 873 """
872 874 data dimension = [nChannels, nHeis]
873 875 """
874 876
875 877 if profileList != None:
876 878
877 879 if self.isThisProfileInList(dataOut.profileIndex, profileList):
878 880
879 881 self.nProfiles = len(profileList)
880 882 dataOut.nProfiles = self.nProfiles
881 883 dataOut.profileIndex = self.profileIndex
882 884 dataOut.flagNoData = False
883 885
884 886 self.incProfileIndex()
885 887 return True
886 888
887 889 if profileRangeList != None:
888 890
889 891 minIndex = profileRangeList[0]
890 892 maxIndex = profileRangeList[1]
891 893
892 894 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
893 895
894 896 self.nProfiles = maxIndex - minIndex + 1
895 897 dataOut.nProfiles = self.nProfiles
896 898 dataOut.profileIndex = self.profileIndex
897 899 dataOut.flagNoData = False
898 900
899 901 self.incProfileIndex()
900 902 return True
901 903
902 904 if rangeList != None:
903 905
904 906 nProfiles = 0
905 907
906 908 for thisRange in rangeList:
907 909 minIndex = thisRange[0]
908 910 maxIndex = thisRange[1]
909 911
910 912 nProfiles += maxIndex - minIndex + 1
911 913
912 914 for thisRange in rangeList:
913 915
914 916 minIndex = thisRange[0]
915 917 maxIndex = thisRange[1]
916 918
917 919 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
918 920
919 921 self.nProfiles = nProfiles
920 922 dataOut.nProfiles = self.nProfiles
921 923 dataOut.profileIndex = self.profileIndex
922 924 dataOut.flagNoData = False
923 925
924 926 self.incProfileIndex()
925 927
926 928 break
927 929
928 930 return True
929 931
930 932
931 933 if beam != None: #beam is only for AMISR data
932 934 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
933 935 dataOut.flagNoData = False
934 936 dataOut.profileIndex = self.profileIndex
935 937
936 938 self.incProfileIndex()
937 939
938 940 return True
939 941
940 942 raise ValueError, "ProfileSelector needs profileList, profileRangeList or rangeList parameter"
941 943
942 944 return False
943 945
944 946 class Reshaper(Operation):
945 947
946 948 def __init__(self, **kwargs):
947 949
948 950 Operation.__init__(self, **kwargs)
949 951
950 952 self.__buffer = None
951 953 self.__nitems = 0
952 954
953 955 def __appendProfile(self, dataOut, nTxs):
954 956
955 957 if self.__buffer is None:
956 958 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
957 959 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
958 960
959 961 ini = dataOut.nHeights * self.__nitems
960 962 end = ini + dataOut.nHeights
961 963
962 964 self.__buffer[:, ini:end] = dataOut.data
963 965
964 966 self.__nitems += 1
965 967
966 968 return int(self.__nitems*nTxs)
967 969
968 970 def __getBuffer(self):
969 971
970 972 if self.__nitems == int(1./self.__nTxs):
971 973
972 974 self.__nitems = 0
973 975
974 976 return self.__buffer.copy()
975 977
976 978 return None
977 979
978 980 def __checkInputs(self, dataOut, shape, nTxs):
979 981
980 982 if shape is None and nTxs is None:
981 983 raise ValueError, "Reshaper: shape of factor should be defined"
982 984
983 985 if nTxs:
984 986 if nTxs < 0:
985 987 raise ValueError, "nTxs should be greater than 0"
986 988
987 989 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
988 990 raise ValueError, "nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs))
989 991
990 992 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
991 993
992 994 return shape, nTxs
993 995
994 996 if len(shape) != 2 and len(shape) != 3:
995 997 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)
996 998
997 999 if len(shape) == 2:
998 1000 shape_tuple = [dataOut.nChannels]
999 1001 shape_tuple.extend(shape)
1000 1002 else:
1001 1003 shape_tuple = list(shape)
1002 1004
1003 1005 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1004 1006
1005 1007 return shape_tuple, nTxs
1006 1008
1007 1009 def run(self, dataOut, shape=None, nTxs=None):
1008 1010
1009 1011 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1010 1012
1011 1013 dataOut.flagNoData = True
1012 1014 profileIndex = None
1013 1015
1014 1016 if dataOut.flagDataAsBlock:
1015 1017
1016 1018 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1017 1019 dataOut.flagNoData = False
1018 1020
1019 1021 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1020 1022
1021 1023 else:
1022 1024
1023 1025 if self.__nTxs < 1:
1024 1026
1025 1027 self.__appendProfile(dataOut, self.__nTxs)
1026 1028 new_data = self.__getBuffer()
1027 1029
1028 1030 if new_data is not None:
1029 1031 dataOut.data = new_data
1030 1032 dataOut.flagNoData = False
1031 1033
1032 1034 profileIndex = dataOut.profileIndex*nTxs
1033 1035
1034 1036 else:
1035 1037 raise ValueError, "nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)"
1036 1038
1037 1039 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1038 1040
1039 1041 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1040 1042
1041 1043 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1042 1044
1043 1045 dataOut.profileIndex = profileIndex
1044 1046
1045 1047 dataOut.ippSeconds /= self.__nTxs
1046 1048
1047 1049 class SplitProfiles(Operation):
1048 1050
1049 1051 def __init__(self, **kwargs):
1050 1052
1051 1053 Operation.__init__(self, **kwargs)
1052 1054
1053 1055 def run(self, dataOut, n):
1054 1056
1055 1057 dataOut.flagNoData = True
1056 1058 profileIndex = None
1057 1059
1058 1060 if dataOut.flagDataAsBlock:
1059 1061
1060 1062 #nchannels, nprofiles, nsamples
1061 1063 shape = dataOut.data.shape
1062 1064
1063 1065 if shape[2] % n != 0:
1064 1066 raise ValueError, "Could not split the data, n=%d has to be multiple of %d" %(n, shape[2])
1065 1067
1066 1068 new_shape = shape[0], shape[1]*n, shape[2]/n
1067 1069
1068 1070 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1069 1071 dataOut.flagNoData = False
1070 1072
1071 1073 profileIndex = int(dataOut.nProfiles/n) - 1
1072 1074
1073 1075 else:
1074 1076
1075 1077 raise ValueError, "Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)"
1076 1078
1077 1079 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1078 1080
1079 1081 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1080 1082
1081 1083 dataOut.nProfiles = int(dataOut.nProfiles*n)
1082 1084
1083 1085 dataOut.profileIndex = profileIndex
1084 1086
1085 1087 dataOut.ippSeconds /= n
1086 1088
1087 1089 class CombineProfiles(Operation):
1088 1090
1089 1091 def __init__(self, **kwargs):
1090 1092
1091 1093 Operation.__init__(self, **kwargs)
1092 1094
1093 1095 self.__remData = None
1094 1096 self.__profileIndex = 0
1095 1097
1096 1098 def run(self, dataOut, n):
1097 1099
1098 1100 dataOut.flagNoData = True
1099 1101 profileIndex = None
1100 1102
1101 1103 if dataOut.flagDataAsBlock:
1102 1104
1103 1105 #nchannels, nprofiles, nsamples
1104 1106 shape = dataOut.data.shape
1105 1107 new_shape = shape[0], shape[1]/n, shape[2]*n
1106 1108
1107 1109 if shape[1] % n != 0:
1108 1110 raise ValueError, "Could not split the data, n=%d has to be multiple of %d" %(n, shape[1])
1109 1111
1110 1112 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1111 1113 dataOut.flagNoData = False
1112 1114
1113 1115 profileIndex = int(dataOut.nProfiles*n) - 1
1114 1116
1115 1117 else:
1116 1118
1117 1119 #nchannels, nsamples
1118 1120 if self.__remData is None:
1119 1121 newData = dataOut.data
1120 1122 else:
1121 1123 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1122 1124
1123 1125 self.__profileIndex += 1
1124 1126
1125 1127 if self.__profileIndex < n:
1126 1128 self.__remData = newData
1127 1129 #continue
1128 1130 return
1129 1131
1130 1132 self.__profileIndex = 0
1131 1133 self.__remData = None
1132 1134
1133 1135 dataOut.data = newData
1134 1136 dataOut.flagNoData = False
1135 1137
1136 1138 profileIndex = dataOut.profileIndex/n
1137 1139
1138 1140
1139 1141 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1140 1142
1141 1143 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1142 1144
1143 1145 dataOut.nProfiles = int(dataOut.nProfiles/n)
1144 1146
1145 1147 dataOut.profileIndex = profileIndex
1146 1148
1147 1149 dataOut.ippSeconds *= n
1148 1150
1149 1151 # import collections
1150 1152 # from scipy.stats import mode
1151 1153 #
1152 1154 # class Synchronize(Operation):
1153 1155 #
1154 1156 # isConfig = False
1155 1157 # __profIndex = 0
1156 1158 #
1157 1159 # def __init__(self, **kwargs):
1158 1160 #
1159 1161 # Operation.__init__(self, **kwargs)
1160 1162 # # self.isConfig = False
1161 1163 # self.__powBuffer = None
1162 1164 # self.__startIndex = 0
1163 1165 # self.__pulseFound = False
1164 1166 #
1165 1167 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1166 1168 #
1167 1169 # #Read data
1168 1170 #
1169 1171 # powerdB = dataOut.getPower(channel = channel)
1170 1172 # noisedB = dataOut.getNoise(channel = channel)[0]
1171 1173 #
1172 1174 # self.__powBuffer.extend(powerdB.flatten())
1173 1175 #
1174 1176 # dataArray = numpy.array(self.__powBuffer)
1175 1177 #
1176 1178 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1177 1179 #
1178 1180 # maxValue = numpy.nanmax(filteredPower)
1179 1181 #
1180 1182 # if maxValue < noisedB + 10:
1181 1183 # #No se encuentra ningun pulso de transmision
1182 1184 # return None
1183 1185 #
1184 1186 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1185 1187 #
1186 1188 # if len(maxValuesIndex) < 2:
1187 1189 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1188 1190 # return None
1189 1191 #
1190 1192 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1191 1193 #
1192 1194 # #Seleccionar solo valores con un espaciamiento de nSamples
1193 1195 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1194 1196 #
1195 1197 # if len(pulseIndex) < 2:
1196 1198 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1197 1199 # return None
1198 1200 #
1199 1201 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1200 1202 #
1201 1203 # #remover senales que se distancien menos de 10 unidades o muestras
1202 1204 # #(No deberian existir IPP menor a 10 unidades)
1203 1205 #
1204 1206 # realIndex = numpy.where(spacing > 10 )[0]
1205 1207 #
1206 1208 # if len(realIndex) < 2:
1207 1209 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1208 1210 # return None
1209 1211 #
1210 1212 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1211 1213 # realPulseIndex = pulseIndex[realIndex]
1212 1214 #
1213 1215 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1214 1216 #
1215 1217 # print "IPP = %d samples" %period
1216 1218 #
1217 1219 # self.__newNSamples = dataOut.nHeights #int(period)
1218 1220 # self.__startIndex = int(realPulseIndex[0])
1219 1221 #
1220 1222 # return 1
1221 1223 #
1222 1224 #
1223 1225 # def setup(self, nSamples, nChannels, buffer_size = 4):
1224 1226 #
1225 1227 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1226 1228 # maxlen = buffer_size*nSamples)
1227 1229 #
1228 1230 # bufferList = []
1229 1231 #
1230 1232 # for i in range(nChannels):
1231 1233 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1232 1234 # maxlen = buffer_size*nSamples)
1233 1235 #
1234 1236 # bufferList.append(bufferByChannel)
1235 1237 #
1236 1238 # self.__nSamples = nSamples
1237 1239 # self.__nChannels = nChannels
1238 1240 # self.__bufferList = bufferList
1239 1241 #
1240 1242 # def run(self, dataOut, channel = 0):
1241 1243 #
1242 1244 # if not self.isConfig:
1243 1245 # nSamples = dataOut.nHeights
1244 1246 # nChannels = dataOut.nChannels
1245 1247 # self.setup(nSamples, nChannels)
1246 1248 # self.isConfig = True
1247 1249 #
1248 1250 # #Append new data to internal buffer
1249 1251 # for thisChannel in range(self.__nChannels):
1250 1252 # bufferByChannel = self.__bufferList[thisChannel]
1251 1253 # bufferByChannel.extend(dataOut.data[thisChannel])
1252 1254 #
1253 1255 # if self.__pulseFound:
1254 1256 # self.__startIndex -= self.__nSamples
1255 1257 #
1256 1258 # #Finding Tx Pulse
1257 1259 # if not self.__pulseFound:
1258 1260 # indexFound = self.__findTxPulse(dataOut, channel)
1259 1261 #
1260 1262 # if indexFound == None:
1261 1263 # dataOut.flagNoData = True
1262 1264 # return
1263 1265 #
1264 1266 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1265 1267 # self.__pulseFound = True
1266 1268 # self.__startIndex = indexFound
1267 1269 #
1268 1270 # #If pulse was found ...
1269 1271 # for thisChannel in range(self.__nChannels):
1270 1272 # bufferByChannel = self.__bufferList[thisChannel]
1271 1273 # #print self.__startIndex
1272 1274 # x = numpy.array(bufferByChannel)
1273 1275 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1274 1276 #
1275 1277 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1276 1278 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1277 1279 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1278 1280 #
1279 1281 # dataOut.data = self.__arrayBuffer
1280 1282 #
1281 1283 # self.__startIndex += self.__newNSamples
1282 1284 #
1283 1285 # return
General Comments 0
You need to be logged in to leave comments. Login now