##// END OF EJS Templates
schain
RSS
fix minHei
Alexander Valdez -
r1676:077ec9b5e257
parent child
Show More
Side by Side Unified
Context file:
r1676:077ec9b5e257 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1629 +1,1627
1 1 import sys
2 2 import numpy,math
3 3 from scipy import interpolate
4 4 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
5 5 from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
6 6 from schainpy.utils import log
7 7 from time import time
8 8 # voltage proc master
9 9
10 10 class VoltageProc(ProcessingUnit):
11 11
12 12 def __init__(self):
13 13
14 14 ProcessingUnit.__init__(self)
15 15
16 16 self.dataOut = Voltage()
17 17 self.flip = 1
18 18 self.setupReq = False
19 19
20 20 def run(self, runNextUnit = 0):
21 21
22 22 if self.dataIn.type == 'AMISR':
23 23 self.__updateObjFromAmisrInput()
24 24
25 25 if self.dataIn.type == 'Voltage':
26 26 self.dataOut.copy(self.dataIn)
27 27 self.dataOut.runNextUnit = runNextUnit
28 28
29 29
30 30 def __updateObjFromAmisrInput(self):
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.flagNoData = self.dataIn.flagNoData
38 38 self.dataOut.data = self.dataIn.data
39 39 self.dataOut.utctime = self.dataIn.utctime
40 40 self.dataOut.channelList = self.dataIn.channelList
41 41 #self.dataOut.timeInterval = self.dataIn.timeInterval
42 42 self.dataOut.heightList = self.dataIn.heightList
43 43 self.dataOut.nProfiles = self.dataIn.nProfiles
44 44
45 45 self.dataOut.nCohInt = self.dataIn.nCohInt
46 46 self.dataOut.ippSeconds = self.dataIn.ippSeconds
47 47 self.dataOut.frequency = self.dataIn.frequency
48 48
49 49 self.dataOut.azimuth = self.dataIn.azimuth
50 50 self.dataOut.zenith = self.dataIn.zenith
51 51
52 52 self.dataOut.beam.codeList = self.dataIn.beam.codeList
53 53 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
54 54 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
55 55
56 56
57 57 class selectChannels(Operation):
58 58
59 59 def run(self, dataOut, channelList):
60 60
61 61 channelIndexList = []
62 62 self.dataOut = dataOut
63 63 for channel in channelList:
64 64 if channel not in self.dataOut.channelList:
65 65 raise ValueError("Channel %d is not in %s" %(channel, str(self.dataOut.channelList)))
66 66
67 67 index = self.dataOut.channelList.index(channel)
68 68 channelIndexList.append(index)
69 69 self.selectChannelsByIndex(channelIndexList)
70 70 return self.dataOut
71 71
72 72 def selectChannelsByIndex(self, channelIndexList):
73 73 """
74 74 Selecciona un bloque de datos en base a canales segun el channelIndexList
75 75
76 76 Input:
77 77 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
78 78
79 79 Affected:
80 80 self.dataOut.data
81 81 self.dataOut.channelIndexList
82 82 self.dataOut.nChannels
83 83 self.dataOut.m_ProcessingHeader.totalSpectra
84 84 self.dataOut.systemHeaderObj.numChannels
85 85 self.dataOut.m_ProcessingHeader.blockSize
86 86
87 87 Return:
88 88 None
89 89 """
90 90
91 91 for channelIndex in channelIndexList:
92 92 if channelIndex not in self.dataOut.channelIndexList:
93 93 raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
94 94
95 95 if self.dataOut.type == 'Voltage':
96 96 if self.dataOut.flagDataAsBlock:
97 97 """
98 98 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
99 99 """
100 100 data = self.dataOut.data[channelIndexList,:,:]
101 101 else:
102 102 data = self.dataOut.data[channelIndexList,:]
103 103
104 104 self.dataOut.data = data
105 105 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
106 106 self.dataOut.channelList = range(len(channelIndexList))
107 107
108 108 elif self.dataOut.type == 'Spectra':
109 109 data_spc = self.dataOut.data_spc[channelIndexList, :]
110 110 data_dc = self.dataOut.data_dc[channelIndexList, :]
111 111
112 112 self.dataOut.data_spc = data_spc
113 113 self.dataOut.data_dc = data_dc
114 114
115 115 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
116 116 self.dataOut.channelList = range(len(channelIndexList))
117 117 self.__selectPairsByChannel(channelIndexList)
118 118
119 119 return 1
120 120
121 121 def __selectPairsByChannel(self, channelList=None):
122 122
123 123 if channelList == None:
124 124 return
125 125
126 126 pairsIndexListSelected = []
127 127 for pairIndex in self.dataOut.pairsIndexList:
128 128 # First pair
129 129 if self.dataOut.pairsList[pairIndex][0] not in channelList:
130 130 continue
131 131 # Second pair
132 132 if self.dataOut.pairsList[pairIndex][1] not in channelList:
133 133 continue
134 134
135 135 pairsIndexListSelected.append(pairIndex)
136 136
137 137 if not pairsIndexListSelected:
138 138 self.dataOut.data_cspc = None
139 139 self.dataOut.pairsList = []
140 140 return
141 141
142 142 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
143 143 self.dataOut.pairsList = [self.dataOut.pairsList[i]
144 144 for i in pairsIndexListSelected]
145 145
146 146 return
147 147
148 148 class selectHeights(Operation):
149 149
150 150 def run(self, dataOut, minHei=None, maxHei=None, minIndex=None, maxIndex=None):
151 151 """
152 152 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
153 153 minHei <= height <= maxHei
154 154
155 155 Input:
156 156 minHei : valor minimo de altura a considerar
157 157 maxHei : valor maximo de altura a considerar
158 158
159 159 Affected:
160 160 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
161 161
162 162 Return:
163 163 1 si el metodo se ejecuto con exito caso contrario devuelve 0
164 164 """
165 165
166 166 self.dataOut = dataOut
167 167
168 168 if type(minHei) == int or type(minHei) == float:
169 169 v_minHei= True
170 170 else:
171 171 v_minHei= False
172 172
173 173 if v_minHei and maxHei:
174 174 if (minHei < self.dataOut.heightList[0]):
175 175 minHei = self.dataOut.heightList[0]
176 176
177 177 if (maxHei > self.dataOut.heightList[-1]):
178 178 maxHei = self.dataOut.heightList[-1]
179 179
180 180 minIndex = 0
181 181 maxIndex = 0
182 182 heights = self.dataOut.heightList
183 183 inda = numpy.where(heights >= minHei)
184 184 indb = numpy.where(heights <= maxHei)
185 185
186 186 try:
187 187 minIndex = inda[0][0]
188 188 except:
189 189 minIndex = 0
190 190
191 191 try:
192 192 maxIndex = indb[0][-1]
193 193 except:
194 194 maxIndex = len(heights)
195 print(minIndex)
196 print(maxIndex)
197 195 self.selectHeightsByIndex(minIndex, maxIndex)
198 196
199 197 return self.dataOut
200 198
201 199 def selectHeightsByIndex(self, minIndex, maxIndex):
202 200 """
203 201 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
204 202 minIndex <= index <= maxIndex
205 203
206 204 Input:
207 205 minIndex : valor de indice minimo de altura a considerar
208 206 maxIndex : valor de indice maximo de altura a considerar
209 207
210 208 Affected:
211 209 self.dataOut.data
212 210 self.dataOut.heightList
213 211
214 212 Return:
215 213 1 si el metodo se ejecuto con exito caso contrario devuelve 0
216 214 """
217 215
218 216 if self.dataOut.type == 'Voltage':
219 217 print(minIndex)
220 218 print(maxIndex)
221 219 if (minIndex < 0) or (minIndex > maxIndex):
222 220 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
223 221
224 222 if (maxIndex >= self.dataOut.nHeights):
225 223 maxIndex = self.dataOut.nHeights
226 224
227 225 #voltage
228 226 if self.dataOut.flagDataAsBlock:
229 227 """
230 228 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
231 229 """
232 230 data = self.dataOut.data[:,:, minIndex:maxIndex]
233 231 else:
234 232 data = self.dataOut.data[:, minIndex:maxIndex]
235 233
236 234 # firstHeight = self.dataOut.heightList[minIndex]
237 235
238 236 self.dataOut.data = data
239 237 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
240 238
241 239 if self.dataOut.nHeights <= 1:
242 240 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
243 241 elif self.dataOut.type == 'Spectra':
244 242 if (minIndex < 0) or (minIndex > maxIndex):
245 243 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
246 244 minIndex, maxIndex))
247 245
248 246 if (maxIndex >= self.dataOut.nHeights):
249 247 maxIndex = self.dataOut.nHeights - 1
250 248
251 249 # Spectra
252 250 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
253 251
254 252 data_cspc = None
255 253 if self.dataOut.data_cspc is not None:
256 254 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
257 255
258 256 data_dc = None
259 257 if self.dataOut.data_dc is not None:
260 258 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
261 259
262 260 self.dataOut.data_spc = data_spc
263 261 self.dataOut.data_cspc = data_cspc
264 262 self.dataOut.data_dc = data_dc
265 263
266 264 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
267 265
268 266 return 1
269 267
270 268
271 269 class filterByHeights(Operation):
272 270
273 271 def run(self, dataOut, window):
274 272
275 273 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
276 274
277 275 if window == None:
278 276 window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
279 277
280 278 newdelta = deltaHeight * window
281 279 r = dataOut.nHeights % window
282 280 newheights = (dataOut.nHeights-r)/window
283 281
284 282 if newheights <= 1:
285 283 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
286 284
287 285 if dataOut.flagDataAsBlock:
288 286 """
289 287 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
290 288 """
291 289 buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)]
292 290 buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window)
293 291 buffer = numpy.sum(buffer,3)
294 292
295 293 else:
296 294 buffer = dataOut.data[:,0:int(dataOut.nHeights-r)]
297 295 buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
298 296 buffer = numpy.sum(buffer,2)
299 297
300 298 dataOut.data = buffer
301 299 dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta
302 300 dataOut.windowOfFilter = window
303 301
304 302 return dataOut
305 303
306 304
307 305 class setH0(Operation):
308 306
309 307 def run(self, dataOut, h0, deltaHeight = None):
310 308
311 309 if not deltaHeight:
312 310 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
313 311
314 312 nHeights = dataOut.nHeights
315 313
316 314 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
317 315
318 316 dataOut.heightList = newHeiRange
319 317
320 318 return dataOut
321 319
322 320
323 321 class deFlip(Operation):
324 322
325 323 def run(self, dataOut, channelList = []):
326 324
327 325 data = dataOut.data.copy()
328 326
329 327 if dataOut.flagDataAsBlock:
330 328 flip = self.flip
331 329 profileList = list(range(dataOut.nProfiles))
332 330
333 331 if not channelList:
334 332 for thisProfile in profileList:
335 333 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
336 334 flip *= -1.0
337 335 else:
338 336 for thisChannel in channelList:
339 337 if thisChannel not in dataOut.channelList:
340 338 continue
341 339
342 340 for thisProfile in profileList:
343 341 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
344 342 flip *= -1.0
345 343
346 344 self.flip = flip
347 345
348 346 else:
349 347 if not channelList:
350 348 data[:,:] = data[:,:]*self.flip
351 349 else:
352 350 for thisChannel in channelList:
353 351 if thisChannel not in dataOut.channelList:
354 352 continue
355 353
356 354 data[thisChannel,:] = data[thisChannel,:]*self.flip
357 355
358 356 self.flip *= -1.
359 357
360 358 dataOut.data = data
361 359
362 360 return dataOut
363 361
364 362
365 363 class setAttribute(Operation):
366 364 '''
367 365 Set an arbitrary attribute(s) to dataOut
368 366 '''
369 367
370 368 def __init__(self):
371 369
372 370 Operation.__init__(self)
373 371 self._ready = False
374 372
375 373 def run(self, dataOut, **kwargs):
376 374
377 375 for key, value in kwargs.items():
378 376 setattr(dataOut, key, value)
379 377
380 378 return dataOut
381 379
382 380
383 381 @MPDecorator
384 382 class printAttribute(Operation):
385 383 '''
386 384 Print an arbitrary attribute of dataOut
387 385 '''
388 386
389 387 def __init__(self):
390 388
391 389 Operation.__init__(self)
392 390
393 391 def run(self, dataOut, attributes):
394 392
395 393 if isinstance(attributes, str):
396 394 attributes = [attributes]
397 395 for attr in attributes:
398 396 if hasattr(dataOut, attr):
399 397 log.log(getattr(dataOut, attr), attr)
400 398
401 399
402 400 class interpolateHeights(Operation):
403 401
404 402 def run(self, dataOut, topLim, botLim):
405 403 #69 al 72 para julia
406 404 #82-84 para meteoros
407 405 if len(numpy.shape(dataOut.data))==2:
408 406 sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
409 407 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
410 408 #dataOut.data[:,botLim:limSup+1] = sampInterp
411 409 dataOut.data[:,botLim:topLim+1] = sampInterp
412 410 else:
413 411 nHeights = dataOut.data.shape[2]
414 412 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
415 413 y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
416 414 f = interpolate.interp1d(x, y, axis = 2)
417 415 xnew = numpy.arange(botLim,topLim+1)
418 416 ynew = f(xnew)
419 417 dataOut.data[:,:,botLim:topLim+1] = ynew
420 418
421 419 return dataOut
422 420
423 421
424 422 class CohInt(Operation):
425 423
426 424 isConfig = False
427 425 __profIndex = 0
428 426 __byTime = False
429 427 __initime = None
430 428 __lastdatatime = None
431 429 __integrationtime = None
432 430 __buffer = None
433 431 __bufferStride = []
434 432 __dataReady = False
435 433 __profIndexStride = 0
436 434 __dataToPutStride = False
437 435 n = None
438 436
439 437 def __init__(self, **kwargs):
440 438
441 439 Operation.__init__(self, **kwargs)
442 440
443 441 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
444 442 """
445 443 Set the parameters of the integration class.
446 444
447 445 Inputs:
448 446
449 447 n : Number of coherent integrations
450 448 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
451 449 overlapping :
452 450 """
453 451
454 452 self.__initime = None
455 453 self.__lastdatatime = 0
456 454 self.__buffer = None
457 455 self.__dataReady = False
458 456 self.byblock = byblock
459 457 self.stride = stride
460 458
461 459 if n == None and timeInterval == None:
462 460 raise ValueError("n or timeInterval should be specified ...")
463 461
464 462 if n != None:
465 463 self.n = n
466 464 self.__byTime = False
467 465 else:
468 466 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
469 467 self.n = 9999
470 468 self.__byTime = True
471 469
472 470 if overlapping:
473 471 self.__withOverlapping = True
474 472 self.__buffer = None
475 473 else:
476 474 self.__withOverlapping = False
477 475 self.__buffer = 0
478 476
479 477 self.__profIndex = 0
480 478
481 479 def putData(self, data):
482 480
483 481 """
484 482 Add a profile to the __buffer and increase in one the __profileIndex
485 483
486 484 """
487 485
488 486 if not self.__withOverlapping:
489 487 self.__buffer += data.copy()
490 488 self.__profIndex += 1
491 489 return
492 490
493 491 #Overlapping data
494 492 nChannels, nHeis = data.shape
495 493 data = numpy.reshape(data, (1, nChannels, nHeis))
496 494
497 495 #If the buffer is empty then it takes the data value
498 496 if self.__buffer is None:
499 497 self.__buffer = data
500 498 self.__profIndex += 1
501 499 return
502 500
503 501 #If the buffer length is lower than n then stakcing the data value
504 502 if self.__profIndex < self.n:
505 503 self.__buffer = numpy.vstack((self.__buffer, data))
506 504 self.__profIndex += 1
507 505 return
508 506
509 507 #If the buffer length is equal to n then replacing the last buffer value with the data value
510 508 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
511 509 self.__buffer[self.n-1] = data
512 510 self.__profIndex = self.n
513 511 return
514 512
515 513
516 514 def pushData(self):
517 515 """
518 516 Return the sum of the last profiles and the profiles used in the sum.
519 517
520 518 Affected:
521 519
522 520 self.__profileIndex
523 521
524 522 """
525 523
526 524 if not self.__withOverlapping:
527 525 data = self.__buffer
528 526 n = self.__profIndex
529 527
530 528 self.__buffer = 0
531 529 self.__profIndex = 0
532 530
533 531 return data, n
534 532
535 533 #Integration with Overlapping
536 534 data = numpy.sum(self.__buffer, axis=0)
537 535 # print data
538 536 # raise
539 537 n = self.__profIndex
540 538
541 539 return data, n
542 540
543 541 def byProfiles(self, data):
544 542
545 543 self.__dataReady = False
546 544 avgdata = None
547 545 # n = None
548 546 # print data
549 547 # raise
550 548 self.putData(data)
551 549
552 550 if self.__profIndex == self.n:
553 551 avgdata, n = self.pushData()
554 552 self.__dataReady = True
555 553
556 554 return avgdata
557 555
558 556 def byTime(self, data, datatime):
559 557
560 558 self.__dataReady = False
561 559 avgdata = None
562 560 n = None
563 561
564 562 self.putData(data)
565 563
566 564 if (datatime - self.__initime) >= self.__integrationtime:
567 565 avgdata, n = self.pushData()
568 566 self.n = n
569 567 self.__dataReady = True
570 568
571 569 return avgdata
572 570
573 571 def integrateByStride(self, data, datatime):
574 572 # print data
575 573 if self.__profIndex == 0:
576 574 self.__buffer = [[data.copy(), datatime]]
577 575 else:
578 576 self.__buffer.append([data.copy(),datatime])
579 577 self.__profIndex += 1
580 578 self.__dataReady = False
581 579
582 580 if self.__profIndex == self.n * self.stride :
583 581 self.__dataToPutStride = True
584 582 self.__profIndexStride = 0
585 583 self.__profIndex = 0
586 584 self.__bufferStride = []
587 585 for i in range(self.stride):
588 586 current = self.__buffer[i::self.stride]
589 587 data = numpy.sum([t[0] for t in current], axis=0)
590 588 avgdatatime = numpy.average([t[1] for t in current])
591 589 # print data
592 590 self.__bufferStride.append((data, avgdatatime))
593 591
594 592 if self.__dataToPutStride:
595 593 self.__dataReady = True
596 594 self.__profIndexStride += 1
597 595 if self.__profIndexStride == self.stride:
598 596 self.__dataToPutStride = False
599 597 # print self.__bufferStride[self.__profIndexStride - 1]
600 598 # raise
601 599 return self.__bufferStride[self.__profIndexStride - 1]
602 600
603 601
604 602 return None, None
605 603
606 604 def integrate(self, data, datatime=None):
607 605
608 606 if self.__initime == None:
609 607 self.__initime = datatime
610 608
611 609 if self.__byTime:
612 610 avgdata = self.byTime(data, datatime)
613 611 else:
614 612 avgdata = self.byProfiles(data)
615 613
616 614
617 615 self.__lastdatatime = datatime
618 616
619 617 if avgdata is None:
620 618 return None, None
621 619
622 620 avgdatatime = self.__initime
623 621
624 622 deltatime = datatime - self.__lastdatatime
625 623
626 624 if not self.__withOverlapping:
627 625 self.__initime = datatime
628 626 else:
629 627 self.__initime += deltatime
630 628
631 629 return avgdata, avgdatatime
632 630
633 631 def integrateByBlock(self, dataOut):
634 632
635 633 times = int(dataOut.data.shape[1]/self.n)
636 634 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
637 635
638 636 id_min = 0
639 637 id_max = self.n
640 638
641 639 for i in range(times):
642 640 junk = dataOut.data[:,id_min:id_max,:]
643 641 avgdata[:,i,:] = junk.sum(axis=1)
644 642 id_min += self.n
645 643 id_max += self.n
646 644
647 645 timeInterval = dataOut.ippSeconds*self.n
648 646 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
649 647 self.__dataReady = True
650 648 return avgdata, avgdatatime
651 649
652 650 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
653 651
654 652 if not self.isConfig:
655 653 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
656 654 self.isConfig = True
657 655 if dataOut.flagDataAsBlock:
658 656 """
659 657 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
660 658 """
661 659 avgdata, avgdatatime = self.integrateByBlock(dataOut)
662 660 dataOut.nProfiles /= self.n
663 661 else:
664 662 if stride is None:
665 663 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
666 664 else:
667 665 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
668 666
669 667
670 668 # dataOut.timeInterval *= n
671 669 dataOut.flagNoData = True
672 670
673 671 if self.__dataReady:
674 672 dataOut.data = avgdata
675 673 if not dataOut.flagCohInt:
676 674 dataOut.nCohInt *= self.n
677 675 dataOut.flagCohInt = True
678 676 dataOut.utctime = avgdatatime
679 677 # print avgdata, avgdatatime
680 678 # raise
681 679 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
682 680 dataOut.flagNoData = False
683 681 return dataOut
684 682
685 683 class Decoder(Operation):
686 684
687 685 isConfig = False
688 686 __profIndex = 0
689 687
690 688 code = None
691 689
692 690 nCode = None
693 691 nBaud = None
694 692
695 693 def __init__(self, **kwargs):
696 694
697 695 Operation.__init__(self, **kwargs)
698 696
699 697 self.times = None
700 698 self.osamp = None
701 699 # self.__setValues = False
702 700 self.isConfig = False
703 701 self.setupReq = False
704 702 def setup(self, code, osamp, dataOut):
705 703
706 704 self.__profIndex = 0
707 705
708 706 self.code = code
709 707
710 708 self.nCode = len(code)
711 709 self.nBaud = len(code[0])
712 710
713 711 if (osamp != None) and (osamp >1):
714 712 self.osamp = osamp
715 713 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
716 714 self.nBaud = self.nBaud*self.osamp
717 715
718 716 self.__nChannels = dataOut.nChannels
719 717 self.__nProfiles = dataOut.nProfiles
720 718 self.__nHeis = dataOut.nHeights
721 719
722 720 if self.__nHeis < self.nBaud:
723 721 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
724 722
725 723 #Frequency
726 724 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
727 725
728 726 __codeBuffer[:,0:self.nBaud] = self.code
729 727
730 728 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
731 729
732 730 if dataOut.flagDataAsBlock:
733 731
734 732 self.ndatadec = self.__nHeis #- self.nBaud + 1
735 733
736 734 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
737 735
738 736 else:
739 737
740 738 #Time
741 739 self.ndatadec = self.__nHeis #- self.nBaud + 1
742 740
743 741 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
744 742
745 743 def __convolutionInFreq(self, data):
746 744
747 745 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
748 746
749 747 fft_data = numpy.fft.fft(data, axis=1)
750 748
751 749 conv = fft_data*fft_code
752 750
753 751 data = numpy.fft.ifft(conv,axis=1)
754 752
755 753 return data
756 754
757 755 def __convolutionInFreqOpt(self, data):
758 756
759 757 raise NotImplementedError
760 758
761 759 def __convolutionInTime(self, data):
762 760
763 761 code = self.code[self.__profIndex]
764 762 for i in range(self.__nChannels):
765 763 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
766 764
767 765 return self.datadecTime
768 766
769 767 def __convolutionByBlockInTime(self, data):
770 768
771 769 repetitions = int(self.__nProfiles / self.nCode)
772 770 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
773 771 junk = junk.flatten()
774 772 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
775 773 profilesList = range(self.__nProfiles)
776 774
777 775 for i in range(self.__nChannels):
778 776 for j in profilesList:
779 777 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
780 778 return self.datadecTime
781 779
782 780 def __convolutionByBlockInFreq(self, data):
783 781
784 782 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
785 783
786 784
787 785 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
788 786
789 787 fft_data = numpy.fft.fft(data, axis=2)
790 788
791 789 conv = fft_data*fft_code
792 790
793 791 data = numpy.fft.ifft(conv,axis=2)
794 792
795 793 return data
796 794
797 795
798 796 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
799 797
800 798 if dataOut.flagDecodeData:
801 799 print("This data is already decoded, recoding again ...")
802 800
803 801 if not self.isConfig:
804 802
805 803 if code is None:
806 804 if dataOut.code is None:
807 805 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
808 806
809 807 code = dataOut.code
810 808 else:
811 809 code = numpy.array(code).reshape(nCode,nBaud)
812 810 self.setup(code, osamp, dataOut)
813 811
814 812 self.isConfig = True
815 813
816 814 if mode == 3:
817 815 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
818 816
819 817 if times != None:
820 818 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
821 819
822 820 if self.code is None:
823 821 print("Fail decoding: Code is not defined.")
824 822 return
825 823
826 824 self.__nProfiles = dataOut.nProfiles
827 825 datadec = None
828 826
829 827 if mode == 3:
830 828 mode = 0
831 829
832 830 if dataOut.flagDataAsBlock:
833 831 """
834 832 Decoding when data have been read as block,
835 833 """
836 834
837 835 if mode == 0:
838 836 datadec = self.__convolutionByBlockInTime(dataOut.data)
839 837 if mode == 1:
840 838 datadec = self.__convolutionByBlockInFreq(dataOut.data)
841 839 else:
842 840 """
843 841 Decoding when data have been read profile by profile
844 842 """
845 843 if mode == 0:
846 844 datadec = self.__convolutionInTime(dataOut.data)
847 845
848 846 if mode == 1:
849 847 datadec = self.__convolutionInFreq(dataOut.data)
850 848
851 849 if mode == 2:
852 850 datadec = self.__convolutionInFreqOpt(dataOut.data)
853 851
854 852 if datadec is None:
855 853 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
856 854
857 855 dataOut.code = self.code
858 856 dataOut.nCode = self.nCode
859 857 dataOut.nBaud = self.nBaud
860 858
861 859 dataOut.data = datadec
862 860 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
863 861
864 862 dataOut.flagDecodeData = True #asumo q la data esta decodificada
865 863
866 864 if self.__profIndex == self.nCode-1:
867 865 self.__profIndex = 0
868 866 return dataOut
869 867
870 868 self.__profIndex += 1
871 869
872 870 return dataOut
873 871
874 872
875 873 class ProfileConcat(Operation):
876 874
877 875 isConfig = False
878 876 buffer = None
879 877
880 878 def __init__(self, **kwargs):
881 879
882 880 Operation.__init__(self, **kwargs)
883 881 self.profileIndex = 0
884 882
885 883 def reset(self):
886 884 self.buffer = numpy.zeros_like(self.buffer)
887 885 self.start_index = 0
888 886 self.times = 1
889 887
890 888 def setup(self, data, m, n=1):
891 889 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
892 890 self.nHeights = data.shape[1]#.nHeights
893 891 self.start_index = 0
894 892 self.times = 1
895 893
896 894 def concat(self, data):
897 895
898 896 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
899 897 self.start_index = self.start_index + self.nHeights
900 898
901 899 def run(self, dataOut, m):
902 900 dataOut.flagNoData = True
903 901
904 902 if not self.isConfig:
905 903 self.setup(dataOut.data, m, 1)
906 904 self.isConfig = True
907 905
908 906 if dataOut.flagDataAsBlock:
909 907 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
910 908
911 909 else:
912 910 self.concat(dataOut.data)
913 911 self.times += 1
914 912 if self.times > m:
915 913 dataOut.data = self.buffer
916 914 self.reset()
917 915 dataOut.flagNoData = False
918 916 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
919 917 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
920 918 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
921 919 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
922 920 dataOut.ippSeconds *= m
923 921 return dataOut
924 922
925 923 class ProfileSelector(Operation):
926 924
927 925 profileIndex = None
928 926 # Tamanho total de los perfiles
929 927 nProfiles = None
930 928
931 929 def __init__(self, **kwargs):
932 930
933 931 Operation.__init__(self, **kwargs)
934 932 self.profileIndex = 0
935 933
936 934 def incProfileIndex(self):
937 935
938 936 self.profileIndex += 1
939 937
940 938 if self.profileIndex >= self.nProfiles:
941 939 self.profileIndex = 0
942 940
943 941 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
944 942
945 943 if profileIndex < minIndex:
946 944 return False
947 945
948 946 if profileIndex > maxIndex:
949 947 return False
950 948
951 949 return True
952 950
953 951 def isThisProfileInList(self, profileIndex, profileList):
954 952
955 953 if profileIndex not in profileList:
956 954 return False
957 955
958 956 return True
959 957
960 958 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
961 959
962 960 """
963 961 ProfileSelector:
964 962
965 963 Inputs:
966 964 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
967 965
968 966 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
969 967
970 968 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
971 969
972 970 """
973 971
974 972 if rangeList is not None:
975 973 if type(rangeList[0]) not in (tuple, list):
976 974 rangeList = [rangeList]
977 975
978 976 dataOut.flagNoData = True
979 977
980 978 if dataOut.flagDataAsBlock:
981 979 """
982 980 data dimension = [nChannels, nProfiles, nHeis]
983 981 """
984 982 if profileList != None:
985 983 dataOut.data = dataOut.data[:,profileList,:]
986 984
987 985 if profileRangeList != None:
988 986 minIndex = profileRangeList[0]
989 987 maxIndex = profileRangeList[1]
990 988 profileList = list(range(minIndex, maxIndex+1))
991 989
992 990 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
993 991
994 992 if rangeList != None:
995 993
996 994 profileList = []
997 995
998 996 for thisRange in rangeList:
999 997 minIndex = thisRange[0]
1000 998 maxIndex = thisRange[1]
1001 999
1002 1000 profileList.extend(list(range(minIndex, maxIndex+1)))
1003 1001
1004 1002 dataOut.data = dataOut.data[:,profileList,:]
1005 1003
1006 1004 dataOut.nProfiles = len(profileList)
1007 1005 dataOut.profileIndex = dataOut.nProfiles - 1
1008 1006 dataOut.flagNoData = False
1009 1007
1010 1008 return dataOut
1011 1009
1012 1010 """
1013 1011 data dimension = [nChannels, nHeis]
1014 1012 """
1015 1013
1016 1014 if profileList != None:
1017 1015
1018 1016 if self.isThisProfileInList(dataOut.profileIndex, profileList):
1019 1017
1020 1018 self.nProfiles = len(profileList)
1021 1019 dataOut.nProfiles = self.nProfiles
1022 1020 dataOut.profileIndex = self.profileIndex
1023 1021 dataOut.flagNoData = False
1024 1022
1025 1023 self.incProfileIndex()
1026 1024 return dataOut
1027 1025
1028 1026 if profileRangeList != None:
1029 1027
1030 1028 minIndex = profileRangeList[0]
1031 1029 maxIndex = profileRangeList[1]
1032 1030
1033 1031 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1034 1032
1035 1033 self.nProfiles = maxIndex - minIndex + 1
1036 1034 dataOut.nProfiles = self.nProfiles
1037 1035 dataOut.profileIndex = self.profileIndex
1038 1036 dataOut.flagNoData = False
1039 1037
1040 1038 self.incProfileIndex()
1041 1039 return dataOut
1042 1040
1043 1041 if rangeList != None:
1044 1042
1045 1043 nProfiles = 0
1046 1044
1047 1045 for thisRange in rangeList:
1048 1046 minIndex = thisRange[0]
1049 1047 maxIndex = thisRange[1]
1050 1048
1051 1049 nProfiles += maxIndex - minIndex + 1
1052 1050
1053 1051 for thisRange in rangeList:
1054 1052
1055 1053 minIndex = thisRange[0]
1056 1054 maxIndex = thisRange[1]
1057 1055
1058 1056 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1059 1057
1060 1058 self.nProfiles = nProfiles
1061 1059 dataOut.nProfiles = self.nProfiles
1062 1060 dataOut.profileIndex = self.profileIndex
1063 1061 dataOut.flagNoData = False
1064 1062
1065 1063 self.incProfileIndex()
1066 1064
1067 1065 break
1068 1066
1069 1067 return dataOut
1070 1068
1071 1069
1072 1070 if beam != None: #beam is only for AMISR data
1073 1071 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
1074 1072 dataOut.flagNoData = False
1075 1073 dataOut.profileIndex = self.profileIndex
1076 1074
1077 1075 self.incProfileIndex()
1078 1076
1079 1077 return dataOut
1080 1078
1081 1079 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
1082 1080
1083 1081
1084 1082 class Reshaper(Operation):
1085 1083
1086 1084 def __init__(self, **kwargs):
1087 1085
1088 1086 Operation.__init__(self, **kwargs)
1089 1087
1090 1088 self.__buffer = None
1091 1089 self.__nitems = 0
1092 1090
1093 1091 def __appendProfile(self, dataOut, nTxs):
1094 1092
1095 1093 if self.__buffer is None:
1096 1094 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1097 1095 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1098 1096
1099 1097 ini = dataOut.nHeights * self.__nitems
1100 1098 end = ini + dataOut.nHeights
1101 1099
1102 1100 self.__buffer[:, ini:end] = dataOut.data
1103 1101
1104 1102 self.__nitems += 1
1105 1103
1106 1104 return int(self.__nitems*nTxs)
1107 1105
1108 1106 def __getBuffer(self):
1109 1107
1110 1108 if self.__nitems == int(1./self.__nTxs):
1111 1109
1112 1110 self.__nitems = 0
1113 1111
1114 1112 return self.__buffer.copy()
1115 1113
1116 1114 return None
1117 1115
1118 1116 def __checkInputs(self, dataOut, shape, nTxs):
1119 1117
1120 1118 if shape is None and nTxs is None:
1121 1119 raise ValueError("Reshaper: shape of factor should be defined")
1122 1120
1123 1121 if nTxs:
1124 1122 if nTxs < 0:
1125 1123 raise ValueError("nTxs should be greater than 0")
1126 1124
1127 1125 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1128 1126 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1129 1127
1130 1128 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1131 1129
1132 1130 return shape, nTxs
1133 1131
1134 1132 if len(shape) != 2 and len(shape) != 3:
1135 1133 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))
1136 1134
1137 1135 if len(shape) == 2:
1138 1136 shape_tuple = [dataOut.nChannels]
1139 1137 shape_tuple.extend(shape)
1140 1138 else:
1141 1139 shape_tuple = list(shape)
1142 1140
1143 1141 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1144 1142
1145 1143 return shape_tuple, nTxs
1146 1144
1147 1145 def run(self, dataOut, shape=None, nTxs=None):
1148 1146
1149 1147 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1150 1148
1151 1149 dataOut.flagNoData = True
1152 1150 profileIndex = None
1153 1151
1154 1152 if dataOut.flagDataAsBlock:
1155 1153
1156 1154 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1157 1155 dataOut.flagNoData = False
1158 1156
1159 1157 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1160 1158
1161 1159 else:
1162 1160
1163 1161 if self.__nTxs < 1:
1164 1162
1165 1163 self.__appendProfile(dataOut, self.__nTxs)
1166 1164 new_data = self.__getBuffer()
1167 1165
1168 1166 if new_data is not None:
1169 1167 dataOut.data = new_data
1170 1168 dataOut.flagNoData = False
1171 1169
1172 1170 profileIndex = dataOut.profileIndex*nTxs
1173 1171
1174 1172 else:
1175 1173 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1176 1174
1177 1175 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1178 1176
1179 1177 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1180 1178
1181 1179 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1182 1180
1183 1181 dataOut.profileIndex = profileIndex
1184 1182
1185 1183 dataOut.ippSeconds /= self.__nTxs
1186 1184
1187 1185 return dataOut
1188 1186
1189 1187 class SplitProfiles(Operation):
1190 1188
1191 1189 def __init__(self, **kwargs):
1192 1190
1193 1191 Operation.__init__(self, **kwargs)
1194 1192
1195 1193 def run(self, dataOut, n):
1196 1194
1197 1195 dataOut.flagNoData = True
1198 1196 profileIndex = None
1199 1197
1200 1198 if dataOut.flagDataAsBlock:
1201 1199
1202 1200 #nchannels, nprofiles, nsamples
1203 1201 shape = dataOut.data.shape
1204 1202
1205 1203 if shape[2] % n != 0:
1206 1204 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1207 1205
1208 1206 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1209 1207
1210 1208 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1211 1209 dataOut.flagNoData = False
1212 1210
1213 1211 profileIndex = int(dataOut.nProfiles/n) - 1
1214 1212
1215 1213 else:
1216 1214
1217 1215 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1218 1216
1219 1217 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1220 1218
1221 1219 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1222 1220
1223 1221 dataOut.nProfiles = int(dataOut.nProfiles*n)
1224 1222
1225 1223 dataOut.profileIndex = profileIndex
1226 1224
1227 1225 dataOut.ippSeconds /= n
1228 1226
1229 1227 return dataOut
1230 1228
1231 1229 class CombineProfiles(Operation):
1232 1230 def __init__(self, **kwargs):
1233 1231
1234 1232 Operation.__init__(self, **kwargs)
1235 1233
1236 1234 self.__remData = None
1237 1235 self.__profileIndex = 0
1238 1236
1239 1237 def run(self, dataOut, n):
1240 1238
1241 1239 dataOut.flagNoData = True
1242 1240 profileIndex = None
1243 1241
1244 1242 if dataOut.flagDataAsBlock:
1245 1243
1246 1244 #nchannels, nprofiles, nsamples
1247 1245 shape = dataOut.data.shape
1248 1246 new_shape = shape[0], shape[1]/n, shape[2]*n
1249 1247
1250 1248 if shape[1] % n != 0:
1251 1249 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1252 1250
1253 1251 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1254 1252 dataOut.flagNoData = False
1255 1253
1256 1254 profileIndex = int(dataOut.nProfiles*n) - 1
1257 1255
1258 1256 else:
1259 1257
1260 1258 #nchannels, nsamples
1261 1259 if self.__remData is None:
1262 1260 newData = dataOut.data
1263 1261 else:
1264 1262 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1265 1263
1266 1264 self.__profileIndex += 1
1267 1265
1268 1266 if self.__profileIndex < n:
1269 1267 self.__remData = newData
1270 1268 #continue
1271 1269 return
1272 1270
1273 1271 self.__profileIndex = 0
1274 1272 self.__remData = None
1275 1273
1276 1274 dataOut.data = newData
1277 1275 dataOut.flagNoData = False
1278 1276
1279 1277 profileIndex = dataOut.profileIndex/n
1280 1278
1281 1279
1282 1280 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1283 1281
1284 1282 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1285 1283
1286 1284 dataOut.nProfiles = int(dataOut.nProfiles/n)
1287 1285
1288 1286 dataOut.profileIndex = profileIndex
1289 1287
1290 1288 dataOut.ippSeconds *= n
1291 1289
1292 1290 return dataOut
1293 1291
1294 1292 class PulsePairVoltage(Operation):
1295 1293 '''
1296 1294 Function PulsePair(Signal Power, Velocity)
1297 1295 The real component of Lag[0] provides Intensity Information
1298 1296 The imag component of Lag[1] Phase provides Velocity Information
1299 1297
1300 1298 Configuration Parameters:
1301 1299 nPRF = Number of Several PRF
1302 1300 theta = Degree Azimuth angel Boundaries
1303 1301
1304 1302 Input:
1305 1303 self.dataOut
1306 1304 lag[N]
1307 1305 Affected:
1308 1306 self.dataOut.spc
1309 1307 '''
1310 1308 isConfig = False
1311 1309 __profIndex = 0
1312 1310 __initime = None
1313 1311 __lastdatatime = None
1314 1312 __buffer = None
1315 1313 noise = None
1316 1314 __dataReady = False
1317 1315 n = None
1318 1316 __nch = 0
1319 1317 __nHeis = 0
1320 1318 removeDC = False
1321 1319 ipp = None
1322 1320 lambda_ = 0
1323 1321
1324 1322 def __init__(self,**kwargs):
1325 1323 Operation.__init__(self,**kwargs)
1326 1324
1327 1325 def setup(self, dataOut, n = None, removeDC=False):
1328 1326 '''
1329 1327 n= Numero de PRF's de entrada
1330 1328 '''
1331 1329 self.__initime = None
1332 1330 self.__lastdatatime = 0
1333 1331 self.__dataReady = False
1334 1332 self.__buffer = 0
1335 1333 self.__profIndex = 0
1336 1334 self.noise = None
1337 1335 self.__nch = dataOut.nChannels
1338 1336 self.__nHeis = dataOut.nHeights
1339 1337 self.removeDC = removeDC
1340 1338 self.lambda_ = 3.0e8/(9345.0e6)
1341 1339 self.ippSec = dataOut.ippSeconds
1342 1340 self.nCohInt = dataOut.nCohInt
1343 1341 print("IPPseconds",dataOut.ippSeconds)
1344 1342
1345 1343 print("ELVALOR DE n es:", n)
1346 1344 if n == None:
1347 1345 raise ValueError("n should be specified.")
1348 1346
1349 1347 if n != None:
1350 1348 if n<2:
1351 1349 raise ValueError("n should be greater than 2")
1352 1350
1353 1351 self.n = n
1354 1352 self.__nProf = n
1355 1353
1356 1354 self.__buffer = numpy.zeros((dataOut.nChannels,
1357 1355 n,
1358 1356 dataOut.nHeights),
1359 1357 dtype='complex')
1360 1358
1361 1359 def putData(self,data):
1362 1360 '''
1363 1361 Add a profile to he __buffer and increase in one the __profiel Index
1364 1362 '''
1365 1363 self.__buffer[:,self.__profIndex,:]= data
1366 1364 self.__profIndex += 1
1367 1365 return
1368 1366
1369 1367 def pushData(self,dataOut):
1370 1368 '''
1371 1369 Return the PULSEPAIR and the profiles used in the operation
1372 1370 Affected : self.__profileIndex
1373 1371 '''
1374 1372 #----------------- Remove DC-----------------------------------
1375 1373 if self.removeDC==True:
1376 1374 mean = numpy.mean(self.__buffer,1)
1377 1375 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1378 1376 dc= numpy.tile(tmp,[1,self.__nProf,1])
1379 1377 self.__buffer = self.__buffer - dc
1380 1378 #------------------Calculo de Potencia ------------------------
1381 1379 pair0 = self.__buffer*numpy.conj(self.__buffer)
1382 1380 pair0 = pair0.real
1383 1381 lag_0 = numpy.sum(pair0,1)
1384 1382 #------------------Calculo de Ruido x canal--------------------
1385 1383 self.noise = numpy.zeros(self.__nch)
1386 1384 for i in range(self.__nch):
1387 1385 daux = numpy.sort(pair0[i,:,:],axis= None)
1388 1386 self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt)
1389 1387
1390 1388 self.noise = self.noise.reshape(self.__nch,1)
1391 1389 self.noise = numpy.tile(self.noise,[1,self.__nHeis])
1392 1390 noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis)
1393 1391 noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1])
1394 1392 #------------------ Potencia recibida= P , Potencia senal = S , Ruido= N--
1395 1393 #------------------ P= S+N ,P=lag_0/N ---------------------------------
1396 1394 #-------------------- Power --------------------------------------------------
1397 1395 data_power = lag_0/(self.n*self.nCohInt)
1398 1396 #------------------ Senal ---------------------------------------------------
1399 1397 data_intensity = pair0 - noise_buffer
1400 1398 data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt)
1401 1399 #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt)
1402 1400 for i in range(self.__nch):
1403 1401 for j in range(self.__nHeis):
1404 1402 if data_intensity[i][j] < 0:
1405 1403 data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j]))
1406 1404
1407 1405 #----------------- Calculo de Frecuencia y Velocidad doppler--------
1408 1406 pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
1409 1407 lag_1 = numpy.sum(pair1,1)
1410 1408 data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
1411 1409 data_velocity = (self.lambda_/2.0)*data_freq
1412 1410
1413 1411 #---------------- Potencia promedio estimada de la Senal-----------
1414 1412 lag_0 = lag_0/self.n
1415 1413 S = lag_0-self.noise
1416 1414
1417 1415 #---------------- Frecuencia Doppler promedio ---------------------
1418 1416 lag_1 = lag_1/(self.n-1)
1419 1417 R1 = numpy.abs(lag_1)
1420 1418
1421 1419 #---------------- Calculo del SNR----------------------------------
1422 1420 data_snrPP = S/self.noise
1423 1421 for i in range(self.__nch):
1424 1422 for j in range(self.__nHeis):
1425 1423 if data_snrPP[i][j] < 1.e-20:
1426 1424 data_snrPP[i][j] = 1.e-20
1427 1425
1428 1426 #----------------- Calculo del ancho espectral ----------------------
1429 1427 L = S/R1
1430 1428 L = numpy.where(L<0,1,L)
1431 1429 L = numpy.log(L)
1432 1430 tmp = numpy.sqrt(numpy.absolute(L))
1433 1431 data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L)
1434 1432 n = self.__profIndex
1435 1433
1436 1434 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1437 1435 self.__profIndex = 0
1438 1436 return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,n
1439 1437
1440 1438
1441 1439 def pulsePairbyProfiles(self,dataOut):
1442 1440
1443 1441 self.__dataReady = False
1444 1442 data_power = None
1445 1443 data_intensity = None
1446 1444 data_velocity = None
1447 1445 data_specwidth = None
1448 1446 data_snrPP = None
1449 1447 self.putData(data=dataOut.data)
1450 1448 if self.__profIndex == self.n:
1451 1449 data_power,data_intensity, data_velocity,data_snrPP,data_specwidth, n = self.pushData(dataOut=dataOut)
1452 1450 self.__dataReady = True
1453 1451
1454 1452 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth
1455 1453
1456 1454
1457 1455 def pulsePairOp(self, dataOut, datatime= None):
1458 1456
1459 1457 if self.__initime == None:
1460 1458 self.__initime = datatime
1461 1459 data_power, data_intensity, data_velocity, data_snrPP, data_specwidth = self.pulsePairbyProfiles(dataOut)
1462 1460 self.__lastdatatime = datatime
1463 1461
1464 1462 if data_power is None:
1465 1463 return None, None, None,None,None,None
1466 1464
1467 1465 avgdatatime = self.__initime
1468 1466 deltatime = datatime - self.__lastdatatime
1469 1467 self.__initime = datatime
1470 1468
1471 1469 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth, avgdatatime
1472 1470
1473 1471 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1474 1472
1475 1473 if not self.isConfig:
1476 1474 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1477 1475 self.isConfig = True
1478 1476 data_power, data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
1479 1477 dataOut.flagNoData = True
1480 1478
1481 1479 if self.__dataReady:
1482 1480 dataOut.nCohInt *= self.n
1483 1481 dataOut.dataPP_POW = data_intensity # S
1484 1482 dataOut.dataPP_POWER = data_power # P
1485 1483 dataOut.dataPP_DOP = data_velocity
1486 1484 dataOut.dataPP_SNR = data_snrPP
1487 1485 dataOut.dataPP_WIDTH = data_specwidth
1488 1486 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1489 1487 dataOut.utctime = avgdatatime
1490 1488 dataOut.flagNoData = False
1491 1489 return dataOut
1492 1490
1493 1491
1494 1492
1495 1493 # import collections
1496 1494 # from scipy.stats import mode
1497 1495 #
1498 1496 # class Synchronize(Operation):
1499 1497 #
1500 1498 # isConfig = False
1501 1499 # __profIndex = 0
1502 1500 #
1503 1501 # def __init__(self, **kwargs):
1504 1502 #
1505 1503 # Operation.__init__(self, **kwargs)
1506 1504 # # self.isConfig = False
1507 1505 # self.__powBuffer = None
1508 1506 # self.__startIndex = 0
1509 1507 # self.__pulseFound = False
1510 1508 #
1511 1509 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1512 1510 #
1513 1511 # #Read data
1514 1512 #
1515 1513 # powerdB = dataOut.getPower(channel = channel)
1516 1514 # noisedB = dataOut.getNoise(channel = channel)[0]
1517 1515 #
1518 1516 # self.__powBuffer.extend(powerdB.flatten())
1519 1517 #
1520 1518 # dataArray = numpy.array(self.__powBuffer)
1521 1519 #
1522 1520 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1523 1521 #
1524 1522 # maxValue = numpy.nanmax(filteredPower)
1525 1523 #
1526 1524 # if maxValue < noisedB + 10:
1527 1525 # #No se encuentra ningun pulso de transmision
1528 1526 # return None
1529 1527 #
1530 1528 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1531 1529 #
1532 1530 # if len(maxValuesIndex) < 2:
1533 1531 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1534 1532 # return None
1535 1533 #
1536 1534 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1537 1535 #
1538 1536 # #Seleccionar solo valores con un espaciamiento de nSamples
1539 1537 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1540 1538 #
1541 1539 # if len(pulseIndex) < 2:
1542 1540 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1543 1541 # return None
1544 1542 #
1545 1543 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1546 1544 #
1547 1545 # #remover senales que se distancien menos de 10 unidades o muestras
1548 1546 # #(No deberian existir IPP menor a 10 unidades)
1549 1547 #
1550 1548 # realIndex = numpy.where(spacing > 10 )[0]
1551 1549 #
1552 1550 # if len(realIndex) < 2:
1553 1551 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1554 1552 # return None
1555 1553 #
1556 1554 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1557 1555 # realPulseIndex = pulseIndex[realIndex]
1558 1556 #
1559 1557 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1560 1558 #
1561 1559 # print "IPP = %d samples" %period
1562 1560 #
1563 1561 # self.__newNSamples = dataOut.nHeights #int(period)
1564 1562 # self.__startIndex = int(realPulseIndex[0])
1565 1563 #
1566 1564 # return 1
1567 1565 #
1568 1566 #
1569 1567 # def setup(self, nSamples, nChannels, buffer_size = 4):
1570 1568 #
1571 1569 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1572 1570 # maxlen = buffer_size*nSamples)
1573 1571 #
1574 1572 # bufferList = []
1575 1573 #
1576 1574 # for i in range(nChannels):
1577 1575 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1578 1576 # maxlen = buffer_size*nSamples)
1579 1577 #
1580 1578 # bufferList.append(bufferByChannel)
1581 1579 #
1582 1580 # self.__nSamples = nSamples
1583 1581 # self.__nChannels = nChannels
1584 1582 # self.__bufferList = bufferList
1585 1583 #
1586 1584 # def run(self, dataOut, channel = 0):
1587 1585 #
1588 1586 # if not self.isConfig:
1589 1587 # nSamples = dataOut.nHeights
1590 1588 # nChannels = dataOut.nChannels
1591 1589 # self.setup(nSamples, nChannels)
1592 1590 # self.isConfig = True
1593 1591 #
1594 1592 # #Append new data to internal buffer
1595 1593 # for thisChannel in range(self.__nChannels):
1596 1594 # bufferByChannel = self.__bufferList[thisChannel]
1597 1595 # bufferByChannel.extend(dataOut.data[thisChannel])
1598 1596 #
1599 1597 # if self.__pulseFound:
1600 1598 # self.__startIndex -= self.__nSamples
1601 1599 #
1602 1600 # #Finding Tx Pulse
1603 1601 # if not self.__pulseFound:
1604 1602 # indexFound = self.__findTxPulse(dataOut, channel)
1605 1603 #
1606 1604 # if indexFound == None:
1607 1605 # dataOut.flagNoData = True
1608 1606 # return
1609 1607 #
1610 1608 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1611 1609 # self.__pulseFound = True
1612 1610 # self.__startIndex = indexFound
1613 1611 #
1614 1612 # #If pulse was found ...
1615 1613 # for thisChannel in range(self.__nChannels):
1616 1614 # bufferByChannel = self.__bufferList[thisChannel]
1617 1615 # #print self.__startIndex
1618 1616 # x = numpy.array(bufferByChannel)
1619 1617 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1620 1618 #
1621 1619 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1622 1620 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1623 1621 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1624 1622 #
1625 1623 # dataOut.data = self.__arrayBuffer
1626 1624 #
1627 1625 # self.__startIndex += self.__newNSamples
1628 1626 #
1629 1627 # return No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now