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